<|start_filename|>geoist/magmod/pymm.c<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#include "common.h" /* common definitions - to be included before Python.h */
#include <Python.h>
#include <numpy/arrayobject.h>

/* module version */
#include "version.h"

/* common python utilities */
#include "py_common.h"

/* Coordinate conversions. */
#include "pymm_cconv.h"

/* spherical harmonic model evaluation */
#include "pymm_sheval.h"

/* evaluation of the associative Legendre functions */
#include "pymm_legendre.h"

/* evaluation of the relative radius powers */
#include "pymm_relradpow.h"

/* evaluation of the series of longitude sine and cosine values */
#include "pymm_lonsincos.h"

/* final spherical-harmonic gradient evaluation */
#include "pymm_sphargrd.h"

/* final spherical-harmonic gradient potential */
#include "pymm_spharpot.h"

/* vector rotations */
#include "pymm_vrot_sph2geod.h"
#include "pymm_vrot_sph2cart.h"
#include "pymm_vrot_cart2sph.h"

/*---------------------------------------------------------------------------*/
/* module's doc string */

#define DOC_PYMM \
"This module provides bindings to the Geomagnetic Model library."

/*---------------------------------------------------------------------------*/
/*define module's methods */
static PyMethodDef pymm_methods[] =
{
    {"vrot_sph2cart", (PyCFunction)vrot_sph2cart, METH_VARARGS|METH_KEYWORDS, DOC_VROT_SPH2CART},
    {"vrot_cart2sph", (PyCFunction)vrot_cart2sph, METH_VARARGS|METH_KEYWORDS, DOC_VROT_CART2SPH},
    {"vrot_sph2geod", (PyCFunction)vrot_sph2geod, METH_VARARGS|METH_KEYWORDS, DOC_VROT_SPH2GEOD},
    {"spharpot", (PyCFunction)spharpot, METH_VARARGS|METH_KEYWORDS, DOC_SPHARPOT},
    {"sphargrd", (PyCFunction)sphargrd, METH_VARARGS|METH_KEYWORDS, DOC_SPHARGRD},
    {"lonsincos", (PyCFunction)lonsincos, METH_VARARGS|METH_KEYWORDS, DOC_LONSINCOS},
    {"relradpow", (PyCFunction)relradpow, METH_VARARGS|METH_KEYWORDS, DOC_RELRADPOW},
    {"legendre", (PyCFunction)legendre, METH_VARARGS|METH_KEYWORDS, DOC_LEGENDRE},
    {"sheval", (PyCFunction)sheval, METH_VARARGS|METH_KEYWORDS, DOC_SHEVAL},
    {"convert", (PyCFunction)convert, METH_VARARGS|METH_KEYWORDS, DOC_CONVERT},
    {NULL, NULL, 0, NULL} /* Sentinel - DO NOT REMOVE! */
} ;

/*---------------------------------------------------------------------------*/
/* module initialization  */

static PyObject* init_module(void)
{
    PyObject *module = init_python_module("_pymm", DOC_PYMM, pymm_methods);
    if (NULL == module)
        goto exit;

    PyObject *dict = PyModule_GetDict(module);
    if (NULL == dict)
        goto exit;

    /* add module specific exception */
    //set_dict_item_str(dict, "MMError", PyErr_NewException("_pymm.MMError", NULL, NULL));

    /* integer constants */
    set_dict_item_str_long(dict, "GEODETIC_ABOVE_WGS84", CT_GEODETIC_ABOVE_WGS84);
    set_dict_item_str_long(dict, "GEOCENTRIC_SPHERICAL", CT_GEOCENTRIC_SPHERICAL);
    set_dict_item_str_long(dict, "GEOCENTRIC_CARTESIAN", CT_GEOCENTRIC_CARTESIAN);
    set_dict_item_str_long(dict, "POTENTIAL", SM_POTENTIAL);
    set_dict_item_str_long(dict, "GRADIENT", SM_GRADIENT);
    set_dict_item_str_long(dict, "POTENTIAL_AND_GRADIENT", SM_POTENTIAL_AND_GRADIENT);
    set_dict_item_str_double(dict, "EARTH_RADIUS", RADIUS);

    /* module metadata */
    set_dict_item_str_str(dict, "__author__", "<NAME> (<EMAIL>)");
    set_dict_item_str_str(dict, "__copyright__", "Copyright (C) 2014 EOX IT Services GmbH");
    set_dict_item_str_str(dict, "__licence__", "EOX licence (MIT style)");
    set_dict_item_str_str(dict, "__version__", VERSION);

  exit:
    return module;
}

#if PY_MAJOR_VERSION == 2

PyMODINIT_FUNC init_pymm(void)
{
    import_array();
    init_module();
}

#else

PyObject* PyInit__pymm(void)
{
    import_array();
    return init_module();
}

#endif


<|start_filename|>geoist/magmod/include/sunpos.h<|end_filename|>
/**
 * @file sunpos.h
 * @author <NAME> <<EMAIL>>
 * @brief Sun Position
 *
 * Calculation of the Solar position.
 *
 * This code is derived from the code kindly provided by <NAME>.
 *
 * Also See:
 * <NAME> (2012), Five new algorithms for the computation of sun
 * position from 2010 to 2110, Solar Energy, vol. 86(5), 1323-1337,
 * doi:10.1016/j.solener.2012.01.024
 * Ph.Blanc (2012), The SG2 algorithm for a fast and accurate computation
 * of the position of the Sun for multi-decadal time period, Solar Energy,
 * vol. 86(5), 3072-3083, doi:10.1016/j.solener.2012.07.018
 *
 * Copyright (C) 2016 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#ifndef SUNPOS_H
#define SUNPOS_H

#include "time_conversion.h"
#include <math.h>

#ifndef SEC2DAYS
#define SEC2DAYS (1.0/86400.0)
#endif

#ifndef PI
#define PI M_PI
#endif

#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923  /* pi/2 */
#endif

#define PI2 (2*M_PI) // 2*PI
#define PIM M_PI_2   // PI/2

#ifndef FMOD_FLOOR
#define FMOD_FLOOR(a,b) ((a)-(b)*floor((a)/(b)))
#endif

/**
 * @brief Convert year, month, day and number of decimal hours to MJD2000
 *
 * Convert year, month, day and number of decimal hours to MJD2000 approximation
 * used by the original sunpos algorithm implementation.
 *
 * This approximation gives correct MJD2000 days from 2000-03-01 to 2400-02-29.

 * The deviations outside this range are caused by:
 *  1) use of truncation instead of rounding towards negative infinity (floor),
 *  2) no compensation for the extra leap year every 400years.
 */

void date_to_sunposmjd2k(double *sunpos_mjd2k, int year, int month, int day, double hours) {
    const int mt = month + 12 * (month <= 2);
    const int yt = year - (month <= 2);
    *sunpos_mjd2k = (
        (int)(365.25*(yt - 2000)) +
        (int)(30.6001*(mt + 1)) -
        (int)(0.01*yt) +
        day
    ) + (1.0/24.0)*hours - 43.0;
}

/**
 * @brief Evaluate solar position in equatorial coordinate system.
 *
 * Algorithm #5 - 1:1 copy of the original implementation
 *
 * <NAME> (2012), Five new algorithms for the computation of sun
 * position from 2010 to 2110, Solar Energy, vol. 86(5), 1323-1337,
 * doi:10.1016/j.solener.2012.01.024
 * Ph.Blanc (2012), The SG2 algorithm for a fast and accurate computation
 * of the position of the Sun for multi-decadal time period, Solar Energy,
 * vol. 86(5), 3072-3083, doi:10.1016/j.solener.2012.07.018
 *
 *  outputs:
 *      decl - declination [rad]
 *      rasc - right ascension [rad] (eastward)
 *      hang - hour angle [rad] (westward!)
 *
 *  inputs:
 *      year - CE year
 *      month - month of the year 1-12
 *      day  - day of the month 1-31
 *      hours  - UT decimal hours
 *      lon - longitude [rad] (eastward)
 *      lat - longitude [rad]
 *      rad - distance from earth centre [km] (set 0 to skip parallax correction)
 *      dtt - UTC offset to TT [seconds]
 *      pres - pressure [atm] (set 0 to skip refraction correction)
 *      temp - temperature [dgC]
 */

void sunpos5original(
    double *decl, double *rasc, double *hang, double *azimuth, double *zenith,
    int year, int month, int day, double hours, double lat, double lon,
    double rad, double dtt, double press, double temp
) {
    const int mt = month + 12 * (month <= 2);
    const int yt = year - (month <= 2);
    const double t = (
        (int)(365.25*(yt - 2000)) +
        (int)(30.6001*(mt + 1)) -
        (int)(0.01*yt) +
        day
    ) + (1.0/24.0)*hours - 21958.0;
    const double te = t + 1.1574e-5*dtt;
    const double wte = 0.0172019715*te;

    const double s1 = sin(wte);
    const double c1 = cos(wte);
    const double s2 = 2.0*s1*c1;
    const double c2 = (c1+s1)*(c1-s1);
    const double s3 = s2*c1 + c2*s1;
    const double c3 = c2*c1 - s2*s1;

    const double l = (
        1.7527901 + 1.7202792159e-2*te + 3.33024e-2*s1 -
        2.0582e-3*c1 + 3.512e-4*s2 - 4.07e-5*c2 + 5.2e-6*s3 -
        9e-7*c3 -8.23e-5*s1*sin(2.92e-5*te) + 1.27e-5*sin(1.49e-3*te - 2.337) +
        1.21e-5*sin(4.31e-3*te + 3.065) + 2.33e-5*sin(1.076e-2*te - 1.533) +
        3.49e-5*sin(1.575e-2*te - 2.358) + 2.67e-5*sin(2.152e-2*te + 0.074) +
        1.28e-5*sin(3.152e-2*te + 1.547) + 3.14e-5*sin(2.1277e-1*te - 0.488)
    );

    const double nu = 9.282e-4*te - 0.8;
    const double dlam = 8.34e-5*sin(nu);
    const double lambda = l + PI + dlam;

    const double epsi = 4.089567e-1 - 6.19e-9*te + 4.46e-5*cos(nu);

    const double sl = sin(lambda);
    const double cl = cos(lambda);
    const double se = sin(epsi);
    const double ce = sqrt(1-se*se);

    *rasc = atan2(sl*ce, cl);
    if (*rasc < 0.0)
        *rasc += PI2;

    *decl = asin(sl*se);

    *hang = 1.7528311 + 6.300388099*t + lon - *rasc + 0.92*dlam;
    *hang = FMOD_FLOOR(*hang + PI, PI2) - PI;

    const double sp = sin(lat);
    const double cp = sqrt((1-sp*sp));
    const double sd = sin(*decl);
    const double cd = sqrt(1-sd*sd);
    const double sh = sin(*hang);
    const double ch = cos(*hang);
    const double se0 = sp*sd + cp*cd*ch;

    // elevation including parallax correction with radius in km
    const double ep = asin(se0) - sqrt(1.0-se0*se0)*rad*(1.0/149597871.0);

    // refraction correction
    double de = 0.0;
    if ((press > 0.0) && (ep > 0.0))
        de = (0.08422*press) / ((273.0+temp)*tan(ep + 0.003138/(ep + 0.08919)));

    *zenith = PIM - ep - de;
    *azimuth = atan2(sh, ch*sp - sd*cp/cd);
    *azimuth = FMOD_FLOOR(*azimuth, PI2) - PI;
}

/**
 * @brief Evaluate solar position in equatorial coordinate system.
 *
 * Algorithm #5 - part 1
 *
 * <NAME> (2012), Five new algorithms for the computation of sun
 * position from 2010 to 2110, Solar Energy, vol. 86(5), 1323-1337,
 * doi:10.1016/j.solener.2012.01.024
 * Ph.Blanc (2012), The SG2 algorithm for a fast and accurate computation
 * of the position of the Sun for multi-decadal time period, Solar Energy,
 * vol. 86(5), 3072-3083, doi:10.1016/j.solener.2012.07.018
 *
 *  outputs:
 *      decl - declination [rad]
 *      rasc - right ascension [rad] (eastward)
 *      hang - hour angle [rad] (westward!)
 *
 *  inputs:
 *      mjd2k  - Modified Julian Data 2000 [days]
 *      dtt - UTC offset to TT [seconds]
 *      lon - longitude [rad] (eastward)
 */

void sunpos5equat(
    double *decl, double *rasc, double *hang,
    double mjd2k, double dtt ,double lon
) {
    // t is decimal number of days since 1940-01-01
    const double t = mjd2k + 43.0 - 21958.0;
    const double te = SEC2DAYS*dtt + t; // Terestrial Time correction
    const double wte = 0.0172019715*te;
    const double s1 = sin(wte);
    const double c1 = cos(wte);
    const double s2 = 2.0*s1*c1;
    const double c2 = (c1+s1)*(c1-s1);
    const double s3 = s2*c1 + c2*s1;
    const double c3 = c2*c1 - s2*s1;

    const double l = (
        1.7527901 + 1.7202792159e-2*te + 3.33024e-2*s1 - 2.0582e-3*c1 +
        3.512e-4*s2 - 4.07e-5*c2 + 5.2e-6*s3 - 9e-7*c3 -
        8.23e-5*s1*sin(2.92e-5*te) + 1.27e-5*sin(1.49e-3*te - 2.337) +
        1.21e-5*sin(4.31e-3*te + 3.065) + 2.33e-5*sin(1.076e-2*te - 1.533) +
        3.49e-5*sin(1.575e-2*te - 2.358) + 2.67e-5*sin(2.152e-2*te + 0.074) +
        1.28e-5*sin(3.152e-2*te + 1.547) + 3.14e-5*sin(2.1277e-1*te - 0.488)
    );

    const double nu = 9.282e-4*te - 0.8;
    const double dlam = 8.34e-5*sin(nu);
    const double lambda = l + PI + dlam;
    const double epsi = 4.089567e-1 - 6.19e-9*te + 4.46e-5*cos(nu);

    const double sl = sin(lambda);
    const double cl = cos(lambda);
    const double se = sin(epsi);
    const double ce = sqrt(1-se*se);

    // equatorial coordinates
    const double _rasc0 = atan2(sl*ce, cl);
    const double _rasc = _rasc0 + (_rasc0 < 0.0 ? PI2 : 0.0);
    const double _decl = asin(sl*se);
    const double _hang0 = (
        // 1.7528311 + 6.300388099*t Greenwich Sideral Time approximation [rad]
        1.7528311 + 6.300388099*t + lon - _rasc + 0.92*dlam
    );
    const double _hang = FMOD_FLOOR(_hang0 + PI, PI2) - PI;

    *rasc = _rasc;
    *decl = _decl;
    *hang = _hang;
}

/**
 * @brief Convert solar equatorial to horizontal coordinates.
 *
 * Algorithm #5 part 2
 *
 * <NAME> (2012), Five new algorithms for the computation of sun
 * position from 2010 to 2110, Solar Energy, vol. 86(5), 1323-1337,
 * doi:10.1016/j.solener.2012.01.024
 * Ph.Blanc (2012), The SG2 algorithm for a fast and accurate computation
 * of the position of the Sun for multi-decadal time period, Solar Energy,
 * vol. 86(5), 3072-3083, doi:10.1016/j.solener.2012.07.018
 *
 *  outputs:
 *      azimuth - azimuth angle [rad]
 *      zenith - zenith angle [rad]
 *
 *  inputs:
 *      decl - declination [rad]
 *      hang - hour angle [rad] (westward!)
 *      lat - latitude [rad]
 *      rad - radius (distance from Earth centre) [km]
 *            (set to zero to disable paralax correction)
 *      pres - pressure [atm] (set 0 to skip refraction correction)
 *      temp - temperature [dgC]
 */

void sunpos5eq2hor(
    double *azimuth, double *zenith,
    double decl, double hang, double lat, double rad, double press, double temp
) {
    // horizontal coordinates
    // move conversion equatorial to horizontal coords to a separate function
    const double sp = sin(lat);
    const double cp = sqrt((1-sp*sp));
    const double sd = sin(decl);
    const double cd = sqrt(1-sd*sd);
    const double sh = sin(hang);
    const double ch = cos(hang);
    const double se0 = sp*sd + cp*cd*ch;

    // elevation including parallax correction with radius in km
    const double ep = asin(se0) - sqrt(1.0-se0*se0)*rad*(1.0/149597871.0);

    // refraction correction
    double de = 0.0;
    if ((press > 0.0) && (ep > 0.0))
        de = (0.08422*press) / ((273.0+temp)*tan(ep + 0.003138/(ep + 0.08919)));

    *zenith = PIM - ep - de;
    *azimuth = atan2(sh, ch*sp - sd*cp/cd);
    *azimuth = FMOD_FLOOR(*azimuth, PI2) - PI;
}

#endif /* SUNPOS_H */


<|start_filename|>geoist/magmod/pymm_relradpow.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings
 * - evaluation of the relative radius power series
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_RELRADPOW_H
#define PYMM_RELRADPOW_H

#include "shc.h"
#include "pymm_aux.h"

/* Earth radius in km */
#define RADIUS  6371.2

#ifndef STR
#define SRINGIFY(x) #x
#define STR(x) SRINGIFY(x)
#endif

/* python function definition */

#define DOC_RELRADPOW "\n"\
"   rrp = relradpow(radius, degree, reference_radius="STR(RADIUS)", is_internal=True)\n"\
"\n"\
"     By default when the 'is_internal' flag is set to True, evaluate\n"\
"     for the given 'radius' and 'reference_radius' relative radius power\n"\
"     series:\n"\
"       (reference_radius/radius)**(i+2) for i in range(0, degree+1) .\n"\
"\n"\
"     When the 'is_internal' flag is set to False, evaluate\n"\
"     for the given 'radius' and 'reference_radius' relative radius power\n"\
"     series:\n"\
"       (radius/reference_radius)**(i+1) for i in range(0, degree+1) .\n"\
"\n"


static PyObject* relradpow(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {
        "radius", "degree", "reference_radius", "is_internal", NULL
    };

    int degree;
    int is_internal;
    double rad, rad0 = RADIUS; // radius and reference radius
    PyObject *obj_is_internal = NULL; // boolean flag
    PyObject *arr_rrp = NULL; // P array
    PyObject *retval = NULL; // output tuple

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "di|dO:relradpow", keywords,
        &rad, &degree, &rad0, &obj_is_internal
    ))
        goto exit;

    is_internal = (obj_is_internal == NULL) || PyObject_IsTrue(obj_is_internal);

    if (degree < 0)
    {
        PyErr_Format(PyExc_ValueError, "%s < 0", keywords[1]);
        goto exit;
    }

    if (rad < 0.0)
    {
        PyErr_Format(PyExc_ValueError, "%s < 0", keywords[0]);
        goto exit;
    }

    if (rad0 <= 0.0)
    {
        PyErr_Format(PyExc_ValueError, "%s <= 0", keywords[2]);
        goto exit;
    }

    // create the output array
    if (NULL == (arr_rrp = _get_new_double_array(1, NULL, degree+1)))
        goto exit;

    // evaluate the relative radius power series
    if (is_internal)
        shc_relradpow_internal(PyArray_DATA(arr_rrp), degree, rad/rad0);
    else
        shc_relradpow_external(PyArray_DATA(arr_rrp), degree, rad/rad0);

    retval = arr_rrp;

  exit:

    // decrease reference counters to the arrays
    if (!retval && arr_rrp){Py_DECREF(arr_rrp);}

    return retval;
}

#endif  /* PYMM_RELRADPOW_H */


<|start_filename|>extensions/TsDeparture_zhangbei/tmpData/debug.json<|end_filename|>
﻿{
  "pars": {
    "Methods": "M"
  },
  "tsdata": "C:\\Users\\chens\\AppData\\Local\\Temp\\tsdata.txt",
  "OutputPath": "C:\\Users\\chens\\AppData\\Local\\Temp\\",
  "ResultFile": "C:\\Users\\chens\\AppData\\Local\\Temp\\result.json",
  "GetResult": "Print",
  "DataMode": "Url"
}

<|start_filename|>geoist/magmod/pytimeconv.c<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Time conversion utilities.
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2018 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#include "common.h" /* common definitions - to be included before Python.h */
#include <Python.h>
#include <numpy/arrayobject.h>

/* module version */
#include "version.h"

/* common python utilities */
#include "py_common.h"

/* time conversion subroutines */
#include "pytimeconv_mjd2000_to_decimal_year.h"
#include "pytimeconv_mjd2000_to_year_fraction.h"
#include "pytimeconv_decimal_year_to_mjd2000.h"

/*---------------------------------------------------------------------------*/
/* module's doc string */

#define DOC_PYTIMECONV \
"Time conversion utilities."

/*---------------------------------------------------------------------------*/
/*define module's methods */
static PyMethodDef pytimeconv_methods[] =
{
    {"mjd2000_to_decimal_year", (PyCFunction)pytimeconv_mjd2000_to_decimal_year, METH_VARARGS|METH_KEYWORDS, DOC_MJD2000_TO_DECIMAL_YEAR},
    {"mjd2000_to_year_fraction", (PyCFunction)pytimeconv_mjd2000_to_year_fraction, METH_VARARGS|METH_KEYWORDS, DOC_MJD2000_TO_YEAR_FRACTION},
    {"decimal_year_to_mjd2000", (PyCFunction)pytimeconv_decimal_year_to_mjd2000, METH_VARARGS|METH_KEYWORDS, DOC_DECIMAL_YEAR_TO_MJD2000},
    {NULL, NULL, 0, NULL} /* Sentinel - DO NOT REMOVE! */
} ;

/*---------------------------------------------------------------------------*/

/* module initialization  */
static PyObject* init_module(void)
{
    PyObject *module = init_python_module("_pytimeconv", DOC_PYTIMECONV, pytimeconv_methods);
    if (NULL == module)
        goto exit;

    PyObject *dict = PyModule_GetDict(module);
    if (NULL == dict)
        goto exit;

    /* metadata */
    set_dict_item_str_str(dict, "__author__", "<NAME> (<EMAIL>)");
    set_dict_item_str_str(dict, "__copyright__", "Copyright (C) 2018 EOX IT Services GmbH");
    set_dict_item_str_str(dict, "__licence__", "EOX licence (MIT style)");
    set_dict_item_str_str(dict, "__version__", VERSION);

  exit:
    return module;
}

#if PY_MAJOR_VERSION == 2

PyMODINIT_FUNC init_pytimeconv(void)
{
    import_array();
    init_module();
}

#else

PyObject* PyInit__pytimeconv(void)
{
    import_array();
    return init_module();
}

#endif


<|start_filename|>geoist/magmod/include/shc.h<|end_filename|>
/**
 * @file sph_harm.h
 * @author <NAME> <<EMAIL>>
 * @brief Spherical Harmonics
 *
 * Various utilities needed by the spherical harmonic model evaluation.
 *
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#ifndef SHC_H
#define SHC_H 1

#define SHC_LEG_HIGH_SCALE 1.0e+280
#define FDIV(a,b) ((double)(a)/(double)(b))

#include <math.h>

/**
 * @brief Allocate and evaluate pre-calculated square-root series
 *
 *   sqrt(i) for i = 0..(2*degree+1)
 *
 *   NOTE: The returned pointer *MUST* be passed to free().
 */

static double * shc_presqrt(int degree)
{
    const int nval = 2*degree + 2;
    int i;
    double * psqrt = malloc(nval*sizeof(double));

    if (NULL == psqrt)
        return NULL;

    for (i = 0; i < nval ; ++i)
    {
        psqrt[i] = sqrt((double) i);
    }

    return psqrt;
}

/**
 * @brief Evaluate the associated Legendre functions
 *
 * Evaluate the Schmidt semi-normalised associated Legendre functions.
 *
 * This version is susceptible to overflow for higher orders (degree > 20)
 * especially close to the poles.
 */

static void shc_legendre_low(double *lp, double *ldp, int degree, double sin_elv, double cos_elv, const double *psqrt)
{
    int i, j;

    // Gauss normalised associated Legendre functions (aLf)
    // Note the ldp sign change because of the derivative being calculated
    // with respect to the elevation (latitude) instead of the
    // inclination (co-latitude).
    lp[0] = 1.0;
    lp[1] = sin_elv ;
    lp[2] = cos_elv;

    ldp[0] = 0.0 ;
    ldp[1] = cos_elv;
    ldp[2] = -sin_elv;

    for (i = 2; i <= degree; ++i)
    {
        const int i_off = (i*(i+1))/2;
        const int i1sq = (i-1)*(i-1);
        const double w_dnm = 1.0/((2*i-1)*(2*i-3));

        for (j = 0; j < (i-1); ++j)
        {
            const int idx0 = i_off + j ;
            const int idx1 = i_off + j - i;
            const int idx2 = i_off + j - i - i + 1;
            const double w = w_dnm*(i1sq - j*j);

            lp [idx0] = sin_elv*lp [idx1] - w*lp [idx2];
            ldp[idx0] = sin_elv*ldp[idx1] - w*ldp[idx2] + cos_elv*lp[idx1];
        }

        {
            const int idx0 = i_off + i;
            const int idx1 = i_off - 1;

            lp [idx0-1] = sin_elv*lp [idx1];
            ldp[idx0-1] = sin_elv*ldp[idx1] + cos_elv*lp[idx1];

            lp [idx0] = cos_elv*lp [idx1];
            ldp[idx0] = cos_elv*ldp[idx1] - sin_elv*lp[idx1];
        }
    }

    // Convert Gauss normalised aLf to Schmidt quasi-normalised aLf
    {
        double sqn1, sqn0 = 1.0 ;

        for (i = 2; i <= degree; ++i)
        {
            const int i_off = (i*(i+1))/2;

            sqn0 *= FDIV(2*i-1,i);
            lp[i_off] *= sqn0;
            ldp[i_off] *= sqn0;

            sqn1 = sqn0*psqrt[2*i]/psqrt[i+1];

            lp[i_off+1] *= sqn1;
            ldp[i_off+1] *= sqn1;

            for (j = 2; j <= i; ++j)
            {
                sqn1 *= psqrt[i-j+1]/psqrt[i+j];
                lp[i_off+j] *= sqn1;
                ldp[i_off+j] *= sqn1;
            }
        }
    }
}


/**
 * @brief Evaluate the associated Legendre functions -
 *
 * Evaluate the Schmidt semi-normalised associated Legendre functions.
 *
 * This subroutine tries to compensate the effect of underflow for high degree
 * values near the poles by scaling of the of the powers of the .
 * [Holmes and Featherstone 2002, <NAME>, 76, 279-299]
 */


static void shc_legendre_high(double *lp, double *ldp, int degree, double sin_elv, double cos_elv, const double *psqrt)
{
    int i, j, idx, idx0;
    const double scale = SHC_LEG_HIGH_SCALE;
    const double rscale = 1.0/SHC_LEG_HIGH_SCALE;
    const double rcos_elv = 1.0/cos_elv ;
    double pm2, pm1, pm0, pmm;
    double pow_cos_elv;

    // Gauss normalised associated Legendre functions (aLf)
    lp[0] = pm2 = 1.0;
    lp[1] = pm1 = sin_elv;

    ldp[0] = 0.0 ;
    ldp[1] = cos_elv;

    for (i = 2, idx = 3; i <= degree; idx += ++i)
    {
        const double f1 = FDIV(2*i-1,i);
        const double f2 = FDIV(i-1,i);

        lp[idx] = pm0 = f1*sin_elv*pm1 - f2*pm2;
        ldp[idx] = i*rcos_elv*(pm1 - sin_elv*pm0);

        pm2 = pm1;
        pm1 = pm0;
    }

    pmm = rscale*psqrt[2];
    pow_cos_elv = scale;

    for (j = 1, idx0 = 2; j < degree; idx0 += ++j + 1)
    {
        const double pow_cos_elv_last = pow_cos_elv;
        pow_cos_elv *= cos_elv;

        // i = j
        pmm *= psqrt[2*j+1]/psqrt[2*j];
        pm2 = pmm/psqrt[2*j+1];
        pm1 = sin_elv*pmm;

        lp[idx0] = pow_cos_elv*pm2;
        ldp[idx0] = -pow_cos_elv_last*(j*sin_elv*pm2);

        // i = j + 1
        idx = idx0 + j + 1;
        lp[idx] = pow_cos_elv*pm1;
        ldp[idx] = -pow_cos_elv_last*((j+1)*sin_elv*pm1 - pmm);

        for(i = j + 2, idx += j + 2; i <= degree; idx += ++i)
        {
            const double f0 = psqrt[i+j]*psqrt[i-j];
            const double f1 = (2*i-1);
            const double f2 = psqrt[i-j-1]*psqrt[i+j-1];
            const double plm = (f1*sin_elv*pm1 - f2*pm2)/f0;

            lp[idx] = pow_cos_elv*plm;
            ldp[idx] = -pow_cos_elv_last*(i*sin_elv*plm - f0*pm1);

            pm2 = pm1;
            pm1 = plm;
        }
    }

    // i = j = degree
    pm2 = pmm/psqrt[2*j];

    lp[idx0] = pow_cos_elv*cos_elv*pm2;
    ldp[idx0] = -pow_cos_elv*(j*sin_elv*pm2);
}


/**
 * @brief Evaluate the associated Legendre functions
 *
 * Evaluate the Schmidt semi-normalised associated Legendre functions and
 * their derivatives with respect to the elevation (latitude) coordinate.
 *
 * parameters:
 *
 *      lp ... output array of the associated Legendre functions
 *      ldp ... output array of the associated Legendre functions' derivatives
 *      degree ... model degree
 *      elv ... elevation (latitude) angle in radians
 *      psqrt ... array of the pre-calculated square roots
 */

static void shc_legendre(double *lp, double *ldp, int degree, double elv, const double *psqrt)
{
    const double sin_elv = sin(elv);
    //const double cos_elv = cos(elv);
    const double cos_elv = sqrt((1.0-sin_elv)*(1.0+sin_elv));
    double *presqrt = shc_presqrt(degree);

    if ((degree <= 16)||(fabs(cos_elv) < 1e-10))
        // low degree model + poles
        shc_legendre_low(lp, ldp, degree, sin_elv, cos_elv, presqrt);
    else
        // high degree model
        shc_legendre_high(lp, ldp, degree, sin_elv, cos_elv, presqrt);

    free((void*)presqrt); // free allocated square roots
}

/**
 * @brief Evaluate the "internal" series of relative radius powers.
 *
 * Evaluate the series of relative radius powers for given relative radius ratio
 * 'relrad' (rad/rad0)
 *
 * (rad0/rad)**(i+2) for i = 0..degree
 *
 */

static void shc_relradpow_internal(double *rrp, int degree, double relrad)
{
    int i;
    double rr = 1.0/relrad;
    double rrp_last = rr;

    for (i = 0; i <= degree; ++i)
    {
        rrp_last *= rr;
        rrp[i] = rrp_last;
    }
}

/**
 * @brief Evaluate the "external" series of relative radius powers.
 *
 * Evaluate the series of relative radius powers for given relative radius ratio
 * 'relrad' (rad/rad0)
 *
 * (rad/rad0)**(i-1) for i = 0..degree
 *
 */

static void shc_relradpow_external(double *rrp, int degree, double relrad)
{
    int i;
    double rrp_val = 1.0;

    if (0 <= degree)
        rrp[0] = rrp_val/relrad;

    if (1 <= degree)
        rrp[1] = rrp_val;

    for (i = 2; i <= degree; ++i)
    {
        rrp_val *= relrad;
        rrp[i] = rrp_val;
    }
}

/**
 * @brief Evaluate the series of azimuth angle (longitude) sines and cosines.
 *
 * Evaluate the series of azimuth angle (longitude - 'lon') sines and cosines.
 *
 *   cos(i*lon) for i = 0...degree
 *   sin(i*lon) for i = 0...degree
 *
 * This subroutine uses a faster evaluation based on pure recurrent
 * addition/substration and multiplication:
 *
 *   sin(i*lon) = cos((i-1)*lon)*sin(lon) + sin((i-1)*lon)*cos(lon)
 *   cos(i*lon) = cos((i-1)*lon)*cos(lon) - sin((i-1)*lon)*sin(lon)
 *
 * The input angle must be in radians.
 */

static void shc_azmsincos(double *lonsin, double *loncos, int degree, double lon)
{
    int i;
    const double sin_lon = sin(lon);
    const double cos_lon = cos(lon);
    double sl, sl_new, cl, cl_new;

    lonsin[0] = 0.0;
    loncos[0] = 1.0;
    lonsin[1] = sl = sin_lon;
    loncos[1] = cl = cos_lon;

    for (i = 2; i <= degree; ++i)
    {
        sl_new = cl*sin_lon + sl*cos_lon;
        cl_new = cl*cos_lon - sl*sin_lon;
        lonsin[i] = sl = sl_new;
        loncos[i] = cl = cl_new;
    }
}

/**
 * @brief Evaluate the series of azimuth angle (longitude) sines and cosines.
 *
 * Evaluate the series of azimuth angle (longitude - 'lon') sines and cosines.
 *
 *   cos(i*lon) for i = 0...degree
 *   sin(i*lon) for i = 0...degree
 *
 * This subroutine contains the reference (slow) implementation evaluation
 * sine and cosine functions for each term of the series.
 *
 * The input angle must be in radians.
 */

static void shc_azmsincos_ref(double *lonsin, double *loncos, int degree, double lon)
{
    int i;
    for (i = 0; i <= degree; ++i)
    {
        const double ilon = i*lon;
        lonsin[i] = sin(ilon);
        loncos[i] = cos(ilon);
    }
}

/**
 * @brief Evaluate the scalar potential and its (spherical) gradient/
 *
 * Spherical harmonic evaluation of the scalar potential and
 * the gradient in the spherical coordinates.
 *
 *  outputs:
 *    vpot - value of the scalar potential
 *    dvel - elevation component of the gradient
 *    dvaz - azimuth component of the gradient
 *    dvrd - radial component of the gradient
 *
 *  inputs:
 *    degree - degree of the model
 *    mode - 3 - evaluate both gradient and potential
 *           1 - evaluate potential only
 *           2 - evaluate gradient only
 *    elv - elevation angle in radians (coordinate - needed in mode 3 and 2)
 *    rad - radius (coordinate needed in mode 3 and 1)
 *    cg, ch - spherical harmonic coefficients [(degree+1)*(degree+2)/2]
 *    lp, ldp - Legendre associative function and their derivatives (with
 *              respect to the elevation coordinate) [(degree+1)*(degree+2)/2]
 *    rrp - relative radius power series [degree+1]
 *    lsin, lcos - series of azimuth angle sines and cosines [degree+1]
 *    is_internal - boolean flag indicating type of the evaluated field
                    set true for the internal or false for the external field.
 */

static void shc_eval(
    double *vpot, double *dvel, double *dvaz, double *dvrd,
    int degree, int mode, double elv, double rad,
    const double *cg, const double *ch, const double *lp, const double *ldp,
    const double *rrp, const double *lsin, const double *lcos, int is_internal
)
// the evaluation
{
    int i, j, dr_scale, dr_offset;
    const double sin_elv = sin(elv);
    const double cos_elv = cos(elv);
    double _vpot, _dvrd, _dvel, _dvaz;

    if (is_internal) {
        dr_scale = 1;
        dr_offset = 1;
    } else {
        dr_scale = -1;
        dr_offset = 0;
    }

    { // i = 0
        const double tmp = cg[0]*rrp[0];
        const int i_dr = dr_scale*dr_offset;

        _vpot = tmp;
        _dvel = 0.0; // north-ward
        _dvaz = 0.0; // east-ward
        _dvrd = -tmp * i_dr ; // up-ward
    }

    for (i = 1; i <= degree; ++i)
    {
        const int i_off = (i*(i+1))/2;
        const int i_dr = dr_scale*(i + dr_offset);

        {
            const double tmp = cg[i_off]*rrp[i];

            _vpot += tmp * lp[i_off];
            _dvel += tmp * ldp[i_off]; // north-ward
            _dvrd -= tmp * lp[i_off] * i_dr ; // up-ward
        }

        for (j = 1; j <= i; ++j)
        {
            const int idx = i_off + j;
            const double tmp0 = (cg[idx]*lcos[j] + ch[idx]*lsin[j])*rrp[i];
            const double tmp1 = (cg[idx]*lsin[j] - ch[idx]*lcos[j])*rrp[i];

            _vpot += tmp0 * lp[idx];
            _dvel += tmp0 * ldp[idx]; // north-ward
            _dvaz -= tmp1 * lp[idx] * j; // east-ward
            _dvrd -= tmp0 * lp[idx] * i_dr ; // up-ward
        }
    }

    if (mode&0x1)
    {
        *vpot = _vpot * rad;
    }

    if (mode&0x2)
    {
        *dvel = _dvel;
        *dvaz = _dvaz / cos_elv;
        *dvrd = _dvrd;

        // handling of the poles
        if ((degree > 0) && (fabs(cos_elv) < 1e-10))
        {
            const double lsin1 = lsin[1], lcos1 = lcos[1];
            double sqn3, sqn1 = 1.0;
            double ps2, ps1 = 1.0, ps0 = 1.0,

            // i = 1
            _dvaz = -(cg[2]*lsin1 - ch[2]*lcos1) * rrp[1];

            for (i = 2; i <= degree; ++i)
            {
                const int idx = 1 + (i*(i+1))/2;
                const double tmp = FDIV((i-1)*(i-1)-1, (2*i-1)*(2*i-3));

                // evaluate ratio between the Gauss-normalised and Schmidt
                // quasi-normalised associated Legendre functions.
                //  Equivalent to: sqrt((j==0?1:2)*(i-j)!/(i+j!))*(2i-1)!!/(i-j)!
                sqn1 = sqn1 * FDIV(2*i-1, i);
                sqn3 = sqn1 * sqrt(FDIV(i*2, i+1));
                ps2 = ps1;
                ps1 = ps0;
                ps0 = sin_elv*ps1 - tmp*ps2;

                _dvaz -= (cg[idx]*lsin1 - ch[idx]*lcos1) * rrp[i] * ps0 * sqn3;
            }

            *dvaz = _dvaz;
        }
    }
}

#undef FDIV
#endif  /*SHC_H*/


<|start_filename|>geoist/magmod/pysunpos.c<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Solar position
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2017 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#include "common.h" /* common definitions - to be included before Python.h */
#include <Python.h>
#include <numpy/arrayobject.h>

/* module version */
#include "version.h"

/* common python utilities */
#include "py_common.h"

/* Sun ephemeris algorithms */
#include "pysunpos.h"
#include "pysunpos_original.h"

/*---------------------------------------------------------------------------*/
/* module's doc string */

#define DOC_PYSUNPOS \
"This module provides bindings to the solar model."

/*---------------------------------------------------------------------------*/
/*define module's methods */
static PyMethodDef pysunpos_methods[] =
{
    {"sunpos", (PyCFunction)pysunpos_sunpos, METH_VARARGS|METH_KEYWORDS, DOC_SUNPOS},
    {"sunpos_original", (PyCFunction)pysunpos_sunpos_original, METH_VARARGS|METH_KEYWORDS, DOC_SUNPOS_ORIGINAL},
    {NULL, NULL, 0, NULL} /* Sentinel - DO NOT REMOVE! */
} ;

/*---------------------------------------------------------------------------*/

static PyObject* init_module(void)
{
    PyObject *module = init_python_module("_pysunpos", DOC_PYSUNPOS, pysunpos_methods);
    if (NULL == module)
        goto exit;

    PyObject *dict = PyModule_GetDict(module);
    if (NULL == dict)
        goto exit;

    /* module metadata */
    set_dict_item_str_str(dict, "__author__", "<NAME> (<EMAIL>)");
    set_dict_item_str_str(dict, "__copyright__", "Copyright (C) 2017 EOX IT Services GmbH");
    set_dict_item_str_str(dict, "__licence__", "EOX licence (MIT style)");
    set_dict_item_str_str(dict, "__version__", VERSION);

  exit:
    return module;
}

#if PY_MAJOR_VERSION == 2

PyMODINIT_FUNC init_pysunpos(void)
{
    import_array();
    init_module();
}

#else

PyObject* PyInit__pysunpos(void)
{
    import_array();
    return init_module();
}

#endif


<|start_filename|>geoist/magmod/pytimeconv_decimal_year_to_mjd2000.h<|end_filename|>
/*-----------------------------------------------------------------------------
#
# MJD2000 to decimal year conversion
# 
#
# Author: <NAME> <<EMAIL>>
# 
# Original Author: <NAME> <<EMAIL>>
#-------------------------------------------------------------------------------
# Copyright (C) 2019 Geoist team
#
#-------------------------------------------------------------------------------*/


#ifndef PYTIMECONV_DECIMAL_YEAR_TO_MJD2000_H
#define PYTIMECONV_DECIMAL_YEAR_TO_MJD2000_H

#include "pymm_aux.h"
#include "time_conversion.h"
#include <math.h>

/*---------------------------------------------------------------------------*/
/* recursive array evaluation */

static void _year_to_mjd2k(ARRAY_DATA arrd_in, ARRAY_DATA arrd_out)
{
    if (arrd_in.ndim > 0)
    {
        npy_intp i;
        for(i = 0; i < arrd_in.dim[0]; ++i)
            _year_to_mjd2k(
                _get_arrd_item(&arrd_in, i),
                _get_arrd_item(&arrd_out, i)
            );
    }
    else
    {
        #define S(a) (*((double*)(a).data))
        S(arrd_out) = decimal_year_to_mjd2k(S(arrd_in));
        #undef S 
    }
}

#define DOC_DECIMAL_YEAR_TO_MJD2000 "\n"\
"   time_mjd2k = decimal_year_to_mjd2000(decimal_year)\n\n"\
"     Output:\n"\
"       time_mjd2k - array of the calculated MJD2000 times.\n"\
"\n"\
"     Parameters:\n"\
"       decimal_year - array of decimal years (up to 15 dimensions).\n"\
"\n"

static PyObject* pytimeconv_decimal_year_to_mjd2000(
    PyObject *self, PyObject *args, PyObject *kwdict
)
{
    static char *keywords[] = {"time_mjd2k", NULL};
    PyObject *obj_in = NULL; // input object
    PyObject *arr_in = NULL; // input array
    PyObject *arr_out = NULL; // return array
    PyObject *retval = NULL; // return object

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "O|:decimal_year_to_mjd2000", keywords, &obj_in
    ))
        goto exit;

    // cast the objects to arrays
    if (NULL == (arr_in = _get_as_double_array(obj_in, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr_in) > (MAX_OUT_ARRAY_NDIM-1))
    {
        PyErr_Format(PyExc_ValueError, "Array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[0],
            PyArray_NDIM(arr_in), (MAX_OUT_ARRAY_NDIM-1));
        goto exit;
    }

    // create output array (adding one extra dimension)
    if (NULL == (arr_out = PyArray_EMPTY(PyArray_NDIM(arr_in), PyArray_DIMS(arr_in), NPY_DOUBLE, 0)))
        goto exit;

    // evaluate solar positions
    _year_to_mjd2k(
        _array_to_arrd(arr_in),
        _array_to_arrd(arr_out)
    );

    // assign return value
    retval = arr_out;

  exit:

    // decrease reference counters to the arrays
    if (arr_in) Py_DECREF(arr_in);
    if (!retval && arr_out) Py_DECREF(arr_out);

    return retval;
}

#endif /* PYTIMECONV_DECIMAL_YEAR_TO_MJD2000_H */


<|start_filename|>geoist/flex/cpwt/cpwt.f90<|end_filename|>
! Copyright 2019 <NAME>
!
! This file is part of PlateFlex.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
!===========================================================================
!
! MODULE conf_cpwt
!
! Configuration module defining global variables used in modules to 
! interface with the Python codes.
!
!===========================================================================

      MODULE conf_cpwt

      IMPLICIT NONE

      REAL, PARAMETER :: pi = 3.141592653589793
      INTEGER, PARAMETER :: na = 11
!
! Wavelet parameter
!
      REAL :: k0

      END MODULE conf_cpwt

!===========================================================================
!
! MODULE cpwt
!
! Contains subroutines to calculate the wavelet transform, scalogram and 
! cross-spectral quantities with errors
!
!===========================================================================

      MODULE cpwt

      CONTAINS

!---------------------------------------------------------------------------
! Subroutine wlet_transform
!
! Subroutine to compute the complex-valued wavelet transform coefficients
! obtained from directional Morlet wavelets
!---------------------------------------------------------------------------

      SUBROUTINE wlet_transform(grid, nx, ny, dx, dy, kf, ns, wt_grid)     

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, nnx, nny, ns, nn(2)
      REAL    :: dx, dy
      REAL    :: grid(nx,ny), grid_mirror(2*nx, 2*ny), kf(ns)
      COMPLEX :: wt_grid(nx,ny,na,ns)

      INTEGER :: is, ia
      REAL    :: da, scales, angle, lam, kk
      CHARACTER(LEN=1) :: trail
      CHARACTER(LEN=30) :: progress

! Allocatable work arrays
      REAL, ALLOCATABLE :: grid_pad(:,:), kx(:), ky(:)
      COMPLEX, ALLOCATABLE :: ft_grid(:,:), ft2_grid(:,:), daughter(:,:) 
!
! Python bindings
!
!f2py REAL, intent(in) :: grid
!f2py INTEGER, intent(hide),depend(grid) :: nx=shape(grid,0), ny=shape(grid,1)
!f2py REAL, intent(in) :: dx, dy
!f2py REAL, intent(in) :: kf
!f2py INTEGER, intent(hide),depend(kf) :: ns=shape(kf,0)
!f2py COMPLEX, intent(out) :: wt_grid

        CALL npow2(2*nx, nnx)
        CALL npow2(2*ny, nny)
        ALLOCATE(grid_pad(nnx,nny))
! 
! Mirror data
!
        grid_mirror = 0.
        grid_mirror(1:nx,1:ny) = grid
        CALL mirror_data(grid_mirror,nx,ny,2*nx,2*ny)
! 
! Taper data
!
        ! CALL taper_data(grid_mirror,2*nx,2*ny,6)
!
! Pad with zeros and remove mean value
!
        grid_pad = 0.
        grid_pad(1:2*nx,1:2*ny) = grid_mirror - sum(grid_mirror)/2./2./nx/ny
!
! Complex form
!
        ALLOCATE(ft_grid(nnx,nny), ft2_grid(nnx,nny), daughter(nnx,nny))
        ft_grid = CMPLX(grid_pad)
!
! Fourier transform
!
        nn(1) = nnx
        nn(2) = nny
        CALL fourn(ft_grid,nn,2,1)
!
! Define wavenumbers
!
        ALLOCATE(kx(nnx), ky(nny))
        CALL defk(nnx,nny,dx,dy,kx,ky)
!
! Define angle parameters as per Kirby (2005)
!
        da = 2.*SQRT(-2.*LOG(0.75))/k0
!_______________________________________________
!     
!     MAIN WAVELET TANSFORM LOOP
!_______________________________________________
!
!
! Loop through scales
!
        trail = '-'
        progress = ''
        WRITE(*, fmt="(1x,a,i2,a)", advance="no")'#loops = ',ns,':'
        DO is = 1,ns
!
! Define wavenumbers and scales
!
          trail = '_'
          kk = kf(is)*1.e3
          lam = 2.*pi/kk
          scales = k0/kk
          ! WRITE(*, fmt="(1x,a)", advance="no")TRIM(progress)//trail
          WRITE(*, fmt="(1x,i2)", advance="no")is
          progress = TRIM(progress)//trail
!
! Loop through angles
!
          DO ia = 1,na
!
! Define angle increments
! 
            angle = REAL(ia-1)*da-pi/2.e0
!
! Calculate daughter wavelet
!
            CALL wave_function(nnx,nny,kx,ky,k0,scales,angle,daughter)
!
! Compute wavelet transform in Fourier space
!
            daughter = daughter*CONJG(ft_grid)
!
! Back to physical space
!
            CALL fourn(daughter,nn,2,-1)
!
! Normalization
!
            ft2_grid = (daughter)/FLOAT(nnx*nny)
!
! Quadrant swap
!
            CALL cfftswap(ft2_grid,nnx,nny)
!
! Store array
!
            wt_grid(:,:,ia,is) = ft2_grid(1:nx,1:ny)

          END DO
        END DO
        WRITE(*,*)

      END SUBROUTINE wlet_transform

!---------------------------------------------------------------------------
! Subroutine wlet_scalogram
!
! Subroutine to compute the real-valued wavelet scalogram obtained from 
! azimuthal-averaging of wavelet transform coefficients
!---------------------------------------------------------------------------

      SUBROUTINE wlet_scalogram(wl_trans, nx, ny, ns, wl_sg, ewl_sg)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, ns
      COMPLEX :: wl_trans(nx,ny,na,ns)
      REAL    :: wl_sg(nx,ny,ns), ewl_sg(nx,ny,ns)
      REAL    :: sg(nx,ny,na,ns)

      INTEGER :: is, ia
!
! Python bindings
!
!f2py COMPLEX, intent(in) :: wl_trans
!f2py INTEGER, intent(hide),depend(wl_trans) :: nx=shape(wl_trans,0)
!f2py INTEGER, intent(hide),depend(wl_trans) :: ny=shape(wl_trans,1)
!f2py INTEGER, intent(hide),depend(wl_trans) :: ns=shape(wl_trans,3)
!f2py REAL, intent(out) :: wl_sg, ewl_sg

        wl_sg = 0.e0
        DO is = 1,ns
          DO ia = 1,na
            sg(:,:,ia,is) = CABS(wl_trans(:,:,ia,is)*CONJG(wl_trans(:,:,ia,is)))
!
! Azimuthal average
!
            wl_sg(:,:,is) = wl_sg(:,:,is) + sg(:,:,ia,is)
          END DO
        END DO
!
! Normalization
!
        wl_sg = wl_sg/FLOAT(na)
!
! Jackknife estimation
!
        CALL jackknife_1d_sgram(sg,nx,ny,na,ns,ewl_sg)

        RETURN
      END SUBROUTINE wlet_scalogram

!---------------------------------------------------------------------------
! Subroutine wlet_xscalogram
!
! Subroutine to compute the complex-valued wavelet cross-scalogram obtained 
! from azimuthal-averaging of wavelet transform coefficients
!---------------------------------------------------------------------------

      SUBROUTINE wlet_xscalogram(wl_trans1, wl_trans2, nx, ny, ns, wl_xsg, ewl_xsg)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, ns
      COMPLEX :: wl_trans1(nx,ny,na,ns), wl_trans2(nx,ny,na,ns)
      COMPLEX :: wl_xsg(nx,ny,ns), xsg(nx,ny,na,ns)
      REAL :: ewl_xsg(nx,ny,ns), ewl_xsg_real(nx,ny,ns), ewl_xsg_imag(nx,ny,ns)

      INTEGER :: is, ia
!
! Python bindings
!
!f2py COMPLEX, intent(in) :: wl_trans1, wl_trans2
!f2py INTEGER, intent(hide),depend(wl_trans1) :: nx=shape(wl_trans1,0)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ny=shape(wl_trans1,1)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ns=shape(wl_trans1,3)
!f2py REAL, intent(out) :: wl_xsg, ewl_xsg

        wl_xsg = CMPLX(0.e0,0.e0)
        DO is = 1,ns
          DO ia = 1, na
            xsg(:,:,ia,is) = wl_trans1(:,:,ia,is)*CONJG(wl_trans2(:,:,ia,is))
!
! Azimuthal average
!
            wl_xsg(:,:,is) = wl_xsg(:,:,is) + xsg(:,:,ia,is)
          END DO
        END DO
!
! Normalization
!
        wl_xsg = wl_xsg/FLOAT(na)
!
! Jackknife estimation
!
        CALL jackknife_1d_sgram(REAL(xsg),nx,ny,na,ns,ewl_xsg_real)
        CALL jackknife_1d_sgram(IMAG(xsg),nx,ny,na,ns,ewl_xsg_imag)
        ewl_xsg = SQRT(ewl_xsg_real**2 + ewl_xsg_imag**2)

        RETURN
      END SUBROUTINE wlet_xscalogram

!---------------------------------------------------------------------------
! Subroutine wlet_admit_coh
!
! Subroutine to compute the complex-valued wavelet admittance and real-valued
! wavelet coherence obtained from azimuthal-averaging of wavelet transform 
! coefficients
!---------------------------------------------------------------------------

      SUBROUTINE wlet_admit_coh(wl_trans1, wl_trans2, nx, ny, ns, &
        wl_admit, ewl_admit, wl_coh, ewl_coh)

      USE conf_cpwt

      IMPLICIT NONE

      INTEGER :: nx, ny, ns
      COMPLEX :: wl_trans1(nx,ny,na,ns), wl_trans2(nx,ny,na,ns)
      REAL    :: wl_admit(nx,ny,ns), wl_coh(nx,ny,ns)
      REAL    :: ewl_admit(nx,ny,ns), ewl_coh(nx,ny,ns)

      REAL    :: sg1(nx,ny,na,ns), sg2(nx,ny,na,ns)
      COMPLEX :: xsg(nx,ny,na,ns)

      REAL    :: wl_sg1(nx,ny,ns), wl_sg2(nx,ny,ns)
      COMPLEX :: wl_xsg(nx,ny,ns)
      INTEGER :: is, ia
!
! Python bindings
!
!f2py COMPLEX, intent(in) :: wl_trans1, wl_trans2
!f2py INTEGER, intent(hide),depend(wl_trans1) :: nx=shape(wl_trans1,0)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ny=shape(wl_trans1,1)
!f2py INTEGER, intent(hide),depend(wl_trans1) :: ns=shape(wl_trans1,3)
!f2py REAL, intent(out) :: wl_admit, ewl_admit, wl_coh, ewl_coh

        DO is = 1,ns
          DO ia = 1, na

            sg1(:,:,ia,is) = CABS(wl_trans1(:,:,ia,is)*CONJG(wl_trans1(:,:,ia,is)))
            sg2(:,:,ia,is) = CABS(wl_trans2(:,:,ia,is)*CONJG(wl_trans2(:,:,ia,is)))
            xsg(:,:,ia,is) = wl_trans1(:,:,ia,is)*CONJG(wl_trans2(:,:,ia,is))
!
! Azimuthal average
!
            wl_sg1(:,:,is) = wl_sg1(:,:,is) + sg1(:,:,ia,is)
            wl_sg2(:,:,is) = wl_sg2(:,:,is) + sg2(:,:,ia,is)
            wl_xsg(:,:,is) = wl_xsg(:,:,is) + xsg(:,:,ia,is)

          END DO
        END DO
!
! Normalization
!
        wl_sg1 = wl_sg1/FLOAT(na)
        wl_sg2 = wl_sg2/FLOAT(na)
        wl_xsg = wl_xsg/FLOAT(na)
!
! Admittance and coherence (squared-real coherency) 
!
        wl_admit = REAL(wl_xsg)/wl_sg1
        wl_coh = REAL(wl_xsg)**2/(wl_sg1*wl_sg2)
        ! wl_coh = wl_xsg/SQRT(wl_sg1*wl_sg2)
!
! Jackknife estimation
!
        PRINT*,'Calculation jackknife error on admittance and coherence'
        CALL jackknife_admit_coh(xsg,sg1,sg2,nx,ny,na,ns,ewl_admit,ewl_coh)

        RETURN
      END SUBROUTINE wlet_admit_coh

      END MODULE cpwt

<|start_filename|>geoist/magmod/pymm_vrot_sph2geod.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - vector rotation - sph2geod
 *  (i.e., vector coordinate system transformation)
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_VROT_SPH2GEOD_H
#define PYMM_VROT_SPH2GEOD_H

#include "math_aux.h"
#include "geo_conv.h"
#include "pymm_aux.h"
#include "pymm_vrot_common.h"

/* recursive vector rotation */

static void _vrot_sph2geod(ARRAY_DATA ad_i, ARRAY_DATA ad_dlat,
                           ARRAY_DATA ad_o)
{
    if (ad_i.ndim > 1)
    {
        npy_intp i;
        for(i = 0; i < ad_i.dim[0]; ++i)
            _vrot_sph2geod(
                _get_arrd_item(&ad_i, i),
                _get_arrd_item(&ad_dlat, i),
                _get_arrd_item(&ad_o, i)
            );
    }
    else
    {
        #define SV(a) (*((double*)((a).data)))
        #define P(a,i) ((double*)((char*)(a).data+(i)*(a).stride[0])) //((double*)((a).data+(i)*(a).stride[0]))
        #define V(a,i) (*P(a,i))

        const double dlat = -DG2RAD*SV(ad_dlat);

        // rotate around the elevation axis
        rot2d(P(ad_o,2), P(ad_o,0), V(ad_i,2), V(ad_i,0), sin(dlat), cos(dlat));
        V(ad_o,1) = V(ad_i,1);

        #undef V
        #undef P
    }
}

/* python function definition */

#define DOC_VROT_SPH2GEOD "\n"\
"   arr_out = vrot_sph2geod(arr_in, arr_dlat)\n"\
"\n"\
"     Rotate vectors from the geocentric spherical (NEC) to the geodetic\n"\
"     (NEC) coordinate frame for a given difference of latitudes\n"\
"     in degrees.\n"\
"     This function can be also used for the inverse rotation from the geodetic\n"\
"     (NEC) to the geocentric spherical (NEC) coordinate frame by setting\n"\
"     the negative difference of latitudes.\n"\
"\n"\
"     The inputs are:\n"\
"         arr_in - array of the input vectors\n"\
"         arr_dlat - array of differences of the latitudes.\n"\
"     A scalar value is also accepted for a single vector rotation.\n"

static PyObject* vrot_sph2geod(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"arr_in", "arr_dlat", NULL};
    PyObject *obj_in = NULL; // input object
    PyObject *obj_dlat = NULL; // input object
    PyObject *arr_in = NULL; // input array
    PyObject *arr_dlat = NULL; // input array
    PyObject *arr_out = NULL; // output array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OO|:vrot_sph2geod", keywords, &obj_in, &obj_dlat
    ))
        goto exit;

    // cast the objects to arrays
    if (NULL == (arr_in=_get_as_double_array(obj_in, 1, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_dlat=_get_as_double_array(obj_dlat, 0, 0, NPY_ALIGNED, keywords[1])))
        goto exit;

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr_in) > MAX_OUT_ARRAY_NDIM)
    {
        PyErr_Format(PyExc_ValueError, "The input array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[0],
            PyArray_NDIM(arr_in), MAX_OUT_ARRAY_NDIM);
        goto exit;
    }

    // check the dimensions
    if (_check_array_dim_eq(arr_in, -1, 3, keywords[0]))
        goto exit;

    if (_vrot_arr_check(arr_in, arr_dlat, keywords[0], keywords[1]))
        goto exit;

    // create the output array
    if (NULL == (arr_out = _get_new_double_array(PyArray_NDIM(arr_in), PyArray_DIMS(arr_in), 3)))
        goto exit;

    // rotate the vector(s)
    _vrot_sph2geod(_array_to_arrd(arr_in), _array_to_arrd(arr_dlat), _array_to_arrd(arr_out));

    retval = arr_out;

  exit:

    // decrease reference counters to the arrays
    if (arr_in){Py_DECREF(arr_in);}
    if (arr_dlat){Py_DECREF(arr_dlat);}
    if (!retval && arr_out){Py_DECREF(arr_out);}

    return retval;
}

#endif  /* PYMM_VROT_SPH2GEOD_H */


<|start_filename|>geoist/magmod/pymm_spharpot.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - spherical-harmonics
 * - evaluation of the potential gradient in the spherical coordinate system
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_SPHARPOT_H
#define PYMM_SPHARPOT_H

#include "shc.h"
#include "pymm_aux.h"
#include "pymm_cconv.h"

/* python function definition */
#define DOC_SPHARPOT "\n"\
"  v_grad = spharpot(radius, degree, coef_g, coef_h, leg_p, rrp, lonsin, loncos)\n"\
"\n"\
"     Spherical harmonic evaluation of the value of the potential\n"\
"     (scalar) field in the (geocentric) spherical coordinates (latitude,\n"\
"     longitude, radius).\n"\
"     The input parameters are:\n"\
"       radius - radius, i.e., distance from the earth centre, at the\n"\
"                evaluated location.\n"\
"       degree - degree of the spherical harmonic model.\n"\
"       coef_g - vector of spherical harmonic model coefficients.\n"\
"       coef_h - vector of spherical harmonic model coefficients.\n"\
"       leg_p - vector the Legendre polynomials.\n"\
"       rrp - vector the relative radius powers.\n"\
"       lonsin - vector of the longitude cosines.\n"\
"       lonsin - vector of the longitude sines.\n"

static PyObject* spharpot(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"radius", "degree", "coef_g", "coef_h",
                                "leg_p", "rpp", "lonsin", "loncos", NULL};

    int degree, nterm;
    double rad;

    PyObject *obj_cg = NULL; // coef_g object
    PyObject *obj_ch = NULL; // coef_h object
    PyObject *obj_lp = NULL; // P object
    PyObject *obj_rrp = NULL; // rel.rad.pow. object
    PyObject *obj_lsin = NULL; // lonsin object
    PyObject *obj_lcos = NULL; // loncos object

    PyObject *arr_out = NULL; // output array
    PyObject *arr_cg = NULL; // coef_g array
    PyObject *arr_ch = NULL; // coef_h array
    PyObject *arr_lp = NULL; // P array
    PyObject *arr_rrp = NULL; // rel.rad.pow. array
    PyObject *arr_lsin = NULL; // lonsin array
    PyObject *arr_lcos = NULL; // loncos array

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "diOOOOOO|:spharpot", keywords,
        &rad, &degree, &obj_cg, &obj_ch, &obj_lp, &obj_rrp, &obj_lsin, &obj_lcos
    ))
        goto exit;

    if (degree < 0)
    {
        PyErr_Format(PyExc_ValueError, "%s < 0", keywords[1]);
        goto exit;
    }

    nterm = ((degree+1)*(degree+2))/2;

    // cast the objects to arrays
    if (NULL == (arr_cg=_get_as_double_array(obj_cg, 1, 1, NPY_IN_ARRAY, keywords[2])))
        goto exit;

    if (NULL == (arr_ch=_get_as_double_array(obj_ch, 1, 1, NPY_IN_ARRAY, keywords[3])))
        goto exit;

    if (NULL == (arr_lp=_get_as_double_array(obj_lp, 1, 1, NPY_IN_ARRAY, keywords[4])))
        goto exit;

    if (NULL == (arr_rrp=_get_as_double_array(obj_rrp, 1, 1, NPY_IN_ARRAY, keywords[5])))
        goto exit;

    if (NULL == (arr_lsin=_get_as_double_array(obj_lsin, 1, 1, NPY_IN_ARRAY, keywords[6])))
        goto exit;

    if (NULL == (arr_lcos=_get_as_double_array(obj_lcos, 1, 1, NPY_IN_ARRAY, keywords[7])))
        goto exit;

    // check the arrays' dimensions
    if (_check_array_dim_le(arr_cg, 0, nterm, keywords[2]))
        goto exit;

    if (_check_array_dim_le(arr_ch, 0, nterm, keywords[3]))
        goto exit;

    if (_check_array_dim_le(arr_lp, 0, nterm, keywords[4]))
        goto exit;

    if (_check_array_dim_le(arr_rrp, 0, degree+1, keywords[5]))
        goto exit;

    if (_check_array_dim_le(arr_lsin, 0, degree+1, keywords[6]))
        goto exit;

    if (_check_array_dim_le(arr_lcos, 0, degree+1, keywords[7]))
        goto exit;

    // allocate the output numpy scalar (zero-dimensional array)
    if (NULL == (arr_out = PyArray_EMPTY(0, NULL, NPY_DOUBLE, 0)))
        goto exit;

    // the evaluation
    {
        double *out = (double*)PyArray_DATA(arr_out);
        shc_eval(
            out, NULL, NULL, NULL, degree, 0x1,
            0.0, rad, PyArray_DATA(arr_cg),
            PyArray_DATA(arr_ch), PyArray_DATA(arr_lp),
            PyArray_DATA(arr_lp), PyArray_DATA(arr_rrp),
            PyArray_DATA(arr_lsin), PyArray_DATA(arr_lcos),
            0 // this flag has no impact on the evaluated potential
        );
    }

  exit:

    // decrease reference counters to the arrays
    if (arr_cg){Py_DECREF(arr_cg);}
    if (arr_ch){Py_DECREF(arr_ch);}
    if (arr_lp){Py_DECREF(arr_lp);}
    if (arr_rrp){Py_DECREF(arr_rrp);}
    if (arr_lsin){Py_DECREF(arr_lsin);}
    if (arr_lcos){Py_DECREF(arr_lcos);}

    return arr_out;
 }

#endif  /* PYMM_SPHARPOT_H */


<|start_filename|>geoist/magmod/pyqd_eval_mlt.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Magnetic Quasi Dipole Coordinates - C python bindings
 * -  Magnetic Local Time evaluation
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2016 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYQD_EVAL_MLT_H
#define PYQD_EVAL_MLT_H

#include "pymm_aux.h"
#include "qdipole/cqdipole.h"

/* python function definition */

#define DOC_EVAL_MLT "\n"\
"   mlt = eval_mlt(qdlon, time, fname)\n"\
"     Inputs:\n"\
"       qdlon - quasi-dipole longitudes(s).\n"\
"       time  - MJD2000 time(s)\n"\
"       fname - file-name of the model text file.\n"\
"     Outputs:\n"\
"       mlt - magnetic local times (s).\n"\
""

static PyObject* eval_mlt(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"qdlon", "time", "fname", NULL};

    const char *model_fname = NULL;
    PyObject *obj_qdlon = NULL; // gclon object
    PyObject *obj_time = NULL; // time object
    PyObject *arr_qdlon = NULL; // qdlon array
    PyObject *arr_time = NULL; // time array
    PyObject *arr_mlt = NULL; // mlt array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OOs:eval_mlt", keywords,
        &obj_qdlon, &obj_time, &model_fname
    ))
        goto exit;

    #define NPY_REQ (NPY_ALIGNED|NPY_CONTIGUOUS)

    // cast the objects to arrays
    if (NULL == (arr_qdlon=_get_as_double_array(obj_qdlon, 0, 1, NPY_REQ, keywords[1])))
        goto exit;

    if (NULL == (arr_time=_get_as_double_array(obj_time, 0, 1, NPY_REQ, keywords[3])))
        goto exit;

    // check the dimensions
    npy_intp ndim = PyArray_NDIM(arr_qdlon);
    npy_intp *dims = PyArray_DIMS(arr_qdlon);

    if(_check_arr_dims_all_eq(arr_time, ndim, dims, keywords[3]))
        goto exit;

    // create the output arrays
    if (NULL == (arr_mlt = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    // evaluate the output values
    c_eval_mlt(
        (double*)PyArray_DATA(arr_mlt),
        (double*)PyArray_DATA(arr_qdlon),
        (double*)PyArray_DATA(arr_time),
        ndim == 0 ? 1 : dims[0],
        model_fname
    );

    retval = arr_mlt;

  exit:

    // decrease reference counters to the arrays
    if (arr_qdlon){Py_DECREF(arr_qdlon);}
    if (arr_time){Py_DECREF(arr_time);}
    if (!retval)
    {
        if (arr_mlt){Py_DECREF(arr_mlt);}
    }

    return retval;
}

#endif  /* PYQD_EVAL_MLT_H */



<|start_filename|>geoist/magmod/pyqd_eval_subsol.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Magnetic Quasi Dipole Coordinates - C python bindings
 * -  sub-solar point evaluation
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2016 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYQD_EVAL_SUBSOL_H
#define PYQD_EVAL_SUBSOL_H

#include "pymm_aux.h"
#include "qdipole/cqdipole.h"

/* python function definition */

#define DOC_EVAL_SUBSOL "\n"\
"   gdlat, gdlon = eval_subsol(time)\n"\
"     Inputs:\n"\
"       time  - MJD2000 time(s)\n"\
"     Outputs:\n"\
"       gdlat - sub-solar point latitude(s).\n"\
"       gdlon - sub-solar point longitudes(s).\n"\
""

static PyObject* eval_subsol(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"time", NULL};

    PyObject *obj_time = NULL; // time object
    PyObject *arr_time = NULL; // time array
    PyObject *arr_gdlat = NULL; // gdlat array
    PyObject *arr_gdlon = NULL; // gdlon array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "O:eval_subsol", keywords, &obj_time
    ))
        goto exit;

    #define NPY_REQ (NPY_ALIGNED|NPY_CONTIGUOUS)

    // cast the objects to arrays
    if (NULL == (arr_time=_get_as_double_array(obj_time, 0, 1, NPY_REQ, keywords[0])))
        goto exit;

    // create the output arrays
    npy_intp ndim = PyArray_NDIM(arr_time);
    npy_intp *dims = PyArray_DIMS(arr_time);

    if (NULL == (arr_gdlat = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_gdlon = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    // evaluate the output values
    c_eval_subsol(
        (double*)PyArray_DATA(arr_gdlat),
        (double*)PyArray_DATA(arr_gdlon),
        (double*)PyArray_DATA(arr_time),
        ndim == 0 ? 1 : dims[0]
    );

    if (NULL == (retval = Py_BuildValue("NN", arr_gdlat, arr_gdlon))) goto exit;

  exit:

    // decrease reference counters to the arrays
    if (arr_time) {Py_DECREF(arr_time);}
    if (!retval)
    {
        if (arr_gdlat) {Py_DECREF(arr_gdlat);}
        if (arr_gdlon) {Py_DECREF(arr_gdlon);}
    }

    return retval;
}

#endif  /* PYQD_EVAL_SUBSOL_H */




<|start_filename|>geoist/magmod/pyqd_eval_qdlatlon.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Magnetic Quasi Dipole Coordinates - C python bindings
 * -  QD coordinates' evaluation
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2016 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYQD_EVAL_QDLATLON_H
#define PYQD_EVAL_QDLATLON_H

#include "pymm_aux.h"
#include "qdipole/cqdipole.h"

/* python function definition */

#define DOC_EVAL_QDLATLON "\n"\
"   qdlat, qdlon = eval_qdlatlon(gclat, gclon, gcrad, time, fname)\n"\
"   qdlat, qdlon, f11, f12, f21, f22, f__ = "\
                       "eval_qdlatlon(gclat, gclon, gcrad, time fname, True)\n"\
"     Inputs:\n"\
"       gclat - geocentric latitude(s).\n"\
"       gclon - geocentric longitudes(s).\n"\
"       gcrad - geocentric radial coordinate(s) in km.\n"\
"       time  - decimal year time(s)\n"\
"       fname - file-name of the model text file.\n"\
"       eval_base - boolean flag indicating whether the base vectors should\n"\
"                   be calculated or not (False by default).\n"\
"     Outputs:\n"\
"       qdlat - quasi-dipole latitude(s).\n"\
"       qdlon - quasi-dipole longitudes(s).\n"\
"       f11, f12, f21, f22 - base vectors.\n"\
"       f__    - F = |F1 x F2|\n"\
""

static PyObject* eval_qdlatlon(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {
        "gclat", "gclon", "gcrad", "time", "fname", "eval_base", NULL
    };

    int eval_base;
    const char *model_fname = NULL;
    PyObject *obj_gclat = NULL; // gclat object
    PyObject *obj_gclon = NULL; // gclon object
    PyObject *obj_gcrad = NULL; // gcrad object
    PyObject *obj_time = NULL; // time object
    PyObject *obj_flag = NULL; //  boolean flag object

    PyObject *arr_gclat = NULL; // gclat array
    PyObject *arr_gclon = NULL; // gclon array
    PyObject *arr_gcrad = NULL; // gcrad array
    PyObject *arr_time = NULL; // time array

    PyObject *arr_qdlat = NULL; // qdlat array
    PyObject *arr_qdlon = NULL; // qdlon array
    PyObject *arr_f11 = NULL; // f11 array
    PyObject *arr_f12 = NULL; // f12 array
    PyObject *arr_f21 = NULL; // f21 array
    PyObject *arr_f22 = NULL; // f22 array
    PyObject *arr_f = NULL; // f array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OOOOs|O:eval_qdlatlon", keywords,
        &obj_gclat, &obj_gclon, &obj_gcrad, &obj_time, &model_fname, &obj_flag
    ))
        goto exit;

    eval_base = obj_flag ? (1 == PyObject_IsTrue(obj_flag)) : 0;

    #define NPY_REQ (NPY_ALIGNED|NPY_CONTIGUOUS)

    // cast the objects to arrays
    if (NULL == (arr_gclat=_get_as_double_array(obj_gclat, 0, 1, NPY_REQ, keywords[0])))
        goto exit;

    if (NULL == (arr_gclon=_get_as_double_array(obj_gclon, 0, 1, NPY_REQ, keywords[1])))
        goto exit;

    if (NULL == (arr_gcrad=_get_as_double_array(obj_gcrad, 0, 1, NPY_REQ, keywords[2])))
        goto exit;

    if (NULL == (arr_time=_get_as_double_array(obj_time, 0, 1, NPY_REQ, keywords[3])))
        goto exit;

    // check the dimensions
    npy_intp ndim = PyArray_NDIM(arr_gclat);
    npy_intp *dims = PyArray_DIMS(arr_gclat);

    if (_check_arr_dims_all_eq(arr_gclon, ndim, dims, keywords[1]))
        goto exit;

    if (_check_arr_dims_all_eq(arr_gcrad, ndim, dims, keywords[2]))
        goto exit;

    if (_check_arr_dims_all_eq(arr_time, ndim, dims, keywords[3]))
        goto exit;

    // create the output arrays
    if (NULL == (arr_qdlat = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_qdlon = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (eval_base)
    {
        if (NULL == (arr_f11 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
            goto exit;

        if (NULL == (arr_f12 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
            goto exit;

        if (NULL == (arr_f21 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
            goto exit;

        if (NULL == (arr_f22 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
            goto exit;

        if (NULL == (arr_f = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
            goto exit;
    }

    // evaluate the output values
    if (eval_base)
    {
        c_eval_qdlatlonvb(
            (double*)PyArray_DATA(arr_qdlat),
            (double*)PyArray_DATA(arr_qdlon),
            (double*)PyArray_DATA(arr_f11),
            (double*)PyArray_DATA(arr_f12),
            (double*)PyArray_DATA(arr_f21),
            (double*)PyArray_DATA(arr_f22),
            (double*)PyArray_DATA(arr_f),
            (double*)PyArray_DATA(arr_time),
            (double*)PyArray_DATA(arr_gcrad),
            (double*)PyArray_DATA(arr_gclat),
            (double*)PyArray_DATA(arr_gclon),
            ndim == 0 ? 1 : dims[0],
            model_fname
        );
    }
    else
    {
        c_eval_qdlatlon(
            (double*)PyArray_DATA(arr_qdlat),
            (double*)PyArray_DATA(arr_qdlon),
            (double*)PyArray_DATA(arr_time),
            (double*)PyArray_DATA(arr_gcrad),
            (double*)PyArray_DATA(arr_gclat),
            (double*)PyArray_DATA(arr_gclon),
            ndim == 0 ? 1 : dims[0],
            model_fname
        );
    }

    if (eval_base)
    {
        retval = Py_BuildValue(
            "NNNNNNN", arr_qdlat, arr_qdlon,
            arr_f11, arr_f12, arr_f21, arr_f22, arr_f
        );
    }
    else
    {
        retval = Py_BuildValue("NN", arr_qdlat, arr_qdlon);
    }

    if (NULL == retval) goto exit;

  exit:

    // decrease reference counters to the arrays
    if (arr_gclat) {Py_DECREF(arr_gclat);}
    if (arr_gclon) {Py_DECREF(arr_gclon);}
    if (arr_gcrad) {Py_DECREF(arr_gcrad);}
    if (arr_time) {Py_DECREF(arr_time);}
    if (!retval)
    {
        if (arr_qdlat) {Py_DECREF(arr_qdlat);}
        if (arr_qdlon) {Py_DECREF(arr_qdlon);}
        if (arr_f11) {Py_DECREF(arr_f11);}
        if (arr_f12) {Py_DECREF(arr_f12);}
        if (arr_f21) {Py_DECREF(arr_f21);}
        if (arr_f22) {Py_DECREF(arr_f22);}
        if (arr_f) {Py_DECREF(arr_f);}
    }

    return retval;
}

#endif  /* PYQD_EVAL_QDLATLON_H */



<|start_filename|>extensions/RandGen_chenshi/tmpData/debug.json<|end_filename|>
﻿{
  "pars": {
    "Chan1": "正态分布",
    "Chan2": "泊松分布",
    "Chan3": "0-1分布",
    "windowsize": 100
  },
  "OutputPath": "C:\\Users\\chens\\AppData\\Local\\Temp\\",
  "ResultFile": "C:\\Users\\chens\\AppData\\Local\\Temp\\result.json",
  "GetResult": "Print",
  "DataMode": "Url"
}

<|start_filename|>geoist/magmod/pyqd_eval_apex.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Magnetic Quasi Dipole Coordinates - C python bindings
 * -  QD coordinates' evaluation
 *
 * NOTE: This subroutine does not return correct MLT values!
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2015 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYQD_EVAL_APEX_H
#define PYQD_EVAL_APEX_H

#include "pymm_aux.h"
#include "qdipole/cqdipole.h"

/* python function definition */

#define DOC_EVAL_APEX "\n"\
"   qdlat, qdlon, mlt, f11, f12, f21, f22 = "\
                                "eval_apex(gclat, gclon, gcrad, time, fname)\n"\
"     Inputs:\n"\
"       gclat - geocentric latitude(s).\n"\
"       gclon - geocentric longitudes(s).\n"\
"       gcrad - geocentric radial coordinate(s) in km.\n"\
"       time  - decimal year time(s)\n"\
"       fname - file-name of the model text file.\n"\
"     Outputs:\n"\
"       qdlat - quasi-dipole latitude(s).\n"\
"       qdlon - quasi-dipole longitudes(s).\n"\
"       mlt - magnetic local times (s).\n"\
"       f11, f12, f21, f22 - base vectors.\n"\
""

static PyObject* eval_apex(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"gclat", "gclon", "gcrad", "time", "fname", NULL};

    const char *model_fname = NULL;
    PyObject *obj_gclat = NULL; // gclat object
    PyObject *obj_gclon = NULL; // gclon object
    PyObject *obj_gcrad = NULL; // gcrad object
    PyObject *obj_time = NULL; // time object
    PyObject *arr_gclat = NULL; // gclat array
    PyObject *arr_gclon = NULL; // gclon array
    PyObject *arr_gcrad = NULL; // gcrad array
    PyObject *arr_time = NULL; // time array
    PyObject *arr_qdlat = NULL; // qdlat array
    PyObject *arr_qdlon = NULL; // qdlon array
    PyObject *arr_mlt = NULL; // mlt array
    PyObject *arr_f11 = NULL; // f11 array
    PyObject *arr_f12 = NULL; // f12 array
    PyObject *arr_f21 = NULL; // f21 array
    PyObject *arr_f22 = NULL; // f22 array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OOOOs:eval_apex", keywords,
        &obj_gclat, &obj_gclon, &obj_gcrad, &obj_time, &model_fname
    ))
        goto exit;

    #define NPY_REQ (NPY_ALIGNED|NPY_CONTIGUOUS)

    // cast the objects to arrays
    if (NULL == (arr_gclat=_get_as_double_array(obj_gclat, 0, 1, NPY_REQ, keywords[0])))
        goto exit;

    if (NULL == (arr_gclon=_get_as_double_array(obj_gclon, 0, 1, NPY_REQ, keywords[1])))
        goto exit;

    if (NULL == (arr_gcrad=_get_as_double_array(obj_gcrad, 0, 1, NPY_REQ, keywords[2])))
        goto exit;

    if (NULL == (arr_time=_get_as_double_array(obj_time, 0, 1, NPY_REQ, keywords[3])))
        goto exit;

    // check the dimensions
    npy_intp ndim = PyArray_NDIM(arr_gclat);
    npy_intp *dims = PyArray_DIMS(arr_gclat);

    if(_check_arr_dims_all_eq(arr_gclon, ndim, dims, keywords[1]))
        goto exit;

    if(_check_arr_dims_all_eq(arr_gcrad, ndim, dims, keywords[2]))
        goto exit;

    if(_check_arr_dims_all_eq(arr_time, ndim, dims, keywords[3]))
        goto exit;

    // create the output arrays
    if (NULL == (arr_qdlat = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_qdlon = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_mlt = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_f11 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_f12 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_f21 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;

    if (NULL == (arr_f22 = PyArray_EMPTY(ndim, dims, NPY_DOUBLE, 0)))
        goto exit;


    // evaluate the output values
    c_make_apex(
        (double*)PyArray_DATA(arr_qdlat),
        (double*)PyArray_DATA(arr_qdlon),
        (double*)PyArray_DATA(arr_mlt),
        (double*)PyArray_DATA(arr_f11),
        (double*)PyArray_DATA(arr_f12),
        (double*)PyArray_DATA(arr_f21),
        (double*)PyArray_DATA(arr_f22),
        (double*)PyArray_DATA(arr_time),
        (double*)PyArray_DATA(arr_gcrad),
        (double*)PyArray_DATA(arr_gclat),
        (double*)PyArray_DATA(arr_gclon),
        ndim == 0 ? 1 : dims[0],
        model_fname
    );

    if (NULL == (retval = Py_BuildValue("NNNNNNN", arr_qdlat, arr_qdlon,
            arr_mlt, arr_f11, arr_f12, arr_f21, arr_f22)))
        goto exit;

  exit:

    // decrease reference counters to the arrays
    if (arr_gclat){Py_DECREF(arr_gclat);}
    if (arr_gclon){Py_DECREF(arr_gclon);}
    if (arr_gcrad){Py_DECREF(arr_gcrad);}
    if (arr_time){Py_DECREF(arr_time);}
    if (!retval)
    {
        if (arr_qdlat){Py_DECREF(arr_qdlat);}
        if (arr_qdlon){Py_DECREF(arr_qdlon);}
        if (arr_mlt){Py_DECREF(arr_mlt);}
        if (arr_f11){Py_DECREF(arr_f11);}
        if (arr_f12){Py_DECREF(arr_f12);}
        if (arr_f21){Py_DECREF(arr_f21);}
        if (arr_f22){Py_DECREF(arr_f22);}
    }

    return retval;
}

#endif  /* PYQD_EVAL_APEX_H */


<|start_filename|>geoist/magmod/pymm_legendre.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings
 * - associative Legendre functions evaluation
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_LEGENDRE_H
#define PYMM_LEGENDRE_H

#include "shc.h"
#include "pymm_aux.h"
#include "pymm_cconv.h"

/* python function definition */

#define DOC_LEGENDRE "\n"\
"   p, dp = legendre(latitude, degree)\n"\
"\n"\
"     For given the 'latitude' in degrees and the model's 'degree' evaluate\n"\
"     the Schmidt semi-normalised associated Legendre functions and its\n"\
"     and its derivatives: \n"\
"         P_n^m(sin(latitude))  and  dP_n^m(sin(latitude))\n"\
"     where n = 0..degree and m = 0..n \n"\
"\n"\
"     The input parameters are:\n"\
"       latitude - spherical latitude in degrees.\n"\
"       degree - degree of the spherical harmonic model.\n"\
"\n"


static PyObject* legendre(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"latitude", "degree", NULL};

    int degree, nterm;
    double lat_sph;
    PyObject *arr_p = NULL; // P array
    PyObject *arr_dp = NULL; // dP array
    PyObject *retval = NULL; // output tuple

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "di:legendre", keywords, &lat_sph, &degree
    ))
        goto exit;

    if (degree < 0)
    {
        PyErr_Format(PyExc_ValueError, "%s < 0", keywords[1]);
        goto exit;
    }

    if ((lat_sph < -90.0) || (lat_sph > 90.0))
    {
        if (lat_sph < 0.0) {
            PyErr_Format(PyExc_ValueError, "%s < -90", keywords[0]);
        } else {
            PyErr_Format(PyExc_ValueError, "%s > 90", keywords[0]);
        }
        goto exit;
    }

    nterm = ((degree+1)*(degree+2))/2;

    // create the output arrays
    if (NULL == (arr_p = _get_new_double_array(1, NULL, nterm)))
        goto exit;

    if (NULL == (arr_dp = _get_new_double_array(1, NULL, nterm)))
        goto exit;

    {
        // allocate and fill the pre-calculated square-roots
        double *psqrt = NULL;
        if (NULL == (psqrt = shc_presqrt(degree)))
        {
            PyErr_Format(PyExc_ValueError, "Memory allocation error!");
            goto exit;
        }

        shc_legendre(PyArray_DATA(arr_p), PyArray_DATA(arr_dp), degree, DG2RAD*lat_sph, psqrt);

        // free the square root array
        free(psqrt);
    }

    if (NULL == (retval = Py_BuildValue("NN", arr_p, arr_dp)))
        goto exit;

  exit:

    // decrease reference counters to the arrays
    if (!retval && arr_p){Py_DECREF(arr_p);}
    if (!retval && arr_dp){Py_DECREF(arr_dp);}

    return retval;
}

#endif  /* PYMM_LEGENDRE_H */



<|start_filename|>parameters.json<|end_filename|>
{
   "resource": "https://graph.microsoft.com/",
   "tenant" : "{your_directoryId}",
   "authorityHostUrl" : "https://login.microsoftonline.com",
   "clientid" : "{your_app_id}",
   "username" : "{username}",
   "password" : <PASSWORD>}"
 }


<|start_filename|>geoist/flex/cpwt/cpwt_sub.f90<|end_filename|>
!===============================================================
! SUBROUTINES cpwt_sub
!
!   Provies subroutines to CPWT program 
!
!   (See MakeFILE for compilation)
!
!===============================================================
!   <NAME>
!   Department of Earth and Planetary Science
!   U! Berkeley
!   <EMAIL>      
!
!   July 2009
!================================================================
   
!_____________________________________________________________
!
!     SUBROUTINE LIBRARY
!_____________________________________________________________
!

!--------------------------------------------
!     Function determines smallest power of 2
!--------------------------------------------
      SUBROUTINE npow2(n, n2)
      n2=1
      DO 
        n2=2*n2
        IF(n2.ge.n)RETURN
      END DO
      RETURN 
      END SUBROUTINE npow2


!------------------------------------------------
!     Subroutine defk
!
!     Computes wavenumbers for given nx,ny,dx,dy
!------------------------------------------------
      SUBROUTINE defk(nx,ny,dx,dy,kx,ky)
      IMPLICIT NONE
      REAL :: pi = 3.141592653589793
      INTEGER :: nx,ny,ix,iy
      REAL :: kx(nx),ky(ny),dx,dy,dkx,dky

      dkx=2*pi/(FLOAT(nx)*dx)
      dky=2*pi/(FLOAT(ny)*dy)
      DO iy=1,ny/2+1
        DO ix=1,nx/2+1
          kx(ix)=FLOAT(ix-1)*dkx
          ky(iy)=FLOAT(iy-1)*dky
        END DO
      END DO
      DO iy=ny/2+2,ny
        DO ix=nx/2+2,nx
          kx(ix)=-kx(nx-ix+2)
          ky(iy)=-ky(ny-iy+2)
        END DO
      END DO
      RETURN
      END SUBROUTINE defk


!---------------------------------------------
!     Subroutine wave_function
!
!     !omputes the daughter Morlet wavelet 
!     in Fourier space (kx,ky) with a 
!     preferred orientation in wavevector space
!     given by 'a', at scale 's'
!---------------------------------------------
      SUBROUTINE wave_function(nx, ny, kx, ky, k0, s, a, daughter)
      IMPLICIT NONE
      REAL :: pi = 3.141592653589793
      INTEGER :: nx, ny, ix, iy
      REAL :: k0, s, a 
      REAL :: norm, expntx, expnty, expnt
      REAL :: kx(nx), ky(ny), kx0, ky0
      COMPLEX :: daughter(nx,ny)

!
! Morlet wavelet in Fourier domain
!
      kx0=k0*cos(a)
      ky0=k0*sin(a)
      norm=s/SQRT(pi)
      DO iy=1,ny
      DO ix=1,nx
        expntx = -0.5*((s*kx(ix) - kx0)**2)
        expnty = -0.5*((s*ky(iy) - ky0)**2)
        daughter(ix,iy) = CMPLX(norm*EXP(expntx)*EXP(expnty))
        expnt = -0.5*((kx(ix)**2 + ky(iy)**2) + (kx0**2 + ky0**2))
        daughter(ix,iy) = daughter(ix,iy) - CMPLX(norm*EXP(expnt))
      END DO
      END DO

      RETURN
      END SUBROUTINE wave_function


!----------------------------------------------------------
!     Subroutine jackknife_1d_admit
!
!     Calculates Jackknife error on admittance.
!----------------------------------------------------------
      SUBROUTINE jackknife_1d_sgram(ss,nx,ny,na,ns,ess)
      IMPLICIT NONE
      INTEGER :: nx,ny,na,ns,ix,iy,ia,ia2,i_k,is
      REAL :: ss(nx,ny,na,ns),p1

      REAL :: mss,dss,ssj(na-1),ssj2(na),ess(nx,ny,ns)

      DO is=1,ns
      DO iy=1,ny
      DO ix=1,nx
      DO ia=1,na

        i_k=0
        p1=0.

        DO ia2=1,na

          IF (ia2==ia) THEN
            CYCLE
          END IF

          i_k=i_k+1

          p1=p1+ss(ix,iy,ia2,is)

          IF (p1.ne.0) THEN
            ssj(i_k)=p1
          ELSE
            ssj(i_k)=0.
          END IF

        END DO

        ssj2(ia)=SUM(ssj)/FLOAT(na-1)

      END DO

      mss=SUM(ssj2)/FLOAT(na)
      dss=0.

      DO ia=1,na
        dss=dss+(ssj2(ia)-mss)**2
      END DO

      dss=(SQRT(dss*FLOAT(na-1)/FLOAT(na)))
      IF (ISNAN(dss)) dss=0.

      ess(ix,iy,is)=dss
      END DO
      END DO
      END DO

      RETURN
      END SUBROUTINE jackknife_1d_sgram


!----------------------------------------------------------
!     Subroutine jackknife_1d_admit
!
!     Calculates Jackknife error on admittance.
!----------------------------------------------------------
      SUBROUTINE jackknife_1d_admit(sxp,ss1,nx,ny,na,ns,sad)
      IMPLICIT NONE
      INTEGER :: nx,ny,na,ns,ix,iy,ia,ia2,i_k,is
      REAL :: ss1(nx,ny,na,ns),p1
      REAL :: sxp(nx,ny,na,ns),xp

      REAL :: mad,dad,ad(na-1),ad2(na),sad(nx,ny,ns)

      DO is=1,ns
      DO iy=1,ny
      DO ix=1,nx
      DO ia=1,na

        i_k=0
        p1=0.
        xp=(0.,0.)

        DO ia2=1,na

          IF (ia2==ia) THEN
            CYCLE
          END IF

          i_k=i_k+1

          p1=p1+ss1(ix,iy,ia2,is)
          xp=xp+sxp(ix,iy,ia2,is)

          IF (p1.ne.0) THEN
            ad(i_k)=xp/p1
          ELSE
            ad(i_k)=0.
          END IF

        END DO

        ad2(ia)=SUM(ad)/FLOAT(na-1)

      END DO

      mad=SUM(ad2)/FLOAT(na)
      dad=0.

      DO ia=1,na
        dad=dad+(ad2(ia)-mad)**2
      END DO

      dad=(SQRT(dad*FLOAT(na-1)/FLOAT(na)))
      IF (ISNAN(dad)) dad=0.

      sad(ix,iy,is)=dad
      END DO
      END DO
      END DO

      RETURN
      END SUBROUTINE jackknife_1d_admit


!----------------------------------------------------------
!     Subroutine jackknife_1d_corr
!
!     !alculates Jackknife error on coherency.
!----------------------------------------------------------
      SUBROUTINE jackknife_1d_corr(sxp,ss1,ss2,nx,ny,na,ns,sco)
      IMPLICIT NONE
      INTEGER :: nx,ny,na,ns,ix,iy,ia,ia2,i_k,is
      REAL :: ss1(nx,ny,na,ns),ss2(nx,ny,na,ns),p1,p2
      REAL :: sxp(nx,ny,na,ns),xp

      REAL :: mco,dco,co(na-1),co2(na),sco(nx,ny,ns)

      DO is=1,ns
      DO iy=1,ny
      DO ix=1,nx
      DO ia=1,na

        i_k=0
        p1=0.
        p2=0.
        xp=(0.,0.)

        DO ia2=1,na

          IF (ia2==ia) THEN
            CYCLE
          END IF

          i_k=i_k+1

          p1=p1+ss1(ix,iy,ia2,is)
          p2=p2+ss2(ix,iy,ia2,is)
          xp=xp+sxp(ix,iy,ia2,is)

          IF ((p1.ne.0).and.(p2.ne.0)) THEN
            co(i_k)=xp/SQRT(p1*p2)
          ELSE
            co(i_k)=0.
          END IF

        END DO

        co2(ia)=SUM(co)/FLOAT(na-1)

      END DO

      mco=SUM(co2)/FLOAT(na)
      dco=0.

      DO ia=1,na
        dco=dco+(co2(ia)-mco)**2
      END DO

      dco=(SQRT(dco*FLOAT(na-1)/FLOAT(na)))
      IF (ISNAN(dco)) dco=0.

      sco(ix,iy,is)=dco
      END DO
      END DO
      END DO

      RETURN
      END SUBROUTINE jackknife_1d_corr


!----------------------------------------------------------
!     Subroutine jackknife_1d_corr
!
!     Calculates Jackknife error on coherence.
!----------------------------------------------------------
      SUBROUTINE jackknife_admit_coh(sxp,ss1,ss2,nx,ny,na,ns,sad,sco)
      IMPLICIT NONE
      INTEGER :: nx,ny,na,ns,ix,iy,ia,ia2,i_k,is
      REAL :: ss1(nx,ny,na,ns),ss2(nx,ny,na,ns),p1,p2
      COMPLEX :: sxp(nx,ny,na,ns),xp

      REAL :: mco,dco,co(na-1),co2(na),sco(nx,ny,ns)
      REAL :: mad,dad,ad(na-1),ad2(na),sad(nx,ny,ns)

      DO is=1,ns
      DO iy=1,ny
      DO ix=1,nx
      DO ia=1,na

        i_k=0
        p1=0.
        p2=0.
        xp=CMPLX(0.,0.)

        DO ia2=1,na

          IF (ia2==ia) THEN
            CYCLE
          END IF

          i_k=i_k+1

          p1=p1+ss1(ix,iy,ia2,is)
          p2=p2+ss2(ix,iy,ia2,is)
          xp=xp+sxp(ix,iy,ia2,is)

          IF ((p1.ne.0).and.(p2.ne.0)) THEN
            co(i_k)=REAL(xp)**2/(p1*p2)
            ad(i_k)=REAL(xp)/p1
          ELSE
            co(i_k)=0.
            ad(i_k)=0.
          END IF

        END DO

        co2(ia)=SUM(co)/FLOAT(na-1)
        ad2(ia)=SUM(ad)/FLOAT(na-1)

      END DO

      mco=SUM(co2)/FLOAT(na)
      mad=SUM(ad2)/FLOAT(na)
      dco=0.
      dad=0.

      DO ia=1,na
        dco=dco+(co2(ia)-mco)**2
        dad=dad+(ad2(ia)-mad)**2
      END DO

      dco=(SQRT(dco*FLOAT(na-1)/FLOAT(na)))
      dad=(SQRT(dad*FLOAT(na-1)/FLOAT(na)))
      IF (ISNAN(dco)) dco=0.
      IF (ISNAN(dad)) dad=0.

      sco(ix,iy,is)=dco
      sad(ix,iy,is)=dad
      END DO
      END DO
      END DO

      RETURN
      END SUBROUTINE jackknife_admit_coh


!----------------------------------------------------------
!     Subroutine mean_std
!
!     !omputes mean and standard deviation from signal
!----------------------------------------------------------
      SUBROUTINE mean_std(hwl,atot,mhwl,shwl)
      IMPLICIT NONE
      INTEGER :: atot,a
      REAL :: hwl(atot),mhwl,shwl

      mhwl=0.
      shwl=0.
      DO a=1,atot
        mhwl=mhwl+hwl(a)
      END DO
      mhwl=mhwl/REAL(atot)
      DO a=1,atot
        shwl=shwl+(hwl(a)-mhwl)**2
      END DO
      shwl=SQRT(shwl/REAL(atot))
      RETURN
      END SUBROUTINE mean_std


!---------------------------------------------
!     Subroutine taper_data
!
!     Puts half taper function on four sides
!     of matrix
!---------------------------------------------
      SUBROUTINE taper_data(dat,maxx,maxy,maxt)
      IMPLICIT NONE
      REAL :: pi = 3.141592653589793
      INTEGER :: maxt,it,ix,iy,maxx,maxy
      REAL :: taperx(maxx/2),tapery(maxy/2)
      REAL :: dat(maxx,maxy)

      taperx=0.
      tapery=0.
      DO it=1,maxt
        tapery(it)=0.5*(1.-COS(pi*(it-1)/maxt))
      END DO
      DO it=1,maxt
        taperx(it)=0.5*(1.-COS(pi*(it-1)/maxt))
      END DO
      DO iy=1,maxy
         DO it=1,maxt
            dat(it,iy)=taperx(it)*dat(it,iy)
            dat(maxx+1-it,iy)=taperx(it)*dat(maxx+1-it,iy)
         END DO
      END DO
      DO ix=1,maxx
         DO it=1,maxt
            dat(ix,it)=tapery(it)*dat(ix,it)
            dat(ix,maxy+1-it)=tapery(it)*dat(ix,maxy+1-it)
         END DO
      END DO
      RETURN
      END SUBROUTINE taper_data


!---------------------------------------------
!     Subroutine mirror_data
!
!     Mirrors 2D data along edges
!---------------------------------------------
      SUBROUTINE mirror_data(dat,nx,ny,nnx,nny)
      IMPLICIT NONE
      INTEGER :: nx,ny,nnx,nny,ix,iy
      REAL :: dat(nnx,nny)

      DO iy=1,ny
      DO ix=1,nx
        dat(ix+nx,iy+ny)=dat(nx-ix+1,ny-iy+1)
        dat(ix+nx,iy)=dat(nx-ix+1,iy)
        dat(ix,iy+ny)=dat(ix,ny-iy+1)
      END DO
      END DO

      RETURN
      END SUBROUTINE mirror_data


!-------------------------------------------------
!     Subroutine cfftswap
!
!     Swaps both directions of a matrix
!-------------------------------------------------
      SUBROUTINE cfftswap(spec,nx,ny)
      IMPLICIT NONE
      COMPLEX :: spec(nx,ny),tmp(nx,ny)
      INTEGER :: nx,ny,ix,iy

      tmp=spec
      DO iy=1,ny
      DO ix=1,nx
        spec(nx-ix+1,ny-iy+1)=tmp(ix,iy)
      END DO
      END DO
      RETURN
      END SUBROUTINE cfftswap


!-------------------------------------------------
!     Subroutine fftswap
!
!     Swaps both directions of a matrix
!-------------------------------------------------
      SUBROUTINE fftswap(spec,nx,ny)
      IMPLICIT NONE
      REAL :: spec(nx,ny),tmp(nx,ny)
      INTEGER :: nx,ny,ix,iy

      tmp=spec
      DO iy=1,ny
      DO ix=1,nx
        spec(nx-ix+1,ny-iy+1)=tmp(ix,iy)
      END DO
      END DO
      RETURN
      END SUBROUTINE fftswap


!-------------------------------------------------
!     Subroutine cfftshift
!
!     Swaps both directions of a matrix
!-------------------------------------------------
      SUBROUTINE cfftshift(spec,nx,ny)
      IMPLICIT NONE
      COMPLEX :: spec(nx,ny),tmp(nx,ny)
      INTEGER :: nx,ny,ix,iy

      tmp=spec
      DO iy=1,ny/2
      DO ix=1,nx/2
        spec(ix,iy)=tmp(nx/2+ix,ny/2+iy)
        spec(nx/2+ix,iy)=tmp(ix,ny/2+iy)
        spec(nx/2+ix,ny/2+iy)=tmp(ix,iy)
        spec(ix,ny/2+iy)=tmp(nx/2+ix,iy)
      END DO
      END DO
      RETURN
      END SUBROUTINE cfftshift


!-------------------------------------------------
!     Subroutine fourn
!
!     !alculates Fast Fourier Transform
!     From Numerical Recipes textbook
!-------------------------------------------------

      SUBROUTINE fourn(data,nn,ndim,isign)
      INTEGER isign,ndim,nn(ndim)
      REAL data(*)
      INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,&
          k2,n,nprev,nrem,ntot
      REAL tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      ntot=1
      do idim=1,ndim
        ntot=ntot*nn(idim)
      end do
      nprev=1
      do idim=1,ndim
        n=nn(idim)
        nrem=ntot/(n*nprev)
        ip1=2*nprev
        ip2=ip1*n
        ip3=ip2*nrem
        i2rev=1
        do i2=1,ip2,ip1
          if(i2.lt.i2rev)then
            do i1=i2,i2+ip1-2,2
              do i3=i1,ip3,ip2
                i3rev=i2rev+i3-i2
                tempr=data(i3)
                tempi=data(i3+1)
                data(i3)=data(i3rev)
                data(i3+1)=data(i3rev+1)
                data(i3rev)=tempr
                data(i3rev+1)=tempi
              end do
            end do
          endif
          ibit=ip2/2
1         if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
            i2rev=i2rev-ibit
            ibit=ibit/2
          goto 1
          endif
          i2rev=i2rev+ibit
        end do
        ifp1=ip1
2       if(ifp1.lt.ip2)then
          ifp2=2*ifp1
          theta=isign*6.28318530717959d0/(ifp2/ip1)
          wpr=-2.d0*sin(0.5d0*theta)**2
          wpi=sin(theta)
          wr=1.d0
          wi=0.d0
          do i3=1,ifp1,ip1
            do i1=i3,i3+ip1-2,2
              do i2=i1,ip3,ifp2
                k1=i2
                k2=k1+ifp1
                tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                data(k2)=data(k1)-tempr
                data(k2+1)=data(k1+1)-tempi
                data(k1)=data(k1)+tempr
                data(k1+1)=data(k1+1)+tempi
              end do
            end do
            wtemp=wr
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wtemp*wpi+wi
          end do
          ifp1=ifp2
        goto 2
        endif
        nprev=n*nprev
      end do
      return
      END SUBROUTINE fourn




<|start_filename|>geoist/magmod/include/time_conversion.h<|end_filename|>
/**
 * @file time_conversion.h
 * @author <NAME> <<EMAIL>>
 * @brief Time Conversion Subroutines
 *
 * Time conversion subroutines.
 *
 * Copyright (C) 2018 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
#ifndef TIME_CONVERSION_H
#define TIME_CONVERSION_H

#include <math.h>

/**
 * @brief Convert MJD2000 to year, month, day and day fraction
 *
 * Convert MJD2000 to year, month, day and number of decimal hours
 *
 * ref: https://en.wikipedia.org/wiki/Julian_day#Julian_or_Gregorian_calendar_from_Julian_day_number
 * Gregorian date formula applied since 1582-10-15
 * Julian date formula applied until 1582-10-04
 */

static void mjd2k_to_date(int *year, int *month, int *day, double *hours, double mjd2k) {
    const int day2k = (int)floor(mjd2k);
    const int d__ = day2k + 2451545;
    const int f0_ = d__ + 1401 ;
    const int f__ = f0_ + (d__ > 2299160 ? (((4*d__ + 274277)/146097)*3)/4 - 38: 0);
    const int e__ = 4*f__ + 3;
    const int h__ = 5*((e__ % 1461)/4) + 2;
    *day = (h__%153)/5 + 1;
    *month = (h__/153 + 2)%12 + 1;
    *year = e__/1461 - 4716 + (14 - *month)/12;
    *hours = 24.0 * (mjd2k - day2k);
}


/**
 * @brief Convert MJD2000 to a whole year.
 *
 * ref: https://en.wikipedia.org/wiki/Julian_day#Julian_or_Gregorian_calendar_from_Julian_day_number
 * Gregorian date formula applied since 1582-10-15
 * Julian date formula applied until 1582-10-04
 */

static int mjd2k_to_year(double mjd2k) {
    const int day2k = (int)floor(mjd2k);
    const int d__ = day2k + 2451545;
    const int f0_ = d__ + 1401 ;
    const int f__ = f0_ + (d__ > 2299160 ? (((4*d__ + 274277)/146097)*3)/4 - 38: 0);
    const int e__ = 4*f__ + 3;
    const int h__ = 5*((e__ % 1461)/4) + 2;
    const int month = (h__/153 + 2)%12 + 1;
    return e__/1461 - 4716 + (14 - month)/12;
}


/**
 * @brief Get MJD2000 of the year start.
 *
 * The function returns number of days between the start of the given year
 * and start of the year 2000.
 *
 * ref: https://en.wikipedia.org/wiki/Julian_day#Converting_Julian_or_Gregorian_calendar_date_to_Julian_day_number
 * ref: http://www.cs.utsa.edu/~cs1063/projects/Spring2011/Project1/jdn-explanation.html
 * Gregorian date formula applied since 1583
 * Julian date formula applied until 1582
 */

static int mjd2k_year_start(int year) {
    const int y__ = year + 4799;
    return (365*y__ + y__/4 - 2483321 + (year > 1582)*(y__/400 - y__/100 + 38));
}


/**
 * @brief Return 1 for a leap year or 0 otherwise.
 *
 * Gregorian date formula applied since 1583
 * Julian date formula applied until 1582
 */
static int is_leap_year(int year) {
    return ((year%4 == 0)&&((year%100 != 0)||(year%400 == 0)||(year <= 1582)));
}

/**
 * @brief Get number of days in the given year.
 *
 * This function returns 356 for a leap year or 366 otherwise.
 *
 * Gregorian date formula applied since 1583
 * Julian date formula applied until 1582
 */
static int days_per_year(int year) {
    return 365 + is_leap_year(year);
}

/**
 * @brief Get fraction of a year passed in one day.
 *
 * This function returns 1/366 for a leap year or 1/355 otherwise.
 *
 * Gregorian date formula applied since 1583
 * Julian date formula applied until 1582
 */
static double years_per_day(int year) {
    return is_leap_year(year) ? 1./366. : 1./365.;
}

/**
 * @brief Convert MJD2000 to decimal year.
 *
 * Gregorian date formula applied since 1582-10-15
 * Julian date formula applied until 1582-10-04
 */

static double mjd2k_to_decimal_year(double mjd2k) {
    const int year = mjd2k_to_year(mjd2k);
    return year + (mjd2k - mjd2k_year_start(year))*years_per_day(year);
}

/**
 * @brief Convert MJD2000 to year fraction.
 *
 * Equivalent to:
 *  decimal_year = mjd2k_to_decimal_year(mjd2k);
 *  year_fraction = decimal_year - floor(decimal_year);
 */

static double mjd2k_to_year_fraction(double mjd2k) {
    const int year = mjd2k_to_year(mjd2k);
    return (mjd2k - mjd2k_year_start(year))*years_per_day(year);
}

/**
 * @brief Convert decimal year to MJD2000.
 *
 * Gregorian date formula applied since 1582-10-15
 * Julian date formula applied until 1582-10-04
 */

static double decimal_year_to_mjd2k(double decimal_year) {
    if (isfinite(decimal_year)) {
        const double whole_year = floor(decimal_year);
        const int year = (int)whole_year;
        return mjd2k_year_start(year) + (decimal_year - whole_year)*days_per_year(year);
    } else {
        return decimal_year;
    }
}

#endif /* TIME_CONVERSION */


<|start_filename|>geoist/magmod/pymm_lonsincos.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings
 * - longitude sin/cos spherical terms' evaluation
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
 */

#ifndef PYMM_LONSINCOS_H
#define PYMM_LONSINCOS_H

#include <math.h>
#include "shc.h"
#include "pymm_aux.h"

/* python function definition */

#define DOC_LONSINCOS "\n"\
"  lonsin, loncos = lonsincos(longitude, degree, fast_alg=True)\n"\
"\n"\
"     For given 'longitude' (geocentric spherical) and 'degree', evaluate\n"\
"     the following sine and cosine series: \n"\
"        cos(i*longitude) for i in range(0, degree+1)\n"\
"        sin(i*longitude) for i in range(0, degree+1)\n"\
"     The longitude has to be entered in dg..\n"\
"     The 'fast_alg' boolean options forces the subroutine to use a faster\n"\
"     but slightly less precise evaluation algorithm.\n"

static PyObject* lonsincos(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"longitude", "degree", "fast_alg", NULL};

    int degree;
    int fast_alg = 1;
    double lon_dg;
    PyObject *arr_lonsin = NULL; // lonsin array
    PyObject *arr_loncos = NULL; // loncos array
    PyObject *retval = NULL; // output tuple

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "di|i:lonsincos", keywords, &lon_dg, &degree, &fast_alg
    ))
        goto exit;

    if (degree < 0)
    {
        PyErr_Format(PyExc_ValueError, "%s < 0", keywords[1]);
        goto exit;
    }

    // create the output arrays
    if (NULL == (arr_lonsin = _get_new_double_array(1, NULL, degree+1)))
        goto exit;

    if (NULL == (arr_loncos = _get_new_double_array(1, NULL, degree+1)))
        goto exit;

    if (NULL == (retval = Py_BuildValue("NN", arr_lonsin, arr_loncos)))
        goto exit;

    // evaluate series
    if (fast_alg)
        shc_azmsincos(PyArray_DATA(arr_lonsin), PyArray_DATA(arr_loncos), degree, DG2RAD*lon_dg);
    else
        shc_azmsincos_ref(PyArray_DATA(arr_lonsin), PyArray_DATA(arr_loncos), degree, DG2RAD*lon_dg);

  exit:

    // decrease reference counters to the arrays
    if (!retval && arr_lonsin){Py_DECREF(arr_lonsin);}
    if (!retval && arr_loncos){Py_DECREF(arr_loncos);}

    return retval;
}

#endif  /* PYMM_LONSINCOS_H */


<|start_filename|>geoist/magmod/pymm_sheval.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - spherical harmonic model evaluation
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_SHEVAL_H
#define PYMM_SHEVAL_H

#include <math.h>
#include <string.h>
#include <stdlib.h>

#include "math_aux.h"
#include "geo_conv.h"
#include "shc.h"
#include "pymm_aux.h"
#include "pymm_coord.h"
#include "pymm_cconv.h"

#ifndef NAN
#define NAN (0.0/0.0)
#endif

#ifndef M_PI
#define M_PI 3.14159265358979323846  /* pi */
#endif

#ifndef M_1_PI
#define M_1_PI 0.31830988618379067154  /* 1/pi */
#endif
/* Earth radius in km */
#ifndef RADIUS
#define RADIUS  6371.2
#endif

#ifndef STR
#define SRINGIFY(x) #x
#define STR(x) SRINGIFY(x)
#endif

/* mode of evaluation enumerate */
typedef enum {
    SM_INVALID = 0x0,
    SM_POTENTIAL = 0x1,
    SM_GRADIENT = 0x2,
    SM_POTENTIAL_AND_GRADIENT = 0x3,
} SHEVAL_MODE;

/* mode of evaluation check */
static SHEVAL_MODE _check_sheval_mode(int mode, const char *label)
{
    switch (mode)
    {
        case SM_POTENTIAL:
            return SM_POTENTIAL;
        case SM_GRADIENT:
            return SM_GRADIENT;
        case SM_POTENTIAL_AND_GRADIENT:
            return SM_POTENTIAL_AND_GRADIENT;
        default:
            PyErr_Format(PyExc_ValueError, "Invalid mode '%s'!", label);
            return SM_INVALID;
    }
}

/* magnetic model - auxiliary structure */
typedef struct {
    int is_internal;
    int degree;
    int nterm;
    int coord_in;
    int coord_out;
    double elps_a;
    double elps_eps2;
    double crad_ref;
    double clat_last;
    double clon_last;
    double crad_last;
    double scale_potential;
    double scale_gradient0;
    double scale_gradient1;
    double scale_gradient2;
    const double *cg;
    const double *ch;
    double *lp;
    double *ldp;
    double *lsin;
    double *lcos;
    double *rrp;
    double *psqrt;
} MODEL;

/* model structure destruction */
static void _sheval_model_destroy(MODEL *model)
{
    if(NULL != model->lp)
        free(model->lp);
    if(NULL != model->ldp)
        free(model->ldp);
    if(NULL != model->lsin)
        free(model->lsin);
    if(NULL != model->lcos)
        free(model->lcos);
    if(NULL != model->rrp)
        free(model->rrp);
    if(NULL != model->psqrt)
        free(model->psqrt);

    memset(model, 0, sizeof(MODEL));
}

/* model structure initialization */
static int _sheval_model_init(MODEL *model, int is_internal, int degree,
    int coord_in, int coord_out, const double *cg, const double *ch,
    const double scale_potential, const double *scale_gradient
    )
{
    const int nterm = ((degree+1)*(degree+2))/2;

    memset(model, 0, sizeof(MODEL));

    model->is_internal = is_internal;
    model->degree = degree;
    model->nterm = nterm;
    model->coord_in = coord_in;
    model->coord_out = coord_out;
    model->elps_a = WGS84_A;
    model->elps_eps2 = WGS84_EPS2;
    model->crad_ref = RADIUS;
    model->clat_last = NAN;
    model->clon_last = NAN;
    model->crad_last = NAN;
    model->scale_potential = scale_potential;
    model->scale_gradient0 = scale_gradient[0];
    model->scale_gradient1 = scale_gradient[1];
    model->scale_gradient2 = scale_gradient[2];
    model->cg = cg;
    model->ch = ch;

    if (NULL == (model->lp = (double*)malloc(nterm*sizeof(double))))
        goto error;

    if (NULL == (model->ldp = (double*)malloc(nterm*sizeof(double))))
        goto error;

    if (NULL == (model->lsin = (double*)malloc((degree+1)*sizeof(double))))
        goto error;

    if (NULL == (model->lcos = (double*)malloc((degree+1)*sizeof(double))))
        goto error;

    if (NULL == (model->rrp = (double*)malloc((degree+1)*sizeof(double))))
        goto error;

    if (NULL == (model->psqrt = shc_presqrt(degree)))
        goto error;

    return 0;

  error:
    _sheval_model_destroy(model);
    return 1;
}

/* single point evaluation */
static void _sheval_model_eval(MODEL *model, int mode, double *fpot,
    double *fx, double *fy, double *fz, double x, double y, double z)
{
    double glat, glon, ghgt;
    double clat, clon, crad;
    double flat, flon, frad;
    double tmp;

    int glat_needed = (model->coord_out == CT_GEODETIC_ABOVE_WGS84);

    void (*shc_relradpow)(double *rrp, int degree, double relrad) = (
        model->is_internal ? shc_relradpow_internal : shc_relradpow_external
    );

    // convert the input coordinates
    switch(model->coord_in)
    {
        case CT_GEODETIC_ABOVE_WGS84:
            glat = x; glon = y; ghgt = z;
            geodetic2geocentric_sph(&crad, &clat, &clon, glat, glon, ghgt, model->elps_a, model->elps_eps2);
            break;
        case CT_GEOCENTRIC_SPHERICAL:
            clat = DG2RAD*x; clon = DG2RAD*y; crad = z;
            if (glat_needed) {
                geocentric_sph2geodetic(&glat, &glon, &ghgt, crad, clat, clon, model->elps_a, model->elps_eps2);
            }
            break;
        case CT_GEOCENTRIC_CARTESIAN:
            cart2sph(&crad, &clat, &clon, x, y, z);
            if (glat_needed) {
                geocentric_sph2geodetic(&glat, &glon, &ghgt, crad, clat, clon, model->elps_a, model->elps_eps2);
            }
            break;
        default:
            return;
    }

    // associative Legendre functions
    if (model->clat_last != clat)
        shc_legendre(model->lp, model->ldp, model->degree, clat, model->psqrt);

    // relative radial powers
    if (model->clon_last != clon)
        shc_azmsincos(model->lsin, model->lcos, model->degree, clon);

    // sines/cosines of the longitude
    if (model->crad_last != crad)
        shc_relradpow(model->rrp, model->degree, crad/model->crad_ref);

    // save the last evaluated coordinate
    model->clat_last = clat;
    model->clon_last = clon;
    model->crad_last = crad;

    // evaluate the model
    shc_eval(
        fpot, &flat, &flon, &frad, model->degree, mode,
        clat, crad, model->cg, model->ch, model->lp, model->ldp,
        model->rrp, model->lsin, model->lcos, model->is_internal
    );

    // project the produced vectors to the desired coordinate system
    if (mode&0x2)
    {
        switch(model->coord_out)
        {
            case CT_GEOCENTRIC_SPHERICAL:
                *fx = flat;
                *fy = flon;
                *fz = frad;
                break;

            case CT_GEODETIC_ABOVE_WGS84:
                tmp = clat - DG2RAD*glat;
                rot2d(fz, fx, frad, flat, sin(tmp), cos(tmp));
                *fy = flon;
                break;

            case CT_GEOCENTRIC_CARTESIAN:
                rot2d(&tmp, fz, frad, flat, sin(clat), cos(clat));
                rot2d(fx, fy, tmp, flon, sin(clon), cos(clon));
                break;
        }

        // final scaling
        *fx *= model->scale_gradient0;
        *fy *= model->scale_gradient1;
        *fz *= model->scale_gradient2;
    }

    if (mode&0x1)
    {
        // final scaling
        *fpot *= model->scale_potential;
    }
}

/* recursive batch model_evaluation for the input numpy arrays of coordinates */
#define S(a) ((double*)(a).data)
#define P(a,i) ((double*)((char*)(a).data+(i)*(a).stride[0])) //((double*)((a).data+(i)*(a).stride[0]))
#define V(a,i) (*P(a,i))

static void _sheval1(ARRAY_DATA arrd_in, ARRAY_DATA arrd_pot, MODEL *model)
{
    if (arrd_in.ndim > 1)
    {
        npy_intp i;
        for(i = 0; i < arrd_in.dim[0]; ++i)
            _sheval1(_get_arrd_item(&arrd_in, i), _get_arrd_item(&arrd_pot, i), model);
        return;
    }

    _sheval_model_eval(model, SM_POTENTIAL, S(arrd_pot), NULL, NULL, NULL,
                    V(arrd_in, 0), V(arrd_in, 1), V(arrd_in, 2));
}

static void _sheval2(ARRAY_DATA arrd_in, ARRAY_DATA arrd_grd, MODEL *model)
{
    if (arrd_in.ndim > 1)
    {
        npy_intp i;
        for(i = 0; i < arrd_in.dim[0]; ++i)
            _sheval2(_get_arrd_item(&arrd_in, i), _get_arrd_item(&arrd_grd, i), model);
        return;
    }

    _sheval_model_eval(model, SM_GRADIENT, NULL, P(arrd_grd, 0), P(arrd_grd, 1),
        P(arrd_grd, 2), V(arrd_in, 0), V(arrd_in, 1), V(arrd_in, 2));

}

static void _sheval3(ARRAY_DATA arrd_in, ARRAY_DATA arrd_pot, ARRAY_DATA arrd_grd, MODEL *model)
{
    if (arrd_in.ndim > 1)
    {
        npy_intp i;
        for(i = 0; i < arrd_in.dim[0]; ++i)
            _sheval3(_get_arrd_item(&arrd_in, i), _get_arrd_item(&arrd_pot, i),
                    _get_arrd_item(&arrd_grd, i), model);
        return;
    }

    _sheval_model_eval(model, SM_POTENTIAL_AND_GRADIENT, S(arrd_pot),
        P(arrd_grd, 0), P(arrd_grd, 1), P(arrd_grd, 2), V(arrd_in, 0),
        V(arrd_in, 1), V(arrd_in, 2));
}

#undef V
#undef P
#undef S

/* python function definition */

#define DOC_SHEVAL "\n"\
"   arr_out = sheval(arr_in, degree, coef_g, coef_h, coord_type_in=GEODETIC_ABOVE_WGS84,\n"\
"                    coord_type_out=GEODETIC_ABOVE_WGS84, mode=GRADIENT, is_internal=True,\n"\
"                    scale_potential=1.0, scale_gradient=1.0)\n"\
"\n"\
"     Parameters:\n"\
"       arr_in - array of 3D coordinates (up to 16 dimensions).\n"\
"       degree - degree of the spherical harmonic model.\n"\
"       coef_g - vector of spherical harmonic model coefficients.\n"\
"       coef_h - vector of spherical harmonic model coefficients.\n"\
"       coord_type_in - type of the input coordinates.\n"\
"       mode - quantity to be evaluated:\n"\
"                  POTENTIAL\n"\
"                  GRADIENT (default)\n"\
"                  POTENTIAL_AND_GRADIENT\n"\
"       is_internal - boolean flag set to True by default. When set to False \n"\
"                     external field evaluation is used.\n"\
"       scale_potential - scalar value multiplied with the result potentials.\n"\
"       scale_gradient - scalar or 3 element array multiplied with the result\n"\
"                        gradient components.\n"\
"\n"


static PyObject* sheval(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {
        "arr_in", "degree", "coef_g", "coef_h", "coord_type_in",
        "coord_type_out", "mode", "is_internal",
        "scale_potential", "scale_gradient", NULL
    };
    int ct_in = CT_GEODETIC_ABOVE_WGS84;
    int ct_out = CT_GEODETIC_ABOVE_WGS84;
    int nterm, degree = 0, mode = 0x2, is_internal;
    double scale_potential = 1.0;
    double scale_gradient[3] = {1.0, 1.0, 1.0};
    PyObject *obj_is_internal = NULL; // boolean flag
    PyObject *obj_in = NULL; // input object
    PyObject *obj_cg = NULL; // coef_g object
    PyObject *obj_ch = NULL; // coef_h object
    PyObject *obj_scale = NULL; // gradient scale object
    PyObject *arr_in = NULL; // input array
    PyObject *arr_cg = NULL; // coef_g array
    PyObject *arr_ch = NULL; // coef_h array
    PyObject *arr_pot = NULL; // output array
    PyObject *arr_grd = NULL; // output array
    PyObject *arr_scale = NULL; // gradient scale array
    PyObject *retval = NULL;
    MODEL model;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OiOO|iiiOdO:sheval", keywords,
        &obj_in, &degree, &obj_cg, &obj_ch, &ct_in, &ct_out, &mode,
        &obj_is_internal, &scale_potential, &obj_scale
    ))
        goto exit;

    is_internal = (obj_is_internal == NULL) || PyObject_IsTrue(obj_is_internal);

    // check the type of the coordinate transformation
    if (CT_INVALID == _check_coord_type(ct_in, keywords[4]))
        goto exit;

    if (CT_INVALID == _check_coord_type(ct_out, keywords[5]))
        goto exit;

    // check the operation mode
    if (SM_INVALID == _check_sheval_mode(mode, keywords[6]))
        goto exit;

    // cast the objects to arrays
    if (NULL == (arr_in=_get_as_double_array(obj_in, 1, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_cg=_get_as_double_array(obj_cg, 1, 1, NPY_IN_ARRAY, keywords[2])))
        goto exit;

    if (NULL == (arr_ch=_get_as_double_array(obj_ch, 1, 1, NPY_IN_ARRAY, keywords[3])))
        goto exit;

    if (degree < 0)
    {
        PyErr_Format(PyExc_ValueError, "Invalid value %d of '%s'!", degree, keywords[1]);
        goto exit;
    }

    // check dimensions of the coefficient arrays
    nterm = ((degree+1)*(degree+2))/2;

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr_in) > MAX_OUT_ARRAY_NDIM)
    {
        PyErr_Format(PyExc_ValueError, "The input array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[0],
            PyArray_NDIM(arr_in), MAX_OUT_ARRAY_NDIM);
        goto exit;
    }

    // check the last dimension (required length of the coordinates vectors)
    if (_check_array_dim_eq(arr_in, -1, 3, keywords[0]))
        goto exit;

    // check the coefficients arrays' dimensions
    if (_check_array_dim_le(arr_cg, 0, nterm, keywords[2]))
        goto exit;

    if (_check_array_dim_le(arr_ch, 0, nterm, keywords[3]))
        goto exit;

    // handle the gradient scale factors
    if (NULL != obj_scale)
    {
        if (NULL == (arr_scale=_get_as_double_array(obj_scale, 0, 1, NPY_IN_ARRAY, keywords[9])))
            goto exit;

        if (_extract_1d_double_array(scale_gradient, 3, arr_scale, keywords[9]))
            goto exit;
    }

    // create the output arrays
    if (mode & SM_POTENTIAL)
        if (NULL == (arr_pot = PyArray_EMPTY(PyArray_NDIM(arr_in)-1, PyArray_DIMS(arr_in), NPY_DOUBLE, 0)))
            goto exit;

    if (mode & SM_GRADIENT)
        if (NULL == (arr_grd = _get_new_double_array(PyArray_NDIM(arr_in), PyArray_DIMS(arr_in), 3)))
            goto exit;

    // evaluate the model

    if(_sheval_model_init(
        &model, is_internal, degree, ct_in, ct_out,
        PyArray_DATA(arr_cg), PyArray_DATA(arr_ch),
        scale_potential, scale_gradient
    ))
        goto exit;

    switch(mode)
    {
        case SM_POTENTIAL:
            _sheval1(_array_to_arrd(arr_in), _array_to_arrd(arr_pot), &model);
            retval = arr_pot;
            break;

        case SM_GRADIENT:
            _sheval2(_array_to_arrd(arr_in), _array_to_arrd(arr_grd), &model);
            retval = arr_grd;
            break;

        case SM_POTENTIAL_AND_GRADIENT:
            _sheval3(_array_to_arrd(arr_in), _array_to_arrd(arr_pot), _array_to_arrd(arr_grd), &model);
            if (NULL == (retval = Py_BuildValue("NN", arr_pot, arr_grd)))
                goto exit;
            break;
    }

    _sheval_model_destroy(&model);

  exit:

    // decrease reference counters to the arrays
    if (arr_in){Py_DECREF(arr_in);}
    if (arr_cg){Py_DECREF(arr_cg);}
    if (arr_ch){Py_DECREF(arr_ch);}
    if (!retval && arr_grd){Py_DECREF(arr_grd);}
    if (!retval && arr_pot){Py_DECREF(arr_pot);}

    return retval;
}

#endif  /* PYMM_SHEVAL_H */


<|start_filename|>extensions/MagMod_chenshi/default.json<|end_filename|>
﻿{"pluginname":"MagMod","copyright":1,"author":"chenshi","Items":[{"name":"地磁参考场","desc":"地磁参考场，度分秒格式的经纬度需换算成度。","icon":"icon.png","language":"py","ui":"index_MagMod.html","pnodefieldname":"inputTable.txt","script":"main_MagMod.py","fieldfile":"outputFields.txt","hasviewer":false,"issourcenode":false,"isbase64":false}],"appendfiles":[],"dnpfile":""}

<|start_filename|>extensions/TsADF_zhangbei/tmpData/default.json<|end_filename|>
﻿{
	"pluginname": "TsADF",
	"copyright": 1,
	"author": "zhangbei",
	"Items": [
		{
			"name": "ADF检验",
			"desc": "Augumented Dickey-Fuller检验",
			"icon": "icon_ADF.png",
			"language": "py",
			"ui": "index_TsADF.html",
			"pnodefieldname": "inputTable.txt",
			"script": "main_TsADF.py",
			"fieldfile": "outputFields.txt",
			"hasviewer": true,
			"issourcenode": false,
			"isbase64": false
		}
	],
	"appendfiles": [],
	"dnpfile": ""
}

<|start_filename|>geoist/magmod/pytimeconv_mjd2000_to_year_fraction.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * MJD2000 to year fraction conversion
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2018 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYTIMECONV_MJD2000_TO_YEAR_FRACTION_H
#define PYTIMECONV_MJD2000_TO_YEAR_FRACTION_H

#include "pymm_aux.h"
#include "time_conversion.h"
#include <math.h>

/*---------------------------------------------------------------------------*/
/* recursive array evaluation */

static void _mjd2k_to_year_fract(ARRAY_DATA arrd_in, ARRAY_DATA arrd_out)
{
    if (arrd_in.ndim > 0)
    {
        npy_intp i;
        for(i = 0; i < arrd_in.dim[0]; ++i)
            _mjd2k_to_year_fract(
                _get_arrd_item(&arrd_in, i),
                _get_arrd_item(&arrd_out, i)
            );
    }
    else
    {
        #define S(a) (*((double*)(a).data))
        S(arrd_out) = mjd2k_to_year_fraction(S(arrd_in));
        #undef S 
    }
}

#define DOC_MJD2000_TO_YEAR_FRACTION "\n"\
"   year_fraction = mjd2000_to_year_fraction(time_mjd2k)\n\n"\
"     Output:\n"\
"       year_fraction - fraction of the year since the latest *-01-01T00:00Z\n"\
"\n"\
"     Parameters:\n"\
"       time_mjd2k - array of MJD2000 times (up to 15 dimensions).\n"\
"\n"

static PyObject* pytimeconv_mjd2000_to_year_fraction(
    PyObject *self, PyObject *args, PyObject *kwdict
)
{
    static char *keywords[] = {"time_mjd2k", NULL};
    PyObject *obj_in = NULL; // input object
    PyObject *arr_in = NULL; // input array
    PyObject *arr_out = NULL; // return array
    PyObject *retval = NULL; // return object

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "O|:mjd2000_to_year_fraction", keywords, &obj_in
    ))
        goto exit;

    // cast the objects to arrays
    if (NULL == (arr_in = _get_as_double_array(obj_in, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr_in) > (MAX_OUT_ARRAY_NDIM-1))
    {
        PyErr_Format(PyExc_ValueError, "Array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[0],
            PyArray_NDIM(arr_in), (MAX_OUT_ARRAY_NDIM-1));
        goto exit;
    }

    // create output array (adding one extra dimension)
    if (NULL == (arr_out = PyArray_EMPTY(PyArray_NDIM(arr_in), PyArray_DIMS(arr_in), NPY_DOUBLE, 0)))
        goto exit;

    // evaluate solar positions
    _mjd2k_to_year_fract(
        _array_to_arrd(arr_in),
        _array_to_arrd(arr_out)
    );

    // assign return value
    retval = arr_out;

  exit:

    // decrease reference counters to the arrays
    if (arr_in) Py_DECREF(arr_in);
    if (!retval && arr_out) Py_DECREF(arr_out);

    return retval;
}

#endif /* PYTIMECONV_MJD2000_TO_YEAR_FRACTION_H */


<|start_filename|>geoist/magmod/py_common.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 *  common shared Python definitions Python 2/3 compatibility
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2018 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PY_COMMON
#define PY_COMMON

#include <Python.h>

/* checking the Python version */

#if PY_MAJOR_VERSION == 2
  #define IS_PYTHON_2
  #if PY_MINOR_VERSION < 6
    #error "Non-supported Python minor version!"
  #endif
#elif PY_MAJOR_VERSION == 3
  #define IS_PYTHON_3
  #if PY_MINOR_VERSION < 4
    #error "Non-supported Python minor version!"
  #endif
#else
  #error "Non-supported Python major version!"
#endif

/* Python 2/3 compatibility */

#ifdef IS_PYTHON_3
#define PyString_FromString PyUnicode_FromString
#endif

/* Python dictionary operations */

static void set_dict_item_str(PyObject *dict, const char * key, PyObject *value)
{
    PyDict_SetItemString(dict, key, value);
    Py_DECREF(value);
}

static void set_dict_item_str_long(PyObject *dict, const char * key, long value)
{
    set_dict_item_str(dict, key, PyLong_FromLong(value));
}

static void set_dict_item_str_double(PyObject *dict, const char * key, double value)
{
    set_dict_item_str(dict, key, PyFloat_FromDouble(value));
}

static void set_dict_item_str_str(PyObject *dict, const char * key, const char *value)
{
    #ifdef IS_PYTHON_2
    set_dict_item_str(dict, key, PyString_FromString(value));
    #else
    set_dict_item_str(dict, key, PyUnicode_FromString(value));
    #endif
}

/* Python module initialization */

#ifdef IS_PYTHON_3
static struct PyModuleDef module_definition = {
    PyModuleDef_HEAD_INIT, NULL, NULL, -1, NULL, NULL, NULL, NULL, NULL
};
#endif

PyObject* init_python_module(
    const char *name,
    const char *doc,
    PyMethodDef *methods
)
{
    PyObject *module = NULL;
#ifdef IS_PYTHON_3
    module_definition.m_name = name;
    module_definition.m_doc = doc;
    module_definition.m_methods = methods;
    module = PyModule_Create(&module_definition);
#else /* IS_PYTHON2 */
    module = Py_InitModule3(name, methods, doc);
#endif
    return module;
}

#endif /* PY_COMMON */


<|start_filename|>geoist/magmod/pymm_vrot_cart2sph.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - vector rotation - cart2sph
 *  (i.e., vector coordinate system transformation)
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_VROT_CART2SPH_H
#define PYMM_VROT_CART2SPH_H

#include <math.h>
#include "math_aux.h"
#include "geo_conv.h"
#include "pymm_aux.h"
#include "pymm_vrot_common.h"

/* recursive vector rotation */

static void _vrot_cart2sph(ARRAY_DATA ad_i, ARRAY_DATA ad_lat,
                           ARRAY_DATA ad_lon, ARRAY_DATA ad_o)
{
    if (ad_i.ndim > 1)
    {
        npy_intp i;
        for(i = 0; i < ad_i.dim[0]; ++i)
            _vrot_cart2sph(
                _get_arrd_item(&ad_i, i),
                _get_arrd_item(&ad_lat, i),
                _get_arrd_item(&ad_lon, i),
                _get_arrd_item(&ad_o, i)
            );
    }
    else
    {
        #define SV(a) (*((double*)((a).data)))
        #define P(a,i) ((double*)((char*)(a).data+(i)*(a).stride[0])) //((double*)((a).data+(i)*(a).stride[0]))
        #define V(a,i) (*P(a,i))

        const double lat = -DG2RAD*SV(ad_lat);
        const double lon = -DG2RAD*SV(ad_lon);
        double tmp;

        // rotate around the azimuth axis
        rot2d(&tmp, P(ad_o,1), V(ad_i,0), V(ad_i,1), sin(lon), cos(lon));

        // rotate around the elevation axis
        rot2d(P(ad_o,2), P(ad_o,0), tmp, V(ad_i,2), sin(lat), cos(lat));

        #undef V
        #undef P
    }
}

/* python function definition */

#define DOC_VROT_CART2SPH "\n"\
"   arr_out = vrot_cart2sph(arr_in, arr_lat, arr_lon)\n"\
"\n"\
"     Rotate vectors from the Cartesian (XYZ) to \n"\
"     the spherical (NEC) coordinate frame for the given latitude\n"\
"     and longitude in degrees.\n"\
"\n"\
"     The inputs are:\n"\
"         arr_in - array of the input vectors\n"\
"         arr_lat - array of latitudes.\n"\
"         arr_lon - array of longitudes.\n"\
"     Scalar lat/lon values are also accepted for a single vector rotation.\n"

static PyObject* vrot_cart2sph(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"arr_in", "arr_lat", "arr_lon", NULL};
    PyObject *obj_in = NULL; // input object
    PyObject *obj_lat = NULL; // input object
    PyObject *obj_lon = NULL; // input object
    PyObject *arr_in = NULL; // input array
    PyObject *arr_lat = NULL; // input array
    PyObject *arr_lon = NULL; // input array
    PyObject *arr_out = NULL; // output array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OOO|:vrot_cart2sph", keywords,
        &obj_in, &obj_lat, &obj_lon
    ))
        goto exit;

    // cast the objects to arrays
    if (NULL == (arr_in=_get_as_double_array(obj_in, 1, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_lat=_get_as_double_array(obj_lat, 0, 0, NPY_ALIGNED, keywords[1])))
        goto exit;

    if (NULL == (arr_lon=_get_as_double_array(obj_lon, 0, 0, NPY_ALIGNED, keywords[2])))
        goto exit;

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr_in) > MAX_OUT_ARRAY_NDIM)
    {
        PyErr_Format(PyExc_ValueError, "The input array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[0],
            PyArray_NDIM(arr_in), MAX_OUT_ARRAY_NDIM);
        goto exit;
    }

    // check the dimensions
    if (_check_array_dim_eq(arr_in, -1, 3, keywords[0]))
        goto exit;

    if (_vrot_arr_check(arr_in, arr_lat, keywords[0], keywords[1]))
        goto exit;

    if (_vrot_arr_check(arr_in, arr_lon, keywords[0], keywords[2]))
        goto exit;

    // create the output array
    if (NULL == (arr_out = _get_new_double_array(PyArray_NDIM(arr_in), PyArray_DIMS(arr_in), 3)))
        goto exit;

    // rotate the vector(s)
    _vrot_cart2sph(_array_to_arrd(arr_in), _array_to_arrd(arr_lat),
                   _array_to_arrd(arr_lon), _array_to_arrd(arr_out));

    retval = arr_out;

  exit:

    // decrease reference counters of the arrays
    if (arr_in){Py_DECREF(arr_in);}
    if (arr_lat){Py_DECREF(arr_lat);}
    if (arr_lon){Py_DECREF(arr_lon);}
    if (!retval && arr_out){Py_DECREF(arr_out);}

    return retval;
}

#endif  /* PYMM_VROT_CART2SPH_H */


<|start_filename|>geoist/magmod/pyqd.c<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Magnetic Quasi Dipole Coordinates - C python bindings
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2015 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#include "common.h" /* common definitions - to be included before Python.h */
#include <Python.h>
#include <numpy/arrayobject.h>

/* module version */
#include "version.h"

/* common python utilities */
#include "py_common.h"

/* quasi dipole coordinates evaluation - old API with wrong MLT */
#include "pyqd_eval_apex.h"

/* Quasi-Dipole coordinates evaluation */
#include "pyqd_eval_qdlatlon.h"

/* Magnetic Local Time coordinates evaluation */
#include "pyqd_eval_mlt.h"

/* sub-solar point coordinates evaluation */
#include "pyqd_eval_subsol.h"

/*---------------------------------------------------------------------------*/
/* module's doc string */

#define DOC_PYQD \
"library evaluation magnetic quasi-dipole coordinates and local magnetic time."

/*---------------------------------------------------------------------------*/
/*define module's methods */
static PyMethodDef pyqd_methods[] =
{
    {"eval_qdlatlon", (PyCFunction)eval_qdlatlon, METH_VARARGS|METH_KEYWORDS, DOC_EVAL_QDLATLON},
    {"eval_mlt", (PyCFunction)eval_mlt, METH_VARARGS|METH_KEYWORDS, DOC_EVAL_MLT},
    {"eval_subsol", (PyCFunction)eval_subsol, METH_VARARGS|METH_KEYWORDS, DOC_EVAL_SUBSOL},
    {"eval_apex", (PyCFunction)eval_apex, METH_VARARGS|METH_KEYWORDS, DOC_EVAL_APEX},
    {NULL, NULL, 0, NULL} /* Sentinel - DO NOT REMOVE! */
} ;

/*---------------------------------------------------------------------------*/
/* module initialization  */

static PyObject* init_module(void)
{
    PyObject *module = init_python_module("_pyqd", DOC_PYQD, pyqd_methods);
    if (NULL == module)
        goto exit;

    PyObject *dict = PyModule_GetDict(module);
    if (NULL == dict)
        goto exit;

    /* add module specific exception */
    //PyDict_SetItemString(dict, "QDError", PyErr_NewException("_pyqd.QDError", NULL, NULL));

    /* metadata */
    set_dict_item_str_str(dict, "__author__", "<NAME> (<EMAIL>)");
    set_dict_item_str_str(dict, "__copyright__", "Copyright (C) 2015 EOX IT Services GmbH");
    set_dict_item_str_str(dict, "__licence__", "EOX licence (MIT style)");
    set_dict_item_str_str(dict, "__version__", VERSION);

  exit:
    return module;
}

/*---------------------------------------------------------------------------*/

#if PY_MAJOR_VERSION == 2

PyMODINIT_FUNC init_pyqd(void)
{
    import_array();
    init_module();
}

#else

PyObject* PyInit__pyqd(void)
{
    import_array();
    return init_module();
}

#endif


<|start_filename|>geoist/magmod/pymm_cconv.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - coordinate conversions
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_CCONV_H
#define PYMM_CCONV_H

#include "geo_conv.h"
#include "pymm_aux.h"
#include "pymm_coord.h"

#ifndef NAN
#define NAN (0.0/0.0)
#endif

#ifndef M_PI
#define M_PI 3.14159265358979323846  /* pi */
#endif

#ifndef M_1_PI
#define M_1_PI 0.31830988618379067154  /* 1/pi */
#endif

/* low-level conversion handlers */

typedef void (*f_conv)(double*, double*, double*, double, double, double);

static void conv_identity(double *x1, double *y1, double *z1,
                        double x0, double y0, double z0)
{
    *x1 = x0; *y1 = y0; *z1 = z0;
}

static void conv_WGS84_to_cartECEF(double *x, double *y, double *z,
                        double lat, double lon, double h)
{
    geodetic2geocentric_cart(x, y, z, lat, lon, h, WGS84_A, WGS84_EPS2);
}

static void conv_WGS84_to_sphECEF(double *lat1, double *lon1, double *r1,
                        double lat0, double lon0, double h0)
{
    geodetic2geocentric_sph(r1, lat1, lon1, lat0, lon0, h0, WGS84_A, WGS84_EPS2);
    *lat1 *= RAD2DG;
    *lon1 *= RAD2DG;
}

static void conv_cartECEF_to_sphECEF(double *lat, double *lon, double *r,
                         double x, double y, double z)
{
    cart2sph(r, lat, lon, x, y, z);
    *lat *= RAD2DG;
    *lon *= RAD2DG;
}

static void conv_cartECEF_to_WGS84(double *lat, double *lon, double *h,
                         double x, double y, double z)
{
    geocentric_cart2geodetic(lat, lon, h, x, y, z, WGS84_A, WGS84_EPS2);
}


static void conv_sphECEF_to_cartECEF(double *x, double *y, double *z,
                         double lat, double lon, double r)
{
    sph2cart(x, y, z, r, DG2RAD*lat, DG2RAD*lon);
}

static void conv_sphECEF_to_WGS84(double *lat1, double *lon1, double *h1,
                           double lat0, double lon0, double r0)
{
    geocentric_sph2geodetic(lat1, lon1, h1, r0, DG2RAD*lat0, DG2RAD*lon0, WGS84_A, WGS84_EPS2);
}

static f_conv _get_fconv(int ct_in, int ct_out)
{
    f_conv mtx[3][3] = {
        {conv_identity, conv_WGS84_to_sphECEF, conv_WGS84_to_cartECEF},
        {conv_sphECEF_to_WGS84, conv_identity, conv_sphECEF_to_cartECEF},
        {conv_cartECEF_to_WGS84, conv_cartECEF_to_sphECEF, conv_identity}
    };

    return mtx[ct_in][ct_out];
}

/* recursive coordinate conversion */

static void _convert(ARRAY_DATA arrd_in, ARRAY_DATA arrd_out, f_conv fconv)
{
    if (arrd_in.ndim > 1)
    {
        npy_intp i;
        for(i = 0; i < arrd_in.dim[0]; ++i)
            _convert(_get_arrd_item(&arrd_in, i), _get_arrd_item(&arrd_out, i), fconv);
    }
    else
    {
        #define P(a,i) ((double*)((char*)(a).data+(i)*(a).stride[0]))   //((double*)((a).data+(i)*(a).stride[0]))
        #define V(a,i) (*P(a,i))

        fconv(P(arrd_out, 0), P(arrd_out, 1), P(arrd_out, 2),
              V(arrd_in, 0), V(arrd_in, 1), V(arrd_in, 2));

        #undef V
        #undef P
    }
}


/* python function definition */

#define DOC_CONVERT "\n"\
"   arr_out = convert(arr_in, coord_type_in=GEODETIC_ABOVE_WGS84, coord_type_out=GEODETIC_ABOVE_WGS84)\n"\
"\n     Convert coordinates between different coordinate systems.\n"


static PyObject* convert(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"arr_in", "coord_type_in", "coord_type_out",NULL};
    int ct_in = CT_GEODETIC_ABOVE_WGS84;
    int ct_out = CT_GEODETIC_ABOVE_WGS84;
    PyObject *obj_in = NULL; // input object
    PyObject *arr_in = NULL; // input array
    PyObject *arr_out = NULL; // output array
    PyObject *retval = NULL;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "O|ii:convert", keywords, &obj_in, &ct_in, &ct_out
    ))
        goto exit;

    // check the type of the coordinate transformation
    if (CT_INVALID == _check_coord_type(ct_in, keywords[0])) goto exit;
    if (CT_INVALID == _check_coord_type(ct_out, keywords[1])) goto exit;

    // cast the object to an array
    if (NULL == (arr_in=_get_as_double_array(obj_in, 1, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr_in) > MAX_OUT_ARRAY_NDIM)
    {
        PyErr_Format(PyExc_ValueError, "The input array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[0],
            PyArray_NDIM(arr_in), MAX_OUT_ARRAY_NDIM);
        goto exit;
    }

    // check the last dimension (required length of the coordinates vectors)
    if (_check_array_dim_eq(arr_in, -1, 3, keywords[0]))
        goto exit;

    // fast-track identity
    if (ct_in == ct_out) return arr_in;

    // create the output array
    if (NULL == (arr_out = _get_new_double_array(PyArray_NDIM(arr_in), PyArray_DIMS(arr_in), 3)))
        goto exit;

    // process
    _convert(_array_to_arrd(arr_in), _array_to_arrd(arr_out), _get_fconv(ct_in, ct_out));

    //retval = Py_None; Py_INCREF(Py_None);
    retval = arr_out;

  exit:

    // decrease reference counters to the arrays
    if (arr_in){Py_DECREF(arr_in);}
    if (!retval && arr_out){Py_DECREF(arr_out);}

    return retval;
}

#endif  /* PYMM_CCONV_H */


<|start_filename|>geoist/magmod/include/geo_conv.h<|end_filename|>
/**
 * @file geo_conv.h
 * @author <NAME> <<EMAIL>>
 * @brief Geo-coordinates conversions.
 *
 * This file contains various geo-coordinates conversion utilities.
 *
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#ifndef GEO_CONV_H
#define GEO_CONV_H 1

#include <math.h>
#include "math_aux.h"

#define DG2RAD (M_PI/180.0)
#define RAD2DG (180.0*M_1_PI)

#define WGS84_A 6378.137
//#define WGS84_B 6356.752314245
#define WGS84_B 6356.7523142
#define WGS84_INVFLAT 298.257223563
#define WGS84_FLAT (1.0/WGS84_INVFLAT)
#define WGS84_EPS2 (1.0 - (WGS84_B*WGS84_B)/(WGS84_A*WGS84_A))
#define WGS84_EPS sqrt(WGS84_EPS2)
#define WGS84_RADIUS 6371.2
#ifndef M_PI
#define M_PI 3.14159265358979323846  /* pi */
#endif

#ifndef M_1_PI
#define M_1_PI 0.31830988618379067154  /* 1/pi */
#endif

/**
 * @brief Low level conversion from geodetic to geocentric coordinates.
 *
 * Required ellipsoid parameters:
 *  elp_a - semi-major axis
 *  elp_e2 - squared first eccentricity
 */

static void _geodetic2geocentric(
    double *rad_xy, double *z,
    double lat_rad, double elv,
    double elp_a, double elp_e2)
{
    const double sin_lat = sin(lat_rad);
    const double cos_lat = cos(lat_rad);

    /* radius of curvature of the ellipsoid */
    const double rc = elp_a / sqrt(1.0 - elp_e2*sin_lat*sin_lat);

    /* radius on the equator plane */
    *rad_xy = (rc + elv)*cos_lat;

    /* z coordinate */
    *z = (rc*(1.0 - elp_e2) + elv)*sin_lat;

}

/**
 * @brief Low level conversion from geocentric to geodetic geocentric coordinates.
 *
 * Uses Ferrari's solution.
 *
 * Required ellipsoid parameters:
 *  elp_a - semi-major axis
 *  elp_e2 - squared first eccentricity
 */

static void _geocentric2geodetic(
    double *lat_rad, double *elv,
    double rad_xy, double z,
    double elp_a, double elp_e2)
{
    const double p = rad_xy;
    /* Ferrari's solution */
    const double pa = p/elp_a;
    const double za = z/elp_a;
    const double pa2 = pa*pa;
    const double za2 = za*za;
    const double ee4 = elp_e2*elp_e2;
    const double rkp0 = (1.0 - elp_e2);
    const double zt = rkp0*za2;
    const double rh = (pa2 + zt - ee4)/6.0;
    const double ss = 0.25*zt*ee4*pa2;
    const double rh3 = rh*rh*rh;
    const double tmp = rh3 + ss + sqrt(ss*(ss+2.0*rh3));
    const double tt = copysign(pow(fabs(tmp), 1.0/3.0), tmp);
    const double uu = rh + tt + rh*rh/tt;
    const double vv = sqrt(uu*uu + ee4*zt);
    const double ww = 0.5*elp_e2*(uu + vv - zt)/vv;
    const double kp = 1.0 + elp_e2*(sqrt(uu + vv + ww*ww) + ww)/(uu + vv);
    const double zkp = z * kp;

    *lat_rad = atan2(zkp, p);
    *elv = norm2d(p, zkp)*(1.0/kp - rkp0)/elp_e2;
}


/**
 * @brief Convert geodetic to geocentric Cartesian coordinates.
 *
 * Covert geodetic coordinates (latitude, longitude, elevation(above ellipsoid)
 * to geocentric Cartesian (ECEF) coordinates.
 *
 * The input geodetic coordinates shall be in degrees. The units of the output
 * Cartesian coordinates are the same as the one of the ellipsoid semi-major axis.
 *
 * Required ellipsoid parameters:
 *  elp_a - semi-major axis
 *  elp_e2 - squared first eccentricity
 */

static void geodetic2geocentric_cart(
    double *x, double *y, double *z,
    double lat, double lon, double elv,
    double elp_a, double elp_e2)
{
    const double lat_rad = DG2RAD*lat;
    const double lon_rad = DG2RAD*lon;
    double rad_xy;

    _geodetic2geocentric(&rad_xy, z, lat_rad, elv, elp_a, elp_e2);

    *x = rad_xy * cos(lon_rad);
    *y = rad_xy * sin(lon_rad);
}


/**
 * @brief Convert geocentric coordinates to geocentric spherical coordinates
 *
 * Covert geodetic coordinates (latitude, longitude, elevation(above ellipsoid)
 * to geocentric spherical coordinates (radius, elevation/theta, azimuth/phi).
 *
 * The input geodetic coordinates shall be in degrees. The units of the output
 * Cartesian coordinates as well as the unit of the radial coordinate are
 * the same as the one of the ellipsoid semi-major axis. The output spherical
 * coordinates are in radians.
 *
 * Required ellipsoid parameters:
 *  elp_a - semi-major axis
 *  elp_e2 - squared first eccentricity
 */

static void geodetic2geocentric_sph(
    double *r, double *th, double *ph,
    double lat, double lon, double elv,
    double elp_a, double elp_e2)
{
    const double lat_rad = DG2RAD*lat;
    const double lon_rad = DG2RAD*lon;
    double rad_xy, z;

    _geodetic2geocentric(&rad_xy, &z, lat_rad, elv, elp_a, elp_e2);

    *r = norm2d(rad_xy, z);
    *th = atan2(z, rad_xy);
    *ph = lon_rad;
}


/**
 * @brief Convert geodetic coordinates to geocentric Cartesian and spherical ones.
 *
 * Covert geodetic coordinates (latitude, longitude, elevation(above ellipsoid))
 * to geocentric coordinates both Cartesian (x, y, z) and spherical (radius,
 * elevation/theta, azimuth/phi). The input geodetic coordinates shall be in
 * degrees. The output spherical coordinates are in radians.
 * The input geodetic coordinates shall be in degrees. The unit of the output
 * the radial coordinate is  the same as the one of the ellipsoid semi-major
 * axis. The output spherical coordinates are in radians.
 *
 * Use when both, spherical and Cartesian, coordinates are needed.
 * Note that this function is more efficient that two calls to
 *  geodetic2geocentric_cart
 *  geodetic2geocentric_sph
 *
 * Required ellipsoid parameters:
 *  elp_a - semi-major axis
 *  elp_e2 - squared first eccentricity
 */

static void geodetic2geocentric(
    double *x, double *y, double *z,
    double *r, double *th, double *ph,
    double lat, double lon, double elv,
    double elp_a, double elp_e2)
{
    const double lat_rad = DG2RAD*lat;
    const double lon_rad = DG2RAD*lon;
    double rad_xy;

    _geodetic2geocentric(&rad_xy, z, lat_rad, elv, elp_a, elp_e2);

    *x = rad_xy * cos(lon_rad);
    *y = rad_xy * sin(lon_rad);

    *r = norm2d(rad_xy, *z);
    *th = atan2(*z, rad_xy);
    *ph = lon_rad;
}

/**
 * @brief Convert spherical to Cartesian coordinates
 *
 * Covert spherical coordinates(radius, elevation/theta, azimuth/phi)
 * to Cartesian coordinates (x, y, z).
 *
 * The spherical coordinates shall be in radians.
 */

static void sph2cart(
    double *x, double *y, double *z,
    double r, double th, double ph)
{
    const double sin_th = sin(th);
    const double cos_th = cos(th);
    const double sin_ph = sin(ph);
    const double cos_ph = cos(ph);

    /* radius on the azimuth plane (z=0)*/
    const double ra = r*cos_th;

    *x = ra*cos_ph;
    *y = ra*sin_ph;
    *z = r*sin_th;
}

/**
 * @brief Convert Cartesian to spherical coordinates
 *
 * Covert Cartesian coordinates (x, y, z) to spherical
 * coordinates (radius, elevation/theta, azimuth/phi)
 *
 * The spherical coordinates are produced in radians.
 */
static void cart2sph(
    double *r, double *th, double *ph,
    double x, double y, double z)
{
    *r = norm3d(x, y, z);
    *th = asin(z/(*r));
    *ph = atan2(y, x);
}


/**
 * @brief Convert Cartesian geocentric to geodetic coordinates
 *
 * Covert Cartesian (ECEF) coordinates (x, y, z) to geodetic coordinates
 * (latitude, longitude, elevation(above ellipsoid).
 *
 * The geodetic coordinates are produced in degrees. The unit of the elevation
 * is the same as the unit of the Cartesian coordinates and of the ellipsoid
 * semi-major axis.
 */

static void geocentric_cart2geodetic(
    double *lat, double *lon, double *elv,
    double x, double y, double z,
    double elp_a, double elp_e2)
{
    double lat_rad;

    _geocentric2geodetic(&lat_rad, elv, norm2d(x, y), z, elp_a, elp_e2);

    *lat = RAD2DG*lat_rad;
    *lon = RAD2DG*atan2(y, x);
}

/**
 * @brief Convert spherical geocentric to geodetic coordinates
 *
 * Convert geocentric spherical coordinates (radius, elevation/theta, azimuth/phi)
 * to geodetic coordinates (latitude, longitude, elevation(above ellipsoid).
 *
 * The spherical coordinates are expect to be in radians.
 * The geodetic coordinates are produced in degrees. The unit of the elevation
 * is the same as the unit of the radial coordinate and of the ellipsoid
 * semi-major axis.
 */

static void geocentric_sph2geodetic(
    double *lat, double *lon, double *elv,
    double r, double th, double ph,
    double elp_a, double elp_e2)
{
    double lat_rad;

    _geocentric2geodetic(&lat_rad, elv, r*cos(th), r*sin(th), elp_a, elp_e2);

    *lat = RAD2DG*lat_rad;
    *lon = RAD2DG*ph;
}

#endif  /*GEO_CONV_H*/


<|start_filename|>geoist/magmod/pymm_sphargrd.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - spherical-harmonics
 * - evaluation of the potential gradient in the spherical coordinate system
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_SPHARGRD_H
#define PYMM_SPHARGRD_H

#include "shc.h"
#include "pymm_aux.h"
#include "pymm_cconv.h"

/* python function definition */
#define DOC_SPHARGRD "\n"\
"  v_grad = sphargrd(latitude, degree, coef_g, coef_h, leg_p, leg_dp, rrp, lonsin, loncos, is_internal=True)\n"\
"\n"\
"     Spherical harmonic evaluation of the gradient of the potential\n"\
"     (scalar) field in the geocentric spherical coordinates (latitude,\n"\
"     longitude, radius).\n"\
"     The input parameters are:\n"\
"       latitude - spherical (or geodetic) latitude in degrees at the evaluated\n"\
"                  location.\n"\
"       degree - degree of the spherical harmonic model.\n"\
"       coef_g - vector of spherical harmonic model coefficients.\n"\
"       coef_h - vector of spherical harmonic model coefficients.\n"\
"       leg_p - vector the Legendre polynomials.\n"\
"       leg_dp - vector the Legendre polynomials' derivations.\n"\
"       rrp - vector the relative radius powers.\n"\
"       lonsin - vector of the longitude cosines.\n"\
"       lonsin - vector of the longitude sines.\n"\
"       is_internal - boolean flag set to True by default. When set to False\n"\
"                     external field evaluation is used.\n"\
"\n"\

static PyObject* sphargrd(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {
        "latitude", "degree", "coef_g", "coef_h", "leg_p", "leg_dp",
        "rpp", "lonsin", "loncos", "is_internal", NULL
    };

    int degree, nterm, is_internal;
    double lat_sph;

    PyObject *obj_is_internal = NULL; // boolean flag
    PyObject *obj_cg = NULL; // coef_g object
    PyObject *obj_ch = NULL; // coef_h object
    PyObject *obj_lp = NULL; // P object
    PyObject *obj_ldp = NULL; // dP object
    PyObject *obj_rrp = NULL; // rel.rad.pow. object
    PyObject *obj_lsin = NULL; // lonsin object
    PyObject *obj_lcos = NULL; // loncos object

    PyObject *arr_out = NULL; // output array
    PyObject *arr_cg = NULL; // coef_g array
    PyObject *arr_ch = NULL; // coef_h array
    PyObject *arr_lp = NULL; // P array
    PyObject *arr_ldp = NULL; // dP array
    PyObject *arr_rrp = NULL; // rel.rad.pow. array
    PyObject *arr_lsin = NULL; // lonsin array
    PyObject *arr_lcos = NULL; // loncos array

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "diOOOOOOO|O:sphargrd", keywords,
        &lat_sph, &degree, &obj_cg, &obj_ch, &obj_lp, &obj_ldp, &obj_rrp,
        &obj_lsin, &obj_lcos, &obj_is_internal
    ))
        goto exit;

    is_internal = (obj_is_internal == NULL) || PyObject_IsTrue(obj_is_internal);

    if (degree < 0)
    {
        PyErr_Format(PyExc_ValueError, "%s < 0", keywords[1]);
        goto exit;
    }

    nterm = ((degree+1)*(degree+2))/2;

    // cast the objects to arrays
    if (NULL == (arr_cg=_get_as_double_array(obj_cg, 1, 1, NPY_IN_ARRAY, keywords[2])))
        goto exit;

    if (NULL == (arr_ch=_get_as_double_array(obj_ch, 1, 1, NPY_IN_ARRAY, keywords[3])))
        goto exit;

    if (NULL == (arr_lp=_get_as_double_array(obj_lp, 1, 1, NPY_IN_ARRAY, keywords[4])))
        goto exit;

    if (NULL == (arr_ldp=_get_as_double_array(obj_ldp, 1, 1, NPY_IN_ARRAY, keywords[5])))
        goto exit;

    if (NULL == (arr_rrp=_get_as_double_array(obj_rrp, 1, 1, NPY_IN_ARRAY, keywords[6])))
        goto exit;

    if (NULL == (arr_lsin=_get_as_double_array(obj_lsin, 1, 1, NPY_IN_ARRAY, keywords[7])))
        goto exit;

    if (NULL == (arr_lcos=_get_as_double_array(obj_lcos, 1, 1, NPY_IN_ARRAY, keywords[8])))
        goto exit;

    // check the arrays' dimensions
    if (_check_array_dim_le(arr_cg, 0, nterm, keywords[2]))
        goto exit;

    if (_check_array_dim_le(arr_ch, 0, nterm, keywords[3]))
        goto exit;

    if (_check_array_dim_le(arr_lp, 0, nterm, keywords[4]))
        goto exit;

    if (_check_array_dim_le(arr_ldp, 0, nterm, keywords[5]))
        goto exit;

    if (_check_array_dim_le(arr_rrp, 0, degree+1, keywords[6]))
        goto exit;

    if (_check_array_dim_le(arr_lsin, 0, degree+1, keywords[7]))
        goto exit;

    if (_check_array_dim_le(arr_lcos, 0, degree+1, keywords[8]))
        goto exit;

    // allocate the output array
    if (NULL == (arr_out = _get_new_double_array(1, NULL, 3)))
        goto exit;

    // the evaluation
    {
        double *out = (double*)PyArray_DATA(arr_out);
        shc_eval(
            NULL, out+0, out+1, out+2, degree, 0x2,
            DG2RAD*lat_sph, 0.0, PyArray_DATA(arr_cg),
            PyArray_DATA(arr_ch), PyArray_DATA(arr_lp),
            PyArray_DATA(arr_ldp), PyArray_DATA(arr_rrp),
            PyArray_DATA(arr_lsin), PyArray_DATA(arr_lcos),
            is_internal
        );
    }

  exit:

    // decrease reference counters to the arrays
    if (arr_cg){Py_DECREF(arr_cg);}
    if (arr_ch){Py_DECREF(arr_ch);}
    if (arr_lp){Py_DECREF(arr_lp);}
    if (arr_ldp){Py_DECREF(arr_ldp);}
    if (arr_rrp){Py_DECREF(arr_rrp);}
    if (arr_lsin){Py_DECREF(arr_lsin);}
    if (arr_lcos){Py_DECREF(arr_lcos);}

    return arr_out;
 }

#endif  /* PYMM_SPHARGRD_H */


<|start_filename|>extensions/TsEDA_chenshi/debug.json<|end_filename|>
﻿{
  "pars": {
    "method": "AN1"
  },
  "inputdata": "C:\\Users\\chens\\AppData\\Local\\Temp\\inputdata.txt",
  "OutputPath": "C:\\Users\\chens\\AppData\\Local\\Temp\\",
  "ResultFile": "C:\\Users\\chens\\AppData\\Local\\Temp\\result.json",
  "GetResult": "Print",
  "DataMode": "FileList"
}

<|start_filename|>geoist/magmod/pymm_vrot_common.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - vector rotation - common
 *  (i.e., vector coordinate system transformation)
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_VROT_COMMON_H
#define PYMM_VROT_COMMON_H

#include "pymm_aux.h"

/* array check*/

static int _vrot_arr_check(
    PyObject *arr_ref, PyObject *arr_checked,
    const char *label_ref, const char *label_checked
)
{
    //if ((PyArray_NDIM(arr_ref) > 1)||(PyArray_NDIM(arr_checked) > 0))
    if (PyArray_NDIM(arr_checked) > 0)
    {
        int d;

        if (PyArray_NDIM(arr_ref) != PyArray_NDIM(arr_checked)+1)
        {
            PyErr_Format(PyExc_ValueError, "Shape mismatch between '%s' and "
                "'%s'!", label_ref, label_checked);
            return 1;
        }

        for (d = 0; d < (PyArray_NDIM(arr_ref)-1); ++d)
        {
            if (PyArray_DIM(arr_ref, d) != PyArray_DIM(arr_checked, d))
            {
                PyErr_Format(PyExc_ValueError, "Shape mismatch between '%s' "
                    "and '%s'!", label_ref, label_checked);
                return 1;
            }
        }
    }

    return 0;
}

#endif  /* PYMM_VROT_COMMON_H */


<|start_filename|>geoist/magmod/pysunpos.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Solar position
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2017 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/
#ifndef PYSUNPOS_H
#define PYSUNPOS_H

#include "pymm_aux.h"
#include "sunpos.h"
#include <math.h>

/*---------------------------------------------------------------------------*/
/* recursive array evaluation */

#ifndef RAD2DEG
#define RAD2DEG (180.0/M_PI)
#endif

#ifndef DEG2RAD
#define DEG2RAD (M_PI/180.0)
#endif

static void _sunpos(
    ARRAY_DATA arrd_dim,
    ARRAY_DATA arrd_mjd,
    ARRAY_DATA arrd_lat,
    ARRAY_DATA arrd_lon,
    ARRAY_DATA arrd_rad,
    ARRAY_DATA arrd_dtt,
    ARRAY_DATA arrd_out
)
{
    if (arrd_dim.ndim > 0)
    {
        npy_intp i;
        for(i = 0; i < arrd_dim.dim[0]; ++i)
            _sunpos(
                _get_arrd_item(&arrd_dim, i),
                _get_arrd_item(&arrd_mjd, i),
                _get_arrd_item(&arrd_lat, i),
                _get_arrd_item(&arrd_lon, i),
                _get_arrd_item(&arrd_rad, i),
                _get_arrd_item(&arrd_dtt, i),
                _get_arrd_item(&arrd_out, i)
            );
    }
    else
    {
        #define P(a,i) ((double*)((char*)(a).data+(i)*(a).stride[0]))
        #define V(a,i) (*P(a,i))

        double rasc, decl, hang, zenith, azimuth;
        double mjd2k, dtt, lat, lon, rad, pres, temp;

        mjd2k = V(arrd_mjd,0); // decimal days since 2000-01-01T00:00:00
        dtt = V(arrd_dtt,0); // seconds (TT to UT offset)
        lat = V(arrd_lat,0) * DEG2RAD; // deg --> rad
        lon = V(arrd_lon,0) * DEG2RAD; // deg --> rad
        rad = V(arrd_rad,0); // km
        pres = 0; // atm
        temp = 0; // dgC

        sunpos5equat(&decl, &rasc, &hang, mjd2k, dtt, lon);
        sunpos5eq2hor(&azimuth, &zenith, decl, hang, lat, rad, pres, temp);

        V(arrd_out,0) = RAD2DEG * decl; // rad --> deg
        V(arrd_out,1) = RAD2DEG * rasc; // rad --> deg
        V(arrd_out,2) = RAD2DEG * hang; // rad --> deg
        V(arrd_out,3) = RAD2DEG * azimuth; // rad --> deg
        V(arrd_out,4) = RAD2DEG * zenith; // rad --> deg

        #undef V
        #undef P
    }
}

#define DOC_SUNPOS "\n"\
"   arr_out = sunpos(time_mjd2k, lat, lon, rad, dtt)\n\n"\
"     Output:\n"\
"       arr_out - array of the Sun equatorial and horizontal coordinates:\n"\
"                 - declination\n"\
"                 - right ascension\n"\
"                 - hour angle\n"\
"                 - azimuth\n"\
"                 - zenith\n"\
"                 All angles are in deg.\n\n"\
"     Parameters:\n"\
"       time_mjd2k - array of MJD2000 times (up to 15 dimensions).\n"\
"       lat - array of latitudes in deg.\n"\
"       lon - array of longitudes in deg.\n"\
"       rad - array of radii in km (0 means no parallax correction).\n"\
"       dtt - array of offsets to TT in sec. \n"

static PyObject* pysunpos_sunpos(PyObject *self, PyObject *args, PyObject *kwdict)
{
    static char *keywords[] = {"time_mjd2k", "lat", "lon", "rad", "dtt", NULL};
    PyObject *obj_mjd = NULL, *arr_mjd = NULL;
    PyObject *obj_lat = NULL, *arr_lat = NULL;
    PyObject *obj_lon = NULL, *arr_lon = NULL;
    PyObject *obj_rad = NULL, *arr_rad = NULL;
    PyObject *obj_dtt = NULL, *arr_dtt = NULL;
    PyObject *arr_out = NULL; // return object
    PyObject *retval = NULL; // return object
    PyObject *arr = NULL;
    int idx = 0;

    // parse input arguments
    if (!PyArg_ParseTupleAndKeywords(
        args, kwdict, "OOOOO|:sunpos", keywords,
        &obj_mjd, &obj_lat, &obj_lon, &obj_rad, &obj_dtt
    ))
        goto exit;

    // cast the objects to arrays
    if (NULL == (arr_mjd = _get_as_double_array(obj_mjd, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_lat = _get_as_double_array(obj_lat, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_lon = _get_as_double_array(obj_lon, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_rad = _get_as_double_array(obj_rad, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    if (NULL == (arr_dtt = _get_as_double_array(obj_dtt, 0, 0, NPY_ALIGNED, keywords[0])))
        goto exit;

    // check the array dimensions and allocate the output array
        // get array with the max.dimension
    arr = arr_mjd;
    idx = 0;
    if (PyArray_NDIM(arr) < PyArray_NDIM(arr_lat)) {idx = 1; arr = arr_lat;}
    if (PyArray_NDIM(arr) < PyArray_NDIM(arr_lon)) {idx = 2; arr = arr_lon;}
    if (PyArray_NDIM(arr) < PyArray_NDIM(arr_rad)) {idx = 3; arr = arr_rad;}
    if (PyArray_NDIM(arr) < PyArray_NDIM(arr_dtt)) {idx = 4; arr = arr_dtt;}

    // check maximum allowed input array dimension
    if (PyArray_NDIM(arr) > (MAX_OUT_ARRAY_NDIM-1))
    {
        PyErr_Format(PyExc_ValueError, "Array dimension of '%s'"\
            " %d exceeds the allowed maximum value %d!", keywords[idx],
            PyArray_NDIM(arr), (MAX_OUT_ARRAY_NDIM-1));
        goto exit;
    }

    // check array dimensions
    if (PyArray_NDIM(arr_mjd) > 0)
        if (_check_equal_shape(arr_mjd, arr, keywords[0], keywords[idx]))
            goto exit;

    if (PyArray_NDIM(arr_lat) > 0)
        if (_check_equal_shape(arr_lat, arr, keywords[1], keywords[idx]))
            goto exit;

    if (PyArray_NDIM(arr_lon) > 0)
        if (_check_equal_shape(arr_lon, arr, keywords[2], keywords[idx]))
            goto exit;

    if (PyArray_NDIM(arr_rad) > 0)
        if (_check_equal_shape(arr_rad, arr, keywords[3], keywords[idx]))
            goto exit;

    if (PyArray_NDIM(arr_dtt) > 0)
        if (_check_equal_shape(arr_dtt, arr, keywords[4], keywords[idx]))
            goto exit;

    // create output array (adding one extra dimension)
    {
        npy_intp *dims_src = PyArray_DIMS(arr);
        npy_intp dims_dst[MAX_OUT_ARRAY_NDIM];
        int i, ndim_src = PyArray_NDIM(arr);

        for (i = 0; i < ndim_src; ++i) dims_dst[i] = dims_src[i];
        dims_dst[ndim_src] = 5;

        if (NULL == (arr_out = PyArray_EMPTY(ndim_src + 1, dims_dst, NPY_DOUBLE, 0)))
            goto exit;
    }

    // evaluate solar positions
    _sunpos(
        _array_to_arrd(arr),
        _array_to_arrd(arr_mjd),
        _array_to_arrd(arr_lat),
        _array_to_arrd(arr_lon),
        _array_to_arrd(arr_rad),
        _array_to_arrd(arr_dtt),
        _array_to_arrd(arr_out)
    );

    // assign return value
    retval = arr_out;

  exit:

    // decrease reference counters to the arrays
    if (arr_mjd) Py_DECREF(arr_mjd);
    if (arr_lat) Py_DECREF(arr_lat);
    if (arr_lon) Py_DECREF(arr_lon);
    if (arr_rad) Py_DECREF(arr_rad);
    if (arr_dtt) Py_DECREF(arr_dtt);
    if (!retval && arr_out) Py_DECREF(arr_out);

    return retval;
}

#endif /* PYSUNPOS_H */


<|start_filename|>geoist/magmod/include/math_aux.h<|end_filename|>
/** 
 * @file math_aux.h
 * @author <NAME> <<EMAIL>>
 * @brief Auxiliary mathematical subroutines.
 *
 * This file contains definitions of auxiliary mathematical subroutines 
 * to be used all over the whole program.
 *
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#ifndef MATH_AUX_H
#define MATH_AUX_H 1

#include<math.h>

/**
 * @brief Evaluate quadratic norm of a 3D vector.
 *  v = (x, y, z)
 *  l = |v| = sqrt(x^2 + y^2 + z^2) 
 */
static double norm3d(double x, double y, double z)
{ 
    return sqrt(x*x + y*y + z*z); 
}

/**
 * @brief Evaluate quadratic norm of a 2D vector.
 *  v = (x, y)
 *  l = |v| = sqrt(x^2 + y^2) 
 */
static double norm2d(double x, double y)
{
    return sqrt(x*x + y*y);
}

/**
 * @brief Evaluate derivative of a quadratic norm of a 3D vector.
 *  l = |v| 
 *  dl = d(|v|)
 *  v = (x, y, z)
 *  dv = (dx, dy, dz) 
 *  if |v| >= 0:
 *    dl = (v/|v|)*dv = (x*dx + y*dy + z*dz)/l
 *  if |v| == 0: 
 *    dl = |dv| = sqrt(dx^2 + dy^2 + dz^2)
 */
static double dnorm3d(double x, double y, double z, double l,
               double dx, double dy, double dz)
{ 
    if (l > 0.0)
        return (x*dx + y*dy + z*dz)/l;
    else if (l == 0.0)
        return norm3d(dx, dy, dz);
    else 
        return NAN;
}

/**
 * @brief Evaluate derivative of a quadratic norm of a 2D vector.
 *  l = |v| 
 *  dl = d(|v|)
 *  v = (x, y)
 *  dv = (dx, dy) 
 *  if |v| >= 0:
 *    dl = (v/|v|)*dv = (x*dx + y*dy)/l
 *  if |v| == 0: 
 *    dl = |dv| = sqrt(dx^2 + dy^2)
 */
static double dnorm2d(double x, double y, double l, double dx, double dy)
{ 
    if (l > 0.0)
        return (x*dx + y*dy)/l;
    else if (l == 0.0)
        return norm2d(dx, dy);
    else 
        return NAN;
}

/**
 * @brief 2D vector rotation.
 *
 * x' = x*cos(a) - y*sin(a)
 * y' = x*sin(a) + y*cos(a)
 */
static void rot2d(double *x_out, double *y_out, double x, double y, double a_sin, double a_cos)
{
    *x_out = x*a_cos - y*a_sin;
    *y_out = x*a_sin + y*a_cos;
}

#endif  /*MATH_AUX_H*/


<|start_filename|>extensions/TsTools_zhangbei/tmpData/debug.json<|end_filename|>
﻿{
  "pars": {
    "Tools": "多项式拟合",
    "WindowSize": 5,
    "Periods": 1,
    "Factor": 0.5,
    "FitNum": "3",
    "Methods": ""
  },
  "tsdata": "C:\\Users\\chens\\AppData\\Local\\Temp\\tsdata.txt",
  "OutputPath": "C:\\Users\\chens\\AppData\\Local\\Temp\\",
  "ResultFile": "C:\\Users\\chens\\AppData\\Local\\Temp\\result.json",
  "GetResult": "Print",
  "DataMode": "Url"
}

<|start_filename|>extensions/TsADF_zhangbei/debug.json<|end_filename|>
﻿{
  "pars": {
    "windowsize": "36"
  },
  "allfields": true,
  "names": {
    "windowsize": "36"
  },
  "tsdata": "C:\\Users\\chens\\AppData\\Local\\Temp\\tsdata.txt",
  "OutputPath": "C:\\Users\\chens\\AppData\\Local\\Temp\\",
  "ResultFile": "C:\\Users\\chens\\AppData\\Local\\Temp\\result.json",
  "GetResult": "Print",
  "DataMode": "Url"
}

<|start_filename|>geoist/magmod/pymm_aux.h<|end_filename|>
/*-----------------------------------------------------------------------------
 *
 * Geomagnetic Model - C python bindings - auxiliary functions
 *
 * Author: <NAME> <<EMAIL>>
 *
 *-----------------------------------------------------------------------------
 * Copyright (C) 2014 EOX IT Services GmbH
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *-----------------------------------------------------------------------------
*/

#ifndef PYMM_AUX_H
#define PYMM_AUX_H
#ifndef M_PI
#define M_PI 3.14159265358979323846  /* pi */
#endif

#ifndef M_1_PI
#define M_1_PI 0.31830988618379067154  /* 1/pi */
#endif
/*
 * Check the input python object and convert it to a double precision NumPy
 * array ensuring the native byte-order.
 */
static PyObject* _get_as_int_array(PyObject *data, int dmin, int dmax,
                int reqs, const char *label)
{
    PyArray_Descr *dtype = PyArray_DescrFromType(NPY_INT32);
    PyObject *arr = PyArray_FromAny(data, dtype, dmin, dmax, reqs, NULL);
    /*
    if (NULL == arr)
        PyErr_Format(PyExc_ValueError, "Failed to cast %s to an array!", label);
    */
    return arr;
}

/*
 * Check the input python object and convert it to a double precision NumPy
 * array ensuring the native byte-order.
 */
static PyObject* _get_as_double_array(PyObject *data, int dmin, int dmax,
                int reqs, const char *label)
{
    PyArray_Descr *dtype = PyArray_DescrFromType(NPY_FLOAT64);
    PyObject *arr = PyArray_FromAny(data, dtype, dmin, dmax, reqs, NULL);
    /*
    if (NULL == arr)
        PyErr_Format(PyExc_ValueError, "Failed to cast %s to an array!", label);
    */
    return arr;
}

/*
 * Get new allocated NumPy array. The first (N-1) dimensions as read from
 * the array of dimensions (allowing easily set the same shape as the input
 * matrix). The last Nth dimension is overridden by the 'dim_last' value.
 */
static PyObject* _get_new_double_array(npy_intp ndim, const npy_intp *dims, npy_intp dim_last)
{
    npy_intp i;
    npy_intp dims_new[MAX_OUT_ARRAY_NDIM];

    if (ndim > MAX_OUT_ARRAY_NDIM)
    {
        PyErr_Format(PyExc_ValueError, "The output array dimension "\
            " %d exceeds the allowed maximum value %d!", (int)ndim, MAX_OUT_ARRAY_NDIM);
        return NULL;
    }

    for (i = 0; i < (ndim-1); ++i)
        dims_new[i] = dims[i];
    if (ndim >= 1)
        dims_new[ndim-1] = dim_last;

    return PyArray_EMPTY(ndim, dims_new, NPY_DOUBLE, 0);
}

/*
 * Check that the dimension of the array has the required value.
 */
static int _check_array_dim_eq(PyObject *arr, int dim, size_t size, const char *label)
{
    if (dim < 0)
        dim += PyArray_NDIM(arr);
    int rv = PyArray_DIM(arr, dim) != size;
    if (rv)
        PyErr_Format(PyExc_ValueError, "The dimension #%d of '%s'"\
            " %ld is not equal the allowed value %ld!", dim, label,
            (size_t)PyArray_DIM(arr, dim), size);
    return rv;
}

/*
 * Check that the dimension of the array is greater then or equal to the required value.
 */
static int _check_array_dim_le(PyObject *arr, int dim, size_t size, const char *label)
{
    if (dim < 0)
        dim += PyArray_NDIM(arr);
    int rv = PyArray_DIM(arr, dim) < size;
    if (rv)
        PyErr_Format(PyExc_ValueError, "The dimension #%d of '%s'"\
            " %ld is lower than the minimum allowed value %ld!", dim, label,
            (size_t)PyArray_DIM(arr, dim), size);
    return rv;
}

/*
 * Check that the array dimensions match the required values
 */

static int _check_arr_dims_all_eq(PyObject *arr, npy_intp ndim, const npy_intp *dims, const char *label)
{
    npy_intp dim;

    if (PyArray_NDIM(arr) != ndim)
    {
        PyErr_Format(PyExc_ValueError, "The number of dimensions of '%s'"\
            " %ld does not match the required value %ld!", label,
            (size_t)(PyArray_NDIM(arr)), (size_t)ndim);
        return 1;
    }

    for (dim = 0; dim < ndim; ++dim)
    {
        if (PyArray_DIM(arr, dim) != dims[dim])
        {
            PyErr_Format(PyExc_ValueError, "The dimensions #%ld of '%s'"\
            " %ld does not match the required value %ld!", (size_t)dim, label,
            (size_t)PyArray_DIM(arr, dim), (size_t)dims[dim]);
            return 1;
        }
    }

    return 0;
}

/*
 * Check that two arrays have the same shape
 */

static int _check_equal_shape(PyObject *arr0, PyObject *arr1, const char *label0, const char *label1)
{
    npy_intp dim;

    if (PyArray_NDIM(arr0) != PyArray_NDIM(arr1))
    {
        PyErr_Format(
            PyExc_ValueError, "Array dimension mismatch between '%s' and '%s'!",
            label0, label1
        );
        return 1;
    }

    for (dim = 0; dim < PyArray_NDIM(arr0); ++dim)
    {
        if (PyArray_DIM(arr0, dim) != PyArray_DIM(arr1, dim))
        {
            PyErr_Format(
                PyExc_ValueError, "Array shape mismatch between '%s' and '%s'!",
                label0, label1
            );
            return 1;
        }
    }

    return 0;
}

/*
 *  Extract 1D Numpy array or broadcast scalar to a 1D C array.
 */
static int _extract_1d_double_array(double *out, size_t size, PyObject *arr_in, const char *label)
{
    npy_intp ndim = PyArray_NDIM(arr_in);
    void *data = PyArray_DATA(arr_in); //void *data = PyArray_DATA(arr_in);

    if (ndim == 0)
    {
        size_t i;
        for (i = 0 ; i < size; ++i)
            out[i] = *((double*)data);
    }
    else if ((ndim == 1) && (PyArray_DIM(arr_in, 0) == size))
    {
        npy_intp stride = PyArray_STRIDE(arr_in, 0);
        size_t i;
        for (i = 0 ; i < size; ++i)
            out[i] = *((double*)((char*)data + i*stride)); //*((double*)(data + i*stride));
    }
    else
    {
        PyErr_Format(PyExc_ValueError, "Invalid %s dimension!", label);
        return 1;
    }

    return 0;
}
/*
 * Extraction of the lower dimensional parts of the arrays.
 */

typedef struct {
    void *data; //void *data; 
    npy_intp ndim;
    const npy_intp *dim;
    const npy_intp *stride;
} ARRAY_DATA;

static ARRAY_DATA _array_to_arrd(PyObject *arr)
{
    ARRAY_DATA arrd = {
        PyArray_DATA(arr),
        PyArray_NDIM(arr),
        PyArray_DIMS(arr),
        PyArray_STRIDES(arr)
    };
    return arrd;
}

static ARRAY_DATA _get_arrd_item(const ARRAY_DATA *arrd, npy_intp idx)
{
    if (arrd->ndim > 0) // extract sub-dimension
    {
        ARRAY_DATA arrd_sub = {
            (char *)arrd->data + idx*arrd->stride[0], //arrd->data + idx*arrd->stride[0],
            arrd->ndim - 1,
            arrd->dim + 1,
            arrd->stride + 1
        };
        return arrd_sub;
    }
    else // treat as scalar
    {
        return *arrd;
    }
}

#endif  /* PYMM_AUX_H */


<|start_filename|>extensions/EQcatlog_chenshi/tmpData/debug.json<|end_filename|>
﻿{
  "pars": {
    "Dataset": "CN"
  },
  "OutputPath": "C:\\Users\\chens\\AppData\\Local\\Temp\\",
  "ResultFile": "C:\\Users\\chens\\AppData\\Local\\Temp\\result.json",
  "GetResult": "Print",
  "DataMode": "Url"
}