centered dipole + MAG frame, weird results

[atmosalex]

Hi, I find that a centered dipole magnetic field object returns strange results when using the MAG frame (sysaxes=6). The following code demonstrates the issue:

import IRBEM as ib
import numpy as np
import datetime

bmodel = 5

mf_GSM = ib.MagFields(options=[1,0,0,0,bmodel], verbose=False, kext='None', sysaxes=2)
mf_GEO = ib.MagFields(options=[1,0,0,0,bmodel], verbose=False, kext='None', sysaxes=1)
mf_GDZ = ib.MagFields(options=[1,0,0,0,bmodel], verbose=False, kext='None', sysaxes=0)
mf_MAG = ib.MagFields(options=[1,0,0,0,bmodel], verbose=False, kext='None', sysaxes=6)
coords = ib.Coords()

X = {}

t_ = datetime.datetime(2015, 1, 1)

pos_GSM = [5, 0, 0] #RE
print("GSM coord.:", pos_GSM, ", [RE]")
pos_GEO = coords.transform(t_, pos_GSM, 'GSM', 'GEO')[0]
print("GEO coord.:", pos_GEO, ", [RE]")
pos_GDZ = coords.transform(t_, pos_GSM, 'GSM', 'GDZ')[0]
print("GDZ coord.:", pos_GDZ, ", [km, deg., deg.]")
pos_MAG = coords.transform(t_, pos_GSM, 'GSM', 'MAG')[0]
print("MAG coord.:", pos_MAG, ", [RE]")

#define some constant common arguments:
X['dateTime'] = t_
maginput = {'Kp': 4 * 10}

#B with GSM coordinates:
X['x1'] = pos_GSM[0]
X['x2'] = pos_GSM[1]
X['x3'] = pos_GSM[2]
#X['sysaxes'] = 2 #ignored
output_dictionary = mf_GSM.get_field_multi(X, maginput)
print("1, get_field_multi results, GSM:")
for key in output_dictionary.keys():
    print("", key, output_dictionary[key])

#B with GEO coordinates:
X['x1'] = pos_GEO[0]
X['x2'] = pos_GEO[1]
X['x3'] = pos_GEO[2]
#X['sysaxes'] = 1 #ignored
output_dictionary = mf_GEO.get_field_multi(X, maginput)
print("2, get_field_multi results, GEO:")
for key in output_dictionary.keys():
    print("", key, output_dictionary[key])

#B with GDZ coordinates:
X['x1'] = pos_GDZ[0]
X['x2'] = pos_GDZ[1]
X['x3'] = pos_GDZ[2]
#X['sysaxes'] = 0 #ignored
output_dictionary = mf_GDZ.get_field_multi(X, maginput)
print("3, get_field_multi results, GDZ:")
for key in output_dictionary.keys():
    print("", key, output_dictionary[key])


#B with MAG coordinates:
X['x1'] = pos_MAG[0]
X['x2'] = pos_MAG[1]
X['x3'] = pos_MAG[2]
#X['sysaxes'] = 6 #ignored
output_dictionary = mf_MAG.get_field_multi(X, maginput)
print("4, get_field_multi results, MAG:")
for key in output_dictionary.keys():
    print("", key, output_dictionary[key])
sys.exit()

This produces the output (note the fourth set of results):

1, get_field_multi results, GSM:
 BxGEO [-258.11145069]
 ByGEO [-3.60580063]
 BzGEO [129.08818197]
 Bl [288.61441666]
2, get_field_multi results, GEO:
 BxGEO [-258.11145069]
 ByGEO [-3.60580063]
 BzGEO [129.08818197]
 Bl [288.61441666]
3, get_field_multi results, GDZ:
 BxGEO [-258.11145066]
 ByGEO [-3.60580062]
 BzGEO [129.08818203]
 Bl [288.61441665]
4, get_field_multi results, MAG:
 BxGEO [-59.52457249]
 ByGEO [-271.74478956]
 BzGEO [106.41965094]
 Bl [297.84819538]

However, when bmodel = 1 (eccentric tilted dipole) the magnetic field vector is the same for each result.

[atmosalex]

Some field line tracing of a "centered dipole" in the MAG frame...

[AntoineBrunet]

Hi, thanks for raising this issue.

Here is what is happening according to me:

• When doing the first conversion from GSM to MAG, the IGRF tilted dipole for the chosen date is used to define the MAG system.
• However, during the call to get_field_multi, the coordinates are converted back to GEO using the centered dipole axis to define the MAG system, hence giving incorrect results.

This is definitely a (somewhat serious) bug! I don't think it's particularly complex to fix, but it could have unintended side effects...

For reference, I think INIT_CD should be removed (and not called anymore in init_fields). The /dipigrf/ commons would then only be used for IGRF dipole informations. It should fix the issue.