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

Back to the Galileo-HIC Homepage

Back to the Caltech- Homepage