#!$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 #include #include #include /*********************************************************************** * * 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, ¢isec)) { 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()