PR #67 amended to retain praxis

.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*

libLanlGeoMag/Lgm/Lgm_QuadPack.h

@@ -5,11 +5,11 @@
 #include <stdio.h>
 #include <math.h>
 
-/*
- *    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 
 
+

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 );

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