/*******************************************************************************
*
* 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
*
* CALCULATED 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
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)
/* 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, 0, 0, 1, 0);
} 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, 0, 0, 1, 0);
} 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, 0, 0, 1, 0);
} else if (geometry && strcmp (geometry, "NULL") && strcmp (geometry, "0")) {
off_display (offdata);
}
%}
END