irfpy.jupiter.kk_2009_port

Krishan Khurana’s magnbetic field model

It is a line-by-line translate of the IDL code in http://lasp.colorado.edu/mop/resources/code/#khuranamodel

If you publish a paper using the code, it is a good idea to cite a paper of Khurana & Schwarzl (2005).

Ported to python by Yoshifumi Futaana, 2014.

Warning

It is always recommended to use irfpy.jupiter.kk2009. This irfpy.juptier.kk_2009_port contains a bug, The keep of this module is kept for the compatibility with the original IDL code.

USAGE:

The kk_2009() returns the magnetic field in the system III coordinates.

>>> results = kk_2009(ctimer(2000, 6, 1, 10, 32, 17), 30.00, np.deg2rad(86.88), np.deg2rad(157.94))
>>> print('{r[0]:.4f} {r[1]:.4f} {r[2]:.4f}'.format(r=results[0]))
28.9079 5.0641 -3.5039
>>> print('{:.4f}'.format(results[1]))
-4.1292
>>> print(results[2])
1

Warning

Do not use “raw” ctime value. Call always ctimer() function where global variables are set for, likely, leapsecond call… YF does not like this feature by obvious reasons, but it is how the original code is implemented.

Warning

It looks that the original IDL code has a serious error. Try the following IDL code in a single session.

.compile kk_2009

; Conversion of coordinate is tried
ctime = ctimer(1994, 8, 21, 3, 0, 0)
vecin = [1., 0., 0.]
vecou = [0., 0., 0.]
jrot, 's3c', 'jss', vecin, vecou, ctime
print, 'First answer', vecou
; First answer     -0.87623097     -0.48189136       0.0000000

; Here, different time is examined.
ctime = ctimer(2000, 6, 1, 10, 32, 17.0)
kk_2009, ctime, 30.0, 86.88 * 3.1416/ 180., 157.94 * 3.1416/ 180., Brm, Btm, Bpm, ianswer

; Conversion of coordinate exactly same as the first one tried.
ctime = ctimer(1994, 8, 21, 3, 0, 0)
vecin = [1., 0., 0.]
vecou = [0., 0., 0.]
jrot, 's3c', 'jss', vecin, vecou, ctime
print, 'Second answer', vecou
; Second answer     -0.92678450     -0.37559352       0.0000000

; The answers differ!!

The module, kk_2009_port is kept the bug as it has been. This is rather to keep the compatibility (for mainly the purpose of debugging/validating the translation) to IDL code.

On the other hand, YF fixed the bug for the module, kk_2009. I eliminate using global variables (system variables in IDL original code). This looks to have been avoided. (Source code lines around 1030).

In practice, using kk_2009 is highly recommended. If you want “bug” version, you simply name kk_2009_port to kk_2009:

> from irfpy.jupiter import kk_2009_port as kk_2009

irfpy.jupiter.kk_2009_port.print_global()

irfpy.jupiter.kk_2009_port.kk_2009(ctime, R, theta, phi)

KK model

Parameters:

• ctime – J2000 time in seconds (can by calculated by the ctimer subroutine located in this code)

• R –

• theta –

• phi – Right-handed system III coordinates. Theta phi in radians.

Returns:

A tuple, ((Br, Bt, Bp), ZNS3, mpCheck). (Br, Bt, Bp) in System III spherical coordinats, ZNS3 for the distance to the current sheet from the jovigraphic equator in system III coordinates, and mpCheck if in magnetopause(=1) or not(=0).

Usage:

First you have to calculate the “ctime”, J2000 time in second.

>>> ctime = ctimer(2000, 6, 1, 10, 32, 17.0)
>>> print(ctime)
1086085937.0

Then, call this function

>>> R, theta, phi = 30.00, np.deg2rad(86.88), np.deg2rad(157.94)
>>> Bvec, icheck = kk_2009(ctime, R, theta, phi)   

Note

Original comments:

####    ;--------------------------------------------------------------------------------------
####    ;---------------Krishan Khurana's Jupiter magnetic field model for IDL-----------------
####    ;--------------------------------------------------------------------------------------
####    ; NAME:
####    ;   kk_2009
####    ;
####    ; PURPOSE:
####    ;   Given right-handed system III spherical coordinates, this program
####    ;    outputs the strength of jupiter's magnetic field in spherical coordinates
####    ;    at the specified coordinates. 
####    ;
####    ; DESCRIPTION:
####    ;   Calculates the magnetic field of a tilted and warped shielded magnetosphere
####    ;    using general deformation technique.
####    ;
####    ; CALLING SEQUENCE:
####    ;   kk_2009, ctime, R, theta, phi, Brm, Btm, Bpm, mpCheck
####    ;
####    ; INPUTS:
####    ;   ctime         - J2000 time in seconds (can by calculated by the
####    ;                   ctimer subroutine located in this code)
####    ;   R, theta, phi - right-handed system III coordinates
####    ;
####    ; OUTPUTS:
####    ;   Brm, Btm, Bpm    - r, theta, phi components of the magnetic field
####    ;                      in nanoteslas [nT]
####    ;   mpCheck          - equals 1 if  inside the magnetopause
####    ;                    - equals 0 if outside the magnetopause
####    ;  
####    ;   note: as a procedure the output to kk_2009 can easily be changed
####    ;         to include more than just the components to the magnetic field
####    ;         - for example, the following can also be included:
####    ;
####    ;   XJSO, YJSO, ZJSO - x, y, z components of the magnetic field in the 
####    ;
####    ;                      juptier-sun-orbital coordinate system (see JROT subroutine),
####    ;   ZNS3             - distance to the current sheet from the
####    ;                      jovigraphic equator in system III coordinates
####    ;
####    ;
####    ; MODIFICATION HISTORY:
####    ;   2003, written by Krishan Khurana in Fortran 77/90
####    ;   2008, ctimer function added by Mariel Desroche
####    ;         (converts year, month, day, hour and second to J2000 seconds)
####    ;   2009, converted to IDL by Adam Shinn
####    ;   2009, leapSecond function added by Adam Shinn and Rob Wilson
####    ;         (replaces eetime function to easily append for future leap seconds)
####    ;
####    ; WEBSITE:
####    ;   This code is available on the website for the Magnetospheres of
####    ;    Outer Planets group at the Univeristy of Colorado at Boulder,
####    ;    associated with the Laboratory for Atmospheric and Space Physics:
####    ;    http://lasp.colorado.edu/MOP/resources/
####    ;--------------------------------------------------------------------------------------

irfpy.jupiter.kk_2009_port.doy(iyr, imon, iday)

DOY calculator

I like simpler calculation using datetime, but keep it as it has been…

>>> print(doy(1996, 3, 1))
61
>>> print(doy(1999, 3, 1))
60
>>> print(doy(2000, 3, 1))
61
>>> print(doy(2100, 3, 1))
60

irfpy.jupiter.kk_2009_port.ctimer(iyr, imon, iday, ihr, imin, sec)

Calculate J2000 time.

If YF understand correctly, this is the epoch from 1966-01-01T00:00:00. In this case, we can make simpler function using dateutil or time modules.

>>> print(ctimer(1966, 1, 1, 0, 0, 0))
0.0
>>> print(ctimer(1966, 1, 1, 0, 0, 1))
1.0
>>> print(ctimer(1967, 1, 1, 0, 0, 0))
31536000.0

The following test is in sample.pro.

>>> print(ctimer(2000, 6, 1, 10, 32, 17))
1086085937.0

irfpy.jupiter.kk_2009_port.leapSecond(ctime)

If leap second identified, consider it.

Indeed, the implementation uses global variables. Very quick way.

What is ctime here? In the end, the returned is ctime - 34year. ctime is 2000-01-01?

Todo

Look for definition of J2000 epoch!

> kk.leapSecond(kk.ctimer(2000, 1, 1, 0, 0, 0))
-43135.81599998474

irfpy.jupiter.kk_2009_port.JSun(ctime)

(To be written)

Parameters:

ctime – J2000 time in seconds (can by calculated by the ctimer subroutine located in this code)

Returns (stheta, sphi, phase):

stheta, sphi is latitude and longitude (in radians) of the Sun in system III (RH). phase is Orbital phase angle of Jupiter (in radians)

Todo

Carefully address the difference of IDL’s mod and python’s %. For IDL, neg mod pos can return neg, but not for python. (Really?)

irfpy.jupiter.kk_2009_port.jrot(from_, to_, vecin, ctime)

Convert by rotation.

Parameters:

• from – Incoming frame

• to – Outgoing frame

• vecin – Vector in the incoming frame.

• ctime – a double precision variable denoting Cline time which can be calculated by using the function program ctimer.

The “from” coordinate system is rotated to system III, then rotated into the “to” coordinate system.

The supported coordinate systems are:

• S3C System III Cartesian (right-handed)

• JSO Jupiter-Sun-Orbital

• JSM Jupite-Sun-Magnetic

• DIP Dipole (cartesian)

• JSS Jupiter-Sun-Spin

irfpy.jupiter.kk_2009_port.cyl2car_pos(rho, phi, Z)

irfpy.jupiter.kk_2009_port.car2cyl_vect(BXcar, BYcar, BZcar, phi)

irfpy.jupiter.kk_2009_port.cyl2car_vect(Brho, Bphi, BZ, phi)

irfpy.jupiter.kk_2009_port.CAR2SPH_pos(X, Y, Z)

irfpy.jupiter.kk_2009_port.SPH2CAR_pos(R, TH, PHI)

irfpy.jupiter.kk_2009_port.CAR2SPH_MAG(BX, BY, BZ, TH, PHI)

irfpy.jupiter.kk_2009_port.PJSO2S3(BXPJSO, BYPJSO, BZPJSO, phi)

irfpy.jupiter.kk_2009_port.S32PJSO(BXS3, BYS3, BZS3, phi)

irfpy.jupiter.kk_2009_port.bessj1(x)

Bessel function.

Returns bessel function with real order of 1.

>>> print('{:.6f}'.format(bessj1(0)))
0.000000
>>> print('{:.6f}'.format(bessj1(0.1)))
0.049938
>>> print('{:.6f}'.format(bessj1(2)))
0.576725
>>> print('{:.6f}'.format(bessj1(-0.1)))
-0.049938

Internally, scipy:scipy.special.jv() is used. The returned values agreed with the original IDL version.

irfpy.jupiter.kk_2009_port.bessj0(x)

bessel function with real order of 0.

>>> print('{:.6f}'.format(bessj0(0)))
1.000000
>>> print('{:.6f}'.format(bessj0(0.1)))
0.997502
>>> print('{:.6f}'.format(bessj0(2)))
0.223891
>>> print('{:.6f}'.format(bessj0(-0.1)))
0.997502

Internally, scipy:scipy.special.jv() is used. The returned values agreed with the original IDL version.

irfpy.jupiter.kk_2009_port.DBSJ2(x)

irfpy.jupiter.kk_2009_port.DBSJ3(x)

irfpy.jupiter.kk_2009_port.B_mp_perp(rho, phi, x, Nmodes)

irfpy.jupiter.kk_2009_port.dipole(B0x, B0y, B0z, x, y, z)

Dipole

Original comments in IDL
; Calculates the field of Jupiter's dipole for shielding in the magnetopause
;
;   (B0x, B0y, B0z) is the dipole moment
;   x, y, z is the position vector
;   Bx, By, Bz is the output field vector

irfpy.jupiter.kk_2009_port.dipole_shielded(PARMOD, x, y, z)

Obtain dipole field?

Original comments in IDL.

; WRITTEN BY K. K. KHURANA   11/2002
; MODIFIED BY H. K. SCHWARZL 11/2003
; PARMOD is an input array (real*8) that contains the model parameters
;   Dimension of PARMOD is: DIMENSION PARMOD(10)
;   currently just PARMOD[0] is used for the dipole tilt angle
; x,y,z input position
; OUTPUT: mag .filed Bx, By, Bz at x, y, z

irfpy.jupiter.kk_2009_port.dipole_shield_cyl_S3(PARMOD, ctime, rmap, pmap, zmap, sphiOut)

irfpy.jupiter.kk_2009_port.B_tail_shield(x, y, z, M, ModeIn)

irfpy.jupiter.kk_2009_port.tail_mag_shield_cyl_S3(ctime, rmap, pmap, zmap, M, sphiOut, Mode)

irfpy.jupiter.kk_2009_port.csheet_struc(rho, phi, XJSO, YJSO, stheta, ctime)

This program calculates the distance of current sheet from the Jovigraphic equator in system III coordinates

>>> zns3 = csheet_struc(30.0, 0.5, 28.0, 1.2, np.deg2rad(20.), ctimer(2000, 1, 1, 0, 0, 0))
>>> print('%.5f' % zns3)
2.54079

The returned value is examined by IDL.

IDL> .compile kk_2009
IDL> csheet_struc, retval, 30.0, 0.5, 28.0, 1.2, 20. * 3.1416 / 180., ctimer(2000, 1, 1, 0, 0, 0) 
IDL> print, retval
2.5407947

irfpy.jupiter.kk_2009_port.csheet_N(rho, phi, zs3, stheta, ctime)

Something related to current sheet…

Parameters:

• rho – Distance from Jupiter (?)

• zs3 – z-value in system III.

• phi – Phase angle (longitude) in system III, radians.

• stheta – Angle (radian) in system III(?)

• ctime – Time in J2000 epoch (?), where anyway ctimer() derives.

Note

In the source code,

\[xs3 = \rho * \cos(\phi)\]

is done. So that phi is the angle in the same coordinate system to xs3 (and likely zs3).

Later, the rho, phi based vector (xs3, ys3, zs3) is converted from System3 to JSO. Therefore, they looks in SYS3 coordinate.

>>> ctime = ctimer(2000, 6, 1, 10, 32, 17.0)
>>> print(ctime)
1086085937.0

>>> rnx, rny, rnz = csheet_N(29.955566, 2.7565556, 1.6321992, 0.05433560, ctime)
>>> print('{:.6f} {:.6f} {:.6f}'.format(rnx, rny, rnz))
-0.051022 0.055858 0.997134

Note

In IDL original program, the results will be (-0.050987206 0.055844548 0.99713675) from the same input.

The difference (only 3-places accuracy) could come from the fact that the calculation in IDL was done in single-precision, where 7-places accuracy is there, but the part of numerical differentiation, similar values are subtracted and only 3-places accuray is left (cancellation of significant digits).

irfpy.jupiter.kk_2009_port.get_mapped_sunangle(sthetaIN, sphiIN, xpx, xpy, xpz, ypx, ypy, ypz, zpx, zpy, zpz)

irfpy.jupiter.kk_2009_port.getBimf(ctime, thetaS3, phiS3)

irfpy.jupiter.kk_2009_port.CheckIfInsideMappedMp(ctime, XS3, YS3, ZS3, ZNS3)

irfpy.jupiter.kk_2009_port.checkIfInsideMagnetop(x, y, z)

irfpy.jupiter.kk_2009_port.tail_mag_notilt_all_modes(rho, phi, zz, RLT, Mode)

According to the IDL code, the unit test is generated.

>>> rho = 29.774504
>>> phi = -2.6928876
>>> zz = 5.7614139
>>> RLT = 2.7842312
>>> Mode = 7
>>> Brcs, Bpcs, Bzcs = tail_mag_notilt_all_modes(rho, phi, zz, RLT, Mode)
>>> print('%.6f' % Brcs)
20.816283
>>> print('%.6f' % Bpcs)
-3.551422
>>> print('%.6f' % Bzcs)
15.714019

irfpy.jupiter.kk_2009_port.JOVIAN_VIP4_no_dipole(NM, r, theta, phi)

irfpy.jupiter.kk_2009_port.doctests()