sun_position/sun_calc.py

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