/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2008, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: FermiChopper. * * %Identification * * Written by: M. Poehlmann, C. Carbogno, H. Schober, E. Farhi * Date: May 2002 * Origin: ILL Grenoble / TU Muenchen * Modified by: K.Lieutenant, June 2005: added phase parameter. Comp validation. * Modified by: EF, Nov 2005: completely rewrote comp. * Modified by: EF, Oct 2008: fix chopper orientation * Modified by: EF, Mar 2009: improved intersection algorithm * * Fermi Chopper with rotating frame. * * %D * Models a fermi chopper with optional supermirror coated blades * supermirror facilities may be disabled by setting m = 0, R0=0 * Slit packages are straight. Chopper slits are separated by an infinitely * thin absorbing material. The effective transmission (resulting from fraction * of the transparent material and its transmission) may be specified. * The chopper slit package width may be specified through the total width 'xwidth' * of the full package or the width 'w' of each single slit. The other parameter * is calculated by: xwidth = nslit*w. The slit package may be made curved and use * super-mirror coating. * * Example: * FermiChopper(phase=-50.0, radius=0.04, nu=100, yheight=0.08, w=0.00022475, nslit=200.0, R0=0.0, Qc=0.02176, alpha=2.33, m=0.0, length=0.012, eff=0.95) * * %VALIDATION * Apr 2005: extensive external test, most problems solved (cf. 'Bugs') * Validated by: K. Lieutenant, E. Farhi * * limitations: * no absorbing blade width used * * %Parameters * INPUT PARAMETERS: * * Geometrical chopper constants: * radius: [m] chopper cylinder radius * nslit: [1] number of chopper slits * length: [m] channel length of the Fermi chopper * w: [m] width of one chopper slit * xwidth: [m] optional total width of slit package * yheight: [m] height of slit package * nu: [Hz] chopper frequency. Omega=2*PI*nu in rad/s, nu*60 in rpm. Positive value corresponds to a counter-clockwise rotation around y. * eff: [1] efficiency = transmission x fraction of transparent material * verbose: [1] set to 1,2 or 3 gives debugging information * curvature: [m-1] Curvature of slits (1/radius of curvature). * * Supermirror constants: * m: [1] m-value of material. Zero means completely absorbing. * alpha: [AA] slope of reflectivity * Qc: [AA-1] critical scattering vector * W: [AA-1] width of supermirror cut-off * R0: [1] low-angle reflectivity * * Constants to reset time of flight: * zero_time: [1] set time to zero: 0=no, 1=once per half cycle, 2=auto adjust phase * phase: [deg] chopper phase at t=0 * delay: [s] sets phase so that transmision is centered on 'delay' * * OUTPUT PARAMETERS: * FCVars : [] structure * * %L * Vitess_ChopperFermi component by * G. Zsigmond, imported from Vitess by K. Lieutenant. * * %End *****************************************************************************/ /* NOTE: * The initial component version (McStas version <= 1.12) was written * with an inverted coordinate frame orientation. This corresponds * to inverting the frequency and phase sign. */ DEFINE COMPONENT FermiChopper DEFINITION PARAMETERS () SETTING PARAMETERS (phase=0, radius=0.04, nu=100, w=0.00022475, nslit=200, R0=0.0, Qc=0.02176, alpha=2.33, m=0.0, W=2e-3, length=0.012, eff=0.95, zero_time=0, xwidth=0, verbose=0, yheight=0.08, curvature=0,delay=0) OUTPUT PARAMETERS() SHARE %{ %include "ref-lib" #ifndef FermiChopper_TimeAccuracy #define FermiChopper_TimeAccuracy 1e-9 #define FermiChopper_MAXITER 100 /* Definition of internal variable structure: all counters */ struct FermiChopper_struct { double omega; /* chopper rotation */ double ph0; /* chopper rotation */ double t0; /* chopper rotation */ double C_slit; /* slit curvature radius in [m] */ double L_slit; /* slit package length [m] */ double sum_t; double sum_v; double sum_N; double sum_N_pass; /* events */ long absorb_alreadyinside; long absorb_topbottom; long absorb_cylentrance; long absorb_sideentrance; long absorb_notreachentrance; long absorb_packentrance; long absorb_slitcoating; long warn_notreachslitwall; long absorb_exitslitpack; long absorb_maxiterations; long absorb_wrongdirection; long absorb_nocontrol; long absorb_cylexit; long warn_notreachslitoutput; char compcurname[256]; }; /***************************************************************************** * FC_zrot: returns Z' in rotating frame, from X,Z and t,omega,ph0 ****************************************************************************/ #pragma acc routine seq double FC_zrot(double X, double Z, double T, struct FermiChopper_struct FCs) { double omega =FCs.omega; double ph0 =FCs.ph0; return( Z*cos(omega*T+ph0)-X*sin(omega*T+ph0) ); } /***************************************************************************** * FC_xrot: returns X' in rotating frame, from X,Z and omega,t,ph0 * additional coordinate shift in case of curved slits ****************************************************************************/ #pragma acc routine seq double FC_xrot(double X, double Z, double T, struct FermiChopper_struct FCs) { double omega =FCs.omega; double ph0 =FCs.ph0; double C_slit=FCs.C_slit; double ret, tmp; ret = X*cos(omega*T+ph0)+Z*sin(omega*T+ph0); if (C_slit) { tmp = fabs(FC_zrot(X, Z, T, FCs)); if (tmp < FCs.L_slit/2) { tmp = (FCs.L_slit/2 - tmp)*C_slit; ret += (1-sqrt(1-tmp*tmp))/C_slit; } } return( ret ); } /***************************************************************************** * FC_xzrot_dt(x,z,vx,vz, t,dt, type='x' or 'z', FCs) * returns X' or Z' in rotating frame, from X,Z and t,omega,ph0 * taking into account propagation with velocity during time dt ****************************************************************************/ #pragma acc routine seq double FC_xzrot_dt(double x, double z, double vx, double vz, double t, double dt, char type, struct FermiChopper_struct FCs) { if (dt) /* with propagation */ return( (type == 'x' ? FC_xrot(x+vx*dt, z+vz*dt, t+dt, FCs) : FC_zrot(x+vx*dt, z+vz*dt, t+dt, FCs)) ); else /* without propagation */ return( (type == 'x' ? FC_xrot(x,z,t,FCs) : FC_zrot(x,z,t,FCs)) ); } /***************************************************************************** * FC_xzbrent(x,z,vx,vz, t,dt, type='x' or 'z', d, FCs) * solves X'=d and Z'=d with Brent algorithm in time interval [0, dt]. * Returns time within [0,dt], from NumRecip in C, chap 9, p360 (zbrent) * ERRORS: return -1 not used * -2 if exceed MAX iteration * -3 no sign change in range ****************************************************************************/ #pragma acc routine seq double FC_xzbrent(double x, double z, double vx, double vz, double t, double dt, char type, double d, struct FermiChopper_struct FCs) { int iter; double a=0,b=dt; double c,dd,e,min1,min2; double tol=FermiChopper_TimeAccuracy; double EPS=FermiChopper_TimeAccuracy; double fa=FC_xzrot_dt(x,z,vx,vz, t,a, type, FCs) - d; double fb=FC_xzrot_dt(x,z,vx,vz, t,b, type, FCs) - d; double fc,p,q,r,s,tol1,xm; if (fb*fa > 0.0) return -3; fc=fb; for (iter=1;iter<=FermiChopper_MAXITER;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=dd=b-a; } if (fabs(fc) < fabs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*fabs(b)+0.5*tol; xm=0.5*(c-b); if (fabs(xm) <= tol1 || fb == 0.0) return b; if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=fabs(p); min1=3.0*xm*q-fabs(tol1*q); min2=fabs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=dd; dd=p/q; } else { dd=xm; e=dd; } } else { dd=xm; e=dd; } a=b; fa=fb; if (fabs(dd) > tol1) b += dd; else b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1)); fb=FC_xzrot_dt(x,z,vx,vz, t,b, type, FCs) - d; } return -2; } /* FC_xzbrent */ /***************************************************************************** * Wrappers to intersection algorithms ****************************************************************************/ #pragma acc routine seq double FC_xintersect(double x, double z, double vx, double vz, double t, double dt, double d, struct FermiChopper_struct FCs) { return(FC_xzbrent(x, z, vx, vz, t, dt, 'x', d, FCs)); } #pragma acc routine seq double FC_zintersect(double x, double z, double vx, double vz, double t, double dt, double d, struct FermiChopper_struct FCs) { return(FC_xzbrent(x, z, vx, vz, t, dt, 'z', d, FCs)); } #endif %} DECLARE %{ struct FermiChopper_struct FCVars; %} INITIALIZE %{ /************************ CALCULATION CONSTANTS *****************************/ strcpy(FCVars.compcurname, NAME_CURRENT_COMP); FCVars.omega = 2*PI*nu; if (!phase && delay) { FCVars.ph0= fmod(-delay*nu*360,360)*DEG2RAD; } else FCVars.ph0 = phase*DEG2RAD; FCVars.sum_t=FCVars.sum_v=FCVars.sum_N=FCVars.sum_N_pass=0; /* check of input parameters */ if (nslit < 1) nslit=1; if (yheight <= 0) exit(printf("FermiChopper: %s: FATAL: unrealistic cylinder yheight =%g [m]\n", NAME_CURRENT_COMP, yheight)); if (m <= 0) { m=0; R0=0; } if (radius <= 0) { printf("FermiChopper: %s: FATAL: Unrealistic cylinder radius radius=%g [m]\n", NAME_CURRENT_COMP, radius); exit(-1); } if (xwidth > 0 && xwidth < radius*2 && nslit > 0) { w = xwidth/nslit; } if (w <= 0) { printf("FermiChopper: %s: FATAL: Slits in the package have unrealistic width w=%g [m]\n", NAME_CURRENT_COMP, w); exit(-1); } if (nslit*w > radius*2) { nslit = floor(radius/w); printf("FermiChopper: %s: Too many slits to fit in the cylinder\n" " Adjusting nslit=%f\n", NAME_CURRENT_COMP, nslit); } if (length > radius*2) { length = 2*sqrt(radius*radius - nslit*w*nslit*w/4); printf("FermiChopper: %s: Slit package is longer than the whole\n" " chopper cylinder. Adjusting length=%g [m]\n", NAME_CURRENT_COMP, length); } if (eff <= 0 || eff > 1) { eff = 0.95; printf("FermiChopper: %s: Efficiency is unrealistic\n" " Adjusting eff=%f\n", NAME_CURRENT_COMP, eff); } if (Qc <= 0) { Qc = 0.02176; m = 0; R0 = 0; } if (W <= 0) W=1e-6; if (curvature) { FCVars.C_slit = curvature; if (1 < fabs(radius*curvature)) exit(printf("FermiChopper: %s: Slit curvature is unrealistic\n", NAME_CURRENT_COMP)); } FCVars.L_slit = length; if (verbose && nu) printf("FermiChopper: %s: Frequency nu=%g [Hz] %g [rpm], time frame=%g [s] phase=%g [deg]\n" , NAME_CURRENT_COMP, nu, nu*60, 2/nu, FCVars.ph0*RAD2DEG); FCVars.absorb_alreadyinside = 0; FCVars.absorb_topbottom = 0; FCVars.absorb_cylentrance = 0; FCVars.absorb_sideentrance = 0; FCVars.absorb_notreachentrance = 0; FCVars.absorb_packentrance = 0; FCVars.absorb_slitcoating = 0; FCVars.warn_notreachslitwall = 0; FCVars.absorb_exitslitpack = 0; FCVars.absorb_maxiterations = 0; FCVars.absorb_wrongdirection = 0; FCVars.absorb_nocontrol = 0; FCVars.absorb_cylexit = 0; FCVars.warn_notreachslitoutput = 0; /* fix for the wrong coordinate frame orientation to come back to McStas XYZ system */ FCVars.omega *= -1; FCVars.ph0 *= -1; FCVars.t0 = -FCVars.ph0/FCVars.omega; %} TRACE %{ /** local CALCULATION VARIABLES**************************************/ /** Interaction with slit package ***************************/ double slit_input; /* length of the slits */ /** Variables for calculating interaction with blades ***************/ double xp1, zp1, vxp1; /** Reflections ***********************************************/ double n1; /** Multiple Reflections ******************************/ int loopcounter=0; /* How many reflections happen? */ /** Time variables *********************************/ double t3=0; /* interaction time at n3 position (slit wall) */ double t1=0,t2=0; /* cylinder intersection time (at entry and exit of slit pack - n1 n2 - or cylinder) */ double dt; /* interaction intervals (e.g. time for exiting slit pack) */ double X[5],Z[5]; /* position of interaction locations in the Fermi chopper rotating frame [0]=cylinder input [1]=slit input [2]=slit wall reflection [3]=slit exit [4]=cylinder exit */ double ref = 0; double par[5] = {R0, Qc, alpha, m, W}; /************************ TIME OF FLIGHT RESET ************************/ if (zero_time == 1 && nu) t -= floor( (t+1/(4*nu))*(2*nu) )/(2*nu); /* zero arrays used to store positions */ for (loopcounter=0; loopcounter < 5; loopcounter++) X[loopcounter] = Z[loopcounter]=0; /************** test, if the neutron interacts with the cylinder ***/ if (cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight)) { if (t1 <= 0) { /* Neutron started inside the cylinder */ if (verbose > 0 && FCVars.absorb_alreadyinside 2) printf("FermiChopper: %s: t1=%8.3g t2=%8.3g xyz=[%8.3g %8.3g %8.3g] v=[%8.3g %8.3g %8.3g] t=%8.3g (init).\n", NAME_CURRENT_COMP, t1, t2, x,y,z,vx,vy,vz,t); dt=t2-t1; /* total time of flight inside the cylinder */ PROP_DT(t1); /* Propagates neutron to entrance of the cylinder */ SCATTER; if (verbose > 2) printf("FermiChopper: %s: PROP_DT t1=%8.3g t2=%8.3g xyz=[%8.3g %8.3g %8.3g] v=[%8.3g %8.3g %8.3g] t=%8.3g (IN cyl).\n", NAME_CURRENT_COMP, t1, t2, x,y,z,vx,vy,vz,t); /* neutron must not enter or leave from top or bottom of cylinder. */ if (fabs(y) >= yheight/2.0 || fabs(y+vy*dt) >= yheight/2.0) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (enter).\n", NAME_CURRENT_COMP, y); #pragma acc atomic FCVars.absorb_topbottom = FCVars.absorb_topbottom + 1; ABSORB; } vxp1 = sqrt(vx*vx+vy*vy+vz*vz); double tmp=p*vxp1; #pragma acc atomic FCVars.sum_v = FCVars.sum_v + tmp; tmp = p*t; #pragma acc atomic FCVars.sum_t = FCVars.sum_t + tmp; #pragma acc atomic FCVars.sum_N = FCVars.sum_N + p; if (zero_time > 1 && FCVars.sum_N) { /* automatic phase adjustment */ double mean_t, mean_phase; mean_t = FCVars.sum_t/FCVars.sum_N; mean_phase = fmod(mean_t*nu*2*PI, 2*PI); /* now we shift the phase so that the neutron pulse is centered on the slit pack */ mean_phase+= radius/vxp1*2*PI*nu; #pragma acc atomic write FCVars.ph0 = mean_phase; } /* neutron must enter the cylinder opening: |X'| < full package width*/ xp1 = FC_xrot(x,z, t,FCVars); /* X'(t) */ if (fabs(xp1) >= nslit*w/2) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside cylinder aperture, x'=%8.3g > %g (enter).\n", NAME_CURRENT_COMP, xp1, nslit*w/2); #pragma acc atomic FCVars.absorb_cylentrance = FCVars.absorb_cylentrance + 1; ABSORB; } /*********************** PROPAGATE TO SLIT PACKAGE **************************/ /* zp1 = Z' at entrance of cylinder Z'(t) */ zp1 = FC_zrot(x,z, t, FCVars); X[0] = xp1; Z[0] = zp1; /* here we should have sqrt(x*x+z*z) == sqrt(xp1*xp1+zp1*zp1) */ t3 = sqrt(x*x+z*z) / sqrt(xp1*xp1+zp1*zp1); if (t3 < 0.99 || 1.01 < t3) { if (verbose > 1 && FCVars.absorb_cylentrance < FermiChopper_MAXITER) printf("FermiChopper: %s: ABSORB Neutron radius on cylinder in static frame r=%g does not match that of rotating frame r'=%g.\n", NAME_CURRENT_COMP, sqrt(x*x+z*z), sqrt(xp1*xp1+zp1*zp1)); #pragma acc atomic FCVars.absorb_cylentrance = FCVars.absorb_cylentrance +1 ; ABSORB; } /* Checking on which side of the Chopper the Neutron enters: sign(Z') */ slit_input = (zp1 > 0 ? length/2 : -length/2); /* time shift to reach slit package in [0,time to exit cylinder]: Z'=slit_input */ /* t3 is used here as a tmp variable, will be redefined in for loop */ t3 = FC_zintersect(x,z,vx,vz, t,dt, slit_input, FCVars); if( (t3 < 0)||(t3 > dt) ) { if (verbose > 2 && FCVars.absorb_notreachentrance < FermiChopper_MAXITER) { printf("FermiChopper: %s: Can not reach entrance of slits. dt=%8.3g t3=%8.3g (intersection:1).\n", NAME_CURRENT_COMP, dt, t3); if (t3 == -3) printf(" No sign change to determine intersection\n"); else if (t3 == -2) printf(" Max iterations reached\n"); else if (t3 < 0) printf(" Error when solving intersection\n"); } #pragma acc atomic FCVars.absorb_notreachentrance = FCVars.absorb_notreachentrance + 1; ABSORB; /* neutron can not reach slit entrance */ } /* Propagating to the slit package entrance */ PROP_DT(t3); /* dt = t2-t1: time in cylinder */ dt -= t3; /* remaining time from slit pack entry to exit of cylinder */ xp1 = FC_xrot(x,z, t, FCVars); /* should be slit_input */ zp1 = FC_zrot(x,z, t, FCVars); X[1] = xp1; Z[1] = zp1; #ifndef OPENACC if (mcdotrace) { /* indicate position of neutron in mcdisplay */ double xp2 = x; double zp2 = z; x = xp1; z=zp1; SCATTER; x=xp2; z=zp2; } else #endif SCATTER; if (verbose > 2) printf("FermiChopper: %s: PROP_DT t=%8.3g dt=%8.3g x'=%8.3g z'=%8.3g length=%g (slit enter).\n", NAME_CURRENT_COMP, t, dt, xp1, zp1, slit_input); /* must have X'< slit package width at package Z'=slit_input */ if (fabs(xp1) >= nslit*w/2) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package, x'=%8.3g > %g (enter).\n", NAME_CURRENT_COMP, xp1, nslit*w/2); #pragma acc atomic FCVars.absorb_packentrance = FCVars.absorb_packentrance + 1; ABSORB; } /* solve Z'=-slit_input for time of exit of slit package */ /* t3 is used here as a tmp variable, will be redefined in for loop */ t3 = FC_zintersect(x,z,vx,vz, t,dt*1.1, -slit_input, FCVars); if((t3 < FermiChopper_TimeAccuracy)||(t3 > dt)) { if (verbose > 1 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) { printf("FermiChopper: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g (intersection:2).\n", NAME_CURRENT_COMP, dt, t3); if (t3 == -3) printf(" No sign change to determine intersection\n"); else if (t3 == -2) printf(" Max iterations reached\n"); else if (t3 < 0) printf(" Error when solving intersection\n"); } #pragma acc atomic FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1; /* we estimate analytically the time to the slit exit */ t3 = length/vxp1; } dt = t3; /* reduce time interval to [0, time of slit exit] */ /* here we should have dt*v = length (slit entrance -> exit) */ t3 = fabs(FC_xzrot_dt(x,z,vx,vz, t, 0 , 'z', FCVars) - FC_xzrot_dt(x,z,vx,vz, t, dt, 'z', FCVars))/length; if (fabs(t3-1) > 0.02) { if (verbose > 0 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) printf("FermiChopper: %s: ABSORB Neutron propagation time v*dt/length=%g in slit does not match its length=%g (slit exit expected).\n", NAME_CURRENT_COMP, t3, length); #pragma acc atomic FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1; ABSORB; } /* dt= time shift to go to exit of slit package (or exit of cylinder in case of error) */ /* |---------| | / | |o * (dt) | | |---------| */ /*********************PROPAGATION INSIDE THE SLIT PACKAGE *******************/ /* index of slit hit at entrance n1=-N/2 (xp1=-) ... N/2 (xp1=+) */ n1 = floor(xp1/w); /******************* BEGIN LOOP FOR MULTIPLE REFLECTIONS ********************/ for (loopcounter=0; loopcounter<=FermiChopper_MAXITER; loopcounter++) { double dt_to_tangent=0; /* time shift to go to tangent intersection */ double n2,n3; /* slit indices */ double xp2, zp2,xp3=0,vxp2=0,vzp1=0; /* position, velocity */ double q; /* used for calculating new velocities after reflection */ int i; /* compute trajectory tangents: m1=Vz'+w.X'(t), m2=Vz'+w.X'(t+dt) */ xp1 = FC_xrot (x,z, t, FCVars); /* X'(t) current position */ xp2 = FC_xzrot_dt(x,z,vx,vz, t, dt, 'x', FCVars); /* X'(t+dt) slit exit */ zp2 = FC_xzrot_dt(x,z,vx,vz, t, dt, 'z', FCVars); /* Z'(t+dt) slit exit */ /* slit index at the end of the slit: */ n2 = floor(xp2/w); /* quick exit for absorbing walls when neutron changes slit index */ if (n2 != n1 && (m <= 0 || R0 <= 0 || Qc <= 0)) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (change slit).\n", NAME_CURRENT_COMP); #pragma acc atomic FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1; ABSORB; } /* compute transversal velocity to determine their intersection */ vxp1= FC_xrot (vx+z*FCVars.omega,vz-x*FCVars.omega, t, FCVars); /* dX'(t)/dt slope at current position*/ vxp2= FC_xrot (vx+(z+vz*dt)*FCVars.omega,vz-(x+vx*dt)*FCVars.omega, t+dt,FCVars); /* dX'(t+dt)/dt slope at slit exit */ /* absolute time at tangent intersection, changed to time shift below */ dt_to_tangent = (vxp1 - vxp2 ? (xp2 - xp1 - dt*vxp2)/(vxp1 - vxp2) : -1); /* If algorithm fails, take the middle of the interval*/ if (dt_to_tangent < 0 || dt_to_tangent > dt) { if (verbose > 1) printf("FermiChopper: %s: WARNING could not determine tangent intersection dt=%g. Using middle.\n", NAME_CURRENT_COMP, dt_to_tangent); dt_to_tangent=dt*0.5; } /* *(dt_to_tangent, xp3) |---------| | / \ | xp1 |o \|(dt) xp2 | | |---------| */ /* point coordinates at tangent intersection/middle point (max deviation from optical axis) */ xp3 = FC_xzrot_dt(x,z,vx,vz, t, dt_to_tangent, 'x', FCVars); /* X'(t+dt_to_tangent) */ /* slit index at the tangent intersection/middle point */ n3 = floor(xp3/w); if (verbose > 2) printf("FermiChopper: %s: t3=%8.3g slit_indices=[%g %g %g] (time at tangent intersection).\n", NAME_CURRENT_COMP, dt_to_tangent, n1, n2, n3); /* change slit means there is a reflection/intersection inside */ if ( n2!=n1 || n3!=n1 ) { double t3a, t3b, distance_Wa, distance_Wb; if (m <= 0 || R0 <= 0 || Qc <= 0) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (change slit).\n", NAME_CURRENT_COMP); #pragma acc atomic FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1; ABSORB; } /* Choosing the first time it isn't in the slit anymore */ if(n2 == n1 && n3 != n1){ n2 = n3; } /* get position of slit wall towards which neutron is propagating */ if (n2 > n1) { /* X' positive side of slit is in principle the first intersection to test*/ distance_Wa = n1*w+w; distance_Wb = n1*w; } else { /* X' negative side of slit */ distance_Wb = n1*w+w; distance_Wa = n1*w; } /* time shift to reach slit wall point in [0,dt_to_tangent]: X'=distance_W slit_wall */ for (i=0; i< 2; i++) { /* first attempt: [0,dt_to_tangent] (to max deviation location) second attempt: [0, dt] (to slit exit) */ double dt_search = (i == 0 ? dt_to_tangent : dt); t3a = FC_xintersect(x,z,vx,vz, t,dt_to_tangent, distance_Wa, FCVars); t3b = FC_xintersect(x,z,vx,vz, t,dt_to_tangent, distance_Wb, FCVars); if (t3b < 0) t3 = t3a; else if (t3a < 0 && t3b >= 0) t3 = t3b; else t3 = (t3a < t3b ? t3a : t3b); if (FermiChopper_TimeAccuracy < t3 && t3 < dt_search) break; /* we found the intersection location */ } /* handle case where intersection search fails */ if (t3 < FermiChopper_TimeAccuracy || t3 >= dt_to_tangent) { if (verbose > 1 && FCVars.warn_notreachslitwall < FermiChopper_MAXITER) { printf("FermiChopper: %s: Can not reach slit wall (iteration %i). dt=%8.3g t3=%8.3g (intersection:3).\n", NAME_CURRENT_COMP, loopcounter, dt_to_tangent, t3); if (t3 == -3) printf(" No sign change to determine intersection\n"); else if (t3 == -2) printf(" Max iterations reached\n"); else if (t3 < 0 || t3 >= dt_to_tangent) printf(" Error when solving intersection\n"); } /* neutron can not reach slit wall. */ #pragma acc atomic FCVars.warn_notreachslitwall = FCVars.warn_notreachslitwall + 1; ABSORB; } /* Propagate to slit wall point (t+t3) on slit n3 wall */ PROP_DT(t3); dt -= t3; /* dt: time remaining to slit exit after propagation */ xp1 = FC_xrot(x,z, t, FCVars); /* X'(t+t3) : on slit wall */ zp1 = FC_zrot(x,z, t, FCVars); /* Z'(t+t3) : on slit wall */ X[2] = xp1; Z[2] = zp1; if (verbose > 2) printf("FermiChopper: %s: PROP_DT t3=%8.3g dt=%8.3g xyz=[%8.3g %8.3g %8.3g] (on wall). z'=%g\n", NAME_CURRENT_COMP, t3, dt, x,y,z, zp1); /* check if intersection point is still in the slit package, else exit loop */ if (fabs(zp1) > length/2 || dt <= 0) { if (verbose > 2) printf("FermiChopper: %s: Neutron is outside slit pack (on slit wall).\n", NAME_CURRENT_COMP); break; } /* here |---o-----| | / \ | |/ \ | | \ |(dt) |---------| */ /* get velocity in rotating frame, on slit wall */ vxp1 = FC_xrot(vx,vz, t, FCVars); vzp1 = FC_zrot(vx,vz, t, FCVars); q = 2*V2Q*(fabs(vxp1)); { /* double ref = 0; double par[] = {R0, Qc, alpha, m, W};*/ StdReflecFunc(q, par, &ref); if (ref>0) p *= ref; else { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (on slit wall).\n", NAME_CURRENT_COMP); #pragma acc atomic FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1; ABSORB; } /* Cutoff ~ 1E-10 */ } #ifndef OPENACC if (mcdotrace) { double xp2 = x; double zp2 = z; /* indicate position of neutron in mcdisplay */ x = FC_xrot(x,z, t,FCVars); z= FC_zrot(x,z, t,FCVars); SCATTER; x=xp2; z=zp2; } else #endif SCATTER; /* reflect perpendicular velocity and compute new velocity in static frame */ vxp1 *= -1; /* apply transposed Transformation matrix */ vx = FC_xrot( vxp1,-vzp1, t,FCVars); vz = FC_zrot(-vxp1, vzp1, t,FCVars); /* recompute time to slit exit */ /* solve Z'=-slit_input for time of exit of slit package */ t3 = FC_zintersect(x,z,vx,vz, t,dt, -slit_input, FCVars); if(t3 < 0 || t3 > dt) { if (verbose > 1 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) { printf("FermiChopper: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g (intersection:4).\n", NAME_CURRENT_COMP, dt, t3); if (t3 == -3) printf(" No sign change to determine intersection\n"); else if (t3 == -2) printf(" Max iterations reached\n"); else if (t3 < 0) printf(" Error when solving intersection\n"); } #pragma acc atomic FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1; ABSORB; /* neutron can not reach slit output. */ } else dt = t3; /* reduce time interval to [0, time of slit exit] */ } /* end if changed slit index */ else { /* neutron remains in the same slit: direct flight */ if (dt > 0) PROP_DT(dt); /* go to slit exit */ break; } } /* end for */ xp1 = FC_xrot(x,z, t,FCVars); zp1 = FC_zrot(x,z, t,FCVars); X[3] = xp1; Z[3] = zp1; /* slit exit */ if (fabs(xp1) >= nslit*w/2) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package, x'=%8.3g (slit exit).\n", NAME_CURRENT_COMP, xp1); #pragma acc atomic FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1; ABSORB; } if (loopcounter >= FermiChopper_MAXITER) { if (verbose > 1 && FCVars.absorb_maxiterations < FermiChopper_MAXITER) printf("FermiChopper: %s: Max iterations %i reached inside slit. Absorb.\n", NAME_CURRENT_COMP, FermiChopper_MAXITER); #pragma acc atomic FCVars.absorb_maxiterations = FCVars.absorb_maxiterations + 1; ABSORB; } /********************* EXIT SLIT PACKAGE ********************************/ /* propagation times to cylinder. Should use t2 to exit */ if (!cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight)) { /* must be inside the cylinder */ if (verbose > 1) printf("FermiChopper: %s: ABSORB Neutron has unexpectidely exited cylinder ! (exiting)\n", NAME_CURRENT_COMP); #pragma acc atomic FCVars.absorb_exitslitpack = FCVars.absorb_exitslitpack + 1; ABSORB; } if (t1 > 0) { if (verbose > 1 && FCVars.absorb_wrongdirection < FermiChopper_MAXITER) printf("FermiChopper: %s: Neutrons are leaving chopper\n" " in the wrong direction. Absorb.\n", NAME_CURRENT_COMP); #pragma acc atomic FCVars.absorb_wrongdirection = FCVars.absorb_wrongdirection + 1; ABSORB; } if (t2 <= 0 && FCVars.absorb_nocontrol < FermiChopper_MAXITER) { if (verbose > 1) printf("FermiChopper: %s: Neutrons are leaving chopper\n" " without any control. Absorb.\n", NAME_CURRENT_COMP); #pragma acc atomic FCVars.absorb_nocontrol = FCVars.absorb_nocontrol + 1; ABSORB; } /* propagate to cylinder surface (exit) */ PROP_DT(t2); SCATTER; xp1 = FC_xrot(x,z, t,FCVars); zp1 = FC_zrot(x,z, t,FCVars); X[4] = xp1; Z[4] = zp1; if (verbose > 2) printf("FermiChopper: %s: t1=%8.3g PROP_DT t2=%8.3g xyz=[%8.3g %8.3g %8.3g] (OUT cyl). z'=%g\n", NAME_CURRENT_COMP, t1, t2, x,y,z, zp1); /* Check if the neutron left the cylinder by its top or bottom */ if (fabs(y) >= yheight/2) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (exiting)\n", NAME_CURRENT_COMP, y); #pragma acc atomic FCVars.absorb_topbottom = FCVars.absorb_topbottom + 1; ABSORB; } /* must have X'< slit package width at package Z'=cylinder output */ if (fabs(xp1) >= nslit*w/2) { if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package cylinder, xp1=%8.3g (exiting).\n", NAME_CURRENT_COMP, xp1); #pragma acc atomic FCVars.absorb_cylexit = FCVars.absorb_cylexit + 1; ABSORB; } /* Transmission coefficent */ p = p*eff; //finite cross section + transmission #pragma acc atomic FCVars.sum_N_pass = FCVars.sum_N_pass + p; } /* end if cylinder_intersect */ %} SAVE %{ double mean_k, mean_v, mean_t, mean_w=0, mean_L=0; if (FCVars.sum_N) { mean_v = FCVars.sum_v/FCVars.sum_N; mean_t = FCVars.sum_t/FCVars.sum_N; mean_k = V2K*mean_v; if (mean_k) mean_L = 2*PI/mean_k; mean_w = VS2E*mean_v*mean_v; /* base opening time */ double div, mean_phase; div = atan2(w,length)/PI*180; mean_phase = fmod(mean_t*nu*360, 360); mean_phase+=radius/mean_v*360*nu; if (mean_phase > 180) mean_phase -= 360; if (!FCVars.sum_N_pass) printf("FermiChopper: %s: No neutron can pass the chopper.\n", NAME_CURRENT_COMP); if (!FCVars.sum_N_pass || verbose) printf("FermiChopper: %s\n" " Mean velocity v = %g [m/s]\n" " Mean wavelength lambda= %g [Angs]\n" " Mean energy omega = %g [meV]\n" " Mean arrival time t = %g [s]\n" " Mean phase = %g [deg] (%s)\n" " Slit pack divergence = %g [deg] (full width)\n" " Opening time dt = %g [s]\n" " Intensity reaching FC = %g [n/s]\n" " Intensity passing FC = %g [n/s]\n" , NAME_CURRENT_COMP, mean_v, mean_L, mean_w, mean_t, mean_phase, (zero_time > 1 ? "set automatically" : "use phase=-(this value) to optimize"), 2*div, (nu ? fabs(div/PI/nu) : 1), FCVars.sum_N, FCVars.sum_N_pass); if (!FCVars.sum_N_pass || verbose) { printf("FermiChopper: %s: Lost events anaylsis\n" " Already inside: %li\n" " By Top/Bottom of cylinder: %li\n" " At cylinder entrance: %li\n" " Hit cyl. entrance sides: %li\n" " Can't prop. to slit pack: %li\n" " At slit pack entrance: %li\n" " On absorbing slit coating: %li\n" " Exiting slit pack: %li\n" " Too many iterations: %li\n" " Prop. in wrong direction : %li\n" " Mad neutron (no control): %li\n" " At cylinder exit: %li\n" , NAME_CURRENT_COMP, FCVars.absorb_alreadyinside, FCVars.absorb_topbottom, FCVars.absorb_cylentrance, FCVars.absorb_sideentrance, FCVars.absorb_notreachentrance, FCVars.absorb_packentrance, FCVars.absorb_slitcoating, FCVars.absorb_exitslitpack, FCVars.absorb_maxiterations, FCVars.absorb_wrongdirection, FCVars.absorb_nocontrol, FCVars.absorb_cylexit); if (FCVars.warn_notreachslitwall || FCVars.warn_notreachslitoutput) printf("Warning: Can not reach slit wall: %li\n" "Warning: Can not reach slit output: %li\n", FCVars.warn_notreachslitwall, FCVars.warn_notreachslitoutput); } } else { printf("FermiChopper: %s: No neutron can reach the chopper.\n", NAME_CURRENT_COMP); } %} MCDISPLAY %{ double index_x=0; double index_z=0; double xpos, zpos; double Nz,Nx; double ymin=-yheight/2.0; double ymax= yheight/2.0; double omega=FCVars.omega; FCVars.omega=0; // FCVars.ph0 =0; Nz = (FCVars.C_slit ? 4 : 1); Nx = (nslit > 11 ? 11 : nslit); FCVars.C_slit *= -1; /* cylinder top/center/bottom */ circle("xz", 0,ymax,0,radius); circle("xz", 0,0 ,0,radius); circle("xz", 0,ymin,0,radius); /* vertical lines to make a kind of volume */ line( 0 ,ymin,-radius, 0 ,ymax,-radius); line( 0 ,ymin, radius, 0 ,ymax, radius); line(-radius,ymin, 0 ,-radius,ymax, 0 ); line( radius,ymin, 0 , radius,ymax, 0 ); /* slit package */ for (index_x = -Nx/2; index_x < Nx/2; index_x++) { for (index_z = -Nz/2; index_z < Nz/2; index_z++) { double xs1, zs1, zs2; double xp1, xp2, zp1, zp2; zs1 = length*index_z/Nz; zs2 = length*(index_z+1)/Nz; xs1 = w*nslit*index_x/Nx; xp1 = FC_xrot(xs1, zs1, 0, FCVars); xp2 = FC_xrot(xs1, zs2, 0, FCVars); zp1 = FC_zrot(xs1, zs1, 0, FCVars); zp2 = FC_zrot(xs1, zs2, 0, FCVars); multiline(5, xp1, ymin, zp1, xp1, ymax, zp1, xp2, ymax, zp2, xp2, ymin, zp2, xp1, ymin, zp1); } } /* cylinder inner sides containing slit package */ double xp1, xp2, zp1, zp2; xpos = nslit*w/2; zpos = sqrt(radius*radius - xpos*xpos); xp1 = FC_xrot(xpos, -zpos, 0, FCVars); xp2 = FC_xrot(xpos, +zpos, 0, FCVars); zp1 = FC_zrot(xpos, -zpos, 0, FCVars); zp2 = FC_zrot(xpos, +zpos, 0, FCVars); multiline(5, xp1, ymin, zp1, xp1, ymax, zp1, xp2, ymax, zp2, xp2, ymin, zp2, xp1, ymin, zp1); xpos *= -1; xp1 = FC_xrot(xpos, -zpos, 0, FCVars); xp2 = FC_xrot(xpos, +zpos, 0, FCVars); zp1 = FC_zrot(xpos, -zpos, 0, FCVars); zp2 = FC_zrot(xpos, +zpos, 0, FCVars); multiline(5, xp1, ymin, zp1, xp1, ymax, zp1, xp2, ymax, zp2, xp2, ymin, zp2, xp1, ymin, zp1); FCVars.omega=omega; %} END