← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.sun2/rsunlib.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.sun2/rsunlib.c
1
 
 
2
 
/*******************************************************************************
3
 
r.sun: rsunlib.c. This program was writen by Jaro Hofierka in Summer 1993 and re-engineered
4
 
in 1996-1999. In cooperation with Marcel Suri and Thomas Huld from JRC in Ispra
5
 
a new version of r.sun was prepared using ESRA solar radiation formulas.
6
 
See manual pages for details.
7
 
(C) 2002 Copyright Jaro Hofierka, Gresaka 22, 085 01 Bardejov, Slovakia, 
8
 
              and GeoModel, s.r.o., Bratislava, Slovakia
9
 
email: hofierka@geomodel.sk, marcel.suri@jrc.it, suri@geomodel.sk
10
 
*******************************************************************************/
11
 
/*
12
 
 * This program is free software; you can redistribute it and/or
13
 
 * modify it under the terms of the GNU General Public License
14
 
 * as published by the Free Software Foundation; either version 2
15
 
 * of the License, or (at your option) any later version.
16
 
 *
17
 
 * This program is distributed in the hope that it will be useful,
18
 
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19
 
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20
 
 * GNU General Public License for more details.
21
 
 *
22
 
 * You should have received a copy of the GNU General Public License
23
 
 * along with this program; if not, write to the Free Software
24
 
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25
 
 */
26
 
 
27
 
/*v. 2.0 July 2002, NULL data handling, JH */
28
 
/*v. 2.1 January 2003, code optimization by Thomas Huld, JH */
29
 
 
30
 
#include <stdio.h>
31
 
#include <math.h>
32
 
#include <grass/gis.h>
33
 
#include <grass/gprojects.h>
34
 
#include <grass/glocale.h>
35
 
#include "sunradstruct.h"
36
 
#include "local_proto.h"
37
 
#include "rsunglobals.h"
38
 
 
39
 
 
40
 
int civilTimeFlag;
41
 
int useCivilTime()
42
 
{
43
 
    return civilTimeFlag;
44
 
}
45
 
void setUseCivilTime(int val)
46
 
{
47
 
    civilTimeFlag = val;
48
 
}
49
 
 
50
 
 
51
 
double angular_loss_denom;
52
 
 
53
 
void setAngularLossDenominator()
54
 
{
55
 
    angular_loss_denom = 1. / (1 - exp(-1. / a_r));
56
 
 
57
 
}
58
 
 
59
 
 
60
 
int useShadowFlag;
61
 
int useShadow()
62
 
{
63
 
    return useShadowFlag;
64
 
}
65
 
void setUseShadow(int val)
66
 
{
67
 
    useShadowFlag = val;
68
 
}
69
 
 
70
 
 
71
 
 
72
 
int useHorizonDataFlag;
73
 
int useHorizonData()
74
 
{
75
 
    return useHorizonDataFlag;
76
 
}
77
 
void setUseHorizonData(int val)
78
 
{
79
 
    useHorizonDataFlag = val;
80
 
}
81
 
 
82
 
double timeOffset;
83
 
double getTimeOffset()
84
 
{
85
 
    return timeOffset;
86
 
}
87
 
void setTimeOffset(double val)
88
 
{
89
 
    timeOffset = val;
90
 
}
91
 
 
92
 
double horizonInterval;
93
 
double getHorizonInterval()
94
 
{
95
 
    return horizonInterval;
96
 
}
97
 
void setHorizonInterval(double val)
98
 
{
99
 
    horizonInterval = val;
100
 
}
101
 
 
102
 
 
103
 
/* com_sol_const(): compute the Solar Constant corrected for the day of the
104
 
   year. The Earth is closest to the Sun (Perigee) on about January 3rd,
105
 
   it is furthest from the sun (Apogee) about July 6th. The 1367 W/m^2 solar
106
 
   constant is at the average 1AU distance, but on Jan 3 it gets up to
107
 
   around 1412.71 W/m^2 and on July 6 it gets down to around 1321 W/m^2.
108
 
   This value is for what hits the top of the atmosphere before any energy
109
 
   is attenuated. */
110
 
double com_sol_const(int no_of_day)
111
 
{
112
 
    double I0, d1;
113
 
 
114
 
    /* Solar constant: 1367.0 W/m^2.
115
 
 
116
 
       Perigee offset: here we call Jan 2 at 8:18pm the Perigee, so day
117
 
        number 2.8408. In angular units that's (2*pi * 2.8408 / 365.25) = 0.048869.
118
 
 
119
 
       Orbital eccentricity: For Earth this is currently about 0.01672,
120
 
        and so the distance to the sun varies by +/- 0.01672 from the
121
 
        mean distance (1AU), so over the year the amplitude of the
122
 
        function is 2*ecc = 0.03344.
123
 
 
124
 
       And 365.25 is of course the number of days in a year.
125
 
    */
126
 
 
127
 
    /*  v W/(m*m) */
128
 
    d1 = pi2 * no_of_day / 365.25;
129
 
    I0 = 1367. * (1 + 0.03344 * cos(d1 - 0.048869));
130
 
 
131
 
    return I0;
132
 
}
133
 
 
134
 
 
135
 
void com_par_const(double longitTime, struct SunGeometryConstDay *sungeom,
136
 
                   struct GridGeometry *gridGeom)
137
 
{
138
 
    double pom;
139
 
    double totOffsetTime;
140
 
 
141
 
    sungeom->lum_C11 = gridGeom->sinlat * sungeom->cosdecl;
142
 
    sungeom->lum_C13 = -gridGeom->coslat * sungeom->sindecl;
143
 
    sungeom->lum_C22 = sungeom->cosdecl;
144
 
    sungeom->lum_C31 = gridGeom->coslat * sungeom->cosdecl;
145
 
    sungeom->lum_C33 = gridGeom->sinlat * sungeom->sindecl;
146
 
 
147
 
    if (fabs(sungeom->lum_C31) >= EPS) {
148
 
        totOffsetTime = timeOffset + longitTime;
149
 
 
150
 
        if (useCivilTime()) {
151
 
            sungeom->timeAngle -= totOffsetTime * HOURANGLE;
152
 
        }
153
 
        pom = -sungeom->lum_C33 / sungeom->lum_C31;
154
 
        if (fabs(pom) <= 1) {
155
 
            pom = acos(pom);
156
 
            pom = (pom * 180) / M_PI;
157
 
            sungeom->sunrise_time = (90 - pom) / 15 + 6;
158
 
            sungeom->sunset_time = (pom - 90) / 15 + 18;
159
 
        }
160
 
        else {
161
 
            if (pom < 0) {
162
 
                /* G_debug(3,"\n Sun is ABOVE the surface during the whole day"); */
163
 
                sungeom->sunrise_time = 0;
164
 
                sungeom->sunset_time = 24;
165
 
            }
166
 
            else {
167
 
                /* G_debug(3,"\n The sun is BELOW the surface during the whole day"); */
168
 
                if (fabs(pom) - 1 <= EPS) {
169
 
                    sungeom->sunrise_time = 12;
170
 
                    sungeom->sunset_time = 12;
171
 
                }
172
 
            }
173
 
        }
174
 
    }
175
 
 
176
 
}
177
 
 
178
 
 
179
 
 
180
 
 
181
 
 
182
 
void com_par(struct SunGeometryConstDay *sungeom,
183
 
             struct SunGeometryVarDay *sunVarGeom,
184
 
             struct GridGeometry *gridGeom, double latitude, double longitude)
185
 
{
186
 
    double pom, xpom, ypom;
187
 
 
188
 
    double costimeAngle;
189
 
    double lum_Lx, lum_Ly;
190
 
 
191
 
    double newLatitude, newLongitude;
192
 
    double inputAngle;
193
 
    double delt_lat, delt_lon;
194
 
    double delt_east, delt_nor;
195
 
    double delt_dist;
196
 
 
197
 
 
198
 
    costimeAngle = cos(sungeom->timeAngle);
199
 
 
200
 
 
201
 
    lum_Lx = -sungeom->lum_C22 * sin(sungeom->timeAngle);
202
 
    lum_Ly = sungeom->lum_C11 * costimeAngle + sungeom->lum_C13;
203
 
    sunVarGeom->sinSolarAltitude =
204
 
        sungeom->lum_C31 * costimeAngle + sungeom->lum_C33;
205
 
 
206
 
    if (fabs(sungeom->lum_C31) < EPS) {
207
 
        if (fabs(sunVarGeom->sinSolarAltitude) >= EPS) {
208
 
            if (sunVarGeom->sinSolarAltitude > 0) {
209
 
                /* G_debug(3,"\tSun is ABOVE area during the whole day"); */
210
 
                sungeom->sunrise_time = 0;
211
 
                sungeom->sunset_time = 24;
212
 
            }
213
 
            else {
214
 
                sunVarGeom->solarAltitude = 0.;
215
 
                sunVarGeom->solarAzimuth = UNDEF;
216
 
                return;
217
 
            }
218
 
        }
219
 
        else {
220
 
            /* G_debug(3,"\tThe Sun is ON HORIZON during the whole day"); */
221
 
            sungeom->sunrise_time = 0;
222
 
            sungeom->sunset_time = 24;
223
 
        }
224
 
    }
225
 
 
226
 
    sunVarGeom->solarAltitude = asin(sunVarGeom->sinSolarAltitude);     /* vertical angle of the sun */
227
 
    /* sinSolarAltitude is sin(solarAltitude) */
228
 
 
229
 
    xpom = lum_Lx * lum_Lx;
230
 
    ypom = lum_Ly * lum_Ly;
231
 
    pom = sqrt(xpom + ypom);
232
 
 
233
 
 
234
 
    if (fabs(pom) > EPS) {
235
 
        sunVarGeom->solarAzimuth = lum_Ly / pom;
236
 
        sunVarGeom->solarAzimuth = acos(sunVarGeom->solarAzimuth);      /* horiz. angle of the Sun */
237
 
        /* solarAzimuth *= RAD; */
238
 
        if (lum_Lx < 0)
239
 
            sunVarGeom->solarAzimuth = pi2 - sunVarGeom->solarAzimuth;
240
 
    }
241
 
    else {
242
 
        sunVarGeom->solarAzimuth = UNDEF;
243
 
    }
244
 
 
245
 
 
246
 
 
247
 
    if (sunVarGeom->solarAzimuth < 0.5 * M_PI)
248
 
        sunVarGeom->sunAzimuthAngle = 0.5 * M_PI - sunVarGeom->solarAzimuth;
249
 
    else
250
 
        sunVarGeom->sunAzimuthAngle = 2.5 * M_PI - sunVarGeom->solarAzimuth;
251
 
 
252
 
 
253
 
    inputAngle = sunVarGeom->sunAzimuthAngle + pihalf;
254
 
    inputAngle = (inputAngle >= pi2) ? inputAngle - pi2 : inputAngle;
255
 
 
256
 
 
257
 
    delt_lat = -0.0001 * cos(inputAngle);  /* Arbitrary small distance in latitude */
258
 
    delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
259
 
 
260
 
    newLatitude = (latitude + delt_lat) * rad2deg;
261
 
    newLongitude = (longitude + delt_lon) * rad2deg;
262
 
 
263
 
 
264
 
    if ((G_projection() != PROJECTION_LL)) {
265
 
        if (pj_do_proj(&newLongitude, &newLatitude, &oproj, &iproj) < 0) {
266
 
            G_fatal_error("Error in pj_do_proj");
267
 
        }
268
 
    }
269
 
 
270
 
    delt_east = newLongitude - gridGeom->xp;
271
 
    delt_nor = newLatitude - gridGeom->yp;
272
 
 
273
 
    delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor);
274
 
 
275
 
 
276
 
    sunVarGeom->stepsinangle = gridGeom->stepxy * delt_nor / delt_dist;
277
 
    sunVarGeom->stepcosangle = gridGeom->stepxy * delt_east / delt_dist;
278
 
 
279
 
    /*
280
 
       sunVarGeom->stepsinangle = stepxy * sin(sunVarGeom->sunAzimuthAngle);
281
 
       sunVarGeom->stepcosangle = stepxy * cos(sunVarGeom->sunAzimuthAngle);
282
 
     */
283
 
    sunVarGeom->tanSolarAltitude = tan(sunVarGeom->solarAltitude);
284
 
 
285
 
    return;
286
 
 
287
 
}
288
 
 
289
 
 
290
 
int searching(double *length, struct SunGeometryVarDay *sunVarGeom,
291
 
              struct GridGeometry *gridGeom)
292
 
{
293
 
    double z2;
294
 
    double curvature_diff;
295
 
    int succes = 0;
296
 
 
297
 
    if (sunVarGeom->zp == UNDEFZ)
298
 
        return 0;
299
 
 
300
 
 
301
 
    gridGeom->yy0 += sunVarGeom->stepsinangle;
302
 
    gridGeom->xx0 += sunVarGeom->stepcosangle;
303
 
    if (((gridGeom->xx0 + (0.5 * gridGeom->stepx)) < 0)
304
 
        || ((gridGeom->xx0 + (0.5 * gridGeom->stepx)) > gridGeom->deltx)
305
 
        || ((gridGeom->yy0 + (0.5 * gridGeom->stepy)) < 0)
306
 
        || ((gridGeom->yy0 + (0.5 * gridGeom->stepy)) > gridGeom->delty))
307
 
        succes = 3;
308
 
    else
309
 
        succes = 1;
310
 
 
311
 
 
312
 
    if (succes == 1) {
313
 
        where_is_point(length, sunVarGeom, gridGeom);
314
 
        if (func == NULL) {
315
 
            gridGeom->xx0 = gridGeom->xg0;
316
 
            gridGeom->yy0 = gridGeom->yg0;
317
 
            return (3);
318
 
        }
319
 
        curvature_diff = EARTHRADIUS * (1. - cos(*length / EARTHRADIUS));
320
 
 
321
 
        z2 = sunVarGeom->z_orig + curvature_diff +
322
 
            *length * sunVarGeom->tanSolarAltitude;
323
 
        if (z2 < sunVarGeom->zp)
324
 
            succes = 2;         /* shadow */
325
 
        if (z2 > sunVarGeom->zmax)
326
 
            succes = 3;         /* no test needed all visible */
327
 
    }
328
 
 
329
 
    if (succes != 1) {
330
 
        gridGeom->xx0 = gridGeom->xg0;
331
 
        gridGeom->yy0 = gridGeom->yg0;
332
 
    }
333
 
    return (succes);
334
 
}
335
 
 
336
 
 
337
 
 
338
 
 
339
 
double lumcline2(struct SunGeometryConstDay *sungeom,
340
 
                 struct SunGeometryVarDay *sunVarGeom,
341
 
                 struct SunGeometryVarSlope *sunSlopeGeom,
342
 
                 struct GridGeometry *gridGeom, unsigned char *horizonpointer)
343
 
{
344
 
    double s = 0;
345
 
    double length;
346
 
    int r = 0;
347
 
 
348
 
    int lowPos, highPos;
349
 
    double timeoffset, horizPos;
350
 
    double horizonHeight;
351
 
 
352
 
    func = cube;
353
 
    sunVarGeom->isShadow = 0;
354
 
 
355
 
    if (useShadow()) {
356
 
        length = 0;
357
 
 
358
 
        if (useHorizonData()) {
359
 
            /* Start is due east, sungeom->timeangle = -pi/2 */
360
 
            /* timeoffset = sungeom->timeAngle+pihalf; */
361
 
            timeoffset = sunVarGeom->sunAzimuthAngle;
362
 
 
363
 
            /*
364
 
               if(timeoffset<0.)
365
 
                  timeoffset+=pi2;
366
 
               else if(timeoffset>pi2)
367
 
                  timeoffset-=pi2;
368
 
               horizPos = arrayNumInt - timeoffset/horizonInterval;
369
 
             */
370
 
 
371
 
            horizPos = timeoffset / getHorizonInterval();
372
 
 
373
 
 
374
 
            lowPos = (int)horizPos;
375
 
            highPos = lowPos + 1;
376
 
            if (highPos == arrayNumInt) {
377
 
                highPos = 0;
378
 
            }
379
 
            horizonHeight = invScale * ((1. -
380
 
                                         (horizPos -
381
 
                                          lowPos)) * horizonpointer[lowPos]
382
 
                                        + (horizPos - lowPos)
383
 
                                        * horizonpointer[highPos]);
384
 
            sunVarGeom->isShadow =
385
 
                (horizonHeight > sunVarGeom->solarAltitude);
386
 
 
387
 
            if (!sunVarGeom->isShadow) {
388
 
                /* if (z_orig != UNDEFZ) {
389
 
                      s = sunSlopeGeom->lum_C31_l
390
 
                        * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
391
 
                        + sunSlopeGeom->lum_C33_l;
392
 
                   } else {
393
 
                     s = sunVarGeom->sinSolarAltitude;
394
 
                   }
395
 
                 */
396
 
                s = sunSlopeGeom->lum_C31_l
397
 
                        * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
398
 
                        + sunSlopeGeom->lum_C33_l;      /* Jenco */
399
 
            }
400
 
 
401
 
        }                       /*  End if useHorizonData() */
402
 
        else {
403
 
            while ((r = searching(&length, sunVarGeom, gridGeom)) == 1) {
404
 
                if (r == 3)
405
 
                    break;      /* no test is needed */
406
 
            }
407
 
 
408
 
 
409
 
            if (r == 2) {
410
 
                sunVarGeom->isShadow = 1;       /* shadow */
411
 
            }
412
 
            else {
413
 
 
414
 
                /* if (z_orig != UNDEFZ) {
415
 
                      s = sunSlopeGeom->lum_C31_l
416
 
                        * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
417
 
                        + sunSlopeGeom->lum_C33_l;
418
 
                   } else {
419
 
                     s = sunVarGeom->sinSolarAltitude;
420
 
                   }
421
 
                 */
422
 
                s = sunSlopeGeom->lum_C31_l
423
 
                        * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
424
 
                        + sunSlopeGeom->lum_C33_l;      /* Jenco */
425
 
            }
426
 
        }
427
 
    }
428
 
    else {
429
 
        /* if (z_orig != UNDEFZ) {
430
 
             s = sunSlopeGeom->lum_C31_l
431
 
                * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
432
 
                + sunSlopeGeom->lum_C33_l;
433
 
           } else {
434
 
             s = sunVarGeom->sinSolarAltitude;
435
 
           }
436
 
         */
437
 
        s = sunSlopeGeom->lum_C31_l
438
 
                * cos(-sungeom->timeAngle - sunSlopeGeom->longit_l)
439
 
                + sunSlopeGeom->lum_C33_l;      /* Jenco */
440
 
    }
441
 
 
442
 
    /* if (s <= 0) return UNDEFZ; ?? */
443
 
    if (s < 0)
444
 
        return 0.;
445
 
 
446
 
    return (s);
447
 
}
448
 
 
449
 
 
450
 
 
451
 
double brad(double sh, double *bh, struct SunGeometryVarDay *sunVarGeom,
452
 
            struct SunGeometryVarSlope *sunSlopeGeom,
453
 
            struct SolarRadVar *sunRadVar)
454
 
{
455
 
    double opticalAirMass, airMass2Linke, rayl, br;
456
 
    double drefract, temp1, temp2, h0refract;
457
 
    double elevationCorr;
458
 
 
459
 
    double locSolarAltitude;
460
 
 
461
 
    locSolarAltitude = sunVarGeom->solarAltitude;
462
 
 
463
 
/* FIXME: please document all coefficients */
464
 
    elevationCorr = exp(-sunVarGeom->z_orig / 8434.5);
465
 
    temp1 = 0.1594 + locSolarAltitude * (1.123 + 0.065656 * locSolarAltitude);
466
 
    temp2 = 1. + locSolarAltitude * (28.9344 + 277.3971 * locSolarAltitude);
467
 
    drefract = 0.061359 * temp1 / temp2;        /* in radians */
468
 
    h0refract = locSolarAltitude + drefract;
469
 
    opticalAirMass = elevationCorr / (sin(h0refract) +
470
 
                                      0.50572 * pow(h0refract * rad2deg +
471
 
                                                    6.07995, -1.6364));
472
 
    airMass2Linke = 0.8662 * sunRadVar->linke;
473
 
    if (opticalAirMass <= 20.) {
474
 
        rayl = 1. / (6.6296 +
475
 
                     opticalAirMass * (1.7513 +
476
 
                                       opticalAirMass * (-0.1202 +
477
 
                                                         opticalAirMass *
478
 
                                                         (0.0065 -
479
 
                                                          opticalAirMass *
480
 
                                                          0.00013))));
481
 
    }
482
 
    else {
483
 
        rayl = 1. / (10.4 + 0.718 * opticalAirMass);
484
 
    }
485
 
    *bh =
486
 
        sunRadVar->cbh * sunRadVar->G_norm_extra *
487
 
        sunVarGeom->sinSolarAltitude * exp(-rayl * opticalAirMass *
488
 
                                           airMass2Linke);
489
 
    if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.)
490
 
        br = *bh * sh / sunVarGeom->sinSolarAltitude;
491
 
    else
492
 
        br = *bh;
493
 
 
494
 
    return (br);
495
 
}
496
 
 
497
 
double brad_angle_loss(double sh, double *bh,
498
 
                       struct SunGeometryVarDay *sunVarGeom,
499
 
                       struct SunGeometryVarSlope *sunSlopeGeom,
500
 
                       struct SolarRadVar *sunRadVar)
501
 
{
502
 
    double p, opticalAirMass, airMass2Linke, rayl, br;
503
 
    double drefract, temp1, temp2, h0refract;
504
 
 
505
 
    double locSolarAltitude;
506
 
 
507
 
    locSolarAltitude = sunVarGeom->solarAltitude;
508
 
 
509
 
/* FIXME: please document all coefficients */
510
 
    p = exp(-sunVarGeom->z_orig / 8434.5);
511
 
    temp1 = 0.1594 + locSolarAltitude * (1.123 + 0.065656 * locSolarAltitude);
512
 
    temp2 = 1. + locSolarAltitude * (28.9344 + 277.3971 * locSolarAltitude);
513
 
    drefract = 0.061359 * temp1 / temp2;        /* in radians */
514
 
    h0refract = locSolarAltitude + drefract;
515
 
    opticalAirMass = p / (sin(h0refract) +
516
 
                          0.50572 * pow(h0refract * rad2deg + 6.07995,
517
 
                                        -1.6364));
518
 
    airMass2Linke = 0.8662 * sunRadVar->linke;
519
 
    if (opticalAirMass <= 20.)
520
 
        rayl =
521
 
            1. / (6.6296 +
522
 
                  opticalAirMass *
523
 
                   (1.7513 + opticalAirMass *
524
 
                    (-0.1202 + opticalAirMass *
525
 
                     (0.0065 - opticalAirMass * 0.00013))));
526
 
    else
527
 
        rayl = 1. / (10.4 + 0.718 * opticalAirMass);
528
 
    *bh =
529
 
        sunRadVar->cbh * sunRadVar->G_norm_extra *
530
 
        sunVarGeom->sinSolarAltitude * exp(-rayl * opticalAirMass *
531
 
                                           airMass2Linke);
532
 
    if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.)
533
 
        br = *bh * sh / sunVarGeom->sinSolarAltitude;
534
 
    else
535
 
        br = *bh;
536
 
 
537
 
    br *= (1 - exp(-sh / a_r)) * angular_loss_denom;
538
 
 
539
 
    return (br);
540
 
}
541
 
 
542
 
 
543
 
 
544
 
double drad(double sh, double bh, double *rr,
545
 
            struct SunGeometryVarDay *sunVarGeom,
546
 
            struct SunGeometryVarSlope *sunSlopeGeom,
547
 
            struct SolarRadVar *sunRadVar)
548
 
{
549
 
    double tn, fd, fx = 0., A1, A2, A3, A1b;
550
 
    double r_sky, kb, dr, gh, a_ln, ln, fg;
551
 
    double dh;
552
 
    double cosslope, sinslope;
553
 
    double locSinSolarAltitude;
554
 
    double locLinke;
555
 
 
556
 
    locLinke = sunRadVar->linke;
557
 
    locSinSolarAltitude = sunVarGeom->sinSolarAltitude;
558
 
 
559
 
    cosslope = cos(sunSlopeGeom->slope);
560
 
    sinslope = sin(sunSlopeGeom->slope);
561
 
 
562
 
/* FIXME: please document all coefficients */
563
 
    tn = -0.015843 + locLinke * (0.030543 + 0.0003797 * locLinke);
564
 
    A1b = 0.26463 + locLinke * (-0.061581 + 0.0031408 * locLinke);
565
 
    if (A1b * tn < 0.0022)
566
 
        A1 = 0.0022 / tn;
567
 
    else
568
 
        A1 = A1b;
569
 
    A2 = 2.04020 + locLinke * (0.018945 - 0.011161 * locLinke);
570
 
    A3 = -1.3025 + locLinke * (0.039231 + 0.0085079 * locLinke);
571
 
 
572
 
    fd = A1 + A2 * locSinSolarAltitude +
573
 
        A3 * locSinSolarAltitude * locSinSolarAltitude;
574
 
    dh = sunRadVar->cdh * sunRadVar->G_norm_extra * fd * tn;
575
 
    gh = bh + dh;
576
 
    if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.) {
577
 
        kb = bh / (sunRadVar->G_norm_extra * locSinSolarAltitude);
578
 
        r_sky = (1. + cosslope) / 2.;
579
 
        a_ln = sunVarGeom->solarAzimuth - sunSlopeGeom->aspect;
580
 
        ln = a_ln;
581
 
        if (a_ln > M_PI)
582
 
            ln = a_ln - pi2;
583
 
        else if (a_ln < -M_PI)
584
 
            ln = a_ln + pi2;
585
 
        a_ln = ln;
586
 
        fg = sinslope - sunSlopeGeom->slope * cosslope -
587
 
            M_PI * sin(0.5 * sunSlopeGeom->slope) * sin(0.5 *
588
 
                                                        sunSlopeGeom->slope);
589
 
        if ((sunVarGeom->isShadow == 1) || sh <= 0.)
590
 
            fx = r_sky + fg * 0.252271;
591
 
        else if (sunVarGeom->solarAltitude >= 0.1) {
592
 
            fx = ((0.00263 - kb * (0.712 + 0.6883 * kb)) * fg + r_sky) *
593
 
                  (1. - kb) + kb * sh / locSinSolarAltitude;
594
 
        }
595
 
        else if (sunVarGeom->solarAltitude < 0.1)
596
 
            fx = ((0.00263 - 0.712 * kb - 0.6883 * kb * kb) * fg +
597
 
                  r_sky) * (1. - kb) + kb *
598
 
                  sinslope * cos(a_ln) /
599
 
                  (0.1 - 0.008 * sunVarGeom->solarAltitude);
600
 
        dr = dh * fx;
601
 
        /* refl. rad */
602
 
        *rr = sunRadVar->alb * gh * (1 - cosslope) / 2.;
603
 
    }
604
 
    else {      /* plane */
605
 
        dr = dh;
606
 
        *rr = 0.;
607
 
    }
608
 
    return (dr);
609
 
}
610
 
 
611
 
#define c2 -0.074
612
 
 
613
 
double drad_angle_loss(double sh, double bh, double *rr,
614
 
                       struct SunGeometryVarDay *sunVarGeom,
615
 
                       struct SunGeometryVarSlope *sunSlopeGeom,
616
 
                       struct SolarRadVar *sunRadVar)
617
 
{
618
 
    double dh;
619
 
    double tn, fd, fx = 0., A1, A2, A3, A1b;
620
 
    double r_sky, kb, dr, gh, a_ln, ln, fg;
621
 
    double cosslope, sinslope;
622
 
 
623
 
    double diff_coeff_angleloss;
624
 
    double refl_coeff_angleloss;
625
 
    double c1;
626
 
    double diff_loss_factor, refl_loss_factor;
627
 
    double locSinSolarAltitude;
628
 
    double locLinke;
629
 
 
630
 
    locSinSolarAltitude = sunVarGeom->sinSolarAltitude;
631
 
    locLinke = sunRadVar->linke;
632
 
    cosslope = cos(sunSlopeGeom->slope);
633
 
    sinslope = sin(sunSlopeGeom->slope);
634
 
 
635
 
/* FIXME: please document all coefficients */
636
 
    tn = -0.015843 + locLinke * (0.030543 + 0.0003797 * locLinke);
637
 
    A1b = 0.26463 + locLinke * (-0.061581 + 0.0031408 * locLinke);
638
 
 
639
 
    if (A1b * tn < 0.0022)
640
 
        A1 = 0.0022 / tn;
641
 
    else
642
 
        A1 = A1b;
643
 
    A2 = 2.04020 + locLinke * (0.018945 - 0.011161 * locLinke);
644
 
    A3 = -1.3025 + locLinke * (0.039231 + 0.0085079 * locLinke);
645
 
 
646
 
    fd = A1 + A2 * locSinSolarAltitude +
647
 
        A3 * locSinSolarAltitude * locSinSolarAltitude;
648
 
    dh = sunRadVar->cdh * sunRadVar->G_norm_extra * fd * tn;
649
 
    gh = bh + dh;
650
 
 
651
 
    if (sunSlopeGeom->aspect != UNDEF && sunSlopeGeom->slope != 0.) {
652
 
 
653
 
        kb = bh / (sunRadVar->G_norm_extra * locSinSolarAltitude);
654
 
        r_sky = (1. + cosslope) / 2.;
655
 
        a_ln = sunVarGeom->solarAzimuth - sunSlopeGeom->aspect;
656
 
        ln = a_ln;
657
 
 
658
 
        if (a_ln > M_PI)
659
 
            ln = a_ln - pi2;
660
 
        else if (a_ln < -M_PI)
661
 
            ln = a_ln + pi2;
662
 
        a_ln = ln;
663
 
        fg = sinslope - sunSlopeGeom->slope * cosslope -
664
 
            M_PI * sin(sunSlopeGeom->slope / 2.) * sin(sunSlopeGeom->slope /
665
 
                                                       2.);
666
 
        if ((sunVarGeom->isShadow) || (sh <= 0.))
667
 
            fx = r_sky + fg * 0.252271;
668
 
        else if (sunVarGeom->solarAltitude >= 0.1) {
669
 
            fx = ((0.00263 - kb * (0.712 + 0.6883 * kb)) * fg + r_sky) * (1. -
670
 
                                                                          kb)
671
 
                + kb * sh / locSinSolarAltitude;
672
 
        }
673
 
        else if (sunVarGeom->solarAltitude < 0.1)
674
 
            fx = ((0.00263 - 0.712 * kb - 0.6883 * kb * kb) * fg +
675
 
                  r_sky) * (1. - kb) + kb * sinslope * cos(a_ln) /
676
 
                  (0.1 - 0.008 * sunVarGeom->solarAltitude);
677
 
 
678
 
        dr = dh * fx;
679
 
        /* refl. rad */
680
 
        *rr = sunRadVar->alb * gh * (1 - cosslope) / 2.;
681
 
    }
682
 
    else {      /* plane */
683
 
        dr = dh;
684
 
        *rr = 0.;
685
 
    }
686
 
 
687
 
    c1 = 4. / (3. * M_PI);
688
 
    diff_coeff_angleloss = sinslope
689
 
        + (M_PI - sunSlopeGeom->slope - sinslope) / (1 + cosslope);
690
 
    refl_coeff_angleloss = sinslope
691
 
        + (sunSlopeGeom->slope - sinslope) / (1 - cosslope);
692
 
 
693
 
    diff_loss_factor
694
 
        = 1. - exp(-(c1 * diff_coeff_angleloss
695
 
                     + c2 * diff_coeff_angleloss * diff_coeff_angleloss)
696
 
                   / a_r);
697
 
    refl_loss_factor
698
 
        = 1. - exp(-(c1 * refl_coeff_angleloss
699
 
                     + c2 * refl_coeff_angleloss * refl_coeff_angleloss)
700
 
                   / a_r);
701
 
 
702
 
    dr *= diff_loss_factor;
703
 
    *rr *= refl_loss_factor;
704
 
 
705
 
 
706
 
 
707
 
    return (dr);
708
 
}