ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/genssky.c
(Generate patch)

Comparing ray/src/gen/genssky.c (file contents):
Revision 2.2 by greg, Fri Jul 19 23:38:28 2024 UTC vs.
Revision 2.3 by greg, Fri Aug 2 18:47:25 2024 UTC

# Line 8 | Line 8 | static const char RCSid[] = "$Id$";
8   #include "copyright.h"
9   #include "resolu.h"
10   #include "rtio.h"
11 – #include "view.h"
11   #include <ctype.h>
12   #ifdef _WIN32
13   #include <windows.h>
# Line 61 | Line 60 | static int make_directory(const char *path) {
60   #endif
61   }
62  
63 + inline static float deg2rad(float deg) { return deg * (PI / 180.); }
64 +
65   static int cvthour(char *hs, int *tsolar, double *hour) {
66    char *cp = hs;
67    int i, j;
# Line 131 | Line 132 | static void basename(const char *path, char *output, s
132    }
133   }
134  
135 < char *join_paths(const char *path1, const char *path2) {
135 > static char *join_paths(const char *path1, const char *path2) {
136    size_t len1 = strlen(path1);
137    size_t len2 = strlen(path2);
138    int need_separator = (path1[len1 - 1] != DIRSEP);
# Line 183 | Line 184 | static void write_header(const int argc, char **argv,
184      printf("%s ", argv[i]);
185    }
186    printf("\n");
187 <  printf("#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
188 <         cloud_cover, grefl, res);
187 >  printf(
188 >      "#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
189 >      cloud_cover, grefl, res);
190   }
191  
192   static void write_rad(const double *sun_radiance, const FVECT sundir,
193 <                      const char skyfile[PATH_MAX],
192 <                      const char grndfile[PATH_MAX]) {
193 >                      const char *ddir, const char *skyfile) {
194    if (sundir[2] > 0) {
195      printf("void spectrum sunrad\n0\n0\n22 380 780 ");
196      /* Normalize to one */
# Line 207 | Line 208 | static void write_rad(const double *sun_radiance, cons
208      printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
209             sundir[2]);
210    }
211 <  printf("void specpict skymap\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
212 <         "-mx\n0\n0\n\n",
211 >  printf("void specpict skyfunc\n5 noop %s . 'atan2(Dy,Dx)/PI+1' "
212 >         "'acos(Dz)/PI'\n0\n0\n\n",
213           skyfile);
213 –
214 –  printf("void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
215 –         "-my\n0\n0\n\n",
216 –         grndfile);
217 –  printf("void mixfunc skyfunc\n4 skymap grndmap if(Dz,1,0) .\n0\n0\n");
214   }
215  
216   static void write_hsr_header(FILE *fp, RESOLU *res) {
217 <  float wvsplit[4] = {380, 480, 588,
222 <                      780}; /* RGB wavelength limits+partitions (nm) */
217 >  float wvsplit[4] = {380, 480, 588, 780};
218    newheader("RADIANCE", fp);
219    fputncomp(NSSAMP, fp);
220    fputwlsplit(wvsplit, fp);
# Line 228 | Line 223 | static void write_hsr_header(FILE *fp, RESOLU *res) {
223    fputsresolu(res, fp);
224   }
225  
226 + static inline float frac(float x) { return x - floor(x); }
227 +
228   int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
229                    DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
230                    const double cloud_cover, const FVECT sundir,
231 <                  const double grefl, const int res, const char *outname) {
231 >                  const double grefl, const int res, const char *outname,
232 >                  const char *ddir) {
233    char skyfile[PATH_MAX];
234    char grndfile[PATH_MAX];
235 <  if (!snprintf(skyfile, sizeof(skyfile), "%s_sky.hsr", outname)) {
235 >  if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP,
236 >                outname)) {
237      fprintf(stderr, "Error setting sky file name\n");
238      return 0;
239    };
240 <  if (!snprintf(grndfile, sizeof(grndfile), "%s_ground.hsr", outname)) {
241 <    fprintf(stderr, "Error setting ground file name\n");
242 <    return 0;
244 <  }
245 <  RESOLU rs = {PIXSTANDARD, res, res};
240 >  int xres = res;
241 >  int yres = xres / 2;
242 >  RESOLU rs = {PIXSTANDARD, xres, yres};
243    FILE *skyfp = fopen(skyfile, "w");
247 –  FILE *grndfp = fopen(grndfile, "w");
248 –  write_hsr_header(grndfp, &rs);
244    write_hsr_header(skyfp, &rs);
250 –  VIEW skyview = {VT_ANG, {0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, 1.,
251 –                  180.,   180.,         0.,           0.,           0.,
252 –                  0.,     {0., 0., 0.}, {0., 0., 0.}, 0.,           0.};
253 –  VIEW grndview = {
254 –      VT_ANG, {0., 0., 0.}, {0., 0., -1.}, {0., 1., 0.}, 1., 180., 180., 0., 0.,
255 –      0.,     0.,           {0., 0., 0.},  {0., 0., 0.}, 0., 0.};
256 –  setview(&skyview);
257 –  setview(&grndview);
245  
246    CNDX[3] = NSSAMP;
247  
248 <  FVECT view_point = {0, 0, ER};
248 >  FVECT view_point = {0, 0, ER + 10};
249    const double radius = VLEN(view_point);
250    const double sun_ct = fdot(view_point, sundir) / radius;
251 <  for (unsigned int j = 0; j < res; ++j) {
252 <    for (unsigned int i = 0; i < res; ++i) {
253 <      RREAL loc[2];
267 <      FVECT rorg = {0};
268 <      FVECT rdir_sky = {0};
269 <      FVECT rdir_grnd = {0};
270 <      SCOLOR sky_radiance = {0};
271 <      SCOLOR ground_radiance = {0};
251 >  for (int j = 0; j < yres; ++j) {
252 >    for (int i = 0; i < xres; ++i) {
253 >      SCOLOR radiance = {0};
254        SCOLR sky_sclr = {0};
273 –      SCOLR ground_sclr = {0};
255  
256 <      pix2loc(loc, &rs, i, j);
257 <      viewray(rorg, rdir_sky, &skyview, loc[0], loc[1]);
258 <      viewray(rorg, rdir_grnd, &grndview, loc[0], loc[1]);
256 >      float px = i / (xres - 1.0);
257 >      float py = j / (yres - 1.0);
258 >      float lambda = ((1 - py) * PI) - (PI / 2.0);
259 >      float phi = (px * 2.0 * PI) - PI;
260  
261 <      const double mu_sky = fdot(view_point, rdir_sky) / radius;
262 <      const double nu_sky = fdot(rdir_sky, sundir);
261 >      FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
262 >                    sin(lambda)};
263  
264 <      const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
265 <      const double nu_grnd = fdot(rdir_grnd, sundir);
264 >      const double mu = fdot(view_point, rdir) / radius;
265 >      const double nu = fdot(rdir, sundir);
266  
267 <      get_sky_radiance(scat_clear, scat1m_clear, radius, mu_sky, sun_ct, nu_sky,
268 <                       sky_radiance);
269 <      get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
270 <                          view_point, rdir_grnd, radius, mu_grnd, sun_ct,
271 <                          nu_grnd, grefl, sundir, ground_radiance);
267 >      /* hit ground */
268 >      if (rdir[2] < 0) {
269 >        get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
270 >                            view_point, rdir, radius, mu, sun_ct, nu, grefl,
271 >                            sundir, radiance);
272 >      } else {
273 >        get_sky_radiance(scat_clear, scat1m_clear, radius, mu, sun_ct, nu,
274 >                         radiance);
275 >      }
276  
277        for (int k = 0; k < NSSAMP; ++k) {
278 <        sky_radiance[k] *= WVLSPAN;
293 <        ground_radiance[k] *= WVLSPAN;
278 >        radiance[k] *= WVLSPAN;
279        }
280  
281        if (cloud_cover > 0) {
282          double zenithbr = get_zenith_brightness(sundir);
283          double grndbr = zenithbr * GNORM;
284 <        double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
285 <        for (int k = 0; k < NSSAMP; ++k) {
286 <          sky_radiance[k] =
287 <              wmean2(sky_radiance[k], skybr * D6415[k], cloud_cover);
288 <          ground_radiance[k] =
289 <              wmean2(ground_radiance[k], grndbr * D6415[k], cloud_cover);
284 >        double skybr = get_overcast_brightness(rdir[2], zenithbr);
285 >        if (rdir[2] < 0) {
286 >          for (int k = 0; k < NSSAMP; ++k) {
287 >            radiance[k] = wmean2(radiance[k], grndbr * D6415[k], cloud_cover);
288 >          }
289 >        } else {
290 >          for (int k = 0; k < NSSAMP; ++k) {
291 >            radiance[k] = wmean2(radiance[k], skybr * D6415[k], cloud_cover);
292 >          }
293          }
294        }
295  
296 <      scolor2scolr(sky_sclr, sky_radiance, 20);
296 >      scolor2scolr(sky_sclr, radiance, 20);
297        putbinary(sky_sclr, LSCOLR, 1, skyfp);
310 –
311 –      scolor2scolr(ground_sclr, ground_radiance, 20);
312 –      putbinary(ground_sclr, LSCOLR, 1, grndfp);
298      }
299    }
300    fclose(skyfp);
316 –  fclose(grndfp);
301  
302    /* Get solar radiance */
303    double sun_radiance[NSSAMP] = {0};
# Line 328 | Line 312 | int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_
312      }
313    }
314  
315 <  write_rad(sun_radiance, sundir, skyfile, grndfile);
315 >  write_rad(sun_radiance, sundir, ddir, skyfile);
316    return 1;
317   }
318  
# Line 351 | Line 335 | static DpPaths get_dppaths(const char *dir, const doub
335   static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
336                                           const int is_summer,
337                                           const double s_latitude) {
354 –  /* Set rayleigh density profile */
338    if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
339      tag[0] = 's';
340      if (is_summer) {
# Line 424 | Line 407 | int main(int argc, char *argv[]) {
407    int got_meridian = 0;
408    double grefl = 0.2;
409    double ccover = 0.0;
410 <  int res = 128;
410 >  int res = 64;
411    double aod = AOD0_CA;
412    char *outname = "out";
413    char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
# Line 432 | Line 415 | int main(int argc, char *argv[]) {
415    char lstag[3];
416    char *ddir = ".";
417  
418 <  if (!strcmp(argv[1], "-defaults")) {
418 >  if (argc == 2 && !strcmp(argv[1], "-defaults")) {
419      printf("-i %d\t\t\t\t#scattering order\n", sorder);
420      printf("-g %f\t\t\t#ground reflectance\n", grefl);
421      printf("-c %f\t\t\t#cloud cover\n", ccover);
# Line 440 | Line 423 | int main(int argc, char *argv[]) {
423      printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
424      printf("-f %s\t\t\t\t#output name (-f)\n", outname);
425      printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
426 <    exit(1);
426 >    exit(0);
427    }
428  
429    if (argc < 4) {
# Line 547 | Line 530 | int main(int argc, char *argv[]) {
530  
531    char gsdir[PATH_MAX];
532    size_t siz = strlen(ddir);
533 <  if (ISDIRSEP(ddir[siz-1]))
534 <    ddir[siz-1] = '\0';
533 >  if (ISDIRSEP(ddir[siz - 1]))
534 >    ddir[siz - 1] = '\0';
535    snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
553 –  printf("gsdir: %s\n", gsdir);
536    if (!make_directory(gsdir)) {
537      fprintf(stderr, "Failed creating atmos_data directory");
538      exit(1);
# Line 576 | Line 558 | int main(int argc, char *argv[]) {
558    write_header(argc, argv, ccover, grefl, res);
559  
560    if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
561 <                     irrad_clear_dp, ccover, sundir, grefl, res, outname)) {
561 >                     irrad_clear_dp, ccover, sundir, grefl, res, outname,
562 >                     ddir)) {
563      fprintf(stderr, "gen_spect_sky failed\n");
564      exit(1);
565    }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines