SOLPOS.C

Go to the documentation of this file.

/*============================================================================
*    Contains:
*        S_date       (computes month/day from day number)
*            INPUTS:     S_year, S_daynum
*            OUTPUTS:    S_year, S_month, S_day
*
*        S_solpos     (computes solar position and intensity
*                      from time and place)
*            INPUTS:     S_month, S_day, S_latitude, S_longitude, S_timezone
*            OPTIONAL:   S_year, S_hour, S_minute, S_second,
*                            S_press, S_temp, S_tilt, S_aspect
*            OUTPUTS:    EVERY variable beginning with "S_"
*                            (defined beneath this comment block)
*
*    Usage:
*         In calling program, just after other 'includes', insert:
*
*              #include "solpos.h"
*
*    Martin Rymes
*    National Renewable Energy Laboratory
*    25 March 1998
*
*    27 April 1999 REVISION:  Corrected leap year in S_date.
*    14 January 2000 REVISION:  Corrected leap year in S_solpos.
*----------------------------------------------------------------------------*/
#include <math.h>
#include <string.h>
#include <stdio.h>



int   S_day       =  -99;     /* Day of month (May 27 = 27, etc.) */
int   S_daynum    = -999;     /* Day number (day of year; Feb 1 = 32 ) */
int   S_hour      =   12;     /* Hour of day, 0 - 23 */
int   S_minute    =    0;     /* Minute of hour, 0 - 59 */
int   S_month     =  -99;     /* Month number (Jan = 1, Feb = 2, etc.) */
int   S_second    =    0;     /* Second of minute, 0 - 59 */
int   S_year      = 2001;     /* 4-digit year (2-digit is assumed 19xx) */

double S_amass;                /* Relative optical airmass */
double S_ampress;              /* Pressure-corrected airmass */
double S_aspect    = 180.0;    /* Azimuth of panel surface (direction it
                                    faces) N=0, E=90, S=180, W=270 */
double S_azim;                 /* Solar azimuth angle:  N=0, E=90, S=180, W=270 */
double S_cosinc;               /* Cosine of solar incidence angle on panel */
double S_dayang;               /* Day angle (daynum*360/year-length) degrees */
double S_declin;               /* Declination--zenith angle of solar noon
                                    at equator, degrees NORTH */
double S_eclong;               /* Ecliptic longitude, degrees */
double S_ecobli;               /* Obliquity of ecliptic */
double S_ectime;               /* Time of ecliptic calculations */
double S_elevref;              /* Solar elevation angle,
                                    deg. from horizon, refracted */
double S_eqntim;               /* Equation of time (TST - LMT), minutes */
double S_erv;                  /* Earth radius vector
                                    (multiplied to solar constant) */
double S_etr;                  /* Extraterrestrial (top-of-atmosphere) W/sq m
                                    global horizontal solar irradiance */
double S_etrn;                 /* Extraterrestrial (top-of-atmosphere) W/sq m
                                    direct normal solar irradiance */
double S_etrtilt;              /* Extraterrestrial (top-of-atmosphere) W/sq m
                                    global irradiance on a tilted surface */
double S_gmst;                 /* Greenwich mean sidereal time, hours */
double S_hrang;                /* Hour angle--hour of sun from solar noon,
                                    degrees WEST */
double S_julday;               /* Julian Day of 1 JAN 2000 minus 2,400,000 days
                                    (in order to regain single precision) */
double S_latitude  =  -99.0;   /* Latitude, degrees north (south negative) */
double S_longitude = -999.0;   /* Longitude, degrees east (west negative) */
double S_lmst;                 /* Local mean sidereal time, degrees */
double S_mnanom;               /* Mean anomaly, degrees */
double S_mnlong;               /* Mean longitude, degrees */
double S_rascen;               /* Right ascension, degrees */
double S_press     = 1013.0;   /* Surface pressure, millibars */
double S_prime;                /* Factor that normalizes Kt, Kn, etc. */
double S_sbcf;                 /* Shadow-band correction factor */
double S_solcon    = 1367.0;   /* Solar constant, 1367 W/sq m */
double S_ssha;                 /* Sunset(/rise) hour angle, degrees */
double S_sunrise;              /* Sunrise time, minutes from midnight,
                                    local, no refraction */
double S_sunset;               /* Sunset time, minutes from midnight,
                                    local, no refraction */
double S_temp      =    0.0;   /* Ambient dry-bulb temperature, degrees C */
double S_tilt      =    0.0;   /* Degrees tilt from horizontal of panel */
double S_timezone  =  -99.0;   /* Time zone, east (west negative).
                                    USA:  Mountain = -7, Central = -6, etc. */
double S_tst;                  /* True solar time, minutes from midnight */
double S_tstfix;               /* True solar time - local standard time */
double S_unprime;              /* Factor that denormalizes Kt', Kn', etc. */
double S_utime;                /* Universal (Greenwich) standard time */
double S_zenetr;               /* Solar zenith angle,
                                    no atmospheric correction (= ETR) */
double S_zenref;               /* Solar zenith angle,
                                    deg. from zenith, refracted */

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
* Temporary variables used only in this file:
*
*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
static int  ndy;   /* Temporary S_daynum */
static int  leap;  /* Leap year counter */
static int  imon;  /* Month (month_days) array counter */
static int  lenyr[2] = { 365, 366 };
                   /* length of year; 0 = non-leap, 1 = leap year */
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 */

static double bottom;      /* denomerator (bottom) of the fraction */
static double c2;          /* cosine of d2 */
static double ca;          /* cosine of the solar azimuth angle */
static double cd;          /* cosine of the day angle or declination */
static double cdcl;        /* ( cd * cl ) */
static double ce;          /* cosine of the solar elevation */
static double cecl;        /* ( ce * cl ) */
static double ch;          /* cosine of the hour angle */
static double cl;          /* cosine of the latitude */
static double cp;          /* cosine of the panel aspect */
static double cssha;       /* cosine of the sunset hour angle */
static double ct;          /* cosine of the panel tilt */
static double cz;          /* cosine of the solar zenith angle */
static double d2;          /* Double S_dayang */
static double degrad = 57.295779513;
                          /* converts from radians to degrees */
static double delta;       /* difference between current year and 1949 */
static double elev;        /* temporary elevation angle */
static double p;           /* used to compute S_sbcf */
static double prestemp;    /* temporary pressure/temperature correction */
static double raddeg = 0.0174532925;
                          /* converts from degrees to radians */
static double refcor;      /* temporary refraction correction */
static double s2;          /* sine of d2 */
static double sa;          /* sine of the solar azimuth angle */
static double sd;          /* sine of the day angle or declination */
static double se;          /* sine of the solar elevation */
static double sl;          /* sine of the latitude */
static double sp;          /* sine of the panel aspect */
static double st;          /* sine of the panel tilt */
static double sz;          /* sine of the solar zenith angle */
static double tanelev;     /* tangent of the solar elevation angle */
static double t1;          /* used to compute S_sbcf */
static double t2;          /* used to compute S_sbcf */
static double top;         /* numerator (top) of the fraction */

/*============================================================================
*    Integer function S_date, adapted from the VAX solar libraries
*
*    This function computes the month/day from the day number.
*
*    Requires:
*        Year and day number:
*            S_year           NOTE:  2-digit year "xx" is assumed "19xx"
*                             DEFAULT S_year = 2001
*            S_daynum         RANGE: -1 to a large positive integer
*
*    Returns:
*        Year, month, day:
*            S_year
*            S_month
*            S_day
*----------------------------------------------------------------------------*/
int S_date (void)
{
    /* We will allow -1 (in order to adjust to the previous year) but no
less */
    if ( S_daynum < -1 ) {
        printf ( "S_date => Please fix the day number: %10d\n", S_daynum );
        return ( -1 );
    }

    /* If a 2-digit year "xx" is given, assume "19xx" */
    if ( S_year < 100 )
        S_year += 1900;

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

    /* If S_daynum is a large number, it is in a future year */
    ndy = S_daynum;
    while ( ndy > lenyr[leap] ) {
        ++S_year;
        ndy -= lenyr[leap];
        if ( ((S_year % 4) == 0) && 
             ( ((S_year % 100) != 0) || ((S_year % 400) == 0) ) )
            leap = 1;
        else
            leap = 0;
    }

    /* A non-positive day number indicates the previous year */
    if ( ndy <= 0 ) {
        --S_year;
        leap  =  0;
        ndy  += 365;
    }

    /* Now that the year is fixed, find the month */
    imon = 12;
    while ( ndy <= month_days [leap][imon] )
        --imon;

    /* Set the month and day of month */
    S_month = imon;
    S_day   = ndy - month_days[leap][imon];

    return ( 0 );
}
/*============================================================================
*    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:
*        Date and time:
*            S_year      DEFAULT 2001
*            S_month
*            S_day
*            S_hour      DEFAULT 12
*            S_minute    DEFAULT    0
*            S_second    DEFAULT 0
*        Location:
*            S_latitude
*            S_longitude
*        Location/time adjuster:
*            S_timezone
*        Atmospheric pressure and temperature:
*            S_press     DEFAULT 1013.0 mb
*            S_temp      DEFAULT 10.0 degrees C
*        Tilt of flat surface that receives solar energy:
*            S_aspect    DEFAULT 180 (South)
*            S_tilt      DEFAULT 0 (Horizontal)
*
*    Returns:
*        everything defined at the top of this listing.
*----------------------------------------------------------------------------*/
int S_solpos (void)
{
     /* No absurd dates, please. */
    if ( (S_month < 1) || (S_month > 12) )
    {
        printf ( "S_solpos => Please fix the month: %10d\n", S_month );
        return ( -1 );
    }
    else if ( (S_day < 1) || (S_day > 31 ) )
    {
        printf ( "S_solpos => Please fix the day of month: %10d\n", S_day );
        return ( -1 );
    }

    /* No absurd times, please. */
    if ( (S_hour < 0) || (S_hour > 23) )
    {
        printf ( "S_solpos => Please fix the hour: %10d\n", S_hour );
        return ( -1 );
    }
    else if ( (S_minute < 0) || (S_minute > 59) )
    {
        printf ( "S_solpos => Please fix the minute: %10d\n", S_minute );
        return ( -1 );
    }
    else if ( (S_second < 0) || (S_second > 59) )
    {
        printf ( "S_solpos => Please fix the second: %10d\n", S_second );
        return ( -1 );
    }
    else if ( (S_timezone < -12.0) || (S_timezone > 12.0) )
    {
        printf ( "S_solpos => Please fix the time zone: %10.4f\n",
S_timezone );
        return ( -1 );
    }
    /* No absurd locations, please. */
    if ( fabs (S_longitude) > 180.0 ) {
        printf ( "S_solpos => Please fix the longitude: %10.4f\n",
S_longitude );
        return ( -1 );
    }
    else if ( fabs (S_latitude) > 90.0 ) {
        printf ( "S_solpos => Please fix the latitude: %10.4f\n",
S_latitude );
        return ( -1 );
    }

    /* No silly temperatures or pressures, please. */
    if ( fabs (S_temp) > 100.0 ) {
        printf ( "S_solpos => Please fix the temperature: %10.4f\n", S_temp );
        return ( -1 );
    }
    else if ( (S_press < 0.0) || (S_press > 2000.0) ) {
        printf ( "S_solpos => Please fix the pressure: %10.4f\n", S_press );
        return ( -1 );
    }

    /* Day number */
    /* (convert day-of-month to day-of-year (non-leap)) */
    S_daynum = S_day + month_days[0][S_month];

    /* (if no year was entered, give it an arbitrary 2001.  Also, turn
        2-digit years into 19xx.  Note that if the year > 1999, it had
        better be 4-digit) */
    if ( S_year < 0 )
        S_year  = 2001;
    else if ( S_year < 100 )
        S_year += 1900;

    /* (now adjust for leap year) */
    if ( ((S_year % 4) == 0) &&
       ( ((S_year % 100) != 0) || ((S_year % 400) == 0) ) &&
             (S_month > 2) )
        S_daynum += 1;

     /* Day angle */
        /*  Iqbal, M.  1983.  An Introduction to Solar Radiation.
            Academic Press, NY., page 3 */
     S_dayang = 360.0 * ( S_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 * S_dayang);
    cd     = cos (raddeg * S_dayang);
    d2     = 2.0 * S_dayang;
    c2     = cos (raddeg * d2);
    s2     = sin (raddeg * d2);

    S_erv  = 1.000110 + 0.034221 * cd + 0.001280 * sd;
    S_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. */
    S_utime = S_hour * 3600.0 + S_minute * 60.0 + S_second;
    S_utime = S_utime / 3600.0 - S_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. */

    delta    = S_year - 1949;
    leap     = (int) ( delta / 4.0 );
    S_julday = 32916.5 + delta * 365.0 + leap + S_daynum + S_utime / 24.0;

    if ( ( (S_year % 100) == 0 ) && ( (S_year % 400) != 0) )
        --S_julday;

    /* 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. */
    S_ectime = S_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. */
    S_mnlong  = 280.460 + 0.9856474 * S_ectime;

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
    S_mnlong -= 360.0 * (int) ( S_mnlong / 360.0 );
    if ( S_mnlong < 0.0 )
        S_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. */
    S_mnanom  = 357.528 + 0.9856003 * S_ectime;

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
    S_mnanom -= 360.0 * (int) ( S_mnanom / 360.0 );
    if ( S_mnanom < 0.0 )
        S_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. */
    S_eclong  = S_mnlong + 1.915 * sin ( S_mnanom * raddeg ) +
                    0.020 * sin ( 2.0 * S_mnanom * raddeg );

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
    S_eclong -= 360.0 * (int) ( S_eclong / 360.0 );
    if ( S_eclong < 0.0 )
        S_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. */
    S_ecobli = 23.439 + 4.0e-07 * S_ectime;

    /* Declination */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
    S_declin = degrad * asin ( sin (S_ecobli * raddeg) *
                               sin (S_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 * S_ecobli ) * sin ( raddeg * S_eclong );
    bottom   =  cos ( raddeg * S_eclong );

    S_rascen =  degrad * atan2 ( top, bottom );

    /* (make it a positive angle) */
    if ( S_rascen < 0.0 )
        S_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. */
    S_gmst  = 6.697375 + 0.0657098242 * S_ectime + S_utime;

    /* (dump the multiples of 24, so the answer is between 0 and 24) */
    S_gmst -= 24.0 * (int) ( S_gmst / 24.0 );
    if ( S_gmst < 0.0 )
        S_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. */
    S_lmst  = S_gmst * 15.0 + S_longitude;

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
    S_lmst -= 360.0 * (int) ( S_lmst / 360.0 );
    if ( S_lmst < 0.)
        S_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. */
    S_hrang = S_lmst - S_rascen;

    /* (force it between -180 and 180 degrees) */
    if ( S_hrang < -180.0 )
        S_hrang += 360.0;
    else if ( S_hrang > 180.0 )
        S_hrang -= 360.0;

    /* ETR solar zenith angle */
        /*  Iqbal, M.  1983.  An Introduction to Solar Radiation.
            Academic Press, NY., page 15 */
    cd = cos ( raddeg * S_declin );
    ch = cos ( raddeg * S_hrang );
    cl = cos ( raddeg * S_latitude );
    sd = sin ( raddeg * S_declin );
    sl = sin ( raddeg * S_latitude );

    cz = sd * sl + cd * cl * ch;

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

    S_zenetr   = acos ( cz ) * degrad;

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

    /* Sunset hour angle */
    /*  Iqbal, M.  1983.  An Introduction to Solar Radiation.
        Academic Press, NY., page 16 */
    S_ssha  = 90.0;
    cdcl    = cd * cl;

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

    /* This keeps the cosine from blowing on roundoff */
    if ( cssha < -1.0  )
        S_ssha = 180.0;
    else if ( cssha > 1.0 )
        S_ssha = 0.0;
    else
        S_ssha = degrad * acos ( cssha );
    }

    /* Shadowband correction factor */
    /*  Drummond, A. J.  1956.  A contribution to absolute pyrheliometry.
        Q. J. R. Meteorol. Soc. 82, pp. 481-493 */
    p      = 0.1526 * pow (cd,3);
    t1     = sl * sd * S_ssha * raddeg;
    t2     = cl * cd * sin ( S_ssha * raddeg );
    S_sbcf = 0.04 + 1.0 / ( 1.0 - p * ( t1 + t2 ) );

    /* TST -> True Solar Time = local standard time + TSTfix */
        /*  Iqbal, M.  1983.  An Introduction to Solar Radiation.
            Academic Press, NY., page 13 */
    S_tst    = ( 180.0 + S_hrang ) * 4.0;
    S_tstfix = S_tst - S_hour * 60.0 - S_minute - S_second / 60.0;
    S_eqntim = S_tstfix + 60.0 * S_timezone - 4.0 * S_longitude;

    /* Sunrise and sunset times */
    if ( S_ssha <= 1.0 ) {
        S_sunrise   =  2999.0;
        S_sunset    = -2999.0;
    }
    else if ( S_ssha >= 179.0 ) {
        S_sunrise   = -2999.0;
        S_sunset    =  2999.0;
    }
    else {
        S_sunrise   = 720.0 - 4.0 * S_ssha - S_tstfix;
        S_sunset    = 720.0 + 4.0 * S_ssha - S_tstfix;
    }

    /* Solar azimuth angle */
        /*  Iqbal, M.  1983.  An Introduction to Solar Radiation.
            Academic Press, NY., page 15 */
    ce         = sin ( raddeg * S_zenetr );
    se         = cos ( raddeg * S_zenetr );

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

        S_azim = 180.0 - acos ( ca ) * degrad;
        if ( S_hrang > 0 )
            S_azim  = 360.0 - S_azim;
    }

    /* Refraction correction */
        /*  Zimmerman, John C.  1981.  Sun-pointing programs and their
accuracy.
            SAND81-0761, Experimental Systems Operation Division 4721,
            Sandia National Laboratories, Albuquerque, NM. */

    elev = 90.0 - S_zenetr;

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

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

        prestemp    = ( S_press * 283.0 ) / ( 1013.0 * ( 273.0 + S_temp ) );
        refcor     *= prestemp / 3600.0;
    }

    /* Refracted solar elevation angle */
    S_elevref = elev + refcor;

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

    /* Refracted solar zenith angle */
    S_zenref  = 90. - S_elevref;

    /* Airmass */
        /*  Kasten, F. and Young, A.  1989.  Revised optical air mass
tables and
            approximation formula.  Applied Optics 28 (22), pp. 4735-4738 */
    cz = cos ( raddeg * S_zenref );
    if ( S_zenref > 93.0 )
        S_amass = -1.0;
    else
        S_amass = 1.0 / ( cz + 0.50572 * pow ((96.07995 - S_zenref),-1.6364) );

    S_ampress   = S_amass * S_press / 1013.0;

    /* 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 */
    S_unprime = 1.031 * exp ( -1.4 / ( 0.9 + 9.4 / S_amass ) ) + 0.1;
    S_prime   = 1.0 / S_unprime;

    /* Cosine of the angle between the sun and a tipped flat surface,
       useful for calculating solar energy on tilted surfaces */
    ca        = cos ( raddeg * S_azim );
    cp        = cos ( raddeg * S_aspect );
    ct        = cos ( raddeg * S_tilt );
    sa        = sin ( raddeg * S_azim );
    sp        = sin ( raddeg * S_aspect );
    st        = sin ( raddeg * S_tilt );
    sz        = sin ( raddeg * S_zenref );
    S_cosinc  = cz * ct + sz * st * ( ca * cp + sa * sp );

    /* Extraterrestrial (top-of-atmosphere) solar irradiance */
    if ( cz > 0.0 ) {
        S_etrn = S_solcon * S_erv;
        S_etr  = S_etrn * cz;
    }
    else {
        S_etrn = 0.0;
        S_etr  = 0.0;
    }
 
    if ( S_cosinc > 0.0 )
        S_etrtilt = S_etrn * S_cosinc;
    else
        S_etrtilt = 0.0;

    return ( 0 );
}