| 1 |
/* Copyright (c) 1994,2006 *Fraunhofer Institut for Solar Energy Systems |
| 2 |
* Heidenhofstr. 2, D-79110 Freiburg, Germany |
| 3 |
* *Agence de l'Environnement et de la Maitrise de l'Energie |
| 4 |
* Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France |
| 5 |
* *BOUYGUES |
| 6 |
* 1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France |
| 7 |
* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 |
| 8 |
* version 2.6 (2021/01/29): dew point dependency added according to Perez publication 1990 (W -> atm_preci_water=exp(0.07*Td-0.075) ). by J. Wienold, EPFL |
| 9 |
*/ |
| 10 |
|
| 11 |
#define _USE_MATH_DEFINES |
| 12 |
#include <stdio.h> |
| 13 |
#include <string.h> |
| 14 |
#include <math.h> |
| 15 |
#include <stdlib.h> |
| 16 |
|
| 17 |
#include "color.h" |
| 18 |
#include "sun.h" |
| 19 |
#include "paths.h" |
| 20 |
|
| 21 |
#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) |
| 22 |
|
| 23 |
double normsc(); |
| 24 |
|
| 25 |
/*static char *rcsid="$Header: /usr/local/cvs/radiance/ray/src/gen/gendaylit.c,v 2.20 2020/11/18 17:18:41 greg Exp $";*/ |
| 26 |
|
| 27 |
float coeff_perez[] = { |
| 28 |
1.3525,-0.2576,-0.2690,-1.4366,-0.7670,0.0007,1.2734,-0.1233,2.8000,0.6004,1.2375,1.000,1.8734,0.6297, |
| 29 |
0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,-1.2219,-0.7730,1.4148,1.1016,-0.2054,0.0367,-3.9128,0.9156, |
| 30 |
6.9750,0.1774,6.4477,-0.1239,-1.5798,-0.5081,-1.7812,0.1080,0.2624,0.0672,-0.2190,-0.4285,-1.1000,-0.2515, |
| 31 |
0.8952,0.0156,0.2782,-0.1812,-4.5000,1.1766,24.7219,-13.0812,-37.7000,34.8438,-5.0000,1.5218,3.9229, |
| 32 |
-2.6204,-0.0156,0.1597,0.4199,-0.5562,-0.5484,-0.6654,-0.2672,0.7117,0.7234,-0.6219,-5.6812,2.6297, |
| 33 |
33.3389,-18.3000,-62.2500,52.0781,-3.5000,0.0016,1.1477,0.1062,0.4659,-0.3296,-0.0876,-0.0329,-0.6000, |
| 34 |
-0.3566,-2.5000,2.3250,0.2937,0.0496,-5.6812,1.8415,21.0000,-4.7656,-21.5906,7.2492,-3.5000,-0.1554, |
| 35 |
1.4062,0.3988,0.0032,0.0766,-0.0656,-0.1294,-1.0156,-0.3670,1.0078,1.4051,0.2875,-0.5328,-3.8500,3.3750, |
| 36 |
14.0000,-0.9999,-7.1406,7.5469,-3.4000,-0.1078,-1.0750,1.5702,-0.0672,0.4016,0.3017,-0.4844,-1.0000, |
| 37 |
0.0211,0.5025,-0.5119,-0.3000,0.1922,0.7023,-1.6317,19.0000,-5.0000,1.2438,-1.9094,-4.0000,0.0250,0.3844, |
| 38 |
0.2656,1.0468,-0.3788,-2.4517,1.4656,-1.0500,0.0289,0.4260,0.3590,-0.3250,0.1156,0.7781,0.0025,31.0625, |
| 39 |
-14.5000,-46.1148,55.3750,-7.2312,0.4050,13.3500,0.6234,1.5000,-0.6426,1.8564,0.5636}; |
| 40 |
|
| 41 |
|
| 42 |
float defangle_theta[] = { |
| 43 |
84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, |
| 44 |
84, 84, 84, 84, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, |
| 45 |
72, 72, 72, 72, 72, 72, 72, 72, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, |
| 46 |
60, 60, 60, 60, 60, 60, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, |
| 47 |
48, 48, 48, 48, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 24, 24, 24, 24, |
| 48 |
24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0}; |
| 49 |
|
| 50 |
float defangle_phi[] = { |
| 51 |
0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, |
| 52 |
276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, |
| 53 |
192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105, |
| 54 |
120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75, |
| 55 |
90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60, |
| 56 |
80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210, |
| 57 |
240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0}; |
| 58 |
/* default values for Berlin */ |
| 59 |
float locus[] = { |
| 60 |
-4.843e9,2.5568e6,0.24282e3,0.23258,-4.843e9,2.5568e6,0.24282e3,0.23258,-1.2848,1.7519,-0.093786}; |
| 61 |
|
| 62 |
|
| 63 |
|
| 64 |
/* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */ |
| 65 |
double sky_brightness(); |
| 66 |
double sky_clearness(); |
| 67 |
|
| 68 |
/* calculation of the direct and diffuse components from the Perez parametrization */ |
| 69 |
double diffuse_irradiance_from_sky_brightness(); |
| 70 |
double direct_irradiance_from_sky_clearness(); |
| 71 |
|
| 72 |
/* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : */ |
| 73 |
/* input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */ |
| 74 |
|
| 75 |
double glob_h_effi_PEREZ(); |
| 76 |
double glob_h_diffuse_effi_PEREZ(); |
| 77 |
double direct_n_effi_PEREZ(); |
| 78 |
|
| 79 |
/*likelihood check of the epsilon, delta, direct and diffuse components*/ |
| 80 |
void check_parametrization(); |
| 81 |
void check_irradiances(); |
| 82 |
void check_illuminances(); |
| 83 |
void illu_to_irra_index(); |
| 84 |
void print_error_sky(); |
| 85 |
|
| 86 |
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]); |
| 87 |
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]); |
| 88 |
double radians(double degres); |
| 89 |
double degres(double radians); |
| 90 |
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z); |
| 91 |
double integ_lv(float *lv,float *theta); |
| 92 |
|
| 93 |
void printdefaults(); |
| 94 |
void check_sun_position(); |
| 95 |
void computesky(); |
| 96 |
void printhead(int ac, char** av); |
| 97 |
void usage_error(char* msg); |
| 98 |
void printsky(); |
| 99 |
|
| 100 |
FILE * frlibopen(char* fname); |
| 101 |
|
| 102 |
/* astronomy and geometry*/ |
| 103 |
double get_eccentricity(); |
| 104 |
double air_mass(); |
| 105 |
|
| 106 |
double solar_sunset(int month, int day); |
| 107 |
double solar_sunrise(int month, int day); |
| 108 |
|
| 109 |
const double AU = 149597890E3; |
| 110 |
const double solar_constant_e = 1367; /* solar constant W/m^2 */ |
| 111 |
const double solar_constant_l = 127500; /* solar constant lux */ |
| 112 |
|
| 113 |
const double half_sun_angle = 0.2665; |
| 114 |
const double half_direct_angle = 2.85; |
| 115 |
|
| 116 |
const double skyclearinf = 1.0; /* limitations for the variation of the Perez parameters */ |
| 117 |
const double skyclearsup = 12.01; |
| 118 |
const double skybriginf = 0.01; |
| 119 |
const double skybrigsup = 0.6; |
| 120 |
|
| 121 |
|
| 122 |
|
| 123 |
/* required values */ |
| 124 |
int year = 0; /* year (optional) */ |
| 125 |
int month, day; /* date */ |
| 126 |
double hour; /* time */ |
| 127 |
int tsolar; /* 0=standard, 1=solar */ |
| 128 |
double altitude, azimuth; /* or solar angles */ |
| 129 |
|
| 130 |
|
| 131 |
|
| 132 |
/* definition of the sky conditions through the Perez parametrization */ |
| 133 |
double skyclearness = 0; |
| 134 |
double skybrightness = 0; |
| 135 |
double solarradiance; |
| 136 |
double diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance; |
| 137 |
double sunzenith, daynumber, atm_preci_water, Td=10.97353115; |
| 138 |
|
| 139 |
/*double sunaltitude_border = 0;*/ |
| 140 |
double diffnormalization = 0; |
| 141 |
double dirnormalization = 0; |
| 142 |
double *c_perez; |
| 143 |
|
| 144 |
int output=0; /* define the unit of the output (sky luminance or radiance): */ |
| 145 |
/* visible watt=0, solar watt=1, lumen=2 */ |
| 146 |
int input=0; /* define the input for the calulation */ |
| 147 |
int color_output=0; |
| 148 |
int suppress_warnings=0; |
| 149 |
|
| 150 |
/* default values */ |
| 151 |
int cloudy = 0; /* 1=standard, 2=uniform */ |
| 152 |
int dosun = 1; |
| 153 |
double zenithbr = -1.0; |
| 154 |
double betaturbidity = 0.1; |
| 155 |
double gprefl = 0.2; |
| 156 |
int S_INTER=0; |
| 157 |
|
| 158 |
|
| 159 |
/* computed values */ |
| 160 |
double sundir[3]; |
| 161 |
double groundbr = 0; |
| 162 |
double F2; |
| 163 |
double solarbr = 0.0; |
| 164 |
int u_solar = 0; /* -1=irradiance, 1=radiance */ |
| 165 |
float timeinterval = 0; |
| 166 |
|
| 167 |
char *progname; |
| 168 |
char errmsg[128]; |
| 169 |
|
| 170 |
double st; |
| 171 |
|
| 172 |
|
| 173 |
int main(int argc, char** argv) |
| 174 |
{ |
| 175 |
int i; |
| 176 |
|
| 177 |
progname = argv[0]; |
| 178 |
if (argc == 2 && !strcmp(argv[1], "-defaults")) { |
| 179 |
printdefaults(); |
| 180 |
return 0; |
| 181 |
} |
| 182 |
if (argc < 4) |
| 183 |
usage_error("arg count"); |
| 184 |
if (!strcmp(argv[1], "-ang")) { |
| 185 |
altitude = atof(argv[2]) * (M_PI/180); |
| 186 |
azimuth = atof(argv[3]) * (M_PI/180); |
| 187 |
month = 0; |
| 188 |
} else { |
| 189 |
month = atoi(argv[1]); |
| 190 |
if (month < 1 || month > 12) |
| 191 |
usage_error("bad month"); |
| 192 |
day = atoi(argv[2]); |
| 193 |
if (day < 1 || day > 31) |
| 194 |
usage_error("bad day"); |
| 195 |
hour = atof(argv[3]); |
| 196 |
if (hour < 0 || hour >= 24) |
| 197 |
usage_error("bad hour"); |
| 198 |
tsolar = argv[3][0] == '+'; |
| 199 |
} |
| 200 |
for (i = 4; i < argc; i++) |
| 201 |
if (argv[i][0] == '-' || argv[i][0] == '+') |
| 202 |
switch (argv[i][1]) { |
| 203 |
case 'd': |
| 204 |
Td = atof(argv[++i]); |
| 205 |
if (Td < -40 || Td > 40) { |
| 206 |
Td=10.97353115; } |
| 207 |
break; |
| 208 |
case 's': |
| 209 |
cloudy = 0; |
| 210 |
dosun = argv[i][0] == '+'; |
| 211 |
break; |
| 212 |
case 'y': |
| 213 |
year = atoi(argv[++i]); |
| 214 |
break; |
| 215 |
case 'R': |
| 216 |
u_solar = argv[i][1] == 'R' ? -1 : 1; |
| 217 |
solarbr = atof(argv[++i]); |
| 218 |
break; |
| 219 |
case 'c': |
| 220 |
cloudy = argv[i][0] == '+' ? 2 : 1; |
| 221 |
dosun = 0; |
| 222 |
break; |
| 223 |
case 'C': |
| 224 |
if (argv[i][2] == 'I' && argv[i][3] == 'E' ) { |
| 225 |
locus[0] = -4.607e9; |
| 226 |
locus[1] = 2.9678e6; |
| 227 |
locus[2] = 0.09911e3; |
| 228 |
locus[3] = 0.244063; |
| 229 |
locus[4] = -2.0064e9; |
| 230 |
locus[5] = 1.9018e6; |
| 231 |
locus[6] = 0.24748e3; |
| 232 |
locus[7] = 0.23704; |
| 233 |
locus[8] = -3.0; |
| 234 |
locus[9] = 2.87; |
| 235 |
locus[10] = -0.275; |
| 236 |
}else{ color_output = 1; |
| 237 |
} |
| 238 |
break; |
| 239 |
case 'l': |
| 240 |
locus[0] = atof(argv[++i]); |
| 241 |
locus[1] = atof(argv[++i]); |
| 242 |
locus[2] = atof(argv[++i]); |
| 243 |
locus[3] = atof(argv[++i]); |
| 244 |
locus[4] = locus[0]; |
| 245 |
locus[5] = locus[1]; |
| 246 |
locus[6] = locus[2]; |
| 247 |
locus[7] = locus[3]; |
| 248 |
locus[8] = atof(argv[++i]); |
| 249 |
locus[9] = atof(argv[++i]); |
| 250 |
locus[10] = atof(argv[++i]); |
| 251 |
break; |
| 252 |
|
| 253 |
case 't': |
| 254 |
betaturbidity = atof(argv[++i]); |
| 255 |
break; |
| 256 |
case 'w': |
| 257 |
suppress_warnings = 1; |
| 258 |
break; |
| 259 |
case 'b': |
| 260 |
zenithbr = atof(argv[++i]); |
| 261 |
break; |
| 262 |
case 'g': |
| 263 |
gprefl = atof(argv[++i]); |
| 264 |
break; |
| 265 |
case 'a': |
| 266 |
s_latitude = atof(argv[++i]) * (M_PI/180); |
| 267 |
break; |
| 268 |
case 'o': |
| 269 |
s_longitude = atof(argv[++i]) * (M_PI/180); |
| 270 |
break; |
| 271 |
case 'm': |
| 272 |
s_meridian = atof(argv[++i]) * (M_PI/180); |
| 273 |
break; |
| 274 |
|
| 275 |
case 'O': |
| 276 |
output = atof(argv[++i]); /*define the unit of the output of the program: |
| 277 |
sky and sun luminance/radiance |
| 278 |
(0==W visible, 1==W solar radiation, 2==lm) */ |
| 279 |
break; |
| 280 |
|
| 281 |
case 'P': |
| 282 |
input = 0; /* Perez parameters: epsilon, delta */ |
| 283 |
skyclearness = atof(argv[++i]); |
| 284 |
skybrightness = atof(argv[++i]); |
| 285 |
break; |
| 286 |
|
| 287 |
case 'W': /* direct normal Irradiance [W/m^2] */ |
| 288 |
input = 1; /* diffuse horizontal Irrad. [W/m^2] */ |
| 289 |
directirradiance = atof(argv[++i]); |
| 290 |
diffuseirradiance = atof(argv[++i]); |
| 291 |
break; |
| 292 |
|
| 293 |
case 'L': /* direct normal Illuminance [Lux] */ |
| 294 |
input = 2; /* diffuse horizontal Ill. [Lux] */ |
| 295 |
directilluminance = atof(argv[++i]); |
| 296 |
diffuseilluminance = atof(argv[++i]); |
| 297 |
break; |
| 298 |
|
| 299 |
case 'G': /* direct horizontal Irradiance [W/m^2] */ |
| 300 |
input = 3; /* diffuse horizontal Irrad. [W/m^2] */ |
| 301 |
directirradiance = atof(argv[++i]); |
| 302 |
diffuseirradiance = atof(argv[++i]); |
| 303 |
break; |
| 304 |
|
| 305 |
case 'E': /* Erbs model based on the */ |
| 306 |
input = 4; /* global-horizontal irradiance [W/m^2] */ |
| 307 |
globalirradiance = atof(argv[++i]); |
| 308 |
break; |
| 309 |
|
| 310 |
case 'i': |
| 311 |
timeinterval = atof(argv[++i]); |
| 312 |
break; |
| 313 |
|
| 314 |
|
| 315 |
default: |
| 316 |
sprintf(errmsg, "unknown option: %s", argv[i]); |
| 317 |
usage_error(errmsg); |
| 318 |
} |
| 319 |
else |
| 320 |
usage_error("bad option"); |
| 321 |
|
| 322 |
if (month && !tsolar && fabs(s_meridian-s_longitude) > 45*M_PI/180) |
| 323 |
fprintf(stderr,"%s: warning: %.1f hours btwn. standard meridian and longitude\n", |
| 324 |
progname, (s_longitude-s_meridian)*12/M_PI); |
| 325 |
|
| 326 |
|
| 327 |
/* dynamic memory allocation for the pointers */ |
| 328 |
if ( (c_perez = calloc(5, sizeof(double))) == NULL ) |
| 329 |
{ fprintf(stderr,"Out of memory error in function main"); return 1; } |
| 330 |
|
| 331 |
|
| 332 |
atm_preci_water=exp(0.07*Td-0.075); |
| 333 |
printhead(argc, argv); |
| 334 |
computesky(); |
| 335 |
printsky(); |
| 336 |
return 0; |
| 337 |
|
| 338 |
} |
| 339 |
|
| 340 |
|
| 341 |
|
| 342 |
|
| 343 |
|
| 344 |
void computesky() |
| 345 |
{ |
| 346 |
|
| 347 |
int j; |
| 348 |
|
| 349 |
float *lv_mod; /* 145 luminance values */ |
| 350 |
float *theta_o, *phi_o; |
| 351 |
double dzeta, gamma; |
| 352 |
double normfactor; |
| 353 |
double erbs_s0, erbs_kt; |
| 354 |
|
| 355 |
|
| 356 |
/* compute solar direction */ |
| 357 |
|
| 358 |
if (month) { /* from date and time */ |
| 359 |
double sd; |
| 360 |
|
| 361 |
st = hour; |
| 362 |
if (year) { /* Michalsky algorithm? */ |
| 363 |
double mjd = mjdate(year, month, day, hour); |
| 364 |
if (tsolar) |
| 365 |
sd = msdec(mjd, NULL); |
| 366 |
else |
| 367 |
sd = msdec(mjd, &st); |
| 368 |
} else { |
| 369 |
int jd = jdate(month, day); /* Julian date */ |
| 370 |
sd = sdec(jd); /* solar declination */ |
| 371 |
if (!tsolar) /* get solar time? */ |
| 372 |
st = hour + stadj(jd); |
| 373 |
} |
| 374 |
|
| 375 |
if(timeinterval) { |
| 376 |
|
| 377 |
if(timeinterval<0) { |
| 378 |
fprintf(stderr, "time interval negative\n"); |
| 379 |
exit(1); |
| 380 |
} |
| 381 |
|
| 382 |
if(fabs(solar_sunrise(month,day)-st)<=timeinterval/120) { |
| 383 |
st= (st+timeinterval/120+solar_sunrise(month,day))/2; |
| 384 |
if(suppress_warnings==0) |
| 385 |
{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } |
| 386 |
} |
| 387 |
|
| 388 |
if(fabs(solar_sunset(month,day)-st)<timeinterval/120) { |
| 389 |
st= (st-timeinterval/120+solar_sunset(month,day))/2; |
| 390 |
if(suppress_warnings==0) |
| 391 |
{ fprintf(stderr, "Solar position corrected at time step %d %d %.3f\n",month,day,hour); } |
| 392 |
} |
| 393 |
|
| 394 |
if((st<solar_sunrise(month,day)-timeinterval/120) || (st>solar_sunset(month,day)+timeinterval/120)) { |
| 395 |
if(suppress_warnings==0) |
| 396 |
{ fprintf(stderr, "Warning: sun position too low, printing error sky at %d %d %.3f\n",month,day,hour); } |
| 397 |
altitude = salt(sd, st); |
| 398 |
azimuth = sazi(sd, st); |
| 399 |
print_error_sky(); |
| 400 |
exit(0); |
| 401 |
} |
| 402 |
} |
| 403 |
else |
| 404 |
|
| 405 |
if(st<solar_sunrise(month,day) || st>solar_sunset(month,day)) { |
| 406 |
if(suppress_warnings==0) |
| 407 |
{ fprintf(stderr, "Warning: sun altitude below zero at time step %i %i %.2f, printing error sky\n",month,day,hour); } |
| 408 |
altitude = salt(sd, st); |
| 409 |
azimuth = sazi(sd, st); |
| 410 |
print_error_sky(); |
| 411 |
exit(0); |
| 412 |
} |
| 413 |
|
| 414 |
altitude = salt(sd, st); |
| 415 |
azimuth = sazi(sd, st); |
| 416 |
|
| 417 |
daynumber = (double)jdate(month, day); |
| 418 |
|
| 419 |
} |
| 420 |
|
| 421 |
|
| 422 |
|
| 423 |
|
| 424 |
|
| 425 |
if (!cloudy && altitude > 87.*M_PI/180.) { |
| 426 |
|
| 427 |
if (suppress_warnings==0) { |
| 428 |
fprintf(stderr, |
| 429 |
"%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n", |
| 430 |
progname); |
| 431 |
} |
| 432 |
altitude = 87.*M_PI/180.; |
| 433 |
} |
| 434 |
|
| 435 |
|
| 436 |
|
| 437 |
sundir[0] = -sin(azimuth)*cos(altitude); |
| 438 |
sundir[1] = -cos(azimuth)*cos(altitude); |
| 439 |
sundir[2] = sin(altitude); |
| 440 |
|
| 441 |
|
| 442 |
/* calculation for the new functions */ |
| 443 |
sunzenith = 90 - altitude*180/M_PI; |
| 444 |
|
| 445 |
|
| 446 |
/* compute the inputs for the calculation of the light distribution over the sky*/ |
| 447 |
if (input==0) /* P */ |
| 448 |
{ |
| 449 |
check_parametrization(); |
| 450 |
diffuseirradiance = diffuse_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/ |
| 451 |
directirradiance = direct_irradiance_from_sky_clearness(); |
| 452 |
check_irradiances(); |
| 453 |
|
| 454 |
if (output==0 || output==2) |
| 455 |
{ |
| 456 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
| 457 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
| 458 |
check_illuminances(); |
| 459 |
} |
| 460 |
} |
| 461 |
|
| 462 |
|
| 463 |
else if (input==1) /* W */ |
| 464 |
{ |
| 465 |
check_irradiances(); |
| 466 |
skybrightness = sky_brightness(); |
| 467 |
skyclearness = sky_clearness(); |
| 468 |
|
| 469 |
check_parametrization(); |
| 470 |
|
| 471 |
if (output==0 || output==2) |
| 472 |
{ |
| 473 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
| 474 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
| 475 |
check_illuminances(); |
| 476 |
} |
| 477 |
|
| 478 |
} |
| 479 |
|
| 480 |
|
| 481 |
else if (input==2) /* L */ |
| 482 |
{ |
| 483 |
check_illuminances(); |
| 484 |
illu_to_irra_index(); |
| 485 |
check_parametrization(); |
| 486 |
} |
| 487 |
|
| 488 |
|
| 489 |
else if (input==3) /* G */ |
| 490 |
{ |
| 491 |
if (altitude<=0) |
| 492 |
{ |
| 493 |
if (suppress_warnings==0) |
| 494 |
fprintf(stderr, "Warning: sun altitude < 0, proceed with irradiance values of zero\n"); |
| 495 |
directirradiance = 0; |
| 496 |
diffuseirradiance = 0; |
| 497 |
} else { |
| 498 |
|
| 499 |
directirradiance=directirradiance/sin(altitude); |
| 500 |
} |
| 501 |
|
| 502 |
check_irradiances(); |
| 503 |
skybrightness = sky_brightness(); |
| 504 |
skyclearness = sky_clearness(); |
| 505 |
check_parametrization(); |
| 506 |
|
| 507 |
if (output==0 || output==2) |
| 508 |
{ |
| 509 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
| 510 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
| 511 |
check_illuminances(); |
| 512 |
} |
| 513 |
|
| 514 |
} |
| 515 |
|
| 516 |
|
| 517 |
else if (input==4) /* E */ /* Implementation of the Erbs model. W.Sprenger (04/13) */ |
| 518 |
{ |
| 519 |
|
| 520 |
if (altitude<=0) |
| 521 |
{ |
| 522 |
if (suppress_warnings==0 && globalirradiance > 50) |
| 523 |
fprintf(stderr, "Warning: global irradiance higher than 50 W/m^2 while the sun altitude is lower than zero\n"); |
| 524 |
globalirradiance = 0; diffuseirradiance = 0; directirradiance = 0; |
| 525 |
|
| 526 |
} else { |
| 527 |
|
| 528 |
erbs_s0 = solar_constant_e*get_eccentricity()*sin(altitude); |
| 529 |
|
| 530 |
if (globalirradiance>erbs_s0) |
| 531 |
{ |
| 532 |
if (suppress_warnings==0) |
| 533 |
fprintf(stderr, "Warning: global irradiance is higher than the time-dependent solar constant s0\n"); |
| 534 |
globalirradiance=erbs_s0*0.999; |
| 535 |
} |
| 536 |
|
| 537 |
erbs_kt=globalirradiance/erbs_s0; |
| 538 |
|
| 539 |
if (erbs_kt<=0.22) diffuseirradiance=globalirradiance*(1-0.09*erbs_kt); |
| 540 |
else if (erbs_kt<=0.8) diffuseirradiance=globalirradiance*(0.9511-0.1604*erbs_kt+4.388*pow(erbs_kt,2)-16.638*pow(erbs_kt,3)+12.336*pow(erbs_kt,4)); |
| 541 |
else if (erbs_kt<1) diffuseirradiance=globalirradiance*(0.165); |
| 542 |
|
| 543 |
directirradiance=globalirradiance-diffuseirradiance; |
| 544 |
|
| 545 |
printf("# erbs_s0, erbs_kt, irr_dir_h, irr_diff: %.3f %.3f %.3f %.3f\n", erbs_s0, erbs_kt, directirradiance, diffuseirradiance); |
| 546 |
printf("# WARNING: the -E option is only recommended for a rough estimation!\n"); |
| 547 |
|
| 548 |
directirradiance=directirradiance/sin(altitude); |
| 549 |
|
| 550 |
} |
| 551 |
|
| 552 |
check_irradiances(); |
| 553 |
skybrightness = sky_brightness(); |
| 554 |
skyclearness = sky_clearness(); |
| 555 |
check_parametrization(); |
| 556 |
|
| 557 |
if (output==0 || output==2) |
| 558 |
{ |
| 559 |
diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/ |
| 560 |
directilluminance = directirradiance*direct_n_effi_PEREZ(); |
| 561 |
check_illuminances(); |
| 562 |
} |
| 563 |
|
| 564 |
} |
| 565 |
|
| 566 |
|
| 567 |
|
| 568 |
|
| 569 |
else { fprintf(stderr,"error at the input arguments"); exit(1); } |
| 570 |
|
| 571 |
|
| 572 |
|
| 573 |
/* normalization factor for the relative sky luminance distribution, diffuse part*/ |
| 574 |
|
| 575 |
if ( (lv_mod = malloc(145*sizeof(float))) == NULL) |
| 576 |
{ |
| 577 |
fprintf(stderr,"Out of memory in function main"); |
| 578 |
exit(1); |
| 579 |
} |
| 580 |
|
| 581 |
/* read the angles */ |
| 582 |
theta_o = defangle_theta; |
| 583 |
phi_o = defangle_phi; |
| 584 |
|
| 585 |
|
| 586 |
/* parameters for the perez model */ |
| 587 |
coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez); |
| 588 |
|
| 589 |
|
| 590 |
|
| 591 |
/*calculation of the modelled luminance */ |
| 592 |
for (j=0;j<145;j++) |
| 593 |
{ |
| 594 |
theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith)); |
| 595 |
|
| 596 |
*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
| 597 |
|
| 598 |
/* fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j)); */ |
| 599 |
} |
| 600 |
|
| 601 |
/* integration of luminance for the normalization factor, diffuse part of the sky*/ |
| 602 |
|
| 603 |
diffnormalization = integ_lv(lv_mod, theta_o); |
| 604 |
|
| 605 |
|
| 606 |
|
| 607 |
/*normalization coefficient in lumen or in watt*/ |
| 608 |
if (output==0) |
| 609 |
{ |
| 610 |
diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY; |
| 611 |
} |
| 612 |
else if (output==1) |
| 613 |
{ |
| 614 |
diffnormalization = diffuseirradiance/diffnormalization; |
| 615 |
} |
| 616 |
else if (output==2) |
| 617 |
{ |
| 618 |
diffnormalization = diffuseilluminance/diffnormalization; |
| 619 |
} |
| 620 |
|
| 621 |
else {fprintf(stderr,"Wrong output specification.\n"); exit(1);} |
| 622 |
|
| 623 |
|
| 624 |
|
| 625 |
|
| 626 |
/* calculation for the solar source */ |
| 627 |
if (output==0) |
| 628 |
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY; |
| 629 |
|
| 630 |
else if (output==1) |
| 631 |
solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
| 632 |
|
| 633 |
else |
| 634 |
solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180))); |
| 635 |
|
| 636 |
|
| 637 |
|
| 638 |
/* Compute the ground radiance */ |
| 639 |
zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez); |
| 640 |
zenithbr*=diffnormalization; |
| 641 |
|
| 642 |
if (skyclearness==1) |
| 643 |
normfactor = 0.777778; |
| 644 |
|
| 645 |
if (skyclearness>=6) |
| 646 |
{ |
| 647 |
F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]); |
| 648 |
normfactor = normsc()/F2/M_PI; |
| 649 |
} |
| 650 |
|
| 651 |
if ( (skyclearness>1) && (skyclearness<6) ) |
| 652 |
{ |
| 653 |
S_INTER=1; |
| 654 |
F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude)); |
| 655 |
normfactor = normsc()/F2/M_PI; |
| 656 |
} |
| 657 |
|
| 658 |
groundbr = zenithbr*normfactor; |
| 659 |
|
| 660 |
if (dosun&&(skyclearness>1)) |
| 661 |
groundbr += 6.8e-5/M_PI*solarradiance*sundir[2]; |
| 662 |
|
| 663 |
groundbr *= gprefl; |
| 664 |
|
| 665 |
|
| 666 |
|
| 667 |
if(*(c_perez+1)>0) |
| 668 |
{ |
| 669 |
if(suppress_warnings==0) |
| 670 |
{ fprintf(stderr, "Warning: positive Perez parameter B (= %lf), printing error sky\n",*(c_perez+1));} |
| 671 |
print_error_sky(); |
| 672 |
exit(0); |
| 673 |
} |
| 674 |
|
| 675 |
|
| 676 |
return; |
| 677 |
} |
| 678 |
|
| 679 |
|
| 680 |
|
| 681 |
|
| 682 |
|
| 683 |
double solar_sunset(int month,int day) |
| 684 |
{ |
| 685 |
float W; |
| 686 |
extern double s_latitude; |
| 687 |
W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); |
| 688 |
return(12+(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); |
| 689 |
} |
| 690 |
|
| 691 |
|
| 692 |
|
| 693 |
|
| 694 |
double solar_sunrise(int month,int day) |
| 695 |
{ |
| 696 |
float W; |
| 697 |
extern double s_latitude; |
| 698 |
W=-1*(tan(s_latitude)*tan(sdec(jdate(month, day)))); |
| 699 |
return(12-(M_PI/2 - atan2(W,sqrt(1-W*W)))*180/(M_PI*15)); |
| 700 |
} |
| 701 |
|
| 702 |
|
| 703 |
|
| 704 |
|
| 705 |
void printsky() |
| 706 |
{ |
| 707 |
|
| 708 |
printf("# Local solar time: %.2f\n", st); |
| 709 |
printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); |
| 710 |
printf("# epsilon, delta, atmospheric precipitable water content : %.4f %.4f %.4f \n", skyclearness, skybrightness,atm_preci_water ); |
| 711 |
|
| 712 |
|
| 713 |
if (dosun&&(skyclearness>1)) |
| 714 |
{ |
| 715 |
printf("\nvoid light solar\n"); |
| 716 |
printf("0\n0\n"); |
| 717 |
printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance); |
| 718 |
printf("\nsolar source sun\n"); |
| 719 |
printf("0\n0\n"); |
| 720 |
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
| 721 |
} else if (dosun) { |
| 722 |
printf("\nvoid light solar\n"); |
| 723 |
printf("0\n0\n"); |
| 724 |
printf("3 0.0 0.0 0.0\n"); |
| 725 |
printf("\nsolar source sun\n"); |
| 726 |
printf("0\n0\n"); |
| 727 |
printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle); |
| 728 |
} |
| 729 |
/* print colored output if activated in command line (-C). Based on model from A. Diakite, TU-Berlin. Implemented by J. Wienold, August 26 2018 */ |
| 730 |
if (color_output==1 && skyclearness < 4.5 && skyclearness >1.065 ) |
| 731 |
{ |
| 732 |
fprintf(stderr, " warning: sky clearness(epsilon)= %f \n",skyclearness); |
| 733 |
fprintf(stderr, " warning: intermediate sky!! \n"); |
| 734 |
fprintf(stderr, " warning: color model for intermediate sky pending \n"); |
| 735 |
fprintf(stderr, " warning: no color output ! \n"); |
| 736 |
color_output=0; |
| 737 |
} |
| 738 |
if (color_output==1) |
| 739 |
{ |
| 740 |
printf("\nvoid colorfunc skyfunc\n"); |
| 741 |
printf("4 skybright_r skybright_g skybright_b perezlum_c.cal\n"); |
| 742 |
printf("0\n"); |
| 743 |
printf("22 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", diffnormalization, groundbr, |
| 744 |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
| 745 |
sundir[0], sundir[1], sundir[2],skyclearness,locus[0],locus[1],locus[2],locus[3],locus[4],locus[5],locus[6],locus[7],locus[8],locus[9],locus[10]); |
| 746 |
}else{ |
| 747 |
printf("\nvoid brightfunc skyfunc\n"); |
| 748 |
printf("2 skybright perezlum.cal\n"); |
| 749 |
printf("0\n"); |
| 750 |
printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr, |
| 751 |
*(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), |
| 752 |
sundir[0], sundir[1], sundir[2]); |
| 753 |
} |
| 754 |
|
| 755 |
} |
| 756 |
|
| 757 |
|
| 758 |
|
| 759 |
void print_error_sky() |
| 760 |
{ |
| 761 |
|
| 762 |
|
| 763 |
sundir[0] = -sin(azimuth)*cos(altitude); |
| 764 |
sundir[1] = -cos(azimuth)*cos(altitude); |
| 765 |
sundir[2] = sin(altitude); |
| 766 |
|
| 767 |
printf("# Local solar time: %.2f\n", st); |
| 768 |
printf("# Solar altitude and azimuth: %.1f %.1f\n", altitude*180/M_PI, azimuth*180/M_PI); |
| 769 |
|
| 770 |
printf("\nvoid brightfunc skyfunc\n"); |
| 771 |
printf("2 skybright perezlum.cal\n"); |
| 772 |
printf("0\n"); |
| 773 |
printf("10 0.00 0.00 0.000 0.000 0.000 0.000 0.000 %f %f %f \n", sundir[0], sundir[1], sundir[2]); |
| 774 |
} |
| 775 |
|
| 776 |
|
| 777 |
|
| 778 |
|
| 779 |
|
| 780 |
void printdefaults() /* print default values */ |
| 781 |
{ |
| 782 |
printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl); |
| 783 |
if (zenithbr > 0.0) |
| 784 |
printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr); |
| 785 |
else |
| 786 |
printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity); |
| 787 |
printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI)); |
| 788 |
printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI)); |
| 789 |
printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI)); |
| 790 |
} |
| 791 |
|
| 792 |
|
| 793 |
|
| 794 |
|
| 795 |
void usage_error(char* msg) /* print usage error and quit */ |
| 796 |
{ |
| 797 |
if (msg != NULL) |
| 798 |
fprintf(stderr, "%s: Use error - %s\n\n", progname, msg); |
| 799 |
fprintf(stderr, "Usage: %s month day hour [-y year] [...]\n", progname); |
| 800 |
fprintf(stderr, " or: %s -ang altitude azimuth [...]\n", progname); |
| 801 |
fprintf(stderr, " followed by: -P epsilon delta [options]\n"); |
| 802 |
fprintf(stderr, " or: [-W|-L|-G] direct_value diffuse_value [options]\n"); |
| 803 |
fprintf(stderr, " or: -E global_irradiance [options]\n\n"); |
| 804 |
fprintf(stderr, " Description:\n"); |
| 805 |
fprintf(stderr, " -P epsilon delta (these are the Perez parameters) \n"); |
| 806 |
fprintf(stderr, " -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
| 807 |
fprintf(stderr, " -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n"); |
| 808 |
fprintf(stderr, " -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n"); |
| 809 |
fprintf(stderr, " -E global-horizontal-irradiance (W/m^2)\n\n"); |
| 810 |
fprintf(stderr, " Output specification with option:\n"); |
| 811 |
fprintf(stderr, " -O [0|1|2] (0=output in W/m^2/sr visible, 1=output in W/m^2/sr solar, 2=output in candela/m^2), default is 0 \n"); |
| 812 |
fprintf(stderr, " gendaylit version 2.5 (2018/04/18) \n\n"); |
| 813 |
exit(1); |
| 814 |
} |
| 815 |
|
| 816 |
|
| 817 |
|
| 818 |
|
| 819 |
double normsc() /* compute normalization factor (E0*F2/L0) */ |
| 820 |
{ |
| 821 |
static double nfc[2][5] = { |
| 822 |
/* clear sky approx. */ |
| 823 |
{2.766521, 0.547665, -0.369832, 0.009237, 0.059229}, |
| 824 |
/* intermediate sky approx. */ |
| 825 |
{3.5556, -2.7152, -1.3081, 1.0660, 0.60227}, |
| 826 |
}; |
| 827 |
register double *nf; |
| 828 |
double x, nsc; |
| 829 |
register int i; |
| 830 |
/* polynomial approximation */ |
| 831 |
nf = nfc[S_INTER]; |
| 832 |
x = (altitude - M_PI/4.0)/(M_PI/4.0); |
| 833 |
nsc = nf[i=4]; |
| 834 |
while (i--) |
| 835 |
nsc = nsc*x + nf[i]; |
| 836 |
|
| 837 |
return(nsc); |
| 838 |
} |
| 839 |
|
| 840 |
|
| 841 |
|
| 842 |
|
| 843 |
|
| 844 |
void printhead(int ac, char** av) /* print command header */ |
| 845 |
{ |
| 846 |
putchar('#'); |
| 847 |
while (ac--) { |
| 848 |
putchar(' '); |
| 849 |
fputs(*av++, stdout); |
| 850 |
} |
| 851 |
putchar('\n'); |
| 852 |
} |
| 853 |
|
| 854 |
|
| 855 |
|
| 856 |
|
| 857 |
|
| 858 |
|
| 859 |
/* Perez models */ |
| 860 |
|
| 861 |
/* Perez global horizontal luminous efficacy model */ |
| 862 |
double glob_h_effi_PEREZ() |
| 863 |
{ |
| 864 |
|
| 865 |
double value; |
| 866 |
double category_bounds[10], a[10], b[10], c[10], d[10]; |
| 867 |
int category_total_number, category_number, i; |
| 868 |
|
| 869 |
check_parametrization(); |
| 870 |
|
| 871 |
|
| 872 |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
| 873 |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_effi_PEREZ \n"); */ |
| 874 |
|
| 875 |
|
| 876 |
/* initialize category bounds (clearness index bounds) */ |
| 877 |
|
| 878 |
category_total_number = 8; |
| 879 |
|
| 880 |
category_bounds[1] = 1; |
| 881 |
category_bounds[2] = 1.065; |
| 882 |
category_bounds[3] = 1.230; |
| 883 |
category_bounds[4] = 1.500; |
| 884 |
category_bounds[5] = 1.950; |
| 885 |
category_bounds[6] = 2.800; |
| 886 |
category_bounds[7] = 4.500; |
| 887 |
category_bounds[8] = 6.200; |
| 888 |
category_bounds[9] = 12.01; |
| 889 |
|
| 890 |
|
| 891 |
/* initialize model coefficients */ |
| 892 |
a[1] = 96.63; |
| 893 |
a[2] = 107.54; |
| 894 |
a[3] = 98.73; |
| 895 |
a[4] = 92.72; |
| 896 |
a[5] = 86.73; |
| 897 |
a[6] = 88.34; |
| 898 |
a[7] = 78.63; |
| 899 |
a[8] = 99.65; |
| 900 |
|
| 901 |
b[1] = -0.47; |
| 902 |
b[2] = 0.79; |
| 903 |
b[3] = 0.70; |
| 904 |
b[4] = 0.56; |
| 905 |
b[5] = 0.98; |
| 906 |
b[6] = 1.39; |
| 907 |
b[7] = 1.47; |
| 908 |
b[8] = 1.86; |
| 909 |
|
| 910 |
c[1] = 11.50; |
| 911 |
c[2] = 1.79; |
| 912 |
c[3] = 4.40; |
| 913 |
c[4] = 8.36; |
| 914 |
c[5] = 7.10; |
| 915 |
c[6] = 6.06; |
| 916 |
c[7] = 4.93; |
| 917 |
c[8] = -4.46; |
| 918 |
|
| 919 |
d[1] = -9.16; |
| 920 |
d[2] = -1.19; |
| 921 |
d[3] = -6.95; |
| 922 |
d[4] = -8.31; |
| 923 |
d[5] = -10.94; |
| 924 |
d[6] = -7.60; |
| 925 |
d[7] = -11.37; |
| 926 |
d[8] = -3.15; |
| 927 |
|
| 928 |
|
| 929 |
|
| 930 |
for (i=1; i<=category_total_number; i++) |
| 931 |
{ |
| 932 |
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
| 933 |
category_number = i; |
| 934 |
} |
| 935 |
|
| 936 |
value = a[category_number] + b[category_number]*atm_preci_water + |
| 937 |
c[category_number]*cos(sunzenith*M_PI/180) + d[category_number]*log(skybrightness); |
| 938 |
|
| 939 |
return(value); |
| 940 |
} |
| 941 |
|
| 942 |
|
| 943 |
|
| 944 |
|
| 945 |
/* global horizontal diffuse efficacy model, according to PEREZ */ |
| 946 |
double glob_h_diffuse_effi_PEREZ() |
| 947 |
{ |
| 948 |
double value; |
| 949 |
double category_bounds[10], a[10], b[10], c[10], d[10]; |
| 950 |
int category_total_number, category_number, i; |
| 951 |
|
| 952 |
check_parametrization(); |
| 953 |
|
| 954 |
|
| 955 |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
| 956 |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */ |
| 957 |
|
| 958 |
/* initialize category bounds (clearness index bounds) */ |
| 959 |
|
| 960 |
category_total_number = 8; |
| 961 |
|
| 962 |
//XXX: category_bounds > 0.1 |
| 963 |
category_bounds[1] = 1; |
| 964 |
category_bounds[2] = 1.065; |
| 965 |
category_bounds[3] = 1.230; |
| 966 |
category_bounds[4] = 1.500; |
| 967 |
category_bounds[5] = 1.950; |
| 968 |
category_bounds[6] = 2.800; |
| 969 |
category_bounds[7] = 4.500; |
| 970 |
category_bounds[8] = 6.200; |
| 971 |
category_bounds[9] = 12.01; |
| 972 |
|
| 973 |
|
| 974 |
/* initialize model coefficients */ |
| 975 |
a[1] = 97.24; |
| 976 |
a[2] = 107.22; |
| 977 |
a[3] = 104.97; |
| 978 |
a[4] = 102.39; |
| 979 |
a[5] = 100.71; |
| 980 |
a[6] = 106.42; |
| 981 |
a[7] = 141.88; |
| 982 |
a[8] = 152.23; |
| 983 |
|
| 984 |
b[1] = -0.46; |
| 985 |
b[2] = 1.15; |
| 986 |
b[3] = 2.96; |
| 987 |
b[4] = 5.59; |
| 988 |
b[5] = 5.94; |
| 989 |
b[6] = 3.83; |
| 990 |
b[7] = 1.90; |
| 991 |
b[8] = 0.35; |
| 992 |
|
| 993 |
c[1] = 12.00; |
| 994 |
c[2] = 0.59; |
| 995 |
c[3] = -5.53; |
| 996 |
c[4] = -13.95; |
| 997 |
c[5] = -22.75; |
| 998 |
c[6] = -36.15; |
| 999 |
c[7] = -53.24; |
| 1000 |
c[8] = -45.27; |
| 1001 |
|
| 1002 |
d[1] = -8.91; |
| 1003 |
d[2] = -3.95; |
| 1004 |
d[3] = -8.77; |
| 1005 |
d[4] = -13.90; |
| 1006 |
d[5] = -23.74; |
| 1007 |
d[6] = -28.83; |
| 1008 |
d[7] = -14.03; |
| 1009 |
d[8] = -7.98; |
| 1010 |
|
| 1011 |
|
| 1012 |
|
| 1013 |
category_number = -1; |
| 1014 |
for (i=1; i<=category_total_number; i++) |
| 1015 |
{ |
| 1016 |
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
| 1017 |
category_number = i; |
| 1018 |
} |
| 1019 |
|
| 1020 |
if (category_number == -1) { |
| 1021 |
if (suppress_warnings==0) |
| 1022 |
fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness); |
| 1023 |
print_error_sky(); |
| 1024 |
exit(0); |
| 1025 |
} |
| 1026 |
|
| 1027 |
|
| 1028 |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*cos(sunzenith*M_PI/180) + |
| 1029 |
d[category_number]*log(skybrightness); |
| 1030 |
|
| 1031 |
return(value); |
| 1032 |
|
| 1033 |
} |
| 1034 |
|
| 1035 |
|
| 1036 |
|
| 1037 |
|
| 1038 |
|
| 1039 |
|
| 1040 |
/* direct normal efficacy model, according to PEREZ */ |
| 1041 |
|
| 1042 |
double direct_n_effi_PEREZ() |
| 1043 |
|
| 1044 |
{ |
| 1045 |
double value; |
| 1046 |
double category_bounds[10], a[10], b[10], c[10], d[10]; |
| 1047 |
int category_total_number, category_number, i; |
| 1048 |
|
| 1049 |
|
| 1050 |
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0) |
| 1051 |
fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/ |
| 1052 |
|
| 1053 |
|
| 1054 |
/* initialize category bounds (clearness index bounds) */ |
| 1055 |
|
| 1056 |
category_total_number = 8; |
| 1057 |
|
| 1058 |
category_bounds[1] = 1; |
| 1059 |
category_bounds[2] = 1.065; |
| 1060 |
category_bounds[3] = 1.230; |
| 1061 |
category_bounds[4] = 1.500; |
| 1062 |
category_bounds[5] = 1.950; |
| 1063 |
category_bounds[6] = 2.800; |
| 1064 |
category_bounds[7] = 4.500; |
| 1065 |
category_bounds[8] = 6.200; |
| 1066 |
category_bounds[9] = 12.1; |
| 1067 |
|
| 1068 |
|
| 1069 |
/* initialize model coefficients */ |
| 1070 |
a[1] = 57.20; |
| 1071 |
a[2] = 98.99; |
| 1072 |
a[3] = 109.83; |
| 1073 |
a[4] = 110.34; |
| 1074 |
a[5] = 106.36; |
| 1075 |
a[6] = 107.19; |
| 1076 |
a[7] = 105.75; |
| 1077 |
a[8] = 101.18; |
| 1078 |
|
| 1079 |
b[1] = -4.55; |
| 1080 |
b[2] = -3.46; |
| 1081 |
b[3] = -4.90; |
| 1082 |
b[4] = -5.84; |
| 1083 |
b[5] = -3.97; |
| 1084 |
b[6] = -1.25; |
| 1085 |
b[7] = 0.77; |
| 1086 |
b[8] = 1.58; |
| 1087 |
|
| 1088 |
c[1] = -2.98; |
| 1089 |
c[2] = -1.21; |
| 1090 |
c[3] = -1.71; |
| 1091 |
c[4] = -1.99; |
| 1092 |
c[5] = -1.75; |
| 1093 |
c[6] = -1.51; |
| 1094 |
c[7] = -1.26; |
| 1095 |
c[8] = -1.10; |
| 1096 |
|
| 1097 |
d[1] = 117.12; |
| 1098 |
d[2] = 12.38; |
| 1099 |
d[3] = -8.81; |
| 1100 |
d[4] = -4.56; |
| 1101 |
d[5] = -6.16; |
| 1102 |
d[6] = -26.73; |
| 1103 |
d[7] = -34.44; |
| 1104 |
d[8] = -8.29; |
| 1105 |
|
| 1106 |
|
| 1107 |
|
| 1108 |
for (i=1; i<=category_total_number; i++) |
| 1109 |
{ |
| 1110 |
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) ) |
| 1111 |
category_number = i; |
| 1112 |
} |
| 1113 |
|
| 1114 |
value = a[category_number] + b[category_number]*atm_preci_water + c[category_number]*exp(5.73*sunzenith*M_PI/180 - 5) + d[category_number]*skybrightness; |
| 1115 |
|
| 1116 |
if (value < 0) value = 0; |
| 1117 |
|
| 1118 |
return(value); |
| 1119 |
} |
| 1120 |
|
| 1121 |
|
| 1122 |
/*check the range of epsilon and delta indexes of the perez parametrization*/ |
| 1123 |
void check_parametrization() |
| 1124 |
{ |
| 1125 |
|
| 1126 |
if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) |
| 1127 |
{ |
| 1128 |
|
| 1129 |
/* limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */ |
| 1130 |
|
| 1131 |
if (skyclearness<skyclearinf){ |
| 1132 |
/* if (suppress_warnings==0) |
| 1133 |
fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */ |
| 1134 |
skyclearness=skyclearinf; |
| 1135 |
} |
| 1136 |
if (skyclearness>skyclearsup){ |
| 1137 |
/* if (suppress_warnings==0) |
| 1138 |
fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */ |
| 1139 |
skyclearness=skyclearsup-0.001; |
| 1140 |
} |
| 1141 |
if (skybrightness<skybriginf){ |
| 1142 |
/* if (suppress_warnings==0) |
| 1143 |
fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */ |
| 1144 |
skybrightness=skybriginf; |
| 1145 |
} |
| 1146 |
if (skybrightness>skybrigsup){ |
| 1147 |
/* if (suppress_warnings==0) |
| 1148 |
fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */ |
| 1149 |
skybrightness=skybrigsup; |
| 1150 |
} |
| 1151 |
|
| 1152 |
return; } |
| 1153 |
else return; |
| 1154 |
} |
| 1155 |
|
| 1156 |
|
| 1157 |
|
| 1158 |
|
| 1159 |
|
| 1160 |
/* validity of the direct and diffuse components */ |
| 1161 |
void check_illuminances() |
| 1162 |
{ |
| 1163 |
if (directilluminance < 0) { |
| 1164 |
if(suppress_warnings==0) |
| 1165 |
{ fprintf(stderr,"Warning: direct illuminance < 0. Using 0.0\n"); } |
| 1166 |
directilluminance = 0.0; |
| 1167 |
} |
| 1168 |
if (diffuseilluminance < 0) { |
| 1169 |
if(suppress_warnings==0) |
| 1170 |
{ fprintf(stderr,"Warning: diffuse illuminance < 0. Using 0.0\n"); } |
| 1171 |
diffuseilluminance = 0.0; |
| 1172 |
} |
| 1173 |
|
| 1174 |
if (directilluminance+diffuseilluminance==0 && altitude > 0) { |
| 1175 |
if(suppress_warnings==0) |
| 1176 |
{ fprintf(stderr,"Warning: zero illuminance at sun altitude > 0, printing error sky\n"); } |
| 1177 |
print_error_sky(); |
| 1178 |
exit(0); |
| 1179 |
} |
| 1180 |
|
| 1181 |
if (directilluminance > solar_constant_l) { |
| 1182 |
if(suppress_warnings==0) |
| 1183 |
{ fprintf(stderr,"Warning: direct illuminance exceeds solar constant\n"); } |
| 1184 |
print_error_sky(); |
| 1185 |
exit(0); |
| 1186 |
} |
| 1187 |
} |
| 1188 |
|
| 1189 |
|
| 1190 |
void check_irradiances() |
| 1191 |
{ |
| 1192 |
if (directirradiance < 0) { |
| 1193 |
if(suppress_warnings==0) |
| 1194 |
{ fprintf(stderr,"Warning: direct irradiance < 0. Using 0.0\n"); } |
| 1195 |
directirradiance = 0.0; |
| 1196 |
} |
| 1197 |
if (diffuseirradiance < 0) { |
| 1198 |
if(suppress_warnings==0) |
| 1199 |
{ fprintf(stderr,"Warning: diffuse irradiance < 0. Using 0.0\n"); } |
| 1200 |
diffuseirradiance = 0.0; |
| 1201 |
} |
| 1202 |
|
| 1203 |
if (directirradiance+diffuseirradiance==0 && altitude > 0) { |
| 1204 |
if(suppress_warnings==0) |
| 1205 |
{ fprintf(stderr,"Warning: zero irradiance at sun altitude > 0, printing error sky\n"); } |
| 1206 |
print_error_sky(); |
| 1207 |
exit(0); |
| 1208 |
} |
| 1209 |
|
| 1210 |
if (directirradiance > solar_constant_e) { |
| 1211 |
if(suppress_warnings==0) |
| 1212 |
{ fprintf(stderr,"Warning: direct irradiance exceeds solar constant\n"); } |
| 1213 |
print_error_sky(); |
| 1214 |
exit(0); |
| 1215 |
} |
| 1216 |
} |
| 1217 |
|
| 1218 |
|
| 1219 |
|
| 1220 |
/* Perez sky's brightness */ |
| 1221 |
double sky_brightness() |
| 1222 |
{ |
| 1223 |
double value; |
| 1224 |
|
| 1225 |
value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity()); |
| 1226 |
|
| 1227 |
return(value); |
| 1228 |
} |
| 1229 |
|
| 1230 |
|
| 1231 |
/* Perez sky's clearness */ |
| 1232 |
double sky_clearness() |
| 1233 |
{ |
| 1234 |
double value; |
| 1235 |
|
| 1236 |
value = ( (diffuseirradiance + directirradiance)/(diffuseirradiance) + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180 ) / (1 + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ; |
| 1237 |
|
| 1238 |
return(value); |
| 1239 |
} |
| 1240 |
|
| 1241 |
|
| 1242 |
|
| 1243 |
/* diffus horizontal irradiance from Perez sky's brightness */ |
| 1244 |
double diffuse_irradiance_from_sky_brightness() |
| 1245 |
{ |
| 1246 |
double value; |
| 1247 |
|
| 1248 |
value = skybrightness / air_mass() * ( solar_constant_e*get_eccentricity()); |
| 1249 |
|
| 1250 |
return(value); |
| 1251 |
} |
| 1252 |
|
| 1253 |
|
| 1254 |
/* direct normal irradiance from Perez sky's clearness */ |
| 1255 |
double direct_irradiance_from_sky_clearness() |
| 1256 |
{ |
| 1257 |
double value; |
| 1258 |
|
| 1259 |
value = diffuse_irradiance_from_sky_brightness(); |
| 1260 |
value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ); |
| 1261 |
|
| 1262 |
return(value); |
| 1263 |
} |
| 1264 |
|
| 1265 |
|
| 1266 |
|
| 1267 |
|
| 1268 |
void illu_to_irra_index() |
| 1269 |
{ |
| 1270 |
double test1=0.1, test2=0.1, d_eff; |
| 1271 |
int counter=0; |
| 1272 |
|
| 1273 |
diffuseirradiance = diffuseilluminance*solar_constant_e/(solar_constant_l); |
| 1274 |
directirradiance = directilluminance*solar_constant_e/(solar_constant_l); |
| 1275 |
skyclearness = sky_clearness(); |
| 1276 |
skybrightness = sky_brightness(); |
| 1277 |
check_parametrization(); |
| 1278 |
|
| 1279 |
|
| 1280 |
while ( ((fabs(diffuseirradiance-test1)>10) || (fabs(directirradiance-test2)>10) |
| 1281 |
|| (!(skyclearness<skyclearinf || skyclearness>skyclearsup)) |
| 1282 |
|| (!(skybrightness<skybriginf || skybrightness>skybrigsup)) ) |
| 1283 |
&& !(counter==9) ) |
| 1284 |
{ |
| 1285 |
|
| 1286 |
test1=diffuseirradiance; |
| 1287 |
test2=directirradiance; |
| 1288 |
counter++; |
| 1289 |
|
| 1290 |
diffuseirradiance = diffuseilluminance/glob_h_diffuse_effi_PEREZ(); |
| 1291 |
d_eff = direct_n_effi_PEREZ(); |
| 1292 |
|
| 1293 |
|
| 1294 |
if (d_eff < 0.1) |
| 1295 |
directirradiance = 0; |
| 1296 |
else |
| 1297 |
directirradiance = directilluminance/d_eff; |
| 1298 |
|
| 1299 |
skybrightness = sky_brightness(); |
| 1300 |
skyclearness = sky_clearness(); |
| 1301 |
check_parametrization(); |
| 1302 |
|
| 1303 |
} |
| 1304 |
|
| 1305 |
|
| 1306 |
return; |
| 1307 |
} |
| 1308 |
|
| 1309 |
static int get_numlin(float epsilon) |
| 1310 |
{ |
| 1311 |
if (epsilon < 1.065) |
| 1312 |
return 0; |
| 1313 |
else if (epsilon < 1.230) |
| 1314 |
return 1; |
| 1315 |
else if (epsilon < 1.500) |
| 1316 |
return 2; |
| 1317 |
else if (epsilon < 1.950) |
| 1318 |
return 3; |
| 1319 |
else if (epsilon < 2.800) |
| 1320 |
return 4; |
| 1321 |
else if (epsilon < 4.500) |
| 1322 |
return 5; |
| 1323 |
else if (epsilon < 6.200) |
| 1324 |
return 6; |
| 1325 |
return 7; |
| 1326 |
} |
| 1327 |
|
| 1328 |
/* sky luminance perez model */ |
| 1329 |
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]) |
| 1330 |
{ |
| 1331 |
|
| 1332 |
float x[5][4]; |
| 1333 |
int i,j,num_lin; |
| 1334 |
double c_perez[5]; |
| 1335 |
|
| 1336 |
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
| 1337 |
{ |
| 1338 |
fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n"); |
| 1339 |
exit(1); |
| 1340 |
} |
| 1341 |
|
| 1342 |
/* correction de modele de Perez solar energy ...*/ |
| 1343 |
if ( (epsilon > 1.065) && (epsilon < 2.8) ) |
| 1344 |
{ |
| 1345 |
if ( Delta < 0.2 ) Delta = 0.2; |
| 1346 |
} |
| 1347 |
|
| 1348 |
|
| 1349 |
num_lin = get_numlin(epsilon); |
| 1350 |
|
| 1351 |
for (i=0;i<5;i++) |
| 1352 |
for (j=0;j<4;j++) |
| 1353 |
{ |
| 1354 |
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
| 1355 |
/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */ |
| 1356 |
} |
| 1357 |
|
| 1358 |
|
| 1359 |
if (num_lin) |
| 1360 |
{ |
| 1361 |
for (i=0;i<5;i++) |
| 1362 |
c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z); |
| 1363 |
} |
| 1364 |
else |
| 1365 |
{ |
| 1366 |
c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z); |
| 1367 |
c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z); |
| 1368 |
c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z); |
| 1369 |
c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3]; |
| 1370 |
c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3]; |
| 1371 |
} |
| 1372 |
|
| 1373 |
|
| 1374 |
return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) * |
| 1375 |
(1 + c_perez[2]*exp(c_perez[3]*gamma) + |
| 1376 |
c_perez[4]*cos(gamma)*cos(gamma) ); |
| 1377 |
} |
| 1378 |
|
| 1379 |
|
| 1380 |
|
| 1381 |
/* coefficients for the sky luminance perez model */ |
| 1382 |
void coeff_lum_perez(double Z, double epsilon, double Delta, float coeff_perez[]) |
| 1383 |
{ |
| 1384 |
float x[5][4]; |
| 1385 |
int i,j,num_lin; |
| 1386 |
|
| 1387 |
if ( (epsilon < skyclearinf) || (epsilon >= skyclearsup) ) |
| 1388 |
{ |
| 1389 |
fprintf(stderr,"Error: epsilon out of range in function coeff_lum_perez!\n"); |
| 1390 |
exit(1); |
| 1391 |
} |
| 1392 |
|
| 1393 |
/* correction du modele de Perez solar energy ...*/ |
| 1394 |
if ( (epsilon > 1.065) && (epsilon < 2.8) ) |
| 1395 |
{ |
| 1396 |
if ( Delta < 0.2 ) Delta = 0.2; |
| 1397 |
} |
| 1398 |
|
| 1399 |
|
| 1400 |
num_lin = get_numlin(epsilon); |
| 1401 |
|
| 1402 |
/*fprintf(stderr,"numlin %d\n", num_lin);*/ |
| 1403 |
|
| 1404 |
for (i=0;i<5;i++) |
| 1405 |
for (j=0;j<4;j++) |
| 1406 |
{ |
| 1407 |
x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j); |
| 1408 |
/* printf("x %d %d vaut %f\n",i,j,x[i][j]); */ |
| 1409 |
} |
| 1410 |
|
| 1411 |
|
| 1412 |
if (num_lin) |
| 1413 |
{ |
| 1414 |
for (i=0;i<5;i++) |
| 1415 |
*(c_perez+i) = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z); |
| 1416 |
|
| 1417 |
} |
| 1418 |
else |
| 1419 |
{ |
| 1420 |
*(c_perez+0) = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z); |
| 1421 |
*(c_perez+1) = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z); |
| 1422 |
*(c_perez+4) = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z); |
| 1423 |
*(c_perez+2) = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3]; |
| 1424 |
*(c_perez+3) = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3]; |
| 1425 |
|
| 1426 |
|
| 1427 |
} |
| 1428 |
|
| 1429 |
|
| 1430 |
return; |
| 1431 |
} |
| 1432 |
|
| 1433 |
|
| 1434 |
|
| 1435 |
/* degrees into radians */ |
| 1436 |
double radians(double degres) |
| 1437 |
{ |
| 1438 |
return degres*(M_PI/180.); |
| 1439 |
} |
| 1440 |
|
| 1441 |
|
| 1442 |
/* radian into degrees */ |
| 1443 |
double degres(double radians) |
| 1444 |
{ |
| 1445 |
return radians*(180./M_PI); |
| 1446 |
} |
| 1447 |
|
| 1448 |
|
| 1449 |
/* calculation of the angles dzeta and gamma */ |
| 1450 |
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z) |
| 1451 |
{ |
| 1452 |
*dzeta = theta; /* dzeta = phi */ |
| 1453 |
if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) ) |
| 1454 |
*gamma = 0; |
| 1455 |
else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 ) |
| 1456 |
{ |
| 1457 |
printf("error in calculation of gamma (angle between point and sun"); |
| 1458 |
exit(1); |
| 1459 |
} |
| 1460 |
else |
| 1461 |
*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)); |
| 1462 |
} |
| 1463 |
|
| 1464 |
|
| 1465 |
|
| 1466 |
double integ_lv(float *lv,float *theta) |
| 1467 |
{ |
| 1468 |
int i; |
| 1469 |
double buffer=0.0; |
| 1470 |
|
| 1471 |
for (i=0;i<145;i++) |
| 1472 |
{ |
| 1473 |
buffer += (*(lv+i))*cos(radians(*(theta+i))); |
| 1474 |
} |
| 1475 |
|
| 1476 |
return buffer*(2.*M_PI/145.); |
| 1477 |
} |
| 1478 |
|
| 1479 |
|
| 1480 |
|
| 1481 |
/* enter day number(double), return E0 = square(R0/R): eccentricity correction factor */ |
| 1482 |
|
| 1483 |
double get_eccentricity() |
| 1484 |
{ |
| 1485 |
double day_angle; |
| 1486 |
double E0; |
| 1487 |
|
| 1488 |
day_angle = 2*M_PI*(daynumber - 1)/365; |
| 1489 |
E0 = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+ |
| 1490 |
0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle); |
| 1491 |
|
| 1492 |
return (E0); |
| 1493 |
} |
| 1494 |
|
| 1495 |
|
| 1496 |
/* enter sunzenith angle (degrees) return relative air mass (double) */ |
| 1497 |
double air_mass() |
| 1498 |
{ |
| 1499 |
double m; |
| 1500 |
if (sunzenith>90) |
| 1501 |
{ |
| 1502 |
if(suppress_warnings==0) |
| 1503 |
{ fprintf(stderr, "Warning: air mass has reached the maximal value\n"); } |
| 1504 |
sunzenith=90; |
| 1505 |
} |
| 1506 |
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) ); |
| 1507 |
return(m); |
| 1508 |
} |
| 1509 |
|
| 1510 |
|
| 1511 |
|
