trend

- 4

Trapped Radiation ENvironment model Development

Example 7 - iso-contour of the shell parameter

[ Description | Source | Results ]

Description

    Note that this exemple has been performed on a Alpha/VMS plateform where the library can be directly linked to IDL (see FAQ G.06 and T.05), and that the IDL routine unilib_init is used to correctly initialized the library.

    In this exemple, the subroutine UL220 is called several time to evaluate the magnetic field intensity Bm and the McIlwain shell parameter L at different geographic positions in the plane including both 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_init
;
;  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 Bm iso-contours, respectively. The red dots indicate geographic location where the subroutine UL220 failed to compute the value of L.