;+
;Function: t04s
;
;Purpose:  generates an array of model magnetic field vectors from
;          a monotonic time series and an array of 3-d position
;          vectors
;
;Keywords:
;         tarray: N array representing the time series in seconds utc
;         rgsm_array: Nx3 array representing the position series in
;         earth radii GSM
;         The following arguments can either be N length arrays or
;         single values
;         pdyn_array: Solar wind pressure (nanoPascals) 
;         dsti_array: DST index(nanoTeslas)
;         yimf_array: y component of the interplanetary magnetic field
;         zimf_array: z component of the interplanetary magnetic field
;         w1_array:  index represents a time integral over a storm
;         w2_array:  index represents a time integral over a storm
;         w3_array:  index represents a time integral over a storm
;         w4_array:  index represents a time integral over a storm
;         w5_array:  index represents a time integral over a storm
;         w6_array:  index represents a time integral over a storm
;
;         period(optional): the amount of time between recalculations of
;             geodipole tilt in seconds(default: 60)  
;             increase this value to decrease run time
;         
;
;Returns: an Nx3 length array or -1L on failure
;
;Example:
;   mag_array = t04s(time_array,pos_array,pdyn_array,dsti_array,yimf_array,zimf_array,w1_array,w2_array,w3_array,w4_array,w5_array,w6_array)
;   mag_array = t04s(time_array,pos_array,pdyn_array,dsti_array,yimf_array,zimf_array,w1_array,w2_array,w3_array,w4_array,w5_array,w6_array,period=10)
;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. Sums the contribution from the internal field model and the
;  external field model.
;  3. 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.
;  4. Position units are in earth radii, be sure to divide your normal
;  units by 6374 km to convert them.
;  5.Find more documentation on the inner workings of the model,
;    any gotchas, and the meaning of the arguments at:
;    http://modelweb.gsfc.nasa.gov/magnetos/data-based/modeling.html
;    -or-
;    http://dysprosium.jhuapl.edu/idl_geopack/
;  6. Definition of W1-W6 can be found at:
;  N. A. Tsyganenko and M. I. Sitnov, Modeling the dynamics of the
;  inner magnetosphere during strong geomagnetic storms, J. Geophys. 
;  Res., v. 110 (A3), A03208, doi: 10.1029/2004JA010798, 2005
;
; $LastChangedBy: pcruce $
; $LastChangedDate: 2008-01-15 16:13:22 -0800 (Tue, 15 Jan 2008) $
; $LastChangedRevision: 2266 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/idl_socware/tags/tdas_5_00/external/IDL_GEOPACK/t04s/t04s.pro $
;-

function t04s, tarray, rgsm_array,pdyn,dsti,yimf,zimf,w1,w2,w3,w4,w5,w6, period = period

  ;sanity tests, setting defaults

  if igp_test() eq 0 then return, -1L

  if not keyword_set(period) then period = 600

  if period le 0. then begin
    message, /contiune, 'period must be positive'
    return, -1L
  endif

  t_size = size(tarray, /dimensions)

  pdyn_size = size(pdyn, /dimensions)

  dsti_size = size(dsti, /dimensions)

  yimf_size = size(yimf, /dimensions)

  zimf_size = size(zimf, /dimensions)

  w1_size = size(w1,/dimensions)

  w2_size = size(w2,/dimensions)

  w3_size = size(w3,/dimensions)

  w4_size = size(w4,/dimensions)

  w5_size = size(w5,/dimensions)

  w6_size = size(w6,/dimensions)

  r_size = size(rgsm_array, /dimensions)

  if n_elements(t_size) ne 1 then begin
    message, /continue, 'tarray has incorrect dimensions'
    return, -1L
  endif

  if n_elements(pdyn_size) ne 1 then begin
    message, /continue, 'pdyn has incorrect dimensions'
    return, -1L
  endif

  if n_elements(dsti_size) ne 1 then begin
    message, /continue, 'dsti has incorrect dimensions'
    return, -1L
  endif

  if n_elements(yimf_size) ne 1 then begin
    message, /continue, 'yimf has incorrect dimensions'
    return, -1L
  endif

  if n_elements(zimf_size) ne 1 then begin
    message, /continue, 'zimf has incorrect dimensions'
    return, -1L
  endif

  if n_elements(w1_size) ne 1 then begin
    message, /continue, 'w1 has incorrect dimensions'
    return, -1L
  endif

  if n_elements(w2_size) ne 1 then begin
    message, /continue, 'w2 has incorrect dimensions'
    return, -1L
  endif

  if n_elements(w3_size) ne 1 then begin
    message, /continue, 'w3 has incorrect dimensions'
    return, -1L
  endif

  if n_elements(w4_size) ne 1 then begin
    message, /continue, 'w4 has incorrect dimensions'
    return, -1L
  endif

  if n_elements(w5_size) ne 1 then begin
    message, /continue, 'w5 has incorrect dimensions'
    return, -1L
  endif

  if n_elements(w6_size) ne 1 then begin
    message, /continue, 'w6 has incorrect dimensions'
    return, -1L
  endif
  
  if n_elements(r_size) ne 2 || r_size[1] ne 3 then begin
    message, /continue, 'rgsm has incorrect dimensions'
    return, -1L
  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 rgsm_array'
    return, -1L
  endif

  if pdyn_size[0] eq 0 then begin
      pdyn_array = replicate(pdyn,t_size)
  endif else if t_size[0] ne pdyn_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in pdyn_array'
      return, -1L
  endif else pdyn_array = pdyn

  if dsti_size[0] eq 0 then begin
      dsti_array = replicate(dsti,t_size)
  endif else if t_size[0] ne dsti_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in dsti_array'
      return, -1L
  endif else dsti_array = dsti

  if yimf_size[0] eq 0 then begin
      yimf_array = replicate(yimf,t_size)
  endif else if t_size[0] ne yimf_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in yimf_array'
      return, -1L
  endif else yimf_array = yimf

  if zimf_size[0] eq 0 then begin
      zimf_array = replicate(zimf,t_size)
  endif else if t_size[0] ne zimf_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in zimf_array'
      return, -1L
  endif else zimf_array = zimf

  if w1_size[0] eq 0 then begin
      w1_array = replicate(w1,t_size)
  endif else if t_size[0] ne w1_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in w1_array'
      return, -1L
  endif else w1_array = w1

  if w2_size[0] eq 0 then begin
      w2_array = replicate(w2,t_size)
  endif else if t_size[0] ne w2_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in w2_array'
      return, -1L
  endif else w2_array = w2

  if w3_size[0] eq 0 then begin
      w3_array = replicate(w3,t_size)
  endif else if t_size[0] ne w3_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in w3_array'
      return, -1L
  endif else w3_array = w3

  if w4_size[0] eq 0 then begin
      w4_array = replicate(w4,t_size)
  endif else if t_size[0] ne w4_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in w4_array'
      return, -1L
  endif else w4_array = w4

  if w5_size[0] eq 0 then begin
      w5_array = replicate(w5,t_size)
  endif else if t_size[0] ne w5_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in w5_array'
      return, -1L
  endif else w5_array = w5

  if w6_size[0] eq 0 then begin
      w6_array = replicate(w6,t_size)
  endif else if t_size[0] ne w6_size[0] then begin
      message, /continue, 'number of times in tarray does not match number of elements in w6_array'
      return, -1L
  endif else w6_array = w6

;  out_array = dindgen(r_size)

  ;defaults to NaN so it will plot properly in tplot and to prevent
  ;insertion of spurious default dindgen values
  out_array = make_array(r_size, /DOUBLE, VALUE = !VALUES.D_NAN)

  tstart = tarray[0]

  tend = tarray[t_size - 1L]

  i = 0L

  ts = time_struct(tarray)

  ct = (tend-tstart)/period

  period = long(period)

  parmod = dblarr(t_size, 10)

  parmod[*, 0] = pdyn_array
  
  parmod[*, 1] = dsti_array

  parmod[*, 2] = yimf_array

  parmod[*, 3] = zimf_array

  parmod[*, 4] = w1_array

  parmod[*, 5] = w2_array
  
  parmod[*, 6] = w3_array

  parmod[*, 7] = w4_array

  parmod[*, 8] = w5_array

  parmod[*, 9] = w6_array

  while i le ceil(ct) do begin

    ;find indices of points to be input this iteration
    idx1 = where(tarray ge tstart + i*period)
    idx2 = where(tarray le tstart + (i+1)*period)

    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 = rgsm_array[idx, 0]
      rgsm_y = rgsm_array[idx, 1]
      rgsm_z = rgsm_array[idx, 2]

      ;calculate internal contribution
      geopack_igrf_gsm,rgsm_x, rgsm_y, rgsm_z, igrf_bx, igrf_by,igrf_bz 

      ;calculate external contribution
      ;iopt = kp+1
      geopack_ts04, parmod[id, *], rgsm_x, rgsm_y, rgsm_z, t04s_bx, t04s_by, t04s_bz, tilt = tilt

      ;total field
      out_array[idx, 0] = igrf_bx + t04s_bx
      out_array[idx, 1] = igrf_by + t04s_by
      out_array[idx, 2] = igrf_bz + t04s_bz

    endif

    i++

  endwhile

return, out_array

end