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