[GRASS-SVN] r70825 - in grass-addons/grass7/raster: . r.sun.mp
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Apr 2 21:54:36 PDT 2017
Author: jhofier
Date: 2017-04-02 21:54:36 -0700 (Sun, 02 Apr 2017)
New Revision: 70825
Added:
grass-addons/grass7/raster/r.sun.mp/
grass-addons/grass7/raster/r.sun.mp/Makefile
grass-addons/grass7/raster/r.sun.mp/README
grass-addons/grass7/raster/r.sun.mp/TODO
grass-addons/grass7/raster/r.sun.mp/dalloc.c
grass-addons/grass7/raster/r.sun.mp/description-pak
grass-addons/grass7/raster/r.sun.mp/local_proto.h
grass-addons/grass7/raster/r.sun.mp/main.c
grass-addons/grass7/raster/r.sun.mp/r.sun.mp.html
grass-addons/grass7/raster/r.sun.mp/rsunglobals.h
grass-addons/grass7/raster/r.sun.mp/rsunlib.c
grass-addons/grass7/raster/r.sun.mp/sunradstruct.h
Log:
r.sun with openmp suppport
Added: grass-addons/grass7/raster/r.sun.mp/Makefile
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/Makefile 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.sun.mp
+
+LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PROJLIB)
+DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP)
+EXTRA_INC = $(PROJINC) $(OCLINCPATH)
+EXTRA_CFLAGS=-fopenmp
+EXTRA_LIBS = $(OCLLIB) $(GISLIB) -lgomp $(MATHLIB)
+#needed? LIBES += $(OCLLIBPATH) or EXTRA_LDFLAGS = $(OCLLIBPATH)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.sun.mp/README
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/README (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/README 2017-04-03 04:54:36 UTC (rev 70825)
@@ -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-addons/grass7/raster/r.sun.mp/TODO
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/TODO (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/TODO 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,34 @@
+TODO
+
+-------
+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
+
+####
+Update
+ http://grasswiki.osgeo.org/wiki/R.sun
+
+####
+Fix http://trac.osgeo.org/grass/ticket/498
+
+pseudo-data test-case
+http://trac.osgeo.org/grass/ticket/498#comment:22
+#spearfish (further north than NC so more defined shadows)
+g.region -d
+r.mapcalc "undulates = (2 + sin( row() * 2 ) + cos( col() * 2 )) * 500"
+r.colors undulates color=bcyr
+nviz undulates
+
+DAY=355
+time r.sun elevin=undulates day=$DAY step=0.05 \
+ beam_rad=rad_test.und.355.beam \
+ diff_rad=rad_test.und.355.diff \
+ refl_rad=rad_test.und.355.refl
+
Added: grass-addons/grass7/raster/r.sun.mp/dalloc.c
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/dalloc.c (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/dalloc.c 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,197 @@
+
+/**
+ * \file dalloc.c
+ *
+ * \brief Matrix memory management functions.
+ *
+ * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ *
+ * \author GRASS GIS Development Team
+ *
+ * \date 2004-2006
+ */
+
+#include <stdlib.h>
+#include <grass/gis.h>
+
+
+/**
+ * \fn double *G_alloc_vector (size_t n)
+ *
+ * \brief Vector matrix memory allocation.
+ *
+ * Allocate a vector (array) of <b>n</b> doubles initialized to zero.
+ *
+ * \param[in] n size of vector to allocate
+ * \return double *
+ */
+
+double *G_alloc_vector(size_t n)
+{
+ return (double *)G_calloc(n, sizeof(double));
+}
+
+
+/**
+ * \fn double **G_alloc_matrix (int rows,int cols)
+ *
+ * \brief Matrix memory allocation.
+ *
+ * Allocate a matrix of <b>rows</b> by <b>cols</b> doubles initialized
+ * to zero.
+ *
+ * \param[in] rows number of rows in matrix
+ * \param[in] cols number of columns in matrix
+ * \return double **
+ */
+
+double **G_alloc_matrix(int rows, int cols)
+{
+ double **m;
+ int i;
+
+ m = (double **)G_calloc(rows, sizeof(double *));
+ m[0] = (double *)G_calloc((size_t) rows * cols, sizeof(double));
+ for (i = 1; i < rows; i++)
+ m[i] = m[i - 1] + cols;
+
+ return m;
+}
+
+
+/**
+ * \fn float *G_alloc_fvector (size_t n)
+ *
+ * \brief Floating point vector memory allocation.
+ *
+ * Allocate a vector (array) of <b>n</b> floats initialized to zero.
+ *
+ * \param[in] n size of vector to allocate
+ * \return float *
+ */
+
+float *G_alloc_fvector(size_t n)
+{
+ return (float *)G_calloc(n, sizeof(float));
+}
+
+
+/**
+ * \fn float **G_alloc_fmatrix (int rows, int cols)
+ *
+ * \brief Floating point matrix memory allocation.
+ *
+ * Allocate a matrix of <b>rows</b> by <b>cols</b> floats initialized
+ * to zero.
+ *
+ * \param[in] rows number of rows in matrix
+ * \param[in] cols number of columns in matrix
+ * \return float **
+ */
+
+float **G_alloc_fmatrix(int rows, int cols)
+{
+ float **m;
+ int i;
+
+ m = (float **)G_calloc(rows, sizeof(float *));
+ m[0] = (float *)G_calloc((size_t) rows * cols, sizeof(float));
+ for (i = 1; i < rows; i++)
+ m[i] = m[i - 1] + cols;
+
+ return m;
+}
+
+
+/**
+ * \fn void G_free_vector (double *v)
+ *
+ * \brief Vector memory deallocation.
+ *
+ * Deallocate a vector (array) of doubles.
+ *
+ * \param[in,out] v vector to free
+ * \return void
+ */
+
+void G_free_vector(double *v)
+{
+ G_free(v);
+ v = NULL;
+
+ return;
+}
+
+
+/**
+ * \fn void G_free_fvector (float *v)
+ *
+ * \brief Vector memory deallocation.
+ *
+ * Deallocate a vector (array) of floats.
+ *
+ * \param[in,out] v vector to free
+ * \return void
+ */
+
+void G_free_fvector(float *v)
+{
+ G_free(v);
+ v = NULL;
+
+ return;
+}
+
+
+/**
+ * \fn void G_free_matrix (double **m)
+ *
+ * \brief Matrix memory deallocation.
+ *
+ * Deallocate a matrix of doubles.
+ *
+ * \param[in,out] m matrix to free
+ * \return void
+ */
+
+void G_free_matrix(double **m)
+{
+ G_free(m[0]);
+ G_free(m);
+ m = NULL;
+
+ return;
+}
+
+
+/**
+ * \fn void G_free_fmatrix (float **m)
+ *
+ * \brief Floating point matrix memory deallocation.
+ *
+ * Deallocate a matrix of floats.
+ *
+ * \param[in,out] m matrix to free
+ * \return void
+ */
+
+void G_free_fmatrix(float **m)
+{
+ G_free(m[0]);
+ G_free(m);
+ m = NULL;
+
+ return;
+}
Added: grass-addons/grass7/raster/r.sun.mp/description-pak
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/description-pak (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/description-pak 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1 @@
+Package created with checkinstall 1.6.2
Added: grass-addons/grass7/raster/r.sun.mp/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/local_proto.h (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/local_proto.h 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,180 @@
+/* main.c */
+
+
+void where_is_point(double *length, struct SunGeometryVarDay *sunVarGeom,
+ struct GridGeometry *gridGeom);
+void where_is_point_parallel(double *length, double *sunVarGeom_zp,
+
+ double *gridGeom_xx0,
+ double *gridGeom_yy0,
+ double *gridGeom_xg0,
+ double *gridGeom_yg0,
+ double *gridGeom_stepx,
+ double *gridGeom_stepy);
+int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
+ struct GridGeometry *gridGeom);
+int searching_parallel(double *length,
+
+ double *sunVarGeom_z_orig,
+ double *sunVarGeom_zmax,
+ double *sunVarGeom_zp,
+ double *sunVarGeom_tanSolarAltitude,
+ double *sunVarGeom_stepsinangle,
+ double *sunVarGeom_stepcosangle,
+
+ double *gridGeom_xx0,
+ double *gridGeom_yy0,
+ double *gridGeom_xg0,
+ double *gridGeom_yg0,
+ double *gridGeom_stepx,
+ double *gridGeom_stepy,
+ double *gridGeom_deltx,
+ double *gridGeom_delty);
+
+
+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 brad_parallel(double sh, double *bh,
+
+ double *sunVarGeom_z_orig,
+ double *sunVarGeom_solarAltitude,
+ double *sunVarGeom_sinSolarAltitude,
+
+ double *sunSlopeGeom_slope,
+ double *sunSlopeGeom_aspect,
+
+ double *sunRadVar_cbh,
+ double *sunRadVar_linke,
+ double *sunRadVar_G_norm_extra);
+double drad(double, double bh, double *rr,
+ struct SunGeometryVarDay *sunVarGeom,
+ struct SunGeometryVarSlope *sunSlopeGeom,
+ struct SolarRadVar *sunRadVar);
+double drad_parallel(double sh, double bh, double *rr,
+
+ int *sunVarGeom_isShadow,
+ double *sunVarGeom_solarAltitude,
+ double *sunVarGeom_sinSolarAltitude,
+ double *sunVarGeom_solarAzimuth,
+
+ double *sunSlopeGeom_aspect,
+ double *sunSlopeGeom_slope,
+
+ double *sunRadVar_cdh,
+ double *sunRadVar_linke,
+ double *sunRadVar_G_norm_extra,
+ double *sunRadVar_alb);
+
+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_parallel(double *sunGeom_lum_C11,
+ double *sunGeom_lum_C13,
+ double *sunGeom_lum_C22,
+ double *sunGeom_lum_C31,
+ double *sunGeom_lum_C33,
+ double *sunGeom_sunrise_time,
+ double *sunGeom_sunset_time,
+ double *sunGeom_timeAngle,
+
+ double *sunVarGeom_solarAltitude,
+ double *sunVarGeom_sinSolarAltitude,
+ double *sunVarGeom_tanSolarAltitude,
+ double *sunVarGeom_solarAzimuth,
+ double *sunVarGeom_sunAzimuthAngle,
+ double *sunVarGeom_stepsinangle,
+ double *sunVarGeom_stepcosangle,
+
+ double *gridGeom_stepxy,
+
+ double latitude);
+void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
+ struct GridGeometry *gridGeom);
+void com_par_const_parallel(double longitTime,
+ double *sunGeom_lum_C11,
+ double *sunGeom_lum_C13,
+ double *sunGeom_lum_C22,
+ double *sunGeom_lum_C31,
+ double *sunGeom_lum_C33,
+ double *sunGeom_sunrise_time,
+ double *sunGeom_sunset_time,
+ double *sunGeom_timeAngle,
+ double *sunGeom_sindecl,
+ double *sunGeom_cosdecl,
+
+ double *gridGeom_sinlat,
+ double *gridGeom_coslat);
+double lumcline2(struct SunGeometryConstDay *sungeom,
+ struct SunGeometryVarDay *sunVarGeom,
+ struct SunGeometryVarSlope *sunSlopeGeom,
+ struct GridGeometry *gridGeom,
+ unsigned char *horizonpointer);
+double lumcline2_parallel(double *sunGeom_timeAngle,
+
+ int *sunVarGeom_isShadow,
+ double *sunVarGeom_solarAltitude,
+ double *sunVarGeom_sunAzimuthAngle,
+ double *sunVarGeom_z_orig,
+ double *sunVarGeom_zmax,
+ double *sunVarGeom_zp,
+ double *sunVarGeom_tanSolarAltitude,
+ double *sunVarGeom_stepsinangle,
+ double *sunVarGeom_stepcosangle,
+
+ double *sunSlopeGeom_longit_l,
+ double *sunSlopeGeom_lum_C31_l,
+ double *sunSlopeGeom_lum_C33_l,
+
+ double *gridGeom_xx0,
+ double *gridGeom_yy0,
+ double *gridGeom_xg0,
+ double *gridGeom_yg0,
+ double *gridGeom_stepx,
+ double *gridGeom_stepy,
+ double *gridGeom_deltx,
+ double *gridGeom_delty,
+
+ 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);
Added: grass-addons/grass7/raster/r.sun.mp/main.c
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/main.c (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/main.c 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,2332 @@
+
+/*******************************************************************************
+r.sun.mp: 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
+(c) 2003-2013 by The GRASS Development Team
+*******************************************************************************/
+/*
+ * 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.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/*v. 2.0 July 2002, NULL data handling, JH */
+/*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
+/*v. 3.0 February 2006, several changes (shadowing algorithm, earth's curvature JH */
+/*v. 4.0 December 2016, OpenMP support Stanislav Zubal, Michal Lacko */
+
+
+#define PARALLEL
+
+#ifdef PARALLEL
+#include <omp.h>
+#define NUM_THREADS "4"
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#ifdef USE_OPENCL
+ #ifdef __APPLE__
+ #include <OpenCL/opencl.h>
+ #else
+ #include <CL/cl.h>
+ #endif
+#endif
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+#include "sunradstruct.h"
+#include "local_proto.h"
+#include "rsunglobals.h"
+
+/* default values */
+#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;
+
+const char *elevin;
+const char *aspin;
+const char *slopein;
+const char *civiltime = NULL;
+const char *linkein = NULL;
+const char *albedo = NULL;
+const char *latin = NULL;
+const char *coefbh = NULL;
+const char *coefdh = NULL;
+const char *incidout = NULL;
+const char *longin = NULL;
+const char *horizon = NULL;
+const char *beam_rad = NULL;
+const char *insol_time = NULL;
+const char *diff_rad = NULL;
+const char *refl_rad = NULL;
+const char *glob_rad = NULL;
+const char *mapset = NULL;
+const char *per;
+
+size_t decimals;
+char *str_step;
+
+
+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);
+
+#ifdef PARALLEL
+
+ int threads;
+
+time_t timer;
+
+void printTimeDiff(const char* message)
+{
+ double seconds;
+ time_t timer2;
+ time(&timer2); //or timer2 = time(NULL);
+ seconds = difftime(timer2,timer);
+ timer = timer2;
+/* G_message(_("TIME_DIFF: <%.0lf>\tmessage <%s>\n"),
+ seconds,
+ message);*/
+}
+
+
+#endif
+
+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 *Pbeam_e,
+ double *Pdiff_e,
+ double *Prefl_e,
+ double *Pinsol_t);
+
+
+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,
+ **latitudeArray = NULL, **longitudeArray, **cbhr = NULL, **cdhr = NULL;
+double op, dp;
+double invstepx, invstepy;
+double sunrise_min = 24., sunrise_max = 0., sunset_min = 24., sunset_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 linke_max = 0., linke_min = 100., albedo_max = 0., albedo_min = 1.0,
+ lat_max = -90., lat_min = 90.;
+double offsetx = 0.5, offsety = 0.5;
+char *ttime;
+
+/*
+ * 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;
+
+/* 1852m/nm * 60nm/degree = 111120 m/deg */
+#define DEGREEINMETERS 111120.
+
+int ll_correction = FALSE;
+double coslatsq;
+
+/* why not use G_distance() here which switches to geodesic/great
+ circle distace as needed? */
+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;
+#ifdef PARALLEL
+ timer = time(NULL);
+ printTimeDiff("Start");
+
+ //int threads;
+
+
+ struct GModule *module;
+ struct
+ {
+ struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
+ *lin, *albedo, *longin, *alb, *latin, *coefbh, *coefdh,
+ *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
+ *glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
+ *horizonstep, *numPartitions, *civilTime, *threads;
+ }
+ parm;
+#endif
+#ifndef PARALLEL
+ struct GModule *module;
+ struct
+ {
+ struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
+ *lin, *albedo, *longin, *alb, *latin, *coefbh, *coefdh,
+ *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
+ *glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
+ *horizonstep, *numPartitions, *civilTime;
+ }
+ parm;
+#endif
+
+ struct
+ {
+ struct Flag *noshade, *saveMemory;
+ }
+ flag;
+
+ struct GridGeometry gridGeom;
+
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("solar"));
+ G_add_keyword(_("sun energy"));
+ G_add_keyword(_("shadow"));
+ module->label = _("Solar irradiance and irradiation model.");
+ 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 the map history 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.");
+
+ parm.elevin = G_define_option();
+ parm.elevin->key = "elevation";
+ 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");
+
+ parm.aspin = G_define_option();
+ parm.aspin->key = "aspect";
+ 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");
+
+ parm.aspect = G_define_option();
+ parm.aspect->key = "aspect_value";
+ 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");
+
+ parm.slopein = G_define_option();
+ parm.slopein->key = "slope";
+ 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");
+
+ parm.slope = G_define_option();
+ parm.slope->key = "slope_value";
+ 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");
+
+ parm.linkein = G_define_option();
+ parm.linkein->key = "linke";
+ 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");
+
+ parm.lin = G_define_option();
+ parm.lin->key = "linke_value";
+ 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");
+
+ 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");
+
+ parm.alb = G_define_option();
+ parm.alb->key = "albedo_value";
+ 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");
+
+ parm.latin = G_define_option();
+ parm.latin->key = "lat";
+ parm.latin->type = TYPE_STRING;
+ parm.latin->required = NO;
+ parm.latin->gisprompt = "old,cell,raster";
+ parm.latin->description =
+ _("Name of input raster map containing latitudes [decimal degrees]");
+ parm.latin->guisection = _("Input");
+
+ parm.longin = G_define_option();
+ parm.longin->key = "long";
+ parm.longin->type = TYPE_STRING;
+ parm.longin->required = NO;
+ parm.longin->gisprompt = "old,cell,raster";
+ parm.longin->description =
+ _("Name of input raster map containing longitudes [decimal degrees]");
+ parm.longin->guisection = _("Input");
+
+ parm.coefbh = G_define_option();
+ parm.coefbh->key = "coeff_bh";
+ 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 (thick cloud) input raster map [0-1]");
+ parm.coefbh->guisection = _("Input");
+
+ parm.coefdh = G_define_option();
+ parm.coefdh->key = "coeff_dh";
+ 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 (haze) input raster map [0-1]");
+ parm.coefdh->guisection = _("Input");
+
+ parm.horizon = G_define_standard_option(G_OPT_R_BASENAME_INPUT);
+ parm.horizon->key = "horizon_basename";
+ parm.horizon->required = NO;
+ parm.horizon->gisprompt = "old,cell,raster";
+ parm.horizon->description = _("The horizon information input map basename");
+ parm.horizon->guisection = _("Input");
+
+ parm.horizonstep = G_define_option();
+ parm.horizonstep->key = "horizon_step";
+ parm.horizonstep->type = TYPE_DOUBLE;
+ parm.horizonstep->required = NO;
+ parm.horizonstep->description =
+ _("Angle step size for multidirectional horizon [degrees]");
+ parm.horizonstep->guisection = _("Input");
+
+ parm.incidout = G_define_option();
+ parm.incidout->key = "incidout";
+ parm.incidout->type = TYPE_STRING;
+ parm.incidout->required = NO;
+ parm.incidout->gisprompt = "new,cell,raster";
+ parm.incidout->description =
+ _("Output incidence angle raster map (mode 1 only)");
+ parm.incidout->guisection = _("Output");
+
+ 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->gisprompt = "new,cell,raster";
+ 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");
+
+ 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->gisprompt = "new,cell,raster";
+ 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");
+
+ 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->gisprompt = "new,cell,raster";
+ 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");
+
+ 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->gisprompt = "new,cell,raster";
+ 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");
+
+ 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->gisprompt = "new,cell,raster";
+ parm.insol_time->description =
+ _("Output insolation time raster map [h] (mode 2 only)");
+ parm.insol_time->guisection = _("Output");
+
+ 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.day->options = "1-365";
+ parm.day->guisection = _("Time");
+
+ 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.step->guisection = _("Time");
+
+ parm.declin = G_define_option();
+ parm.declin->key = "declination";
+ 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.ltime->options = "0-24";
+ parm.ltime->guisection = _("Time");
+
+#ifdef PARALLEL
+ parm.threads = G_define_option();
+ parm.threads->key = "threads";
+ parm.threads->type = TYPE_INTEGER;
+ parm.threads->answer = NUM_THREADS;
+ parm.threads->required = NO;
+ parm.threads->description =
+ _("Number of threads which will be used for parallel compute");
+ parm.threads->guisection = _("Parameters");
+#endif
+ /*
+ * parm.startTime = G_define_option();
+ * parm.startTime->key = "start_time";
+ * 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 = "end_time";
+ * 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 = "distance_step";
+ 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 = "npartitions";
+ 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 = "civil_time";
+ 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");
+ parm.civilTime->guisection = _("Time");
+
+ flag.noshade = G_define_flag();
+ flag.noshade->key = 'p';
+ flag.noshade->description =
+ _("Do not 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);
+
+
+ G_get_set_window(&cellhd);
+
+ 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);
+
+ setUseShadow(!flag.noshade->answer);
+
+ saveMemory = flag.saveMemory->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;
+ longin = parm.longin->answer;
+
+ civiltime = parm.civilTime->answer;
+#ifdef PARALLEL
+ sscanf(parm.threads->answer, "%d", &threads);
+ if (threads < 1)
+ {
+ G_warning(_("<%d> is not valid number of threads. Number of threads will be set on <%d>"),
+ threads, abs(threads));
+ threads = abs(threads);
+ }
+ omp_set_num_threads(threads);
+ G_message(_("Number of threads <%d>"),threads);
+ int nthreads, tid;
+#endif
+#ifdef PARALLEL
+ printTimeDiff("M1");
+#endif
+ if (civiltime != NULL) {
+ setUseCivilTime(TRUE);
+
+ if (longin == NULL)
+ G_fatal_error(
+ _("You must give the longitude raster if you use civil time"));
+
+ if (sscanf(parm.civilTime->answer, "%lf", &civilTime) != 1)
+ G_fatal_error(_("Error reading civil time zone value"));
+
+ if (civilTime < -24. || civilTime > 24.)
+ G_fatal_error(_("Invalid civil time zone value"));
+
+ /* Normalize if somebody should be weird enough to give more than +- 12
+ * hours offset. */
+ if (civilTime < -12.)
+ civilTime += 24.;
+ else if (civilTime > 12.)
+ civilTime -= 24;
+ }
+ else {
+ setUseCivilTime(FALSE);
+ }
+
+ 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);
+
+ if (sscanf(parm.step->answer, "%lf", &step) != 1)
+ G_fatal_error(_("Error reading time step size"));
+ if (step <= 0.0 || step > 24.0)
+ G_fatal_error(_("Invalid time step size"));
+
+ if (parm.horizonstep->answer != NULL) {
+ if (sscanf(parm.horizonstep->answer, "%lf", &horizonStep) != 1)
+ G_fatal_error(_("Error reading horizon step size"));
+ str_step = parm.horizonstep->answer;
+ 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."));
+ }
+
+ ttime = 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 requires time parameter to be set"));
+
+ G_message(_("Mode 2: integrated daily irradiation for a given day of the year"));
+ }
+
+ /*
+ * 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);
+ G_message(_("Using Linke constant: %lf"), singleLinke);
+ } else {
+ G_message(_("Using Linke map <%s>"), parm.linkein->answer);
+ }
+
+ if (parm.albedo->answer == NULL){
+ sscanf(parm.alb->answer, "%lf", &singleAlbedo);
+ G_message(_("Using albedo constant: %lf"), singleAlbedo);
+ } else {
+ G_message(_("Using albedo map <%s>"), parm.albedo->answer);
+ }
+
+ if (parm.slopein->answer == NULL){
+ sscanf(parm.slope->answer, "%lf", &singleSlope);
+ G_message(_("Using slope constant: %lf"), singleSlope);
+ singleSlope *= deg2rad;
+ } else {
+ G_message(_("Using slope map <%s>"), parm.slopein->answer);
+ }
+
+ if (parm.aspin->answer == NULL){
+ sscanf(parm.aspect->answer, "%lf", &singleAspect);
+ G_message(_("Using aspect constant: %lf"), singleAspect);
+ singleAspect *= deg2rad;
+ } else {
+ G_message(_("Using aspect map <%s>"), parm.aspin->answer);
+ }
+
+ 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 (ttime != 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 (ttime != 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 */
+ /* we need to do this even if latin= and longin= were given because
+ rsunlib's com_par() needs iproj and oproj */
+ if ((G_projection() != PROJECTION_LL)) {
+ 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"));
+
+ 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"));
+ }
+
+ if ((latin != NULL || longin != NULL) && (G_projection() == PROJECTION_LL))
+ G_warning(_("latin and longin raster maps have no effect when in a Lat/Lon location"));
+ /* true about longin= when civiltime is used? */
+ /* civiltime needs longin= but not latin= for non-LL projections -
+ better would be it just use pj_proj() if it needs those?? */
+ if (latin != NULL && longin == NULL)
+ G_fatal_error(_("Both latin and longin raster maps must be given, or neither"));
+
+/**********end of parser - ******************************/
+
+
+ if ((G_projection() == PROJECTION_LL))
+ ll_correction = TRUE;
+
+ G_debug(3, "calculate() starts...");
+ calculate(singleSlope, singleAspect, singleAlbedo, singleLinke, gridGeom);
+ G_debug(3, "OUTGR() starts...");
+ OUTGR();
+#ifdef PARALLEL
+ printTimeDiff("M4");
+#endif
+ exit(EXIT_SUCCESS);
+}
+
+
+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;
+ char formatString[10];
+ double angle_deg = 0.;
+
+ finalRow = m - offset - m / numPartitions;
+ if (finalRow < 0) {
+ finalRow = 0;
+ }
+
+ numRows = m / numPartitions;
+
+ cell1 = Rast_allocate_f_buf();
+
+ if (z == NULL) {
+ if (!(z = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //z = (float **)G_malloc(sizeof(float *) * (numRows));
+
+
+ //for (l = 0; l < numRows; l++) {
+ // z[l] = (float *)G_malloc(sizeof(float) * (n));
+ //}
+ }
+
+ fd1 = Rast_open_old(elevin, "");
+
+ if (slopein != NULL) {
+ cell3 = Rast_allocate_f_buf();
+ if (s == NULL) {
+ if (!(s = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //s = (float **)G_malloc(sizeof(float *) * (numRows));
+
+ //for (l = 0; l < numRows; l++) {
+ // s[l] = (float *)G_malloc(sizeof(float) * (n));
+ //}
+ }
+ fd3 = Rast_open_old(slopein, "");
+
+ }
+
+ if (aspin != NULL) {
+ cell2 = Rast_allocate_f_buf();
+
+ if (o == NULL) {
+ if (!(o = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //o = (float **)G_malloc(sizeof(float *) * (numRows));
+
+ //for (l = 0; l < numRows; l++) {
+ // o[l] = (float *)G_malloc(sizeof(float) * (n));
+ //}
+ }
+ fd2 = Rast_open_old(aspin, "");
+
+ }
+
+ if (linkein != NULL) {
+ cell4 = Rast_allocate_f_buf();
+ if (li == NULL) {
+ if (!(li = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //li = (float **)G_malloc(sizeof(float *) * (numRows));
+ //for (l = 0; l < numRows; l++)
+ // li[l] = (float *)G_malloc(sizeof(float) * (n));
+
+ }
+ fd4 = Rast_open_old(linkein, "");
+ }
+
+ if (albedo != NULL) {
+ cell5 = Rast_allocate_f_buf();
+ if (a == NULL) {
+ if (!(a = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //a = (float **)G_malloc(sizeof(float *) * (numRows));
+ //for (l = 0; l < numRows; l++)
+ // a[l] = (float *)G_malloc(sizeof(float) * (n));
+ }
+ fd5 = Rast_open_old(albedo, "");
+ }
+
+ if (latin != NULL) {
+ cell6 = Rast_allocate_f_buf();
+ if (latitudeArray == NULL) {
+ if (!(latitudeArray = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //latitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
+ //for (l = 0; l < numRows; l++)
+ // latitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
+ }
+ fd6 = Rast_open_old(latin, "");
+ }
+
+ if (longin != NULL) {
+ cell7 = Rast_allocate_f_buf();
+ if (latitudeArray == NULL) {
+ if (!(longitudeArray = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //longitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
+ //for (l = 0; l < numRows; l++)
+ // longitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
+ }
+
+ fd7 = Rast_open_old(longin, "");
+ }
+
+ if (coefbh != NULL) {
+ rast1 = Rast_allocate_f_buf();
+ if (cbhr == NULL) {
+ if (!(cbhr = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //cbhr = (float **)G_malloc(sizeof(float *) * (numRows));
+ //for (l = 0; l < numRows; l++)
+ // cbhr[l] = (float *)G_malloc(sizeof(float) * (n));
+ }
+ fr1 = Rast_open_old(coefbh, "");
+ }
+
+ if (coefdh != NULL) {
+ rast2 = Rast_allocate_f_buf();
+ if (cdhr == NULL) {
+ if (!(cdhr = G_alloc_fmatrix(numRows, n)))
+ {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ //cdhr = (float **)G_malloc(sizeof(float *) * (numRows));
+ //for (l = 0; l < numRows; l++)
+ // cdhr[l] = (float *)G_malloc(sizeof(float) * (n));
+ }
+ fr2 = Rast_open_old(coefdh, "");
+ }
+
+ if (useHorizonData()) {
+ if (horizonarray == NULL) {
+ horizonarray =
+ (unsigned char *)G_calloc(arrayNumInt *
+ numRows * n, sizeof(char));
+
+ horizonbuf = (FCELL **) G_calloc(arrayNumInt, sizeof(FCELL *) );
+ fd_shad = (int *)G_calloc(arrayNumInt, sizeof(int));
+// 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(ttime != NULL)
+ * {
+ *
+ * horizonbuf[0]=Rast_allocate_f_buf();
+ * sprintf(shad_filename, "%s_%02d", horizon, arrayNumInt);
+ * if((mapset=G_find_raster2(shad_filename,""))==NULL)
+ * G_message("Horizon file no. %d not found\n", arrayNumInt);
+ *
+ * fd_shad[0] = Rast_open_old(shad_filename,mapset);
+ * }
+ * else
+ * {
+ */
+ decimals = G_get_num_decimals(str_step);
+ angle_deg = 0;
+ for (i = 0; i < arrayNumInt; i++) {
+ horizonbuf[i] = Rast_allocate_f_buf();
+ shad_filename = G_generate_basename(horizon, angle_deg,
+ 3, decimals);
+ fd_shad[i] = Rast_open_old(shad_filename, "");
+ angle_deg += horizonStep;
+ G_free(shad_filename);
+ }
+ /*
+ numDigits = (int)(log10(1. * arrayNumInt)) + 1;
+ sprintf(formatString, "%%s_%%0%dd", numDigits);
+ for (i = 0; i < arrayNumInt; i++) {
+ horizonbuf[i] = Rast_allocate_f_buf();
+ sprintf(shad_filename, formatString, horizon, i);
+ fd_shad[i] = Rast_open_old(shad_filename, "");
+ } */
+ }
+ /*
+ * }
+ */
+
+ 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;
+ Rast_get_f_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--) {
+ Rast_get_f_row(fd1, cell1, row);
+ if (aspin != NULL)
+ Rast_get_f_row(fd2, cell2, row);
+ if (slopein != NULL)
+ Rast_get_f_row(fd3, cell3, row);
+ if (linkein != NULL)
+ Rast_get_f_row(fd4, cell4, row);
+ if (albedo != NULL)
+ Rast_get_f_row(fd5, cell5, row);
+ if (latin != NULL)
+ Rast_get_f_row(fd6, cell6, row);
+ if (longin != NULL)
+ Rast_get_f_row(fd7, cell7, row);
+ if (coefbh != NULL)
+ Rast_get_f_row(fr1, rast1, row);
+ if (coefdh != NULL)
+ Rast_get_f_row(fr2, rast2, row);
+
+ row_rev = m - row - 1;
+ rowrevoffset = row_rev - offset;
+
+ for (j = 0; j < n; j++) {
+ if (!Rast_is_f_null_value(cell1 + j))
+ z[rowrevoffset][j] = (float)cell1[j];
+ else
+ z[rowrevoffset][j] = UNDEFZ;
+
+ if (aspin != NULL) {
+ if (!Rast_is_f_null_value(cell2 + j))
+ o[rowrevoffset][j] = (float)cell2[j];
+ else
+ o[rowrevoffset][j] = UNDEFZ;
+ }
+ if (slopein != NULL) {
+ if (!Rast_is_f_null_value(cell3 + j))
+ s[rowrevoffset][j] = (float)cell3[j];
+ else
+ s[rowrevoffset][j] = UNDEFZ;
+ }
+
+ if (linkein != NULL) {
+ if (!Rast_is_f_null_value(cell4 + j))
+ li[rowrevoffset][j] = (float)cell4[j];
+ else
+ li[rowrevoffset][j] = UNDEFZ;
+ }
+
+ if (albedo != NULL) {
+ if (!Rast_is_f_null_value(cell5 + j))
+ a[rowrevoffset][j] = (float)cell5[j];
+ else
+ a[rowrevoffset][j] = UNDEFZ;
+ }
+
+ if (latin != NULL) {
+ if (!Rast_is_f_null_value(cell6 + j))
+ latitudeArray[rowrevoffset][j] = (float)cell6[j];
+ else
+ latitudeArray[rowrevoffset][j] = UNDEFZ;
+ }
+
+ if (longin != NULL) {
+ if (!Rast_is_f_null_value(cell7 + j))
+ longitudeArray[rowrevoffset][j] = (float)cell7[j];
+ else
+ longitudeArray[rowrevoffset][j] = UNDEFZ;
+ }
+
+ if (coefbh != NULL) {
+ if (!Rast_is_f_null_value(rast1 + j))
+ cbhr[rowrevoffset][j] = (float)rast1[j];
+ else
+ cbhr[rowrevoffset][j] = UNDEFZ;
+ }
+
+ if (coefdh != NULL) {
+ if (!Rast_is_f_null_value(rast2 + j))
+ cdhr[rowrevoffset][j] = (float)rast2[j];
+ else
+ cdhr[rowrevoffset][j] = UNDEFZ;
+ }
+ }
+ }
+
+ Rast_close(fd1);
+ G_free(cell1);
+
+ if (aspin != NULL) {
+ G_free(cell2);
+ Rast_close(fd2);
+ }
+ if (slopein != NULL) {
+ G_free(cell3);
+ Rast_close(fd3);
+ }
+ if (linkein != NULL) {
+ G_free(cell4);
+ Rast_close(fd4);
+ }
+ if (albedo != NULL) {
+ G_free(cell5);
+ Rast_close(fd5);
+ }
+ if (latin != NULL) {
+ G_free(cell6);
+ Rast_close(fd6);
+ }
+ if (longin != NULL) {
+ G_free(cell7);
+ Rast_close(fd7);
+ }
+ if (coefbh != NULL) {
+ G_free(rast1);
+ Rast_close(fr1);
+ }
+ if (coefdh != NULL) {
+ G_free(rast2);
+ Rast_close(fr2);
+ }
+
+
+ if (useHorizonData()) {
+ for (i = 0; i < arrayNumInt; i++) {
+ Rast_close(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 && latitudeArray[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 = Rast_allocate_f_buf();
+ fd7 = Rast_open_fp_new(incidout);
+ }
+
+ if (beam_rad != NULL) {
+ cell8 = Rast_allocate_f_buf();
+ fd8 = Rast_open_fp_new(beam_rad);
+ }
+
+ if (insol_time != NULL) {
+ cell11 = Rast_allocate_f_buf();
+ fd11 = Rast_open_fp_new(insol_time);
+ }
+
+ if (diff_rad != NULL) {
+ cell9 = Rast_allocate_f_buf();
+ fd9 = Rast_open_fp_new(diff_rad);
+ }
+
+ if (refl_rad != NULL) {
+ cell10 = Rast_allocate_f_buf();
+ fd10 = Rast_open_fp_new(refl_rad);
+ }
+
+ if (glob_rad != NULL) {
+ cell12 = Rast_allocate_f_buf();
+ fd12 = Rast_open_fp_new(glob_rad);
+ }
+
+ if (m != Rast_window_rows())
+ G_fatal_error("OOPS: rows changed from %d to %d", m, Rast_window_rows());
+ if (n != Rast_window_cols())
+ G_fatal_error("OOPS: cols changed from %d to %d", n, Rast_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)
+ Rast_set_f_null_value(cell7 + j, 1);
+ else
+ cell7[j] = (FCELL) lumcl[i][j];
+ }
+ Rast_put_f_row(fd7, cell7);
+ }
+
+ if (beam_rad != NULL) {
+ for (j = 0; j < n; j++) {
+ if (beam[i][j] == UNDEFZ)
+ Rast_set_f_null_value(cell8 + j, 1);
+ else
+ cell8[j] = (FCELL) beam[i][j];
+ }
+ Rast_put_f_row(fd8, cell8);
+ }
+
+ if (glob_rad != NULL) {
+ for (j = 0; j < n; j++) {
+ if (globrad[i][j] == UNDEFZ)
+ Rast_set_f_null_value(cell12 + j, 1);
+ else
+ cell12[j] = (FCELL) globrad[i][j];
+ }
+ Rast_put_f_row(fd12, cell12);
+ }
+
+
+ if (insol_time != NULL) {
+ for (j = 0; j < n; j++) {
+ if (insol[i][j] == UNDEFZ)
+ Rast_set_f_null_value(cell11 + j, 1);
+ else
+ cell11[j] = (FCELL) insol[i][j];
+ }
+ Rast_put_f_row(fd11, cell11);
+ }
+
+
+ if (diff_rad != NULL) {
+ for (j = 0; j < n; j++) {
+ if (diff[i][j] == UNDEFZ)
+ Rast_set_f_null_value(cell9 + j, 1);
+ else
+ cell9[j] = (FCELL) diff[i][j];
+ }
+ Rast_put_f_row(fd9, cell9);
+ }
+
+ if (refl_rad != NULL) {
+ for (j = 0; j < n; j++) {
+ if (refl[i][j] == UNDEFZ)
+ Rast_set_f_null_value(cell10 + j, 1);
+ else
+ cell10[j] = (FCELL) refl[i][j];
+ }
+ Rast_put_f_row(fd10, cell10);
+ }
+
+ }
+
+ if (incidout != NULL) {
+ Rast_close(fd7);
+ Rast_write_history(incidout, &hist);
+ }
+ if (beam_rad != NULL) {
+ Rast_close(fd8);
+ Rast_write_history(beam_rad, &hist);
+ }
+ if (diff_rad != NULL) {
+ Rast_close(fd9);
+ Rast_write_history(diff_rad, &hist);
+ }
+ if (refl_rad != NULL) {
+ Rast_close(fd10);
+ Rast_write_history(refl_rad, &hist);
+ }
+ if (insol_time != NULL) {
+ Rast_close(fd11);
+ Rast_write_history(insol_time, &hist);
+ }
+ if (glob_rad != NULL) {
+ Rast_close(fd12);
+ Rast_write_history(glob_rad, &hist);
+ }
+
+ return 1;
+}
+
+
+/**********************************************************/
+
+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 *Pbeam_e,
+ double *Pdiff_e,
+ double *Prefl_e,
+ double *Pinsol_t)
+{
+
+ 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 (ttime != 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 */
+ *Pbeam_e += ra;
+ }
+ else {
+ *Pbeam_e = 0.;
+ bh = 0.;
+ }
+
+ if ((diff_rad != NULL) || (glob_rad != NULL)) {
+ dra = drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar); /* diffuse rad. */
+ *Pdiff_e += dra;
+ }
+ if ((refl_rad != NULL) || (glob_rad != NULL)) {
+ if ((diff_rad == NULL) && (glob_rad == NULL))
+ drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
+ *Prefl_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.)) {
+ *Pinsol_t += dfr;
+ ra = brad(s0, &bh, sunVarGeom, sunSlopeGeom, sunRadVar);
+ *Pbeam_e += dfr * ra;
+ ra = 0.;
+ }
+ else {
+ bh = 0.;
+ }
+ if ((diff_rad != NULL) || (glob_rad != NULL)) {
+ dra =
+ drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom,
+ sunRadVar);
+ *Pdiff_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);
+ }
+ *Prefl_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;
+// }
+// }
+ if (incidout != NULL) {
+ lumcl = (float **)G_calloc((m), sizeof(float *));
+ for (l = 0; l < m; l++) {
+ lumcl[l] = (float *)G_calloc((n), sizeof(float *));
+ }
+ for (j = 0; j < m; j++) {
+ for (i = 0; i < n; i++)
+ lumcl[j][i] = UNDEFZ;
+ }
+ }
+
+ if (beam_rad != NULL) {
+ beam = (float **)G_calloc((m), sizeof(float *));
+ for (l = 0; l < m; l++) {
+ beam[l] = (float *)G_calloc((n), sizeof(float *));
+ }
+
+ for (j = 0; j < m; j++) {
+ for (i = 0; i < n; i++)
+ beam[j][i] = UNDEFZ;
+ }
+ }
+
+ if (insol_time != NULL) {
+ insol = (float **)G_calloc((m), sizeof(float *));
+ for (l = 0; l < m; l++) {
+ insol[l] = (float *)G_calloc((n), sizeof(float *));
+ }
+
+ for (j = 0; j < m; j++) {
+ for (i = 0; i < n; i++)
+ insol[j][i] = UNDEFZ;
+ }
+ }
+
+ if (diff_rad != NULL) {
+ diff = (float **)G_calloc((m), sizeof(float *));
+ for (l = 0; l < m; l++) {
+ diff[l] = (float *)G_calloc((n), sizeof(float *));
+ }
+
+ for (j = 0; j < m; j++) {
+ for (i = 0; i < n; i++)
+ diff[j][i] = UNDEFZ;
+ }
+ }
+
+ if (refl_rad != NULL) {
+ refl = (float **)G_calloc((m), sizeof(float *));
+ for (l = 0; l < m; l++) {
+ refl[l] = (float *)G_calloc((n), sizeof(float *));
+ }
+
+ for (j = 0; j < m; j++) {
+ for (i = 0; i < n; i++)
+ refl[j][i] = UNDEFZ;
+ }
+ }
+
+ if (glob_rad != NULL) {
+ globrad = (float **)G_calloc((m), sizeof(float *));
+ for (l = 0; l < m; l++) {
+ globrad[l] = (float *)G_calloc((n), sizeof(float *));
+ }
+
+ 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.);
+ }
+#ifdef PARALLEL
+ printTimeDiff("M2");
+#endif
+ int shadowoffset_base = shadowoffset;
+ 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;
+ shadowoffset_base = (j % (numRows)) * n * arrayNumInt;
+ #pragma omp parallel firstprivate(q1,tan_lam_l,z1,i,shadowoffset,longitTime,coslat,coslatsq,func,latitude,longitude,sin_phi_l,latid_l,sin_u,cos_u,sin_v,cos_v,lum, gridGeom,sunGeom,sunVarGeom,sunSlopeGeom,sunRadVar, elevin,aspin,slopein,civiltime,linkein,albedo,latin,coefbh,coefdh,incidout,longin,horizon,beam_rad,insol_time,diff_rad,refl_rad,glob_rad,mapset,per,decimals,str_step )
+{
+ #pragma omp for schedule(dynamic)
+ for (i = 0; i < n; i++) {
+ shadowoffset = shadowoffset_base + (arrayNumInt * i);
+ //G_message("\n tid: %d", omp_get_thread_num());
+ if (useCivilTime()) {
+ /* sun travels 15deg per hour, so 1 TZ every 15 deg and 15 TZs * 24hrs = 360deg */
+ longitTime = -longitudeArray[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;
+if (i >= n || arrayOffset >= m)
+{
+ G_message("i: %d, n: %d, arrayOffset: %d, m: %d", i,n,arrayOffset,m);
+}
+
+ 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];
+ linke_max = AMAX1(linke_max, sunRadVar.linke);
+ linke_min = AMIN1(linke_min, sunRadVar.linke);
+ }
+ if (albedo != NULL) {
+ sunRadVar.alb = a[arrayOffset][i];
+ albedo_max = AMAX1(albedo_max, sunRadVar.alb);
+ albedo_min = AMIN1(albedo_min, sunRadVar.alb);
+ }
+ if (latin != NULL) {
+ latitude = latitudeArray[arrayOffset][i];
+ lat_max = AMAX1(lat_max, latitude);
+ lat_min = AMIN1(lat_min, latitude);
+ latitude *= deg2rad;
+ }
+ if (longin != NULL) {
+ longitude = longitudeArray[arrayOffset][i];
+ /* lon_max = AMAX1(lon_max, longitude); */
+ /* lon_min = AMIN1(lon_min, longitude); */
+ longitude *= deg2rad;
+ }
+
+ if ((G_projection() != PROJECTION_LL)) {
+
+ if (latin == NULL || longin == NULL) {
+ /* if either is missing we have to calc both from current projection */
+ longitude = gridGeom.xp;
+ latitude = gridGeom.yp;
+
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+ G_fatal_error("Error in pj_do_proj");
+ }
+
+ lat_max = AMAX1(lat_max, latitude);
+ lat_min = AMIN1(lat_min, latitude);
+ latitude *= deg2rad;
+ longitude *= deg2rad;
+ }
+ }
+ else { /* ll projection */
+ latitude = gridGeom.yp;
+ longitude = gridGeom.xp;
+ lat_max = AMAX1(lat_max, latitude);
+ lat_min = AMIN1(lat_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 (ttime != 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);
+ sunrise_min = AMIN1(sunrise_min, sunGeom.sunrise_time);
+ sunrise_max = AMAX1(sunrise_max, sunGeom.sunrise_time);
+ sunset_min = AMIN1(sunset_min, sunGeom.sunset_time);
+ sunset_max = AMAX1(sunset_max, sunGeom.sunset_time);
+
+ }
+
+ if (incidout != NULL) {
+ com_par(&sunGeom, &sunVarGeom, &gridGeom, latitude,
+ longitude);
+ lum =
+ lumcline2(&sunGeom, &sunVarGeom, &sunSlopeGeom,
+ &gridGeom, horizonarray + shadowoffset);
+ if (lum > 0.) {
+ lum = rad2deg * asin(lum);
+ lumcl[j][i] = (float)lum;
+ }
+ else lumcl[j][i] = UNDEFZ;
+ }
+
+ double Pbeam_e = 0.;
+ double Pdiff_e = 0.;
+ double Prefl_e = 0.;
+ double Pinsol_t = 0.;
+
+ if (someRadiation) {
+ joules2(&sunGeom, &sunVarGeom, &sunSlopeGeom, &sunRadVar,
+ &gridGeom, horizonarray + shadowoffset, latitude,
+ longitude,
+ &Pbeam_e,
+ &Pdiff_e,
+ &Prefl_e,
+ &Pinsol_t
+ );
+ if (beam_rad != NULL)
+ beam[j][i] = (float)Pbeam_e;
+ if (insol_time != NULL)
+ insol[j][i] = (float)Pinsol_t;
+ /* G_debug(3,"\n %f",insol[j][i]); */
+ if (diff_rad != NULL)
+ diff[j][i] = (float)Pdiff_e;
+ if (refl_rad != NULL)
+ refl[j][i] = (float)Prefl_e;
+ if (glob_rad != NULL)
+ globrad[j][i] = (float)(Pbeam_e + Pdiff_e + Prefl_e);
+ }
+
+ } /* undefs */
+ //shadowoffset += arrayNumInt;
+ }
+ }
+ arrayOffset++;
+}
+#ifdef PARALLEL
+ printTimeDiff("M3");
+#endif
+ /* re-use &hist, but try all to initiate it for any case */
+ /* note this will result in incorrect map titles */
+ if (incidout != NULL) {
+ Rast_short_history(incidout, "raster", &hist);
+ }
+ else if (beam_rad != NULL) {
+ Rast_short_history(beam_rad, "raster", &hist);
+ }
+ else if (diff_rad != NULL) {
+ Rast_short_history(diff_rad, "raster", &hist);
+ }
+ else if (refl_rad != NULL) {
+ Rast_short_history(refl_rad, "raster", &hist);
+ }
+ else if (insol_time != NULL) {
+ Rast_short_history(insol_time, "raster", &hist);
+ }
+ else if (glob_rad != NULL) {
+ Rast_short_history(glob_rad, "raster", &hist);
+ }
+ else
+ G_fatal_error
+ ("Failed to init map history: no output maps requested!");
+
+ Rast_append_format_history(
+ &hist,
+ " ----------------------------------------------------------------");
+ Rast_append_format_history(
+ &hist,
+ " Day [1-365]: %d",
+ day);
+
+ if (ttime != NULL)
+ Rast_append_format_history(
+ &hist,
+ " Local (solar) time (decimal hr.): %.4f", timo);
+
+ Rast_append_format_history(
+ &hist,
+ " Solar constant (W/m^2): 1367");
+ Rast_append_format_history(
+ &hist,
+ " Extraterrestrial irradiance (W/m^2): %f",
+ sunRadVar.G_norm_extra);
+ Rast_append_format_history(
+ &hist,
+ " Declination (rad): %f", -declination);
+
+ Rast_append_format_history(
+ &hist,
+ " Latitude min-max(deg): %.4f - %.4f",
+ lat_min, lat_max);
+
+ if (ttime != NULL) {
+ Rast_append_format_history(
+ &hist,
+ " Sunrise time (hr.): %.2f",
+ sunGeom.sunrise_time);
+ Rast_append_format_history(
+ &hist,
+ " Sunset time (hr.): %.2f",
+ sunGeom.sunset_time);
+ Rast_append_format_history(
+ &hist,
+ " Daylight time (hr.): %.2f",
+ sunGeom.sunset_time - sunGeom.sunrise_time);
+ }
+ else {
+ Rast_append_format_history(
+ &hist,
+ " Sunrise time min-max (hr.): %.2f - %.2f",
+ sunrise_min, sunrise_max);
+ Rast_append_format_history(
+ &hist,
+ " Sunset time min-max (hr.): %.2f - %.2f",
+ sunset_min, sunset_max);
+ Rast_append_format_history(
+ &hist,
+ " Time step (hr.): %.4f", step);
+ }
+
+ if (incidout != NULL || ttime != NULL) {
+ Rast_append_format_history(
+ &hist,
+ " Solar altitude (deg): %.4f",
+ sunVarGeom.solarAltitude * rad2deg);
+ Rast_append_format_history(
+ &hist,
+ " Solar azimuth (deg): %.4f",
+ sunVarGeom.solarAzimuth * rad2deg);
+ }
+
+ if (linkein == NULL)
+ Rast_append_format_history(
+ &hist,
+ " Linke turbidity factor: %.1f",
+ sunRadVar.linke);
+ else
+ Rast_append_format_history(
+ &hist,
+ " Linke turbidity factor min-max: %.1f-%.1f",
+ linke_min, linke_max);
+
+ if (albedo == NULL)
+ Rast_append_format_history(
+ &hist,
+ " Ground albedo: %.3f",
+ sunRadVar.alb);
+ else
+ Rast_append_format_history(
+ &hist,
+ " Ground albedo min-max: %.3f-%.3f",
+ albedo_min, albedo_max);
+
+ Rast_append_format_history(
+ &hist,
+ " -----------------------------------------------------------------");
+
+ Rast_command_history(&hist);
+ /* don't call Rast_write_history() until after Rast_close() or it just gets overwritten */
+
+} /* End of ) function */
+
+
+
+double com_declin(int no_of_day)
+{
+ double d1, decl;
+
+ /* stretch day number in the following calculation for siderial effect? */
+ /* ? double siderial_day = no_of_day + ((no_of_day * 0.25) / 365.) ? */
+ /* or just change d1 to : d1 = pi2 * no_of_day / 365.0; ? */
+ 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;
+}
Added: grass-addons/grass7/raster/r.sun.mp/r.sun.mp.html
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/r.sun.mp.html (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/r.sun.mp.html 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,384 @@
+<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 incorporated by default. This can be done either internally by
+calculatoion of the shadowing effect directly from the digital elevation model
+or by specifying raster maps of the horizon height which is much faster. These
+horizon raster maps can be calculated 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>The model considers a shadowing effect of the local topography unless switched
+off with the <em>-p</em> flag.
+<b>r.sun</b> 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>
+The solar incidence angle raster map <i>incidout</i> is computed specifying
+elevation raster map <i>elevation</i>, aspect raster map <i>aspect</i>, slope
+steepness raster map <i>slope,</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 explicitly for large areas by input raster map <i>lat_in</i>
+ with interpolated latitude values. 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>
+The solar <i>declination</i> parameter 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>linke</i>, <i>lin</i>
+) and ground albedo (<i>albedo</i>, <i>alb</i>) can be set.
+<p>Besides clear-sky radiations, the user can compute a real-sky radiation (beam,
+diffuse) using <i>coeff_bh</i> and <i>coeff_dh</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
+provided as raster maps with spatial distribution of these coefficients separately
+for beam and diffuse radiation (see Suri and Hofierka, 2004, section 3.2).
+
+<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_in</i>,
+elevation <i>elevation</i>, slope <i>slope</i> and aspect <i>aspect</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.
+
+<p>
+The solar radiation maps for a given day are computed by integrating the
+relevant irradiance between sunrise and sunset times for that day. The
+user can set a finer or coarser time step used for all-day radiation
+calculations with the <i>step</i> option. The 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 (and more jagged) results. As the sun moves through approx.
+15° of the sky in an hour, the default <i>step</i> of half an hour will
+produce 7.5° steps in the data. For relatively smooth output with the
+sun placed for every degree of movement in the sky you should set the
+<i>step</i> to 4 minutes or less. <i>step</i><tt>=0.05</tt> is equivalent
+to every 3 minutes. Of course setting the time step to be very fine
+proportionally increases the module's running time.
+<p>The output units are in Wh per squared meter per given
+day [Wh/(m*m)/day]. The incidence angle and irradiance/irradiation maps are
+computed with the shadowing influence of relief by default. It is also possible
+for them to be computed without this influence using the planar flag (<i>-p</i>).
+In mountainous areas this can lead to very different results! The user should be
+aware that taking into account the shadowing effect of relief can slow
+down the speed of computation, especially when the sun altitude is low.
+
+<p>
+When considering the shadowing effect, speed and precision of computation
+can be controlled by the <i>distance_step</i> parameter, 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 a <i>distance_step</i> less than 1.0 (i.e. sampling
+points will be computed at <i>distance_step</i> * cellsize distance), <em>r.sun</em> takes
+the 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>distance_step</i>=1.0 should give reasonable results for most cases (e.g.
+on DEM). The <i>distance_step</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 output
+maps as zero values. Areas with NULL data are considered as no barrier with
+shadowing effect.
+
+<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>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.
+
+
+<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.
+
+<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>
+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>
+Average monthly values of the Linke turbidity coefficient for a mild climate
+in the northern hemisphere (see reference literature for your study area):
+
+<table border="1">
+<tr><th>Month</th><th>Jan</th><th>Feb</th><th>Mar</th><th>Apr</th><th>May</th><th>Jun</th><th>Jul</th><th>Aug</th><th>Sep</th><th>Oct</th><th>Nov</th><th>Dec</th><th>annual</th></tr>
+<tr><td>mountains</td><td>1.5</td><td>1.6</td><td>1.8</td><td>1.9</td><td>2.0</td><td>2.3</td><td>2.3</td><td>2.3</td><td>2.1</td><td>1.8</td><td>1.6</td><td>1.5</td><td>1.90</td></tr>
+<tr><td>rural</td><td>2.1</td><td>2.2</td><td>2.5</td><td>2.9</td><td>3.2</td><td>3.4</td><td>3.5</td><td>3.3</td><td>2.9</td><td>2.6</td><td>2.3</td><td>2.2</td><td>2.75</td></tr>
+<tr><td>city</td><td>3.1</td><td>3.2</td><td>3.5</td><td>4.0</td><td>4.2</td><td>4.3</td><td>4.4</td><td>4.3</td><td>4.0</td><td>3.6</td><td>3.3</td><td>3.1</td><td>3.75</td></tr>
+<tr><td>industrial</td><td>4.1</td><td>4.3</td><td>4.7</td><td>5.3</td><td>5.5</td><td>5.7</td><td>5.8</td><td>5.7</td><td>5.3</td><td>4.9</td><td>4.5</td><td>4.2</td><td>5.00</td></tr>
+</table>
+
+<p>Planned improvements include the use of the SOLPOS algorithm for solar
+geometry calculations and internal computation of aspect and slope.
+
+<h3>Solar time</h3>
+
+By default r.sun calculates times as true solar time, whereby solar noon is
+always exactly 12 o'clock everywhere in the current region. Depending on where
+the zone of interest is located in the related time zone, this may cause
+differences of up to an hour, in some cases (like Western Spain) even more.
+On top of this, the offset varies during the year according to the Equation
+of Time.
+<p>
+To overcome this problem, the user can use the option <em>civil_time=<timezone_offset></em>
+in r.sun to make it use real-world (wall clock) time. For example, for Central
+Europe the timezone offset is +1, +2 when daylight saving time is in effect.
+<p>
+<!-- WE DON'T KNOW, check source code:
+If the user use the <em>civil_time</em> parameter, also the longitude needs to
+be supplied as a raster map with the <em>long_in</em> parameter. Within a
+latlon location, such a map can be easily made with:
+
+<div class="code"><pre>
+r.mapcalc "lon_raster = x()"
+</pre></div>
+
+END OF WE DON'T KNOW
+-->
+
+<h3>Extraction of shadow maps</h3>
+A map of shadows can be extracted from the solar incidence angle map
+(incidout). Areas with zero values are shadowed. This will not work
+if the <em>-p</em> flag has been used.
+
+<h3>Large maps and out of memory problems</h3>
+
+With a large number or columns and rows, <b>r.sun</b> can consume
+significant amount of memory. While output raster maps are not
+partitionable, the input raster maps are using the <em>npartitions</em>
+parameter.
+
+In case of out of memory error (<tt>ERROR: G_malloc: out of memory</tt>), the
+<em>npartitions</em> parameter can be used to run a segmented calculation
+which consumes less memory during the computations.
+
+The amount of memory by <b>r.sun</b> is estimated as follows:
+
+<div class="code"><pre>
+# without input raster map partitioning:
+# memory requirements: 4 bytes per raster cell
+# rows,cols: rows and columns of current region (find out with g.region)
+# IR: number of input raster maps without horizon maps
+# OR: number of output raster maps
+memory_bytes = rows*cols*(IR*4 + horizon_steps + OR*4)
+
+# with input raster map partitioning:
+memory_bytes = rows*cols*((IR*4+horizon_steps)/npartitions + OR*4)
+</pre></div>
+
+<h2>EXAMPLES</h2>
+
+North Carolina example (considering also cast shadows):
+<div class="code"><pre>
+g.region raster=elevation -p
+
+# calculate horizon angles (to speed up the subsequent r.sun calculation)
+r.horizon elevation=elevation step=30 bufferzone=200 basename=horangle \
+ maxdistance=5000
+
+# slope + aspect
+r.slope.aspect elevation=elevation aspect=aspect.dem slope=slope.dem
+
+# calculate global radiation for day 180 at 2p.m., using r.horizon output
+r.sun elevation=elevation horizon_basename=horangle horizon_step=30 \
+ aspect=aspect.dem slope=slope.dem glob_rad=global_rad day=180 time=14
+# result: output global (total) irradiance/irradiation [W.m-2] for given day/time
+r.univar global_rad
+</pre></div>
+
+<p>
+Calculation of the integrated daily irradiation for a region in North-Carolina
+for a given day of the year at 30m resolution. Here day 172 (i.e., 21 June
+in non-leap years):
+
+<div class="code"><pre>
+g.region raster=elev_ned_30m -p
+
+# considering cast shadows
+r.sun elevation=elev_ned_30m linke_value=2.5 albedo_value=0.2 day=172 \
+ beam_rad=b172 diff_rad=d172 \
+ refl_rad=r172 insol_time=it172
+
+d.mon wx0
+# show irradiation raster map [Wh.m-2.day-1]
+d.rast.leg b172
+# show insolation time raster map [h]
+d.rast.leg it172
+</pre></div>
+
+We can compute the day of year from a specific date in Python:
+<div class="code"><pre>
+>>> import datetime
+>>> datetime.datetime(2014, 6, 21).timetuple().tm_yday
+172
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.horizon.html">r.horizon</a>,
+<a href="r.slope.aspect.html">r.slope.aspect</a>,
+<a href="r.sunhours.html">r.sunhours</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>
+</em>
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li> Hofierka, J., Suri, M. (2002): The solar radiation model for Open source
+GIS: implementation and applications. International
+GRASS users conference in Trento, Italy, September 2002.
+(<a href="http://skagit.meas.ncsu.edu/~jaroslav/trento/Hofierka_Jaroslav.pdf">PDF</a>)
+<li>
+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.
+<li>
+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.
+<li>
+Kasten, F. (1996). The Linke turbidity factor based on improved values of
+the integral Rayleigh optical thickness. Solar Energy, 56 (3), 239-244.
+<li>
+Kasten, F., Young, A. T. (1989). Revised optical air mass tables and approximation
+formula. Applied Optics, 28, 4735-4738.
+<li>
+Kittler, R., Mikler, J. (1986): Basis of the utilization of solar radiation
+(in Slovak). VEDA, Bratislava, p. 150.
+<li>
+Krcho, J. (1990). Morfometrická analza a digitálne modely georeliéfu
+(Morphometric analysis and digital models of georelief, in Slovak).
+VEDA, Bratislava.
+<li>
+Muneer, T. (1990). Solar radiation model for Europe. Building services engineering
+research and technology, 11, 4, 153-163.
+<li>
+Neteler, M., Mitasova, H. (2002): Open Source GIS: A GRASS GIS Approach, Kluwer
+Academic Publishers. (Appendix explains formula;
+<a href="http://www.grassbook.org/">r.sun script download</a>)
+<li>
+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.
+<li>
+Page, J., Albuisson, M., Wald, L. (2001). The European solar radiation atlas:
+a valuable digital tool. Solar Energy, 71, 81-83.
+<li>
+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.
+<li>
+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).
+<li>
+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>
+</ul>
+
+<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. This program is free software under the GNU General Public License (>=v2)
+<address>
+<a href="MAILTO:hofierka 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: 2016-02-19 02:42:11 +0100 (Pá, 19 úno 2016) $</i>
Added: grass-addons/grass7/raster/r.sun.mp/rsunglobals.h
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/rsunglobals.h (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/rsunglobals.h 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,61 @@
+
+/*******************************************************************************
+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.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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);
Added: grass-addons/grass7/raster/r.sun.mp/rsunlib.c
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/rsunlib.c (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/rsunlib.c 2017-04-03 04:54:36 UTC (rev 70825)
@@ -0,0 +1,697 @@
+/****************************************************************************
+ 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 the manual page 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
+
+ (C) 2011 by Hamish Bowman, and the GRASS Development Team
+****************************************************************************/
+/*
+ * 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.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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/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;
+}
+
+
+/* com_sol_const(): compute the Solar Constant corrected for the day of the
+ year. The Earth is closest to the Sun (Perigee) on about January 3rd,
+ it is furthest from the sun (Apogee) about July 6th. The 1367 W/m^2 solar
+ constant is at the average 1AU distance, but on Jan 3 it gets up to
+ around 1412.71 W/m^2 and on July 6 it gets down to around 1321 W/m^2.
+ This value is for what hits the top of the atmosphere before any energy
+ is attenuated. */
+double com_sol_const(int no_of_day)
+{
+ double I0, d1;
+
+ /* Solar constant: 1367.0 W/m^2.
+
+ Perigee offset: here we call Jan 2 at 8:18pm the Perigee, so day
+ number 2.8408. In angular units that's (2*pi * 2.8408 / 365.25) = 0.048869.
+
+ Orbital eccentricity: For Earth this is currently about 0.01672,
+ and so the distance to the sun varies by +/- 0.01672 from the
+ mean distance (1AU), so over the year the amplitude of the
+ function is 2*ecc = 0.03344.
+
+ And 365.25 is of course the number of days in a year.
+ */
+
+ /* 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 inputAngle;
+ double delt_lat, delt_lon;
+ double delt_lat_m, delt_lon_m;
+ 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;
+
+ /* 1852m * 60 * 0.0001rad * 180/pi= 636.67m */
+ delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
+ delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
+
+ delt_lat_m = delt_lat * (180/M_PI) * 1852*60;
+ delt_lon_m = delt_lon * (180/M_PI) * 1852*60 * cos(latitude);
+ delt_dist = sqrt(delt_lat_m * delt_lat_m + delt_lon_m * delt_lon_m);
+
+/*
+ sunVarGeom->stepsinangle = gridGeom->stepxy * sin(sunVarGeom->sunAzimuthAngle);
+ sunVarGeom->stepcosangle = gridGeom->stepxy * cos(sunVarGeom->sunAzimuthAngle);
+*/
+ sunVarGeom->stepsinangle = gridGeom->stepxy * delt_lat_m / delt_dist;
+ sunVarGeom->stepcosangle = gridGeom->stepxy * delt_lon_m / delt_dist;
+
+
+ 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 UNDEFZ; ?? */
+ 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;
+
+/* FIXME: please document all coefficients */
+ 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;
+
+/* FIXME: please document all coefficients */
+ 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);
+
+/* FIXME: please document all coefficients */
+ 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);
+
+/* FIXME: please document all coefficients */
+ 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);
+}
Added: grass-addons/grass7/raster/r.sun.mp/sunradstruct.h
===================================================================
--- grass-addons/grass7/raster/r.sun.mp/sunradstruct.h (rev 0)
+++ grass-addons/grass7/raster/r.sun.mp/sunradstruct.h 2017-04-03 04:54:36 UTC (rev 70825)
@@ -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.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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;
+
+};
More information about the grass-commit
mailing list