[GRASS-SVN] r34554 - in grass/branches/develbranch_6/raster: . r.horizon r.sun r.sun2

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Nov 27 16:45:30 EST 2008


Author: neteler
Date: 2008-11-27 16:45:30 -0500 (Thu, 27 Nov 2008)
New Revision: 34554

Added:
   grass/branches/develbranch_6/raster/r.horizon/
   grass/branches/develbranch_6/raster/r.horizon/Makefile
   grass/branches/develbranch_6/raster/r.horizon/TODO
   grass/branches/develbranch_6/raster/r.horizon/description.html
   grass/branches/develbranch_6/raster/r.horizon/local_proto.h
   grass/branches/develbranch_6/raster/r.horizon/main.c
   grass/branches/develbranch_6/raster/r.sun/DEPRECATED
   grass/branches/develbranch_6/raster/r.sun2/
   grass/branches/develbranch_6/raster/r.sun2/Makefile
   grass/branches/develbranch_6/raster/r.sun2/README
   grass/branches/develbranch_6/raster/r.sun2/TODO
   grass/branches/develbranch_6/raster/r.sun2/description.html
   grass/branches/develbranch_6/raster/r.sun2/local_proto.h
   grass/branches/develbranch_6/raster/r.sun2/main.c
   grass/branches/develbranch_6/raster/r.sun2/rsunglobals.h
   grass/branches/develbranch_6/raster/r.sun2/rsunlib.c
   grass/branches/develbranch_6/raster/r.sun2/sunradstruct.h
Log:
moved here from GRASS Addons; deactivated r.sun

Added: grass/branches/develbranch_6/raster/r.horizon/Makefile
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/Makefile	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.horizon/Makefile	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.horizon
+
+LIBES =  $(GPROJLIB) $(GISLIB)
+DEPENDENCIES = $(GPROJDEP) $(GISDEP)
+EXTRA_INC = $(PROJINC)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd


Property changes on: grass/branches/develbranch_6/raster/r.horizon/Makefile
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.horizon/TODO
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/TODO	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.horizon/TODO	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,9 @@
+Probably the sun position calculation shouldbe replaced
+with 
+
+ G_calc_solar_position()
+
+currently used in r.sunmask. G_calc_solar_position() is based on
+solpos.c from NREL and very precise.
+
+MN 4/2001

Added: grass/branches/develbranch_6/raster/r.horizon/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/description.html	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.horizon/description.html	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,177 @@
+<H2>DESCRIPTION</H2>
+
+<B>r.horizon</B> computes the angular height of terrain horizon in
+radians. It reads a raster of elevation data and outputs the horizon
+outline in one of two modes:
+
+<ul>
+<li> single point: as a series of horizon 
+heights in the specified directions from the given point. The results are
+written to the stdout.
+<li> raster: in this case the output is
+one or more rasters, with each point in a raster giving the horizon
+height in a specific direction. One raster is created for each direction.
+</ul>
+
+<P>
+The directions are given as azimuthal angles (in degrees), with
+the angle starting with 0 towards East and moving counterclockwise
+(North is 90, etc.). The calculation takes into account the actual
+projection, so the angles are corrected for direction distortions
+imposed by it. The directions are thus aligned to those of the
+geographic projection and not the coordinate system given by the rows
+and columns of the raster map. This correction implies that the
+resulting cardinal directions represent true orientation towards the
+East, North, West and South. The only exception of this feature is
+LOCATION with x,y coordinate system, where this correction is
+not applied. 
+
+<H3>Flags:</H3>
+<dl>
+  <dt><B>-d</B>
+  <dd>Output horizon height in degrees (the default is radians)</dl>
+</dl>
+
+<H3>Input parameters:</H3>
+<P>The <I>elevin </I>parameter is an input elevation raster map. If
+the buffer options are used (see below), this raster should extend
+over the area that accommodate the presently defined region plus
+defined buffer zones. 
+</P>
+<P>The <I>step </I>parameter gives the angle step (in degrees)
+between successive azimuthal directions for the calculation of the
+horizon. Thus, a value of 5 for the <I>step </I>will give a total of
+360/5=72 directions (72 rasters if used in the raster mode). 
+</P>
+<P>The <I>direction</I> parameter gives the initial direction of the
+first output. This parameter acts as an direction angle offset. For
+example, if you want to get horizon angles for directions 45 and 225
+degrees, the <I>direction</I> should be set to 45 and <I>step</I> to
+180. If you only want one single direction, use this parameter to
+specify desired direction of horizon angle, and set the <I>step</I>
+size to 0 degrees. 
+</P>
+<P>The <I>dist </I>controls the sampling distance step size for the
+search for horizon along the line of sight. The default value is 1.0
+meaning that the step size will be taken from the raster resolution.
+Setting the value below 1.0 might slightly improve results for
+directions apart from the cardinal ones, but increasing the
+processing load of the search algorithm. 
+</P>
+<P>The <I>maxdistance</I> value gives a maximum distance to move away
+from the origin along the line of sight in order to search for the
+horizon height. The smaller this value the faster the calculation but
+the higher the risk that you may miss a terrain feature that can
+contribute significantly to the horizon outline.</P>
+<P>The <I>coord</I> parameter takes a pair of easting-northing values
+in the current coordinate system and calculates the values of angular
+height of the horizon around this point. To achieve the
+consistency of the results, the point coordinate is
+aligned to the midpoint of the closest elevation raster cell. 
+</P>
+<P>If an analyzed point (or raster cell) lies close to the edge of
+the defined region, the horizon calculation may not be realistic,
+since it may not see some significant terrain features which could
+have contributed to the horizon, because these features are outside
+the region. There are to options how to set the size of the buffer
+that is used to increase the area of the horizon analysis. The
+<I>bufferzone</I> parameter allows you to specify the same size of
+buffer for all cardinal directions and the parameters <I>e_buff</I>,
+<I>n_buff</I>, <I>s_buff</I>, and <I>w_buff</I> allow you to specify
+a buffer size individually for each of the four directions. The
+buffer parameters influence only size of the read elevation map,
+while the analysis in the raster mode will be done only for the
+area specified by the current region definition.</P>
+<P>The <I>horizon </I>parameter gives the prefix of the output
+horizon raster maps. The raster name of each horizon direction
+raster will be constructed as <I>horizon_</I>NNN , where NNN counts
+upwards from 0 to total number of directions. If you use <B>r.horizon</B>
+in the single point mode this option will be ignored. 
+</P>
+
+<p>
+At the moment the elevation and maximum distance must be measured in meters, 
+even if you use geographical coordinates (longitude/latitude). If your 
+projection is based on distance (easting and northing), these too must 
+be in meters. The buffer parameters must be in the same units as the 
+raster coordinates.
+</p>
+
+<H2>METHOD</H2>
+<P>The calculation method is based on the method used in <B>r.sun</B>
+to calculate shadows. It starts at a very shallow angle and walks
+along the line of sight and asks at each step whether the line of
+sight &quot;hits&quot; the terrain. If so, the angle is increased to
+allow the line of sight to pass just above the terrain at that point.
+This is continued until the line of sight reaches a height that is
+higher than any point in the region or until it reaches the border of
+the region (see also the <I>bufferzone,e_buff</I>, <I>n_buff</I>,
+<I>s_buff</I>, and <I>w_buff</I>). The the number of lines of sight (azimuth 
+directions) is determined from the <I>direction</I> and
+<I>step </I>parameters. The method takes into account the curvature
+of the Earth whereby remote features will seem to be lower than they
+actually are. It also accounts for the changes of angles towards
+cardinal directions caused by the projection (see above). 
+</P>
+
+<H2>EXAMPLE</H2>
+
+Single point mode:
+<div class="code"><pre>
+r.horizon elevin=DEM step=30 direction=15 bufferzone=200 coord=47.302,7.365 dist=0.7 &gt; horizon.out
+</pre></div>
+
+
+Raster map mode (for r.sun2):
+<div class="code"><pre>
+# we put a bufferzone of 10% of maxdistance around the study area
+r.horizon elevin=DEM step=30 bufferzone=200 horizon=horangle dist=0.7 maxdistance=2000
+</pre></div>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="r.sun2.html">r.sun2</A>,
+<A HREF="r.los.html">r.los</A></em>
+
+<H2>REFERENCES</H2>
+<P>Hofierka J., 1997. Direct solar radiation modelling within an
+open GIS environment. Proceedings of JEC-GI'97 conference in Vienna,
+Austria, IOS Press Amsterdam, 575-584 
+</P>
+<P>Hofierka J., Huld T., Cebecauer T., Suri M., 2007. Open Source Solar 
+Radiation Tools for Environmental and Renewable Energy Applications,
+<a href="http://www.isess.org/papers.asp?year=2007">International Symposium on 
+Environmental Software Systems</a>, Prague, 2007</P>
+<P>Neteler M., Mitasova H., 2004. Open Source GIS: A GRASS GIS
+Approach, <A HREF="http://www.springer.com/geography/gis+cartography/book/978-0-387-35767-6">Springer</A>, New York.
+ISBN: 1-4020-8064-6, 2nd Edition 2004 (reprinted 2005), 424 pages 
+</P>
+<P>Project <A HREF="http://re.jrc.ec.europa.eu/pvgis/">PVGIS</A>, European 
+Commission, DG Joint Research Centre 2001-2007</P>
+<P>Suri M., Hofierka J., 2004.
+A New GIS-based Solar Radiation Model and Its Application for
+Photovoltaic Assessments. <A HREF="http://www.blackwellpublishing.com/toc.asp?ref=1361-1682">Transactions
+in GIS</A>, 8(2), 175-190</P>
+
+<H2>AUTHORS</H2>
+Thomas Huld, Joint Research Centre of
+the European Commission, Ispra, Italy 
+<br>
+Tomas Cebecauer, Joint Research Centre
+of the European Commission, Ispra, Italy 
+<br>
+Jaroslav Hofierka, GeoModel s.r.o.,
+Bratislava, Slovakia <BR>Marcel Suri, Joint Research Centre of the
+European Commission, Ispra, Italy</P>
+&copy; 2007, Thomas Huld, Tomas Cebecauer, Jaroslav Hofierka, Marcel Suri 
+</P>
+
+<ADDRESS STYLE="margin-bottom: 0.2in"><A HREF="mailto:Thomas.Huld at jrc.it">Thomas.Huld at jrc.it</A>
+<A HREF="mailto:Tomas.Cebecauer at jrc.it">Tomas.Cebecauer at jrc.it</A>
+<A HREF="mailto:hofi at geomodel.sk">hofierka at geomodel.sk</A>
+<A HREF="mailto:Marcel.Suri at jrc.it">Marcel.Suri at jrc.it</A> 
+</ADDRESS>
+
+<P><I>Last changed: $Date: 2007/05/16 16:22:04 $</I> 

Added: grass/branches/develbranch_6/raster/r.horizon/local_proto.h
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/local_proto.h	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.horizon/local_proto.h	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,23 @@
+/* sun11.c */
+int INPUT(void);
+int OUTGR(void);
+double amax1(double, double);
+double amin1(double, double);
+int min(int, int);
+int max(int, int);
+void com_par(void);
+double lumcline2(void);
+double joules2(void);
+int quadrant(void);
+double coef_of_line(void);
+void new_point_x(int, double *, double *, double *);
+void new_point_y(int, double *, double *, double *);
+void which_one(double, double, double, double, double, double);
+int combine_x(int, int, int, int);
+int combine_y(int, int, int, int);
+int vertex(int, int);
+int mesh_vertex(void);
+int mesh_line(void);
+void calculate(void);
+double com_sol_const(void);
+double com_declin(int);


Property changes on: grass/branches/develbranch_6/raster/r.horizon/local_proto.h
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.horizon/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/main.c	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.horizon/main.c	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,1270 @@
+/*******************************************************************************
+r.horizon: This module does one of two things:
+
+1) Using a map of the terrain elevation it calculates for a set of points 
+the angle height of the horizon for each point, using an angle interval given
+by the user.
+
+2) For a given minimum angle it calculates one or more raster map giving the mazimum
+distance to a point on the horizon.  
+
+This program was written in 2006 by Tfomas Huld and Tomas Cebecauer, 
+Joint Research Centre of the European Commission, based on bits of the r.sun module by Jaro Hofierka
+
+*******************************************************************************/
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the
+ *   Free Software Foundation, Inc.,
+ *   59 Temple Place - Suite 330,
+ *   Boston, MA  02111-1307, USA.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+
+#define WHOLE_RASTER 1
+#define SINGLE_POINT 0
+#define RAD      (180. / M_PI)
+#define DEG      ((M_PI)/180.)
+#define EARTHRADIUS 6371000.
+#define UNDEF    0.		/* undefined value for terrain aspect */
+#define UNDEFZ   -9999.		/* undefined value for elevation */
+#define BIG      1.e20
+#define SMALL    1.e-20
+#define EPS      1.e-4
+#define DIST     "1.0"
+#define DEGREEINMETERS 111120.
+#define TANMINANGLE 0.008727	/* tan of minimum horizon angle (0.5 deg) */
+
+#define AMAX1(arg1, arg2) ((arg1) >= (arg2) ? (arg1) : (arg2))
+#define DISTANCE1(x1, x2, y1, y2) (sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2)))
+
+
+FILE *fw;
+
+const double rad2deg = 180. / M_PI;
+const double deg2rad = M_PI / 180.;
+const double pihalf = M_PI * 0.5;
+const double twopi = 2. * M_PI;
+const double invEarth = 1. / EARTHRADIUS;
+
+const double minAngle = DEG;
+
+char *elevin;
+char *latin = NULL;
+char *horizon = NULL;
+char *mapset = NULL;
+char *per;
+char shad_filename[256];
+
+struct Cell_head cellhd;
+struct Key_value *in_proj_info, *in_unit_info;
+struct pj_info iproj;
+struct pj_info oproj;
+
+struct Cell_head new_cellhd;
+double bufferZone = 0., ebufferZone = 0., wbufferZone = 0.,
+       nbufferZone = 0., sbufferZone = 0.;
+
+int INPUT(void);
+int OUTGR(int numrows, int numcols);
+double amax1(double, double);
+double amin1(double, double);
+int min(int, int);
+int max(int, int);
+void com_par(double angle);
+int is_shadow(void);
+double horizon_height(void);
+void calculate_shadow();
+double calculate_shadow_onedirection(double shadow_angle);
+
+int new_point();
+double searching();
+int test_low_res();
+
+/*void where_is_point();
+   void cube(int, int);
+ */
+
+void calculate(double xcoord, double ycoord, int buffer_e, int buffer_w,
+	       int buffer_s, int buffer_n);
+
+
+int ip, jp, ip100, jp100;
+int n, m, m100, n100;
+int degreeOutput = 0;
+float **z, **z100, **horizon_raster;
+double stepx, stepy, stepxhalf, stepyhalf, stepxy, xp, yp, op, dp, xg0, xx0,
+    yg0, yy0, deltx, delty;
+double invstepx, invstepy, distxy;
+double offsetx, offsety;
+double single_direction;
+
+/*int arrayNumInt; */
+double xmin, xmax, ymin, ymax, zmax = 0.;
+int d, day, tien = 0;
+
+double length, maxlength = BIG, zmult = 1.0, step = 0.0, dist;
+double fixedMaxLength = BIG;
+char *tt, *lt;
+double z_orig, zp;
+double h0, tanh0, angle;
+double stepsinangle, stepcosangle, sinangle, cosangle, distsinangle,
+    distcosangle;
+double TOLER;
+
+int mode;
+int isMode()
+{
+    return mode;
+}
+void setMode(int val)
+{
+    mode = val;
+}
+
+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 xcoord, ycoord;
+
+
+    struct GModule *module;
+    struct
+    {
+	struct Option *elevin, *dist, *coord, *direction, *horizon, *step,
+	    *bufferzone, *e_buff, *w_buff, *n_buff, *s_buff, *maxdistance;
+    } parm;
+
+    struct
+    {
+	struct Flag *degreeOutput;
+    }
+    flag;
+
+    G_gisinit(argv[0]);
+    module = G_define_module();
+
+    module->keywords = _("raster");
+    module->label =
+	_("Horizon angle computation from a digital elevation model.");
+    module->description =
+	_("Computes horizon angle height from a digital elevation model. The module has two"
+	 " different modes of operation:  "
+	 "1. Computes the entire horizon around a single point whose coordinates are"
+	 " given with the 'coord' option. The horizon height (in radians). "
+	 "2. Computes one or more raster maps of the horizon height in a single direction. "
+	 " The input for this is the angle (in degrees), which is measured "
+	 " counterclockwise with east=0, north=90 etc. The output is the horizon height in radians.");
+
+    if (G_get_set_window(&cellhd) == -1)
+	G_fatal_error("G_get_set_window() failed");
+    stepx = cellhd.ew_res;
+    stepy = cellhd.ns_res;
+    stepxhalf = stepx / 2.;
+    stepyhalf = stepy / 2.;
+    invstepx = 1. / stepx;
+    invstepy = 1. / stepy;
+    /*
+       offsetx = 2. *  invstepx;
+       offsety = 2. * invstepy;
+       offsetx = 0.5*stepx;
+       offsety = 0.5*stepy;
+     */
+    offsetx = 0.5;
+    offsety = 0.5;
+
+    n /*n_cols */  = cellhd.cols;
+    m /*n_rows */  = cellhd.rows;
+
+    n100 = ceil(n / 100.);
+    m100 = ceil(m / 100.);
+
+    xmin = cellhd.west;
+    ymin = cellhd.south;
+    xmax = cellhd.east;
+    ymax = cellhd.north;
+    deltx = fabs(cellhd.east - cellhd.west);
+    delty = fabs(cellhd.north - cellhd.south);
+
+    parm.elevin = G_define_option();
+    parm.elevin->key = "elevin";
+    parm.elevin->type = TYPE_STRING;
+    parm.elevin->required = YES;
+    parm.elevin->gisprompt = "old,cell,raster";
+    parm.elevin->description =
+	_("Name of the input elevation raster map [meters]");
+    parm.elevin->guisection = _("Input_options");
+
+
+    parm.direction = G_define_option();
+    parm.direction->key = "direction";
+    parm.direction->type = TYPE_DOUBLE;
+    parm.direction->required = NO;
+    parm.direction->description =
+	_("Direction in which you want to know the horizon height");
+    parm.direction->guisection = _("Input_options");
+
+    parm.step = G_define_option();
+    parm.step->key = "step";
+    parm.step->type = TYPE_DOUBLE;
+    parm.step->required = NO;
+    parm.step->description =
+	_("For multidirectional horizon: the step size in degrees");
+    parm.step->guisection = _("Input_options");
+
+    parm.bufferzone = G_define_option();
+    parm.bufferzone->key = "bufferzone";
+    parm.bufferzone->type = TYPE_DOUBLE;
+    parm.bufferzone->required = NO;
+    parm.bufferzone->description =
+	_("For horizon rasters, read from the DEM an extra buffer around the present region");
+    parm.bufferzone->guisection = _("Input_options");
+
+    parm.e_buff = G_define_option();
+    parm.e_buff->key = "e_buff";
+    parm.e_buff->type = TYPE_DOUBLE;
+    parm.e_buff->required = NO;
+    parm.e_buff->description =
+	_("For horizon rasters, read from the DEM an extra buffer eastward the present region");
+    parm.e_buff->guisection = _("Input_options");
+
+    parm.w_buff = G_define_option();
+    parm.w_buff->key = "w_buff";
+    parm.w_buff->type = TYPE_DOUBLE;
+    parm.w_buff->required = NO;
+    parm.w_buff->description =
+	_("For horizon rasters, read from the DEM an extra buffer westward the present region");
+    parm.w_buff->guisection = _("Input_options");
+
+    parm.n_buff = G_define_option();
+    parm.n_buff->key = "n_buff";
+    parm.n_buff->type = TYPE_DOUBLE;
+    parm.n_buff->required = NO;
+    parm.n_buff->description =
+	_("For horizon rasters, read from the DEM an extra buffer northward the present region");
+    parm.n_buff->guisection = _("Input_options");
+
+    parm.s_buff = G_define_option();
+    parm.s_buff->key = "s_buff";
+    parm.s_buff->type = TYPE_DOUBLE;
+    parm.s_buff->required = NO;
+    parm.s_buff->description =
+	_("For horizon rasters, read from the DEM an extra buffer southward the present region");
+    parm.s_buff->guisection = _("Input_options");
+
+    parm.maxdistance = G_define_option();
+    parm.maxdistance->key = "maxdistance";
+    parm.maxdistance->type = TYPE_DOUBLE;
+    parm.maxdistance->required = NO;
+    parm.maxdistance->description =
+	_("The maximum distance to consider when finding the horizon height");
+    parm.maxdistance->guisection = _("Input_options");
+
+
+    parm.horizon = G_define_option();
+    parm.horizon->key = "horizon";
+    parm.horizon->type = TYPE_STRING;
+    parm.horizon->required = NO;
+    parm.horizon->gisprompt = "old,cell,raster";
+    parm.horizon->description = _("Prefix of the horizon raster output maps");
+    parm.horizon->guisection = _("Output_options");
+
+
+    parm.coord = G_define_option();
+    parm.coord->key = "coord";
+    parm.coord->type = TYPE_DOUBLE;
+    parm.coord->key_desc = "east,north";
+    parm.coord->required = NO;
+    parm.coord->description =
+	_("Coordinate for which you want to calculate the horizon");
+    parm.coord->guisection = _("Output_options");
+
+    parm.dist = G_define_option();
+    parm.dist->key = "dist";
+    parm.dist->type = TYPE_DOUBLE;
+    parm.dist->answer = DIST;
+    parm.dist->required = NO;
+    parm.dist->description = _("Sampling distance step coefficient (0.5-1.5)");
+    parm.dist->guisection = _("Output_options");
+
+
+    flag.degreeOutput = G_define_flag();
+    flag.degreeOutput->key = 'd';
+    flag.degreeOutput->description =
+	_("Write output in degrees (default is radians)");
+
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    degreeOutput = flag.degreeOutput->answer;
+
+
+    elevin = parm.elevin->answer;
+
+    if (parm.coord->answer == NULL) {
+	setMode(WHOLE_RASTER);
+    }
+    else {
+	setMode(SINGLE_POINT);
+	if (sscanf(parm.coord->answer, "%lf,%lf", &xcoord, &ycoord) != 2) {
+	    G_fatal_error
+		("Can't read the coordinates from the \"coord\" option.");
+
+	}
+
+	/* Transform the coordinates to row/column */
+
+	/*
+	   xcoord = (int) ((xcoord-xmin)/stepx);
+	   ycoord = (int) ((ycoord-ymin)/stepy);
+	 */
+    }
+
+    if (parm.direction->answer != NULL)
+	sscanf(parm.direction->answer, "%lf", &single_direction);
+
+
+    if (isMode(WHOLE_RASTER)) {
+	if ((parm.direction->answer == NULL) && (parm.step->answer == NULL)) {
+	    G_fatal_error
+		(_("You didn't specify a direction value or step size. Aborting."));
+	}
+
+
+	if (parm.horizon->answer == NULL) {
+	    G_fatal_error
+		(_("You didn't specify a horizon raster name. Aborting."));
+
+	}
+	horizon = parm.horizon->answer;
+	if (parm.step->answer != NULL)
+	    sscanf(parm.step->answer, "%lf", &step);
+    }
+    else {
+
+	if (parm.step->answer == NULL) {
+	    G_fatal_error
+		(_("You didn't specify an angle step size. Aborting."));
+
+	}
+	sscanf(parm.step->answer, "%lf", &step);
+
+
+    }
+
+    if (step == 0.0) {
+	step = 360.;
+    }
+
+    if (parm.bufferzone->answer != NULL) {
+	if (sscanf(parm.bufferzone->answer, "%lf", &bufferZone) != 1) {
+	    G_fatal_error(_("Could not read bufferzone size. Aborting."));
+	}
+    }
+
+    if (parm.e_buff->answer != NULL) {
+	if (sscanf(parm.e_buff->answer, "%lf", &ebufferZone) != 1) {
+	    G_fatal_error(_("Could not read east bufferzone size. Aborting."));
+	}
+    }
+
+    if (parm.w_buff->answer != NULL) {
+	if (sscanf(parm.w_buff->answer, "%lf", &wbufferZone) != 1) {
+	    G_fatal_error(_("Could not read west bufferzone size. Aborting."));
+	}
+    }
+
+    if (parm.s_buff->answer != NULL) {
+	if (sscanf(parm.s_buff->answer, "%lf", &sbufferZone) != 1) {
+	    G_fatal_error
+		(_("Could not read south bufferzone size. Aborting."));
+	}
+    }
+
+
+
+    if (parm.n_buff->answer != NULL) {
+	if (sscanf(parm.n_buff->answer, "%lf", &nbufferZone) != 1) {
+	    G_fatal_error
+		(_("Could not read north bufferzone size. Aborting."));
+	}
+    }
+
+    if (parm.maxdistance->answer != NULL) {
+	if (sscanf(parm.maxdistance->answer, "%lf", &fixedMaxLength) != 1) {
+	    G_fatal_error(_("Could not read maximum distance. Aborting."));
+	}
+    }
+
+
+    sscanf(parm.dist->answer, "%lf", &dist);
+
+    stepxy = dist * 0.5 * (stepx + stepy);
+    distxy = dist;
+    TOLER = stepxy * EPS;
+
+    if (bufferZone > 0. || ebufferZone > 0. || wbufferZone > 0. ||
+	sbufferZone > 0. || nbufferZone > 0.) {
+	new_cellhd = cellhd;
+
+	if (ebufferZone == 0.)
+	    ebufferZone = bufferZone;
+	if (wbufferZone == 0.)
+	    wbufferZone = bufferZone;
+	if (sbufferZone == 0.)
+	    sbufferZone = bufferZone;
+	if (nbufferZone == 0.)
+	    nbufferZone = bufferZone;
+
+	new_cellhd.rows += (int)((nbufferZone + sbufferZone) / stepy);
+	new_cellhd.cols += (int)((ebufferZone + wbufferZone) / stepx);
+
+	new_cellhd.north += nbufferZone;
+	new_cellhd.south -= sbufferZone;
+	new_cellhd.east += ebufferZone;
+	new_cellhd.west -= wbufferZone;
+
+	xmin = new_cellhd.west;
+	ymin = new_cellhd.south;
+	xmax = new_cellhd.east;
+	ymax = new_cellhd.north;
+	deltx = fabs(new_cellhd.east - new_cellhd.west);
+	delty = fabs(new_cellhd.north - new_cellhd.south);
+
+	n /*n_cols */  = new_cellhd.cols;
+	m /*n_rows */  = new_cellhd.rows;
+	/*G_debug(3,"%lf %lf %lf %lf \n",ymax, ymin, xmin,xmax);
+	 */
+	n100 = ceil(n / 100.);
+	m100 = ceil(m / 100.);
+
+	if (G_set_window(&new_cellhd) == -1)
+	    exit(EXIT_FAILURE);
+    }
+
+    struct Key_Value *in_proj_info, *in_unit_info;
+
+    if ((in_proj_info = G_get_projinfo()) == NULL)
+	G_fatal_error
+	    (_("Can't get projection info of current location: please set latitude via 'lat' or 'latin' option!"));
+
+    if ((in_unit_info = G_get_projunits()) == NULL)
+	G_fatal_error(_("Can't get projection units of current location"));
+
+    if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
+	G_fatal_error
+	    (_("Can't get projection key values of current location"));
+
+    G_free_key_value(in_proj_info);
+    G_free_key_value(in_unit_info);
+
+    /* Set output projection to latlong w/ same ellipsoid */
+    oproj.zone = 0;
+    oproj.meters = 1.;
+    sprintf(oproj.proj, "ll");
+    if ((oproj.pj = pj_latlong_from_proj(iproj.pj)) == NULL)
+	G_fatal_error(_("Unable to set up lat/long projection parameters"));
+
+
+
+
+/**********end of parser - ******************************/
+
+
+
+    INPUT();
+    calculate(xcoord, ycoord, (int)(ebufferZone / stepx),
+	      (int)(wbufferZone / stepx), (int)(sbufferZone / stepy),
+	      (int)(nbufferZone / stepy));
+
+    if (bufferZone > 0.) {
+	/* Set the region window back to the original */
+	if (G_set_window(&cellhd) == -1)
+	    exit(EXIT_FAILURE);
+    }
+
+    /*sorry, I've moved OUTGR to calculate() - into the loop */
+    /*      if(isMode(WHOLE_RASTER))
+       {
+       OUTGR(cellhd.rows,cellhd.cols);
+       }
+     */
+    if (G_set_window(&cellhd) == -1)
+	exit(EXIT_FAILURE);
+
+    exit(EXIT_SUCCESS);
+}
+
+/**********************end of main.c*****************/
+
+int INPUT(void)
+{
+    FCELL *cell1;
+    int fd1, row, row_rev;
+    int l, i, j, k;
+    int lmax, kmax;
+
+    cell1 = G_allocate_f_raster_buf();
+
+    z = (float **)G_malloc(sizeof(float *) * (m));
+    z100 = (float **)G_malloc(sizeof(float *) * (m100));
+
+    for (l = 0; l < m; l++) {
+	z[l] = (float *)G_malloc(sizeof(float) * (n));
+    }
+    for (l = 0; l < m100; l++) {
+	z100[l] = (float *)G_malloc(sizeof(float) * (n100));
+    }
+    /*read Z raster */
+
+    if ((mapset = G_find_cell(elevin, "")) == NULL)
+	G_fatal_error(_("Raster map <%s> not found"), elevin);
+
+    fd1 = G_open_cell_old(elevin, mapset);
+
+    for (row = 0; row < m; row++) {
+	G_get_f_raster_row(fd1, cell1, row);
+
+	for (j = 0; j < n; j++) {
+	    row_rev = m - row - 1;
+
+	    if (!G_is_f_null_value(cell1 + j))
+		z[row_rev][j] = (float)cell1[j];
+	    else
+		z[row_rev][j] = UNDEFZ;
+
+	}
+    }
+    G_close_cell(fd1);
+
+    /*create low resolution array 100 */
+    for (i = 0; i < m100; i++) {
+	lmax = (i + 1) * 100;
+	if (lmax > m)
+	    lmax = m;
+
+	for (j = 0; j < n100; j++) {
+	    zmax = SMALL;
+	    kmax = (j + 1) * 100;
+	    if (kmax > n)
+		kmax = n;
+	    for (l = (i * 100); l < lmax; l++) {
+		for (k = (j * 100); k < kmax; k++) {
+		    zmax = amax1(zmax, z[l][k]);
+		}
+	    }
+	    z100[i][j] = zmax;
+	    /* G_debug(3,"%d %d %lf\n", i, j, z100[i][j]); */
+	}
+    }
+
+
+    /*find max Z */
+    for (i = 0; i < m; i++) {
+	for (j = 0; j < n; j++) {
+	    zmax = amax1(zmax, z[i][j]);
+
+	}
+    }
+
+    return 1;
+}
+
+
+
+
+int OUTGR(int numrows, int numcols)
+{
+    FCELL *cell1 = NULL;
+
+    int fd1 = 0;
+    int i, iarc, j;
+
+    if (G_set_window(&cellhd) < 0)
+	exit(EXIT_FAILURE);
+
+    if (horizon != NULL) {
+	cell1 = G_allocate_f_raster_buf();
+	fd1 = G_open_fp_cell_new(shad_filename);
+	if (fd1 < 0)
+	    G_fatal_error(_("Unable to create raster map %s"), shad_filename);
+    }
+
+
+    if (numrows != G_window_rows())
+	G_fatal_error(_("OOPS: rows changed from %d to %d"), numrows,
+		      G_window_rows());
+
+    if (numcols != G_window_cols())
+	G_fatal_error(_("OOPS: cols changed from %d to %d"), numcols,
+		      G_window_cols());
+
+    for (iarc = 0; iarc < numrows; iarc++) {
+	i = numrows - iarc - 1;
+
+	if (horizon != NULL) {
+	    for (j = 0; j < numcols; j++) {
+		if (horizon_raster[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell1 + j, 1);
+		else
+		    cell1[j] = (FCELL) horizon_raster[i][j];
+	    }
+	    G_put_f_raster_row(fd1, cell1);
+	}
+
+    }				/* End loop over rows. */
+
+    G_close_cell(fd1);
+
+
+
+    return 1;
+}
+
+
+double amax1(arg1, arg2)
+     double arg1;
+     double arg2;
+{
+    double res;
+
+    if (arg1 >= arg2) {
+	res = arg1;
+    }
+    else {
+	res = arg2;
+    }
+    return res;
+}
+
+double amin1(arg1, arg2)
+     double arg1;
+     double arg2;
+{
+    double res;
+
+    if (arg1 <= arg2) {
+	res = arg1;
+    }
+    else {
+	res = arg2;
+    }
+    return res;
+}
+
+
+int min(arg1, arg2)
+     int arg1;
+     int arg2;
+{
+    int res;
+
+    if (arg1 <= arg2) {
+	res = arg1;
+    }
+    else {
+	res = arg2;
+    }
+    return res;
+}
+
+int max(arg1, arg2)
+     int arg1;
+     int arg2;
+{
+    int res;
+
+    if (arg1 >= arg2) {
+	res = arg1;
+    }
+    else {
+	res = arg2;
+    }
+    return res;
+}
+
+
+
+/**********************************************************/
+
+void com_par(double angle)
+{
+    sinangle = sin(angle);
+    if (fabs(sinangle) < 0.0000001) {
+	sinangle = 0.;
+    }
+    cosangle = cos(angle);
+    if (fabs(cosangle) < 0.0000001) {
+	cosangle = 0.;
+    }
+    distsinangle = 32000;
+    distcosangle = 32000;
+
+    if (sinangle != 0.) {
+	distsinangle = 100. / (distxy * sinangle);
+    }
+    if (cosangle != 0.) {
+	distcosangle = 100. / (distxy * cosangle);
+    }
+
+    stepsinangle = stepxy * sinangle;
+    stepcosangle = stepxy * cosangle;
+}
+
+double horizon_height(void)
+{
+    double height;
+
+    tanh0 = 0.;
+    length = 0;
+
+    height = searching();
+
+    xx0 = xg0;
+    yy0 = yg0;
+
+    return height;
+}
+
+
+double calculate_shadow_onedirection(double shadow_angle)
+{
+
+    shadow_angle = horizon_height();
+
+    return shadow_angle;
+}
+
+
+
+void calculate_shadow()
+{
+    double dfr_rad;
+
+    int i;
+    int printCount;
+    double shadow_angle;
+    double printangle;
+    double sx, sy;
+    double xp, yp;
+    double latitude, longitude;
+    double delt_lat, delt_lon;
+    double delt_east, delt_nor;
+    double delt_dist;
+
+    double angle;
+
+    printCount = 360. / fabs(step);
+
+
+    if (printCount < 1)
+	printCount = 1;
+
+
+    dfr_rad = step * deg2rad;
+
+    xp = xmin + xx0;
+    yp = ymin + yy0;
+
+    angle = (single_direction * deg2rad) + pihalf;
+
+
+    maxlength = fixedMaxLength;
+
+    for (i = 0; i < printCount; i++) {
+
+	ip = jp = 0;
+
+
+	sx = xx0 * invstepx;
+	sy = yy0 * invstepy;
+	ip100 = floor(sx / 100.);
+	jp100 = floor(sy / 100.);
+
+
+	if ((G_projection() != PROJECTION_LL)) {
+
+	    longitude = xp;
+	    latitude = yp;
+
+	    if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+		G_fatal_error(_("Error in pj_do_proj"));
+	    }
+	}
+	else {			/* ll projection */
+	    latitude = yp;
+	    longitude = xp;
+	}
+	latitude *= deg2rad;
+	longitude *= deg2rad;
+
+
+	delt_lat = -0.0001 * cos(angle);
+	delt_lon = 0.0001 * sin(angle) / cos(latitude);
+
+	latitude = (latitude + delt_lat) * rad2deg;
+	longitude = (longitude + delt_lon) * rad2deg;
+
+	if (pj_do_proj(&longitude, &latitude, &oproj, &iproj) < 0) {
+	    G_fatal_error(_("Error in pj_do_proj"));
+	}
+
+	delt_east = longitude - xp;
+	delt_nor = latitude - yp;
+
+	delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor);
+
+	stepsinangle = stepxy * delt_nor / delt_dist;
+	stepcosangle = stepxy * delt_east / delt_dist;
+
+
+	shadow_angle = horizon_height();
+
+
+	if (degreeOutput) {
+	    shadow_angle *= rad2deg;
+	}
+	printangle = angle * rad2deg - 90.;
+	if (printangle < 0.)
+	    printangle += 360;
+	else if (printangle >= 360.)
+	    printangle -= 360;
+
+	G_message("%lf, %lf", printangle, shadow_angle);
+
+	angle += dfr_rad;
+
+	if (angle < 0.)
+	    angle += twopi;
+	else if (angle > twopi)
+	    angle -= twopi;
+    }				/* end of for loop over angles */
+}
+
+/*////////////////////////////////////////////////////////////////////// */
+
+
+int new_point()
+{
+    int iold, jold;
+    int succes = 1, succes2 = 1;
+    double sx, sy;
+    double dx, dy;
+
+    iold = ip;
+    jold = jp;
+
+    while (succes) {
+	yy0 += stepsinangle;
+	xx0 += stepcosangle;
+
+
+	/* offset 0.5 cell size to get the right cell i, j */
+	sx = xx0 * invstepx + offsetx;
+	sy = yy0 * invstepy + offsety;
+	ip = (int)sx;
+	jp = (int)sy;
+
+	/* test outside of raster */
+	if ((ip < 0) || (ip >= n) || (jp < 0) || (jp >= m))
+	    return (3);
+
+	if ((ip != iold) || (jp != jold)) {
+	    dx = (double)ip *stepx;
+	    dy = (double)jp *stepy;
+
+	    length = distance(xg0, dx, yg0, dy);  /* dist from orig. grid point to the current grid point */
+	    succes2 = test_low_res();
+	    if (succes2 == 1) {
+		zp = z[jp][ip];
+		return (1);
+	    }
+	}
+    }
+    return -1;
+}
+
+
+int test_low_res()
+{
+    int iold100, jold100;
+    double sx, sy;
+    int delx, dely, mindel;
+    double zp100, z2, curvature_diff;
+
+    iold100 = ip100;
+    jold100 = jp100;
+    ip100 = floor(ip / 100.);
+    jp100 = floor(jp / 100.);
+    /*test the new position with low resolution */
+    if ((ip100 != iold100) || (jp100 != jold100)) {
+	/* G_debug(3,"%d %d %d %d\n",ip,jp, iold100,jold100); */
+	/*  replace with approximate version
+	   curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
+	 */
+	curvature_diff = 0.5 * length * length * invEarth;
+	z2 = z_orig + curvature_diff + length * tanh0;
+	zp100 = z100[jp100][ip100];
+	/*G_debug(3,"%d %d %lf %lf \n",ip,jp,z2,zp100); */
+
+	if (zp100 <= z2)
+	    /*skip to the next lowres cell */
+	{
+	    delx = 32000;
+	    dely = 32000;
+	    if (cosangle > 0.) {
+		sx = xx0 * invstepx + offsetx;
+		delx =
+		    floor(fabs
+			  ((ceil(sx / 100.) - (sx / 100.)) * distcosangle));
+	    }
+	    if (cosangle < 0.) {
+		sx = xx0 * invstepx + offsetx;
+		delx =
+		    floor(fabs
+			  ((floor(sx / 100.) - (sx / 100.)) * distcosangle));
+	    }
+	    if (sinangle > 0.) {
+		sy = yy0 * invstepy + offsety;
+		dely =
+		    floor(fabs
+			  ((ceil(sy / 100.) - (sy / 100.)) * distsinangle));
+	    }
+	    else if (sinangle < 0.) {
+		sy = yy0 * invstepy + offsety;
+		dely =
+		    floor(fabs
+			  ((floor(jp / 100.) - (sy / 100.)) * distsinangle));
+	    }
+
+	    mindel = min(delx, dely);
+	    /*G_debug(3,"%d %d %d %lf %lf\n",ip, jp, mindel,xg0, yg0);*/
+
+	    yy0 = yy0 + (mindel * stepsinangle);
+	    xx0 = xx0 + (mindel * stepcosangle);
+	    /*G_debug(3,"  %lf %lf\n",xx0,yy0);*/
+
+	    return (3);
+	}
+	else {
+	    return (1);		/*change of low res array - new cell is reaching limit for high resolution processing */
+	}
+    }
+    else {
+	return (1);		/*no change of low res array */
+    }
+}
+
+
+double searching()
+{
+    double z2;
+    double curvature_diff;
+    int succes = 1;
+
+    if (zp == UNDEFZ)
+	return 0;
+
+    while (1) {
+	succes = new_point();
+
+	if (succes != 1) {
+	    break;
+	}
+	/*
+	   curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
+	 */
+	curvature_diff = 0.5 * length * length * invEarth;
+
+	z2 = z_orig + curvature_diff + length * tanh0;
+
+	if (z2 < zp) {
+	    tanh0 = (zp - z_orig - curvature_diff) / length;
+	}
+
+
+	if (z2 >= zmax) {
+	    break;
+	}
+
+	if (length >= maxlength) {
+	    break;
+	}
+
+    }
+
+    return atan(tanh0);
+}
+
+
+
+/*////////////////////////////////////////////////////////////////////// */
+
+void calculate(double xcoord, double ycoord, int buffer_e, int buffer_w,
+	       int buffer_s, int buffer_n)
+{
+    int i, j, l, k = 0;
+    int numDigits;
+
+    int xindex, yindex;
+    double shadow_angle;
+    double coslat;
+
+    double latitude, longitude;
+    double xp, yp;
+    double inputAngle;
+    double delt_lat, delt_lon;
+    double delt_east, delt_nor;
+    double delt_dist;
+
+    char formatString[10];
+
+
+    int hor_row_start = buffer_s;
+    int hor_row_end = m - buffer_n;
+
+    int hor_col_start = buffer_w;
+    int hor_col_end = n - buffer_e;
+
+    int hor_numrows = m - (buffer_s + buffer_n);
+    int hor_numcols = n - (buffer_e + buffer_w);
+
+    /*      char shad_filename[256]; */
+    int arrayNumInt;
+    double dfr_rad;
+
+    xindex = (int)((xcoord - xmin) / stepx);
+    yindex = (int)((ycoord - ymin) / stepy);
+
+    if ((G_projection() == PROJECTION_LL)) {
+	ll_correction = 1;
+    }
+
+
+    if (isMode() == SINGLE_POINT) {
+	/* Calculate the horizon for one single point */
+
+	/* 
+	   xg0 = xx0 = (double)xcoord * stepx;
+	   yg0 = yy0 = (double)ycoord * stepy;
+	   xg0 = xx0 = xcoord -0.5*stepx -xmin;
+	   yg0 = yy0 = ycoord -0.5*stepy-ymin;
+	   xg0 = xx0 = xindex*stepx -0.5*stepx;
+	   yg0 = yy0 = yindex*stepy -0.5*stepy;
+	 */
+	xg0 = xx0 = xindex * stepx;
+	yg0 = yy0 = yindex * stepy;
+
+
+	if (ll_correction) {
+	    coslat = cos(deg2rad * (ymin + yy0));
+	    coslatsq = coslat * coslat;
+	}
+
+	G_debug(3, "yindex: %d, xindex %d", yindex, xindex);
+	z_orig = zp = z[yindex][xindex];
+
+	calculate_shadow();
+
+    }
+    else {
+
+	/****************************************************************/
+	/*  The loop over raster points starts here!                    */
+	/****************************************************************/
+
+	if (horizon != NULL) {
+	    horizon_raster = (float **)G_malloc(sizeof(float *) * (hor_numrows));
+	    for (l = 0; l < hor_numrows; l++) {
+		horizon_raster[l] =
+		    (float *)G_malloc(sizeof(float) * (hor_numcols));
+	    }
+
+	    for (j = 0; j < hor_numrows; j++) {
+		for (i = 0; i < hor_numcols; i++)
+		    horizon_raster[j][i] = 0.;
+	    }
+	}
+	/*
+	   definition of horizon angle in loop
+	 */
+	if (step == 0.0) {
+	    dfr_rad = 0;
+	    arrayNumInt = 1;
+	    sprintf(shad_filename, "%s", horizon);
+	}
+	else {
+	    dfr_rad = step * deg2rad;
+	    arrayNumInt = (int)(360. / fabs(step));
+	}
+
+	numDigits = (int)(log10(1. * arrayNumInt)) + 1;
+	sprintf(formatString, "%%s_%%0%dd", numDigits);
+
+	for (k = 0; k < arrayNumInt; k++) {
+	    if (step != 0.0)
+		sprintf(shad_filename, formatString, horizon, k);
+	    angle = (single_direction * deg2rad) + (dfr_rad * k);
+	    /*              
+	       com_par(angle);
+	     */
+	    G_message(_("Calculating map %01d of %01d (angle %lf, raster map <%s>)"), (k + 1), arrayNumInt,
+		   angle * rad2deg, shad_filename);
+
+	    for (j = hor_row_start; j < hor_row_end; j++) {
+		G_percent(j - hor_row_start, hor_numrows - 1, 2);
+		shadow_angle = 15 * deg2rad;
+		for (i = hor_col_start; i < hor_col_end; i++) {
+		    ip100 = floor(i / 100.);
+		    jp100 = floor(j / 100.);
+		    ip = jp = 0;
+		    xg0 = xx0 = (double)i *stepx;
+
+		    xp = xmin + xx0;
+		    yg0 = yy0 = (double)j *stepy;
+
+		    yp = ymin + yy0;
+		    length = 0;
+		    if (ll_correction) {
+			coslat = cos(deg2rad * yp);
+			coslatsq = coslat * coslat;
+		    }
+
+		    longitude = xp;
+		    latitude = yp;
+
+
+		    if ((G_projection() != PROJECTION_LL)) {
+
+			if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) <	0) {
+			    G_fatal_error("Error in pj_do_proj");
+			}
+		    }
+		    latitude *= deg2rad;
+		    longitude *= deg2rad;
+
+		    inputAngle = angle + pihalf;
+		    inputAngle =
+			(inputAngle >=
+			 twopi) ? inputAngle - twopi : inputAngle;
+
+
+		    delt_lat = -0.0001 * cos(inputAngle);	/* Arbitrary small distance in latitude */
+		    delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
+
+		    latitude = (latitude + delt_lat) * rad2deg;
+		    longitude = (longitude + delt_lon) * rad2deg;
+
+		    if ((G_projection() != PROJECTION_LL)) {
+			if (pj_do_proj(&longitude, &latitude, &oproj, &iproj) < 0) {
+			    G_fatal_error("Error in pj_do_proj");
+			}
+		    }
+
+		    delt_east = longitude - xp;
+		    delt_nor = latitude - yp;
+
+		    delt_dist =
+			sqrt(delt_east * delt_east + delt_nor * delt_nor);
+
+		    sinangle = delt_nor / delt_dist;
+		    if (fabs(sinangle) < 0.0000001) {
+			sinangle = 0.;
+		    }
+		    cosangle = delt_east / delt_dist;
+		    if (fabs(cosangle) < 0.0000001) {
+			cosangle = 0.;
+		    }
+		    distsinangle = 32000;
+		    distcosangle = 32000;
+
+		    if (sinangle != 0.) {
+			distsinangle = 100. / (distxy * sinangle);
+		    }
+		    if (cosangle != 0.) {
+			distcosangle = 100. / (distxy * cosangle);
+		    }
+
+		    stepsinangle = stepxy * sinangle;
+		    stepcosangle = stepxy * cosangle;
+
+
+		    z_orig = zp = z[j][i];
+		    maxlength = (zmax - z_orig) / TANMINANGLE;
+		    maxlength =
+			(maxlength <
+			 fixedMaxLength) ? maxlength : fixedMaxLength;
+		    if (z_orig != UNDEFZ) {
+
+			/*G_debug(3,"**************new line %d %d\n", i, j);
+			 */
+			shadow_angle = horizon_height();
+			if (degreeOutput) {
+			    shadow_angle *= rad2deg;
+			}
+
+			/*
+			   if((j==1400)&&(i==1400))
+			   {
+			   G_debug(3, "coordinates=%f,%f, raster no. %d, horizon=%f\n",
+			   xp, yp, k, shadow_angle);
+			   }
+			 */
+			horizon_raster[j - buffer_s][i - buffer_w] =
+			    shadow_angle;
+
+		    }		/* undefs */
+		}
+	    }
+
+	    OUTGR(cellhd.rows, cellhd.cols);
+
+	    /* empty array */
+	    for (j = 0; j < hor_numrows; j++) {
+		for (i = 0; i < hor_numcols; i++)
+		    horizon_raster[j][i] = 0.;
+	    }
+
+	    /*return back the buffered region */
+	    if (bufferZone > 0.) {
+		if (G_set_window(&new_cellhd) == -1)
+		    exit(0);
+	    }
+	}
+
+    }
+
+}


Property changes on: grass/branches/develbranch_6/raster/r.horizon/main.c
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.sun/DEPRECATED
===================================================================
--- grass/branches/develbranch_6/raster/r.sun/DEPRECATED	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun/DEPRECATED	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,2 @@
+no longer maintained.
+Continued in ../r.sun2/

Added: grass/branches/develbranch_6/raster/r.sun2/Makefile
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/Makefile	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/Makefile	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.sun
+
+LIBES =  $(GPROJLIB) $(GISLIB)
+DEPENDENCIES = $(GPROJDEP) $(GISDEP)
+EXTRA_INC = $(PROJINC)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd


Property changes on: grass/branches/develbranch_6/raster/r.sun2/Makefile
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.sun2/README
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/README	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/README	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,81 @@
+From Thomas.Huld jrc.it Thu Jul 26 11:01:44 2007
+To: Markus Neteler <neteler itc.it>
+CC: Jaro Hofierka <hofierka geomodel.sk>, "marcel.suri jrc.it"
+        <marcel.suri jrc.it>, Tomas Cebecauer <tomas.cebecauer jrc.it>
+Date: Thu, 26 Jul 2007 11:01:44 +0200
+Subject: Re: r.horizon source code
+Lines: 1533
+
+I'm sending the r.sun sources to you here, even though the description is not 
+yet up to date. I think it works OK now, though I haven't tried all possible 
+combinations of options. For instance, I haven't ever used the incidence 
+angle or insolation time outputs. There's no reason why they should have 
+stopped working but there are always changes that break something unexpected.
+
+The main changes  are:
+
+- possibility to read in horizon rasters to use with the shadow calculation. 
+
+- Calculation of the shadowing effect should now also work in lat/lon 
+projection.
+
+- There was a bug in r.sun whereby the shadowing effect was calculated wrong 
+when the left-right and up-down directions in the projection are different 
+from east-west and north-south, as is the case in the projection we use 
+(Lambert Azimuthal) when moving away from the point of origin. This has been 
+fixed now though only checked with Lambert Azimuthal.
+
+- I'm not sure if this went in already in a previous version, but the 
+algorithm now also takes into account the curvature of the earth when 
+calculating shadows.
+
+- new output possibility: glob_rad, which is the sum of the three radiation 
+components
+
+- new input parameter: <i>civiltime</i>. When this parameter is given, the 
+single time calculation will calculate the irradiance at the *same* time for 
+the entire raster instead of using the local solar time. The value of 
+<i>civiltime</i> is the timezone, relative to GMT. (+1 for central Europe)
+
+- new input parameter <i>numpartitions</i> No change in the calculation 
+methods, but the program will read the input rasters and do the calculations 
+in a number of chunks, instead of reading in everything at the start. Output 
+is still only at the very end. This is only to save memory. It will not work 
+if you try to calculate shadows without the horizon information, and the 
+program will tell you so.
+
+- I reorganized the source code quite a lot. In particular, I spun off several 
+functions into a separate source file called rsunlib.c. I also organized many 
+of the global variables into <i>struct</i>'s for an easier overview. I did 
+this because we work here with several derivatives of r.sun:
+
+	* r.sunyear which calculates the optimum inclination angle for maximum 
+insolation
+	* r.pv which calculates the PV output taking into account also temperature 
+effects
+	* r.sun.2axis which is r.sun for a moving plane, whose normal always points 
+to the sun (two-axis tracking solar energy systems).
+
+Splitting the source code made it easier to reuse bits for these r.sun 
+derivatives. I haven't updated the Makefile to account for these changes yet. 
+Look at the little script  "compdebug" in the source directory.
+
+
+There may be other changes that I've forgotten about but nothing dramatic.
+
+Have fun!
+
+
+Thomas
+
+
+
+-- 
+--------------------------------------------------
+Thomas Huld   
+Joint Research Centre of the European Commission
+T.P. 450
+I-21020 Ispra, Italy
+phone: +39 0332785273
+e-mail: Thomas.Huld at jrc.it
+--------------------------------------------------

Added: grass/branches/develbranch_6/raster/r.sun2/TODO
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/TODO	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/TODO	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,10 @@
+######################################################################
+Probably the sun position calculation shouldbe replaced
+with 
+
+ G_calc_solar_position()
+
+currently used in r.sunmask. G_calc_solar_position() is based on
+solpos.c from NREL and very precise.
+
+MN 4/2001

Added: grass/branches/develbranch_6/raster/r.sun2/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/description.html	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/description.html	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,282 @@
+<h2>DESCRIPTION</h2>
+
+<b>r.sun</b> computes beam (direct), diffuse and ground reflected solar
+irradiation raster maps for given day, latitude, surface and atmospheric
+conditions. Solar parameters (e.g. time of sunrise and sunset, declination,
+extraterrestrial irradiance, daylight length) are stored in the resultant maps'
+history files. Alternatively, the local time can be specified to compute solar
+incidence angle and/or irradiance raster maps. The shadowing effect of the
+topography is optionally incorporated. This can be done either by calculating
+the shadowing effect directly from the digital elevation model or using rasters
+of the horizon height which is much faster. The horizon rasters can be
+constructed using  <a href="r.horizon.html">r.horizon</a>.
+<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>
+
+The solar geometry of the model is based on the works of Krcho (1990), later
+improved by Jenco (1992). The equations describing Sun &#8211; Earth position as
+well as an interaction of the solar radiation with atmosphere were originally 
+based on the formulas suggested by Kitler and Mikler (1986). This component 
+was considerably updated by the results and suggestions of the working group
+co-ordinated by Scharmer and Greif (2000) (this algorithm might be replaced
+by SOLPOS algorithm-library included in GRASS within <a href="r.sunmask.html">
+r.sunmask</a>
+ command). The model computes all three components of global radiation (beam, 
+diffuse and reflected) for the clear sky conditions, i.e. not taking into 
+consideration the spatial and temporal variation of clouds. The extent and
+spatial resolution of the modelled area, as well as integration over time,
+are limited only by the memory and data storage resources. The model is built
+to fulfil user needs in various fields of science (hydrology, climatology,
+ecology and environmental sciences, photovoltaics, engineering, etc.) for
+continental, regional up to the landscape scales. 
+<p>As an option the model considers a shadowing effect of the local topography. 
+The r.sun program works in two modes. In the first mode it calculates for the set 
+local time a solar incidence angle [degrees] and solar irradiance values [W.m-2].
+In the second mode daily sums of solar radiation [Wh.m-2.day-1] are computed
+within a set day. By a scripting the two modes can be used separately or
+in a combination to provide estimates for any desired time interval. The
+model accounts for sky obstruction by local relief features. Several solar
+parameters are saved in the resultant maps' history files, which may be viewed
+with the <a href="r.info.html">r.info</a> command.
+</p>
+<p>The solar incidence angle raster map <i>incidout</i> is computed specifying 
+elevation raster map <i>elevin</i>, aspect raster map <i>aspin</i>, slope 
+steepness raster map <i>slopin,</i> given the day <i>day</i> and local time
+<i>time</i>. There is no need to define latitude for locations with known
+and defined projection/coordinate system (check it with the <a href="g.proj.html">
+g.proj</a>
+ command). If you have undefined projection, (x,y) system, etc. then the latitude
+can be defined explicitely for large areas by input raster map <i>latin</i>
+ with interpolated latitude values or, for smaller areas, a single region 
+latitude value <i>lat</i> can be used instead. All input raster maps must
+be floating point (FCELL) raster maps. Null data in maps are excluded from
+the computation (and also speeding-up the computation), so each output raster
+map will contain null data in cells according to all input raster maps. The
+user can use <a href="r.null.html">r.null</a>
+ command to create/reset null file for your input raster maps. <br>
+The specified day <i>day</i> is the number of the day of the general year
+where January 1 is day no.1 and December 31 is 365. Time <i>time</i> must
+be a local (solar) time (i.e. NOT a zone time, e.g. GMT, CET) in decimal system,
+e.g. 7.5 (= 7h 30m A.M.), 16.1 = 4h 6m P.M.. </p>
+<p>Setting the solar declination <i>declin</i> by user is an option to override
+the value computed by the internal routine for the day of the year. The value
+of geographical latitude can be set as a constant for the whole computed
+region or, as an option, a grid representing spatially distributed values
+over a large region. The geographical latitude must be also in decimal system
+with positive values for northern hemisphere and negative for southern one.
+In similar principle the Linke turbidity factor (<i>linkein</i>, <i>lin</i>
+) and ground albedo (<i>albedo</i>, <i>alb</i>) can be set. </p>
+<p>Besides clear-sky radiations, the user can compute a real-sky radiation (beam,
+diffuse) using <i>coefbh</i> and <i>coefdh </i>input raster maps defining
+the fraction of the respective clear-sky radiations reduced by atmospheric
+factors (e.g. cloudiness). The value is between 0-1. Usually these
+coefficients can be obtained from a long-terms meteorological measurements.</p>
+<p>The solar irradiation or irradiance raster maps <i>beam_rad</i>, <i>diff_rad</i>
+, <i>refl_rad</i> are computed for a given day <i>day,</i> latitude <i>lat
+(latin), </i>elevation <i>elevin</i>, slope <i>slopein</i> and aspect <i>
+aspin</i> raster maps. For convenience, the output raster given as <i>glob_rad</i>
+will output the sum of the three radiation components. The program uses 
+the Linke atmosphere turbidity factor and ground albedo coefficient. 
+A default, single value of Linke factor is <i>lin</i>=3.0 and 
+is near the annual average for rural-city areas. The Linke
+factor for an absolutely clear atmosphere is <i>lin</i>=1.0. See notes below
+to learn more about this factor. The incidence solar angle is the angle between
+horizon and solar beam vector. The solar radiation maps for given day are
+computed integrating the relevant irradiance between sunrise and sunset times
+for given day. The user can set finer or coarser time step <i>step</i> used 
+for all-day radiation calculations. A default value of <i>step</i> is 0.5 
+hour. Larger steps (e.g. 1.0-2.0) can speed-up calculations but produce less
+reliable results. The output units are in Wh per squared meter per given
+day [Wh/(m*m)/day]. The incidence angle and irradiance/irradiation maps can
+be computed without shadowing influence of relief by default or they can
+be computed with this influence using the flag <i>-s</i>. In mountainous areas
+this can lead to very different results! The user should be aware that taken
+into account the shadowing effect of relief can slow
+down the speed of computing especially when the sun altitude is low.
+ When considering shadowing effect (flag <i>-s</i>) speed and precision computing
+can be controlled by a parameter <i>dist</i> which defines the sampling density
+at which the visibility of a grid cell is computed in the direction of a
+path of the solar flow. It also defines the method by which the obstacle's
+altitude is computed. When choosing <i>dist</i> less than 1.0 (i.e. sampling
+points will be computed at <i>dist</i> * cellsize distance), r.sun takes
+altitude from the nearest grid point. Values above 1.0 will use the maximum
+altitude value found in the nearest 4 surrounding grid points. The default
+value <i>dist</i>=1.0 should give reasonable results for most cases (e.g.
+on DEM). <i>Dist</i> value defines a multiplying coefficient for sampling
+distance. This basic sampling distance equals to the arithmetic average of
+both cell sizes. The reasonable values are in the range 0.5-1.5.  The values
+below 0.5 will decrease and values above 1.0 will increase the computing
+speed. Values greater than 2.0 may produce estimates with lower accuracy
+in highly dissected relief. The fully shadowed areas are written to the ouput
+maps as zero values. Areas with NULL data are considered as no barrier with
+shadowing effect .</p>
+<p>The maps' history files are generated containing the following listed 
+parameters used in the computation: <br>
+- Solar constant 1367 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>
+- Decimal hour (Alternative 1 only) <br>
+- Sunrise and sunset (min-max) over a horizontal plane <br>
+- Daylight lengths <br>
+- Geographical latitude (min-max) <br>
+- Linke turbidity factor (min-max) <br>
+- Ground albedo (min-max) </p>
+<p>The user can use a nice shellcript with variable
+day to compute radiation for some time interval within the year (e.g. vegetation
+or winter period). Elevation, aspect and slope input values should not be
+reclassified into coarser categories. This could lead to incorrect results. 
+</p>
+
+<h2> OPTIONS</h2>
+<P>Currently, there are two modes of r.sun.
+In the first mode it calculates solar incidence angle and solar irradiance
+raster maps using the set local time. In the second mode daily sums of solar
+irradiation [Wh.m-2.day-1] are computed for a specified day. </p>
+
+<h2>
+NOTES</h2>
+
+Solar energy is an important input parameter in different models concerning 
+energy industry, landscape, vegetation, evapotranspiration, snowmelt or remote
+sensing. Solar rays incidence angle maps can be effectively used in radiometric
+and topographic corrections in mountainous and hilly terrain where very accurate
+investigations should be performed. 
+<p>
+The clear-sky solar radiation model applied in the r.sun is based on the
+work undertaken for development of European Solar Radiation Atlas (Scharmer 
+and Greif 2000, Page et al. 2001, Rigollier 2001). The clear sky model estimates
+the global radiation from the sum of its beam, diffuse and reflected components.
+The main difference between solar radiation models for inclined surfaces
+in Europe is the treatment of the diffuse component. In the European climate
+this component is often the largest source of estimation error. Taking into
+consideration the existing models and their limitation the European Solar
+Radiation Atlas team selected the Muneer (1990) model as it has a sound theoretical
+basis and thus more potential for later improvement. </p>
+<p>
+Details of underlying equations used in this program can be found in the
+reference literature cited below or book published by Neteler and Mitasova: 
+Open Source GIS: A GRASS GIS Approach (published in Kluwer Academic Publishers 
+in 2002). </p>
+<p>
+Average monthly values of the Linke turbidity coefficient for a mild climate
+(see reference literature for your study area): </p>
+<pre>
+Month&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Jan&nbsp; Feb&nbsp; Mar&nbsp; Apr&nbsp; May&nbsp; Jun&nbsp; Jul&nbsp; Aug&nbsp; Sep&nbsp; Oct&nbsp; Nov&nbsp; Dec&nbsp; annual<br>
+mountains&nbsp; 1.5&nbsp; 1.6&nbsp; 1.8&nbsp; 1.9&nbsp; 2.0&nbsp; 2.3&nbsp; 2.3&nbsp; 2.3&nbsp; 2.1&nbsp; 1.8&nbsp; 1.6&nbsp; 1.5&nbsp; 1.90&nbsp;&nbsp;<br>
+rural&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 2.1&nbsp; 2.2&nbsp; 2.5&nbsp; 2.9&nbsp; 3.2&nbsp; 3.4&nbsp; 3.5&nbsp; 3.3&nbsp; 2.9&nbsp; 2.6&nbsp; 2.3&nbsp; 2.2&nbsp; 2.75&nbsp;&nbsp;<br>
+city&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 3.1&nbsp; 3.2&nbsp; 3.5&nbsp; 4.0&nbsp; 4.2&nbsp; 4.3&nbsp; 4.4&nbsp; 4.3&nbsp; 4.0&nbsp; 3.6&nbsp; 3.3&nbsp; 3.1&nbsp; 3.75&nbsp;&nbsp;<br>
+industrial 4.1&nbsp; 4.3&nbsp; 4.7&nbsp; 5.3&nbsp; 5.5&nbsp; 5.7&nbsp; 5.8&nbsp; 5.7&nbsp; 5.3&nbsp; 4.9&nbsp; 4.5&nbsp; 4.2&nbsp; 5.00
+</pre>
+<P>
+Planned improvements include the use of the SOLPOS algorithm for solar 
+geometry calculations and internal computation of aspect and slope.
+
+<h3>Shadow maps</h3>
+A map of shadows can be extracted from the solar incidence angle map
+(incidout). Areas with zero values are shadowed. The <em>-s</em> flag
+has to be used.
+
+<h3>Large maps and out of memory problems</h3>
+
+In case of out of memory error (<tt>ERROR: G_malloc: out of memory</tt>), the
+<em>numpartitions</em> parameter can be used to run a segmented calculation
+which consumes less memory during the computations.
+
+<h3>Example</h3>
+
+Spearfish example (considering also cast shadows):
+<div class="code"><pre>
+g.region rast=elevation.dem -p
+
+# calculate horizons
+# (we put a bufferzone of 10% of maxdistance around the study area)
+r.horizon elevin=elevation.dem step=30 bufferzone=200 horizon=horangle dist=0.7 maxdistance=2000
+
+# slope + aspect
+r.slope.aspect elevation=elevation.dem aspect=aspect.dem slope=slope.dem
+
+# calculate global radiation for day 180 at 14:00hs
+r.sun2 -s elevation.dem horizon=horangle horizonstep=30 aspin=aspect.dem \
+          slopein=slope.dem glob_rad=global_rad day=180 time=14
+</pre></div>
+
+<h2>SEE ALSO</h2>
+<a href="r.slope.aspect.html">r.slope.aspect</a>,
+<a href="r.sunmask.html">r.sunmask</a>,
+<a href="g.proj.html">g.proj</a>,
+<a href="r.null.html">r.null</a>,
+<a href="v.surf.rst.html">v.surf.rst</a>
+
+
+<h2>REFERENCES</h2>
+
+Hofierka, J., Suri, M. (2002): The solar radiation model for Open source
+GIS: implementation and applications. Manuscript submitted to the International
+GRASS users conference in Trento, Italy, September 2002. 
+<p>
+Hofierka, J. (1997). Direct solar radiation modelling within an open GIS
+environment. Proceedings of JEC-GI'97 conference in Vienna, Austria, IOS
+Press Amsterdam, 575-584. </p>
+<p>
+Jenco, M. (1992). Distribution of direct solar radiation on georelief and
+its modelling by means of complex digital model of terrain (in Slovak). Geograficky
+casopis, 44, 342-355. </p>
+<p>
+Kasten, F. (1996). The Linke turbidity factor based on improved values of
+the integral Rayleigh optical thickness. Solar Energy, 56 (3), 239-244. </p>
+<p>
+Kasten, F., Young, A. T. (1989). Revised optical air mass tables and approximation
+formula. Applied Optics, 28, 4735-4738. </p>
+<p>
+Kittler, R., Mikler, J. (1986): Basis of the utilization of solar radiation 
+(in Slovak). VEDA, Bratislava, p. 150. </p>
+<p>
+Krcho, J. (1990). Morfometrick&aacute; analza a digit&aacute;lne modely georeli&eacute;fu
+(Morphometric analysis and digital models of georelief). VEDA,
+Bratislava (in Slovak). </p>
+<p>
+Muneer, T. (1990). Solar radiation model for Europe. Building services engineering
+research and technology, 11, 4, 153-163. </p>
+<p>
+Neteler, M., Mitasova, H. (2002): Open Source GIS: A GRASS GIS Approach, Kluwer
+Academic Publishers/Springer. </p>
+<p>
+Page, J. ed. (1986). Prediction of solar radiation on inclined surfaces. Solar
+energy R&amp;D in the European Community, series F &#8211; Solar radiation data,
+Dordrecht (D. Reidel), 3, 71, 81-83. </p>
+<p>
+Page, J., Albuisson, M., Wald, L. (2001). The European solar radiation atlas:
+a valuable digital tool. Solar Energy, 71, 81-83. </p>
+<p>
+Rigollier, Ch., Bauer, O., Wald, L. (2000). On the clear sky model of the
+ESRA - European Solar radiation Atlas - with respect to the Heliosat method.
+Solar energy, 68, 33-48. </p>
+<p>
+Scharmer, K., Greif, J., eds., (2000). The European solar radiation atlas,
+Vol. 2: Database and exploitation software. Paris (Les Presses de l&#8217; &Eacute;cole
+des Mines). </p>
+
+<p>Joint Research Centre: <a href="http://re.jrc.ec.europa.eu/pvgis/">GIS solar radiation database for Europe</a> and
+<a href="http://re.jrc.ec.europa.eu/pvgis/solres/solmod3.htm">Solar radiation and GIS</a>
+
+<h2>AUTHORS</h2>
+
+Jaroslav Hofierka, GeoModel, s.r.o. Bratislava, Slovakia <br>
+                                                                        
+Marcel Suri, GeoModel, s.r.o. Bratislava, Slovakia <br>
+
+Thomas Huld, JRC, Italy <br>
+
+&copy; 2007, Jaroslav Hofierka, Marcel Suri 
+<address>
+<a href="MAILTO:hofi at geomodel.sk">hofierka at geomodel.sk</a>
+<a href="MAILTO:suri at geomodel.sk">suri at geomodel.sk</a>
+</address>
+
+<p><i>Last changed: $Date: 2008-07-21 10:07:01 +0200 (Mon, 21 Jul 2008) $</i> </p>

Added: grass/branches/develbranch_6/raster/r.sun2/local_proto.h
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/local_proto.h	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/local_proto.h	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,65 @@
+/* main.c */
+
+
+void where_is_point(double *length, struct SunGeometryVarDay *sunVarGeom,	
+	struct GridGeometry *gridGeom);
+int searching(double *length, struct SunGeometryVarDay *sunVarGeom,	
+	struct GridGeometry *gridGeom);
+
+int useCivilTime();
+void setUseCivilTime(int val);
+int useShadow();
+void setUseShadow(int val);
+int useShadowData();
+void setUseShadowData(int val);
+int useHorizonData();
+void setUseHorizonData(int val);
+double getTimeOffset();
+void setTimeOffset(double val);
+double getHorizonInterval();
+void setHorizonInterval(double val);
+void setAngularLossDenominator();
+
+
+void cube(int, int);
+
+double com_sol_const(int no_of_day);
+
+
+double brad(double, double *bh, struct SunGeometryVarDay *sunVarGeom,
+		struct SunGeometryVarSlope *sunSlopeGeom, 
+		struct SolarRadVar *sunRadVar);
+double drad(double, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
+		struct SunGeometryVarSlope *sunSlopeGeom, 
+		struct SolarRadVar *sunRadVar);
+
+
+double brad_angle_loss(double, double *bh, struct SunGeometryVarDay *sunVarGeom,
+	    struct SunGeometryVarSlope *sunSlopeGeom, 
+	    struct SolarRadVar *sunRadVar);
+double drad_angle_loss(double, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
+	    struct SunGeometryVarSlope *sunSlopeGeom, 
+	    struct SolarRadVar *sunRadVar);
+
+
+void com_par( struct SunGeometryConstDay *sungeom, 
+		struct SunGeometryVarDay *sunVarGeom, 
+		struct GridGeometry *gridGeom,
+		double latitude, double longitude);
+void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
+		struct GridGeometry *gridGeom);
+double lumcline2(struct SunGeometryConstDay *sungeom,
+	struct SunGeometryVarDay *sunVarGeom,
+		struct SunGeometryVarSlope *sunSlopeGeom, 
+		struct GridGeometry *gridGeom,
+		unsigned char *horizonpointer);
+
+
+typedef	double (*BeamRadFunc)(double sh, double *bh, struct SunGeometryVarDay *sunVarGeom,
+	struct SunGeometryVarSlope *sunSlopeGeom, 
+	struct SolarRadVar *sunRadVar);
+	
+typedef	double (*DiffRadFunc)(double sh, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
+	struct SunGeometryVarSlope *sunSlopeGeom, 
+	struct SolarRadVar *sunRadVar);
+


Property changes on: grass/branches/develbranch_6/raster/r.sun2/local_proto.h
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.sun2/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/main.c	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/main.c	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,2136 @@
+/*******************************************************************************
+r.sun: This program was writen by Jaro Hofierka in Summer 1993 and re-engineered
+in 1996-1999. In cooperation with Marcel Suri and Thomas Huld from JRC in Ispra
+a new version of r.sun was prepared using ESRA solar radiation formulas.
+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 Thomas.Huld at jrc.it
+*******************************************************************************/
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the
+ *   Free Software Foundation, Inc.,
+ *   59 Temple Place - Suite 330,
+ *   Boston, MA  02111-1307, USA.
+ */
+
+/*v. 2.0 July 2002, NULL data handling, JH */
+/*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+#include "sunradstruct.h"
+#include "local_proto.h"
+#include "rsunglobals.h"
+
+#define NUM_PARTITIONS "1"
+#define SKIP    "1"
+#define BIG      1.e20
+#define LINKE    "3.0"
+#define SLOPE    "0.0"
+#define ASPECT   "270"
+#define ALB      "0.2"
+#define STEP     "0.5"
+#define BSKY      1.0
+#define DSKY      1.0
+#define DIST     "1.0"
+
+#define SCALING_FACTOR 150.
+const double invScale = 1. / SCALING_FACTOR;
+
+#define AMAX1(arg1, arg2) ((arg1) >= (arg2) ? (arg1) : (arg2))
+#define AMIN1(arg1, arg2) ((arg1) <= (arg2) ? (arg1) : (arg2))
+#define DISTANCE1(x1, x2, y1, y2) (sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2)))
+#define DISTANCE2(x00, y00) ((xx0 - x00)*(xx0 - x00) + (yy0 - y00)*(yy0 - y00))
+
+const double pihalf = M_PI * 0.5;
+const double pi2 = M_PI * 2.;
+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;
+FILE *fw;
+
+
+char *elevin;
+char *aspin;
+char *slopein;
+char *civiltime = NULL;
+char *linkein = NULL;
+char *albedo = NULL;
+char *latin = NULL;
+char *coefbh = NULL;
+char *coefdh = NULL;
+char *incidout = NULL;
+char *longin = NULL;
+char *horizon = NULL;
+char *beam_rad = NULL;
+char *insol_time = NULL;
+char *diff_rad = NULL;
+char *refl_rad = NULL;
+char *glob_rad = NULL;
+char *mapset = NULL;
+char *per;
+char *shade;
+char mapname[1024];
+
+struct Cell_head cellhd;
+struct pj_info iproj;
+struct pj_info oproj;
+struct History hist;
+
+
+int INPUT_part(int offset, double *zmax);
+int OUTGR(void);
+int min(int, int);
+int max(int, int);
+
+void cube(int, int);
+void (*func) (int, int);
+
+void joules2(struct SunGeometryConstDay *sunGeom,
+	     struct SunGeometryVarDay *sunVarGeom,
+	     struct SunGeometryVarSlope *sunSlopeGeom,
+	     struct SolarRadVar *sunRadVar,
+	     struct GridGeometry *gridGeom,
+	     unsigned char *horizonpointer, double latitude, double longitude);
+
+
+void calculate(double singleSlope, double singleAspect,
+	       double singleAlbedo, double singleLinke,
+	       struct GridGeometry gridGeom);
+double com_declin(int);
+
+int n, m, ip, jp;
+int d, day;
+int saveMemory, numPartitions = 1;
+long int shadowoffset = 0;
+int varCount_global = 0;
+int bitCount_global = 0;
+int arrayNumInt = 1;
+float **z = NULL, **o = NULL, **s = NULL, **li = NULL, **a = NULL, **la =
+    NULL, **longitArray, **cbhr = NULL, **cdhr = NULL;
+double op, dp;
+double invstepx, invstepy;
+double sr_min = 24., sr_max = 0., ss_min = 24., ss_max = 0.;
+
+float **lumcl, **beam, **insol, **diff, **refl, **globrad;
+unsigned char *horizonarray = NULL;
+double civilTime;
+/*
+ * double startTime, endTime;
+ */
+double xmin, xmax, ymin, ymax;
+double declin, step, dist;
+double li_max = 0., li_min = 100., al_max = 0., al_min = 1.0, la_max = -90.,
+    la_min = 90.;
+double offsetx = 0.5, offsety = 0.5;
+char *tt, *lt;
+/*
+ * double slope;
+ */
+double o_orig, z1;
+/*
+ * double lum_C11, lum_C13, lum_C22, lum_C31, lum_C33;
+ * double sinSolarAltitude; */
+
+/* Sine of the solar altitude (angle above horizon) 
+ */
+/*
+ * double sunrise_time, sunset_time;
+ * double solarAltitude;
+ * double tanSolarAlt, solarAzimuth;
+ * double stepsinangle, stepcosangle;
+ * double angle;
+ */
+double horizonStep;
+double ltime, tim, timo;
+double declination;		/* Contains the negative of the declination at the chosen day */
+/*
+ * double lum_C31_l, lum_C33_l;
+ */
+double beam_e, diff_e, refl_e, rr, insol_t;
+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;
+    double singleAspect;
+    double singleAlbedo;
+    double singleLinke;
+
+    struct GModule *module;
+    struct
+    {
+	struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
+	    *lin, *albedo, *longin, *alb, *latin, *lat, *coefbh, *coefdh,
+	    *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad, *glob_rad,
+	    *day, *step, *declin, *ltime, *dist, *horizon, *horizonstep,
+	    *numPartitions, *civilTime;
+    }
+    parm;
+
+    struct
+    {
+	struct Flag *shade, *saveMemory;
+    }
+    flag;
+
+    struct GridGeometry gridGeom;
+
+
+    G_gisinit(argv[0]);
+    module = G_define_module();
+
+    module->description =
+	_("Computes direct (beam), diffuse and reflected solar irradiation raster "
+	 "maps for given day, latitude, surface and atmospheric conditions. Solar "
+	 "parameters (e.g. sunrise, sunset times, declination, extraterrestrial "
+	 "irradiance, daylight length) are saved in a local text file. "
+	 "Alternatively, a local time can be specified to compute solar "
+	 "incidence angle and/or irradiance raster maps. The shadowing effect of "
+	 "the topography is optionally incorporated.");
+
+    if (G_get_set_window(&cellhd) == -1)
+	G_fatal_error("G_get_set_window() failed");
+    gridGeom.stepx = cellhd.ew_res;
+    gridGeom.stepy = cellhd.ns_res;
+    invstepx = 1. / gridGeom.stepx;
+    invstepy = 1. / gridGeom.stepy;
+    n /*n_cols */  = cellhd.cols;
+    m /*n_rows */  = cellhd.rows;
+    xmin = cellhd.west;
+    ymin = cellhd.south;
+    xmax = cellhd.east;
+    ymax = cellhd.north;
+    gridGeom.deltx = fabs(cellhd.east - cellhd.west);
+    gridGeom.delty = fabs(cellhd.north - cellhd.south);
+
+    parm.elevin = G_define_option();
+    parm.elevin->key = "elevin";
+    parm.elevin->type = TYPE_STRING;
+    parm.elevin->required = YES;
+    parm.elevin->gisprompt = "old,cell,raster";
+    parm.elevin->description = _("Name of the input elevation raster map [meters]");
+    parm.elevin->guisection = _("Input_options");
+
+    parm.aspin = G_define_option();
+    parm.aspin->key = "aspin";
+    parm.aspin->type = TYPE_STRING;
+    parm.aspin->required = NO;
+    parm.aspin->gisprompt = "old,cell,raster";
+    parm.aspin->description = _("Name of the input aspect map (terrain aspect or azimuth of the solar panel) [decimal degrees]");
+    parm.aspin->guisection = _("Input_options");
+
+    parm.aspect = G_define_option();
+    parm.aspect->key = "aspect";
+    parm.aspect->type = TYPE_DOUBLE;
+    parm.aspect->answer = ASPECT;
+    parm.aspect->required = NO;
+    parm.aspect->description =
+	_("A single value of the orientation (aspect), 270 is south");
+    parm.aspect->guisection = _("Input_options");
+
+    parm.slopein = G_define_option();
+    parm.slopein->key = "slopein";
+    parm.slopein->type = TYPE_STRING;
+    parm.slopein->required = NO;
+    parm.slopein->gisprompt = "old,cell,raster";
+    parm.slopein->description = _("Name of the input slope raster map (terrain slope or solar panel inclination) [decimal degrees]");
+    parm.slopein->guisection = _("Input_options");
+
+    parm.slope = G_define_option();
+    parm.slope->key = "slope";
+    parm.slope->type = TYPE_DOUBLE;
+    parm.slope->answer = SLOPE;
+    parm.slope->required = NO;
+    parm.slope->description = _("A single value of inclination (slope)");
+    parm.slope->guisection = _("Input_options");
+
+    parm.linkein = G_define_option();
+    parm.linkein->key = "linkein";
+    parm.linkein->type = TYPE_STRING;
+    parm.linkein->required = NO;
+    parm.linkein->gisprompt = "old,cell,raster";
+    parm.linkein->description =
+	_("Name of the Linke atmospheric turbidity coefficient input raster map [-]");
+    parm.linkein->guisection = _("Input_options");
+
+    if (parm.linkein->answer == NULL) {
+	parm.lin = G_define_option();
+	parm.lin->key = "lin";
+	parm.lin->type = TYPE_DOUBLE;
+	parm.lin->answer = LINKE;
+	parm.lin->required = NO;
+	parm.lin->description =
+	    _("A single value of the Linke atmospheric turbidity coefficient [-]");
+	parm.lin->guisection = _("Input_options");
+    }
+
+    parm.albedo = G_define_option();
+    parm.albedo->key = "albedo";
+    parm.albedo->type = TYPE_STRING;
+    parm.albedo->required = NO;
+    parm.albedo->gisprompt = "old,cell,raster";
+    parm.albedo->description = _("Name of the ground albedo coefficient input raster map [-]");
+    parm.albedo->guisection = _("Input_options");
+
+    if (parm.albedo->answer == NULL) {
+	parm.alb = G_define_option();
+	parm.alb->key = "alb";
+	parm.alb->type = TYPE_DOUBLE;
+	parm.alb->answer = ALB;
+	parm.alb->required = NO;
+	parm.alb->description = _("A single value of the ground albedo coefficient [-]");
+	parm.alb->guisection = _("Input_options");
+    }
+
+    parm.latin = G_define_option();
+    parm.latin->key = "latin";
+    parm.latin->type = TYPE_STRING;
+    parm.latin->required = NO;
+    parm.latin->gisprompt = "old,cell,raster";
+    parm.latin->description = _("Name of the latitudes input raster map [decimal degrees]");
+    parm.latin->guisection = _("Input_options");
+
+    if (parm.latin->answer == NULL) {
+	parm.lat = G_define_option();
+	parm.lat->key = "lat";
+	parm.lat->type = TYPE_DOUBLE;
+	parm.lat->required = NO;
+	parm.lat->description = _("A single value of latitude [decimal degrees]");
+	parm.lat->guisection = _("Input_options");
+    }
+
+    parm.longin = G_define_option();
+    parm.longin->key = "longin";
+    parm.longin->type = TYPE_STRING;
+    parm.longin->required = NO;
+    parm.longin->gisprompt = "old,cell,raster";
+    parm.longin->description = _("Name of the longitude input raster map [decimal degrees]");
+    parm.longin->guisection = _("Input_options");
+
+    parm.coefbh = G_define_option();
+    parm.coefbh->key = "coefbh";
+    parm.coefbh->type = TYPE_STRING;
+    parm.coefbh->required = NO;
+    parm.coefbh->gisprompt = "old,cell,raster";
+    parm.coefbh->description =
+	_("Name of real-sky beam radiation coefficient raster map [-]");
+    parm.coefbh->guisection = _("Input_options");
+
+    parm.coefdh = G_define_option();
+    parm.coefdh->key = "coefdh";
+    parm.coefdh->type = TYPE_STRING;
+    parm.coefdh->required = NO;
+    parm.coefdh->gisprompt = "old,cell,raster";
+    parm.coefdh->description =
+	_("Name of real-sky diffuse radiation coefficient raster map [-]");
+    parm.coefdh->guisection = _("Input_options");
+
+    parm.horizon = G_define_option();
+    parm.horizon->key = "horizon";
+    parm.horizon->type = TYPE_STRING;
+    parm.horizon->required = NO;
+    parm.horizon->gisprompt = "old,cell,raster";
+    parm.horizon->description = _("The horizon information file prefix");
+    parm.horizon->guisection = _("Input_options");
+
+    parm.horizonstep = G_define_option();
+    parm.horizonstep->key = "horizonstep";
+    parm.horizonstep->type = TYPE_DOUBLE;
+    parm.horizonstep->required = NO;
+    parm.horizonstep->description =
+	_("Angle step size for the horizon information [degrees]");
+    parm.horizonstep->guisection = _("Input_options");
+
+    parm.incidout = G_define_option();
+    parm.incidout->key = "incidout";
+    parm.incidout->type = TYPE_STRING;
+    parm.incidout->required = NO;
+    parm.incidout->description = _("Output incidence angle raster map (mode 1 only)");
+    parm.incidout->guisection = _("Output_options");
+
+    parm.beam_rad = G_define_option();
+    parm.beam_rad->key = "beam_rad";
+    parm.beam_rad->type = TYPE_STRING;
+    parm.beam_rad->required = NO;
+    parm.beam_rad->description =
+	_("Output beam irradiance [W.m-2] (mode 1) or irradiation raster map [Wh.m-2.day-1] (mode 2)");
+    parm.beam_rad->guisection = _("Output_options");
+
+    parm.insol_time = G_define_option();
+    parm.insol_time->key = "insol_time";
+    parm.insol_time->type = TYPE_STRING;
+    parm.insol_time->required = NO;
+    parm.insol_time->description = _("Output insolation time raster map [h] (mode 2 only)");
+    parm.insol_time->guisection = _("Output_options");
+
+    parm.diff_rad = G_define_option();
+    parm.diff_rad->key = "diff_rad";
+    parm.diff_rad->type = TYPE_STRING;
+    parm.diff_rad->required = NO;
+    parm.diff_rad->description =
+	_("Output diffuse irradiance [W.m-2] (mode 1) or irradiation raster map [Wh.m-2.day-1] (mode 2)");
+    parm.diff_rad->guisection = _("Output_options");
+
+    parm.refl_rad = G_define_option();
+    parm.refl_rad->key = "refl_rad";
+    parm.refl_rad->type = TYPE_STRING;
+    parm.refl_rad->required = NO;
+    parm.refl_rad->description =
+	_("Output ground reflected irradiance [W.m-2] (mode 1) or irradiation raster map [Wh.m-2.day-1] (mode 2)");
+    parm.refl_rad->guisection = _("Output_options");
+
+    parm.glob_rad = G_define_option();
+    parm.glob_rad->key = "glob_rad";
+    parm.glob_rad->type = TYPE_STRING;
+    parm.glob_rad->required = NO;
+    parm.glob_rad->description =
+	_("Output global (total) irradiance/irradiation [W.m-2] (mode 1) or irradiance/irradiation raster map [Wh.m-2.day-1] (mode 2)");
+    parm.glob_rad->guisection = _("Output_options");
+
+
+    parm.day = G_define_option();
+    parm.day->key = "day";
+    parm.day->type = TYPE_INTEGER;
+    parm.day->required = YES;
+    parm.day->description = _("No. of day of the year (1-365)");
+
+    parm.step = G_define_option();
+    parm.step->key = "step";
+    parm.step->type = TYPE_DOUBLE;
+    parm.step->answer = STEP;
+    parm.step->required = NO;
+    parm.step->description = _("Time step when computing all-day radiation sums [decimal hours]");
+
+    parm.declin = G_define_option();
+    parm.declin->key = "declin";
+    parm.declin->type = TYPE_DOUBLE;
+    parm.declin->required = NO;
+    parm.declin->description =
+	_("Declination value (overriding the internally computed value) [radians]");
+
+    parm.ltime = G_define_option();
+    parm.ltime->key = "time";
+    parm.ltime->type = TYPE_DOUBLE;
+    /*          parm.ltime->answer = TIME; */
+    parm.ltime->required = NO;
+    parm.ltime->description = _("Local (solar) time (to be set for mode 1 only) [decimal hours]");
+
+    /*
+     * parm.startTime = G_define_option();
+     * parm.startTime->key = "starttime";
+     * parm.startTime->type = TYPE_DOUBLE;
+     * parm.startTime->required = NO;
+     * parm.startTime->description = _("Starting time for calculating results for several different times.");
+     * 
+     * parm.endTime = G_define_option();
+     * parm.endTime->key = "endtime";
+     * parm.endTime->type = TYPE_DOUBLE;
+     * parm.endTime->required = NO;
+     * parm.endTime->description = _("End time for calculating results for several different times.)";
+     */
+
+    parm.dist = G_define_option();
+    parm.dist->key = "dist";
+    parm.dist->type = TYPE_DOUBLE;
+    parm.dist->answer = DIST;
+    parm.dist->required = NO;
+    parm.dist->description = _("Sampling distance step coefficient (0.5-1.5)");
+
+    parm.numPartitions = G_define_option();
+    parm.numPartitions->key = "numpartitions";
+    parm.numPartitions->type = TYPE_INTEGER;
+    parm.numPartitions->answer = NUM_PARTITIONS;
+    parm.numPartitions->required = NO;
+    parm.numPartitions->description =
+	_("Read the input files in this number of chunks");
+
+    parm.civilTime = G_define_option();
+    parm.civilTime->key = "civiltime";
+    parm.civilTime->type = TYPE_DOUBLE;
+    parm.civilTime->required = NO;
+    parm.civilTime->description =
+	_("Civil time zone value, if none, the time will be local solar time");
+
+
+    flag.shade = G_define_flag();
+    flag.shade->key = 's';
+    flag.shade->description =
+	_("Incorporate the shadowing effect of terrain");
+
+    flag.saveMemory = G_define_flag();
+    flag.saveMemory->key = 'm';
+    flag.saveMemory->description =
+	_("Use the low-memory version of the program");
+
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+
+    setUseShadow(flag.shade->answer);
+
+    /*
+     * if(shd)
+     * {
+     * 
+     * }
+     */
+    saveMemory = flag.saveMemory->answer;
+    civiltime = parm.civilTime->answer;
+
+
+    elevin = parm.elevin->answer;
+    aspin = parm.aspin->answer;
+    slopein = parm.slopein->answer;
+    linkein = parm.linkein->answer;
+    albedo = parm.albedo->answer;
+    latin = parm.latin->answer;
+
+
+    if (civiltime != NULL) {
+	setUseCivilTime(1);
+	longin = parm.longin->answer;
+
+	if (longin == NULL) {
+	    G_fatal_error(_("You must give the longitude raster if you use civil time"));
+
+	}
+	sscanf(parm.civilTime->answer, "%lf", &civilTime);
+
+	/* Normalize if somebody should be weird enough to give more than +- 12 
+	 * hours offset. */
+
+	civilTime -= 24 * ((int)(civilTime / 24.));
+	if (civilTime < -12.) {
+	    civilTime += 24.;
+	}
+	else if (civilTime > 12.) {
+	    civilTime -= 24;
+	}
+
+    }
+    else {
+	setUseCivilTime(0);
+    }
+    coefbh = parm.coefbh->answer;
+    coefdh = parm.coefdh->answer;
+    incidout = parm.incidout->answer;
+    horizon = parm.horizon->answer;
+    setUseHorizonData(horizon != NULL);
+    beam_rad = parm.beam_rad->answer;
+    insol_time = parm.insol_time->answer;
+    diff_rad = parm.diff_rad->answer;
+    refl_rad = parm.refl_rad->answer;
+    glob_rad = parm.glob_rad->answer;
+
+    if ((insol_time != NULL) && (incidout != NULL))
+	G_fatal_error(_("insol_time and incidout are incompatible options"));
+
+    sscanf(parm.day->answer, "%d", &day);
+    sscanf(parm.step->answer, "%lf", &step);
+
+    if (parm.step->answer != NULL) {
+	if (sscanf(parm.step->answer, "%lf", &step) != 1)
+	    G_fatal_error(_("Error reading time step size"));
+    }
+    if (parm.horizonstep->answer != NULL) {
+	if (sscanf(parm.horizonstep->answer, "%lf", &horizonStep) != 1)
+	    G_fatal_error(_("Error reading horizon step size"));
+	if(horizonStep > 0.)
+	    setHorizonInterval(deg2rad * horizonStep);
+	else
+	    G_fatal_error(_("The horizon step size must be greater than 0."));
+
+    }
+	else if(useHorizonData()) {
+		G_fatal_error(_("If you use the horizon option you must also set the 'horizonstep' parameter."));
+	     }
+
+
+    tt = parm.ltime->answer;
+    if (parm.ltime->answer != NULL) {
+	if (insol_time != NULL)
+	    G_fatal_error(_("Time and insol_time are incompatible options"));
+	G_message(_("Mode 1: instantaneous solar incidence angle & irradiance using a set local time"));
+	sscanf(parm.ltime->answer, "%lf", &timo);
+    }
+    else {
+	if (incidout != NULL)
+	    G_fatal_error(_("incidout requres time parameter to be set"));
+	G_message(_("Mode 2: integrated daily irradiation"));
+    }
+
+    /*      
+     * if (parm.startTime->answer != NULL) sscanf(parm.startTime->answer, "%lf", &startTime);
+     * if (parm.endTime->answer != NULL) sscanf(parm.endTime->answer, "%lf", &endTime);
+     * 
+     * if((parm.startTime->answer != NULL) ||(parm.endTime->answer != NULL))
+     * {
+     */
+    /* The start and end times should only be defined if the single
+     * time is not defined, and if the step size *is* defined. */
+    /*      
+     * if(parm.step->answer==NULL)
+     *    G_fatal_error("If you want to use a time interval you must also define a step size.");
+     * if(parm.ltime->answer!=NULL)
+     *    G_fatal_error("If you want to use a time interval you can't define a single time value.\n");
+     * if((parm.startTime->answer==NULL)||(parm.endTime->answer==NULL))
+     *    G_fatal_error("If you want to use a time interval both the start and end times must be defined.\n");
+     * }
+     */
+    if (parm.linkein->answer == NULL)
+	sscanf(parm.lin->answer, "%lf", &singleLinke);
+    if (parm.albedo->answer == NULL)
+	sscanf(parm.alb->answer, "%lf", &singleAlbedo);
+    lt = parm.lat->answer;
+    /*
+     * if (parm.lat->answer != NULL)
+     * sscanf(parm.lat->answer, "%lf", &latitude);
+     */
+/* HB 6/2008: why is the above commented out? instead of sscanf, maybe nicer to use:
+    G_scan_lat(parm.lat->answer, &latitude);
+*/
+    if (parm.slopein->answer == NULL)
+	sscanf(parm.slope->answer, "%lf", &singleSlope);
+    singleSlope *= deg2rad;
+
+    if (parm.aspin->answer == NULL)
+	sscanf(parm.aspect->answer, "%lf", &singleAspect);
+    singleAspect *= deg2rad;
+
+    if (parm.coefbh->answer == NULL)
+	cbh = BSKY;
+    if (parm.coefdh->answer == NULL)
+	cdh = DSKY;
+    sscanf(parm.dist->answer, "%lf", &dist);
+
+    if (parm.numPartitions->answer != NULL) {
+	sscanf(parm.numPartitions->answer, "%d", &numPartitions);
+	if (useShadow() && (!useHorizonData()) && (numPartitions != 1)) {
+	    /* If you calculate shadows on the fly, the number of partitions
+	     * must be one.
+	     */
+	    G_fatal_error(_("If you use -s and no horizon rasters, numpartitions must be =1"));
+
+	}
+
+    }
+
+    gridGeom.stepxy = dist * 0.5 * (gridGeom.stepx + gridGeom.stepy);
+    TOLER = gridGeom.stepxy * EPS;
+
+    /* The save memory scheme will not work if you want to calculate shadows
+     * on the fly. If you calculate without shadow effects or if you have the
+     * shadows pre-calculated, there is no problem. */
+
+    if (saveMemory && useShadow() && (!useHorizonData()))
+	G_fatal_error(_("If you want to save memory and to use shadows, you must use pre-calculated horizons."));
+
+    if (parm.declin->answer == NULL)
+	declination = com_declin(day);
+    else {
+	sscanf(parm.declin->answer, "%lf", &declin);
+	declination = -declin;
+    }
+
+
+    /*
+     * if (lt != NULL)
+     * latitude = -latitude * deg2rad;
+     */
+    if (tt != 0) {
+	/* Shadow for just one time during the day */
+	if (horizon == NULL) {
+	    arrayNumInt = 1;
+	}
+	else if (useHorizonData()) {
+	    arrayNumInt = (int)(360. / horizonStep);
+
+	}
+    }
+    else {
+	if (useHorizonData()) {
+	    /*        Number of bytes holding the horizon information */
+	    arrayNumInt = (int)(360. / horizonStep);
+	}
+    }
+
+    if (tt != NULL) {
+
+	tim = (timo - 12) * 15;
+	/* converting to degrees */
+	/* Jenco (12-timeAngle) * 15 */
+	if (tim < 0)
+	    tim += 360;
+	tim = M_PI * tim / 180;
+	/* conv. to radians */
+    }
+
+    /* Set up parameters for projection to lat/long if necessary */
+
+
+    struct Key_Value *in_proj_info, *in_unit_info;
+
+    if ((in_proj_info = G_get_projinfo()) == NULL)
+	G_fatal_error
+	    (_("Can't get projection info of current location: please set latitude via 'lat' or 'latin' option!"));
+
+    if ((in_unit_info = G_get_projunits()) == NULL)
+	G_fatal_error(_("Can't get projection units of current location"));
+
+    if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
+	G_fatal_error(_("Can't get projection key values of current location"));
+
+    G_free_key_value(in_proj_info);
+    G_free_key_value(in_unit_info);
+
+    /* Set output projection to latlong w/ same ellipsoid */
+    oproj.zone = 0;
+    oproj.meters = 1.;
+    sprintf(oproj.proj, "ll");
+    if ((oproj.pj = pj_latlong_from_proj(iproj.pj)) == NULL)
+	G_fatal_error(_("Unable to set up lat/long projection parameters"));
+
+
+/**********end of parser - ******************************/
+
+    if ((G_projection() == PROJECTION_LL)) {
+	ll_correction = 1;
+    }
+
+
+    calculate(singleSlope, singleAspect, singleAlbedo, singleLinke, gridGeom);
+    OUTGR();
+
+    return 1;
+}
+
+
+int INPUT_part(int offset, double *zmax)
+{
+    int finalRow, rowrevoffset;
+    int numRows;
+    int numDigits;
+    FCELL *cell1 = NULL, *cell2 = NULL;
+    FCELL *cell3 = NULL, *cell4 = NULL, *cell5 = NULL, *cell6 = NULL, *cell7 =
+	NULL;
+    FCELL *rast1 = NULL, *rast2 = NULL;
+    static FCELL **horizonbuf;
+    unsigned char *horizonpointer;
+    int fd1 = -1, fd2 = -1, fd3 = -1, fd4 = -1, fd5 = -1, fd6 = -1, fd7 =
+	-1, row, row_rev;
+    static int *fd_shad;
+    int fr1 = -1, fr2 = -1;
+    int l, i, j;
+    char shad_filename[256];
+    char formatString[10];
+
+    finalRow = m - offset - m / numPartitions;
+    if (finalRow < 0) {
+	finalRow = 0;
+    }
+
+    numRows = m / numPartitions;
+
+    cell1 = G_allocate_f_raster_buf();
+
+    if (z == NULL) {
+	z = (float **)G_malloc(sizeof(float *) * (numRows));
+
+
+	for (l = 0; l < numRows; l++) {
+	    z[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+    }
+
+
+    if ((mapset = G_find_cell(elevin, "")) == NULL)
+	G_fatal_error("Elevation raster file not found");
+
+
+
+    fd1 = G_open_cell_old(elevin, mapset);
+
+    if (slopein != NULL) {
+	cell3 = G_allocate_f_raster_buf();
+	if (s == NULL) {
+	    s = (float **)G_malloc(sizeof(float *) * (numRows));
+
+	    for (l = 0; l < numRows; l++) {
+		s[l] = (float *)G_malloc(sizeof(float) * (n));
+	    }
+
+	}
+	if ((mapset = G_find_cell(slopein, "")) == NULL)
+	       G_fatal_error(_("Raster map <%s> not found"),slopein);
+	fd3 = G_open_cell_old(slopein, mapset);
+
+    }
+
+    if (aspin != NULL) {
+	cell2 = G_allocate_f_raster_buf();
+
+	if (o == NULL) {
+	    o = (float **)G_malloc(sizeof(float *) * (numRows));
+
+	    for (l = 0; l < numRows; l++) {
+		o[l] = (float *)G_malloc(sizeof(float) * (n));
+	    }
+	}
+
+	if ((mapset = G_find_cell(aspin, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),aspin);
+	fd2 = G_open_cell_old(aspin, mapset);
+
+    }
+
+
+    if (linkein != NULL) {
+	cell4 = G_allocate_f_raster_buf();
+	if (li == NULL) {
+	    li = (float **)G_malloc(sizeof(float *) * (numRows));
+	    for (l = 0; l < numRows; l++)
+		li[l] = (float *)G_malloc(sizeof(float) * (n));
+
+	}
+	if ((mapset = G_find_cell(linkein, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),linkein);
+
+	fd4 = G_open_cell_old(linkein, mapset);
+    }
+
+    if (albedo != NULL) {
+	cell5 = G_allocate_f_raster_buf();
+	if (a == NULL) {
+	    a = (float **)G_malloc(sizeof(float *) * (numRows));
+	    for (l = 0; l < numRows; l++)
+		a[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+	if ((mapset = G_find_cell(albedo, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),albedo);
+
+	fd5 = G_open_cell_old(albedo, mapset);
+    }
+
+    if (latin != NULL) {
+	cell6 = G_allocate_f_raster_buf();
+	if (la == NULL) {
+	    la = (float **)G_malloc(sizeof(float *) * (numRows));
+	    for (l = 0; l < numRows; l++)
+		la[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+	if ((mapset = G_find_cell(latin, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),latin);
+
+	fd6 = G_open_cell_old(latin, mapset);
+    }
+
+    if (longin != NULL) {
+	cell7 = G_allocate_f_raster_buf();
+	longitArray = (float **)G_malloc(sizeof(float *) * (numRows));
+	for (l = 0; l < numRows; l++)
+	    longitArray[l] = (float *)G_malloc(sizeof(float) * (n));
+
+	if ((mapset = G_find_cell(longin, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),longin);
+
+	fd7 = G_open_cell_old(longin, mapset);
+    }
+
+
+
+    if (coefbh != NULL) {
+	rast1 = G_allocate_f_raster_buf();
+
+	if (cbhr == NULL) {
+	    cbhr = (float **)G_malloc(sizeof(float *) * (numRows));
+	    for (l = 0; l < numRows; l++)
+		cbhr[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+	if ((mapset = G_find_cell(coefbh, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),coefbh);
+
+	fr1 = G_open_cell_old(coefbh, mapset);
+    }
+
+    if (coefdh != NULL) {
+	rast2 = G_allocate_f_raster_buf();
+	if (cdhr == NULL) {
+	    cdhr = (float **)G_malloc(sizeof(float *) * (numRows));
+	    for (l = 0; l < numRows; l++)
+		cdhr[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+	if ((mapset = G_find_cell(coefdh, "")) == NULL)
+	    G_fatal_error(_("Raster map <%s> not found"),coefdh);
+
+	fr2 = G_open_cell_old(coefdh, mapset);
+    }
+
+
+
+
+    if (useHorizonData()) {
+	if (horizonarray == NULL) {
+	    horizonarray =
+		(unsigned char *)G_malloc(sizeof(char) * arrayNumInt * numRows *
+					n);
+
+	    horizonbuf = (FCELL **) G_malloc(sizeof(FCELL *) * arrayNumInt);
+	    fd_shad = (int *)G_malloc(sizeof(int) * arrayNumInt);
+	}
+	/*
+	 * if(tt != NULL)
+	 * {
+	 * 
+	 * horizonbuf[0]=G_allocate_f_raster_buf();
+	 * sprintf(shad_filename, "%s_%02d", horizon, arrayNumInt);
+	 * if((mapset=G_find_cell(shad_filename,""))==NULL)
+	 * G_message("Horizon file no. %d not found\n", arrayNumInt);
+	 * 
+	 * fd_shad[0] = G_open_cell_old(shad_filename,mapset);
+	 * }
+	 * else
+	 * {
+	 */
+	numDigits = (int)(log10(1. * arrayNumInt)) + 1;
+	sprintf(formatString, "%%s_%%0%dd", numDigits);
+	for (i = 0; i < arrayNumInt; i++) {
+		horizonbuf[i] = G_allocate_f_raster_buf();
+	    sprintf(shad_filename, formatString, horizon, i);
+	    if ((mapset = G_find_cell(shad_filename, "")) == NULL)
+		G_fatal_error(_("Horizon file no. %d <%s> not found"), i, shad_filename);
+
+	    fd_shad[i] = G_open_cell_old(shad_filename, mapset);
+	}
+    }
+    /*
+     * }
+     */
+
+    if (useHorizonData()) {
+
+
+	for (i = 0; i < arrayNumInt; i++) {
+	    for (row = m - offset - 1; row >= finalRow; row--) {
+
+		row_rev = m - row - 1;
+		rowrevoffset = row_rev - offset;
+		G_get_f_raster_row(fd_shad[i], horizonbuf[i], row);
+		horizonpointer = horizonarray + (ssize_t) arrayNumInt * n * rowrevoffset;
+		for (j = 0; j < n; j++) {
+
+		    horizonpointer[i] = (char)(rint(SCALING_FACTOR *
+						    AMIN1(horizonbuf[i][j],
+							  256 * invScale)));
+		    horizonpointer += arrayNumInt;
+
+		}
+	    }
+	}
+
+
+    }
+
+
+
+    for (row = m - offset - 1; row >= finalRow; row--) {
+	G_get_f_raster_row(fd1, cell1, row);
+	if (aspin != NULL)
+	    G_get_f_raster_row(fd2, cell2, row);
+	if (slopein != NULL)
+	    G_get_f_raster_row(fd3, cell3, row);
+	if (linkein != NULL)
+	    G_get_f_raster_row(fd4, cell4, row);
+	if (albedo != NULL)
+	    G_get_f_raster_row(fd5, cell5, row);
+	if (latin != NULL)
+	    G_get_f_raster_row(fd6, cell6, row);
+	if (longin != NULL)
+	    G_get_f_raster_row(fd7, cell7, row);
+	if (coefbh != NULL)
+	    G_get_f_raster_row(fr1, rast1, row);
+	if (coefdh != NULL)
+	    G_get_f_raster_row(fr2, rast2, row);
+
+
+
+	row_rev = m - row - 1;
+	rowrevoffset = row_rev - offset;
+
+	for (j = 0; j < n; j++) {
+	    if (!G_is_f_null_value(cell1 + j))
+		z[rowrevoffset][j] = (float)cell1[j];
+	    else
+		z[rowrevoffset][j] = UNDEFZ;
+
+	    if (aspin != NULL) {
+		if (!G_is_f_null_value(cell2 + j))
+		    o[rowrevoffset][j] = (float)cell2[j];
+		else
+		    o[rowrevoffset][j] = UNDEFZ;
+	    }
+	    if (slopein != NULL) {
+		if (!G_is_f_null_value(cell3 + j))
+		    s[rowrevoffset][j] = (float)cell3[j];
+		else
+		    s[rowrevoffset][j] = UNDEFZ;
+	    }
+
+	    if (linkein != NULL) {
+		if (!G_is_f_null_value(cell4 + j))
+		    li[rowrevoffset][j] = (float)cell4[j];
+		else
+		    li[rowrevoffset][j] = UNDEFZ;
+	    }
+
+	    if (albedo != NULL) {
+		if (!G_is_f_null_value(cell5 + j))
+		    a[rowrevoffset][j] = (float)cell5[j];
+		else
+		    a[rowrevoffset][j] = UNDEFZ;
+	    }
+
+	    if (latin != NULL) {
+		if (!G_is_f_null_value(cell6 + j))
+		    la[rowrevoffset][j] = (float)cell6[j];
+		else
+		    la[rowrevoffset][j] = UNDEFZ;
+	    }
+
+	    if (longin != NULL) {
+		if (!G_is_f_null_value(cell7 + j))
+		    longitArray[rowrevoffset][j] = (float)cell7[j];
+		else
+		    longitArray[rowrevoffset][j] = UNDEFZ;
+	    }
+
+	    if (coefbh != NULL) {
+		if (!G_is_f_null_value(rast1 + j))
+		    cbhr[rowrevoffset][j] = (float)rast1[j];
+		else
+		    cbhr[rowrevoffset][j] = UNDEFZ;
+	    }
+
+	    if (coefdh != NULL) {
+		if (!G_is_f_null_value(rast2 + j))
+		    cdhr[rowrevoffset][j] = (float)rast2[j];
+		else
+		    cdhr[rowrevoffset][j] = UNDEFZ;
+	    }
+
+
+
+
+
+	}
+    }
+
+    G_close_cell(fd1);
+    G_free(cell1);
+
+    if (aspin != NULL) {
+	G_free(cell2);
+	G_close_cell(fd2);
+    }
+    if (slopein != NULL) {
+	G_free(cell3);
+	G_close_cell(fd3);
+    }
+    if (linkein != NULL) {
+	G_free(cell4);
+	G_close_cell(fd4);
+    }
+    if (albedo != NULL) {
+	G_free(cell5);
+	G_close_cell(fd5);
+    }
+    if (latin != NULL) {
+	G_free(cell6);
+	G_close_cell(fd6);
+    }
+    if (longin != NULL) {
+	G_free(cell7);
+	G_close_cell(fd7);
+    }
+    if (coefbh != NULL) {
+	G_free(rast1);
+	G_close_cell(fr1);
+    }
+    if (coefdh != NULL) {
+	G_free(rast2);
+	G_close_cell(fr2);
+    }
+
+
+    if (useHorizonData()) {
+	for (i = 0; i < arrayNumInt; i++) {
+	    G_close_cell(fd_shad[i]);
+	    G_free(horizonbuf[i]);
+	}
+    }
+
+
+/*******transformation of angles from 0 to east counterclock
+        to 0 to north clocwise, for ori=0 upslope flowlines
+        turn the orientation 2*M_PI ************/
+
+    /* needs to be eliminated */
+
+    /*for (i = 0; i < m; ++i) */
+    for (i = 0; i < numRows; i++) {
+	for (j = 0; j < n; j++) {
+	    *zmax = AMAX1(*zmax, z[i][j]);
+	    if (aspin != NULL) {
+		if (o[i][j] != 0.) {
+		    if (o[i][j] < 90.)
+			o[i][j] = 90. - o[i][j];
+		    else
+			o[i][j] = 450. - o[i][j];
+		}
+		/*   G_debug(3,"o,z = %d  %d i,j, %d %d \n", o[i][j],z[i][j],i,j); */
+
+		if ((aspin != NULL) && o[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+		if ((slopein != NULL) && s[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+		if (linkein != NULL && li[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+		if (albedo != NULL && a[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+		if (latin != NULL && la[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+		if (coefbh != NULL && cbhr[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+		if (coefdh != NULL && cdhr[i][j] == UNDEFZ)
+		    z[i][j] = UNDEFZ;
+	    }
+
+	}
+    }
+
+    return 1;
+}
+
+int OUTGR(void)
+{
+    FCELL *cell7 = NULL, *cell8 = NULL, *cell9 = NULL, *cell10 = NULL, *cell11 =
+	NULL, *cell12 = NULL;
+    int fd7 = -1, fd8 = -1, fd9 = -1, fd10 = -1, fd11 = -1, fd12 = -1;
+    int i, iarc, j;
+
+    if (incidout != NULL) {
+	cell7 = G_allocate_f_raster_buf();
+	fd7 = G_open_fp_cell_new(incidout);
+	if (fd7 < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), incidout);
+    }
+
+    if (beam_rad != NULL) {
+	cell8 = G_allocate_f_raster_buf();
+	fd8 = G_open_fp_cell_new(beam_rad);
+	if (fd8 < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), beam_rad);
+    }
+
+    if (insol_time != NULL) {
+	cell11 = G_allocate_f_raster_buf();
+	fd11 = G_open_fp_cell_new(insol_time);
+	if (fd11 < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), insol_time);
+    }
+
+    if (diff_rad != NULL) {
+	cell9 = G_allocate_f_raster_buf();
+	fd9 = G_open_fp_cell_new(diff_rad);
+	if (fd9 < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), diff_rad);
+    }
+
+    if (refl_rad != NULL) {
+	cell10 = G_allocate_f_raster_buf();
+	fd10 = G_open_fp_cell_new(refl_rad);
+	if (fd10 < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), refl_rad);
+    }
+
+    if (glob_rad != NULL) {
+	cell12 = G_allocate_f_raster_buf();
+	fd12 = G_open_fp_cell_new(glob_rad);
+	if (fd12 < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), glob_rad);
+    }
+
+
+    if (G_set_window(&cellhd) < 0)
+	G_fatal_error("Cannot set region to output region!");
+
+    if (m != G_window_rows())
+	G_fatal_error("OOPS: rows changed from %d to %d", m, G_window_rows());
+    if (n != G_window_cols())
+	G_fatal_error("OOPS: cols changed from %d to %d", n, G_window_cols());
+
+    for (iarc = 0; iarc < m; iarc++) {
+	i = m - iarc - 1;
+	if (incidout != NULL) {
+	    for (j = 0; j < n; j++) {
+		if (lumcl[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell7 + j, 1);
+		else
+		    cell7[j] = (FCELL) lumcl[i][j];
+	    }
+	    G_put_f_raster_row(fd7, cell7);
+	}
+
+	if (beam_rad != NULL) {
+	    for (j = 0; j < n; j++) {
+		if (beam[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell8 + j, 1);
+		else
+		    cell8[j] = (FCELL) beam[i][j];
+
+	    }
+	    G_put_f_raster_row(fd8, cell8);
+	}
+
+	if (glob_rad != NULL) {
+	    for (j = 0; j < n; j++) {
+		if (globrad[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell12 + j, 1);
+		else
+		    cell12[j] = (FCELL) globrad[i][j];
+
+	    }
+	    G_put_f_raster_row(fd12, cell12);
+	}
+
+
+	if (insol_time != NULL) {
+	    for (j = 0; j < n; j++) {
+		if (insol[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell11 + j, 1);
+		else
+		    cell11[j] = (FCELL) insol[i][j];
+	    }
+	    G_put_f_raster_row(fd11, cell11);
+	}
+
+
+	if (diff_rad != NULL) {
+	    for (j = 0; j < n; j++) {
+		if (diff[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell9 + j, 1);
+		else
+		    cell9[j] = (FCELL) diff[i][j];
+	    }
+	    G_put_f_raster_row(fd9, cell9);
+	}
+
+	if (refl_rad != NULL) {
+	    for (j = 0; j < n; j++) {
+		if (refl[i][j] == UNDEFZ)
+		    G_set_f_null_value(cell10 + j, 1);
+		else
+		    cell10[j] = (FCELL) refl[i][j];
+	    }
+	    G_put_f_raster_row(fd10, cell10);
+	}
+
+    }
+
+    if (incidout != NULL) {
+	G_close_cell(fd7);
+	G_write_history(incidout, &hist);
+    }
+    if (beam_rad != NULL) {
+	G_close_cell(fd8);
+	G_write_history(beam_rad, &hist);
+    }
+    if (diff_rad != NULL) {
+	G_close_cell(fd9);
+	G_write_history(diff_rad, &hist);
+    }
+    if (refl_rad != NULL) {
+	G_close_cell(fd10);
+	G_write_history(refl_rad, &hist);
+    }
+    if (insol_time != NULL) {
+	G_close_cell(fd11);
+	G_write_history(insol_time, &hist);
+    }
+    if (glob_rad != NULL) {
+	G_close_cell(fd12);
+	G_write_history(glob_rad, &hist);
+    }
+
+    return 1;
+}
+
+/*  min(), max() are unused
+ * int min(arg1, arg2)
+ * int arg1;
+ * int arg2;
+ * {
+ * int res;
+ * if (arg1 <= arg2) {
+ * res = arg1;
+ * }
+ * else {
+ * res = arg2;
+ * }
+ * return res;
+ * }
+ * 
+ * int max(arg1, arg2)
+ * int arg1;
+ * int arg2;
+ * {
+ * int res;
+ * if (arg1 >= arg2) {
+ * res = arg1;
+ * }
+ * else {
+ * res = arg2;
+ * }
+ * return res;
+ * }
+ */
+
+
+
+/**********************************************************/
+
+
+void joules2(struct SunGeometryConstDay *sunGeom,
+	     struct SunGeometryVarDay *sunVarGeom,
+	     struct SunGeometryVarSlope *sunSlopeGeom,
+	     struct SolarRadVar *sunRadVar,
+	     struct GridGeometry *gridGeom, unsigned char *horizonpointer,
+	     double latitude, double longitude)
+{
+
+    double s0, dfr, dfr_rad;
+    double ra, dra;
+    int ss = 1;
+    double firstTime;
+    double firstAngle, lastAngle;
+    int srStepNo;
+    double bh;
+    double rr;
+
+    beam_e = 0.;
+    diff_e = 0.;
+    refl_e = 0.;
+    insol_t = 0.;
+
+
+    com_par(sunGeom, sunVarGeom, gridGeom, latitude, longitude);
+
+    if (tt != NULL) {		/*irradiance */
+
+	s0 = lumcline2(sunGeom, sunVarGeom, sunSlopeGeom, gridGeom,
+		       horizonpointer);
+
+	if (sunVarGeom->solarAltitude > 0.) {
+	    if ((!sunVarGeom->isShadow) && (s0 > 0.)) {
+		ra = brad(s0, &bh, sunVarGeom, sunSlopeGeom, sunRadVar);	/* beam radiation */
+		beam_e += ra;
+	    }
+	    else {
+		beam_e = 0.;
+		bh = 0.;
+	    }
+
+	    if ((diff_rad != NULL) || (glob_rad != NULL)) {
+		dra = drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);	/* diffuse rad. */
+		diff_e += dra;
+	    }
+	    if ((refl_rad != NULL) || (glob_rad != NULL)) {
+		if ((diff_rad == NULL) && (glob_rad == NULL))
+		    drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
+		refl_e += rr;	/* reflected rad. */
+	    }
+	}			/* solarAltitude */
+    }
+    else {
+	/* all-day radiation */
+
+	srStepNo = (int)(sunGeom->sunrise_time / step);
+
+	if ((sunGeom->sunrise_time - srStepNo * step) > 0.5 * step) {
+	    firstTime = (srStepNo + 1.5) * step;
+	}
+	else {
+	    firstTime = (srStepNo + 0.5) * step;
+	}
+
+
+	firstAngle = (firstTime - 12) * HOURANGLE;
+	lastAngle = (sunGeom->sunset_time - 12) * HOURANGLE;
+
+
+
+
+	dfr_rad = step * HOURANGLE;
+
+	sunGeom->timeAngle = firstAngle;
+
+	varCount_global = 0;
+
+
+	dfr = step;
+
+	while (ss == 1) {
+
+	    com_par(sunGeom, sunVarGeom, gridGeom, latitude, longitude);
+	    s0 = lumcline2(sunGeom, sunVarGeom, sunSlopeGeom, gridGeom,
+			   horizonpointer);
+
+	    if (sunVarGeom->solarAltitude > 0.) {
+
+		if ((!sunVarGeom->isShadow) && (s0 > 0.)) {
+		    insol_t += dfr;
+		    ra = brad(s0, &bh, sunVarGeom, sunSlopeGeom, sunRadVar);
+		    beam_e += dfr * ra;
+		    ra = 0.;
+		}
+		else {
+		    bh = 0.;
+		}
+		if ((diff_rad != NULL) || (glob_rad != NULL)) {
+		    dra =
+			drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
+		    diff_e += dfr * dra;
+		    dra = 0.;
+		}
+		if ((refl_rad != NULL) || (glob_rad != NULL)) {
+		    if ((diff_rad == NULL) && (glob_rad == NULL)) {
+			drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
+		    }
+		    refl_e += dfr * rr;
+		    rr = 0.;
+		}
+	    }			/* illuminated */
+
+
+	    sunGeom->timeAngle = sunGeom->timeAngle + dfr_rad;
+
+	    if (sunGeom->timeAngle > lastAngle) {
+		ss = 0;		/* we've got the sunset */
+	    }
+	}			/* end of while */
+    }				/* all-day radiation */
+
+}
+
+/*////////////////////////////////////////////////////////////////////// */
+
+
+
+
+/*
+ * void where_is_point(void)
+ * {
+ * double sx, sy;
+ * double dx, dy;
+ * double adx, ady;
+ * int i, j;
+ * 
+ * sx = xx0 * invstepx + TOLER;
+ * sy = yy0 * invstepy + TOLER;
+ * 
+ * i = (int)sx;
+ * j = (int)sy;
+ * if (i < n - 1 && j < m - 1) {
+ * 
+ * dx = xx0 - (double)i *stepx;
+ * dy = yy0 - (double)j *stepy;
+ * 
+ * adx = fabs(dx);
+ * ady = fabs(dy);
+ * 
+ * if ((adx > TOLER) && (ady > TOLER)) {
+ * cube(j, i);
+ * return;
+ * }
+ * else if ((adx > TOLER) && (ady < TOLER)) {
+ * line_x(j, i);
+ * return;
+ * }
+ * else if ((adx < TOLER) && (ady > TOLER)) {
+ * line_y(j, i);
+ * return;
+ * }
+ * else if ((adx < TOLER) && (ady < TOLER)) {
+ * vertex(j, i);
+ * return;
+ * }
+ * 
+ * 
+ * }
+ * else {
+ * func = NULL;
+ * }
+ * }
+ * 
+ */
+
+void where_is_point(double *length, struct SunGeometryVarDay *sunVarGeom,
+		    struct GridGeometry *gridGeom)
+{
+    double sx, sy;
+    double dx, dy;
+    /*              double adx, ady; */
+    int i, j;
+
+    sx = gridGeom->xx0 * invstepx + offsetx;	/* offset 0.5 cell size to get the right cell i, j */
+    sy = gridGeom->yy0 * invstepy + offsety;
+
+    i = (int)sx;
+    j = (int)sy;
+
+    /*      if (i < n-1  && j < m-1) {    to include last row/col */
+    if (i <= n - 1 && j <= m - 1) {
+
+	dx = (double)i *gridGeom->stepx;
+	dy = (double)j *gridGeom->stepy;
+
+	*length = distance(gridGeom->xg0, dx, gridGeom->yg0, dy);	/* dist from orig. grid point to the current grid point */
+
+	sunVarGeom->zp = z[j][i];
+
+	/*
+	 * cube(j, i);
+	 */
+    }
+}
+
+/*
+ * void vertex(jmin, imin)
+ * int jmin, imin;
+ * {
+ * zp = z[jmin][imin];
+ * if ((zp == UNDEFZ))
+ * func = NULL;
+ * }
+ * void line_x(jmin, imin)
+ * int jmin, imin;
+ * {
+ * double c1, c2;
+ * double d1, d2, e1, e2;
+ * e1 = (double)imin *stepx;
+ * e2 = (double)(imin + 1) * stepx;
+ * 
+ * c1 = z[jmin][imin];
+ * c2 = z[jmin][imin + 1];
+ * if (!((c1 == UNDEFZ) || (c2 == UNDEFZ))) {
+ * 
+ * if (dist <= 1.0) {
+ * d1 = (xx0 - e1) / (e2 - e1);
+ * d2 = 1 - d1;
+ * if (d1 < d2)
+ * zp = c1;
+ * else
+ * zp = c2;
+ * }
+ * 
+ * if (dist > 1.0)
+ * zp = AMAX1(c1, c2);
+ * }
+ * else
+ * func = NULL;
+ * }
+ * 
+ * 
+ * void line_y(jmin, imin)
+ * int jmin, imin;
+ * {
+ * double c1, c2;
+ * double d1, d2, e1, e2;
+ * e1 = (double)jmin *stepy;
+ * e2 = (double)(jmin + 1) * stepy;
+ * 
+ * c1 = z[jmin][imin];
+ * c2 = z[jmin + 1][imin];
+ * if (!((c1 == UNDEFZ) || (c2 == UNDEFZ))) {
+ * 
+ * if (dist <= 1.0) {
+ * d1 = (yy0 - e1) / (e2 - e1);
+ * d2 = 1 - d1;
+ * if (d1 < d2)
+ * zp = c1;
+ * else
+ * zp = c2;
+ * }
+ * 
+ * if (dist > 1.0)
+ * zp = AMAX1(c1, c2);
+ * 
+ * }
+ * else
+ * func = NULL;
+ * 
+ * }
+ * 
+ * void cube(jmin, imin)
+ * int jmin, imin;
+ * {
+ * int i, ig = 0;
+ * double x1, x2, y1, y2;
+ * double v[4], vmin = BIG;
+ * double c[4], cmax = -BIG;
+ * 
+ * x1 = (double)imin *stepx;
+ * x2 = x1 + stepx;
+ * 
+ * y1 = (double)jmin *stepy;
+ * y2 = y1 + stepy;
+ * 
+ * v[0] = DISTANCE2(x1, y1);
+ * 
+ * if (v[0] < vmin) {
+ * ig = 0;
+ * vmin = v[0];
+ * }
+ * v[1] = DISTANCE2(x2, y1);
+ * 
+ * if (v[1] < vmin) {
+ * ig = 1;
+ * vmin = v[1];
+ * }
+ * 
+ * v[2] = DISTANCE2(x2, y2);
+ * if (v[2] < vmin) {
+ * ig = 2;
+ * vmin = v[2];
+ * }
+ * 
+ * v[3] = DISTANCE2(x1, y2);
+ * if (v[3] < vmin) {
+ * ig = 3;
+ * vmin = v[3];
+ * }
+ * 
+ * c[0] = z[jmin][imin];
+ * c[1] = z[jmin][imin + 1];
+ * c[2] = z[jmin + 1][imin + 1];
+ * c[3] = z[jmin + 1][imin];
+ * 
+ * 
+ * if (dist <= 1.0) {
+ * 
+ * if (c[ig] != UNDEFZ)
+ * zp = c[ig];
+ * else
+ * func = NULL;
+ * return;
+ * }
+ * 
+ * if (dist > 1.0) {
+ * for (i = 0; i < 4; i++) {
+ * if (c[i] != UNDEFZ) {
+ * cmax = AMAX1(cmax, c[i]);
+ * zp = cmax;
+ * }
+ * else
+ * func = NULL;
+ * }
+ * }
+ * }
+ */
+
+/*
+ * void cube(jmin, imin)
+ * int jmin, imin;
+ * {
+ * zp = z[jmin][imin];
+ * 
+ * }
+ */
+
+void cube(jmin, imin)
+{
+}
+
+
+/*////////////////////////////////////////////////////////////////////// */
+
+void calculate(double singleSlope, double singleAspect, double singleAlbedo,
+	       double singleLinke, struct GridGeometry gridGeom)
+{
+    int i, j, l;
+    /*                      double energy; */
+    int someRadiation;
+    int numRows;
+    int arrayOffset;
+    double lum, q1;
+    double dayRad;
+    double latid_l, cos_u, cos_v, sin_u, sin_v;
+    double sin_phi_l, tan_lam_l;
+    double zmax;
+    double longitTime = 0.;
+    double locTimeOffset;
+    double latitude, longitude;
+    double coslat;
+
+
+    struct SunGeometryConstDay sunGeom;
+    struct SunGeometryVarDay sunVarGeom;
+    struct SunGeometryVarSlope sunSlopeGeom;
+    struct SolarRadVar sunRadVar;
+
+    sunSlopeGeom.slope = singleSlope;
+    sunSlopeGeom.aspect = singleAspect;
+    sunRadVar.alb = singleAlbedo;
+    sunRadVar.linke = singleLinke;
+    sunRadVar.cbh = 1.0;
+    sunRadVar.cdh = 1.0;
+
+    sunGeom.sindecl = sin(declination);
+    sunGeom.cosdecl = cos(declination);
+
+
+    someRadiation = (beam_rad != NULL) || (insol_time != NULL) ||
+	(diff_rad != NULL) || (refl_rad != NULL) || (glob_rad != NULL);
+
+
+    if (incidout != NULL) {
+	lumcl = (float **)G_malloc(sizeof(float *) * (m));
+	for (l = 0; l < m; l++) {
+	    lumcl[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+	for (j = 0; j < m; j++) {
+	    for (i = 0; i < n; i++)
+		lumcl[j][i] = UNDEFZ;
+	}
+    }
+
+    if (beam_rad != NULL) {
+	beam = (float **)G_malloc(sizeof(float *) * (m));
+	for (l = 0; l < m; l++) {
+	    beam[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+
+	for (j = 0; j < m; j++) {
+	    for (i = 0; i < n; i++)
+		beam[j][i] = UNDEFZ;
+	}
+    }
+
+    if (insol_time != NULL) {
+	insol = (float **)G_malloc(sizeof(float *) * (m));
+	for (l = 0; l < m; l++) {
+	    insol[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+
+	for (j = 0; j < m; j++) {
+	    for (i = 0; i < n; i++)
+		insol[j][i] = UNDEFZ;
+	}
+    }
+
+    if (diff_rad != NULL) {
+	diff = (float **)G_malloc(sizeof(float *) * (m));
+	for (l = 0; l < m; l++) {
+	    diff[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+
+	for (j = 0; j < m; j++) {
+	    for (i = 0; i < n; i++)
+		diff[j][i] = UNDEFZ;
+	}
+    }
+
+    if (refl_rad != NULL) {
+	refl = (float **)G_malloc(sizeof(float *) * (m));
+	for (l = 0; l < m; l++) {
+	    refl[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+
+	for (j = 0; j < m; j++) {
+	    for (i = 0; i < n; i++)
+		refl[j][i] = UNDEFZ;
+	}
+    }
+
+    if (glob_rad != NULL) {
+	globrad = (float **)G_malloc(sizeof(float *) * (m));
+	for (l = 0; l < m; l++) {
+	    globrad[l] = (float *)G_malloc(sizeof(float) * (n));
+	}
+
+	for (j = 0; j < m; j++) {
+	    for (i = 0; i < n; i++)
+		globrad[j][i] = UNDEFZ;
+	}
+    }
+
+
+
+    sunRadVar.G_norm_extra = com_sol_const(day);
+
+    numRows = m / numPartitions;
+
+
+    if (useCivilTime()) {
+	/* We need to calculate the deviation of the local solar time from the 
+	 * "local clock time". */
+	dayRad = 2. * M_PI * day / 365.25;
+	locTimeOffset =
+	    -0.128 * sin(dayRad - 0.04887) - 0.165 * sin(2 * dayRad + 0.34383);
+
+	/* Time offset due to timezone as input by user */
+
+	locTimeOffset += civilTime;
+	setTimeOffset(locTimeOffset);
+    }
+    else {
+	setTimeOffset(0.);
+
+    }
+
+
+
+    for (j = 0; j < m; j++) {
+	G_percent(j, m - 1, 2);
+
+	if (j % (numRows) == 0) {
+	    INPUT_part(j, &zmax);
+	    arrayOffset = 0;
+	    shadowoffset = 0;
+
+	}
+	sunVarGeom.zmax = zmax;
+
+
+	for (i = 0; i < n; i++) {
+
+	    if (useCivilTime()) {
+		longitTime = -longitArray[arrayOffset][i] / 15.;
+	    }
+
+
+	    gridGeom.xg0 = gridGeom.xx0 = (double)i *gridGeom.stepx;
+	    gridGeom.yg0 = gridGeom.yy0 = (double)j *gridGeom.stepy;
+
+
+	    gridGeom.xp = xmin + gridGeom.xx0;
+	    gridGeom.yp = ymin + gridGeom.yy0;
+
+	    if (ll_correction) {
+		coslat = cos(deg2rad * gridGeom.yp);
+		coslatsq = coslat * coslat;
+	    }
+
+	    func = NULL;
+
+	    sunVarGeom.z_orig = z1 = sunVarGeom.zp = z[arrayOffset][i];
+
+	    if (sunVarGeom.z_orig != UNDEFZ) {
+		if (aspin != NULL) {
+		    o_orig = o[arrayOffset][i];
+		    if (o[arrayOffset][i] != 0.)
+			sunSlopeGeom.aspect = o[arrayOffset][i] * deg2rad;
+		    else
+			sunSlopeGeom.aspect = UNDEF;
+		}
+		if (slopein != NULL) {
+		    sunSlopeGeom.slope = s[arrayOffset][i] * deg2rad;
+		}
+		if (linkein != NULL) {
+		    sunRadVar.linke = li[arrayOffset][i];
+		    li_max = AMAX1(li_max, sunRadVar.linke);
+		    li_min = AMIN1(li_min, sunRadVar.linke);
+		}
+		if (albedo != NULL) {
+		    sunRadVar.alb = a[arrayOffset][i];
+		    al_max = AMAX1(al_max, sunRadVar.alb);
+		    al_min = AMIN1(al_min, sunRadVar.alb);
+		}
+		if (latin != NULL) {
+		    latitude = la[arrayOffset][i];
+		    la_max = AMAX1(la_max, latitude);
+		    la_min = AMIN1(la_min, latitude);
+		    latitude *= deg2rad;
+		}
+		if ((G_projection() != PROJECTION_LL)) {
+
+			longitude = gridGeom.xp;
+			latitude = gridGeom.yp;
+
+			if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+			    G_fatal_error("Error in pj_do_proj");
+			}
+
+			la_max = AMAX1(la_max, latitude);
+			la_min = AMIN1(la_min, latitude);
+			latitude *= deg2rad;
+			longitude *= deg2rad;
+		}
+		else {	/* ll projection */
+			latitude = gridGeom.yp;
+			longitude = gridGeom.xp;
+			la_max = AMAX1(la_max, latitude);
+			la_min = AMIN1(la_min, latitude);
+			latitude *= deg2rad;
+			longitude *= deg2rad;
+		}
+
+		if (coefbh != NULL) {
+		    sunRadVar.cbh = cbhr[arrayOffset][i];
+		}
+		if (coefdh != NULL) {
+		    sunRadVar.cdh = cdhr[arrayOffset][i];
+		}
+		cos_u = cos(M_PI / 2 - sunSlopeGeom.slope);	/* = sin(slope) */
+		sin_u = sin(M_PI / 2 - sunSlopeGeom.slope);	/* = cos(slope) */
+		cos_v = cos(M_PI / 2 + sunSlopeGeom.aspect);
+		sin_v = sin(M_PI / 2 + sunSlopeGeom.aspect);
+
+		if (tt != NULL)
+		    sunGeom.timeAngle = tim;
+
+		gridGeom.sinlat = sin(-latitude);
+		gridGeom.coslat = cos(-latitude);
+
+		sin_phi_l =
+		    -gridGeom.coslat * cos_u * sin_v + gridGeom.sinlat * sin_u;
+		latid_l = asin(sin_phi_l);
+
+		q1 = gridGeom.sinlat * cos_u * sin_v + gridGeom.coslat * sin_u;
+		tan_lam_l = -cos_u * cos_v / q1;
+		sunSlopeGeom.longit_l = atan(tan_lam_l);
+		sunSlopeGeom.lum_C31_l = cos(latid_l) * sunGeom.cosdecl;
+		sunSlopeGeom.lum_C33_l = sin_phi_l * sunGeom.sindecl;
+
+		if ((incidout != NULL) || someRadiation) {
+		    com_par_const(longitTime, &sunGeom, &gridGeom);
+		    sr_min = AMIN1(sr_min, sunGeom.sunrise_time);
+		    sr_max = AMAX1(sr_max, sunGeom.sunrise_time);
+		    ss_min = AMIN1(ss_min, sunGeom.sunset_time);
+		    ss_max = AMAX1(ss_max, sunGeom.sunset_time);
+
+		}
+
+		if (incidout != NULL) {
+		    com_par(&sunGeom, &sunVarGeom, &gridGeom, latitude,
+			    longitude);
+		    lum =
+			lumcline2(&sunGeom, &sunVarGeom, &sunSlopeGeom,
+				  &gridGeom, horizonarray + shadowoffset);
+		    lum = rad2deg * asin(lum);
+		    lumcl[j][i] = (float)lum;
+		}
+		if (someRadiation) {
+		    joules2(&sunGeom, &sunVarGeom, &sunSlopeGeom, &sunRadVar,
+			    &gridGeom, horizonarray + shadowoffset, latitude,
+			    longitude);
+		    if (beam_rad != NULL)
+			beam[j][i] = (float)beam_e;
+		    if (insol_time != NULL)
+			insol[j][i] = (float)insol_t;
+		    /*  G_debug(3,"\n %f",insol[j][i]); */
+		    if (diff_rad != NULL)
+			diff[j][i] = (float)diff_e;
+		    if (refl_rad != NULL)
+			refl[j][i] = (float)refl_e;
+		    if (glob_rad != NULL)
+			globrad[j][i] = (float)(beam_e + diff_e + refl_e);
+		}
+
+	    }			/* undefs */
+	    shadowoffset += arrayNumInt;
+	}
+	arrayOffset++;
+
+    }
+
+    /* re-use &hist, but try all to initiate it for any case */
+    /*   note this will result in incorrect map titles       */
+    if (incidout != NULL) {
+	G_short_history(incidout, "raster", &hist);
+    }
+    else if (beam_rad != NULL) {
+	G_short_history(beam_rad, "raster", &hist);
+    }
+    else if (diff_rad != NULL) {
+	G_short_history(diff_rad, "raster", &hist);
+    }
+    else if (refl_rad != NULL) {
+	G_short_history(refl_rad, "raster", &hist);
+    }
+    else if (insol_time != NULL) {
+	G_short_history(insol_time, "raster", &hist);
+    }
+    else if (glob_rad != NULL) {
+	G_short_history(glob_rad, "raster", &hist);
+    }
+    else
+	G_fatal_error("Failed to init map history: no output maps requested!");
+
+    sprintf(hist.edhist[0],
+	    " ----------------------------------------------------------------");
+    sprintf(hist.edhist[1], " Day [1-365]:                              %d",
+	    day);
+    hist.edlinecnt = 2;
+
+    if (tt != NULL) {
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Local (solar) time (decimal hr.):         %.4f", timo);
+	hist.edlinecnt++;
+    }
+
+    sprintf(hist.edhist[hist.edlinecnt],
+	    " Solar constant (W/m^2):                   1367");
+    sprintf(hist.edhist[hist.edlinecnt + 1],
+	    " Extraterrestrial irradiance (W/m^2):      %f",
+	    sunRadVar.G_norm_extra);
+    sprintf(hist.edhist[hist.edlinecnt + 2],
+	    " Declination (rad):                        %f", -declination);
+    hist.edlinecnt += 3;
+
+    if (lt != NULL)
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Latitude (deg):                           %.4f",
+		-latitude * rad2deg);
+    else
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Latitude min-max(deg):                    %.4f - %.4f",
+		la_min, la_max);
+    hist.edlinecnt++;
+
+    if (tt != NULL) {
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Sunrise time (hr.):                       %.2f",
+		sunGeom.sunrise_time);
+	sprintf(hist.edhist[hist.edlinecnt + 1],
+		" Sunset time (hr.):                        %.2f",
+		sunGeom.sunset_time);
+	sprintf(hist.edhist[hist.edlinecnt + 2],
+		" Daylight time (hr.):                      %.2f",
+		sunGeom.sunset_time - sunGeom.sunrise_time);
+    }
+    else {
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Sunrise time min-max (hr.):               %.2f - %.2f",
+		sr_min, sr_max);
+	sprintf(hist.edhist[hist.edlinecnt + 1],
+		" Sunset time min-max (hr.):                %.2f - %.2f",
+		ss_min, ss_max);
+	sprintf(hist.edhist[hist.edlinecnt + 2],
+		" Time step (hr.):                          %.4f", step);
+    }
+    hist.edlinecnt += 3;
+
+    if (incidout != NULL || tt != NULL) {
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Solar altitude (deg):                     %.4f",
+		sunVarGeom.solarAltitude * rad2deg);
+	sprintf(hist.edhist[hist.edlinecnt + 1],
+		" Solar azimuth (deg):                      %.4f",
+		sunVarGeom.solarAzimuth * rad2deg);
+	hist.edlinecnt += 2;
+    }
+
+    if (linkein == NULL)
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Linke turbidity factor:                   %.1f",
+		sunRadVar.linke);
+    else
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Linke turbidity factor min-max:           %.1f-%.1f", li_min,
+		li_max);
+    hist.edlinecnt++;
+
+    if (albedo == NULL)
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Ground albedo:                            %.3f",
+		sunRadVar.alb);
+    else
+	sprintf(hist.edhist[hist.edlinecnt],
+		" Ground albedo min-max:                    %.3f-%.3f", al_min,
+		al_max);
+    hist.edlinecnt++;
+
+    sprintf(hist.edhist[hist.edlinecnt],
+	    " -----------------------------------------------------------------");
+    hist.edlinecnt++;
+
+    /* don't call G_write_history() until after G_close_cell() or it just gets overwritten */
+
+}				/* End of ) function */
+
+
+
+double com_declin(int no_of_day)
+{
+    double d1, decl;
+
+    d1 = pi2 * no_of_day / 365.25;
+    decl = asin(0.3978 * sin(d1 - 1.4 + 0.0355 * sin(d1 - 0.0489)));
+    decl = -decl;
+    /*      G_debug(3," declination : %lf\n", decl); */
+
+    return (decl);
+}
+
+
+
+int test(void)
+{
+    /* not finshed yet */
+    int dej;
+
+    G_message("\n ddd: %f", declin);
+    dej = asin(-declin / 0.4093) * 365. / pi2 + 81;
+    /*        dej = asin(-declin/23.35 * deg2rad) / 0.9856 - 284; */
+    /*      dej = dej - 365; */
+    G_message("\n d: %d ", dej);
+    if (dej < day - 5 || dej > day + 5)
+	return 0;
+    else
+	return 1;
+}


Property changes on: grass/branches/develbranch_6/raster/r.sun2/main.c
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.sun2/rsunglobals.h
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/rsunglobals.h	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/rsunglobals.h	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,62 @@
+/*******************************************************************************
+r.sun: rsunglobals.h. This program was writen by Jaro Hofierka in Summer 1993 and re-engineered
+in 1996-1999. In cooperation with Marcel Suri and Thomas Huld from JRC in Ispra
+a new version of r.sun was prepared using ESRA solar radiation formulas.
+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
+*******************************************************************************/
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the
+ *   Free Software Foundation, Inc.,
+ *   59 Temple Place - Suite 330,
+ *   Boston, MA  02111-1307, USA.
+ */
+
+/*v. 2.0 July 2002, NULL data handling, JH */
+/*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
+
+#define EARTHRADIUS 6371000.
+/* undefined value for terrain aspect */
+#define UNDEF    0.        
+/* internal undefined value for NULL */
+#define UNDEFZ   -9999.        
+
+/* Constant for calculating angular loss */
+#define a_r 0.155
+
+extern int varCount_global;
+extern int bitCount_global;
+extern int arrayNumInt;
+
+/*
+extern double xp;
+extern double yp;
+*/
+
+extern double angular_loss_denom;
+
+extern const double invScale;
+extern const double pihalf;
+extern const double pi2;
+extern const double deg2rad;
+extern const double rad2deg;
+
+extern struct pj_info iproj;
+extern struct pj_info oproj;
+
+
+extern void (*func) (int, int);
+


Property changes on: grass/branches/develbranch_6/raster/r.sun2/rsunglobals.h
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.sun2/rsunlib.c
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/rsunlib.c	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/rsunlib.c	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,699 @@
+/*******************************************************************************
+r.sun: rsunlib.c. This program was writen by Jaro Hofierka in Summer 1993 and re-engineered
+in 1996-1999. In cooperation with Marcel Suri and Thomas Huld from JRC in Ispra
+a new version of r.sun was prepared using ESRA solar radiation formulas.
+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
+*******************************************************************************/
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the
+ *   Free Software Foundation, Inc.,
+ *   59 Temple Place - Suite 330,
+ *   Boston, MA  02111-1307, USA.
+ */
+
+/*v. 2.0 July 2002, NULL data handling, JH */
+/*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
+
+#include <stdio.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+#include "sunradstruct.h"
+#include "local_proto.h"
+#include "rsunglobals.h"
+
+
+int civilTimeFlag;
+int useCivilTime()
+	{
+	return civilTimeFlag;
+	}
+void setUseCivilTime(int val)
+	{
+	civilTimeFlag=val;
+	}
+
+
+double angular_loss_denom;
+	
+void setAngularLossDenominator()
+{
+	angular_loss_denom=1./(1-exp(-1./a_r));
+	
+}
+
+
+int useShadowFlag;
+int useShadow()
+	{
+	return useShadowFlag;
+	}
+void setUseShadow(int val)
+	{
+	useShadowFlag=val;
+	}
+
+
+
+int useHorizonDataFlag;
+int useHorizonData()
+	{
+	return useHorizonDataFlag;
+	}
+void setUseHorizonData(int val)
+	{
+	useHorizonDataFlag=val;
+	}
+
+double timeOffset;
+double getTimeOffset()
+	{
+	return timeOffset;
+	}
+void setTimeOffset(double val)
+	{
+	timeOffset=val;
+	}
+
+double horizonInterval;
+double getHorizonInterval()
+	{
+	return horizonInterval;
+	}
+void setHorizonInterval(double val)
+	{
+	horizonInterval=val;
+	}
+
+
+
+double com_sol_const(int no_of_day)
+	{
+		double I0, d1;
+
+		/*  v W/(m*m) */
+		d1 = pi2 * no_of_day / 365.25;
+		I0 = 1367. * (1 + 0.03344 * cos(d1 - 0.048869));
+
+		return I0;
+	}
+
+	
+	
+
+void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
+		   struct GridGeometry *gridGeom)
+	{
+	double pom;
+	double totOffsetTime;
+
+	sungeom->lum_C11 = gridGeom->sinlat * sungeom->cosdecl;
+	sungeom->lum_C13 = -gridGeom->coslat * sungeom->sindecl;
+	sungeom->lum_C22 = sungeom->cosdecl;
+	sungeom->lum_C31 = gridGeom->coslat * sungeom->cosdecl;
+	sungeom->lum_C33 = gridGeom->sinlat * sungeom->sindecl;
+
+    if (fabs(sungeom->lum_C31) >= EPS) {
+        totOffsetTime = timeOffset + longitTime;
+        
+        if(useCivilTime())
+            {
+            sungeom->timeAngle -= totOffsetTime*HOURANGLE;
+            }
+    	pom = -sungeom->lum_C33 / sungeom->lum_C31;
+    	if (fabs(pom) <= 1) {
+        	pom = acos(pom);
+        	pom = (pom * 180) / M_PI;
+        	sungeom->sunrise_time = (90 - pom) / 15 + 6;
+        	sungeom->sunset_time = (pom - 90) / 15 + 18;
+    	    }
+        else {
+            if (pom < 0) {
+                /*        G_debug(3,"\n Sun is ABOVE the surface during the whole day"); */
+                sungeom->sunrise_time = 0;
+                sungeom->sunset_time = 24;
+                }
+            else {
+            /*                G_debug(3,"\n The sun is BELOW the surface during the whole day"); */
+                if (fabs(pom) - 1 <= EPS) {
+                    sungeom->sunrise_time = 12;
+                    sungeom->sunset_time = 12;
+                }
+            }
+        }
+    }
+
+}
+
+
+
+
+
+void com_par( struct SunGeometryConstDay *sungeom, 
+	      struct SunGeometryVarDay *sunVarGeom, struct GridGeometry *gridGeom,
+		double latitude, double longitude)
+{
+	double pom, xpom, ypom;
+
+	double costimeAngle;
+	double lum_Lx, lum_Ly;
+
+	double newLatitude, newLongitude;
+	double inputAngle;
+   	double delt_lat, delt_lon;
+   	double delt_east, delt_nor;
+	double delt_dist;
+
+
+	costimeAngle = cos(sungeom->timeAngle);
+
+
+
+	lum_Lx = -sungeom->lum_C22 * sin(sungeom->timeAngle);
+	lum_Ly = sungeom->lum_C11 * costimeAngle + sungeom->lum_C13;
+	sunVarGeom->sinSolarAltitude = sungeom->lum_C31 * costimeAngle + sungeom->lum_C33;
+
+	if (fabs(sungeom->lum_C31) < EPS) 
+		{
+		if (fabs(sunVarGeom->sinSolarAltitude) >= EPS) 
+			{
+			if (sunVarGeom->sinSolarAltitude > 0) 
+				{
+        		/*                        G_debug(3,"\tSun is ABOVE area during the whole day"); */
+        			sungeom->sunrise_time = 0;
+        			sungeom->sunset_time = 24;
+        			}
+            		else 
+				{
+        			sunVarGeom->solarAltitude = 0.;
+        			sunVarGeom->solarAzimuth = UNDEF;
+        			return;
+                		}
+            		}
+        	else 
+			{
+            /*                      G_debug(3,"\tThe Sun is ON HORIZON during the whole day"); */
+            		sungeom->sunrise_time = 0;
+            		sungeom->sunset_time = 24;
+            		}
+        	}
+
+	sunVarGeom->solarAltitude = asin(sunVarGeom->sinSolarAltitude);        /* vertical angle of the sun */
+	/* sinSolarAltitude is sin(solarAltitude) */
+
+	xpom = lum_Lx * lum_Lx;
+	ypom = lum_Ly * lum_Ly;
+	pom = sqrt(xpom + ypom);
+
+
+	if (fabs(pom) > EPS) 
+		{
+        	sunVarGeom->solarAzimuth = lum_Ly / pom;
+        	sunVarGeom->solarAzimuth = acos(sunVarGeom->solarAzimuth);        /* horiz. angle of the Sun */
+        /*                      solarAzimuth *= RAD; */
+        	if (lum_Lx < 0)
+            		sunVarGeom->solarAzimuth = pi2 - sunVarGeom->solarAzimuth;
+        	}
+    	else 
+		{
+	        sunVarGeom->solarAzimuth = UNDEF;
+        	}
+
+
+
+    	if (sunVarGeom->solarAzimuth < 0.5 * M_PI)
+        	sunVarGeom->sunAzimuthAngle = 0.5 * M_PI - sunVarGeom->solarAzimuth;
+    	else
+        	sunVarGeom->sunAzimuthAngle = 2.5 * M_PI - sunVarGeom->solarAzimuth;
+
+
+	inputAngle=sunVarGeom->sunAzimuthAngle+pihalf;
+	inputAngle=(inputAngle>=pi2)?inputAngle-pi2:inputAngle;
+
+
+	delt_lat = -0.0001*cos(inputAngle); /* Arbitrary small distance in latitude */
+	delt_lon = 0.0001*sin(inputAngle)/cos(latitude);
+
+	newLatitude = (latitude+delt_lat)*rad2deg;
+	newLongitude = (longitude+delt_lon)*rad2deg;
+
+
+	if ((G_projection() != PROJECTION_LL))
+             	{
+	        if (pj_do_proj(&newLongitude, &newLatitude, &oproj, &iproj) <0)
+			{
+               		G_fatal_error("Error in pj_do_proj");
+        		}
+		}
+
+	delt_east=newLongitude-gridGeom->xp;
+	delt_nor=newLatitude-gridGeom->yp;
+
+	delt_dist=sqrt(delt_east*delt_east+delt_nor*delt_nor);
+
+
+
+	sunVarGeom->stepsinangle = gridGeom->stepxy*delt_nor/delt_dist;
+	sunVarGeom->stepcosangle = gridGeom->stepxy*delt_east/delt_dist;
+
+/*
+	sunVarGeom->stepsinangle = stepxy * sin(sunVarGeom->sunAzimuthAngle);
+	sunVarGeom->stepcosangle = stepxy * cos(sunVarGeom->sunAzimuthAngle);
+*/	
+	sunVarGeom->tanSolarAltitude = tan(sunVarGeom->solarAltitude);
+
+	return;
+
+}
+
+
+int searching(double *length, struct SunGeometryVarDay *sunVarGeom, 
+	      struct GridGeometry *gridGeom)
+{
+	double z2;
+	double curvature_diff;
+	int succes = 0;
+
+	if (sunVarGeom->zp == UNDEFZ)
+		return 0;
+
+
+	gridGeom->yy0 += sunVarGeom->stepsinangle;
+	gridGeom->xx0 += sunVarGeom->stepcosangle;
+	if (((gridGeom->xx0+(0.5*gridGeom->stepx)) < 0) 
+		|| ((gridGeom->xx0+(0.5*gridGeom->stepx)) > gridGeom->deltx) 
+		|| ((gridGeom->yy0+(0.5*gridGeom->stepy)) < 0) 
+		|| ((gridGeom->yy0+(0.5*gridGeom->stepy)) > gridGeom->delty))
+		succes=3;
+	else
+		succes=1;
+
+
+	if (succes == 1) 
+	{
+		where_is_point(length, sunVarGeom,gridGeom);
+		if (func == NULL)
+		{
+			gridGeom->xx0 = gridGeom->xg0;
+			gridGeom->yy0 = gridGeom->yg0;
+			return (3);
+		}
+		curvature_diff = EARTHRADIUS*(1.-cos(*length/EARTHRADIUS));
+        
+		z2 = sunVarGeom->z_orig + curvature_diff + *length * sunVarGeom->tanSolarAltitude;
+		if (z2 < sunVarGeom->zp)
+			succes = 2;        /* shadow */
+		if (z2 > sunVarGeom->zmax)
+			succes = 3;        /* no test needed all visible */
+	}
+
+	if(succes!=1)
+	{
+		gridGeom->xx0 = gridGeom->xg0;
+		gridGeom->yy0 = gridGeom->yg0;
+	}
+	return (succes);
+}
+
+
+
+
+double lumcline2( 
+	struct SunGeometryConstDay *sungeom,
+	struct SunGeometryVarDay *sunVarGeom,
+	struct SunGeometryVarSlope *sunSlopeGeom, 
+	struct GridGeometry *gridGeom,
+		unsigned char *horizonpointer)
+	{
+	double s = 0;
+	double length;
+	int r = 0;
+
+	int lowPos, highPos;
+	double timeoffset, horizPos;
+	double horizonHeight;
+
+	func = cube;
+	sunVarGeom->isShadow = 0;
+
+	if (useShadow()) 
+		{
+        	length = 0;
+        
+		if(useHorizonData())
+			{
+			/* Start is due east, sungeom->timeangle = -pi/2 */		
+/*
+		timeoffset = sungeom->timeAngle+pihalf;
+*/
+		timeoffset = sunVarGeom->sunAzimuthAngle;
+
+/*
+		if(timeoffset<0.)
+			timeoffset+=pi2;
+		else if(timeoffset>pi2)
+			timeoffset-=pi2;
+		horizPos = arrayNumInt - timeoffset/horizonInterval;
+*/			
+		
+			horizPos = timeoffset/getHorizonInterval();
+
+
+			lowPos = (int)horizPos;
+			highPos=lowPos+1;
+			if(highPos==arrayNumInt)
+				{
+				highPos=0;
+				}
+			horizonHeight = invScale*(
+				(1.-(horizPos-lowPos))*horizonpointer[lowPos]
+				+(horizPos-lowPos)
+				*horizonpointer[highPos]);
+			sunVarGeom->isShadow = (horizonHeight>sunVarGeom->solarAltitude);			
+
+              		if (!sunVarGeom->isShadow)
+                		{
+/*
+                		if (z_orig != UNDEFZ) 
+					{
+  		                	s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; 
+
+					}
+                		else
+                			{
+		                	s = sunVarGeom->sinSolarAltitude;
+                			}
+*/
+	                	s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; /* Jenco */
+               			}
+
+
+			}  /*  End if useHorizonData() */
+	        else
+	        	{
+	       		while ((r = searching(&length, sunVarGeom, gridGeom)) == 1) 
+	       			{
+				if (r == 3)
+            				break;        /* no test is needed */
+        			}
+
+
+
+
+	                  if (r == 2) 
+		  		{
+				sunVarGeom->isShadow = 1; /* shadow */
+                  		} 
+			   else
+                		{
+
+/*
+		               if (z_orig != UNDEFZ) 
+					{
+			
+			                s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; 
+                			}
+		                else
+                			{
+                			s = sunVarGeom->sinSolarAltitude;
+                			}
+*/
+		                s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; /* Jenco */
+               			}
+       			}
+        	}
+        else
+            	{
+		/*
+            	if (z_orig != UNDEFZ) 
+			{
+              		s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; 
+            		}
+            	else
+			{
+			s = sunVarGeom->sinSolarAltitude;
+			}
+*/
+      		s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; /* Jenco */
+            
+            
+            	}
+
+      	if (s < 0) return 0.;
+      	return (s);
+	}
+
+
+
+double brad(double sh, double *bh, struct SunGeometryVarDay *sunVarGeom,
+		struct SunGeometryVarSlope *sunSlopeGeom, 
+		struct SolarRadVar *sunRadVar)
+{
+	double opticalAirMass, airMass2Linke, rayl, br;
+	double drefract, temp1, temp2, h0refract;
+	double elevationCorr;
+
+	double locSolarAltitude;
+	
+	locSolarAltitude=sunVarGeom->solarAltitude;
+
+
+	elevationCorr = exp(-sunVarGeom->z_orig / 8434.5);
+	temp1 = 0.1594 + locSolarAltitude * (1.123 + 0.065656 * locSolarAltitude);
+	temp2 = 1. + locSolarAltitude * (28.9344 + 277.3971 * locSolarAltitude);
+	drefract = 0.061359 * temp1 / temp2;    /* in radians */
+	h0refract = locSolarAltitude + drefract;
+	opticalAirMass = elevationCorr / (sin(h0refract) +
+	  0.50572 * pow(h0refract * rad2deg + 6.07995, -1.6364));
+	airMass2Linke = 0.8662 * sunRadVar->linke;
+	if (opticalAirMass <= 20.)
+		{
+		rayl = 1. / (6.6296 +
+	  		opticalAirMass * (1.7513 +
+	    		opticalAirMass * (-0.1202 + opticalAirMass * (0.0065 - opticalAirMass * 0.00013))));
+		}
+	else
+		{
+		rayl = 1. / (10.4 + 0.718 * opticalAirMass);
+		}
+	*bh = sunRadVar->cbh * sunRadVar->G_norm_extra * sunVarGeom->sinSolarAltitude * exp(-rayl * opticalAirMass * airMass2Linke);
+	if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.)
+		br = *bh * sh / sunVarGeom->sinSolarAltitude;
+	else
+		br = *bh;
+
+	return (br);
+}
+
+double brad_angle_loss(double sh, double *bh, struct SunGeometryVarDay *sunVarGeom,
+		       struct SunGeometryVarSlope *sunSlopeGeom, 
+		       struct SolarRadVar *sunRadVar)
+{
+	double p, opticalAirMass, airMass2Linke, rayl, br;
+	double drefract, temp1, temp2, h0refract;
+
+	double locSolarAltitude;
+	
+	locSolarAltitude=sunVarGeom->solarAltitude;
+
+
+	p = exp(-sunVarGeom->z_orig / 8434.5);
+	temp1 = 0.1594 + locSolarAltitude * (1.123 + 0.065656 * locSolarAltitude);
+	temp2 = 1. + locSolarAltitude * (28.9344 + 277.3971 * locSolarAltitude);
+	drefract = 0.061359 * temp1 / temp2;    /* in radians */
+	h0refract = locSolarAltitude + drefract;
+	opticalAirMass = p / (sin(h0refract) +
+			0.50572 * pow(h0refract * rad2deg + 6.07995, -1.6364));
+	airMass2Linke = 0.8662 * sunRadVar->linke;
+	if (opticalAirMass <= 20.)
+		rayl =
+				1. / (6.6296 +
+				opticalAirMass * (1.7513 +
+				opticalAirMass * (-0.1202 + opticalAirMass * (0.0065 - opticalAirMass * 0.00013))));
+	else
+		rayl = 1. / (10.4 + 0.718 * opticalAirMass);
+	*bh = sunRadVar->cbh * sunRadVar->G_norm_extra * sunVarGeom->sinSolarAltitude * exp(-rayl * opticalAirMass * airMass2Linke);
+	if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.)
+		br = *bh * sh / sunVarGeom->sinSolarAltitude;
+	else
+		br = *bh;
+
+	br *= (1-exp(-sh/a_r))*angular_loss_denom;
+    
+	return (br);
+}
+
+
+
+double drad(double sh, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
+		struct SunGeometryVarSlope *sunSlopeGeom, 
+		struct SolarRadVar *sunRadVar)
+{
+    double tn, fd, fx = 0., A1, A2, A3, A1b;
+    double r_sky, kb, dr, gh, a_ln, ln, fg;
+    double dh;
+    double cosslope, sinslope;
+	double locSinSolarAltitude;
+	double locLinke;
+	
+	locLinke = sunRadVar->linke;
+	locSinSolarAltitude=sunVarGeom->sinSolarAltitude;
+
+    cosslope = cos(sunSlopeGeom->slope);
+    sinslope = sin(sunSlopeGeom->slope);
+
+    tn = -0.015843 + locLinke * (0.030543 + 0.0003797 * locLinke);
+    A1b = 0.26463 + locLinke * (-0.061581 + 0.0031408 * locLinke);
+    if (A1b * tn < 0.0022)
+        A1 = 0.0022 / tn;
+    else
+        A1 = A1b;
+    A2 = 2.04020 + locLinke * (0.018945 - 0.011161 * locLinke);
+    A3 = -1.3025 + locLinke * (0.039231 + 0.0085079 * locLinke);
+
+    fd = A1 + A2 * locSinSolarAltitude + A3 * locSinSolarAltitude * locSinSolarAltitude;
+    dh = sunRadVar->cdh * sunRadVar->G_norm_extra * fd * tn;
+    gh = bh + dh;
+    if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.) {
+        kb = bh / (sunRadVar->G_norm_extra * locSinSolarAltitude);
+    	r_sky = (1. + cosslope) / 2.;
+    	a_ln = sunVarGeom->solarAzimuth - sunSlopeGeom->aspect;
+    	ln = a_ln;
+    	if (a_ln > M_PI)
+        	ln = a_ln - pi2;
+    	else if (a_ln < -M_PI)
+        	ln = a_ln + pi2;
+    	a_ln = ln;
+    	fg = sinslope - sunSlopeGeom->slope * cosslope -
+        	M_PI * sin(0.5*sunSlopeGeom->slope) * sin(0.5*sunSlopeGeom->slope);
+    	if ((sunVarGeom->isShadow == 1) || sh <= 0.)
+        	fx = r_sky + fg * 0.252271;
+    	else if (sunVarGeom->solarAltitude >= 0.1) {
+        	fx = ((0.00263 - kb * (0.712 + 0.6883 * kb)) * fg + r_sky) * (1. -
+                                    	  kb) +
+        	kb * sh / locSinSolarAltitude;
+    	}
+    	else if (sunVarGeom->solarAltitude < 0.1)
+        	fx = ((0.00263 - 0.712 * kb - 0.6883 * kb * kb) * fg +
+        	  r_sky) * (1. - kb) + kb * sinslope * cos(a_ln) / (0.1 -
+                                    	  0.008 *
+                                    	  sunVarGeom->solarAltitude);
+    	dr = dh * fx;
+    	/* refl. rad */
+    	*rr = sunRadVar->alb * gh * (1 - cosslope) / 2.;
+    }
+    else {            /* plane */
+        dr = dh;
+        *rr = 0.;
+    }
+    return (dr);
+}
+
+#define c2 -0.074
+
+double drad_angle_loss(double sh, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
+		       struct SunGeometryVarSlope *sunSlopeGeom, 
+		       struct SolarRadVar *sunRadVar)
+{
+	double dh;
+	double tn, fd, fx = 0., A1, A2, A3, A1b;
+	double r_sky, kb, dr, gh, a_ln, ln, fg;
+	double cosslope, sinslope;
+
+	double diff_coeff_angleloss;
+	double refl_coeff_angleloss;
+	double c1;
+	double diff_loss_factor, refl_loss_factor;
+	double locSinSolarAltitude;
+	double locLinke;
+	
+	locSinSolarAltitude=sunVarGeom->sinSolarAltitude;
+	locLinke = sunRadVar->linke;
+	cosslope = cos(sunSlopeGeom->slope);
+	sinslope = sin(sunSlopeGeom->slope);
+
+
+	tn = -0.015843 + locLinke * (0.030543 + 0.0003797 * locLinke);
+	A1b = 0.26463 + locLinke * (-0.061581 + 0.0031408 * locLinke);
+	if (A1b * tn < 0.0022)
+		A1 = 0.0022 / tn;
+	else
+		A1 = A1b;
+	A2 = 2.04020 + locLinke * (0.018945 - 0.011161 * locLinke);
+	A3 = -1.3025 + locLinke * (0.039231 + 0.0085079 * locLinke);
+
+	fd = A1 + A2 * locSinSolarAltitude + A3 * locSinSolarAltitude * locSinSolarAltitude;
+	dh = sunRadVar->cdh * sunRadVar->G_norm_extra * fd * tn;
+	gh = bh + dh;
+	if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.) {
+		kb = bh / (sunRadVar->G_norm_extra * locSinSolarAltitude);
+		r_sky = (1. + cosslope) / 2.;
+		a_ln = sunVarGeom->solarAzimuth - sunSlopeGeom->aspect;
+		ln = a_ln;
+		if (a_ln > M_PI)
+			ln = a_ln - pi2;
+		else if (a_ln < -M_PI)
+			ln = a_ln + pi2;
+		a_ln = ln;
+		fg = sinslope - sunSlopeGeom->slope * cosslope -
+				M_PI * sin(sunSlopeGeom->slope / 2.) * sin(sunSlopeGeom->slope / 2.);
+		if ((sunVarGeom->isShadow) || (sh <= 0.))
+			fx = r_sky + fg * 0.252271;
+		else if (sunVarGeom->solarAltitude >= 0.1) {
+			fx = ((0.00263 - kb * (0.712 + 0.6883 * kb)) * fg + r_sky) * (1. -
+					kb) +
+					kb * sh / locSinSolarAltitude;
+		}
+		else if (sunVarGeom->solarAltitude < 0.1)
+			fx = ((0.00263 - 0.712 * kb - 0.6883 * kb * kb) * fg +
+					r_sky) * (1. - kb) + kb * sinslope * cos(a_ln) / (0.1 -
+					0.008 *
+					sunVarGeom->solarAltitude);
+		dr = dh * fx;
+		/* refl. rad */
+		*rr = sunRadVar->alb * gh * (1 - cosslope) / 2.;
+	}
+	else {            /* plane */
+		dr = dh;
+		*rr = 0.;
+	}
+
+	c1 = 4./(3.*M_PI);
+	diff_coeff_angleloss = sinslope 
+			+ (M_PI-sunSlopeGeom->slope-sinslope)/(1+cosslope);
+	refl_coeff_angleloss = sinslope 
+			+ (sunSlopeGeom->slope-sinslope)/(1-cosslope);
+	
+	diff_loss_factor 
+			= 1. - exp(-(c1*diff_coeff_angleloss
+			+c2*diff_coeff_angleloss*diff_coeff_angleloss)
+			/a_r);
+	refl_loss_factor 
+			= 1. - exp(-(c1*refl_coeff_angleloss
+			+c2*refl_coeff_angleloss*refl_coeff_angleloss)
+			/a_r);
+	
+	dr *= diff_loss_factor;
+	*rr *= refl_loss_factor;
+
+
+
+	return (dr);
+}
+
+


Property changes on: grass/branches/develbranch_6/raster/r.sun2/rsunlib.c
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/raster/r.sun2/sunradstruct.h
===================================================================
--- grass/branches/develbranch_6/raster/r.sun2/sunradstruct.h	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.sun2/sunradstruct.h	2008-11-27 21:45:30 UTC (rev 34554)
@@ -0,0 +1,109 @@
+/*******************************************************************************
+r.sun: sunradstruct.h. This program was writen by Jaro Hofierka in Summer 1993 and re-engineered
+in 1996-1999. In cooperation with Marcel Suri and Thomas Huld from JRC in Ispra
+a new version of r.sun was prepared using ESRA solar radiation formulas.
+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
+*******************************************************************************/
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the
+ *   Free Software Foundation, Inc.,
+ *   59 Temple Place - Suite 330,
+ *   Boston, MA  02111-1307, USA.
+ */
+
+/*v. 2.0 July 2002, NULL data handling, JH */
+/*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
+
+#define EPS      1.e-4
+#define HOURANGLE M_PI/12.
+
+
+struct SunGeometryConstDay 
+	{
+	double lum_C11;
+	double lum_C13;
+	double lum_C22;
+	double lum_C31;
+	double lum_C33;
+	double sunrise_time;
+	double sunset_time;
+	double timeAngle;
+	double sindecl;
+	double cosdecl;
+	
+	};
+
+
+struct SunGeometryVarDay 
+	{
+	int isShadow;
+	double z_orig;
+	double zmax;
+	double zp;
+	double solarAltitude;
+	double sinSolarAltitude;
+	double tanSolarAltitude;
+	double solarAzimuth;
+	double sunAzimuthAngle;
+	double stepsinangle;
+	double stepcosangle;
+	
+	};
+
+
+struct SunGeometryVarSlope 
+	{
+	double longit_l; /* The "longitude" difference between the inclined */
+			/* and orientated plane and the instantaneous solar position */
+	double lum_C31_l;
+	double lum_C33_l;
+	double slope;
+	double aspect;
+	
+	};
+
+
+
+struct SolarRadVar 
+	{
+	double cbh;
+	double cdh;
+	double linke;
+	double G_norm_extra;
+	double alb;
+	
+	};
+
+	
+	
+struct GridGeometry
+	{
+	double xp;
+	double yp;
+	double xx0;
+	double yy0;
+	double xg0;
+	double yg0;
+	double stepx;
+	double stepy;
+	double deltx;
+	double delty;
+	double stepxy;
+	double sinlat;
+	double coslat;
+	
+	};


Property changes on: grass/branches/develbranch_6/raster/r.sun2/sunradstruct.h
___________________________________________________________________
Name: svn:eol-style
   + native



More information about the grass-commit mailing list