[docs]defdaily_insolation(lat,day,orb=const.orb_present,S0=const.S0,day_type=1):"""Compute daily average insolation given latitude, time of year and orbital parameters. Orbital parameters can be interpolated to any time in the last 5 Myears with ``climlab.solar.orbital.OrbitalTable`` (see example above). Longer orbital tables are available with ``climlab.solar.orbital.LongOrbitalTable`` Inputs can be scalar, ``numpy.ndarray``, or ``xarray.DataArray``. The return value will be ``numpy.ndarray`` if **all** the inputs are ``numpy``. Otherwise ``xarray.DataArray``. **Function-call argument** \n :param array lat: Latitude in degrees (-90 to 90). :param array day: Indicator of time of year. See argument ``day_type`` for details about format. :param dict orb: a dictionary with three members (as provided by ``climlab.solar.orbital.OrbitalTable``) * ``'ecc'`` - eccentricity * unit: dimensionless * default value: ``0.017236`` * ``'long_peri'`` - longitude of perihelion (precession angle) * unit: degrees * default value: ``281.37`` * ``'obliquity'`` - obliquity angle * unit: degrees * default value: ``23.446`` :param float S0: solar constant \n - unit: :math:`\\textrm{W}/\\textrm{m}^2` \n - default value: ``1365.2`` :param int day_type: Convention for specifying time of year (+/- 1,2) [optional]. *day_type=1* (default): day input is calendar day (1-365.24), where day 1 is January first. The calendar is referenced to the vernal equinox which always occurs at day 80. *day_type=2:* day input is solar longitude (0-360 degrees). Solar longitude is the angle of the Earth's orbit measured from spring equinox (21 March). Note that calendar days and solar longitude are not linearly related because, by Kepler's Second Law, Earth's angular velocity varies according to its distance from the sun. :raises: :exc:`ValueError` if day_type is neither 1 nor 2 :returns: Daily average solar radiation in unit :math:`\\textrm{W}/\\textrm{m}^2`. Dimensions of output are ``(lat.size, day.size, ecc.size)`` :rtype: array Code is fully vectorized to handle array input for all arguments. \n Orbital arguments should all have the same sizes. This is automatic if computed from :func:`~climlab.solar.orbital.OrbitalTable.lookup_parameters` For more information about computation of solar insolation see the :ref:`Tutorial` chapter. """# Inputs can be scalar, numpy vector, or xarray.DataArray.# If numpy, convert to xarray so that it will broadcast correctlylat_is_xarray=Trueday_is_xarray=Trueiftype(lat)isnp.ndarray:lat_is_xarray=Falselat=xr.DataArray(lat,coords=[lat],dims=['lat'])iftype(day)isnp.ndarray:day_is_xarray=Falseday=xr.DataArray(day,coords=[day],dims=['day'])ecc=orb['ecc']long_peri=orb['long_peri']obliquity=orb['obliquity']# Convert precession angle and latitude to radiansphi=deg2rad(lat)# lambda_long (solar longitude) is the angular distance along Earth's orbit measured from spring equinox (21 March)ifday_type==1:# calendar dayslambda_long=solar_longitude(day,orb)elifday_type==2:#solar longitude (1-360) is specified in input, no need to convert days to longitudelambda_long=deg2rad(day)else:raiseValueError('Invalid day_type.')# Compute declination angle of the sundelta=arcsin(sin(deg2rad(obliquity))*sin(lambda_long))# suppress warning message generated by arccos here!oldsettings=np.seterr(invalid='ignore')# Compute Ho, the hour angle at sunrise / sunset# Check for no sunrise or no sunset: Berger 1978 eqn (8),(9)Ho=xr.where(abs(delta)-pi/2+abs(phi)<0.,# there is sunset/sunrisearccos(-tan(phi)*tan(delta)),# otherwise figure out if it's all night or all dayxr.where(phi*delta>0.,pi,0.))# this is not really the daily average cosine of the zenith angle...# it's the integral from sunrise to sunset of that quantity...coszen=Ho*sin(phi)*sin(delta)+cos(phi)*cos(delta)*sin(Ho)# Compute insolation: Berger 1978 eq (10)Fsw=S0/pi*((1+ecc*cos(lambda_long-deg2rad(long_peri)))**2/(1-ecc**2)**2*coszen)ifnot(lat_is_xarrayorday_is_xarray):# Dimensional ordering consistent with previous numpy codereturnFsw.transpose().valueselse:returnFsw