;+
; NAME:
;  vty16_temp_terms_bare_iter
; 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_iter(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:
;  temp_t: terms of approximate temperture, K, s.t.
;     temp = temp_t[0] + temp_t[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
;                      + temp_t[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-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_iter, sol_terms,flux_int,emis,freq,therminertia, thermcond, debug=debug, i_loc=i_loc

    _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_termsis 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
    dtemp_dzeta = flux_int/thermcond
    if t_2d then begin
       if size(dtemp_dzeta,/n_dim) eq 0 then dtemp_dzeta = replicate(dtemp_dzeta, n_loc)
       if size(emis,/n_dim) eq 0 then em = replicate(emis, n_loc) else em=emis
    endif

    in = re(sol_0)
    n_t = 12 * (n_term + 1)
    temp_zeta0 = vty16_thermwave(temp_terms, findgen(n_t)*2*!pi/n_t, flux_int/thermcond, 1., 0.)
    if n_loc eq 1 then begin
       out = mean(emis * !phys.sigma * temp_zeta0^4)

       i_iter = 0
       max_iter = 10
       while abs(out-in)/in gt 0.0001 and (i_iter lt max_iter) do begin
          if temp_0 eq 0. then stop
          temp_0 = temp_0 + (in-out)/(4*emis * !phys.sigma * temp_0^3)
          if t_2d then temp_terms[*,0] = temp_0 else temp_terms[0] = temp_0
          if (1) then begin
             for m = 1, n_term do begin
                ; sol_m is the the m-th term of sol_terms
                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
             dtemp_dzeta = flux_int/thermcond
             if t_2d then begin
                if size(dtemp_dzeta,/n_dim) eq 0 then dtemp_dzeta = replicate(dtemp_dzeta, n_loc)
                if size(emis,/n_dim) eq 0 then em = replicate(emis, n_loc) else em=emis
             endif
             
          endif
          temp_zeta0 = vty16_thermwave(temp_terms, findgen(n_t)*2*!pi/n_t, flux_int/thermcond, 1., 0.)
          out = mean(emis * !phys.sigma * temp_zeta0^4)
          i_iter += 1
       endwhile    
    endif else begin
       out = em * !phys.sigma * total(temp_zeta0^4,1)/n_t
       in = re(sol_0)
       indx = where(re(sol_0) gt 0., n_indx)
       for ii = 0, n_indx-1 do begin
          i = indx[ii]
          while abs(out[i]-in[i])/in[i] gt 0.0001 do begin
             if temp_0[i] eq 0. then stop
             temp_0[i] = temp_0[i] + (in[i]-out[i])/(4*em[i] * !phys.sigma * temp_0[i]^3)
             if t_2d then temp_terms[*,0] = temp_0 else temp_terms[i,0] = temp_0[i]
             temp_zeta0_i = vty16_thermwave(temp_terms[i,*], findgen(n_t)*2*!pi/n_t, dtemp_dzeta[i], 1., 0.)
             out[i] = mean(em[i] * !phys.sigma * temp_zeta0_i^4)
             if size(out,/n_dim) eq 0 then stop
          endwhile    
       endfor
    endelse

    return, temp_terms

end