/******************************************************************************* * * 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 * Version: $Revision 1.0$ * Origin: ILL/LLB * Release: McStas 2.0 * * 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) * or parabolic (use focus1 and focus2 as optical focal length). * The lens cross-section is seen as a 2*radius disk from the beam Z axis. * 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. * * Common material: * 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 * perfluoropolymer * * 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 * %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) OUTPUT PARAMETERS (rho, bc, my_s, my_a_v, mean_n, events) STATE PARAMETERS (x,y,z,vx,vy,vz,t,sx,sy,p) SHARE %{ /* Lens_roughness: function to rotate normal vector around axis for roughness * with specified tilt angle (deg) * RETURNS: rotated nx,ny,nz coordinates */ void Lens_roughness(double *nx, double *ny, double *nz, double tilt) { 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 */ 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; int ret = solve_2nd_order(t0,A,B,C); t1 = NULL; return(ret); } /* 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. */ double Lens_intersect(double x,double y,double z, double vx,double vy,double vz, double shift, double radius, char *type, double arg) { double t0, t1; if (!vz || !type || !radius) return -1; if (strstr(type, "spher") && 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 (strstr(type, "parabol") && 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 (strstr(type, "plan")) { /* 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, bc; double my_s, my_a_v; double mean_n; double events; %} 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(sigma_coh*100/4/PI)*1e-5; /* bound coherent cross section */ 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, ny, nz; /* normal vector to surface */ double ax, ay, az; /* temporary vector for surface refraction */ double alpha1, beta1, theta1, theta2; /* angles for refraction computation */ /* ======================== First face of the lens ========================== */ /* determine intersection time */ if (r1 && !focus1) dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, "sphere", r1); else if (!r1 && focus1) dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, "parabola", -focus1); else dt = Lens_intersect(x,y,z, vx, vy, vz, thickness/2, radius, "plane", phiy1); if (dt < 0) ABSORB; /* propagate to surface */ PROP_DT(dt); /* check for lens cross section (radius) */ if (x*x+y*y > radius*radius) ABSORB; 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)); mean_n += n*p; events += p; theta_RMS=atan(2*RMS/lambda); /* cone angle from RMS roughtness */ /* compute normal vector */ 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); /* compute incoming angles w.r.t surface */ NORM(nx,ny,nz); vec_prod(ax,ay,az,nx,ny,nz,-vx,-vy,-vz); 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 */ if (r2 && !focus2) dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, "sphere", -r2); else if (!r2 && focus2) dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, "parabola", focus2); else dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, "plane", phiy2); /* =============== 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 */ p *= p_mult; /* apply absorption correction */ if (dt<0) 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; p *= p_mult; SCATTER; /* recompute new intersection time to go to second surface (velocity changed during scattering event) */ if (r2 && !focus2) dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, "sphere", -r2); else if (!r2 && focus2) dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, "parabola", -focus2); else dt = Lens_intersect(x,y,z, vx, vy, vz, -thickness/2, radius, "plane", phiy2); } /* end scattering handling in material*/ /* propagate to surface */ PROP_DT(dt); /* check for lens cross section (radius) */ if (x*x+y*y > radius*radius) ABSORB; SCATTER; v=sqrt(vx*vx+vy*vy+vz*vz); /* compute normal vector */ 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); /* compute incoming angles w.r.t surface */ NORM(nx,ny,nz); vec_prod(ax,ay,az,nx,ny,nz,-vx,-vy,-vz); 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; %} 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 (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/(16*focus1); } else { dist_parab1=0.0; } distance=-(radius*radius/4)/(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 { circle("xy",0,0,(radius*radius/4)/(4*focus1)+distance-thickness/2,sqrt((4*fabs(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+(radius*radius/4)/(4*focus1)-eps-thickness/2,sqrt(4*fabs(focus1)*fabs(distance))*cos(theta_line),sqrt(4*fabs(focus1)*fabs(distance))*sin(theta_line),distance+(radius*radius/4)/(4*focus1)-thickness/2); theta_line+=eps2; } } } } else { /* plane1 */ line(height2,-height2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),-height2,-height2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1)); line(height2,height2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),-height2,height2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1)); line(-height2,-height2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),-height2,height2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1)); line(height2,-height2,(-thickness/2-radius/2*sin(phiy1))/cos(phiy1),height2,height2,(-thickness/2+radius/2*sin(phiy1))/cos(phiy1)); } 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/(16*focus2); }else{ dist_parab2=0.0; } distance=(radius*radius/4)/(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*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 { circle("xy",0,0,-(radius*radius/4)/(4*focus2)+distance+thickness/2,sqrt((4*fabs(focus2)*fabs(distance)))); distance+=eps; theta_line=0; while (theta_line<2*PI){ line(sqrt(4*fabs(focus2)*(fabs(distance-eps)))*cos(theta_line),sqrt(4*fabs(focus2)*(fabs(distance-eps)))*sin(theta_line),distance-eps+thickness/2-(radius*radius/4)/(4*focus2),sqrt(4*fabs(focus2)*fabs(distance))*cos(theta_line),sqrt(4*fabs(focus2)*fabs(distance))*sin(theta_line),distance+thickness/2-(radius*radius/4)/(4*focus2)); theta_line+=eps2; } } } } else { /* plane 2 */ line(height,-height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2),-height,-height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2)); line(height,height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2),-height,height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2)); line(-height,-height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2),-height,height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2)); line(height,-height,(thickness/2-radius/2*sin(phiy2))/cos(phiy2),height,height,(thickness/2+radius/2*sin(phiy2))/cos(phiy2)); } /* 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/2,0,-thickness/2+dist_parab1,-radius/2,0,thickness/2+dist_parab2); line(radius/2,0,-thickness/2+dist_parab1,radius/2,0,thickness/2+dist_parab2); line(0,-radius/2,-thickness/2+dist_parab1,0,-radius/2,thickness/2+dist_parab2); line(0,radius/2,-thickness/2+dist_parab1,0,radius/2,thickness/2+dist_parab2); } } %} END