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_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

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

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;

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;