sun_position/sun_calc.py

   1 ### BEGIN GPL LICENSE BLOCK #####
   2 #
   3 #  This program is free software; you can redistribute it and/or
   4 #  modify it under the terms of the GNU General Public License
   5 #  as published by the Free Software Foundation; either version 2
   6 #  of the License, or (at your option) any later version.
   7 #
   8 #  This program is distributed in the hope that it will be useful,
   9 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 #  GNU General Public License for more details.
  12 #
  13 #  You should have received a copy of the GNU General Public License
  14 #  along with this program; if not, write to the Free Software Foundation,
  15 #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16 #
  17 # ##### END GPL LICENSE BLOCK #####
  18
  19 import bpy
  20 from bpy.app.handlers import persistent
  21 from mathutils import Euler
  22 import math
  23 from math import degrees, radians, pi
  24 import datetime
  25 from .geo import parse_position
  26
  27
  28 ############################################################################
  29 #
  30 # SunClass is used for storing intermediate sun calculations.
  31 #
  32 ############################################################################
  33
  34 class SunClass:
  35
  36     class TazEl:
  37         time = 0.0
  38         azimuth = 0.0
  39         elevation = 0.0
  40
  41     class CLAMP:
  42         elevation = 0.0
  43         azimuth = 0.0
  44         az_start_sun = 0.0
  45         az_start_env = 0.0
  46
  47     sunrise = TazEl()
  48     sunset = TazEl()
  49     solar_noon = TazEl()
  50     rise_set_ok = False
  51
  52     bind = CLAMP()
  53     bind_to_sun = False
  54
  55     latitude = 0.0
  56     longitude = 0.0
  57     elevation = 0.0
  58     azimuth = 0.0
  59
  60     month = 0
  61     day = 0
  62     year = 0
  63     day_of_year = 0
  64     time = 0.0
  65
  66     UTC_zone = 0
  67     sun_distance = 0.0
  68     use_daylight_savings = False
  69
  70
  71 sun = SunClass()
  72
  73
  74 def sun_update(self, context):
  75     update_time(context)
  76     move_sun(context)
  77
  78 def parse_coordinates(self, context):
  79     error_message = "ERROR: Could not parse coordinates"
  80     sun_props = context.scene.sun_pos_properties
  81
  82     if sun_props.co_parser:
  83         parsed_co = parse_position(sun_props.co_parser)
  84
  85         if parsed_co is not None and len(parsed_co) == 2:
  86             sun_props.latitude, sun_props.longitude = parsed_co
  87         elif sun_props.co_parser != error_message:
  88             sun_props.co_parser = error_message
  89
  90         # Clear prop
  91     if sun_props.co_parser not in {'', error_message}:
  92         sun_props.co_parser = ''
  93
  94 @persistent
  95 def sun_handler(scene):
  96     update_time(bpy.context)
  97     move_sun(bpy.context)
  98
  99
 100 ############################################################################
 101 #
 102 # move_sun() will cycle through all the selected objects
 103 # and call set_sun_position and set_sun_rotations
 104 # to place them in the sky.
 105 #
 106 ############################################################################
 107
 108
 109 def move_sun(context):
 110     addon_prefs = context.preferences.addons[__package__].preferences
 111     sun_props = context.scene.sun_pos_properties
 112
 113     if sun_props.usage_mode == "HDR":
 114         nt = context.scene.world.node_tree.nodes
 115         env_tex = nt.get(sun_props.hdr_texture)
 116
 117         if sun.bind_to_sun != sun_props.bind_to_sun:
 118             # bind_to_sun was just toggled
 119             sun.bind_to_sun = sun_props.bind_to_sun
 120             sun.bind.az_start_sun = sun_props.hdr_azimuth
 121             if env_tex:
 122                 sun.bind.az_start_env = env_tex.texture_mapping.rotation.z
 123
 124         if env_tex and sun_props.bind_to_sun:
 125             az = sun_props.hdr_azimuth - sun.bind.az_start_sun + sun.bind.az_start_env
 126             env_tex.texture_mapping.rotation.z = az
 127
 128         if sun_props.sun_object:
 129             sun.theta = math.pi / 2 - sun_props.hdr_elevation
 130             sun.phi = -sun_props.hdr_azimuth
 131
 132             obj = sun_props.sun_object
 133             set_sun_position(obj, sun_props.sun_distance)
 134             rotation_euler = Euler((sun_props.hdr_elevation - pi/2,
 135                                     0, -sun_props.hdr_azimuth))
 136
 137             set_sun_rotations(obj, rotation_euler)
 138         return
 139
 140     local_time = sun_props.time
 141     zone = -sun_props.UTC_zone
 142     sun.use_daylight_savings = sun_props.use_daylight_savings
 143     if sun.use_daylight_savings:
 144         zone -= 1
 145
 146     north_offset = degrees(sun_props.north_offset)
 147
 148     if addon_prefs.show_rise_set:
 149         calc_sunrise_sunset(rise=True)
 150         calc_sunrise_sunset(rise=False)
 151
 152     get_sun_position(local_time, sun_props.latitude, sun_props.longitude,
 153             north_offset, zone, sun_props.month, sun_props.day, sun_props.year,
 154             sun_props.sun_distance)
 155
 156     if sun_props.use_sky_texture and sun_props.sky_texture:
 157         sky_node = bpy.context.scene.world.node_tree.nodes.get(sun_props.sky_texture)
 158         if sky_node is not None and sky_node.type == "TEX_SKY":
 159             locX = math.sin(sun.phi) * math.sin(-sun.theta)
 160             locY = math.sin(sun.theta) * math.cos(sun.phi)
 161             locZ = math.cos(sun.theta)
 162             sky_node.texture_mapping.rotation.z = 0.0
 163             sky_node.sun_direction = locX, locY, locZ
 164
 165     # Sun object
 166     if ((sun_props.use_sun_object or sun_props.usage_mode == 'HDR')
 167             and sun_props.sun_object
 168             and sun_props.sun_object.name in context.view_layer.objects):
 169         obj = sun_props.sun_object
 170         set_sun_position(obj, sun_props.sun_distance)
 171         rotation_euler = Euler((math.radians(sun.elevation - 90), 0,
 172                                 math.radians(-sun.az_north)))
 173         set_sun_rotations(obj, rotation_euler)
 174
 175     # Sun collection
 176     if (addon_prefs.show_object_collection
 177             and sun_props.use_object_collection
 178             and sun_props.object_collection):
 179         sun_objects = sun_props.object_collection.objects
 180         object_count = len(sun_objects)
 181         if sun_props.object_collection_type == 'DIURNAL':
 182             # Diurnal motion
 183             if object_count > 1:
 184                 time_increment = sun_props.time_spread / (object_count - 1)
 185                 local_time = local_time + time_increment * (object_count - 1)
 186             else:
 187                 time_increment = sun_props.time_spread
 188
 189             for obj in sun_objects:
 190                 get_sun_position(local_time, sun_props.latitude,
 191                                  sun_props.longitude, north_offset, zone,
 192                                  sun_props.month, sun_props.day,
 193                                  sun_props.year, sun_props.sun_distance)
 194                 set_sun_position(obj, sun_props.sun_distance)
 195                 local_time -= time_increment
 196                 obj.rotation_euler = (
 197                     (math.radians(sun.elevation - 90), 0,
 198                      math.radians(-sun.az_north)))
 199         else:
 200             # Analemma
 201             day_increment = 365 / object_count
 202             day = sun_props.day_of_year + day_increment * (object_count - 1)
 203             for obj in sun_objects:
 204                 dt = (datetime.date(sun_props.year, 1, 1) +
 205                       datetime.timedelta(day - 1))
 206                 get_sun_position(local_time, sun_props.latitude,
 207                                  sun_props.longitude, north_offset, zone,
 208                                  dt.month, dt.day, sun_props.year,
 209                                  sun_props.sun_distance)
 210                 set_sun_position(obj, sun_props.sun_distance)
 211                 day -= day_increment
 212                 obj.rotation_euler = (
 213                     (math.radians(sun.elevation - 90), 0,
 214                      math.radians(-sun.az_north)))
 215
 216 def update_time(context):
 217     sun_props = context.scene.sun_pos_properties
 218
 219     if sun_props.use_day_of_year:
 220         dt = (datetime.date(sun_props.year, 1, 1) +
 221               datetime.timedelta(sun_props.day_of_year - 1))
 222         sun.day = dt.day
 223         sun.month = dt.month
 224         sun.day_of_year = sun_props.day_of_year
 225         if sun_props.day != dt.day:
 226             sun_props.day = dt.day
 227         if sun_props.month != dt.month:
 228             sun_props.month = dt.month
 229     else:
 230         dt = datetime.date(sun_props.year, sun_props.month, sun_props.day)
 231         day_of_year = dt.timetuple().tm_yday
 232         if sun_props.day_of_year != day_of_year:
 233             sun_props.day_of_year = day_of_year
 234         sun.day = sun_props.day
 235         sun.month = sun_props.month
 236         sun.day_of_year = day_of_year
 237     sun.year = sun_props.year
 238     sun.longitude = sun_props.longitude
 239     sun.latitude = sun_props.latitude
 240     sun.UTC_zone = sun_props.UTC_zone
 241
 242
 243 def format_time(the_time, daylight_savings, longitude, UTC_zone=None):
 244     if UTC_zone is not None:
 245         if daylight_savings:
 246             UTC_zone += 1
 247         the_time -= UTC_zone
 248
 249     the_time %= 24
 250
 251     hh = int(the_time)
 252     mm = (the_time - int(the_time)) * 60
 253     ss = int((mm - int(mm)) * 60)
 254
 255     return ("%02i:%02i:%02i" % (hh, mm, ss))
 256
 257
 258 def format_hms(the_time):
 259     hh = str(int(the_time))
 260     min = (the_time - int(the_time)) * 60
 261     sec = int((min - int(min)) * 60)
 262     mm = "0" + str(int(min)) if min < 10 else str(int(min))
 263     ss = "0" + str(sec) if sec < 10 else str(sec)
 264
 265     return (hh + ":" + mm + ":" + ss)
 266
 267
 268 def format_lat_long(lat_long, is_latitude):
 269     hh = str(abs(int(lat_long)))
 270     min = abs((lat_long - int(lat_long)) * 60)
 271     sec = abs(int((min - int(min)) * 60))
 272     mm = "0" + str(int(min)) if min < 10 else str(int(min))
 273     ss = "0" + str(sec) if sec < 10 else str(sec)
 274     if lat_long == 0:
 275         coord_tag = " "
 276     else:
 277         if is_latitude:
 278             coord_tag = " N" if lat_long > 0 else " S"
 279         else:
 280             coord_tag = " E" if lat_long > 0 else " W"
 281
 282     return hh + "° " + mm + "' " + ss + '"' + coord_tag
 283
 284
 285
 286
 287 ############################################################################
 288 #
 289 # Calculate the actual position of the sun based on input parameters.
 290 #
 291 # The sun positioning algorithms below are based on the National Oceanic
 292 # and Atmospheric Administration's (NOAA) Solar Position Calculator
 293 # which rely on calculations of Jean Meeus' book "Astronomical Algorithms."
 294 # Use of NOAA data and products are in the public domain and may be used
 295 # freely by the public as outlined in their policies at
 296 #               www.nws.noaa.gov/disclaimer.php
 297 #
 298 # The calculations of this script can be verified with those of NOAA's
 299 # using the Azimuth and Solar Elevation displayed in the SunPos_Panel.
 300 # NOAA's web site is:
 301 #               http://www.esrl.noaa.gov/gmd/grad/solcalc
 302 ############################################################################
 303
 304
 305 def get_sun_position(local_time, latitude, longitude, north_offset,
 306                    utc_zone, month, day, year, distance):
 307
 308     addon_prefs = bpy.context.preferences.addons[__package__].preferences
 309     sun_props = bpy.context.scene.sun_pos_properties
 310
 311     longitude *= -1                 # for internal calculations
 312     utc_time = local_time + utc_zone   # Set Greenwich Meridian Time
 313
 314     if latitude > 89.93:            # Latitude 90 and -90 gives
 315         latitude = radians(89.93)  # erroneous results so nudge it
 316     elif latitude < -89.93:
 317         latitude = radians(-89.93)
 318     else:
 319         latitude = radians(latitude)
 320
 321     t = julian_time_from_y2k(utc_time, year, month, day)
 322
 323     e = radians(obliquity_correction(t))
 324     L = apparent_longitude_of_sun(t)
 325     solar_dec = sun_declination(e, L)
 326     eqtime = calc_equation_of_time(t)
 327
 328     time_correction = (eqtime - 4 * longitude) + 60 * utc_zone
 329     true_solar_time = ((utc_time - utc_zone) * 60.0 + time_correction) % 1440
 330
 331     hour_angle = true_solar_time / 4.0 - 180.0
 332     if hour_angle < -180.0:
 333         hour_angle += 360.0
 334
 335     csz = (math.sin(latitude) * math.sin(solar_dec) +
 336            math.cos(latitude) * math.cos(solar_dec) *
 337            math.cos(radians(hour_angle)))
 338     if csz > 1.0:
 339         csz = 1.0
 340     elif csz < -1.0:
 341         csz = -1.0
 342
 343     zenith = math.acos(csz)
 344
 345     az_denom = math.cos(latitude) * math.sin(zenith)
 346
 347     if abs(az_denom) > 0.001:
 348         az_rad = ((math.sin(latitude) *
 349                   math.cos(zenith)) - math.sin(solar_dec)) / az_denom
 350         if abs(az_rad) > 1.0:
 351             az_rad = -1.0 if (az_rad < 0.0) else 1.0
 352         azimuth = 180.0 - degrees(math.acos(az_rad))
 353         if hour_angle > 0.0:
 354             azimuth = -azimuth
 355     else:
 356         azimuth = 180.0 if (latitude > 0.0) else 0.0
 357
 358     if azimuth < 0.0:
 359         azimuth = azimuth + 360.0
 360
 361     exoatm_elevation = 90.0 - degrees(zenith)
 362
 363     if sun_props.use_refraction:
 364         if exoatm_elevation > 85.0:
 365             refraction_correction = 0.0
 366         else:
 367             te = math.tan(radians(exoatm_elevation))
 368             if exoatm_elevation > 5.0:
 369                 refraction_correction = (
 370                     58.1 / te - 0.07 / (te ** 3) + 0.000086 / (te ** 5))
 371             elif (exoatm_elevation > -0.575):
 372                 s1 = (-12.79 + exoatm_elevation * 0.711)
 373                 s2 = (103.4 + exoatm_elevation * (s1))
 374                 s3 = (-518.2 + exoatm_elevation * (s2))
 375                 refraction_correction = 1735.0 + exoatm_elevation * (s3)
 376             else:
 377                 refraction_correction = -20.774 / te
 378
 379         refraction_correction = refraction_correction / 3600
 380         solar_elevation = 90.0 - (degrees(zenith) - refraction_correction)
 381
 382     else:
 383         solar_elevation = 90.0 - degrees(zenith)
 384
 385     solar_azimuth = azimuth
 386     solar_azimuth += north_offset
 387
 388     sun.az_north = solar_azimuth
 389
 390     sun.theta = math.pi / 2 - radians(solar_elevation)
 391     sun.phi = radians(solar_azimuth) * -1
 392     sun.azimuth = azimuth
 393     sun.elevation = solar_elevation
 394
 395
 396 def set_sun_position(obj, distance):
 397     locX = math.sin(sun.phi) * math.sin(-sun.theta) * distance
 398     locY = math.sin(sun.theta) * math.cos(sun.phi) * distance
 399     locZ = math.cos(sun.theta) * distance
 400
 401     #----------------------------------------------
 402     # Update selected object in viewport
 403     #----------------------------------------------
 404     obj.location = locX, locY, locZ
 405
 406
 407 def set_sun_rotations(obj, rotation_euler):
 408     rotation_quaternion = rotation_euler.to_quaternion()
 409     obj.rotation_quaternion = rotation_quaternion
 410
 411     if obj.rotation_mode in {'XZY', 'YXZ', 'YZX', 'ZXY','ZYX'}:
 412         obj.rotation_euler = rotation_quaternion.to_euler(obj.rotation_mode)
 413     else:
 414         obj.rotation_euler = rotation_euler
 415
 416     rotation_axis_angle = obj.rotation_quaternion.to_axis_angle()
 417     obj.rotation_axis_angle = (rotation_axis_angle[1],
 418                                *rotation_axis_angle[0])
 419
 420
 421 def calc_sunrise_set_UTC(rise, jd, latitude, longitude):
 422     t = calc_time_julian_cent(jd)
 423     eq_time = calc_equation_of_time(t)
 424     solar_dec = calc_sun_declination(t)
 425     hour_angle = calc_hour_angle_sunrise(latitude, solar_dec)
 426     if not rise:
 427         hour_angle = -hour_angle
 428     delta = longitude + degrees(hour_angle)
 429     time_UTC = 720 - (4.0 * delta) - eq_time
 430     return time_UTC
 431
 432
 433 def calc_sun_declination(t):
 434     e = radians(obliquity_correction(t))
 435     L = apparent_longitude_of_sun(t)
 436     solar_dec = sun_declination(e, L)
 437     return solar_dec
 438
 439
 440 def calc_hour_angle_sunrise(lat, solar_dec):
 441     lat_rad = radians(lat)
 442     HAarg = (math.cos(radians(90.833)) /
 443             (math.cos(lat_rad) * math.cos(solar_dec))
 444             - math.tan(lat_rad) * math.tan(solar_dec))
 445     if HAarg < -1.0:
 446         HAarg = -1.0
 447     elif HAarg > 1.0:
 448         HAarg = 1.0
 449     HA = math.acos(HAarg)
 450     return HA
 451
 452
 453 def calc_solar_noon(jd, longitude, timezone, dst):
 454     t = calc_time_julian_cent(jd - longitude / 360.0)
 455     eq_time = calc_equation_of_time(t)
 456     noon_offset = 720.0 - (longitude * 4.0) - eq_time
 457     newt = calc_time_julian_cent(jd + noon_offset / 1440.0)
 458     eq_time = calc_equation_of_time(newt)
 459
 460     nv = 780.0 if dst else 720.0
 461     noon_local = (nv- (longitude * 4.0) - eq_time + (timezone * 60.0)) % 1440
 462     sun.solar_noon.time = noon_local / 60.0
 463
 464
 465 def calc_sunrise_sunset(rise):
 466     zone = -sun.UTC_zone
 467
 468     jd = get_julian_day(sun.year, sun.month, sun.day)
 469     time_UTC = calc_sunrise_set_UTC(rise, jd, sun.latitude, sun.longitude)
 470     new_time_UTC = calc_sunrise_set_UTC(rise, jd + time_UTC / 1440.0,
 471                      sun.latitude, sun.longitude)
 472     time_local = new_time_UTC + (-zone * 60.0)
 473     tl = time_local / 60.0
 474     get_sun_position(tl, sun.latitude, sun.longitude, 0.0,
 475             zone, sun.month, sun.day, sun.year,
 476             sun.sun_distance)
 477     if sun.use_daylight_savings:
 478         time_local += 60.0
 479         tl = time_local / 60.0
 480     tl %= 24.0
 481     if rise:
 482         sun.sunrise.time = tl
 483         sun.sunrise.azimuth = sun.azimuth
 484         sun.sunrise.elevation = sun.elevation
 485         calc_solar_noon(jd, sun.longitude, -zone, sun.use_daylight_savings)
 486         get_sun_position(sun.solar_noon.time, sun.latitude, sun.longitude,
 487             0.0, zone, sun.month, sun.day, sun.year,
 488             sun.sun_distance)
 489         sun.solar_noon.elevation = sun.elevation
 490     else:
 491         sun.sunset.time = tl
 492         sun.sunset.azimuth = sun.azimuth
 493         sun.sunset.elevation = sun.elevation
 494
 495 ##########################################################################
 496 ## Get the elapsed julian time since 1/1/2000 12:00 gmt
 497 ## Y2k epoch (1/1/2000 12:00 gmt) is Julian day 2451545.0
 498 ##########################################################################
 499
 500
 501 def julian_time_from_y2k(utc_time, year, month, day):
 502     century = 36525.0  # Days in Julian Century
 503     epoch = 2451545.0  # Julian Day for 1/1/2000 12:00 gmt
 504     jd = get_julian_day(year, month, day)
 505     return ((jd + (utc_time / 24)) - epoch) / century
 506
 507
 508 def get_julian_day(year, month, day):
 509     if month <= 2:
 510         year -= 1
 511         month += 12
 512     A = math.floor(year / 100)
 513     B = 2 - A + math.floor(A / 4.0)
 514     jd = (math.floor((365.25 * (year + 4716.0))) +
 515          math.floor(30.6001 * (month + 1)) + day + B - 1524.5)
 516     return jd
 517
 518
 519 def calc_time_julian_cent(jd):
 520     t = (jd - 2451545.0) / 36525.0
 521     return t
 522
 523
 524 def sun_declination(e, L):
 525     return (math.asin(math.sin(e) * math.sin(L)))
 526
 527
 528 def calc_equation_of_time(t):
 529     epsilon = obliquity_correction(t)
 530     ml = radians(mean_longitude_sun(t))
 531     e = eccentricity_earth_orbit(t)
 532     m = radians(mean_anomaly_sun(t))
 533     y = math.tan(radians(epsilon) / 2.0)
 534     y = y * y
 535     sin2ml = math.sin(2.0 * ml)
 536     cos2ml = math.cos(2.0 * ml)
 537     sin4ml = math.sin(4.0 * ml)
 538     sinm = math.sin(m)
 539     sin2m = math.sin(2.0 * m)
 540     etime = (y * sin2ml - 2.0 * e * sinm + 4.0 * e * y *
 541              sinm * cos2ml - 0.5 * y ** 2 * sin4ml - 1.25 * e ** 2 * sin2m)
 542     return (degrees(etime) * 4)
 543
 544
 545 def obliquity_correction(t):
 546     ec = obliquity_of_ecliptic(t)
 547     omega = 125.04 - 1934.136 * t
 548     return (ec + 0.00256 * math.cos(radians(omega)))
 549
 550
 551 def obliquity_of_ecliptic(t):
 552     return ((23.0 + 26.0 / 60 + (21.4480 - 46.8150) / 3600 * t -
 553             (0.00059 / 3600) * t ** 2 + (0.001813 / 3600) * t ** 3))
 554
 555
 556 def true_longitude_of_sun(t):
 557     return (mean_longitude_sun(t) + equation_of_sun_center(t))
 558
 559
 560 def calc_sun_apparent_long(t):
 561     o = true_longitude_of_sun(t)
 562     omega = 125.04 - 1934.136 * t
 563     lamb = o - 0.00569 - 0.00478 * math.sin(radians(omega))
 564     return lamb
 565
 566
 567 def apparent_longitude_of_sun(t):
 568     return (radians(true_longitude_of_sun(t) - 0.00569 - 0.00478 *
 569             math.sin(radians(125.04 - 1934.136 * t))))
 570
 571
 572 def mean_longitude_sun(t):
 573     return (280.46646 + 36000.76983 * t + 0.0003032 * t ** 2) % 360
 574
 575
 576 def equation_of_sun_center(t):
 577     m = radians(mean_anomaly_sun(t))
 578     c = ((1.914602 - 0.004817 * t - 0.000014 * t ** 2) * math.sin(m) +
 579         (0.019993 - 0.000101 * t) * math.sin(m * 2) +
 580          0.000289 * math.sin(m * 3))
 581     return c
 582
 583
 584 def mean_anomaly_sun(t):
 585     return (357.52911 + t * (35999.05029 - 0.0001537 * t))
 586
 587
 588 def eccentricity_earth_orbit(t):
 589     return (0.016708634 - 0.000042037 * t - 0.0000001267 * t ** 2)