#include "invert.h" /********************************************************************/ /* */ /* array access */ /* */ /* We have 4 types of arrays: */ /* Tabulated at midpoint (aM[i] = aM(i)) or end (aE[i] = aE(i+1/2)) */ /* Linear ( aMl or aEl ~ r ) or exponential (aMe or aEe ~ exp(r/H) */ /* */ /* There are 4 types of operations: */ /* Find the midpoint (a[i]), the endpoint (a[i+1/2]), */ /* the start (a[i-1/2]) and the change (a[i+1/2]-a[i-1/2]) */ /* */ /* There are thus 16 functions: */ /* midMl, midMe, midEl, midEe, endMl ..., begMl ..., delMl ... */ /* */ /********************************************************************/ double midMl(double *a, int i) { return a[i]; } double midMe(double *a, int i) { return a[i]; } double midEl(double *a, int i) { return (i>0 ? (a[i] + a[i-1])/2. : a[i] - (a[i+1]-a[i])/2. ); } double midEe(double *a, int i) { return (i>0 ? a[i] * sqrt(a[i-1] / a[i]) : a[i] * sqrt(a[i]/a[i+1]) ); } double endMl(double *a, int i) { return (a[i] + a[i+1])/2.; } double endMe(double *a, int i) { return a[i] * sqrt(a[i+1] / a[i]); } double endEl(double *a, int i) { return a[i]; } double endEe(double *a, int i) { return a[i]; } double begMl(double *a, int i) { return (i>0 ? (a[i] + a[i-1])/2. : a[i] - (a[i+1]-a[i])/2. ); } double begMe(double *a, int i) { return (i>0 ? a[i] * sqrt(a[i-1]/a[i]) : a[i] * sqrt(a[i]/a[i+1]) ); } double begEl(double *a, int i) { return (i>0 ? a[i-1] : a[i] - (a[i+1]-a[i]) ); } double begEe(double *a, int i) { return (i>0 ? a[i-1] : a[i] * (a[i]/a[i+1]) ); } double delMl(double *a, int i) { return ( endMl(a,i) - begMl(a,i) ); } double delMe(double *a, int i) { return ( endMe(a,i) - begMe(a,i) ); } double delEl(double *a, int i) { return ( endEl(a,i) - begEl(a,i) ); } double delEe(double *a, int i) { return ( endEe(a,i) - begEe(a,i) ); } /********************************************************************/ /* */ /* alloc_mat */ /* */ /* Given an array of maxi unassigned pointers to double, allocates */ /* an n-long array of doubles for each. Return (1), or (0) of FAIL */ /* */ /********************************************************************/ int alloc_mat(double **a, int maxi, int n) { int i; for(i = 0; i < maxi; i++) { a[i] = (double *) malloc( n * sizeof (double) ); if (a[i] == (double *) NULL) return (0); } return (1); } int free_mat(double **a, int maxi) { int i; for(i = 0; i < maxi; i++) { free(a[i]); } return (1); } /********************************************************************/ /* */ /* various transidentals */ /* */ /* erf erfc erf2 - error function, complementary erf, general erf */ /* cosh sinh arccosh arcsinh - hyperbolic trig. functions */ /* All definitions are from Abramowitz & Stegun. */ /* */ /********************************************************************/ /* From Abramowitz & Stegun, 7.1.25 WARNING - THIS IS OFF BY 5% AT Y = 5.7 double erfc(double x) { double a1 = 0.3480242; double a2 = -0.0958798; double a3 = 0.7478556; double p = 0.47047; double t = 1.0 / (1.0 + p * x); return ( (a1 + (a2 + a3 *t) *t ) *t ) * exp(-x * x); } double erf(double x) { return (1 - erfc(x)); } double erf2(double xlow, double xhigh) { return( erfc(xlow) - erfc(xhigh) ); } */ /* double cosh(double x) { return ( (exp(x) + exp(-x))/2. ); } double sinh(double x) { return ( (exp(x) - exp(-x))/2. ); } */ double arccosh(double x) { return ( log(x + sqrt(x*x-1)) ); } double arcsinh(double x) { return ( log(x + sqrt(x*x+1)) ); } /* integrals of 2 integral( t^n exp(-t^2), t = y, infinity) */ double yintegral(double y, int n) { double sqrtpi, erfcy, y2, y4, res; switch(n) { case 0: return ( sqrt(PI) * erfc(y) ); case 2: return ( sqrt(PI) * erfc(y)/2. + y * exp(-y*y) ); case 4: return ( sqrt(PI) * erfc(y) * 3./4. + y * exp(-y*y) * (y*y + 3./2.) ); default: return (0.); } } /********************************************************************/ /* */ /* output to outfile */ /* */ /********************************************************************/ void out1(char *fn, char *fmt, int n, double *d) { FILE * fp; int i; fp = fopen(fn, "w"); for (i = 0; i < n; i++) fprintf(fp, fmt, d[i]); fclose(fp); } void out2(char *fn, char *fmt, int n, double *d1, double *d2) { FILE * fp; int i; fp = fopen(fn, "w"); for (i = 0; i < n; i++) fprintf(fp, fmt, d1[i], d2[i]); fclose(fp); } void out3(char *fn, char *fmt, int n, double *d1, double *d2, double *d3) { FILE * fp; int i; fp = fopen(fn, "w"); for (i = 0; i < n; i++) fprintf(fp, fmt, d1[i], d2[i], d3[i]); fclose(fp); } void out4(char *fn, char *fmt, int n, double *d1, double *d2, double *d3, double *d4) { FILE * fp; int i; if (strcmp(fn, "-") == 0) fp = stdout; else fp = fopen(fn, "w"); for (i = 0; i < n; i++) fprintf(fp, fmt, d1[i], d2[i], d3[i], d4[i]); if (strcmp(fn, "-") != 0) fclose(fp); } void out_states(char *fn, double *in[], int n) { FILE * fp; int i, j; char *lbl[] = { "indx\t", "r\t","theta\t", "dtheta\n" }; char *fmt[] = { "%3d\t", "%8.3f\t", "%11.5e\t", "%11.5e\n"}; if (strcmp(fn, "-") == 0) fp = stdout; else fp = fopen(fn, "w"); for (j = 0; j <= 3; j++) fprintf(fp,lbl[j]); for (i = 0; i < n; i++) { j = 0; fprintf(fp,fmt[j], i); j++; /* indx */ fprintf(fp,fmt[j], in[R_I][i]/KM); j++; /* r */ fprintf(fp,fmt[j], in[THETA_I][i]); j++; /* theta */ fprintf(fp,fmt[j], in[DTHETA_I][i]); j++; /* dtheta */ } if (strcmp(fn, "-") != 0) fclose(fp); } /***************************************************************/ /* */ /* dump_states - write the state array to a file */