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