;+ ; NAME: ; vt3d_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 = vt3d_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 vt3d library. ;- function vt3d_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