;+
; NAME:
;  vty16_solar_mu
; PURPOSE: (one line)
;  Return cosine of the solar incidence angle, mu >= 0
; DESCRIPTION:
;  Return cosine of the solar incidence angle, mu >= 0
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  mu0 = vty16_solar_mu(lat, lon, lat0, lon0)
; INPUTS:
;  lat - scalar or array of latitudes (radian)
;  lon - scalar or array of longitudes (radian)
;  lat0 - scalar or array of sub-solar latitude (radian)
;  lon0 - scalar or array of sub-solar longitude (radian)
; KEYWORD INPUT PARAMETERS:
;  lonavg - return the longitudinal average for the given lat and lat0
;  coslat, sinlat, coslon, sinlon : can pass these if they've
;  been already computed
; OUTPUTS:
;  mu0 = cos(theta0), where theta0 is angle of the sun with respect to the zenith
;  If 
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  lat and lon must be scalars or both 1-D arrays of length n_loc
;  lat0 and lon0 must be scalars or both 1-D arrays of length n_time
; PROCEDURE:
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;  Resubmitted to Icarus.
;  Eq. 3.1-5 (Eq 3.2-4a for longitudinal averaged mu0)
; MODIFICATION HISTORY:
;  Written 2011 March 6, by Leslie Young, SwRI
;  2016 Mar 22 LAY. Modified for inclusion in vty16 library.
;-

function vty16_solar_mu, lat, lon, lat0, lon0, lonavg=lonavg

    coslat = cos(lat)
    sinlat = sin(lat)
    coslon = cos(lon)
    sinlon = sin(lon)

    sinlat0 = sin(lat0)
    coslat0 = cos(lat0)

    n_loc = n_elements(lat)
    n_time = n_elements(lat0)
    mu0 = fltarr(n_loc, n_time)
    ; do an outer product of lat and lat0
    ; so that h_max has the right dimension
    h_max = acos((-1 > (-tan(lat)#tan(lat0)) < 1))
    indx = where(h_max gt 0., nindx)
    if nindx eq 0 then return, mu0
    if keyword_set(lonavg) then begin
       mu0 = ( 0. > ( (sinlat#sinlat0)*h_max + (coslat#coslat0)*sin(h_max)) / !pi )
    endif else begin
       if n_time eq 1 then begin
          mu0[indx] =( sinlat[indx]*sinlat0 + coslat[indx]*coslat0*cos(lon0 - lon[indx]) > 0.)
       endif else begin
          for i_time = 0, n_time-1 do begin
             mu0[indx, i_time] =( sinlat[indx]*sinlat0[i_time] + $
                                  coslat[indx]*coslat0[i_time]*cos(lon0[i_time] - lon[indx]) $
                                  > 0.)
          endfor
       endelse
    endelse
    return, mu0

end