/*============================================================================ * 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 "solpos00.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) *----------------------------------------------------------------------------*/ #include #include #include #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 */ 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) ) fprintf(stderr, "S_decode ==> Please fix the year: %d [1950-2050]\n", pdat->year); if ( code & (1L << S_MONTH_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the month: %d\n", pdat->month); if ( code & (1L << S_DAY_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the day-of-month: %d\n", pdat->day); if ( code & (1L << S_DOY_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the day-of-year: %d\n", pdat->daynum); if ( code & (1L << S_HOUR_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the hour: %d\n", pdat->hour); if ( code & (1L << S_MINUTE_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the minute: %d\n", pdat->minute); if ( code & (1L << S_SECOND_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the second: %d\n", pdat->second); if ( code & (1L << S_TZONE_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the time zone: %f\n", pdat->timezone); if ( code & (1L << S_INTRVL_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the interval: %d\n", pdat->interval); if ( code & (1L << S_LAT_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the latitude: %f\n", pdat->latitude); if ( code & (1L << S_LON_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the longitude: %f\n", pdat->longitude); if ( code & (1L << S_TEMP_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the temperature: %f\n", pdat->temp); if ( code & (1L << S_PRESS_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the pressure: %f\n", pdat->press); if ( code & (1L << S_TILT_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the tilt: %f\n", pdat->tilt); if ( code & (1L << S_ASPECT_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the aspect: %f\n", pdat->aspect); if ( code & (1L << S_SBWID_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the shadowband width: %f\n", pdat->sbwid); if ( code & (1L << S_SBRAD_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the shadowband radius: %f\n", pdat->sbrad); if ( code & (1L << S_SBSKY_ERROR) ) fprintf(stderr, "S_decode ==> Please fix the shadowband sky factor: %f\n", pdat->sbsky); }