European Space Agency

Royal Belgian Institute for Aeronomy

Ex. #7. iso-contour of the shell parameter

Description

Note that this example has been performed using IDL, where the library can be directly linked to IDL (see FAQ G.06, and FAQ T.05), and that the IDL routine unilib is used to correctly initialized the library.

In this example, the subroutine ul220() is called several times to evaluate the magnetic field intensity \(B_m\) and the McIlwain shell parameter L at different geographic positions in the plane that includes the Earth rotation axis and the geomagnetic dipole axis.

The IDL routine CONTOUR is then used to display the results. The geomagnetic field model IGRF, epoch 1995, is used without external magnetic field for this computation. It is selected by a call to the subroutines um510() and um520(). The longitude of the geomagnetic dipole axis is obtained by a call to subroutine ut986().

Be aware that the iso-contours of the shell parameter do not necessarily correspond to magnetic field lines.

Source

;
;  Initialisation of the library
;
   unilib
;
;  Selection of the magnetic field model
;
   ifail   =  0L
   label   = '12345678901234567890123456789012'
   status  = CALL_EXTERNAL( 'unilib', 'UM510', 0L, 1995.0d0, label, 6L, ifail )
   IF ifail LT 0 THEN STOP
   status  = CALL_EXTERNAL( 'unilib', 'UM520', 0L, 0.0d0, DBLARR(10), $
                            label, 6L, ifail )
   IF ifail LT 0 THEN STOP
;
;  Get the geomagnetic dipole axis
;
   msun    = {zxyz}
   gm      = 0.0d0
   clatdip = 0.0d0
   elngdip = 0.0d0
   re      = 0.0d0
   status  = CALL_EXTERNAL( 'unilib', 'UT986', gm, clatdip, elngdip, msun, re )
;
;  Define a grid of points
;
   colat   = [INDGEN(21)*150./20.+15.,165.-INDGEN(21)*150./20.,15.] # (FLTARR(21)+1)
   elong   = [FLTARR(21,21),FLTARR(21,21)+180.,FLTARR(1,21)] + elngdip
   radius  = (FLTARR(43)+1) # (INDGEN(21)/2.+1)
   z       = radius * COS( colat*!dtor )
   x       = radius * SIN( colat*!dtor ) * COS( (elong-elngdip)*!dtor )
   vallm   = FLTARR (43,21)
   valbm   = FLTARR (43,21)
;
;  Intermediate variables
;
   fbm     = 0.d0
   flm     = 0.d0
   fkm     = 0.d0
   fsm     = 0.d0
   fbeq    = 0.d0
   fs      = 0.d0
   mpos    = {zgeo}
;
;  Loop on all the points
;
   FOR i = 0, N_ELEMENTS( radius ) - 1 DO BEGIN
     mpos.radius = radius(i)*re
     mpos.colat  = colat(i)
     mpos.elong  = elong(i)
     ifail       = 0L
     status      = CALL_EXTERNAL( 'unilib', 'UL220', mpos, 90.0d0, 1L, $
                                  fbm, flm, fkm, fsm, fbeq, fs, ifail)
       vallm(i)  = flm
       valbm(i)  = fbm
   ENDFOR
;
;  Plot of the results
;
   range   = 5 * [-1,1]
   rapp    = 1.3
   PLOT, [0], /nodata, xsty=3, ysty=3, xra=range*rapp, yra=range, chars=1.5
;
   levlm   = [ 1.2, 1.5, 2., 3., 4., 6., 8. ]
   levbm   = [ 0.002, 0.005, 0.01, 0.03, 0.1, 0.2, 0.4 ]
;
   i = WHERE( vallm LT 0 )
   IF i(0) GE 0 THEN BEGIN
     OPLOT, x(i), z(i), psym=3, color=!p.color*0.2+!p.background*0.8
     valbm(i) = -1
   ENDIF
;
   CONTOUR, vallm, x, z, min_val=0., /overplot, levels=levlm, $
            c_labels=levlm*0, c_colors=!p.color*0.8+!p.background*0.2
   CONTOUR, valbm, x, z, min_val=0., /overplot, levels=levbm, $
            c_labels=levbm*0, c_colors=!p.color*0.5+!p.background*0.5
;
   END

Results

The green and blue curves correspond to the L and \(B_m\) iso-contours, respectively. The red dots indicate geographic location where the subroutine ul220() failed to compute the value of L.

---

UNILIB/tags/v3.02