tsssss / geopack

Python version of geopack and Tsyganenko models
MIT License
30 stars 12 forks source link

Difference in output of t96 and t01 model between IDL and Python geopacks! #9

Closed qudsiramiz closed 2 years ago

qudsiramiz commented 2 years ago

Hello, I am testing and comparing the outputs of this geopack for two different models with those of IDL and I am getting very different results.

I have attached the figures which I plotted for this purpose. For each panel in the figure I have plotted the difference between the output from IDL and Python (b_out_idl - b_out_py). Difference from T96 are plotted in the first row whereas the difference between T01 model are plotted in the second row. The value of x_gsm is at the top left position whereas the value of DST used for the code is at top right

For both T96 and T01 model far from the earth, into the tail, the difference between two outputs is minimal, almost always less than 2nT, and thus can be ignored in most cases (Figure 1, for x_gsm = 9.98 R_E).

However, there seem to be significant difference in output when we move close to the Earth or when we are on the dayside. Though the difference is smaller for T01 outputs compared to those of T96, they are still significant, specially far from the Earth on the dayside (Figure 2).

I wonder if anyone else have had similar issues.

The parameters I used for the codes are as follows:

time = 2015-01-01 00:00:00 GMT par = [5, 0, 1, 1, 0, 0, 0, 0, 0, 0] x_gsm = np.linspace(-15.1, 15, 61) y_gsm = np.linspace(-15.1, 15, 61) z_gsm = np.linspace(-15.1, 15, 61)

I can provide the full python and IDL codes if that will help with the reproduction of these differences.

Figure 1 Figure_3

Figure 2 Figure_2

Figure 3 Figure_1

tsssss commented 2 years ago

Hello,

It'll be helpful if you could attach the IDL and Python codes. Also, just want to check if you can replicate the examples in https://github.com/tsssss/geopack/blob/master/README.md.

Thanks, Sheng

On Mon, 18 Oct 2021 at 14:35, Ramiz Qudsi @.***> wrote:

Hello, I am testing and comparing the outputs of this geopack for two different models with those of IDL and I am getting very different results.

I have attached the figures which I plotted for this purpose. For each panel in the figure I have plotted the difference between the output from IDL and Python (b_out_idl - b_out_py). Difference from T96 are plotted in the first row whereas the difference between T01 model are plotted in the second row. The value of x_gsm is at the top left position whereas the value of DST used for the code is at top right

For both T96 and T01 model far from the earth, into the tail, the difference between two outputs is minimal, almost always less than 2nT, and thus can be ignored in most cases (Figure 1, for x_gsm = 9.98 R_E).

However, there seem to be significant difference in output when we move close to the Earth or when we are on the dayside. Though the difference is smaller for T01 outputs compared to those of T96, they are still significant, specially far from the Earth on the dayside (Figure 2).

I wonder if anyone else have had similar issues.

The parameters I used for the codes are as follows:

time = 2015-01-01 00:00:00 GMT par = [5, 0, 1, 1, 0, 0, 0, 0, 0, 0] x_gsm = np.linspace(-15.1, 15, 61) y_gsm = np.linspace(-15.1, 15, 61) z_gsm = np.linspace(-15.1, 15, 61)

I can provide the full python and IDL codes if that will help with the reproduction of these differences.

Figure 1 [image: Figure_3] https://user-images.githubusercontent.com/25327981/137784185-4058c4c8-93fc-42ee-8520-dc452245a665.PNG

Figure 2 [image: Figure_2] https://user-images.githubusercontent.com/25327981/137787058-989c43f6-9cc4-4dd5-95bc-bcf39a657d75.PNG

Figure 3 [image: Figure_1] https://user-images.githubusercontent.com/25327981/137784239-fb801129-daa7-41a8-a991-f3a2824a95c6.PNG

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/tsssss/geopack/issues/9, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADBMLBWEE4HDX5AVDBKBVELUHRSIDANCNFSM5GHJU4CQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

qudsiramiz commented 2 years ago

Hey Sheng, Yes I was able to reproduce the example which are in Readme file. They are exactly same.

I have uploaded the codes in my repository and have linked them here:

idl_code python_code idl_data_edit plotting_code

I have tried to make the code consistent so that you can run them on your system with minimal edits. Also, please note that the codes are listed in the order they are run!

tsssss commented 2 years ago

Hello,

The T96 coef was recently updated but I haven't uploaded the change to pip. This may explain the T96 difference. Have you tried the latest codes at https://github.com/tsssss/geopack/tree/master/geopack?

Below are my tests at 2 selected points (using the updated T96 coef). I think the values from Python are in general agreement with those from IDL.

I noticed your par is different in Python [5,-100,1,1,...] from IDL [5,0,1,1,...] Maybe this is why the T01 results are so different? As shown below for my test Point #2, the T01 difference between Python and IDL is <1 nT. But it reads <-20 nT in Figure 2 bx-t01.

Let me know if this helps.

Thanks, Sheng

Point #1. for point in your IDL code at Line 5

import geopack.geopack as gp

Load IGRF coefficients ...

import datetime

from dateutil import parser

t1 = datetime.datetime(2015,1,1,0,0,0)

t0 = datetime.datetime(1970,1,1)

ut = (t1-t0).total_seconds()

ps = gp.recalc(ut)

xgsm,ygsm,zgsm = [0.95,0.5,5]

par = [5,0,1,1,0,0,0,0,0,0]

bx_t96,by_t96,bz_t96 = gp.t96.t96(par, ps, xgsm,ygsm,zgsm)

bx_t96,by_t96,bz_t96

(16.664905592587036, 0.5368630998346164, -7.766438571515584)

bx_t01,by_t01,bz_t01 = gp.t01.t01(par, ps, xgsm,ygsm,zgsm)

bx_t01,by_t01,bz_t01

(11.747186665975818, 0.9160073121038579, 7.7570418087141615)

bx_igrf,by_igrf,bz_igrf = gp.igrf_gsm(xgsm,ygsm,zgsm)

bx_igrf,by_igrf,bz_igrf

(-198.86095806029613, -58.91994471219549, -329.89745764187376)

IDL> print, [bx_t96, by_t96, bz_t96]

   16.665817      0.53700743      -7.7665478

IDL> print, bx_t01, by_t01, bz_t01

   11.944837      0.77270680       7.3268356

IDL> print, bx_igrf, by_igrf, bz_igrf

  -198.85908      -58.922195      -329.91322

Point #2, for a point in Figure 2

xgsm,ygsm,zgsm = [-10,0,0]

bx_igrf,by_igrf,bz_igrf = gp.igrf_gsm(xgsm,ygsm,zgsm)

bx_igrf,by_igrf,bz_igrf

(24.538720457834742, -0.3195681706868667, 26.965386128130422)

bx_t01,by_t01,bz_t01 = gp.t01.t01(par, ps, xgsm,ygsm,zgsm)

bx_t01,by_t01,bz_t01

(37.939421306182794, 0.19814634751828214, -13.693608165120146)

bx_t96,by_t96,bz_t96 = gp.t96.t96(par, ps, xgsm,ygsm,zgsm)

bx_t96,by_t96,bz_t96

(35.226061425661655, 0.24487109090381676, -15.711727726040502)

IDL> print, [bx_t96, by_t96, bz_t96]

   35.195693      0.24487109      -15.529020

IDL> print, bx_t01, by_t01, bz_t01

   38.210491      0.18678968      -13.792622

IDL> print, bx_igrf, by_igrf, bz_igrf

   24.537546     -0.31953970       26.966250

On Mon, 18 Oct 2021 at 14:58, Ramiz Qudsi @.***> wrote:

Hey Sheng, Yes I was able to reproduce the example which are in Readme file. They are exactly same.

I have uploaded the codes in my repository and have linked them here:

idl_code https://github.com/qudsiramiz/rcn/blob/main/codes/idl_code_gh.pro python_code https://github.com/qudsiramiz/rcn/blob/main/codes/python_code_gh.py idl_data_edit https://github.com/qudsiramiz/rcn/blob/main/codes/write_edited_file.py plotting_code https://github.com/qudsiramiz/rcn/blob/main/codes/idl_ipy_comparison.py

I have tried to make the code consistent so that you can run them on your system with minimal edits. Also, please note that the codes are listed in the order they are run!

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/tsssss/geopack/issues/9#issuecomment-946116782, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADBMLBQ6ALOL4PT6L3DIRBTUHR37TANCNFSM5GHJU4CQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

qudsiramiz commented 2 years ago

Hey Sheng, Sorry for the delay in response. But yes it worked after I updated it with the latest version which is here on GitHub. Thank you so much for your help in debugging this!