/**************************************************************************** * * 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 supermirror 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 * reflectivity (rPar) should probably be set to 0. * * Reflectivity values can either come from datafiles or from * sets of parameters to the standard analytic reflectivity function commonly * used for neutron guides: * R=R0; q 0 print out some internal runtime parameters * * %L * * %E *******************************************************************************/ DEFINE COMPONENT Pol_guide_mirror DEFINITION PARAMETERS () SETTING PARAMETERS ( vector rUpPar={1.0, 0.0219, 4.07, 3.2, 0.003}, vector rDownPar={1.0, 0.0219, 4.07, 3.2, 0.003}, vector rPar={1.0, 0.0219, 4.07, 3.2, 0.003}, string rData="",string rUpData="",string rDownData="", xwidth, yheight, length, int debug=0) OUTPUT PARAMETERS () /* 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; Coords normalBot; Coords normalLeft; Coords normalRight; Coords normalInOut; Coords pointTop; Coords pointBot; Coords pointLeft; Coords pointRight; Coords pointIn; Coords pointOut; t_Table rTable; t_Table rUpTable; t_Table rDownTable; int rTableFlag; int rUpTableFlag; int rDownTableFlag; %} INITIALIZE %{ if (strlen(rData) && strcmp(rData,"NULL")){ if (Table_Read(&rTable, rData, 1) <= 0) { fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rPar); exit(1); } rTableFlag=1; }else{ rTableFlag=0; } if (strlen(rUpData) && strcmp(rUpData,"NULL")){ if (Table_Read(&rUpTable, rUpData, 1) <= 0) { fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rUpPar); exit(1); } rUpTableFlag=1; }else { rUpTableFlag=0; } if (strlen(rDownData) && strcmp(rDownData,"NULL")){ if (Table_Read(&rDownTable, rDownData, 1) <= 0) { fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n", NAME_CURRENT_COMP, rDownPar); exit(1); } rDownTableFlag=1; }else{ rDownTableFlag=0; } 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 %{ /* 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); 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)); if(rTableFlag){ refWeight=Table_Value(rTable, fabs(Q), 1); }else{ StdReflecFunc(fabs(Q), rPar, &refWeight); } p *= refWeight; SCATTER; } else { /* we have hit the mirror */ if(rUpTableFlag){ Rup=Table_Value(rUpTable,fabs(Q),1); }else{ StdReflecFunc(fabs(Q), rUpPar, &Rup); } if(rDownTableFlag){ Rdown=Table_Value(rDownTable,fabs(Q),1); }else{ StdReflecFunc(fabs(Q), rDownPar, &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); } } %} 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