/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Refractor * * %I * Written by: E. Farhi, B. Cubitt * Date: Oct 2014 * Origin: ILL * * A refractor material/shape, which can be used to model e.g. lenses and prisms. * * %D * Single bulk material shape that can be used as a prism or lens. * * NEUTRON INTERACTION PROCESSES: * The bulk material can reflect, refract, scatter and absorb neutrons, depending * on the material cross sections and incident angles. * * The refracting material is specified from its molar weight, density, coherent * scattering cross section. The refractive index is computed as: * n = sqrt(1-(lambda*lambda*rho*bc/PI)) is the refraction index * * The surface can be coated when specifying a critical wavevector Qc, with e.g. * Qc=m*0.0219 for a super mirror coating. The mirror coating can be suppressed * by setting Qc=0. The critical wavevector is then set to * Qc = 4*sqrt(PI*rho*bc); with * rho= density*6.02214179*1e23*1e-24/weight; * bc = sqrt(fabs(sigma_coh)*100/4/PI)*1e-5; * * COMPONENT GEOMETRY: * The component shape can be a sphere, box, cylinder, biconcave spherical lens * or any other shape defined from an external OFF/PLY file. * sphere: radius * cylinder: radius, yheight * box: xwidth, yheight, zdepth * OFF/PLY: geometry="filename.off or ply", xwidth, yheight, zdepth (bounding box) * lens_sphere: geometry="lens_sphere", radius, zdepth (thickness) * * The lens_sphere geometry is composed of two concave half spheres of same radius, * separated with a minimal thickness zdepth along the Z axis. * * Optionally, you can specify the 'geometry' parameter as a OFF/PLY file name. * The complex geometry option handles any closed non-convex polyhedra. * It computes the intersection points of the neutron ray with the object * transparently, so that it can be used like a regular sample object. * It supports the PLY, OFF and NOFF file format but not COFF (colored faces). * Such files may be generated from XYZ data using: * qhull < coordinates.xyz Qx Qv Tv o > geomview.off * or * powercrust coordinates.xyz * and viewed with geomview or java -jar jroff.jar (see below). * * All geometries are centred. The bulk material fills the shape, but can be * set 'outside' when density is given a negative value. In this case, the material * outside the bulk is void (vacuum). * * Usually, you should stack more than one of these to get a significant effect * on the neutron beam, so-called 'compound refractive lens'. * The focal length for N lenses with focal 'f' is f/N, where f=R/(1-n) * and R = r/2 for a spherical lens with curvature radius 'r' * * COMMON MATERIALS: * Should have high coherent, and low incoherent and absorption cross sections * Be: density=1.85, weight=9.0121, sigma_coh=7.63, sigma_inc=0.0018,sigma_abs=0.0076 * Pb: density=11.115,weight=207.2, sigma_coh=11.115,sigma_inc=0.003, sigma_abs=0.171 * Pb206: sigma_coh=10.68, sigma_inc=0 , sigma_abs=0.03 * Pb208: sigma_coh=11.34, sigma_inc=0 , sigma_abs=0.00048 * Zr: density=6.52, weight=91.224, sigma_coh=6.44, sigma_inc=0.02, sigma_abs=0.185 * Zr90: sigma_coh=5.1, sigma_inc=0 , sigma_abs=0.011 * Zr94: sigma_coh=8.4, sigma_inc=0 , sigma_abs=0.05 * Bi: density=9.78, weight=208.98, sigma_coh=9.148, sigma_inc=0.0084,sigma_abs=0.0338 * Mg: density=1.738, weight=24.3, sigma_coh=3.631, sigma_inc=0.08, sigma_abs=0.063 * MgF2: density=3.148, weight=62.3018,sigma_coh=11.74, sigma_inc=0.0816,sigma_abs=0.0822 * diamond: density=3.52, weight=12.01, sigma_coh=5.551, sigma_inc=0.001, sigma_abs=0.0035 * Quartz/silica: density=2.53, weight=60.08, sigma_coh=10.625,sigma_inc=0.0056,sigma_abs=0.1714 * Si: density=2.329, weight=28.0855,sigma_coh=2.1633,sigma_inc=0.004, sigma_abs=0.171 * Al: density=2.7, weight=26.98, sigma_coh=1.495, sigma_inc=0.0082,sigma_abs=0.231 * Ni: density=8.908, weight=58.69, sigma_coh=13.3, sigma_inc=5.2, sigma_abs=4.49 * Mn: (bc < 0) density=7.21, weight=54.94, sigma_coh=-1.75, sigma_inc=0.4, sigma_abs=13.3 * perfluoropolymer(PTFE/Teflon/CF2): * density=2.2, weight=50.007, sigma_coh=13.584,sigma_inc=0.0026,sigma_abs=0.0227 * Organic molecules with C,O,H,F * * Among the most commonly available and machinable materials are MgF2, SiO2, Si, and Al. * %P * INPUT PARAMETERS: * * xwidth: [m] width of box * yheight: [m] height of box/cylinder * zdepth: [m] depth of box * radius: [m] radius of sphere/cylinder * R0: [1] Low-angle reflectivity * Qc: [Angs-1] critical scattering vector, e.g. Qc=0.0219 for Ni coating. Set Qc=0 to use the bulk critical grazing angles. * sigma_coh: [barn] coherent cross section of refracting material. Use negative value to indicate a negative coherent scattering length * sigma_inc: [barn] incoherent cross section * sigma_abs: [barn] thermal absorption cross section * density: [g/cm3] density of the refracting material. density < 0 means the material is outside/before the shape. * weight: [g/mol] molar mass of the refracting material * geometry: [str] OFF/PLY geometry file name, or NULL to use simple shape A spherical bi-concave lens can be obtained with geometry="lens_sphere" and given radius and zdepth * p_interact: [1] MC Probability for scattering the ray; otherwise transmit. Use 0 to compute true probability, or specify it as e.g. 0.05 * focus_scatter: [deg] angle in which to scatter in bulk, with probability 'p_interact' * RMS: [Angs] root mean square wavyness of the surface * verbose: [1] flag to display detailed component behaviour * p_scatter: [1] flag to allow scattering in the refractor bulk * p_reflect: [1] flag to allow reflection (grazing angle) at the surface * p_refract: [1] flag to allow refraction at the refractor surface * * OUTPUT PARAMETERS: * theta1: [deg] incoming angle to the surface * theta2: [deg] outgoing angle to the surface * SCATTERED: [] 0=transmitted, 1=scattered in bulk, 2=refracted, 3=reflected * * %L * M. L. Goldberger et al, Phys. Rev. 71, 294 - 310 (1947) * %L * Sears V.F. Neutron optics. An introduction to the theory of neutron optical phenomena and their applications. Oxford University Press, 1989. * %L * H. Park et al. Measured operational neutron energies of compound refractive lenses. Nuclear Instruments and Methods B, 251:507-511, 2006. * %L * J. Feinstein and R. H. Pantell. Characteristics of the thick, compound refractive lens. Applied Optics, 42 No. 4:719-723, 2001. * %L * Geomview and Object File Format (OFF) * %L * Java version of Geomview (display only) jroff.jar * %L * Meshlab can view OFF and PLY files * %L * qhull for points to OFF conversion * %L * powercrust for points to OFF conversion * %E *******************************************************************************/ DEFINE COMPONENT Refractor DEFINITION PARAMETERS () SETTING PARAMETERS ( xwidth=0, yheight=0, zdepth=0, radius=0, string geometry="NULL", R0=0.99, sigma_coh=11.74, density=3.148, weight=62.3018, sigma_inc=0, sigma_abs=0, Qc=0, p_interact=0.05, RMS=0, focus_scatter=10, verbose=0, p_scatter=1, p_reflect=1, p_refract=1) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "read_table-lib" %include "interoff-lib" %include "ref-lib" /* wrappers to intersection routines which also return the normal vector */ int box_intersect_n(double *t0, double *t1, double x, double y, double z, double vx, double vy, double vz, double dx, double dy, double dz, double *nx, double *ny, double *nz) { double dt=0; int intersect = box_intersect(t0, t1, x,y,z, vx,vy,vz, dx,dy,dz); /* determine normal vector depending on hit face */ if (intersect) { dt = *t0 <= 0 || *t0 > *t1 ? *t1 : *t0; if (dt < 0) intersect = 0; } if (intersect && dx && dy && dz) { x += vx*dt; y += vy*dt; z += vz*dt; /* determine hit face: difference to plane is closest to 0 */ *nx = trunc(2.002*x/dx); *ny = trunc(2.002*y/dy); *nz = trunc(2.002*z/dz); } return(intersect); } // box_intersect_n /* sphere: when radius < 0 we always use the further intersection point to determine the normal vector */ int sphere_intersect_n(double *t0, double *t1, double x, double y, double z, double vx, double vy, double vz, double r, double *nx, double *ny, double *nz) { double dt=0; int intersect = sphere_intersect(t0, t1, x,y,z, vx,vy,vz, r); if (intersect) { if (r >0) /* First intersection in positive time */ dt = *t0 <= 0 || *t0 > *t1 ? *t1 : *t0; else /* always second intersection */ dt = *t1; if (dt < 0) intersect = 0; } if (intersect && r) { /* must propagate locally to determine normal vector. */ x += vx*dt; y += vy*dt; z += vz*dt; *nx = x; *ny = y; *nz = z; } return(intersect); } // sphere_intersect_n int cylinder_intersect_n(double *t0, double *t1, double x, double y, double z, double vx, double vy, double vz, double r, double h, double *nx, double *ny, double *nz) { double dt=0; int intersect = cylinder_intersect(t0, t1, x,y,z, vx,vy,vz, r,h); if (intersect) { if (r >0) /* First intersection in positive time */ dt = *t0 <= 0 || *t0 > *t1 ? *t1 : *t0; else /* always second intersection */ dt = *t1; if (dt < 0) intersect = 0; } if (intersect && r) { /* must propagate locally to determine normal vector */ x += vx*dt; z += vz*dt; *nx = x; *ny = 0; /* cylinder is vertical */ *nz = z; } return(intersect); } // cylinder_intersect_n // two sphere separated with a given thickness int lens_sphere_intersect_n(double *t0, double *t1, double x, double y, double z, double vx, double vy, double vz, double r, double dz, double *nx, double *ny, double *nz) { /* two spherical lenses, concave, with given radius, thickness=zdepth at meniscus, plus cross section=xwidth*yheight */ int intersect1 = 0, intersect2 = 0; double nx1=0,nx2=0,ny1=0,ny2=0,nz1=0,nz2=0; double t2=0,t3=0; *t0=*t1=0; intersect1 = sphere_intersect_n(t0, t1, x,y,z+r+dz/2, vx,vy,vz, -r, &nx1,&ny1,&nz1); /* -r: use the further intersection with sphere */ intersect2 = sphere_intersect_n(&t2, &t3, x,y,z-r-dz/2, vx,vy,vz, r, &nx2,&ny2,&nz2); /* usual case with 4 intersections: return t1 and t2 */ if (intersect1 || intersect2) { if (intersect1) *t0=*t1; if (intersect2) *t1= t2; } if (!intersect1 && !intersect2) return(0); if (!intersect1 || !intersect2) { if (!intersect1) { /* intersection with the 2st sphere: return t2 and some other external intersection (container box entrance) t0 */ intersect1 = box_intersect_n(t0,&t3, x,y,z, vx,vy,vz, 2*r, 2*r, 2*r+dz, &nx1,&ny1,&nz1); } else { /* intersection with the 1st sphere: return t1 and some other external intersection (container box exit) t3 */ intersect2 = box_intersect_n(&t2,t1, x,y,z, vx,vy,vz, 2*r, 2*r, 2*r+dz, &nx2,&ny2,&nz2); } } /* the normal vector corresponds to the first forward intersection */ if (intersect1 || intersect2) { if (*t0 < 0) { *nx = nx2; *ny = ny2; *nz = nz2; } else { *nx = nx1; *ny = ny1; *nz = nz1; } } /* no intersection: use intersection with container box */ return(intersect1+intersect2); } // lens_sphere_intersect_n /* Surface_wavyness: function to rotate normal vector around axis for roughness * with specified tilt angle (rad) * RETURNS: rotated nx,ny,nz coordinates */ #pragma acc routine void Surface_wavyness(double *nx, double *ny, double *nz, double tilt, _class_particle *_particle) { double nt_x, nt_y, nt_z; /* transverse vector */ double n1_x, n1_y, n1_z; /* normal vector (tmp) */ /* normal vector n_z = [ 0,1,0], n_t = n x n_z; */ vec_prod(nt_x,nt_y,nt_z, *nx,*ny,*nz, 0,1,0); /* rotate n with angle wav_z around n_t -> n1 */ tilt /= (sqrt(8*log(2)))*randnorm(); rotate(n1_x,n1_y,n1_z, *nx,*ny,*nz, tilt, nt_x,nt_y,nt_z); /* rotate n1 with angle phi around n -> nt */ rotate(nt_x,nt_y,nt_z, n1_x,n1_y,n1_z, 2*PI*rand01(), *nx,*ny,*nz); *nx=nt_x; *ny=nt_y; *nz=nt_z; } %} DECLARE %{ double rho; double bc; off_struct offdata; double mean_n; double events; double theta1; double theta2; %} INITIALIZE %{ mean_n=0; events=0; theta1=0; theta2=0; /* set geometry from input geometry parameters */ if (!geometry || !strlen(geometry) || !strcmp(geometry, "NULL") || !strcmp(geometry, "0")) { if (radius && yheight && !xwidth) strcpy(geometry, "cylinder"); else if (xwidth && yheight && zdepth && !radius) strcpy(geometry, "box"); else if (radius && !yheight && !xwidth && !zdepth) strcpy(geometry, "sphere"); } else if (radius && zdepth && !strcmp(geometry, "lens_sphere")) { /* NOP: OK */ } else if (strcmp(geometry, "NULL") && strcmp(geometry, "0") && xwidth && yheight && zdepth) { #ifndef USE_OFF fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n"); exit(-1); #else /* off/ply case */ if (!off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) exit(printf("Refractor: %s: FATAL: invalid OFF/PLY geometry specification for file '%s'.\n", NAME_CURRENT_COMP, geometry)); #endif } else exit(printf("Refractor: %s: FATAL: invalid geometry specification.\n" " Check geometry,xwidth,yheight,zdepth,radius\n", NAME_CURRENT_COMP)); /* compute refraction parameters, if needed */ if (density==0 || weight<=0) exit(printf("Refractor: %s: FATAL: invalid material density or molar weight: density=%g weight=%g\n", NAME_CURRENT_COMP, density, weight)); rho= fabs(density)*6.02214179*1e23*1e-24/weight; /* per at/Angs^3 */ if (sigma_coh==0) exit(printf("Refractor: %s: FATAL: invalid material coherent cross section: sigma_coh=%g\n", NAME_CURRENT_COMP, sigma_coh)); bc=sqrt(fabs(sigma_coh)*100/4/PI)*1e-5; /* bound coherent scattering length */ if (sigma_coh<0) bc *= -1; /* use bulk to compute reflectivity critical angle/wavevector */ if (Qc <=0) Qc = 4*sqrt(PI*rho*fabs(bc)); MPI_MASTER( printf("Refractor: %s: rho.bc=%g [10-6 Angs-2] Qc=%g [Angs-1]. geometry=%s\n", NAME_CURRENT_COMP, rho*bc*1e6, Qc, geometry); ); %} TRACE %{ int intersect = 0; int iterations= 0; char event[256]; theta1=theta2=0; #ifdef OPENACC #ifdef USE_OFF off_struct thread_offdata = offdata; #endif #else #define thread_offdata offdata #endif do { double nx=0, ny=0, nz=0; double t0=0, t1=0, dt=0; intersect = 0; iterations++; #ifndef OPENACC strcpy(event, " "); #endif /* determine intersection times and normal vector with geometry */ if (!strcmp(geometry, "sphere")) { intersect = sphere_intersect_n(&t0, &t1, x,y,z, vx,vy,vz, radius, &nx,&ny,&nz); } else if (!strcmp(geometry, "lens_sphere")) { intersect = lens_sphere_intersect_n(&t0, &t1, x,y,z, vx,vy,vz, radius, zdepth, &nx,&ny,&nz); } else if (!strcmp(geometry, "cylinder")) { intersect = cylinder_intersect_n(&t0, &t1, x,y,z, vx,vy,vz, radius, yheight, &nx,&ny,&nz); } else if (!strcmp(geometry, "box")) { intersect = box_intersect_n(&t0, &t1, x,y,z, vx,vy,vz, xwidth, yheight, zdepth, &nx,&ny,&nz); } #ifdef USE_OFF else if (geometry && strcmp(geometry, "NULL") && strcmp(geometry, "0")) { Coords n0, n1; intersect = off_intersect(&t0, &t1, &n0, &n1, x, y, z, vx, vy, vz, 0, 0, 0, thread_offdata ); if (intersect) { coords_get(t0 <= 0 || t0 > t1 ? n1 : n0, &nx, &ny, &nz); } } #endif if (intersect) { if (t1 < t0) { double tmp=t0; t0=t1; t1=tmp; } dt = t0 <= 1e-10 ? t1 : t0; /* full propagation time inside geometry */ if (dt < 0) dt=0; } if (dt<=0) break; #ifndef OPENACC strcpy(event, "none"); #endif if (verbose) printf("%s:%i: [neutron=%li iteration=%i] intersect %s=%i t0=%g t1=%g dt=%g event=(intersect) %s\n", NAME_CURRENT_COMP, __LINE__, mcget_run_num(), iterations, geometry, intersect, t0,t1, dt, event); /* scattering/absorption in bulk, refraction/reflection at interface */ if (intersect) { char scatter_me=0, refract_me=0; double my_t=0, d_path=0; double v=0,lambda=0; char inbulk = t0 < 0 && t1 > 0; /* when density is given negative, we assume bulk refractive material is 'outside' or 'before' the shape */ if (density < 0) inbulk = !inbulk; v = sqrt(vx*vx+vy*vy+vz*vz); if (!v) ABSORB; lambda = 3956.0032/v; d_path = v*dt; /* full length in material */ /* compute probability to scatter/absorb inside bulk */ if (inbulk) { /* scatter inside bulk */ double ws = 0; double p_trans= 0, p_scatt=0, mc_trans=0, mc_scatt=0; int flag = 0; double my_a_v = (rho * 100 * sigma_abs); double my_a = my_a_v*2200/v; double my_s = rho * 100 *(sigma_inc+sigma_coh); my_t = my_a + my_s; ws = my_s/my_t; /* (inc+coh)/(inc+coh+abs) */ /* Proba of transmission along length d_path */ p_trans = exp(-my_t*d_path); p_scatt = 1 - p_trans; flag = 0; /* flag used for propagation to exit point before ending */ /* are we next to the exit ? probably no scattering (avoid rounding errors) */ if (my_s*d_path <= 4e-7) { flag = 1; /* No interaction before the exit */ } /* force a given fraction of the beam to scatter */ if (p_interact>0 && p_interact<=1) { /* we force a portion of the beam to interact */ /* This is used to improve statistics on single scattering (and multiple) */ mc_trans = 1-p_interact; } else { mc_trans = p_trans; /* 1 - p_scatt */ } mc_scatt = 1 - mc_trans; /* portion of beam to scatter (or force to) */ if (mc_scatt <= 0 || mc_scatt>1) flag=1; /* MC choice: Interaction or transmission ? */ scatter_me = !flag && ws && p_scatt && mc_scatt > 0 && (mc_scatt >= 1 || rand01() < mc_scatt); /* account for absorption and retain scattered fraction */ /* we have chosen portion mc_scatt of beam instead of p_scatt, so we compensate */ if (scatter_me) p *= ws * fabs(p_scatt/mc_scatt); else if (p_trans && mc_trans) p *= fabs(p_trans/mc_trans); /* attenuate beam by portion which is scattered (and left along) */ if (verbose) printf("%s:%i: [neutron=%li iteration=%i] intersect %s=%i t0=%g t1=%g dt=%g event=(inbulk) %s scatter_me=%i\n", NAME_CURRENT_COMP, __LINE__, mcget_run_num(), iterations, geometry, intersect, t0,t1, dt, event, scatter_me); } /* if inbulk */ /* handle absorption and scattering in material */ /* req: my_t, d_path, focus_ah, focus_aw */ if (p_scatter && scatter_me) { double dl=0, solid_angle=0; double vx2=0, vy2=0, vz2=0; if (my_t*d_path < 1e-6){ /* For very weak scattering, use simple uniform sampling of scattering point to avoid rounding errors. */ dl = rand0max(d_path); /* length */ } else { double a = rand0max((1 - exp(-my_t*d_path))); dl = -log(1 - a) / my_t; /* length */ } PROP_DT(dl/v); dt = 0; // we propagate in the bulk SCATTER; /* 1 */ /* scatter randomly in cone to take into account material sigma */ randvec_target_circle(&vx2, &vy2, &vz2, &solid_angle, vx,vy,vz, focus_scatter*DEG2RAD); vx=vx2; vy=vy2; vz=vz2; if (solid_angle) { p *= solid_angle/4/PI; NORM(vx, vy, vz); vx*=v; vy*=v; vz*=v; /* scattering is elastic */ } /* force to recompute intersection with new neutron direction/position */ refract_me = 0; #ifndef OPENACC strcpy(event, "scatter"); #endif if (verbose) printf("%s:%i: [neutron=%li iteration=%i] intersect %s=%i t0=%g t1=%g dt=%g event=(scatter) %s\n", NAME_CURRENT_COMP, __LINE__, mcget_run_num(), iterations, geometry, intersect, t0,t1, dl/v, event); } else { refract_me = 1; } /* if scatter */ if ((p_refract||p_reflect) && refract_me) { // refract or reflect on the surface double n=0; /* refractive index */ double n1=0, n2=0; double q=0, R=0; double par[] = {R0, Qc, 6.0, 1.0, 1.0/300.0}; /* propagate to surface */ PROP_DT(dt); dt=0; SCATTER; SCATTER; /* 2 */ /* compute refraction index */ n = sqrt(1-(lambda*lambda*rho*bc/PI)); mean_n += n; events++; /* compute incoming angle */ if (inbulk) { n1=n; n2=1; } /* from bulk to void */ else { n1=1; n2=n; } /* from void to bulk */ /* compute reflectivity */ q = fabs(2*vz*V2Q); /* Reflectivity (see component Guide). */ StdReflecFunc(q, par, &R); /* tilt normal vector for roughness, in cone theta_RMS */ NORM(nx,ny,nz); /* cone angle from RMS roughness = atan(2*RMS/lambda) */ if (RMS>0) Surface_wavyness(&nx, &ny, &nz, atan(2*RMS/lambda), _particle); /* Snell-Descartes formula for refraction n1 sin(theta1) = n2 sin(theta2) */ /* https://en.wikipedia.org/wiki/Snell's_law */ Coords N = coords_set(nx,ny,nz); // normal vector to surface Coords V = coords_set(vx,vy,vz); // incoming velocity Coords I = coords_scale(V, 1/v); // normalised ray = v/|v| // theta1: incident angle to the surface normal double cos_theta1 = -coords_sp(N,I); // cos(theta1) = -N.I if (fabs(cos_theta1) > 1) ABSORB; // should never occur... theta1 = acos(cos_theta1)*RAD2DEG; // reflected ray: probability R // reflected beam: I + 2cos(theta1).N Coords I_reflect = coords_add(I, coords_scale(N, 2*cos_theta1)); // reflected velocity: I_reflect.v Coords V_reflect = coords_scale(I_reflect, v); // compute refracted angle theta2... double sqr_cos_theta2 = 1-(n1/n2)*(n1/n2)*(1-cos_theta1*cos_theta1); // now choose which one to use, and compute outgoing velocity if (0 < sqr_cos_theta2 && sqr_cos_theta2 < 1) { // refraction is possible // theta2: refracted angle to the surface normal double cos_theta2= sqrt(sqr_cos_theta2); // select reflection (or refraction) from Monte-Carlo choice with probability R // in this case we expect R to be small (q > Qc) if (p_reflect && 0 < R && R < 1 && rand01() < R) { // choose reflection from MC theta2 = theta1; coords_get(V_reflect, &vx, &vy, &vz); #ifndef OPENACC strcpy(event, "reflect"); #endif } else if (p_refract) { // compute refracted ray theta2 = acos(cos_theta2)*RAD2DEG; Coords I_refract = coords_add(coords_scale(I, n1/n2), coords_scale(N, n1/n2*cos_theta1 + (cos_theta1 < 0 ? cos_theta2 : -cos_theta2) )); Coords V_refract = coords_scale(I_refract, v); coords_get(V_refract, &vx, &vy, &vz); #ifndef OPENACC strcpy(event, "refract"); #endif SCATTER; /* 3 */ } } else if (p_reflect) { // only reflection: below total reflection theta2 = theta1; if (0 < R && R < 1) p *= R; // should be R0 coords_get(V_reflect, &vx, &vy, &vz); #ifndef OPENACC strcpy(event, "reflect"); #endif } /* propagate by a small time so that we leave the surface */ PROP_DT(1e-9); if (verbose) printf("%s:%i: [neutron=%li iteration=%i] intersect %s=%i t0=%g t1=%g dt=%g event=(reflect/refract) %s theta1=%g theta2=%g 1-n=%g |xy|=%g nz=%g\n", NAME_CURRENT_COMP, __LINE__, mcget_run_num(), iterations, geometry, intersect, t0,t1, dt, event, theta1, theta2, 1-n, sqrt(x*x+y*y), nz); } /* if reflect/refract */ } /* if intersect */ else break; if (verbose) printf("%s:%i: [neutron=%li iteration=%i] intersect %s=%i t0=%g t1=%g dt=%g event=%s\n", NAME_CURRENT_COMP, __LINE__, mcget_run_num(), iterations, geometry, intersect, t0,t1, dt, event); } while (intersect && iterations<100); %} FINALLY %{ /* * The focal length for N lenses with focal 'f' is f/N, where f=R/(1-n) * and R = r/2 for a spherical lens with curvature radius 'r' */ if (radius && !strcmp(geometry, "lens_sphere") ) { double focus, R; mean_n /= events; R = radius/2; focus = R/(1-mean_n); MPI_MASTER( printf("Refractor: %s: %s focal length f=%g [m]. Focal length for N lenses is f/N.\n" "mean n=%g, events=%g\n", NAME_CURRENT_COMP, geometry, focus, mean_n , events); ); } %} MCDISPLAY %{ /* show geometry */ if (!strcmp(geometry, "sphere")) { circle("xy", 0, 0, 0, radius); circle("xz", 0, 0, 0, radius); circle("yz", 0, 0, 0, radius); } else if (!strcmp(geometry, "cylinder")) { circle("xz", 0, yheight/2.0, 0, radius); circle("xz", 0, -yheight/2.0, 0, radius); line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0); line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0); line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius); line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius); } else if (!strcmp(geometry, "box")) { box(0,0,0, xwidth, yheight, zdepth); } else if (!strcmp(geometry, "lens_sphere")) { // sphere: x^2+y^2+z^2=radius. In 2D: radius^2=z^2+r^2 int index; for (index=0; index<=3; index++) { double z=radius*(index == 3 ? 0.95 : index/3.0); double r=sqrt(radius*radius-z*z); circle("xy", 0, 0, z-(radius+zdepth/2), r); circle("xy", 0, 0,-z+(radius+zdepth/2), r); } box(0,0,0, 2*radius, 2*radius, 2*radius+zdepth); } else if (!strcmp(geometry, "lens_parabola")) { // parabola: z = (x^2+y^2)/4/radius. In 2D: z = r^2/4/radius int index; for (index=0; index<=3; index++) { double z=radius*(index == 3 ? 0.95 : index/3.0); double r=sqrt(z*4*radius); circle("xy", 0, 0,-z-(zdepth/2), r); circle("xy", 0, 0, z+(zdepth/2), r); } box(0,0,0, 4*radius, 4*radius, 2*radius+zdepth); } else if (geometry && strcmp(geometry, "NULL") && strcmp(geometry, "0")) { off_display(offdata); } %} END