##// END OF EJS Templates
DATABASES » madrigal
RSS
Test fix upload
Test fix upload
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r16:a7d6f530f4cf
Show More
updateMadDerivations.py
1807 lines | 54.5 KiB | text/x-python | PythonLexer
/ source / madpy / scripts / bin / updateMadDerivations.py
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 #!$MADROOT/bin/python
"""updateMadDerivations.py is a script that autogenerates ../../madrigal/_derive_new.c
based on the latest methods in source/madpy/madrigal/derivation.py.
Run after derivation.py has been updated. Must be run from madpy/scripts/bin directory
$Id: updateMadDerivations.py 7046 2019-10-07 19:57:14Z brideout $
"""
# standard python imports
import os
# madrigal imports
import madrigal.derivation
# verify we are in the right directory
pwd = os.getcwd()
if not pwd.endswith('source/madpy/scripts/bin'):
raise IOError('This script must be run from source/madpy/scripts/bin, not %s' % (pwd))
codeTemplate1 = """ if (strcmp(methodName, "%s")==0)
%s(inLen, inValues, outLen, outValues, stderr);
"""
codeTemplate2 = """ else if (strcmp(methodName, "%s")==0)
%s(inLen, inValues, outLen, outValues, stderr);
"""
# first template - start of code
startTemplate = """/*
* This code exists to expose all the calculation methods in madDeriveMethods via
* a single python call: madDispatch. madDispatch takes three arguments:
* 1. a string name of the function to call
* 2. a list of floats as inputs to the call
* 3. a list of floats to be set by the call
*
* This file is autogenerated by updateMadDerivations.py. Please
* do not directly edit this file. updateMadDerivations.py takes
* as input the derivation.py file, and creates C code that
* matches all calls in it. The only files that are manually edited
* when a new method is added are derivation.py and
* madDeriveMethods.[c,h]
*
* This file also contains some help methods to expose fuctionality of
* the geolib. However, again, this file should not be directly modified
* in Subversion - instead, the new helper code must be added to
* updateMadDerivations.py, and this file reautogenerated.
*
* $Id: updateMadDerivations.py 7046 2019-10-07 19:57:14Z brideout $
*/
#include <Python.h>
#include <numpy/ndarrayobject.h>
#include <madDeriveMethods.h>
#include <date.h>
/***********************************************************************
*
* madDispatch runs an underlying C derivation method
*
* arguments:
* 1. Python string representing name of method to call
* 2. Python double 1D numpy vector representing the input arguments
* 3. Python double 1D numpy vector representing the output - will be modified
*
* returns:
* Python int 0
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_madDispatch(PyObject * self, PyObject * args)
{
/* local variables */
char errBuf[1000] = "";
int inLen = 0, outLen = 0;
double * inValues;
double * outValues;
/* input variables */
char * methodName = NULL;
PyArrayObject *inArr;
PyArrayObject *outArr;
if (!PyArg_ParseTuple(args, "sO!O!",
&methodName,
&PyArray_Type, &inArr,
&PyArray_Type, &outArr))
return NULL;
/* check that one-D arrays */
if (inArr->nd != 1 || outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions in in and out arrays must be 1, not %i and %i",
inArr->nd, outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
inLen = inArr->dimensions[0];
outLen = outArr->dimensions[0];
inValues = (double *)inArr->data;
outValues = (double *)outArr->data;
/* call the right C method - the following code is autogenerated */
"""
# end template
endTemplate = """ /* end autogenerated code */
else
{
sprintf(errBuf, "methodName %s not found", methodName);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* success */
return(Py_BuildValue("i", 0));
}
/* methods to call fortran/C methods from python */
/***********************************************************************
*
* madIsrim runs Shunrong's underlying Fortran method isrim
*
* arguments:
* 1. Python double representing UT1
* 2. Python int representing kinst
* 3. Python double represent solar local time
* 4. Python double representing GDALT in km
* 5. Python double representing GDLAT in degrees
* 6. Python double representing F107 24 hours before UT1 in W/m2/Hz
* 7. Python double representing AP3 3 hours before
* 8. Python int representing IPAR flag
* 1, FOR NE
* 2, FOR TE
* 3, FOR TI
* 4, FOR PARALLEL ION DRIFT, + UPWARD
* 5, for HMAX and NMAX
* 9. Python double 1D numpy vector size (4,) representing the output - will be modified
*
* returns:
* Python int 0
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_madIsrim(PyObject * self, PyObject * args)
{
/* local variables */
char errBuf[1000] = "";
double * outValues;
int iyr=0, imd=0, ihm=0, ics=0;
int month=0, day=0, i=0;
double doy = 0.0; /* day of year (1-366) */
int ipar2;
/* outputs of Shunrong's model */
double output[4];
/* input variables */
double ut1;
int kinst;
double slt;
double gdalt;
double gdlat;
double f107;
double ap3;
int ipar;
PyArrayObject *outArr;
/* init */
for (i=0; i<4; i++)
output[i] = 0.0; /* also used as an input by Shunrong's code */
if (!PyArg_ParseTuple(args, "didddddiO!",
&ut1, &kinst, &slt, &gdalt,
&gdlat, &f107, &ap3, &ipar,
&PyArray_Type, &outArr))
return NULL;
/* check that one-D array */
if (outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of out array must be 1, not %i",
outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
if (ipar == 5)
{
/* set to get hmax and nmax */
output[3] = 1.0;
ipar2 = 1; /* ne required to be set */
}
else
ipar2 = ipar;
outValues = (double *)outArr->data;
/* get doy */
dinvmadptr(ut1, &iyr, &imd, &ihm, &ics);
month = imd/100;
day = imd - (100*month);
doy = (double)madGetDayno(iyr, month, day);
ISRIM_F77(&kinst, &doy, &slt, &gdalt, &gdlat, &f107, &ap3, &ipar2, output);
if (ipar == 1)
{
/* handle Ne case */
if ((output[3] < 0.0) || (output[0] <= 0.0))
{
outValues[0] = NAN;
outValues[1] = NAN;
}
else
{
outValues[0] = output[0];
outValues[1] = log10(output[0]);
}
}
else if (ipar == 5)
{
/* handle Hmax, Nmax case */
if (output[1] < 1.0)
{
outValues[0] = NAN;
outValues[1] = NAN;
}
else
{
outValues[0] = output[1];
outValues[1] = output[2];
}
}
else /* all others just return one value */
{
if (output[3] < 0.0)
outValues[0] = NAN;
else
outValues[0] = output[0];
}
/* success */
return(Py_BuildValue("i", 0));
}
/***********************************************************************
*
* madRunIri runs underlying IRI methods via c method run_iri_3
*
* arguments:
* 1-6. Python ints representing year, month, day, hour, min, sec
* 7. Python double representing GDLAT in degrees
* 8. Python double representing GLON in degrees
* 9. Python double representing GDALT in km
* 10. Numpy int32 array of length 13 with 13 previous values of ap3 (last is present value)
* 11. Python double representing F107 in W/m2/Hz
* 12. Python double 1D numpy vector representing the 11 outputs:
* NE_IRI Electron density in #/meter3. units: m-3
* NEL_IRI Log of electron density in #/meter3. units:log10(m-3)
* TN_IRI IRI Neutral temperature
* TI_IRI IRI Ion temperature
* TE_IRI IRI Electron temperature
* PO+_IRI IRI Composition - [O+]/Ne
* PNO+_IRI IRI Composition - [NO+]/Ne
* PO2+_IRI IRI Composition - [O2+]/Ne
* PHE+_IRI IRI Composition - [HE+]/Ne
* PH+_IRI IRI Composition - [H+]/Ne
* PN+_IRI IRI Composition - [N+]/Ne
*
* returns:
* Python int 0
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_madRunIri(PyObject * self, PyObject * args)
{
/* local variables */
char errBuf[1000] = "";
int * ap3Values;
double * outValues;
/* input variables */
int year, month, day, hour, minute, second;
double gdlat, glon, gdalt;
double f107;
PyArrayObject *ap3Arr, *outArr;
if (!PyArg_ParseTuple(args, "iiiiiidddO!dO!",
&year, &month, &day, &hour, &minute, &second,
&gdlat, &glon, &gdalt,
&PyArray_Type, &ap3Arr,
&f107,
&PyArray_Type, &outArr))
return NULL;
/* check that one-D array */
if (ap3Arr->nd != 1 || outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of ap3 or out array must be 1, not %i and %i",
ap3Arr->nd, outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
ap3Values = (int *)ap3Arr->data;
outValues = (double *)outArr->data;
run_iri_3(year, month, day, hour, minute, second,
gdlat, glon, gdalt,
ap3Values, f107, outValues);
/* success */
return(Py_BuildValue("i", 0));
}
/***********************************************************************
*
* madRunMsis runs underlying MSIS methods via fortran method gtd7d
*
* arguments:
* 1-6. Python ints representing year, month, day, hour, min, sec
* 7. Python double representing GDLAT in degrees
* 8. Python double representing GLON in degrees
* 9. Python double representing GDALT in km
* 10. Numpy double array of length 7 with AP data from present and up to 57 hours prior to
* present. (see gtd7d for algorithm)
* 11. Python double representing FBAR in 1.0e-22 W/m2/Hz
* 12. Python double representing F107 24 hours previous in 1.0e-22 W/m2/Hz
* 13. Python double 1D numpy vector representing the 26 outputs:
* TNM - Model neutral atmosphere temperature
* TINFM - Model Exospheric temperature
* MOL - Log10 (nutrl atm mass density in kg/m3)
* NTOTL - Log10 (nutrl atm number density in m-3)
* NN2L - Nutrl atm compositn-log10([N2] in m-3)
* NO2L - Nutrl atm compositn-log10([O2] in m-3)
* NOL - Nutrl atm composition-log10([O] in m-3)
* NARL - Nutrl atm compositn-log10([AR] in m-3)
* NHEL - Nutrl atm compositn-log10([HE] in M-3)
* NHL - Nutrl atm composition-log10([H] in m-3)
* NN4SL - Nutrl atm compstn-log10([N(4S)] in m-3)
* NPRESL - Neutral atmospher log10(pressure in Pa)
* PSH - Pressure scale height (m)
* DTNM - Error in Model neutral atmosphere temperature
* DTINFM - Error in Model Exospheric temperature
* DMOL - Error in Log10 (nutrl atm mass density in kg/m3)
* DNTOTL - Error in Log10 (nutrl atm number density in m-3)
* DNN2L - Error in Nutrl atm compositn-log10([N2] in m-3)
* DNO2L - Error in Nutrl atm compositn-log10([O2] in m-3)
* DNOL - Error in Nutrl atm composition-log10([O] in m-3)
* DNARL - Error in Nutrl atm compositn-log10([AR] in m-3)
* DNHEL - Error in Nutrl atm compositn-log10([HE] in M-3)
* DNHL - Error in Nutrl atm composition-log10([H] in m-3)
* DNN4SL - Error in Nutrl atm compstn-log10([N(4S)] in m-3)
* DNPRESL - Error in Neutral atmospher log10(pressure in Pa)
* DPSH - Error in Pressure scale height (m)
*
* returns:
* Python int 0
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_madRunMsis(PyObject * self, PyObject * args)
{
/* constants */
const double k = 1.3807e-23;
const double g = 9.8;
const double errPercent = 0.2; /* assumed error percentage */
/* local variables */
char errBuf[1000] = "";
double * apValues;
double * outValues;
double ave_mass = 0.0;
int i = 0;
/* declaration of data to be passed into gtd7d */
int iyd = 0; /* date in form YYDDD */
double sec = 0.0; /* seconds since beginning of UT day */
double slt = 0.0; /* slt=sec/3600+glon/15 */
int mass = 48; /* tells gtd7d to calculate all */
/* declaration of data returned from gtd7d */
double d[9];
double t[2];
/* input variables */
int year, month, day, hour, minute, second;
double gdlat, glon, gdalt;
double f107, fbar;
PyArrayObject *apArr, *outArr;
/* init arrays to keep purify happy */
for (i=0; i<9; i++)
d[i] = 0.0;
for (i=0; i<2; i++)
t[i] = 0.0;
if (!PyArg_ParseTuple(args, "iiiiiidddO!ddO!",
&year, &month, &day, &hour, &minute, &second,
&gdlat, &glon, &gdalt,
&PyArray_Type, &apArr,
&fbar, &f107,
&PyArray_Type, &outArr))
return NULL;
/* check that one-D array */
if (apArr->nd != 1 || outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of apArr or out array must be 1, not %i and %i",
apArr->nd, outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* force glon to between -180 and 180 */
while (glon < -180.0) glon += 360.0;
while (glon > +180.0) glon -= 360.0;
/* weird time input calculations */
iyd = (year%100)*1000; /* set up year part of YYDDD */
iyd += madGetDayno(year, month, day); /* set up day part of YYDDD */
sec = hour*3600.0;
sec += minute*60.0;
slt = sec/3600.0 + glon/15.0;
apValues = (double *)apArr->data;
outValues = (double *)outArr->data;
/* call main MSIS method */
GTD7D_F77(&iyd,
&sec,
&gdalt,
&gdlat,
&glon,
&slt,
&fbar,
&f107,
apValues,
&mass,
d,
t);
/* get average mass in kg */
if (d[0]+d[1]+d[2]+d[3]+d[4]+d[6]+d[7] > 0.0 && d[5] > 0.0)
{
ave_mass = d[5]/(d[0]+d[1]+d[2]+d[3]+d[4]+d[6]+d[7]);
/* convert from grams to kg */
ave_mass = ave_mass*1.e-3;
/* get pressure */
outValues[11]=(d[0]+d[1]+d[2]+d[3]+d[4]+d[6]+d[7])*1.e6*k*t[1];
/* get scale height */
outValues[12]=k*t[1]/(ave_mass*g);
}
else
{
outValues[11]=NAN;
outValues[12]=NAN;
}
/* copy results */
outValues[0] = t[1]; /* tnm */
outValues[1] = t[0]; /* tinfm */
if (d[5] > 0.0)
outValues[2] = log10(d[5]*1.e3); /* mol */
else
outValues[2] = NAN;
if (d[0]+d[1]+d[2]+d[3]+d[4]+d[6]+d[7] > 0.0)
outValues[3] = log10(1.e6*(d[0]+d[1]+d[2]+d[3]+d[4]+d[6]+d[7])); /* neutral density */
else
outValues[3] = NAN;
if (d[2] > 0.0)
outValues[4] = log10(1.e6*d[2]); /* N2 density */
else
outValues[4] = NAN;
if (d[3] > 0.0)
outValues[5] = log10(1.e6*d[3]); /* O2 density */
else
outValues[5] = NAN;
if (d[1] > 0.0)
outValues[6] = log10(1.e6*d[1]); /* O density */
else
outValues[6] = NAN;
if (d[4] > 0.0)
outValues[7] = log10(1.e6*d[4]); /* Ar density */
else
outValues[7] = NAN;
if (d[0] > 0.0)
outValues[8] = log10(1.e6*d[0]); /* He density */
else
outValues[8] = NAN;
if (d[6] > 0.0)
outValues[9] = log10(1.e6*d[6]); /* H density */
else
outValues[9] = NAN;
if (d[7] > 0.0)
outValues[10] = log10(1.e6*d[7]); /* N - stable density */
else
outValues[10] = NAN;
/* simple error calculation */
for (i=0; i<13; i++)
{
if (isnan(outValues[i]))
outValues[i+13] = NAN;
else if (i==0 || i==1 || i==12)
outValues[i+13] = outValues[i]*errPercent; /* linear parm */
else
outValues[i+13] = log10( pow(10.0,outValues[i]) * errPercent ); /* log parm */
}
/* success */
return(Py_BuildValue("i", 0));
}
/***********************************************************************
*
* madRunTsygan runs underlying Tsyganenko method getTsyganenkoField
*
* arguments:
* 1. Python double representing mid_time in seconds since 1950-01-01
* 2. Python double representing GDLAT in degrees
* 3. Python double representing GLON in degrees
* 4. Python double representing GDALT in km
* 5. Numpy double array of length 24 representing last 24 hours of SWSPD in m/s (last is present value)
* 6. Python double representing imf in y gsm direction in nTesla at mid_time
* 7. Numpy double array of length 24 representing last 24 hours of z gsm values in nT (last is present value)
* 8. Numpy double array of length 24 representing last 24 hours of SWDEN in m^-3 (last is present value)
* 9. Python double representing dst in nTesla at mid_time
* 10. Python double 1D numpy vector representing the 4 outputs:
* "TSYG_EQ_XGSM","TSYG_EQ_YGSM","TSYG_EQ_XGSE","TSYG_EQ_YGSE"
*
* returns:
* Python int 0
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_madRunTsygan(PyObject * self, PyObject * args)
{
/* local variables */
char errBuf[1000] = "";
double * swspdValues;
double * zgsmValues;
double * swdenValues;
double * outValues;
/* input variables */
double mid_time;
double gdlat, glon, gdalt;
double dst, ygsm;
PyArrayObject *swspdArr, *zgsmArr;
PyArrayObject *swdenArr, *outArr;
if (!PyArg_ParseTuple(args, "ddddO!dO!O!dO!",
&mid_time,
&gdlat, &glon, &gdalt,
&PyArray_Type, &swspdArr,
&ygsm,
&PyArray_Type, &zgsmArr,
&PyArray_Type, &swdenArr,
&dst,
&PyArray_Type, &outArr))
return NULL;
/* check that one-D arrays */
if (swspdArr->nd != 1 || zgsmArr->nd != 1 || swdenArr->nd != 1 || outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of all arrays must be 1, not %i, %i %i, and %i",
swspdArr->nd, zgsmArr->nd, swdenArr->nd, outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* force glon to between -180 and 180 */
while (glon < -180.0) glon += 360.0;
while (glon > +180.0) glon -= 360.0;
swspdValues = (double *)swspdArr->data;
zgsmValues = (double *)zgsmArr->data;
swdenValues = (double *)swdenArr->data;
outValues = (double *)outArr->data;
getTsyganenkoField(mid_time,
gdlat, glon, gdalt,
swspdValues, ygsm,
zgsmValues, swdenValues,
dst,
outValues,
outValues + 1,
outValues + 2,
outValues + 3);
/* success */
return(Py_BuildValue("i", 0));
}
/***********************************************************************
*
* traceTsyganenkoField returns the endpoint of a magnetic line using Tsyganenko model
*
* Traces to conjugate point, intersection with a given altitude in the
* northern or southern hemisphere, to the apex, or to GSM XY plane
* Uses TsyganenkoF fields.
* Input arguments are either GSM or Geodetic, depending on inputType argument.
* Output arguments are either GSM or Geodetic, depending on outputType
* argument.
*
* arguments:
* Python int - year
* Python int - month
* Python int - day
* Python int - hour
* Python int - min
* Python int - sec
* Python int - inputType (0 for geodetic, 1 for GSM)
* Python int - outputType (0 for geodetic, 1 for GSM)
* (the following parameter depend on inputType)
* Python double in1 - geodetic altitude or ZGSM of starting point
* Python double in2 - geodetic latitude or XGSM of starting point
* Python double in3 - longitude or YGSM of starting point
* Numpy double array of swspd_array - array of last 24 hourly solar wind speeds in m/s
* Python double imf_ygsm_now - imf in y gsm direction in nTesla now
* Numpy double array of imf_zgsm_array - array of last 24 hourly imf in z gsm direction in nTesla measured now
* Numpy double array of swden_array - array of last 24 hourly solar wind density in m^-3
* Python double dst - dst in nTesla
* Python int qualifier - 0 for conjugate, 1 for north_alt, 2 for south_alt, 3 for apex, 4 for GSM XY plane
* Python double stopAlt - altitude to stop trace at, if qualifier is north_alt or south_alt.
* If other qualifier, this parameter is ignored
* Python double 1D numpy vector representing the 3 outputs, whose meaning depends on outputType
* end_1 - geodetic altitude or ZGSM of starting point
* end_2 - geodetic latitude or XGSM of starting point
* end_3 - longitude or YGSM of starting point
*
* If trace fails, all three return values will be set to nan
*
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_traceTsyganenkoField(PyObject * self, PyObject * args)
{
/* python inputs */
int year = 0;
int month = 0;
int day = 0;
int hour = 0;
int min = 0;
int sec = 0;
int inputType = 0;
int outputType = 0;
double in1 = 0.0;
double in2 = 0.0;
double in3 = 0.0;
PyArrayObject *swspdArr, *zgsmArr;
double imf_ygsm_now;
PyArrayObject *swdenArr, *outArr;
double dst;
int qualifier = 0;
double stopAlt = 0.0;
/* local variables */
char errBuf[1000] = "";
int dayOfYear = 0;
int imod = 0; /* flag to set conversion direction */
int result = 0, i=0;
double * swspdValues;
double * zgsmValues;
double * swdenValues;
double * outValues;
double mid_time;
double out1, out2, out3;
/* inputs required for the latest Tsygenenko code - hard coded solar wind speeds */
double vx = 0.0;
double vy = 29.8;
double vz = 0.0;
/* temporary variable for holding position during conversion */
double dim1 = 0.0;
double dim2 = 0.0;
double dim3 = 0.0;
// get arguments
if (!PyArg_ParseTuple(args, "iiiiiiiidddO!dO!O!didO!",
&year, &month, &day, &hour, &min, &sec,
&inputType, &outputType,
&in1, &in2, &in3,
&PyArray_Type, &swspdArr,
&imf_ygsm_now,
&PyArray_Type, &zgsmArr,
&PyArray_Type, &swdenArr,
&dst,
&qualifier, &stopAlt,
&PyArray_Type, &outArr))
{
return NULL;
}
/* check that one-D arrays */
if (swspdArr->nd != 1 || zgsmArr->nd != 1 || swdenArr->nd != 1 || outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of all arrays must be 1, not %i, %i %i, and %i",
swspdArr->nd, zgsmArr->nd, swdenArr->nd, outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* force glon to between -180 and 180 if geodetic */
if (inputType == 0)
{
while (in3 < -180.0) in3 += 360.0;
while (in3 > +180.0) in3 -= 360.0;
}
swspdValues = (double *)swspdArr->data;
zgsmValues = (double *)zgsmArr->data;
swdenValues = (double *)swdenArr->data;
outValues = (double *)outArr->data;
/* convert time to seconds since 1/1/1950 */
mid_time = getKey(year, month, day, hour, min, sec);
/* get day of year as required by Tsyganenko */
dayOfYear = madGetDayno(year, month, day);
/* set Tsyganenko common block data to the right time - not thread safe!!!! */
TS_RECALC_F77(&year, &dayOfYear, &hour, &min, &sec, &vx, &vy, &vz);
/* convert in* to geodetic if needed */
if (inputType == 1)
{
imod = -1; /* converst gsm to cartesian */
GEOGSM_F77(&dim1, &dim2, &dim3, &in2, &in3, &in1, &imod);
/* now dim1 XGEO, dim2 YGEO, dim3 ZGEO */
/* convert cartesian to spherical */
TS_SPHCAR_F77(&in1, &in2, &in3, &dim1, &dim2, &dim3, &imod);
/* now in1 R, in2 theta, in3 phi */
/* convert spherical to geocentric */
dim1 = in1*6371.2;
dim2 = 90.0 - (in2/0.01745329);
dim3 = in3/0.01745329;
imod = 2; /* convert geocentric coordinates to geodetic */
in3 = dim3;
CONVRT_F77(&imod, &in2, &in1, &dim2, &dim1);
/* now in1 is gdalt, in2 is gdlat, in3 is glon */
}
result = traceTsyganenkoField(mid_time,
in2,
in3,
in1,
swspdValues,
imf_ygsm_now,
zgsmValues,
swdenValues,
dst,
qualifier,
stopAlt,
&out2,
&out3,
&out1);
/* out1 gdalt or ZGSM, out2 gdlat or XGSM, out3 glon or YGSM */
/* convert to GSM if result okay and outputType == 1 and qualifier != 4 (since then already GSM output) */
if (outputType == 1 && result == 0 && qualifier != 4)
{
imod = 1;
/* get geocentric coordinates from geodetic */
CONVRT_F77(&imod, &out2, &out1, &dim2, &dim1);
dim3 = out3;
/* dim2 now gclat, dim1 = rkm, dim3 = glon */
/* convert from geocentric to spherical */
out2 = (90.0 - dim2)*0.01745329;
out3 = dim3*0.01745329;
out1 = dim1/6371.2;
/* now out1 = R, out2 = theta, out3 = phi */
/* get cartesian coordinates from spherical */
TS_SPHCAR_F77(&out1, &out2, &out3, &dim1, &dim2, &dim3, &imod);
/* now dim1-XGEO, dim2=YGEO, dim3=ZGEO */
/* convert cartesian into solar magnetospheric ones */
GEOGSM_F77(&dim1, &dim2, &dim3, &out2, &out3, &out1, &imod);
/* now out1 = ZGSM, out2=XGSM, out3 = YGSM */
}
/* convert to Geodetic if result okay and outputType == 0 and qualifier == 4 (since then GSM output) */
if (outputType == 0 && result == 0 && qualifier == 4)
{
imod = -1; /* converst gsm to cartesian */
GEOGSM_F77(&dim1, &dim2, &dim3, &out2, &out3, &out1, &imod);
/* now dim1 XGEO, dim2 YGEO, dim3 ZGEO */
/* convert cartesian to spherical */
TS_SPHCAR_F77(&out1, &out2, &out3, &dim1, &dim2, &dim3, &imod);
/* now out1 R, out2 theta, out3 phi */
/* convert spherical to geocentric */
dim1 = out1*6371.2;
dim2 = 90.0 - (out2/0.01745329);
dim3 = out3/0.01745329;
imod = 2; /* convert geocentric coordinates to geodetic */
out3 = dim3;
CONVRT_F77(&imod, &out2, &out1, &dim2, &dim1);
/* now out1 is gdalt, out2 is gdlat, out3 is glon */
}
/* force out3 to between -180 and 180 if geodetic */
if (outputType == 0)
{
while (out3 < -180.0) out3 += 360.0;
while (out3 > +180.0) out3 -= 360.0;
}
if (result != 0)
for (i=0; i<3; i++)
outValues[i] = NAN;
else
{
outValues[0] = out1;
outValues[1] = out2;
outValues[2] = out3;
}
/* success */
return(Py_BuildValue("i", 0));
}
/***********************************************************************
*
* traceIGRFField returns the endpoint of a magnetic line using IGRF model via lintra
*
* Traces to conjugate point, intersection with a given altitude in the
* northern or southern hemisphere, to the apex, or to GSM XY plane
* Uses IGRF fields.
* Input arguments are either GSM or Geodetic, depending on inputType argument.
* Output arguments are either GSM or Geodetic, depending on outputType
* argument.
*
* arguments:
* Python int - year
* Python int - inputType (0 for geodetic, 1 for GSM)
* Python int - outputType (0 for geodetic, 1 for GSM)
* (the following parameter depend on inputType)
* Python double in1 - geodetic altitude or ZGSM of starting point
* Python double in2 - geodetic latitude or XGSM of starting point
* Python double in3 - longitude or YGSM of starting point
* Python int qualifier - 0 for conjugate, 1 for north_alt, 2 for south_alt, 3 for apex
* 4 for GSM XY plane - but not possible for lintra, so raises an error
* Python double stopAlt - altitude to stop trace at, if qualifier is north_alt or south_alt.
* If other qualifier, this parameter is ignored
* Python double 1D numpy vector representing the 3 outputs, whose meaning depends on outputType
* end_1 - geodetic altitude or ZGSM of starting point
* end_2 - geodetic latitude or XGSM of starting point
* end_3 - longitude or YGSM of starting point
*
* If trace fails, all three return values will be set to nan
*
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_traceIGRFField(PyObject * self, PyObject * args)
{
/* python inputs */
int year = 0;
int inputType = 0;
int outputType = 0;
double in1 = 0.0;
double in2 = 0.0;
double in3 = 0.0;
int qualifier = 0;
double stopAlt = 0.0;
PyArrayObject *outArr;
/* local variables */
char errBuf[1000] = "";
int result = 0;
double * outValues;
double fyear;
/* temporary variable for holding position during conversion */
double dim1 = 0.0;
double dim2 = 0.0;
double dim3 = 0.0;
/* geocentric data */
double gclat = 0.0;
double rkm = 0.0;
double plat = 0.0; /* end point geocentric latitude (deg) */
double plon = 0.0; /* end point geocentric longitude (deg) */
double prkm = 0.0; /* end point radial distance (km) */
double arc = 0.0; /* arc length of field line traced (km) */
double arad = 0.0; /* apex radius of field line (earth radii) */
double alat = 0.0; /* apex latitude of field line (deg) */
double alon = 0.0; /* apex longitude of field line (deg) */
double end_gdlat;
double end_glon;
double end_gdalt;
/* conversion direction */
int imod = 0;
int ier = 0; /* error indicator */
int npr = 0; /* used by lintra */
int init = 2; /* used by lintra */
int istop = 0; /* sets direction of lintra trace */
// get arguments
if (!PyArg_ParseTuple(args, "iiidddidO!",
&year,
&inputType, &outputType,
&in1, &in2, &in3,
&qualifier, &stopAlt,
&PyArray_Type, &outArr))
{
return NULL;
}
if (qualifier == 4)
{
sprintf(errBuf, "Qualifier for IGRF model cannot be 4");
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* check that one-D arrays */
if (outArr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of output array must be 1, not %i",
outArr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* force glon/in3 to between -180 and 180 if geodetic */
if (inputType == 0)
{
while (in3 < -180.0) in3 += 360.0;
while (in3 > +180.0) in3 -= 360.0;
}
outValues = (double *)outArr->data;
fyear = (double)year;
/* convert in* to geodetic if needed */
if (inputType == 1)
{
imod = -1; /* converst gsm to cartesian */
GEOGSM_F77(&dim1, &dim2, &dim3, &in2, &in3, &in1, &imod);
/* now dim1 XGEO, dim2 YGEO, dim3 ZGEO */
/* convert cartesian to spherical */
TS_SPHCAR_F77(&in1, &in2, &in3, &dim1, &dim2, &dim3, &imod);
/* now in1 R, in2 theta, in3 phi */
/* convert spherical to geocentric */
dim1 = in1*6371.2;
dim2 = 90.0 - (in2/0.01745329);
dim3 = in3/0.01745329;
imod = 2; /* convert geocentric coordinates to geodetic */
in3 = dim3;
CONVRT_F77(&imod, &in2, &in1, &dim2, &dim1);
/* now in1 is gdalt, in2 is gdlat, in3 is glon */
}
imod = 1; /* convert geodetic to geocentric */
CONVRT_F77(&imod, &in2, &in1, &gclat, &rkm);
switch (qualifier)
{
case 0: /* conjugate - try the opposite hemisphere first, but its possible
we'll need to try this hemisphere */
istop = 1;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&in1,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert back to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
/* check for success - more than 50 km traced, and nearly the same altitude returned */
if (arc > 50.0)
{
if (fabs(end_gdalt - in1) < 50.0)
{
end_glon = plon;
break;
}
}
/* opposite hemisphere failed - try the same hemisphere */
istop = -1;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&in1,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert back to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
/* check for success - more than 50 km traced, and nearly the same altitude returned */
if (arc > 50.0)
{
if (fabs(end_gdalt - in1) < 50.0)
{
end_glon = plon;
break;
}
}
/* both hemispheres failed */
end_gdalt = NAN;
end_gdlat = NAN;
end_glon = NAN;
break;
case 1: /* north_alt - randomly try opposite hemisphere first (50/50 guess)*/
istop = 1;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&stopAlt,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert back to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
/* check for success - more than 50 km traced, and north of starting point */
if (arc > 50.0)
{
if (end_gdlat - in2 > 0.1)
{
end_glon = plon;
break;
}
/* or we may have started below the intercept */
if (in1 < stopAlt && end_gdlat - in2 < -0.1)
{
end_glon = plon;
break;
}
}
/* opposite hemisphere failed - try the same hemisphere */
istop = -1;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&stopAlt,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert back to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
/* check for success */
if (ier == 0)
{
end_glon = plon;
break;
}
/* both hemispheres failed */
end_gdalt = NAN;
end_gdlat = NAN;
end_glon = NAN;
break;
case 2: /* south_alt - randomly try opposite hemisphere first (50/50 guess)*/
istop = 1;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&stopAlt,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert back to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
/* check for success - more than 50 km traced, and south of starting point */
if (arc > 50.0)
{
if (end_gdlat - in2 < -0.1)
{
end_glon = plon;
break;
}
/* or we may have started below the intercept */
if (in1 < stopAlt && end_gdlat - in2 > 0.1)
{
end_glon = plon;
break;
}
}
/* opposite hemisphere failed - try the same hemisphere */
istop = -1;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&stopAlt,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert back to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
/* check for success - more than 50 km traced, and south of starting point */
if (ier == 0)
{
end_glon = plon;
break;
}
/* both hemispheres failed */
end_gdalt = NAN;
end_gdlat = NAN;
end_glon = NAN;
break;
case 3: /* apex - set istop = 0 */
istop = 0;
LINTRA_F77(&fyear,
&gclat,
&in3,
&rkm,
&in1,
&in1,
&plat,
&plon,
&prkm,
&arc,
&arad,
&alat,
&alon,
&istop,
&npr,
&init,
&ier);
/* convert apex to geodetic */
imod = 2;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &plat, &prkm);
end_glon = plon;
break;
default:
sprintf(errBuf, "qualifier must be 0-3, not %i",
qualifier);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
/* end_gdalt gdalt or ZGSM, end_gdlat gdlat or XGSM, end_glon glon or YGSM */
/* convert to GSM if result okay and outputType == 1 and qualifier != 4 (since then already GSM output) */
if (outputType == 1 && result == 0 && qualifier != 4)
{
imod = 1;
/* get geocentric coordinates from geodetic */
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &dim2, &dim1);
dim3 = end_glon;
/* dim2 now gclat, dim1 = rkm, dim3 = glon */
/* convert from geocentric to spherical */
end_gdlat = (90.0 - dim2)*0.01745329;
end_glon = dim3*0.01745329;
end_gdalt = dim1/6371.2;
/* now end_gdalt = R, end_gdlat = theta, end_glon = phi */
/* get cartesian coordinates from spherical */
TS_SPHCAR_F77(&end_gdalt, &end_gdlat, &end_glon, &dim1, &dim2, &dim3, &imod);
/* now dim1-XGEO, dim2=YGEO, dim3=ZGEO */
/* convert cartesian into solar magnetospheric ones */
GEOGSM_F77(&dim1, &dim2, &dim3, &end_gdlat, &end_glon, &end_gdalt, &imod);
/* now end_gdalt = ZGSM, end_gdlat=XGSM, end_glon = YGSM */
}
/* convert to Geodetic if result okay and outputType == 0 and qualifier == 4 (since then GSM output) */
if (outputType == 0 && result == 0 && qualifier == 4)
{
imod = -1; /* converst gsm to cartesian */
GEOGSM_F77(&dim1, &dim2, &dim3, &end_gdlat, &end_glon, &end_gdalt, &imod);
/* now dim1 XGEO, dim2 YGEO, dim3 ZGEO */
/* convert cartesian to spherical */
TS_SPHCAR_F77(&end_gdalt, &end_gdlat, &end_glon, &dim1, &dim2, &dim3, &imod);
/* now end_gdalt R, end_gdlat theta, end_glon phi */
/* convert spherical to geocentric */
dim1 = end_gdalt*6371.2;
dim2 = 90.0 - (end_gdlat/0.01745329);
dim3 = end_glon/0.01745329;
imod = 2; /* convert geocentric coordinates to geodetic */
end_glon = dim3;
CONVRT_F77(&imod, &end_gdlat, &end_gdalt, &dim2, &dim1);
/* now end_gdalt is gdalt, end_gdlat is gdlat, end_glon is glon */
}
outValues[0] = end_gdalt;
outValues[1] = end_gdlat;
outValues[2] = end_glon;
/* success */
return(Py_BuildValue("i", 0));
}
/***********************************************************************
*
* faradayRotation calculates total phase shift due to a single pass of an
* electromagnetic wave from a ground instrument to a point
* in space. Also calculates the total tec along that pass,
* plus the total tec is the path were extended to 22000 km.
*
* Uses quasi-longitudinal approximation, so not appropriate if perp
* to magnetic field.
* Uses coord and IRI from geolib
*
* Inputs:
* Python ints - year month day hour min sec
* Python floats - sgdlat slon sgdalt gdlat glon gdalt freq
* Numpy int32 array of length 13 with 13 previous values of ap3 (last is present value)
* Python double representing F107 in W/m2/Hz
*
* where sgdlat slon sgdalt are geodetic station location,
* gdlat glon gdalt are geodetic point locations, and freq is
* wave frequency in Hz
*
* Returns: a tuple with three items:
* 1. one way faraday rotation in radians, NAN if error
* 2. total tec from station to point in electrons/m^2
* 3. total tec from station to 22000 km along line through
* point in electrons/m^2
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_faradayRotation(PyObject * self, PyObject * args)
{
char errBuf[1000] = "";
int year = 0;
int month = 0;
int day = 0;
int hour = 0;
int min = 0;
int sec = 0;
double sgdlat = 0.0;
double slon = 0.0;
double sgdalt = 0.0;
double gdlat = 0.0;
double glon = 0.0;
double gdalt = 0.0;
double freq = 0.0;
double tec = 0.0;
double total_tec = 0.0;
double result = 0;
int * ap3Values;
double f107;
PyArrayObject *ap3Arr;
PyObject * pyRotation;
PyObject * pyTec;
PyObject * pyTecTotal;
// get arguments
if (!PyArg_ParseTuple(args, "iiiiiidddddddO!d",
&year, &month, &day,
&hour, &min, &sec,
&sgdlat, &slon, &sgdalt,
&gdlat, &glon, &gdalt,
&freq,
&PyArray_Type, &ap3Arr,
&f107))
{
return NULL;
}
/* check that one-D array */
if (ap3Arr->nd != 1)
{
sprintf(errBuf, "Number of dimensions of ap3 array must be 1, not %i\\n",
ap3Arr->nd);
PyErr_SetString(PyExc_TypeError, errBuf);
return NULL;
}
ap3Values = (int *)ap3Arr->data;
result = faraday_rotation_3(year, month, day, hour, min, sec,
sgdlat, slon, sgdalt,
gdlat, glon, gdalt, ap3Values, f107, freq,
&tec, &total_tec);
// create min longitude
pyRotation = PyFloat_FromDouble(result);
pyTec = PyFloat_FromDouble(tec);
pyTecTotal = PyFloat_FromDouble(total_tec);
//return all info;
return Py_BuildValue("(OOO)", pyRotation, pyTec, pyTecTotal);
}
/***********************************************************************
*
* radarToGeodetic returns gdalt,gdlat,glon given radar location and az,el,range
*
*
* arguments:
* 1. radar geodetic latitude
* 2. radar longitude
* 3. radar geodetic altitude
* 4. azimuth
* 5. elevation
* 6. range in km
*
* returns:
* 1. Python tuple with (gdlat,glon,gdalt) of point
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_radarToGeodetic(PyObject * self, PyObject * args)
{
PyObject *retObj;
/* input arguments */
double slatgd = 0.0;
double slon = 0.0;
double saltgd = 0.0;
double az = 0.0;
double el = 0.0;
double range = 0.0;
/* return values */
double gdalt = 0.0;
double gdlat = 0.0;
double glon = 0.0;
/* arguments needed by point */
double sr = 0.0;
double slatgc = 0.0;
int direction = 1;
double pr = 0.0;
double gclat = 0.0;
// get input arguments
if (!PyArg_ParseTuple(args, "dddddd",
&slatgd,
&slon,
&saltgd,
&az,
&el,
&range))
{
return NULL;
}
if (range < 0.0)
{
PyErr_SetString(PyExc_TypeError, "Range must be > 0 for radarToGeodetic");
return NULL;
}
/* force slon to between -180 and 180 */
while (slon < -180.0) slon += 360.0;
while (slon > +180.0) slon -= 360.0;
/* get geocentric coordinates needed by point */
CONVRT_F77(&direction,&slatgd,&saltgd,&slatgc,&sr);
POINT_F77(&sr,&slatgc,&slon,&az,&el,&range,
&pr,&gclat,&glon);
direction = 2;
CONVRT_F77(&direction,&gdlat,&gdalt,&gclat,&pr);
/* success */
retObj = Py_BuildValue("(ddd)", gdlat,glon,gdalt);
return(retObj);
}
/***********************************************************************
*
* geodeticToRadar given gdalt,gdlat,glon and radar location, returns az,el,range
*
*
* arguments:
* 1. radar geodetic latitude
* 2. radar longitude
* 3. radar geodetic altitude
* 4. point geodetic latitude
* 5. point longitude
* 6. point geodetic altitude
*
* returns:
* 1. Python tuple with (az,el,range) of radar view of point
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_geodeticToRadar(PyObject * self, PyObject * args)
{
PyObject *retObj;
/* input arguments */
double slatgd = 0.0;
double slon = 0.0;
double saltgd = 0.0;
double gdalt = 0.0;
double gdlat = 0.0;
double glon = 0.0;
/* return values */
double az = 0.0;
double el = 0.0;
double range = 0.0;
/* arguments needed by look */
double sr = 0.0;
double slatgc = 0.0;
double pr = 0.0;
double gclat = 0.0;
int direction = 1;
// get input arguments
if (!PyArg_ParseTuple(args, "dddddd",
&slatgd,
&slon,
&saltgd,
&gdlat,
&glon,
&gdalt))
{
return NULL;
}
/* force slon to between -180 and 180 */
while (slon < -180.0) slon += 360.0;
while (slon > +180.0) slon -= 360.0;
/* force glon to between -180 and 180 */
while (glon < -180.0) glon += 360.0;
while (glon > +180.0) glon -= 360.0;
/* get geocentric coordinates needed by look */
CONVRT_F77(&direction,&slatgd,&saltgd,&slatgc,&sr);
CONVRT_F77(&direction,&gdlat,&gdalt,&gclat,&pr);
LOOK_F77(&sr,&slatgc,&slon,&pr,&gclat,&glon,&az,&el,&range);
/* force az to between -180 and 180 */
while (az < -180.0) az += 360.0;
while (az > +180.0) az -= 360.0;
/* success */
retObj = Py_BuildValue("(ddd)", az,el,range);
return(retObj);
}
/***********************************************************************
*
* getUtFromDate returns the number of seconds as a float since 1/1/1950
*
* arguments:
* Year as Python int
* Month as Python int
* Day as Python int
* Hour as Python int
* Min as Python int
* Sec as Python int
* Centisec as Python int
*
* returns:
* Python double giving the number of seconds as a float since 1/1/1950
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_getUtFromDate(PyObject * self, PyObject * args)
{
int year = 0;
int month = 0;
int day = 0;
int hour = 0;
int min = 0;
int sec = 0;
int centisec = 0;
double ut = 0.0;
PyObject *retObj;
// get date arguments
if (!PyArg_ParseTuple(args, "iiiiiii", &year, &month, &day, &hour, &min, &sec, &centisec))
{
return NULL;
}
ut = getKey(year, month, day, hour, min, sec);
/* add centiseconds */
ut += centisec/100.0;
retObj = Py_BuildValue("d", ut);
return(retObj);
}
/***********************************************************************
*
* getDateFromUt returns a date as a list of integers given ut as seconds since 1/1/1950
*
* arguments:
* Python double giving the number of seconds as a float since 1/1/1950
*
* returns: list containing
* Year as Python int
* Month as Python int
* Day as Python int
* Hour as Python int
* Min as Python int
* Sec as Python int
* Centisec as Python int
* Day of year as Python int
*
* throws:
* PyExc_TypeError if illegal argument passed
*
*/
static PyObject * _derive_getDateFromUt(PyObject * self, PyObject * args)
{
int year = 0;
int month = 0;
int day = 0;
int imd = 0;
int hour = 0;
int min = 0;
int ihm = 0;
int sec = 0;
int centisec = 0;
int ics = 0;
int dayOfYear = 0;
double ut = 0.0;
PyObject * retDateList;
PyObject * pyTempValue; /* python object to add to lists, then release */
// get ut
if (!PyArg_ParseTuple(args, "d", &ut))
{
return NULL;
}
dinvmadptr(ut, &year, &imd, &ihm, &ics);
month = imd/100;
day = imd - month*100;
hour = ihm/100;
min = ihm - hour*100;
sec = ics/100;
centisec = ics - sec*100;
dayOfYear = madGetDayno(year, month, day);
retDateList = PyList_New(0);
pyTempValue = PyInt_FromLong(year);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(month);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(day);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(hour);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(min);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(sec);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(centisec);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
pyTempValue = PyInt_FromLong(dayOfYear);
PyList_Append(retDateList, pyTempValue);
Py_DECREF(pyTempValue);
return(retDateList);
}
/********** Initialization code for module ******************************/
static PyMethodDef _deriveMethods[] =
{
{"madDispatch", _derive_madDispatch, METH_VARARGS, "Call a derivation method"},
{"madIsrim", _derive_madIsrim, METH_VARARGS, "Call Shunrongs model"},
{"madRunIri", _derive_madRunIri, METH_VARARGS, "Call IRI"},
{"madRunMsis", _derive_madRunMsis, METH_VARARGS, "Call MSIS"},
{"madRunTsygan", _derive_madRunTsygan, METH_VARARGS, "Call Tsyganenko"},
{"traceTsyganenkoField", _derive_traceTsyganenkoField, METH_VARARGS, "Trace Tsyganenko field"},
{"traceIGRFField", _derive_traceIGRFField, METH_VARARGS, "Trace IGRF field"},
{"faradayRotation", _derive_faradayRotation, METH_VARARGS, "Get Faraday Rotation"},
{"radarToGeodetic", _derive_radarToGeodetic, METH_VARARGS, "Convert radar coords to geodetic"},
{"geodeticToRadar", _derive_geodeticToRadar, METH_VARARGS, "Convert geodetic to radar coords"},
{"getUtFromDate", _derive_getUtFromDate, METH_VARARGS, "get secs from 1950 UT from date"},
{"getDateFromUt", _derive_getDateFromUt, METH_VARARGS, "get date as secs from 1950 UT"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
void init_derive(void)
{
PyImport_AddModule("_derive");
Py_InitModule("_derive", _deriveMethods);
import_array();
}
"""
###### main code starts here ##########
f = open('../../madrigal/_derive_new.c', 'w')
f.write(startTemplate)
madMethods = madrigal.derivation.MadrigalDerivedMethods
methodNames = list(madMethods.keys())
for methodName in methodNames:
if methodName == methodNames[1]: # first method is python - so skip it
f.write(codeTemplate1 % (methodName, methodName))
else:
items = madMethods[methodName]
if len(items) >= 3:
if items[2] == 'python':
continue
f.write(codeTemplate2 % (methodName, methodName))
f.write(endTemplate)
f.close()