Back to index

radiance  4R0+20100331
gendaylit.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: gendaylit.c,v 2.3 2009/06/20 21:34:34 greg Exp $";
00003 #endif
00004 /*        Copyright (c) 1994       *Fraunhofer Institut for Solar Energy Systems
00005  *                          Oltmannstr 5, D-79100 Freiburg, Germany
00006  *                          *Agence de l'Environnement et de la Maitrise de l'Energie
00007  *                          Centre de Valbonne, 500 route des Lucioles, 06565 Sophia Antipolis Cedex, France
00008  *                          *BOUYGUES 
00009  *                          1 Avenue Eugene Freyssinet, Saint-Quentin-Yvelines, France
00010 */
00011 
00012 
00013 
00014 /*
00015  *  gendaylit.c      program to generate the angular distribution of the daylight.
00016  *                   Our zenith is along the Z-axis, the X-axis
00017  *                   points east, and the Y-axis points north.
00018 */
00019 
00020 #include  <stdio.h>
00021 #include  <string.h>
00022 #include  <math.h>
00023 #include  <stdlib.h>
00024 #include  <ctype.h>
00025 
00026 #include  "rtio.h"
00027 #include  "fvect.h"
00028 #include  "color.h"
00029 #include  "paths.h"
00030 
00031 extern int jdate(int month, int day);
00032 extern double stadj(int  jd);
00033 extern double sdec(int  jd);
00034 extern double salt(double sd, double st);
00035 extern double sazi(double sd, double st);
00036 
00037 double  normsc();
00038 
00039 #define       DATFILE              "coeff_perez.dat"
00040 
00041 
00042 
00043 /* Perez sky parametrization : epsilon and delta calculations from the direct and diffuse irradiances */
00044 double sky_brightness();
00045 double sky_clearness();
00046 
00047 /* calculation of the direct and diffuse components from the Perez parametrization */
00048 double diffus_irradiance_from_sky_brightness();
00049 double        direct_irradiance_from_sky_clearness();
00050 
00051 
00052 /* Perez global horizontal, diffuse horizontal and direct normal luminous efficacy models : input w(cm)=2cm, solar zenith angle(degrees); output efficacy(lm/W) */
00053 double        glob_h_effi_PEREZ();
00054 double        glob_h_diffuse_effi_PEREZ();
00055 double        direct_n_effi_PEREZ();
00056 /*likelihood check of the epsilon, delta, direct and diffuse components*/
00057 void   check_parametrization();
00058 void   check_irradiances();
00059 void   check_illuminances();
00060 void   illu_to_irra_index();
00061 
00062 
00063 /* Perez sky luminance model */
00064 int    lect_coeff_perez(char *filename,float **coeff_perez);
00065 double  calc_rel_lum_perez(double dzeta,double gamma,double Z,
00066               double epsilon,double Delta,float *coeff_perez);
00067 /* coefficients for the sky luminance perez model */
00068 void   coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez);
00069 double radians(double degres);
00070 double degres(double radians);
00071 void   theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);
00072 double        integ_lv(float *lv,float *theta);
00073 float  *theta_ordered(char *filename);
00074 float  *phi_ordered(char *filename);
00075 void   skip_comments(FILE *fp);
00076 
00077 
00078 
00079 /* astronomy and geometry*/
00080 double        get_eccentricity();
00081 double        air_mass();
00082 double        get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut);
00083 
00084 
00085 /* date*/
00086 int     jdate(int month, int day);
00087 
00088 
00089 
00090 
00091 
00092 /* sun calculation constants */
00093 extern double  s_latitude;
00094 extern double  s_longitude;
00095 extern double  s_meridian;
00096 
00097 const double  AU = 149597890E3;
00098 const double  solar_constant_e = 1367;    /* solar constant W/m^2 */
00099 const double         solar_constant_l = 127.5;   /* solar constant klux */
00100 
00101 const double  half_sun_angle = 0.2665;
00102 const double  half_direct_angle = 2.85;
00103 
00104 const double  skyclearinf = 1.000; /* limitations for the variation of the Perez parameters */
00105 const double  skyclearsup = 12.1;
00106 const double  skybriginf = 0.01;
00107 const double  skybrigsup = 0.6;
00108 
00109 
00110 
00111 /* required values */
00112 int    month, day;                        /* date */
00113 double  hour;                             /* time */
00114 int    tsolar;                                   /* 0=standard, 1=solar */
00115 double  altitude, azimuth;                /* or solar angles */
00116 
00117 
00118 
00119 /* definition of the sky conditions through the Perez parametrization */
00120 double skyclearness, skybrightness;
00121 double solarradiance;       /*radiance of the sun disk and of the circumsolar area*/
00122 double diffusilluminance, directilluminance, diffusirradiance, directirradiance;
00123 double sunzenith, daynumber=150, atm_preci_water=2;
00124 
00125 double        diffnormalization, dirnormalization;
00126 double        *c_perez;
00127 
00128 int    output=0;     /*define the unit of the output (sky luminance or radiance): visible watt=0, solar watt=1, lumen=2*/
00129 int    input=0;      /*define the input for the calulation*/
00130 
00131        /* default values */
00132 int  cloudy = 0;                          /* 1=standard, 2=uniform */
00133 int  dosun = 1;
00134 double  zenithbr = -1.0;
00135 double  betaturbidity = 0.1;
00136 double  gprefl = 0.2;
00137 int    S_INTER=0;
00138 
00139        /* computed values */
00140 double  sundir[3];
00141 double  groundbr;
00142 double  F2;
00143 double  solarbr = 0.0;
00144 int    u_solar = 0;                       /* -1=irradiance, 1=radiance */
00145 
00146 char  *progname;
00147 char  errmsg[128];
00148 
00149 
00150 main(argc, argv)
00151 int  argc;
00152 char  *argv[];
00153 {
00154        int  i;
00155 
00156        progname = argv[0];
00157        if (argc == 2 && !strcmp(argv[1], "-defaults")) {
00158               printdefaults();
00159               exit(0);
00160        }
00161        if (argc < 4)
00162               userror("arg count");
00163        if (!strcmp(argv[1], "-ang")) {
00164               altitude = atof(argv[2]) * (M_PI/180);
00165               azimuth = atof(argv[3]) * (M_PI/180);
00166               month = 0;
00167        } else {
00168               month = atoi(argv[1]);
00169               if (month < 1 || month > 12)
00170                      userror("bad month");
00171               day = atoi(argv[2]);
00172               if (day < 1 || day > 31)
00173                      userror("bad day");
00174               hour = atof(argv[3]);
00175               if (hour < 0 || hour >= 24)
00176                      userror("bad hour");
00177               tsolar = argv[3][0] == '+';
00178        }
00179        for (i = 4; i < argc; i++)
00180               if (argv[i][0] == '-' || argv[i][0] == '+')
00181                      switch (argv[i][1]) {
00182                      case 's':
00183                             cloudy = 0;
00184                             dosun = argv[i][0] == '+';
00185                             break;
00186                      case 'r':
00187                      case 'R':
00188                             u_solar = argv[i][1] == 'R' ? -1 : 1;
00189                             solarbr = atof(argv[++i]);
00190                             break;
00191                      case 'c':
00192                             cloudy = argv[i][0] == '+' ? 2 : 1;
00193                             dosun = 0;
00194                             break;
00195                      case 't':
00196                             betaturbidity = atof(argv[++i]);
00197                             break;
00198                      case 'b':
00199                             zenithbr = atof(argv[++i]);
00200                             break;
00201                      case 'g':
00202                             gprefl = atof(argv[++i]);
00203                             break;
00204                      case 'a':
00205                             s_latitude = atof(argv[++i]) * (M_PI/180);
00206                             break;
00207                      case 'o':
00208                             s_longitude = atof(argv[++i]) * (M_PI/180);
00209                             break;
00210                      case 'm':
00211                             s_meridian = atof(argv[++i]) * (M_PI/180);
00212                             break;
00213 
00214                      
00215                      case 'O':
00216                             output = atof(argv[++i]);   /*define the unit of the output of the program : 
00217                                                         sky and sun luminance/radiance (0==W visible, 1==W solar radiation, 2==lm) 
00218                                                         default is set to 0*/
00219                             break;
00220                             
00221                      case 'P':
00222                             input = 0;                         /* Perez parameters: epsilon, delta */
00223                             skyclearness = atof(argv[++i]);
00224                             skybrightness = atof(argv[++i]);
00225                             break;
00226 
00227                      case 'W':                                 /* direct normal Irradiance [W/m^2] */
00228                             input = 1;                         /* diffuse horizontal Irrad. [W/m^2] */
00229                             directirradiance = atof(argv[++i]);
00230                             diffusirradiance = atof(argv[++i]);
00231                             break;
00232                             
00233                      case 'L':                                 /* direct normal Illuminance [Lux] */
00234                             input = 2;                         /* diffuse horizontal Ill. [Lux] */
00235                             directilluminance = atof(argv[++i]);
00236                             diffusilluminance = atof(argv[++i]);
00237                             break;
00238                      
00239                      case 'G':                                 /* direct horizontal Irradiance [W/m^2] */
00240                             input = 3;                         /* diffuse horizontal Irrad. [W/m^2] */
00241                             directirradiance = atof(argv[++i]);
00242                             diffusirradiance = atof(argv[++i]);
00243                             break;
00244                             
00245                      
00246                      default:
00247                             sprintf(errmsg, "unknown option: %s", argv[i]);
00248                             userror(errmsg);
00249                      }
00250               else
00251                      userror("bad option");
00252 
00253        if (fabs(s_meridian-s_longitude) > 30*M_PI/180)
00254               fprintf(stderr,
00255                   "%s: warning: %.1f hours btwn. standard meridian and longitude\n",
00256                   progname, (s_longitude-s_meridian)*12/M_PI);
00257 
00258 
00259        /* allocation dynamique de memoire pour les pointeurs */
00260        if ( (c_perez = malloc(5*sizeof(double))) == NULL )
00261        {
00262               fprintf(stderr,"Out of memory error in function main !");
00263               exit(1);
00264        }
00265 
00266 
00267        printhead(argc, argv);
00268 
00269        computesky();
00270        printsky();
00271 
00272        exit(0);
00273 }
00274 
00275 
00276 computesky()                /* compute sky parameters */
00277 {
00278 
00279        /* new variables */
00280        int    j, i;
00281        float  *lv_mod;  /* 145 luminance values*/
00282          /* 145 directions for the calculation of the normalization coefficient, coefficient Perez model */
00283        float  *theta_o, *phi_o, *coeff_perez;
00284        double        dzeta, gamma;
00285        double diffusion;
00286        double normfactor;
00287 
00288 
00289 
00290        /* compute solar direction */
00291 
00292        if (month) {                /* from date and time */
00293               int  jd;
00294               double  sd, st;
00295 
00296               jd = jdate(month, day);            /* Julian date */
00297               sd = sdec(jd);                     /* solar declination */
00298               if (tsolar)                 /* solar time */
00299                      st = hour;
00300               else
00301                      st = hour + stadj(jd);
00302               altitude = salt(sd, st);
00303               azimuth = sazi(sd, st);
00304               
00305               daynumber = (double)jdate(month, day);
00306 
00307        }
00308        if (!cloudy && altitude > 87.*M_PI/180.) {
00309               fprintf(stderr,
00310                   "%s: warning - sun too close to zenith, reducing altitude to 87 degrees\n",
00311                   progname);
00312               printf(
00313                   "# warning - sun too close to zenith, reducing altitude to 87 degrees\n");
00314               altitude = 87.*M_PI/180.;
00315        }
00316        sundir[0] = -sin(azimuth)*cos(altitude);
00317        sundir[1] = -cos(azimuth)*cos(altitude);
00318        sundir[2] = sin(altitude);
00319        
00320               
00321        /* calculation for the new functions */
00322        sunzenith = 90 - altitude*180/M_PI;
00323        
00324        
00325 
00326 /* compute the inputs for the calculation of the light distribution over the sky*/
00327        if (input==0)
00328               {
00329               check_parametrization();
00330               diffusirradiance = diffus_irradiance_from_sky_brightness(); /*diffuse horizontal irradiance*/
00331               directirradiance = direct_irradiance_from_sky_clearness();
00332               check_irradiances();
00333               
00334               if (output==0 || output==2)
00335                      {
00336                      diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
00337                      directilluminance = directirradiance*direct_n_effi_PEREZ();
00338                      check_illuminances();
00339                      }
00340               }
00341        
00342 
00343        else if (input==1)
00344               {
00345               check_irradiances();
00346               skybrightness = sky_brightness();
00347               skyclearness =  sky_clearness();
00348               check_parametrization();
00349 
00350               if (output==0 || output==2)
00351                      {
00352                      diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
00353                      directilluminance = directirradiance*direct_n_effi_PEREZ();
00354                      check_illuminances();
00355                      }
00356 
00357               }
00358                      
00359        
00360        else if (input==2)
00361               {             
00362               check_illuminances();
00363               illu_to_irra_index();
00364               check_parametrization();
00365               }
00366               
00367 
00368        else if (input==3)
00369               {
00370                      if (altitude<=0)
00371                      {
00372                      fprintf(stderr, "solar zenith angle larger than 90 \n the models used are not more valid\n");
00373                      exit(1);
00374                      }
00375 
00376               directirradiance=directirradiance/sin(altitude);
00377               check_irradiances();
00378               skybrightness = sky_brightness();
00379               skyclearness =  sky_clearness();
00380               check_parametrization();
00381 
00382               if (output==0 || output==2)
00383                      {
00384                      diffusilluminance = diffusirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
00385                      directilluminance = directirradiance*direct_n_effi_PEREZ();
00386                      check_illuminances();
00387                      }
00388 
00389               }
00390 
00391        
00392        else   {fprintf(stderr,"error in giving the input arguments"); exit(1);}
00393 
00394 
00395        
00396 /* normalization factor for the relative sky luminance distribution, diffuse part*/
00397 
00398        /* allocation dynamique de memoire pour les pointeurs */
00399        if ( (coeff_perez = malloc(8*20*sizeof(float))) == NULL )
00400        {
00401               fprintf(stderr,"Out of memory error in function main !");
00402               exit(1);
00403        }
00404 
00405 /* read the coefficients for the Perez sky luminance model */
00406        if (lect_coeff_perez(DATFILE, &coeff_perez) > 0)
00407        {
00408               fprintf(stderr,"lect_coeff_perez does not work\n");
00409               exit(2);
00410        }
00411 
00412        if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
00413        {
00414               fprintf(stderr,"Out of memory in function main");
00415               exit(1);
00416        }
00417 
00418        /* read the angles */
00419        theta_o = theta_ordered("defangle.dat");
00420        phi_o = phi_ordered("defangle.dat");
00421 
00422 /* parameters for the perez model */
00423        coeff_lum_perez(radians(sunzenith), skyclearness, skybrightness, coeff_perez);
00424 
00425 /*calculation of the modelled luminance */
00426        for (j=0;j<145;j++)
00427        {
00428               theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
00429               *(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
00430               /*printf("theta, phi, lv_mod %lf\t %lf\t %lf\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));*/
00431        }
00432 
00433        /* integration of luminance for the normalization factor, diffuse part of the sky*/
00434        diffnormalization = integ_lv(lv_mod, theta_o);
00435        /*printf("perez integration %lf\n", diffnormalization);*/
00436        
00437        
00438 
00439 
00440 /*normalization coefficient in lumen or in watt*/
00441        if (output==0)
00442               {
00443               diffnormalization = diffusilluminance/diffnormalization/WHTEFFICACY;
00444               }
00445        else if (output==1)
00446               {
00447               diffnormalization = diffusirradiance/diffnormalization;
00448               }
00449        else if (output==2)
00450               {
00451               diffnormalization = diffusilluminance/diffnormalization;
00452               }
00453 
00454        else   {fprintf(stderr,"output argument : wrong number"); exit(1);}
00455 
00456 
00457 
00458 
00459 /* calculation for the solar source */    
00460        if (output==0) 
00461               solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
00462               
00463        else if (output==1)
00464               solarradiance = directirradiance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
00465        
00466        else
00467               solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)));
00468        
00469                      
00470 
00471 
00472 /* Compute the ground radiance */
00473 zenithbr=calc_rel_lum_perez(0.0,radians(sunzenith),radians(sunzenith),skyclearness,skybrightness,coeff_perez);
00474 zenithbr*=diffnormalization;
00475 /*
00476 fprintf(stderr, "gendaylit : the actual zenith radiance(W/m^2/sr) or luminance(cd/m^2) is : %.0lf\n", zenithbr);
00477 */
00478 
00479 if (skyclearness==1)
00480        normfactor = 0.777778;
00481               
00482 if (skyclearness>=6)
00483        {             
00484        F2 = 0.274*(0.91 + 10.0*exp(-3.0*(M_PI/2.0-altitude)) + 0.45*sundir[2]*sundir[2]);
00485        normfactor = normsc()/F2/M_PI;
00486        }
00487 
00488 if ( (skyclearness>1) && (skyclearness<6) )
00489        {
00490        S_INTER=1;
00491        F2 = (2.739 + .9891*sin(.3119+2.6*altitude)) * exp(-(M_PI/2.0-altitude)*(.4441+1.48*altitude));
00492        normfactor = normsc()/F2/M_PI;
00493        }
00494 
00495 groundbr = zenithbr*normfactor;
00496 printf("# Ground ambient level: %.1f\n", groundbr);
00497 
00498 if (dosun&&(skyclearness>1)) 
00499 groundbr += 6.8e-5/M_PI*solarradiance*sundir[2];        
00500 
00501 groundbr *= gprefl;
00502 
00503 
00504 
00505 return;
00506 }
00507 
00508 
00509 
00510 
00511 
00512 
00513 
00514 printsky()                  /* print out sky */
00515 {
00516        if (dosun&&(skyclearness>1)) 
00517               {             
00518               printf("\nvoid light solar\n");
00519               printf("0\n0\n");
00520               printf("3 %.3e %.3e %.3e\n", solarradiance, solarradiance, solarradiance);
00521               printf("\nsolar source sun\n");
00522               printf("0\n0\n");
00523               printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
00524               }
00525               
00526        if (dosun&&(skyclearness==1)) 
00527               {             
00528               printf("\nvoid light solar\n");
00529               printf("0\n0\n");
00530               printf("3 0.0 0.0 0.0\n");
00531               printf("\nsolar source sun\n");
00532               printf("0\n0\n");
00533               printf("4 %f %f %f %f\n", sundir[0], sundir[1], sundir[2], 2*half_sun_angle);
00534               }
00535        
00536 
00537        printf("\nvoid brightfunc skyfunc\n");
00538        printf("2 skybright perezlum.cal\n");
00539        printf("0\n");
00540        printf("10 %.3e %.3e %lf %lf %lf %lf %lf %f %f %f \n", diffnormalization, groundbr,
00541            *(c_perez+0),*(c_perez+1),*(c_perez+2),*(c_perez+3),*(c_perez+4), 
00542            sundir[0], sundir[1], sundir[2]);
00543 }
00544 
00545 
00546 printdefaults()                    /* print default values */
00547 {
00548        printf("-g %f\t\t\t# Ground plane reflectance\n", gprefl);
00549        if (zenithbr > 0.0)
00550               printf("-b %f\t\t\t# Zenith radiance (watts/ster/m^2\n", zenithbr);
00551        else
00552               printf("-t %f\t\t\t# Atmospheric betaturbidity\n", betaturbidity);
00553        printf("-a %f\t\t\t# Site latitude (degrees)\n", s_latitude*(180/M_PI));
00554        printf("-o %f\t\t\t# Site longitude (degrees)\n", s_longitude*(180/M_PI));
00555        printf("-m %f\t\t\t# Standard meridian (degrees)\n", s_meridian*(180/M_PI));
00556 }
00557 
00558 
00559 userror(msg)                /* print usage error and quit */
00560 char  *msg;
00561 {
00562        if (msg != NULL)
00563               fprintf(stderr, "%s: Use error - %s\n", progname, msg);
00564        fprintf(stderr, "Usage: %s month day hour [-P|-W|-L] direct_value diffus_value [options]\n", progname);
00565        fprintf(stderr, "or   : %s -ang altitude azimuth [-P|-W|-L] direct_value diffus_value [options]\n", progname);
00566        fprintf(stderr, "    -P epsilon delta  (these are the Perez parameters) \n");
00567        fprintf(stderr, "    -W direct-normal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
00568        fprintf(stderr, "    -L direct-normal-illuminance diffuse-horizontal-illuminance (lux)\n");
00569        fprintf(stderr, "    -G direct-horizontal-irradiance diffuse-horizontal-irradiance (W/m^2)\n");
00570        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");
00571        exit(1);
00572 }
00573 
00574 
00575 
00576 double
00577 normsc()                    /* compute normalization factor (E0*F2/L0) */
00578 {
00579        static double  nfc[2][5] = {
00580                             /* clear sky approx. */
00581               {2.766521, 0.547665, -0.369832, 0.009237, 0.059229},
00582                             /* intermediate sky approx. */
00583               {3.5556, -2.7152, -1.3081, 1.0660, 0.60227},
00584        };
00585        register double  *nf;
00586        double  x, nsc;
00587        register int  i;
00588                                    /* polynomial approximation */
00589        nf = nfc[S_INTER];
00590        x = (altitude - M_PI/4.0)/(M_PI/4.0);
00591        nsc = nf[i=4];
00592        while (i--)
00593               nsc = nsc*x + nf[i];
00594 
00595        return(nsc);
00596 }
00597 
00598 
00599 
00600 printhead(ac, av)           /* print command header */
00601 register int  ac;
00602 register char  **av;
00603 {
00604        putchar('#');
00605        while (ac--) {
00606               putchar(' ');
00607               fputs(*av++, stdout);
00608        }
00609        putchar('\n');
00610 }
00611 
00612 
00613 
00614 
00615 void
00616 skip_comments(FILE *fp)            /* skip comments in file */
00617 {
00618        int    c;
00619        
00620        while ((c = getc(fp)) != EOF)
00621               if (c == '#') {
00622                      while ((c = getc(fp)) != EOF)
00623                             if (c == '\n')
00624                                    break;
00625               } else if (!isspace(c)) {
00626                      ungetc(c, fp);
00627                      break;
00628               }
00629 }
00630 
00631 
00632 
00633 /* Perez models */
00634 
00635 /* Perez global horizontal luminous efficacy model */
00636 double glob_h_effi_PEREZ()
00637 {
00638 
00639        double        value;
00640        double        category_bounds[10], a[10], b[10], c[10], d[10];
00641        int    category_total_number, category_number, i;
00642 
00643 
00644 if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
00645 fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
00646 
00647        /* initialize category bounds (clearness index bounds) */
00648 
00649        category_total_number = 8;
00650 
00651        category_bounds[1] = 1;
00652        category_bounds[2] = 1.065;
00653        category_bounds[3] = 1.230;
00654        category_bounds[4] = 1.500;
00655        category_bounds[5] = 1.950;
00656        category_bounds[6] = 2.800;
00657        category_bounds[7] = 4.500;
00658        category_bounds[8] = 6.200;
00659        category_bounds[9] = 12.01;
00660 
00661 
00662        /* initialize model coefficients */
00663        a[1] = 96.63;
00664        a[2] = 107.54;
00665        a[3] = 98.73;
00666        a[4] = 92.72;
00667        a[5] = 86.73;
00668        a[6] = 88.34;
00669        a[7] = 78.63;
00670        a[8] = 99.65;
00671 
00672        b[1] = -0.47;
00673        b[2] = 0.79;
00674        b[3] = 0.70;
00675        b[4] = 0.56;
00676        b[5] = 0.98;
00677        b[6] = 1.39;
00678        b[7] = 1.47;
00679        b[8] = 1.86;
00680 
00681        c[1] = 11.50;
00682        c[2] = 1.79;
00683        c[3] = 4.40;
00684        c[4] = 8.36;
00685        c[5] = 7.10;
00686        c[6] = 6.06;
00687        c[7] = 4.93;
00688        c[8] = -4.46;
00689 
00690        d[1] = -9.16;
00691        d[2] = -1.19;
00692        d[3] = -6.95;
00693        d[4] = -8.31;
00694        d[5] = -10.94;
00695        d[6] = -7.60;
00696        d[7] = -11.37;
00697        d[8] = -3.15;
00698 
00699 
00700 
00701 
00702        for (i=1; i<=category_total_number; i++)
00703        {
00704               if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
00705                      category_number = i;
00706        }
00707 
00708        value = a[category_number] + b[category_number]*atm_preci_water  + 
00709            c[category_number]*cos(sunzenith*M_PI/180) +  d[category_number]*log(skybrightness);
00710 
00711        return(value);
00712 }
00713 
00714 
00715 /* global horizontal diffuse efficacy model, according to PEREZ */
00716 double glob_h_diffuse_effi_PEREZ()
00717 {
00718        double        value;
00719        double        category_bounds[10], a[10], b[10], c[10], d[10];
00720        int    category_total_number, category_number, i;
00721 
00722        
00723 if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
00724 fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
00725 
00726 /* initialize category bounds (clearness index bounds) */
00727 
00728        category_total_number = 8;
00729 
00730        category_bounds[1] = 1;
00731        category_bounds[2] = 1.065;
00732        category_bounds[3] = 1.230;
00733        category_bounds[4] = 1.500;
00734        category_bounds[5] = 1.950;
00735        category_bounds[6] = 2.800;
00736        category_bounds[7] = 4.500;
00737        category_bounds[8] = 6.200;
00738        category_bounds[9] = 12.01;
00739 
00740 
00741        /* initialize model coefficients */
00742        a[1] = 97.24;
00743        a[2] = 107.22;
00744        a[3] = 104.97;
00745        a[4] = 102.39;
00746        a[5] = 100.71;
00747        a[6] = 106.42;
00748        a[7] = 141.88;
00749        a[8] = 152.23;
00750 
00751        b[1] = -0.46;
00752        b[2] = 1.15;
00753        b[3] = 2.96;
00754        b[4] = 5.59;
00755        b[5] = 5.94;
00756        b[6] = 3.83;
00757        b[7] = 1.90;
00758        b[8] = 0.35;
00759 
00760        c[1] = 12.00;
00761        c[2] = 0.59;
00762        c[3] = -5.53;
00763        c[4] = -13.95;
00764        c[5] = -22.75;
00765        c[6] = -36.15;
00766        c[7] = -53.24;
00767        c[8] = -45.27;
00768 
00769        d[1] = -8.91;
00770        d[2] = -3.95;
00771        d[3] = -8.77;
00772        d[4] = -13.90;
00773        d[5] = -23.74;
00774        d[6] = -28.83;
00775        d[7] = -14.03;
00776        d[8] = -7.98;
00777 
00778 
00779 
00780 
00781        for (i=1; i<=category_total_number; i++)
00782        {
00783               if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
00784                      category_number = i;
00785        }
00786 
00787        value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) + 
00788            d[category_number]*log(skybrightness);
00789 
00790        return(value);
00791 }
00792 
00793 
00794 /* direct normal efficacy model, according to PEREZ */
00795 
00796 double direct_n_effi_PEREZ()
00797 
00798 {
00799 double        value;
00800 double        category_bounds[10], a[10], b[10], c[10], d[10];
00801 int    category_total_number, category_number, i;
00802 
00803 
00804 if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
00805 fprintf(stderr, "Warning : skyclearness or skybrightness out of range ; \n Check your input parameters\n");
00806 
00807 
00808 /* initialize category bounds (clearness index bounds) */
00809 
00810 category_total_number = 8;
00811 
00812 category_bounds[1] = 1;
00813 category_bounds[2] = 1.065;
00814 category_bounds[3] = 1.230;
00815 category_bounds[4] = 1.500;
00816 category_bounds[5] = 1.950;
00817 category_bounds[6] = 2.800;
00818 category_bounds[7] = 4.500;
00819 category_bounds[8] = 6.200;
00820 category_bounds[9] = 12.1;
00821 
00822 
00823 /* initialize model coefficients */
00824 a[1] = 57.20;
00825 a[2] = 98.99;
00826 a[3] = 109.83;
00827 a[4] = 110.34;
00828 a[5] = 106.36;
00829 a[6] = 107.19;
00830 a[7] = 105.75;
00831 a[8] = 101.18;
00832 
00833 b[1] = -4.55;
00834 b[2] = -3.46;
00835 b[3] = -4.90;
00836 b[4] = -5.84;
00837 b[5] = -3.97;
00838 b[6] = -1.25;
00839 b[7] = 0.77;
00840 b[8] = 1.58;
00841 
00842 c[1] = -2.98;
00843 c[2] = -1.21;
00844 c[3] = -1.71;
00845 c[4] = -1.99;
00846 c[5] = -1.75;
00847 c[6] = -1.51;
00848 c[7] = -1.26;
00849 c[8] = -1.10;
00850 
00851 d[1] = 117.12;
00852 d[2] = 12.38;
00853 d[3] = -8.81;
00854 d[4] = -4.56;
00855 d[5] = -6.16;
00856 d[6] = -26.73;
00857 d[7] = -34.44;
00858 d[8] = -8.29;
00859 
00860 
00861 
00862 for (i=1; i<=category_total_number; i++)
00863 {
00864 if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
00865 category_number = i;
00866 }
00867 
00868 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;
00869 
00870 if (value < 0) value = 0;
00871 
00872 return(value);
00873 }
00874 
00875 
00876 /*check the range of epsilon and delta indexes of the perez parametrization*/
00877 void check_parametrization()
00878 {
00879 if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<=skybriginf || skybrightness>skybrigsup)
00880               {
00881               fprintf(stderr,"sky clearness or sky brightness out of range %lf\t %lf\n", skyclearness, skybrightness);
00882               exit(1);      
00883               }
00884        else return;
00885 }
00886 
00887 
00888 /* likelihood of the direct and diffuse components */
00889 void   check_illuminances()
00890 {
00891        if (!( (directilluminance>=0) && (directilluminance<=solar_constant_l*1000) && (diffusilluminance>0) ))
00892        {
00893        fprintf(stderr,"direct or diffuse illuminances out of range\n");
00894        exit(1);
00895        }
00896 return;
00897 }
00898 
00899 
00900 void   check_irradiances()
00901 {
00902        if (!( (directirradiance>=0) && (directirradiance<=solar_constant_e) && (diffusirradiance>0) ))
00903        {
00904        fprintf(stderr,"direct or diffuse irradiances out of range\n");
00905        exit(1);
00906        }      
00907 return;
00908 }
00909        
00910 
00911 
00912 /* Perez sky's brightness */
00913 double sky_brightness()
00914 {
00915 double value;
00916 
00917 value = diffusirradiance * air_mass() / ( solar_constant_e*get_eccentricity());
00918 
00919 return(value);
00920 }
00921 
00922 
00923 /* Perez sky's clearness */
00924 double sky_clearness()
00925 {
00926 double value;
00927 
00928 value = ( (diffusirradiance + directirradiance)/(diffusirradiance) + 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) ;
00929 
00930 return(value);
00931 }
00932 
00933 
00934 
00935 /* diffus horizontal irradiance from Perez sky's brightness */
00936 double diffus_irradiance_from_sky_brightness()
00937 {
00938        double value;
00939 
00940        value = skybrightness / air_mass() * ( solar_constant_e*get_eccentricity());
00941 
00942        return(value);
00943 }
00944 
00945 
00946 /* direct normal irradiance from Perez sky's clearness */
00947 double direct_irradiance_from_sky_clearness()
00948 {
00949        double value;
00950 
00951        value = diffus_irradiance_from_sky_brightness();
00952        value = value * ( (skyclearness-1) * (1+1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) );
00953 
00954        return(value);
00955 }
00956 
00957 
00958 void illu_to_irra_index(void)
00959 {
00960 double test1=0.1, test2=0.1;
00961 int    counter=0;    
00962 
00963 diffusirradiance = diffusilluminance*solar_constant_e/(solar_constant_l*1000);
00964 directirradiance = directilluminance*solar_constant_e/(solar_constant_l*1000);
00965 skyclearness =  sky_clearness(); 
00966 skybrightness = sky_brightness();
00967 if (skyclearness>12) skyclearness=12;
00968 if (skybrightness<0.05) skybrightness=0.01; 
00969        
00970        
00971 while ( ((fabs(diffusirradiance-test1)>10) || (fabs(directirradiance-test2)>10)
00972               || skyclearness>skyclearinf || skyclearness<skyclearsup 
00973               || skybrightness>skybriginf || skybrightness<skybrigsup )
00974                && !(counter==5) )
00975        {
00976               /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
00977 
00978        test1=diffusirradiance;
00979        test2=directirradiance;     
00980        counter++;
00981        
00982        diffusirradiance = diffusilluminance/glob_h_diffuse_effi_PEREZ();
00983        directirradiance = directilluminance/direct_n_effi_PEREZ();
00984        /*fprintf(stderr, "conversion illuminance into irradiance %lf\t %lf\n", diffusirradiance, directirradiance);*/
00985        
00986        skybrightness = sky_brightness();
00987        skyclearness =  sky_clearness();
00988        if (skyclearness>12) skyclearness=12;
00989        if (skybrightness<0.05) skybrightness=0.01; 
00990 
00991               /*fprintf(stderr, "%lf\t %lf\n", skybrightness, skyclearness);*/
00992 
00993        }
00994 
00995 
00996 return;
00997 }             
00998 
00999 
01000 int lect_coeff_perez(char *filename,float **coeff_perez)
01001 {
01002        FILE *fcoeff_perez;
01003        float temp;
01004        int i,j;
01005 
01006        if ((fcoeff_perez = frlibopen(filename)) == NULL)
01007        {
01008               fprintf(stderr,"file %s cannot be opened\n", filename);
01009               return 1; /* il y a un probleme de fichier */
01010        }
01011        else
01012        {
01013               /*printf("file %s  open\n", filename);*/
01014        }
01015        
01016        skip_comments(fcoeff_perez);
01017 
01018        for (i=0;i<8;i++)
01019               for (j=0;j<20;j++)
01020               {
01021                      fscanf(fcoeff_perez,"%f",&temp);
01022                      *(*coeff_perez+i*20+j) = temp;
01023               }
01024        fclose(fcoeff_perez);
01025 
01026        return 0; /* tout est OK */
01027 }
01028 
01029 
01030 
01031 /* sky luminance perez model */
01032 double calc_rel_lum_perez(double dzeta,double gamma,double Z,
01033 double epsilon,double Delta,float *coeff_perez)
01034 {
01035        float x[5][4];
01036        int i,j,num_lin;
01037        double c_perez[5];
01038 
01039        if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
01040        {
01041               fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
01042               exit(1);
01043        }
01044 
01045        /* correction de modele de Perez solar energy ...*/
01046        if ( (epsilon > 1.065) && (epsilon < 2.8) )
01047        {
01048               if ( Delta < 0.2 ) Delta = 0.2;
01049        }
01050 
01051        if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
01052        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
01053        if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
01054        if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
01055        if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
01056        if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
01057        if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
01058        if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
01059 
01060        for (i=0;i<5;i++)
01061               for (j=0;j<4;j++)
01062               {
01063                      x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
01064                      /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
01065               }
01066 
01067 
01068        if (num_lin)
01069        {
01070               for (i=0;i<5;i++)
01071                      c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
01072        }
01073        else
01074        {
01075               c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
01076               c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
01077               c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
01078               c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
01079               c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
01080        }
01081 
01082 
01083        return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) *
01084            (1 + c_perez[2]*exp(c_perez[3]*gamma) +
01085            c_perez[4]*cos(gamma)*cos(gamma) );
01086 }
01087 
01088 
01089 
01090 /* coefficients for the sky luminance perez model */
01091 void coeff_lum_perez(double Z, double epsilon, double Delta, float *coeff_perez)
01092 {
01093        float x[5][4];
01094        int i,j,num_lin;
01095 
01096        if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
01097        {
01098               fprintf(stderr,"Epsilon out of range in function calc_rel_lum_perez !\n");
01099               exit(1);
01100        }
01101 
01102        /* correction du modele de Perez solar energy ...*/
01103        if ( (epsilon > 1.065) && (epsilon < 2.8) )
01104        {
01105               if ( Delta < 0.2 ) Delta = 0.2;
01106        }
01107 
01108        if ( (epsilon >= 1.000) && (epsilon < 1.065) ) num_lin = 0;
01109        if ( (epsilon >= 1.065) && (epsilon < 1.230) ) num_lin = 1;
01110        if ( (epsilon >= 1.230) && (epsilon < 1.500) ) num_lin = 2;
01111        if ( (epsilon >= 1.500) && (epsilon < 1.950) ) num_lin = 3;
01112        if ( (epsilon >= 1.950) && (epsilon < 2.800) ) num_lin = 4;
01113        if ( (epsilon >= 2.800) && (epsilon < 4.500) ) num_lin = 5;
01114        if ( (epsilon >= 4.500) && (epsilon < 6.200) ) num_lin = 6;
01115        if ( (epsilon >= 6.200) && (epsilon < 14.00) ) num_lin = 7;
01116 
01117        for (i=0;i<5;i++)
01118               for (j=0;j<4;j++)
01119               {
01120                      x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
01121                      /* printf("x %d %d vaut %f\n",i,j,x[i][j]); */
01122               }
01123 
01124 
01125        if (num_lin)
01126        {
01127               for (i=0;i<5;i++)
01128                      *(c_perez+i) = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
01129 
01130        }
01131        else
01132        {
01133               *(c_perez+0) = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
01134               *(c_perez+1) = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
01135               *(c_perez+4) = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
01136               *(c_perez+2) = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
01137               *(c_perez+3) = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
01138 
01139 
01140        }
01141 
01142 
01143        return;
01144 }
01145 
01146 
01147 /* degrees into radians */
01148 double radians(double degres)
01149 {
01150        return degres*M_PI/180.0;
01151 }
01152 
01153 /* radian into degrees */
01154 double degres(double radians)
01155 {
01156        return radians/M_PI*180.0;
01157 }
01158 
01159 /* calculation of the angles dzeta and gamma */
01160 void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
01161 {
01162        *dzeta = theta; /* dzeta = phi */
01163        if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) )
01164               *gamma = 0;
01165        else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
01166        {
01167               printf("error in calculation of gamma (angle between point and sun");
01168               exit(3);
01169        }
01170        else
01171               *gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
01172 }
01173 
01174 
01175 /********************************************************************************/
01176 /*     Fonction: theta_ordered                                               */
01177 /*                                                                    */
01178 /*     In: char *filename                                             */
01179 /*                                                                    */
01180 /*     Out: float *                                                   */
01181 /*                                                                    */
01182 /*     Update: 29/08/93                                               */
01183 /*                                                                    */
01184 /*     Rem: theta en degres                                           */
01185 /*                                                                    */
01186 /*     But: fournit les valeurs de theta du fichier d'entree a la memoire    */
01187 /*                                                                    */
01188 /********************************************************************************/
01189 float *theta_ordered(char *filename)
01190 {
01191        int i;
01192        float buffer,*ptr;
01193        FILE *file_in;
01194 
01195        if ( (file_in = frlibopen(filename)) == NULL )
01196        {
01197               fprintf(stderr,"Cannot open file %s in function theta_ordered\n",filename);
01198               exit(1);
01199        }
01200        
01201        skip_comments(file_in);
01202 
01203        if ( (ptr = malloc(145*sizeof(float))) == NULL )
01204        {
01205               fprintf(stderr,"Out of memory in function theta_ordered\n");
01206               exit(1);
01207        }
01208 
01209        for (i=0;i<145;i++)
01210        {
01211               fscanf(file_in,"%f",&buffer);
01212               *(ptr+i) = buffer;
01213               fscanf(file_in,"%f",&buffer);
01214        }
01215 
01216        fclose(file_in);
01217        return ptr;
01218 }
01219 
01220 
01221 /********************************************************************************/
01222 /*     Fonction: phi_ordered                                                 */
01223 /*                                                                    */
01224 /*     In: char *filename                                             */
01225 /*                                                                    */
01226 /*     Out: float *                                                   */
01227 /*                                                                    */
01228 /*     Update: 29/08/93                                               */
01229 /*                                                                    */
01230 /*     Rem: valeurs de Phi en DEGRES                                         */
01231 /*                                                                    */
01232 /*     But: mettre les angles contenus dans le fichier d'entree dans la memoire */
01233 /*                                                                    */
01234 /********************************************************************************/
01235 float *phi_ordered(char *filename)
01236 {
01237        int i;
01238        float buffer,*ptr;
01239        FILE *file_in;
01240 
01241        if ( (file_in = frlibopen(filename)) == NULL )
01242        {
01243               fprintf(stderr,"Cannot open file %s in function phi_ordered\n",filename);
01244               exit(1);
01245        }
01246        
01247        skip_comments(file_in);
01248 
01249        if ( (ptr = malloc(145*sizeof(float))) == NULL )
01250        {
01251               fprintf(stderr,"Out of memory in function phi_ordered");
01252               exit(1);
01253        }
01254 
01255        for (i=0;i<145;i++)
01256        {
01257               fscanf(file_in,"%f",&buffer);
01258               fscanf(file_in,"%f",&buffer);
01259               *(ptr+i) = buffer;
01260        }
01261 
01262        fclose(file_in);
01263        return ptr;
01264 }
01265 
01266 
01267 /********************************************************************************/
01268 /*     Fonction: integ_lv                                             */
01269 /*                                                                    */
01270 /*     In: float *lv,*theta                                           */
01271 /*         int sun_pos                                                       */
01272 /*                                                                    */
01273 /*     Out: double                                                    */
01274 /*                                                                    */
01275 /*     Update: 29/08/93                                               */
01276 /*                                                                    */
01277 /*     Rem:                                                           */
01278 /*                                                                    */
01279 /*     But: calcul l'integrale de luminance relative sans la dir. du soleil  */
01280 /*                                                                    */
01281 /********************************************************************************/
01282 double integ_lv(float *lv,float *theta)
01283 {
01284        int i;
01285        double buffer=0.0;
01286 
01287        for (i=0;i<145;i++)
01288               buffer += (*(lv+i))*cos(radians(*(theta+i)));
01289 
01290        return buffer*2*M_PI/144;
01291 
01292 }
01293 
01294 
01295 
01296 
01297 
01298 
01299 /* enter day number(double), return E0 = square(R0/R):  eccentricity correction factor  */
01300 
01301 double get_eccentricity()
01302 {
01303        double day_angle;
01304        double E0;
01305 
01306        day_angle  = 2*M_PI*(daynumber - 1)/365;
01307        E0         = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
01308            0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);
01309 
01310        return (E0);
01311 
01312 }
01313 
01314 
01315 /* enter sunzenith angle (degrees) return relative air mass (double) */
01316 double        air_mass()
01317 {
01318 double m;
01319 
01320 if (sunzenith>90)
01321        {
01322        fprintf(stderr, "solar zenith angle larger than 90 in fuction air_mass():\n the models used are not more valid\n");
01323        exit(1);
01324        }
01325        
01326 m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
01327 return(m);
01328 }
01329 
01330 
01331 double get_angle_sun_direction(double sun_zenith, double sun_azimut, double direction_zenith, double direction_azimut)
01332 
01333 {
01334 
01335 double angle;
01336 
01337 
01338 if (sun_zenith == 0)
01339         puts("WARNING: zenith_angle = 0 in function get_angle_sun_vert_plan");
01340 
01341 angle = acos( cos(sun_zenith*M_PI/180)*cos(direction_zenith*M_PI/180) + sin(sun_zenith*M_PI/180)*sin(direction_zenith*M_PI/180)*cos((sun_azimut-direction_azimut)*M_PI/180) );
01342 angle = angle*180/M_PI;
01343 return(angle);
01344 }
01345 
01346 
01347 
01348 
01349 
01350 
01351 
01352