[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 "hits" 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 > 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>
+© 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 – 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 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec annual<br>
+mountains 1.5 1.6 1.8 1.9 2.0 2.3 2.3 2.3 2.1 1.8 1.6 1.5 1.90 <br>
+rural 2.1 2.2 2.5 2.9 3.2 3.4 3.5 3.3 2.9 2.6 2.3 2.2 2.75 <br>
+city 3.1 3.2 3.5 4.0 4.2 4.3 4.4 4.3 4.0 3.6 3.3 3.1 3.75 <br>
+industrial 4.1 4.3 4.7 5.3 5.5 5.7 5.8 5.7 5.3 4.9 4.5 4.2 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á analza a digitálne modely georelié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&D in the European Community, series F – 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’ É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>
+
+© 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