;+
; NAME:
;  vt3d_eflux_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:
;  flux_terms = vt3d_eflux_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]
;     sol = sol_terms[0] + sol_terms[1]*exp(_i * freq * t) 
;                    + sol_terms[2]*exp(2 *_i freq t) + ..
;     
;  flux_int : internal hear 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:
;  flux_terms: terms of approximate emitted flux, erg cm^-2 s^-1, s.t.
;     flux = flux_terms[0] + flux_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
;                      + flux_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
;  flux_terms is complex[n_term+1] or complex[n_loc,n_term+1]
; OPTIONAL OUTPUTS:
;  temp_terms: terms of approximate temperture, K, s.t.
;     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:
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   Eq 3.2-17
; MODIFICATION HISTORY:
;  Written 2011 Dec 31, by Leslie Young, SwRI
;-
function vt3d_eflux_terms_bare, sol_terms,flux_int,emis,freq,therminertia, temp_terms=temp_terms

    _i = complex(0,1)

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

    ; first compose into flux terms
    flux_0 = (sol_0 + flux_int)    
    temp_0 = vt3d_temp_term0_bare(sol_0, flux_int, emis)
    phi_e = vt3d_dfluxdtemp_emit(emis, temp_0)
    phi_s = vt3d_dfluxdtemp_substrate(freq, therminertia)
    theta_s_over4 =  phi_s/phi_e ; thermal parameter / 4

    ; 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)
    n_phase = 5000L
    phase = findgen(n_phase) * 2.*!pi/n_phase

    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
    if t_2d then flux_terms = complexarr(n_loc,n_term+1) else flux_terms = complexarr(n_term+1)
    if t_2d then flux_terms[*,0] = flux_0 else flux_terms[0] = flux_0
    for m = 1, n_term do begin
       sol_m = extract_last(sol_terms,m)
       if t_2d then begin
          flux_terms[*,m] = sol_m * 1. / (1. + sqrt(_i * m) * theta_s_over4)
       endif else begin
          flux_terms[m] = sol_m * 1./ (1. + sqrt(_i * m) * theta_s_over4)
       endelse
    endfor

    flux = vt3d_solwave(flux_terms, phase)
    temp = ( (flux / (emis * !phys.sigma)) > 0.)^(0.25)
    if t_2d then begin
       for i_loc = 0, n_loc-1 do begin
          temp_terms[i_loc,*] = sft(phase, temp[i_loc,*], n_term)
       endfor
    endif else begin
          temp_terms = sft(phase, temp, n_term)
    endelse

    return, flux_terms

end