James Cobb
/
solpos
Important changes to repositories hosted on mbed.com
Mbed hosted mercurial repositories are deprecated and are due to be permanently deleted in July 2026.
To keep a copy of this software download the repository Zip archive or clone locally using Mercurial.
It is also possible to export all your personal repositories from the account settings page.
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