Cotrans

Compare IDL code to python code using cotrans functions, including the IGRF-13 model.

IDL code: SPEDAS

; Compile contrans library
cotrans_lib
; Define some data
d = [[245.0, -102.0, 251.0], [775.0, 10.0, -10], [121.0, 545.0, -1.0], [304.65, -205.3, 856.1], [464.34, -561.55, -356.22]]
; Define times
t = [1577112800, 1577308800, 1577598800, 1577608800, 1577998800]
t0 = time_string(t)
t1 = time_struct(t0)
print, t0
; Compute direction of Earth's magnetic axis in GEO, using the IGRF model. 
cdipdir_vect,transpose(t1.year[*]),transpose(t1.doy[*]),gd1,gd2,gd3
print, gd1,gd2,gd3
; Compute sun direction in GEI system
csundir_vect,transpose(t1.year[*]),transpose(t1.doy[*]),transpose(t1.hour[*]),transpose(t1.min[*]),transpose(t1.sec[*]),gst,slong,sra,sdec,obliq
print,gst,slong,sra,sdec,obliq
; Compute GEI to GSE transformation.
tgeigse_vect,transpose(t1.year[*]),transpose(t1.doy[*]),transpose(t1.hour[*]),transpose(t1.min[*]),transpose(t1.sec[*]),transpose(d[0, *]),transpose(d[1, *]),transpose(d[2, *]),xgse,ygse,zgse
print,xgse,ygse,zgse
; Compute GSE to GSM transformation.
tgsegsm_vect,transpose(t1.year[*]),transpose(t1.doy[*]),transpose(t1.hour[*]),transpose(t1.min[*]),transpose(t1.sec[*]),transpose(d[0, *]),transpose(d[1, *]),transpose(d[2, *]),xgsm,ygsm,zgsm
print,xgsm,ygsm,zgsm

IDL results:

IDL> print, t0
2019-12-23/14:53:20
2019-12-25/21:20:00
2019-12-29/05:53:20
2019-12-29/08:40:00
2020-01-02/21:00:00
IDL> print, gd1,gd2,gd3
    0.0486864    0.0486846    0.0486811    0.0486811    0.0486776
    -0.156116    -0.156111    -0.156101    -0.156101    -0.156091
     0.986538     0.986539     0.986541     0.986541     0.986542
IDL> print,gst,slong,sra,sdec,obliq
      5.50119     0.944179      3.24176      3.97097     0.994296
      4.73823      4.77856      4.83825      4.84031      4.92061
      4.74045      4.78439      4.84933      4.85157      4.93861
    -0.408904    -0.408101    -0.405626    -0.405513    -0.399712
     0.409047     0.409047     0.409047     0.409047     0.409047
IDL> print,xgse,ygse,zgse
    0.0618591      45.9846     -480.516     -112.061      738.670
      245.080      773.652      182.717      321.559      318.591
      270.862     -13.1524     -217.683      867.127     -103.484
IDL> print,xgsm,ygsm,zgsm
      245.000      775.000      121.000      304.650      464.340
     -83.8585      11.7033      545.000     -118.216     -460.356
      257.629     -7.93938     0.929342      872.399     -479.899

Python code: pySPEDAS

from cotrans_lib import *
d = [[245.0, -102.0, 251.0], [775.0, 10.0, -10], [121.0, 545.0, -1.0], [304.65, -205.3, 856.1], [464.34, -561.55, -356.22]]
t = [1577112800, 1577308800, 1577598800, 1577608800, 1577998800]
a = cdipdir_vect(t)
b = csundir_vect(t)
gse = tgeigse_vect(t, d)
gsm = tgsegsm_vect(t, d)

Results:

a
Out[2]: 
(array([0.04868638, 0.04868462, 0.0486811 , 0.0486811 , 0.04867761]),
 array([-0.15611589, -0.15611097, -0.15610113, -0.15610113, -0.15609089]),
 array([0.98653812, 0.98653899, 0.98654072, 0.98654072, 0.98654251]))

b
Out[3]: 
(array([5.50118781, 0.94417896, 3.24175902, 3.9709706 , 0.9942959 ]),
 array([4.73823063, 4.77856341, 4.83825495, 4.84031354, 4.92060693]),
 array([4.74044499, 4.78438591, 4.84932932, 4.85156629, 4.93861416]),
 array([-0.40890358, -0.40810141, -0.40562582, -0.40551314, -0.39971157]),
 array([0.4090472 , 0.40904719, 0.40904716, 0.40904716, 0.40904714]))

gse
Out[4]: 
(array([ 6.17328386e-02,  4.59847513e+01, -4.80515919e+02, -1.12061112e+02, 7.38670168e+02]),
 array([245.07961052, 773.65200071, 182.71689404, 321.55872916, 318.59101371]),
 array([ 270.86155264,  -13.15235508, -217.68322895,  867.12698798, -103.48369801]))

gsm
Out[5]: 
(array([245.  , 775.  , 121.  , 304.65, 464.34]),
 array([ -83.85842685,   11.70325822,  545.00012517, -118.21614302, -460.35592829]),
 array([ 257.6291215 ,   -7.93937951,    0.92928139,  872.39913086, -479.89947926]))