mattkjames7 / PyGeopack

Wrapper for geopack-08 used for the Tsyganenko magnetic field models
GNU General Public License v3.0
11 stars 2 forks source link

T96 model gives different results when compiled with pygeopack vs Geopack #20

Closed EvaKraemer closed 9 months ago

EvaKraemer commented 10 months ago

I have been comparing the pygeopack package with the 'original' Geopack_08 package. Unfortunately, I am not able to reconstruct the same result using the same input parameters in both versions of the T96 model.

In the following I have an example of a set of input parameters that lie inside the magnetopause when compiled with Fortran (double precission), but not with python. I can also provide the Fortran code if needed.

import PyGeopack as gp import numpy as np

model = 'T96' r = [11.46, 0, 0] v_sw = [-400, 10, 0] date = 20160101 ut = 1.0166666666 coord = 'GSE' pdyn_sw = 1.3455 dst = -10 b_sw = [5,5,5]

b_field = gp.ModelField(r[0],r[1],r[2], date,ut,Model=model, CoordIn=coord,CoordOut=coord, WithinMPonly=False, Vx = v_sw[0], Vy = v_sw[1], Vz = v_sw[2], Pdyn = pdyn_sw, SymH = dst, By = b_sw[1], Bz = b_sw[2]) b_field = np.array(b_field).flatten() print('Inside magnetopause: {}'.format(~np.all(np.isnan(b_field))))

Are there some fundamental differences that I am missing?

mattkjames7 commented 10 months ago

Hi EvaKraemer,

Thanks for using my code and bringing this to my attention.

It's possible that I have made some sort of mistake somewhere - would you mind providing me with the field vector that you get with these parameters using Fortran please?

When I run my code like so:

bx,by,bz = gp.ModelField(11.46,0.0,0.0,20160101,1.016666666,Model="T96",CoordIn="GSE",CoordOut="GSE",WithinMPOnly=False,Vx=-400.0,Vy=10.0,Vz=0.0,Pdyn=1.3455,SymH=-10.0,By=5.0,Bz=5.0)

I get the vector [2.12870323, 12.74425398, 32.29218772] if I set WithinMPOnly=False. I get nan when I set WithinMPOnly=True which, if this is supposed to be a position within the MP, suggests that I have made a mistake with my function to check if the position is within the magnetosphere.

If the field vector is wrong, then it could be a mistake somewhere in the coordinate transforms or elsewhere, which would be pretty annoying!

Cheers, Matt.

EvaKraemer commented 9 months ago

Hi again, I noticed that I was not working with the latest release where you fixed the Real*8 in the T96 function in python. I found the same mistake in the fortran code last week and after fixing that and now updating the PyGeopack I get very similar results. So it seems like you did everything correct with the coordinate transformation.

Sorry for that mistake from my side.

Best, Eva

mattkjames7 commented 9 months ago

Oh brilliant, thanks for finding this!