;+
; NAME:
;  vty16_solwave
; PURPOSE: (one line)
;  Return solar wave
; DESCRIPTION:
;  Return solar wave
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  sol = vty16_solwave(sol_terms, phase)
; INPUTS:
;  sol_terms = [sol_0, sol_1, ...], erg cm^-2 s^-1 : complex[n_term+1] or complex[n_loc,n_term+1]
;   s.t., at the surface, sol(phase) = sol_0 + sol_1 exp(i*phase) + temp_2 exp(2*i*phase) + ...
;  phase: float or float[n_t], where n_t is the number of time steps. frequency * time, radian
; OUTPUTS:
;  sol: solar flux, erg cm^-2 s^-1. float, float[n_t], float[n_loc], or float[n_t,n_loc]
;
;   sol_t           phase  sol  
;   n_term+1        scalar scaler
;   n_loc,n_term+1  scalar n_loc
;   n_term+1        n_t    n_t
;   n_loc,n_term+1  n_t    n_t,n_loc
;   
; RESTRICTIONS:
; PROCEDURE:
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   Eq 3.2-1
; MODIFICATION HISTORY:
;  Written 2012 Jan 11, by Leslie Young, SwRI
;  2016 Mar 22 LAY. Modified for inclusion in vty16 library.
;-
function vty16_solwave, sol_terms, phase

    _i = complex(0,1)

    ; get the dimensions of sol_t, phase
    sol_0 = extract_last(sol_terms,0,t_dim,n_term) ; scalar or n_loc
    n_term -= 1

    p_dim = size(phase,/n_dim)  ; 0 or 1
    n_t = n_elements(phase) 

    n_loc = n_elements(sol_0)

    ;increase dimension of each into local copy
    if t_dim lt 2 then sol_t = replicate(1.,n_loc)##sol_terms else sol_t=sol_terms
    if p_dim lt 1 then p = [phase] else p=phase

    ; calculate sol at phase = 0
    sol_0 = extract_last(sol_t,0); n_loc

    sol = replicate(1,n_t)#complex(sol_0,0)
    for m = 1,n_term do begin
       sol_m = extract_last(sol_t,m); n_loc
       sol += exp(_i * m * p) # (sol_m)
    endfor

    case 1 of 

       n_t eq 1 and n_loc eq 1: return, re(sol[0])
       n_t eq 1 and n_loc gt 1: return, reform(re(sol),n_loc)

       n_t gt 1 and n_loc eq 1: return, reform(re(sol),n_t)
       n_t gt 1 and n_loc gt 1: return, reform(re(sol),n_t,n_loc)

    endcase

end