diff options
| -rw-r--r-- | src/solar.py | 320 |
1 files changed, 213 insertions, 107 deletions
diff --git a/src/solar.py b/src/solar.py index ead1653..b14ed47 100644 --- a/src/solar.py +++ b/src/solar.py @@ -23,42 +23,50 @@ import time SOLAR_ELEVATION_SUNSET_SUNRISE = 0.0 ''' -:float The Sun's elevation at sunset and sunrise, measured in degrees +:float The Sun's elevation at sunset and sunrise, + measured in degrees ''' SOLAR_ELEVATION_CIVIL_DUSK_DAWN = -6.0 ''' -:float The Sun's elevation at civil dusk and civil dawn, measured in degrees +:float The Sun's elevation at civil dusk and civil + dawn, measured in degrees ''' SOLAR_ELEVATION_NAUTICAL_DUSK_DAWN = -12.0 ''' -:float The Sun's elevation at nautical dusk and nautical dawn, measured in degrees +:float The Sun's elevation at nautical dusk and + nautical dawn, measured in degrees ''' SOLAR_ELEVATION_ASTRONOMICAL_DUSK_DAWN = -18.0 ''' -:float The Sun's elevation at astronomical dusk and astronomical dawn, measured in degrees +:float The Sun's elevation at astronomical dusk + and astronomical dawn, measured in degrees ''' SOLAR_ELEVATION_RANGE_TWILIGHT = (-18.0, 0.0) ''' -:(float, float) The Sun's lowest and highest elevation during all periods of twilight, measured in degrees +:(float, float) The Sun's lowest and highest elevation during + all periods of twilight, measured in degrees ''' SOLAR_ELEVATION_RANGE_CIVIL_TWILIGHT = (-6.0, 0.0) ''' -:(float, float) The Sun's lowest and highest elevation during civil twilight, measured in degrees +:(float, float) The Sun's lowest and highest elevation + during civil twilight, measured in degrees ''' SOLAR_ELEVATION_RANGE_NAUTICAL_TWILIGHT = (-12.0, -6.0) ''' -:(float, float) The Sun's lowest and highest elevation during nautical twilight, measured in degrees +:(float, float) The Sun's lowest and highest elevation + during nautical twilight, measured in degrees ''' SOLAR_ELEVATION_RANGE_ASTRONOMICAL_TWILIGHT = (-18.0, -12.0) ''' -:(float, float) The Sun's lowest and highest elevation during astronomical twilight, measured in degrees +:(float, float) The Sun's lowest and highest elevation during + astronomical twilight, measured in degrees ''' @@ -73,7 +81,7 @@ def sun(latitude, longitude, t = None, low = -6.0, high = 3.0): @param low:float The 100 % night limit elevation of the Sun (highest when not visible) @param high:float The 100 % day limit elevation of the Sun (lowest while fully visible) @return :float The visibilty of the Sun, 0 during the night, 1 during the day, - between 0 and 1 during twilight. Other values will not occur. + between 0 and 1 during twilight. Other values will not occur ''' t = julian_centuries() if t is None else t e = solar_elevation(latitude, longitude, t) @@ -85,7 +93,7 @@ def sun(latitude, longitude, t = None, low = -6.0, high = 3.0): # The following functions are used to calculate the result for `sun` # (most of them) but could be used for anything else. There name is # should tell you enough, `t` (and `noon`) is in Julian centuries -# except for in the convertion methods +# except for in the convertion methods. def julian_day_to_epoch(t): @@ -279,7 +287,8 @@ def sun_apparent_longitude(t): def mean_ecliptic_obliquity(t): ''' - Calculates the mean ecliptic obliquity of the Sun's apparent motion without variation correction + Calculates the mean ecliptic obliquity of the Sun's + apparent motion without variation correction @param t:float The time in Julian Centuries @return :float The uncorrected mean obliquity, in radians @@ -292,7 +301,8 @@ def mean_ecliptic_obliquity(t): def corrected_mean_ecliptic_obliquity(t): ''' - Calculates the mean ecliptic obliquity of the Sun's apparent motion with variation correction + Calculates the mean ecliptic obliquity of the Sun's + apparent motion with variation correction @param t:float The time in Julian Centuries @return :float The mean obliquity, in radians @@ -317,7 +327,8 @@ def solar_declination(t): def equation_of_time(t): ''' - Calculates the equation of time, the discrepancy between apparent and mean solar time + Calculates the equation of time, the discrepancy + between apparent and mean solar time @param t:float The time in Julian Centuries @return :float The equation of time, in degrees @@ -338,7 +349,8 @@ def hour_angle_from_elevation(latitude, declination, elevation): ''' Calculates the solar hour angle from the Sun's elevation - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param longitude:float The longitude in degrees eastwards + from Greenwich, negative for westwards @param declination:float The declination, in degrees @param hour_angle:float The Sun's elevation, in degrees @return :float The solar hour angle, in degrees @@ -356,7 +368,8 @@ def elevation_from_hour_angle(latitude, declination, hour_angle): ''' Calculates the Sun's elevation from the solar hour angle - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param longitude:float The longitude in degrees eastwards + from Greenwich, negative for westwards @param declination:float The declination, in degrees @param hour_angle:float The solar hour angle, in degrees @return :float The Sun's elevation, in degrees @@ -372,7 +385,8 @@ def time_of_solar_noon(t, longitude): Calculates the time of the closest solar noon @param t:float A time close to the seeked time, in Julian Centuries - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards @return :float The time, in Julian Centuries, of the closest solar noon ''' t, rc = julian_centuries_to_julian_day(t), longitude @@ -384,12 +398,15 @@ def time_of_solar_noon(t, longitude): def time_of_solar_elevation(t, noon, latitude, longitude, elevation): ''' - Calculates the time the Sun has a specified apparent elevation at a geographical position + Calculates the time the Sun has a specified apparent + elevation at a geographical position @param t:float A time close to the seeked time, in Julian Centuries @param noon:float The time of the closest solar noon - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards @param elevation:float The solar elevation, in degrees @return :float The time, in Julian Centuries, of the specified elevation ''' @@ -407,12 +424,16 @@ def time_of_solar_elevation(t, noon, latitude, longitude, elevation): def solar_elevation_from_time(t, latitude, longitude): ''' - Calculates the Sun's elevation as apparent from a geographical position + Calculates the Sun's elevation as apparent + from a geographical position @param t:float The time in Julian Centuries - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards - @return :float The Sun's apparent at the specified time as seen from the specified position, + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards + @return :float The Sun's apparent at the specified time + as seen from the specified position, measured in degrees ''' rc = julian_centuries_to_julian_day(t) @@ -424,12 +445,17 @@ def solar_elevation_from_time(t, latitude, longitude): def solar_elevation(latitude, longitude, t = None): ''' - Calculates the Sun's elevation as apparent from a geographical position - - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float The Sun's apparent at the specified time as seen from the specified position, + Calculates the Sun's elevation as apparent + from a geographical position + + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float The Sun's apparent at the specified time + as seen from the specified position, measured in degrees ''' rc = julian_centuries() if t is None else t @@ -441,12 +467,16 @@ def solar_elevation(latitude, longitude, t = None): # TODO document def have_sunrise_and_sunset(latitude, t = None): ''' - Determine whether solar declination currently is so that there can be sunrises and sunsets. - If not, you either have 24-hour daytime or 24-hour nighttime. + Determine whether solar declination currently is + so that there can be sunrises and sunsets. If not, + you either have 24-hour daytime or 24-hour nighttime. - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param t:float? The time in Julian Centuries, `None` for the current time - @return Whether there can be sunrises and sunsets where you are located + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param t:float? The time in Julian Centuries, `None` + for the current time + @return Whether there can be sunrises and + sunsets where you are located ''' t = julian_centuries() if t is None else t d = degrees(solar_declination(t)) @@ -465,9 +495,12 @@ def is_summer(latitude, t = None): ''' Determine whether it is summer - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param t:float? The time in Julian Centuries, `None` for the current time - @return Whether it is summer on the hemisphere you are located on + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param t:float? The time in Julian Centuries, `None` + for the current time + @return Whether it is summer on the hemisphere + you are located on ''' t = julian_centuries() if t is None else t d = solar_declination(t) @@ -479,9 +512,12 @@ def is_winter(latitude, t = None): ''' Determine whether it is winter - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param t:float? The time in Julian Centuries, `None` for the current time - @return Whether it is winter on the hemisphere you are located on + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param t:float? The time in Julian Centuries, `None` + for the current time + @return Whether it is winter on the hemisphere + you are located on ''' t = julian_centuries() if t is None else t d = solar_declination(t) @@ -492,15 +528,21 @@ def is_winter(latitude, t = None): # TODO document def solar_prediction(delta, requested, fun, epsilon = 0.000001, span = 0.01, t = None): ''' - Predict the time point of the next or previous time an arbitrary condition is meet + Predict the time point of the next or previous + time an arbitrary condition is meet - @param delta:float Iteration step size, negative for past event, positive for future event - @param requested:float The value returned by `fun` for which to calculate the time point of occurrence + @param delta:float Iteration step size, negative for past + event, positive for future event + @param requested:float The value returned by `fun` for which to + calculate the time point of occurrence @param fun:(t:float)→float Function that calculate the data of interest @param epsilon:float Error tolerance for `requested` - @param span:float The number of Julian centuries (0,01 for one year) to restrict the search to - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within the specified time span + @param span:float The number of Julian centuries (0,01 for + one year) to restrict the search to + @param t:float? The time in Julian Centuries, `None` for + the current time + @return :float? The calculated time point, `None` if none + were found within the specified time span ''' t = julian_centuries() if t is None else t t1 = t2 = t @@ -542,8 +584,10 @@ def future_past_equinox(delta, t = None): ''' Predict the time point of the next or previous equinox - @param delta:float Iteration step size, negative for past event, positive for future event - @param t:float? The time in Julian Centuries, `None` for the current time + @param delta:float Iteration step size, negative for + past event, positive for future event + @param t:float? The time in Julian Centuries, `None` + for the current time @return :float The calculated time point ''' return solar_prediction(delta, 0, solar_declination, t = t) @@ -554,8 +598,10 @@ def future_equinox(t = None): ''' Predict the time point of the next equinox - @param delta:float Iteration step size, negative for past event, positive for future event - @param t:float? The time in Julian Centuries, `None` for the current time + @param delta:float Iteration step size, negative for + past event, positive for future event + @param t:float? The time in Julian Centuries, `None` + for the current time @return :float The calculated time point ''' return future_past_equinox(0.01 / 2000, t) @@ -566,8 +612,10 @@ def past_equinox(t = None): ''' Predict the time point of the previous equinox - @param delta:float Iteration step size, negative for past event, positive for future event - @param t:float? The time in Julian Centuries, `None` for the current time + @param delta:float Iteration step size, negative for + past event, positive for future event + @param t:float? The time in Julian Centuries, `None` + for the current time @return :float The calculated time point ''' return future_past_equinox(0.01 / -2000, t) @@ -579,13 +627,16 @@ def future_past_solstice(delta, t = None): ''' Predict the time point of the next or previous solstice - @param delta:float Iteration step size, negative for past event, positive for future event - @param t:float? The time in Julian Centuries, `None` for the current time + @param delta:float Iteration step size, negative for + past event, positive for future event + @param t:float? The time in Julian Centuries, `None` + for the current time @return :float The calculated time point ''' e = 0.00001 - fun = lambda t : (solar_declination(t + e) - solar_declination(t - e)) / 2 - return solar_prediction(delta, 0, fun, t = t) + fun = solar_declination + dfun = lambda t : (fun(t + e) - fun(t - e)) / 2 + return solar_prediction(delta, 0, dfun, t = t) # TODO document @@ -593,7 +644,8 @@ def future_solstice(t = None): ''' Predict the time point of the next solstice - @param t:float? The time in Julian Centuries, `None` for the current time + @param t:float? The time in Julian Centuries, + `None` for the current time @return :float The calculated time point ''' return future_past_solstice(0.01 / 2000, t) @@ -604,7 +656,8 @@ def past_solstice(t = None): ''' Predict the time point of the previous solstice - @param t:float? The time in Julian Centuries, `None` for the current time + @param t:float? The time in Julian Centuries, + `None` for the current time @return :float The calculated time point ''' return future_past_solstice(0.01 / -2000, t) @@ -613,111 +666,164 @@ def past_solstice(t = None): def future_past_elevation(delta, latitude, longitude, elevation, t = None): ''' - Predict the time point of the next or previous time the Sun reaches or reached a specific elevation + Predict the time point of the next or previous time + the Sun reaches or reached a specific elevation - @param delta:float Iteration step size, negative for past event, positive for future event - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param delta:float Iteration step size, negative for past + event, positive for future event + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards @param elevation:float The elevation of interest - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within a year + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float? The calculated time point, `None` if + none were found within a year ''' - return solar_prediction(delta, elevation, lambda t : solar_elevation(latitude, longitude, t), t = t) + fun = lambda t : solar_elevation(latitude, longitude, t) + return solar_prediction(delta, elevation, fun, t = t) def future_elevation(latitude, longitude, elevation, t = None): ''' - Predict the time point of the next time the Sun reaches a specific elevation + Predict the time point of the next time the Sun + reaches a specific elevation - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards @param elevation:float The elevation of interest - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within a year + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float? The calculated time point, `None` if + none were found within a year ''' return future_past_elevation(0.01 / 2000, latitude, longitude, elevation, t) def past_elevation(latitude, longitude, elevation, t = None): ''' - Predict the time point of the previous time the Sun reached a specific elevation + Predict the time point of the previous time the Sun + reached a specific elevation - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards @param elevation:float The elevation of interest - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within a year + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float? The calculated time point, `None` if + none were found within a year ''' return future_past_elevation(0.01 / -2000, latitude, longitude, elevation, t) # TODO document -def future_past_elevation_derivative(delta, latitude, longitude, elevation_derivative, t = None): - ''' - Predict the time point of the next or previous time the Sun reaches or reached a specific elevation derivative - - @param delta:float Iteration step size, negative for past event, positive for future event - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards - @param elevation_derivative:float The elevation derivative value of interest - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within a year +def future_past_elevation_derivative(delta, latitude, longitude, derivative, t = None): + ''' + Predict the time point of the next or previous time the + Sun reaches or reached a specific elevation derivative + + @param delta:float Iteration step size, negative for past + event, positive for future event + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards + @param derivative:float The elevation derivative value of interest + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float? The calculated time point, `None` if + none were found within a year ''' fun = lambda t : solar_elevation(latitude, longitude, t) dfun = lambda t : (fun(t + e) - fun(t - e)) / 2 - return solar_prediction(delta, elevation_derivative, dfun, t = t) + return solar_prediction(delta, derivative, dfun, t = t) # TODO document -def future_elevation_derivative(latitude, longitude, elevation_derivative, t = None): +def future_elevation_derivative(latitude, longitude, derivative, t = None): ''' - Predict the time point of the next time the Sun reaches a specific elevation derivative + Predict the time point of the next time the + Sun reaches a specific elevation derivative - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards - @param elevation_derivative:float The elevation derivative value of interest - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within a year + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards + @param derivative:float The elevation derivative value of interest + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float? The calculated time point, `None` if + none were found within a year ''' - return future_past_elevation_derivative(0.01 / 2000, latitude, longitude, elevation_derivative, t) + return future_past_elevation_derivative(0.01 / 2000, latitude, longitude, derivative, t) # TODO document -def past_elevation_derivative(latitude, longitude, elevation_derivative, t = None): +def past_elevation_derivative(latitude, longitude, derivative, t = None): ''' - Predict the time point of the previous time the Sun reached a specific elevation derivative + Predict the time point of the previous time + the Sun reached a specific elevation derivative - @param latitude:float The latitude in degrees northwards from the equator, negative for southwards - @param longitude:float The longitude in degrees eastwards from Greenwich, negative for westwards - @param elevation_derivative:float The elevation derivative value of interest - @param t:float? The time in Julian Centuries, `None` for the current time - @return :float? The calculated time point, `None` if none were found within a year + @param latitude:float The latitude in degrees northwards from + the equator, negative for southwards + @param longitude:float The longitude in degrees eastwards from + Greenwich, negative for westwards + @param derivative:float The elevation derivative value of interest + @param t:float? The time in Julian Centuries, `None` + for the current time + @return :float? The calculated time point, `None` + if none were found within a year ''' - return future_past_elevation_derivative(0.01 / -2000, latitude, longitude, elevation_derivative, t) + return future_past_elevation_derivative(0.01 / -2000, latitude, longitude, derivative, t) # TODO: This algorithm is imprecise, gives an incorrent sunrise and I do not fully know its behaviour def sunrise_equation(latitude, longitude, t = None): + # Calculate Julian Cycle j_cent = julian_centuries() if t is None else t j_date = julian_centuries_to_julian_day(j_cent) j_cycle = int(j_date - 2451545.0009 - longitude / 360 + 0.5) + + # Calculate approximate solar noon and solar man anomaly approx_solar_noon = 451545.0009 + longitude / 360 + j_cycle solar_mean_anomaly = int(357.5291 + 0.98560028 * (j_cycle - 2451545)) % 360 - equation_of_centre = 1.9148 * sin(1 * solar_mean_anomaly) + + # Calculate solar equation of centre + equation_of_centre = 1.9148 * sin(1 * solar_mean_anomaly) equation_of_centre += 0.0200 * sin(2 * solar_mean_anomaly) equation_of_centre += 0.0003 * sin(3 * solar_mean_anomaly) + + # Calculate solar ecliptic longitude ecliptic_longitude = (solar_mean_anomaly + 102.9372 + equation_of_centre + 180) % 360 - solar_transit = approx_solar_noon + 0.0053 * sin(solar_mean_anomaly) - 0.0069 * sin(2 * ecliptic_longitude) + + # Calculate solar transit + solar_transit = approx_solar_noon + 0.0053 * sin(solar_mean_anomaly) + solar_transit -= 0.0069 * sin(2 * ecliptic_longitude) + + # Calculate solar declination declination = asin(sin(ecliptic_longitude) * sin(radians(23.45))) - hour_angle = sin(radians(-0.83)) + + # Calculate solar hour angle + hour_angle = sin(radians(-0.83)) hour_angle -= sin(latitude) * sin(declination) hour_angle /= cos(latitude) * cos(declination) hour_angle = degrees(acos(hour_angle)) - sunset = 2451545.0009 + (hour_angle + longitude) / 360 + j_cycle + solar_transit - approx_solar_noon + + # Calculate time of sunset and sunrise + sunset = 2451545.0009 + (hour_angle + longitude) / 360 + sunset += j_cycle + solar_transit - approx_solar_noon sunrise = 2 * solar_transit - sunset - return (julian_day_to_julian_centuries(sunset), julian_day_to_julian_centuries(sunrise)) + + # Convert to Julian Centuries + return (julian_day_to_julian_centuries(sunset), + julian_day_to_julian_centuries(sunrise)) |