/**************************************************************************** * * McStas, neutron ray-tracing package * Copyright 1997-2003, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Pol_guide_mirror * * %I * Written by: Erik B Knudsen * Date: July 2018 * Origin: DTU Physics * * Polarising guide with a supermirror along its diagonal. * * %D * Based on Pol_guide_vmirror by P. Christiansen. * Models a rectangular guide with entrance centered on the Z axis and * with one supermirros sitting on the diagonal inside. * The entrance lies in the X-Y plane. Draws a true depiction * of the guide with mirror and trajectories. * The polarisation is handled similar to in Monochromator_pol. * The reflec functions are handled similar to Pol_mirror. * The up direction is hardcoded to be along the y-axis (0, 1, 0) * * Note that this component can also be used as a frame overlap-mirror * if the up and down reflectivities are set equal. In this case the wall * refletivity (rPar) should probably be set to absorb. * * The parameters can either be * double pointer initializations (e.g. {R0, Qc, alpha, m, W}) * or table names (e.g."supermirror_m2.rfl" AND useTables=1). * NB! This might cause warnings by the compiler that can be ignored. * * GRAVITY: YES * * %BUGS * No absorption by mirror. * * %P * INPUT PARAMETERS: * * xwidth: [m] Width at the guide entry * yheight: [m] Height at the guide entry * length: [m] length of guide * rFunc: [1] Guide Reflection function * rPar: [1] Guide Parameters for rFunc * rUpFunc: [1] Mirror Reflection function for spin up * rUpPar: [1] Mirror Parameters for rUpFunc * rDownFunc: [1] Mirror Reflection function for spin down * rDownPar: [1] Mirror Parameters for rDownFunc * useTables: [1] Parameters are 0: Values, 1: Table names * debug: [1] if debug > 0 print out some internal runtime parameters * * OUTPUT PARAMETERS: * * localG: [m/s/s] Gravity vector in guide reference system * normalTop: [1] One of several normal vectors used for defining the geometry * pointTop: [1] One of several points used for defining the geometry * rParPtr: One of several pointers to reflection parameters used with the ref. functions. [] * SCATTERED: [] is 1 for reflected, and 2 for transmitted neutrons * * * %L * * %E *******************************************************************************/ DEFINE COMPONENT Pol_guide_mirror DEFINITION PARAMETERS (xwidth, yheight, length, rFunc=StdReflecFunc, rUpFunc=StdReflecFunc, rDownFunc=StdReflecFunc, rPar ={1.0, 0.0219, 4.07, 3.2, 0.003}, rUpPar ={1.0, 0.0219, 4.07, 3.2, 0.003}, rDownPar={0.1, 0.0219, 4.07, 3.2, 0.003}, useTables=0) SETTING PARAMETERS (int debug=0) OUTPUT PARAMETERS (localG, normalTop, normalBot, normalLeft, normalRight, normalInOut, pointTop, pointBot, pointLeft, pointRight, pointIn, pointOut, rParPtr, rUpParPtr, rDownParPtr) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "pol-lib" %include "ref-lib" %} DECLARE %{ Coords localG; Coords normalTop, normalBot, normalLeft, normalRight, normalInOut; Coords pointTop, pointBot, pointLeft, pointRight, pointIn, pointOut; #if (useTables) t_Table *rParPtr = 0; t_Table *rUpParPtr = 0; t_Table *rDownParPtr = 0; #else double rParPtr[] = (double [])rPar; double rUpParPtr[] = (double [])rUpPar; double rDownParPtr[] = (double [])rDownPar; #endif %} INITIALIZE %{ #if (useTables) rParPtr = (t_Table*) malloc(sizeof(t_Table)); rUpParPtr = (t_Table*) malloc(sizeof(t_Table)); rDownParPtr = (t_Table*) malloc(sizeof(t_Table)); if (Table_Read(rParPtr, rPar, 1) <= 0) { fprintf(stderr,"Pol_guide_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rPar); exit(1); } if (Table_Read(rUpParPtr, rUpPar, 1) <= 0) { fprintf(stderr,"Pol_guide_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rUpPar); exit(1); } if (Table_Read(rDownParPtr, rDownPar, 1) <= 0) { fprintf(stderr,"Pol_guide_mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rDownPar); exit(1); } #endif if ((xwidth<=0) || (yheight<= 0) || (length<=0)) { fprintf(stderr, "Pol_guide_mirror: %s: NULL or negative length scale!\n" "ERROR (xwidth,yheight,length). Exiting\n", NAME_CURRENT_COMP); exit(1); } if (mcgravitation) { localG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0)); fprintf(stdout,"Pol_guide_mirror: %s: Gravity is on!\n", NAME_CURRENT_COMP); } else localG = coords_set(0, 0, 0); // To be able to handle the situation properly where a component of // the gravity is along the z-axis we also define entrance (in) and // exit (out) planes // The entrance and exit plane are defined by the normal vector // (0, 0, 1) // and the two points pointIn=(0, 0, 0) and pointOut=(0, 0, length) normalInOut = coords_set(0, 0, 1); pointIn = coords_set(0, 0, 0); pointOut = coords_set(0, 0, length); // Top plane (+y dir) can be spanned by (1, 0, 0) & (0, 0, 1) // and the point (0, yheight/2, 0) // A normal vector is (0, 1, 0) normalTop = coords_set(0, 1, 0); pointTop = coords_set(0, yheight/2, 0); // Bottom plane (-y dir) can be spanned by (1, 0, 0) & (0, 0, 1) // and the point (0, -yheight/2, 0) // A normal vector is (0, 1, 0) normalBot = coords_set(0, 1, 0); pointBot = coords_set(0, -yheight/2, 0); // Left plane (+x dir) can be spanned by (0, 1, 0) & (0, 0, 1) // and the point (xwidth/2, 0, 0) // A normal vector is (1, 0, 0) normalLeft = coords_set(1, 0, 0); pointLeft = coords_set(xwidth/2, 0, 0); // Right plane (-x dir) can be spanned by (0, 1, 0) & (0, 0, 1) // and the point (-xwidth/2, 0, 0) // A normal vector is (1, 0, 0) normalRight = coords_set(1, 0, 0); pointRight = coords_set(-xwidth/2, 0, 0); %} TRACE %{ #if (useTables) #undef rFunc #undef rUpFunc #undef rDownFunc #define rFunc TableReflecFunc #define rUpFunc TableReflecFunc #define rDownFunc TableReflecFunc #endif /* time threshold */ const double tThreshold = 1e-10/sqrt(vx*vx + vy*vy + vz*vz); const double xwhalf = xwidth/2; const double norm = 1.0/sqrt(xwidth*xwidth + length*length); double R; Coords normalMirror, pointMirror; Coords* normalPointer = 0; // Pol variables double FN, FM, Rup, Rdown, refWeight; /* Propagate neutron to guide entrance. */ PROP_Z0; if (!inside_rectangle(x, y, xwidth, yheight)) ABSORB; SCATTER; normalMirror = coords_set(norm*length, 0, -norm*xwidth); pointMirror = coords_set(-xwhalf, 0, 0); for(;;) { double tLeft, tRight, tTop, tBot, tIn, tOut, tMirror; double tUp, tSide, time, endtime; double Q; //, dummy1, dummy2, dummy3; Coords vVec, xVec; int mirrorReflect; mirrorReflect = 0; xVec = coords_set(x, y, z); vVec = coords_set(vx, vy, vz); solve_2nd_order(&tTop, NULL, 0.5*coords_sp(normalTop,localG), coords_sp(normalTop, vVec), coords_sp(normalTop, coords_sub(xVec, pointTop))); solve_2nd_order(&tBot, NULL, 0.5*coords_sp(normalBot,localG), coords_sp(normalBot, vVec), coords_sp(normalBot, coords_sub(xVec, pointBot))); solve_2nd_order(&tRight, NULL, 0.5*coords_sp(normalRight,localG), coords_sp(normalRight, vVec), coords_sp(normalRight, coords_sub(xVec, pointRight))); solve_2nd_order(&tLeft, NULL, 0.5*coords_sp(normalLeft,localG), coords_sp(normalLeft, vVec), coords_sp(normalLeft, coords_sub(xVec, pointLeft))); solve_2nd_order(&tIn, NULL, 0.5*coords_sp(normalInOut,localG), coords_sp(normalInOut, vVec), coords_sp(normalInOut, coords_sub(xVec, pointIn))); solve_2nd_order(&tOut, NULL, 0.5*coords_sp(normalInOut,localG), coords_sp(normalInOut, vVec), coords_sp(normalInOut, coords_sub(xVec, pointOut))); solve_2nd_order(&tMirror, NULL, 0.5*coords_sp(normalMirror,localG), coords_sp(normalMirror, vVec), coords_sp(normalMirror, coords_sub(xVec, pointMirror))); double nx,ny,nz,px,py,pz; coords_get(normalMirror,&nx,&ny,&nz); coords_get(pointMirror,&px,&py,&pz); plane_intersect(&tMirror, x,y,z,vx,vy,vz, nx, ny, nz, px, py, pz); /* Choose appropriate reflection time */ if (tTop>tThreshold && (tToptThreshold && (tLefttThreshold && (tUptThreshold && tMirrortThreshold && (tOut endtime) break; if(time <= tThreshold) { printf("Time below threshold!\n"); fprintf(stdout, "Pol_guide_mirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n" "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP, tTop, tBot, tRight, tLeft, tUp, tSide, time); break; } if(debug>0 && time==tLeft) { fprintf(stdout, "\nPol_guide_mirror: %s: Left side hit: x, v, normal, point, gravity\n", NAME_CURRENT_COMP); coords_print(xVec); coords_print(vVec); coords_print(normalLeft); coords_print(pointLeft); coords_print(localG); fprintf(stdout, "\nA: %f, B: %f, C: %f, tLeft: %f\n", 0.5*coords_sp(normalLeft,localG),coords_sp(normalLeft, vVec), coords_sp(normalLeft, coords_sub(xVec, pointLeft)), tLeft); } if(debug>0) fprintf(stdout, "Pol_guide_mirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n" "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP, tTop, tBot, tRight, tLeft, tUp, tSide, time); if(debug>0) fprintf(stdout, "Pol_guide_mirror: %s: Start v: (%f, %f, %f)\n", NAME_CURRENT_COMP, vx, vy, vz); PROP_DT(time); if (mcgravitation) vVec = coords_set(vx, vy, vz); //SCATTER; if(time==tTop) normalPointer = &normalTop; else if(time==tBot) normalPointer = &normalBot; else if(time==tRight) normalPointer = &normalRight; else if(time==tLeft) normalPointer = &normalLeft; else if(time==tMirror) normalPointer = &normalMirror; else fprintf(stderr, "Pol_guide_mirror: %s: This should never happen!!!!\n", NAME_CURRENT_COMP); Q = 2*coords_sp(vVec, *normalPointer)*V2K; if(!mirrorReflect) { // we have hit one of the sides. Always reflect. vVec = coords_add(vVec, coords_scale(*normalPointer, -Q*K2V)); rFunc(fabs(Q), rParPtr, &refWeight); p *= refWeight; SCATTER; } else { // we have hit the mirror rUpFunc(fabs(Q), rUpParPtr, &Rup); rDownFunc(fabs(Q), rDownParPtr, &Rdown); if (Rup < 0) ABSORB; if (Rup > 1) Rup =1 ; if (Rdown < 0) ABSORB; if (Rdown > 1) Rdown =1 ; GetMonoPolFNFM(Rup, Rdown, &FN, &FM); GetMonoPolRefProb(FN, FM, sy, &refWeight); // check that refWeight is meaningful if (refWeight < 0) ABSORB; if (refWeight > 1) refWeight =1 ; if ( rand01()1.000001) { // check that polarisation is meaningfull fprintf(stderr, "Pol_guide_mirror: %s: polarisation |s|=%g > 1 s=[%g,%g,%g]\n", NAME_CURRENT_COMP, sx*sx+sy*sy+sz*sz, sx, sy, sz); } } if(p==0) { ABSORB; break; } // set new velocity vector coords_get(vVec, &vx, &vy, &vz); if(debug>0){ fprintf(stdout, "Pol_guide_mirror: %s: End v: (%f, %f, %f)\n", NAME_CURRENT_COMP, vx, vy, vz); } } #if (useTables) #undef rFunc #undef rUpFunc #undef rDownFunc #define rFunc StdReflecFunc #define rUpFunc StdReflecFunc #define rDownFunc StdReflecFunc #endif %} MCDISPLAY %{ int i, j; // draw box box(0, 0, length/2.0, xwidth, yheight, length); for(j = -1; j<=1; j+=2) line(-xwidth/2.0, j*yheight/2, 0, xwidth/2.0, j*yheight/2, length); %} END