Source Code for rngfun.c
Source Code for rngfun.c
The following gives the source code listing for rngfun.c. This page is currently under construction. Please check back soon for more detail.
/* 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; i0.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;i21.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; i21.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;i0;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: Monday 28 April, 1997