1 |
+ |
#include "color.h" |
2 |
|
#ifndef lint |
3 |
< |
static const char RCSid[] = "$Id$"; |
3 |
> |
static const char RCSid[] = |
4 |
> |
"$Id$"; |
5 |
|
#endif |
6 |
|
/* Main function for generating spectral sky */ |
7 |
|
/* Cloudy sky computed as weight average of clear and cie overcast sky */ |
162 |
|
return (a * x + b * y) / (a + b); |
163 |
|
} |
164 |
|
|
165 |
< |
static double get_zenith_brightness(const double sundir[3]) { |
165 |
> |
static double get_overcast_zenith_brightness(const double sundir[3]) { |
166 |
|
double zenithbr; |
167 |
|
if (sundir[2] < 0) { |
168 |
|
zenithbr = 0; |
192 |
|
cloud_cover, grefl, res); |
193 |
|
} |
194 |
|
|
195 |
< |
static void write_rad(const double *sun_radiance, const FVECT sundir, |
196 |
< |
const char *ddir, const char *skyfile) { |
195 |
> |
static void write_rad(const double *sun_radiance, const double intensity, |
196 |
> |
const FVECT sundir, const char *ddir, |
197 |
> |
const char *skyfile) { |
198 |
|
if (sundir[2] > 0) { |
199 |
|
printf("void spectrum sunrad\n0\n0\n22 380 780 "); |
197 |
– |
/* Normalize to one */ |
198 |
– |
double sum = 0.0; |
200 |
|
int i; |
201 |
|
for (i = 0; i < NSSAMP; ++i) { |
202 |
< |
sum += sun_radiance[i]; |
202 |
> |
printf("%.3f ", sun_radiance[i]); |
203 |
|
} |
203 |
– |
double mean = sum / NSSAMP; |
204 |
– |
for (i = 0; i < NSSAMP; ++i) { |
205 |
– |
printf("%.3f ", sun_radiance[i] / mean); |
206 |
– |
} |
207 |
– |
double intensity = mean * WVLSPAN; |
204 |
|
printf("\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n", intensity, |
205 |
|
intensity, intensity); |
206 |
|
printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1], |
227 |
|
DATARRAY *scat1m_clear, DATARRAY *irrad_clear, |
228 |
|
const double cloud_cover, const FVECT sundir, |
229 |
|
const double grefl, const int res, const char *outname, |
230 |
< |
const char *ddir) { |
230 |
> |
const char *ddir, const double dirnorm, const double difhor) { |
231 |
|
char skyfile[PATH_MAX]; |
236 |
– |
char grndfile[PATH_MAX]; |
232 |
|
if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP, |
233 |
|
outname)) { |
234 |
|
fprintf(stderr, "Error setting sky file name\n"); |
245 |
|
FVECT view_point = {0, 0, ER + 10}; |
246 |
|
const double radius = VLEN(view_point); |
247 |
|
const double sun_ct = fdot(view_point, sundir) / radius; |
248 |
+ |
|
249 |
+ |
double overcast_zenithbr = get_overcast_zenith_brightness(sundir); |
250 |
+ |
double overcast_grndbr = overcast_zenithbr * GNORM; |
251 |
+ |
|
252 |
+ |
double dif_ratio = 1; |
253 |
+ |
if (difhor > 0) { |
254 |
+ |
DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad_clear, radius, sun_ct); |
255 |
+ |
double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0; |
256 |
+ |
double diffuse_irradiance = 0; |
257 |
+ |
int l; |
258 |
+ |
for (l = 0; l < NSSAMP; ++l) { |
259 |
+ |
diffuse_irradiance += indirect_irradiance_clear->arr.d[l] * 20; /* 20nm interval */ |
260 |
+ |
} |
261 |
+ |
free(indirect_irradiance_clear); |
262 |
+ |
diffuse_irradiance = wmean2(diffuse_irradiance, overcast_ghi, cloud_cover); |
263 |
+ |
dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15; /* fudge */ |
264 |
+ |
} |
265 |
|
int i, j, k; |
266 |
|
for (j = 0; j < yres; ++j) { |
267 |
|
for (i = 0; i < xres; ++i) { |
294 |
|
} |
295 |
|
|
296 |
|
if (cloud_cover > 0) { |
297 |
< |
double zenithbr = get_zenith_brightness(sundir); |
286 |
< |
double grndbr = zenithbr * GNORM; |
287 |
< |
double skybr = get_overcast_brightness(rdir[2], zenithbr); |
297 |
> |
double skybr = get_overcast_brightness(rdir[2], overcast_zenithbr); |
298 |
|
if (rdir[2] < 0) { |
299 |
|
for (k = 0; k < NSSAMP; ++k) { |
300 |
< |
radiance[k] = wmean2(radiance[k], grndbr * D6415[k], cloud_cover); |
300 |
> |
radiance[k] = wmean2(radiance[k], overcast_grndbr * D6415[k], cloud_cover); |
301 |
|
} |
302 |
|
} else { |
303 |
|
for (k = 0; k < NSSAMP; ++k) { |
306 |
|
} |
307 |
|
} |
308 |
|
|
309 |
< |
scolor2scolr(sky_sclr, radiance, 20); |
309 |
> |
for (k = 0; k < NSSAMP; ++k) { |
310 |
> |
radiance[k] *= dif_ratio; |
311 |
> |
} |
312 |
> |
|
313 |
> |
scolor2scolr(sky_sclr, radiance, NSSAMP); |
314 |
|
putbinary(sky_sclr, LSCOLR, 1, skyfp); |
315 |
|
} |
316 |
|
} |
321 |
|
get_solar_radiance(tau_clear, scat_clear, scat1m_clear, sundir, radius, |
322 |
|
sun_ct, sun_radiance); |
323 |
|
if (cloud_cover > 0) { |
324 |
< |
double zenithbr = get_zenith_brightness(sundir); |
311 |
< |
double skybr = get_overcast_brightness(sundir[2], zenithbr); |
324 |
> |
double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr); |
325 |
|
int i; |
326 |
|
for (i = 0; i < NSSAMP; ++i) { |
327 |
|
sun_radiance[i] = |
329 |
|
} |
330 |
|
} |
331 |
|
|
332 |
< |
write_rad(sun_radiance, sundir, ddir, skyfile); |
332 |
> |
/* Normalize */ |
333 |
> |
double sum = 0.0; |
334 |
> |
for (i = 0; i < NSSAMP; ++i) { |
335 |
> |
sum += sun_radiance[i]; |
336 |
> |
} |
337 |
> |
double mean = sum / NSSAMP; |
338 |
> |
for (i = 0; i < NSSAMP; ++i) { |
339 |
> |
sun_radiance[i] /= mean; |
340 |
> |
} |
341 |
> |
double intensity = mean * WVLSPAN; |
342 |
> |
if (dirnorm > 0) { |
343 |
> |
intensity = dirnorm / SOLOMG / WHTEFFICACY; |
344 |
> |
} |
345 |
> |
|
346 |
> |
write_rad(sun_radiance, intensity, sundir, ddir, skyfile); |
347 |
|
return 1; |
348 |
|
} |
349 |
|
|
446 |
|
char lstag[3]; |
447 |
|
char *ddir = "."; |
448 |
|
int i; |
449 |
+ |
double dirnorm = 0; /* direct normal illuminance */ |
450 |
+ |
double difhor = 0; /* diffuse horizontal illuminance */ |
451 |
|
|
452 |
|
if (argc == 2 && !strcmp(argv[1], "-defaults")) { |
453 |
|
printf("-i %d\t\t\t\t#scattering order\n", sorder); |
463 |
|
if (argc < 4) { |
464 |
|
fprintf(stderr, |
465 |
|
"Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r " |
466 |
< |
"res -n nproc -c ccover -l mie -g grefl -f outpath\n", |
466 |
> |
"res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum " |
467 |
> |
"-g grefl -f outpath\n", |
468 |
|
argv[0]); |
469 |
|
return 0; |
470 |
|
} |
524 |
|
case 'o': |
525 |
|
s_longitude = atof(argv[++i]) * (PI / 180.0); |
526 |
|
break; |
527 |
+ |
case 'L': |
528 |
+ |
dirnorm = atof(argv[++i]); |
529 |
+ |
difhor = atof(argv[++i]); |
530 |
+ |
break; |
531 |
|
case 'p': |
532 |
|
ddir = argv[++i]; |
533 |
|
break; |
597 |
|
write_header(argc, argv, ccover, grefl, res); |
598 |
|
|
599 |
|
if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, |
600 |
< |
irrad_clear_dp, ccover, sundir, grefl, res, outname, |
601 |
< |
ddir)) { |
600 |
> |
irrad_clear_dp, ccover, sundir, grefl, res, outname, ddir, |
601 |
> |
dirnorm, difhor)) { |
602 |
|
fprintf(stderr, "gen_spect_sky failed\n"); |
603 |
|
exit(1); |
604 |
|
} |