r.sun: allow user to set solar constant (#482)

OSGeo · Apr 4, 2020 · 9b0d34b · 9b0d34b
1 parent 9de8682
commit 9b0d34b
Show file tree

Hide file tree

Showing 4 changed files with 21 additions and 5 deletions.
diff --git a/raster/r.sun/main.c b/raster/r.sun/main.c
@@ -187,6 +187,7 @@ double o_orig, z1;
 double horizonStep;
 double ltime, tim, timo;
 double declination;	/* Contains the negative of the declination at the chosen day */
+double solar_constant;
 
 /*
  * double lum_C31_l, lum_C33_l;
@@ -230,7 +231,7 @@ int main(int argc, char *argv[])
 	struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
 	    *lin, *albedo, *longin, *alb, *latin, *coefbh, *coefdh,
 	    *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
-	    *glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
+	    *glob_rad, *day, *step, *declin, *solar_cnst, *ltime, *dist, *horizon,
 	    *horizonstep, *numPartitions, *civilTime, *threads;
     }
     parm;
@@ -470,6 +471,14 @@ int main(int argc, char *argv[])
     parm.declin->description =
 	_("Declination value (overriding the internally computed value) [radians]");
 
+    parm.solar_cnst = G_define_option();
+    parm.solar_cnst->key = "solar_constant";
+    parm.solar_cnst->type = TYPE_DOUBLE;
+    parm.solar_cnst->required = NO;
+    parm.solar_cnst->answer = "1367";
+    parm.solar_cnst->description =
+	_("Solar constant [W/m^2]");
+
     parm.ltime = G_define_option();
     parm.ltime->key = "time";
     parm.ltime->type = TYPE_DOUBLE;
@@ -741,6 +750,11 @@ int main(int argc, char *argv[])
 	declination = -declin;
     }
 
+    if (parm.solar_cnst->answer)
+        sscanf(parm.solar_cnst->answer, "%lf", &solar_constant);
+    else
+        solar_constant = 1367;
+
     if (ttime != 0) {
 	/* Shadow for just one time during the day */
 	if (horizon == NULL) {
@@ -2026,7 +2040,7 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 
     Rast_append_format_history(
 	&hist,
-	" Solar constant (W/m^2):                   1367");
+	" Solar constant (W/m^2):                   %f", solar_constant);
     Rast_append_format_history(
 	&hist,
 	" Extraterrestrial irradiance (W/m^2):      %f",

diff --git a/raster/r.sun/r.sun.html b/raster/r.sun/r.sun.html
@@ -137,7 +137,7 @@ <h2>DESCRIPTION</h2>
 <p>
 The maps' history files are generated containing the following listed 
 parameters used in the computation: <br>
-- Solar constant 1367 W.m-2 <br>
+- Solar constant used W.m-2 <br>
 - Extraterrestrial irradiance on a plane perpendicular to the solar beam [W.m-2] <br>
 - Day of the year <br>
 - Declination [radians] <br>

diff --git a/raster/r.sun/rsunglobals.h b/raster/r.sun/rsunglobals.h
@@ -54,6 +54,8 @@ extern const double pi2;
 extern const double deg2rad;
 extern const double rad2deg;
 
+extern double solar_constant;
+
 extern struct pj_info iproj;
 extern struct pj_info oproj;
 

diff --git a/raster/r.sun/rsunlib.c b/raster/r.sun/rsunlib.c
@@ -113,7 +113,7 @@ double com_sol_const(int no_of_day)
 {
     double I0, d1;
 
-    /* Solar constant: 1367.0 W/m^2.
+    /* Solar constant: 1367.0 W/m^2. Note: solar constant is parameter.
 
        Perigee offset: here we call Jan 2 at 8:18pm the Perigee, so day
 	number 2.8408. In angular units that's (2*pi * 2.8408 / 365.25) = 0.048869.
@@ -128,7 +128,7 @@ double com_sol_const(int no_of_day)
 
     /*  v W/(m*m) */
     d1 = pi2 * no_of_day / 365.25;
-    I0 = 1367. * (1 + 0.03344 * cos(d1 - 0.048869));
+    I0 = solar_constant * (1 + 0.03344 * cos(d1 - 0.048869));
 
     return I0;
 }