[GRASS-SVN] r30670 - grass/branches/releasebranch_6_3/raster/r.sun

[svn_grass at osgeo.org]

Author: neteler
Date: 2008-03-21 06:18:05 -0400 (Fri, 21 Mar 2008)
New Revision: 30670

Modified:
   grass/branches/releasebranch_6_3/raster/r.sun/description.html
   grass/branches/releasebranch_6_3/raster/r.sun/main.c
Log:
Thomas Huld <Thomas.Huld jrc.it>: bugfix: the code did not calculate the distance along the line of sight to the sun correctly, it calculated them in lat/lon distances, while the height difference was in meters! This was never correct in r.sun.
This is a workaround for lat/lon coordinates and it then REQUIRES that the elevation map is in meters. The rule is: - lat/lon coordinates: elevation in meters - Other coordinates: elevation in the same unit as the easting-northing coordinates. (merge from HEAD)


Modified: grass/branches/releasebranch_6_3/raster/r.sun/description.html
===================================================================
--- grass/branches/releasebranch_6_3/raster/r.sun/description.html	2008-03-21 10:16:52 UTC (rev 30669)
+++ grass/branches/releasebranch_6_3/raster/r.sun/description.html	2008-03-21 10:18:05 UTC (rev 30670)
@@ -10,6 +10,13 @@
 for the earth curvature is internally calucated.<br>
 The units of the parameters are specified in brackets, a hyphen in the brackets
 explains that the parameter has no units.
+<P>
+For latitude-longitude coordinates it requires that the elevation map is in meters.
+The rules are:
+<ul>
+<li> lat/lon coordinates: elevation in meters;
+<li> Other coordinates: elevation in the same unit as the easting-northing coordinates.
+</ul>
 
 <P>
 The solar geometry of the model is based on the works of Krcho (1990), later
@@ -254,6 +261,8 @@
                                                                         
 Marcel Suri, GeoModel, s.r.o. Bratislava, Slovakia <br>
 
+Thomas Huld, JRC, Italy <br>
+
 &copy; 2002, Jaroslav Hofierka, Marcel Suri 
 <address>
 <a href="MAILTO:hofi at geomodel.sk">hofierka at geomodel.sk</a>

Modified: grass/branches/releasebranch_6_3/raster/r.sun/main.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/r.sun/main.c	2008-03-21 10:16:52 UTC (rev 30669)
+++ grass/branches/releasebranch_6_3/raster/r.sun/main.c	2008-03-21 10:18:05 UTC (rev 30670)
@@ -6,7 +6,8 @@
 See manual pages for details.
 (C) 2002 Copyright Jaro Hofierka, Gresaka 22, 085 01 Bardejov, Slovakia, 
               and GeoModel, s.r.o., Bratislava, Slovakia
-email: hofierka at geomodel.sk,marcel.suri at jrc.it,suri at geomodel.sk
+email: hofierka at geomodel.sk,marcel.suri at jrc.it,suri at geomodel.sk,
+       Thomas Huld <Thomas.Huld at jrc.it> (lat/long fix)
 *******************************************************************************/
 /*
  * This program is free software; you can redistribute it and/or
@@ -58,6 +59,9 @@
 #include <grass/gprojects.h>
 #include <grass/glocale.h>
 
+const double deg2rad = M_PI / 180.;
+const double rad2deg = 180. / M_PI;
+
 FILE *felevin, *faspin, *fslopein, *flinkein, *falbedo, *flatin;
 FILE *fincidout, *fbeam_rad, *finsol_time, *fdiff_rad, *frefl_rad;
 
@@ -130,6 +134,23 @@
 double beam_e, diff_e, refl_e, bh, dh, rr, insol_t;
 double cbh, cdh;
 
+#define DEGREEINMETERS 111120.
+
+int ll_correction = 0;
+double coslatsq;
+
+double distance(double x1, double x2, double y1, double y2)
+{
+    if (ll_correction) {
+	return DEGREEINMETERS * sqrt(coslatsq * (x1 - x2) * (x1 - x2)
+				     + (y1 - y2) * (y1 - y2));
+    }
+    else {
+	return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+    }
+}
+
+
 int main(int argc, char *argv[])
 {
 
@@ -389,6 +410,10 @@
     diff_rad = parm.diff_rad->answer;
     refl_rad = parm.refl_rad->answer;
 
+    if ((G_projection() == PROJECTION_LL)) {
+	ll_correction = 1;
+    }
+
     if ((insol_time != NULL) && (incidout != NULL))
 	G_fatal_error(_("insol_time and incidout are incompatible options"));
 
@@ -1159,7 +1184,7 @@
 	dx = (double)i *stepx;
 	dy = (double)j *stepy;
 
-	length = DISTANCE1(xg0, dx, yg0, dy);	/* dist from orig. grid point to the current grid point */
+	length = distance(xg0, dx, yg0, dy);	/* dist from orig. grid point to the current grid point */
 
 	cube(j, i);
 	return;
@@ -1288,6 +1313,11 @@
 	    func = NULL;
 	    length = 0;
 
+	    if (ll_correction) {
+		coslat = cos(deg2rad * yp);
+		coslatsq = coslat * coslat;
+	    }
+
 	    z_orig = z1 = zp = z[j][i];
 	    o_orig = o[j][i];