Selaa lähdekoodia

move r.sunhours to trunk

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@56152 15284696-431f-4ddb-bdfa-cd5b030d7da7
Markus Metz 12 vuotta sitten
vanhempi
commit
8cf7246f7a

+ 12 - 0
raster/r.sunhours/Makefile

@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.sunhours
+
+LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(PROJLIB)
+DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP)
+
+EXTRA_INC = $(PROJINC)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

+ 523 - 0
raster/r.sunhours/main.c

@@ -0,0 +1,523 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.sunhours
+ * AUTHOR(S):    Markus Metz
+ * PURPOSE:      Calculates solar azimuth and angle, and 
+ *               sunshine hours (also called daytime period)
+ * COPYRIGHT:    (C) 2010-2013 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.
+ *
+ *****************************************************************************/
+
+/* TODO: always use solpos if time is Greenwich standard time */
+
+#define MAIN
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+#include "solpos00.h"
+
+void set_solpos_time(struct posdata *, int, int, int, int, int, int);
+void set_solpos_longitude(struct posdata *, double );
+int roundoff(double *);
+
+int main(int argc, char *argv[])
+{
+    struct GModule *module;
+    struct
+    {
+	struct Option *elev, *azimuth, *sunhours, *year,
+	    *month, *day, *hour, *minutes, *seconds;
+	struct Flag *lst_time, *no_solpos;
+    } parm;
+    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;
+    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();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("solar"));
+    module->label = _("Calculates solar elevation, solar azimuth, and sun hours.");
+    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.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.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.year = G_define_option();
+    parm.year->key = "year";
+    parm.year->type = TYPE_INTEGER;
+    parm.year->required = YES;
+    parm.year->description = _("Year");
+    parm.year->options = "1950-2050";
+    parm.year->guisection = _("Time");
+
+    parm.month = G_define_option();
+    parm.month->key = "month";
+    parm.month->type = TYPE_INTEGER;
+    parm.month->required = NO;
+    parm.month->label = _("Month");
+    parm.month->description = _("If not given, day is interpreted as day of the year");
+    parm.month->options = "1-12";
+    parm.month->guisection = _("Time");
+
+    parm.day = G_define_option();
+    parm.day->key = "day";
+    parm.day->type = TYPE_INTEGER;
+    parm.day->required = YES;
+    parm.day->description = _("Day");
+    parm.day->options = "1-366";
+    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");
+    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");
+    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");
+    parm.seconds->options = "0-60";
+    parm.seconds->answer = "0";
+    parm.seconds->guisection = _("Time");
+
+    parm.lst_time = G_define_flag();
+    parm.lst_time->key = 't';
+    parm.lst_time->description = _("Time is local sidereal time, not Greenwich standard time");
+
+    parm.no_solpos = G_define_flag();
+    parm.no_solpos->key = 's';
+    parm.no_solpos->description = _("Do not use solpos algorithm of NREL");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    G_get_window(&window);
+
+    /* 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."));
+
+    year = atoi(parm.year->answer);
+    if (parm.month->answer)
+	month = atoi(parm.month->answer);
+    else
+	month = -1;
+
+    day = atoi(parm.day->answer);
+    hour = atoi(parm.hour->answer);
+    minutes = atoi(parm.minutes->answer);
+    seconds = atoi(parm.seconds->answer);
+
+    lst_time = (parm.lst_time->answer != 0);
+    use_solpos = (parm.no_solpos->answer == 0);
+
+    /* init variables */
+    ha = 180;
+    ha_cos = 0;
+    sd_cos = 0;
+    sd_sin = 1;
+    north_gc_cos = 0;
+    north_gc_sin = 1;
+    se_sin = 0;
+
+    if (use_solpos && lst_time) {
+	G_warning(_("NREL solpos algorithm uses Greenwich standard time."));
+	G_warning(_("Time will be interpreted as Greenwich standard time."));
+
+	lst_time = 0;
+    }
+    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 NREL solpos."));
+    }
+
+    if ((G_projection() != PROJECTION_LL)) {
+	if (window.proj == 0)
+	    G_fatal_error(_("Current projection is x,y (undefined)."));
+
+	do_reproj = 1;
+
+	/* read current projection info */
+	if ((in_proj_info = G_get_projinfo()) == NULL)
+	    G_fatal_error(_("Cannot get projection info of current location"));
+
+	if ((in_unit_info = G_get_projunits()) == NULL)
+	    G_fatal_error(_("Cannot get projection units of current location"));
+
+	if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
+	    G_fatal_error(_("Cannot get projection key values of current location"));
+
+	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"));
+    }
+
+    /* 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;
+    }
+    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);
+
+    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
+	 ***************************************************************/
+
+	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) {
+	/* sun declination */
+	/***************************************************************
+	 * The declination of the sun is the angle between
+	 * the rays of the sun and the plane of the Earth's equator.
+	 ***************************************************************/
+
+	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));
+
+	G_debug(1, "sun declination: %.5f", s_declination * RAD2DEG);
+	G_debug(1, "sun declination (solpos): %.5f", pd.declin);
+
+	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)"));
+	    }
+	    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);
+	}
+
+	/* 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);
+
+	if (0 && lst_time) {
+	    pd.timezone = 0;
+	    pd.time_updated = pd.longitude_updated = 1;
+	    S_solpos(&pd);
+	}
+    }
+
+    if (elev_name) {
+	if ((elev_fd = Rast_open_new(elev_name, FCELL_TYPE)) < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), elev_name);
+
+	elevbuf = Rast_allocate_f_buf();
+    }
+    else {
+	elevbuf = NULL;
+	elev_fd = -1;
+    }
+
+    if (azimuth_name) {
+	if ((azimuth_fd = Rast_open_new(azimuth_name, FCELL_TYPE)) < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), azimuth_name);
+
+	azimuthbuf = Rast_allocate_f_buf();
+    }
+    else {
+	azimuthbuf = NULL;
+	azimuth_fd = -1;
+    }
+
+    if (sunhour_name) {
+	if ((sunhour_fd = Rast_open_new(sunhour_name, FCELL_TYPE)) < 0)
+	    G_fatal_error(_("Unable to create raster map <%s>"), sunhour_name);
+
+	sunhourbuf = Rast_allocate_f_buf();
+    }
+    else {
+	sunhourbuf = NULL;
+	sunhour_fd = -1;
+    }
+
+    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..."));
+    }
+
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    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);
+	    }
+	    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.;
+		if (sunhourbuf[col] > 24.)
+		    sunhourbuf[col] = 24.;
+		if (sunhourbuf[col] < 0.)
+		    sunhourbuf[col] = 0.;
+	    }
+
+	}
+	if (elev_name)
+	    Rast_put_f_row(elev_fd, elevbuf);
+	if (azimuth_name)
+	    Rast_put_f_row(azimuth_fd, azimuthbuf);
+	if (sunhour_name)
+	    Rast_put_f_row(sunhour_fd, sunhourbuf);
+    }
+    G_percent(1, 1, 2);
+
+    if (elev_name) {
+	Rast_close(elev_fd);
+	/* writing history file */
+	Rast_short_history(elev_name, "raster", &hist);
+	Rast_command_history(&hist);
+	Rast_write_history(elev_name, &hist);
+    }
+    if (azimuth_name) {
+	Rast_close(azimuth_fd);
+	/* writing history file */
+	Rast_short_history(azimuth_name, "raster", &hist);
+	Rast_command_history(&hist);
+	Rast_write_history(azimuth_name, &hist);
+    }
+    if (sunhour_name) {
+	Rast_close(sunhour_fd);
+	/* writing history file */
+	Rast_short_history(sunhour_name, "raster", &hist);
+	Rast_command_history(&hist);
+	Rast_write_history(sunhour_name, &hist);
+    }
+
+    G_done_msg(" ");
+
+    exit(EXIT_SUCCESS);
+}
+
+void set_solpos_time(struct posdata *pdat, int year, int month, int day,
+                    int hour, int minute, int second)
+{
+    pdat->year = year; 
+    pdat->month = month; 
+    pdat->day = day; 
+    pdat->daynum = day; 
+    pdat->hour = hour; 
+    pdat->minute = minute; 
+    pdat->second = second;
+    pdat->timezone = 0;
+
+    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 0;
+}

+ 42 - 0
raster/r.sunhours/r.sunhours.html

@@ -0,0 +1,42 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.sunhours</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>
+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>
+The sun azimuth angle is here defined as the azimuth angle in degrees
+of the sun from due north in a clockwise direction.
+<p>
+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).
+
+<h2>EXAMPLES</h2>
+
+Calculate sun's elevation angle map for 2010-10-11 at 14:00h solar time:
+<div class="code"><pre>
+r.sunhours elevation=sun_elev year=2010 month=10 day=11 hour=14 minute=00
+</pre></div>
+<p>
+
+Calculate photoperiod of day of year 1 (1st January) 2012
+<div class="code"><pre>
+r.sunhours -s sunhour=photoperiod_doy_001 year=2012 day=1
+</pre></div>
+
+<h2>AUTHOR</h2>
+
+Markus Metz
+
+<p><i>Last changed: $Date$</i>

Tiedoston diff-näkymää rajattu, sillä se on liian suuri
+ 1040 - 0
raster/r.sunhours/solpos00.c


+ 376 - 0
raster/r.sunhours/solpos00.h

@@ -0,0 +1,376 @@
+
+#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
+*
+*    Contains:
+*        S_solpos     (computes the solar position and intensity
+*                      from time and place)
+*            INPUTS:     (from posdata)
+*                          year, month, day, hour, minute, second,
+*                          latitude, longitude, timezone, interval
+*            OPTIONAL:   (from posdata; defaults from S_init function)
+*                            press   DEFAULT 1013.0 (standard pressure)
+*                            temp    DEFAULT   10.0 (standard temperature)
+*                            tilt    DEFAULT    0.0 (horizontal panel)
+*                            aspect  DEFAULT  180.0 (South-facing panel)
+*                            sbwid   DEFAULT    7.6 (shadowband width)
+*                            sbrad   DEFAULT   31.7 (shadowband radius)
+*                            sbsky   DEFAULT   0.04 (shadowband sky factor)
+*
+*            OUTPUTS:    (posdata) daynum, amass, ampress, azim, cosinc,
+*                        elevref, etr, etrn, etrtilt, prime,
+*                        sbcf, sretr, ssetr, unprime, zenref
+*
+*            RETURNS:   Long int status code (defined in solpos.h)
+*
+*    Usage:
+*         In calling program, along with other 'includes', insert:
+*
+*              #include "solpos.h"
+*
+*    Martin Rymes
+*    National Renewable Energy Laboratory
+*    25 March 1998
+*
+*----------------------------------------------------------------------------*/
+
+/*============================================================================
+*
+*     Define the function codes
+*
+*----------------------------------------------------------------------------*/
+#define L_DOY    0x0001
+#define L_GEOM   0x0002
+#define L_ZENETR 0x0004
+#define L_SSHA   0x0008
+#define L_SBCF   0x0010
+#define L_TST    0x0020
+#define L_SRSS   0x0040
+#define L_SOLAZM 0x0080
+#define L_REFRAC 0x0100
+#define L_AMASS  0x0200
+#define L_PRIME  0x0400
+#define L_TILT   0x0800
+#define L_ETR    0x1000
+#define L_ALL    0xFFFF
+
+/*============================================================================
+*
+*     Define the bit-wise masks for each function
+*
+*----------------------------------------------------------------------------*/
+#define S_DOY    ( L_DOY                          )
+#define S_GEOM   ( L_GEOM   | S_DOY               )
+#define S_ZENETR ( L_ZENETR | S_GEOM              )
+#define S_SSHA   ( L_SSHA   | S_GEOM              )
+#define S_SBCF   ( L_SBCF   | S_SSHA              )
+#define S_TST    ( L_TST    | S_GEOM              )
+#define S_SRSS   ( L_SRSS   | S_SSHA   | S_TST    )
+#define S_SOLAZM ( L_SOLAZM | S_ZENETR            )
+#define S_REFRAC ( L_REFRAC | S_ZENETR            )
+#define S_AMASS  ( L_AMASS  | S_REFRAC            )
+#define S_PRIME  ( L_PRIME  | S_AMASS             )
+#define S_TILT   ( L_TILT   | S_SOLAZM | S_REFRAC )
+#define S_ETR    ( L_ETR    | S_REFRAC            )
+#define S_ALL    ( L_ALL                          )
+
+
+/*============================================================================
+*
+*     Enumerate the error codes
+*     (Bit positions are from least significant to most significant)
+*
+*----------------------------------------------------------------------------*/
+/*          Code          Bit       Parameter            Range
+   ===============     ===  ===================  =============   */
+enum
+{ S_YEAR_ERROR,			/*  0   year                  1950 -  2050   */
+    S_MONTH_ERROR,		/*  1   month                    1 -    12   */
+    S_DAY_ERROR,		/*  2   day-of-month             1 -    31   */
+    S_DOY_ERROR,		/*  3   day-of-year              1 -   366   */
+    S_HOUR_ERROR,		/*  4   hour                     0 -    24   */
+    S_MINUTE_ERROR,		/*  5   minute                   0 -    59   */
+    S_SECOND_ERROR,		/*  6   second                   0 -    59   */
+    S_TZONE_ERROR,		/*  7   time zone              -12 -    12   */
+    S_INTRVL_ERROR,		/*  8   interval (seconds)       0 - 28800   */
+    S_LAT_ERROR,		/*  9   latitude               -90 -    90   */
+    S_LON_ERROR,		/* 10   longitude             -180 -   180   */
+    S_TEMP_ERROR,		/* 11   temperature (deg. C)  -100 -   100   */
+    S_PRESS_ERROR,		/* 12   pressure (millibars)     0 -  2000   */
+    S_TILT_ERROR,		/* 13   tilt                   -90 -    90   */
+    S_ASPECT_ERROR,		/* 14   aspect                -360 -   360   */
+    S_SBWID_ERROR,		/* 15   shadow band width (cm)   1 -   100   */
+    S_SBRAD_ERROR,		/* 16   shadow band radius (cm)  1 -   100   */
+    S_SBSKY_ERROR
+};				/* 17   shadow band sky factor  -1 -     1   */
+
+struct posdata
+{
+
+  /***** ALPHABETICAL LIST OF COMMON VARIABLES *****/
+    /* Each comment begins with a 1-column letter code:
+       I:  INPUT variable
+       O:  OUTPUT variable
+       T:  TRANSITIONAL variable used in the algorithm,
+       of interest only to the solar radiation
+       modelers, and available to you because you
+       may be one of them.
+
+       The FUNCTION column indicates which sub-function
+       within solpos must be switched on using the
+       "function" parameter to calculate the desired
+       output variable.  All function codes are
+       defined in the solpos.h file.  The default
+       S_ALL switch calculates all output variables.
+       Multiple functions may be or'd to create a
+       composite function switch.  For example,
+       (S_TST | S_SBCF). Specifying only the functions
+       for required output variables may allow solpos
+       to execute more quickly.
+
+       The S_DOY mask works as a toggle between the
+       input date represented as a day number (daynum)
+       or as month and day.  To set the switch (to
+       use daynum input), the function is or'd; to
+       clear the switch (to use month and day input),
+       the function is inverted and and'd.
+
+       For example:
+       pdat->function |= S_DOY (sets daynum input)
+       pdat->function &= ~S_DOY (sets month and day input)
+
+       Whichever date form is used, S_solpos will
+       calculate and return the variables(s) of the
+       other form.  See the soltest.c program for
+       other examples. */
+
+    /* VARIABLE        I/O  Function    Description */
+    /* -------------  ----  ----------  --------------------------------------- */
+
+    int day;			/* I/O: S_DOY      Day of month (May 27 = 27, etc.)
+				   solpos will CALCULATE this by default,
+				   or will optionally require it as input
+				   depending on the setting of the S_DOY
+				   function switch. */
+    int daynum;			/* I/O: S_DOY      Day number (day of year; Feb 1 = 32 )
+				   solpos REQUIRES this by default, but
+				   will optionally calculate it from
+				   month and day depending on the setting
+				   of the S_DOY function switch. */
+    int function;		/* I:              Switch to choose functions for desired
+				   output. */
+    int hour;			/* I:              Hour of day, 0 - 23, DEFAULT = 12 */
+    int interval;		/* I:              Interval of a measurement period in
+				   seconds.  Forces solpos to use the
+				   time and date from the interval
+				   midpoint. The INPUT time (hour,
+				   minute, and second) is assumed to
+				   be the END of the measurement
+				   interval. */
+    int minute;			/* I:              Minute of hour, 0 - 59, DEFAULT = 0 */
+    int month;			/* I/O: S_DOY      Month number (Jan = 1, Feb = 2, etc.)
+				   solpos will CALCULATE this by default,
+				   or will optionally require it as input
+				   depending on the setting of the S_DOY
+				   function switch. */
+    int second;			/* I:              Second of minute, 0 - 59, DEFAULT = 0 */
+    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 */
+    float ampress;		/* O:  S_AMASS    Pressure-corrected airmass */
+    float aspect;		/* I:             Azimuth of panel surface (direction it
+				   faces) N=0, E=90, S=180, W=270,
+				   DEFAULT = 180 */
+    float azim;			/* O:  S_SOLAZM   Solar azimuth angle:  N=0, E=90, S=180,
+				   W=270 */
+    float cosinc;		/* O:  S_TILT     Cosine of solar incidence angle on
+				   panel */
+    float coszen;		/* O:  S_REFRAC   Cosine of refraction corrected solar
+				   zenith angle */
+    float dayang;		/* T:  S_GEOM     Day angle (daynum*360/year-length)
+				   degrees */
+    float declin;		/* T:  S_GEOM     Declination--zenith angle of solar noon
+				   at equator, degrees NORTH */
+    float eclong;		/* T:  S_GEOM     Ecliptic longitude, degrees */
+    float ecobli;		/* T:  S_GEOM     Obliquity of ecliptic */
+    float ectime;		/* T:  S_GEOM     Time of ecliptic calculations */
+    float elevetr;		/* O:  S_ZENETR   Solar elevation, no atmospheric
+				   correction (= ETR) */
+    float elevref;		/* O:  S_REFRAC   Solar elevation angle,
+				   deg. from horizon, refracted */
+    float eqntim;		/* T:  S_TST      Equation of time (TST - LMT), minutes */
+    float erv;			/* T:  S_GEOM     Earth radius vector
+				   (multiplied to solar constant) */
+    float etr;			/* O:  S_ETR      Extraterrestrial (top-of-atmosphere)
+				   W/sq m global horizontal solar
+				   irradiance */
+    float etrn;			/* O:  S_ETR      Extraterrestrial (top-of-atmosphere)
+				   W/sq m direct normal solar
+				   irradiance */
+    float etrtilt;		/* O:  S_TILT     Extraterrestrial (top-of-atmosphere)
+				   W/sq m global irradiance on a tilted
+				   surface */
+    float gmst;			/* T:  S_GEOM     Greenwich mean sidereal time, hours */
+    float hrang;		/* T:  S_GEOM     Hour angle--hour of sun from solar noon,
+				   degrees WEST */
+    float julday;		/* T:  S_GEOM     Julian Day of 1 JAN 2000 minus
+				   2,400,000 days (in order to regain
+				   single precision) */
+    float latitude;		/* I:             Latitude, degrees north (south negative) */
+    float longitude;		/* I:             Longitude, degrees east (west negative) */
+    float lmst;			/* T:  S_GEOM     Local mean sidereal time, degrees */
+    float mnanom;		/* T:  S_GEOM     Mean anomaly, degrees */
+    float mnlong;		/* T:  S_GEOM     Mean longitude, degrees */
+    float rascen;		/* T:  S_GEOM     Right ascension, degrees */
+    float press;		/* I:             Surface pressure, millibars, used for
+				   refraction correction and ampress */
+    float prime;		/* O:  S_PRIME    Factor that normalizes Kt, Kn, etc. */
+    float sbcf;			/* O:  S_SBCF     Shadow-band correction factor */
+    float sbwid;		/* I:             Shadow-band width (cm) */
+    float sbrad;		/* I:             Shadow-band radius (cm) */
+    float sbsky;		/* I:             Shadow-band sky factor */
+    float solcon;		/* I:             Solar constant (NREL uses 1367 W/sq m) */
+    float ssha;			/* T:  S_SRHA     Sunset(/rise) hour angle, degrees */
+    float sretr;		/* O:  S_SRSS     Sunrise time, minutes from midnight,
+				   local, WITHOUT refraction */
+    float ssetr;		/* O:  S_SRSS     Sunset time, minutes from midnight,
+				   local, WITHOUT refraction */
+    float temp;			/* I:             Ambient dry-bulb temperature, degrees C,
+				   used for refraction correction */
+    float tilt;			/* I:             Degrees tilt from horizontal of panel */
+    float timezone;		/* I:             Time zone, east (west negative).
+				   USA:  Mountain = -7, Central = -6, etc. */
+    float tst;			/* T:  S_TST      True solar time, minutes from midnight */
+    float tstfix;		/* T:  S_TST      True solar time - local standard time */
+    float unprime;		/* O:  S_PRIME    Factor that denormalizes Kt', Kn', etc. */
+    float utime;		/* T:  S_GEOM     Universal (Greenwich) standard time */
+    float zenetr;		/* T:  S_ZENETR   Solar zenith angle, no atmospheric
+				   correction (= ETR) */
+    float zenref;		/* O:  S_REFRAC   Solar zenith angle, deg. from zenith,
+				   refracted */
+};
+
+/* For users that wish to access individual functions, the following table
+   lists all output and transition variables, the L_ mask for the function
+   that calculates them, and all the input variables required by that function.
+   The function variable is set to the L_ mask, which will force S_solpos to
+   only call the required function.  L_ masks may be ORed as desired.
+
+   VARIABLE      Mask       Required Variables
+   ---------  ----------  ---------------------------------------
+   amass      L_AMASS    zenref, press
+   ampress    L_AMASS    zenref, press
+   azim       L_SOLAZM   elevetr, declin, latitude, hrang
+   cosinc     L_TILT     azim, aspect, tilt, zenref, coszen,etrn
+   coszen     L_REFRAC   elevetr, press, temp
+   dayang     L_GEOM     All date, time, and location inputs
+   declin     L_GEOM     All date, time, and location inputs
+   eclong     L_GEOM     All date, time, and location inputs
+   ecobli     L_GEOM     All date, time, and location inputs
+   ectime     L_GEOM     All date, time, and location inputs
+   elevetr    L_ZENETR   declin, latitude, hrang
+   elevref    L_REFRAC   elevetr, press, temp
+   eqntim     L_TST      hrang, hour, minute, second, interval
+   erv        L_GEOM     All date, time, and location inputs
+   etr        L_ETR      coszen, solcon, erv
+   etrn       L_ETR      coszen, solcon, erv
+   etrtilt    L_TILT     azim, aspect, tilt, zenref, coszen, etrn
+   gmst       L_GEOM     All date, time, and location inputs
+   hrang      L_GEOM     All date, time, and location inputs
+   julday     L_GEOM     All date, time, and location inputs
+   lmst       L_GEOM     All date, time, and location inputs
+   mnanom     L_GEOM     All date, time, and location inputs
+   mnlong     L_GEOM     All date, time, and location inputs
+   rascen     L_GEOM     All date, time, and location inputs
+   prime      L_PRIME    amass
+   sbcf       L_SBCF     latitude, declin, ssha, sbwid, sbrad, sbsky
+   ssha       L_SRHA     latitude, declin
+   sretr      L_SRSS     ssha, tstfix
+   ssetr      L_SRSS     ssha, tstfix
+   tst        L_TST      hrang, hour, minute, second, interval
+   tstfix     L_TST      hrang, hour, minute, second, interval
+   unprime    L_PRIME    amass
+   utime      L_GEOM     All date, time, and location inputs
+   zenetr     L_ZENETR   declination, latitude, hrang
+   zenref     L_REFRAC   elevetr, press, temp
+
+ */
+
+/*============================================================================
+*    Long int function S_solpos, adapted from the NREL VAX solar libraries
+*
+*    This function calculates the apparent solar position and intensity
+*    (theoretical maximum solar energy) based on the date, time, and
+*    location on Earth. (DEFAULT values are from the optional S_posinit
+*    function.)
+*
+*    Requires:
+*        Date and time:
+*            year
+*            month  (optional without daynum)
+*            day    (optional without daynum)
+*            daynum
+*            hour
+*            minute
+*            second
+*        Location:
+*            latitude
+*            longitude
+*        Location/time adjuster:
+*            timezone
+*        Atmospheric pressure and temperature:
+*            press     DEFAULT 1013.0 mb
+*            temp      DEFAULT 10.0 degrees C
+*        Tilt of flat surface that receives solar energy:
+*            aspect    DEFAULT 180 (South)
+*            tilt      DEFAULT 0 (Horizontal)
+*        Shadow band parameters:
+*            sbwid     DEFAULT 7.6 cm
+*            sbrad     DEFAULT 31.7 cm
+*            sbsky     DEFAULT 0.04
+*        Functionality
+*            function  DEFAULT S_ALL (all output parameters computed)
+*
+*    Returns:
+*        everything defined at the top of this listing.
+*----------------------------------------------------------------------------*/
+long S_solpos(struct posdata *pdat);
+
+/*============================================================================
+*    Void function S_init
+*
+*    This function initiates all of the input functions to S_Solpos().
+*    NOTE: This function is optional if you initialize all input parameters
+*          in your calling code.
+*
+*    Requires: Pointer to a posdata structure, members of which are
+*           initialized.
+*
+*    Returns: Void
+*
+*----------------------------------------------------------------------------*/
+void S_init(struct posdata *pdat);
+
+
+/*============================================================================
+*    Void function S_decode
+*
+*    This function decodes the error codes from S_solpos return value
+*
+*    INPUTS: Long integer S_solpos return value, struct posdata*
+*
+*    OUTPUTS: Descriptive text of errors to stderr
+*----------------------------------------------------------------------------*/
+void S_decode(long code, struct posdata *pdat);