/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2009, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Lens
*
* %Identification
*
* Written by: C. Monzat/E. Farhi/S. Desert/G. Euzen
* Date: 2009
* Origin: ILL/LLB
*
* Refractive lens with absorption, incoherent scattering and surface imperfection.
*
* %Description
* Refractive Lens with absorption, incoherent scattering and surface imperfection.
* Geometry may be:
* spherical (use r1 and r2 to specify radius of curvature),
* planar (use phiy1 and phiy2 to specify rotation angle of plane w.r.t beam)
* parabolic (use focus1 and focus2 as optical focal length).
*
* 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).
*
* The lens cross-section is seen as a 2*radius disk from the beam Z axis, except
* when a 'geometry' file is given.
* 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'
* R = focus for a parabolic lens with focal of the parabola
* and n = sqrt(1-(lambda*lambda*rho*bc/PI)) with
* bc = sqrt(fabs(sigma_coh)*100/4/PI)*1e-5
* rho = density*6.02214179*1e23*1e-24/weight
*
* 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
* 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
*
* Example: Lens(r1=0.025,r2=0.025, thickness=0.001,radius=0.0150)
*
* %BUGS
* parabolic shape is not fully validated yet, but should do still...
*
* %Parameters
* INPUT PARAMETERS:
* r1: [m] radius of the first circle describing the lens. r>0 means concave face, r<0 means convex, r=0 means plane
* r2: [m] radius of the second circle describing the lens r>0 means concave face, r<0 means convex, r=0 means plane
* focus1: [m] focal of the first parabola describing the lens
* focus2: [m] focal of the second parabola describing the lens
* phiy1: [deg] angle of plane1 (r1=0) around y vertical axis
* phiy2: [deg] angle of plane2 (r2=0) around y vertical axis
* thickness: [m] thickness of the lens between its two surfaces
* radius: [m] radius of the lens section, e.g. beam size.
* density: [g/cm3] density of the material for the lens
* weight: [g/mol] molar mass of the material
* focus_aw: [deg] vertical angle to forward focus after scattering
* focus_ah: [deg] horizontal angle to forward focus after scattering
* sigma_coh: [barn] coherent cross section
* sigma_inc: [barn] incoherent cross section
* sigma_abs: [barn] thermal absorption cross section
* p_interact: [1] MC Probability for scattering the ray; otherwise transmit
* RMS: [Angs] root mean square roughness of the surface
* geometry: [str] Name of an Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust
*
* %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 operationnal 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
* qhull
* %L
* powercrust
* %E
*******************************************************************************/
DEFINE COMPONENT Lens
DEFINITION PARAMETERS ()
SETTING PARAMETERS (r1=0,r2=0,focus1=0,focus2=0,phiy1=0,phiy2=0,
thickness=0.001,radius=0.015,
sigma_coh=11.74,sigma_inc=0.0816,sigma_abs=0.0822,
density=3.148,weight=62.3018,
p_interact=0.1,focus_aw=10,focus_ah=10,RMS=0, string geometry=0)
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
/* support for OFF/PLY geometry */
%include "read_table-lib"
%include "interoff-lib"
/* Lens_roughness: function to rotate normal vector around axis for roughness
* with specified tilt angle (deg)
* RETURNS: rotated nx,ny,nz coordinates
*/
#pragma acc routine seq
void Lens_roughness(double *nx, double *ny, double *nz, double tilt
#ifdef OPENACC
, _class_particle* _particle
#endif
)
{
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 *= DEG2RAD/(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;
}
/* parabola_intersect: Calculate intersection between a line and a parabola with
* axis along z, focal length f, centered at (0,0,0)
* RETURNS 0 when no intersection is found
* or 1 in case of intersection with resulting times t0 and t1 */
#pragma acc routine seq
int parabola_intersect(double *t0, double *t1, double x, double y, double z,
double vx, double vy, double vz, double f)
{
/* equation of line: (x,y,z) = (vx,vy,vz)*t+(x,y,z)_0 */
/* equation of parabola: z = (x^2+y^2)/4/f
that is: 4*f*(vz*t+z0) = (vx*t+x0)^2 + (vy*t+y0)^2 */
double A = vx*vx+vy*vy;
double B = 2*vx*x+2*vy*y-4*vz*f;
double C = x*x+y*y-4*f*z;
return(solve_2nd_order(t0, t1, A,B,C));
}
/* Lens_intersect: function to compute intersection with lens surface
* of section 2*radius (beam size)
* which can be
* type="spherical" (arg=curvature radius),
* type="parabolic" (arg=focus) or
* type="planar" (arg=phi around y).
* The Surface intersection along the Z axis is 'shift' (e.g. +/-thickness/2)
* RETURNS: intersection time to surface, or negative for no intersection
* neutron must then be propagated on the surface with PROP_DT,
* and checked for actual intersection.
*/
#pragma acc routine seq
double Lens_intersect(double x,double y,double z,
double vx,double vy,double vz,
double shift, double radius,
int type, double arg)
{
double t0, t1;
if (!vz || !type || !radius) return -1;
if (type==1 && arg) { /* spherical */
if (sphere_intersect(&t0,&t1, x,y,z+shift+arg, vx,vy,vz, fabs(arg))) {
return ( arg > 0 ? t1 : t0 ); /* concave : convex surface */
} else
return(-1);
} else if (type==2 && arg) { /* parabolic */
if (parabola_intersect(&t0, &t1, x,y,z+shift, vx,vy,vz, arg)) {
if (t0 < 0 && 0 < t1) return(t1);
if (t1 < 0 && 0 < t0) return(t0);
return(t1 < t0 ? t1 : t0);
} else
return(-1);
/* end parabola */
} else if (type==3) { /* planar */
if (plane_intersect(&t0, x,y,z+shift, vx,vy,vz, -sin(arg),0,-cos(arg), 0,0,0)>0)
return(t0);
else
return(-1);
} else /* unsupported geometry */
return(-1);
}
%}
DECLARE
%{
double rho;
double bc;
double my_s;
double my_a_v;
double mean_n;
double events;
off_struct offdata;
%}
INITIALIZE
%{
/* density: Number of atoms by AA^3, pow(10,-24) stands for conversion from cm^3 to AA^3 */
if (density<=0 || weight<=0)
exit(printf("Lens: %s: FATAL: invalid material density or molar weight: density=%g weight=%g\n",
NAME_CURRENT_COMP, density, weight));
rho= density*6.02214179*1e23*1e-24/weight;
if (sigma_coh==0)
exit(printf("Lens: %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 cross section */
if (sigma_coh<0) bc *= -1;
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
#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
/* init the OFF/PLY without centering and scaling */
if (!off_init(geometry, 0, 0, 0, 1, &offdata))
exit(printf("Lens: %s: FATAL: could not initialize geometry file %s [OFF/PLY]\n",
NAME_CURRENT_COMP, geometry));
#endif
}
phiy1 *= DEG2RAD; /* planes orientation */
phiy2 *= DEG2RAD;
focus_aw *= DEG2RAD;
focus_ah *= DEG2RAD;
my_s = rho * 100 *(sigma_inc+sigma_coh);
my_a_v= rho * 100 * sigma_abs;
%}
TRACE
%{
double n; /* refractive index */
double v, lambda; /* neutron speed and wavelength */
double dt; /* time to intersection */
double theta_RMS;
double nx=0, ny=0, nz=0; /* normal vector to surface */
double ax, ay, az; /* temporary vector for surface refraction */
double alpha1, beta1, theta1, theta2; /* angles for refraction computation */
#ifdef OPENACC
#ifdef USE_OFF
off_struct thread_offdata = offdata;
#endif
#else
#define thread_offdata offdata
#endif
/* iterate two faces intersection until no more is possible */
do {
/* ======================== First face of the lens ========================== */
/* determine intersection time */
nx=ny=nz=0;
#ifdef USE_OFF
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
double t1=0,t0=0;
Coords n, n0, n1;
int intersect=off_intersect (&t0, &t1, &n0, &n1, x,y,z, vx,vy,vz, 0, 0, 0, thread_offdata);
if (!intersect) dt=-1;
else if (t0 < 0 && 0 < t1) { dt = t1; n=n1; }
else if (t1 < 0 && 0 < t0) { dt = t0; n=n0; }
else {
if (t1 < t0) { dt=t1; n=n1; }
else { dt=t0; n=n0; }
}
coords_get(n, &nx, &ny, &nz);
} else
#endif
if (r1 && !focus1) {
dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, 1, r1);
} else if (!r1 && focus1){
if (focus1>0){
dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, 2, -focus1);
}else{
dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2 + radius*radius/(4*fabs(focus1)) , radius, 2, -focus1);
}
} else {
dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, 3, phiy1);
}
if (dt <= 0) break; /* no intersection: exit from main loop */
/* propagate to surface (1) */
PROP_DT(dt);
/* check for lens cross section (radius) */
if (radius && x*x+y*y > radius*radius) break; /* outside from cross section: exit from main loop */
SCATTER;
v=sqrt(vx*vx+vy*vy+vz*vz);
if (v) lambda=3956.0032/v; else ABSORB;
/* compute refractive index */
/* without inc nor abs: n = sqrt(1-(lambda*lambda*rho*bc/PI)); */
/* M. L. Goldberger et al, Phys. Rev. 71, 294 - 310 (1947) */
/* alpha1 = (sigma_inc+sigma_abs*2200/v)/2/lambda;
n = 1-lambda*lambda*rho/2/PI*sqrt(bc*bc-alpha1*alpha1); */
n = sqrt(1-(lambda*lambda*rho*bc/PI));
#pragma acc atomic
mean_n += n*p;
#pragma acc atomic
events += p;
theta_RMS=atan(2*RMS/lambda); /* cone angle from RMS roughness */
/* compute normal vector (1) when not given by OFF/PLY */
if (!nx && !ny && !nz) {
dt = 0;
if (r1 && !focus1) {
double sign=r1/fabs(r1);
nx=-sign*x;
ny=-sign*y;
nz=sign*(-z-thickness/2-r1);
} else if (!r1 && focus1) {
nx=-x/2/focus1;
ny=-y/2/focus1;
nz=-1;
} else {
nx=-sin(phiy1);
ny=0;
nz=-cos(phiy1);
}
}
/* tilt normal vector for roughness, in cone theta_RMS */
if (RMS>0) Lens_roughness(&nx, &ny, &nz, theta_RMS
#ifdef OPENACC
, _particle
#endif
);
/* compute incoming angles w.r.t surface */
NORM(nx,ny,nz);
vec_prod(ax,ay,az,nx,ny,nz,-vx,-vy,-vz);
/* if n and v are parallel - no refraction.*/
if (ax!=0 || ay!=0 || az!=0){
theta1 = atan2(sqrt(ax*ax+ay*ay+az*az),scalar_prod(nx,ny,nz,-vx,-vy,-vz));
/* Fresnel formula for refraction */
theta2 = asin(sin(theta1)/n);
/* compute new velocity after refraction */
alpha1 = (sin(theta2))/(sin(theta1));
beta1 = (alpha1*v*cos(theta1))-(v*cos(theta2));
vx = beta1*nx + alpha1*vx;
vy = beta1*ny + alpha1*vy;
vz = beta1*nz + alpha1*vz;
}
/* ======================== Second face of the lens ========================= */
/* determine intersection time to go to second surface */
nx=ny=nz=0;
#ifdef USE_OFF
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
double t1=0,t0=0;
Coords n, n0, n1;
int intersect=off_intersect (&t0, &t1, &n0, &n1, x,y,z, vx,vy,vz, 0, 0, 0, thread_offdata);
if (!intersect) dt=-1;
else if (t0 < 0 && 0 < t1) { dt = t1; n=n1; }
else if (t1 < 0 && 0 < t0) { dt = t0; n=n0; }
else {
if (t1 < t0) { dt=t1; n=n1; }
else { dt=t0; n=n0; }
}
coords_get(n, &nx, &ny, &nz);
} else
#endif
if (r2 && !focus2)
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 1, -r2);
else if (!r2 && focus2){
if (focus2>0){
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 2, focus2);
}else{
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2 - radius*radius/(4*fabs(focus2)) , radius, 2, focus2);
}
}else{
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 3, phiy2);
}
if (dt <= 0) break; /* no intersection: exit from main loop */
/* =============== Absorption and scattering between the two faces =========== */
double my_a,my_t,p_trans,p_scatt,mc_trans,mc_scatt,ws,d_path;
double p_mult=1;
int flag=0;
double l_i=0;
double solid_angle=0;
my_a = my_a_v*(2200/v);
my_t = my_a + my_s;
ws = my_s/my_t; /* (inc+coh)/(inc+coh+abs) */
d_path = v * dt;
/* 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 ? */
if (!flag && mc_scatt > 0 && (mc_scatt >= 1 || (rand01()) < mc_scatt)) { /* Interaction neutron/sample */
p_mult *= ws; /* Update weight ; account for absorption and retain scattered fraction */
/* we have chosen portion mc_scatt of beam instead of p_scatt, so we compensate */
p_mult *= fabs(p_scatt/mc_scatt); /* lower than 1 */
} else {
flag = 1; /* Transmission : no interaction neutron/sample */
p_mult *= fabs(p_trans/mc_trans); /* attenuate beam by portion which is scattered (and left along) */
}
if (flag) { /* propagate directly to secound surface */
if (!isnan(p_mult)) {
p *= p_mult; /* apply absorption correction */
} else {
ABSORB;
}
}
else
{
double a;
if (my_t*d_path < 1e-6){
/* For very weak scattering, use simple uniform sampling of scattering
point to avoid rounding errors. */
dt = rand0max(d_path); /* length */
} else {
a = rand0max((1 - exp(-my_t*d_path)));
dt = -log(1 - a) / my_t; /* length */
}
l_i = dt;/* Penetration in sample: scattering+abs */
dt /= v; /* Time from present position to scattering point */
PROP_DT(dt); /* Point of scattering */
randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,0,0,1, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
NORM(vx, vy, vz);
vx*=v;
vy*=v;
vz*=v;
p_mult *= solid_angle/4/PI;
if (!isnan(p_mult)) {
p *= p_mult;
} else {
ABSORB;
}
SCATTER;
/* recompute new intersection time to go to second surface (velocity changed during scattering event) */
nx=ny=nz=0;
#ifdef USE_OFF
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
double t1=0,t0=0;
Coords n, n0, n1;
int intersect=off_intersect (&t0, &t1, &n0, &n1, x,y,z, vx,vy,vz, 0, 0, 0, thread_offdata);
if (!intersect) dt=-1;
else if (t0 < 0 && 0 < t1) { dt = t1; n=n1; }
else if (t1 < 0 && 0 < t0) { dt = t0; n=n0; }
else {
if (t1 < t0) { dt=t1; n=n1; }
else { dt=t0; n=n0; }
}
coords_get(n, &nx, &ny, &nz);
} else
#endif
if (r2 && !focus2)
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 1, -r2);
else if (!r2 && focus2)
if (focus2>0){
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 2, focus2);
}else{
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2 - radius*radius/(4*fabs(focus2)) , radius, 2, focus2);
}
else
dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, 3, phiy2);
if (dt <= 0) break; /* no intersection: exit from main loop */
} /* end scattering handling in material*/
/* propagate to surface */
PROP_DT(dt);
/* check for lens cross section (radius) */
if (radius && x*x+y*y > radius*radius) break; /* outside from cross section: exit from main loop */
SCATTER;
v=sqrt(vx*vx+vy*vy+vz*vz);
/* compute normal vector (2) when not given by OFF/PLY */
if (!nx && !ny && !nz) {
dt = 0;
if (r2 && !focus2) {
double sign=r2/fabs(r2);
nx=sign*x;
ny=sign*y;
nz=sign*(z-thickness/2-r2);
} else if (!r2 && focus2){
nx=x/2/focus2;
ny=y/2/focus2;
nz=-1;
} else {
nx=-sin(phiy2);
ny=0;
nz=-cos(phiy2);
}
}
/* tilt normal vector for roughness, in cone theta_RMS */
if (RMS>0) Lens_roughness(&nx, &ny, &nz, theta_RMS
#ifdef OPENACC
, _particle
#endif
);
/* compute incoming angles w.r.t surface */
NORM(nx,ny,nz);
vec_prod(ax,ay,az,nx,ny,nz,-vx,-vy,-vz);
/* if n and v are parallel - no refraction.*/
if (ax!=0 || ay!=0 || az!=0){
theta1 = atan2(sqrt(ax*ax+ay*ay+az*az),scalar_prod(nx,ny,nz,-vx,-vy,-vz));
/* Fresnel formula for refraction */
theta2=asin(sin(theta1)*n);
/* compute new velocity after refraction */
alpha1 = (sin(theta2))/(sin(theta1));
beta1 = (alpha1*v*cos(theta1))-(v*cos(theta2));
vx = beta1*nx + alpha1*vx;
vy = beta1*ny + alpha1*vy;
vz = beta1*nz + alpha1*vz;
}
} while (dt>0); /* loop until no more intersection */
%}
FINALLY %{
/* compute focal length */
double f1=0, f2=0, f=0;
double n = mean_n / events;
if (n != 1) {
if (r1 && !focus1) f1 = r1/2/(1-n);
else if (!r1 && focus1) f1 = focus1/(1-n);
if (r2 && !focus2) f2 = r2/2/(1-n);
else if (!r2 && focus2) f2 = focus2/(1-n);
}
if (f1) f+= 1/f1;
if (f2) f+= 1/f2;
if (f) {
f = 1/f;
printf("Lens: %s: focal length f=%g [m]. Focal length for N lenses is f/N.\n", NAME_CURRENT_COMP, f);
}
%}
MCDISPLAY
%{
magnify("xy");
double theta1,theta00=0,theta01=0,eps,eps2,theta_line,distance,height,height2,dist_parab1,dist_parab2;
height=radius/2;
height2=radius/2;
dist_parab1=0.0;
dist_parab2=0.0;
/* draw circles to describe the 2 surfaces */
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) { /* OFF file */
off_display(offdata);
}
else {
if (r1 && !focus1) { /* sphere1 */
theta00=asin(fabs((radius/2)/r1));
if (r1<0 && r2<=0 && (-r1+r1*cos(theta00)>thickness/2) ){
theta00=acos((fabs(r1)-thickness/2)/fabs(r1));
height=fabs(r1)*sin(theta00);
}
theta1=-theta00;
eps=theta00/4;
eps2=2*PI/4;
while (theta1<0) {
circle("xy",0,0,-thickness/2-r1+r1*cos(theta1),fabs(r1)*sin(theta1));
theta1+=eps;
theta_line=0;
while (theta_line<2*PI){
line(fabs(r1)*sin(theta1-eps)*cos(theta_line),fabs(r1)*sin(theta1-eps)*sin(theta_line),-thickness/2-r1+r1*cos(theta1-eps),fabs(r1)*sin(theta1)*cos(theta_line),fabs(r1)*sin(theta1)*sin(theta_line),-thickness/2-r1+r1*cos(theta1));
theta_line+=eps2;
}
}
} else if (!r1 && focus1) { /* parabola1 */
if (focus1>0){
dist_parab1=-radius*radius/(4*focus1);
} else {
dist_parab1=0.0;
}
distance=-(radius*radius)/(4*focus1);
eps=-distance/4;
eps2=2*PI/4;
while (focus1*distance<0) {
if (focus1>0){
circle("xy",0,0,distance-thickness/2,sqrt((4*focus1*fabs(distance))));
distance+=eps;
theta_line=0;
while (theta_line<2*PI){
line(sqrt(4*fabs(focus1)*(fabs(distance-eps)))*cos(theta_line),sqrt(4*fabs(focus1)*(fabs(distance-eps)))*sin(theta_line),distance-eps-thickness/2,sqrt(4*fabs(focus1)*fabs(distance))*cos(theta_line),sqrt(4*fabs(focus1)*fabs(distance))*sin(theta_line),distance-thickness/2);
theta_line+=eps2;
}
} else {
theta_line=0;
double rp,rpe,zz;
rp=sqrt(4*fabs(focus1)*fabs(distance));
zz=-thickness/2-radius*radius/(4*fabs(focus1))+distance;
circle("xy",0,0,zz,rp);
distance+=eps;
zz=-thickness/2-radius*radius/(4*fabs(focus1))+distance;
rpe=sqrt(4*fabs(focus1)*fabs(distance));
theta_line=0;
while (theta_line<2*PI){
line(rpe*cos(theta_line),rpe*sin(theta_line),zz, rp*cos(theta_line), rp*sin(theta_line), zz-eps);
theta_line+=eps2;
}
}
}
} else { /* plane 1 */
/* phiy1 is rotation angle around 'y' at z=-thickness. width/height is radius. */
double x1= height2*cos(phiy1);
double y1= height2;
double z1=-height2*sin(phiy1)-thickness/2;
multiline(5, x1, y1, z1,
x1, -y1, z1,
-x1, -y1,-z1,
-x1, y1,-z1,
x1, y1, z1);
}
if (r2 && !focus2) { /* sphere 2 */
theta01=asin(fabs((radius/2)/r2));
if (r1<=0&&r2<0&&(r2-r2*cos(theta01)<-thickness/2)){
theta01=acos((fabs(r2)-thickness/2)/fabs(r2));
height2=fabs(r2)*sin(theta01);
}
theta1=PI-theta01;
eps=(PI-theta1)/4;
eps2=2*PI/4;
while (theta10){
dist_parab2=radius*radius/(4*focus2);
}else{
dist_parab2=0.0;
}
distance=(radius*radius)/(4*focus2);
height2=sqrt((4*focus2*fabs(distance)));
eps=-distance/4;
eps2=2*PI/4;
while (focus2*distance>0){
if (focus2>0) {
circle("xy",0,0,distance+thickness/2,sqrt((4*fabs(focus2)*fabs(distance))));
distance+=eps;
theta_line=0;
while (theta_line<2*PI) {
line(sqrt(4*focus2*(fabs(distance-eps)))*cos(theta_line),sqrt(4*focus2*(fabs(distance-eps)))*sin(theta_line),distance-eps+thickness/2,sqrt(4*focus2*fabs(distance))*cos(theta_line),sqrt(4*focus2*fabs(distance))*sin(theta_line),distance+thickness/2);
theta_line+=eps2;
}
} else {
double rp,rpe,zz;
rp=sqrt(4*fabs(focus2)*fabs(distance));
zz=thickness/2+(radius*radius)/(4*fabs(focus2))+distance;
circle("xy",0,0,zz,rp);
distance+=eps;
zz=thickness/2+(radius*radius)/(4*fabs(focus2))+distance;
rpe=sqrt(4*fabs(focus2)*fabs(distance));
theta_line=0;
while (theta_line<2*PI){
line(rpe*cos(theta_line),rpe*sin(theta_line),zz, rp*cos(theta_line), rp*sin(theta_line), zz-eps);
theta_line+=eps2;
}
}
}
} else { /* plane 2 */
/* phiy2 is rotation angle around 'y' at z=+thickness. width/height is radius. */
double x1= height2*cos(phiy2);
double y1= height2;
double z1=-height2*sin(phiy2)+thickness/2;
multiline(5, x1, y1, z1,
x1, -y1, z1,
-x1, -y1,-z1,
-x1, y1,-z1,
x1, y1, z1);
}
/* draw outer containing cylinder */
if (r1 && r2 && !focus1 && !focus2){
line(-fabs(r1)*sin(theta00),-fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),-fabs(r2)*sin(theta01),fabs(r2)*sin(theta01),thickness/2+r2-r2*cos(theta01));
line(fabs(r1)*sin(theta00),-fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),fabs(r2)*sin(theta01),fabs(r1)*sin(theta00),thickness/2+r2-r2*cos(theta01));
line(-fabs(r1)*sin(theta00),-fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),fabs(r2)*sin(theta01),-fabs(r1)*sin(theta00),thickness/2+r2-r2*cos(theta01));
line(-fabs(r1)*sin(theta00),fabs(r1)*sin(theta00),-thickness/2-r1+r1*cos(theta00),fabs(r2)*sin(theta01),fabs(r1)*sin(theta00),thickness/2+r2-r2*cos(theta01));
}
else if (r1 && !r2) {
if (!focus2){
line(-height,0,-thickness/2-r1+r1*cos(theta00),-height,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(height,0,-thickness/2-r1+r1*cos(theta00),height,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
line(0,-height,-thickness/2-r1+r1*cos(theta00),0,-height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(0,height,-thickness/2-r1+r1*cos(theta00),0,height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
}
else
{
line(-fabs(r1)*sin(theta00),0,-thickness/2-r1+r1*cos(theta00),-radius/2,0,thickness/2+radius*radius/(16*focus2));
line(fabs(r1)*sin(theta00),0,-thickness/2-r1+r1*cos(theta00),radius/2,0,thickness/2+radius*radius/(16*focus2));
}
}
else if (r2 && !r1) {
if (!focus1){
line(-height2,0,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),-height2,0,thickness/2+r2-r2*cos(theta01));
line(height2,0,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),height2,0,thickness/2+r2-r2*cos(theta01));
line(0,height2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),0,height2,thickness/2+r2-r2*cos(theta01));
line(0,height2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),0,height2,thickness/2+r2-r2*cos(theta01));
}
else
{
line(-height2,0,-thickness/2-radius*radius/(16*focus1),-height2,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(height2,0,-thickness/2-radius*radius/(16*focus1),height2,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
line(0,height2,-thickness/2-radius*radius/(16*focus1),0,height2,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(0,height2,-thickness/2-radius*radius/(16*focus1),0,height2,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
}
}
else {
if (!focus1 && !focus2){
line(-radius/2,0,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),-radius/2,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(radius/2,0,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),radius/2,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
line(0,-radius/2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),0,-radius/2,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(0,radius/2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),0,radius/2,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
}
else if (!focus1){
line(-radius/2,0,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),-radius/2,0,thickness/2+dist_parab2);
line(radius/2,0,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),radius/2,0,thickness/2+dist_parab2);
line(0,-radius/2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1),0,-radius/2,thickness/2+dist_parab2);
line(0,radius/2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),0,radius/2,thickness/2+dist_parab2);
}
else if (!focus2){
line(-height,0,-thickness/2+dist_parab1,-height,0,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(height,0,-thickness/2+dist_parab1,height,0,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
line(0,-height,-thickness/2+dist_parab1,0,-height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2));
line(0,height,-thickness/2+dist_parab1,0,height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2));
}
else {
line(-radius,0,-thickness/2+dist_parab1,-radius,0,thickness/2+dist_parab2);
line(radius,0,-thickness/2+dist_parab1,radius,0,thickness/2+dist_parab2);
line(0,-radius,-thickness/2+dist_parab1,0,-radius,thickness/2+dist_parab2);
line(0,radius,-thickness/2+dist_parab1,0,radius,thickness/2+dist_parab2);
}
}
} /* end of not a OFF/PLY */
%}
END