ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gendaylit.c
Revision: 2.14
Committed: Fri Sep 6 16:54:06 2013 UTC (10 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.13: +181 -155 lines
Log Message:
Latest updates and bug fixesx from Wendelin

File Contents

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