;+
; NAME:
;  insol_fac
; PURPOSE: (one line)
;  return the insolation factor: insolation / normal insolation
; DESCRIPTION:
;
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  fac = insolfac(sslat, lat, dlon) or insolfac(sslat, lat)
; INPUTS:
;  sslat = sub solar latitude in radians
;  lat = latitude in radians
; OPTIONAL INPUT PARAMETERS:
;  dlon = longitude difference in radians.
;     This cann be thought of as either local time of day
;     with noon = 0, or as the difference between the subsolar
;     longitude and the surface element longitude
;  if dlon is not provided, then this returns the average over a rotation
; OUTPUTS:
;  cos(solar angle) or 
;  1/(2 pi) * integral_over_daylight( cos(solar angle) dlon)
; PROCEDURE:
;  normal to sun is
;    s = [cos(sslat),        0,                 sin(sslat)]
;  normal to point on surface is
;    p = [cos(lat) cos(dlon), cos(lat) sin(dlon), sin(lat) ]
;  and so the dot product is
;    mu = cos(sslat) cos(lat) cos(dlon) + sin(sslat) sin(lat)
; MODIFICATION HISTORY:
;  Written 2006 Dec 26, Leslie Young SwRI
;-

function insolfac, sslat, lat, dlon

    if n_params() eq 3 then begin

        ; insolation at a single point

        mu = cos(sslat) * cos(lat) * cos(dlon) + sin(sslat) * sin(lat)
        return, 0 > mu
    endif else begin

        ; insolation averaged over a day

        ; find the limit of integration

        ; it's easier to make sure we're not on the pole by
        ; adding a hair to the angles than it is to check
        ; division by zero, since this way we're general and
        ; can take lists or scalars for sslat or lat
        mach = machar()
        maxlat = !pi/2. * (1 - mach.eps)
        sins = sin(sslat)
        coss = cos((-maxlat) > sslat < maxlat)
        sinp = sin(lat)
        cosp = cos((-maxlat) > lat < maxlat)

        coslonlim = -1.d > (- (sins * sinp) / (coss * cosp)) < 1.d
        lonlim = acos(coslonlim)
        mu = (coss * cosp * sin(lonlim) + sins * sinp * lonlim)/!dpi

        return, mu

    endelse

end