[GRASS-SVN] r35965 - grass/trunk/raster/r.sun2
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 19 19:07:48 EST 2009
Author: hamish
Date: 2009-02-19 19:07:47 -0500 (Thu, 19 Feb 2009)
New Revision: 35965
Modified:
grass/trunk/raster/r.sun2/local_proto.h
grass/trunk/raster/r.sun2/main.c
grass/trunk/raster/r.sun2/rsunglobals.h
grass/trunk/raster/r.sun2/rsunlib.c
grass/trunk/raster/r.sun2/sunradstruct.h
Log:
run tools/grass_indent.sh
Modified: grass/trunk/raster/r.sun2/local_proto.h
===================================================================
--- grass/trunk/raster/r.sun2/local_proto.h 2009-02-20 00:04:28 UTC (rev 35964)
+++ grass/trunk/raster/r.sun2/local_proto.h 2009-02-20 00:07:47 UTC (rev 35965)
@@ -1,10 +1,10 @@
/* main.c */
-void where_is_point(double *length, struct SunGeometryVarDay *sunVarGeom,
- struct GridGeometry *gridGeom);
-int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
- struct GridGeometry *gridGeom);
+void where_is_point(double *length, struct SunGeometryVarDay *sunVarGeom,
+ struct GridGeometry *gridGeom);
+int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
+ struct GridGeometry *gridGeom);
int useCivilTime();
void setUseCivilTime(int val);
@@ -27,39 +27,43 @@
double brad(double, double *bh, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
- struct SolarRadVar *sunRadVar);
-double drad(double, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
- struct SolarRadVar *sunRadVar);
-
-
-double brad_angle_loss(double, double *bh, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
+ struct SunGeometryVarSlope *sunSlopeGeom,
struct SolarRadVar *sunRadVar);
-double drad_angle_loss(double, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
+double drad(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);
+double brad_angle_loss(double, double *bh,
+ struct SunGeometryVarDay *sunVarGeom,
+ struct SunGeometryVarSlope *sunSlopeGeom,
+ struct SolarRadVar *sunRadVar);
+double drad_angle_loss(double, double bh, double *rr,
+ struct SunGeometryVarDay *sunVarGeom,
+ struct SunGeometryVarSlope *sunSlopeGeom,
+ struct SolarRadVar *sunRadVar);
+
+
+void com_par(struct SunGeometryConstDay *sungeom,
+ struct SunGeometryVarDay *sunVarGeom,
+ struct GridGeometry *gridGeom,
+ double latitude, double longitude);
void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
- struct GridGeometry *gridGeom);
+ struct GridGeometry *gridGeom);
double lumcline2(struct SunGeometryConstDay *sungeom,
- struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
- struct GridGeometry *gridGeom,
- unsigned char *horizonpointer);
+ struct SunGeometryVarDay *sunVarGeom,
+ struct SunGeometryVarSlope *sunSlopeGeom,
+ struct GridGeometry *gridGeom,
+ unsigned char *horizonpointer);
-typedef double (*BeamRadFunc)(double sh, double *bh, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
- struct SolarRadVar *sunRadVar);
-
-typedef double (*DiffRadFunc)(double sh, double bh, double *rr, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
- struct SolarRadVar *sunRadVar);
+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);
Modified: grass/trunk/raster/r.sun2/main.c
===================================================================
--- grass/trunk/raster/r.sun2/main.c 2009-02-20 00:04:28 UTC (rev 35964)
+++ grass/trunk/raster/r.sun2/main.c 2009-02-20 00:07:47 UTC (rev 35965)
@@ -1,3 +1,4 @@
+
/*******************************************************************************
r.sun: This program was writen by Jaro Hofierka in Summer 1993 and re-engineered
in 1996-1999. In cooperation with Marcel Suri and Thomas Huld from JRC in Ispra
@@ -107,7 +108,8 @@
struct SunGeometryVarSlope *sunSlopeGeom,
struct SolarRadVar *sunRadVar,
struct GridGeometry *gridGeom,
- unsigned char *horizonpointer, double latitude, double longitude);
+ unsigned char *horizonpointer, double latitude,
+ double longitude);
void calculate(double singleSlope, double singleAspect,
@@ -131,6 +133,7 @@
float **lumcl, **beam, **insol, **diff, **refl, **globrad;
unsigned char *horizonarray = NULL;
double civilTime;
+
/*
* double startTime, endTime;
*/
@@ -140,10 +143,12 @@
la_min = 90.;
double offsetx = 0.5, offsety = 0.5;
char *tt, *lt;
+
/*
* double slope;
*/
double o_orig, z1;
+
/*
* double lum_C11, lum_C13, lum_C22, lum_C31, lum_C33;
* double sinSolarAltitude; */
@@ -160,6 +165,7 @@
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;
*/
@@ -199,9 +205,9 @@
{
struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
*lin, *albedo, *longin, *alb, *latin, *lat, *coefbh, *coefdh,
- *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad, *glob_rad,
- *day, *step, *declin, *ltime, *dist, *horizon, *horizonstep,
- *numPartitions, *civilTime;
+ *incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
+ *glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
+ *horizonstep, *numPartitions, *civilTime;
}
parm;
@@ -233,7 +239,8 @@
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->description =
+ _("Name of the input elevation raster map [meters]");
parm.elevin->guisection = _("Input_options");
parm.aspin = G_define_option();
@@ -241,7 +248,8 @@
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->description =
+ _("Name of the input aspect map (terrain aspect or azimuth of the solar panel) [decimal degrees]");
parm.aspin->guisection = _("Input_options");
parm.aspect = G_define_option();
@@ -258,7 +266,8 @@
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->description =
+ _("Name of the input slope raster map (terrain slope or solar panel inclination) [decimal degrees]");
parm.slopein->guisection = _("Input_options");
parm.slope = G_define_option();
@@ -294,7 +303,8 @@
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->description =
+ _("Name of the ground albedo coefficient input raster map [-]");
parm.albedo->guisection = _("Input_options");
if (parm.albedo->answer == NULL) {
@@ -303,7 +313,8 @@
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->description =
+ _("A single value of the ground albedo coefficient [-]");
parm.alb->guisection = _("Input_options");
}
@@ -312,7 +323,8 @@
parm.latin->type = TYPE_STRING;
parm.latin->required = NO;
parm.latin->gisprompt = "old,cell,raster";
- parm.latin->description = _("Name of the latitudes input raster map [decimal degrees]");
+ parm.latin->description =
+ _("Name of the latitudes input raster map [decimal degrees]");
parm.latin->guisection = _("Input_options");
if (parm.latin->answer == NULL) {
@@ -320,7 +332,8 @@
parm.lat->key = "lat";
parm.lat->type = TYPE_DOUBLE;
parm.lat->required = NO;
- parm.lat->description = _("A single value of latitude [decimal degrees]");
+ parm.lat->description =
+ _("A single value of latitude [decimal degrees]");
parm.lat->guisection = _("Input_options");
}
@@ -329,7 +342,8 @@
parm.longin->type = TYPE_STRING;
parm.longin->required = NO;
parm.longin->gisprompt = "old,cell,raster";
- parm.longin->description = _("Name of the longitude input raster map [decimal degrees]");
+ parm.longin->description =
+ _("Name of the longitude input raster map [decimal degrees]");
parm.longin->guisection = _("Input_options");
parm.coefbh = G_define_option();
@@ -371,7 +385,8 @@
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->description =
+ _("Output incidence angle raster map (mode 1 only)");
parm.incidout->guisection = _("Output_options");
parm.beam_rad = G_define_option();
@@ -388,7 +403,8 @@
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->description =
+ _("Output insolation time raster map [h] (mode 2 only)");
parm.insol_time->guisection = _("Output_options");
parm.diff_rad = G_define_option();
@@ -430,7 +446,8 @@
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->description =
+ _("Time step when computing all-day radiation sums [decimal hours]");
parm.declin = G_define_option();
parm.declin->key = "declin";
@@ -444,7 +461,8 @@
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->description =
+ _("Local (solar) time (to be set for mode 1 only) [decimal hours]");
/*
* parm.startTime = G_define_option();
@@ -465,7 +483,8 @@
parm.dist->type = TYPE_DOUBLE;
parm.dist->answer = DIST;
parm.dist->required = NO;
- parm.dist->description = _("Sampling distance step coefficient (0.5-1.5)");
+ parm.dist->description =
+ _("Sampling distance step coefficient (0.5-1.5)");
parm.numPartitions = G_define_option();
parm.numPartitions->key = "numpartitions";
@@ -582,14 +601,14 @@
if (parm.horizonstep->answer != NULL) {
if (sscanf(parm.horizonstep->answer, "%lf", &horizonStep) != 1)
G_fatal_error(_("Error reading horizon step size"));
- if(horizonStep > 0.)
+ 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."));
- }
+ else if (useHorizonData()) {
+ G_fatal_error(_("If you use the horizon option you must also set the 'horizonstep' parameter."));
+ }
tt = parm.ltime->answer;
@@ -632,10 +651,10 @@
* if (parm.lat->answer != NULL)
* sscanf(parm.lat->answer, "%lf", &latitude);
*/
-/* HB 6/2008: why is the above commented out? instead of sscanf, maybe nicer to use:
- G_scan_lat(parm.lat->answer, &latitude);
- MN 2/2009: latitude doesn't exist! also G_scan_lat() does not exist in GRASS
-*/
+ /* HB 6/2008: why is the above commented out? instead of sscanf, maybe nicer to use:
+ G_scan_lat(parm.lat->answer, &latitude);
+ MN 2/2009: latitude doesn't exist! also G_scan_lat() does not exist in GRASS
+ */
if (parm.slopein->answer == NULL)
sscanf(parm.slope->answer, "%lf", &singleSlope);
@@ -713,27 +732,27 @@
/* Set up parameters for projection to lat/long if necessary */
if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {
- struct Key_Value *in_proj_info, *in_unit_info;
+ struct Key_Value *in_proj_info, *in_unit_info;
- if ((in_proj_info = G_get_projinfo()) == NULL)
- G_fatal_error
- (_("Can't get projection info of current location: please set latitude via 'lat' or 'latin' option!"));
+ if ((in_proj_info = G_get_projinfo()) == NULL)
+ G_fatal_error
+ (_("Can't get projection info of current location: please set latitude via 'lat' or 'latin' option!"));
- if ((in_unit_info = G_get_projunits()) == NULL)
- G_fatal_error(_("Can't get projection units of current location"));
+ if ((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"));
+ 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);
+ 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"));
+ /* Set output projection to latlong w/ same ellipsoid */
+ oproj.zone = 0;
+ oproj.meters = 1.;
+ sprintf(oproj.proj, "ll");
+ if ((oproj.pj = pj_latlong_from_proj(iproj.pj)) == NULL)
+ G_fatal_error(_("Unable to set up lat/long projection parameters"));
}
/**********end of parser - ******************************/
@@ -790,7 +809,7 @@
if ((mapset = G_find_cell2(elevin, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),elevin);
+ G_fatal_error(_("Raster map <%s> not found"), elevin);
fd1 = G_open_cell_old(elevin, mapset);
@@ -805,7 +824,7 @@
}
}
if ((mapset = G_find_cell2(slopein, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),slopein);
+ G_fatal_error(_("Raster map <%s> not found"), slopein);
fd3 = G_open_cell_old(slopein, mapset);
}
@@ -822,7 +841,7 @@
}
if ((mapset = G_find_cell2(aspin, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),aspin);
+ G_fatal_error(_("Raster map <%s> not found"), aspin);
fd2 = G_open_cell_old(aspin, mapset);
}
@@ -836,7 +855,7 @@
}
if ((mapset = G_find_cell2(linkein, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),linkein);
+ G_fatal_error(_("Raster map <%s> not found"), linkein);
fd4 = G_open_cell_old(linkein, mapset);
}
@@ -848,7 +867,7 @@
a[l] = (float *)G_malloc(sizeof(float) * (n));
}
if ((mapset = G_find_cell2(albedo, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),albedo);
+ G_fatal_error(_("Raster map <%s> not found"), albedo);
fd5 = G_open_cell_old(albedo, mapset);
}
@@ -860,7 +879,7 @@
la[l] = (float *)G_malloc(sizeof(float) * (n));
}
if ((mapset = G_find_cell2(latin, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),latin);
+ G_fatal_error(_("Raster map <%s> not found"), latin);
fd6 = G_open_cell_old(latin, mapset);
}
@@ -871,7 +890,7 @@
longitArray[l] = (float *)G_malloc(sizeof(float) * (n));
if ((mapset = G_find_cell2(longin, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),longin);
+ G_fatal_error(_("Raster map <%s> not found"), longin);
fd7 = G_open_cell_old(longin, mapset);
}
@@ -883,7 +902,7 @@
cbhr[l] = (float *)G_malloc(sizeof(float) * (n));
}
if ((mapset = G_find_cell2(coefbh, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),coefbh);
+ G_fatal_error(_("Raster map <%s> not found"), coefbh);
fr1 = G_open_cell_old(coefbh, mapset);
}
@@ -895,15 +914,15 @@
cdhr[l] = (float *)G_malloc(sizeof(float) * (n));
}
if ((mapset = G_find_cell2(coefdh, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"),coefdh);
+ G_fatal_error(_("Raster map <%s> not found"), coefdh);
fr2 = G_open_cell_old(coefdh, mapset);
}
if (useHorizonData()) {
if (horizonarray == NULL) {
horizonarray =
- (unsigned char *)G_malloc(sizeof(char) * arrayNumInt * numRows *
- n);
+ (unsigned char *)G_malloc(sizeof(char) * arrayNumInt *
+ numRows * n);
horizonbuf = (FCELL **) G_malloc(sizeof(FCELL *) * arrayNumInt);
fd_shad = (int *)G_malloc(sizeof(int) * arrayNumInt);
@@ -925,10 +944,11 @@
numDigits = (int)(log10(1. * arrayNumInt)) + 1;
sprintf(formatString, "%%s_%%0%dd", numDigits);
for (i = 0; i < arrayNumInt; i++) {
- horizonbuf[i] = G_allocate_f_raster_buf();
+ horizonbuf[i] = G_allocate_f_raster_buf();
sprintf(shad_filename, formatString, horizon, i);
if ((mapset = G_find_cell2(shad_filename, "")) == NULL)
- G_fatal_error(_("Horizon file no. %d <%s> not found"), i, shad_filename);
+ G_fatal_error(_("Horizon file no. %d <%s> not found"), i,
+ shad_filename);
fd_shad[i] = G_open_cell_old(shad_filename, mapset);
}
}
@@ -944,7 +964,8 @@
row_rev = m - row - 1;
rowrevoffset = row_rev - offset;
G_get_f_raster_row(fd_shad[i], horizonbuf[i], row);
- horizonpointer = horizonarray + (ssize_t) arrayNumInt * n * rowrevoffset;
+ horizonpointer =
+ horizonarray + (ssize_t) arrayNumInt *n * rowrevoffset;
for (j = 0; j < n; j++) {
horizonpointer[i] = (char)(rint(SCALING_FACTOR *
@@ -1132,8 +1153,8 @@
int OUTGR(void)
{
- FCELL *cell7 = NULL, *cell8 = NULL, *cell9 = NULL, *cell10 = NULL, *cell11 =
- NULL, *cell12 = NULL;
+ 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;
@@ -1379,13 +1400,15 @@
}
if ((diff_rad != NULL) || (glob_rad != NULL)) {
dra =
- drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
+ drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom,
+ sunRadVar);
diff_e += dfr * dra;
dra = 0.;
}
if ((refl_rad != NULL) || (glob_rad != NULL)) {
if ((diff_rad == NULL) && (glob_rad == NULL)) {
- drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom, sunRadVar);
+ drad(s0, bh, &rr, sunVarGeom, sunSlopeGeom,
+ sunRadVar);
}
refl_e += dfr * rr;
rr = 0.;
@@ -1458,6 +1481,7 @@
{
double sx, sy;
double dx, dy;
+
/* double adx, ady; */
int i, j;
@@ -1637,6 +1661,7 @@
double singleLinke, struct GridGeometry gridGeom)
{
int i, j, l;
+
/* double energy; */
int someRadiation;
int numRows;
@@ -1754,7 +1779,8 @@
* "local clock time". */
dayRad = 2. * M_PI * day / 365.25;
locTimeOffset =
- -0.128 * sin(dayRad - 0.04887) - 0.165 * sin(2 * dayRad + 0.34383);
+ -0.128 * sin(dayRad - 0.04887) - 0.165 * sin(2 * dayRad +
+ 0.34383);
/* Time offset due to timezone as input by user */
@@ -1826,29 +1852,29 @@
latitude *= deg2rad;
}
/* MN 2/2009: should it be??
- if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {
- */
+ if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {
+ */
if ((G_projection() != PROJECTION_LL)) {
- longitude = gridGeom.xp;
- latitude = gridGeom.yp;
+ longitude = gridGeom.xp;
+ latitude = gridGeom.yp;
- if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
- G_fatal_error("Error in pj_do_proj");
- }
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+ G_fatal_error("Error in pj_do_proj");
+ }
- la_max = AMAX1(la_max, latitude);
- la_min = AMIN1(la_min, latitude);
- latitude *= deg2rad;
- longitude *= deg2rad;
+ la_max = AMAX1(la_max, latitude);
+ la_min = AMIN1(la_min, latitude);
+ latitude *= deg2rad;
+ longitude *= deg2rad;
}
- else { /* ll projection */
- latitude = gridGeom.yp;
- longitude = gridGeom.xp;
- la_max = AMAX1(la_max, latitude);
- la_min = AMIN1(la_min, latitude);
- latitude *= deg2rad;
- longitude *= deg2rad;
+ else { /* ll projection */
+ latitude = gridGeom.yp;
+ longitude = gridGeom.xp;
+ la_max = AMAX1(la_max, latitude);
+ la_min = AMIN1(la_min, latitude);
+ latitude *= deg2rad;
+ longitude *= deg2rad;
}
if (coefbh != NULL) {
@@ -1869,10 +1895,12 @@
gridGeom.coslat = cos(-latitude);
sin_phi_l =
- -gridGeom.coslat * cos_u * sin_v + gridGeom.sinlat * sin_u;
+ -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;
+ 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;
@@ -1941,7 +1969,8 @@
G_short_history(glob_rad, "raster", &hist);
}
else
- G_fatal_error("Failed to init map history: no output maps requested!");
+ G_fatal_error
+ ("Failed to init map history: no output maps requested!");
sprintf(hist.edhist[0],
" ----------------------------------------------------------------");
@@ -2013,8 +2042,8 @@
sunRadVar.linke);
else
sprintf(hist.edhist[hist.edlinecnt],
- " Linke turbidity factor min-max: %.1f-%.1f", li_min,
- li_max);
+ " Linke turbidity factor min-max: %.1f-%.1f",
+ li_min, li_max);
hist.edlinecnt++;
if (albedo == NULL)
@@ -2023,8 +2052,8 @@
sunRadVar.alb);
else
sprintf(hist.edhist[hist.edlinecnt],
- " Ground albedo min-max: %.3f-%.3f", al_min,
- al_max);
+ " Ground albedo min-max: %.3f-%.3f",
+ al_min, al_max);
hist.edlinecnt++;
sprintf(hist.edhist[hist.edlinecnt],
Modified: grass/trunk/raster/r.sun2/rsunglobals.h
===================================================================
--- grass/trunk/raster/r.sun2/rsunglobals.h 2009-02-20 00:04:28 UTC (rev 35964)
+++ grass/trunk/raster/r.sun2/rsunglobals.h 2009-02-20 00:07:47 UTC (rev 35965)
@@ -1,3 +1,4 @@
+
/*******************************************************************************
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
@@ -30,9 +31,9 @@
#define EARTHRADIUS 6371000.
/* undefined value for terrain aspect */
-#define UNDEF 0.
+#define UNDEF 0.
/* internal undefined value for NULL */
-#define UNDEFZ -9999.
+#define UNDEFZ -9999.
/* Constant for calculating angular loss */
#define a_r 0.155
@@ -42,9 +43,9 @@
extern int arrayNumInt;
/*
-extern double xp;
-extern double yp;
-*/
+ extern double xp;
+ extern double yp;
+ */
extern double angular_loss_denom;
@@ -59,4 +60,3 @@
extern void (*func) (int, int);
-
Modified: grass/trunk/raster/r.sun2/rsunlib.c
===================================================================
--- grass/trunk/raster/r.sun2/rsunlib.c 2009-02-20 00:04:28 UTC (rev 35964)
+++ grass/trunk/raster/r.sun2/rsunlib.c 2009-02-20 00:07:47 UTC (rev 35965)
@@ -1,3 +1,4 @@
+
/*******************************************************************************
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
@@ -40,122 +41,121 @@
int civilTimeFlag;
int useCivilTime()
- {
- return civilTimeFlag;
- }
+{
+ return civilTimeFlag;
+}
void setUseCivilTime(int val)
- {
- civilTimeFlag=val;
- }
+{
+ civilTimeFlag = val;
+}
double angular_loss_denom;
-
+
void setAngularLossDenominator()
{
- angular_loss_denom=1./(1-exp(-1./a_r));
-
+ angular_loss_denom = 1. / (1 - exp(-1. / a_r));
+
}
int useShadowFlag;
int useShadow()
- {
- return useShadowFlag;
- }
+{
+ return useShadowFlag;
+}
void setUseShadow(int val)
- {
- useShadowFlag=val;
- }
+{
+ useShadowFlag = val;
+}
int useHorizonDataFlag;
int useHorizonData()
- {
- return useHorizonDataFlag;
- }
+{
+ return useHorizonDataFlag;
+}
void setUseHorizonData(int val)
- {
- useHorizonDataFlag=val;
- }
+{
+ useHorizonDataFlag = val;
+}
double timeOffset;
double getTimeOffset()
- {
- return timeOffset;
- }
+{
+ return timeOffset;
+}
void setTimeOffset(double val)
- {
- timeOffset=val;
- }
+{
+ timeOffset = val;
+}
double horizonInterval;
double getHorizonInterval()
- {
- return horizonInterval;
- }
+{
+ return horizonInterval;
+}
void setHorizonInterval(double val)
- {
- horizonInterval=val;
- }
+{
+ horizonInterval = val;
+}
double com_sol_const(int no_of_day)
- {
- double I0, d1;
+{
+ double I0, d1;
- /* v W/(m*m) */
- d1 = pi2 * no_of_day / 365.25;
- I0 = 1367. * (1 + 0.03344 * cos(d1 - 0.048869));
+ /* v W/(m*m) */
+ d1 = pi2 * no_of_day / 365.25;
+ I0 = 1367. * (1 + 0.03344 * cos(d1 - 0.048869));
- return I0;
- }
+ return I0;
+}
-
-
+
+
void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
struct GridGeometry *gridGeom)
- {
- double pom;
- double totOffsetTime;
+{
+ 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;
+ 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;
- }
- }
- }
+ 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;
+ }
+ }
+ }
}
}
@@ -164,536 +164,539 @@
-void com_par( struct SunGeometryConstDay *sungeom,
- struct SunGeometryVarDay *sunVarGeom, struct GridGeometry *gridGeom,
- double latitude, double longitude)
+void com_par(struct SunGeometryConstDay *sungeom,
+ struct SunGeometryVarDay *sunVarGeom,
+ struct GridGeometry *gridGeom, double latitude, double longitude)
{
- double pom, xpom, ypom;
+ double pom, xpom, ypom;
- double costimeAngle;
- double lum_Lx, lum_Ly;
+ double costimeAngle;
+ double lum_Lx, lum_Ly;
- double newLatitude, newLongitude;
- double inputAngle;
- double delt_lat, delt_lon;
- double delt_east, delt_nor;
- double delt_dist;
+ double newLatitude, newLongitude;
+ double inputAngle;
+ double delt_lat, delt_lon;
+ double delt_east, delt_nor;
+ double delt_dist;
- costimeAngle = cos(sungeom->timeAngle);
+ 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;
+ 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;
- }
- }
+ 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) */
+ 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);
+ 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 (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;
+ 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;
+ inputAngle = sunVarGeom->sunAzimuthAngle + pihalf;
+ inputAngle = (inputAngle >= pi2) ? inputAngle - pi2 : inputAngle;
- delt_lat = -0.0001*cos(inputAngle); /* Arbitrary small distance in latitude */
- delt_lon = 0.0001*sin(inputAngle)/cos(latitude);
+ delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
+ delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
- newLatitude = (latitude+delt_lat)*rad2deg;
- newLongitude = (longitude+delt_lon)*rad2deg;
+ newLatitude = (latitude + delt_lat) * rad2deg;
+ newLongitude = (longitude + delt_lon) * rad2deg;
- if ((G_projection() != PROJECTION_LL))
- {
- if (pj_do_proj(&newLongitude, &newLatitude, &oproj, &iproj) <0)
- {
- G_fatal_error("Error in pj_do_proj");
- }
- }
+ if ((G_projection() != PROJECTION_LL)) {
+ if (pj_do_proj(&newLongitude, &newLatitude, &oproj, &iproj) < 0) {
+ G_fatal_error("Error in pj_do_proj");
+ }
+ }
- delt_east=newLongitude-gridGeom->xp;
- delt_nor=newLatitude-gridGeom->yp;
+ delt_east = newLongitude - gridGeom->xp;
+ delt_nor = newLatitude - gridGeom->yp;
- delt_dist=sqrt(delt_east*delt_east+delt_nor*delt_nor);
+ delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor);
- sunVarGeom->stepsinangle = gridGeom->stepxy*delt_nor/delt_dist;
- sunVarGeom->stepcosangle = gridGeom->stepxy*delt_east/delt_dist;
+ sunVarGeom->stepsinangle = gridGeom->stepxy * delt_nor / delt_dist;
+ sunVarGeom->stepcosangle = gridGeom->stepxy * delt_east / delt_dist;
-/*
- sunVarGeom->stepsinangle = stepxy * sin(sunVarGeom->sunAzimuthAngle);
- sunVarGeom->stepcosangle = stepxy * cos(sunVarGeom->sunAzimuthAngle);
-*/
- sunVarGeom->tanSolarAltitude = tan(sunVarGeom->solarAltitude);
+ /*
+ sunVarGeom->stepsinangle = stepxy * sin(sunVarGeom->sunAzimuthAngle);
+ sunVarGeom->stepcosangle = stepxy * cos(sunVarGeom->sunAzimuthAngle);
+ */
+ sunVarGeom->tanSolarAltitude = tan(sunVarGeom->solarAltitude);
- return;
+ return;
}
-int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
+int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
struct GridGeometry *gridGeom)
{
- double z2;
- double curvature_diff;
- int succes = 0;
+ double z2;
+ double curvature_diff;
+ int succes = 0;
- if (sunVarGeom->zp == UNDEFZ)
- return 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;
+ 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) {
+ 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));
- if(succes!=1)
- {
- gridGeom->xx0 = gridGeom->xg0;
- gridGeom->yy0 = gridGeom->yg0;
- }
- return (succes);
+ 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;
+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;
+ int lowPos, highPos;
+ double timeoffset, horizPos;
+ double horizonHeight;
- func = cube;
- sunVarGeom->isShadow = 0;
+ 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 (useShadow()) {
+ length = 0;
-/*
- if(timeoffset<0.)
- timeoffset+=pi2;
- else if(timeoffset>pi2)
- timeoffset-=pi2;
- horizPos = arrayNumInt - timeoffset/horizonInterval;
-*/
-
- horizPos = timeoffset/getHorizonInterval();
+ 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;
+ */
- 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);
+ horizPos = timeoffset / getHorizonInterval();
- 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 */
- }
+ 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;
- } /* End if useHorizonData() */
- else
- {
- while ((r = searching(&length, sunVarGeom, gridGeom)) == 1)
- {
- if (r == 3)
- break; /* no test is needed */
- }
+ }
+ 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 (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 */
-
-
- }
+ if (z_orig != UNDEFZ)
+ {
- if (s < 0) return 0.;
- return (s);
+ s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l;
+ }
+ else
+ {
+ s = sunVarGeom->sinSolarAltitude;
+ }
+ */
+ s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; /* Jenco */
+ }
}
+ }
+ else {
+ /*
+ if (z_orig != UNDEFZ)
+ {
+ s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l;
+ }
+ else
+ {
+ s = sunVarGeom->sinSolarAltitude;
+ }
+ */
+ s = sunSlopeGeom->lum_C31_l * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) + sunSlopeGeom->lum_C33_l; /* Jenco */
+ }
+ if (s < 0)
+ return 0.;
+ return (s);
+}
+
+
+
double brad(double sh, double *bh, struct SunGeometryVarDay *sunVarGeom,
- struct SunGeometryVarSlope *sunSlopeGeom,
- struct SolarRadVar *sunRadVar)
+ struct SunGeometryVarSlope *sunSlopeGeom,
+ struct SolarRadVar *sunRadVar)
{
- double opticalAirMass, airMass2Linke, rayl, br;
- double drefract, temp1, temp2, h0refract;
- double elevationCorr;
+ double opticalAirMass, airMass2Linke, rayl, br;
+ double drefract, temp1, temp2, h0refract;
+ double elevationCorr;
- double locSolarAltitude;
-
- locSolarAltitude=sunVarGeom->solarAltitude;
+ double locSolarAltitude;
+ locSolarAltitude = sunVarGeom->solarAltitude;
- elevationCorr = exp(-sunVarGeom->z_orig / 8434.5);
- temp1 = 0.1594 + locSolarAltitude * (1.123 + 0.065656 * locSolarAltitude);
- temp2 = 1. + locSolarAltitude * (28.9344 + 277.3971 * locSolarAltitude);
- drefract = 0.061359 * temp1 / temp2; /* in radians */
- h0refract = locSolarAltitude + drefract;
- opticalAirMass = elevationCorr / (sin(h0refract) +
- 0.50572 * pow(h0refract * rad2deg + 6.07995, -1.6364));
- airMass2Linke = 0.8662 * sunRadVar->linke;
- if (opticalAirMass <= 20.)
- {
- rayl = 1. / (6.6296 +
- opticalAirMass * (1.7513 +
- opticalAirMass * (-0.1202 + opticalAirMass * (0.0065 - opticalAirMass * 0.00013))));
- }
- else
- {
- rayl = 1. / (10.4 + 0.718 * opticalAirMass);
- }
- *bh = sunRadVar->cbh * sunRadVar->G_norm_extra * sunVarGeom->sinSolarAltitude * exp(-rayl * opticalAirMass * airMass2Linke);
- if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.)
- br = *bh * sh / sunVarGeom->sinSolarAltitude;
- else
- br = *bh;
- return (br);
+ 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,
+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 p, opticalAirMass, airMass2Linke, rayl, br;
+ double drefract, temp1, temp2, h0refract;
- double locSolarAltitude;
-
- locSolarAltitude=sunVarGeom->solarAltitude;
+ double locSolarAltitude;
+ locSolarAltitude = sunVarGeom->solarAltitude;
- p = exp(-sunVarGeom->z_orig / 8434.5);
- temp1 = 0.1594 + locSolarAltitude * (1.123 + 0.065656 * locSolarAltitude);
- temp2 = 1. + locSolarAltitude * (28.9344 + 277.3971 * locSolarAltitude);
- drefract = 0.061359 * temp1 / temp2; /* in radians */
- h0refract = locSolarAltitude + drefract;
- opticalAirMass = p / (sin(h0refract) +
- 0.50572 * pow(h0refract * rad2deg + 6.07995, -1.6364));
- airMass2Linke = 0.8662 * sunRadVar->linke;
- if (opticalAirMass <= 20.)
- rayl =
- 1. / (6.6296 +
- opticalAirMass * (1.7513 +
- opticalAirMass * (-0.1202 + opticalAirMass * (0.0065 - opticalAirMass * 0.00013))));
- else
- rayl = 1. / (10.4 + 0.718 * opticalAirMass);
- *bh = sunRadVar->cbh * sunRadVar->G_norm_extra * sunVarGeom->sinSolarAltitude * exp(-rayl * opticalAirMass * airMass2Linke);
- if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.)
- br = *bh * sh / sunVarGeom->sinSolarAltitude;
- else
- br = *bh;
- br *= (1-exp(-sh/a_r))*angular_loss_denom;
-
- return (br);
+ 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 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;
+ double locSinSolarAltitude;
+ double locLinke;
+ locLinke = sunRadVar->linke;
+ locSinSolarAltitude = sunVarGeom->sinSolarAltitude;
+
cosslope = cos(sunSlopeGeom->slope);
sinslope = sin(sunSlopeGeom->slope);
tn = -0.015843 + locLinke * (0.030543 + 0.0003797 * locLinke);
A1b = 0.26463 + locLinke * (-0.061581 + 0.0031408 * locLinke);
if (A1b * tn < 0.0022)
- A1 = 0.0022 / tn;
+ A1 = 0.0022 / tn;
else
- A1 = A1b;
+ 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;
+ 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.;
+ 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.;
+ 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,
+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 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);
+ double diff_coeff_angleloss;
+ double refl_coeff_angleloss;
+ double c1;
+ double diff_loss_factor, refl_loss_factor;
+ double locSinSolarAltitude;
+ double locLinke;
+ locSinSolarAltitude = sunVarGeom->sinSolarAltitude;
+ locLinke = sunRadVar->linke;
+ cosslope = cos(sunSlopeGeom->slope);
+ sinslope = sin(sunSlopeGeom->slope);
- tn = -0.015843 + locLinke * (0.030543 + 0.0003797 * locLinke);
- A1b = 0.26463 + locLinke * (-0.061581 + 0.0031408 * locLinke);
- if (A1b * tn < 0.0022)
- A1 = 0.0022 / tn;
- else
- A1 = A1b;
- A2 = 2.04020 + locLinke * (0.018945 - 0.011161 * locLinke);
- A3 = -1.3025 + locLinke * (0.039231 + 0.0085079 * locLinke);
- fd = A1 + A2 * locSinSolarAltitude + A3 * locSinSolarAltitude * locSinSolarAltitude;
- dh = sunRadVar->cdh * sunRadVar->G_norm_extra * fd * tn;
- gh = bh + dh;
- if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.) {
- kb = bh / (sunRadVar->G_norm_extra * locSinSolarAltitude);
- r_sky = (1. + cosslope) / 2.;
- a_ln = sunVarGeom->solarAzimuth - sunSlopeGeom->aspect;
- ln = a_ln;
- if (a_ln > M_PI)
- ln = a_ln - pi2;
- else if (a_ln < -M_PI)
- ln = a_ln + pi2;
- a_ln = ln;
- fg = sinslope - sunSlopeGeom->slope * cosslope -
- M_PI * sin(sunSlopeGeom->slope / 2.) * sin(sunSlopeGeom->slope / 2.);
- if ((sunVarGeom->isShadow) || (sh <= 0.))
- fx = r_sky + fg * 0.252271;
- else if (sunVarGeom->solarAltitude >= 0.1) {
- fx = ((0.00263 - kb * (0.712 + 0.6883 * kb)) * fg + r_sky) * (1. -
- kb) +
- kb * sh / locSinSolarAltitude;
- }
- else if (sunVarGeom->solarAltitude < 0.1)
- fx = ((0.00263 - 0.712 * kb - 0.6883 * kb * kb) * fg +
- r_sky) * (1. - kb) + kb * sinslope * cos(a_ln) / (0.1 -
- 0.008 *
- sunVarGeom->solarAltitude);
- dr = dh * fx;
- /* refl. rad */
- *rr = sunRadVar->alb * gh * (1 - cosslope) / 2.;
+ 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 { /* plane */
- dr = dh;
- *rr = 0.;
- }
+ 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;
+ 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);
-}
+ return (dr);
+}
Modified: grass/trunk/raster/r.sun2/sunradstruct.h
===================================================================
--- grass/trunk/raster/r.sun2/sunradstruct.h 2009-02-20 00:04:28 UTC (rev 35964)
+++ grass/trunk/raster/r.sun2/sunradstruct.h 2009-02-20 00:07:47 UTC (rev 35965)
@@ -1,3 +1,4 @@
+
/*******************************************************************************
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
@@ -32,78 +33,78 @@
#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 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 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 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 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;
-
- };
+{
+ 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