[GRASS-SVN] r52877 - grass-addons/grass6/raster/r.sun.angle
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Aug 24 11:13:55 PDT 2012
Author: mmetz
Date: 2012-08-24 11:13:55 -0700 (Fri, 24 Aug 2012)
New Revision: 52877
Removed:
grass-addons/grass6/raster/r.sun.angle/README
grass-addons/grass6/raster/r.sun.angle/TODO
grass-addons/grass6/raster/r.sun.angle/g_solposition.c
grass-addons/grass6/raster/r.sun.angle/global.h
Modified:
grass-addons/grass6/raster/r.sun.angle/Makefile
grass-addons/grass6/raster/r.sun.angle/description.html
grass-addons/grass6/raster/r.sun.angle/main.c
grass-addons/grass6/raster/r.sun.angle/solpos00.c
grass-addons/grass6/raster/r.sun.angle/solpos00.h
Log:
r.sun.angle: sun elevation, azimuth, photoperiod
Modified: grass-addons/grass6/raster/r.sun.angle/Makefile
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/Makefile 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/Makefile 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = r.sunmask
+PGM = r.sun.angle
LIBES = $(GPROJLIB) $(GISLIB)
DEPENDENCIES = $(GPROJDEP) $(GISDEP)
Deleted: grass-addons/grass6/raster/r.sun.angle/README
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/README 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/README 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,45 +0,0 @@
-
-r.sunmask can be used in two ways:
-
-a) define sun position:
-r.sunmask elev=dgm25 output=dgm25.shadow alt=47.394203 azim=138.604904
-
-b) define location, date, time
-r.sunmask -v el=dgm5 out=shadow year=1997 month=9 day=25 hour=11 minute=15 \
- sec=0 timezone=1 east=3570720 north=5766095
-
-or (using map center coordinates)
-r.sunmask -v el=dgm5 out=shadow year=1997 month=9 day=25 hour=11 minute=15 \
- sec=0 timezone=1
-
-Calculating sun position... (using solpos (V. 11 April 2001) from NREL)
- 1997.09.25, daynum 268, time: 11:15:00
- long: 10.030458, lat: 52.025555, timezone: 2.000000
- Solar position: sun azimuth 145.121445,
- sun angle above horz.(refraction corrected) 31.662413
- Sunrise time (without refraction): 07:16
- Sunset time (without refraction): 19:07
-Calculating shadows from DEM...
-
---------------------------
-
-Since the solpos algorithm expects lat/long coordinates, the module
-transforms the local coordinates to lat/long on-the-fly
-
---------------------------
-daylight savings: Rather than convert time to GMT, the solpos
-algorithm uses what is called Local Standard Time, which is generally
-defined politically as an offset from GMT. So the key is the offset from
-GMT, which the solpos Time Zone parameter. If the user specifies clock time
-(different for winter and summer), he/she would have to change the Time Zone
-parameter seasonally in r.sunmask (timezone parameter).
-
-
---------------------------
-
-solpos questions can be directed to
-> Steve Wilcox (Stephen_Wilcox at nrel.gov)
-> National Renewable Energy Laboratory
-> 1617 Cole Boulevard, Golden, CO 80401-3393
-> voice: 1-303-275-4061; fax: 1-303-275-4675
- http://rredc.nrel.gov/solar/codes_algs/solpos/
Deleted: grass-addons/grass6/raster/r.sun.angle/TODO
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/TODO 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/TODO 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,2 +0,0 @@
-improve speed for FP maps
-
Modified: grass-addons/grass6/raster/r.sun.angle/description.html
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/description.html 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/description.html 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,114 +1,36 @@
<h2>DESCRIPTION</h2>
-<em>r.sunmask</em> creates an output map layer based on an input elevation
-raster map layer and the sun position. The output map layer contains the
-cast shadow areas arising from sun shine and elevations. The user can define
-the sun position either directly or the module calculates it from given
-location and date/time parameters using the
-<a href=http://rredc.nrel.gov/>NREL</a> sun position algorithm. So either
-"A:"-parameters to specify the exact known sun position or "B:-parameters"
-to specify date/time for sun position calculation by r.sunmask itself have
-to be used.
+<em>r.sun.angle</em> calculates sun elevation and sun azimuth angles for
+the given time of day and each grid cell in the current region.
+Additionally, the photoperiod (sunshine hours) can be calculated.
<p>
-The module performs sunset/sunrise checks and refraction correction for sun
-position calculation. Local coordinate systems are internally transformed to
-latitude/longitude for the SOLPOS algorithm. The elevation is not considered
-in the sunset/sunrise calculations.
-
-<h2>NOTES</h2>
-
-r.sunmask and daylight savings: Rather than convert time to GMT, the solpos
-algorithm uses what is called Local Standard Time, which is generally
-defined politically as an offset from GMT. So the key is the offset from
-GMT, which the solpos Time Zone parameter. If the user specifies clock time
-(different for winter and summer), he/she would have to change the Time Zone
-parameter seasonally in r.sunmask (timezone parameter). See also
-<a href="http://en.wikipedia.org/wiki/Daylight_saving_time_by_country">Daylight saving time by region and country</a>.
+Sun elevation, height, height angle, or solar altitude angle is the
+angle in degrees between the horizon and a line that points from the
+site towards the centre of the sun.
<p>
-
-Note: In latitude/longitude locations the position coordinates pair
-(east/west) has to be specified in decimal degree (not D:M:S). If
-not specified, the map center's coordinates will be used.
-Also <em>g.region -l</em> displays the map center's coordinates.
-
+The sun azimuth angle is here defined as the azimuth angle in degrees
+of the sun from due north in a clockwise direction.
<p>
-Note for module usage with <em>-g</em> flag and calculations
-close to sunset/sunrise:
+The time used here is defined such that 12:00 (high noon) is the time
+when the sun has reached its highest point in the sky at the current site,
+unless the <em>-t</em> flag is used in which case time is interpreted as
+Greenwich standard time.
+<p>
+If a <em>sunhour</em> output map is specified, the module calculates
+sunshine hours for the given day. This option requires both Greenwhich
+standard time (<em>-t</em> flag) and the solpos algorithm of NREL
+(<em>-s</em> flag).
-<pre>
- [...]
- sunangleabovehorizont=0.434240
- sunrise=07:59:19
- sunset=16:25:17
- Time (07:59:02) is before sunrise (07:59:19)!
- WARNING: Nothing to calculate. Please verify settings.
- No map calculation requested. Finished.
-</pre>
+<h2>EXAMPLES</h2>
-In above calculation appears to be a mistake as
-the program indicates that we are before sunrise while
-the <i>sun angle above horizon</i> is already positive.
-The reason is that <i>sun angle above horizon</i> is
-calculated with correction for atmosphere refraction while
-<i>sunrise</i> and <i>sunset</i> are calculated <b>without</b>
-correction for atmosphere refraction. The output without
-<em>-g</em> flag contains related indications.
-
-<h2>EXAMPLE</H2>
-
-Example for North Carolina sample data set for the calculation
-of sun position angles and more:
-
+Calculate sun's elevation angle map for 2010-10-11 at 14:00h solar time:
<div class="code"><pre>
-# set the region to a place near Raleigh (NC)
-g.region rast=elev_lid792_1m -p
-
-# compute the sun position
-r.sunmask -s --v elev_lid792_1m out=dummy year=2012 month=2 \
- day=22 hour=10 minute=30 sec=0 timezone=-5
-Using map center coordinates: 638650.000000 220375.000000
-Calculating sun position... (using solpos (V. 11 April 2001) from NREL)
-2012/02/22, daynum: 53, time: 10:30:00 (decimal time: 10.500000)
-long: -78.678856, lat: 35.736160, timezone: -5.000000
-Solar position: sun azimuth: 143.006424, sun angle above horz. (refraction corrected): 36.233883
-Sunrise time (without refraction): 06:58:11
-Sunset time (without refraction): 17:58:47
-
-# with -g flag
-r.sunmask -s -g --v elev_lid792_1m out=dummy year=2012 month=2 \
- day=22 hour=10 minute=30 sec=0 timezone=-5
-Using map center coordinates: 638650.000000 220375.000000
-Calculating sun position... (using solpos (V. 11 April 2001) from NREL)
-date=2012/02/22
-daynum=53
-time=10:30:00
-decimaltime=10.500000
-longitudine=-78.678856
-latitude=35.736160
-timezone=-5.000000
-sunazimuth=143.006424
-sunangleabovehorizon=36.233883
-sunrise=06:58:11
-sunset=17:58:47
+r.sun.angle elevation=sun_elev year=2010 month=10 day=11 hour=14 minute=00
</pre></div>
-
-<h2>Acknowledgements</h2>
-Acknowledgements: National Renewable Energy Laboratory for
-their <a href=http://rredc.nrel.gov/solar/codes_algs/solpos/>SOLPOS 2.0</a> sun position
-algorithm.
-
-<h2>SEE ALSO</h2>
-<em>
-<a href="g.region.html">g.region</a>,
-<a href="r.sun.html">r.sun</a>,
-<a href="r.slope.aspect.html">r.slope.aspect</a>
-</em>
-
<h2>AUTHOR</h2>
-Janne Soimasuo, Finland 1994<br>
-update to FP by Huidae Cho 2001<br>
-added solpos algorithm feature by Markus Neteler 2001
+Markus Metz
+
<p><i>Last changed: $Date$</i>
Deleted: grass-addons/grass6/raster/r.sun.angle/g_solposition.c
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/g_solposition.c 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/g_solposition.c 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,173 +0,0 @@
-/*
- * G_calc_solar_position() calculates solar position parameters from
- * given position, date and time
- *
- * Written by Markus Neteler <neteler at geog.uni-hannover.de>
- * with kind help from Morten Hulden <morten untamo.net>
- *
- *----------------------------------------------------------------------------
- * using solpos.c with permission from
- * From rredc at nrel.gov Wed Mar 21 18:37:25 2001
- * Message-Id: <v04220805b6de9b1ad6ff@[192.174.39.30]>
- * Mary Anderberg
- * http://rredc.nrel.gov
- * National Renewable Energy Laboratory
- * 1617 Cole Boulevard
- * Golden, Colorado, USA 80401
- *
- * http://rredc.nrel.gov/solar/codes_algs/solpos/
- *
- * G_calc_solar_position is based on: soltest.c
- * by
- * Martin Rymes
- * National Renewable Energy Laboratory
- * 25 March 1998
- *----------------------------------------------------------------------------*/
-
-/* uncomment to get debug output */
-
-#include <math.h>
-#include <string.h>
-#include <stdio.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include <grass/gprojects.h>
-#include "solpos00.h"
-
-struct posdata pd, *pdat; /* declare solpos data struct and a pointer for it */
-
-
-long calc_solar_position(double longitude, double latitude, double timezone,
- int year, int month, int day, int hour, int minute,
- int second)
-{
-
- /* Note: this code is valid from year 1950 to 2050 (solpos restriction)
-
- - the algorithm will compensate for leap year.
- - longitude, latitude: decimal degree
- - timezone: DO NOT ADJUST FOR DAYLIGHT SAVINGS TIME.
- - timezone: negative for zones west of Greenwich
- - lat/long: east and north positive
- - the atmospheric refraction is calculated for 1013hPa, 15°C
- - time: local time from your watch
- */
-
- long retval; /* to capture S_solpos return codes */
- struct Key_Value *in_proj_info, *in_unit_info; /* projection information of input map */
- struct pj_info iproj; /* input map proj parameters */
- struct pj_info oproj; /* output map proj parameters */
- extern struct Cell_head window;
- int inside;
-
-
- /* we don't like to run G_calc_solar_position in xy locations */
- if (window.proj == 0)
- G_fatal_error(_("Unable to calculate sun position in un-projected locations. "
- "Specify sunposition directly."));
-
- pdat = &pd; /* point to the structure for convenience */
-
- /* Initialize structure to default values. (Optional only if ALL input
- parameters are initialized in the calling code, which they are not
- in this example.) */
-
- S_init(pdat);
-
- /* check if given point is in current window */
- G_debug(1, "window.north: %f, window.south: %f\n", window.north,
- window.south);
- G_debug(1, "window.west: %f, window.east : %f\n", window.west,
- window.east);
-
- inside = 0;
- if (latitude >= window.south && latitude <= window.north &&
- longitude >= window.west && longitude <= window.east)
- inside = 1;
- if (!inside)
- G_warning(_("Specified point %f, %f outside of current region, "
- "is that intended? Anyway, it will be used."),
- longitude, latitude);
-
- /* if coordinates are not in lat/long format, transform them: */
- if ((G_projection() != PROJECTION_LL) && window.proj != 0) {
- G_debug(1, "Transforming input coordinates to lat/long (req. for solar position)");
-
- /* read current projection info */
- if ((in_proj_info = G_get_projinfo()) == NULL)
- G_fatal_error(_("Unable to get projection info of current location"));
-
- if ((in_unit_info = G_get_projunits()) == NULL)
- G_fatal_error(_("Unable to get projection units of current location"));
-
- if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
- G_fatal_error(_("Unable to get projection key values of current location"));
-
- G_free_key_value(in_proj_info);
- G_free_key_value(in_unit_info);
-
- /* Try using pj_print_proj_params() instead of all this */
- G_debug(1, "Projection found in location:");
- G_debug(1, "IN: meter: %f zone: %i proj: %s (iproj struct)",
- iproj.meters, iproj.zone, iproj.proj);
- G_debug(1, "IN coord: longitude: %f, latitude: %f", longitude,
- latitude);
- /* see src/include/projects.h, struct PJconsts */
-
- /* set output projection to lat/long for solpos */
- 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");
-
- /* XX do the transform
- * outx outy in_info out_info */
- if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
- G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
- }
-
- G_debug(1, "Transformation to lat/long:");
- G_debug(1, "OUT: longitude: %f, latitude: %f", longitude,
- latitude);
-
- } /* transform if not LL */
-
- pdat->longitude = longitude; /* Note that latitude and longitude are */
- pdat->latitude = latitude; /* in DECIMAL DEGREES, not Deg/Min/Sec */
- pdat->timezone = timezone; /* DO NOT ADJUST FOR DAYLIGHT SAVINGS TIME. */
-
- pdat->year = year; /* The year */
- pdat->function &= ~S_DOY;
- pdat->month = month;
- pdat->day = day; /* the algorithm will compensate for leap year, so
- you just count days). S_solpos can be
- configured to accept day-of-the year */
-
- /* The time of day (STANDARD (GMT) time) */
-
- pdat->hour = hour;
- pdat->minute = minute;
- pdat->second = second;
-
- /* Let's assume that the temperature is 20 degrees C and that
- the pressure is 1013 millibars. The temperature is used for the
- atmospheric refraction correction, and the pressure is used for the
- refraction correction and the pressure-corrected airmass. */
-
- pdat->temp = 20.0;
- pdat->press = 1013.0;
-
- /* Finally, we will assume that you have a flat surface
- facing nowhere, tilted at latitude. */
-
- pdat->tilt = pdat->latitude; /* tilted at latitude */
- pdat->aspect = 180.0;
-
- /* perform the calculation */
- retval = S_solpos(pdat); /* S_solpos function call: returns long value */
- S_decode(retval, pdat); /* prints an error in case of problems */
-
- return retval;
-
-}
Deleted: grass-addons/grass6/raster/r.sun.angle/global.h
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/global.h 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/global.h 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,12 +0,0 @@
-#ifdef MAIN
-#define GLOBAL
-#else
-#define GLOBAL extern
-#endif
-
-GLOBAL float asol, phi0, sun_zenith, sun_azimuth; /* from nadir, from north */
-GLOBAL int sunset;
-
-/* proto */
-long calc_solar_position(double, double, double, int, int, int, int, int,
- int);
Modified: grass-addons/grass6/raster/r.sun.angle/main.c
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/main.c 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/main.c 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,44 +1,19 @@
/****************************************************************************
*
- * MODULE: r.sunmask
- * AUTHOR(S): Janne Soimasuo, Finland 1994 (original contributor)
- * update to FP by Huidae Cho <grass4u gmail.com> 2001
- * added solpos algorithm feature by Markus Neteler 2001
- * Brad Douglas <rez touchofmadness.com>, Glynn Clements <glynn gclements.plus.com>,
- * Hamish Bowman <hamish_b yahoo.com>, Paul Kelly <paul-grass stjohnspoint.co.uk>
- * PURPOSE:
- * COPYRIGHT: (C) 1999-2006 by the GRASS Development Team
+ * MODULE: r.sun.angle
+ * AUTHOR(S): Markus Metz
+ * PURPOSE: Calculates solar azimuth and angle
+ * COPYRIGHT: (C) 2010 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
-/*
- * r.sunmask:
- * Calculates the cast shadow areas from a DEM
- *
- * input: DEM (int, float, double)
- * output: binary shadow map
- * no shadow: null()
- * shadow: 1
- *
- * Algorithm source: unknown (Janne Soimasuo?)
- * Original module author: Janne Soimasuo, Finland 1994
- *
- * GPL >= 2
- *
- **********************
- * Added solpol sun position calculation:
- * Markus Neteler 4/2001
- *
- **********************
- * MN 2/2001: attempt to update to FP
- * Huidae Cho 3/2001: FP update done
- * but it's somewhat slow with non-CELL maps
- *********************************************************************/
+/* TODO: always use solpos if time is Greenwich standard time */
+
#define MAIN
#include <stdio.h>
@@ -46,105 +21,78 @@
#include <string.h>
#include <math.h>
#include <grass/gis.h>
+#include <grass/gprojects.h>
#include <grass/glocale.h>
-
-#include "global.h"
#include "solpos00.h"
-/* to be displayed in r.sunmask */
-static char *SOLPOSVERSION = "11 April 2001";
+void set_solpos_time(struct posdata *, int, int, int, int, int, int);
+void set_solpos_longitude(struct posdata *, double );
+int roundoff(double *);
-extern struct posdata pd, *pdat; /* declare a posdata struct and a pointer for
- it (if desired, the structure could be
- allocated dynamically with G_malloc) */
-struct Cell_head window;
-
-union RASTER_PTR
-{
- void *v;
- CELL *c;
- FCELL *f;
- DCELL *d;
-};
-
-#ifdef RASTER_VALUE_FUNC
-double raster_value(union RASTER_PTR buf, int data_type, int col);
-#else
-#define raster_value(buf, data_type, col) ((double)(data_type == CELL_TYPE ? buf.c[col] : (data_type == FCELL_TYPE ? buf.f[col] : buf.d[col])))
-#endif
-
int main(int argc, char *argv[])
{
- char *mapset;
- extern struct Cell_head window;
- union RASTER_PTR elevbuf, tmpbuf, outbuf;
- CELL min, max;
- DCELL dvalue, dvalue2, dmin, dmax;
- struct History hist;
- RASTER_MAP_TYPE data_type;
- struct Range range;
- struct FPRange fprange;
- double drow, dcol;
- int elev_fd, output_fd, zeros;
+ struct GModule *module;
struct
{
- struct Option *opt1, *opt2, *opt3, *opt4, *north, *east, *year,
- *month, *day, *hour, *minutes, *seconds, *timezone;
+ struct Option *elev, *azimuth, *sunhours, *year,
+ *month, *day, *hour, *minutes, *seconds;
+ struct Flag *gst_time, *solpos;
} parm;
- struct Flag *flag1, *flag2, *flag3, *flag4;
- struct GModule *module;
- char *name, *outname;
- double dazi, dalti;
- double azi, alti;
- double nstep, estep;
- double maxh;
- double east, east1, north, north1;
- int row1, col1;
- char OK;
- double timezone;
+ struct Cell_head window;
+ FCELL *elevbuf, *azimuthbuf, *sunhourbuf;
+ struct History hist;
+
+ /* projection information of input map */
+ struct Key_Value *in_proj_info, *in_unit_info;
+ struct pj_info iproj; /* input map proj parameters */
+ struct pj_info oproj; /* output map proj parameters */
+
+ char *elev_name, *azimuth_name, *sunhour_name;
+ int elev_fd, azimuth_fd, sunhour_fd;
+ double ha, ha_cos, s_gamma, s_elevation, s_azimuth;
+ double s_declination, sd_sin, sd_cos;
+ double se_sin, sa_cos;
+ double east, east_ll, north, north_ll;
+ double north_gc, north_gc_sin, north_gc_cos; /* geocentric latitude */
+ double ba2;
int year, month, day, hour, minutes, seconds;
- long retval;
- int solparms, locparms, use_solpos;
- double sunrise, sunset, current_time;
- int sretr = 0, ssetr = 0, sretr_sec = 0, ssetr_sec = 0;
- double dsretr, dssetr;
+ int doy; /* day of year */
+ int row, col, nrows, ncols;
+ int do_reproj = 0;
+ int lst_time = 1;
+ int use_solpos = 0;
+ struct posdata pd;
G_gisinit(argv[0]);
module = G_define_module();
- module->keywords = _("raster, sun position");
- module->label = _("Calculates cast shadow areas from sun position and elevation raster map.");
- module->description = _("Either exact sun position (A) is specified, or date/time to calculate "
- "the sun position (B) by r.sunmask itself.");
+ module->keywords = _("raster");
+ module->label = _("Calculates solar elevation and azimuth.");
+ module->description = _("Solar elevation: the angle between the direction of the geometric center "
+ "of the sun's apparent disk and the (idealized) horizon. "
+ "Solar azimuth: the angle from due north in clockwise direction.");
- parm.opt1 = G_define_standard_option(G_OPT_R_ELEV);
+ parm.elev = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.elev->key = "elevation";
+ parm.elev->label = _("Output raster map with solar elevation angle");
+ parm.elev->required = NO;
+
+ parm.azimuth = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.azimuth->key = "azimuth";
+ parm.azimuth->label = _("Output raster map with solar azimuth angle");
+ parm.azimuth->required = NO;
- parm.opt2 = G_define_standard_option(G_OPT_R_OUTPUT);
- parm.opt2->required = NO;
-
- parm.opt3 = G_define_option();
- parm.opt3->key = "altitude";
- parm.opt3->type = TYPE_DOUBLE;
- parm.opt3->required = NO;
- parm.opt3->options = "0-89.999";
- parm.opt3->description =
- _("Altitude of the sun above horizon, degrees (A)");
- parm.opt3->guisection = _("Position");
+ parm.sunhours = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.sunhours->key = "sunhour";
+ parm.sunhours->label = _("Output raster map with sunshine hours");
+ parm.sunhours->description = _("Sunshine hours require solpos and Greenwich standard time");
+ parm.sunhours->required = NO;
- parm.opt4 = G_define_option();
- parm.opt4->key = "azimuth";
- parm.opt4->type = TYPE_DOUBLE;
- parm.opt4->required = NO;
- parm.opt4->options = "0-360";
- parm.opt4->description =
- _("Azimuth of the sun from the north, degrees (A)");
- parm.opt4->guisection = _("Position");
-
parm.year = G_define_option();
parm.year->key = "year";
parm.year->type = TYPE_INTEGER;
- parm.year->required = NO;
- parm.year->description = _("Year (B)");
+ parm.year->required = YES;
+ parm.year->description = _("Year");
parm.year->options = "1950-2050";
parm.year->guisection = _("Time");
@@ -152,411 +100,415 @@
parm.month->key = "month";
parm.month->type = TYPE_INTEGER;
parm.month->required = NO;
- parm.month->description = _("Month (B)");
+ parm.month->label = _("Month");
+ parm.month->description = _("If not given, day is interpreted as day of the year");
parm.month->options = "0-12";
parm.month->guisection = _("Time");
parm.day = G_define_option();
parm.day->key = "day";
parm.day->type = TYPE_INTEGER;
- parm.day->required = NO;
- parm.day->description = _("Day (B)");
- parm.day->options = "0-31";
+ parm.day->required = YES;
+ parm.day->description = _("Day");
+ parm.day->options = "0-365";
parm.day->guisection = _("Time");
parm.hour = G_define_option();
parm.hour->key = "hour";
parm.hour->type = TYPE_INTEGER;
- parm.hour->required = NO;
- parm.hour->description = _("Hour (B)");
+ parm.hour->required = YES;
+ parm.hour->description = _("Hour");
parm.hour->options = "0-24";
+ parm.hour->answer = "12";
parm.hour->guisection = _("Time");
parm.minutes = G_define_option();
parm.minutes->key = "minute";
parm.minutes->type = TYPE_INTEGER;
- parm.minutes->required = NO;
- parm.minutes->description = _("Minutes (B)");
+ parm.minutes->required = YES;
+ parm.minutes->description = _("Minutes");
parm.minutes->options = "0-60";
+ parm.minutes->answer = "0";
parm.minutes->guisection = _("Time");
parm.seconds = G_define_option();
parm.seconds->key = "second";
parm.seconds->type = TYPE_INTEGER;
- parm.seconds->required = NO;
- parm.seconds->description = _("Seconds (B)");
+ parm.seconds->required = YES;
+ parm.seconds->description = _("Seconds");
parm.seconds->options = "0-60";
+ parm.seconds->answer = "0";
parm.seconds->guisection = _("Time");
- parm.timezone = G_define_option();
- parm.timezone->key = "timezone";
- parm.timezone->type = TYPE_INTEGER;
- parm.timezone->required = NO;
- parm.timezone->label =
- _("Timezone");
- parm.timezone->description = _("East positive, offset from GMT, also use to adjust daylight savings");
- parm.timezone->guisection = _("Time");
+ parm.gst_time = G_define_flag();
+ parm.gst_time->key = 't';
+ parm.gst_time->description = _("Time is Greenwich standard time, not local sidereal time");
- parm.east = G_define_option();
- parm.east->key = "east";
- parm.east->key_desc = "value";
- parm.east->type = TYPE_STRING;
- parm.east->required = NO;
- parm.east->label =
- _("Easting coordinate (point of interest)");
- parm.east->description = _("Default: map center");
- parm.east->guisection = _("Position");
+ parm.solpos = G_define_flag();
+ parm.solpos->key = 's';
+ parm.solpos->description = _("Use solpos algorithm of NREL");
- parm.north = G_define_option();
- parm.north->key = "north";
- parm.north->key_desc = "value";
- parm.north->type = TYPE_STRING;
- parm.north->required = NO;
- parm.north->label =
- _("Northing coordinate (point of interest)");
- parm.north->description = _("Default: map center");
- parm.north->guisection = _("Position");
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
- flag1 = G_define_flag();
- flag1->key = 'z';
- flag1->description = _("Don't ignore zero elevation");
+ G_get_window(&window);
- flag2 = G_define_flag();
- flag2->key = 'v';
- flag2->description =
- _("Verbose output (also print out sun position etc.)");
+ /* require at least one output */
+ elev_name = parm.elev->answer;
+ azimuth_name = parm.azimuth->answer;
+ sunhour_name = parm.sunhours->answer;
+ if (!elev_name && !azimuth_name && !sunhour_name)
+ G_fatal_error(_("No output requested, exiting."));
- flag3 = G_define_flag();
- flag3->key = 's';
- flag3->description = _("Calculate sun position only and exit");
- flag3->guisection = _("Print");
-
- flag4 = G_define_flag();
- flag4->key = 'g';
- flag4->description =
- _("Print the sun position output in shell script style");
- flag4->guisection = _("Print");
+ sscanf(parm.year->answer, "%i", &year);
+ if (parm.month->answer)
+ sscanf(parm.month->answer, "%i", &month);
+ else
+ month = -1;
+ sscanf(parm.day->answer, "%i", &day);
+ sscanf(parm.hour->answer, "%i", &hour);
+ sscanf(parm.minutes->answer, "%i", &minutes);
+ sscanf(parm.seconds->answer, "%i", &seconds);
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ lst_time = (parm.gst_time->answer == 0);
+ use_solpos = parm.solpos->answer;
- zeros = flag1->answer;
+ /* init variables */
+ ha = 180;
+ ha_cos = 0;
+ sd_cos = 0;
+ sd_sin = 1;
+ north_gc_cos = 0;
+ north_gc_sin = 1;
+ se_sin = 0;
- G_get_window(&window);
+ if (use_solpos && lst_time) {
+ G_warning(_("NREL solpos algorithm uses Greenwich standard time."));
+ G_warning(_("Time will be interpreted as Greenwich standard time."));
- /* if not given, get east and north: XX */
- if (!parm.north->answer || !parm.east->answer) {
- north = (window.north - window.south) / 2. + window.south;
- east = (window.west - window.east) / 2. + window.east;
- G_verbose_message(_("Using map center coordinates: %f %f"), east, north);
+ lst_time = 0;
}
- else { /* user defined east, north: */
-
- sscanf(parm.north->answer, "%lf", &north);
- sscanf(parm.east->answer, "%lf", &east);
- if (strlen(parm.east->answer) == 0)
- G_fatal_error(_("Empty east coordinate specified"));
- if (strlen(parm.north->answer) == 0)
- G_fatal_error(_("Empty north coordinate specified"));
+ if (!use_solpos) {
+ if (lst_time)
+ G_message(_("Time will be interpreted as local sidereal time."));
+ else
+ G_message(_("Time will be interpreted as Greenwich standard time."));
+
+ if (sunhour_name)
+ G_fatal_error(_("Sunshine hours require solpos"));
}
- /* check which method to use for sun position:
- either user defines directly sun position or it is calculated */
+ if ((G_projection() != PROJECTION_LL)) {
+ if (window.proj == 0)
+ G_fatal_error(_("Current projection is x,y (undefined)."));
- if (parm.opt3->answer && parm.opt4->answer)
- solparms = 1; /* opt3 & opt4 complete */
- else
- solparms = 0; /* calculate sun position */
+ do_reproj = 1;
- if (parm.year->answer && parm.month->answer && parm.day->answer &&
- parm.hour->answer && parm.minutes->answer && parm.seconds->answer &&
- parm.timezone->answer)
- locparms = 1; /* complete */
- else
- locparms = 0;
+ /* read current projection info */
+ if ((in_proj_info = G_get_projinfo()) == NULL)
+ G_fatal_error(_("Can't get projection info of current location"));
- if (solparms && locparms) /* both defined */
- G_fatal_error(_("Either define sun position or location/date/time parameters"));
+ if ((in_unit_info = G_get_projunits()) == NULL)
+ G_fatal_error(_("Can't get projection units of current location"));
- if (!solparms && !locparms) /* nothing defined */
- G_fatal_error(_("Neither sun position nor east/north, date/time/timezone definition are complete"));
+ 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 here, one definition was complete */
- if (locparms) {
- G_message(_("Calculating sun position... (using solpos (V. %s) from NREL)"),
- SOLPOSVERSION);
- use_solpos = 1;
+ G_free_key_value(in_proj_info);
+ G_free_key_value(in_unit_info);
+
+ /* output projection to lat/long w/ same ellipsoid as input */
+ 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 update lat/long projection parameters"));
}
- else {
- G_message(_("Using user defined sun azimuth, altitude settings (ignoring eventual other values)"));
- use_solpos = 0;
- }
- name = parm.opt1->answer;
- outname = parm.opt2->answer;
- if (!use_solpos) {
- sscanf(parm.opt3->answer, "%lf", &dalti);
- sscanf(parm.opt4->answer, "%lf", &dazi);
+ /* always init pd */
+ S_init(&pd);
+ pd.function = S_GEOM;
+ if (use_solpos) {
+ pd.function = S_ZENETR;
+ if (azimuth_name)
+ pd.function = S_SOLAZM;
+ if (sunhour_name)
+ pd.function |= S_SRSS;
}
- else {
- sscanf(parm.year->answer, "%i", &year);
- sscanf(parm.month->answer, "%i", &month);
- sscanf(parm.day->answer, "%i", &day);
- sscanf(parm.hour->answer, "%i", &hour);
- sscanf(parm.minutes->answer, "%i", &minutes);
- sscanf(parm.seconds->answer, "%i", &seconds);
- sscanf(parm.timezone->answer, "%lf", &timezone);
- }
+ if (month == -1)
+ doy = day;
+ else
+ doy = dom2doy2(year, month, day);
+
+ set_solpos_time(&pd, year, 1, doy, hour, minutes, seconds);
+ set_solpos_longitude(&pd, 0);
+ pd.latitude = 0;
+ S_solpos(&pd);
- /* NOTES: G_calc_solar_position ()
- - the algorithm will compensate for leap year.
- - longitude, latitude: decimal degree
- - timezone: DO NOT ADJUST FOR DAYLIGHT SAVINGS TIME.
- - timezone: negative for zones west of Greenwich
- - lat/long: east and north positive
- - the atmospheric refraction is calculated for 1013hPa, 15�C
- - time: local time from your watch
+ if (lst_time) {
+ /* hour angle */
+ /***************************************************************
+ * The hour angle of a point on the Earth's surface is the angle
+ * through which the earth would turn to bring the meridian of
+ * the point directly under the sun. This angular displacement
+ * represents time (1 hour = 15 degrees).
+ * The hour angle is negative in the morning, zero at 12:00,
+ * and positive in the afternoon
+ ***************************************************************/
- Order of parameters:
- long, lat, timezone, year, month, day, hour, minutes, seconds
- */
+ ha = 15.0 * (hour + minutes / 60.0 + seconds / 3600.0) - 180.;
+ G_debug(1, "Solar hour angle, degrees: %.2f", ha);
+ ha *= DEG2RAD;
+ ha_cos = cos(ha);
+ roundoff(&ha_cos);
+ }
- if (use_solpos) {
- G_debug(3, "\nlat:%f long:%f", north, east);
- retval =
- calc_solar_position(east, north, timezone, year, month, day,
- hour, minutes, seconds);
+ if (!use_solpos) {
+ /* sun declination */
+ /***************************************************************
+ * The declination of the sun is the angle between
+ * the rays of the sun and the plane of the Earth's equator.
+ ***************************************************************/
- /* Remove +0.5 above if you want round-down instead of round-to-nearest */
- sretr = (int)floor(pdat->sretr); /* sunrise */
- dsretr = pdat->sretr;
- sretr_sec = (int)
- floor(((dsretr - floor(dsretr)) * 60 -
- floor((dsretr - floor(dsretr)) * 60)) * 60);
- ssetr = (int)floor(pdat->ssetr); /* sunset */
- dssetr = pdat->ssetr;
- ssetr_sec = (int)
- floor(((dssetr - floor(dssetr)) * 60 -
- floor((dssetr - floor(dssetr)) * 60)) * 60);
+ s_gamma = (2 * M_PI * (doy - 1)) / 365;
+ G_debug(1, "fractional year in radians: %.2f", s_gamma);
+ /* sun declination for day of year with Fourier series representation
+ * NOTE: based on 1950, override with solpos */
+ s_declination = (0.006918 - 0.399912 * cos(s_gamma) + 0.070257 * sin(s_gamma) -
+ 0.006758 * cos(2 * s_gamma) + 0.000907 * sin(2 * s_gamma) -
+ 0.002697 * cos(3 * s_gamma) + 0.00148 * sin(3 * s_gamma));
- /* print the results */
- if (retval == 0) { /* error check */
- if (flag2->answer || (flag3->answer && !flag2->answer)) {
- if (flag4->answer) {
- fprintf(stdout, "date=%d/%02d/%02d\n", pdat->year,
- pdat->month, pdat->day);
- fprintf(stdout, "daynum=%d\n", pdat->daynum);
- fprintf(stdout, "time=%02i:%02i:%02i\n", pdat->hour,
- pdat->minute, pdat->second);
- fprintf(stdout, "decimaltime=%f\n",
- pdat->hour + (pdat->minute * 100.0 / 60.0 +
- pdat->second * 100.0 / 3600.0) /
- 100.);
- fprintf(stdout, "longitudine=%f\n", pdat->longitude);
- fprintf(stdout, "latitude=%f\n", pdat->latitude);
- fprintf(stdout, "timezone=%f\n", pdat->timezone);
- fprintf(stdout, "sunazimuth=%f\n", pdat->azim);
- fprintf(stdout, "sunangleabovehorizon=%f\n",
- pdat->elevref);
+ G_debug(1, "sun declination: %.5f", s_declination * RAD2DEG);
+ G_debug(1, "sun declination (solpos): %.5f", pd.declin);
- if (sretr / 60 <= 24.0) {
- fprintf(stdout, "sunrise=%02d:%02d:%02d\n",
- sretr / 60, sretr % 60, sretr_sec);
- fprintf(stdout, "sunset=%02d:%02d:%02d\n", ssetr / 60,
- ssetr % 60, ssetr_sec);
- }
- }
- else {
- fprintf(stdout, "%d/%02d/%02d, daynum: %d, time: %02i:%02i:%02i (decimal time: %f)\n",
- pdat->year, pdat->month, pdat->day,
- pdat->daynum, pdat->hour, pdat->minute,
- pdat->second,
- pdat->hour + (pdat->minute * 100.0 / 60.0 +
- pdat->second * 100.0 / 3600.0) /
- 100.);
- fprintf(stdout, "long: %f, lat: %f, timezone: %f\n",
- pdat->longitude, pdat->latitude,
- pdat->timezone);
- fprintf(stdout, "Solar position: sun azimuth: %f, sun angle above horz. (refraction corrected): %f\n",
- pdat->azim, pdat->elevref);
-
- if (sretr / 60 <= 24.0) {
- fprintf(stdout, "Sunrise time (without refraction): %02d:%02d:%02d\n",
- sretr / 60, sretr % 60, sretr_sec);
- fprintf(stdout, "Sunset time (without refraction): %02d:%02d:%02d\n",
- ssetr / 60, ssetr % 60, ssetr_sec);
- }
- }
+ if (lst_time) {
+ north_ll = (window.north + window.south) / 2;
+ east_ll = (window.east + window.west) / 2;
+ if (do_reproj) {
+ if (pj_do_proj(&east_ll, &north_ll, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
}
- sunrise = pdat->sretr / 60.; /* decimal minutes */
- sunset = pdat->ssetr / 60.;
- current_time =
- pdat->hour + (pdat->minute / 60.) + (pdat->second / 3600.);
+ pd.timezone = east_ll / 15.;
+ pd.time_updated = 1;
+ set_solpos_longitude(&pd, east_ll);
+ G_debug(1, "fake timezone: %.2f", pd.timezone);
+ S_solpos(&pd);
+ G_debug(1, "Solar hour angle (solpos), degrees: %.2f", pd.hrang);
}
- else /* fatal error in G_calc_solar_position() */
- G_fatal_error(_("Please correct settings"));
- }
- if (use_solpos) {
- dalti = pdat->elevref;
- dazi = pdat->azim;
- } /* otherwise already defined */
+ /* always use solpos sun declination */
+ s_declination = pd.declin * DEG2RAD;
+ sd_sin = sin(s_declination);
+ roundoff(&sd_sin);
+ sd_cos = cos(s_declination);
+ roundoff(&sd_cos);
+ G_debug(1, "sun declination (solpos): %.5f", s_declination * RAD2DEG);
- /* check sunrise */
- if (use_solpos) {
- G_debug(3, "current_time:%f sunrise:%f", current_time, sunrise);
- if ((current_time < sunrise)) {
- if (sretr / 60 <= 24.0)
- G_message(_("Time (%02i:%02i:%02i) is before sunrise (%02d:%02d:%02d)"),
- pdat->hour, pdat->minute, pdat->second, sretr / 60,
- sretr % 60, sretr_sec);
- else
- G_message(_("Time (%02i:%02i:%02i) is before sunrise"),
- pdat->hour, pdat->minute, pdat->second);
-
- G_warning(_("Nothing to calculate. Please verify settings."));
+ if (0 && lst_time) {
+ pd.timezone = 0;
+ pd.time_updated = pd.longitude_updated = 1;
+ S_solpos(&pd);
}
- if ((current_time > sunset)) {
- if (sretr / 60 <= 24.0)
- G_message(_("Time (%02i:%02i:%02i) is after sunset (%02d:%02d:%02d)"),
- pdat->hour, pdat->minute, pdat->second, ssetr / 60,
- ssetr % 60, ssetr_sec);
- else
- G_message(_("Time (%02i:%02i:%02i) is after sunset"),
- pdat->hour, pdat->minute, pdat->second);
- G_warning(_("Nothing to calculate. Please verify settings."));
- }
}
- if (flag3->answer && (use_solpos == 1)) { /* we only want the sun position */
- exit(EXIT_SUCCESS);
+ if (elev_name) {
+ if ((elev_fd = G_open_raster_new(elev_name, FCELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), elev_name);
+
+ elevbuf = G_allocate_f_raster_buf();
}
- else if (flag3->answer && (use_solpos == 0)) {
- /* are you joking ? */
- G_message(_("You already know the sun position"));
- exit(EXIT_SUCCESS);
+ else {
+ elevbuf = NULL;
+ elev_fd = -1;
}
- if (!outname)
- G_fatal_error(_("Option <%s> required"), parm.opt2->key);
-
- /* Search for output layer in all mapsets ? yes. */
- mapset = G_find_cell2(name, "");
+ if (azimuth_name) {
+ if ((azimuth_fd = G_open_raster_new(azimuth_name, FCELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), azimuth_name);
- if ((elev_fd = G_open_cell_old(name, mapset)) < 0)
- G_fatal_error(_("Unable to open raster map <%s>"), name);
- if ((output_fd = G_open_cell_new(outname)) < 0)
- G_fatal_error(_("Unable to create raster map <%s>"), outname);
+ azimuthbuf = G_allocate_f_raster_buf();
+ }
+ else {
+ azimuthbuf = NULL;
+ azimuth_fd = -1;
+ }
- data_type = G_get_raster_map_type(elev_fd);
- elevbuf.v = G_allocate_raster_buf(data_type);
- tmpbuf.v = G_allocate_raster_buf(data_type);
- outbuf.v = G_allocate_raster_buf(CELL_TYPE); /* binary map */
+ if (sunhour_name) {
+ if ((sunhour_fd = G_open_raster_new(sunhour_name, FCELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), sunhour_name);
- if (data_type == CELL_TYPE) {
- if ((G_read_range(name, mapset, &range)) < 0)
- G_fatal_error(_("Unable to open range file for raster map <%s>"), name);
- G_get_range_min_max(&range, &min, &max);
- dmin = (double)min;
- dmax = (double)max;
+ sunhourbuf = G_allocate_f_raster_buf();
}
else {
- G_read_fp_range(name, mapset, &fprange);
- G_get_fp_range_min_max(&fprange, &dmin, &dmax);
+ sunhourbuf = NULL;
+ sunhour_fd = -1;
}
- azi = 2 * M_PI * dazi / 360;
- alti = 2 * M_PI * dalti / 360;
- nstep = cos(azi) * window.ns_res;
- estep = sin(azi) * window.ew_res;
- row1 = 0;
+ if (elev_name && azimuth_name) {
+ G_message(_("Calculating solar elevation and azimuth..."));
+ }
+ else if (elev_name) {
+ G_message(_("Calculating solar elevation..."));
+ }
+ else if (azimuth_name) {
+ G_message(_("Calculating solar azimuth..."));
+ }
- G_message(_("Calculating shadows from DEM..."));
- while (row1 < window.rows) {
- G_percent(row1, window.rows, 2);
- col1 = 0;
- drow = -1;
- if (G_get_raster_row(elev_fd, elevbuf.v, row1, data_type) < 0)
- G_fatal_error(_("Unable to read raster map <%s> row %d"),
- name, row1);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
- while (col1 < window.cols) {
- dvalue = raster_value(elevbuf, data_type, col1);
- /* outbuf.c[col1]=1; */
- G_set_null_value(&outbuf.c[col1], 1, CELL_TYPE);
- OK = 1;
- east = G_col_to_easting(col1 + 0.5, &window);
- north = G_row_to_northing(row1 + 0.5, &window);
- east1 = east;
- north1 = north;
- if (dvalue == 0.0 && !zeros)
- OK = 0;
- while (OK == 1) {
- east += estep;
- north += nstep;
- if (north > window.north || north < window.south
- || east > window.east || east < window.west)
- OK = 0;
- else {
- maxh = tan(alti) *
- sqrt((north1 - north) * (north1 - north) +
- (east1 - east) * (east1 - east));
- if ((maxh) > (dmax - dvalue))
- OK = 0;
- else {
- dcol = G_easting_to_col(east, &window);
- if (drow != G_northing_to_row(north, &window)) {
- drow = G_northing_to_row(north, &window);
- G_get_raster_row(elev_fd, tmpbuf.v, (int)drow,
- data_type);
- }
- dvalue2 = raster_value(tmpbuf, data_type, (int)dcol);
- if ((dvalue2 - dvalue) > (maxh)) {
- OK = 0;
- outbuf.c[col1] = 1;
- }
- }
+ ba2 = 6356752.3142 / 6378137.0;
+ ba2 = ba2 * ba2;
+
+ for (row = 0; row < nrows; row++) {
+
+ G_percent(row, nrows, 2);
+
+ /* get cell center northing */
+ north = window.north - (row + 0.5) * window.ns_res;
+ north_ll = north;
+
+ for (col = 0; col < ncols; col++) {
+ long int retval;
+ s_elevation = s_azimuth = -1.;
+
+ /* get cell center easting */
+ east = window.west + (col + 0.5) * window.ew_res;
+ east_ll = east;
+
+ if (do_reproj) {
+ north_ll = north;
+
+ if (pj_do_proj(&east_ll, &north_ll, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+ }
+
+ /* geocentric latitude */
+ north_gc = atan(ba2 * tan(DEG2RAD * north_ll));
+ north_gc_sin = sin(north_gc);
+ roundoff(&north_gc_sin);
+ north_gc_cos = cos(north_gc);
+ roundoff(&north_gc_cos);
+
+ if (!lst_time) {
+ set_solpos_longitude(&pd, east_ll);
+ pd.latitude = north_gc * RAD2DEG;
+ retval = S_solpos(&pd);
+ S_decode(retval, &pd);
+ G_debug(3, "solpos hour angle: %.5f", pd.hrang);
+ }
+
+ /* solar elevation angle */
+ if (!use_solpos) {
+ if (!lst_time) {
+ ha = pd.hrang;
+ ha_cos = cos(ha * DEG2RAD);
+ roundoff(&ha_cos);
}
+ se_sin = ha_cos * sd_cos * north_gc_cos + sd_sin * north_gc_sin;
+ roundoff(&se_sin);
+ s_elevation = RAD2DEG * asin(se_sin);
}
- G_debug(3, "Analysing col %i", col1);
- col1 += 1;
+ else /* use_solpos && !lst_time */
+ s_elevation = pd.elevetr;
+
+ if (elev_name)
+ elevbuf[col] = s_elevation;
+
+ if (azimuth_name) {
+ /* solar azimuth angle */
+ if (!use_solpos) {
+ sa_cos = (se_sin * north_gc_sin - sd_sin) /
+ (cos(DEG2RAD * s_elevation) * north_gc_cos);
+ roundoff(&sa_cos);
+ s_azimuth = RAD2DEG * acos(sa_cos);
+
+ /* morning */
+ s_azimuth = 180. - RAD2DEG * acos(sa_cos);
+ if (ha > 0) /* afternoon */
+ s_azimuth = 360.0 - s_azimuth;
+ }
+ else
+ s_azimuth = pd.azim;
+
+ azimuthbuf[col] = s_azimuth;
+ }
+
+ if (sunhour_name)
+ sunhourbuf[col] = (pd.ssetr - pd.sretr) / 60.;
+
}
- G_debug(3, "Writing result row %i of %i", row1, window.rows);
- G_put_raster_row(output_fd, outbuf.c, CELL_TYPE);
- row1 += 1;
+ if (elev_name)
+ G_put_f_raster_row(elev_fd, elevbuf);
+ if (azimuth_name)
+ G_put_f_raster_row(azimuth_fd, azimuthbuf);
+ if (sunhour_name)
+ G_put_f_raster_row(sunhour_fd, sunhourbuf);
}
- G_percent(1, 1, 1);
+ G_percent(1, 1, 2);
- G_close_cell(output_fd);
- G_close_cell(elev_fd);
+ if (elev_name) {
+ G_close_cell(elev_fd);
+ /* writing history file */
+ G_short_history(elev_name, "raster", &hist);
+ G_command_history(&hist);
+ G_write_history(elev_name, &hist);
+ }
+ if (azimuth_name) {
+ G_close_cell(azimuth_fd);
+ /* writing history file */
+ G_short_history(azimuth_name, "raster", &hist);
+ G_command_history(&hist);
+ G_write_history(azimuth_name, &hist);
+ }
+ if (sunhour_name) {
+ G_close_cell(sunhour_fd);
+ /* writing history file */
+ G_short_history(sunhour_name, "raster", &hist);
+ G_command_history(&hist);
+ G_write_history(sunhour_name, &hist);
+ }
- /* writing history file */
- G_short_history(outname, "raster", &hist);
- G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation map: %s", name);
- G_command_history(&hist);
- G_write_history(outname, &hist);
+ G_done_msg(" ");
exit(EXIT_SUCCESS);
}
-#ifdef RASTER_VALUE_FUNC
-double raster_value(union RASTER_PTR buf, int data_type, int col)
+void set_solpos_time(struct posdata *pdat, int year, int month, int day,
+ int hour, int minute, int second)
{
- double retval;
+ pdat->year = year;
+ pdat->month = month;
+ pdat->day = day;
+ pdat->daynum = day;
+ pdat->hour = hour;
+ pdat->minute = minute;
+ pdat->second = second;
+ pdat->timezone = 0;
- switch (data_type) {
- case CELL_TYPE:
- retval = (double)buf.c[col];
- break;
- case FCELL_TYPE:
- retval = (double)buf.f[col];
- break;
- case DCELL_TYPE:
- retval = (double)buf.d[col];
- break;
+ pdat->time_updated = 1;
+ pdat->longitude_updated = 1;
+}
+
+void set_solpos_longitude(struct posdata *pdat, double longitude)
+{
+ pdat->longitude = longitude;
+
+ pdat->longitude_updated = 1;
+}
+
+int roundoff(double *x)
+{
+ /* watch out for the roundoff errors */
+ if (fabs(*x) > 1.0) {
+ if (*x > 0.0)
+ *x = 1.0;
+ else
+ *x = -1.0;
+
+ return 1;
}
- return retval;
+ return 0;
}
-#endif
Modified: grass-addons/grass6/raster/r.sun.angle/solpos00.c
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/solpos00.c 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/solpos00.c 2012-08-24 18:13:55 UTC (rev 52877)
@@ -127,8 +127,6 @@
/* cumulative number of days prior to beginning of month */
-static float degrad = 57.295779513; /* converts from radians to degrees */
-static float raddeg = 0.0174532925; /* converts from degrees to radians */
/*============================================================================
* Local function prototypes
@@ -289,6 +287,9 @@
pdat->sbrad = 31.7; /* Eppley shadow band radius */
pdat->sbsky = 0.04; /* Drummond factor for partly cloudy skies */
pdat->function = S_ALL; /* compute all parameters */
+
+ pdat->time_updated = 0;
+ pdat->longitude_updated = 0;
}
@@ -392,6 +393,19 @@
}
+int dom2doy2(int year, int month, int day)
+{
+ int doy = day + month_days[0][month];
+
+ /* (adjust for leap year) */
+ if (((year % 4) == 0) &&
+ (((year % 100) != 0) || ((year % 400) == 0)) &&
+ (month > 2))
+ doy += 1;
+
+ return doy;
+}
+
/*============================================================================
* Local void function doy2dom
*
@@ -447,147 +461,178 @@
float top; /* numerator (top) of the fraction */
int leap; /* leap year counter */
- /* Day angle */
- /* Iqbal, M. 1983. An Introduction to Solar Radiation.
- Academic Press, NY., page 3 */
- pdat->dayang = 360.0 * (pdat->daynum - 1) / 365.0;
+ /** time-dependent variables **/
+ if (pdat->time_updated) {
+ /* Day angle */
+ /* Iqbal, M. 1983. An Introduction to Solar Radiation.
+ Academic Press, NY., page 3 */
+ pdat->dayang = 360.0 * (pdat->daynum - 1) / 365.0;
- /* Earth radius vector * solar constant = solar energy */
- /* Spencer, J. W. 1971. Fourier series representation of the
- position of the sun. Search 2 (5), page 172 */
- sd = sin(raddeg * pdat->dayang);
- cd = cos(raddeg * pdat->dayang);
- d2 = 2.0 * pdat->dayang;
- c2 = cos(raddeg * d2);
- s2 = sin(raddeg * d2);
+ /* Earth radius vector * solar constant = solar energy */
+ /* Spencer, J. W. 1971. Fourier series representation of the
+ position of the sun. Search 2 (5), page 172 */
+ sd = sin(DEG2RAD * pdat->dayang);
+ cd = cos(DEG2RAD * pdat->dayang);
+ d2 = 2.0 * pdat->dayang;
+ c2 = cos(DEG2RAD * d2);
+ s2 = sin(DEG2RAD * d2);
- pdat->erv = 1.000110 + 0.034221 * cd + 0.001280 * sd;
- pdat->erv += 0.000719 * c2 + 0.000077 * s2;
+ pdat->erv = 1.000110 + 0.034221 * cd + 0.001280 * sd;
+ pdat->erv += 0.000719 * c2 + 0.000077 * s2;
- /* Universal Coordinated (Greenwich standard) time */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->utime =
- pdat->hour * 3600.0 +
- pdat->minute * 60.0 + pdat->second - (float)pdat->interval / 2.0;
- pdat->utime = pdat->utime / 3600.0 - pdat->timezone;
+ /* Universal Coordinated (Greenwich standard) time */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->utime =
+ pdat->hour * 3600.0 +
+ pdat->minute * 60.0 + pdat->second - (float)pdat->interval / 2.0;
+ pdat->utime = pdat->utime / 3600.0 - pdat->timezone;
- /* Julian Day minus 2,400,000 days (to eliminate roundoff errors) */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
+ /* Julian Day minus 2,400,000 days (to eliminate roundoff errors) */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
- /* No adjustment for century non-leap years since this function is
- bounded by 1950 - 2050 */
- delta = pdat->year - 1949;
- leap = (int)(delta / 4.0);
- pdat->julday =
- 32916.5 + delta * 365.0 + leap + pdat->daynum + pdat->utime / 24.0;
+ /* No adjustment for century non-leap years since this function is
+ bounded by 1950 - 2050 */
+ delta = pdat->year - 1949;
+ leap = (int)(delta / 4.0);
+ pdat->julday =
+ 32916.5 + delta * 365.0 + leap + pdat->daynum + pdat->utime / 24.0;
- /* Time used in the calculation of ecliptic coordinates */
- /* Noon 1 JAN 2000 = 2,400,000 + 51,545 days Julian Date */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->ectime = pdat->julday - 51545.0;
+ /* Time used in the calculation of ecliptic coordinates */
+ /* Noon 1 JAN 2000 = 2,400,000 + 51,545 days Julian Date */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->ectime = pdat->julday - 51545.0;
- /* Mean longitude */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->mnlong = 280.460 + 0.9856474 * pdat->ectime;
+ /* Mean longitude */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->mnlong = 280.460 + 0.9856474 * pdat->ectime;
- /* (dump the multiples of 360, so the answer is between 0 and 360) */
- pdat->mnlong -= 360.0 * (int)(pdat->mnlong / 360.0);
- if (pdat->mnlong < 0.0)
- pdat->mnlong += 360.0;
+ /* (dump the multiples of 360, so the answer is between 0 and 360) */
+ pdat->mnlong -= 360.0 * (int)(pdat->mnlong / 360.0);
+ if (pdat->mnlong < 0.0)
+ pdat->mnlong += 360.0;
- /* Mean anomaly */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->mnanom = 357.528 + 0.9856003 * pdat->ectime;
+ /* Mean anomaly */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->mnanom = 357.528 + 0.9856003 * pdat->ectime;
- /* (dump the multiples of 360, so the answer is between 0 and 360) */
- pdat->mnanom -= 360.0 * (int)(pdat->mnanom / 360.0);
- if (pdat->mnanom < 0.0)
- pdat->mnanom += 360.0;
+ /* (dump the multiples of 360, so the answer is between 0 and 360) */
+ pdat->mnanom -= 360.0 * (int)(pdat->mnanom / 360.0);
+ if (pdat->mnanom < 0.0)
+ pdat->mnanom += 360.0;
- /* Ecliptic longitude */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->eclong = pdat->mnlong + 1.915 * sin(pdat->mnanom * raddeg) +
- 0.020 * sin(2.0 * pdat->mnanom * raddeg);
+ /* Ecliptic longitude */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->eclong = pdat->mnlong + 1.915 * sin(pdat->mnanom * DEG2RAD) +
+ 0.020 * sin(2.0 * pdat->mnanom * DEG2RAD);
- /* (dump the multiples of 360, so the answer is between 0 and 360) */
- pdat->eclong -= 360.0 * (int)(pdat->eclong / 360.0);
- if (pdat->eclong < 0.0)
- pdat->eclong += 360.0;
+ /* (dump the multiples of 360, so the answer is between 0 and 360) */
+ pdat->eclong -= 360.0 * (int)(pdat->eclong / 360.0);
+ if (pdat->eclong < 0.0)
+ pdat->eclong += 360.0;
- /* Obliquity of the ecliptic */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
+ /* Obliquity of the ecliptic */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
- /* 02 Feb 2001 SMW corrected sign in the following line */
- /* pdat->ecobli = 23.439 + 4.0e-07 * pdat->ectime; */
- pdat->ecobli = 23.439 - 4.0e-07 * pdat->ectime;
+ /* 02 Feb 2001 SMW corrected sign in the following line */
+ /* pdat->ecobli = 23.439 + 4.0e-07 * pdat->ectime; */
+ /* more precisely: */
+ /* pdat->ecobli = 23.43929 - 3.56237e-07 * pdat->ectime; */
+ pdat->ecobli = 23.439 - 4.0e-07 * pdat->ectime;
+ G_debug(1, "axial tilt (solpos): %.5f", pdat->ecobli);
- /* Declination */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->declin = degrad * asin(sin(pdat->ecobli * raddeg) *
- sin(pdat->eclong * raddeg));
+ if (pdat->ectime != 0) {
+ /* convert ectime from days to years */
+ double ectime = pdat->ectime / 36525.0;
+ double ectime2 = ectime * ectime;
+ double ectime3 = ectime2 * ectime;
- /* Right ascension */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- top = cos(raddeg * pdat->ecobli) * sin(raddeg * pdat->eclong);
- bottom = cos(raddeg * pdat->eclong);
+ /* recommended by the International Astronomical Union in 2000: */
+ /* accurate for the year 2000 +- 5000 years
+ * result is in arc seconds */
+ pdat->ecobli = 84381.448 - 46.84024 * ectime - 5.9e-5 * ectime2 + 1.813e-3 * ectime3;
+ /* convert to degrees */
+ pdat->ecobli /= 3600.0;
- pdat->rascen = degrad * atan2(top, bottom);
+ G_debug(1, "axial tilt (IAU 2000): %.5f", pdat->ecobli);
+ }
+ else
+ pdat->ecobli = 23.43929;
- /* (make it a positive angle) */
- if (pdat->rascen < 0.0)
- pdat->rascen += 360.0;
+ /* Declination */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->declin = RAD2DEG * asin(sin(pdat->ecobli * DEG2RAD) *
+ sin(pdat->eclong * DEG2RAD));
- /* Greenwich mean sidereal time */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->gmst = 6.697375 + 0.0657098242 * pdat->ectime + pdat->utime;
+ /* Right ascension */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ top = cos(DEG2RAD * pdat->ecobli) * sin(DEG2RAD * pdat->eclong);
+ bottom = cos(DEG2RAD * pdat->eclong);
- /* (dump the multiples of 24, so the answer is between 0 and 24) */
- pdat->gmst -= 24.0 * (int)(pdat->gmst / 24.0);
- if (pdat->gmst < 0.0)
- pdat->gmst += 24.0;
+ pdat->rascen = RAD2DEG * atan2(top, bottom);
- /* Local mean sidereal time */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->lmst = pdat->gmst * 15.0 + pdat->longitude;
+ /* (make it a positive angle) */
+ if (pdat->rascen < 0.0)
+ pdat->rascen += 360.0;
- /* (dump the multiples of 360, so the answer is between 0 and 360) */
- pdat->lmst -= 360.0 * (int)(pdat->lmst / 360.0);
- if (pdat->lmst < 0.)
- pdat->lmst += 360.0;
+ /* Greenwich mean sidereal time */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->gmst = 6.697375 + 0.0657098242 * pdat->ectime + pdat->utime;
- /* Hour angle */
- /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
- approximate solar position (1950-2050). Solar Energy 40 (3),
- pp. 227-235. */
- pdat->hrang = pdat->lmst - pdat->rascen;
+ /* (dump the multiples of 24, so the answer is between 0 and 24) */
+ pdat->gmst -= 24.0 * (int)(pdat->gmst / 24.0);
+ if (pdat->gmst < 0.0)
+ pdat->gmst += 24.0;
- /* (force it between -180 and 180 degrees) */
- if (pdat->hrang < -180.0)
- pdat->hrang += 360.0;
- else if (pdat->hrang > 180.0)
- pdat->hrang -= 360.0;
+ pdat->time_updated = 0;
+ }
+
+ /** longitude-dependent variables **/
+ if (pdat->longitude_updated) {
+ /* Local mean sidereal time */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->lmst = pdat->gmst * 15.0 + pdat->longitude;
+
+ /* (dump the multiples of 360, so the answer is between 0 and 360) */
+ pdat->lmst -= 360.0 * (int)(pdat->lmst / 360.0);
+ if (pdat->lmst < 0.)
+ pdat->lmst += 360.0;
+
+ /* Hour angle */
+ /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for
+ approximate solar position (1950-2050). Solar Energy 40 (3),
+ pp. 227-235. */
+ pdat->hrang = pdat->lmst - pdat->rascen;
+
+ /* (force it between -180 and 180 degrees) */
+ if (pdat->hrang < -180.0)
+ pdat->hrang += 360.0;
+ else if (pdat->hrang > 180.0)
+ pdat->hrang -= 360.0;
+
+ pdat->longitude_updated = 0;
+ }
}
@@ -613,7 +658,7 @@
cz = -1.0;
}
- pdat->zenetr = acos(cz) * degrad;
+ pdat->zenetr = acos(cz) * RAD2DEG;
/* (limit the degrees below the horizon to 9 [+90 -> 99]) */
if (pdat->zenetr > 99.0)
@@ -647,7 +692,7 @@
else if (cssha > 1.0)
pdat->ssha = 0.0;
else
- pdat->ssha = degrad * acos(cssha);
+ pdat->ssha = RAD2DEG * acos(cssha);
}
else if (((pdat->declin >= 0.0) && (pdat->latitude > 0.0)) ||
((pdat->declin < 0.0) && (pdat->latitude < 0.0)))
@@ -670,8 +715,8 @@
localtrig(pdat, tdat);
p = 0.6366198 * pdat->sbwid / pdat->sbrad * pow(tdat->cd, 3);
- t1 = tdat->sl * tdat->sd * pdat->ssha * raddeg;
- t2 = tdat->cl * tdat->cd * sin(pdat->ssha * raddeg);
+ t1 = tdat->sl * tdat->sd * pdat->ssha * DEG2RAD;
+ t2 = tdat->cl * tdat->cd * sin(pdat->ssha * DEG2RAD);
pdat->sbcf = pdat->sbsky + 1.0 / (1.0 - p * (t1 + t2));
}
@@ -739,8 +784,8 @@
float se; /* sine of the solar elevation */
localtrig(pdat, tdat);
- ce = cos(raddeg * pdat->elevetr);
- se = sin(raddeg * pdat->elevetr);
+ ce = cos(DEG2RAD * pdat->elevetr);
+ se = sin(DEG2RAD * pdat->elevetr);
pdat->azim = 180.0;
cecl = ce * tdat->cl;
@@ -751,7 +796,7 @@
else if (ca < -1.0)
ca = -1.0;
- pdat->azim = 180.0 - acos(ca) * degrad;
+ pdat->azim = 180.0 - acos(ca) * RAD2DEG;
if (pdat->hrang > 0)
pdat->azim = 360.0 - pdat->azim;
}
@@ -779,7 +824,7 @@
/* Otherwise, we have refraction */
else {
- tanelev = tan(raddeg * pdat->elevetr);
+ tanelev = tan(DEG2RAD * pdat->elevetr);
if (pdat->elevetr >= 5.0)
refcor = 58.1 / tanelev -
0.07 / (pow(tanelev, 3)) + 0.000086 / (pow(tanelev, 5));
@@ -806,7 +851,7 @@
/* Refracted solar zenith angle */
pdat->zenref = 90.0 - pdat->elevref;
- pdat->coszen = cos(raddeg * pdat->zenref);
+ pdat->coszen = cos(DEG2RAD * pdat->zenref);
}
@@ -826,7 +871,7 @@
}
else {
pdat->amass =
- 1.0 / (cos(raddeg * pdat->zenref) + 0.50572 *
+ 1.0 / (cos(DEG2RAD * pdat->zenref) + 0.50572 *
pow((96.07995 - pdat->zenref), -1.6364));
pdat->ampress = pdat->amass * pdat->press / 1013.0;
@@ -886,15 +931,15 @@
if (tdat->sd < -900.0) { /* sd was initialized -999 as flag */
tdat->sd = 1.0; /* reflag as having completed calculations */
if (pdat->function | CD_MASK)
- tdat->cd = cos(raddeg * pdat->declin);
+ tdat->cd = cos(DEG2RAD * pdat->declin);
if (pdat->function | CH_MASK)
- tdat->ch = cos(raddeg * pdat->hrang);
+ tdat->ch = cos(DEG2RAD * pdat->hrang);
if (pdat->function | CL_MASK)
- tdat->cl = cos(raddeg * pdat->latitude);
+ tdat->cl = cos(DEG2RAD * pdat->latitude);
if (pdat->function | SD_MASK)
- tdat->sd = sin(raddeg * pdat->declin);
+ tdat->sd = sin(DEG2RAD * pdat->declin);
if (pdat->function | SL_MASK)
- tdat->sl = sin(raddeg * pdat->latitude);
+ tdat->sl = sin(DEG2RAD * pdat->latitude);
}
}
@@ -917,13 +962,13 @@
/* Cosine of the angle between the sun and a tipped flat surface,
useful for calculating solar energy on tilted surfaces */
- ca = cos(raddeg * pdat->azim);
- cp = cos(raddeg * pdat->aspect);
- ct = cos(raddeg * pdat->tilt);
- sa = sin(raddeg * pdat->azim);
- sp = sin(raddeg * pdat->aspect);
- st = sin(raddeg * pdat->tilt);
- sz = sin(raddeg * pdat->zenref);
+ ca = cos(DEG2RAD * pdat->azim);
+ cp = cos(DEG2RAD * pdat->aspect);
+ ct = cos(DEG2RAD * pdat->tilt);
+ sa = sin(DEG2RAD * pdat->azim);
+ sp = sin(DEG2RAD * pdat->aspect);
+ st = sin(DEG2RAD * pdat->tilt);
+ sz = sin(DEG2RAD * pdat->zenref);
pdat->cosinc = pdat->coszen * ct + sz * st * (ca * cp + sa * sp);
if (pdat->cosinc > 0.0)
Modified: grass-addons/grass6/raster/r.sun.angle/solpos00.h
===================================================================
--- grass-addons/grass6/raster/r.sun.angle/solpos00.h 2012-08-24 17:55:38 UTC (rev 52876)
+++ grass-addons/grass6/raster/r.sun.angle/solpos00.h 2012-08-24 18:13:55 UTC (rev 52877)
@@ -1,4 +1,9 @@
+#define RAD2DEG 57.295779513 /* converts from radians to degrees */
+#define DEG2RAD 0.0174532925 /* converts from degrees to radians */
+
+int dom2doy2(int year, int month, int day);
+
/*============================================================================
*
* NAME: solpos.h
@@ -178,6 +183,9 @@
int year; /* I: 4-digit year (2-digit year is NOT
allowed */
+ int time_updated; /* recalculate time-dependent variables */
+ int longitude_updated; /* recalculate longitude-dependent variables */
+
/***** FLOATS *****/
float amass; /* O: S_AMASS Relative optical airmass */
More information about the grass-commit
mailing list