diff --git a/.gitignore b/.gitignore index 93fc562a..794faf9a 100644 --- a/.gitignore +++ b/.gitignore @@ -29,12 +29,14 @@ libLanlGeoMag/.deps/ *.o *.lo *.png +*.gcno /Makefile Tools/Makefile tests/Makefile libLanlGeoMag/Makefile libLanlGeoMag/Lgm/Makefile libLanlGeoMag/EopData/Makefile +libLanlGeoMag/*.gcno /build-aux/ /autom4te.cache/ configure @@ -63,3 +65,4 @@ Python/Makefile /tests/check_*.log /tests/check_*.trs /tests/test-suite.log +_build* diff --git a/libLanlGeoMag/Lgm/Lgm_QuadPack.h b/libLanlGeoMag/Lgm/Lgm_QuadPack.h index a5b7248c..bf268224 100644 --- a/libLanlGeoMag/Lgm/Lgm_QuadPack.h +++ b/libLanlGeoMag/Lgm/Lgm_QuadPack.h @@ -5,11 +5,11 @@ #include #include -/* - * QuadPack.h +/** + * @file QuadPack.h + * @brief Cified FORTRAN QUADPACK routines */ - #define TRUE 1 #define FALSE 0 #define dmax1(a, b) ( ((a) > (b)) ? (a) : (b) ) @@ -22,38 +22,174 @@ typedef int _qpInfo; double d1mach( int i ); -int dqags(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, - double epsabs, double epsrel, double *result, double *abserr, int *neval, - int *ier, int limit, int lenw, int *last, int *iwork, double *work, int verbosity ); - -int dqagse(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, - double epsabs, double epsrel, int limit, double *result, double *abserr, - int *neval, int *ier, double *alist, double *blist, double *rlist, - double *elist, int *iord, int *last); - - -int dqagp(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, - int npts2, double *points, double epsabs, double epsrel, double *result, - double *abserr, int *neval, int *ier, int leniw, int lenw, int *last, - int *iwork, double *work, int verbosity ); - -int dqagpe(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, - int npts2, double *points, double epsabs, double epsrel, int limit, - double *result, double *abserr, int *neval, int *ier, double *alist, - double *blist, double *rlist, double *elist, double *pts, int *iord, - int *level, int *ndin, int *last); - - -int dqk21(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, - double *result, double *abserr, double *resabs, double *resasc); - -int dqelg(int n, double epstab[], double *result, double *abserr, double res3la[], - int *nres); - -int dqpsrt(int limit, int last, int *maxerr, double *ermax, double elist[], - int iord[], int *nrmax); - -void PrintQuadpackError( int ); +/*! \brief FORTRAN QUADPACK routine DQAGS + * Definite integral, double-precision, general purpose + */ +int dqags(double (*f)( double, _qpInfo *), //!< integrand function f(x, param) + _qpInfo *qpInfo, //!< auxiliary information to pass to function + double a, //!< integral lower limit + double b, //!< integral upper limit + double epsabs, //!< absolute accuracy requested + double epsrel, //!< relative accuracy requested + double *result, //!< calculated integral of f from a to b + double *abserr, //!< estimated modulus of absolute error in result + int *neval, //!< number of function evaluations + int *ier, //!< flag indicating error if positive + int limit, //!< positive number determining maximum #subintervals used + int lenw, //!< must be at least limit*4 + int *last, //!< on return indicates #subintervals produced + int *iwork, //!< must be allocated with size >= limit + double *work, //!< must be allocated with size >= lenw + int verbosity //!< detail level of messaging to stdout + ); + +/*! \brief FORTRAN QUADPACK routine DQAGSE + * Definite integral, double-precision, general purpose, finer control + * and more information returned than DQAGS + */ +int dqagse(double (*f)( double, _qpInfo *), //!< integrand function f(x, param) + _qpInfo *qpInfo, //!< auxiliary information to pass to function + double a, //!< integral lower limit + double b, //!< integral upper limit + double epsabs, //!< absolute accuracy requested + double epsrel, //!< relative accuracy requested + int limit, //!< positive number determining maximum #subintervals used + double *result, //!< calculated integral of f from a to b + double *abserr, //!< estimated modulus of absolute error in result + int *neval, //!< number of function evaluations + int *ier, //!< flag indicating error if positive + double *alist, //!< must be allocated with size >= limit, left endpoints of + //!< subintervals in the partition of (a, b) + double *blist, //!< must be allocated with size >= limit, right endpoints of + //!< subintervals in the partition of (a, b) + double *rlist, //!< must be allocated with size >= limit, approximations of + //!< integral on the subintervals + double *elist, //!< must be allocated with size >= limit, moduli of absolute + //!< error estimates on the subintervals + int *iord, //!< must be allocated with size>= limit, pointers to largest + //!< error estimates + int *last //!< number of subintervals actually produced in the subdivision + ); + +/*! \brief FORTRAN QUADPACK routine DQAGP + * Definite integral, double-precision, singularities or discontinuities + */ +int dqagp(double (*f)( double, _qpInfo *), //!< integrand function f(x, param) + _qpInfo *qpInfo, //!< auxuliary information to pass to function + double a, //!< integral lower limit + double b, //!< integral upper limit + int npts2, //!< number of user-supplied break points + 2 + double *points, //!< user-supplied breakpoints + double epsabs, //!< absolute accuracy requested + double epsrel, //!< relative accuracy requested + double *result, //!< calculated integral of f from a to b + double *abserr, //!< estimated modulus of absolute error + int *neval, //!< number of integrand evaluations + int *ier, //!< flag indicating error if positive + int leniw, //!< dimensioning parameter for iwork + int lenw, //!< dimensioning parameter for work + int *last, //!< on return, #subintervals produced in the subdivision + int *iwork, //!< must be allocated with size >= leniw, pointers to largest + //!< error estimates + double *work, //!< must be allocated with size at least lenw, pointers to + //!< interval endpoints, integral approximations, corresponding + //!< error estimates, and integration limits/break points sorted + int verbosity //!< detail level of messaging to stdout +); + + +/*! \brief FORTRAN QUADPACK routine DQAGPE + * Definite integral, double-precision, singularities or discontinuities + * More control and information returned than DQAGP + */ +int dqagpe(double (*f)( double, _qpInfo *), //!< integrand function f(x, param) + _qpInfo *qpInfo, //!< auxuliary information to pass to function + double a, //!< integral lower limit + double b, //!< integral upper limit + int npts2, //!< number of user-supplied break points + 2 + double *points, //!< user-supplied breakpoints + double epsabs, //!< absolute accuracy requested + double epsrel, //!< relative accuracy requested + int limit, //!< upper bound on number of subintervals of (a, b) to use + double *result, //!< calculated integral of f from a to b + double *abserr, //!< estimated modulus of absolute error + int *neval, //!< number of integrand evaluations + int *ier, //!< flag indicating error if positive + double *alist, //!< must be allocated with size at least limit, + //!< left endpoints of subintervals + double *blist, //!< must be allocated with size at least limit, + //!< right endpoints of subintervals + double *rlist, //!< must be allocated with size at least limit, + //!< integral approximations on the subintervals + double *elist, //!< must be allocated with size at least limit, + //!< moduli of absolute error estimates on subintervals + double *pts, //!< must be allocated with size at least npts2, + //!< integration limits and breakpoints in ascending sequence + int *iord, //!< must be allocated with size at least limit, + //!< pointers to error estimates in decreasing order + int *level, //!< must be allocated with size at least limit, + //!< subdivision levels of subintervals + int *ndin, //!< must be allocated with size at least limit, + //!< indicates whether an interval's error estimate was increased + //!< artificially after first integral step in order to push + //!< a subdivision + int *last //!< number of subintervals actually produced in subdivision + ); + +/*! \brief FORTRAN QUADPACK routine DQK21 + * Definite integral using 21-point Gauss-Kronrod rule, with error + * estimate j = integral of |f| over (a, b) + */ +int dqk21(double (*f)( double, _qpInfo *), //!< integrand function f(x, param) + _qpInfo *qpInfo, //!< auxiliary information to pass to function + double a, //!< integral lower limit + double b, //!< integral upper limit + double *result, //!< calculated integral of f from a to b + double *abserr, //!< estimate of modulus of absolute error + double *resabs, //!< approximation to error estimate j + double *resasc //!< approximation to the integral of |f-i/(b-a)| over (a, b) + ); + +/*! \brief FORTRAN QUADPACK routine DQELG + * estimates limit of a sequence of approximations and error + * Uses P. Wynn's epsilon algorithm + */ +int dqelg(int n, //!< new element in first column of epsilon table + double epstab[], //!< vector of size 52 containing elements of the two lower + //!< diagonals of the triangular epsilon table, numbered + //!< starting at the right-hand corner of the triangle + double *result, //!< resulting approximation + double *abserr, //!< estimate of absolute error from RESULT and the 3 previous + double res3la[], //!< vector of dimension 3 containing the last 3 results + int *nres //!< number of calls to the routine + ); + +/** + * \brief FORTRAN QUADPACK routine DQPSRT + * Maintains the descending ordering in the list of the local error + * estimated resulting from the interval subdivision process. At each + * call two error estimates are inserted using the sequential search + * method, top-down for the largest error estimate and bottom-up for the + * smallest error estimate. + */ +int dqpsrt(int limit, //!< maximum number of error estimates the list can have + int last, //!< number of error estimates currently in the list + int *maxerr, //!< points to NRMAX-th largest error estimate currently in list + double *ermax, //!< NRMAX-th largest error estimate + double elist[], //!< dimension LAST vector containing error estimates + int iord[], //!< dimension LAST vector whose first K elements contain + //!< pointers to error estimates in decreasing sequence + int *nrmax //!< MAXERR = IORD(NRMAX), in accordance with the prophecy + ); + +/** + * \brief PrintQuadpackError Prints plain-English description of ier flag output + * Will print "Unknown error" if given something that isn't a known + * error, even if that something is "no error happened" + */ +void PrintQuadpackError( int ier //!< error flag given by QUADPACK routine + ); #endif + diff --git a/libLanlGeoMag/Lgm/Lgm_SummersDiffCoeff.h b/libLanlGeoMag/Lgm/Lgm_SummersDiffCoeff.h index 911928c3..a0d70003 100644 --- a/libLanlGeoMag/Lgm/Lgm_SummersDiffCoeff.h +++ b/libLanlGeoMag/Lgm/Lgm_SummersDiffCoeff.h @@ -62,7 +62,7 @@ typedef struct Lgm_SummersInfo { double aStarEq; //!< Summer's \f$ \alpha^* \f$ value which is \f$ \Omega_e/\omega^2_{pe} \f$. // double dB; //!< Value of wave amplitude [nT]. void *BwFuncData; //!< Pointer to data that may be needed by BwFunc() - double (*BwFunc)(); //!< Function to return Bw as a function of latitude. + double (*BwFunc)(double, void *); //!< Function to return Bw as a function of latitude. double Omega_eEq; //!< Equatorial gyrofrequency of electrons [Hz]. double Omega_SigEq; //!< Equatorial gyrofrequency of species [Hz]. double w1; //!< Lower frequency cutoff [Hz]. @@ -93,7 +93,7 @@ typedef struct Lgm_SummersInfo { int Lgm_SummersDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)( double, void * ), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *Daa_ba, double *Dap_ba, double *Dpp_ba); //int Lgm_GlauertAndHorneDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)( double, void * ), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, double x1, double x2, int numberOfWaveNormalAngleDistributions, double *xm, double *dx, double *weightsOnWaveNormalAngleDistributions,int WaveMode, int Species, double MaxWaveLat, double *Daa_ba, double *Dap_ba, double *Dpp_ba); int Lgm_GlauertAndHorneDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)( double, void * ), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, double x1, double x2, int numberOfWaveNormalAngleDistributions, double *xm, double *dx, double *weightsOnWaveNormalAngleDistributions,int WaveMode, int Species, double MaxWaveLat, int nNw, int nPlasmaParameters, double aStarMin, double aStarMax, double *Nw, double *Daa_ba, double *Dap_ba, double *Dpp_ba); -int Lgm_SummersDxxDerivsBounceAvg( int DerivScheme, double ha, int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *dDaa, double *dDap); +int Lgm_SummersDxxDerivsBounceAvg( int DerivScheme, double ha, int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(double, void *), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *dDaa, double *dDap); double Lgm_ePlasmaFreq( double Density ); double Lgm_GyroFreq( double q, double B, double m ); double CdipIntegrand_Sb( double Lat, _qpInfo *qpInfo ); diff --git a/libLanlGeoMag/Lgm/praxis.h b/libLanlGeoMag/Lgm/praxis.h new file mode 100644 index 00000000..fc8e960c --- /dev/null +++ b/libLanlGeoMag/Lgm/praxis.h @@ -0,0 +1,43 @@ +#ifndef LGM_PRAXIS_H +#define LGM_PRAXIS_H + +double *allocate_real_vector(int l, int u); +double **allocate_real_matrix(int lr, int ur, int lc, int uc); +void free_real_vector(double *v, int l); +void free_real_matrix(double **m, int lr, int ur, int lc); +void inivec(int l, int u, double *a, double x); +void inimat(int lr, int ur, int lc, int uc, double **a, double x); +void dupvec(int l, int u, int shift, double *a, double *b); +void dupcolvec(int l, int u, int j, double **a, double *b); +void dupmat(int l, int u, int i, int j, double **a, double **b); +double mattam(int l, int u, int i, int j, double **a, double **b); +double matmat(int l, int u, int i, int j, double **a, double **b); +void hshreabid(double **a, int m, int n, double *d, double *b, double *em); +void mulcol(int l, int u, int i, int j, double **a, double **b, double x); +void mulrow(int l, int u, int i, int j, double **a, double **b, double x); +double tammat(int l, int u, int i, int j, double **a, double **b); +double vecvec(int l, int u, int shift, double *a, double *b); +void elmrow(int l, int u, int i, int j, double **a, double **b, double x); +void elmcol(int l, int u, int i, int j, double **a, double **b, double x); +void ichrowcol(int l, int u, int i, int j, double **a); +void elmveccol(int l, int u, int i, double *a, double **b, double x); +int qrisngvaldec(double **a, int m, int n, double *val, double **v, double *em); +int qrisngvaldecbid(double *d, double *b, int m, int n, double **u, double **v, + double *em); +void psttfmmat(double **a, int n, double **v, double *b); +void pretfmmat(double **a, int m, int n, double *d); +void praxismin(int j, int nits, double *d2, double *x1, double *f1, int fk, + int n, double *x, double **v, double *qa, double *qb, double *qc, + double qd0, double qd1, double *q0, double *q1, int *nf, int *nl, + double *fx, double m2, double m4, double dmin, double ldt, + double reltol, double abstol, double small, double h, + double (*funct)(double *, void *data), int *data); +void rotcol(int l, int u, int i, int j, double **a, double c, double s); +double praxisflin(double l, int j, int n, double *x, double **v, double *qa, + double *qb, double *qc, double qd0, double qd1, double *q0, + double *q1, int *nf, double (*f)(double *, void *), + int *data); +void praxis(int n, double *x, int *data, double (*funct)(double *, void *data), + double *in, double *out); + +#endif diff --git a/libLanlGeoMag/Lgm_QuadPack.c b/libLanlGeoMag/Lgm_QuadPack.c index a97a5689..3c2e9f82 100644 --- a/libLanlGeoMag/Lgm_QuadPack.c +++ b/libLanlGeoMag/Lgm_QuadPack.c @@ -3,27 +3,13 @@ /* * QUADPACK DQAGS Routine converted to C */ -int dqags(f, qpInfo, a, b, epsabs, epsrel, result, abserr, neval, ier, limit, lenw, last, iwork, work, verbosity ) -double (*f)( double, _qpInfo *); /* The integrand function -- I.e. the function to integrate */ -_qpInfo *qpInfo; /* Auxilliary information to pass to function (to avoid making globals) */ -double a; /* Lower Limit of integration. */ -double b; /* Upper limit of integration. */ -double epsabs; /* Absolute accuracy requested. */ -double epsrel; /* Relative accuracy requested. */ -double *result; /* The desired result. I.e. integral of f() from a to b */ -double *abserr; /* Estimate of the modulus of the absolute error in the result */ -int *neval; /* The number of integrand evaluations performed. */ -int *ier; /* Error flag. An error occurred if ier > 0. See below. */ -int limit; -int lenw; -int *last; -int *iwork; -double *work; -int verbosity; -{ +int dqags(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, + double epsabs, double epsrel, double *result, double *abserr, int *neval, + int *ier, int limit, int lenw, int *last, int *iwork, double *work, int verbosity ) { + /* * * Begin prologue: dqags @@ -248,34 +234,13 @@ int TMPiwork[600]; } - - - /* * QUADPACK DQAGSE Routine converted to C */ -int dqagse(f, qpInfo, a, b, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) -double (*f)( double, _qpInfo *); /* The integrand function -- I.e. the function to integrate */ -_qpInfo *qpInfo; /* Auxilliary information to pass to function (to avoid making globals) */ -double a; /* Lower Limit of integration. */ -double b; /* Upper limit of integration. */ -double epsabs; /* Absolute accuracy requested. */ -double epsrel; /* Relative accuracy requested. */ -int limit; -double *result; /* The desired result. I.e. integral of f() from a to b */ -double *abserr; /* Estimate of the modulus of the absolute error in the result */ -int *neval; /* The number of integrand evaluations performed. */ -int *ier; /* Error flag. An error occurred if ier > 0. See below. */ -double *alist; -double *blist; -double *rlist; -double *elist; -int *iord; -int *last; -{ - - - +int dqagse(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, + double epsabs, double epsrel, int limit, double *result, double *abserr, + int *neval, int *ier, double *alist, double *blist, double *rlist, + double *elist, int *iord, int *last) { /* * @@ -429,7 +394,7 @@ int *last; double area, abseps, area1, area12, area2, a1; - double a2, b1, b2, correc=0.0, defabs, defab1, defab2, d1mach(); + double a2, b1, b2, correc=0.0, defabs, defab1, defab2; double dres, epmach, erlarg=0.0, erlast, errbnd, errmax; double error1, error2, erro12, errsum, ertest=0.0, oflow, resabs, reseps; double res3la[4], rlist2[53], small=0.0, uflow; @@ -863,10 +828,6 @@ int *last; } - - - - /* * QUADPACK DQELG Routine converted to C */ @@ -928,7 +889,7 @@ int dqelg(int n, double epstab[], double *result, double *abserr, double res3la[ - double delta1, delta2, delta3, d1mach(); + double delta1, delta2, delta3; double epmach, epsinf, error, err1, err2, err3, e0, e1, e1abs, e2, e3; double oflow, res, ss, tol1, tol2, tol3; int i, ib, ib2, ie, indx, k1, k2, k3, limexp, newelm, num; @@ -1125,26 +1086,11 @@ int dqelg(int n, double epstab[], double *result, double *abserr, double res3la[ } - - - - /* * QUADPACK DQK21 Routine converted to C */ -int dqk21(f, qpInfo, a, b, result, abserr, resabs, resasc) -double (*f)( double, _qpInfo *); /* The integrand function -- I.e. the function to integrate */ -_qpInfo *qpInfo; /* Auxilliary information to pass to function (to avoid making globals) */ -double a; /* Lower Limit of integration. */ -double b; /* Upper limit of integration. */ -double *result; /* The desired result. I.e. integral of f() from a to b */ -double *abserr; /* Estimate of the modulus of the absolute error in the result */ -double *resabs; /* */ -double *resasc; /* */ -{ - - - +int dqk21(double (*f)( double, _qpInfo *), _qpInfo* qpInfo, double a, double b, + double* result, double* abserr, double* resabs, double* resasc) { /* * @@ -1200,18 +1146,11 @@ double *resasc; /* */ * end prologue dqk21 */ - - - - double absc, centr, dhlgth, d1mach(); + double absc, centr, dhlgth; double epmach, fc, fsum, fval1, fval2, fv1[11], fv2[11], hlgth; double resg, resk, reskh, uflow; int j, jtw, jtwm1; - - - - /* * the abscissae and weights are given for the interval (-1,1). * because of symmetry only the positive abscissae and their @@ -1233,8 +1172,6 @@ double *resasc; /* */ * bell labs, nov. 1981. */ - - double wg[] = { 0.0, 0.066671344308688137593568809893332, 0.149451349150580593145776339657697, @@ -1268,8 +1205,6 @@ double *resasc; /* */ 0.147739104901338491374841515972068, 0.149445554002916905664936468389821 }; - - /* * * list of major variables @@ -1292,9 +1227,6 @@ double *resasc; /* */ * uflow is the smallest positive magnitude. */ - - - /* * first executable statement dqk21 */ @@ -1305,7 +1237,6 @@ double *resasc; /* */ hlgth = 0.5*(b-a); dhlgth = fabs(hlgth); - /* * compute the 21-point kronrod approximation to * the integral, and estimate the absolute error. @@ -1328,8 +1259,6 @@ double *resasc; /* */ *resabs += wgk[jtw]*(fabs(fval1)+fabs(fval2)); } - - for (j = 1; j<=5; ++j) { jtwm1 = 2*j-1; absc = hlgth*xgk[jtwm1]; @@ -1342,8 +1271,6 @@ double *resasc; /* */ *resabs += wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); } - - reskh = resk*0.5; *resasc = wgk[11]*fabs(fc-reskh); @@ -1358,14 +1285,7 @@ double *resasc; /* */ if ( *resabs > uflow/(50.0*epmach) ) *abserr = dmax1( (epmach*50.0)*(*resabs), *abserr ); - - - return(1); - - - - } diff --git a/libLanlGeoMag/Lgm_QuadPack2.c b/libLanlGeoMag/Lgm_QuadPack2.c index f7b71261..c8d2f347 100644 --- a/libLanlGeoMag/Lgm_QuadPack2.c +++ b/libLanlGeoMag/Lgm_QuadPack2.c @@ -446,7 +446,7 @@ int *last; double area, abseps, area1, area12, area2, a1; - double a2, b1, b2, correc=0.0, defabs, defab1, defab2, d1mach(); + double a2, b1, b2, correc=0.0, defabs, defab1, defab2; double dres, epmach, erlarg=0.0, erlast, errbnd, errmax; double error1, error2, erro12, errsum, ertest=0.0, oflow, resabs, reseps; double res3la[4], rlist2[53], small=0.0, uflow; diff --git a/libLanlGeoMag/Lgm_QuadPack3.c b/libLanlGeoMag/Lgm_QuadPack3.c index 1ea4605b..a7420bc7 100644 --- a/libLanlGeoMag/Lgm_QuadPack3.c +++ b/libLanlGeoMag/Lgm_QuadPack3.c @@ -3,28 +3,13 @@ /* * QUADPACK DQAGP Routine converted to C */ -int dqagp(f, qpInfo, a, b, npts2, points, epsabs, epsrel, result, abserr, neval, ier, leniw, lenw, last, iwork, work, verbosity ) -double (*f)( double, _qpInfo *); /* The integrand function -- I.e. the function to integrate */ -_qpInfo *qpInfo; /* Auxilliary information to pass to function (to avoid making globals) */ -double a; /* Lower Limit of integration. */ -double b; /* Upper limit of integration. */ -int npts2; -double *points; -double epsabs; /* Absolute accuracy requested. */ -double epsrel; /* Relative accuracy requested. */ -double *result; /* The desired result. I.e. integral of f() from a to b */ -double *abserr; /* Estimate of the modulus of the absolute error in the result */ -int *neval; /* The number of integrand evaluations performed. */ -int *ier; /* Error flag. An error occurred if ier > 0. See below. */ -int leniw; -int lenw; -int *last; -int *iwork; -double *work; -int verbosity; -{ +int dqagp(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, + int npts2, double *points, double epsabs, double epsrel, double *result, + double *abserr, int *neval, int *ier, int leniw, int lenw, int *last, + int *iwork, double *work, int verbosity) { + /* BEGIN PROLOGUE DQAGP * PURPOSE The routine calculates an approximation result to a given * definite integral I = Integral of F over (A,B), @@ -292,32 +277,14 @@ int verbosity; /* * QUADPACK DQAGPE Routine converted to C */ -int dqagpe(f, qpInfo, a, b, npts2, points, epsabs, epsrel, limit, result, abserr, neval, ier, alist, blist, rlist, elist, pts, iord, level, ndin, last) -double (*f)( double, _qpInfo *); /* The integrand function -- I.e. the function to integrate */ -_qpInfo *qpInfo; /* Auxilliary information to pass to function (to avoid making globals) */ -double a; /* Lower Limit of integration. */ -double b; /* Upper limit of integration. */ -int npts2; -double *points; -double epsabs; /* Absolute accuracy requested. */ -double epsrel; /* Relative accuracy requested. */ -int limit; -double *result; /* The desired result. I.e. integral of f() from a to b */ -double *abserr; /* Estimate of the modulus of the absolute error in the result */ -int *neval; /* The number of integrand evaluations performed. */ -int *ier; /* Error flag. An error occurred if ier > 0. See below. */ -double *alist; -double *blist; -double *rlist; -double *elist; -double *pts; -int *iord; -int *level; -int *ndin; -int *last; -{ +int dqagpe(double (*f)( double, _qpInfo *), _qpInfo *qpInfo, double a, double b, + int npts2, double *points, double epsabs, double epsrel, int limit, + double *result, double *abserr, int *neval, int *ier, double *alist, + double *blist, double *rlist, double *elist, double *pts, int *iord, + int *level, int *ndin, int *last) { + /* BEGIN PROLOGUE DQAGPE * PURPOSE Approximate a given definite integral I = Integral of F * over (A,B), hopefully satisfying the accuracy claim: @@ -522,7 +489,7 @@ int *last; double area, abseps, area1, area12, area2, a1; - double a2, b1, b2, correc=0.0, defabs, defab1, defab2, d1mach(); + double a2, b1, b2, correc=0.0, defabs, defab1, defab2; double dres, epmach, erlarg=0.0, erlast, errbnd, errmax; double error1, error2, erro12, errsum, ertest=0.0, oflow, resabs, reseps; double res3la[4], rlist2[53], uflow; diff --git a/libLanlGeoMag/Lgm_SummersDiffCoeff.c b/libLanlGeoMag/Lgm_SummersDiffCoeff.c index 99e91ef8..7b937d35 100644 --- a/libLanlGeoMag/Lgm_SummersDiffCoeff.c +++ b/libLanlGeoMag/Lgm_SummersDiffCoeff.c @@ -172,7 +172,7 @@ double Lgm_GyroFreq( double q, double B, double m ) { * \date 2011 * */ -int Lgm_SummersDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *Daa_ba, double *Dap_ba, double *Dpp_ba) { +int Lgm_SummersDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(double, void *), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *Daa_ba, double *Dap_ba, double *Dpp_ba) { double T, a, b, E0, Omega_eEq, Omega_SigEq, Beq, Rho; double epsabs, epsrel, abserr, work[2002], points[10]; @@ -434,7 +434,7 @@ int Lgm_SummersDxxBounceAvg( int Version, double Alpha0, double Ek, double L, * spread in wave normal angles for that component */ -int Lgm_GlauertAndHorneDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, double x1, double x2, int numberOfWaveNormalAngleDistributions, double *xmArray, double *dxArray, double *weightsOnWaveNormalAngleDistributions, int WaveMode, int Species, double MaxWaveLat, int nNw, int nPlasmaParameters, double aStarMin, double aStarMax, double *Nw, double *Daa_ba, double *Dap_ba, double *Dpp_ba) { +int Lgm_GlauertAndHorneDxxBounceAvg( int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(double, void *), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, double x1, double x2, int numberOfWaveNormalAngleDistributions, double *xmArray, double *dxArray, double *weightsOnWaveNormalAngleDistributions, int WaveMode, int Species, double MaxWaveLat, int nNw, int nPlasmaParameters, double aStarMin, double aStarMax, double *Nw, double *Daa_ba, double *Dap_ba, double *Dpp_ba) { double T, a, b, E0, Omega_eEq, Omega_SigEq, Beq, Rho; @@ -684,7 +684,7 @@ printf("done finding integral\n"); * \date 2011 * */ -int Lgm_SummersDxxDerivsBounceAvg( int DerivScheme, double ha, int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *dDaa, double *dDap) { +int Lgm_SummersDxxDerivsBounceAvg( int DerivScheme, double ha, int Version, double Alpha0, double Ek, double L, void *BwFuncData, double (*BwFunc)(double, void *), double n1, double n2, double n3, double aStarEq, int Directions, double w1, double w2, double wm, double dw, int WaveMode, int Species, double MaxWaveLat, double *dDaa, double *dDap) { double a, h, H, faa[7], fap[7], Daa_ba, Dap_ba, Dpp_ba; int i, N; diff --git a/libLanlGeoMag/TA16.c b/libLanlGeoMag/TA16.c index 4edd8c0d..0678d025 100755 --- a/libLanlGeoMag/TA16.c +++ b/libLanlGeoMag/TA16.c @@ -272,7 +272,7 @@ int Lgm_SetCoeffs_TA16(long int Date, double UTC, LgmTA16_Info *ta) { if ( (fp = fopen( datafile, "r" )) != NULL ) { // to start with, just loop over... should we actually // be interpolating linearly between values? - while ( fgets( &tmpstr, 1300, fp ) != NULL ) { + while ( fgets( tmpstr, 1300, fp ) != NULL ) { ncols = sscanf( tmpstr, "%d %d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf %lf", &year, &doy, &hour, &minute, &bx_av, &by_av, &bz_av, &vx, &vy, &vz, &nden, &temp, &symh, &IMFflag, &SWflag, diff --git a/libLanlGeoMag/Tsyg2007.c b/libLanlGeoMag/Tsyg2007.c index 59e7903b..49e422a6 100644 --- a/libLanlGeoMag/Tsyg2007.c +++ b/libLanlGeoMag/Tsyg2007.c @@ -192,15 +192,15 @@ void Lgm_SetCoeffs_TS07( long int Date, double UTC, LgmTsyg2007_Info *t ){ if ( (fp = fopen( Filename, "r" )) != NULL ) { for ( k=1; k<=101; k++ ) { - fgets( &tmpstr, 512, fp); - sscanf( &tmpstr, "%lf", &t->A[k] ); + fgets( tmpstr, 512, fp); + sscanf( tmpstr, "%lf", &t->A[k] ); //fscanf( fp, "%lf%*[\n]\n", &t->A[k] ); //printf("t->A[%d] = %g\n", k, t->A[k]); } while ((!foundP) && (!feof(fp))) { - fgets( &tmpstr, 512, fp); - if ( strstr( &tmpstr, p_str) != NULL ) { //check line for Pdyn, if present read value - sscanf( &tmpstr, "%*s %lf", &t->Pdyn); + fgets( tmpstr, 512, fp); + if ( strstr( tmpstr, p_str) != NULL ) { //check line for Pdyn, if present read value + sscanf( tmpstr, "%*s %lf", &t->Pdyn); foundP = TRUE; } } diff --git a/libLanlGeoMag/praxis.c b/libLanlGeoMag/praxis.c index 0978659f..9454b9da 100644 --- a/libLanlGeoMag/praxis.c +++ b/libLanlGeoMag/praxis.c @@ -3,9 +3,9 @@ * praxis.c -- not sure where the original source of this code is. * */ +#include #include #include -#include // #include // commented out BAL 3Mar2011 for float.h instead #include @@ -19,1021 +19,826 @@ #endif #ifdef MSWIN -#include #include +#include #endif - - - /* * Some defines that may not be known * by all gcc compilers... */ #ifndef RAND_MAX - /* 2^15 - 1 */ +/* 2^15 - 1 */ #define RAND_MAX (32767.0) #endif +#include "Lgm/praxis.h" +void praxis(int n, double *x, int *data, double (*funct)(double *, void *data), + double *in, double *out) { + int illc, i, j, k, k2, nl, maxf, nf, kl, kt, ktm, emergency; + double s, sl, dn, dmin, fx, f1, lds, ldt, sf, df, qf1, qd0, qd1, qa, qb, qc, + m2, m4, small, vsmall, large, vlarge, scbd, ldfac, t2, macheps, reltol, + abstol, h, **v, *d, *y, *z, *q0, *q1, **a, em[8], l; - double *allocate_real_vector(); - double **allocate_real_matrix(); - void free_real_vector(); - void free_real_matrix(); - void inivec(); - void inimat(); - void dupvec(); - void dupmat(); - void dupcolvec(); - void mulrow(); - void mulcol(); - double vecvec(); - double tammat(); - double mattam(); - void ichrowcol(); - void elmveccol(); - int qrisngvaldec(); - void praxismin(); - - - -void praxis( int n, double *x, int *data, double (*funct)(double *, void *data), double *in, double *out) { - - int illc,i,j,k,k2,nl,maxf,nf,kl,kt,ktm,emergency; - double s,sl,dn,dmin,fx,f1,lds,ldt,sf,df,qf1,qd0,qd1,qa,qb,qc,m2,m4, - small,vsmall,large,vlarge,scbd,ldfac,t2,macheps,reltol, - abstol,h,**v,*d,*y,*z,*q0,*q1,**a,em[8],l; - - /* - * Seed random number generator - */ + /* + * Seed random number generator + */ #ifdef MSWIN - srand(34084320); + srand(34084320); #else - srand48(34084320); + srand48(34084320); #endif -// for (i=0; i<8; ++i) x[i+1] = (double)data->x[i]; - d=allocate_real_vector(1,n); - y=allocate_real_vector(1,n); - z=allocate_real_vector(1,n); - q0=allocate_real_vector(1,n); - q1=allocate_real_vector(1,n); - v=allocate_real_matrix(1,n,1,n); - a=allocate_real_matrix(1,n,1,n); - - // heuristic numbers: - // - // If the axes may be badly scaled (which is to be avoided if - // possible), then set scbd = 10. otherwise set scbd=1. - // - // If the problem is known to be ill-conditioned, set ILLC = true. - // - // KTM is the number of iterations without improvement before the - // algorithm terminates. KTM = 4 is very cautious; usually KTM = 1 - // is satisfactory. - // - - macheps=in[0]; - reltol=in[1]; - abstol=in[2]; - maxf=in[5]; - h=in[6]; - scbd=in[7]; - ktm=in[8]; - illc = in[9] < 0.0; - small=macheps*macheps; - vsmall=small*small; - large=1.0/small; - vlarge=1.0/vsmall; - m2=reltol; - m4=sqrt(m2); - srand(1); - ldfac = (illc ? 0.1 : 0.01); - kt=nl=0; - nf=1; - out[3]=qf1=fx=(*funct)(x, data); - abstol=t2=small+fabs(abstol); - dmin=small; - if (h < abstol*100.0) h=abstol*100; - ldt=h; - inimat(1,n,1,n,v,0.0); - for (i=1; i<=n; i++) v[i][i]=1.0; - d[1]=qd0=qd1=0.0; - dupvec(1,n,0,q1,x); - inivec(1,n,q0,0.0); - emergency=0; - - while (1) { - sf=d[1]; - d[1]=s=0.0; - praxismin(1,2,&(d[1]),&s,&fx,0, - n,x,v,&qa,&qb,&qc,qd0,qd1,q0,q1,&nf, - &nl,&fx,m2,m4,dmin,ldt,reltol,abstol,small,h,funct, data); - if (s <= 0.0) mulcol(1,n,1,1,v,v,-1.0); - if (sf <= 0.9*d[1] || 0.9*sf >= d[1]) inivec(2,n,d,0.0); - for (k=2; k<=n; k++) { - dupvec(1,n,0,y,x); - sf=fx; - illc = (illc || kt > 0); - while (1) { - kl=k; - df=0.0; - if (illc) { - /* random stop to get off resulting valley */ - for (i=1; i<=n; i++) { - s=z[i]=(0.1*ldt+t2*pow(10.0,kt))* + // for (i=0; i<8; ++i) x[i+1] = (double)data->x[i]; + d = allocate_real_vector(1, n); + y = allocate_real_vector(1, n); + z = allocate_real_vector(1, n); + q0 = allocate_real_vector(1, n); + q1 = allocate_real_vector(1, n); + v = allocate_real_matrix(1, n, 1, n); + a = allocate_real_matrix(1, n, 1, n); + + // heuristic numbers: + // + // If the axes may be badly scaled (which is to be avoided if + // possible), then set scbd = 10. otherwise set scbd=1. + // + // If the problem is known to be ill-conditioned, set ILLC = true. + // + // KTM is the number of iterations without improvement before the + // algorithm terminates. KTM = 4 is very cautious; usually KTM = 1 + // is satisfactory. + // + + macheps = in[0]; + reltol = in[1]; + abstol = in[2]; + maxf = in[5]; + h = in[6]; + scbd = in[7]; + ktm = in[8]; + illc = in[9] < 0.0; + small = macheps * macheps; + vsmall = small * small; + large = 1.0 / small; + vlarge = 1.0 / vsmall; + m2 = reltol; + m4 = sqrt(m2); + srand(1); + ldfac = (illc ? 0.1 : 0.01); + kt = nl = 0; + nf = 1; + out[3] = qf1 = fx = (*funct)(x, data); + abstol = t2 = small + fabs(abstol); + dmin = small; + if (h < abstol * 100.0) { + h = abstol * 100; + } + ldt = h; + inimat(1, n, 1, n, v, 0.0); + for (i = 1; i <= n; i++) { + v[i][i] = 1.0; + } + d[1] = qd0 = qd1 = 0.0; + dupvec(1, n, 0, q1, x); + inivec(1, n, q0, 0.0); + emergency = 0; + + while (1) { + sf = d[1]; + d[1] = s = 0.0; + praxismin(1, 2, &(d[1]), &s, &fx, 0, n, x, v, &qa, &qb, &qc, qd0, qd1, q0, + q1, &nf, &nl, &fx, m2, m4, dmin, ldt, reltol, abstol, small, h, + funct, data); + if (s <= 0.0) { + mulcol(1, n, 1, 1, v, v, -1.0); + } + if (sf <= 0.9 * d[1] || 0.9 * sf >= d[1]) { + inivec(2, n, d, 0.0); + } + for (k = 2; k <= n; k++) { + dupvec(1, n, 0, y, x); + sf = fx; + illc = (illc || kt > 0); + while (1) { + kl = k; + df = 0.0; + if (illc) { + /* random stop to get off resulting valley */ + for (i = 1; i <= n; i++) { + s = z[i] = (0.1 * ldt + t2 * pow(10.0, kt)) * #ifdef MSWIN - ((double)(rand())/RAND_MAX-0.5); + ((double)(rand()) / RAND_MAX - 0.5); #else - (drand48()-0.5); + (drand48() - 0.5); #endif - elmveccol(1,n,i,x,v,s); - } - fx=(*funct)(x, data); - nf++; - } - for (k2=k; k2<=n; k2++) { - sl=fx; - s=0.0; - praxismin(k2,2,&(d[k2]),&s,&fx,0, - n,x,v,&qa,&qb,&qc,qd0,qd1,q0,q1,&nf, - &nl,&fx,m2,m4,dmin,ldt,reltol,abstol,small,h,funct, data); - s = illc ? d[k2]*(s+z[k2])*(s+z[k2]) : sl-fx; - if (df < s) { - df=s; - kl=k2; - } - } - if (!illc && df < fabs(100.0*macheps*fx)) - illc=1; - else - break; - } - for (k2=1; k2<=k-1; k2++) { - s=0.0; - praxismin(k2,2,&(d[k2]),&s,&fx,0, - n,x,v,&qa,&qb,&qc,qd0,qd1,q0,q1,&nf, - &nl,&fx,m2,m4,dmin,ldt,reltol,abstol,small,h,funct, data); - } - f1=fx; - fx=sf; - lds=0.0; - for (i=1; i<=n; i++) { - sl=x[i]; - x[i]=y[i]; - y[i] = sl -= y[i]; - lds += sl*sl; - } - lds=sqrt(lds); - if (lds > small) { - for (i=kl-1; i>=k; i--) { - for (j=1; j<=n; j++) v[j][i+1]=v[j][i]; - d[i+1]=d[i]; - } - d[k]=0.0; - dupcolvec(1,n,k,v,y); - mulcol(1,n,k,k,v,v,1.0/lds); - praxismin(k,4,&(d[k]),&lds,&f1,1, - n,x,v,&qa,&qb,&qc,qd0,qd1,q0,q1,&nf, - &nl,&fx,m2,m4,dmin,ldt,reltol,abstol,small,h,funct, data); - if (lds <= 0.0) { - lds = -lds; - mulcol(1,n,k,k,v,v,-1.0); - } - } - ldt *= ldfac; - if (ldt < lds) ldt=lds; - t2=m2*sqrt(vecvec(1,n,0,x,x))+abstol; - kt = (ldt > 0.5*t2) ? 0 : kt+1; - if (kt > ktm) { - out[1]=0.0; - emergency=1; - } - } - if (emergency) break; - /* quad */ - s=fx; - fx=qf1; - qf1=s; - qd1=0.0; - for (i=1; i<=n; i++) { - s=x[i]; - x[i]=l=q1[i]; - q1[i]=s; - qd1 += (s-l)*(s-l); - } - l=qd1=sqrt(qd1); - s=0.0; - if ((qd0*qd1 > DBL_MIN) && (nl >=3*n*n)) { - praxismin(0,2,&s,&l,&qf1,1, - n,x,v,&qa,&qb,&qc,qd0,qd1,q0,q1,&nf, - &nl,&fx,m2,m4,dmin,ldt,reltol,abstol,small,h,funct, data); - qa=l*(l-qd1)/(qd0*(qd0+qd1)); - qb=(l+qd0)*(qd1-l)/(qd0*qd1); - qc=l*(l+qd0)/(qd1*(qd0+qd1)); - } else { - fx=qf1; - qa=qb=0.0; - qc=1.0; - } - qd0=qd1; - for (i=1; i<=n; i++) { - s=q0[i]; - q0[i]=x[i]; - x[i]=qa*s+qb*x[i]+qc*q1[i]; - } - /* end of quad */ - dn=0.0; - for (i=1; i<=n; i++) { - d[i]=1.0/sqrt(d[i]); - if (dn < d[i]) dn=d[i]; - } - for (j=1; j<=n; j++) { - s=d[j]/dn; - mulcol(1,n,j,j,v,v,s); - } - if (scbd > 1.0) { - s=vlarge; - for (i=1; i<=n; i++) { - sl=z[i]=sqrt(mattam(1,n,i,i,v,v)); - if (sl < m4) z[i]=m4; - if (s > sl) s=sl; - } - for (i=1; i<=n; i++) { - sl=s/z[i]; - z[i]=1.0/sl; - if (z[i] > scbd) { - sl=1.0/scbd; - z[i]=scbd; - } - mulrow(1,n,i,i,v,v,sl); - } - } - for (i=1; i<=n; i++) ichrowcol(i+1,n,i,i,v); - em[0]=em[2]=macheps; - em[4]=10*n; - em[6]=vsmall; - dupmat(1,n,1,n,a,v); - if (qrisngvaldec(a,n,n,d,v,em) != 0) { - out[1]=2.0; - emergency=1; - } - if (emergency) break; - if (scbd > 1.0) { - for (i=1; i<=n; i++) mulrow(1,n,i,i,v,v,z[i]); - for (i=1; i<=n; i++) { - s=sqrt(tammat(1,n,i,i,v,v)); - d[i] *= s; - s=1.0/s; - mulcol(1,n,i,i,v,v,s); - } - } - for (i=1; i<=n; i++) { - s=dn*d[i]; - d[i] = (s > large) ? vsmall : - ((s < small) ? vlarge : 1.0/(s*s)); - } - /* sort */ - for (i=1; i<=n-1; i++) { - k=i; - s=d[i]; - for (j=i+1; j<=n; j++) - if (d[j] > s) { - k=j; - s=d[j]; - } - if (k > i) { - d[k]=d[i]; - d[i]=s; - for (j=1; j<=n; j++) { - s=v[j][i]; - v[j][i]=v[j][k]; - v[j][k]=s; - } - } - } - /* end of sort */ - dmin=d[n]; - if (dmin < small) dmin=small; - illc = (m2*d[1]) > dmin; - if (nf >= maxf) { - out[1]=1.0; - break; - } - } - out[2]=fx; - out[4]=nf; - out[5]=nl; - out[6]=ldt; - free_real_vector(d,1); - free_real_vector(y,1); - free_real_vector(z,1); - free_real_vector(q0,1); - free_real_vector(q1,1); - free_real_matrix(v,1,n,1); - free_real_matrix(a,1,n,1); - -// for (i=0; i<40; ++i) data->x[i] = (double)x[i+1]; - + elmveccol(1, n, i, x, v, s); + } + fx = (*funct)(x, data); + nf++; + } + for (k2 = k; k2 <= n; k2++) { + sl = fx; + s = 0.0; + praxismin(k2, 2, &(d[k2]), &s, &fx, 0, n, x, v, &qa, &qb, &qc, qd0, + qd1, q0, q1, &nf, &nl, &fx, m2, m4, dmin, ldt, reltol, + abstol, small, h, funct, data); + s = illc ? d[k2] * (s + z[k2]) * (s + z[k2]) : sl - fx; + if (df < s) { + df = s; + kl = k2; + } + } + if (!illc && df < fabs(100.0 * macheps * fx)) { + illc = 1; + } else { + break; + } + } + for (k2 = 1; k2 <= k - 1; k2++) { + s = 0.0; + praxismin(k2, 2, &(d[k2]), &s, &fx, 0, n, x, v, &qa, &qb, &qc, qd0, qd1, + q0, q1, &nf, &nl, &fx, m2, m4, dmin, ldt, reltol, abstol, + small, h, funct, data); + } + f1 = fx; + fx = sf; + lds = 0.0; + for (i = 1; i <= n; i++) { + sl = x[i]; + x[i] = y[i]; + y[i] = sl -= y[i]; + lds += sl * sl; + } + lds = sqrt(lds); + if (lds > small) { + for (i = kl - 1; i >= k; i--) { + for (j = 1; j <= n; j++) { + v[j][i + 1] = v[j][i]; + } + d[i + 1] = d[i]; + } + d[k] = 0.0; + dupcolvec(1, n, k, v, y); + mulcol(1, n, k, k, v, v, 1.0 / lds); + praxismin(k, 4, &(d[k]), &lds, &f1, 1, n, x, v, &qa, &qb, &qc, qd0, qd1, + q0, q1, &nf, &nl, &fx, m2, m4, dmin, ldt, reltol, abstol, + small, h, funct, data); + if (lds <= 0.0) { + lds = -lds; + mulcol(1, n, k, k, v, v, -1.0); + } + } + ldt *= ldfac; + if (ldt < lds) { + ldt = lds; + } + t2 = m2 * sqrt(vecvec(1, n, 0, x, x)) + abstol; + kt = (ldt > 0.5 * t2) ? 0 : kt + 1; + if (kt > ktm) { + out[1] = 0.0; + emergency = 1; + } + } + if (emergency) { + break; + } + /* quad */ + s = fx; + fx = qf1; + qf1 = s; + qd1 = 0.0; + for (i = 1; i <= n; i++) { + s = x[i]; + x[i] = l = q1[i]; + q1[i] = s; + qd1 += (s - l) * (s - l); + } + l = qd1 = sqrt(qd1); + s = 0.0; + if ((qd0 * qd1 > DBL_MIN) && (nl >= 3 * n * n)) { + praxismin(0, 2, &s, &l, &qf1, 1, n, x, v, &qa, &qb, &qc, qd0, qd1, q0, q1, + &nf, &nl, &fx, m2, m4, dmin, ldt, reltol, abstol, small, h, + funct, data); + qa = l * (l - qd1) / (qd0 * (qd0 + qd1)); + qb = (l + qd0) * (qd1 - l) / (qd0 * qd1); + qc = l * (l + qd0) / (qd1 * (qd0 + qd1)); + } else { + fx = qf1; + qa = qb = 0.0; + qc = 1.0; + } + qd0 = qd1; + for (i = 1; i <= n; i++) { + s = q0[i]; + q0[i] = x[i]; + x[i] = qa * s + qb * x[i] + qc * q1[i]; + } + /* end of quad */ + dn = 0.0; + for (i = 1; i <= n; i++) { + d[i] = 1.0 / sqrt(d[i]); + if (dn < d[i]) { + dn = d[i]; + } + } + for (j = 1; j <= n; j++) { + s = d[j] / dn; + mulcol(1, n, j, j, v, v, s); + } + if (scbd > 1.0) { + s = vlarge; + for (i = 1; i <= n; i++) { + sl = z[i] = sqrt(mattam(1, n, i, i, v, v)); + if (sl < m4) { + z[i] = m4; + } + if (s > sl) { + s = sl; + } + } + for (i = 1; i <= n; i++) { + sl = s / z[i]; + z[i] = 1.0 / sl; + if (z[i] > scbd) { + sl = 1.0 / scbd; + z[i] = scbd; + } + mulrow(1, n, i, i, v, v, sl); + } + } + for (i = 1; i <= n; i++) { + ichrowcol(i + 1, n, i, i, v); + } + em[0] = em[2] = macheps; + em[4] = 10 * n; + em[6] = vsmall; + dupmat(1, n, 1, n, a, v); + if (qrisngvaldec(a, n, n, d, v, em) != 0) { + out[1] = 2.0; + emergency = 1; + } + if (emergency) { + break; + } + if (scbd > 1.0) { + for (i = 1; i <= n; i++) { + mulrow(1, n, i, i, v, v, z[i]); + } + for (i = 1; i <= n; i++) { + s = sqrt(tammat(1, n, i, i, v, v)); + d[i] *= s; + s = 1.0 / s; + mulcol(1, n, i, i, v, v, s); + } + } + for (i = 1; i <= n; i++) { + s = dn * d[i]; + d[i] = (s > large) ? vsmall : ((s < small) ? vlarge : 1.0 / (s * s)); + } + /* sort */ + for (i = 1; i <= n - 1; i++) { + k = i; + s = d[i]; + for (j = i + 1; j <= n; j++) { + if (d[j] > s) { + k = j; + s = d[j]; + } + } + if (k > i) { + d[k] = d[i]; + d[i] = s; + for (j = 1; j <= n; j++) { + s = v[j][i]; + v[j][i] = v[j][k]; + v[j][k] = s; + } + } + } + /* end of sort */ + dmin = d[n]; + if (dmin < small) { + dmin = small; + } + illc = (m2 * d[1]) > dmin; + if (nf >= maxf) { + out[1] = 1.0; + break; + } + } + out[2] = fx; + out[4] = nf; + out[5] = nl; + out[6] = ldt; + free_real_vector(d, 1); + free_real_vector(y, 1); + free_real_vector(z, 1); + free_real_vector(q0, 1); + free_real_vector(q1, 1); + free_real_matrix(v, 1, n, 1); + free_real_matrix(a, 1, n, 1); + + // for (i=0; i<40; ++i) data->x[i] = (double)x[i+1]; } -void praxismin(j, nits, d2, x1, f1, fk, n, x, v, qa, qb, qc, qd0, qd1, q0, q1, nf, nl, - fx, m2, m4, dmin, ldt, reltol, abstol, small, h, funct, data) -int j; -int nits; -double *d2; -double *x1; -double *f1; -int fk; -int n; -double x[]; -double **v; -double *qa; -double *qb; -double *qc; -double qd0; -double qd1; -double q0[]; -double q1[]; -int *nf; -int *nl; -double *fx; -double m2; -double m4; -double dmin; -double ldt; -double reltol; -double abstol; -double small; -double h; -double (*funct)(double *, void *data); -int *data; - - - -{ - /* this function is internally used by PRAXIS */ - - double praxisflin(); - int k,dz,loop; - double x2,xm,f0,f2,fm,d1,t2,s,sf1,sx1; - - sf1 = *f1; - sx1 = *x1; - k=0; - xm=0.0; - f0 = fm = *fx; - dz = *d2 < reltol; - s=sqrt(vecvec(1,n,0,x,x)); - t2=m4*sqrt(fabs(*fx)/(dz ? dmin : *d2)+s*ldt)+m2*ldt; - s=s*m4+abstol; - if (dz && (t2 > s)) t2=s; - if (t2 < small) t2=small; - if (t2 > 0.01*h) t2=0.01*h; - if (fk && (*f1 <= fm)) { - xm = *x1; - fm = *f1; - } - if (!fk || (fabs(*x1) < t2)) { - *x1 = (*x1 > 0.0) ? t2 : -t2; - *f1=praxisflin(*x1,j,n,x,v,qa,qb,qc,qd0,qd1,q0,q1,nf,funct, data); - } - if (*f1 <= fm) { - xm = *x1; - fm = *f1; - } - loop=1; - while (loop) { - if (dz) { - /* evaluate praxisflin at another point and - estimate the second derivative */ - x2 = (f0 < *f1) ? -(*x1) : (*x1)*2.0; - f2=praxisflin(x2,j,n,x,v,qa,qb,qc,qd0,qd1,q0,q1,nf,funct, data); - if (f2 <= fm) { - xm=x2; - fm=f2; - } - *d2=(x2*((*f1)-f0)-(*x1)*(f2-f0))/((*x1)*x2*((*x1)-x2)); - } - /* estimate first derivative at 0 */ - d1=((*f1)-f0)/(*x1)-(*x1)*(*d2); - dz=1; - x2 = (*d2 <= small) ? ((d1 < 0.0) ? h : -h) : -0.5*d1/(*d2); - if (fabs(x2) > h) x2 = (x2 > 0.0) ? h : -h; - while (1) { - f2=praxisflin(x2,j,n,x,v,qa,qb,qc,qd0,qd1,q0,q1,nf,funct, data); - if (k < nits && f2 > f0) { - k++; - if (f0 < *f1 && (*x1)*x2 > 0.0) break; - x2=0.5*x2; - } else { - loop=0; - break; - } - } - } - (*nl)++; - if (f2 > fm) - x2=xm; - else - fm=f2; - *d2 = (fabs(x2*(x2-(*x1))) > small) ? - ((x2*((*f1)-f0)-(*x1)*(fm-f0))/((*x1)*x2*((*x1)-x2))) : - ((k > 0) ? 0.0 : *d2); - if (*d2 <= small) *d2=small; - *x1=x2; - *fx=fm; - if (sf1 < *fx) { - *fx=sf1; - *x1=sx1; - } - if (j > 0) elmveccol(1,n,j,x,v,*x1); +void praxismin(int j, int nits, double *d2, double *x1, double *f1, int fk, + int n, double *x, double **v, double *qa, double *qb, double *qc, + double qd0, double qd1, double *q0, double *q1, int *nf, int *nl, + double *fx, double m2, double m4, double dmin, double ldt, + double reltol, double abstol, double small, double h, + double (*funct)(double *, void *data), int *data) { + /* this function is internally used by PRAXIS */ + + int k, dz, loop; + double x2, xm, f0, f2, fm, d1, t2, s, sf1, sx1; + + sf1 = *f1; + sx1 = *x1; + k = 0; + xm = 0.0; + f0 = fm = *fx; + dz = *d2 < reltol; + s = sqrt(vecvec(1, n, 0, x, x)); + t2 = m4 * sqrt(fabs(*fx) / (dz ? dmin : *d2) + s * ldt) + m2 * ldt; + s = s * m4 + abstol; + if (dz && (t2 > s)) { + t2 = s; + } + if (t2 < small) { + t2 = small; + } + if (t2 > 0.01 * h) { + t2 = 0.01 * h; + } + if (fk && (*f1 <= fm)) { + xm = *x1; + fm = *f1; + } + if (!fk || (fabs(*x1) < t2)) { + *x1 = (*x1 > 0.0) ? t2 : -t2; + *f1 = praxisflin(*x1, j, n, x, v, qa, qb, qc, qd0, qd1, q0, q1, nf, funct, + data); + } + if (*f1 <= fm) { + xm = *x1; + fm = *f1; + } + loop = 1; + while (loop) { + if (dz) { + /* evaluat e praxisflin at another point and + estimate the second derivative */ + x2 = (f0 < *f1) ? -(*x1) : (*x1) * 2.0; + f2 = praxisflin(x2, j, n, x, v, qa, qb, qc, qd0, qd1, q0, q1, nf, funct, + data); + if (f2 <= fm) { + xm = x2; + fm = f2; + } + *d2 = + (x2 * ((*f1) - f0) - (*x1) * (f2 - f0)) / ((*x1) * x2 * ((*x1) - x2)); + } + /* estimate first derivative at 0 */ + d1 = ((*f1) - f0) / (*x1) - (*x1) * (*d2); + dz = 1; + x2 = (*d2 <= small) ? ((d1 < 0.0) ? h : -h) : -0.5 * d1 / (*d2); + if (fabs(x2) > h) { + x2 = (x2 > 0.0) ? h : -h; + } + while (1) { + f2 = praxisflin(x2, j, n, x, v, qa, qb, qc, qd0, qd1, q0, q1, nf, funct, + data); + if (k < nits && f2 > f0) { + k++; + if (f0 < *f1 && (*x1) * x2 > 0.0) { + break; + } + x2 = 0.5 * x2; + } else { + loop = 0; + break; + } + } + } + (*nl)++; + if (f2 > fm) { + x2 = xm; + } else { + fm = f2; + } + *d2 = (fabs(x2 * (x2 - (*x1))) > small) + ? ((x2 * ((*f1) - f0) - (*x1) * (fm - f0)) / + ((*x1) * x2 * ((*x1) - x2))) + : ((k > 0) ? 0.0 : *d2); + if (*d2 <= small) { + *d2 = small; + } + *x1 = x2; + *fx = fm; + if (sf1 < *fx) { + *fx = sf1; + *x1 = sx1; + } + if (j > 0) { + elmveccol(1, n, j, x, v, *x1); + } } -double praxisflin(l, j, n, x, v, qa, qb, qc, qd0, qd1, q0, q1, nf, funct, data) - -double l; -int j; -int n; -double x[]; -double **v; -double *qa; -double *qb; -double *qc; -double qd0; -double qd1; -double q0[]; -double q1[]; -int *nf; -double (*funct)(double *, void *); -int *data; +double praxisflin(double l, int j, int n, double *x, double **v, double *qa, + double *qb, double *qc, double qd0, double qd1, double *q0, + double *q1, int *nf, double (*f)(double *, void *), + int *data) { + /* this function is internally used by PRAXISMIN */ -{ - /* this function is internally used by PRAXISMIN */ - - int i; - double *t,result; - - t=allocate_real_vector(1,n); - if (j > 0) - for (i=1; i<=n; i++) t[i]=x[i]+l*v[i][j]; - else { - /* search along parabolic space curve */ - *qa=l*(l-qd1)/(qd0*(qd0+qd1)); - *qb=(l+qd0)*(qd1-l)/(qd0*qd1); - *qc=l*(l+qd0)/(qd1*(qd0+qd1)); - for (i=1; i<=n; i++) t[i]=(*qa)*q0[i]+(*qb)*x[i]+(*qc)*q1[i]; - } - (*nf)++; - result=(*funct)(t, data); - free_real_vector(t,1); - return result; + double* t = allocate_real_vector(1, n); + if (j > 0) { + for (int i = 1; i <= n; i++) { + t[i] = x[i] + l * v[i][j]; + } + } else { + /* search along parabolic space curve */ + *qa = l * (l - qd1) / (qd0 * (qd0 + qd1)); + *qb = (l + qd0) * (qd1 - l) / (qd0 * qd1); + *qc = l * (l + qd0) / (qd1 * (qd0 + qd1)); + for (int i = 1; i <= n; i++) { + t[i] = (*qa) * q0[i] + (*qb) * x[i] + (*qc) * q1[i]; + } + } + (*nf)++; + double result = (*f)(t, data); + free_real_vector(t, 1); + return result; } -void dupcolvec(l, u, j, a, b) -int l; -int u; -int j; -double **a; -double b[]; -{ - for (; l<=u; l++) a[l][j]=b[l]; +void dupcolvec(int l, int u, int j, double **a, double *b) { + for (; l <= u; l++) { + a[l][j] = b[l]; + } } -void dupmat(l, u, i, j, a, b) -int l; -int u; -int i; -int j; -double **a; -double **b; -{ - int k; - - for (; l<=u; l++) - for (k=i; k<=j; k++) a[l][k]=b[l][k]; +void dupmat(int l, int u, int i, int j, double **a, double **b) { + for (; l <= u; l++) { + for (int k = i; k <= j; k++) { + a[l][k] = b[l][k]; + } + } } - - -void dupvec(l, u, shift, a, b) -int l; -int u; -int shift; -double a[]; -double b[]; - -{ - for (; l<=u; l++) a[l]=b[l+shift]; +void dupvec(int l, int u, int shift, double *a, double *b) { + for (; l <= u; l++) { + a[l] = b[l + shift]; + } } - - -void elmveccol(l, u, i, a, b, x) -int l; -int u; -int i; -double a[]; -double **b; -double x; -{ - for (; l<=u; l++) a[l] += b[l][i]*x; +void elmveccol(int l, int u, int i, double *a, double **b, double x) { + for (; l <= u; l++) { + a[l] += b[l][i] * x; + } } - - -void ichrowcol(l, u, i, j, a) -int l; -int u; -int i; -int j; -double **a; -{ - double r; - - for (; l<=u; l++) { - r=a[i][l]; - a[i][l]=a[l][j]; - a[l][j]=r; - } +void ichrowcol(int l, int u, int i, int j, double **a) { + double r; + for (; l <= u; l++) { + r = a[i][l]; + a[i][l] = a[l][j]; + a[l][j] = r; + } } +void inimat(int lr, int ur, int lc, int uc, double **a, double x) { + int j; - -void inimat(lr, ur, lc, uc, a, x) -int lr; -int ur; -int lc; -int uc; -double **a; -double x; -{ - int j; - - for (; lr<=ur; lr++) - for (j=lc; j<=uc; j++) a[lr][j]=x; + for (; lr <= ur; lr++) { + for (j = lc; j <= uc; j++) { + a[lr][j] = x; + } + } } - -void inivec(l, u, a, x) -int l; -int u; -double a[]; -double x; -{ - for (; l<=u; l++) a[l]=x; +void inivec(int l, int u, double *a, double x) { + for (; l <= u; l++) { + a[l] = x; + } } - - - -double mattam(l, u, i, j, a, b) -int l; -int u; -int i; -int j; -double **a; -double **b; -{ - int k; - double s; - - s=0.0; - for (k=l; k<=u; k++) s += a[i][k]*b[j][k]; - return (s); +double mattam(int l, int u, int i, int j, double **a, double **b) { + double s = 0; + for (int k = l; k <= u; k++) { + s += a[i][k] * b[j][k]; + } + return (s); } - - -void mulcol(l, u, i, j, a, b, x) -int l; -int u; -int i; -int j; -double **a; -double **b; -double x; -{ - for (; l<=u; l++) a[l][i]=b[l][j]*x; +void mulcol(int l, int u, int i, int j, double **a, double **b, double x) { + for (; l <= u; l++) { + a[l][i] = b[l][j] * x; + } } - -void mulrow(l, u, i, j, a, b, x) -int l; -int u; -int i; -int j; -double **a; -double **b; -double x; -{ - for (; l<=u; l++) a[i][l]=b[j][l]*x; +void mulrow(int l, int u, int i, int j, double **a, double **b, double x) { + for (; l <= u; l++) { + a[i][l] = b[j][l] * x; + } } - - - -double tammat(l, u, i, j, a, b) -int l; -int u; -int i; -int j; -double **a; -double **b; -{ - int k; - double s; - - s=0.0; - for (k=l; k<=u; k++) s += a[k][i]*b[k][j]; - return (s); +double tammat(int l, int u, int i, int j, double **a, double **b) { + double s = 0; + for (int k = l; k <= u; k++) { + s += a[k][i] * b[k][j]; + } + return (s); } - - -double vecvec(l, u, shift, a, b) -int l; -int u; -int shift; -double a[]; -double b[]; -{ - int k; - double s; - - s=0.0; - for (k=l; k<=u; k++) s += a[k]*b[k+shift]; - return (s); +double vecvec(int l, int u, int shift, double *a, double *b) { + double s = 0; + for (int k = l; k <= u; k++) { + s += a[k] * b[k + shift]; + } + return (s); } - - -int qrisngvaldec(a, m, n, val, v, em) -double **a; -int m; -int n; -double val[]; -double **v; -double em[]; -{ - double *allocate_real_vector(); - void free_real_vector(); - void hshreabid(); - void psttfmmat(); - void pretfmmat(); - int qrisngvaldecbid(); - int i; - double *b; - - b=allocate_real_vector(1,n); - hshreabid(a,m,n,val,b,em); - psttfmmat(a,n,v,b); - pretfmmat(a,m,n,val); - i=qrisngvaldecbid(val,b,m,n,a,v,em); - free_real_vector(b,1); - return i; +int qrisngvaldec(double **a, int m, int n, double *val, double **v, + double *em) { + double *b = allocate_real_vector(1, n); + hshreabid(a, m, n, val, b, em); + psttfmmat(a, n, v, b); + pretfmmat(a, m, n, val); + int i = qrisngvaldecbid(val, b, m, n, a, v, em); + free_real_vector(b, 1); + return i; } - -int qrisngvaldecbid(d, b, m, n, u, v, em) -double d[]; -double b[]; -int m; -int n; -double **u; -double **v; -double em[]; -{ - void rotcol(); - int n0,n1,k,k1,i,i1,count,max,rnk; - double tol,bmax,z,x,y,g,h,f,c,s,min; - - tol=em[2]*em[1]; - count=0; - bmax=0.0; - max=em[4]; - min=em[6]; - rnk=n0=n; - do { - k=n; - n1=n-1; - while (1) { - k--; - if (k <= 0) break; - if (fabs(b[k]) >= tol) { - if (fabs(d[k]) < tol) { - c=0.0; - s=1.0; - for (i=k; i<=n1; i++) { - f=s*b[i]; - b[i] *= c; - i1=i+1; - if (fabs(f) < tol) break; - g=d[i1]; - d[i1]=h=sqrt(f*f+g*g); - c=g/h; - s = -f/h; - rotcol(1,m,k,i1,u,c,s); - } - break; - } - } else { - if (fabs(b[k]) > bmax) bmax=fabs(b[k]); - break; - } - } - if (k == n1) { - if (d[n] < 0.0) { - d[n] = -d[n]; - for (i=1; i<=n0; i++) v[i][n] = -v[i][n]; - } - if (d[n] <= min) rnk--; - n=n1; - } else { - count++; - if (count > max) break; - k1=k+1; - z=d[n]; - x=d[k1]; - y=d[n1]; - g = (n1 == 1) ? 0.0 : b[n1-1]; - h=b[n1]; - f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); - g=sqrt(f*f+1.0); - f=((x-z)*(x+z)+h*(y/((f < 0.0) ? f-g : f+g)-h))/x; - c=s=1.0; - for (i=k1+1; i<=n; i++) { - i1=i-1; - g=b[i1]; - y=d[i]; - h=s*g; - g *= c; - z=sqrt(f*f+h*h); - c=f/z; - s=h/z; - if (i1 != k1) b[i1-1]=z; - f=x*c+g*s; - g=g*c-x*s; - h=y*s; - y *= c; - rotcol(1,n0,i1,i,v,c,s); - d[i1]=z=sqrt(f*f+h*h); - c=f/z; - s=h/z; - f=c*g+s*y; - x=c*y-s*g; - rotcol(1,m,i1,i,u,c,s); - } - b[n1]=f; - d[n]=x; - } - } while (n > 0); - em[3]=bmax; - em[5]=count; - em[7]=rnk; - return n; +int qrisngvaldecbid(double *d, double *b, int m, int n, double **u, double **v, + double *em) { + int n0, n1, k, k1, i, i1, count, max, rnk; + double tol, bmax, z, x, y, g, h, f, c, s, min; + + tol = em[2] * em[1]; + count = 0; + bmax = 0.0; + max = em[4]; + min = em[6]; + rnk = n0 = n; + do { + k = n; + n1 = n - 1; + while (1) { + k--; + if (k <= 0) + break; + if (fabs(b[k]) >= tol) { + if (fabs(d[k]) < tol) { + c = 0.0; + s = 1.0; + for (i = k; i <= n1; i++) { + f = s * b[i]; + b[i] *= c; + i1 = i + 1; + if (fabs(f) < tol) { + break; + } + g = d[i1]; + d[i1] = h = sqrt(f * f + g * g); + c = g / h; + s = -f / h; + rotcol(1, m, k, i1, u, c, s); + } + break; + } + } else { + if (fabs(b[k]) > bmax) { + bmax = fabs(b[k]); + } + break; + } + } + if (k == n1) { + if (d[n] < 0.0) { + d[n] = -d[n]; + for (i = 1; i <= n0; i++) { + v[i][n] = -v[i][n]; + } + } + if (d[n] <= min) { + rnk--; + } + n = n1; + } else { + count++; + if (count > max) { + break; + } + k1 = k + 1; + z = d[n]; + x = d[k1]; + y = d[n1]; + g = (n1 == 1) ? 0.0 : b[n1 - 1]; + h = b[n1]; + f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); + g = sqrt(f * f + 1.0); + f = ((x - z) * (x + z) + h * (y / ((f < 0.0) ? f - g : f + g) - h)) / x; + c = s = 1.0; + for (i = k1 + 1; i <= n; i++) { + i1 = i - 1; + g = b[i1]; + y = d[i]; + h = s * g; + g *= c; + z = sqrt(f * f + h * h); + c = f / z; + s = h / z; + if (i1 != k1) { + b[i1 - 1] = z; + } + f = x * c + g * s; + g = g * c - x * s; + h = y * s; + y *= c; + rotcol(1, n0, i1, i, v, c, s); + d[i1] = z = sqrt(f * f + h * h); + c = f / z; + s = h / z; + f = c * g + s * y; + x = c * y - s * g; + rotcol(1, m, i1, i, u, c, s); + } + b[n1] = f; + d[n] = x; + } + } while (n > 0); + em[3] = bmax; + em[5] = count; + em[7] = rnk; + return n; } -void free_real_vector(v, l) -double *v; -int l; -{ - free((char *) (v+l)); + +void free_real_vector(double *v, int l) { + free((char *)(v + l)); } +void free_real_matrix(double **m, int lr, int ur, int lc) { + for (int i = ur; i >= lr; i--) { + free((char *)(m[i] + lc)); + } + free((char *)(m + lr)); +} +double *allocate_real_vector(int l, int u) { + double* p = (double *)malloc((unsigned)(u - l + 1) * sizeof(double)); + if (!p) { + fprintf(stderr, "Memory allocation failure in allocate_real_vector\n"); + exit(1); + } -void free_real_matrix(m, lr, ur, lc) -double **m; -int lr; -int ur; -int lc; -{ - int i; - for (i=ur; i>=lr; i--) free((char *) (m[i]+lc)); - free((char *) (m+lr)); + return (p - l); } +double **allocate_real_matrix(int lr, int ur, int lc, int uc) { + double **p = (double **)malloc((unsigned)(ur - lr + 1) * sizeof(double *)); + if (!p) { + fprintf(stderr, "Memory allocation failure in allocate_real_matrix\n"); + exit(1); + } + p -= lr; -double *allocate_real_vector(l, u) -int l; -int u; -{ - double *p; - - p = (double *)malloc((unsigned) (u-l+1)*sizeof(double)); + for (int i = lr; i <= ur; i++) { + p[i] = (double *)malloc((unsigned)(uc - lc + 1) * sizeof(double)); if (!p) { - fprintf(stderr, "Memory allocation failure in allocate_real_vector\n"); - exit(1); + fprintf(stderr, "Memory allocation failure in allocate_real_matrix\n"); + exit(1); } + p[i] -= lc; + } - return(p-l); + return p; } +void rotcol(int l, int u, int i, int j, double **a, double c, double s) { + for (; l <= u; l++) { + double x = a[l][i]; + double y = a[l][j]; + a[l][i] = x * c + y * s; + a[l][j] = y * c - x * s; + } +} - -double **allocate_real_matrix(lr, ur, lc, uc) -int lr; -int ur; -int lc; -int uc; -{ - int i; - double **p; - - - p = (double **)malloc((unsigned) (ur-lr+1)*sizeof(double *)); - if (!p) { - fprintf(stderr, "Memory allocation failure in allocate_real_matrix\n"); - exit(1); +void psttfmmat(double **a, int n, double **v, double *b) { + int i, i1, j; + double h; + + i1 = n; + v[n][n] = 1.0; + for (i = n - 1; i >= 1; i--) { + h = b[i] * a[i][i1]; + if (h < 0.0) { + for (j = i1; j <= n; j++) { + v[j][i] = a[i][j] / h; + } + for (j = i1; j <= n; j++) { + elmcol(i1, n, j, i, v, v, matmat(i1, n, i, j, a, v)); + } } - - p -= lr; - - for (i=lr; i<=ur; i++){ - p[i] = (double *)malloc((unsigned) (uc-lc+1)*sizeof(double)); - if (!p) { - fprintf(stderr, "Memory allocation failure in allocate_real_matrix\n"); - exit(1); - } - p[i] -= lc; + for (j = i1; j <= n; j++) { + v[i][j] = v[j][i] = 0.0; } - - return p; -} - - -void rotcol(l, u, i, j, a, c, s) -int l; -int u; -int i; -int j; -double **a; -double c; -double s; -{ - double x, y; - - for (; l<=u; l++) { - x=a[l][i]; - y=a[l][j]; - a[l][i]=x*c+y*s; - a[l][j]=y*c-x*s; - } + v[i][i] = 1.0; + i1 = i; + } } +void pretfmmat(double **a, int m, int n, double *d) { -void psttfmmat(a, n, v, b) -double **a; -int n; -double **v; -double b[]; -{ - double matmat(); - void elmcol(); - int i,i1,j; - double h; - - i1=n; - v[n][n]=1.0; - for (i=n-1; i>=1; i--) { - h=b[i]*a[i][i1]; - if (h < 0.0) { - for (j=i1; j<=n; j++) v[j][i]=a[i][j]/h; - for (j=i1; j<=n; j++) - elmcol(i1,n,j,i,v,v,matmat(i1,n,i,j,a,v)); - } - for (j=i1; j<=n; j++) v[i][j]=v[j][i]=0.0; - v[i][i]=1.0; - i1=i; - } -} -void pretfmmat(a, m, n, d) -double **a; -int m; -int n; -double d[]; -{ - double tammat(); - void elmcol(); - int i,i1,j; - double g,h; - - for (i=n; i>=1; i--) { - i1=i+1; - g=d[i]; - h=g*a[i][i]; - for (j=i1; j<=n; j++) a[i][j]=0.0; - if (h < 0.0) { - for (j=i1; j<=n; j++) - elmcol(i,m,j,i,a,a,tammat(i1,m,i,j,a,a)/h); - for (j=i; j<=m; j++) a[j][i] /= g; - } else - for (j=i; j<=m; j++) a[j][i]=0.0; - a[i][i] += 1.0; - } + for (int i = n; i >= 1; i--) { + int i1 = i + 1; + double g = d[i]; + double h = g * a[i][i]; + for (int j = i1; j <= n; j++) { + a[i][j] = 0.0; + } + if (h < 0.0) { + for (int j = i1; j <= n; j++) { + elmcol(i, m, j, i, a, a, tammat(i1, m, i, j, a, a) / h); + } + for (int j = i; j <= m; j++) { + a[j][i] /= g; + } + } else { + for (int j = i; j <= m; j++) { + a[j][i] = 0.0; + } + } + a[i][i] += 1.0; + } } -void hshreabid(a, m, n, d, b, em) -double **a; -int m; -int n; -double d[]; -double b[]; -double em[]; -{ - double tammat(); - double mattam(); - void elmcol(); - void elmrow(); - int i,j,i1; - double norm,machtol,w,s,f,g,h; - - norm=0.0; - for (i=1; i<=m; i++) { - w=0.0; - for (j=1; j<=n; j++) w += fabs(a[i][j]); - if (w > norm) norm=w; - } - machtol=em[0]*norm; - em[1]=norm; - for (i=1; i<=n; i++) { - i1=i+1; - s=tammat(i1,m,i,i,a,a); - if (s < machtol) - d[i]=a[i][i]; - else { - f=a[i][i]; - s += f*f; - d[i] = g = (f < 0.0) ? sqrt(s) : -sqrt(s); - h=f*g-s; - a[i][i]=f-g; - for (j=i1; j<=n; j++) - elmcol(i,m,j,i,a,a,tammat(i,m,i,j,a,a)/h); - } - if (i < n) { - s=mattam(i1+1,n,i,i,a,a); - if (s < machtol) - b[i]=a[i][i1]; - else { - f=a[i][i1]; - s += f*f; - b[i] = g = (f < 0.0) ? sqrt(s) : -sqrt(s); - h=f*g-s; - a[i][i1]=f-g; - for (j=i1; j<=m; j++) - elmrow(i1,n,j,i,a,a,mattam(i1,n,i,j,a,a)/h); - } - } - } +void hshreabid(double **a, int m, int n, double *d, double *b, double *em) { + double norm = 0; + for (int i = 1; i <= m; i++) { + double w = 0.0; + for (int j = 1; j <= n; j++) { + w += fabs(a[i][j]); + } + if (w > norm) { + norm = w; + } + } + double machtol = em[0] * norm; + em[1] = norm; + for (int i = 1; i <= n; i++) { + int i1 = i + 1; + double s = tammat(i1, m, i, i, a, a); + if (s < machtol) { + d[i] = a[i][i]; + } else { + double f = a[i][i]; + s += f * f; + double g = (f < 0.0) ? sqrt(s) : -sqrt(s); + d[i] = g; + double h = f * g - s; + a[i][i] = f - g; + for (int j = i1; j <= n; j++) { + elmcol(i, m, j, i, a, a, tammat(i, m, i, j, a, a) / h); + } + } + if (i < n) { + s = mattam(i1 + 1, n, i, i, a, a); + if (s < machtol) { + b[i] = a[i][i1]; + } else { + double f = a[i][i1]; + s += f * f; + double g = (f < 0.0) ? sqrt(s) : -sqrt(s); + b[i] = g; + double h = f * g - s; + a[i][i1] = f - g; + for (int j = i1; j <= m; j++) { + elmrow(i1, n, j, i, a, a, mattam(i1, n, i, j, a, a) / h); + } + } + } + } } - - -void elmcol(l, u, i, j, a, b, x) -int l; -int u; -int i; -int j; -double **a; -double **b; -double x; -{ - for (; l<=u; l++) a[l][i] += b[l][j]*x; +void elmcol(int l, int u, int i, int j, double **a, double **b, double x) { + for (; l <= u; l++) { + a[l][i] += b[l][j] * x; + } } - - - -void elmrow(l, u, i, j, a, b, x) -int l; -int u; -int i; -int j; -double **a; -double **b; -double x; -{ - for (; l<=u; l++) a[i][l] += b[j][l]*x; +void elmrow(int l, int u, int i, int j, double **a, double **b, double x) { + for (; l <= u; l++) { + a[i][l] += b[j][l] * x; + } } - - -double matmat(l, u, i, j, a, b) -int l; -int u; -int i; -int j; -double **a; -double **b; -{ - int k; - double s; - - s=0.0; - for (k=l; k<=u; k++) s += a[i][k]*b[k][j]; - return (s); +double matmat(int l, int u, int i, int j, double **a, double **b) { + double s = 0; + for (int k = l; k <= u; k++) { + s += a[i][k] * b[k][j]; + } + return (s); } diff --git a/libLanlGeoMag/xvgifwr2.c b/libLanlGeoMag/xvgifwr2.c index 7ce6c67a..0ad4ac77 100644 --- a/libLanlGeoMag/xvgifwr2.c +++ b/libLanlGeoMag/xvgifwr2.c @@ -11,8 +11,6 @@ * */ - - /***************************************************************** * Portions of this code Copyright (C) 1989 by Michael Mauldin. * Permission is granted to use this file in whole or in @@ -43,10 +41,6 @@ #include #include -/* -#define PARM(a) a -*/ -#define PARM(a) () #define CONV24_8BIT 0 #define CONV24_24BIT 1 #define PIC8 CONV24_8BIT @@ -59,7 +53,7 @@ typedef long int count_int; typedef unsigned long u_long; */ -void xvbzero PARM((char *, size_t)); +void xvbzero(char *, size_t); static int Width, Height; static int curx, cury; @@ -67,20 +61,18 @@ static long CountDown; static int Interlace; /* static byte bw[2] = {0, 0xff}; */ -static void putword PARM((int, FILE *)); -static void compress PARM((int, FILE *, unsigned char *, int)); -static void output PARM((int)); -static void cl_block PARM((void)); -static void cl_hash PARM((count_int)); -static void char_init PARM((void)); -static void char_out PARM((int)); -static void flush_char PARM((void)); +static void putword(int, FILE *); +static void compress(int, FILE *, unsigned char *, int); +static void output(int); +static void cl_block(void); +static void cl_hash(count_int); +static void char_init(void); +static void char_out(int); +static void flush_char(void); static unsigned char pc2nc[256],r1[256],g1[256],b1[256]; - - /* * Added by MGH Jan 29, 2006 */ @@ -108,26 +100,14 @@ void SwapIntBytes( int *a ){ memcpy( a, c, si); } - - - - - - /*************************************************************/ -int WriteGIF(fp, pic, ptype, w, h, rmap, gmap, bmap, numcols, colorstyle, - comment) - FILE *fp; - unsigned char *pic; - int ptype, w,h; - unsigned char *rmap, *gmap, *bmap; - int numcols, colorstyle; - char *comment; -{ +int WriteGIF(FILE *fp, unsigned char *pic, int ptype, int w, int h, + unsigned char *rmap, unsigned char *gmap, unsigned char *bmap, + int numcols, int colorstyle, char *comment) { // int DEBUG=1; int RWidth, RHeight; int LeftOfs, TopOfs; - int ColorMapSize, InitCodeSize, Background, BitsPerPixel; + int ColorMapSize, Background, BitsPerPixel; int i,j,nc; unsigned char *pic8; @@ -141,18 +121,15 @@ int WriteGIF(fp, pic, ptype, w, h, rmap, gmap, bmap, numcols, colorstyle, } else pic8 = pic; - /* * If we are bigendian, we need to swap order of multi-byte words * Added by MGH Jan 29, 2006. */ SwapBytes = BigEndian(); - Interlace = 0; Background = 0; - for (i=0; i<256; i++) { pc2nc[i] = r1[i] = g1[i] = b1[i] = 0; } /* compute number of unique colors */ @@ -190,8 +167,7 @@ int WriteGIF(fp, pic, ptype, w, h, rmap, gmap, bmap, numcols, colorstyle, CountDown = w * h; /* # of pixels we'll be doing */ - if (BitsPerPixel <= 1) InitCodeSize = 2; - else InitCodeSize = BitsPerPixel; + int InitCodeSize = (BitsPerPixel <= 1) ? 2 : BitsPerPixel; curx = cury = 0; @@ -235,7 +211,6 @@ int WriteGIF(fp, pic, ptype, w, h, rmap, gmap, bmap, numcols, colorstyle, fputc(0, fp); /* future expansion byte */ - if (colorstyle == 1) { /* greyscale */ for (i=0; i>8)&0xff, fp); } - - - /***********************************************************************/ - static unsigned long cur_accum = 0; static int cur_bits = 0; - - - #define min(a,b) ((a>b) ? b : a) #define XV_BITS 12 /* BITS was already defined on some systems */ @@ -327,7 +288,6 @@ static int cur_bits = 0; typedef unsigned char char_type; - static int n_bits; /* number of bits/code */ static int maxbits = XV_BITS; /* user settable max # bits/code */ static int maxcode; /* maximum code, given n_bits */ @@ -390,12 +350,8 @@ static int EOFCode; /********************************************************/ -static void compress(init_bits, outfile, data, len) -int init_bits; -FILE *outfile; -unsigned char *data; -int len; -{ +static void compress(int init_bits, FILE* outfile, unsigned char* data, + int len) { register long fcode; register int i = 0; register int c; @@ -524,9 +480,7 @@ unsigned long masks[] = { 0x0000, 0x0001, 0x0003, 0x0007, 0x000F, 0x01FF, 0x03FF, 0x07FF, 0x0FFF, 0x1FFF, 0x3FFF, 0x7FFF, 0xFFFF }; -static void output(code) -int code; -{ +static void output(int code) { cur_accum &= masks[cur_bits]; if (cur_bits > 0) @@ -583,8 +537,7 @@ int code; /********************************/ -static void cl_block () /* table clear for block compress */ -{ +static void cl_block () { /* table clear for block compress */ /* Clear out the hash table */ cl_hash ( (count_int) hsize ); @@ -596,9 +549,7 @@ static void cl_block () /* table clear for block compress */ /********************************/ -static void cl_hash(hsize) /* reset code table */ -register count_int hsize; -{ +static void cl_hash(register count_int hsize) { /* reset code table */ register count_int *htab_p = htab+hsize; register long i; register long m1 = -1; @@ -657,9 +608,7 @@ static char accum[ 256 ]; * Add a character to the end of the current packet, and if it is 254 * characters, flush the packet to disk. */ -static void char_out(c) -int c; -{ +static void char_out(int c) { accum[ a_count++ ] = c; if( a_count >= 254 ) flush_char(); @@ -668,17 +617,14 @@ int c; /* * Flush the packet to disk, and reset the accumulator */ -static void flush_char() -{ +static void flush_char() { if( a_count > 0 ) { fputc(a_count, g_outfile ); fwrite(accum, (size_t) 1, (size_t) a_count, g_outfile ); a_count = 0; } } -void xvbzero(s, len) - char *s; - size_t len; -{ + +void xvbzero(char* s, size_t len) { for ( ; len>0; len--) *s++ = 0; }