Unit test and change on daily global insolation on inclined surface.

functions/G_al_i.R → functions/G_ali.R

@@ -9,7 +9,7 @@ Gh_eq = dget(here("functions", "G_h.R"))
 source(here("functions", "albedo.R"))
 
 function(Ls, phi, longitude, T_s, Z=Z_eq(Ls=Ls, T_s=T_s, phi=phi, nfft=nfft), tau, al=albedo(latitude=phi, longitude=longitude, tau=tau), beta, nfft){
-  Gal_i = al * Gh_eq(Ls=Ls, phi=phi, longitude=longitude, Z=Z, tau=tau, al=al, nfft=nfft) * sin((beta*pi/180) / 2)^2
-  return(Gal_i)
+  Gali = al * Gh_eq(Ls=Ls, phi=phi, longitude=longitude, Z=Z, tau=tau, al=al, nfft=nfft) * sin((beta*pi/180) / 2)^2
+  return(Gali)
 }
 

functions/G_bi.R

@@ -40,70 +40,74 @@ function(Ls, phi, T_s, Z=Z_eq(Ls=Ls, T_s=T_s, phi=phi, nfft=nfft), tau, beta, ga
   #   between the Mars solar time T and the hour angle as for the Earth.
   omega_deg = 15 * T_s - 180
 
-  # Equation 6 (1993): Solar Azimuth Angle [deg]
-  solar_azimuth_angle = function(){
-    x = sin(phi*pi/180) * cos(delta) * cos(omega_deg*pi/180)
-    y = cos(phi*pi/180) * sin(delta)
-    z = sin(Z*pi/180)
-
-    gamma_s = acos((x - y) / z) * 180/pi # [deg]
-
-    # Alternatively, Equation 7 (1993): Solar Azimuth Angle [deg]
-    # x = sin(phi*pi/180) * cos(Z*pi/180) - sin(delta)
-    # y = cos(phi*pi/180) * sin(Z*pi/180)
-    # 
-    # gamma_s = acos(x / y) * 180/pi # [deg]
-
-    return(gamma_s)
-  } 
-
-  # Set to omega to exactly zero when omega near zero instead because sometimes it is essentially zero with values like -2.8421709430404e-14 deg.
-  for(i in 1:length(omega_deg)){
-    if(omega_deg[i] > -0.1 && omega_deg[i] < 0.1){
-      omega_deg[i] = 0
-    }
+  # Set to omega to exactly zero when omega near zero because sometimes it is essentially zero with values like -2.8421709430404e-14 deg.
+  zero = function(x){
+    return(
+      ifelse(x > -1e-5 && x < 1e-5, 0, x)
+    )
   }
 
+  # Apply the function that will set omega to zero if it is close to zero.
+  omega_deg = sapply(omega_deg, zero)
+
+  # Equation 6 (1993): Solar Azimuth Angle [deg]
+  x = sin(phi*pi/180) * cos(delta) * cos(omega_deg*pi/180)
+  y = cos(phi*pi/180) * sin(delta)
+  z = sin(Z*pi/180)
+
   # From (32) in (1993): It is solar noon when omega is 0 deg. This translates to gamma_s = 0 deg.
-  gamma_s = ifelse(omega_deg == 0, 0, solar_azimuth_angle())
+  # The following operation will result in a value very close to 1 but not exactly 1 if it is
+  # solar noon, i.e. when omega is 0 deg. When op = 1 -> acos(op) = 0, i.e. (32) in (1993).
+  op = ((x - y) / z)
 
+  # The following function is to make sure that op is exactly 1 when omega is 0 deg
+  one = function(x){
+    return(
+      ifelse(x > 0,
+             ifelse(x > 1 && x < 1+1e-5, 1, x),  # If x is positive.
+             ifelse(x < -1 && x > -1-1e-5, 1, x)) # If x is negative.
+    )
+  }
+
+  # Apply the function.
+  op = sapply(op, one)
 
+  # Calculate gamma_s.
+  gamma_s = acos(op) * 180/pi # [deg]
+
+  # Alternatively, Equation 7 (1993): Solar Azimuth Angle [deg]
+  # x = sin(phi*pi/180) * cos(Z*pi/180) - sin(delta)
+  # y = cos(phi*pi/180) * sin(Z*pi/180)
+  # 
+  # gamma_s = acos(x / y) * 180/pi # [deg]
+
   # Sun Angle of Incidence [rad].
-  #   Equation 13 (1993): Inclined surface.
-  # FIXME: Equation 23 (1993): Vertical surface.
-  #        Equation 27 (1993): Facing the equator.
-  #        Equation 29 (1993): Beta = 0.
   sun_angle_of_incidence = function(){
-    # if(beta == 0){ # Equation 29 (1993): Beta = 0.
-    #   return(NULL)
-    #   
-    # }else if(beta == 90){ #Equation 29 (1993): Beta = 0.
-    #   return(NULL)
-    #   
-    # }else if(gamma_c == 0){ # Equation 27 (1993): Facing the equator.
-    #   return(NULL)
-    #   
-    # }else{
-    #   #  Inclined surface.
-    #   i = cos(beta * pi/180) * cos(Z * pi/180)
-    #   j = sin(beta * pi/180) * sin(Z * pi/180) * cos((gamma_s - gamma_c) * pi/180) # Does not matter when beta = 0 because it leads to j = 0.
-    #   teta = acos(i + j) # [rad]
-    # }
-
-    # Inclined surface.
-    i = cos(beta * pi/180) * cos(Z * pi/180)
-    j = sin(beta * pi/180) * sin(Z * pi/180) * cos((gamma_s - gamma_c) * pi/180) # Does not matter when beta = 0 because it leads to j = 0.
-    teta = acos(i + j) # [rad]
+    teta = NULL
 
+    # (23) in (1993) Vertical surface.
+    if(beta == 90){
+      i = cos(gamma_s * pi/180) * cos(gamma_c * pi/180)
+      #TODO:  Double check this, why wouldn't the paper use squared instead of multiply the by the same value? 
+      #       Check with (24) in (1993). 
+      j = sin(gamma_s * pi/180) * sin(gamma_s * pi/180) 
+      teta = acos(sin(Z * pi/180) * (i + j))# [rad]
+
+    }else{
+      # (13) in (1993): Inclined surface.
+      # (27) in (1993): Surface facing the equator.
+      i = cos(beta * pi/180) * cos(Z * pi/180)
+      j = sin(beta * pi/180) * sin(Z * pi/180) * cos((gamma_s - gamma_c) * pi/180) # Does not matter when beta = 0 because it leads to j = 0.
+      teta = acos(i + j) # [rad]
+    }
+
     return(teta)
   }
 
-
   # Sun Angle of Incidence [rad] on an inclined surface.
   teta = sun_angle_of_incidence()
 
   Gbi = Gb_eq(Ls=Ls, Z=Z, tau=tau) * cos(teta)
 
   return(Gbi)
-}
-
+}

functions/G_i.R

@@ -15,7 +15,7 @@ Z_eq = dget(here("functions", "Z.R"))
 
 Gbi_eq = dget(here("functions", "G_bi.R"))
 Gdi_eq = dget(here("functions", "G_di.R"))
-Gali_eq = dget(here("functions", "G_al_i.R"))#
+Gali_eq = dget(here("functions", "G_ali.R"))#
 
 source(here("functions", "albedo.R"))
 

functions/H_al_i.R → functions/H_ali.R

@@ -1,6 +1,6 @@
 library(here)
 
-Ial_i_eq = dget(here("functions", "I_al_i.R"))
+Iali_eq = dget(here("functions", "I_ali.R"))
 
 source(here("functions", "albedo.R"))
 
@@ -9,9 +9,9 @@ function(Ls, phi, longitude, tau, al=albedo(latitude=phi, longitude=longitude, t
     stop("Surface azimuth angle gamma_c must between -180 and 180 degress with zero south, east negative, and west positive.")
   }
 
-  # H_i is obtained by integrating Ial_i over the period from sunrise to sunset.
-  H_i = Ial_i_eq(Ls=Ls, phi=phi, longitude=longitude, tau=tau, T_start=0, T_end=24, al=al, beta=beta, gamma_c=gamma_c, nfft=nfft)
+  # H_alii is obtained by integrating I_ali over the period from sunrise to sunset.
+  H_ali = Iali_eq(Ls=Ls, phi=phi, longitude=longitude, tau=tau, T_start=0, T_end=24, al=al, beta=beta, gamma_c=gamma_c, nfft=nfft)
 
   # Return result.
-  return(H_i)
+  return(H_ali)
 }

functions/H_bi.R

@@ -1,15 +1,15 @@
 library(here)
 
-Ib_beta_eq = dget(here("functions", "I_b_beta.R"))
+Ibi_eq = dget(here("functions", "I_bi.R"))
 
-function(Ls, phi, tau, al=0.1, beta, gamma_c, nfft){
+function(Ls, phi, tau, beta, gamma_c, nfft){
   if(gamma_c > 180 || gamma_c < -180){
     stop("Surface azimuth angle gamma_c must between -180 and 180 degress with zero south, east negative, and west positive.")
   }
 
-  # H_b_beta is obtained by integrating I_b_beta over the period from sunrise to sunset.
-  H_b_beta = Ib_beta_eq(Ls=Ls, phi=phi, tau=tau, T_start=0, T_end=24, al=al, beta=beta, gamma_c=gamma_c, nfft=nfft)
+  # H_bi is obtained by integrating I_bi over the period from sunrise to sunset.
+  H_bi = Ibi_eq(Ls=Ls, phi=phi, tau=tau, T_start=0, T_end=24, beta=beta, gamma_c=gamma_c, nfft=nfft)
 
   # Return result.
-  return(H_b_beta)
+  return(H_bi)
 }

functions/I_al_i.R → functions/I_ali.R

@@ -1,6 +1,6 @@
 library(here)
 
-Gal_i_eq = dget(here("functions", "G_al_i.R"))
+Gali_eq = dget(here("functions", "G_ali.R"))
 
 source(here("functions", "albedo.R"))
 
@@ -26,12 +26,12 @@ function(Ls, phi, longitude, tau, T_start, T_end, al=albedo(latitude=phi, longit
 
   # Step 2: Calculate insolation.
   interand = function(T_s){
-    G_al_i = Gal_i_eq(Ls=Ls, phi=phi, longitude=longitude, T_s=T_s, tau=tau, al=al, beta=beta, nfft=nfft)
-    return(G_al_i)
+    Gali = Gali_eq(Ls=Ls, phi=phi, longitude=longitude, T_s=T_s, tau=tau, al=al, beta=beta, nfft=nfft)
+    return(Gali)
   }
 
-  I_al_i = integrate(interand, T_start, T_end)
+  Iali = integrate(interand, T_start, T_end)
 
   # Return integration result.
-  return(I_al_i$value)
+  return(Iali$value)
 }

functions/I_bi.R

@@ -1,11 +1,11 @@
 library(here)
 
-Gb_beta_eq = dget(here("functions", "G_b_beta.R"))
+Gbi_eq = dget(here("functions", "G_bi.R"))
 
 # Constrain T_start and T_end based on sunrise and sunset times.
 constrain_solar_time_range = dget(here("utils", "constrain_solar_time_range.R")) 
 
-function(Ls, phi, tau, T_start, T_end, al=0.1, beta, gamma_c, nfft){
+function(Ls, phi, tau, T_start, T_end, beta, gamma_c, nfft){
 
   if(gamma_c > 180 || gamma_c < -180){
     stop("Surface azimuth angle gamma_c must between -180 and 180 degres with zero south, east negative, and west positive.")
@@ -29,12 +29,12 @@ function(Ls, phi, tau, T_start, T_end, al=0.1, beta, gamma_c, nfft){
   # Step 2: Calculate insolation.
 
   interand = function(T_s){
-    G_b_beta = Gb_beta_eq(Ls=Ls, phi=phi, T_s=T_s, tau=tau, al=al, beta=beta, gamma_c=gamma_c, nfft=nfft)
-    return(G_b_beta)
+    G_bi = Gbi_eq(Ls=Ls, phi=phi, T_s=T_s, tau=tau, beta=beta, gamma_c=gamma_c, nfft=nfft)
+    return(G_bi)
   }
 
-  I_b_beta = integrate(interand, T_start, T_end)
+  I_bi = integrate(interand, T_start, T_end)
 
   # Return integration result.
-  return(I_b_beta$value)
+  return(I_bi$value)
 }

functions/f.R

@@ -139,7 +139,8 @@ f_90 = function(Z, tau, al=0.1){
 f_analytical = function(Z, tau, al=0.1){
   # Check for and warn against parameters that would result in lagest errors (max. 7%).
   if(tau > 5){
-    warning(paste("Large error encountered with τ = ", tau, " greater than 5 (maximum error is 7%). ", 
+    # TODO: Make warning.
+    message(paste("Large error encountered with τ = ", tau, " greater than 5 (maximum error is 7%). ", 
                   "Consider using the f_89 and f_90 table lookup implementation of the normalized net flux function instead of its analytical expression.",
                   sep="")
     )
@@ -155,7 +156,8 @@ f_analytical = function(Z, tau, al=0.1){
   # handle warning message.
   for(w_msg in warning_msg){
     if(w_msg != ""){
-      warning(w_msg)
+      # TODO: Make warning.
+      message(w_msg)
     }
   }
 

test/data.R

@@ -0,0 +1,37 @@
+library(here)
+Hbi_eq = dget(here("functions", "H_bi.R"))
+
+test_data_dir = "test/data/"
+
+spacecrafts = list(
+  "VL1" = list(
+    "latitude" = 22.3,
+    "longitude" = -49.97,
+    "beta_optimal" = 6.5
+  ),
+  "VL2" = list(
+    "latitude" = 47.7,
+    "longitude" = 134.29,
+    "beta_optimal" = 22
+  )
+)
+
+# Expected data
+expected_data = list(
+  "VL1" = list(
+    "tau" = read.csv(here(test_data_dir, "discretized_tau_at_vl1_fig_3_1991_update.csv")),
+    "insolation" =
+      list(
+        "beta_equals_phi" = read.csv(here(test_data_dir, "daily_insolation_on_an_inclined_surface_for_beta_equals_phi_at_vl1_table_ii_1993.csv")),
+        "beta_optimal" = read.csv(here(test_data_dir, "daily_insolation_on_optimal_inclined_angle_beta_at_vl1_table_iii_1993.csv"))
+      )
+  ),
+  "VL2" = list(
+    "tau" = read.csv(here(test_data_dir, "discretized_tau_at_vl2_fig_3_1991_update.csv")),
+    "insolation" =
+      list(
+        "beta_equals_phi" = read.csv(here(test_data_dir, "daily_insolation_on_an_inclined_surface_for_beta_equals_phi_at_vl2_table_ii_1993.csv")),
+        "beta_optimal" = read.csv(here(test_data_dir, "daily_insolation_on_optimal_inclined_angle_beta_at_vl2_table_iii_1993.csv"))
+      )
+  )
+)

test/data/daily_insolation_on_an_inclined_surface_for_beta_equals_phi_at_vl1_table_ii_1993.csv

@@ -46,7 +46,7 @@ Ls,H,Hb,Hd,Hal
 220,2030.8,204.1,1798,28.7
 225,2007.7,218.6,1762.1,27
 230,2043.8,273.3,1745.9,24.6
-235,2027.8,207.5,1707.4,22.9
+235,2027.8,297.5,1707.4,22.9
 240,2025.7,326.9,1677.3,21.5
 245,2068.8,393.1,1656,19.7
 250,2125.1,479.8,1627.4,17.9

test/data/daily_insolation_on_an_inclined_surface_for_beta_equals_phi_at_vl2_table_ii_1993.csv

@@ -9,7 +9,7 @@ Ls,H,Hb,Hd,Hal
 35,3260.4,1978.6,1155.6,126.2
 40,3361.8,2127.0,1102.4,132.4
 45,3312.9,1979.0,1200.1,133.8
-50,2919.9,12313,1566.2,122.2
+50,2919.9,1231.5,1566.2,122.2
 55,3353.0,1971.8,1240.6,140.6
 60,3370.6,1964.1,1259.0,143.5
 65,3387.8,1965.2,1276.4,146.2
@@ -30,7 +30,7 @@ Ls,H,Hb,Hd,Hal
 140,3852.2,2559.9,1141.2,151.1
 145,3955.8,2782.9,1022.5,150.4
 150,3861.6,2643.5,1114.3,143.8
-155,1850.2,2613.5,1097.3,139.4
+155,3850.2,2613.5,1097.3,139.4
 160,3950.3,2834.4,978.6,137.3
 165,3657.3,2395.0,1135.7,126.6
 170,3593.8,2365.8,1107.0,121.0