;-------------------------------------------------------------------------------------- ; 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