API tools faq
paste
martineau

Sun.py

martineau
Apr 19th, 2021 (edited)
324
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 19.18 KB | None | 0 0
raw download clone embed print report
  1. # -*- coding: iso-8859-1 -*-
  2. """
  3. SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
  4.             the length of the day at any date and latitude
  5.  
  6. Written as DAYLEN.C, 1989-08-16
  7.  
  8. Modified to SUNRISET.C, 1992-12-01
  9.  
  10. (c) Paul Schlyter, 1989, 1992
  11.  
  12. Released to the public domain by Paul Schlyter, December 1992
  13.  
  14. Direct conversion to Java
  15. Sean Russell <ser@germane-software.com>
  16.  
  17. Conversion to Python Class, 2002-03-21
  18. Henrik Härkönen <radix@kortis.to>
  19.  
  20. Solar Altitude added by Miguel Tremblay 2005-01-16
  21. Solar flux, equation of time and import of python library
  22.  added by Miguel Tremblay 2007-11-22
  23.  
  24.  
  25. 2007-12-12 - v1.5 by Miguel Tremblay: bug fix to solar flux calculation
  26.  
  27.  
  28. """
  29.  
  30. SUN_PY_VERSION = 1.5
  31.  
  32. import math
  33. from math import pi
  34.  
  35. import calendar
  36.  
  37. class Sun:
  38.  
  39.     def __init__(self):
  40.         """"""
  41.  
  42.         # Some conversion factors between radians and degrees
  43.         self.RADEG = 180.0 / pi
  44.         self.DEGRAD = pi / 180.0
  45.         self.INV360 = 1.0 / 360.0
  46.  
  47.  
  48.     def daysSince2000Jan0(self, y, m, d):
  49.         """A macro to compute the number of days elapsed since 2000 Jan 0.0
  50.           (which is equal to 1999 Dec 31, 0h UT)"""
  51.         return (367*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530)
  52.  
  53.  
  54.     # The trigonometric functions in degrees
  55.     def sind(self, x):
  56.         """Returns the sin in degrees"""
  57.         return math.sin(x * self.DEGRAD)
  58.  
  59.     def cosd(self, x):
  60.         """Returns the cos in degrees"""
  61.         return math.cos(x * self.DEGRAD)
  62.  
  63.     def tand(self, x):
  64.         """Returns the tan in degrees"""
  65.         return math.tan(x * self.DEGRAD)
  66.  
  67.     def atand(self, x):
  68.         """Returns the arc tan in degrees"""
  69.         return math.atan(x) * self.RADEG
  70.  
  71.     def asind(self, x):
  72.         """Returns the arc sin in degrees"""
  73.         return math.asin(x) * self.RADEG
  74.  
  75.     def acosd(self, x):
  76.         """Returns the arc cos in degrees"""
  77.         return math.acos(x) * self.RADEG
  78.  
  79.     def atan2d(self, y, x):
  80.         """Returns the atan2 in degrees"""
  81.         return math.atan2(y, x) * self.RADEG
  82.  
  83.     # Following are some macros around the "workhorse" function __daylen__
  84.     # They mainly fill in the desired values for the reference altitude
  85.     # below the horizon, and also selects whether this altitude should
  86.     # refer to the Sun's center or its upper limb.
  87.  
  88.     def dayLength(self, year, month, day, lon, lat):
  89.         """
  90.        This macro computes the length of the day, from sunrise to sunset.
  91.        Sunrise/set is considered to occur when the Sun's upper limb is
  92.        35 arc minutes below the horizon (this accounts for the refraction
  93.        of the Earth's atmosphere).
  94.        """
  95.         return self.__daylen__(year, month, day, lon, lat, -35.0/60.0, 1)
  96.  
  97.     def dayCivilTwilightLength(self, year, month, day, lon, lat):
  98.         """
  99.        This macro computes the length of the day, including civil twilight.
  100.        Civil twilight starts/ends when the Sun's center is 6 degrees below
  101.        the horizon.
  102.        """
  103.         return self.__daylen__(year, month, day, lon, lat, -6.0, 0)
  104.  
  105.     def dayNauticalTwilightLength(self, year, month, day, lon, lat):
  106.         """
  107.        This macro computes the length of the day, incl. nautical twilight.
  108.        Nautical twilight starts/ends when the Sun's center is 12 degrees
  109.        below the horizon.
  110.        """
  111.         return self.__daylen__(year, month, day, lon, lat, -12.0, 0)
  112.  
  113.     def dayAstronomicalTwilightLength(self, year, month, day, lon, lat):
  114.         """
  115.        This macro computes the length of the day, incl. astronomical twilight.
  116.        Astronomical twilight starts/ends when the Sun's center is 18 degrees
  117.        below the horizon.
  118.        """
  119.         return self.__daylen__(year, month, day, lon, lat, -18.0, 0)
  120.  
  121.     def sunRiseSet(self, year, month, day, lon, lat):
  122.         """
  123.        This macro computes times for sunrise/sunset.
  124.        Sunrise/set is considered to occur when the Sun's upper limb is
  125.        35 arc minutes below the horizon (this accounts for the refraction
  126.        of the Earth's atmosphere).
  127.        """
  128.         return self.__sunriset__(year, month, day, lon, lat, -35.0/60.0, 1)
  129.  
  130.  
  131.     def civilTwilight(self, year, month, day, lon, lat):
  132.         """
  133.        This macro computes the start and end times of civil twilight.
  134.        Civil twilight starts/ends when the Sun's center is 6 degrees below
  135.        the horizon.
  136.        """
  137.         return self.__sunriset__(year, month, day, lon, lat, -6.0, 0)
  138.  
  139.     def nauticalTwilight(self, year, month, day, lon, lat):
  140.         """
  141.        This macro computes the start and end times of nautical twilight.
  142.        Nautical twilight starts/ends when the Sun's center is 12 degrees
  143.        below the horizon.
  144.        """
  145.         return self.__sunriset__(year, month, day, lon, lat, -12.0, 0)
  146.  
  147.     def astronomicalTwilight(self, year, month, day, lon, lat):
  148.         """
  149.        This macro computes the start and end times of astronomical twilight.
  150.        Astronomical twilight starts/ends when the Sun's center is 18 degrees
  151.        below the horizon.
  152.        """
  153.         return self.__sunriset__(year, month, day, lon, lat, -18.0, 0)
  154.  
  155.     # The "workhorse" function for sun rise/set times
  156.     def __sunriset__(self, year, month, day, lon, lat, altit, upper_limb):
  157.         """
  158.        Note: year,month,date = calendar date, 1801-2099 only.
  159.              Eastern longitude positive, Western longitude negative
  160.              Northern latitude positive, Southern latitude negative
  161.              The longitude value IS critical in this function!
  162.              altit = the altitude which the Sun should cross
  163.                      Set to -35/60 degrees for rise/set, -6 degrees
  164.                      for civil, -12 degrees for nautical and -18
  165.                      degrees for astronomical twilight.
  166.                upper_limb: non-zero -> upper limb, zero -> center
  167.                      Set to non-zero (e.g. 1) when computing rise/set
  168.                      times, and to zero when computing start/end of
  169.                      twilight.
  170.              *rise = where to store the rise time
  171.              *set  = where to store the set  time
  172.                      Both times are relative to the specified altitude,
  173.                      and thus this function can be used to compute
  174.                      various twilight times, as well as rise/set times
  175.        Return value:  0 = sun rises/sets this day, times stored at
  176.                           *trise and *tset.
  177.                      +1 = sun above the specified 'horizon' 24 hours.
  178.                           *trise set to time when the sun is at south,
  179.                           minus 12 hours while *tset is set to the south
  180.                           time plus 12 hours. 'Day' length = 24 hours
  181.                      -1 = sun is below the specified 'horizon' 24 hours
  182.                           'Day' length = 0 hours, *trise and *tset are
  183.                            both set to the time when the sun is at south.
  184.        """
  185.         # Compute d of 12h local mean solar time
  186.         d = self.daysSince2000Jan0(year,month,day) + 0.5 - (lon/360.0)
  187.  
  188.         # Compute local sidereal time of this moment
  189.         sidtime = self.revolution(self.GMST0(d) + 180.0 + lon)
  190.  
  191.         # Compute Sun's RA + Decl at this moment
  192.         res = self.sunRADec(d)
  193.         sRA = res[0]
  194.         sdec = res[1]
  195.         sr = res[2]
  196.  
  197.         # Compute time when Sun is at south - in hours UT
  198.         tsouth = 12.0 - self.rev180(sidtime - sRA)/15.0;
  199.  
  200.         # Compute the Sun's apparent radius, degrees
  201.         sradius = 0.2666 / sr;
  202.  
  203.         # Do correction to upper limb, if necessary
  204.         if upper_limb:
  205.             altit = altit - sradius
  206.  
  207.         # Compute the diurnal arc that the Sun traverses to reach
  208.         # the specified altitude altit:
  209.  
  210.         cost = (self.sind(altit) - self.sind(lat) * self.sind(sdec))/\
  211.                (self.cosd(lat) * self.cosd(sdec))
  212.  
  213.         if cost >= 1.0:
  214.             rc = -1
  215.             t = 0.0           # Sun always below altit
  216.  
  217.         elif cost <= -1.0:
  218.             rc = +1
  219.             t = 12.0;         # Sun always above altit
  220.  
  221.         else:
  222.             t = self.acosd(cost)/15.0   # The diurnal arc, hours
  223.  
  224.  
  225.         # Store rise and set times - in hours UT
  226.         return (tsouth-t, tsouth+t)
  227.  
  228.  
  229.     def __daylen__(self, year, month, day, lon, lat, altit, upper_limb):
  230.         """
  231.        Note: year,month,date = calendar date, 1801-2099 only.
  232.              Eastern longitude positive, Western longitude negative
  233.              Northern latitude positive, Southern latitude negative
  234.              The longitude value is not critical. Set it to the correct
  235.              longitude if you're picky, otherwise set to, say, 0.0
  236.              The latitude however IS critical - be sure to get it correct
  237.              altit = the altitude which the Sun should cross
  238.                      Set to -35/60 degrees for rise/set, -6 degrees
  239.                      for civil, -12 degrees for nautical and -18
  240.                      degrees for astronomical twilight.
  241.                upper_limb: non-zero -> upper limb, zero -> center
  242.                      Set to non-zero (e.g. 1) when computing day length
  243.                      and to zero when computing day+twilight length.
  244.  
  245.        """
  246.  
  247.         # Compute d of 12h local mean solar time
  248.         d = self.daysSince2000Jan0(year,month,day) + 0.5 - (lon/360.0)
  249.  
  250.         # Compute obliquity of ecliptic (inclination of Earth's axis)
  251.         obl_ecl = 23.4393 - 3.563E-7 * d
  252.  
  253.         # Compute Sun's position
  254.         res = self.sunpos(d)
  255.         slon = res[0]
  256.         sr = res[1]
  257.  
  258.         # Compute sine and cosine of Sun's declination
  259.         sin_sdecl = self.sind(obl_ecl) * self.sind(slon)
  260.         cos_sdecl = math.sqrt(1.0 - sin_sdecl * sin_sdecl)
  261.  
  262.         # Compute the Sun's apparent radius, degrees
  263.         sradius = 0.2666 / sr
  264.  
  265.         # Do correction to upper limb, if necessary
  266.         if upper_limb:
  267.             altit = altit - sradius
  268.  
  269.  
  270.         cost = (self.sind(altit) - self.sind(lat) * sin_sdecl) / \
  271.                (self.cosd(lat) * cos_sdecl)
  272.         if cost >= 1.0:
  273.             t = 0.0             # Sun always below altit
  274.  
  275.         elif cost <= -1.0:
  276.             t = 24.0      # Sun always above altit
  277.  
  278.         else:
  279.             t = (2.0/15.0) * self.acosd(cost);     # The diurnal arc, hours
  280.  
  281.         return t
  282.  
  283.  
  284.     def sunpos(self, d):
  285.         """
  286.        Computes the Sun's ecliptic longitude and distance
  287.        at an instant given in d, number of days since
  288.        2000 Jan 0.0.  The Sun's ecliptic latitude is not
  289.        computed, since it's always very near 0.
  290.        """
  291.  
  292.         # Compute mean elements
  293.         M = self.revolution(356.0470 + 0.9856002585 * d)
  294.         w = 282.9404 + 4.70935E-5 * d
  295.         e = 0.016709 - 1.151E-9 * d
  296.  
  297.         # Compute true longitude and radius vector
  298.         E = M + e * self.RADEG * self.sind(M) * (1.0 + e * self.cosd(M))
  299.         x = self.cosd(E) - e
  300.         y = math.sqrt(1.0 - e*e) * self.sind(E)
  301.         r = math.sqrt(x*x + y*y)              #Solar distance
  302.         v = self.atan2d(y, x)                 # True anomaly
  303.         lon = v + w                        # True solar longitude
  304.         if lon >= 360.0:
  305.             lon = lon - 360.0   # Make it 0..360 degrees
  306.  
  307.         return (lon,r)
  308.  
  309.  
  310.     def sunRADec(self, d):
  311.         """
  312.        Returns the angle of the Sun (RA)
  313.        the declination (dec) and the distance of the Sun (r)
  314.        for a given day d.
  315.        """
  316.  
  317.         # Compute Sun's ecliptical coordinates
  318.         res = self.sunpos(d)
  319.         lon = res[0]  # True solar longitude
  320.         r = res[1]    # Solar distance
  321.  
  322.         # Compute ecliptic rectangular coordinates (z=0)
  323.         x = r * self.cosd(lon)
  324.         y = r * self.sind(lon)
  325.  
  326.         # Compute obliquity of ecliptic (inclination of Earth's axis)
  327.         obl_ecl = 23.4393 - 3.563E-7 * d
  328.  
  329.         # Convert to equatorial rectangular coordinates - x is unchanged
  330.         z = y * self.sind(obl_ecl)
  331.         y = y * self.cosd(obl_ecl)
  332.  
  333.         # Convert to spherical coordinates
  334.         RA = self.atan2d(y, x)
  335.         dec = self.atan2d(z, math.sqrt(x*x + y*y))
  336.  
  337.         return (RA, dec, r)
  338.  
  339.  
  340.     def revolution(self, x):
  341.         """
  342.        This function reduces any angle to within the first revolution
  343.        by subtracting or adding even multiples of 360.0 until the
  344.        result is >= 0.0 and < 360.0
  345.  
  346.        Reduce angle to within 0..360 degrees
  347.        """
  348.         return (x - 360.0 * math.floor(x * self.INV360))
  349.  
  350.     def rev180(self, x):
  351.         """
  352.        Reduce angle to within +180..+180 degrees
  353.        """
  354.         return (x - 360.0 * math.floor(x * self.INV360 + 0.5))
  355.  
  356.     def GMST0(self, d):
  357.         """
  358.        This function computes GMST0, the Greenwich Mean Sidereal Time
  359.        at 0h UT (i.e. the sidereal time at the Greenwhich meridian at
  360.        0h UT).  GMST is then the sidereal time at Greenwich at any
  361.        time of the day.  I've generalized GMST0 as well, and define it
  362.        as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at
  363.        other times than 0h UT as well.  While this sounds somewhat
  364.        contradictory, it is very practical:  instead of computing
  365.        GMST like:
  366.  
  367.         GMST = (GMST0) + UT * (366.2422/365.2422)
  368.  
  369.        where (GMST0) is the GMST last time UT was 0 hours, one simply
  370.        computes:
  371.  
  372.         GMST = GMST0 + UT
  373.  
  374.        where GMST0 is the GMST "at 0h UT" but at the current moment!
  375.        Defined in this way, GMST0 will increase with about 4 min a
  376.        day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)
  377.        is equal to the Sun's mean longitude plus/minus 180 degrees!
  378.        (if we neglect aberration, which amounts to 20 seconds of arc
  379.        or 1.33 seconds of time)
  380.        """
  381.         # Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr
  382.         # L = M + w, as defined in sunpos().  Since I'm too lazy to
  383.         # add these numbers, I'll let the C compiler do it for me.
  384.         # Any decent C compiler will add the constants at compile
  385.         # time, imposing no runtime or code overhead.
  386.  
  387.         sidtim0 = self.revolution((180.0 + 356.0470 + 282.9404) +
  388.                                      (0.9856002585 + 4.70935E-5) * d)
  389.         return sidtim0;
  390.  
  391.     def solar_altitude(self, latitude, year, month, day):
  392.         """
  393.        Compute the altitude of the sun. No atmospherical refraction taken
  394.        in account.
  395.        Altitude of the southern hemisphere are given relative to
  396.        true north.
  397.        Altitude of the northern hemisphere are given relative to
  398.        true south.
  399.        Declination is between 23.5° North and 23.5° South depending
  400.        on the period of the year.
  401.        Source of formula for altitude is PhysicalGeography.net
  402.        http://www.physicalgeography.net/fundamentals/6h.html
  403.        """
  404.         # Compute declination
  405.         N = self.daysSince2000Jan0(year, month, day)
  406.         res =  self.sunRADec(N)
  407.         declination = res[1]
  408.         sr = res[2]
  409.  
  410.         # Compute the altitude
  411.         altitude = 90.0 - latitude  + declination
  412.  
  413.         # In the tropical and  in extreme latitude, values over 90 may occurs.
  414.         if altitude > 90:
  415.             altitude = 90 - (altitude-90)
  416.  
  417.         if altitude < 0:
  418.             altitude = 0
  419.  
  420.         return altitude
  421.  
  422.     def get_max_solar_flux(self, latitude, year, month, day):
  423.         """
  424.        Compute the maximal solar flux to reach the ground for this date and
  425.        latitude.
  426.        Originaly comes from Environment Canada weather forecast model.
  427.        Information was of the public domain before release by Environment Canada
  428.        Output is in W/M^2.
  429.        """
  430.  
  431.         (fEot, fR0r, tDeclsc) = self.equation_of_time(year, month, day, latitude)
  432.         fSF = (tDeclsc[0]+tDeclsc[1])*fR0r
  433.  
  434.         # In the case of a negative declinaison, solar flux is null
  435.         if fSF < 0:
  436.             fCoeff = 0
  437.         else:
  438.             fCoeff =  -1.56e-12*fSF**4 + 5.972e-9*fSF**3 -\
  439.                      8.364e-6*fSF**2  + 5.183e-3*fSF - 0.435
  440.  
  441.         fSFT = fSF * fCoeff
  442.  
  443.         if fSFT < 0:
  444.             fSFT=0
  445.  
  446.         return fSFT
  447.  
  448.     def equation_of_time(self, year, month, day, latitude):
  449.         """
  450.        Description: Subroutine computing the part of the equation of time
  451.                     needed in the computing of the theoritical solar flux
  452.                     Correction originating of the CMC GEM model.
  453.  
  454.        Parameters:  int nTime : cTime for the correction of the time.
  455.  
  456.        Returns: tuple (double fEot, double fR0r, tuple tDeclsc)
  457.                 dEot: Correction for the equation of time
  458.                 dR0r: Corrected solar constant for the equation of time
  459.                 tDeclsc: Declinaison
  460.        """
  461.         # Julian date
  462.         nJulianDate = self.Julian(year, month, day)
  463.         # Check if it is a leap year
  464.         if(calendar.isleap(year)):
  465.             fDivide = 366.0
  466.         else:
  467.             fDivide = 365.0
  468.         # Correction for "equation of time"
  469.         fA = nJulianDate/fDivide*2*pi
  470.         fR0r = self.__Solcons(fA)*0.1367e4
  471.         fRdecl = 0.412*math.cos((nJulianDate+10.0)*2.0*pi/fDivide-pi)
  472.         fDeclsc1 = self.sind(latitude)*math.sin(fRdecl)
  473.         fDeclsc2 = self.cosd(latitude)*math.cos(fRdecl)
  474.         tDeclsc = (fDeclsc1, fDeclsc2)
  475.         # in minutes
  476.         fEot = 0.002733 -7.343*math.sin(fA)+ .5519*math.cos(fA) \
  477.                - 9.47*math.sin(2.0*fA) - 3.02*math.cos(2.0*fA) \
  478.                - 0.3289*math.sin(3.*fA) -0.07581*math.cos(3.0*fA) \
  479.                -0.1935*math.sin(4.0*fA) -0.1245*math.cos(4.0*fA)
  480.         # Express in fraction of hour
  481.         fEot = fEot/60.0
  482.         # Express in radians
  483.         fEot = fEot*15*pi/180.0
  484.  
  485.         return (fEot, fR0r, tDeclsc)
  486.  
  487.     def __Solcons(self, dAlf):
  488.         """
  489.        Name: __Solcons
  490.  
  491.        Parameters: [I] double dAlf : Solar constant to correct the excentricity
  492.  
  493.        Returns: double dVar : Variation of the solar constant
  494.  
  495.        Functions Called: cos, sin
  496.  
  497.        Description:  Statement function that calculates the variation of the
  498.          solar constant as a function of the julian day. (dAlf, in radians)
  499.  
  500.        Notes: Comes from the
  501.  
  502.        Revision History:
  503.        Author          Date            Reason
  504.        Miguel Tremblay      June 30th 2004
  505.        """
  506.  
  507.         dVar = 1.0/(1.0-9.464e-4*math.sin(dAlf)-0.01671*math.cos(dAlf)- \
  508.                     + 1.489e-4*math.cos(2.0*dAlf)-2.917e-5*math.sin(3.0*dAlf)- \
  509.                     + 3.438e-4*math.cos(4.0*dAlf))**2
  510.         return dVar
  511.  
  512.  
  513.     def Julian(self, year, month, day):
  514.         """
  515.        Return julian day.
  516.        """
  517.         if calendar.isleap(year): # Bissextil year, 366 days
  518.             lMonth = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366]
  519.         else: # Normal year, 365 days
  520.             lMonth = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365]
  521.  
  522.         nJulian = lMonth[month-1] + day
  523.         return nJulian
  524.  
  525.  
  526.  
  527. if __name__ == "__main__":
  528.  
  529.     k = Sun()
  530.     print k.get_max_solar_flux(46.2, 2004, 01, 30)
  531. #    print k.sunRiseSet(2002, 3, 22, 25.42, 62.15)
  532.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!