Azimuth calculations

README.md

@@ -45,7 +45,7 @@ using SkyDomes
 lat = 52.0*π/180.0 # latitude in radians
 DOY = 182
 f = 0.5 # solar noon
-Ig, Idir, Idif = clear_sky(lat = lat, DOY = DOY, f = f) # W/m2
+Ig, Idir, Idif, theta, phi = clear_sky(lat = lat, DOY = DOY, f = f) # W/m2
 ```
 
 The values `Ig`, `Idir` and `Idif` are the total, direct and diffuse solar
@@ -95,12 +95,15 @@ that purpose:
 sources = sky(scene,
              Idir = Idir_PAR, # Direct solar radiation from above
              nrays_dir = 1_000_000, # Number of rays for direct solar radiation
+             theta_dir = theta, # Solar zenith angle for direct solar radiation
+             phi_dir = phi, # Solar azimuth angle for direct solar radiation
              Idif = Idif_PAR, # Diffuse solar radiation from above
              nrays_dif = 10_000_000, # Total number of rays for diffuse solar radiation
              sky_model = StandardSky, # Angular distribution of solar radiation
              dome_method = equal_solid_angles, # Discretization of the SkyDomes dome
              ntheta = 9, # Number of discretization steps in the zenith angle
              nphi = 12) # Number of discretization steps in the azimuth angle
+render!(sources) # To visualize the dome of light sources
 ```
 
 The function takes the scene as input to ensure that light sources scale with

src/SolarIrradiance.jl

@@ -13,12 +13,16 @@ This code is based on:
 # Calculate the declination angle
 declination(DOY) = 23.45 * π / 180.0 * sin(2π / 365 * (DOY - 81))
 
-# Calculate the solar zenith angle and its cosine at a given location and time.
-function solar_zenith_angle(; lat, t, dec)
+# Calculate the solar zenith and azimuth angle and its cosine at a given location and time.
+function solar_angles(; lat, t, dec)
     h = (t - 12.0) * 15.0 * π / 180.0 # hour angle of the sun
     cos_theta = cos(lat) * cos(dec) * cos(h) + sin(lat) * sin(dec) # cosine of zenith angle
     theta = acos(cos_theta) # zenith angle
-    return (cos_theta, theta) # return both values as a tuple
+    # cosine of azimuth angle for each half of the day
+    cos_phi = (sin(dec)*cos(lat) - cos(h)*cos(dec)*sin(lat))/sin(theta)
+    phi = acos(clamp(cos_phi, -1, 1)) # azimuth angle without correction
+    phi = h < 0 ? phi : 2π - phi # compass azimuth angle with correction
+    return (cos_theta, theta, phi) # return both values as a tuple
 end
 
 #=
@@ -72,6 +76,7 @@ A named tuple with fields:
 - `Idir`: direct solar radiation on the horizontal plane in W/m^2
 - `Idif`: diffuse solar radiation on the horizontal plane in W/m^2
 - `theta`: solar zenith angle in radians
+- `phi`: solar azimuth angle in radians
 
 # References
 Ineichen P., Perez R., A new airmass independent formulation for
@@ -87,7 +92,7 @@ function clear_sky(; lat, DOY, f, altitude = 0.0, TL = 4.0)
     dec = declination(DOY) # declination angle of the sun
     DL = day_length(lat, dec)
     t = timeOfDay(f, DL)
-    cos_theta, theta = solar_zenith_angle(; lat = lat, dec = dec, t = t)
+    cos_theta, theta, phi = solar_angles(; lat = lat, dec = dec, t = t)
     # Clear sky model by Ineichen and Perez (2002)
     am = air_mass(cos_theta, theta)
     fh1 = exp(-altitude / 8000)
@@ -99,7 +104,7 @@ function clear_sky(; lat, DOY, f, altitude = 0.0, TL = 4.0)
     # Direct solar radiation on the horizontal plane
     b = 0.664 + 0.163 / fh1
     Idir = Io * cos_theta * b * exp(-0.09 * am * (TL - 1))
-    return (Ig = Ig, Idir = Idir, Idif = Ig - Idir)
+    return (Ig = Ig, Idir = Idir, Idif = Ig - Idir, theta = theta, phi = phi)
 end
 
 # Convert solar irradiance to a particular waveband based on a reference solar

test/test_raytracing.jl

@@ -20,7 +20,7 @@ let
     lat = 52.0 * π / 180.0
     DOY = 182
     f = 0.5
-    Ig, Idir, Idif = SkyDomes.clear_sky(lat = lat, DOY = DOY, f = f)
+    Ig, Idir, Idif, theta = SkyDomes.clear_sky(lat = lat, DOY = DOY, f = f)
 
     # 2. Convert to different wavebands (PAR, NIR, UV, red, green, blue)
     Idir_PAR = SkyDomes.waveband_conversion(Itype = :direct,
@@ -34,7 +34,7 @@ let
 
     # 3. Create a sky dome of directional light sources using the solar irradiance
     sources = SkyDomes.sky(scene, Idir = (Idir_PAR, Idir_PAR), Idif = (Idif_PAR, Idif_PAR),
-        nrays_dif = 10_000_000, nrays_dir = 1_000_000,
+        nrays_dif = 10_000_000, nrays_dir = 1_000_000, theta_dir = theta,
         sky_model = StandardSky,
         dome_method = equal_solid_angles, ntheta = 9, nphi = 12)
 

test/test_solar_irradiance.jl

@@ -9,13 +9,15 @@ let
     DOYs = rand(1:365, n)
     decs = SkyDomes.declination.(DOYs)
     ts = 24.0 .* rand(n)
-    temp = [SkyDomes.solar_zenith_angle(lat = lats[i], dec = decs[i], t = ts[i])
+    temp = [SkyDomes.solar_angles(lat = lats[i], dec = decs[i], t = ts[i])
             for i in 1:n]
-    cos_theta, theta = Tuple(getindex.(temp, i) for i in 1:2)
+    cos_theta, theta, phi = Tuple(getindex.(temp, i) for i in 1:3)
     @test maximum(theta) <= π
     @test minimum(theta) >= 0
     @test maximum(cos_theta) <= 1.0
     @test minimum(cos_theta) >= -1.0
+    @test maximum(phi) <= 2π
+    @test minimum(phi) >= 0
 
     # Compute air mass for the different cos_theta and theta from above
     thetas = 0.0:0.01:(π / 2)
@@ -46,4 +48,17 @@ let
             end
         end
     end
+
+    # Check that zenith angles are calculated correctly
+    lat = 0.0*pi/180
+    DOY = 182
+    dec = SkyDomes.declination(DOY)
+    ts = 0.1:0.1:24.0
+    temp = [SkyDomes.solar_angles(lat = lat, dec = dec, t = t)
+            for t in ts]
+    cos_thetas, thetas, phis = Tuple(getindex.(temp, i) for i in 1:3)
+    day = cos_thetas .> 0
+    # plot(ts[day], theta[day].*180.0./pi)
+    # plot!(ts[day], phis[day].*180.0./pi)
+
 end