docs: def solar_longitude( day, orb=const.orb_present, days_per_year = None ):
      """Estimates solar longitude from calendar day.
  
      Method is using an approximation from :cite:`Berger_1978` section 3
      (lambda = 0 at spring equinox).
  
      **Function-call arguments** \n
  
      :param array day:           Indicator of time of year.
      :param dict orb:            a dictionary with three members (as provided by
                                  :class:`~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 days_per_year: number of days in a year (optional)
                                  (default: 365.2422)
                                  Reads the length of the year from
                                  :mod:`~climlab.utils.constants` if available.
      :returns:                   solar longitude ``lambda_long``
                                  in dimension``( day.size, ecc.size )``
      :rtype:                     array
  
      Works for both scalar and vector orbital parameters.
  
  
  
      """
      if days_per_year is None:
          days_per_year = const.days_per_year
      ecc = orb['ecc']
      long_peri_rad = deg2rad( orb['long_peri'])
      delta_lambda = (day - 80.) * 2*pi/days_per_year
      beta = sqrt(1 - ecc**2)
      lambda_long_m = -2*((ecc/2 + (ecc**3)/8 ) * (1+beta) * sin(-long_peri_rad) -
          (ecc**2)/4 * (1/2 + beta) * sin(-2*long_peri_rad) + (ecc**3)/8 *
          (1/3 + beta) * sin(-3*long_peri_rad)) + delta_lambda
      lambda_long = ( lambda_long_m + (2*ecc - (ecc**3)/4)*sin(lambda_long_m - long_peri_rad) +
          (5/4)*(ecc**2) * sin(2*(lambda_long_m - long_peri_rad)) + (13/12)*(ecc**3)
          * sin(3*(lambda_long_m - long_peri_rad)) )
      return lambda_long