text: /*-----------------------------------------------------------------------------
  ; AACGM library
  ;
  ; a collection of C routines intended to fully exploit the functionality of
  ; the AACGM coordinates [Shepherd, 2014] including use of the AACGM
  ; coefficients and field line tracing
  ;
  ; 20140611 SGS v0.0  Modification to existing AACGM C and IDL software, but
  ;                    includes additional features: linear interpolation and
  ;                    fieldline tracing.
  ; 20140702 SGS v0.1  Made operational, moved globals to library, fixed bug
  ;                    in malloc.
  ; 20140703 SGS v0.2  Added AACGM_GetDateTime() function to provide capability
  ;                    of getting date and time used in calculations (Jeff).
  ;                    Changed behavior to abort (with message) if not data/time
  ;                    is set.
  ;                    Variable fact changed to double from unsigned long which
  ;                    is not big enought on some 32-bit systems.
  ;                    Switched from NAN to HUGE_VAL for undefined result
  ; 20140918 SGS v1.0  change function names to _v2 for wider distribution
  ; 20150810 SGS v1.1  added code to default to geodetic coordinates for inverse
  ;                    transformation. This code was left out in the C version.
  ; 20170308 SGS v1.2  Added static to global variables in order to work with RST
  ;                    library.
  ;
  ; Functions:
  ;
  ; AACGM_v2_Rylm
  ; AACGM_v2_Alt2CGM  - not used
  ; AACGM_v2_CGM2Alt
  ; AACGM_v2_Sgn
  ; convert_geo_coord_v2
  ; AACGM_v2_LoadCoefFP
  ; AACGM_v2_LoadCoef
  ; AACGM_v2_LoadCoefs
  ; AACGM_v2_Convert
  ; AACGM_v2_SetDateTime
  ; AACGM_v2_GetDateTime
  ; AACGM_v2_SetNow
  ; AACGM_v2_errmsg
  ;
  
  ; AACGM_v2_Dayno
  ; AACGM_v2_Trace
  ; AACGM_v2_Trace_inv
  ;
  ;------------------------------------------------------------------------------
  */
  
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  #include <math.h>
  #include <time.h>
  #include "aacgmlib_v2.h"
  #include "igrflib.h"
  #include "genmag.h"
  
  #define DEBUG 0
  
  /* put these in the library header file when you figure out how to do so... */
  static struct {
    int year;
    int month;
    int day;
    int hour;
    int minute;
    int second;
    int dayno;
    int daysinyear;
  } aacgm_date = {-1,-1,-1,-1,-1,-1,-1,-1};
  
  static int myear = 0;       /* model year: 5-year epoch */
  static double fyear = 0.;   /* floating point year */
  
  static int myear_old = -1;
  static double fyear_old = -1.;
  
  static double height_old[2] = {-1,-1};
  
  static struct {
    double coef[AACGM_KMAX][NCOORD][POLYORD][NFLAG];      /* interpolated coefs */
    double coefs[AACGM_KMAX][NCOORD][POLYORD][NFLAG][2];  /* bracketing coefs */
  } sph_harm_model;
  
  /* SGS added for MSC compatibility */
  #ifndef complex
  struct complex {
    double x;
    double y;
  };
  #endif
  
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_Rylm
  ;
  ; PURPOSE:
  ;       Computes an array of real spherical harmonic function values
  ;       Y_lm(phi,theta) for a given colatitiude (phi) and longitude (theta)
  ;       for all the values up to l = order, which is typically 10. The
  ;       values are stored in a 1D array of dimension (order+1)^2. The
  ;       indexing scheme used is:
  ;
  ;        l    0  1  1  1  2  2  2  2  2  3  3  3  3  3  3  3  4  4  4  4  4 ...
  ;        m    0 -1  0  1 -2 -1  0  1  2 -3 -2 -1  0  1  2  3 -4 -3 -2 -1  0 ...
  ;C & IDL j    0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 ...
  ;FORTRAN j    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 ...
  ; 
  ; CALLING SEQUENCE:
  ;       AACGM_v2_Rylm, colat,lon,order, ylmval
  ;
  ;     Input Arguments:  
  ;       colat         - The colatitude of the point for which the spherical
  ;                       harmonic Y_lm is to be calculated
  ;
  ;       lon           - The longitude of the point for which the spherical
  ;                       harmonic Y_lm is to be calculated
  ;
  ;       order         - The order of the spherical harmonic function expansion.
  ;                       The total number of terms computed will be (order+1)^2
  ;
  ;     Output Argument:
  ;       ylmval        - 1D array of spherical harmonic functions at the point
  ;                       (colat,lon)
  ;
  ; HISTORY:
  ;
  ; Revision 1.1  94/10/12  15:24:21  15:24:21  baker (Kile Baker S1G)
  ; Initial revision
  ;
  ; subsequent revisions, porting to C and IDL by baker, wing and barnes.
  ;
  ; NOTES by SGS:
  ;
  ; It is likely that the original version was taken from FORTRAN and used array
  ; indexing that begins with 1. Indexing is somewhat more natural using the
  ; zeros-based indexing of C/IDL. Indices have thus been changed from the
  ; original version.
  ;
  ; It appears that the original version used unnormalized spherical harmonic
  ; functions. I suspect this might be better, but realized it too late. The
  ; coefficients I derived are for orthonormal spherical harmonic functions
  ; which then require the same for evaluation. I believe that the original
  ; authors used orthogonal spherical harmonic functions which eliminate the
  ; need for computing the normalization factors. I suspect this is just fine,
  ; but have not tested it.
  ; 
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_Rylm(double colat, double lon, int order, double *ylmval)
  {
    struct complex z1, z2;
    struct complex q_fac, q_val;
    int k, l, m, kk;
    int ia,ib,ic,id;
    double cos_theta, sin_theta;  
    double cos_lon, sin_lon;
    double l2, tl, fac;
    double ca, cb, d1;
    double *ffff;
    /* long unsigned *fact;    not big enough for 32-bit machines */
    double *fact;
  
    cos_theta = cos(colat);
    sin_theta = sin(colat);
  
    cos_lon = cos(lon);
    sin_lon = sin(lon);
  
    d1 = -sin_theta;
    z2.x = cos_lon;
    z2.y = sin_lon;
  
    z1.x = d1 * z2.x;
    z1.y = d1 * z2.y;
    q_fac = z1;
  
    /*
     * Generate Zonal Harmonics (P_l^(m=0) for l = 1,order) using recursion
     * relation (6.8.7), p. 252, Numerical Recipes in C, 2nd. ed., Press. W.
     * et al. Cambridge University Press, 1992) for case where m = 0.
     *
     * l Pl = cos(theta) (2l-1) Pl-1 - (l-1) Pl-2          (6.8.7)
     *
     * where Pl = P_l^(m=0) are the associated Legendre polynomials
     *
     */
  
    ylmval[0] = 1;            /* l = 0, m = 0 */
    ylmval[2] = cos_theta;    /* l = 1, m = 0 */
  
    for (l=2; l<=order; l++) {
      /* indices for previous two values: k = l * (l+1) + m with m=0 */
      ia = (l-2)*(l-1);
      ib = (l-1)*l;
      ic = l * (l+1);
  
      ylmval[ic] = (cos_theta * (2*l-1) * ylmval[ib] - (l-1)*ylmval[ia])/l;
    }
  
    /*
     * Generate P_l^l for l = 1 to (order+1)^2 using algorithm based upon (6.8.8)
     * in Press et al., but incorporate longitude dependence, i.e., sin/cos (phi)
     *
     * Pll = (-1)^l (2l-1)!! (sin^2(theta))^(l/2)
     *
     * where Plm = P_l^m are the associated Legendre polynomials
     *
     */
  
    q_val = q_fac;
    ylmval[3] = q_val.x;    /* l = 1, m = +1 */
    ylmval[1] = -q_val.y;   /* l = 1, m = -1 */
    for (l=2; l<=order; l++) {
      d1    = l*2 - 1.;
      z2.x  = d1 * q_fac.x;
      z2.y  = d1 * q_fac.y;
      z1.x  = z2.x * q_val.x - z2.y * q_val.y;
      z1.y  = z2.x * q_val.y + z2.y * q_val.x;
      q_val = z1;
  
      /* indices for previous two values: k = l * (l+1) + m */
      ia = l*(l+2);   /* m = +l */
      ib = l*l;       /* m = -l */
  
      ylmval[ia] =  q_val.x;
      ylmval[ib] = -q_val.y;
    }
  
    /*
     * Generate P_l,l-1 to P_(order+1)^2,l-1 using algorithm based upon (6.8.9)
     * in Press et al., but incorporate longitude dependence, i.e., sin/cos (phi)
     *
     * Pl,l-1 = cos(theta) (2l-1) Pl-1,l-1
     *
     */
  
    for (l=2; l<=order; l++) {
  
      l2 = l*l;
      tl = 2*l;
      /* indices for Pl,l-1; Pl-1,l-1; Pl,-(l-1); Pl-1,-(l-1) */
      ia = l2 - 1;
      ib = l2 - tl + 1;
      ic = l2 + tl - 1;
      id = l2 + 1;
  
      fac = tl - 1;
      ylmval[ic] = fac * cos_theta * ylmval[ia];      /* Pl,l-1   */
      ylmval[id] = fac * cos_theta * ylmval[ib];      /* Pl,-(l-1) */
    }
  
    /*
     * Generate remaining P_l+2,m to P_(order+1)^2,m for each m = 1 to order-2
     * using algorithm based upon (6.8.7) in Press et al., but incorporate
     * longitude dependence, i.e., sin/cos (phi).
     *
     * for each m value 1 to order-2 we have P_mm and P_m+1,m so we can compute
     * P_m+2,m; P_m+3,m; etc.
     *
     */
  
    for (m=1; m<=order-2; m++) {
      for (l=m+2; l<=order; l++) {
        ca = ((double) (2*l-1))/(l-m);
        cb = ((double) (l+m-1))/(l-m);
  
        l2 = l*l;
        ic = l2 + l + m;
        ib = l2 - l + m;
        ia = l2 - l - l - l + 2 + m;
        /* positive m */
        ylmval[ic] = ca * cos_theta * ylmval[ib] - cb * ylmval[ia];
  
        ic -= (m+m);
        ib -= (m+m);
        ia -= (m+m);
        /* negative m */
        ylmval[ic] = ca * cos_theta * ylmval[ib] - cb * ylmval[ia];
      }
    }
  
    /*
     * Normalization added here (SGS)
     * 
     * Note that this is NOT the standard spherical harmonic normalization factors
     * 
     * The recursive algorithms above treat positive and negative values of m in
     * the same manner. In order to use these algorithms the normalization must
     * also be modified to reflect the symmetry.
     * 
     * Output values have been checked against those obtained using the internal
     * IDL legendre() function to obtain the various associated legendre
     * polynomials.
     * 
     * As stated above, I think that this normalization may be unnecessary. The
     * important thing is that the various spherical harmonics are orthogonal,
     * rather than orthonormal.
     * 
     */
  
    /* determine array of factorials */
    /*fact = (long unsigned *)malloc(sizeof(long unsigned)*(2*order+2));*/
    fact = (double *)malloc(sizeof(double)*(2*order+2));
    if (fact == NULL) return (-1);
    fact[0] = fact[1] = 1;
    for (k=2; k <= 2*order+1; k++) fact[k] = k*fact[k-1];
  
    ffff = (double *)malloc(sizeof(double)*(order+1)*(order+1));
    if (ffff == NULL) return (-1);
  
    /* determine normalization factors */
    for (l=0; l<=order; l++) {
      for (m=0; m<=l; m++) {
        k = l * (l+1) + m;      /* 1D index for l,m */
        ffff[k] = sqrt((2*l+1)/(4*M_PI) * fact[l-m]/fact[l+m]);
        ylmval[k] *= ffff[k];
      }
      for (m=-l; m<0; m++) {
        k  = l * (l+1) + m;     /* 1D index for l,m */
        kk = l * (l+1) - m;
        ylmval[k] *= ffff[kk] * ((-m % 2) ? -1 : 1);
      }
    }
  
    #if DEBUG > 1
    for (k=0; k<2*order+1; k++)
      /*printf("%03d %lu\n", k, fact[k]);*/
      printf("%03d %lf\n", k, fact[k]);
    #endif
  
    free(fact);
    free(ffff);
  
    #if DEBUG > 10
    for (l=0; l<=order; l++) {
      for (m=0; m<=l; m++) {
        k = l * (l+1) + m;      /* 1D index for l,m */
        printf("%03d %lf\n", k, ylmval[k]);
      }
    }
    #endif
  
    return (0);
  } 
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_Alt2CGM
  ;
  ; PURPOSE:
  ;       Transformation from so-called 'at-altitude' coordinates to AACGM.
  ;       The purpose of this function is to scale the latitudes in such a
  ;       way so that there is no gap. The problem is that for non-zero
  ;       altitudes (h) are range of latitudes near the equator lie on dipole
  ;       field lines that near reach the altitude h, and are therefore not
  ;       accessible. This is the inverse transformation.
  ;
  ;       cos (lat_aacgm) = sqrt( Re/(Re + h) ) cos (lat_at-alt)
  ;       
  ;
  ; CALLING SEQUENCE:
  ;       AACGM_v2_Alt2CGM,r_height_in,r_lat_alt,r_lat_adj
  ;
  ;     Input Arguments:  
  ;       r_height_in   - The altitude (h)
  ;       r_lat_alt     - The 'at-altitude' latitude
  ;
  ;     Output Arguments:  
  ;       r_lat_adj     - The corrected latitude, i.e., AACGM latitude
  ;
  ; HISTORY:
  ;
  ; This function is unchanged from its original version (Baker ?)
  ;     
  ;+-----------------------------------------------------------------------------
  */
  
  void AACGM_v2_Alt2CGM(double r_height_in, double r_lat_alt, double *r_lat_adj)
  {
    double eps =1e-9;
    double unim =0.9999999;
  
    double r1;
    double r0, ra;
  
    #if DEBUG > 0
    printf("AACGM_v2_Alt2CGM\n");
    #endif
  
    /* Computing 2nd power */
    r1 = cos(r_lat_alt*DTOR);
    ra = r1 * r1;
    if (ra < eps) ra = eps;
  
    r0 = (r_height_in/RE + 1.) / ra;
    if (r0 < unim) r0 = unim;
    
    r1 = acos(sqrt(1/r0));
    *r_lat_adj = AACGM_v2_Sgn(r1, r_lat_alt)/DTOR;
  } 
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_CGM2Alt
  ;
  ; PURPOSE:
  ;       Transformation from AACGM to so-called 'at-altitude' coordinates.
  ;       The purpose of this function is to scale the latitudes in such a
  ;       way so that there is no gap. The problem is that for non-zero
  ;       altitudes (h) are range of latitudes near the equator lie on dipole
  ;       field lines that near reach the altitude h, and are therefore not
  ;       accessible. This mapping closes the gap.
  ;
  ;       cos (lat_at-alt) = sqrt( (Re + h)/Re ) cos (lat_aacgm)
  ;       
  ;
  ; CALLING SEQUENCE:
  ;       AACGM_v2_CGM2Alt,r_height_in,r_lat_in,r_lat_adj, error
  ;
  ;     Input Arguments:  
  ;       r_height_in   - The altitude (h)
  ;       r_lat_in      - The AACGM latitude
  ;
  ;     Output Arguments:  
  ;       r_lat_adj     - The 'at-altitude' latitude
  ;       error         - variable is set if latitude is below the value that
  ;                       is mapped to the origin
  ;
  ; HISTORY:
  ;
  ; This function is unchanged from its original version (Baker ?)
  ;     
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_CGM2Alt(double r_height_in, double r_lat_in, double *r_lat_adj)
  {
    double unim=1;
    double r1;
    double ra;
    int error=0;
  
    #if DEBUG > 0
    printf("AACGM_v2_CGM2Alt\n");
    #endif
  
    /* convert from AACGM to at-altitude coordinates */
    r1 = cos(r_lat_in*DTOR);
    ra = (r_height_in/RE + 1)*(r1*r1);
    if (ra > unim) {
      ra = unim;
      error = 1;
    }
  
    r1 = acos(sqrt(ra));
    *r_lat_adj = AACGM_v2_Sgn(r1,r_lat_in)/DTOR;
  
    return (error);
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_Sgn
  ;
  ; PURPOSE:
  ;       return the signed quantity of a variable where the magnitude is given
  ;       by the first argument and the sign is given by the second argument.
  ;
  ; CALLING SEQUENCE:
  ;       AACGM_v2_Sgn, a, b
  ;     
  ;     Input Arguments:  
  ;       a             - magnitude
  ;       b             - sign
  ;
  ;     Return Value:
  ;       signed quantity
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  double AACGM_v2_Sgn(double a, double b)
  {
    double x=0;
    x = (a >= 0) ? a : -a;
    return (b >= 0) ? x: -x;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       convert_geo_coord_v2
  ;
  ; PURPOSE:
  ;       Second-level function used to determine the lat/lon of the input
  ;       coordinates.
  ;
  ; CALLING SEQUENCE:
  ;       err = convert_geo_coord_v2(in_lat,in_lon,height, out_lat,out_lon,
  ;                               code,order);
  ;     
  ;     Input Arguments:  
  ;       in_lat        - latitude in degrees
  ;       in_lon        - longitude in degrees
  ;       height        - height above Earth in km
  ;       code          - bitwise code for passing options into converter
  ;                       G2A         - geographic (geodetic) to AACGM-v2
  ;                       A2G         - AACGM-v2 to geographic (geodetic)
  ;                       TRACE       - use field-line tracing, not coefficients
  ;                       ALLOWTRACE  - use trace only above 2000 km
  ;                       BADIDEA     - use coefficients above 2000 km
  ;       order         - integer order of spherical harmonic expansion
  ;
  ;     Output Arguments:
  ;       out_lat       - pointer to output latitude in degrees
  ;       out_lon       - pointer to output longitude in degrees
  ;
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int convert_geo_coord_v2(double lat_in, double lon_in, double height_in,
                        double *lat_out, double *lon_out, int code, int order) {
  
  /*  int i,j,k,l,m,f,a,t,flag;*/
    int i,j,k,l,m,flag;
    int i_err64, err;
  
  /*  extern int rylm(); */
    double ylmval[AACGM_KMAX];
    double colat_temp;
    double lon_output;
      
    double lat_adj=0;
  /*  double lat_alt=0; */
    double colat_input; 
     
    double alt_var_cu=0, lon_temp=0, alt_var_sq=0, alt_var_qu=0;
    double colat_output=0, r=0, x=0, y=0, z=0;
    double ztmp, fac;
    double alt_var=0;
    double lon_input=0;
  
    static double cint[AACGM_KMAX][NCOORD][NFLAG];
  
    #if DEBUG > 0
    printf("convert_geo_coord_v2\n");
    #endif
  
    /* no date/time set so use current time */
    if (aacgm_date.year < 0) {    /* AACGM_v2_SetNow();*/
      AACGM_v2_errmsg(0);
      return -128;
    }
  
    /* TRACE */ /* call tracing functions here and return */
    if ((code & TRACE) || (height_in > MAXALT && (code & ALLOWTRACE))) {
      if (A2G & code) {   /* AACGM-v2 to geographic */
        err = AACGM_v2_Trace_inv(lat_in,lon_in,height_in, lat_out,lon_out);
  
        /* v2.3 moved to AACGM_v2_Convert
        if ((code & GEOCENTRIC) == 0) {
          geoc2geod(*lat_out,*lon_out,(RE+height_in)/RE, llh);
          *lat_out = llh[0];
        } */
      } else {
        err = AACGM_v2_Trace(lat_in,lon_in,height_in, lat_out,lon_out);
      }
  
      return (err);
    }
  
    /* determine the altitude dependence of the coefficients */
    flag = (A2G & code);    /* 0 for G2A; 1 for A2G */
    if (height_in != height_old[flag]) {
      alt_var = height_in/(double)MAXALT;
      alt_var_sq = alt_var * alt_var;
      alt_var_cu = alt_var * alt_var_sq;
      alt_var_qu = alt_var * alt_var_cu;
  
      #if DEBUG > 1
      printf("alt_var = %lf\n", alt_var);
      printf("alt_var_qu = %lf\n", alt_var_qu);
      #endif
  
      #if DEBUG > 0
      printf("** HEIGHT INTERPOLATION **\n");
      #endif
  
      for (i=0; i<NCOORD; i++) {
        for (j=0; j<AACGM_KMAX;j++) {
          /* change to allow general polynomial approximation */
          cint[j][i][flag] =  sph_harm_model.coef[j][i][0][flag] +
                              sph_harm_model.coef[j][i][1][flag]*alt_var+
                              sph_harm_model.coef[j][i][2][flag]*alt_var_sq+
                              sph_harm_model.coef[j][i][3][flag]*alt_var_cu+
                              sph_harm_model.coef[j][i][4][flag]*alt_var_qu;
          #if DEBUG > 10
          printf("%35.30lf %35.30lf\n", cint[j][i][flag],
                              sph_harm_model.coef[j][i][0][flag]);
          #endif
  
        }
      }
  
      height_old[flag] = height_in;
    }
    #if DEBUG > 1
    printf("cint[0][0][%d] = %lf\n", flag, cint[0][0][flag]);
    printf("cint[%d][0][%d] = %lf\n",
              AACGM_KMAX-1, flag, cint[AACGM_KMAX-1][0][flag]);
    printf("cint[%d][%d][%d] = %lf\n", AACGM_KMAX-1, NCOORD-1, flag,
                                        cint[AACGM_KMAX-1][NCOORD-1][flag]);
    #endif
  
    x = y = z = 0;
  
    lon_input = lon_in*DTOR;
  
    if (flag == 0) {
      colat_input = (90.-lat_in)*DTOR;
    } else {
      /* use intermediate "at-altitude" coordinates for inverse trans. */
      i_err64 = AACGM_v2_CGM2Alt(height_in, lat_in, &lat_adj);
  
      if (i_err64 != 0) return -64;
      colat_input= (90. - lat_adj)*DTOR;
    }
  
    /* Compute the values of the spherical harmonic functions.
     * NOTE: this function was adapted to use orthonormal SH functions */
    AACGM_v2_Rylm(colat_input, lon_input, order, ylmval);
  
    for (l=0; l<=order; l++) {
      for (m=-l; m<=l; m++) {
        k = l * (l+1) + m;      /* SGS: changes indexing */
  
        x += cint[k][0][flag]*ylmval[k];
        y += cint[k][1][flag]*ylmval[k];
        z += cint[k][2][flag]*ylmval[k];
      }
    }
   
    /* COMMENT: SGS
     * 
     * This answers one of my questions about how the coordinates for AACGM are
     * guaranteed to be on the unit sphere. Here they compute xyz indpendently
     * using the SH coefficients for each coordinate. They reject anything that
     * is +/- .1 Re from the surface of the Earth. They then scale each xyz
     * coordinate by the computed radial distance. This is a TERRIBLE way to do
     * things... but necessary for the inverse transformation.
     */
  
    /* SGS - new method that ensures position is on unit sphere and results in a
     *       much better fit. Uses z coordinate only for sign, i.e., hemisphere.
     */
    if (flag == 0) {
      fac = x*x + y*y;
      if (fac > 1.) {
        /* we are in the forbidden region and the solution is undefined */
        *lat_out = HUGE_VAL;
        *lon_out = HUGE_VAL;
  
        return -64;
      }
  
      ztmp = sqrt(1. - fac);
      z = (z < 0) ? -ztmp : ztmp;
  
      colat_temp = acos(z);
    } else {
    /* SGS - for inverse the old normalization produces lower overall errors...*/
      r = sqrt(x*x + y*y + z*z);
      if ((r< 0.9) || (r > 1.1)) return -32;
  
      z /= r;
      x /= r;
      y /= r;
  
      if (z > 1.) colat_temp = 0;
      else if (z < -1.) colat_temp = M_PI;
      else colat_temp = acos(z);
    }
    
    if ((fabs(x) < 1e-8) && (fabs(y) < 1e-8)) lon_temp = 0;
    else lon_temp = atan2(y,x);
  
    lon_output = lon_temp;
  
    /* SGS - not used for forward transformation */
    /*
    if (flag == 0) {
      lat_alt =90 - colat_temp*180/PI;
      altitude_to_cgm(height_in, lat_alt,&lat_adj);
      colat_output = (90. - lat_adj) * PI/180;
    } else colat_output = colat_temp;
    */
    colat_output = colat_temp;
  
    *lat_out = (double) (90. - colat_output/DTOR);
    *lon_out = (double) (lon_output/DTOR);
  
    /* v2.3 moved to AACGM_v2_Convert
    if ((code & GEOCENTRIC) == 0 && (code & A2G)) {
      geoc2geod(*lat_out,*lon_out,(RE+height_in)/RE, llh);
      *lat_out = llh[0];
    } */
  
    return 0;
  } 
  
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_LoadCoefFP
  ;
  ; PURPOSE:
  ;       Load a set of spherical harmonic coefficients.
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_LoadCoefFP(fp, code);
  ;     
  ;     Input Arguments:  
  ;       fp            - FILE pointer to open coefficient file
  ;       code          - 0 for 1st set; 1 for 2nd set
  ;
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_LoadCoefFP(FILE *fp, int code)
  {
    /*  char tmp[64]; */
    double tmp;
    int f,l,a,t;
  
    #if DEBUG > 0
    printf("AACGM_v2_LoadCoefFP\n");
    #endif
  
    if (fp==NULL) {
      #if DEBUG > 0
      printf("NULL file pointer in AACGM_v2_LoadCoefFP\n");
      #endif
      return -1;
    }
  
    for (f=0;f<NFLAG;f++) { 
      for (l=0;l<POLYORD;l++) {
        for (a=0;a<NCOORD;a++) { 
          for (t=0;t<AACGM_KMAX;t++) {
            if (fscanf(fp, "%lf", &tmp) != 1) {
              #if DEBUG > 0
              printf("FILE error, aborting\n");
              #endif
              fclose(fp);
              return -1;
            }
  
            sph_harm_model.coefs[t][a][l][f][code] = tmp;
          }
        }
      }
    }
  
    #if DEBUG > 10
    for (f=0;f<NFLAG;f++) { 
      for (l=0;l<POLYORD;l++) {
        for (a=0;a<NCOORD;a++) { 
          for (t=0;t<AACGM_KMAX;t++) {
            printf("%lf ", sph_harm_model.coefs[t][a][l][f][code]);
          }
          printf("\n");
        }
        printf("\n");
      }
      printf("\n");
    }
    #endif
  
  return 0;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_LoadCoef
  ;
  ; PURPOSE:
  ;       Load a set of spherical harmonic coefficients.
  ;
  ; CALLING SEQUENCE:
  ;       ret = AACGM_v2_LoadCoef(fname,code);
  ;     
  ;     Input Arguments:  
  ;       fname         - filename containing the AACGM coefficient set
  ;       code          - 0 for 1st set; 1 for 2nd set
  ;
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_LoadCoef(char *fname, int code)
  {
    FILE *fp=NULL;
    int err=-2;
  
    #if DEBUG > 0
    printf("loading %s\n", fname);
    #endif
  
    fp = fopen(fname,"r");
    if (fp==NULL) {
      #if DEBUG > 0
      printf("NULL file pointer in AACGM_v2_LoadCoef\n");
      #endif
      return -1;
    }
  
    err = AACGM_v2_LoadCoefFP(fp, code);
    if (err != 0) {
      #if DEBUG > 0
      printf("error in AACGM_v2_LoadCoefFP\n");
      #endif
      return -2;
    }
    fclose(fp);
  
    return 0;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_LoadCoefs
  ;
  ; PURPOSE:
  ;       Load two sets of spherical harmonic coefficients.
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_LoadCoefs(myear);
  ;     
  ;     Input Arguments:  
  ;       myear         - 5-year epoch year prior to desired time; bracketing
  ;                       set is +5 years.
  ;
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_LoadCoefs(int year)
  {
    char fname[256];
    char root[256];
    char yrstr[5];  
    int ret=0;
  
    #if DEBUG > 0
    printf("AACGM_v2_LoadCoefs\n");
    #endif
    /* default location of coefficient files */
    strcpy(root,getenv("AACGM_v2_DAT_PREFIX"));  
    if (strlen(root)==0) {
      AACGM_v2_errmsg(2);
      return -1;
    }
  
    if (year <= 0) return -1;
    sprintf(yrstr,"%4.4d",year);  
  
    strcpy(fname,root);
    strcat(fname,yrstr);
    strcat(fname,".asc");
    #if DEBUG > 0
    printf("AACGM_v2_LoadCoefs: %s\n", fname);
    #endif
    ret = AACGM_v2_LoadCoef(fname,G2A);   /* forward coefficients */
    if (ret != 0) return ret;
  
    sprintf(yrstr,"%4.4d",year+5);  
    strcpy(fname,root);
    strcat(fname,yrstr);
    strcat(fname,".asc");
    #if DEBUG > 0
    printf("AACGM_v2_LoadCoefs: %s\n", fname);
    #endif
    ret += AACGM_v2_LoadCoef(fname,A2G);  /* inverse coefficients */
  
    myear_old = year;
  
    return ret;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_Convert
  ;
  ; PURPOSE:
  ;       Main function called by many SD plotting routines that are written
  ;       in C.
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_Convert(in_lat, in_lon, height,
  ;                 out_lat, out_lon, r, code);
  ;     
  ;     Input Arguments:  
  ;       int_lat       - double precision input latitude in degrees
  ;       int_lon       - double precision input longitude in degrees
  ;       height        - altitude in km
  ;       code          - bitwise code for passing options into converter
  ;                       G2A         - geographic (geodetic) to AACGM-v2
  ;                       A2G         - AACGM-v2 to geographic (geodetic)
  ;                       TRACE       - use field-line tracing, not coefficients
  ;                       ALLOWTRACE  - use trace only above 2000 km
  ;                       BADIDEA     - use coefficients above 2000 km
  ;                       GEOCENTRIC  - assume inputs are geocentric w/ RE=6371.2
  ;
  ;     Output Arguments:  
  ;       out_lat       - output latitude in degrees
  ;       out_lon       - output longitude in degrees
  ;       r             - geocentric radial distance in Re
  ;
  ;     Return Value:
  ;       error code
  ;
  ;
  ; NOTES:
  ;
  ;       All AACGM-v2 conversions are done in geocentric coordinates using a
  ;           value of 6371.2 km for the Earth radius.
  ;
  ;       For G2A conversion inputs are geographic latitude, longitude and
  ;           height (glat,glon,height), specified as either geocentric or
  ;           geodetic (default). For geodetic inputs a conversion to geocentric
  ;           coordinates is performed, which changes the values of
  ;           glat,glon,height. The output is AACGM-v2 latitude, longitude and
  ;           the geocentric radius (mlat,mlon,r) using the geocentric height
  ;           in units of RE.
  ;
  ;        For A2G conversion inputs are AACGM-v2 latitude, longitude and the
  ;            geocentric height (mlat,mlon,height). The latter can be obtained
  ;            from the r output of the G2A conversion. The output is geographic
  ;            latitude, longitude and height (glat,glon,height). If the
  ;            gedodetic option is desired (default) a conversion of the outputs
  ;            is performed, which changes the values of glat,glon,height.
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_Convert(double in_lat, double in_lon, double height,
                    double *out_lat, double *out_lon, double *r, int code)
  {
    int err;
    int order=10;   /* pass in so a lower order would be allowed? */
    double rtp[3];
    double llh[3];
  
    #if DEBUG > 0
    printf("AACGM_v2_Convert\n");
    #endif
  
    /* latitude out of bounds */
    if (fabs(in_lat) > 90.) {
      fprintf(stderr, "ERROR: latitude must be in the range -90 to +90 degrees: "
                      "%lf\n", in_lat);
      return -8;
    }
  
    /* longitude out of bounds */
  /* SGS v2.3 removing requirement that longitude be -180 to 180. Does not seems
   *          to matter and is inconsistent with IDL version: -180 to 180.
  
    if ((in_lon < -180) || (in_lon > 180)) {
      fprintf(stderr, "ERROR: longitude must be in the range -180 to 180 "
                      "degrees: %lf\n", in_lon);
      return -16;
    }
   */
  
    /* if forward calculation (G2A) and input coordinates are given in geodetic
       coordinates (default) then must first convert to geocentric coordinates */
    if ((code & GEOCENTRIC) == 0 && (code & A2G) == 0) {
      geod2geoc(in_lat,in_lon,height, rtp);
  
      /* modify lat/lon/alt to geocentric values */
      in_lat = 90. - rtp[1]/DTOR;
      in_lon = rtp[2]/DTOR;
      height = (rtp[0]-1.)*RE;
    }
  
    /* height < 0 km */
    if ((height < 0) && (code&VERBOSE)) {
      fprintf(stderr, "WARNING: coordinate transformations are not intended "
                      "for altitudes < 0 km: %lf\n", height);
    /*  return -2; */
    }
  
    /* height > 2000 km not allowed for coefficients */
    if (height > MAXALT && !(code & (TRACE|ALLOWTRACE|BADIDEA))) {
      fprintf(stderr, "ERROR: coefficients are not valid for altitudes "
                      "above %d km: %lf.\n", MAXALT, height);
      fprintf(stderr, "       You must either use field-line tracing "
                      "(TRACE or ALLOWTRACE) or\n"
                      "       indicate that you know this is a very bad idea "
                      "(BADIDEA)\n\n");
      return -4;
    }
  
    /* all inputs are geocentric */
    err = convert_geo_coord_v2(in_lat,in_lon,height, out_lat,out_lon, code,order);
    /* all outputs are geocentric */
  
    if ((code & A2G) == 0) {    /* forward: G2A */
      *r = (height + RE)/RE;    /* geocentric radial distance in RE */
    } else {                    /* inverse: A2G */
      if ((code & GEOCENTRIC) == 0) { /* geodetic outputs */
        geoc2geod(*out_lat,*out_lon,(RE+height)/RE, llh);
        *out_lat = llh[0];
        height = llh[2];
      }
      *r = height;              /* height in km */
    }
  
    if (err !=0) return -1;
    return 0;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_SetDateTime
  ;
  ; PURPOSE:
  ;       Function to set date and time. MUST be called at least once BEFORE
  ;       any calls to AACGM_v2_Convert.
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_SetDateTime(year, month, day, hour, minute, second);
  ;     
  ;     Input Arguments:  
  ;       year          - year [1900-2025)
  ;       month         - month of year [01-12]
  ;       day           - day of month [01-31]
  ;       hour          - hour of day [00-24]
  ;       minute        - minute of hour [00-60]
  ;       second        - second of minute [00-60]
  ;
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_SetDateTime(int year, int month, int day,
                        int hour, int minute, int second)
  {
    int err, doy, ndays;
    double fyear;
  
    doy = dayno(year,month,day,&ndays);
    fyear = year + ((doy-1) + (hour + (minute + second/60.)/60.)/24.) / ndays;
  
    if ((fyear < IGRF_FIRST_EPOCH) || (fyear >= IGRF_LAST_EPOCH + 5.)) {
      AACGM_v2_errmsg(1);
      return (-1);
    }
  
    aacgm_date.year   = year;
    aacgm_date.month  = month;
    aacgm_date.day    = day;
    aacgm_date.hour   = hour;
    aacgm_date.minute = minute;
    aacgm_date.second = second;
    aacgm_date.dayno  = doy;
    aacgm_date.daysinyear = ndays;
  
    #if DEBUG > 0
    printf("AACGM_v2_SetDateTime\n");
    printf("%03d: %04d%02d%02d %02d%02d:%02d\n",
          aacgm_date.dayno, aacgm_date.year, aacgm_date.month, aacgm_date.day,
          aacgm_date.hour, aacgm_date.minute, aacgm_date.second);
    #endif
  
    err = AACGM_v2_TimeInterp();
  
    return err;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_GetDateTime
  ;
  ; PURPOSE:
  ;       Function to get date and time.
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_GetDateTime(year, month, day, hour, minute, second, dayno);
  ;     
  ;     Output Arguments (integer pointers):  
  ;       year          - year [1965-2014]
  ;       month         - month of year [01-12]
  ;       day           - day of month [01-31]
  ;       hour          - hour of day [00-24]
  ;       minute        - minute of hour [00-60]
  ;       second        - second of minute [00-60]
  ;       dayno         - day of year [01-366]
  ;
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_GetDateTime(int *year, int *month, int *day,
                        int *hour, int *minute, int *second, int *dayno)
  {
    *year   = aacgm_date.year;
    *month  = aacgm_date.month;
    *day    = aacgm_date.day;
    *hour   = aacgm_date.hour;
    *minute = aacgm_date.minute;
    *second = aacgm_date.second;
    *dayno  = aacgm_date.dayno;
  
    return 0;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_SetNow
  ;
  ; PURPOSE:
  ;       Function to set date and time to current computer time in UT.
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_SetNow();
  ;     
  ;     Return Value:
  ;       error code
  ;
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_SetNow(void)
  {
    /* current time */
    int err, doy,ndays;
    double fyear;
    time_t now;
    struct tm *tm_now;
  
    time(&now);
    tm_now = gmtime(&now);    /* right now in UT */
  
    doy = dayno(1900 + tm_now->tm_year,tm_now->tm_mon,tm_now->tm_mday,&ndays);
    fyear = 1900 + tm_now->tm_year + ((doy-1) + (tm_now->tm_hour +
                  (tm_now->tm_min + tm_now->tm_sec/60.)/60.)/24.) / ndays;
  
    if ((fyear < IGRF_FIRST_EPOCH) || (fyear >= IGRF_LAST_EPOCH + 5.)) {
      fprintf(stderr, "\nDate range for AACGM-v2 is [%4d - %4d)\n\n",
                        IGRF_FIRST_EPOCH, IGRF_LAST_EPOCH + 5);
      fprintf(stderr, "%04d%02d%02d %02d%02d:%02d\n", tm_now->tm_year,
                       tm_now->tm_mon, tm_now->tm_mday, tm_now->tm_hour,
                       tm_now->tm_min, tm_now->tm_sec);
      return (-1);
    }
  
    aacgm_date.year   = (*tm_now).tm_year + 1900;
    aacgm_date.month  = (*tm_now).tm_mon  + 1;
    aacgm_date.day    = (*tm_now).tm_mday;
    aacgm_date.hour   = (*tm_now).tm_hour;
    aacgm_date.minute = (*tm_now).tm_min;
    aacgm_date.second = (*tm_now).tm_sec;
    aacgm_date.dayno  = (*tm_now).tm_yday + 1;
    aacgm_date.daysinyear = ndays;
  
    #if DEBUG > 0
    printf("AACGM_v2_SetNow\n");
    printf("%03d: %04d%02d%02d %02d%02d:%02d\n",
          aacgm_date.dayno, aacgm_date.year, aacgm_date.month, aacgm_date.day,
          aacgm_date.hour, aacgm_date.minute, aacgm_date.second);
    #endif
  
    err = AACGM_v2_TimeInterp();
  
    return err;
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_errmsg
  ;
  ; PURPOSE:
  ;       Display error message because no date and time have been set.
  ;
  ; CALLING SEQUENCE:
  ;       AACGM_v2_errmsg(code);
  ;     
  ;+-----------------------------------------------------------------------------
  */
  
  void AACGM_v2_errmsg(int ecode)
  {
  
    fprintf(stderr, "\n"
    "**************************************************************************"
    "\n");
  
    switch (ecode) {
      case 0: /* no Date/Time set */
    fprintf(stderr,
    "* AACGM-v2 ERROR: No Date/Time Set                                       *\n"
    "*                                                                        *\n"
    "* You must specifiy the date and time in order to use AACGM coordinates, *\n"
    "* which depend on the internal (IGRF) magnetic field. Before calling     *\n"
    "* AACGM_v2_Convert() you must set the date and time to the integer values*\n"
    "* using the function:                                                    *\n"
    "*                                                                        *\n"
    "*   AACGM_v2_SetDateTime(year,month,day,hour,minute,second);             *\n"
    "*                                                                        *\n"
    "* or to the current computer time in UT using the function:              *\n"
    "*                                                                        *\n"
    "*   AACGM_v2_SetNow();                                                   *\n"
    "*                                                                        *\n"
    "* subsequent calls to AACGM_v2_Convert() will use the last date and time *\n"
    "* that was set, so update to the actual date and time that is desired.   *"
    "\n");
      break;
  
      case 1: /* Date/Time out of bounds */
    fprintf(stderr,
    "* AACGM-v2 ERROR: Date out of bounds                                     *\n"
    "*                                                                        *\n"
    "* The current date range for AACGM-v2 coordinates is [1990-2025), which  *\n"
    "* corresponds to the date range for the IGRF12 model, including the      *\n"
    "* 5-year secular variation.                                              *"
    "\n");
      break;
  
      case 2: /* COEF Path not set */
    fprintf(stderr,
    "* AACGM-v2 ERROR: AACGM_v2_DAT_PREFIX path not set *\n"
    "*                                                                        *\n"
    "* You must set the environment variable AACGM_v2_DAT_PREFIX to the       *\n"
    "\n");
      break;
    }
    fprintf(stderr,
    "**************************************************************************"
    "\n\n");
  }
  
  /*-----------------------------------------------------------------------------
  ;
  ; NAME:
  ;       AACGM_v2_TimeInterp
  ;
  ; PURPOSE:
  ;       Interpolate coefficients between adjacent 5-year epochs
  ;
  ; CALLING SEQUENCE:
  ;       err = AACGM_v2_TimeInterp();
  ;     
  ;+-----------------------------------------------------------------------------
  */
  
  int AACGM_v2_TimeInterp(void)
  {
    int myear,f,l,a,t,err;
    double fyear;
  
    /* myear is the epoch model year */
    myear = aacgm_date.year/5*5;
    if (myear != myear_old) {   /* load the new coefficients, if needed */
      err = AACGM_v2_LoadCoefs(myear);
      if (err != 0) return err;
      fyear_old = -1;           /* force time interpolation */
      height_old[0] = -1.;      /* force height interpolation */
      height_old[1] = -1.;
    }
  
    /* fyear is the floating point time */
    fyear = aacgm_date.year + ((aacgm_date.dayno-1) + (aacgm_date.hour +
                      (aacgm_date.minute + aacgm_date.second/60.)/60.)/24.)/
                      aacgm_date.daysinyear;
  
    /* time interpolation right here */
    if (fyear != fyear_old) {
      #if DEBUG > 0
      printf("** TIME INTERPOLATION **\n");
      #endif
  
      for (f=0;f<NFLAG;f++)
      for (l=0;l<POLYORD;l++)
      for (a=0;a<NCOORD;a++)
      for (t=0;t<AACGM_KMAX;t++)
        sph_harm_model.coef[t][a][l][f] = sph_harm_model.coefs[t][a][l][f][0] +
            (fyear - myear) * (sph_harm_model.coefs[t][a][l][f][1] -
                              sph_harm_model.coefs[t][a][l][f][0])/5;
  
      height_old[0] = -1.;      /* force height interpolation because coeffs */
      height_old[1] = -1.;      /* have changed */
  
      fyear_old = fyear;
    }
  
    return (0);
  }
  
  
  int AACGM_v2_Trace(double lat_in, double lon_in, double alt,
                      double *lat_out, double *lon_out)
  {
    int err, kk, idir, below;
    unsigned long k,niter;
    double ds, dsRE, dsRE0, eps, Lshell;
    double rtp[3],xyzg[3],xyzm[3],xyzc[3],xyzp[3];
  
    /* set date for IGRF model */
    IGRF_SetDateTime(aacgm_date.year, aacgm_date.month, aacgm_date.day,
                      aacgm_date.hour, aacgm_date.minute, aacgm_date.second);
  
    /* Q: these could eventually be command-line options */
    ds    = 1.;
    dsRE  = ds/RE;
    dsRE0 = dsRE;
    eps   = 1.e-4/RE;
  
    /* for the model we are doing the tracing in geocentric coordinates */
    rtp[0] = (RE+alt)/RE;       /* distance in RE; 1.0 is surface of sphere */
    rtp[1] = (90.-lat_in)*DTOR;   /* colatitude in radians */
    rtp[2] = lon_in*DTOR;         /* longitude in radians */
  
    k = 0L;
    /* convert position to Cartesian coords */
    sph2car(rtp, xyzg);
  
    /* convert to magnetic Dipole coordinates */
    geo2mag(xyzg, xyzm);
  
    idir = (xyzm[2] > 0.) ? -1 : 1;   /* N or S hemisphere */
  
    dsRE = dsRE0;
  
    /*
    ; trace to magnetic equator
    ;
    ; Note that there is the possibility that the magnetic equator lies
    ; at an altitude above the surface of the Earth but below the starting
    ; altitude. I am not certain of the definition of CGM, but these
    ; fieldlines map to very different locations than the solutions that
    ; lie above the starting altitude. I am considering the solution for
    ; this set of fieldlines as undefined; just like those that lie below
    ; the surface of the Earth.
    ;
    ; Added a check for when tracing goes below altitude so as not to contiue
    ; tracing beyond what is necessary.
    ;
    ; Also making sure that stepsize does not go to zero
    */
    below = 0;
    while (!below && idir*xyzm[2] < 0.) {
  
      for (kk=0;kk<3;kk++) xyzp[kk] = xyzg[kk]; /* save as previous */
  
      AACGM_v2_RK45(xyzg, idir, &dsRE, eps, 1); /* set to 0 for RK4: /noadapt) */
  
      /* make sure that stepsize does not go to zero */
      if (dsRE*RE < 1e-2) dsRE = 1e-2/RE;
  
      /* convert to magnetic Dipole coordinates */
      geo2mag(xyzg, xyzm);
  
      below = ((xyzg[0]*xyzg[0]+xyzg[1]*xyzg[1]+xyzg[2]*xyzg[2]) <
               (RE+alt)*(RE+alt)/(RE*RE));
      k++;
    }
    niter = k;
  
    if (!below && niter > 1) {
      /* now bisect stepsize (fixed) to land on magnetic equator w/in 1 m */
      for (k=0;k<3;k++) xyzc[k] = xyzp[k];
      kk = 0L;
      while (dsRE > 1e-3/RE) {
        dsRE *= .5;
        for (k=0;k<3;k++) xyzp[k] = xyzc[k];
        AACGM_v2_RK45(xyzc, idir, &dsRE, eps, 0);  /* using RK4 */
        geo2mag(xyzc,xyzm);
  
        /* Is it possible that resetting here causes a doubling of the tol? */
        if (idir * xyzm[2] > 0)
        for (k=0;k<3;k++) xyzc[k] = xyzp[k];
  
        kk++;
      }
      niter += kk;
    } else {
      /*if (below) printf("BELOW\n");*/
      for (k=0;k<3;k++) xyzc[k] = xyzg[k];    /* just use last value */
    }
  
    /* 'trace' back to reference surface along Dipole field lines */
    Lshell = sqrt(xyzc[0]*xyzc[0] + xyzc[1]*xyzc[1] + xyzc[2]*xyzc[2]);
    if (Lshell < (RE+alt)/RE) { /* magnetic equator is below ... */
      *lat_out = NAN;
      *lon_out = NAN;
  
      err = -1;
    } else {
      geo2mag(xyzc, xyzm);  /* geographic to magnetic */
      car2sph(xyzm, rtp);
  
      *lat_out = -idir*acos(sqrt(1./Lshell))/DTOR;
      *lon_out = rtp[2]/DTOR;
      if (*lon_out > 180) *lon_out -= 360.;
  
      err = 0;
    }
  
    return (err);
  }
  
  int AACGM_v2_Trace_inv(double lat_in, double lon_in, double alt,
                      double *lat_out, double *lon_out)
  {
    int err, kk, idir;
    unsigned long k,niter;
    double ds, dsRE, dsRE0, eps, Lshell;
    double rtp[3],xyzg[3],xyzm[3],xyzc[3],xyzp[3];
  
    /* set date for IGRF model */
    IGRF_SetDateTime(aacgm_date.year, aacgm_date.month, aacgm_date.day,
                      aacgm_date.hour, aacgm_date.minute, aacgm_date.second);
  
    /* Q: these could eventually be command-line options */
    ds    = 1.;
    dsRE  = ds/RE;
    dsRE0 = dsRE;
    eps   = 1.e-4/RE;
  
    /* Q: Test this */
    /* poles map to infinity */
    if (fabs(fabs(lat_in) - 90.) < 1e-6)
      lat_in += (lat_in > 0) ? -1e-6 : 1e-6;
  
    Lshell = 1./(cos(lat_in*DTOR)*cos(lat_in*DTOR));
    if (Lshell <(RE+alt)/RE) { /* solution does not exist; the starting
                                * position at the magnetic equator is below
                                * the altitude of interest */
      *lat_out = NAN;
      *lon_out = NAN;
      err = -1;
    } else {
      /* magnetic Cartesian coordinates of fieldline trace starting point */
      xyzm[0] = Lshell*cos(lon_in*DTOR);
      xyzm[1] = Lshell*sin(lon_in*DTOR);
      xyzm[2] = 0.;
  
      /* geographic Cartesian coordinates of starting point */
      mag2geo(xyzm, xyzg);
  
      /* geographic spherical coordinates of starting point */
      car2sph(xyzg,rtp);
  
      k = 0L;
  
      /* direction of trace is determined by the starting hemisphere? */
      idir = (lat_in > 0) ? 1 : -1;
  
      dsRE = dsRE0;
  
      /* trace back to altitude above Earth */
      while (rtp[0] >  (RE + alt)/RE) {
        for (kk=0;kk<3;kk++) xyzp[kk] = xyzg[kk]; /* save as previous */
  
        AACGM_v2_RK45(xyzg, idir, &dsRE, eps, 1); /* set to 0 for RK4: /noadapt)*/
  
        /* make sure that stepsize does not go to zero */
        if (dsRE*RE < 5e-1) dsRE = 5e-1/RE;
  
        car2sph(xyzg, rtp);
        k++;
      }
      niter = k;
  
      if (niter > 1) {
        /* now bisect stepsize (fixed) to land on magnetic equator w/in 1 m */
        for (k=0;k<3;k++) xyzc[k] = xyzp[k];
        kk = 0L;
        while (dsRE > 1e-3/RE) {
          dsRE *= .5;
          for (k=0;k<3;k++) xyzp[k] = xyzc[k];
          AACGM_v2_RK45(xyzc, idir, &dsRE, eps, 0);  /* using RK4 */
  
          car2sph(xyzc, rtp);
  
          if (rtp[0] < (RE + alt)/RE)
            for (k=0;k<3;k++) xyzc[k] = xyzp[k];
  
          kk++;
        }
        niter += kk;
      }
  
      *lat_out = 90. - rtp[1]/DTOR;
      *lon_out = rtp[2]/DTOR;
      if (*lon_out > 180) *lon_out -= 360.;
      err = 0;
    }
  
    return (err);
  }
  
  