[GRASS-SVN] r30653 - grass-addons/raster/r.sun_horizon/r.sun2

[svn_grass at osgeo.org]

Author: neteler
Date: 2008-03-20 09:20:53 -0400 (Thu, 20 Mar 2008)
New Revision: 30653

Modified:
   grass-addons/raster/r.sun_horizon/r.sun2/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.


Modified: grass-addons/raster/r.sun_horizon/r.sun2/main.c
===================================================================
--- grass-addons/raster/r.sun_horizon/r.sun2/main.c	2008-03-19 22:39:54 UTC (rev 30652)
+++ grass-addons/raster/r.sun_horizon/r.sun2/main.c	2008-03-20 13:20:53 UTC (rev 30653)
@@ -174,6 +174,30 @@
 double cbh, cdh;
 double TOLER;
 
+
+#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[])
 {
 	double singleSlope;
@@ -710,6 +734,12 @@
 
 /**********end of parser - ******************************/
 
+	if ((G_projection() == PROJECTION_LL))
+		{
+		ll_correction=1;	
+		}
+
+
     calculate(singleSlope, singleAspect, singleAlbedo,
 	      singleLinke,gridGeom);
     OUTGR();
@@ -1557,7 +1587,7 @@
         dx = (double)i *gridGeom->stepx;
         dy = (double)j *gridGeom->stepy;
 
-	*length = DISTANCE1(gridGeom->xg0, dx, gridGeom->yg0, dy); /* dist from orig. grid point to the current grid point */
+	*length = distance(gridGeom->xg0, dx, gridGeom->yg0, dy); /* dist from orig. grid point to the current grid point */
 
 	sunVarGeom->zp = z[j][i];
 
@@ -1732,6 +1762,7 @@
 	double longitTime=0.;
 	double locTimeOffset;
 	double latitude, longitude;
+	double coslat;
 	
 	
 	struct SunGeometryConstDay sunGeom;
@@ -1885,6 +1916,11 @@
 			gridGeom.xp=xmin+gridGeom.xx0;
 			gridGeom.yp=ymin+gridGeom.yy0;
 
+			if(ll_correction)
+			{
+				coslat=cos(deg2rad*gridGeom.yp);
+				coslatsq = coslat*coslat;
+			}
 
         		func = NULL;