;+
;Procedure: tt16
;
;Purpose:  tplot wrapper for the functional interface to the IDL Geopack
;          implementation of the TA16 field model.
;
;Input:
;          pos_gsm_tvar: the tplot variable storing the position in
;              gsm coordinates
;
;
;Keywords:
;         pdyn(optional): Solar wind pressure(nanoPascals) should either be a
;              string naming a tplot variable or an array or a single
;              value. If a tplot input is used it will be interpolated to
;              match the time inputs from the position var. Non-tplot array values
;              must match the number of times in the tplot input for pos_gsm_tvar
;
;         yimf(optional): y component of the interplanetary magnetic field
;             should either be a string naming a tplot variable or an
;             array or a single value.  If a tplot input is used it will
;             be interpolated to match the time inputs from the position
;             var. Non-tplot array values must match the number of times in the
;             tplot input for pos_gsm_tvar
;
;         symh(optional): Sym-H index (nT),
;             longitudinal symmetric component of the ring current;
;             should either be a string naming a tplot variable or an
;             array or a single value.  If a tplot input is used it will
;             be interpolated to match the time inputs from the position
;             var. Non-tplot array values must match the number of times in the
;             tplot input for pos_gsm_tvar
;
;         symc(optional): sliding average of Sym-H over 30-min interval,
;             centered on the current time moment.
;             Either symh or symc must be provided.
;             symc can be computed from symh using GEOPACK_GETSYMHC
;             should either be a string naming a tplot variable or an
;             array or a single value.  If a tplot input is used it will
;             be interpolated to match the time inputs from the position
;             var. Non-tplot array values must match the number of times in the
;             tplot input for pos_gsm_tvar
;
;         xind(optional)
;
;         parmod(optional): can input the Nx10 parmod array used by the
;             fortran Tsyganenko model instead of inputing parameters as
;             separate arrays. If passed as a raw array it will not be
;             modified or interpolated so be sure its has the correct
;             number of entries. It can also be passed as a tplot variable
;             name in which case it will be interpolated. If values are
;             passed individually and as par, the par values will be overwritten.
;
;
;         period(optional): the amount of time between recalculations of
;             geodipole tilt in seconds(default: 60)  increase this
;             value to decrease run time
;
;         get_nperiod(optional): Return the number of periods used in the time interval
;
;         newname(optional):the name of the output variable.
;              (default: pos_gsm_tvar+'_bta16') This option is ignored if
;              globbing is used.
;
;         error(optional): named variable in which to return the
;              error state of this procedure call. 1 = success, 0 = failure
;
;         get_tilt(optional):  Set this value to a tplot variable name in which the geodipole tilt for each period will be returned
;              One sample will be returned for each period with time at the center of the period.
;
;         set_tilt(optional): Set this to a tplot variable name or an array of values containing the dipole tilt that should be used.
;              If a tplot input is used it will be interpolated to match the time inputs from the position
;              var. Non-tplot array values must match the number of times in the tplot input for pos_gsm_tvar
;              Notes:
;                  1) set_tilt will cause add_tilt to be ignored
;                  2) Due to this routine adding IGRF to the returned field, you cannot use set_tilt = 0 and give input
;                      position values in SM coordinates; input position values are required to be in GSM coordinates due to the
;                      IGRF calculation
;
;         add_tilt(optional): Set this to a tplot variable name or an array of values containing the values to be added to the dipole tilt
;              that should be used for each period. If a tplot input is used it will be interpolated to match the time inputs from the position
;              var. Non-tplot array values must match the number of times in the tplot input for pos_gsm_tvar
;
;         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
;
;         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.
;

; Output: Stores the result of the field model calculations in tplot variables
;
; Notes:
;        1. converts from normal gsm to rgsm by dividing vectors by earth's
;            radius(6371.2 km) ie inputs should be in km
;            6371.2 = the value used in the GEOPACK FORTRAN code for Re
;        2. Input must be in GSM coordinates
;        3. Haje Korth's IDL/Geopack DLM must be installed for this
;            procedure to work
;        4. either the variables setting parmod or the variables
;            setting the individual parameter arrays should be set because
;            the defaults aren't scientifically accurate
;        5. model parameters that are input as tplot variables they
;            will be interpolated to match the time values on the input
;            position
;
;  See Newell 2007 for details:
;  https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2006JA012015
;
;  TA16 model description:
;  https://geo.phys.spbu.ru/~tsyganenko/empirical-models/magnetic_field/ta16
;  https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016JA023217
;
;  The N-index calculation is implemented in omni2nindex.pro
;
; $LastChangedBy: jwl $
; $LastChangedDate: 2022-09-16 16:30:39 -0700 (Fri, 16 Sep 2022) $
; $LastChangedRevision: 31098 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_6_1/external/IDL_GEOPACK/ta16/tta16.pro $
;-

pro tta16, pos_gsm_tvar, pdyn=pdyn, yimf=yimf, symh=symh, symc=symc, $
  xind=xind, parmod=parmod, period=period, $
  get_nperiod=get_nperiod, newname=newname, error=error, get_tilt=get_tilt, $
  set_tilt=set_tilt, add_tilt=add_tilt, geopack_2008=geopack_2008, $
  exact_tilt_times=exact_tilt_times, skip_ta16_load=skip_ta16_load

  COMPILE_OPT idl2
  error = 0
  if ta16_supported() eq 0 then return

  if not keyword_set(pos_gsm_tvar) then begin
    message, /continue, 'pos_gsm_tvar must be set'
    return
  endif

  var_names = tnames(pos_gsm_tvar)

  if(var_names[0] eq '') then begin
    message, /continue, 'No valid tplot_variables match pos_gsm_tvar'
    return
  endif else if(n_elements(var_names) ne 1) then begin
    message, /continue, 'tta16 only accepts one position tplot variable as input'
    return
  end

  var_name = var_names[0];just in case

  if n_elements(parmod) gt 0 then begin
    ; parmod: (1) Pdyn [nPa], (2) SymHc, (3) XIND, (4) IMF By [nT].
    if size(parmod,/type) eq 7 then begin
      if tnames(parmod) eq '' then message,'parmod variable not valid tplot variable'
      tinterpol_mxn,parmod,var_name,newname='par_out'
      get_data,'par_out',data=dat
      par_temp = dat.y
    endif else par_temp = parmod

    if n_elements(pdyn) gt 0 then begin
      pdyn_dat = tsy_valid_param(pdyn, var_name)
      if(size(pdyn_dat, /n_dim) eq 0 && pdyn_dat[0] eq -1L) then return
    endif else begin
      pdyn_dat = par_temp[*,0]
    endelse

    if n_elements(yimf) gt 0 then begin
      yimf_dat = tsy_valid_param(yimf, var_name)
      if(size(yimf_dat, /n_dim) eq 0 && yimf_dat[0] eq -1L) then return
    endif else begin
      yimf_dat = par_temp[*,3]
    endelse

    if n_elements(xind) gt 0 then begin
      xind_dat = tsy_valid_param(xind, var_name)
      if(size(xind_dat, /n_dim) eq 0 && xind_dat[0] eq -1L) then return
    endif else begin
      xind_dat = par_temp[*,2]
    endelse

    ; If symc is provided, use it, otherwise compute it from symh
    if n_elements(symc) gt 0 then begin
      symc_dat = tsy_valid_param(symc, var_name)
      if(size(symc_dat, /n_dim) eq 0 && symc_dat[0] eq -1L) then return
    endif else if n_elements(symh) gt 0 then begin
      symc=symh+'_c'
      symh2symc, symh=symh, pdyn=pdyn, trange=trange, newname=symc
      symc_dat = tsy_valid_param(symc, var_name)
      if(size(symc_dat, /n_dim) eq 0 && symc_dat[0] eq -1L) then return
    endif else begin
      symc_dat = par_temp[*,1]
    endelse

  endif else begin
    ; the user didn't provide a parameter array (parmod), need to check
    ; the individual parameter keywords
    if undefined(pdyn) || undefined(yimf) || (undefined(symc) && undefined(symh) ) || undefined(xind) then begin
      dprint, dlevel = 1, 'Error, missing one or more of the required model parameters'
      return
    endif
    pdyn_dat = tsy_valid_param(pdyn, var_name)
    if(size(pdyn_dat, /n_dim) eq 0 && pdyn_dat[0] eq -1L) then return

    yimf_dat = tsy_valid_param(yimf, var_name)
    if(size(yimf_dat, /n_dim) eq 0 && yimf_dat[0] eq -1L) then return

    xind_dat = tsy_valid_param(xind, var_name)
    if(size(xind_dat, /n_dim) eq 0 && xind_dat[0] eq -1L) then return

    ; If symc is provided, use it, otherwise compute it from symh
    if n_elements(symc) gt 0 then begin
      symc_dat = tsy_valid_param(symc, var_name)
    endif else if n_elements(symh) gt 0 then begin
      symc=symh+'_c'
      symh2symc, symh=symh, pdyn=pdyn, trange=trange, newname=symc
      symc_dat = tsy_valid_param(symc, var_name)
      if(size(symc_dat, /n_dim) eq 0 && symc_dat[0] eq -1L) then return
    endif
    if(size(symc_dat, /n_dim) eq 0 && symc_dat[0] eq -1L) then return

  endelse

  get_data, var_name, data = d, dlimits = dl, limits = l

  ; check the variable's coordinate system
  if tsy_valid_coords(dl, geopack_2008 = geopack_2008) eq -1 then return

  if n_elements(add_tilt) gt 0 then begin
    add_tilt_dat = tsy_valid_param(add_tilt, var_name)
    if(size(add_tilt, /n_dim) eq 0 && add_tilt_dat[0] eq -1L) then return
  endif

  if n_elements(set_tilt) gt 0 then begin
    set_tilt_dat = tsy_valid_param(set_tilt, var_name)
    if(size(set_tilt, /n_dim) eq 0 && set_tilt_dat[0] eq -1L) then return
  endif

  ;do the calculation, division converts position into earth radii units
  if n_elements(set_tilt) gt 0 then begin
    mag_array = ta16(d.x, d.y/6371.2, pdyn_dat, yimf_dat, symc_dat, xind_dat, $
      period=period, get_nperiod=get_nperiod, $
      get_period_times=period_times_dat, get_tilt=tilt_dat, set_tilt=set_tilt_dat, $
      geopack_2008=geopack_2008, exact_tilt_times=exact_tilt_times,skip_ta16_load=skip_ta16_load)
  endif else if n_elements(add_tilt) gt 0 then begin
    mag_array = ta16(d.x, d.y/6371.2, pdyn_dat, yimf_dat, symc_dat, xind_dat, $
      period=period, get_nperiod=get_nperiod, $
      get_period_times=period_times_dat, get_tilt=tilt_dat, add_tilt=add_tilt_dat, $
      exact_tilt_times=exact_tilt_times,skip_ta16_load=skip_ta16_load)
  endif else begin
    mag_array = ta16(d.x, d.y/6371.2, pdyn_dat, yimf_dat, symc_dat, xind_dat, $
      period=period, get_nperiod=get_nperiod, $
      get_period_times=period_times_dat, get_tilt=tilt_dat,geopack_2008=geopack_2008, $
      exact_tilt_times=exact_tilt_times, skip_ta16_load=skip_ta16_load)
  endelse

  if size(mag_array, /n_dim) eq 0 && mag_array[0] eq -1L then begin
    message, /continue, 'Tsyganenko model query failed, returning'
    return
  endif

  if is_string(get_tilt) then begin
    store_data,get_tilt,data={x:period_times_dat,y:tilt_dat}
  endif

  ;sometimes v element is present, sometimes not
  ;if it is around it is stored in output so information is not lost
  str_element, d, 'v', success = s

  if s eq 1 then $
    d_out = {x:d.x, y:mag_array, v:d.v} $
  else $
    d_out = {x:d.x, y:mag_array}

  if ~keyword_set(geopack_2008) then geopack_2008 = 0
  pd = 'External magnetic field in nT'
  desc = 'B Field using Geopack model TA16' 
  g_att = {input_var: pos_gsm_tvar, input_coord_sys: 'gsm', input_units: 'km', geopack_2008: geopack_2008}
  data_att = {project:'GEOPACK_DLM', observatory:'B Field', instrument:'TA16', units:'nT', coord_sys:'gsm', description:pd}
  dlg = {geopack:g_att, data_att: data_att, spec:0, log:0, colors: [2,4,6], labels: ['b_x','b_y','b_z'], ysubtitle: '[nT]', description: desc}

  if keyword_set(newname) then begin
    store_data, newname, data = d_out, dlimits = dlg
  endif else begin
    store_data, var_name +'_bta16', data = d_out, dlimits = dlg
  endelse

  ;signal success
  error = 1

  return
end