Source Code to wrcspni.c
Source Code to wrcspni.c
The following listing gives the source code to the C function wrcspni.c. This page is currently undergoing construction. Please check back for more detail soon.
/* Uses Rick Cook's range-extension parameters and formula 29 January, 1982 */
/* Fixed problem in zcaldl 11/9/84 */
/* janni table and functions to calculate r, e, z, m 25 August, 1980 */
#define MAXPTS 140
#define MICRON 10000.0 /* range units per cm. 10000 for microns */
#define MCRN 0.1 /* === MICRON * 1.0E-5 */
static float janni[MAXPTS][2] = {
{0.10,0.000122},{0.15,0.000174},{0.20,0.000233},
{0.30,0.000364},{0.40,0.000513},{0.50,0.000678},
{0.60,0.000860},{0.70,0.001059},{0.80,0.001273},
{0.90,0.001500},{1.00,0.001737},{1.20,0.002250},
{1.40,0.002820},{1.60,0.003445},{1.80,0.004125},
{2.00,0.004854},{2.20,0.005637},{2.40,0.006471},
{2.60,0.007355},{2.80,0.008287},{3.00,0.009267},
{3.20,0.010296},{3.40,0.011371},{3.60,0.012494},
{3.80,0.013664},{4.00,0.014877},{4.20,0.016137},
{4.40,0.017441},{4.60,0.018788},{4.80,0.020183},
{5.00,0.02162},{5.50,0.02540},{6.00,0.02944},
{6.50,0.03374},{7.00,0.03831},{7.50,0.04312},
{8.00,0.04817},{8.50,0.05347},{9.00,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.00,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;
#include
#include
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 ;
set3spl(MAXPTS,rang,enrg,spx,sppx);
spy = janni[0];
sppy = janni[MAXPTS/2];
set3spl(MAXPTS,enrg,rang,spy,sppy);
apm[8] = apm[1] * apm[2];
}
double rofej(e)
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)
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)
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)
double e,z,m;
{
double epn,u,rcor,rnge,rofej(),cfun(),sqrt(),pow();
epn = e/m;
if (epn>=0.) u = 137.0*sqrt(epn/469.10)/z;
else { /*printf ("err%lf",u);*/ u = 0.0; }
rcor = cfun(u)*m*pow(z,0.666666666667)*MCRN;
rnge = rofej(epn)*m/(z*z) + rcor;
return(rnge);
}
double eorc(r,z,m)
double r,z,m;
{
double q,e1,e2,r1,r2,eofrj(),roec(),fabs();
double efr;
int i;
/*printf ("eorc of %lf",r);*/
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)
double z,de,e,t;
{
double fabs(),roec(),pow();
double q,f1,f2,m1,m2,ei;
double 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;
/*
fprintf(stderr,"\n e3,z1,m1 = %10.2f %10.2f %10.2f\n",e3,z1,m1);
*/
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;
/*
fprintf(stderr,"\n%10.2lf %10.2lf %16.7le\n",z1,z2,f2);
*/
for (i=0;i 10000.) {
return(-3.);
}
f1=f2;
/*
fprintf(stderr,"\n e3,z2,m2 = %10.2f %10.2f %10.2f\n",e3,z2,m2);
*/
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;
/*
fprintf(stderr,"\n%10.2lf %10.2lf %16.7le\n",z1,z2,f2);
*/
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;
}
/*
fprintf(stderr,"\n%10.2f %5.0f %5.0f\n",z2,e1,e3);
*/
return(-z2);
}
Last Updated: Monday 28 April, 1997