James Cobb
This code is imported direct from http://rredc.nrel.gov/solar/codesandalgorithms/solpos/ Copyrights are retained in the code for each algorithm used. The library calculates the apparent solar position and intensity (theoretical maximum solar energy) based on the date, time, and location on Earth. The calculated data can be used to predict solar radiation, to be used in meteorological, solar energy and irrigation applications.
Embed:
(wiki syntax)
Show/hide line numbers
solpos.cpp
00001 /*============================================================================ 00002 * Contains: 00003 * S_solpos (computes solar position and intensity 00004 * from time and place) 00005 * 00006 * INPUTS: (via posdata struct) year, daynum, hour, 00007 * minute, second, latitude, longitude, timezone, 00008 * intervl 00009 * OPTIONAL: (via posdata struct) month, day, press, temp, tilt, 00010 * aspect, function 00011 * OUTPUTS: EVERY variable in the struct posdata 00012 * (defined in solpos.h) 00013 * 00014 * NOTE: Certain conditions exist during which some of 00015 * the output variables are undefined or cannot be 00016 * calculated. In these cases, the variables are 00017 * returned with flag values indicating such. In other 00018 * cases, the variables may return a realistic, though 00019 * invalid, value. These variables and the flag values 00020 * or invalid conditions are listed below: 00021 * 00022 * amass -1.0 at zenetr angles greater than 93.0 00023 * degrees 00024 * ampress -1.0 at zenetr angles greater than 93.0 00025 * degrees 00026 * azim invalid at zenetr angle 0.0 or latitude 00027 * +/-90.0 or at night 00028 * elevetr limited to -9 degrees at night 00029 * etr 0.0 at night 00030 * etrn 0.0 at night 00031 * etrtilt 0.0 when cosinc is less than 0 00032 * prime invalid at zenetr angles greater than 93.0 00033 * degrees 00034 * sretr +/- 2999.0 during periods of 24 hour sunup or 00035 * sundown 00036 * ssetr +/- 2999.0 during periods of 24 hour sunup or 00037 * sundown 00038 * ssha invalid at the North and South Poles 00039 * unprime invalid at zenetr angles greater than 93.0 00040 * degrees 00041 * zenetr limited to 99.0 degrees at night 00042 * 00043 * S_init (optional initialization for all input parameters in 00044 * the posdata struct) 00045 * INPUTS: struct posdata* 00046 * OUTPUTS: struct posdata* 00047 * 00048 * (Note: initializes the required S_solpos INPUTS above 00049 * to out-of-bounds conditions, forcing the user to 00050 * supply the parameters; initializes the OPTIONAL 00051 * S_solpos inputs above to nominal values.) 00052 * 00053 * S_decode (optional utility for decoding the S_solpos return code) 00054 * INPUTS: long integer S_solpos return value, struct posdata* 00055 * OUTPUTS: text to stderr 00056 * 00057 * Usage: 00058 * In calling program, just after other 'includes', insert: 00059 * 00060 * #include "solpos00.h" 00061 * 00062 * Function calls: 00063 * S_init(struct posdata*) [optional] 00064 * . 00065 * . 00066 * [set time and location parameters before S_solpos call] 00067 * . 00068 * . 00069 * int retval = S_solpos(struct posdata*) 00070 * S_decode(int retval, struct posdata*) [optional] 00071 * (Note: you should always look at the S_solpos return 00072 * value, which contains error codes. S_decode is one option 00073 * for examining these codes. It can also serve as a 00074 * template for building your own application-specific 00075 * decoder.) 00076 * 00077 * Martin Rymes 00078 * National Renewable Energy Laboratory 00079 * 25 March 1998 00080 * 00081 * 27 April 1999 REVISION: Corrected leap year in S_date. 00082 * 13 January 2000 REVISION: SMW converted to structure posdata parameter 00083 * and subdivided into functions. 00084 * 01 February 2001 REVISION: SMW corrected ecobli calculation 00085 * (changed sign). Error is small (max 0.015 deg 00086 * in calculation of declination angle) 00087 *----------------------------------------------------------------------------*/ 00088 #include <math.h> 00089 #include <string.h> 00090 #include <stdio.h> 00091 #include "solpos00.h" 00092 00093 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00094 * 00095 * Structures defined for this module 00096 * 00097 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 00098 struct trigdata /* used to pass calculated values locally */ 00099 { 00100 float cd; /* cosine of the declination */ 00101 float ch; /* cosine of the hour angle */ 00102 float cl; /* cosine of the latitude */ 00103 float sd; /* sine of the declination */ 00104 float sl; /* sine of the latitude */ 00105 }; 00106 00107 00108 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00109 * 00110 * Temporary global variables used only in this file: 00111 * 00112 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 00113 static int month_days[2][13] = { { 0, 0, 31, 59, 90, 120, 151, 00114 181, 212, 243, 273, 304, 334 }, 00115 { 0, 0, 31, 60, 91, 121, 152, 00116 182, 213, 244, 274, 305, 335 } }; 00117 /* cumulative number of days prior to beginning of month */ 00118 00119 static float degrad = 57.295779513; /* converts from radians to degrees */ 00120 static float raddeg = 0.0174532925; /* converts from degrees to radians */ 00121 00122 /*============================================================================ 00123 * Local function prototypes 00124 ============================================================================*/ 00125 static long int validate ( struct posdata *pdat); 00126 static void dom2doy( struct posdata *pdat ); 00127 static void doy2dom( struct posdata *pdat ); 00128 static void geometry ( struct posdata *pdat ); 00129 static void zen_no_ref ( struct posdata *pdat, struct trigdata *tdat ); 00130 static void ssha( struct posdata *pdat, struct trigdata *tdat ); 00131 static void sbcf( struct posdata *pdat, struct trigdata *tdat ); 00132 static void tst( struct posdata *pdat ); 00133 static void srss( struct posdata *pdat ); 00134 static void sazm( struct posdata *pdat, struct trigdata *tdat ); 00135 static void refrac( struct posdata *pdat ); 00136 static void amass( struct posdata *pdat ); 00137 static void prime( struct posdata *pdat ); 00138 static void etr( struct posdata *pdat ); 00139 static void tilt( struct posdata *pdat ); 00140 static void localtrig( struct posdata *pdat, struct trigdata *tdat ); 00141 00142 /*============================================================================ 00143 * Long integer function S_solpos, adapted from the VAX solar libraries 00144 * 00145 * This function calculates the apparent solar position and the 00146 * intensity of the sun (theoretical maximum solar energy) from 00147 * time and place on Earth. 00148 * 00149 * Requires (from the struct posdata parameter): 00150 * Date and time: 00151 * year 00152 * daynum (requirement depends on the S_DOY switch) 00153 * month (requirement depends on the S_DOY switch) 00154 * day (requirement depends on the S_DOY switch) 00155 * hour 00156 * minute 00157 * second 00158 * interval DEFAULT 0 00159 * Location: 00160 * latitude 00161 * longitude 00162 * Location/time adjuster: 00163 * timezone 00164 * Atmospheric pressure and temperature: 00165 * press DEFAULT 1013.0 mb 00166 * temp DEFAULT 10.0 degrees C 00167 * Tilt of flat surface that receives solar energy: 00168 * aspect DEFAULT 180 (South) 00169 * tilt DEFAULT 0 (Horizontal) 00170 * Function Switch (codes defined in solpos.h) 00171 * function DEFAULT S_ALL 00172 * 00173 * Returns (via the struct posdata parameter): 00174 * everything defined in the struct posdata in solpos.h. 00175 *----------------------------------------------------------------------------*/ 00176 long S_solpos (struct posdata *pdat) 00177 { 00178 long int retval; 00179 00180 struct trigdata trigdat, *tdat; 00181 00182 tdat = &trigdat; /* point to the structure */ 00183 00184 /* initialize the trig structure */ 00185 tdat->sd = -999.0; /* flag to force calculation of trig data */ 00186 tdat->cd = 1.0; 00187 tdat->ch = 1.0; /* set the rest of these to something safe */ 00188 tdat->cl = 1.0; 00189 tdat->sl = 1.0; 00190 00191 if ((retval = validate ( pdat )) != 0) /* validate the inputs */ 00192 return retval; 00193 00194 00195 if ( pdat->function & L_DOY ) 00196 doy2dom( pdat ); /* convert input doy to month-day */ 00197 else 00198 dom2doy( pdat ); /* convert input month-day to doy */ 00199 00200 if ( pdat->function & L_GEOM ) 00201 geometry( pdat ); /* do basic geometry calculations */ 00202 00203 if ( pdat->function & L_ZENETR ) /* etr at non-refracted zenith angle */ 00204 zen_no_ref( pdat, tdat ); 00205 00206 if ( pdat->function & L_SSHA ) /* Sunset hour calculation */ 00207 ssha( pdat, tdat ); 00208 00209 if ( pdat->function & L_SBCF ) /* Shadowband correction factor */ 00210 sbcf( pdat, tdat ); 00211 00212 if ( pdat->function & L_TST ) /* true solar time */ 00213 tst( pdat ); 00214 00215 if ( pdat->function & L_SRSS ) /* sunrise/sunset calculations */ 00216 srss( pdat ); 00217 00218 if ( pdat->function & L_SOLAZM ) /* solar azimuth calculations */ 00219 sazm( pdat, tdat ); 00220 00221 if ( pdat->function & L_REFRAC ) /* atmospheric refraction calculations */ 00222 refrac( pdat ); 00223 00224 if ( pdat->function & L_AMASS ) /* airmass calculations */ 00225 amass( pdat ); 00226 00227 if ( pdat->function & L_PRIME ) /* kt-prime/unprime calculations */ 00228 prime( pdat ); 00229 00230 if ( pdat->function & L_ETR ) /* ETR and ETRN (refracted) */ 00231 etr( pdat ); 00232 00233 if ( pdat->function & L_TILT ) /* tilt calculations */ 00234 tilt( pdat ); 00235 00236 return 0; 00237 } 00238 00239 00240 /*============================================================================ 00241 * Void function S_init 00242 * 00243 * This function initiates all of the input parameters in the struct 00244 * posdata passed to S_solpos(). Initialization is either to nominal 00245 * values or to out of range values, which forces the calling program to 00246 * specify parameters. 00247 * 00248 * NOTE: This function is optional if you initialize ALL input parameters 00249 * in your calling code. Note that the required parameters of date 00250 * and location are deliberately initialized out of bounds to force 00251 * the user to enter real-world values. 00252 * 00253 * Requires: Pointer to a posdata structure, members of which are 00254 * initialized. 00255 * 00256 * Returns: Void 00257 *----------------------------------------------------------------------------*/ 00258 void S_init(struct posdata *pdat) 00259 { 00260 pdat->day = -99; /* Day of month (May 27 = 27, etc.) */ 00261 pdat->daynum = -999; /* Day number (day of year; Feb 1 = 32 ) */ 00262 pdat->hour = -99; /* Hour of day, 0 - 23 */ 00263 pdat->minute = -99; /* Minute of hour, 0 - 59 */ 00264 pdat->month = -99; /* Month number (Jan = 1, Feb = 2, etc.) */ 00265 pdat->second = -99; /* Second of minute, 0 - 59 */ 00266 pdat->year = -99; /* 4-digit year */ 00267 pdat->interval = 0; /* instantaneous measurement interval */ 00268 pdat->aspect = 180.0; /* Azimuth of panel surface (direction it 00269 faces) N=0, E=90, S=180, W=270 */ 00270 pdat->latitude = -99.0; /* Latitude, degrees north (south negative) */ 00271 pdat->longitude = -999.0; /* Longitude, degrees east (west negative) */ 00272 pdat->press = 1013.0; /* Surface pressure, millibars */ 00273 pdat->solcon = 1367.0; /* Solar constant, 1367 W/sq m */ 00274 pdat->temp = 15.0; /* Ambient dry-bulb temperature, degrees C */ 00275 pdat->tilt = 0.0; /* Degrees tilt from horizontal of panel */ 00276 pdat->timezone = -99.0; /* Time zone, east (west negative). */ 00277 pdat->sbwid = 7.6; /* Eppley shadow band width */ 00278 pdat->sbrad = 31.7; /* Eppley shadow band radius */ 00279 pdat->sbsky = 0.04; /* Drummond factor for partly cloudy skies */ 00280 pdat->function = S_ALL; /* compute all parameters */ 00281 } 00282 00283 00284 /*============================================================================ 00285 * Local long int function validate 00286 * 00287 * Validates the input parameters 00288 *----------------------------------------------------------------------------*/ 00289 static long int validate ( struct posdata *pdat) 00290 { 00291 00292 long int retval = 0; /* start with no errors */ 00293 00294 /* No absurd dates, please. */ 00295 if ( pdat->function & L_GEOM ) 00296 { 00297 if ( (pdat->year < 1950) || (pdat->year > 2050) ) /* limits of algoritm */ 00298 retval |= (1L << S_YEAR_ERROR); 00299 if ( !(pdat->function & S_DOY) && ((pdat->month < 1) || (pdat->month > 12))) 00300 retval |= (1L << S_MONTH_ERROR); 00301 if ( !(pdat->function & S_DOY) && ((pdat->day < 1) || (pdat->day > 31)) ) 00302 retval |= (1L << S_DAY_ERROR); 00303 if ( (pdat->function & S_DOY) && ((pdat->daynum < 1) || (pdat->daynum > 366)) ) 00304 retval |= (1L << S_DOY_ERROR); 00305 00306 /* No absurd times, please. */ 00307 if ( (pdat->hour < 0) || (pdat->hour > 24) ) 00308 retval |= (1L << S_HOUR_ERROR); 00309 if ( (pdat->minute < 0) || (pdat->minute > 59) ) 00310 retval |= (1L << S_MINUTE_ERROR); 00311 if ( (pdat->second < 0) || (pdat->second > 59) ) 00312 retval |= (1L << S_SECOND_ERROR); 00313 if ( (pdat->hour == 24) && (pdat->minute > 0) ) /* no more than 24 hrs */ 00314 retval |= ( (1L << S_HOUR_ERROR) | (1L << S_MINUTE_ERROR) ); 00315 if ( (pdat->hour == 24) && (pdat->second > 0) ) /* no more than 24 hrs */ 00316 retval |= ( (1L << S_HOUR_ERROR) | (1L << S_SECOND_ERROR) ); 00317 if ( fabs (pdat->timezone) > 12.0 ) 00318 retval |= (1L << S_TZONE_ERROR); 00319 if ( (pdat->interval < 0) || (pdat->interval > 28800) ) 00320 retval |= (1L << S_INTRVL_ERROR); 00321 00322 /* No absurd locations, please. */ 00323 if ( fabs (pdat->longitude) > 180.0 ) 00324 retval |= (1L << S_LON_ERROR); 00325 if ( fabs (pdat->latitude) > 90.0 ) 00326 retval |= (1L << S_LAT_ERROR); 00327 } 00328 00329 /* No silly temperatures or pressures, please. */ 00330 if ( (pdat->function & L_REFRAC) && (fabs (pdat->temp) > 100.0) ) 00331 retval |= (1L << S_TEMP_ERROR); 00332 if ( (pdat->function & L_REFRAC) && 00333 (pdat->press < 0.0) || (pdat->press > 2000.0) ) 00334 retval |= (1L << S_PRESS_ERROR); 00335 00336 /* No out of bounds tilts, please */ 00337 if ( (pdat->function & L_TILT) && (fabs (pdat->tilt) > 180.0) ) 00338 retval |= (1L << S_TILT_ERROR); 00339 if ( (pdat->function & L_TILT) && (fabs (pdat->aspect) > 360.0) ) 00340 retval |= (1L << S_ASPECT_ERROR); 00341 00342 /* No oddball shadowbands, please */ 00343 if ( (pdat->function & L_SBCF) && 00344 (pdat->sbwid < 1.0) || (pdat->sbwid > 100.0) ) 00345 retval |= (1L << S_SBWID_ERROR); 00346 if ( (pdat->function & L_SBCF) && 00347 (pdat->sbrad < 1.0) || (pdat->sbrad > 100.0) ) 00348 retval |= (1L << S_SBRAD_ERROR); 00349 if ( (pdat->function & L_SBCF) && ( fabs (pdat->sbsky) > 1.0) ) 00350 retval |= (1L << S_SBSKY_ERROR); 00351 00352 return retval; 00353 } 00354 00355 00356 /*============================================================================ 00357 * Local Void function dom2doy 00358 * 00359 * Converts day-of-month to day-of-year 00360 * 00361 * Requires (from struct posdata parameter): 00362 * year 00363 * month 00364 * day 00365 * 00366 * Returns (via the struct posdata parameter): 00367 * year 00368 * daynum 00369 *----------------------------------------------------------------------------*/ 00370 static void dom2doy( struct posdata *pdat ) 00371 { 00372 pdat->daynum = pdat->day + month_days[0][pdat->month]; 00373 00374 /* (adjust for leap year) */ 00375 if ( ((pdat->year % 4) == 0) && 00376 ( ((pdat->year % 100) != 0) || ((pdat->year % 400) == 0) ) && 00377 (pdat->month > 2) ) 00378 pdat->daynum += 1; 00379 } 00380 00381 00382 /*============================================================================ 00383 * Local void function doy2dom 00384 * 00385 * This function computes the month/day from the day number. 00386 * 00387 * Requires (from struct posdata parameter): 00388 * Year and day number: 00389 * year 00390 * daynum 00391 * 00392 * Returns (via the struct posdata parameter): 00393 * year 00394 * month 00395 * day 00396 *----------------------------------------------------------------------------*/ 00397 static void doy2dom(struct posdata *pdat) 00398 { 00399 int imon; /* Month (month_days) array counter */ 00400 int leap; /* leap year switch */ 00401 00402 /* Set the leap year switch */ 00403 if ( ((pdat->year % 4) == 0) && 00404 ( ((pdat->year % 100) != 0) || ((pdat->year % 400) == 0) ) ) 00405 leap = 1; 00406 else 00407 leap = 0; 00408 00409 /* Find the month */ 00410 imon = 12; 00411 while ( pdat->daynum <= month_days [leap][imon] ) 00412 --imon; 00413 00414 /* Set the month and day of month */ 00415 pdat->month = imon; 00416 pdat->day = pdat->daynum - month_days[leap][imon]; 00417 } 00418 00419 00420 /*============================================================================ 00421 * Local Void function geometry 00422 * 00423 * Does the underlying geometry for a given time and location 00424 *----------------------------------------------------------------------------*/ 00425 static void geometry ( struct posdata *pdat ) 00426 { 00427 float bottom; /* denominator (bottom) of the fraction */ 00428 float c2; /* cosine of d2 */ 00429 float cd; /* cosine of the day angle or delination */ 00430 float d2; /* pdat->dayang times two */ 00431 float delta; /* difference between current year and 1949 */ 00432 float s2; /* sine of d2 */ 00433 float sd; /* sine of the day angle */ 00434 float top; /* numerator (top) of the fraction */ 00435 int leap; /* leap year counter */ 00436 00437 /* Day angle */ 00438 /* Iqbal, M. 1983. An Introduction to Solar Radiation. 00439 Academic Press, NY., page 3 */ 00440 pdat->dayang = 360.0 * ( pdat->daynum - 1 ) / 365.0; 00441 00442 /* Earth radius vector * solar constant = solar energy */ 00443 /* Spencer, J. W. 1971. Fourier series representation of the 00444 position of the sun. Search 2 (5), page 172 */ 00445 sd = sin (raddeg * pdat->dayang); 00446 cd = cos (raddeg * pdat->dayang); 00447 d2 = 2.0 * pdat->dayang; 00448 c2 = cos (raddeg * d2); 00449 s2 = sin (raddeg * d2); 00450 00451 pdat->erv = 1.000110 + 0.034221 * cd + 0.001280 * sd; 00452 pdat->erv += 0.000719 * c2 + 0.000077 * s2; 00453 00454 /* Universal Coordinated (Greenwich standard) time */ 00455 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00456 approximate solar position (1950-2050). Solar Energy 40 (3), 00457 pp. 227-235. */ 00458 pdat->utime = 00459 pdat->hour * 3600.0 + 00460 pdat->minute * 60.0 + 00461 pdat->second - 00462 (float)pdat->interval / 2.0; 00463 pdat->utime = pdat->utime / 3600.0 - pdat->timezone; 00464 00465 /* Julian Day minus 2,400,000 days (to eliminate roundoff errors) */ 00466 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00467 approximate solar position (1950-2050). Solar Energy 40 (3), 00468 pp. 227-235. */ 00469 00470 /* No adjustment for century non-leap years since this function is 00471 bounded by 1950 - 2050 */ 00472 delta = pdat->year - 1949; 00473 leap = (int) ( delta / 4.0 ); 00474 pdat->julday = 00475 32916.5 + delta * 365.0 + leap + pdat->daynum + pdat->utime / 24.0; 00476 00477 /* Time used in the calculation of ecliptic coordinates */ 00478 /* Noon 1 JAN 2000 = 2,400,000 + 51,545 days Julian Date */ 00479 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00480 approximate solar position (1950-2050). Solar Energy 40 (3), 00481 pp. 227-235. */ 00482 pdat->ectime = pdat->julday - 51545.0; 00483 00484 /* Mean longitude */ 00485 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00486 approximate solar position (1950-2050). Solar Energy 40 (3), 00487 pp. 227-235. */ 00488 pdat->mnlong = 280.460 + 0.9856474 * pdat->ectime; 00489 00490 /* (dump the multiples of 360, so the answer is between 0 and 360) */ 00491 pdat->mnlong -= 360.0 * (int) ( pdat->mnlong / 360.0 ); 00492 if ( pdat->mnlong < 0.0 ) 00493 pdat->mnlong += 360.0; 00494 00495 /* Mean anomaly */ 00496 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00497 approximate solar position (1950-2050). Solar Energy 40 (3), 00498 pp. 227-235. */ 00499 pdat->mnanom = 357.528 + 0.9856003 * pdat->ectime; 00500 00501 /* (dump the multiples of 360, so the answer is between 0 and 360) */ 00502 pdat->mnanom -= 360.0 * (int) ( pdat->mnanom / 360.0 ); 00503 if ( pdat->mnanom < 0.0 ) 00504 pdat->mnanom += 360.0; 00505 00506 /* Ecliptic longitude */ 00507 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00508 approximate solar position (1950-2050). Solar Energy 40 (3), 00509 pp. 227-235. */ 00510 pdat->eclong = pdat->mnlong + 1.915 * sin ( pdat->mnanom * raddeg ) + 00511 0.020 * sin ( 2.0 * pdat->mnanom * raddeg ); 00512 00513 /* (dump the multiples of 360, so the answer is between 0 and 360) */ 00514 pdat->eclong -= 360.0 * (int) ( pdat->eclong / 360.0 ); 00515 if ( pdat->eclong < 0.0 ) 00516 pdat->eclong += 360.0; 00517 00518 /* Obliquity of the ecliptic */ 00519 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00520 approximate solar position (1950-2050). Solar Energy 40 (3), 00521 pp. 227-235. */ 00522 00523 /* 02 Feb 2001 SMW corrected sign in the following line */ 00524 /* pdat->ecobli = 23.439 + 4.0e-07 * pdat->ectime; */ 00525 pdat->ecobli = 23.439 - 4.0e-07 * pdat->ectime; 00526 00527 /* Declination */ 00528 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00529 approximate solar position (1950-2050). Solar Energy 40 (3), 00530 pp. 227-235. */ 00531 pdat->declin = degrad * asin ( sin (pdat->ecobli * raddeg) * 00532 sin (pdat->eclong * raddeg) ); 00533 00534 /* Right ascension */ 00535 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00536 approximate solar position (1950-2050). Solar Energy 40 (3), 00537 pp. 227-235. */ 00538 top = cos ( raddeg * pdat->ecobli ) * sin ( raddeg * pdat->eclong ); 00539 bottom = cos ( raddeg * pdat->eclong ); 00540 00541 pdat->rascen = degrad * atan2 ( top, bottom ); 00542 00543 /* (make it a positive angle) */ 00544 if ( pdat->rascen < 0.0 ) 00545 pdat->rascen += 360.0; 00546 00547 /* Greenwich mean sidereal time */ 00548 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00549 approximate solar position (1950-2050). Solar Energy 40 (3), 00550 pp. 227-235. */ 00551 pdat->gmst = 6.697375 + 0.0657098242 * pdat->ectime + pdat->utime; 00552 00553 /* (dump the multiples of 24, so the answer is between 0 and 24) */ 00554 pdat->gmst -= 24.0 * (int) ( pdat->gmst / 24.0 ); 00555 if ( pdat->gmst < 0.0 ) 00556 pdat->gmst += 24.0; 00557 00558 /* Local mean sidereal time */ 00559 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00560 approximate solar position (1950-2050). Solar Energy 40 (3), 00561 pp. 227-235. */ 00562 pdat->lmst = pdat->gmst * 15.0 + pdat->longitude; 00563 00564 /* (dump the multiples of 360, so the answer is between 0 and 360) */ 00565 pdat->lmst -= 360.0 * (int) ( pdat->lmst / 360.0 ); 00566 if ( pdat->lmst < 0.) 00567 pdat->lmst += 360.0; 00568 00569 /* Hour angle */ 00570 /* Michalsky, J. 1988. The Astronomical Almanac's algorithm for 00571 approximate solar position (1950-2050). Solar Energy 40 (3), 00572 pp. 227-235. */ 00573 pdat->hrang = pdat->lmst - pdat->rascen; 00574 00575 /* (force it between -180 and 180 degrees) */ 00576 if ( pdat->hrang < -180.0 ) 00577 pdat->hrang += 360.0; 00578 else if ( pdat->hrang > 180.0 ) 00579 pdat->hrang -= 360.0; 00580 } 00581 00582 00583 /*============================================================================ 00584 * Local Void function zen_no_ref 00585 * 00586 * ETR solar zenith angle 00587 * Iqbal, M. 1983. An Introduction to Solar Radiation. 00588 * Academic Press, NY., page 15 00589 *----------------------------------------------------------------------------*/ 00590 static void zen_no_ref ( struct posdata *pdat, struct trigdata *tdat ) 00591 { 00592 float cz; /* cosine of the solar zenith angle */ 00593 00594 localtrig( pdat, tdat ); 00595 cz = tdat->sd * tdat->sl + tdat->cd * tdat->cl * tdat->ch; 00596 00597 /* (watch out for the roundoff errors) */ 00598 if ( fabs (cz) > 1.0 ) { 00599 if ( cz >= 0.0 ) 00600 cz = 1.0; 00601 else 00602 cz = -1.0; 00603 } 00604 00605 pdat->zenetr = acos ( cz ) * degrad; 00606 00607 /* (limit the degrees below the horizon to 9 [+90 -> 99]) */ 00608 if ( pdat->zenetr > 99.0 ) 00609 pdat->zenetr = 99.0; 00610 00611 pdat->elevetr = 90.0 - pdat->zenetr; 00612 } 00613 00614 00615 /*============================================================================ 00616 * Local Void function ssha 00617 * 00618 * Sunset hour angle, degrees 00619 * Iqbal, M. 1983. An Introduction to Solar Radiation. 00620 * Academic Press, NY., page 16 00621 *----------------------------------------------------------------------------*/ 00622 static void ssha( struct posdata *pdat, struct trigdata *tdat ) 00623 { 00624 float cssha; /* cosine of the sunset hour angle */ 00625 float cdcl; /* ( cd * cl ) */ 00626 00627 localtrig( pdat, tdat ); 00628 cdcl = tdat->cd * tdat->cl; 00629 00630 if ( fabs ( cdcl ) >= 0.001 ) { 00631 cssha = -tdat->sl * tdat->sd / cdcl; 00632 00633 /* This keeps the cosine from blowing on roundoff */ 00634 if ( cssha < -1.0 ) 00635 pdat->ssha = 180.0; 00636 else if ( cssha > 1.0 ) 00637 pdat->ssha = 0.0; 00638 else 00639 pdat->ssha = degrad * acos ( cssha ); 00640 } 00641 else if ( ((pdat->declin >= 0.0) && (pdat->latitude > 0.0 )) || 00642 ((pdat->declin < 0.0) && (pdat->latitude < 0.0 )) ) 00643 pdat->ssha = 180.0; 00644 else 00645 pdat->ssha = 0.0; 00646 } 00647 00648 00649 /*============================================================================ 00650 * Local Void function sbcf 00651 * 00652 * Shadowband correction factor 00653 * Drummond, A. J. 1956. A contribution to absolute pyrheliometry. 00654 * Q. J. R. Meteorol. Soc. 82, pp. 481-493 00655 *----------------------------------------------------------------------------*/ 00656 static void sbcf( struct posdata *pdat, struct trigdata *tdat ) 00657 { 00658 float p, t1, t2; /* used to compute sbcf */ 00659 00660 localtrig( pdat, tdat ); 00661 p = 0.6366198 * pdat->sbwid / pdat->sbrad * pow (tdat->cd,3); 00662 t1 = tdat->sl * tdat->sd * pdat->ssha * raddeg; 00663 t2 = tdat->cl * tdat->cd * sin ( pdat->ssha * raddeg ); 00664 pdat->sbcf = pdat->sbsky + 1.0 / ( 1.0 - p * ( t1 + t2 ) ); 00665 00666 } 00667 00668 00669 /*============================================================================ 00670 * Local Void function tst 00671 * 00672 * TST -> True Solar Time = local standard time + TSTfix, time 00673 * in minutes from midnight. 00674 * Iqbal, M. 1983. An Introduction to Solar Radiation. 00675 * Academic Press, NY., page 13 00676 *----------------------------------------------------------------------------*/ 00677 static void tst( struct posdata *pdat ) 00678 { 00679 pdat->tst = ( 180.0 + pdat->hrang ) * 4.0; 00680 pdat->tstfix = 00681 pdat->tst - 00682 (float)pdat->hour * 60.0 - 00683 pdat->minute - 00684 (float)pdat->second / 60.0 + 00685 (float)pdat->interval / 120.0; /* add back half of the interval */ 00686 00687 /* bound tstfix to this day */ 00688 while ( pdat->tstfix > 720.0 ) 00689 pdat->tstfix -= 1440.0; 00690 while ( pdat->tstfix < -720.0 ) 00691 pdat->tstfix += 1440.0; 00692 00693 pdat->eqntim = 00694 pdat->tstfix + 60.0 * pdat->timezone - 4.0 * pdat->longitude; 00695 00696 } 00697 00698 00699 /*============================================================================ 00700 * Local Void function srss 00701 * 00702 * Sunrise and sunset times (minutes from midnight) 00703 *----------------------------------------------------------------------------*/ 00704 static void srss( struct posdata *pdat ) 00705 { 00706 if ( pdat->ssha <= 1.0 ) { 00707 pdat->sretr = 2999.0; 00708 pdat->ssetr = -2999.0; 00709 } 00710 else if ( pdat->ssha >= 179.0 ) { 00711 pdat->sretr = -2999.0; 00712 pdat->ssetr = 2999.0; 00713 } 00714 else { 00715 pdat->sretr = 720.0 - 4.0 * pdat->ssha - pdat->tstfix; 00716 pdat->ssetr = 720.0 + 4.0 * pdat->ssha - pdat->tstfix; 00717 } 00718 } 00719 00720 00721 /*============================================================================ 00722 * Local Void function sazm 00723 * 00724 * Solar azimuth angle 00725 * Iqbal, M. 1983. An Introduction to Solar Radiation. 00726 * Academic Press, NY., page 15 00727 *----------------------------------------------------------------------------*/ 00728 static void sazm( struct posdata *pdat, struct trigdata *tdat ) 00729 { 00730 float ca; /* cosine of the solar azimuth angle */ 00731 float ce; /* cosine of the solar elevation */ 00732 float cecl; /* ( ce * cl ) */ 00733 float se; /* sine of the solar elevation */ 00734 00735 localtrig( pdat, tdat ); 00736 ce = cos ( raddeg * pdat->elevetr ); 00737 se = sin ( raddeg * pdat->elevetr ); 00738 00739 pdat->azim = 180.0; 00740 cecl = ce * tdat->cl; 00741 if ( fabs ( cecl ) >= 0.001 ) { 00742 ca = ( se * tdat->sl - tdat->sd ) / cecl; 00743 if ( ca > 1.0 ) 00744 ca = 1.0; 00745 else if ( ca < -1.0 ) 00746 ca = -1.0; 00747 00748 pdat->azim = 180.0 - acos ( ca ) * degrad; 00749 if ( pdat->hrang > 0 ) 00750 pdat->azim = 360.0 - pdat->azim; 00751 } 00752 } 00753 00754 00755 /*============================================================================ 00756 * Local Int function refrac 00757 * 00758 * Refraction correction, degrees 00759 * Zimmerman, John C. 1981. Sun-pointing programs and their 00760 * accuracy. 00761 * SAND81-0761, Experimental Systems Operation Division 4721, 00762 * Sandia National Laboratories, Albuquerque, NM. 00763 *----------------------------------------------------------------------------*/ 00764 static void refrac( struct posdata *pdat ) 00765 { 00766 float prestemp; /* temporary pressure/temperature correction */ 00767 float refcor; /* temporary refraction correction */ 00768 float tanelev; /* tangent of the solar elevation angle */ 00769 00770 /* If the sun is near zenith, the algorithm bombs; refraction near 0 */ 00771 if ( pdat->elevetr > 85.0 ) 00772 refcor = 0.0; 00773 00774 /* Otherwise, we have refraction */ 00775 else { 00776 tanelev = tan ( raddeg * pdat->elevetr ); 00777 if ( pdat->elevetr >= 5.0 ) 00778 refcor = 58.1 / tanelev - 00779 0.07 / ( pow (tanelev,3) ) + 00780 0.000086 / ( pow (tanelev,5) ); 00781 else if ( pdat->elevetr >= -0.575 ) 00782 refcor = 1735.0 + 00783 pdat->elevetr * ( -518.2 + pdat->elevetr * ( 103.4 + 00784 pdat->elevetr * ( -12.79 + pdat->elevetr * 0.711 ) ) ); 00785 else 00786 refcor = -20.774 / tanelev; 00787 00788 prestemp = 00789 ( pdat->press * 283.0 ) / ( 1013.0 * ( 273.0 + pdat->temp ) ); 00790 refcor *= prestemp / 3600.0; 00791 } 00792 00793 /* Refracted solar elevation angle */ 00794 pdat->elevref = pdat->elevetr + refcor; 00795 00796 /* (limit the degrees below the horizon to 9) */ 00797 if ( pdat->elevref < -9.0 ) 00798 pdat->elevref = -9.0; 00799 00800 /* Refracted solar zenith angle */ 00801 pdat->zenref = 90.0 - pdat->elevref; 00802 pdat->coszen = cos( raddeg * pdat->zenref ); 00803 } 00804 00805 00806 /*============================================================================ 00807 * Local Void function amass 00808 * 00809 * Airmass 00810 * Kasten, F. and Young, A. 1989. Revised optical air mass 00811 * tables and approximation formula. Applied Optics 28 (22), 00812 * pp. 4735-4738 00813 *----------------------------------------------------------------------------*/ 00814 static void amass( struct posdata *pdat ) 00815 { 00816 if ( pdat->zenref > 93.0 ) 00817 { 00818 pdat->amass = -1.0; 00819 pdat->ampress = -1.0; 00820 } 00821 else 00822 { 00823 pdat->amass = 00824 1.0 / ( cos (raddeg * pdat->zenref) + 0.50572 * 00825 pow ((96.07995 - pdat->zenref),-1.6364) ); 00826 00827 pdat->ampress = pdat->amass * pdat->press / 1013.0; 00828 } 00829 } 00830 00831 00832 /*============================================================================ 00833 * Local Void function prime 00834 * 00835 * Prime and Unprime 00836 * Prime converts Kt to normalized Kt', etc. 00837 * Unprime deconverts Kt' to Kt, etc. 00838 * Perez, R., P. Ineichen, Seals, R., & Zelenka, A. 1990. Making 00839 * full use of the clearness index for parameterizing hourly 00840 * insolation conditions. Solar Energy 45 (2), pp. 111-114 00841 *----------------------------------------------------------------------------*/ 00842 static void prime( struct posdata *pdat ) 00843 { 00844 pdat->unprime = 1.031 * exp ( -1.4 / ( 0.9 + 9.4 / pdat->amass ) ) + 0.1; 00845 pdat->prime = 1.0 / pdat->unprime; 00846 } 00847 00848 00849 /*============================================================================ 00850 * Local Void function etr 00851 * 00852 * Extraterrestrial (top-of-atmosphere) solar irradiance 00853 *----------------------------------------------------------------------------*/ 00854 static void etr( struct posdata *pdat ) 00855 { 00856 if ( pdat->coszen > 0.0 ) { 00857 pdat->etrn = pdat->solcon * pdat->erv; 00858 pdat->etr = pdat->etrn * pdat->coszen; 00859 } 00860 else { 00861 pdat->etrn = 0.0; 00862 pdat->etr = 0.0; 00863 } 00864 } 00865 00866 00867 /*============================================================================ 00868 * Local Void function localtrig 00869 * 00870 * Does trig on internal variable used by several functions 00871 *----------------------------------------------------------------------------*/ 00872 static void localtrig( struct posdata *pdat, struct trigdata *tdat ) 00873 { 00874 /* define masks to prevent calculation of uninitialized variables */ 00875 #define SD_MASK ( L_ZENETR | L_SSHA | S_SBCF | S_SOLAZM ) 00876 #define SL_MASK ( L_ZENETR | L_SSHA | S_SBCF | S_SOLAZM ) 00877 #define CL_MASK ( L_ZENETR | L_SSHA | S_SBCF | S_SOLAZM ) 00878 #define CD_MASK ( L_ZENETR | L_SSHA | S_SBCF ) 00879 #define CH_MASK ( L_ZENETR ) 00880 00881 if ( tdat->sd < -900.0 ) /* sd was initialized -999 as flag */ 00882 { 00883 tdat->sd = 1.0; /* reflag as having completed calculations */ 00884 if ( pdat->function | CD_MASK ) 00885 tdat->cd = cos ( raddeg * pdat->declin ); 00886 if ( pdat->function | CH_MASK ) 00887 tdat->ch = cos ( raddeg * pdat->hrang ); 00888 if ( pdat->function | CL_MASK ) 00889 tdat->cl = cos ( raddeg * pdat->latitude ); 00890 if ( pdat->function | SD_MASK ) 00891 tdat->sd = sin ( raddeg * pdat->declin ); 00892 if ( pdat->function | SL_MASK ) 00893 tdat->sl = sin ( raddeg * pdat->latitude ); 00894 } 00895 } 00896 00897 00898 /*============================================================================ 00899 * Local Void function tilt 00900 * 00901 * ETR on a tilted surface 00902 *----------------------------------------------------------------------------*/ 00903 static void tilt( struct posdata *pdat ) 00904 { 00905 float ca; /* cosine of the solar azimuth angle */ 00906 float cp; /* cosine of the panel aspect */ 00907 float ct; /* cosine of the panel tilt */ 00908 float sa; /* sine of the solar azimuth angle */ 00909 float sp; /* sine of the panel aspect */ 00910 float st; /* sine of the panel tilt */ 00911 float sz; /* sine of the refraction corrected solar zenith angle */ 00912 00913 00914 /* Cosine of the angle between the sun and a tipped flat surface, 00915 useful for calculating solar energy on tilted surfaces */ 00916 ca = cos ( raddeg * pdat->azim ); 00917 cp = cos ( raddeg * pdat->aspect ); 00918 ct = cos ( raddeg * pdat->tilt ); 00919 sa = sin ( raddeg * pdat->azim ); 00920 sp = sin ( raddeg * pdat->aspect ); 00921 st = sin ( raddeg * pdat->tilt ); 00922 sz = sin ( raddeg * pdat->zenref ); 00923 pdat->cosinc = pdat->coszen * ct + sz * st * ( ca * cp + sa * sp ); 00924 00925 if ( pdat->cosinc > 0.0 ) 00926 pdat->etrtilt = pdat->etrn * pdat->cosinc; 00927 else 00928 pdat->etrtilt = 0.0; 00929 00930 } 00931 00932 00933 /*============================================================================ 00934 * Void function S_decode 00935 * 00936 * This function decodes the error codes from S_solpos return value 00937 * 00938 * Requires the long integer return value from S_solpos 00939 * 00940 * Returns descriptive text to stderr 00941 *----------------------------------------------------------------------------*/ 00942 void S_decode(long code, struct posdata *pdat) 00943 { 00944 if ( code & (1L << S_YEAR_ERROR) ) 00945 fprintf(stderr, "S_decode ==> Please fix the year: %d [1950-2050]\n", 00946 pdat->year); 00947 if ( code & (1L << S_MONTH_ERROR) ) 00948 fprintf(stderr, "S_decode ==> Please fix the month: %d\n", 00949 pdat->month); 00950 if ( code & (1L << S_DAY_ERROR) ) 00951 fprintf(stderr, "S_decode ==> Please fix the day-of-month: %d\n", 00952 pdat->day); 00953 if ( code & (1L << S_DOY_ERROR) ) 00954 fprintf(stderr, "S_decode ==> Please fix the day-of-year: %d\n", 00955 pdat->daynum); 00956 if ( code & (1L << S_HOUR_ERROR) ) 00957 fprintf(stderr, "S_decode ==> Please fix the hour: %d\n", 00958 pdat->hour); 00959 if ( code & (1L << S_MINUTE_ERROR) ) 00960 fprintf(stderr, "S_decode ==> Please fix the minute: %d\n", 00961 pdat->minute); 00962 if ( code & (1L << S_SECOND_ERROR) ) 00963 fprintf(stderr, "S_decode ==> Please fix the second: %d\n", 00964 pdat->second); 00965 if ( code & (1L << S_TZONE_ERROR) ) 00966 fprintf(stderr, "S_decode ==> Please fix the time zone: %f\n", 00967 pdat->timezone); 00968 if ( code & (1L << S_INTRVL_ERROR) ) 00969 fprintf(stderr, "S_decode ==> Please fix the interval: %d\n", 00970 pdat->interval); 00971 if ( code & (1L << S_LAT_ERROR) ) 00972 fprintf(stderr, "S_decode ==> Please fix the latitude: %f\n", 00973 pdat->latitude); 00974 if ( code & (1L << S_LON_ERROR) ) 00975 fprintf(stderr, "S_decode ==> Please fix the longitude: %f\n", 00976 pdat->longitude); 00977 if ( code & (1L << S_TEMP_ERROR) ) 00978 fprintf(stderr, "S_decode ==> Please fix the temperature: %f\n", 00979 pdat->temp); 00980 if ( code & (1L << S_PRESS_ERROR) ) 00981 fprintf(stderr, "S_decode ==> Please fix the pressure: %f\n", 00982 pdat->press); 00983 if ( code & (1L << S_TILT_ERROR) ) 00984 fprintf(stderr, "S_decode ==> Please fix the tilt: %f\n", 00985 pdat->tilt); 00986 if ( code & (1L << S_ASPECT_ERROR) ) 00987 fprintf(stderr, "S_decode ==> Please fix the aspect: %f\n", 00988 pdat->aspect); 00989 if ( code & (1L << S_SBWID_ERROR) ) 00990 fprintf(stderr, "S_decode ==> Please fix the shadowband width: %f\n", 00991 pdat->sbwid); 00992 if ( code & (1L << S_SBRAD_ERROR) ) 00993 fprintf(stderr, "S_decode ==> Please fix the shadowband radius: %f\n", 00994 pdat->sbrad); 00995 if ( code & (1L << S_SBSKY_ERROR) ) 00996 fprintf(stderr, "S_decode ==> Please fix the shadowband sky factor: %f\n", 00997 pdat->sbsky); 00998 } 00999
Generated on Tue Jul 12 2022 17:04:07 by
1.7.2