snippet_jupiter.sample_traceΒΆ
A sample for the tracing (ported from IDL code)
''' A sample for the tracing (ported from IDL code)
'''
import numpy as np
from irfpy.jupiter.kk_2009 import * # For usual-use purpose
## from irfpy.jupiter.kk_2009_port import * # For debug/port purpose
def trace():
#### ;--------------------------------------------------------------------------------------
#### ; NAME:
#### ; sampleTrace
#### ;
#### ; PURPOSE:
#### ; This procedure shows how to use the IDL Khurana Jupiter magnetic
#### ; field model to trace field lines and draw them in IDL's iplot tool.
#### ;
#### ; DESCRIPTION:
#### ; At a given time, the trace procedure starts at a given latitude and longitude
#### ; on the planet and uses steps scaled by the magnetic field magnitude to trace along
#### ; fieldlines that eventually end up back at the planet. The fieldlines are traced in
#### ; a while loop, and then plotted in IDL's iplot program. In addition to the fieldlines,
#### ; a wire-frame oblate spheroid is also plotted to compare the size and orientation of
#### ; the fieldlines with respect to Jupiter.
#### ;
#### ; 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 trace
#### ; 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
#### imon = 6
#### iday = 1
#### ihr = 10
#### imin = 32
#### sec = 17.0
iyr = 2000
imon = 6
iday = 1
ihr = 10
imin = 32
sec = 17.0
####
#### ; define the output variable of ctimer to be a double precision number
#### ctime = 0d
ctime = 0.
####
#### ; 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 = ctimer(iyr, imon, iday, ihr, imin, sec)
####
#### ; instead of defining a fixed latitude, we can use a for loop to step through
#### ; several latitudes, and therefore trace several lines
#### for itheta = 130, 155, 5 do begin
for itheta in range(130, 156, 5):
print('#THETA: ', itheta)
#### ; note: field line tracing from Jupiter's southern hemisphere is usually easier
#### ; because there is a magnetic anomaly on the northern hemisphere
#### ; convert theta from degrees to radians
#### theta = itheta/!radeg ; system variable !radeg converts between radians and degrees
theta = np.deg2rad(itheta)
#### ; note: in right-handed system 3 coordinates, 158 degrees longitude is in the
#### ; direction of tilt of the magnetic pole
#### phi = 158./!radeg
phi = np.deg2rad(158.)
####
#### ; now to define the starting point for the radial coordinate -
#### ; if we start from the planet's surface (Jupiter's surface is
#### ; defined at the pressure of 1 bar) then we need to account
#### ; for the oblateness of Jupiter which is 1/15.4 less at the
#### ; poles of the planet than at the equator
#### rpole = (1. - 1./15.4)
rpole = (1. - 1./15.4)
#### rpole2 = Rpole*Rpole
rpole2 = rpole*rpole
#### sin2 = sin(theta)*sin(theta)
sin2 = np.sin(theta)**2
#### cos2 = cos(theta)*cos(theta)
cos2 = np.cos(theta)**2
#### rhob = sqrt(sin2 + cos2*rpole2)
rhob = np.sqrt(sin2 + cos2*rpole2)
#### r = rhob
r = rhob
####
#### ; transform the starting points from spherical to cartesian
#### ; because we eventually trace the line in cartesian coordinates
#### xs3 = rhob*sin(theta)*cos(phi)
xs3 = rhob*np.sin(theta)*np.cos(phi)
#### ys3 = rhob*sin(theta)*sin(phi)
ys3 = rhob*np.sin(theta)*np.sin(phi)
#### zs3 = rhob*cos(theta)
zs3 = rhob*np.cos(theta)
####
#### ; these points start off the x, y, and z arrays which are used
#### ; to store the points on the fieldlines
#### x = [xs3]
#### y = [ys3]
#### z = [zs3]
x = [xs3]
y = [ys3]
z = [zs3]
####
#### ; set the magnetopause check to 1 at the start of tracing each line -
#### ; this makes it so the while loop starts if the last line hit the magnetopuase (mpCheck = 0)
#### mpCheck = 1
mpCheck = 1
####
#### ; the while loop makes it so the line traces until it hits the planet (r < rhob)
#### ; or until the lines traces up to the magnetopause (mpCheck = 0)
#### while (r ge rhob) && (mpCheck eq 1) do begin
while (r >= rhob) and (mpCheck == 1):
#### ; call the Khurana magnetic field model, the only inputs are
#### kk_2009, ctime, r, theta, phi, Brm, Btm, Bpm, mpCheck
(Brm, Btm, Bpm), zns3, mpCheck = kk_2009(ctime, r, theta, phi)
####
#### ; compute the magnitude of the field
#### BnT = sqrt(Brm*Brm + Btm*Btm + Bpm*Bpm)
BnT = np.sqrt(Brm*Brm + Btm*Btm + Bpm*Bpm)
####
#### ; transform the spherical magnetic field values into cartesian
#### Bx = Brm*sin(theta)*cos(phi) + Btm*cos(theta)*cos(phi) - Bpm*sin(phi)
#### By = Brm*sin(theta)*sin(phi) + Btm*cos(theta)*sin(phi) + Bpm*cos(phi)
#### Bz = Brm*cos(theta) - Btm*sin(theta)
Bx = Brm*np.sin(theta)*np.cos(phi) + Btm*np.cos(theta)*np.cos(phi) - Bpm*np.sin(phi)
By = Brm*np.sin(theta)*np.sin(phi) + Btm*np.cos(theta)*np.sin(phi) + Bpm*np.cos(phi)
Bz = Brm*np.cos(theta) - Btm*np.sin(theta)
print('#RESULT: ', r, theta, phi, Bx, By, Bz)
####
#### ; ds controls the step size (it is porportional to the field magnitude)
#### ; note: a negative step means you are tracing from the south to the north
#### ds = -.05/alog10(BnT + 2.)
#### delx = ds*Bx/(BnT + 2.)
#### dely = ds*By/(BnT + 2.)
#### delz = ds*Bz/(BnT + 2.)
ds = -.05/np.log10(BnT + 2.)
delx = ds*Bx/(BnT + 2.)
dely = ds*By/(BnT + 2.)
delz = ds*Bz/(BnT + 2.)
####
#### ; step the field in cartesian
#### x_new = xs3 + delx
#### y_new = ys3 + dely
#### z_new = zs3 + delz
x_new = xs3 + delx
y_new = ys3 + dely
z_new = zs3 + delz
####
#### ; transform from cartesian to spherical coordinates to feed back into kk_2009
#### R = sqrt(x_new*x_new + y_new*y_new + z_new*z_new)
#### theta = atan(sqrt(x_new*x_new + y_new*y_new), z_new)
#### phi = atan(y_new, x_new)
r = np.sqrt(x_new*x_new + y_new*y_new + z_new*z_new)
theta = np.arctan2(np.sqrt(x_new*x_new + y_new*y_new), z_new)
phi = np.arctan2(y_new, x_new)
####
#### ; recalculate the oblateness at this point's latitude
#### sin2 = sin(theta)*sin(theta)
#### cos2 = cos(theta)*cos(theta)
#### rhob = sqrt(sin2 + cos2*rpole2)
sin2 = np.sin(theta)*np.sin(theta)
cos2 = np.cos(theta)*np.cos(theta)
rhob = np.sqrt(sin2 + cos2*rpole2)
####
#### xs3 = x_new
#### ys3 = y_new
#### zs3 = z_new
xs3 = x_new
ys3 = y_new
zs3 = z_new
####
#### ; store the values into the arrays
#### ; note: this is a good way to store the values of the fieldline without having
#### ; to define the size of the array beforehand, you're just placing the new
#### ; values at the end of the array
#### x = [x, x_new]
#### y = [y, y_new]
#### z = [z, z_new]
x.append(x_new)
y.append(y_new)
z.append(z_new)
#### endwhile
#### ; iplot the x, y, and z arrays to create a 3D plot that you can visualize the fieldlines with
#### ; you can just as easily use the normal plot procedure, or even plot3D - but iplot provides
#### ; a quick and easy 3D plot visualization that you can rotate and move around
#### ; see iplot in the IDL help manual for keyword help
#### iplot, x, y, z, /overplot, /scale_isotropic, thick = 2.
#### endfor
#### ;now that the fieldlines are plotted, the spheroid program overplots an oblate spheroid which
#### ; represents Jupiter so that you can compare size of the fieldlines to the planet
#### spheroid
####
#### ; to explore on your own:
#### ; step through the parameters of theta or phi using for loops
#### ; to create more fieldlines
#### ; step through time to see how the fieldlines vary with time
#### ; (but keep in mind, the planet rotates with time as well!)
#### end
####
#### ; the spheroid procedure overplots a wire-frame oblate spheroid
#### ; note: 1. corresponds to 1. Juptier radius
#### pro spheroid
#### theta = (180.)*(findgen(101)/100.)/!radeg
####
#### sin2 = sin(theta)*sin(theta)
#### cos2 = cos(theta)*cos(theta)
#### ; 0.874346375465 is the square of the distance to the pole = (1. - 1./15.4)^2.
#### r = sqrt(sin2 + cos2*0.874346375465)
####
#### ; draw longitudinal lines every 10 degrees
#### for iphi = 0, 350, 10 do begin
#### phi = iphi/!radeg
#### x = r*sin(theta)*cos(phi)
#### y = r*sin(theta)*sin(phi)
#### z = r*cos(theta)
#### iplot, x, y, z, /overplot, /scale_isotropic, color = [255, 190, 72], thick = 1.
#### endfor
####
#### ; draw latitudinal lines every 15 degrees
#### phi = (360.)*(findgen(101)/100.)/!radeg
#### for itheta = 45, 135, 15 do begin
#### theta = itheta/!radeg
#### sin2 = sin(theta)*sin(theta)
#### cos2 = cos(theta)*cos(theta)
#### ; 0.874346375465 is the square of the distance to the pole = (1. - 1./15.4)^2.
#### r = sqrt(sin2 + cos2*0.874346375465)
#### x = r*sin(theta)*cos(phi)
#### y = r*sin(theta)*sin(phi)
#### z = replicate(r*cos(theta), 101)
#### iplot, x, y, z, /overplot, /scale_isotropic, color = [240, 155, 0], thick = 1.
#### endfor
#### end
def main():
trace()
if __name__ == '__main__':
main()