Module: Solar

Defined in:
lib/solar.rb,
lib/solar/lambert.rb,
lib/solar/support.rb,
lib/solar/version.rb,
lib/solar/passages.rb,
lib/solar/position.rb,
lib/solar/day_night.rb,
lib/solar/radiation.rb

Overview

Calculation of solar position, rise & set times for a given position & time. Algorithms are taken from Jean Meeus, Astronomical Algorithms Some code & ideas taken from John P. Power’s astro-algo: astro-algo.rubyforge.org/astro-algo/

Constant Summary collapse

ALTITUDES = 
{
  official:     -50/60.0,
  civil:            -6.0,
  nautical:        -12.0,
  astronomical:    -18.0
}
VERSION = 
"0.1.2"
DEBUG_RADIATION = 
ENV['DEBUG_SOLAR_RADIATION']

Class Method Summary collapse

Class Method Details

.day_or_night(t, longitude, latitude, options = {}) ⇒ Object

Day-night (or twilight) status at a given position and time. Returns :night, :day or :twilight

Options:

  • :twilight_zenith is the zenith for the sun at dawn (at the beginning of twilight) and at dusk (end of twilight). Default value is :civil

  • :day_zenith is the zenith for the san at sunrise and sun set. Default value is :official (sun aparently under the horizon, tangent to it)

These parameters can be assigned zenith values in degrees or these symbols: :official, :civil, :nautical or :astronomical.

A simple day/night result (returning either :day or :night) can be requested by setting the :simple option to true (which usses the official day definition), or by setting a :zenith parameter to define the kind of day-night distinction.

By setting the :detailed option to true, the result will be one of: :night, :astronomical_twilight, :nautical_twilight, :civil_twilight, :day



28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 
# File 'lib/solar/day_night.rb', line 28

def day_or_night(t, longitude, latitude, options = {})
  h, _ = position(t, longitude, latitude)
  options = { zenith: :official } if options[:simple]
  if options[:detailed]
    if h < ALTITUDES[:astronomical]
      :night
    elsif  h < ALTITUDES[:nautical]
      :astronomical_twilight
    elsif h < ALTITUDES[:civil]
      :nautical_twilight
    elsif h < ALTITUDES[:official]
      :civil_twilight
    else
      :day
    end
  else
    # Determined :night / :twilight / :day state;
    # twilight/day definition can be changed with options :zenith or :twilight_zenith, :day_zenith
    if options[:zenith]
      # only :day / :night distinction as defined by :zenith
      twilight_altitude = day_altitude = altitude_from_options(options)
    else
      twilight_altitude = altitude_from_options(
        zenith: options[:twilight_zenith] || :civil
      )
      day_altitude = altitude_from_options(
        zenith: options[:day_zenith] || :official
      )
    end
    if h > day_altitude
      :day
    else
      h <= twilight_altitude ? :night : :twilight
    end
  end
end

.illumination_factor(sun_azimuth, sun_elevation, slope, aspect) ⇒ Object

slope, sun_azimuth: degrees (0: horizontal; 90: vertical) aspect: degrees from North towards East sun_elevation: degrees (0 on horizon, 90 on zenith)



7
8
9
10
11
12
13
14
 
# File 'lib/solar/lambert.rb', line 7

def illumination_factor(sun_azimuth, sun_elevation, slope, aspect)
  s = sun_vector(sun_azimuth, sun_elevation)
  nh = vertical_vector
  n  = normal_from_slope_aspect(slope, aspect)
  # Problem: when sun near horizon...
  f = dot(s, n) / dot(s, nh)
  f < 0 ? 0.0 : f
end

.illumination_factor_at(t, longitude, latitude, slope, aspect) ⇒ Object 



16
17
18
19
20
 
# File 'lib/solar/lambert.rb', line 16

def illumination_factor_at(t, longitude, latitude, slope, aspect)
  sun_elevation, sun_azimuth = Solar.position(t, longitude, latitude)
  return 0.0 if sun_elevation < 0 # should return 1.0, really
  illumination_factor(sun_azimuth, sun_elevation, slope, aspect)
end

.mean_illumination_factor_on(date, latitude, slope, aspect, options = {}) ⇒ Object 



22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 
# File 'lib/solar/lambert.rb', line 22

def mean_illumination_factor_on(date, latitude, slope, aspect, options = {})
  t1, t2, t3 = Solar.passages(date, 0, latitude, altitude: 24.0)
  if options[:noon]
    return factor_at(t2, 0.0, latitude, slope, aspect)
  end
  if n = options[:n]
    dt = (t3 - t1).to_f*24*3600/n
  else
    dt = (options[:dt] || 60*30.0).to_f
  end
  f = 0.0
  n = 0
  # max_f = 0.0
  (t1..t3).step(dt).each do |t|
    f += illumination_factor_at(Time.at(t), 0.0, latitude, slope, aspect)
    n += 1
  end
  f / n
end

.passages(date, longitude, latitude, options = {}) ⇒ Object

Solar passages [rising, transit, setting] for a given date (UTC) and position.

The :zenith or :altitude of the sun can be passed as an argument, which can be numeric (in degrees) or symbolic: :official, :civil, :nautical or :astronomical.

In circumpolar case:

  • If Sun never rises, returns 00:00:00 on Date for all passages.

  • If Sun never sets, returns 00:00:00 (rising), 12:00:00 (transit), 24:00:00 (setting) on Date for all passages.



58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
 
# File 'lib/solar/passages.rb', line 58

def passages(date, longitude, latitude, options = {})
  ho = altitude_from_options(options)
  t = to_jc(jd_r(date.to_datetime))
  theta0 = (100.46061837 + (36000.770053608 + (0.000387933 - t/38710000)*t)*t) % 360

  # Calculate apparent right ascention and declination for 0 hr Dynamical time for three days (degrees)
  ra = []
  decl = []
  -1.upto(1) do |i|
    declination, right_ascention = equatorial_position_rad((date+i).to_datetime)
    ra << to_deg(right_ascention)
    decl << to_deg(declination)
  end
  # tweak right ascention around 180 degrees (autumnal equinox)
  if ra[0] > ra[1]
    ra[0] -= 360
  end
  if ra[2] < ra[1]
    ra[2] += 360
  end

  ho_rad, latitude_rad = [ho, latitude].map{|x| to_rad(x)}
  decl_rad = decl.map { |x| to_rad(x) }

  # approximate Hour Angle (degrees)
  ha = Math.sin(ho_rad) /
       Math.cos(latitude_rad) * Math.cos(decl_rad[1]) -
       Math.tan(latitude_rad) * Math.tan(decl_rad[1])
  # handle circumpolar. see note 2 at end of chapter
  if ha.abs <= 1
    ha = to_deg(Math.acos(ha))
  elsif ha > 1  # circumpolar - sun never rises
    # format sunrise, sunset & solar noon as DateTime
    sunrise = date.to_datetime
    transit = date.to_datetime
    sunset = date.to_datetime
    return [sunrise, transit, sunset]
  else  # cirumpolar - sun never sets
    # format sunrise, sunset & solar noon as DateTime
    sunrise = date.to_datetime
    transit = date.to_datetime + 0.5
    sunset = date.to_datetime + 1
    return [sunrise, transit, sunset]
  end
  # approximate m (fraction of 1 day)
  # store days added or subtracted to add in later
  m = []
  days = [0]*3
  (0..2).each do |i|
    case i
    when 0
      m[i] = (ra[1] - longitude - theta0) / 360 # transit
      day_offset = +1
    when 1
      m[i] = m[0] - ha / 360 # rising
      day_offset = -1
    when 2
      m[i] = m[0] + ha / 360 # setting
      day_offset = -1
    end

    until m[i] >= 0
      m[i] += 1
      days[i] += day_offset
    end
    until m[i] <= 1
      m[i] -= 1
      days[i] -= day_offset
    end
  end
  theta = [] # apparent sidereal time (degrees)
  ra2 = []   # apparent right ascension (degrees)
  decl2 = [] # apparent declination (degrees)
  h = []   # local hour angle (degrees)
  alt = []   # altitude (degrees)
  delta_m = [1]*3
  while delta_m[0] >= 0.01 || delta_m[1] >= 0.01 || delta_m[2] >= 0.01
    0.upto(2) do |i|
      theta[i] = theta0 + 360.985647 * m[i]
      n = m[i] + delta_t(date.to_datetime).to_r / 86_400
      a = ra[1] - ra[0]
      b = ra[2] - ra[1]
      c = b - a
      ra2[i] = ra[1] + n / 2 * (a + b + n * c)

      n = m[i] + delta_t(date.to_datetime).to_r / 86_400
      a = decl[1] - decl[0]
      b = decl[2] - decl[1]
      c = b - a
      decl2[i] = decl[1] + n / 2 * (a + b + n * c)

      h[i] = theta[i] + longitude - ra2[i]

      alt[i] = to_deg Math.asin(
        Math.sin(latitude_rad) * Math.sin(to_rad(decl2[i])) +
        Math.cos(latitude_rad) * Math.cos(to_rad(decl2[i])) *
        Math.cos(to_rad(h[i]))
      )
    end
    # adjust m
    delta_m[0] = -h[0] / 360
    1.upto(2) do |i|
      delta_m[i] = (alt[i] - ho) /
      360 * Math.cos(to_rad(decl2[i])) * Math.cos(latitude_rad) * Math.sin(to_rad(h[i]))
    end
    0.upto(2) do |i|
      m[i] += delta_m[i]
    end
  end
  # format sunrise, sunset & solar noon as DateTime
  sunrise = date.to_datetime + m[1] + days[1]
  transit = date.to_datetime + m[0] + days[0]
  sunset = date.to_datetime + m[2] + days[2]
  [sunrise, transit, sunset]
end

.position(t, longitude, latitude) ⇒ Object

Sun horizontal coordinates (relative position) in degrees:

  • elevation (altitude over horizon) in degrees; positive upwards

  • azimuth in degrees measured clockwise (towards East) from North direction



8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
 
# File 'lib/solar/position.rb', line 8

def position(t, longitude, latitude)

  delta_rad, alpha_rad = equatorial_position_rad(t)
  alpha_deg = to_deg(alpha_rad)
  # alpha_h += 360 if alpha_h < 0

  # t as Julian centuries of 36525 ephemeris days form the epoch J2000.0
  if false
    # Float
    jd = jd_f(t)
  else
    # Rational
    jd = jd_r(t)
  end
  t = to_jc(jd)

  # Sidereal time at Greenwich
  theta = 280.46061837 +
          360.98564736629*(jd - 2_451_545) +
          (0.000387933 - t/38_710_000)*t*t

  # Reduce magnitude to minimize errors
  theta %= 360

  # Local hour angle
  h = theta + longitude - alpha_deg
  h %= 360

  latitude_rad = to_rad(latitude)
  h_rad = to_rad(h)

  # Local horizontal coordinates : Meeus pg 89
  altitude_rad = Math.asin(Math.sin(latitude_rad)*Math.sin(delta_rad) +
                 Math.cos(latitude_rad)*Math.cos(delta_rad)*Math.cos(h_rad))
  azimuth_rad = Math.atan2(
    Math.sin(h_rad),
    Math.cos(h_rad) * Math.sin(latitude_rad) -
    Math.tan(delta_rad) * Math.cos(latitude_rad)
  )

  [to_deg(altitude_rad), (180 + to_deg(azimuth_rad)) % 360]
end

.radiation(t, longitude, latitude, options = {}) ⇒ Object

Radiation in W/m^2 at time t at a position given by longitude, latitud on terrain of given slope and aspect, given the global radiation on a horizontal plane and optionally the diffuse radiation on a horizontal plane. Computes mean radiation in a short period of length givien by :t_step (10 minute by default) centered at t. Basically the method used is the oned presented in:

  • Solar Engineering of Thermal Processes 1991 Duffie, Beckman

Options:

  • :slope slope angle in degrees (0-horizontal to 90-vertical)

  • :aspect clockwise from North in degrees

  • :t_step time step in hours; if radiation values are given they are considered as the mean values for the time step.

  • :global_radiation Global radation measured on a horizontal plane in W/m^2 (mean value over the time step)

  • :diffuse_radiation Diffuse radation measured on a horizontal plane in W/m^2 (mean value over the time step)

  • :albedo of the terrain, used to compute reflected radiation.

This can be used with a measured :global_radiation and optional :diffuse_radiation, both measured on horizontal to compute the estimated global radiation on a sloped plane.

It can be also used by giving the :clearness_index to compute estimated radiation.



32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
 
# File 'lib/solar/radiation.rb', line 32

def radiation(t, longitude, latitude, options = {})
  # TODO: parameterize on the .... algorithms
  #       consider method of [Neteler-2008] Appendix A pg. 377-381

  slope = options[:slope] ||  0.0
  aspect = options[:aspect] || 0.0

  # time step in hours
  t_step = options[:t_step] || 1/6.0
  # global measured radiation in W/m^2 (mean over time step) (horizontal)
  g = options[:global_radiation]
  # optional: diffuse radiation (horizontal)
  g_d = options[:diffuse_radiation]
  k_t = options[:clearness_index]
  k_t ||= 1.0 if g.nil?

  # ground reflectance (albedo) as %
  rg = options[:albedo] || 0.2

  # t is assumed at half the time step

  t_utc = t.utc
  n = t_utc.yday # day of the year (a number between 1 and 365)
  lambda = to_rad(longitude)
  phi = to_rad(latitude)

  d = solar_declination(n)

  # utc time in hours
  tu = t_utc.hour + t_utc.min/60.0 + t_utc.sec/3600.0
  b = to_rad(360.0*(n-1)/365.0)
  # equation of time
  e = 3.82*(0.000075+0.001868*Math.cos(b) - 0.032077*Math.sin(b) - 0.014615*Math.cos(2*b) - 0.04089*Math.sin(2*b))
  # solar time in hours
  ts = tu + longitude/15.0 + e

  # hour angle (omega) in radians
  w = to_rad((ts - 12.0)*15.0)
  cos_w = Math.cos(w)
  sin_w = Math.sin(w)

  # extraterrestrial_normal_radiation in W/m^2
  g_on = ext_radiation(t)

  cos_phi = Math.cos(phi)
  sin_phi = Math.sin(phi)
  cos_d = Math.cos(d)
  sin_d = Math.sin(d)

  # zenith angle in radians
  # eq. 1.6.5 (pg 15) [1991-Duffie]
  cos_phi_z = cos_phi*cos_d*cos_w + sin_phi*sin_d

  # extraterrestrial horizontal radiation in W/m^2
  g_o = g_on*cos_phi_z

  # hour_angle at beginning of time step
  w1 = to_rad((ts - t_step/2 - 12.0)*15.0)
  # hour_angle at end of time step
  w2 = to_rad((ts + t_step/2 - 12.0)*15.0)

  # extraterrestrial horizontal radiation in W/m^2  averaged over the time step
  # [1991-Duffie pg 37] eq. 1.10.1, 1.10.3
  # TODO: for long time steps we should average as:
  #  g_o_a = 12/Math::PI * g_on * ( cos_phi*cos_d*(Math.sin(w2)-Math.sin(w1)) + (w2-w1)*sin_phi*sin_d )
  g_o_a = g_o
  g_o_a = 0 if g_o_a < 0

  # clearness index
  k_t ||= g/g_o_a

  # either k_t or g must be defined by the user
  g ||= (g_o_a == 0 ? 0.0 : k_t*g_o_a)

  # diffuse fraction
  if k_t.infinite?
    df = 1.0 # actual df may be around 0.5; for the purpose of computing g_d it is 1
  else
    solar_elevation = 90 - to_deg(Math.acos(cos_phi_z))
    df = diffuse_fraction(k_t, solar_elevation)
  end

  # diffuse radiation W/m^2
  g_d ||= df*g

  # beam radiation
  g_b = g - g_d

  # slope
  beta = to_rad(slope)
  # azimuth
  gamma = to_rad(aspect-180)

  cos_beta = Math.cos(beta)
  sin_beta = Math.sin(beta)
  cos_gamma = Math.cos(gamma)
  sin_gamma = Math.sin(gamma)

  # angle of incidence
  # eq (1.6.2) pg 14 [1991-Duffie]
  # eq (3) "Analytical integrated functions for daily solar radiation on slopes" - Allen, Trezza, Tasumi
  cos_phi = sin_d*sin_phi*cos_beta \
            - sin_d*cos_phi*sin_beta*cos_gamma \
            + cos_d*cos_phi*cos_beta*cos_w \
            + cos_d*sin_phi*sin_beta*cos_gamma*cos_w \
            + cos_d*sin_beta*sin_gamma*sin_w

  # ratio of beam radiation on tilted surface to beam radiation on horizontal
  # [1991-Duffie pg 23-24] eq. 1.8.1
  # rb = illumination_factor_at(t, longitude, latitude, slope, aspect)
  rb = cos_phi / cos_phi_z
  rb = 0.0 if rb < 0.0

  # anisotropy index
  if k_t.infinite? || g_o_a == 0
    ai = 0.0
  else
    ai = g_b / g_o_a
  end

  # horizontal brightening factor
  if g != 0
    f = Math.sqrt(g_b / g)
  else
    f = 1.0
  end

  if DEBUG_RADIATION
    sun_elevation, sun_azimuth = Solar.position(t_utc, longitude, latitude)
    rb2 = illumination_factor_at(t_utc, longitude, latitude, slope, aspect)
    puts ""
    puts "  @#{t_utc} #{sun_elevation}  [#{90-sun_elevation}] az: #{sun_azimuth} <#{Math.acos(cos_phi_z)*180/Math::PI}>"
    puts "  kt:#{k_t} df:#{df}   g: #{g} gon: #{g_on}"
    puts "  rb: #{rb} (#{rb2}) g_b:#{g_b} g_o_a:#{g_o_a} g_d:#{g_d}"
    puts "  -> #{(g_b + g_d*ai)*rb}+#{g_d*(1-ai)*((1 + cos_beta)/2)*(1 + f*Math.sin(beta/2)**3)}+#{g*rg*(1 - cos_beta)/2}"
  end

  # global radiation on slope according to HDKR model
  # eq. 2.16.7 (pg.92) [1991-Duffie]
  # three terms:
  # * direct and circumsolar: (g_b + g_d*ai)*rb
  # * diffuse: g_d*(1-ai)*((1 + cos_beta)/2)*(1 + f*Math.sin(beta/2)**3)
  # * reflected: g*rg*(1 - cos_beta)/2
  (g_b + g_d*ai)*rb + g_d*(1-ai)*((1 + cos_beta)/2)*(1 + f*Math.sin(beta/2)**3) + g*rg*(1 - cos_beta)/2
end

.rise(date, longitude, latitude, options = {}) ⇒ Object

Sun rise time for a given date (UTC) and position.

The :zenith or :altitude of the sun can be passed as an argument, which can be numeric (in degrees) or symbolic: :official, :civil, :nautical or :astronomical.

nil is returned if the sun doesn’t rise at the date and position.



11
12
13
14
15
16
17
18
 
# File 'lib/solar/passages.rb', line 11

def rise(date, longitude, latitude, options = {})
  rising, _, setting = passages(date, longitude, latitude, options)
  if rising == setting || (setting - rising) == 1
    nil # rising==setting => no rise; setting-rising == 1 => no set
  else
    rising
  end
end

.rise_and_set(date, longitude, latitude, options = {}) ⇒ Object

Rise and set times as given by rise() and set()



38
39
40
41
42
43
44
45
 
# File 'lib/solar/passages.rb', line 38

def rise_and_set(date, longitude, latitude, options = {})
  rising, _, setting = passages(date, longitude, latitude, options)
  if rising==setting || (setting-rising)==1
    nil # rising==setting => no rise; setting-rising == 1 => no set
  else
    [rising, setting]
  end
end

.set(date, longitude, latitude, options = {}) ⇒ Object

Sun set time for a given date (UTC) and position.

The :zenith or :altitude of the sun can be passed as an argument, which can be numeric (in degrees) or symbolic: :official, :civil, :nautical or :astronomical.

nil is returned if the sun doesn’t set at the date and position.



28
29
30
31
32
33
34
35
 
# File 'lib/solar/passages.rb', line 28

def set(date, longitude, latitude, options = {})
  rising, _, setting = passages(date, longitude, latitude, options)
  if rising==setting || (setting-rising)==1
    nil # rising==setting => no rise; setting-rising == 1 => no set
  else
    setting
  end
end