Added the Boberg extension to the TS89 model for Kp > 6

drflei commented 3 years ago

Added the Boberg extension to TS89 model for Kp > 6 cases with iopt fixed at 7, but set dst = -200 nT for Kp=7, dst = -225 nT for Kp=8 and dst = -250 nT for Kp=9.

Test:

Slightly changed version of the testLStarOutput() in python/IRBEM/test_IRBEM.py ... def testLStarOutput(test_datetime = True): """ This test function will test is the make_lstar1() function works correctly. If you run this, the output should be the follwing.

{'MLT': [8.34044753112316], 'xj': [9.898414822276834], 'lstar': [-1e+31],
'Lm': [4.631806704496794], 'bmin': [268.5087756309121], 
'blocal': [39730.828875776126]}
"""
model = MagFields(options = [0,0,0,0,0], kext=4, verbose = True)
LLA = {}
LLA['x1'] = [450, 450]
LLA['x2'] = [55, 55]
LLA['x3'] = [255, 260]
if test_datetime:
    LLA['dateTime'] = [datetime.datetime(1990, 1, 1, 0,0,0), 
    datetime.datetime(1990, 1, 1, 0, 0, 0)]
else:
    LLA['dateTime'] = ['2015-02-02T06:12:43', '2015-02-02T06:12:43']
maginput = {'Kp':[0, 0]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[10, 10]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[20, 20]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[30, 30]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[40, 40]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[50, 50]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[60, 60]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[70, 70]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[80, 80]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)
maginput = {'Kp':[90, 90]}
model.make_lstar(LLA, maginput)
print(model.make_lstar_output)

Here are the execution output:

Running test: testLStarOutput() Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [5.726767773795413, -6.08881441322539], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48498.58885311038, 48665.86598627261], 'bmin': [150.96834279788308, 124.79880818678406], 'Lstar': [-1e+31, -1e+31], 'xj': [13.074669052872109, 14.073675533005904]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [5.7635479857400185, 6.138350342299773], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48503.610693187664, 48671.07403890945], 'bmin': [147.29752674756213, 120.80145719869908], 'Lstar': [-1e+31, -1e+31], 'xj': [13.17591726279675, 14.210069718562913]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [5.809872671789241, -6.203015666878452], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48509.2913033562, 48677.00748679707], 'bmin': [142.46568199157505, 115.71729195564826], 'Lstar': [-1e+31, -1e+31], 'xj': [13.303434483426177, 14.38811678674148]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [5.8338249661986525, 6.237881159894027], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48512.56951403513, 48680.511895893746], 'bmin': [140.0767020028152, 112.48991401343247], 'Lstar': [-1e+31, -1e+31], 'xj': [13.3693768234096, 14.484124716882935]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [5.941125828452119, -6.372435754886049], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48524.502641178085, 48692.68095773457], 'bmin': [130.53154442819468, 104.20227116714513], 'Lstar': [-1e+31, -1e+31], 'xj': [13.66475703768701, 14.85464981195754]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [6.174932760948292, 6.659065148639292], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48547.72054225265, 48716.41410984769], 'bmin': [108.96338641171882, 85.62662045330771], 'Lstar': [-1e+31, -1e+31], 'xj': [14.308462717202248, 15.644060658150499]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [6.5142096349690135, 7.060671006585493], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48571.29471376352, 48739.970185952676], 'bmin': [95.59812314466937, 77.07028738219243], 'Lstar': [-1e+31, -1e+31], 'xj': [15.242634824095866, 16.75026439855387]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [6.9925651410983605, 7.676748522787302], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48594.90379924422, 48764.11036171313], 'bmin': [72.490159586074, 56.87836871038282], 'Lstar': [-1e+31, -1e+31], 'xj': [16.560043970157178, 18.447583943962748]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [-7.424891434245514, 8.230030913956039], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48611.279379472326, 48780.85404333654], 'bmin': [58.723620448094735, 44.938714502644544], 'Lstar': [-1e+31, -1e+31], 'xj': [17.75098861165504, 19.972294720181697]} Prepping magnetic field inputs. Done prepping magnetic field inputs. Running IRBEM-LIB make_lstar {'Lm': [-7.957711921697585, 8.914272679952461], 'MLT': [15.854697006914213, 16.27533825441823], 'blocal': [48627.65803611577, 48797.60064965626], 'bmin': [46.68116415809574, 34.56608067397459], 'Lstar': [-1e+31, -1e+31], 'xj': [19.21911603059771, 21.858324625253953]}

Pull Request Checklist

New code is easier to review, integrate and maintain if it's clear, commented, and consistent with the style of the rest of the IRBEM code. Please try to follow as much of this checklist as you can for your PR. If you can't hit everything, or don't know how to, then submit the PR and the rest can be discussed during review. Some additional suggestions are given below.

Please also see our Code of Conduct so that the PRBEM community remains welcoming and inclusive.

Go ahead and replace the text above (and this line!) with your pull request description.

[ ] Pull request has descriptive title
[ ] Pull request gives overview of changes
[ ] New code has inline comments where necessary
[ ] Appropriate documentation has been written
[ ] Relevant issues are linked to (e.g. See issue # or Closes #)

AntoineBrunet commented 3 years ago

Hello @drflei , thank you for this PR. If I understand correctly, it implements Boberg et al. 1995, but modeling Dst as a function of Kp. I have the following remarks, and I'd like other people to also give their opinions:

I'm unsure about modifying an existing model (here TS89). I think it would be OK if this extension is very standard. In this case it should at least be properly documented in the manual.
If this extension is not a standard one, but of interest, we could add ad TS89+Boberg to the field list.
Since Dst is in maginput, why not drive this extension directly with Dst? It would require passing it to the T89C routine, but it should not be too complicated.

drflei commented 3 years ago

The TS89 is valid for kp up to 6. In IRBEM, one can use the model for higher kps but the field model is same for all kp>6 cases, i.e., it is the kp=6 field. The motivation for this is to extend the TS89 to the higher kps using the Boberg version of the model, but without impact to the kp<=6 cases. So what's implemented is a hybrid of the TS89 and the Boberg version:

For kp< 7, it is the TS89 model
for kp>=7, it is the Boberg model with fixed Dst(nT) for each kp, kp7:-200, kp8:-225, kp9:-250

This is an extension to the TS89, i.e., it is still a single parameter, kp, model and there is no change in the kp<7 cases.
This extension is widely used but there is no standard or consensus on the kp-dst combinations. For this reason, I think I better withdraw the pull request and not imposing others to follow my choice of the kp/dst combinations.

AntoineBrunet commented 3 years ago

I'll wait for others to pitch in, but my gut feeling is that we should integrate this extension, but drive it using Dst from maginput. This way you can still set the Dst values based on Kp as you are doing, but manually in the maginput parameter. This seems to me like the most versatile option.

VinniDoddu commented 3 years ago

Hello All, I agree with @AntoineBrunet. The already implemented version of Tsyganenko 89 model (TSY89c) is only Kp-dependent and is a reference. Many studies have relied on it, with its well-known limitations. Thus I would also suggest to add a new magnetic field model (in kext) to take into this extension, with the C5(Dst) function implemented in. However, it seems more interesting to let the user choose the Dst value for high Kp values, only adding a kind of disclaimer to take care of the results he/she could obtained. If there is not yet published work on that, maybe your implementation @drflei may help some people analyze carefully its benefits, as Tsy89c is one of the fastest magnetic field model.

mshumko commented 3 years ago

@drflei I second @AntoineBrunet that we should add the extension into a new kext option "TS89+Boberg" and not overwrite the "T89" option.

Once this is implemented on the Fortran side, to make this work on Python, all we have to do it add the option to the extModels lookup list.

drflei commented 3 years ago

I agree it is better to have a separate TS89BOB model rather than the current hybrid approach. It will be the full Boberg modification, i.e., the model requires two inputs, Kp and Dst. We can recommend Dst/Kp combinations for usage as an extension to the TS89C model for higher Kps, but ultimately that's the user's decision.

I am happy to make the necessary changes in onera_desp_lib.f and created a Ts89bob.f from the Tsyganenko89.f, but I am reluctant to made the wider changes required for the integration of a new model. Hope one of you can take the lead for this, starting with the decision on which kext value is assigned to the new TS89BOB model..

PRBEM / IRBEM

Added the Boberg extension to the TS89 model for Kp > 6 #11

Pull Request Checklist