1 ### BEGIN GPL LICENSE BLOCK #####
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.
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.
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.
17 # ##### END GPL LICENSE BLOCK #####
20 from bpy
.app
.handlers
import persistent
21 from mathutils
import Euler
23 from math
import degrees
, radians
, pi
25 from .geo
import parse_position
28 ############################################################################
30 # SunClass is used for storing intermediate sun calculations.
32 ############################################################################
68 use_daylight_savings
= False
74 def sun_update(self
, context
):
78 def parse_coordinates(self
, context
):
79 error_message
= "ERROR: Could not parse coordinates"
80 sun_props
= context
.scene
.sun_pos_properties
82 if sun_props
.co_parser
:
83 parsed_co
= parse_position(sun_props
.co_parser
)
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
91 if sun_props
.co_parser
not in {'', error_message
}:
92 sun_props
.co_parser
= ''
95 def sun_handler(scene
):
96 update_time(bpy
.context
)
100 ############################################################################
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.
106 ############################################################################
109 def move_sun(context
):
110 addon_prefs
= context
.preferences
.addons
[__package__
].preferences
111 sun_props
= context
.scene
.sun_pos_properties
113 if sun_props
.usage_mode
== "HDR":
114 nt
= context
.scene
.world
.node_tree
.nodes
115 env_tex
= nt
.get(sun_props
.hdr_texture
)
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
122 sun
.bind
.az_start_env
= env_tex
.texture_mapping
.rotation
.z
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
128 if sun_props
.sun_object
:
129 sun
.theta
= math
.pi
/ 2 - sun_props
.hdr_elevation
130 sun
.phi
= -sun_props
.hdr_azimuth
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
))
137 set_sun_rotations(obj
, rotation_euler
)
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
:
146 north_offset
= degrees(sun_props
.north_offset
)
148 if addon_prefs
.show_rise_set
:
149 calc_sunrise_sunset(rise
=True)
150 calc_sunrise_sunset(rise
=False)
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
)
156 if 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 sky_node
.sun_elevation
= math
.radians(sun
.elevation
)
165 sky_node
.sun_rotation
= math
.radians(sun
.az_north
)
168 if (sun_props
.sun_object
is not None
169 and sun_props
.sun_object
.name
in context
.view_layer
.objects
):
170 obj
= sun_props
.sun_object
171 set_sun_position(obj
, sun_props
.sun_distance
)
172 rotation_euler
= Euler((math
.radians(sun
.elevation
- 90), 0,
173 math
.radians(-sun
.az_north
)))
174 set_sun_rotations(obj
, rotation_euler
)
177 if sun_props
.object_collection
is not None:
178 sun_objects
= sun_props
.object_collection
.objects
179 object_count
= len(sun_objects
)
180 if sun_props
.object_collection_type
== 'DIURNAL':
183 time_increment
= sun_props
.time_spread
/ (object_count
- 1)
184 local_time
= local_time
+ time_increment
* (object_count
- 1)
186 time_increment
= sun_props
.time_spread
188 for obj
in sun_objects
:
189 get_sun_position(local_time
, sun_props
.latitude
,
190 sun_props
.longitude
, north_offset
, zone
,
191 sun_props
.month
, sun_props
.day
,
192 sun_props
.year
, sun_props
.sun_distance
)
193 set_sun_position(obj
, sun_props
.sun_distance
)
194 local_time
-= time_increment
195 obj
.rotation_euler
= (
196 (math
.radians(sun
.elevation
- 90), 0,
197 math
.radians(-sun
.az_north
)))
200 day_increment
= 365 / object_count
201 day
= sun_props
.day_of_year
+ day_increment
* (object_count
- 1)
202 for obj
in sun_objects
:
203 dt
= (datetime
.date(sun_props
.year
, 1, 1) +
204 datetime
.timedelta(day
- 1))
205 get_sun_position(local_time
, sun_props
.latitude
,
206 sun_props
.longitude
, north_offset
, zone
,
207 dt
.month
, dt
.day
, sun_props
.year
,
208 sun_props
.sun_distance
)
209 set_sun_position(obj
, sun_props
.sun_distance
)
211 obj
.rotation_euler
= (
212 (math
.radians(sun
.elevation
- 90), 0,
213 math
.radians(-sun
.az_north
)))
215 def update_time(context
):
216 sun_props
= context
.scene
.sun_pos_properties
218 if sun_props
.use_day_of_year
:
219 dt
= (datetime
.date(sun_props
.year
, 1, 1) +
220 datetime
.timedelta(sun_props
.day_of_year
- 1))
223 sun
.day_of_year
= sun_props
.day_of_year
224 if sun_props
.day
!= dt
.day
:
225 sun_props
.day
= dt
.day
226 if sun_props
.month
!= dt
.month
:
227 sun_props
.month
= dt
.month
229 dt
= datetime
.date(sun_props
.year
, sun_props
.month
, sun_props
.day
)
230 day_of_year
= dt
.timetuple().tm_yday
231 if sun_props
.day_of_year
!= day_of_year
:
232 sun_props
.day_of_year
= day_of_year
233 sun
.day
= sun_props
.day
234 sun
.month
= sun_props
.month
235 sun
.day_of_year
= day_of_year
236 sun
.year
= sun_props
.year
237 sun
.longitude
= sun_props
.longitude
238 sun
.latitude
= sun_props
.latitude
239 sun
.UTC_zone
= sun_props
.UTC_zone
242 def format_time(the_time
, daylight_savings
, longitude
, UTC_zone
=None):
243 if UTC_zone
is not None:
251 mm
= (the_time
- int(the_time
)) * 60
252 ss
= int((mm
- int(mm
)) * 60)
254 return ("%02i:%02i:%02i" % (hh
, mm
, ss
))
257 def format_hms(the_time
):
258 hh
= str(int(the_time
))
259 min = (the_time
- int(the_time
)) * 60
260 sec
= int((min - int(min)) * 60)
261 mm
= "0" + str(int(min)) if min < 10 else str(int(min))
262 ss
= "0" + str(sec
) if sec
< 10 else str(sec
)
264 return (hh
+ ":" + mm
+ ":" + ss
)
267 def format_lat_long(lat_long
, is_latitude
):
268 hh
= str(abs(int(lat_long
)))
269 min = abs((lat_long
- int(lat_long
)) * 60)
270 sec
= abs(int((min - int(min)) * 60))
271 mm
= "0" + str(int(min)) if min < 10 else str(int(min))
272 ss
= "0" + str(sec
) if sec
< 10 else str(sec
)
277 coord_tag
= " N" if lat_long
> 0 else " S"
279 coord_tag
= " E" if lat_long
> 0 else " W"
281 return hh
+ "° " + mm
+ "' " + ss
+ '"' + coord_tag
286 ############################################################################
288 # Calculate the actual position of the sun based on input parameters.
290 # The sun positioning algorithms below are based on the National Oceanic
291 # and Atmospheric Administration's (NOAA) Solar Position Calculator
292 # which rely on calculations of Jean Meeus' book "Astronomical Algorithms."
293 # Use of NOAA data and products are in the public domain and may be used
294 # freely by the public as outlined in their policies at
295 # www.nws.noaa.gov/disclaimer.php
297 # The calculations of this script can be verified with those of NOAA's
298 # using the Azimuth and Solar Elevation displayed in the SunPos_Panel.
299 # NOAA's web site is:
300 # http://www.esrl.noaa.gov/gmd/grad/solcalc
301 ############################################################################
304 def get_sun_position(local_time
, latitude
, longitude
, north_offset
,
305 utc_zone
, month
, day
, year
, distance
):
307 addon_prefs
= bpy
.context
.preferences
.addons
[__package__
].preferences
308 sun_props
= bpy
.context
.scene
.sun_pos_properties
310 longitude
*= -1 # for internal calculations
311 utc_time
= local_time
+ utc_zone
# Set Greenwich Meridian Time
313 if latitude
> 89.93: # Latitude 90 and -90 gives
314 latitude
= radians(89.93) # erroneous results so nudge it
315 elif latitude
< -89.93:
316 latitude
= radians(-89.93)
318 latitude
= radians(latitude
)
320 t
= julian_time_from_y2k(utc_time
, year
, month
, day
)
322 e
= radians(obliquity_correction(t
))
323 L
= apparent_longitude_of_sun(t
)
324 solar_dec
= sun_declination(e
, L
)
325 eqtime
= calc_equation_of_time(t
)
327 time_correction
= (eqtime
- 4 * longitude
) + 60 * utc_zone
328 true_solar_time
= ((utc_time
- utc_zone
) * 60.0 + time_correction
) % 1440
330 hour_angle
= true_solar_time
/ 4.0 - 180.0
331 if hour_angle
< -180.0:
334 csz
= (math
.sin(latitude
) * math
.sin(solar_dec
) +
335 math
.cos(latitude
) * math
.cos(solar_dec
) *
336 math
.cos(radians(hour_angle
)))
342 zenith
= math
.acos(csz
)
344 az_denom
= math
.cos(latitude
) * math
.sin(zenith
)
346 if abs(az_denom
) > 0.001:
347 az_rad
= ((math
.sin(latitude
) *
348 math
.cos(zenith
)) - math
.sin(solar_dec
)) / az_denom
349 if abs(az_rad
) > 1.0:
350 az_rad
= -1.0 if (az_rad
< 0.0) else 1.0
351 azimuth
= 180.0 - degrees(math
.acos(az_rad
))
355 azimuth
= 180.0 if (latitude
> 0.0) else 0.0
358 azimuth
= azimuth
+ 360.0
360 exoatm_elevation
= 90.0 - degrees(zenith
)
362 if sun_props
.use_refraction
:
363 if exoatm_elevation
> 85.0:
364 refraction_correction
= 0.0
366 te
= math
.tan(radians(exoatm_elevation
))
367 if exoatm_elevation
> 5.0:
368 refraction_correction
= (
369 58.1 / te
- 0.07 / (te
** 3) + 0.000086 / (te
** 5))
370 elif (exoatm_elevation
> -0.575):
371 s1
= (-12.79 + exoatm_elevation
* 0.711)
372 s2
= (103.4 + exoatm_elevation
* (s1
))
373 s3
= (-518.2 + exoatm_elevation
* (s2
))
374 refraction_correction
= 1735.0 + exoatm_elevation
* (s3
)
376 refraction_correction
= -20.774 / te
378 refraction_correction
= refraction_correction
/ 3600
379 solar_elevation
= 90.0 - (degrees(zenith
) - refraction_correction
)
382 solar_elevation
= 90.0 - degrees(zenith
)
384 solar_azimuth
= azimuth
385 solar_azimuth
+= north_offset
387 sun
.az_north
= solar_azimuth
389 sun
.theta
= math
.pi
/ 2 - radians(solar_elevation
)
390 sun
.phi
= radians(solar_azimuth
) * -1
391 sun
.azimuth
= azimuth
392 sun
.elevation
= solar_elevation
395 def set_sun_position(obj
, distance
):
396 locX
= math
.sin(sun
.phi
) * math
.sin(-sun
.theta
) * distance
397 locY
= math
.sin(sun
.theta
) * math
.cos(sun
.phi
) * distance
398 locZ
= math
.cos(sun
.theta
) * distance
400 #----------------------------------------------
401 # Update selected object in viewport
402 #----------------------------------------------
403 obj
.location
= locX
, locY
, locZ
406 def set_sun_rotations(obj
, rotation_euler
):
407 rotation_quaternion
= rotation_euler
.to_quaternion()
408 obj
.rotation_quaternion
= rotation_quaternion
410 if obj
.rotation_mode
in {'XZY', 'YXZ', 'YZX', 'ZXY','ZYX'}:
411 obj
.rotation_euler
= rotation_quaternion
.to_euler(obj
.rotation_mode
)
413 obj
.rotation_euler
= rotation_euler
415 rotation_axis_angle
= obj
.rotation_quaternion
.to_axis_angle()
416 obj
.rotation_axis_angle
= (rotation_axis_angle
[1],
417 *rotation_axis_angle
[0])
420 def calc_sunrise_set_UTC(rise
, jd
, latitude
, longitude
):
421 t
= calc_time_julian_cent(jd
)
422 eq_time
= calc_equation_of_time(t
)
423 solar_dec
= calc_sun_declination(t
)
424 hour_angle
= calc_hour_angle_sunrise(latitude
, solar_dec
)
426 hour_angle
= -hour_angle
427 delta
= longitude
+ degrees(hour_angle
)
428 time_UTC
= 720 - (4.0 * delta
) - eq_time
432 def calc_sun_declination(t
):
433 e
= radians(obliquity_correction(t
))
434 L
= apparent_longitude_of_sun(t
)
435 solar_dec
= sun_declination(e
, L
)
439 def calc_hour_angle_sunrise(lat
, solar_dec
):
440 lat_rad
= radians(lat
)
441 HAarg
= (math
.cos(radians(90.833)) /
442 (math
.cos(lat_rad
) * math
.cos(solar_dec
))
443 - math
.tan(lat_rad
) * math
.tan(solar_dec
))
448 HA
= math
.acos(HAarg
)
452 def calc_solar_noon(jd
, longitude
, timezone
, dst
):
453 t
= calc_time_julian_cent(jd
- longitude
/ 360.0)
454 eq_time
= calc_equation_of_time(t
)
455 noon_offset
= 720.0 - (longitude
* 4.0) - eq_time
456 newt
= calc_time_julian_cent(jd
+ noon_offset
/ 1440.0)
457 eq_time
= calc_equation_of_time(newt
)
459 nv
= 780.0 if dst
else 720.0
460 noon_local
= (nv
- (longitude
* 4.0) - eq_time
+ (timezone
* 60.0)) % 1440
461 sun
.solar_noon
.time
= noon_local
/ 60.0
464 def calc_sunrise_sunset(rise
):
467 jd
= get_julian_day(sun
.year
, sun
.month
, sun
.day
)
468 time_UTC
= calc_sunrise_set_UTC(rise
, jd
, sun
.latitude
, sun
.longitude
)
469 new_time_UTC
= calc_sunrise_set_UTC(rise
, jd
+ time_UTC
/ 1440.0,
470 sun
.latitude
, sun
.longitude
)
471 time_local
= new_time_UTC
+ (-zone
* 60.0)
472 tl
= time_local
/ 60.0
473 get_sun_position(tl
, sun
.latitude
, sun
.longitude
, 0.0,
474 zone
, sun
.month
, sun
.day
, sun
.year
,
476 if sun
.use_daylight_savings
:
478 tl
= time_local
/ 60.0
481 sun
.sunrise
.time
= tl
482 sun
.sunrise
.azimuth
= sun
.azimuth
483 sun
.sunrise
.elevation
= sun
.elevation
484 calc_solar_noon(jd
, sun
.longitude
, -zone
, sun
.use_daylight_savings
)
485 get_sun_position(sun
.solar_noon
.time
, sun
.latitude
, sun
.longitude
,
486 0.0, zone
, sun
.month
, sun
.day
, sun
.year
,
488 sun
.solar_noon
.elevation
= sun
.elevation
491 sun
.sunset
.azimuth
= sun
.azimuth
492 sun
.sunset
.elevation
= sun
.elevation
494 ##########################################################################
495 ## Get the elapsed julian time since 1/1/2000 12:00 gmt
496 ## Y2k epoch (1/1/2000 12:00 gmt) is Julian day 2451545.0
497 ##########################################################################
500 def julian_time_from_y2k(utc_time
, year
, month
, day
):
501 century
= 36525.0 # Days in Julian Century
502 epoch
= 2451545.0 # Julian Day for 1/1/2000 12:00 gmt
503 jd
= get_julian_day(year
, month
, day
)
504 return ((jd
+ (utc_time
/ 24)) - epoch
) / century
507 def get_julian_day(year
, month
, day
):
511 A
= math
.floor(year
/ 100)
512 B
= 2 - A
+ math
.floor(A
/ 4.0)
513 jd
= (math
.floor((365.25 * (year
+ 4716.0))) +
514 math
.floor(30.6001 * (month
+ 1)) + day
+ B
- 1524.5)
518 def calc_time_julian_cent(jd
):
519 t
= (jd
- 2451545.0) / 36525.0
523 def sun_declination(e
, L
):
524 return (math
.asin(math
.sin(e
) * math
.sin(L
)))
527 def calc_equation_of_time(t
):
528 epsilon
= obliquity_correction(t
)
529 ml
= radians(mean_longitude_sun(t
))
530 e
= eccentricity_earth_orbit(t
)
531 m
= radians(mean_anomaly_sun(t
))
532 y
= math
.tan(radians(epsilon
) / 2.0)
534 sin2ml
= math
.sin(2.0 * ml
)
535 cos2ml
= math
.cos(2.0 * ml
)
536 sin4ml
= math
.sin(4.0 * ml
)
538 sin2m
= math
.sin(2.0 * m
)
539 etime
= (y
* sin2ml
- 2.0 * e
* sinm
+ 4.0 * e
* y
*
540 sinm
* cos2ml
- 0.5 * y
** 2 * sin4ml
- 1.25 * e
** 2 * sin2m
)
541 return (degrees(etime
) * 4)
544 def obliquity_correction(t
):
545 ec
= obliquity_of_ecliptic(t
)
546 omega
= 125.04 - 1934.136 * t
547 return (ec
+ 0.00256 * math
.cos(radians(omega
)))
550 def obliquity_of_ecliptic(t
):
551 return ((23.0 + 26.0 / 60 + (21.4480 - 46.8150) / 3600 * t
-
552 (0.00059 / 3600) * t
** 2 + (0.001813 / 3600) * t
** 3))
555 def true_longitude_of_sun(t
):
556 return (mean_longitude_sun(t
) + equation_of_sun_center(t
))
559 def calc_sun_apparent_long(t
):
560 o
= true_longitude_of_sun(t
)
561 omega
= 125.04 - 1934.136 * t
562 lamb
= o
- 0.00569 - 0.00478 * math
.sin(radians(omega
))
566 def apparent_longitude_of_sun(t
):
567 return (radians(true_longitude_of_sun(t
) - 0.00569 - 0.00478 *
568 math
.sin(radians(125.04 - 1934.136 * t
))))
571 def mean_longitude_sun(t
):
572 return (280.46646 + 36000.76983 * t
+ 0.0003032 * t
** 2) % 360
575 def equation_of_sun_center(t
):
576 m
= radians(mean_anomaly_sun(t
))
577 c
= ((1.914602 - 0.004817 * t
- 0.000014 * t
** 2) * math
.sin(m
) +
578 (0.019993 - 0.000101 * t
) * math
.sin(m
* 2) +
579 0.000289 * math
.sin(m
* 3))
583 def mean_anomaly_sun(t
):
584 return (357.52911 + t
* (35999.05029 - 0.0001537 * t
))
587 def eccentricity_earth_orbit(t
):
588 return (0.016708634 - 0.000042037 * t
- 0.0000001267 * t
** 2)