Geografic-Information-System/raster/r.sunmask/solpos00.c

/*============================================================================
*    Contains:
*        S_solpos     (computes solar position and intensity
*                      from time and place)
*
*            INPUTS:     (via posdata struct) year, daynum, hour,
*                        minute, second, latitude, longitude, timezone,
*                        intervl
*            OPTIONAL:   (via posdata struct) month, day, press, temp, tilt,
*                        aspect, function
*            OUTPUTS:    EVERY variable in the struct posdata
*                            (defined in solpos.h)
*
*                       NOTE: Certain conditions exist during which some of
*                       the output variables are undefined or cannot be
*                       calculated.  In these cases, the variables are
*                       returned with flag values indicating such.  In other
*                       cases, the variables may return a realistic, though
*                       invalid, value. These variables and the flag values
*                       or invalid conditions are listed below:
*
*                       amass     -1.0 at zenetr angles greater than 93.0
*                                 degrees
*                       ampress   -1.0 at zenetr angles greater than 93.0
*                                 degrees
*                       azim      invalid at zenetr angle 0.0 or latitude
*                                 +/-90.0 or at night
*                       elevetr   limited to -9 degrees at night
*                       etr       0.0 at night
*                       etrn      0.0 at night
*                       etrtilt   0.0 when cosinc is less than 0
*                       prime     invalid at zenetr angles greater than 93.0
*                                 degrees
*                       sretr     +/- 2999.0 during periods of 24 hour sunup or
*                                 sundown
*                       ssetr     +/- 2999.0 during periods of 24 hour sunup or
*                                 sundown
*                       ssha      invalid at the North and South Poles
*                       unprime   invalid at zenetr angles greater than 93.0
*                                 degrees
*                       zenetr    limited to 99.0 degrees at night
*
*        S_init       (optional initialization for all input parameters in
*                      the posdata struct)
*           INPUTS:     struct posdata*
*           OUTPUTS:    struct posdata*
*
*                     (Note: initializes the required S_solpos INPUTS above
*                      to out-of-bounds conditions, forcing the user to
*                      supply the parameters; initializes the OPTIONAL
*                      S_solpos inputs above to nominal values.)
*
*       S_decode      (optional utility for decoding the S_solpos return code)
*           INPUTS:     long integer S_solpos return value, struct posdata*
*           OUTPUTS:    text to stderr
*
*    Usage:
*         In calling program, just after other 'includes', insert:
*
*              #include "solpos.h"
*
*         Function calls:
*              S_init(struct posdata*)  [optional]
*              .
*              .
*              [set time and location parameters before S_solpos call]
*              .
*              .
*              int retval = S_solpos(struct posdata*)
*              S_decode(int retval, struct posdata*) [optional]
*                  (Note: you should always look at the S_solpos return
*                   value, which contains error codes. S_decode is one option
*                   for examining these codes.  It can also serve as a
*                   template for building your own application-specific
*                   decoder.)
*
*    Martin Rymes
*    National Renewable Energy Laboratory
*    25 March 1998
*
*    27 April 1999 REVISION:  Corrected leap year in S_date.
*    13 January 2000 REVISION:  SMW converted to structure posdata parameter
*                               and subdivided into functions.
*    01 February 2001 REVISION: SMW corrected ecobli calculation
*                               (changed sign). Error is small (max 0.015 deg
*                               in calculation of declination angle)
*    11 April 2001 REVISION:    SMW corrected parenthesis grouping in
*                               validation() function for pressure and
*                               shadowband.  As it was, the validation
*                               routine would have checked for and reported
*                               high limit errors when the functions had not
*                               been selected.
*----------------------------------------------------------------------------*/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include "solpos00.h"

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 *
 * Structures defined for this module
 *
 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
struct trigdata			/* used to pass calculated values locally */
{
    float cd;			/* cosine of the declination */
    float ch;			/* cosine of the hour angle */
    float cl;			/* cosine of the latitude */
    float sd;			/* sine of the declination */
    float sl;			/* sine of the latitude */
};


/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 *
 * Temporary global variables used only in this file:
 *
 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
static int month_days[2][13] = { {0, 0, 31, 59, 90, 120, 151,
				  181, 212, 243, 273, 304, 334},
{0, 0, 31, 60, 91, 121, 152,
 182, 213, 244, 274, 305, 335}
};

		   /* cumulative number of days prior to beginning of month */

/* TODO: the "degrad" - "raddeg" variable names are flipped - typo or severe bug?? */
static float degrad = 57.295779513;	/* converts from radians to degrees */
static float raddeg = 0.0174532925;	/* converts from degrees to radians */

/*============================================================================
*    Local function prototypes
============================================================================*/
static long int validate(struct posdata *pdat);
static void dom2doy(struct posdata *pdat);
static void doy2dom(struct posdata *pdat);
static void geometry(struct posdata *pdat);
static void zen_no_ref(struct posdata *pdat, struct trigdata *tdat);
static void ssha(struct posdata *pdat, struct trigdata *tdat);
static void sbcf(struct posdata *pdat, struct trigdata *tdat);
static void tst(struct posdata *pdat);
static void srss(struct posdata *pdat);
static void sazm(struct posdata *pdat, struct trigdata *tdat);
static void refrac(struct posdata *pdat);
static void amass(struct posdata *pdat);
static void prime(struct posdata *pdat);
static void etr(struct posdata *pdat);
static void tilt(struct posdata *pdat);
static void localtrig(struct posdata *pdat, struct trigdata *tdat);

/*============================================================================
*    Long integer function S_solpos, adapted from the VAX solar libraries
*
*    This function calculates the apparent solar position and the
*    intensity of the sun (theoretical maximum solar energy) from
*    time and place on Earth.
*
*    Requires (from the struct posdata parameter):
*        Date and time:
*            year
*            daynum   (requirement depends on the S_DOY switch)
*            month    (requirement depends on the S_DOY switch)
*            day      (requirement depends on the S_DOY switch)
*            hour
*            minute
*            second
*            interval  DEFAULT 0
*        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)
*        Function Switch (codes defined in solpos.h)
*            function  DEFAULT S_ALL
*
*    Returns (via the struct posdata parameter):
*        everything defined in the struct posdata in solpos.h.
*----------------------------------------------------------------------------*/
long S_solpos(struct posdata *pdat)
{
    long int retval;

    struct trigdata trigdat, *tdat;

    tdat = &trigdat;		/* point to the structure */

    /* initialize the trig structure */
    tdat->sd = -999.0;		/* flag to force calculation of trig data */
    tdat->cd = 1.0;
    tdat->ch = 1.0;		/* set the rest of these to something safe */
    tdat->cl = 1.0;
    tdat->sl = 1.0;

    if ((retval = validate(pdat)) != 0)	/* validate the inputs */
	return retval;


    if (pdat->function & L_DOY)
	doy2dom(pdat);		/* convert input doy to month-day */
    else
	dom2doy(pdat);		/* convert input month-day to doy */

    if (pdat->function & L_GEOM)
	geometry(pdat);		/* do basic geometry calculations */

    if (pdat->function & L_ZENETR)	/* etr at non-refracted zenith angle */
	zen_no_ref(pdat, tdat);

    if (pdat->function & L_SSHA)	/* Sunset hour calculation */
	ssha(pdat, tdat);

    if (pdat->function & L_SBCF)	/* Shadowband correction factor */
	sbcf(pdat, tdat);

    if (pdat->function & L_TST)	/* true solar time */
	tst(pdat);

    if (pdat->function & L_SRSS)	/* sunrise/sunset calculations */
	srss(pdat);

    if (pdat->function & L_SOLAZM)	/* solar azimuth calculations */
	sazm(pdat, tdat);

    if (pdat->function & L_REFRAC)	/* atmospheric refraction calculations */
	refrac(pdat);

    if (pdat->function & L_AMASS)	/* airmass calculations */
	amass(pdat);

    if (pdat->function & L_PRIME)	/* kt-prime/unprime calculations */
	prime(pdat);

    if (pdat->function & L_ETR)	/* ETR and ETRN (refracted) */
	etr(pdat);

    if (pdat->function & L_TILT)	/* tilt calculations */
	tilt(pdat);

    return 0;
}


/*============================================================================
*    Void function S_init
*
*    This function initiates all of the input parameters in the struct
*    posdata passed to S_solpos().  Initialization is either to nominal
*    values or to out of range values, which forces the calling program to
*    specify parameters.
*
*    NOTE: This function is optional if you initialize ALL input parameters
*          in your calling code.  Note that the required parameters of date
*          and location are deliberately initialized out of bounds to force
*          the user to enter real-world values.
*
*    Requires: Pointer to a posdata structure, members of which are
*           initialized.
*
*    Returns: Void
*----------------------------------------------------------------------------*/
void S_init(struct posdata *pdat)
{
    pdat->day = -99;		/* Day of month (May 27 = 27, etc.) */
    pdat->daynum = -999;	/* Day number (day of year; Feb 1 = 32 ) */
    pdat->hour = -99;		/* Hour of day, 0 - 23 */
    pdat->minute = -99;		/* Minute of hour, 0 - 59 */
    pdat->month = -99;		/* Month number (Jan = 1, Feb = 2, etc.) */
    pdat->second = -99;		/* Second of minute, 0 - 59 */
    pdat->year = -99;		/* 4-digit year */
    pdat->interval = 0;		/* instantaneous measurement interval */
    pdat->aspect = 180.0;	/* Azimuth of panel surface (direction it
				   faces) N=0, E=90, S=180, W=270 */
    pdat->latitude = -99.0;	/* Latitude, degrees north (south negative) */
    pdat->longitude = -999.0;	/* Longitude, degrees east (west negative) */
    pdat->press = 1013.0;	/* Surface pressure, millibars */
    pdat->solcon = 1367.0;	/* Solar constant, 1367 W/sq m */
    pdat->temp = 15.0;		/* Ambient dry-bulb temperature, degrees C */
    pdat->tilt = 0.0;		/* Degrees tilt from horizontal of panel */
    pdat->timezone = -99.0;	/* Time zone, east (west negative). */
    pdat->sbwid = 7.6;		/* Eppley shadow band width */
    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 */
}


/*============================================================================
*    Local long int function validate
*
*    Validates the input parameters
*----------------------------------------------------------------------------*/
static long int validate(struct posdata *pdat)
{

    long int retval = 0;	/* start with no errors */

    /* No absurd dates, please. */
    if (pdat->function & L_GEOM) {
	if ((pdat->year < 1950) || (pdat->year > 2050))	/* limits of algoritm */
	    retval |= (1L << S_YEAR_ERROR);
	if (!(pdat->function & S_DOY) &&
	    ((pdat->month < 1) || (pdat->month > 12)))
	    retval |= (1L << S_MONTH_ERROR);
	if (!(pdat->function & S_DOY) &&
	    ((pdat->day < 1) || (pdat->day > 31)))
	    retval |= (1L << S_DAY_ERROR);
	if ((pdat->function & S_DOY) &&
	    ((pdat->daynum < 1) || (pdat->daynum > 366)))
	    retval |= (1L << S_DOY_ERROR);

	/* No absurd times, please. */
	if ((pdat->hour < 0) || (pdat->hour > 24))
	    retval |= (1L << S_HOUR_ERROR);
	if ((pdat->minute < 0) || (pdat->minute > 59))
	    retval |= (1L << S_MINUTE_ERROR);
	if ((pdat->second < 0) || (pdat->second > 59))
	    retval |= (1L << S_SECOND_ERROR);
	if ((pdat->hour == 24) && (pdat->minute > 0))	/* no more than 24 hrs */
	    retval |= ((1L << S_HOUR_ERROR) | (1L << S_MINUTE_ERROR));
	if ((pdat->hour == 24) && (pdat->second > 0))	/* no more than 24 hrs */
	    retval |= ((1L << S_HOUR_ERROR) | (1L << S_SECOND_ERROR));
	if (fabs(pdat->timezone) > 12.0)
	    retval |= (1L << S_TZONE_ERROR);
	if ((pdat->interval < 0) || (pdat->interval > 28800))
	    retval |= (1L << S_INTRVL_ERROR);

	/* No absurd locations, please. */
	if (fabs(pdat->longitude) > 180.0)
	    retval |= (1L << S_LON_ERROR);
	if (fabs(pdat->latitude) > 90.0)
	    retval |= (1L << S_LAT_ERROR);
    }

    /* No silly temperatures or pressures, please. */
    if ((pdat->function & L_REFRAC) && (fabs(pdat->temp) > 100.0))
	retval |= (1L << S_TEMP_ERROR);
    if ((pdat->function & L_REFRAC) &&
	((pdat->press < 0.0) || (pdat->press > 2000.0)))
	retval |= (1L << S_PRESS_ERROR);

    /* No out of bounds tilts, please */
    if ((pdat->function & L_TILT) && (fabs(pdat->tilt) > 180.0))
	retval |= (1L << S_TILT_ERROR);
    if ((pdat->function & L_TILT) && (fabs(pdat->aspect) > 360.0))
	retval |= (1L << S_ASPECT_ERROR);

    /* No oddball shadowbands, please */
    if ((pdat->function & L_SBCF) &&
	((pdat->sbwid < 1.0) || (pdat->sbwid > 100.0)))
	retval |= (1L << S_SBWID_ERROR);
    if ((pdat->function & L_SBCF) &&
	((pdat->sbrad < 1.0) || (pdat->sbrad > 100.0)))
	retval |= (1L << S_SBRAD_ERROR);
    if ((pdat->function & L_SBCF) && (fabs(pdat->sbsky) > 1.0))
	retval |= (1L << S_SBSKY_ERROR);

    return retval;
}


/*============================================================================
*    Local Void function dom2doy
*
*    Converts day-of-month to day-of-year
*
*    Requires (from struct posdata parameter):
*            year
*            month
*            day
*
*    Returns (via the struct posdata parameter):
*            year
*            daynum
*----------------------------------------------------------------------------*/
static void dom2doy(struct posdata *pdat)
{
    pdat->daynum = pdat->day + month_days[0][pdat->month];

    /* (adjust for leap year) */
    if (((pdat->year % 4) == 0) &&
	(((pdat->year % 100) != 0) || ((pdat->year % 400) == 0)) &&
	(pdat->month > 2))
	pdat->daynum += 1;
}


/*============================================================================
*    Local void function doy2dom
*
*    This function computes the month/day from the day number.
*
*    Requires (from struct posdata parameter):
*        Year and day number:
*            year
*            daynum
*
*    Returns (via the struct posdata parameter):
*            year
*            month
*            day
*----------------------------------------------------------------------------*/
static void doy2dom(struct posdata *pdat)
{
    int imon;			/* Month (month_days) array counter */
    int leap;			/* leap year switch */

    /* Set the leap year switch */
    if (((pdat->year % 4) == 0) &&
	(((pdat->year % 100) != 0) || ((pdat->year % 400) == 0)))
	leap = 1;
    else
	leap = 0;

    /* Find the month */
    imon = 12;
    while (pdat->daynum <= month_days[leap][imon])
	--imon;

    /* Set the month and day of month */
    pdat->month = imon;
    pdat->day = pdat->daynum - month_days[leap][imon];
}


/*============================================================================
*    Local Void function geometry
*
*    Does the underlying geometry for a given time and location
*----------------------------------------------------------------------------*/
static void geometry(struct posdata *pdat)
{
    float bottom;		/* denominator (bottom) of the fraction */
    float c2;			/* cosine of d2 */
    float cd;			/* cosine of the day angle or delination */
    float d2;			/* pdat->dayang times two */
    float delta;		/* difference between current year and 1949 */
    float s2;			/* sine of d2 */
    float sd;			/* sine of the day angle */
    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;

    /* 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);

    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;

    /* 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;

    /* 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;

    /* (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;

    /* (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);

    /* (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. */

    /* 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;

    /* 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));

    /* 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);

    pdat->rascen = degrad * atan2(top, bottom);

    /* (make it a positive angle) */
    if (pdat->rascen < 0.0)
	pdat->rascen += 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;

    /* (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;

    /* 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;
}


/*============================================================================
*    Local Void function zen_no_ref
*
*    ETR solar zenith angle
*       Iqbal, M.  1983.  An Introduction to Solar Radiation.
*            Academic Press, NY., page 15
*----------------------------------------------------------------------------*/
static void zen_no_ref(struct posdata *pdat, struct trigdata *tdat)
{
    float cz;			/* cosine of the solar zenith angle */

    localtrig(pdat, tdat);
    cz = tdat->sd * tdat->sl + tdat->cd * tdat->cl * tdat->ch;

    /* (watch out for the roundoff errors) */
    if (fabs(cz) > 1.0) {
	if (cz >= 0.0)
	    cz = 1.0;
	else
	    cz = -1.0;
    }

    pdat->zenetr = acos(cz) * degrad;

    /* (limit the degrees below the horizon to 9 [+90 -> 99]) */
    if (pdat->zenetr > 99.0)
	pdat->zenetr = 99.0;

    pdat->elevetr = 90.0 - pdat->zenetr;
}


/*============================================================================
*    Local Void function ssha
*
*    Sunset hour angle, degrees
*       Iqbal, M.  1983.  An Introduction to Solar Radiation.
*            Academic Press, NY., page 16
*----------------------------------------------------------------------------*/
static void ssha(struct posdata *pdat, struct trigdata *tdat)
{
    float cssha;		/* cosine of the sunset hour angle */
    float cdcl;			/* ( cd * cl ) */

    localtrig(pdat, tdat);
    cdcl = tdat->cd * tdat->cl;

    if (fabs(cdcl) >= 0.001) {
	cssha = -tdat->sl * tdat->sd / cdcl;

	/* This keeps the cosine from blowing on roundoff */
	if (cssha < -1.0)
	    pdat->ssha = 180.0;
	else if (cssha > 1.0)
	    pdat->ssha = 0.0;
	else
	    pdat->ssha = degrad * acos(cssha);
    }
    else if (((pdat->declin >= 0.0) && (pdat->latitude > 0.0)) ||
	     ((pdat->declin < 0.0) && (pdat->latitude < 0.0)))
	pdat->ssha = 180.0;
    else
	pdat->ssha = 0.0;
}


/*============================================================================
*    Local Void function sbcf
*
*    Shadowband correction factor
*       Drummond, A. J.  1956.  A contribution to absolute pyrheliometry.
*            Q. J. R. Meteorol. Soc. 82, pp. 481-493
*----------------------------------------------------------------------------*/
static void sbcf(struct posdata *pdat, struct trigdata *tdat)
{
    float p, t1, t2;		/* used to compute sbcf */

    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);
    pdat->sbcf = pdat->sbsky + 1.0 / (1.0 - p * (t1 + t2));

}


/*============================================================================
*    Local Void function tst
*
*    TST -> True Solar Time = local standard time + TSTfix, time
*      in minutes from midnight.
*        Iqbal, M.  1983.  An Introduction to Solar Radiation.
*            Academic Press, NY., page 13
*----------------------------------------------------------------------------*/
static void tst(struct posdata *pdat)
{
    pdat->tst = (180.0 + pdat->hrang) * 4.0;
    pdat->tstfix = pdat->tst - (float)pdat->hour * 60.0 - pdat->minute - (float)pdat->second / 60.0 + (float)pdat->interval / 120.0;	/* add back half of the interval */

    /* bound tstfix to this day */
    while (pdat->tstfix > 720.0)
	pdat->tstfix -= 1440.0;
    while (pdat->tstfix < -720.0)
	pdat->tstfix += 1440.0;

    pdat->eqntim =
	pdat->tstfix + 60.0 * pdat->timezone - 4.0 * pdat->longitude;

}


/*============================================================================
*    Local Void function srss
*
*    Sunrise and sunset times (minutes from midnight)
*----------------------------------------------------------------------------*/
static void srss(struct posdata *pdat)
{
    if (pdat->ssha <= 1.0) {
	pdat->sretr = 2999.0;
	pdat->ssetr = -2999.0;
    }
    else if (pdat->ssha >= 179.0) {
	pdat->sretr = -2999.0;
	pdat->ssetr = 2999.0;
    }
    else {
	pdat->sretr = 720.0 - 4.0 * pdat->ssha - pdat->tstfix;
	pdat->ssetr = 720.0 + 4.0 * pdat->ssha - pdat->tstfix;
    }
}


/*============================================================================
*    Local Void function sazm
*
*    Solar azimuth angle
*       Iqbal, M.  1983.  An Introduction to Solar Radiation.
*            Academic Press, NY., page 15
*----------------------------------------------------------------------------*/
static void sazm(struct posdata *pdat, struct trigdata *tdat)
{
    float ca;			/* cosine of the solar azimuth angle */
    float ce;			/* cosine of the solar elevation */
    float cecl;			/* ( ce * cl ) */
    float se;			/* sine of the solar elevation */

    localtrig(pdat, tdat);
    ce = cos(raddeg * pdat->elevetr);
    se = sin(raddeg * pdat->elevetr);

    pdat->azim = 180.0;
    cecl = ce * tdat->cl;
    if (fabs(cecl) >= 0.001) {
	ca = (se * tdat->sl - tdat->sd) / cecl;
	if (ca > 1.0)
	    ca = 1.0;
	else if (ca < -1.0)
	    ca = -1.0;

	pdat->azim = 180.0 - acos(ca) * degrad;
	if (pdat->hrang > 0)
	    pdat->azim = 360.0 - pdat->azim;
    }
}


/*============================================================================
*    Local Int function refrac
*
*    Refraction correction, degrees
*        Zimmerman, John C.  1981.  Sun-pointing programs and their
*            accuracy.
*            SAND81-0761, Experimental Systems Operation Division 4721,
*            Sandia National Laboratories, Albuquerque, NM.
*----------------------------------------------------------------------------*/
static void refrac(struct posdata *pdat)
{
    float prestemp;		/* temporary pressure/temperature correction */
    float refcor;		/* temporary refraction correction */
    float tanelev;		/* tangent of the solar elevation angle */

    /* If the sun is near zenith, the algorithm bombs; refraction near 0 */
    if (pdat->elevetr > 85.0)
	refcor = 0.0;

    /* Otherwise, we have refraction */
    else {
	tanelev = tan(raddeg * pdat->elevetr);
	if (pdat->elevetr >= 5.0)
	    refcor = 58.1 / tanelev -
		0.07 / (pow(tanelev, 3)) + 0.000086 / (pow(tanelev, 5));
	else if (pdat->elevetr >= -0.575)
	    refcor = 1735.0 +
		pdat->elevetr * (-518.2 + pdat->elevetr * (103.4 +
							   pdat->elevetr *
							   (-12.79 +
							    pdat->elevetr *
							    0.711)));
	else
	    refcor = -20.774 / tanelev;

	prestemp = (pdat->press * 283.0) / (1013.0 * (273.0 + pdat->temp));
	refcor *= prestemp / 3600.0;
    }

    /* Refracted solar elevation angle */
    pdat->elevref = pdat->elevetr + refcor;

    /* (limit the degrees below the horizon to 9) */
    if (pdat->elevref < -9.0)
	pdat->elevref = -9.0;

    /* Refracted solar zenith angle */
    pdat->zenref = 90.0 - pdat->elevref;
    pdat->coszen = cos(raddeg * pdat->zenref);
}


/*============================================================================
*    Local Void function  amass
*
*    Airmass
*       Kasten, F. and Young, A.  1989.  Revised optical air mass
*            tables and approximation formula.  Applied Optics 28 (22),
*            pp. 4735-4738
*----------------------------------------------------------------------------*/
static void amass(struct posdata *pdat)
{
    if (pdat->zenref > 93.0) {
	pdat->amass = -1.0;
	pdat->ampress = -1.0;
    }
    else {
	pdat->amass =
	    1.0 / (cos(raddeg * pdat->zenref) + 0.50572 *
		   pow((96.07995 - pdat->zenref), -1.6364));

	pdat->ampress = pdat->amass * pdat->press / 1013.0;
    }
}


/*============================================================================
*    Local Void function prime
*
*    Prime and Unprime
*    Prime  converts Kt to normalized Kt', etc.
*       Unprime deconverts Kt' to Kt, etc.
*            Perez, R., P. Ineichen, Seals, R., & Zelenka, A.  1990.  Making
*            full use of the clearness index for parameterizing hourly
*            insolation conditions. Solar Energy 45 (2), pp. 111-114
*----------------------------------------------------------------------------*/
static void prime(struct posdata *pdat)
{
    pdat->unprime = 1.031 * exp(-1.4 / (0.9 + 9.4 / pdat->amass)) + 0.1;
    pdat->prime = 1.0 / pdat->unprime;
}


/*============================================================================
*    Local Void function etr
*
*    Extraterrestrial (top-of-atmosphere) solar irradiance
*----------------------------------------------------------------------------*/
static void etr(struct posdata *pdat)
{
    if (pdat->coszen > 0.0) {
	pdat->etrn = pdat->solcon * pdat->erv;
	pdat->etr = pdat->etrn * pdat->coszen;
    }
    else {
	pdat->etrn = 0.0;
	pdat->etr = 0.0;
    }
}


/*============================================================================
*    Local Void function localtrig
*
*    Does trig on internal variable used by several functions
*----------------------------------------------------------------------------*/
static void localtrig(struct posdata *pdat, struct trigdata *tdat)
{
    /* define masks to prevent calculation of uninitialized variables */
#define SD_MASK ( L_ZENETR | L_SSHA | S_SBCF | S_SOLAZM )
#define SL_MASK ( L_ZENETR | L_SSHA | S_SBCF | S_SOLAZM )
#define CL_MASK ( L_ZENETR | L_SSHA | S_SBCF | S_SOLAZM )
#define CD_MASK ( L_ZENETR | L_SSHA | S_SBCF )
#define CH_MASK ( L_ZENETR )

    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);
	if (pdat->function | CH_MASK)
	    tdat->ch = cos(raddeg * pdat->hrang);
	if (pdat->function | CL_MASK)
	    tdat->cl = cos(raddeg * pdat->latitude);
	if (pdat->function | SD_MASK)
	    tdat->sd = sin(raddeg * pdat->declin);
	if (pdat->function | SL_MASK)
	    tdat->sl = sin(raddeg * pdat->latitude);
    }
}


/*============================================================================
*    Local Void function tilt
*
*    ETR on a tilted surface
*----------------------------------------------------------------------------*/
static void tilt(struct posdata *pdat)
{
    float ca;			/* cosine of the solar azimuth angle */
    float cp;			/* cosine of the panel aspect */
    float ct;			/* cosine of the panel tilt */
    float sa;			/* sine of the solar azimuth angle */
    float sp;			/* sine of the panel aspect */
    float st;			/* sine of the panel tilt */
    float sz;			/* sine of the refraction corrected solar zenith angle */


    /* 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);
    pdat->cosinc = pdat->coszen * ct + sz * st * (ca * cp + sa * sp);

    if (pdat->cosinc > 0.0)
	pdat->etrtilt = pdat->etrn * pdat->cosinc;
    else
	pdat->etrtilt = 0.0;

}


/*============================================================================
*    Void function S_decode
*
*    This function decodes the error codes from S_solpos return value
*
*    Requires the long integer return value from S_solpos
*
*    Returns descriptive text to stderr
*----------------------------------------------------------------------------*/
void S_decode(long code, struct posdata *pdat)
{
    if (code & (1L << S_YEAR_ERROR))
	G_warning(_("S_decode ==> Please fix the year: %d [1950-2050]"),
		  pdat->year);
    if (code & (1L << S_MONTH_ERROR))
	G_warning(_("S_decode ==> Please fix the month: %d"), pdat->month);
    if (code & (1L << S_DAY_ERROR))
	G_warning(_("S_decode ==> Please fix the day-of-month: %d"),
		  pdat->day);
    if (code & (1L << S_DOY_ERROR))
	G_warning(_("S_decode ==> Please fix the day-of-year: %d"),
		  pdat->daynum);
    if (code & (1L << S_HOUR_ERROR))
	G_warning(_("S_decode ==> Please fix the hour: %d"), pdat->hour);
    if (code & (1L << S_MINUTE_ERROR))
	G_warning(_("S_decode ==> Please fix the minute: %d"), pdat->minute);
    if (code & (1L << S_SECOND_ERROR))
	G_warning(_("S_decode ==> Please fix the second: %d"), pdat->second);
    if (code & (1L << S_TZONE_ERROR))
	G_warning(_("S_decode ==> Please fix the time zone: %f"),
		  pdat->timezone);
    if (code & (1L << S_INTRVL_ERROR))
	G_warning(_("S_decode ==> Please fix the interval: %d"),
		  pdat->interval);
    if (code & (1L << S_LAT_ERROR))
	G_warning(_("S_decode ==> Please fix the latitude: %f"),
		  pdat->latitude);
    if (code & (1L << S_LON_ERROR))
	G_warning(_("S_decode ==> Please fix the longitude: %f"),
		  pdat->longitude);
    if (code & (1L << S_TEMP_ERROR))
	G_warning(_("S_decode ==> Please fix the temperature: %f"),
		  pdat->temp);
    if (code & (1L << S_PRESS_ERROR))
	G_warning(_("S_decode ==> Please fix the pressure: %f"), pdat->press);
    if (code & (1L << S_TILT_ERROR))
	G_warning(_("S_decode ==> Please fix the tilt: %f"), pdat->tilt);
    if (code & (1L << S_ASPECT_ERROR))
	G_warning(_("S_decode ==> Please fix the aspect: %f"), pdat->aspect);
    if (code & (1L << S_SBWID_ERROR))
	G_warning(_("S_decode ==> Please fix the shadowband width: %f"),
		  pdat->sbwid);
    if (code & (1L << S_SBRAD_ERROR))
	G_warning(_("S_decode ==> Please fix the shadowband radius: %f"),
		  pdat->sbrad);
    if (code & (1L << S_SBSKY_ERROR))
	G_warning(_("S_decode ==> Please fix the shadowband sky factor: %f"),
		  pdat->sbsky);
}