Add tool for calculating L-shell

[ericthewizard]

In IDL, we have a tool for calculating the L-shell value given spacecraft position data (calculate_lshell.pro); we still need to add this

[nickssl]

In IDL SPEDAS, when calculating the L-shell value, we use the following call to geopack:

geopack_trace, gsm_re[1,i], gsm_re[2,i], gsm_re[3,i], 1, 0, out_foot_x, out_foot_y, out_foot_z,/igrf,/equator,/refine,rlim=60.0D

From the geopack manual, the two keywords /refine and /equator are the following:

REFINE: Refined foot point locations using trace routine by Vassilis Angelopoulos.

EQUATOR (REFINE only): Trace to equatorial plane. The equator is defined as the location on the field line where the radial magnetic field component reverses sign.

The python function trace() does not contain any similar keywords, perhaps we should add it.

https://github.com/tsssss/geopack/blob/master/geopack/geopack.py

[jameswilburlewis]

Related issue (implement trace routines): https://github.com/spedas/pyspedas/issues/609

I'm not sure if the 'refine' option is present in Tsyganenko's original Fortran code (which is the basis of the geopack python package). So we might need help from Haje or Vassilis if we want to have this option in Python.

[xandrd]

@jameswilburlewis, is Tsyganenko's original Fortran code traces the field lines?

I think what @nickssl is referring to is related to field tracing not the magnetic field model. As far as I know everyone uses there own tracing method e.g., IRBEM, Geopack, SpacePy (LANL), etc. It make sense adding a keyword which will define the method.

[xandrd]

I have created a prototype of the function in a separate branch: https://github.com/spedas/pyspedas/blob/calculate_lshell/pyspedas/geopack/calculate_lshell.py

[nickssl]

The main issue here is to make sure that IDL geopack_trace and python geopack trace give the same results.

[nickssl]

To calculate l-shell we need the /equator keyword in geopack trace() function, currently this is missing from the python implementation.

[xandrd]

There is something I don't understand. In SPEDAS we call geopack_trace without any keywords for external magnetic field. And we also set option par=0. According to documentation this option is ignored for dipole. However, we do set flag for internal magnetic field to be /igrf:

    ;trace field line from position to equator
    geopack_trace, gsm_re[1,i], gsm_re[2,i], gsm_re[3,i], 1, 0, out_foot_x, $
        out_foot_y, out_foot_z,/igrf,/equator,/refine,rlim=60.0D

I could not find any source code for IDL Geopack DLM. It is only available in a form of compiled downloadable library. The Geopack-08 FOTRAN code must explicitly specify the external field, e.g., parameter EXNAME:

C==============================================================================
C
      SUBROUTINE TRACE_08 (XI,YI,ZI,DIR,DSMAX,ERR,RLIM,R0,IOPT,PARMOD,
     * EXNAME,INNAME,XF,YF,ZF,XX,YY,ZZ,L,LMAX)
C
C  TRACES A FIELD LINE FROM AN ARBITRARY POINT OF SPACE TO THE EARTH'S
C  SURFACE OR TO A MODEL LIMITING BOUNDARY.

Since in IDL we do not explicitly specify the external field, I am not sure how does it work. The IDL Geopack DLM documentation does not explain what is the default option for geopack_trace without external magnetic field keywords.

Note, that in Python version of geopack has default external field specifyed as T89, e.g., exname='t89'. And python version of geopack essentially mimics the Geopack-08 FOTRAN code, except for the default call options.

def trace(xi,yi,zi,dir,rlim=10,r0=1,parmod=2,exname='t89',inname='igrf',maxloop=1000):
    """
    Traces a field line from an arbitrary point of space to the earth's surface or
    to a model limiting boundary.
    """

Here are my questions:

• What is default external magnetic field in geopack_trace function of IDL Geopack DLM?
• Or, is there an source code for IDL Geopack DLM?

[jameswilburlewis]

Haje hasn't released the source code he's using for the DLM. My understanding is that it's based mostly on Tsyganenko's Fortran code, with some extra functionality, contributed by Vassilis, to improve the accuracy of field line tracing and permit tracing to the equator. I guess there must be some default arguments in the DLM code for geopack_trace. The external model for geopack_trace to use is specified by setting one of the keywords T89, T96, T01, etc. to 1. The documentation doesn't mention which external model it defaults to, or what its model parameters default to. In the calculate_lshell code, the 'par' argument is set to the scalar 0. So I'm guessing the default model is "dipole", since 0 is not a valid value for T89's iopt parameter, and the other models need a a 10-element parameter array. It looks like that routine has only had minor updates since it was added in 2011....at some point, someone might have started to add geopack_2008 support, but didn't finish the job. All I see is the feature test, with no logic to implement it. It also seems to assume that the input points are north of the equator, because the step direction is always anti-parallel to the B vector.

In IDL SPEDAS we have ttrace2iono and ttrace2equator wrapper routines (which convert tplot variables to arrays, and do some other housekeeping), which call lower level routines trace2iono and trace2equator, which eventually call geopack_trace. These routines let you specify the external model, the model parameters, whether to use geopack_2008, etc., so I would consider those the preferred routines to do field line tracing, rather than calling geopack_trace directly.

The Python geopack library doesn't yet support tracing to the equator, so the first order of business will be to get that working. I hope I'll have something this week, based on scipy.integrate.RK45.