sun_position/sun_calc.py

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