/* Uses Rick Cook's range-extension. See WRC thesis. * janni table basis */ #include#include #define MAXPTS 140 #define MICRON 10000.0 /* microns per cm (Janni in cm Si */ #define MCRN 0.1 /* === MICRON * 1.0E-5 */ static float janni[MAXPTS][2] = { {0.1,0.000122}, {0.15,.000174}, {0.2,0.000233}, {0.3,0.000364}, {0.4,0.000513}, {0.5,0.000678}, {0.6,0.000860}, {0.7,0.001059}, {0.8,0.001273}, {0.9,0.001500}, {1.0,0.001737}, {1.2,0.002250}, {1.4,0.002820}, {1.6,0.003445}, {1.8,0.004125}, {2.0,0.004854}, {2.2,0.005637}, {2.4,0.006471}, {2.6,0.007355}, {2.8,0.008287}, {3.0,0.009267}, {3.2,0.010296}, {3.4,0.011371}, {3.6,0.012494}, {3.8,0.013664}, {4.0,0.014877}, {4.2,0.016137}, {4.4,0.017441}, {4.6,0.018788}, {4.8,0.020183}, {5.00,0.02162}, {5.50,0.02540}, {6.0,0.02944}, {6.50,0.03374}, {7.00,0.03831}, {7.5,0.04312}, {8.00,0.04817}, {8.50,0.05347}, {9.0,0.05902}, {9.50,0.06482}, {10.0,0.07084}, {11.0,0.08360}, {12.0,0.09730}, {13.0,0.11190}, {14.0,0.12741}, {15.0,0.14380}, {16.0,0.16107}, {17.0,0.17921}, {18.0,0.19820}, {19.0,0.21803}, {20.0,0.23871}, {22.0,0.28253}, {24.0,0.32962}, {26.0,0.37989}, {28.0,0.43332}, {30.0,0.48982}, {32.0,0.54936}, {34.0,0.61191}, {36.0,0.67742}, {38.0,0.74583}, {40.0,0.81713}, {45.0,1.00770}, {50.0,1.21560}, {55.0,1.4403}, {60.0,1.6813}, {65.0,1.9382}, {70.0,2.2107}, {75.0,2.4983}, {80.0,2.8008}, {90.0,3.4489}, {100.,4.1527}, {110.,4.9097}, {120.,5.7180}, {130.,6.5756}, {140.,7.4806}, {150.,8.4314}, {160.,9.4263}, {170.,10.464}, {180.,11.542}, {190.,12.661}, {200.,13.817}, {210.,15.011}, {220.,16.241}, {230.,17.505}, {240.,18.804}, {250.,20.135}, {260.,21.497}, {270.,22.890}, {280.,24.313}, {290.,25.765}, {300.,27.245}, {310.,28.752}, {320.,30.285}, {330.,31.844}, {340.,33.427}, {350.,35.035}, {360.,36.666}, {370.,38.320}, {380.,39.997}, {390.,41.694}, {400.,43.413}, {410.,45.152}, {420.,46.911}, {430.,48.689}, {440.,50.485}, {450.,52.300}, {460.,54.133}, {470.,55.982}, {480.,57.850}, {490.,59.734}, {500.,61.633}, {510.,63.548}, {520.,65.478}, {530.,67.423}, {540.,69.383}, {550.,71.356}, {560.,73.343}, {570.,75.344}, {580.,77.358}, {590.,79.384}, {600.,81.423}, {620.,85.536}, {640.,89.696}, {660.,93.899}, {680.,98.145}, {700.,102.43}, {720.,106.75}, {740.,111.11}, {760.,115.51}, {780.,119.93}, {800.,124.39}, {820.,128.88}, {840.,133.40}, {860.,137.95}, {880.,142.52}, {900.,147.12}, {920.,151.74}, {940.,156.39}, {960.,161.07}, {1000.,170.52} }; static double apm[10] = {3.2209,2.0040,0.92383,1.1242,0.28267,0.98002,0.0,0.0,0.0,0.0}; static float error,enrg[MAXPTS],rang[MAXPTS],spx[MAXPTS],sppx[MAXPTS]; static float *spy,*sppy; static int numits; setre(err,num) /* takes logs of e, r and converts range to microns */ int num; double err; { int i; double log(),exp(); for (i=0; i 0.0) ? err : 1.0e-4 ; numits = (num>0 ) ? num : 10 ; spy = janni[0]; sppy = janni[MAXPTS/2]; set3spl(MAXPTS,rang,enrg,spx,sppx); set3spl(MAXPTS,enrg,rang,spy,sppy); apm[8] = apm[1] * apm[2]; } double rofej(e) /* energy to range based on Janni */ double e; { double loge,logr,log(),exp(),get3spl(); if (e < 0.1) e = 0.1; loge=log(e); logr=get3spl(loge,MAXPTS,rang,enrg,spx,sppx); return(exp(logr)); } double eofrj(r) /* range to energy based on Janni */ double r; { double loge,logr,log(),exp(),get3spl(); if (r < 1.0) r = 1.0; logr=log(r); loge=get3spl(logr,MAXPTS,enrg,rang,spy,sppy); return(exp(loge)); } double cfun(u) /* correction function */ double u; { double exp(),log(),sqrt(); if (u > 3.0) u = 3.0; return(apm[8]*(exp(-u*(1.0/apm[2]+u/apm[5]))-1.0) + apm[0]*(apm[3]-apm[4]*log(1.0+exp((apm[3]-u)/apm[4])))); } double roec(e,z,m) /* energy to range, corrected */ double e,z,m; { double epn,u,rcor,rnge,rofej(),cfun(),sqrt(),pow(); epn = e/m; u = 137.0*sqrt(epn/469.10)/z; rcor = cfun(u)*m*pow(z,0.666666666667)*MCRN; rnge = rofej(epn)*m/(z*z) + rcor; return(rnge); } double eorc(r,z,m) /* range to energy, corrected */ double r,z,m; { double efr,q,e1,e2,r1,r2,eofrj(),roec(),fabs(); int i; efr = error/1000.0; e1 = eofrj(r*z*z/m)*m; e2 = e1 * 0.90; r2 = r - roec(e1,z,m); for (i=0;i 21.0) m1 = 2.132 * z1; else m1 = 40.0 + (z1 - 20.0) * 4.772; if (z2<20.0) m2 = 2.0 * z2; else if (z2>21.0) m2 = 2.132 * z2; else m2 = 40.0 + (z2 - 20.0) * 4.772; f2 = t - roec(ei,z1,m1) + roec(e,z1,m1); for (i=0; i 21.0) m2 = 2.132 * z2; else m2 = 40.0 + (z2 - 20.0) * 4.772; } return(-z2); } double mcal(z,de,e,t) /* mass calculation */ double z,de,e,t; { double fabs(),roec(),pow(); double q,f1,f2,m1,m2,ei,efr; int i; if (z<0.1 ||de<0.01 ||e<0.01 ||t<1.0) return(-1.0); efr = error/100.0; ei = e + de; m1 = 2.0 * z; m2 = pow(ei,1.77) - pow(e,1.77); m2 = pow(t*z*z/(m2*11.9),-1.3); f2 = t - roec(ei,z,m1) + roec(e,z,m1); for (i=0; i = xmax. Dave Chenette 25 February, 1981 */ set3spl(n,y,x,sp,spp) int n; float *y,*x,*sp,*spp; { register int i; double a,b,c,d,e,f; *sp = 0.0; *spp = 0.0; for (i=1;i 0;i--) *(spp+i) = *(sp+i) * *(spp+i+1) + *(spp+i); *(spp) = 0.0; for (i=0;i = *(x+ix)) return(*(y+ix) + *(sp+ix) * (arg - *(x+ix))); in = 0; while ((ix-in) > 1) { i = (ix+in+1)/2; if (arg < *(x+i)) ix = i; else if (arg > *(x+i)) in = i; else return(*(y+i)); } i = in; a = arg - *(x+i); b = (*(spp+i+1) - *(spp+i))/(*(x+i+1) - *(x+i)); a = *(y+i)+a*(*(sp+i)+a*(*(spp+i)/2.0+a*b/6.0)); return(a); } /***********************************************************/ double zcaldl(e1,e3,t1,t2,e2) double e1, e3, t1, t2, *e2; { double eorc(), roec(),fabs(),pow(); double q,f1,f2,z1,z2,m1,m2,ei,efr, rng3; int i; if (e1 < .01 || e3 < .01 || t1 < 1. || t2 < 1.) return(-1.); efr=error/100.; ei=e3+2.*e1; m1=pow(ei,1.77)-pow(e3+e1,1.77); z1=pow(6.98*m1/t1,.36); z2=1.2*z1; if (z1 < 20.) m1=2.*z1; else if (z1 > 21.) m1=2.132*z1; else m1=40.+(z1-20.)*4.772; if (z2 < 20.) m2=2.*z2; else if (z2 > 21.) m2=2.132*z2; else m2=40.+(z2-20.)*4.772; rng3 = roec(e3,z1,m1); *e2 = eorc(t1+t2+rng3,z1,m1) - e1 - e3; f2 = (e1 - eorc(t1+t2+rng3,z1,m1) + eorc(t2+rng3,z1,m1))/e1; for (i=0;i 10000.) { return(-3.); } f1=f2; rng3 = roec(e3,z2,m2); *e2 = eorc(t1+t2+rng3,z2,m2) - e1 - e3; f2 = (e1 - eorc(t1+t2+rng3,z2,m2) + eorc(t2+rng3,z2,m2))/e1; if (fabs(f2) < .00003) return(z2); if (fabs(q = 1. - f1/f2) < efr) return(-2.); if (fabs(q = (z2 - z1)/q) < error) return(z2); z1=z2; m1=m2; z2 -= q; if (z2 < 20.) m2=2.*z2; else if (z2 > 21.) m2=2.132*z2; else m2=40.+(z2-20.)*4.772; } return(-z2); }

Last Updated: