;+
; NAME:
;  vty16_temp_terms_bare
; PURPOSE: (one line)
;  Return temperature in balance with mean absorbed solar flux
; DESCRIPTION:
;  Return temperature in balance with mean absorbed solar flux
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  temp_terms = vty16_temp_terms_bare(sol_terms,flux_int,emis,freq,therminertia)
; INPUTS:
;  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
;   complex[n_term+1] or complex[n_loc,n_term+1]
;     defined so that (if _i = sqrt(-1) )
;     sol = sol_terms[0] + sol_terms[1]*exp(   _i * freq * t) 
;                        + sol_terms[2]*exp(2 *_i * freq * t) + ..
;     
;  flux_int : internal heat flux, erg/cm^2/s: float or float[n_loc]
;  emis : thermal emissivity, unitless : float or float[n_loc]
;  freq : frequency of forcing
;  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
; OUTPUTS:
;  temp_terms: terms of approximate temperture, K, such that
;     temp = temp_terms[0] + temp_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
;                      + temp_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
;  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
; RESTRICTIONS:
; PROCEDURE:
; PROCEDURE:
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   Eq 3.2-10
; MODIFICATION HISTORY:
;  Written 2011 Dec 31, by Leslie Young, SwRI
;  2016 Mar 24 LAY. Modified for inclusion in vty16 library.
;-
function vty16_temp_terms_bare, sol_terms,flux_int,emis,freq,therminertia

    _i = complex(0,1)

    sol_0 = extract_last(sol_terms,0,s_dim,n_term)
    n_term -= 1

    temp_0 = vty16_temp_term0_bare(sol_0, flux_int, emis)
    phi_e = vty16_dfluxdtemp_emit(emis, temp_0)
    phi_s = vty16_dfluxdtemp_substrate(freq, therminertia)
    ; get dimensionality
    foo = temp_0 * phi_s
    t_2d = (size(foo, /n_dim) gt 0) or (s_dim eq 2) ; 1 if temp_terms is a 2-D array
    n_loc = n_elements(foo)

    if t_2d then temp_terms = complexarr(n_loc,n_term+1) else temp_terms = complexarr(n_term+1)
    if t_2d then temp_terms[*,0] = temp_0 else temp_terms[0] = temp_0
    for m = 1, n_term do begin
       sol_m = extract_last(sol_terms,m)
       if t_2d then begin
          temp_terms[*,m] = sol_m / (phi_e + sqrt(_i * m) * phi_s)
       endif else begin
          temp_terms[m] = sol_m / (phi_e + sqrt(_i * m) * phi_s)
       endelse
    endfor

    return, temp_terms

end