[GRASS-SVN] r46638 - grass/branches/releasebranch_6_4/raster/r.sun2
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 9 05:57:13 EDT 2011
Author: neteler
Date: 2011-06-09 02:57:13 -0700 (Thu, 09 Jun 2011)
New Revision: 46638
Modified:
grass/branches/releasebranch_6_4/raster/r.sun2/rsunlib.c
Log:
partial code layout/comments sync with 6.5
Modified: grass/branches/releasebranch_6_4/raster/r.sun2/rsunlib.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.sun2/rsunlib.c 2011-06-09 08:35:04 UTC (rev 46637)
+++ grass/branches/releasebranch_6_4/raster/r.sun2/rsunlib.c 2011-06-09 09:57:13 UTC (rev 46638)
@@ -6,7 +6,7 @@
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
+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
@@ -102,11 +102,30 @@
}
-
+/* 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));
@@ -115,8 +134,6 @@
}
-
-
void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
struct GridGeometry *gridGeom)
{
@@ -144,12 +161,12 @@
}
else {
if (pom < 0) {
- /* G_debug(3,"\n Sun is ABOVE the surface during the whole day"); */
+ /* 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"); */
+ /* 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;
@@ -183,7 +200,6 @@
costimeAngle = cos(sungeom->timeAngle);
-
lum_Lx = -sungeom->lum_C22 * sin(sungeom->timeAngle);
lum_Ly = sungeom->lum_C11 * costimeAngle + sungeom->lum_C13;
sunVarGeom->sinSolarAltitude =
@@ -192,7 +208,7 @@
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"); */
+ /* G_debug(3,"\tSun is ABOVE area during the whole day"); */
sungeom->sunrise_time = 0;
sungeom->sunset_time = 24;
}
@@ -203,7 +219,7 @@
}
}
else {
- /* G_debug(3,"\tThe Sun is ON HORIZON during the whole day"); */
+ /* G_debug(3,"\tThe Sun is ON HORIZON during the whole day"); */
sungeom->sunrise_time = 0;
sungeom->sunset_time = 24;
}
@@ -220,7 +236,7 @@
if (fabs(pom) > EPS) {
sunVarGeom->solarAzimuth = lum_Ly / pom;
sunVarGeom->solarAzimuth = acos(sunVarGeom->solarAzimuth); /* horiz. angle of the Sun */
- /* solarAzimuth *= RAD; */
+ /* solarAzimuth *= RAD; */
if (lum_Lx < 0)
sunVarGeom->solarAzimuth = pi2 - sunVarGeom->solarAzimuth;
}
@@ -240,7 +256,7 @@
inputAngle = (inputAngle >= pi2) ? inputAngle - pi2 : inputAngle;
- delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in 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;
@@ -259,7 +275,6 @@
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;
@@ -344,16 +359,14 @@
if (useHorizonData()) {
/* Start is due east, sungeom->timeangle = -pi/2 */
- /*
- timeoffset = sungeom->timeAngle+pihalf;
- */
+ /* timeoffset = sungeom->timeAngle+pihalf; */
timeoffset = sunVarGeom->sunAzimuthAngle;
/*
if(timeoffset<0.)
- timeoffset+=pi2;
+ timeoffset+=pi2;
else if(timeoffset>pi2)
- timeoffset-=pi2;
+ timeoffset-=pi2;
horizPos = arrayNumInt - timeoffset/horizonInterval;
*/
@@ -388,7 +401,6 @@
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) {
@@ -397,8 +409,6 @@
}
-
-
if (r == 2) {
sunVarGeom->isShadow = 1; /* shadow */
}
@@ -435,8 +445,10 @@
}
+ /* if (s <= 0) return UNDEFZ; ?? */
if (s < 0)
return 0.;
+
return (s);
}
@@ -454,7 +466,7 @@
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);
@@ -500,7 +512,7 @@
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);
@@ -513,12 +525,10 @@
if (opticalAirMass <= 20.)
rayl =
1. / (6.6296 +
- opticalAirMass * (1.7513 +
- opticalAirMass * (-0.1202 +
- opticalAirMass *
- (0.0065 -
- opticalAirMass *
- 0.00013))));
+ opticalAirMass *
+ (1.7513 + opticalAirMass *
+ (-0.1202 + opticalAirMass *
+ (0.0065 - opticalAirMass * 0.00013))));
else
rayl = 1. / (10.4 + 0.718 * opticalAirMass);
*bh =
@@ -555,6 +565,7 @@
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)
@@ -584,21 +595,19 @@
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;
+ 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);
+ 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 */
+ else { /* plane */
dr = dh;
*rr = 0.;
}
@@ -629,9 +638,10 @@
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
@@ -643,11 +653,14 @@
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)
@@ -665,15 +678,14 @@
}
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);
+ 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 */
+ else { /* plane */
dr = dh;
*rr = 0.;
}
More information about the grass-commit
mailing list