;+
;Procedure: trace2equator
;
;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) as there is no Tsyganenko code to 
;         perform equatorial traces.
;
;Input:
;         tarray: N length array storing the position times in seconds utc since 1970
;
;         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)
;
;Keywords:
;         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','gse','geo','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','gse','geo','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
;             
;         R0(optional):  radius of a sphere (in re), defining the inner boundary of the tracing region
;         (usually, earth's surface or the ionosphere, where r0~1.0)
;         if the field line reaches that sphere from outside, its inbound tracing is
;         terminated and the crossing point coordinates xf,yf,zf  are calculated.
;         (units are km if /km is set)
;
;         RLIM(optional) - radius of a sphere (in re), defining the outer boundary of the tracing region;
;         if the field line reaches that boundary from inside, its outbound tracing is
;         terminated and the crossing point coordinates xf,yf,zf are calculated.(default 60 RE)
;         (units are km if /km is set)   
;
;         NOBOUNDARY(optional): Override boundary limits.
;
;         STORM(optional): Specify storm-time version of T01 external 
;             magnetic field model use together with /T01.
;
;         add_tilt:  Increment the default dipole tilt used by the model with
;             a user provided tilt in degrees.  Result will be produced with TSY_DEFAULT_TILT+ADD_TILT
;             Value can be set to an N length array an M length array or a single element array. 
;             N is the number of time elements for the data.  M is the number of periods in the time interval.(determined by the period keyword)
;             If single element is provided the same correction will be applied to all periods.   
;             If an N length array is provided, the data will be re-sampled to an M length array. Consequently, if
;             the values change quickly, the period may need to be shortened. 
;         
;         get_tilt: Returns the dipole_tilt parameter used for each period. 
;             Returned value has a number of elements equal to the value returned by get_nperiod
;         
;         set_tilt: Alternative alternative dipole_tilt value rather than the geopack tilt.
;             This input can be an M length array, and N length array or a single elemnt.
;             Value can be set to an N length array an M length array or a single element array. 
;             N is the number of time elements for the data.  M is the number of periods in the time interval.(determined by the period keyword)
;             If an N length array is provided, the data will be re-sampled to an M length array. Consequently, if
;             the values change quickly, the period may need to be shortened. 
;             set_tilt will cause add_tilt to be ignored. 
;
;         exact_tilt_times (optional):  Set this keyword to avoid grouping similar times (default 10 minutes) and instead
;              recalculate the dipole tilt at each input time
;                   
;         get_nperiod: Returns the number of periods used for the time interval=  ceil((end_time-start_time)/period)
;     
;         geopack_2008 (optional): Set this keyword to use the latest version (2008) of the Geopack
;              library. Version 9.2 of the IDL Geopack DLM is required for this keyword to work.
;              
;          ts07_param_dir (optional): Location of TS07 parameter directory
;          
;          ts07_param_file (optional): Location of TS07 parameter file
;          
;          skip_ts07_load (optional): Do not reset TS07 parameter directory or reload parameter file
;          
;Example: 
;        trace2equator,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: jwl $
; $LastChangedDate: 2021-10-19 17:05:44 -0700 (Tue, 19 Oct 2021) $
; $LastChangedRevision: 30381 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_5_0/external/IDL_GEOPACK/trace/trace2equator.pro $
;-

pro trace2equator, 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, r0=r0, rlim=rlim, add_tilt=add_tilt, $
    get_tilt=get_tilt, set_tilt=set_tilt, get_nperiod=get_nperiod, get_period_times=get_period_times, $
    geopack_2008=geopack_2008, exact_tilt_times=exact_tilt_times, $
    ts07_param_dir=ts07_param_dir, ts07_param_file=ts07_param_file, skip_ts07_load=skip_ts07_load, _extra=_extra

    if undefined(geopack_2008) then geopack_2008=0
    if undefined(exact_tilt_times) then exact_tilt_times=0
    if undefined(skip_ts07_load) then skip_ts07_load=0

    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', 'ts07','ta15b','ta15n']
    
    ; 6371.2 = the value used in the GEOPACK FORTRAN code for Re
    km_in_re = 6371.2D
    
    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(geopack_2008=geopack_2008) 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 n_elements(r0) gt 0 then r02 = double(r0)
    if n_elements(rlim) gt 0 then rlim2 = double(rlim)
    
    if keyword_set(km) then begin
       in_pos_array2 = in_pos_array2/km_in_re
       if n_elements(r02) gt 0 then r02 = r02/km_in_re
       if n_elements(rlim2) gt 0 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,"trace2equator".'
       endif
    endelse
    
    if internal_model2 eq 'igrf' then IGRF = 1 else IGRF = 0
    
    if external_model2 ne 'none' and (n_elements(par) eq 0) 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
      TS07 = 0
      TA15B = 0
      TA15N = 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
    
      ; Range check for Kp
      ; Valid IOPT values are in the range [1.0,7.0]
      ; kp2iopt will replace any out-of-range iopt values with 1.0 or 7.0 as appropriate,
      ; so we only warn here instead of throwing an error.
      
      par_idx_low = where(par_array lt 0.0)
    
      if par_idx_low[0] ne -1L then begin
          message, /continue, 'par (Kp) has value less than 0 (will be treated as 0)'
      endif
    
      par_idx_high = where(par_array gt 6.0)
    
      if par_idx_high[0] ne -1L then begin
          message, /continue, 'par (Kp) has value greater than 6 (will be treated as 6)'
      endif
      
      par_array=kp2iopt(par_array)  ; Convert Kp values to valid iopt values
      
    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
       if external_model2 eq 'ts07' then TS07 = 1 else TS07 = 0
       if external_model2 eq 'ta15n' then TA15N=1 else TA15N = 0
       if external_model2 eq 'ta15b' then TA15B=1 else TA15B = 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 if ' + external_model2 + ' is set'
          return
       endif else par_array = double(par)
    endif else begin
       T89 = 0
       T96 = 0
       T01 = 0
       TS04 = 0
       TS07 = 0
       TA15B = 0
       TA15N = 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
    
    ts = time_struct(tarray2)
    
    tstart = tarray2[0]
    
    ;all calculations will be done in GSM (or GSW, for Geopack 2008) internally
    if geopack_2008 then begin
        if in_coord2 eq 'gei' then begin
           cotrans,in_pos_array2,in_pos_array2,tarray2,/gei2gse
        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
        endif else if in_coord2 eq 'sm' then begin 
           cotrans,in_pos_array2,in_pos_array2,tarray2,/sm2gse
        endif else if in_coord2 eq 'gsm' then begin
           cotrans,in_pos_array2,in_pos_array2,tarray2, /gsm2gse
        endif 
        ; cotrans transformed the coordinates to GSE, use Geopack to transform to GSW
        geopack_recalc_08, ts[0].year, ts[0].doy, ts[0].hour, ts[0].min, ts[0].sec, tilt = tilt
        geopack_conv_coord_08, in_pos_array2[*,0], in_pos_array2[*,1], in_pos_array2[*,2], x_out_gsw, y_out_gsw, z_out_gsw, /from_gse, /to_gsw
        in_pos_array2 = [[x_out_gsw], [y_out_gsw], [z_out_gsw]]
    endif else begin
        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
    endelse
    
    ;intialize and check period
    if not keyword_set(period) then period2 = 60.0D  $
    else period2 = double(period)
    
    if period2 le 0. then begin
      message, /continue, '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)
    
    
    idx = where((tarray2[1:n_elements(tarray2)-1] - tarray2[0:n_elements(tarray2)-2]) lt 0,nonmonotone_times)
      
    if nonmonotone_times gt 0 then begin
      dprint,'Warning some times are non monotonic, this may cause unreliable results'
    endif
    
    ;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
       if ~exact_tilt_times then begin
        tend = tarray2[t_size[0] - 1L]
        ;number of iterations is interval length divided by period length
        ct = ceil((tend-tstart)/period2)
       endif else begin
        ct = t_size[0] - 1L
       endelse           
    endelse
    
    nperiod = ct+1
    
    if arg_present(get_nperiod) then begin
      get_nperiod = nperiod
    endif
    
    if arg_present(get_tilt) then begin
      get_tilt = dblarr(nperiod)
    endif
    
    if arg_present(get_period_times) then begin
      if arg_present(out_trace_array) or (exact_tilt_times) then begin
        ; If traces are requested, or if exact_tilt_times keyword specified, tilt times are identical to input times
        get_period_times = tarray
      endif else begin
        ; Otherwise, divide interval up into periods and calculate the tilt at the start of each period.
        get_period_times = tstart + dindgen(nperiod)*period2
      endelse
    endif
    
    if n_elements(add_tilt) gt 0 then begin
      if n_elements(add_tilt) eq 1 then begin
        tilt_value = replicate(add_tilt[0],nperiod)
      endif else if n_elements(add_tilt) eq nperiod then begin
        tilt_value = add_tilt
      endif else if n_elements(add_tilt) eq t_size[0] then begin
      if ~exact_tilt_times then begin
        period_abcissas = tstart + dindgen(nperiod)*period
      endif else begin
        period_abcissas = tarray
      endelse
        tilt_value = interpol(add_tilt,tarray,period_abcissas)
      endif else begin
        dprint,'Error: add_tilt values do not match data values or period values'
        return
      endelse
    endif
    
    if n_elements(set_tilt) gt 0 then begin
      if n_elements(set_tilt) eq 1 then begin
        tilt_value = replicate(set_tilt[0],nperiod)
      endif else if n_elements(set_tilt) eq nperiod then begin
        tilt_value = set_tilt
      endif else if n_elements(set_tilt) eq t_size[0] then begin
        ;resample tilt values to period intervals, using middle of sample
      if ~exact_tilt_times then begin
        period_abcissas = tstart + dindgen(nperiod)*period
      endif else begin
        period_abcissas = tarray
      endelse
         tilt_value = interpol(set_tilt,tarray,period_abcissas)
      endif else begin
        dprint,'Error: set_tilt values do not match data values or period values'
        return
      endelse
    endif
 
    if TS07 eq 1 then begin

      if ~keyword_set(skip_ts07_load) then begin
        
         ; find start, end years and download parameter files
         time_start = strmid(time_string(tarray[0]), 0, 4)
         time_end = strmid(time_string(tarray[n_elements(tarray)-1]), 0, 4)
   
         if time_end gt time_start then begin
           tdiff = 0 + time_end - time_start
           years = [time_start]
           for i=1, tdiff do begin
             years = [years, time_start + i]
           endfor
         endif else ts07_years = [time_start]
   
         ts07_download, years=ts07_years
   
         ; Directory for model parameters.
         if not keyword_set(ts07_param_dir) then begin
        ts07_param_dir = !spedas.geopack_param_dir
         endif else if ~file_test(ts07_param_dir, /READ, /directory) then begin
           ts07_param_dir = !spedas.geopack_param_dir
         endif
         ; Directory that contains the coeficient files
         GEOPACK_TS07_SETPATH, ts07_param_dir[0]
         if keyword_set(ts07_param_file) then begin
           ; Coeficcient filename only, without the full path
           GEOPACK_TS07_LOADCOEF, file_basename(ts07_param_file[0])
         endif else begin
           dprint,'Error: ts07_param_file must be specified if using model ts07'
           return
         endelse
    endif else begin
        dprint,'Skipping loading of TS07 parameters'
    endelse
  endif

   
    i = 0L
 
    tilt = 0.0D   ; Ensure tilt is always defined
       
    while i le ct do begin
 
      ; Default to last calculated value, in case no points lie in this interval     
      if n_elements(get_tilt) gt 0 then begin
        get_tilt[i] = tilt
      endif

       ;call for each point individually if field line traces are requrested
       if arg_present(out_trace_array) then begin
    
          ;recalculate magnetic dipole
          if geopack_2008 then begin
            geopack_recalc_08, ts[i].year,ts[i].doy, ts[i].hour, ts[i].min, ts[i].sec, tilt = tilt
          endif else begin
            geopack_recalc, ts[i].year,ts[i].doy, ts[i].hour, ts[i].min, ts[i].sec, tilt = tilt
          endelse
    
          ;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 || TS07 eq 1 || TA15B eq 1 || TA15N eq 1 then par_iter = par_array[i,*] $
          else par_iter = ''
    
          ;account for user tilt.
          if n_elements(tilt_value) gt 0 then begin
            if n_elements(set_tilt) gt 0 then begin
              tilt = tilt_value[i]
            endif else if n_elements(add_tilt) gt 0 then begin
              tilt = tilt+tilt_value[i]
            endif
          endif
          
          if n_elements(get_tilt) gt 0 then begin
            get_tilt[i] = tilt
          endif
    
      ;    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,_extra=_extra,/REFINE,/EQUATOR
          if geopack_2008 then begin
            geopack_trace_08,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, TS07 = TS07, TA15B = TA15B, TA15N = TA15N, /refine, /equator, _extra = _extra
          endif else begin
            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, TS07 = TS07,TA15B = TA15B, TA15N = TA15N, /refine, /equator, _extra = _extra
          endelse
          
          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
          if ~exact_tilt_times then begin
            ;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)
          endif else begin
            idx = i
          endelse
     
          if idx[0] ne -1L then begin 
          
             id = idx[0]
    
             ;recalculate geomagnetic dipole
             if geopack_2008 then begin
                geopack_recalc_08, ts[id].year,ts[id].doy, ts[id].hour, ts[id].min, ts[id].sec, tilt = tilt
             endif else begin
                geopack_recalc, ts[id].year,ts[id].doy, ts[id].hour, ts[id].min, ts[id].sec, tilt = tilt
             endelse
             
             ;account for user tilt.
             if n_elements(tilt_value) gt 0 then begin
               if n_elements(set_tilt) gt 0 then begin
                 tilt = tilt_value[i]
               endif else if n_elements(add_tilt) gt 0 then begin
                 tilt = tilt+tilt_value[i]
               endif
             endif
              
             if n_elements(get_tilt) gt 0 then begin
               get_tilt[i] = tilt
             endif
    
             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 || TS07 eq 1 then par_iter = par_array[id,*] else par_iter = ''
             
             if geopack_2008 then begin
                geopack_trace_08,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,TS07=TS07,_extra=_extra,/REFINE,/EQUATOR
             endif else begin
                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,TS07=TS07,_extra=_extra,/REFINE,/EQUATOR
             endelse 
             
             ;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 geopack_2008 then begin
              ; if geopack 2008 is being used, need to convert back to GSM
              geopack_conv_coord_08, tr_temp[*,0], tr_temp[*,1], tr_temp[*,2], x_out_gse, y_out_gse, z_out_gse, /from_gsw, /to_gse
              tr_temp = [[x_out_gse], [y_out_gse], [z_out_gse]]
              ; convert from GSE to GSM
              cotrans, tr_temp, tr_temp, t_temp, /gse2gsm
          endif
          
          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
    
    ; if geopack 2008 is being used, need to convert back to GSM
    if geopack_2008 then begin
        geopack_conv_coord_08, out_foot_array[*,0], out_foot_array[*,1], out_foot_array[*,2], x_footout_gse, y_footout_gse, z_footout_gse, /from_gsw, /to_gse
        out_foot_array = [[x_footout_gse], [y_footout_gse], [z_footout_gse]]
        ; convert from GSE to GSM
        cotrans, out_foot_array, out_foot_array, tarray2, /gse2gsm
    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