;+
;Procedure: trace2iono
;
;Purpose: Generates a model field line footprint from an array of given
; positions and times,will also trace field lines at the user's request
; This program will always use the refined foot point mappings(a
; modification made by Vassilis Angelopoulos to provide more accurate
; mappings of foot points) unless you request the standard mapping;
;
;Keywords:
; tarray: N length array storing the position times in seconds utc
;
; in_pos_array: Nx3 array representing the position series (in
; km gsm by default)
;
; out_foot_array: named variable in which to store the footprints
; calculated by this function(will be an Nx3 array) will be returned
; (in RE gsm by default)
;
; out_trace_array(optional): named variable in which to store the traces of
; field lines leading to footprints. Because traces are of variable
; length, the returned array will be of dimensions NxDx3
; D is the maximum number of vectors in any of the traces.
; Shorter traces will have NaNs filling the space at the end
; of the array
;
; in_coord(optional): set this keyword to a string indicating the
; coordinate system input position data is in.
; (can be 'gei','geo','gse','gsm',or 'sm' default: gsm)
;
; out_coord(optional): set this keyword to a string indicating the
; coordinate system output data should be in.
; (can be 'gei','geo','gse','gsm',or 'sm' default: gsm)
;
; internal_model(optional): set this keyword to a string
; indicating the internal model that should be used in tracing
; (can be 'dip' or 'igrf' default:igrf)
;
; external_model(optional): set this keyword to a string
; indicating the external model that should be used in tracing
; (can be 'none','t89','t96','t01', or 't04s' default: none)
;
; /SOUTH(optional): set this keyword to indicate that fields
; should be traced towards the southern hemisphere. By default
; they trace north.
;
; /KM(optional): set this keyword to indicate that input
; and output will be in KM not RE
;
; par(optional): parameter input for the external field model
; if using t89 then it should be an N element array containing
; kp values or a single kp value, if using t96,t01,t04s it
; should be an N x 10 element array of parmod values or a 10
; element array or a 1x10 element array. At the moment if an
; external model is set and this is not set an error will be thrown.
;
; period(optional): the amount of time between recalculations of
; geodipole tilt and input of new model parameters in
; seconds (default: 60) increase this value to decrease
; run time
; if field line traces are requested this parameter is
; ignored and new model parameters are input on each
; iteration
;
; error(optional): named variable in which to return the error state
; of the procedure. 1 for success, 0 for failure
;
; standard_mapping(optional): Set to use Tsyganenko's
; unmodified version instead of the Angelopoulos's
; refined version
;
; R0(optional): Minimum trace distance in RE,unless /km is set
;
; RLIM(optional): Maximum trace distance in RE,unless /km is
; set (default: 60 RE)
;
; /NOBOUNDARY(optional): Override boundary limits.
;
; /STORM(optional): Specify storm-time version of T01 external
; magnetic field model use together with /T01.
;
;
;Example: trace2iono,in_time,in_pos,out_foot
;
;Notes:
; 1. Relies on the IDL/Geopack Module provided by Haje Korth JHU/APL
; and N.A. Tsyganenko NASA/GSFC, if the module is not installed
; this function will fail.
; 2. Has a loop with number of iterations =
; (tarray[n_elements(t_array)]-tarray[0])/period
; This means that as period becomes smaller the amount time of this
; function should take will grow quickly.
; 3. If the trace_array variable is set
; the period variable will be ignored. The program will
; recalculate for each value, this will cause the program to
; run very slowly.
; 4. All calculations are done internally in double precision
;
;
; $LastChangedBy: jimm $
; $LastChangedDate: 2008-05-13 18:17:46 -0700 (Tue, 13 May 2008) $
; $LastChangedRevision: 3082 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/idl_socware/tags/tdas_5_02/external/IDL_GEOPACK/trace/trace2iono.pro $
;-
pro trace2iono,tarray,in_pos_array,out_foot_array, out_trace_array = out_trace_array, in_coord = in_coord, out_coord = out_coord, internal_model = internal_model, external_model = external_model, south = south, km = km, par=par,period=period,error = error, standard_mapping=standard_mapping,r0=r0,rlim=rlim,_extra=_extra
error = 0
;constant arrays used for input validation
valid_coords = ['gei', 'gse', 'geo','gsm', 'sm']
valid_internals = ['dip', 'igrf']
valid_externals = ['none', 't89', 't96', 't01', 't04s']
km_in_re = 6374.4D
if not keyword_set(rlim) then $
if keyword_set(km) then $
rlim = 60D*km_in_re $
else $
rlim = 60D
;test to make sure idl/geopack library is installed
if igp_test() eq 0 then return
if not keyword_set(tarray) then begin
message, /continue, 'tarray must be set'
return
endif
if not keyword_set(in_pos_array) then begin
message, /continue, 'in_pos_array must be set'
return
endif
;be sure to test this predicate
if not arg_present(out_foot_array) then begin
message, /continue, 'out_foot_array must be set'
return
endif
if keyword_set(in_coord) then begin
in_coord2 = strlowcase(in_coord)
if(strfilter(valid_coords, in_coord2) eq '') then begin
message, /continue, 'in_coord not a valid coordinate name'
return
endif
endif else in_coord2 = 'gsm'
if keyword_set(out_coord) then begin
out_coord2 = strlowcase(out_coord)
if(strfilter(valid_coords, out_coord2) eq '') then begin
message, /continue, 'out_coord not a valid coordinate name'
return
endif
endif else out_coord2 = 'gsm'
if keyword_set(internal_model) then begin
internal_model2 = strlowcase(internal_model)
if(strfilter(valid_internals, internal_model2) eq '') then begin
message, /continue, 'internal_model not a valid internal model name'
return
endif
endif else internal_model2 = 'igrf'
if keyword_set(external_model) then begin
external_model2 = strlowcase(external_model)
if(strfilter(valid_externals, external_model2) eq '') then begin
message, /continue, 'external_model not a valid external model name'
return
endif
endif else external_model2 = 'none'
if keyword_set(south) then dir = 1.0D $
else dir = -1.0D
;convert inputs into double precision to ensure consistency of calculations
tarray2 = double(tarray)
in_pos_array2 = double(in_pos_array)
if keyword_set(r0) then r02 = double(r0)
if keyword_set(rlim) then rlim2 = double(rlim)
if keyword_set(km) then begin
in_pos_array2 = in_pos_array2/km_in_re
if keyword_set(r02) then r02 = r02/km_in_re
if keyword_set(rlim2) then rlim2 = rlim2/km_in_re
endif else begin
idx_temp = where(abs(in_pos_array2) gt km_in_re)
if idx_temp[0] ne -1L then begin
message,/continue,'!!!! WARNING !!!! the magnitude of your rgsm values suggests your data may be in km'
message,/continue,'Default is Re, please set keyword "km" or see calling sequence by typing doc_library,"trace2iono".'
endif
endelse
if internal_model2 eq 'igrf' then IGRF = 1 else IGRF = 0
if external_model2 ne 'none' and not keyword_set(par) then begin
message,/continue,'par must be set if external model is set'
return
endif
if external_model2 eq 't89' then begin
;these switches used to tell the
;IDL/GEOPACK trace function which
;model to use
T89 = 1
T96 = 0
T01 = 0
TS04 = 0
;intialize and check par array
if size(par,/n_dim) eq 0 then par_array = make_array(n_elements(tarray2),/DOUBLE,value=par)
if size(par,/n_dim) eq 1 then begin
if n_elements(par) ne n_elements(tarray2) then begin
message,/continue,'par must have the same number of elements as tarray '
return
endif else par_array = double(par)
endif
if size(par,/n_dim) gt 1 then begin
message,/continue,'par must have 0 or 1 dimensions when used with model t89'
return
endif
par_idx_low = where(par_array lt 1)
if par_idx_low[0] ne -1L then begin
message, /continue, 'par has value less than 1'
return
endif
par_idx_high = where(par_array gt 7)
if par_idx_high[0] ne -1L then begin
message, /continue, 'par has value greater than 7'
return
endif
endif else if external_model2 ne 'none' then begin
;these switches used to tell the
;IDL/GEOPACK trace function which
;model to use
T89 = 0
if external_model2 eq 't96' then T96 = 1 else T96 = 0
if external_model2 eq 't01' then T01 = 1 else T01 = 0
if external_model2 eq 't04s' then TS04 = 1 else TS04 = 0
;check par array
s = size(par,/dimensions)
if n_elements(s) eq 1 && s[0] eq 10 then par_array = transpose(rebin(par,10,n_elements(tarray2))) else $
if n_elements(s) eq 2 && s[0] eq 1 && s[1] eq 10 then par_array = rebin(par,n_elements(tarray2),10) else $
if n_elements(s) ne 2 || s[1] ne 10 || s[0] ne n_elements(tarray2) then begin
message,/continue,'par must be an N x 10 array or 10 element array if ' + external_model2 + ' is set'
return
endif else par_array = double(par)
endif else begin
T89 = 0
T96 = 0
T01 = 0
TS04 = 0
endelse
;check input array dimensions
t_size = size(tarray2, /dimensions)
r_size = size(in_pos_array2, /dimensions)
if n_elements(t_size) ne 1 then begin
message, /continue, 'tarray has incorrect dimensions'
return
endif
if n_elements(r_size) ne 2 || r_size[1] ne 3 then begin
message, /continue, 'in_pos_array has incorrect dimensions'
return
endif
if t_size[0] ne r_size[0] then begin
message, /continue, 'number of times in tarray does not match number of positions in in_pos_array'
return
endif
;all calculations will be done in gsm internally
if in_coord2 eq 'gei' then begin
cotrans,in_pos_array2,in_pos_array2,tarray2,/gei2gse
cotrans,in_pos_array2,in_pos_array2,tarray2,/gse2gsm
endif else if in_coord2 eq 'geo' then begin
cotrans,in_pos_array2,in_pos_array2,tarray2,/geo2gei
cotrans,in_pos_array2,in_pos_array2,tarray2,/gei2gse
cotrans,in_pos_array2,in_pos_array2,tarray2,/gse2gsm
endif else if in_coord2 eq 'gse' then $
cotrans,in_pos_array2,in_pos_array2,tarray2,/gse2gsm $
else if in_coord2 eq 'sm' then $
cotrans,in_pos_array2,in_pos_array2,tarray2,/sm2gsm
;intialize and check period
if not keyword_set(period) then period2 = 60.0D $
else period2 = double(period)
if period2 le 0. then begin
message, /contiune, 'period must be positive'
return
endif
;defaults to NaN so it will plot properly in tplot and to prevent
;insertion of spurious default dindgen values
out_foot_array = make_array(r_size, /DOUBLE, VALUE = !VALUES.D_NAN)
ts = time_struct(tarray2)
;loop boundaries and storage if traces requested
if arg_present(out_trace_array) then begin
;if traces are requested every point is processed individually
ct = t_size[0] - 1L
;allocate space for traces, since traces may have different lengths
;pointers are allocated
tr_ptr_arr = ptrarr(t_size[0])
;the maximum trace size will be stored so we'll know how large to
;make the array that is ultimately output
max_trace_size = 0
endif else begin ;loop boundaries if traces are not requested
tstart = tarray2[0]
tend = tarray2[t_size[0] - 1L]
;number of iterations is interval length divided by period length
ct = ceil((tend-tstart)/period2)
endelse
i = 0L
while i le ct do begin
;call for each point individually if field line traces are requrested
if arg_present(out_trace_array) then begin
;recalculate magnetic dipole
geopack_recalc, ts[i].year,ts[i].doy, ts[i].hour, ts[i].min, ts[i].sec, tilt = tilt
;calculate which par values should be used on this iteration
if T89 eq 1 then par_iter = par_array[i] $
else if T96 eq 1 || T01 eq 1 || TS04 eq 1 then par_iter = par_array[i,*] $
else par_iter = ''
; geopack_trace,in_pos_array2[i,0],in_pos_array2[i,1],in_pos_array2[i,2],dir,par_iter,out_foot_array[i,0],out_foot_array[i,1],out_foot_array[i,2],R0=R02,RLIM=RLIM2,fline = trgsm_out,tilt=tilt,IGRF=IGRF,T89=T89,T96=T96,T01=T01,TS04=TS04,/refine,/ionosphere,_extra=_extra
if keyword_set(standard_mapping) then $
geopack_trace, in_pos_array2[i, 0], in_pos_array2[i, 1], in_pos_array2[i, 2], dir, par_iter, out_foot_x, out_foot_y, out_foot_z, R0 = R02, RLIM = RLIM2, fline = trgsm_out, tilt = tilt, IGRF = IGRF, T89 = T89, T96 = T96, T01 = T01, TS04 = TS04, _extra = _extra $
else $
geopack_trace, in_pos_array2[i, 0], in_pos_array2[i, 1], in_pos_array2[i, 2], dir, par_iter, out_foot_x, out_foot_y, out_foot_z, R0 = R02, RLIM = RLIM2, fline = trgsm_out, tilt = tilt, IGRF = IGRF, T89 = T89, T96 = T96, T01 = T01, TS04 = TS04, /refine, /ionosphere, _extra = _extra
out_foot_array[i, 0] = out_foot_x
out_foot_array[i, 1] = out_foot_y
out_foot_array[i, 2] = out_foot_z
;store pointer to field trace
tr_ptr_arr[i] = ptr_new(trgsm_out)
tr_size = size(trgsm_out,/dimensions)
;store the maximum trace size so it
;can be used to calculate output storage
if tr_size[0] gt max_trace_size then max_trace_size = tr_size[0]
endif else begin ;calculate over an interval if traces are not requested
;find indices of points in the interval for this iteration
idx1 = where(tarray2 ge tstart + i*period2)
idx2 = where(tarray2 le tstart + (i+1)*period2)
idx = ssl_set_intersection(idx1, idx2)
if idx[0] ne -1L then begin
id = idx[0]
;recalculate geomagnetic dipole
geopack_recalc, ts[id].year,ts[id].doy, ts[id].hour, ts[id].min, ts[id].sec, tilt = tilt
rgsm_x = in_pos_array2[idx, 0]
rgsm_y = in_pos_array2[idx, 1]
rgsm_z = in_pos_array2[idx, 2]
;calculate which par values should be used on this iteration
if T89 eq 1 then par_iter = par_array[id] $
else if T96 eq 1 || T01 eq 1 || TS04 eq 1 then par_iter = par_array[id,*] else par_iter = ''
if keyword_set(standard_mapping) then $
geopack_trace,rgsm_x,rgsm_y,rgsm_z,dir,par_iter,foot_x,foot_y,foot_z,R0=R02,RLIM=RLIM2,tilt=tilt,IGRF=IGRF,T89=T89,T96=T96,T01=T01,TS04=TS04,_extra=_extra $
else $
geopack_trace,rgsm_x,rgsm_y,rgsm_z,dir,par_iter,foot_x,foot_y,foot_z,R0=R02,RLIM=RLIM2,tilt=tilt,IGRF=IGRF,T89=T89,T96=T96,T01=T01,TS04=TS04,/refine,/ionosphere,_extra=_extra
;output foot
out_foot_array[idx, 0] = foot_x
out_foot_array[idx, 1] = foot_y
out_foot_array[idx, 2] = foot_z
endif
endelse
i++
endwhile
;copy pointer version of trace into output array
if arg_present(out_trace_array) then begin
out_trace_array = make_array([t_size[0],max_trace_size,3],/DOUBLE, VALUE = !VALUES.D_NAN)
for i = 0L,t_size[0]-1L do begin
tr_temp = *tr_ptr_arr[i]
ptr_free,tr_ptr_arr[i]
s_temp = size(tr_temp,/dimensions)
;all points within each trace have the same time
t_temp = replicate(tarray2[i],s_temp[0])
;convert trace into the output coordinate system
if out_coord2 eq 'gei' then begin
cotrans,tr_temp,tr_temp,t_temp,/gsm2gse
cotrans,tr_temp,tr_temp,t_temp,/gse2gei
endif else if out_coord2 eq 'geo' then begin
cotrans,tr_temp,tr_temp,t_temp,/gsm2gse
cotrans,tr_temp,tr_temp,t_temp,/gse2gei
cotrans,tr_temp,tr_temp,t_temp,/gei2geo
endif else if out_coord2 eq 'gse' then $
cotrans,tr_temp,tr_temp,t_temp,/gsm2gse $
else if out_coord2 eq 'sm' then $
cotrans,tr_temp,tr_temp,t_temp,/gsm2sm
out_trace_array[i,0:(s_temp[0]-1L),*] = tr_temp
endfor
;convert from re to km
if keyword_set(km) then out_trace_array *= km_in_re
endif
;convert footprint into the output coordinate system
if out_coord2 eq 'gei' then begin
cotrans,out_foot_array,out_foot_array,tarray2,/gsm2gse
cotrans,out_foot_array,out_foot_array,tarray2,/gse2gei
endif else if out_coord2 eq 'geo' then begin
cotrans,out_foot_array,out_foot_array,tarray2,/gsm2gse
cotrans,out_foot_array,out_foot_array,tarray2,/gse2gei
cotrans,out_foot_array,out_foot_array,tarray2,/gei2geo
endif else if out_coord2 eq 'gse' then $
cotrans,out_foot_array,out_foot_array,tarray2,/gsm2gse $
else if out_coord2 eq 'sm' then $
cotrans,out_foot_array,out_foot_array,tarray2,/gsm2sm
;convert from re to km
if keyword_set(km) then out_foot_array *= km_in_re
;signal success
error = 1
end