snippet_jupiter.sample_kkmodelΒΆ
A port of KK model sample.pro
''' A port of KK model sample.pro
'''
import numpy as np
from irfpy.jupiter import kk_2009
def main():
# ;--------------------------------------------------------------------------------------
# ; NAME:
# ; sample
# ;
# ; PURPOSE:
# ; This procedure shows how to use the IDL Khurana Jupiter magnetic
# ; field model in the simplest way.
# ;
# ; DESCRIPTION:
# ; We start out by defining a time, and a point in space to feed into
# ; the magnetic field model to obtain an output of magnetic field
# ; values and finally, compute the magnitude of the magnetic field at
# ; the point given.
# ; Because the position coordinates that are fed into the model
# ; need to be in spherical coordinates, there will be an example of
# ; coordinate transformation from cartesian to spherical, as well as
# ; an example of switching from one coordinate system to another
# ; (see the JROT subroutine for more information, as well as a description
# ; of supported jovian coordinate systems).
# ;
# ; INPUTS (to the kk_2009.pro field model):
# ; ctime - J2000 time in seconds (can by calculated by the
# ; ctimer subroutine located at the bottom of this code)
# ; R, theta, phi - right-handed system III coordinates
# ;
# ; OUTPUTS (to the kk_2009.pro field model):
# ; Brm, Btm, Bpm - r, theta, phi components of the magnetic field
# ; mpCheck - equals 1 if inside the magnetopause
# ; - equals 0 if outside the magnetopause
# ;
# ; MODIFICATION HISTORY:
# ; 2010, written by Adam Shinn
# ;
# ; WEBSITE:
# ; This code is available on the website for the Magnetospheres of
# ; Outer Planets group for the Univeristy of Colorado at Boulder,
# ; associated with the Laboratory for Atmospheric and Space Physics.
# ; http://lasp.colorado.edu/MOP/resources/
# ;--------------------------------------------------------------------------------------
# pro sample
# ; define a point in time to feed into the magnetic field model
# ; define year, month, day, hour, minute and second
# ; note: iyr, imon, iday, ihr, and imin are integers while sec is a
# ; floating point variable
# ; note: default time is 1-JUN-2000 10:32:17.0 (dipole is pointing
# ; in the direction of noon)
# iyr = 2000
iyr = 2000
# imon = 6
imon = 6
# iday = 1
iday = 1
# ihr = 10
ihr = 10
# imin = 32
imin = 32
# sec = 17.0
sec = 17.0
#
# ; define the output variable of ctimer to be a double precision number
# ctime = 0d
#
# ; call the ctimer function (located at the bottom of kk_2009.pro)
# ; the output is a double precision number that stands for J2000 time in seconds
# ctime = ctimer(iyr, imon, iday, ihr, imin, sec)
ctime = kk_2009.ctimer(iyr, imon, iday, ihr, imin, sec)
print(iyr, imon, iday, ihr, imin, sec)
print(ctime)
#
# ; define a position in space to feed into the magnetic field model
# ; the input coordinates to the model need to be in the form:
# ; R, theta, phi in a right-handed system 3 coordinate system
# ; aside: Jupiter right-handed system 3 (S3) coordinate system
# ; is defined with Jupiter's rotation axis as the azimuthal Z-axis,
# ; the X-axis and Y-axis form a plane normal to the Z-axis,
# ; and the tilt of the northern magnetic "dipole" is pointed
# ; in the direction of 158 degrees longitude
# ; (see the JROT subroutine in kk_2009.pro for more information)
# ; in the jupiter-sun-orbital coordinate system, the longitude phi = 0
# ; refers to the direction of the sun (the direction of noon)
# ; so let's define the point in space to evaluate to be at coordinates(default):
# ; R = 30 [jupiter radii]
# ; theta = 90 [degrees] (equatorial region)
# ; note: theta = 0 [degrees] is defined to be at the north pole
# ; phi = 0 [degrees] (sunward in JSO coordinates)
# RJSO = 30. ; [juipter radii]
# thetaJSO = 90./!radeg ; [radians] the IDL system variable !radeg converts between radians and degrees
# phiJSO = 0. ; [radians]
RJSO = 30. # [juipter radii]
thetaJSO = np.deg2rad(90.) # [radians]
phiJSO = 0. # [radians]
#
# ; since the input coorinates need to be in S3 coordinates, we must transform
# ; the JSO coordinates into cartesian, then use the JROT subroutine to convert
# ; to S3 cartesian coordinates, then transform back into spherical coordinates
# xJSO = RJSO*sin(thetaJSO)*cos(phiJSO)
# yJSO = RJSO*sin(thetaJSO)*sin(phiJSO)
# zJSO = RJSO*cos(thetaJSO)
xJSO = RJSO*np.sin(thetaJSO)*np.cos(phiJSO)
yJSO = RJSO*np.sin(thetaJSO)*np.sin(phiJSO)
zJSO = RJSO*np.cos(thetaJSO)
print(xJSO, yJSO, zJSO)
#
# ; the JROT subroutine takes the form of:
# ; jrot, "input coordinates", "output coordinates", [xin, yin, zin], vecout, ctime
# ; (vecout is a 3 element cartesian vector just like the input vector)
# ; see the JROT subroutine in kk_2009.pro for more information
# jrot, 'jso', 's3c', [xJSO, yJSO, zJSO], vecout, ctime
# xS3 = vecout[0]
# yS3 = vecout[1]
# zS3 = vecout[2]
vecout = kk_2009.jrot('jso', 's3c', np.array([xJSO, yJSO, zJSO]), ctime)
xS3 = vecout[0]
yS3 = vecout[1]
zS3 = vecout[2]
#
# ; convert from cartesian to spherical
# R = sqrt(xS3*xS3 + yS3*yS3 + zS3*zS3)
# theta = atan(sqrt(xS3*xS3 + yS3*yS3), zS3)
# phi = atan(yS3, xS3)
R = np.sqrt(xS3*xS3 + yS3*yS3 + zS3*zS3)
theta = np.arctan2(np.sqrt(xS3*xS3 + yS3*yS3), zS3)
phi = np.arctan2(yS3, xS3)
print(R, theta, phi)
#
# ; print input
# print, '--------------------model input--------------------'
# print, "time :", imon, "/", iday, "/", iyr, format = '(A8, I2, A1, I2, A1, I4)'
# print, "position in right-handed S3 spherical coordinates:"
# print, "r :", R, " [jupiter radii]", format = '(A8, F10.2, A16)'
# print, "theta :", theta*!radeg, " [degrees]", format = '(A8, F10.2, A10)'
# print, "phi :", phi*!radeg, " [degrees]", format = '(A8, F10.2, A10)'
#
# ; call the khurana code, supplying time in J2000 seconds, and
# ; position in spherical right-handed system 3 coordinates
# kk_2009, ctime, R, theta, phi, Brm, Btm, Bpm, ianswer
(Brm, Btm, Bpm), zns3, ipcheck = kk_2009.kk_2009(ctime, R, theta, phi)
# ; print output
# print, '-------------------model output--------------------'
# print, "magnetic field strength in nanotesla"
# print, "Br :", brm, " [nT]", format = '(A9, F10.2, A5)'
# print, "Btheta :", btm, " [nT]", format = '(A9, F10.2, A5)'
# print, "Bphi :", bpm, " [nT]", format = '(A9, F10.2, A5)'
# print, "magnetic field magnitude"
# print, " ", sqrt(brm*brm + btm*btm + bpm*bpm), " [nT]", format = '(A9, F10.2, A5)'
#
# ; for practical uses of the magnetic field model:
# ; use a for loop to step through time to observe the time variation
# ; of jupiter's magnetic field
# ; use two for loops to step through x and y to create a grid of
# ; field values for contour plots of latitudinal planes of the space
# ; within the magnetopause
# ; use a whle loop and clever spacial stepping to trace magnetic fieldlines
# ; (see the tracing sample also available on our group website)
# end ; sample.pro
#### Output is
#--------------------model input--------------------
# time : 6/ 1/2000
#position in right-handed S3 spherical coordinates:
# r : 30.00 [jupiter radii]
# theta : 86.88 [degrees]
# phi : 157.94 [degrees]
#-------------------model output--------------------
#magnetic field strength in nanotesla
# Br : 28.91 [nT]
# Btheta : 5.06 [nT]
# Bphi : -3.50 [nT]
#magnetic field magnitude
# 29.56 [nT]
print (Brm, Btm, Bpm)
if __name__ == "__main__":
main()