;--------------------------------------------------------------------------------------
; 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
; 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)
; 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
; 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
; note: in right-handed system 3 coordinates, 158 degrees longitude is in the
; direction of tilt of the magnetic pole
phi = 158./!radeg
; 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)
rpole2 = Rpole*Rpole
sin2 = sin(theta)*sin(theta)
cos2 = cos(theta)*cos(theta)
rhob = sqrt(sin2 + cos2*rpole2)
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)
ys3 = rhob*sin(theta)*sin(phi)
zs3 = rhob*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]
; 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
; 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
; call the Khurana magnetic field model, the only inputs are
kk_2009, ctime, r, theta, phi, Brm, Btm, Bpm, mpCheck
; compute the magnitude of the field
BnT = 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)
; 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.)
; step the field in cartesian
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)
; recalculate the oblateness at this point's latitude
sin2 = sin(theta)*sin(theta)
cos2 = cos(theta)*cos(theta)
rhob = sqrt(sin2 + cos2*rpole2)
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]
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