;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; ; *** WARNING *** These models have been initially verified but require ; more testing. Use with caution!!! ; ; Notes: ; 1. The THM model returns the closest results at larger distances and ; the LM model - for smaller distances. ; 2. While the model scripts work, there is no such thing as the best ; or most accurate model. ; ;------------------------------------------------------------------------------ ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;NAME: ; sm_ns_model ; ;PURPOSE: ; This routine calculates the position along the zaxis at a specific ; x and y location. The distance to the neutral sheet is the z ; component of the position data where gsm has been converted to sm, ; dz2ns=sm[*,2]. ; ;CALLING SEQUENCE: ; zthemis=neutral_sheet(time, position) ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in GSM coordinates in re (pos[*,3]) ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;NOTES: ; For the nominal mission, THEMIS used this model for the inner probes ; ;HISTORY: ; ;----------------------------------------------------------------------------- FUNCTION sm_ns_model, time, gsm_pos, z2ns=z2ns ; convert gsm to sm coordinates cotrans,gsm_pos,sm_pos,time,/GSM2SM z2ns = -sm_pos[*,2] dz2ns = z2ns RETURN, dz2ns END ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;NAME: ; themis_NS_model ; ;PURPOSE: ; This routine calculates the position along the zaxis at a specific ; x and y location. The themis model is used for this calculation. ; The themis model uses z-sm (converted from z-gsm) for the inner probes ; and the Hammond model for the outer probes. ; ;CALLING SEQUENCE: ; dz2ns=neutral_sheet(time, position) ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in GSM coordinates in re (pos[*,3]) ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;NOTES; ; Reference: ; The themis model uses z-sm (converted from z-gsm) for the inner probes ; and the Hammond model (default) for the outer probes. The algorithm can be found ; in ssllib neutralsheet.pro. ; ;HISTORY: ; ;----------------------------------------------------------------------------- FUNCTION themis_ns_model, time, gsm_pos, z2ns=z2ns ; initialize constants and variables re = 6378. h0 = 8.6 ;10.5 ; hinge point of the neutral sheet rad = !pi/180. z2NS = make_array(n_elements(time), /double) dz2NS = make_array(n_elements(time), /double) ; constants used in hammond model H1=8.6 Y0=20.2 D=12.2 ; calculate the radial distance rdist = sqrt(gsm_pos[*,0]^2 + gsm_pos[*,1]^2 + gsm_pos[*,2]^2) ; Use the sm coordinates for radial distances <= h0 (8.6) sm_ind = where(rdist LE h0, ncnt) IF ncnt GT 0 THEN BEGIN cotrans,gsm_pos[sm_ind,*],sm_pos,time[sm_ind],/GSM2SM z2ns[sm_ind] = sm_pos[*,2] dz2ns[sm_ind] = -sm_pos[*,2] ENDIF ; Use the Hammond model for radial distances > h0 (8.6) lr_ind = where(rdist GT h0, ncnt) IF ncnt GT 0 THEN BEGIN ; initialize variables x = gsm_pos[lr_ind,0] y = gsm_pos[lr_ind,1] z = gsm_pos[lr_ind,2] tilt = make_array(n_elements(x), /double) zns = make_array(n_elements(x), /double) ; check input time format and convert to doy, hr, min FOR i=0,n_elements(x)-1 DO BEGIN IF size(time, /type) EQ 5 or size(time, /type) EQ 7 THEN BEGIN time_struc = time_struct(time[lr_ind[i]]) yr = time_struc.year doy = time_struc.doy hr = time_struc.hour mm = time_struc.min sc = time_struc.sec ENDIF ELSE BEGIN print, 'Invalid time format. Format: YYYY-MM-DD/hh:mm:ss or seconds since 1970' RETURN, -1 ENDELSE ; calculate the tilt in degrees geopack_recalc, yr, doy, hr, mm, sc, tilt=tt tilt[i]=tt*rad ;convert to radians ENDFOR ; hammond model iless=where(abs(y) LT Y0, jless) IF (jless GT 0) THEN zns[iless]=((H1+D)*sqrt(1-y(iless)^2/Y0^2)-D)*sin(tilt(iless)) imore=where(abs(y) GE Y0, jmore) IF (jmore GT 0) THEN zns[imore]=-D*sin(tilt(imore)) z2ns[lr_ind] = zns dz2ns[lr_ind] = zns - z ENDIF RETURN, dz2ns END ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;NAME: ; aen_ns_model ; ;PURPOSE: This program is to find the AEN(Analytical Equatorial Neutral) sheet in the ; magnetopause in different time and position ; ;CALLING SEQUENCE: ; zaen=neutral_sheet(time, position) ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in GSM coordinates in re (pos[*,3] ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;NOTES: ; ;C References: ;(1) AEN(Analytical Equatorial Neutral): ; Zhu, M. and R.-L. Xu, 1994, A continuous neutral sheet model and a normal ; curved coordinate system in the magnetotail, Chinese J. Space Science, 14, ; (4)269, (in Chinese). ; Wang, Z.-D. and R.-L. Xu, Neutral Sheet Observed on ISEE Satellite, ; Geophysical Research Letter, 21, (19)2087, 1994. ;(2) Magnetopause model: ; Sibeck, D. G., R. E. Lopez, and R. C. Roelof, Solar wind control of the ; magnetopause shape, location, and motion, J. Grophys. Res., 96, 5489, 1991 ; ;HISTORY: ; ;----------------------------------------------------------------------------- FUNCTION aen_ns_model, time, gsm_pos, z2ns=z2ns ; initialize constants rad = !pi/180. h0 = 12.6/!pi z2ns = make_array(n_elements(time), /double) tilt = make_array(n_elements(time), /double) FOR i=0,n_elements(time)-1 DO BEGIN ; convert time into year, doy, .... IF size(time, /type) EQ 5 or size(time, /type) EQ 7 THEN BEGIN time_struc = time_struct(time[i]) yr = time_struc.year doy = time_struc.doy hr = time_struc.hour mm = time_struc.min sc = time_struc.sec ENDIF ELSE BEGIN print, 'Invalid time format. Format: YYYY-MM-DD/hh:mm:ss or seconds since 1970' RETURN, -1 ENDELSE ; calculate the tilt angle geopack_recalc, yr, doy, hr, mm, sc, tilt=tt ; calcuate the position of the neutral sheet z2NS[i] = -h0 * sin(tt*rad) * atan(gsm_pos[i,0]/5) * (2*cos(gsm_pos[i,1]/6)) ENDFOR ; calculate the spacecraft distance to the neutral sheet along the x-axes dz2ns = z2NS - gsm_pos[*,2] RETURN, dz2NS END ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;NAME: ; den_ns_model ; ;PURPOSE: ; This program finds the DEN(Displaced Equatorial Neutral) sheet inside ; the magnetopause in different tine and positions. The routine calculates ; the position along the zaxis at a specific location. ; ;CALLING SEQUENCE: ; dz2ns=neutral_sheet(time, position) ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in GSM coordinates in re (pos[*,3]) ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;NOTES: ;References: ;(1) DEN(Displaced Equatorial Neutral): ; Xu, R.-L., A Displaced Equatorial Neutral Sheet Surface Observed on ISEE-2 ; Satellite, J. Atmospheric and Terrestrial Phys., 58, 1085, 1991 ;(2) Magnetopause model: ; Sibeck, D. G., R. E. Lopez, and R. C. Roelof, Solar wind control of the ; magnetopause shape, location, and motion, J. Grophys. Res., 96, 5489, 1991 ;Original Authors of the FORTRAN source code: ;Ronglan XU and Lei LI, Center for Space Sci. and Applied Res., ;Chinese Academy of Sciences, PO Box 8701, Beijing 100080, China ;E-mail: XURL@SUN.IHEP.AC.CN, XURL@SUN20.CSSAR.AC.CN ; ;This source code was ported from the original FORTRAN source code into IDL ;The original source code only calculated to 10.05 RE. In this IDL version ;that restriction was increased to 25. ; ;HISTORY: ; ;----------------------------------------------------------------------------- ; helper function PRO SMPF, xgsm, rmp, im im = 0 rmp2 = -0.14*xgsm^2 - 18.2*xgsm+217.4 IF rmp2 LT 0 THEN BEGIN im = 9999.99 ENDIF ELSE BEGIN rmp = sqrt(rmp2) IF xgsm LE -65 THEN rmp = 28.5 ENDELSE END ; helper function PRO SFA4, aa, bb, cc, dd, x ndx = 0 xmin = 0 xmax = 50 ndxmax = 3 dx = 1 x = xmin yy = x^4 + aa*x^3 + bb*x^2 + cc*x + dd WHILE ndx LE ndxmax DO BEGIN x = x+dx IF x GE xmax THEN BEGIN ndx = 0 RETURN ENDIF y = x^4 + aa*x^3 + bb*x^2 + cc*x + dd ry = y/yy IF ry LT 0. THEN BEGIN x = x-dx dx = dx/10. ndx = ndx+1 ENDIF ELSE BEGIN yy = y ENDELSE ENDWHILE END ; helper function PRO SD1, til, H, H1, xgsm, d ct = cos(til) xx = xgsm xh = -H*ct IF xgsm GE xh THEN xx = xh ; calculate the radius of the cross section IF xx LE -5. THEN rm = 9*(10-3*xx)/(10-xx)+3 IF xx GT -5. THEN rm = sqrt(18^2-(xx+5)^2) rm2 = RM^2 ; in cross_section areas above and below the neutral ; sheet aa = 4*H-(32*rm2*H^2)/(!pi^2*H1^2*(H-xx/ct)) bb = 2*H^2*(3.-8.*rm2/(!pi^2*H1^2)) cc = 4*(H^3) dd = H^4 sfa4, aa, bb, cc, dd, x d = x IF xgsm GE xh THEN BEGIN fk = -x/sqrt(-xh) d = -fk*sqrt(-xgsm) ENDIF END FUNCTION den_ns_model, time, gsm_pos, z2ns=z2ns ; calcuate the position of the neutral sheet along z axis H = 25.5 H1 = 25.05 rad = !pi/180. done = 0 z2ns = make_array(n_elements(time), /float) FOR i = 0, n_elements(time)-1 DO BEGIN xgsm = gsm_pos[i,0] ygsm = gsm_pos[i,1] IF size(time, /type) EQ 5 or size(time, /type) EQ 7 THEN BEGIN time_struc = time_struct(time[i]) yr = time_struc.year doy = time_struc.doy hr = time_struc.hour mm = time_struc.min sc = time_struc.sec ENDIF ELSE BEGIN print, 'Invalid time format. Format: YYYY-MM-DD/hh:mm:ss or seconds since 1970' RETURN, -1 ENDELSE ; get tilt angle of magnetic pole geopack_recalc, yr, doy, hr, mm, sc, tilt=tt tilt = tt * rad IF xgsm GT -100. THEN BEGIN sd1, tilt, H, H1, xgsm, d ym21 = ((H1*(H+d))^2) * (1-(xgsm/(H*cos(tilt)))^2) ym22 = (H+d)^2 - (d-xgsm/cos(tilt))^2 ym2 = ym21/ym22 IF ym2 LT 0 THEN BEGIN ie[i] = 2 CONTINUE ENDIF ym = sqrt(ym2) xd2 = ((H*cos(tilt))^2) * (1-(ygsm/H1)^2) IF abs(ygsm) GT H1 THEN xd2 = 0 ; find the equatorial region xd = sqrt(xd2) rd = sqrt(xd^2+ygsm^2) rsm = sqrt(xgsm^2+ygsm^2) IF xgsm GT 0 OR rsm le rd THEN BEGIN z2ns[i] = -xgsm*sin(tilt)/cos(tilt) done = 1 ENDIF IF abs(ygsm) GT ym && done NE 1 THEN BEGIN z2ns[i] = -d*sin(tilt) done = 1 ENDIF IF done NE 1 THEN BEGIN z2ns[i] = ((H+d)*sqrt(1-(ygsm^2)/ym2)-d)*sin(tilt) ENDIF ENDIF ; end of xgsm LE 0. ENDFOR dz2ns = z2NS - gsm_pos[*,2] RETURN, dz2ns END ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;NAME: ; fairfield_NS_model ; ;PURPOSE: ; This routine calculates the position along the zaxis at a specific ; x and y location. The Fairfield model is used to this calculation. ; ;CALLING SEQUENCE: ; zaen=neutral_sheet(time, position, [rmp=rmp], [ie=ie] ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in GSM coordinates in re (pos[*,3]) ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;NOTES: ;Reference: ; A statistical determination of the shape and position of the ; geomagnetic neutral sheet, Journal of Geophysical Research, ; Vol. 85, No A2, pages 775-780, February 1, 1980 ; Author - D. Fairfield ; ;HISTORY: ; ;----------------------------------------------------------------------------- FUNCTION fairfield_ns_model, time, gsm_pos, z2ns=z2ns ; constants (in re) h0 = 10.5 y0 = 22.5 d = 14. rad = !pi/180. tilt = make_array(n_elements(time), /double) FOR i=0,n_elements(time)-1 DO BEGIN ; check input time format and convert to doy, hr, min IF size(time, /type) EQ 5 OR size(time, /type) EQ 7 THEN BEGIN time_struc = time_struct(time[i]) yr = time_struc.year doy = time_struc.doy hr = time_struc.hour mm = time_struc.min sc = time_struc.sec ENDIF ELSE BEGIN print, 'Invalid time format. Format: YYYY-MM-DD/hh:mm:ss or seconds since 1970' RETURN, -1 ENDELSE ; calculate tilt angle of geomagnetic axis geopack_recalc, yr, doy, hr, mm, sc, tilt=tt tilt[i] = tt*rad ENDFOR ; calcuate the position of the neutral sheet along z axis z2ns = ((h0 + d) * sqrt(1 - gsm_pos[*,1]^2/y0^2) - d)*sin(tilt) dz2ns = z2ns - gsm_pos[*,2] RETURN, dz2NS END ;+ ;NAME: ; den_fairfield_ns_model ; ;PURPOSE: ; This routine calculates the position along the zaxis at a specific ; x and y location. ; ;CALLING SEQUENCE: ; dz2ns=den_fairfield_ns_model(time, position) ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in GSM coordinates in re (pos[*,3] ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;HISTORY: ; ;----------------------------------------------------------------------------- FUNCTION den_fairfield_ns_model, time, gsm_pos, z2ns=z2ns ; initialize constants rad = !pi/180. z2ns = make_array(n_elements(time), /double) ; Use the den model for radial distances <12.re rdist = sqrt(gsm_pos[*,0]^2 + gsm_pos[*,1]^2 + gsm_pos[*,2]^2) sm_ind = where(rdist LE 10., ncnt) IF ncnt GT 0 THEN z2ns[sm_ind] = (den_ns_model(time[sm_ind], gsm_pos[sm_ind,*])) ; use the fairfield model for radial distances >12.re lr_ind = where(rdist GT 10., ncnt) IF ncnt GT 0 THEN z2ns[lr_ind] = (fairfield_ns_model(time[lr_ind], gsm_pos[lr_ind,*])) dz2ns = z2ns RETURN, dz2NS END ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;NAME: ; lopez_NS_model ; ;PURPOSE: ; This routine calculates the position along the zaxis at a specific ; x and y location. The Lopez model is used for this calculation. ; ;CALLING SEQUENCE: ; dz2NS=lopez_NS_model(time, position, kp, mlt) ; ;INPUT: ; time - string or double format ; double(s) seconds since 1970 ; string(s) format: YYYY-MM-DD/hh:mm:ss ; gsm_pos - position vector in re, GSM Coordinates (pos[*,3]) ; kp - kp index value ; mlt - magnetic local time in degrees (0=midnight) ; ;OUTPUT: ; dz2NS - distance of the spacecraft to the neutral sheet along z axis in re ; ;KEYWORD [Optional]: ; z2NS - position of the neutral sheet along the z-axis in re ; ;NOTES: ;Reference: ; The position of the magnetotail neutral sheet in the near-Earth Region, ; Geophysical Research Letters, Vol. 17, No 10, pages 1617-1620, 1990 ; Author - Ramon E. Lopez ; ; The lopez model is best used for distances <8.8 RE ; ;HISTORY: ; ;----------------------------------------------------------------------------- FUNCTION lopez_ns_model, time, gsm_pos, kp=kp, mlt=mlt, z2ns=z2ns ; constants re = 6378. rad = !pi/180. IF ~keyword_set(kp) THEN kp = 0 IF ~keyword_set(mlt) THEN mlt = 0.0 tilt = make_array(n_elements(time), /double) FOR i=0,n_elements(time)-1 DO BEGIN ; check input time format and convert to doy, hr, min IF size(time, /type) EQ 5 or size(time, /type) EQ 7 THEN BEGIN time_struc = time_struct(time[i]) yr = time_struc.year doy = time_struc.doy hr = time_struc.hour mm = time_struc.min sc = time_struc.sec ENDIF ELSE BEGIN print, 'Invalid time format. Format: YYYY-MM-DD/hh:mm:ss or seconds since 1970' RETURN, -1 ENDELSE ; calculate tilt angle of geomagnetic axis geopack_recalc, yr, doy, hr, mm, sc, tilt=tt tilt[i] = tt * rad ENDFOR ; calcuate the position of the neutral sheet along z axis rdist = sqrt(gsm_pos[*,0]^2 + gsm_pos[*,1]^2 + gsm_pos[*,2]^2) mlat = -(0.14*Kp + 0.69) * ((cos(rad*mlt))^.3333333) * (0.065*(rdist^0.8) - 0.16) * tilt ; convert magnetic latitude to position rthph2xyz, rdist, mlat/rad, mlt, x, y, z cotrans, gsm_pos, sm_pos, time, /gsm2sm, /ignore_dlimits sm_pos[*,2]=z cotrans, sm_pos, new_pos, time, /sm2gsm, /ignore_dlimits z2NS = new_pos[*,2] dz2NS = z2NS - gsm_pos[*,2] RETURN, dz2NS END ; Helper function for the lopez model ; converts spherical to cartesian coordinates ; NOTE: th,ph in degrees, and th is latitude (not colatitude) (i.e. [-90->90]) PRO rthph2xyz,r,th,ph,x,y,z FLAG=6.8792E+28 FLAG98=0.98*FLAG PI=3.1415926535898 thrad=th*PI/180. phrad=ph*PI/180. sth=sin(thrad) cth=cos(thrad) sph=sin(phrad) cph=cos(phrad) x=r*cth*cph y=r*cth*sph z=r*sth iflags=where((r GT FLAG98) OR (th GT FLAG98) OR (ph GT FLAG98), iany) IF (iany GT 0) THEN BEGIN x(iflags)=FLAG y(iflags)=FLAG z(iflags)=FLAG ENDIF RETURN END PRO neutral_sheet_models ; does nothing. ; call cotrans_lib at the beginning of any routine ; that needs to use any cotrans_lib routines, to ensure ; that they are compiled. END