CrIS Code Collection
libcris.c
Go to the documentation of this file.
1#include "libcris.h"
2
3/*****************************************************************************/
4
5void add_att(
6 const int ncid,
7 const int varid,
8 const char *unit,
9 const char *long_name) {
10
11 /* Set long name... */
12 NC(nc_put_att_text(ncid, varid, "long_name", strlen(long_name), long_name));
13
14 /* Set units... */
15 NC(nc_put_att_text(ncid, varid, "units", strlen(unit), unit));
16}
17
18/*****************************************************************************/
19
20void add_var(
21 const int ncid,
22 const char *varname,
23 const char *unit,
24 const char *longname,
25 int type,
26 int dimid[],
27 int *varid,
28 int ndims) {
29
30 /* Check if variable exists... */
31 if (nc_inq_varid(ncid, varname, varid) != NC_NOERR) {
32
33 /* Define variable... */
34 NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
35
36 /* Set long name... */
37 NC(nc_put_att_text
38 (ncid, *varid, "long_name", strlen(longname), longname));
39
40 /* Set units... */
41 NC(nc_put_att_text(ncid, *varid, "units", strlen(unit), unit));
42 }
43}
44
45/*****************************************************************************/
46
47void background_poly_help(
48 const double *xx,
49 double *yy,
50 const int n,
51 const int dim) {
52
53 double chisq, xx2[WX > WY ? WX : WY], yy2[WX > WY ? WX : WY];
54
55 size_t n2 = 0;
56
57 /* Check for nan... */
58 for (size_t i = 0; i < (size_t) n; i++)
59 if (gsl_finite(yy[i])) {
60 xx2[n2] = xx[i];
61 yy2[n2] = yy[i];
62 n2++;
63 }
64 if ((int) n2 < dim || n2 < 0.9 * n) {
65 for (size_t i = 0; i < (size_t) n; i++)
66 yy[i] = GSL_NAN;
67 return;
68 }
69
70 /* Allocate... */
71 gsl_multifit_linear_workspace *work =
72 gsl_multifit_linear_alloc((size_t) n2, (size_t) dim);
73 gsl_matrix *cov = gsl_matrix_alloc((size_t) dim, (size_t) dim);
74 gsl_matrix *X = gsl_matrix_alloc((size_t) n2, (size_t) dim);
75 gsl_vector *c = gsl_vector_alloc((size_t) dim);
76 gsl_vector *x = gsl_vector_alloc((size_t) n2);
77 gsl_vector *y = gsl_vector_alloc((size_t) n2);
78
79 /* Compute polynomial fit... */
80 for (size_t i = 0; i < (size_t) n2; i++) {
81 gsl_vector_set(x, i, xx2[i]);
82 gsl_vector_set(y, i, yy2[i]);
83 for (size_t i2 = 0; i2 < (size_t) dim; i2++)
84 gsl_matrix_set(X, i, i2, pow(gsl_vector_get(x, i), (double) i2));
85 }
86 gsl_multifit_linear(X, y, c, cov, &chisq, work);
87 for (size_t i = 0; i < (size_t) n; i++)
88 yy[i] = gsl_poly_eval(c->data, (int) dim, xx[i]);
89
90 /* Free... */
91 gsl_multifit_linear_free(work);
92 gsl_matrix_free(cov);
93 gsl_matrix_free(X);
94 gsl_vector_free(c);
95 gsl_vector_free(x);
96 gsl_vector_free(y);
97}
98
99/*****************************************************************************/
100
101void background_poly(
102 wave_t *wave,
103 int dim_x,
104 int dim_y) {
105
106 /* Check parameters... */
107 if (dim_x <= 0 && dim_y <= 0)
108 return;
109
110 /* Copy temperatures to background... */
111 for (int ix = 0; ix < wave->nx; ix++)
112 for (int iy = 0; iy < wave->ny; iy++) {
113 wave->bg[ix][iy] = wave->temp[ix][iy];
114 wave->pt[ix][iy] = 0;
115 }
116
117 /* Compute fit in x-direction... */
118 if (dim_x > 0)
119 for (int iy = 0; iy < wave->ny; iy++) {
120 double x[WX], y[WX];
121 for (int ix = 0; ix < wave->nx; ix++) {
122 x[ix] = (double) ix;
123 y[ix] = wave->bg[ix][iy];
124 }
125 background_poly_help(x, y, wave->nx, dim_x);
126 for (int ix = 0; ix < wave->nx; ix++)
127 wave->bg[ix][iy] = y[ix];
128 }
129
130 /* Compute fit in y-direction... */
131 if (dim_y > 0)
132 for (int ix = 0; ix < wave->nx; ix++) {
133 double x[WX], y[WX];
134 for (int iy = 0; iy < wave->ny; iy++) {
135 x[iy] = (int) iy;
136 y[iy] = wave->bg[ix][iy];
137 }
138 background_poly_help(x, y, wave->ny, dim_y);
139 for (int iy = 0; iy < wave->ny; iy++)
140 wave->bg[ix][iy] = y[iy];
141 }
142
143 /* Recompute perturbations... */
144 for (int ix = 0; ix < wave->nx; ix++)
145 for (int iy = 0; iy < wave->ny; iy++)
146 wave->pt[ix][iy] = wave->temp[ix][iy] - wave->bg[ix][iy];
147}
148
149/*****************************************************************************/
150
151void background_smooth(
152 wave_t *wave,
153 int npts_x,
154 int npts_y) {
155
156 const double dmax = 2500.;
157
158 static double help[WX][WY];
159
160 int n;
161
162 /* Check parameters... */
163 if (npts_x <= 0 && npts_y <= 0)
164 return;
165
166 /* Smooth background... */
167 for (int ix = 0; ix < wave->nx; ix++)
168 for (int iy = 0; iy < wave->ny; iy++) {
169
170 /* Init... */
171 n = 0;
172 help[ix][iy] = 0;
173
174 /* Set maximum range... */
175 int dx = GSL_MIN(GSL_MIN(npts_x, ix), wave->nx - 1 - ix);
176 int dy = GSL_MIN(GSL_MIN(npts_y, iy), wave->ny - 1 - iy);
177
178 /* Average... */
179 for (int i = ix - dx; i <= ix + dx; i++)
180 for (int j = iy - dy; j <= iy + dy; j++)
181 if (fabs(wave->x[ix] - wave->x[i]) < dmax &&
182 fabs(wave->y[iy] - wave->y[j]) < dmax) {
183 help[ix][iy] += wave->bg[i][j];
184 n++;
185 }
186
187 /* Normalize... */
188 if (n > 0)
189 help[ix][iy] /= n;
190 else
191 help[ix][iy] = GSL_NAN;
192 }
193
194 /* Recalculate perturbations... */
195 for (int ix = 0; ix < wave->nx; ix++)
196 for (int iy = 0; iy < wave->ny; iy++) {
197 wave->bg[ix][iy] = help[ix][iy];
198 wave->pt[ix][iy] = wave->temp[ix][iy] - wave->bg[ix][iy];
199 }
200}
201
202/*****************************************************************************/
203
204void create_background(
205 wave_t *wave) {
206
207 /* Loop over grid points... */
208 for (int ix = 0; ix < wave->nx; ix++)
209 for (int iy = 0; iy < wave->ny; iy++) {
210
211 /* Set background for 4.3 micron BT measurements... */
212 wave->bg[ix][iy] =
213 235.626 + 5.38165e-6 * gsl_pow_2(wave->x[ix] -
214 0.5 * (wave->x[0] +
215 wave->x[wave->nx - 1])) -
216 1.78519e-12 * gsl_pow_4(wave->x[ix] -
217 0.5 * (wave->x[0] + wave->x[wave->nx - 1]));
218
219 /* Set temperature perturbation... */
220 wave->pt[ix][iy] = 0;
221
222 /* Set temperature... */
223 wave->temp[ix][iy] = wave->bg[ix][iy];
224 }
225}
226
227/*****************************************************************************/
228
229void create_noise(
230 wave_t *wave,
231 double nedt) {
232
233 /* Initialize random number generator... */
234 gsl_rng_env_setup();
235 gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
236 gsl_rng_set(r, (unsigned long int) time(NULL));
237
238 /* Add noise to temperature... */
239 if (nedt > 0)
240 for (int ix = 0; ix < wave->nx; ix++)
241 for (int iy = 0; iy < wave->ny; iy++)
242 wave->temp[ix][iy] += gsl_ran_gaussian(r, nedt);
243
244 /* Free... */
245 gsl_rng_free(r);
246}
247
248/*****************************************************************************/
249
250void create_wave(
251 wave_t *wave,
252 double amp,
253 double lx,
254 double ly,
255 double phi,
256 double fwhm) {
257
258 /* Loop over grid points... */
259 for (int ix = 0; ix < wave->nx; ix++)
260 for (int iy = 0; iy < wave->ny; iy++) {
261
262 /* Set wave perturbation... */
263 wave->pt[ix][iy] = amp * cos((lx != 0 ? 2 * M_PI / lx : 0) * wave->x[ix]
264 + (ly !=
265 0 ? 2 * M_PI / ly : 0) * wave->y[iy]
266 - phi * M_PI / 180.)
267 * (fwhm > 0 ? exp(-0.5 * gsl_pow_2((wave->x[ix]) / (lx * fwhm) * 2.35)
268 -
269 0.5 * gsl_pow_2((wave->y[iy]) / (ly * fwhm) *
270 2.35)) : 1.0);
271
272 /* Add perturbation to temperature... */
273 wave->temp[ix][iy] += wave->pt[ix][iy];
274 }
275}
276
277/*****************************************************************************/
278
279void day2doy(
280 const int year,
281 const int mon,
282 const int day,
283 int *doy) {
284
285 const int
286 d0[12] = { 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 },
287 d0l[12] = { 1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336 };
288
289 /* Get day of year... */
290 if (year % 400 == 0 || (year % 100 != 0 && year % 4 == 0))
291 *doy = d0l[mon - 1] + day - 1;
292 else
293 *doy = d0[mon - 1] + day - 1;
294}
295
296/*****************************************************************************/
297
298void doy2day(
299 const int year,
300 const int doy,
301 int *mon,
302 int *day) {
303
304 const int
305 d0[12] = { 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 },
306 d0l[12] = { 1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336 };
307
308 int i;
309
310 /* Get month and day... */
311 if (year % 400 == 0 || (year % 100 != 0 && year % 4 == 0)) {
312 for (i = 11; i > 0; i--)
313 if (d0l[i] <= doy)
314 break;
315 *mon = i + 1;
316 *day = doy - d0l[i] + 1;
317 } else {
318 for (i = 11; i > 0; i--)
319 if (d0[i] <= doy)
320 break;
321 *mon = i + 1;
322 *day = doy - d0[i] + 1;
323 }
324}
325
326/*****************************************************************************/
327
328void fit_wave(
329 wave_t *wave,
330 double amp,
331 double phi,
332 double kx,
333 double ky,
334 double *chisq) {
335
336 /* Init... */
337 *chisq = 0;
338
339 /* Calculate fit... */
340 for (int ix = 0; ix < wave->nx; ix++)
341 for (int iy = 0; iy < wave->ny; iy++) {
342 wave->fit[ix][iy]
343 = amp * cos(kx * wave->x[ix] + ky * wave->y[iy]
344 - phi * M_PI / 180.);
345 *chisq += POW2(wave->fit[ix][iy] - wave->pt[ix][iy]);
346 }
347
348 /* Calculate chisq... */
349 *chisq /= (double) (wave->nx * wave->ny);
350}
351
352/*****************************************************************************/
353
354void fft_help(
355 double *fcReal,
356 double *fcImag,
357 int n) {
358
359 double data[2 * PMAX];
360
361 /* Check size... */
362 if (n > PMAX)
363 ERRMSG("Too many data points!");
364
365 /* Allocate... */
366 gsl_fft_complex_wavetable *wavetable =
367 gsl_fft_complex_wavetable_alloc((size_t) n);
368 gsl_fft_complex_workspace *workspace =
369 gsl_fft_complex_workspace_alloc((size_t) n);
370
371 /* Set data (real, complex)... */
372 for (int i = 0; i < n; i++) {
373 data[2 * i] = fcReal[i];
374 data[2 * i + 1] = fcImag[i];
375 }
376
377 /* Calculate FFT... */
378 gsl_fft_complex_forward(data, 1, (size_t) n, wavetable, workspace);
379
380 /* Copy data... */
381 for (int i = 0; i < n; i++) {
382 fcReal[i] = data[2 * i];
383 fcImag[i] = data[2 * i + 1];
384 }
385
386 /* Free... */
387 gsl_fft_complex_wavetable_free(wavetable);
388 gsl_fft_complex_workspace_free(workspace);
389}
390
391/*****************************************************************************/
392
393void fft(
394 wave_t *wave,
395 double *Amax,
396 double *phimax,
397 double *lhmax,
398 double *kxmax,
399 double *kymax,
400 double *alphamax,
401 double *betamax,
402 char *filename) {
403
404 static double A[PMAX][PMAX], phi[PMAX][PMAX], kx[PMAX], ky[PMAX],
405 cutReal[PMAX], cutImag[PMAX], boxImag[PMAX][PMAX], boxReal[PMAX][PMAX];
406
407 int imin, imax, jmin, jmax;
408
409 /* Find box... */
410 imin = jmin = 9999;
411 imax = jmax = -9999;
412 for (int i = 0; i < wave->nx; i++)
413 for (int j = 0; j < wave->ny; j++)
414 if (gsl_finite(wave->var[i][j])) {
415 imin = GSL_MIN(imin, i);
416 imax = GSL_MAX(imax, i);
417 jmin = GSL_MIN(jmin, j);
418 jmax = GSL_MAX(jmax, j);
419 }
420 const int nx = imax - imin + 1;
421 const int ny = jmax - jmin + 1;
422
423 /* Copy data... */
424 for (int i = imin; i <= imax; i++)
425 for (int j = jmin; j <= jmax; j++) {
426 if (gsl_finite(wave->pt[i][j]))
427 boxReal[i - imin][j - jmin] = wave->pt[i][j];
428 else
429 boxReal[i - imin][j - jmin] = 0.0;
430 boxImag[i - imin][j - jmin] = 0.0;
431 }
432
433 /* FFT of the rows... */
434 for (int i = 0; i < nx; i++) {
435 for (int j = 0; j < ny; j++) {
436 cutReal[j] = boxReal[i][j];
437 cutImag[j] = boxImag[i][j];
438 }
439 fft_help(cutReal, cutImag, ny);
440 for (int j = 0; j < ny; j++) {
441 boxReal[i][j] = cutReal[j];
442 boxImag[i][j] = cutImag[j];
443 }
444 }
445
446 /* FFT of the columns... */
447 for (int j = 0; j < ny; j++) {
448 for (int i = 0; i < nx; i++) {
449 cutReal[i] = boxReal[i][j];
450 cutImag[i] = boxImag[i][j];
451 }
452 fft_help(cutReal, cutImag, nx);
453 for (int i = 0; i < nx; i++) {
454 boxReal[i][j] = cutReal[i];
455 boxImag[i][j] = cutImag[i];
456 }
457 }
458
459 /* Get frequencies, amplitude, and phase... */
460 for (int i = 0; i < nx; i++)
461 kx[i] = 2. * M_PI * ((i < nx / 2) ? (double) i : -(double) (nx - i))
462 / (nx * fabs(wave->x[imax] - wave->x[imin]) / (nx - 1.0));
463 for (int j = 0; j < ny; j++)
464 ky[j] = 2. * M_PI * ((j < ny / 2) ? (double) j : -(double) (ny - j))
465 / (ny * fabs(wave->y[jmax] - wave->y[jmin]) / (ny - 1.0));
466 for (int i = 0; i < nx; i++)
467 for (int j = 0; j < ny; j++) {
468 A[i][j]
469 = (i == 0 && j == 0 ? 1.0 : 2.0) / (nx * ny)
470 * sqrt(gsl_pow_2(boxReal[i][j]) + gsl_pow_2(boxImag[i][j]));
471 phi[i][j]
472 = 180. / M_PI * atan2(boxImag[i][j], boxReal[i][j]);
473 }
474
475 /* Check frequencies... */
476 for (int i = 0; i < nx; i++)
477 for (int j = 0; j < ny; j++)
478 if (kx[i] == 0 || ky[j] == 0) {
479 A[i][j] = GSL_NAN;
480 phi[i][j] = GSL_NAN;
481 }
482
483 /* Find maximum... */
484 *Amax = 0;
485 for (int i = 0; i < nx; i++)
486 for (int j = 0; j < ny / 2; j++)
487 if (gsl_finite(A[i][j]) && A[i][j] > *Amax) {
488 *Amax = A[i][j];
489 *phimax = phi[i][j];
490 *kxmax = kx[i];
491 *kymax = ky[j];
492 }
493
494 /* Get horizontal wavelength... */
495 *lhmax = 2 * M_PI / sqrt(gsl_pow_2(*kxmax) + gsl_pow_2(*kymax));
496
497 /* Get propagation direction in xy-plane... */
498 *alphamax = 90. - 180. / M_PI * atan2(*kxmax, *kymax);
499
500 /* Get propagation direction in lon,lat-plane... */
501 *betamax = *alphamax
502 +
503 180. / M_PI *
504 atan2(wave->lat[wave->nx / 2 >
505 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
506 - wave->lat[wave->nx / 2 <
507 wave->nx - 1 ? wave->nx / 2 +
508 1 : wave->nx / 2][wave->ny / 2],
509 wave->lon[wave->nx / 2 >
510 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
511 - wave->lon[wave->nx / 2 <
512 wave->nx - 1 ? wave->nx / 2 +
513 1 : wave->nx / 2][wave->ny / 2]);
514
515 /* Save FFT data... */
516 if (filename != NULL) {
517
518 /* Write info... */
519 LOG(1, "Write FFT data: %s", filename);
520
521 /* Create file... */
522 FILE *out;
523 if (!(out = fopen(filename, "w")))
524 ERRMSG("Cannot create file!");
525
526 /* Write header... */
527 fprintf(out,
528 "# $1 = altitude [km]\n"
529 "# $2 = wavelength in x-direction [km]\n"
530 "# $3 = wavelength in y-direction [km]\n"
531 "# $4 = wavenumber in x-direction [1/km]\n"
532 "# $5 = wavenumber in y-direction [1/km]\n"
533 "# $6 = amplitude [K]\n" "# $7 = phase [rad]\n");
534
535 /* Write data... */
536 for (int i = nx - 1; i > 0; i--) {
537 fprintf(out, "\n");
538 for (int j = ny / 2; j > 0; j--) {
539 int i2 = (i == nx / 2 ? 0 : i);
540 int j2 = (j == ny / 2 ? 0 : j);
541 fprintf(out, "%g %g %g %g %g %g %g\n", wave->z,
542 (kx[i2] != 0 ? 2 * M_PI / kx[i2] : 0),
543 (ky[j2] != 0 ? 2 * M_PI / ky[j2] : 0),
544 kx[i2], ky[j2], A[i2][j2], phi[i2][j2]);
545 }
546 }
547
548 /* Close file... */
549 fclose(out);
550 }
551}
552
553/*****************************************************************************/
554
555void gauss(
556 wave_t *wave,
557 double fwhm) {
558
559 static double help[WX][WY];
560
561 /* Check parameters... */
562 if (fwhm <= 0)
563 return;
564
565 /* Compute sigma^2... */
566 const double sigma2 = gsl_pow_2(fwhm / 2.3548);
567
568 /* Loop over data points... */
569 for (int ix = 0; ix < wave->nx; ix++)
570 for (int iy = 0; iy < wave->ny; iy++) {
571
572 /* Init... */
573 double wsum = 0;
574 help[ix][iy] = 0;
575
576 /* Average... */
577 for (int ix2 = 0; ix2 < wave->nx; ix2++)
578 for (int iy2 = 0; iy2 < wave->ny; iy2++) {
579 double d2 = gsl_pow_2(wave->x[ix] - wave->x[ix2])
580 + gsl_pow_2(wave->y[iy] - wave->y[iy2]);
581 if (d2 <= 9 * sigma2) {
582 double w = exp(-d2 / (2 * sigma2));
583 wsum += w;
584 help[ix][iy] += w * wave->pt[ix2][iy2];
585 }
586 }
587
588 /* Normalize... */
589 wave->pt[ix][iy] = help[ix][iy] / wsum;
590 }
591}
592
593/*****************************************************************************/
594
595void hamming(
596 wave_t *wave,
597 int niter) {
598
599 static double help[WX][WY];
600
601 /* Iterations... */
602 for (int iter = 0; iter < niter; iter++) {
603
604 /* Filter in x direction... */
605 for (int ix = 0; ix < wave->nx; ix++)
606 for (int iy = 0; iy < wave->ny; iy++)
607 help[ix][iy]
608 = 0.23 * wave->pt[ix > 0 ? ix - 1 : ix][iy]
609 + 0.54 * wave->pt[ix][iy]
610 + 0.23 * wave->pt[ix < wave->nx - 1 ? ix + 1 : ix][iy];
611
612 /* Filter in y direction... */
613 for (int ix = 0; ix < wave->nx; ix++)
614 for (int iy = 0; iy < wave->ny; iy++)
615 wave->pt[ix][iy]
616 = 0.23 * help[ix][iy > 0 ? iy - 1 : iy]
617 + 0.54 * help[ix][iy]
618 + 0.23 * help[ix][iy < wave->ny - 1 ? iy + 1 : iy];
619 }
620}
621
622/*****************************************************************************/
623
624void intpol_x(
625 wave_t *wave,
626 int n) {
627
628 double dummy, x[WX], xc[WX][3], xc2[WX][3], y[WX];
629
630 /* Check parameters... */
631 if (n <= 0)
632 return;
633 if (n > WX)
634 ERRMSG("Too many data points!");
635
636 /* Set new x-coordinates... */
637 for (int i = 0; i < n; i++)
638 x[i] = LIN(0.0, wave->x[0], n - 1.0, wave->x[wave->nx - 1], i);
639
640 /* Allocate... */
641 gsl_interp_accel *acc = gsl_interp_accel_alloc();
642 gsl_spline *spline =
643 gsl_spline_alloc(gsl_interp_cspline, (size_t) wave->nx);
644
645 /* Loop over scans... */
646 for (int iy = 0; iy < wave->ny; iy++) {
647
648 /* Interpolate Cartesian coordinates... */
649 for (int ix = 0; ix < wave->nx; ix++)
650 geo2cart(0, wave->lon[ix][iy], wave->lat[ix][iy], xc[ix]);
651 for (int ic = 0; ic < 3; ic++) {
652 for (int ix = 0; ix < wave->nx; ix++)
653 y[ix] = xc[ix][ic];
654 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
655 for (int i = 0; i < n; i++)
656 xc2[i][ic] = gsl_spline_eval(spline, x[i], acc);
657 }
658 for (int i = 0; i < n; i++)
659 cart2geo(xc2[i], &dummy, &wave->lon[i][iy], &wave->lat[i][iy]);
660
661 /* Interpolate temperature... */
662 for (int ix = 0; ix < wave->nx; ix++)
663 y[ix] = wave->temp[ix][iy];
664 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
665 for (int i = 0; i < n; i++)
666 wave->temp[i][iy] = gsl_spline_eval(spline, x[i], acc);
667
668 /* Interpolate background... */
669 for (int ix = 0; ix < wave->nx; ix++)
670 y[ix] = wave->bg[ix][iy];
671 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
672 for (int i = 0; i < n; i++)
673 wave->bg[i][iy] = gsl_spline_eval(spline, x[i], acc);
674
675 /* Interpolate perturbations... */
676 for (int ix = 0; ix < wave->nx; ix++)
677 y[ix] = wave->pt[ix][iy];
678 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
679 for (int i = 0; i < n; i++)
680 wave->pt[i][iy] = gsl_spline_eval(spline, x[i], acc);
681
682 /* Interpolate variance... */
683 for (int ix = 0; ix < wave->nx; ix++)
684 y[ix] = wave->var[ix][iy];
685 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
686 for (int i = 0; i < n; i++)
687 wave->var[i][iy] = gsl_spline_eval(spline, x[i], acc);
688 }
689
690 /* Free... */
691 gsl_spline_free(spline);
692 gsl_interp_accel_free(acc);
693
694 /* Set new x-coordinates... */
695 for (int i = 0; i < n; i++)
696 wave->x[i] = x[i];
697 wave->nx = n;
698}
699
700/*****************************************************************************/
701
702void median(
703 wave_t *wave,
704 int dx) {
705
706 static double data[WX * WY], help[WX][WY];
707
708 size_t n;
709
710 /* Check parameters... */
711 if (dx <= 0)
712 return;
713
714 /* Loop over data points... */
715 for (int ix = 0; ix < wave->nx; ix++)
716 for (int iy = 0; iy < wave->ny; iy++) {
717
718 /* Init... */
719 n = 0;
720
721 /* Get data... */
722 for (int ix2 = GSL_MAX(ix - dx, 0);
723 ix2 < GSL_MIN(ix + dx, wave->nx - 1); ix2++)
724 for (int iy2 = GSL_MAX(iy - dx, 0);
725 iy2 < GSL_MIN(iy + dx, wave->ny - 1); iy2++) {
726 data[n] = wave->pt[ix2][iy2];
727 n++;
728 }
729
730 /* Normalize... */
731 gsl_sort(data, 1, n);
732 help[ix][iy] = gsl_stats_median_from_sorted_data(data, 1, n);
733 }
734
735 /* Loop over data points... */
736 for (int ix = 0; ix < wave->nx; ix++)
737 for (int iy = 0; iy < wave->ny; iy++)
738 wave->pt[ix][iy] = help[ix][iy];
739}
740
741/*****************************************************************************/
742
743void merge_y(
744 wave_t *wave1,
745 wave_t *wave2) {
746
747 /* Check data... */
748 if (wave1->nx != wave2->nx)
749 ERRMSG("Across-track sizes do not match!");
750 if (wave1->ny + wave2->ny > WY)
751 ERRMSG("Too many data points!");
752
753 /* Get offset in y direction... */
754 const double y =
755 wave1->y[wave1->ny - 1] + (wave1->y[wave1->ny - 1] -
756 wave1->y[0]) / (wave1->ny - 1);
757
758 /* Merge data... */
759 for (int ix = 0; ix < wave2->nx; ix++)
760 for (int iy = 0; iy < wave2->ny; iy++) {
761 wave1->y[wave1->ny + iy] = y + wave2->y[iy];
762 wave1->lon[ix][wave1->ny + iy] = wave2->lon[ix][iy];
763 wave1->lat[ix][wave1->ny + iy] = wave2->lat[ix][iy];
764 wave1->temp[ix][wave1->ny + iy] = wave2->temp[ix][iy];
765 wave1->bg[ix][wave1->ny + iy] = wave2->bg[ix][iy];
766 wave1->pt[ix][wave1->ny + iy] = wave2->pt[ix][iy];
767 wave1->var[ix][wave1->ny + iy] = wave2->var[ix][iy];
768 }
769
770 /* Increment counter... */
771 wave1->ny += wave2->ny;
772}
773
774/*****************************************************************************/
775
776void noise(
777 wave_t *wave,
778 double *mu,
779 double *sig) {
780
781 int n = 0;
782
783 /* Init... */
784 *mu = 0;
785 *sig = 0;
786
787 /* Estimate noise (Immerkaer, 1996)... */
788 for (int ix = 1; ix < wave->nx - 1; ix++)
789 for (int iy = 1; iy < wave->ny - 1; iy++) {
790
791 /* Check data... */
792 int okay = 1;
793 for (int ix2 = ix - 1; ix2 <= ix + 1; ix2++)
794 for (int iy2 = iy - 1; iy2 <= iy + 1; iy2++)
795 if (!gsl_finite(wave->temp[ix2][iy2]))
796 okay = 0;
797 if (!okay)
798 continue;
799
800 /* Get mean noise... */
801 n++;
802 *mu += wave->temp[ix][iy];
803 *sig += gsl_pow_2(+4. / 6. * wave->temp[ix][iy]
804 - 2. / 6. * (wave->temp[ix - 1][iy]
805 + wave->temp[ix + 1][iy]
806 + wave->temp[ix][iy - 1]
807 + wave->temp[ix][iy + 1])
808 + 1. / 6. * (wave->temp[ix - 1][iy - 1]
809 + wave->temp[ix + 1][iy - 1]
810 + wave->temp[ix - 1][iy + 1]
811 + wave->temp[ix + 1][iy + 1]));
812 }
813
814 /* Normalize... */
815 *mu /= (double) n;
816 *sig = sqrt(*sig / (double) n);
817}
818
819/*****************************************************************************/
820
821void noise_pert(
822 pert_t *pert,
823 int track0,
824 int track1,
825 double *mu,
826 double *sig) {
827
828 int n = 0;
829
830 /* Init... */
831 *mu = 0;
832 *sig = 0;
833
834 /* Estimate noise (Immerkaer, 1996)... */
835 for (int track = track0; track < track1; track++)
836 for (int xtrack = 0; xtrack < pert->nxtrack; xtrack++) {
837
838 /* Check data... */
839 int okay = 1;
840 for (int ifov = 0; ifov < pert->nfov; ifov++)
841 if (!(gsl_finite(pert->bt[track][xtrack][ifov])
842 && gsl_finite(pert->pt[track][xtrack][ifov])))
843 okay = 0;
844 if (!okay)
845 continue;
846
847 /* Get mean noise... */
848 n++;
849 *mu += pert->bt[track][xtrack][4];
850 *sig += gsl_pow_2(+4. / 6. * pert->pt[track][xtrack][4]
851 - 2. / 6. * (pert->pt[track][xtrack][1]
852 + pert->pt[track][xtrack][3]
853 + pert->pt[track][xtrack][5]
854 + pert->pt[track][xtrack][7])
855 + 1. / 6. * (pert->pt[track][xtrack][0]
856 + pert->pt[track][xtrack][2]
857 + pert->pt[track][xtrack][6]
858 + pert->pt[track][xtrack][8]));
859 }
860
861 /* Normalize... */
862 *mu /= (double) n;
863 *sig = sqrt(*sig / (double) n);
864}
865
866/*****************************************************************************/
867
868void period(
869 wave_t *wave,
870 double lxymax,
871 double dlxy,
872 double *Amax,
873 double *phimax,
874 double *lhmax,
875 double *kxmax,
876 double *kymax,
877 double *alphamax,
878 double *betamax,
879 char *filename) {
880
881 static double kx[PMAX], ky[PMAX], A[PMAX][PMAX], phi[PMAX][PMAX],
882 cx[PMAX][WX], cy[PMAX][WY], sx[PMAX][WX], sy[PMAX][WY], a, b, c;
883
884 int imin, imax, jmin, jmax, lmax = 0, mmax = 0;
885
886 /* Compute wavenumbers and periodogram coefficients... */
887 for (double lx = -lxymax; lx <= lxymax; lx += dlxy) {
888 kx[lmax] = (lx != 0 ? 2 * M_PI / lx : 0);
889 for (int i = 0; i < wave->nx; i++) {
890 cx[lmax][i] = cos(kx[lmax] * wave->x[i]);
891 sx[lmax][i] = sin(kx[lmax] * wave->x[i]);
892 }
893 if ((++lmax) > PMAX)
894 ERRMSG("Too many wavenumbers for periodogram!");
895 }
896 for (double ly = 0; ly <= lxymax; ly += dlxy) {
897 ky[mmax] = (ly != 0 ? 2 * M_PI / ly : 0);
898 for (int j = 0; j < wave->ny; j++) {
899 cy[mmax][j] = cos(ky[mmax] * wave->y[j]);
900 sy[mmax][j] = sin(ky[mmax] * wave->y[j]);
901 }
902 if ((++mmax) > PMAX)
903 ERRMSG("Too many wavenumbers for periodogram!");
904 }
905
906 /* Find area... */
907 imin = jmin = 9999;
908 imax = jmax = -9999;
909 for (int i = 0; i < wave->nx; i++)
910 for (int j = 0; j < wave->ny; j++)
911 if (gsl_finite(wave->var[i][j])) {
912 imin = GSL_MIN(imin, i);
913 imax = GSL_MAX(imax, i);
914 jmin = GSL_MIN(jmin, j);
915 jmax = GSL_MAX(jmax, j);
916 }
917
918 /* Get Nyquist frequencies... */
919 const double kx_ny =
920 M_PI / fabs((wave->x[imax] - wave->x[imin]) /
921 ((double) imax - (double) imin));
922 const double ky_ny =
923 M_PI / fabs((wave->y[jmax] - wave->y[jmin]) /
924 ((double) jmax - (double) jmin));
925
926 /* Loop over wavelengths... */
927 for (int l = 0; l < lmax; l++)
928 for (int m = 0; m < mmax; m++) {
929
930 /* Check frequencies... */
931 if (kx[l] == 0 || fabs(kx[l]) > kx_ny ||
932 ky[m] == 0 || fabs(ky[m]) > ky_ny) {
933 A[l][m] = GSL_NAN;
934 phi[l][m] = GSL_NAN;
935 continue;
936 }
937
938 /* Compute periodogram... */
939 a = b = c = 0;
940 for (int i = imin; i <= imax; i++)
941 for (int j = jmin; j <= jmax; j++)
942 if (gsl_finite(wave->var[i][j])) {
943 a += wave->pt[i][j] * (cx[l][i] * cy[m][j] - sx[l][i] * sy[m][j]);
944 b += wave->pt[i][j] * (sx[l][i] * cy[m][j] + cx[l][i] * sy[m][j]);
945 c++;
946 }
947 a *= 2. / c;
948 b *= 2. / c;
949
950 /* Get amplitude and phase... */
951 A[l][m] = sqrt(gsl_pow_2(a) + gsl_pow_2(b));
952 phi[l][m] = atan2(b, a) * 180. / M_PI;
953 }
954
955 /* Find maximum... */
956 *Amax = 0;
957 for (int l = 0; l < lmax; l++)
958 for (int m = 0; m < mmax; m++)
959 if (gsl_finite(A[l][m]) && A[l][m] > *Amax) {
960 *Amax = A[l][m];
961 *phimax = phi[l][m];
962 *kxmax = kx[l];
963 *kymax = ky[m];
964 }
965
966 /* Get horizontal wavelength... */
967 *lhmax = 2 * M_PI / sqrt(gsl_pow_2(*kxmax) + gsl_pow_2(*kymax));
968
969 /* Get propagation direction in xy-plane... */
970 *alphamax = 90. - 180. / M_PI * atan2(*kxmax, *kymax);
971
972 /* Get propagation direction in lon,lat-plane... */
973 *betamax = *alphamax
974 +
975 180. / M_PI *
976 atan2(wave->lat[wave->nx / 2 >
977 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
978 - wave->lat[wave->nx / 2 <
979 wave->nx - 1 ? wave->nx / 2 +
980 1 : wave->nx / 2][wave->ny / 2],
981 wave->lon[wave->nx / 2 >
982 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
983 - wave->lon[wave->nx / 2 <
984 wave->nx - 1 ? wave->nx / 2 +
985 1 : wave->nx / 2][wave->ny / 2]);
986
987 /* Save periodogram data... */
988 if (filename != NULL) {
989
990 /* Write info... */
991 LOG(1, "Write periodogram data: %s", filename);
992
993 /* Create file... */
994 FILE *out;
995 if (!(out = fopen(filename, "w")))
996 ERRMSG("Cannot create file!");
997
998 /* Write header... */
999 fprintf(out,
1000 "# $1 = altitude [km]\n"
1001 "# $2 = wavelength in x-direction [km]\n"
1002 "# $3 = wavelength in y-direction [km]\n"
1003 "# $4 = wavenumber in x-direction [1/km]\n"
1004 "# $5 = wavenumber in y-direction [1/km]\n"
1005 "# $6 = amplitude [K]\n" "# $7 = phase [rad]\n");
1006
1007 /* Write data... */
1008 for (int l = 0; l < lmax; l++) {
1009 fprintf(out, "\n");
1010 for (int m = 0; m < mmax; m++)
1011 fprintf(out, "%g %g %g %g %g %g %g\n", wave->z,
1012 (kx[l] != 0 ? 2 * M_PI / kx[l] : 0),
1013 (ky[m] != 0 ? 2 * M_PI / ky[m] : 0),
1014 kx[l], ky[m], A[l][m], phi[l][m]);
1015 }
1016
1017 /* Close file... */
1018 fclose(out);
1019 }
1020}
1021
1022/*****************************************************************************/
1023
1024#if 0
1025void pert2wave(
1026 pert_t *pert,
1027 wave_t *wave,
1028 int track0,
1029 int track1,
1030 int xtrack0,
1031 int xtrack1) {
1032
1033 double x0[3], x1[3];
1034
1035 /* Check ranges... */
1036 track0 = GSL_MIN(GSL_MAX(track0, 0), pert->ntrack - 1);
1037 track1 = GSL_MIN(GSL_MAX(track1, 0), pert->ntrack - 1);
1038 xtrack0 = GSL_MIN(GSL_MAX(xtrack0, 0), pert->nxtrack - 1);
1039 xtrack1 = GSL_MIN(GSL_MAX(xtrack1, 0), pert->nxtrack - 1);
1040
1041 /* Set size... */
1042 wave->nx = xtrack1 - xtrack0 + 1;
1043 if (wave->nx > WX)
1044 ERRMSG("Too many across-track values!");
1045 wave->ny = track1 - track0 + 1;
1046 if (wave->ny > WY)
1047 ERRMSG("Too many along-track values!");
1048
1049 /* Loop over footprints... */
1050 for (int itrack = track0; itrack <= track1; itrack++)
1051 for (int ixtrack = xtrack0; ixtrack <= xtrack1; ixtrack++) {
1052
1053 /* Get distances... */
1054 if (itrack == track0) {
1055 wave->x[0] = 0;
1056 if (ixtrack > xtrack0) {
1057 geo2cart(0, pert->lon[itrack][ixtrack - 1],
1058 pert->lat[itrack][ixtrack - 1], x0);
1059 geo2cart(0, pert->lon[itrack][ixtrack],
1060 pert->lat[itrack][ixtrack], x1);
1061 wave->x[ixtrack - xtrack0] =
1062 wave->x[ixtrack - xtrack0 - 1] + DIST(x0, x1);
1063 }
1064 }
1065 if (ixtrack == xtrack0) {
1066 wave->y[0] = 0;
1067 if (itrack > track0) {
1068 geo2cart(0, pert->lon[itrack - 1][ixtrack],
1069 pert->lat[itrack - 1][ixtrack], x0);
1070 geo2cart(0, pert->lon[itrack][ixtrack],
1071 pert->lat[itrack][ixtrack], x1);
1072 wave->y[itrack - track0] =
1073 wave->y[itrack - track0 - 1] + DIST(x0, x1);
1074 }
1075 }
1076
1077 /* Save geolocation... */
1078 wave->time = pert->time[(track0 + track1) / 2][(xtrack0 + xtrack1) / 2];
1079 wave->z = 0;
1080 wave->lon[ixtrack - xtrack0][itrack - track0] =
1081 pert->lon[itrack][ixtrack];
1082 wave->lat[ixtrack - xtrack0][itrack - track0] =
1083 pert->lat[itrack][ixtrack];
1084
1085 /* Save temperature data... */
1086 wave->temp[ixtrack - xtrack0][itrack - track0]
1087 = pert->bt[itrack][ixtrack];
1088 wave->bg[ixtrack - xtrack0][itrack - track0]
1089 = pert->bt[itrack][ixtrack] - pert->pt[itrack][ixtrack];
1090 wave->pt[ixtrack - xtrack0][itrack - track0]
1091 = pert->pt[itrack][ixtrack];
1092 wave->var[ixtrack - xtrack0][itrack - track0]
1093 = pert->var[itrack][ixtrack];
1094 }
1095}
1096#endif
1097
1098/*****************************************************************************/
1099
1100int read_cris_l1(
1101 char *filename,
1102 cris_l1_t *l1,
1103 int apo) {
1104
1105 int ncid, varid, dimid;
1106
1107 size_t n;
1108
1109 /* Write info... */
1110 LOG(1, "Read CrIS Level-1B file: %s", filename);
1111
1112 /* Open netCDF file... */
1113 if (nc_open(filename, NC_NOWRITE, &ncid) != NC_NOERR) {
1114 WARN("Cannot open file!");
1115 return 0;
1116 }
1117
1118 /* Check dimensions... */
1119 NC(nc_inq_dimid(ncid, "atrack", &dimid));
1120 NC(nc_inq_dimlen(ncid, dimid, &n));
1121 if (n != L1_NTRACK)
1122 ERRMSG("Level-1B file dimension mismatch!");
1123
1124 NC(nc_inq_dimid(ncid, "xtrack", &dimid));
1125 NC(nc_inq_dimlen(ncid, dimid, &n));
1126 if (n != L1_NXTRACK)
1127 ERRMSG("Level-1B file dimension mismatch!");
1128
1129 NC(nc_inq_dimid(ncid, "fov", &dimid));
1130 NC(nc_inq_dimlen(ncid, dimid, &n));
1131 if (n != L1_NFOV)
1132 ERRMSG("Level-1B file dimension mismatch!");
1133
1134 NC(nc_inq_dimid(ncid, "wnum_lw", &dimid));
1135 NC(nc_inq_dimlen(ncid, dimid, &n));
1136 if (n != L1_NCHAN_LW)
1137 ERRMSG("Level-1B file dimension mismatch!");
1138
1139 NC(nc_inq_dimid(ncid, "wnum_mw", &dimid));
1140 NC(nc_inq_dimlen(ncid, dimid, &n));
1141 if (n != L1_NCHAN_MW)
1142 ERRMSG("Level-1B file dimension mismatch!");
1143
1144 NC(nc_inq_dimid(ncid, "wnum_sw", &dimid));
1145 NC(nc_inq_dimlen(ncid, dimid, &n));
1146 if (n != L1_NCHAN_SW)
1147 ERRMSG("Level-1B file dimension mismatch!");
1148
1149 /* Read satellite position... */
1150 NC(nc_inq_varid(ncid, "obs_time_tai93", &varid));
1151 NC(nc_get_var_double(ncid, varid, l1->time[0]));
1152 NC(nc_inq_varid(ncid, "lon", &varid));
1153 NC(nc_get_var_double(ncid, varid, l1->lon[0][0]));
1154 NC(nc_inq_varid(ncid, "lat", &varid));
1155 NC(nc_get_var_double(ncid, varid, l1->lat[0][0]));
1156 NC(nc_inq_varid(ncid, "sat_alt", &varid));
1157 NC(nc_get_var_double(ncid, varid, l1->sat_z));
1158 NC(nc_inq_varid(ncid, "subsat_lon", &varid));
1159 NC(nc_get_var_double(ncid, varid, l1->sat_lon));
1160 NC(nc_inq_varid(ncid, "subsat_lat", &varid));
1161 NC(nc_get_var_double(ncid, varid, l1->sat_lat));
1162
1163 /* Convert units... */
1164 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1165 l1->sat_z[itrack] /= 1e3;
1166
1167 /* Read spectra... */
1168 NC(nc_inq_varid(ncid, "wnum_lw", &varid));
1169 NC(nc_get_var_double(ncid, varid, l1->nu_lw));
1170 NC(nc_inq_varid(ncid, "rad_lw", &varid));
1171 NC(nc_get_var_float(ncid, varid, l1->rad_lw[0][0][0]));
1172 NC(nc_inq_varid(ncid, "wnum_mw", &varid));
1173 NC(nc_get_var_double(ncid, varid, l1->nu_mw));
1174 NC(nc_inq_varid(ncid, "rad_mw", &varid));
1175 NC(nc_get_var_float(ncid, varid, l1->rad_mw[0][0][0]));
1176 NC(nc_inq_varid(ncid, "wnum_sw", &varid));
1177 NC(nc_get_var_double(ncid, varid, l1->nu_sw));
1178 NC(nc_inq_varid(ncid, "rad_sw", &varid));
1179 NC(nc_get_var_float(ncid, varid, l1->rad_sw[0][0][0]));
1180
1181 /* Read noise... */
1182 NC(nc_inq_varid(ncid, "nedn_lw", &varid));
1183 NC(nc_get_var_float(ncid, varid, l1->nedn_lw[0]));
1184 NC(nc_inq_varid(ncid, "nedn_mw", &varid));
1185 NC(nc_get_var_float(ncid, varid, l1->nedn_mw[0]));
1186 NC(nc_inq_varid(ncid, "nedn_sw", &varid));
1187 NC(nc_get_var_float(ncid, varid, l1->nedn_sw[0]));
1188
1189 /* Check quality flags... */
1190 NC(nc_inq_varid(ncid, "rad_lw_qc", &varid));
1191 NC(nc_get_var_short(ncid, varid, l1->qual_lw[0][0]));
1192 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1193 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1194 for (int ifov = 0; ifov < L1_NFOV; ifov++)
1195 if (l1->qual_lw[itrack][ixtrack][ifov] > 1)
1196 for (int ichan = 0; ichan < L1_NCHAN_LW; ichan++)
1197 l1->rad_lw[itrack][ixtrack][ifov][ichan] = GSL_NAN;
1198
1199 NC(nc_inq_varid(ncid, "rad_mw_qc", &varid));
1200 NC(nc_get_var_short(ncid, varid, l1->qual_mw[0][0]));
1201 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1202 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1203 for (int ifov = 0; ifov < L1_NFOV; ifov++)
1204 if (l1->qual_mw[itrack][ixtrack][ifov] > 1)
1205 for (int ichan = 0; ichan < L1_NCHAN_MW; ichan++)
1206 l1->rad_mw[itrack][ixtrack][ifov][ichan] = GSL_NAN;
1207
1208 NC(nc_inq_varid(ncid, "rad_sw_qc", &varid));
1209 NC(nc_get_var_short(ncid, varid, l1->qual_sw[0][0]));
1210 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1211 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1212 for (int ifov = 0; ifov < L1_NFOV; ifov++)
1213 if (l1->qual_sw[itrack][ixtrack][ifov] > 1)
1214 for (int ichan = 0; ichan < L1_NCHAN_SW; ichan++)
1215 l1->rad_sw[itrack][ixtrack][ifov][ichan] = GSL_NAN;
1216
1217 /* Apodization... */
1218 if (apo)
1219 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1220 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1221 for (int ifov = 0; ifov < L1_NFOV; ifov++) {
1222
1223 double help[L1_NCHAN_MW];
1224 for (int ichan = 0; ichan < L1_NCHAN_LW; ichan++)
1225 help[ichan] = 0;
1226 for (int ichan = 1; ichan < L1_NCHAN_LW - 1; ichan++)
1227 help[ichan]
1228 = 0.23 * l1->rad_lw[itrack][ixtrack][ifov][ichan - 1]
1229 + 0.54 * l1->rad_lw[itrack][ixtrack][ifov][ichan]
1230 + 0.23 * l1->rad_lw[itrack][ixtrack][ifov][ichan + 1];
1231 for (int ichan = 1; ichan < L1_NCHAN_LW - 1; ichan++)
1232 l1->rad_lw[itrack][ixtrack][ifov][ichan] = (float) help[ichan];
1233 l1->rad_lw[itrack][ixtrack][ifov][0] =
1234 l1->rad_lw[itrack][ixtrack][ifov][L1_NCHAN_LW - 1] = GSL_NAN;
1235
1236 for (int ichan = 0; ichan < L1_NCHAN_MW; ichan++)
1237 help[ichan] = 0;
1238 for (int ichan = 1; ichan < L1_NCHAN_MW - 1; ichan++)
1239 help[ichan]
1240 = 0.23 * l1->rad_mw[itrack][ixtrack][ifov][ichan - 1]
1241 + 0.54 * l1->rad_mw[itrack][ixtrack][ifov][ichan]
1242 + 0.23 * l1->rad_mw[itrack][ixtrack][ifov][ichan + 1];
1243 for (int ichan = 1; ichan < L1_NCHAN_MW - 1; ichan++)
1244 l1->rad_mw[itrack][ixtrack][ifov][ichan] = (float) help[ichan];
1245 l1->rad_mw[itrack][ixtrack][ifov][0] =
1246 l1->rad_mw[itrack][ixtrack][ifov][L1_NCHAN_MW - 1] = GSL_NAN;
1247
1248 for (int ichan = 0; ichan < L1_NCHAN_SW; ichan++)
1249 help[ichan] = 0;
1250 for (int ichan = 1; ichan < L1_NCHAN_SW - 1; ichan++)
1251 help[ichan]
1252 = 0.23 * l1->rad_sw[itrack][ixtrack][ifov][ichan - 1]
1253 + 0.54 * l1->rad_sw[itrack][ixtrack][ifov][ichan]
1254 + 0.23 * l1->rad_sw[itrack][ixtrack][ifov][ichan + 1];
1255 for (int ichan = 1; ichan < L1_NCHAN_SW - 1; ichan++)
1256 l1->rad_sw[itrack][ixtrack][ifov][ichan] = (float) help[ichan];
1257 l1->rad_sw[itrack][ixtrack][ifov][0] =
1258 l1->rad_sw[itrack][ixtrack][ifov][L1_NCHAN_SW - 1] = GSL_NAN;
1259 }
1260
1261 /* Close file... */
1262 NC(nc_close(ncid));
1263
1264 /* Return success... */
1265 return 1;
1266}
1267
1268/*****************************************************************************/
1269
1270void read_pert(
1271 char *filename,
1272 char *pertname,
1273 int dc,
1274 pert_t *pert) {
1275
1276 static char varname[LEN];
1277
1278 static double help[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV];
1279
1280 static int dimid[3], ncid, varid;
1281
1282 static size_t ntrack, nxtrack, nfov, start[3] = { 0, 0, 0 }, count[3] =
1283 { 1, 1, 1 };
1284
1285 /* Write info... */
1286 LOG(1, "Read perturbation data: %s", filename);
1287
1288 /* Open netCDF file... */
1289 NC(nc_open(filename, NC_NOWRITE, &ncid));
1290
1291 /* Get dimensions... */
1292 NC(nc_inq_dimid(ncid, "NTRACK", &dimid[0]));
1293 NC(nc_inq_dimlen(ncid, dimid[0], &ntrack));
1294 if (ntrack > PERT_NTRACK)
1295 ERRMSG("Too many scans!");
1296
1297 NC(nc_inq_dimid(ncid, "NXTRACK", &dimid[1]));
1298 NC(nc_inq_dimlen(ncid, dimid[1], &nxtrack));
1299 if (nxtrack > PERT_NXTRACK)
1300 ERRMSG("Too many tracks!");
1301
1302 NC(nc_inq_dimid(ncid, "NFOV", &dimid[2]));
1303 NC(nc_inq_dimlen(ncid, dimid[2], &nfov));
1304 if (nfov > PERT_NFOV)
1305 ERRMSG("Too many field of views!");
1306
1307 pert->ntrack = (int) ntrack;
1308 pert->nxtrack = (int) nxtrack;
1309 pert->nfov = (int) nfov;
1310 count[2] = nfov;
1311
1312 /* Read data... */
1313 NC(nc_inq_varid(ncid, "time", &varid));
1314 for (size_t itrack = 0; itrack < ntrack; itrack++)
1315 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1316 start[0] = itrack;
1317 start[1] = ixtrack;
1318 NC(nc_get_vara_double
1319 (ncid, varid, start, count, pert->time[itrack][ixtrack]));
1320 }
1321
1322 NC(nc_inq_varid(ncid, "lon", &varid));
1323 for (size_t itrack = 0; itrack < ntrack; itrack++)
1324 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1325 start[0] = itrack;
1326 start[1] = ixtrack;
1327 NC(nc_get_vara_double
1328 (ncid, varid, start, count, pert->lon[itrack][ixtrack]));
1329 }
1330
1331 NC(nc_inq_varid(ncid, "lat", &varid));
1332 for (size_t itrack = 0; itrack < ntrack; itrack++)
1333 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1334 start[0] = itrack;
1335 start[1] = ixtrack;
1336 NC(nc_get_vara_double
1337 (ncid, varid, start, count, pert->lat[itrack][ixtrack]));
1338 }
1339
1340 NC(nc_inq_varid(ncid, "bt_8mu", &varid));
1341 for (size_t itrack = 0; itrack < ntrack; itrack++)
1342 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1343 start[0] = itrack;
1344 start[1] = ixtrack;
1345 NC(nc_get_vara_double
1346 (ncid, varid, start, count, pert->dc[itrack][ixtrack]));
1347 }
1348
1349 NC(nc_inq_varid(ncid, "bt_10mu", &varid));
1350 for (size_t itrack = 0; itrack < ntrack; itrack++)
1351 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1352 start[0] = itrack;
1353 start[1] = ixtrack;
1354 NC(nc_get_vara_double
1355 (ncid, varid, start, count, help[itrack][ixtrack]));
1356 }
1357
1358 for (size_t itrack = 0; itrack < ntrack; itrack++)
1359 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++)
1360 for (size_t ifov = 0; ifov < nfov; ifov++)
1361 if (dc == 2
1362 || (dc == 0 && !gsl_finite(pert->dc[itrack][ixtrack][ifov])))
1363 pert->dc[itrack][ixtrack][ifov] = help[itrack][ixtrack][ifov];
1364
1365 sprintf(varname, "bt_%s", pertname);
1366 NC(nc_inq_varid(ncid, varname, &varid));
1367 for (size_t itrack = 0; itrack < ntrack; itrack++)
1368 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1369 start[0] = itrack;
1370 start[1] = ixtrack;
1371 NC(nc_get_vara_double
1372 (ncid, varid, start, count, pert->bt[itrack][ixtrack]));
1373 }
1374
1375 sprintf(varname, "bt_%s_pt", pertname);
1376 NC(nc_inq_varid(ncid, varname, &varid));
1377 for (size_t itrack = 0; itrack < ntrack; itrack++)
1378 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1379 start[0] = itrack;
1380 start[1] = ixtrack;
1381 NC(nc_get_vara_double
1382 (ncid, varid, start, count, pert->pt[itrack][ixtrack]));
1383 }
1384
1385 sprintf(varname, "bt_%s_var", pertname);
1386 NC(nc_inq_varid(ncid, varname, &varid));
1387 for (size_t itrack = 0; itrack < ntrack; itrack++)
1388 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1389 start[0] = itrack;
1390 start[1] = ixtrack;
1391 NC(nc_get_vara_double
1392 (ncid, varid, start, count, pert->var[itrack][ixtrack]));
1393 }
1394
1395 /* Close file... */
1396 NC(nc_close(ncid));
1397}
1398
1399/*****************************************************************************/
1400
1401void read_retr(
1402 char *filename,
1403 ret_t *ret) {
1404
1405 static double help[NDS * NPG];
1406
1407 int dimid, ids, ncid, varid;
1408
1409 size_t nds, np, ntrack, nxtrack;
1410
1411 /* Write info... */
1412 LOG(1, "Read retrieval data: %s", filename);
1413
1414 /* Open netCDF file... */
1415 NC(nc_open(filename, NC_NOWRITE, &ncid));
1416
1417 /* Read new retrieval file format... */
1418 if (nc_inq_dimid(ncid, "L1_NTRACK", &dimid) == NC_NOERR) {
1419
1420 /* Get dimensions... */
1421 NC(nc_inq_dimid(ncid, "RET_NP", &dimid));
1422 NC(nc_inq_dimlen(ncid, dimid, &np));
1423 ret->np = (int) np;
1424 if (ret->np > NPG)
1425 ERRMSG("Too many data points!");
1426
1427 NC(nc_inq_dimid(ncid, "L1_NTRACK", &dimid));
1428 NC(nc_inq_dimlen(ncid, dimid, &ntrack));
1429 NC(nc_inq_dimid(ncid, "L1_NXTRACK", &dimid));
1430 NC(nc_inq_dimlen(ncid, dimid, &nxtrack));
1431 ret->nds = (int) (ntrack * nxtrack);
1432 if (ret->nds > NDS)
1433 ERRMSG("Too many data sets!");
1434
1435 /* Read time... */
1436 NC(nc_inq_varid(ncid, "l1_time", &varid));
1437 NC(nc_get_var_double(ncid, varid, help));
1438 ids = 0;
1439 for (size_t itrack = 0; itrack < ntrack; itrack++)
1440 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1441 for (int ip = 0; ip < ret->np; ip++)
1442 ret->time[ids][ip] = help[ids];
1443 ids++;
1444 }
1445
1446 /* Read altitudes... */
1447 NC(nc_inq_varid(ncid, "ret_z", &varid));
1448 NC(nc_get_var_double(ncid, varid, help));
1449 ids = 0;
1450 for (size_t itrack = 0; itrack < ntrack; itrack++)
1451 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1452 for (int ip = 0; ip < ret->np; ip++)
1453 ret->z[ids][ip] = help[ip];
1454 ids++;
1455 }
1456
1457 /* Read longitudes... */
1458 NC(nc_inq_varid(ncid, "l1_lon", &varid));
1459 NC(nc_get_var_double(ncid, varid, help));
1460 ids = 0;
1461 for (size_t itrack = 0; itrack < ntrack; itrack++)
1462 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1463 for (int ip = 0; ip < ret->np; ip++)
1464 ret->lon[ids][ip] = help[ids];
1465 ids++;
1466 }
1467
1468 /* Read latitudes... */
1469 NC(nc_inq_varid(ncid, "l1_lat", &varid));
1470 NC(nc_get_var_double(ncid, varid, help));
1471 ids = 0;
1472 for (size_t itrack = 0; itrack < ntrack; itrack++)
1473 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1474 for (int ip = 0; ip < ret->np; ip++)
1475 ret->lat[ids][ip] = help[ids];
1476 ids++;
1477 }
1478
1479 /* Read temperatures... */
1480 NC(nc_inq_varid(ncid, "ret_temp", &varid));
1481 NC(nc_get_var_double(ncid, varid, help));
1482 ids = 0;
1483 for (size_t itrack = 0; itrack < ntrack; itrack++)
1484 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1485 for (int ip = 0; ip < ret->np; ip++)
1486 ret->t[ids][ip] =
1487 help[(itrack * nxtrack + ixtrack) * (size_t) np + (size_t) ip];
1488 ids++;
1489 }
1490 }
1491
1492 /* Read old retrieval file format... */
1493 if (nc_inq_dimid(ncid, "np", &dimid) == NC_NOERR) {
1494
1495 /* Get dimensions... */
1496 NC(nc_inq_dimid(ncid, "np", &dimid));
1497 NC(nc_inq_dimlen(ncid, dimid, &np));
1498 ret->np = (int) np;
1499 if (ret->np > NPG)
1500 ERRMSG("Too many data points!");
1501
1502 NC(nc_inq_dimid(ncid, "nds", &dimid));
1503 NC(nc_inq_dimlen(ncid, dimid, &nds));
1504 ret->nds = (int) nds;
1505 if (ret->nds > NDS)
1506 ERRMSG("Too many data sets!");
1507
1508 /* Read data... */
1509 NC(nc_inq_varid(ncid, "time", &varid));
1510 NC(nc_get_var_double(ncid, varid, help));
1511 read_retr_help(help, ret->nds, ret->np, ret->time);
1512
1513 NC(nc_inq_varid(ncid, "z", &varid));
1514 NC(nc_get_var_double(ncid, varid, help));
1515 read_retr_help(help, ret->nds, ret->np, ret->z);
1516
1517 NC(nc_inq_varid(ncid, "lon", &varid));
1518 NC(nc_get_var_double(ncid, varid, help));
1519 read_retr_help(help, ret->nds, ret->np, ret->lon);
1520
1521 NC(nc_inq_varid(ncid, "lat", &varid));
1522 NC(nc_get_var_double(ncid, varid, help));
1523 read_retr_help(help, ret->nds, ret->np, ret->lat);
1524
1525 NC(nc_inq_varid(ncid, "press", &varid));
1526 NC(nc_get_var_double(ncid, varid, help));
1527 read_retr_help(help, ret->nds, ret->np, ret->p);
1528
1529 NC(nc_inq_varid(ncid, "temp", &varid));
1530 NC(nc_get_var_double(ncid, varid, help));
1531 read_retr_help(help, ret->nds, ret->np, ret->t);
1532
1533 NC(nc_inq_varid(ncid, "temp_apr", &varid));
1534 NC(nc_get_var_double(ncid, varid, help));
1535 read_retr_help(help, ret->nds, ret->np, ret->t_apr);
1536
1537 NC(nc_inq_varid(ncid, "temp_total", &varid));
1538 NC(nc_get_var_double(ncid, varid, help));
1539 read_retr_help(help, ret->nds, ret->np, ret->t_tot);
1540
1541 NC(nc_inq_varid(ncid, "temp_noise", &varid));
1542 NC(nc_get_var_double(ncid, varid, help));
1543 read_retr_help(help, ret->nds, ret->np, ret->t_noise);
1544
1545 NC(nc_inq_varid(ncid, "temp_formod", &varid));
1546 NC(nc_get_var_double(ncid, varid, help));
1547 read_retr_help(help, ret->nds, ret->np, ret->t_fm);
1548
1549 NC(nc_inq_varid(ncid, "temp_cont", &varid));
1550 NC(nc_get_var_double(ncid, varid, help));
1551 read_retr_help(help, ret->nds, ret->np, ret->t_cont);
1552
1553 NC(nc_inq_varid(ncid, "temp_res", &varid));
1554 NC(nc_get_var_double(ncid, varid, help));
1555 read_retr_help(help, ret->nds, ret->np, ret->t_res);
1556
1557 NC(nc_inq_varid(ncid, "chisq", &varid));
1558 NC(nc_get_var_double(ncid, varid, ret->chisq));
1559 }
1560
1561 /* Close file... */
1562 NC(nc_close(ncid));
1563}
1564
1565/*****************************************************************************/
1566
1567void read_retr_help(
1568 double *help,
1569 int nds,
1570 int np,
1571 double mat[NDS][NPG]) {
1572
1573 int n = 0;
1574
1575 for (int ip = 0; ip < np; ip++)
1576 for (int ids = 0; ids < nds; ids++)
1577 mat[ids][ip] = help[n++];
1578}
1579
1580/*****************************************************************************/
1581
1582void read_wave(
1583 char *filename,
1584 wave_t *wave) {
1585
1586 FILE *in;
1587
1588 char line[LEN];
1589
1590 double rtime, rz, rlon, rlat, rx, ry, ryold = -1e10,
1591 rtemp, rbg, rpt, rvar, rfit;
1592
1593 /* Init... */
1594 wave->nx = 0;
1595 wave->ny = 0;
1596
1597 /* Write info... */
1598 LOG(1, "Read wave data: %s", filename);
1599
1600 /* Open file... */
1601 if (!(in = fopen(filename, "r")))
1602 ERRMSG("Cannot open file!");
1603
1604 /* Read data... */
1605 while (fgets(line, LEN, in))
1606 if (sscanf(line, "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", &rtime,
1607 &rz, &rlon, &rlat, &rx, &ry, &rtemp, &rbg, &rpt,
1608 &rvar, &rfit) == 11) {
1609
1610 /* Set index... */
1611 if (ry != ryold) {
1612 if ((++wave->ny >= WY))
1613 ERRMSG("Too many y-values!");
1614 wave->nx = 0;
1615 } else if ((++wave->nx) >= WX)
1616 ERRMSG("Too many x-values!");
1617 ryold = ry;
1618
1619 /* Save data... */
1620 wave->time = rtime;
1621 wave->z = rz;
1622 wave->lon[wave->nx][wave->ny] = rlon;
1623 wave->lat[wave->nx][wave->ny] = rlat;
1624 wave->x[wave->nx] = rx;
1625 wave->y[wave->ny] = ry;
1626 wave->temp[wave->nx][wave->ny] = rtemp;
1627 wave->bg[wave->nx][wave->ny] = rbg;
1628 wave->pt[wave->nx][wave->ny] = rpt;
1629 wave->var[wave->nx][wave->ny] = rvar;
1630 wave->var[wave->nx][wave->ny] = rfit;
1631 }
1632
1633 /* Increment counters... */
1634 wave->nx++;
1635 wave->ny++;
1636
1637 /* Close file... */
1638 fclose(in);
1639}
1640
1641/*****************************************************************************/
1642
1643#if 0
1644void rad2wave(
1645 cris_l1_t *gran,
1646 double *nu,
1647 int nd,
1648 wave_t *wave) {
1649
1650 double x0[3], x1[3];
1651
1652 int ichan[AIRS_RAD_CHANNEL];
1653
1654 /* Get channel numbers... */
1655 for (int id = 0; id < nd; id++) {
1656 for (ichan[id] = 0; ichan[id] < AIRS_RAD_CHANNEL; ichan[id]++)
1657 if (fabs(gran->nominal_freq[ichan[id]] - nu[id]) < 0.1)
1658 break;
1659 if (ichan[id] >= AIRS_RAD_CHANNEL)
1660 ERRMSG("Could not find channel!");
1661 }
1662
1663 /* Set size... */
1664 wave->nx = AIRS_RAD_GEOXTRACK;
1665 wave->ny = AIRS_RAD_GEOTRACK;
1666 if (wave->nx > WX || wave->ny > WY)
1667 ERRMSG("Wave struct too small!");
1668
1669 /* Set Cartesian coordinates... */
1670 geo2cart(0, gran->Longitude[0][0], gran->Latitude[0][0], x0);
1671 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
1672 geo2cart(0, gran->Longitude[0][xtrack], gran->Latitude[0][xtrack], x1);
1673 wave->x[xtrack] = DIST(x0, x1);
1674 }
1675 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++) {
1676 geo2cart(0, gran->Longitude[track][0], gran->Latitude[track][0], x1);
1677 wave->y[track] = DIST(x0, x1);
1678 }
1679
1680 /* Set geolocation... */
1681 wave->time =
1682 gran->Time[AIRS_RAD_GEOTRACK / 2][AIRS_RAD_GEOXTRACK / 2] - 220838400;
1683 wave->z = 0;
1684 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
1685 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
1686 wave->lon[xtrack][track] = gran->Longitude[track][xtrack];
1687 wave->lat[xtrack][track] = gran->Latitude[track][xtrack];
1688 }
1689
1690 /* Set brightness temperature... */
1691 for (int track = 0; track < AIRS_RAD_GEOTRACK; track++)
1692 for (int xtrack = 0; xtrack < AIRS_RAD_GEOXTRACK; xtrack++) {
1693 wave->temp[xtrack][track] = 0;
1694 wave->bg[xtrack][track] = 0;
1695 wave->pt[xtrack][track] = 0;
1696 wave->var[xtrack][track] = 0;
1697 for (int id = 0; id < nd; id++) {
1698 if ((gran->state[track][xtrack] != 0)
1699 || (gran->ExcludedChans[ichan[id]] > 2)
1700 || (gran->CalChanSummary[ichan[id]] & 8)
1701 || (gran->CalChanSummary[ichan[id]] & (32 + 64))
1702 || (gran->CalFlag[track][ichan[id]] & 16))
1703 wave->temp[xtrack][track] = GSL_NAN;
1704 else
1705 wave->temp[xtrack][track]
1706 += brightness(gran->radiances[track][xtrack][ichan[id]] * 1e-3,
1707 gran->nominal_freq[ichan[id]]) / nd;
1708 }
1709 }
1710}
1711#endif
1712
1713/*****************************************************************************/
1714
1715void ret2wave(
1716 ret_t *ret,
1717 wave_t *wave,
1718 int dataset,
1719 int ip) {
1720
1721 double x0[3], x1[3];
1722
1723 /* Initialize... */
1724 wave->nx = 90;
1725 if (wave->nx > WX)
1726 ERRMSG("Too many across-track values!");
1727 wave->ny = 135;
1728 if (wave->ny > WY)
1729 ERRMSG("Too many along-track values!");
1730 if (ip < 0 || ip >= ret->np)
1731 ERRMSG("Altitude index out of range!");
1732
1733 /* Loop over data sets and data points... */
1734 for (int ids = 0; ids < ret->nds; ids++) {
1735
1736 /* Get horizontal indices... */
1737 int ix = ids % 90;
1738 int iy = ids / 90;
1739
1740 /* Get distances... */
1741 if (iy == 0) {
1742 geo2cart(0.0, ret->lon[0][0], ret->lat[0][0], x0);
1743 geo2cart(0.0, ret->lon[ids][ip], ret->lat[ids][ip], x1);
1744 wave->x[ix] = DIST(x0, x1);
1745 }
1746 if (ix == 0) {
1747 geo2cart(0.0, ret->lon[0][0], ret->lat[0][0], x0);
1748 geo2cart(0.0, ret->lon[ids][ip], ret->lat[ids][ip], x1);
1749 wave->y[iy] = DIST(x0, x1);
1750 }
1751
1752 /* Save geolocation... */
1753 wave->time = ret->time[0][0];
1754 if (ix == 0 && iy == 0)
1755 wave->z = ret->z[ids][ip];
1756 wave->lon[ix][iy] = ret->lon[ids][ip];
1757 wave->lat[ix][iy] = ret->lat[ids][ip];
1758
1759 /* Save temperature... */
1760 if (dataset == 1)
1761 wave->temp[ix][iy] = ret->t[ids][ip];
1762 else if (dataset == 2)
1763 wave->temp[ix][iy] = ret->t_apr[ids][ip];
1764 }
1765}
1766
1767/*****************************************************************************/
1768
1769void variance(
1770 wave_t *wave,
1771 double dh) {
1772
1773 /* Check parameters... */
1774 if (dh <= 0)
1775 return;
1776
1777 /* Compute squared radius... */
1778 const double dh2 = gsl_pow_2(dh);
1779
1780 /* Get sampling distances... */
1781 const int dx =
1782 (int) (dh / fabs(wave->x[wave->nx - 1] - wave->x[0]) * (wave->nx - 1.0) +
1783 1);
1784 const int dy =
1785 (int) (dh / fabs(wave->y[wave->ny - 1] - wave->y[0]) * (wave->ny - 1.0) +
1786 1);
1787
1788 /* Loop over data points... */
1789 for (int ix = 0; ix < wave->nx; ix++)
1790 for (int iy = 0; iy < wave->ny; iy++) {
1791
1792 /* Init... */
1793 double mu = 0, help = 0;
1794 int n = 0;
1795
1796 /* Get data... */
1797 for (int ix2 = GSL_MAX(ix - dx, 0);
1798 ix2 <= GSL_MIN(ix + dx, wave->nx - 1); ix2++)
1799 for (int iy2 = GSL_MAX(iy - dy, 0);
1800 iy2 <= GSL_MIN(iy + dy, wave->ny - 1); iy2++)
1801 if ((gsl_pow_2(wave->x[ix] - wave->x[ix2])
1802 + gsl_pow_2(wave->y[iy] - wave->y[iy2])) <= dh2)
1803 if (gsl_finite(wave->pt[ix2][iy2])) {
1804 mu += wave->pt[ix2][iy2];
1805 help += gsl_pow_2(wave->pt[ix2][iy2]);
1806 n++;
1807 }
1808
1809 /* Compute local variance... */
1810 if (n > 1)
1811 wave->var[ix][iy] = help / n - gsl_pow_2(mu / n);
1812 else
1813 wave->var[ix][iy] = GSL_NAN;
1814 }
1815}
1816
1817/*****************************************************************************/
1818
1819void write_wave(
1820 char *filename,
1821 wave_t *wave) {
1822
1823 FILE *out;
1824
1825 /* Write info... */
1826 LOG(1, "Write wave data: %s", filename);
1827
1828 /* Create file... */
1829 if (!(out = fopen(filename, "w")))
1830 ERRMSG("Cannot create file!");
1831
1832 /* Write header... */
1833 fprintf(out,
1834 "# $1 = time (seconds since 2000-01-01T00:00Z)\n"
1835 "# $2 = altitude [km]\n"
1836 "# $3 = longitude [deg]\n"
1837 "# $4 = latitude [deg]\n"
1838 "# $5 = across-track distance [km]\n"
1839 "# $6 = along-track distance [km]\n"
1840 "# $7 = temperature [K]\n"
1841 "# $8 = background [K]\n"
1842 "# $9 = perturbation [K]\n"
1843 "# $10 = variance [K^2]\n" "# $11 = fitting model [K]\n");
1844
1845 /* Write data... */
1846 for (int j = 0; j < wave->ny; j++) {
1847 fprintf(out, "\n");
1848 for (int i = 0; i < wave->nx; i++)
1849 fprintf(out, "%.2f %g %g %g %g %g %g %g %g %g %g\n",
1850 wave->time, wave->z, wave->lon[i][j], wave->lat[i][j],
1851 wave->x[i], wave->y[j], wave->temp[i][j], wave->bg[i][j],
1852 wave->pt[i][j], wave->var[i][j], wave->fit[i][j]);
1853 }
1854
1855 /* Close file... */
1856 fclose(out);
1857}
cart2geo
void cart2geo(const double *x, double *z, double *lon, double *lat)
Convert Cartesian coordinates to geolocation.
Definition: jurassic.c:108
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Convert geolocation to Cartesian coordinates.
Definition: jurassic.c:3500
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
POW2
#define POW2(x)
Compute x^2.
Definition: jurassic.h:187
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
WARN
#define WARN(...)
Print warning message.
Definition: jurassic.h:231
LOG
#define LOG(level,...)
Print log message.
Definition: jurassic.h:221
DIST
#define DIST(a, b)
Compute Cartesian distance between two vectors.
Definition: jurassic.h:134
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:164
day2doy
void day2doy(const int year, const int mon, const int day, int *doy)
Get day of year from date.
Definition: libcris.c:279
merge_y
void merge_y(wave_t *wave1, wave_t *wave2)
Merge wave structs in y-direction.
Definition: libcris.c:743
fit_wave
void fit_wave(wave_t *wave, double amp, double phi, double kx, double ky, double *chisq)
Evaluate wave fit...
Definition: libcris.c:328
create_wave
void create_wave(wave_t *wave, double amp, double lx, double ly, double phi, double fwhm)
Add linear wave pattern...
Definition: libcris.c:250
noise_pert
void noise_pert(pert_t *pert, int track0, int track1, double *mu, double *sig)
Estimate noise from perurbations.
Definition: libcris.c:821
background_poly_help
void background_poly_help(const double *xx, double *yy, const int n, const int dim)
Get background based on polynomial fits.
Definition: libcris.c:47
fft_help
void fft_help(double *fcReal, double *fcImag, int n)
Calculate 1-D FFT...
Definition: libcris.c:354
background_smooth
void background_smooth(wave_t *wave, int npts_x, int npts_y)
Smooth background.
Definition: libcris.c:151
intpol_x
void intpol_x(wave_t *wave, int n)
Interpolate to regular grid in x-direction.
Definition: libcris.c:624
fft
void fft(wave_t *wave, double *Amax, double *phimax, double *lhmax, double *kxmax, double *kymax, double *alphamax, double *betamax, char *filename)
Calculate 2-D FFT...
Definition: libcris.c:393
noise
void noise(wave_t *wave, double *mu, double *sig)
Estimate noise.
Definition: libcris.c:776
read_retr_help
void read_retr_help(double *help, int nds, int np, double mat[NDS][NPG])
Convert array.
Definition: libcris.c:1567
hamming
void hamming(wave_t *wave, int niter)
Apply Hamming filter to perturbations...
Definition: libcris.c:595
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libcris.c:101
create_noise
void create_noise(wave_t *wave, double nedt)
Add noise to perturbations and temperatures...
Definition: libcris.c:229
read_cris_l1
int read_cris_l1(char *filename, cris_l1_t *l1, int apo)
Read CrIS Level-1 data.
Definition: libcris.c:1100
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libcris.c:1769
ret2wave
void ret2wave(ret_t *ret, wave_t *wave, int dataset, int ip)
Convert CrIS retrieval results to wave analysis struct.
Definition: libcris.c:1715
add_att
void add_att(const int ncid, const int varid, const char *unit, const char *long_name)
Add variable attributes to netCDF file.
Definition: libcris.c:5
read_retr
void read_retr(char *filename, ret_t *ret)
Read CrIS retrieval data.
Definition: libcris.c:1401
doy2day
void doy2day(const int year, const int doy, int *mon, int *day)
Get date from day of year.
Definition: libcris.c:298
create_background
void create_background(wave_t *wave)
Set background...
Definition: libcris.c:204
read_pert
void read_pert(char *filename, char *pertname, int dc, pert_t *pert)
Read radiance perturbation data.
Definition: libcris.c:1270
read_wave
void read_wave(char *filename, wave_t *wave)
Read wave analysis data.
Definition: libcris.c:1582
gauss
void gauss(wave_t *wave, double fwhm)
Apply Gaussian filter to perturbations...
Definition: libcris.c:555
add_var
void add_var(const int ncid, const char *varname, const char *unit, const char *longname, int type, int dimid[], int *varid, int ndims)
Add variable to netCDF file.
Definition: libcris.c:20
median
void median(wave_t *wave, int dx)
Apply median filter to perturbations...
Definition: libcris.c:702
write_wave
void write_wave(char *filename, wave_t *wave)
Write wave analysis data.
Definition: libcris.c:1819
period
void period(wave_t *wave, double lxymax, double dlxy, double *Amax, double *phimax, double *lhmax, double *kxmax, double *kymax, double *alphamax, double *betamax, char *filename)
Compute periodogram.
Definition: libcris.c:868
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libcris.h:61
rad2wave
void rad2wave(cris_l1_t *cris_l1, double *nu, int nd, wave_t *wave)
Convert CrIS radiance data to wave analysis struct.
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libcris.h:42
PMAX
#define PMAX
Maximum number of data points for spectral analysis.
Definition: libcris.h:54
L1_NCHAN_LW
#define L1_NCHAN_LW
Number of CrIS longwave radiance channels.
Definition: libcris.h:30
L1_NXTRACK
#define L1_NXTRACK
Across-track size of CrIS radiance granule.
Definition: libcris.h:24
L1_NTRACK
#define L1_NTRACK
Along-track size of CrIS radiance granule.
Definition: libcris.h:21
L1_NFOV
#define L1_NFOV
Number of field of views of CrIS radiance granule.
Definition: libcris.h:27
WX
#define WX
Across-track size of wave analysis data.
Definition: libcris.h:48
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libcris.h:39
NDS
#define NDS
Maximum number of data sets per granule.
Definition: libcris.h:15
pert2wave
void pert2wave(pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
Convert radiance perturbation data to wave analysis struct.
L1_NCHAN_MW
#define L1_NCHAN_MW
Number of CrIS midwave radiance channels.
Definition: libcris.h:33
L1_NCHAN_SW
#define L1_NCHAN_SW
Number of CrIS shortwave radiance channels.
Definition: libcris.h:36
NPG
#define NPG
Maximum number of data points per granule.
Definition: libcris.h:18
PERT_NFOV
#define PERT_NFOV
Number of field of views of perturbation data.
Definition: libcris.h:45
WY
#define WY
Along-track size of wave analysis data.
Definition: libcris.h:51
cris_l1_t
CrIS Level-1 data.
Definition: libcris.h:72
cris_l1_t::nu_sw
double nu_sw[L1_NCHAN_SW]
Shortwave channel frequencies [cm^-1].
Definition: libcris.h:117
cris_l1_t::sat_z
double sat_z[L1_NTRACK]
Satellite altitude [km].
Definition: libcris.h:84
cris_l1_t::nedn_sw
float nedn_sw[L1_NFOV][L1_NCHAN_SW]
Longwave radiance noise [W/(m^2 sr cm^-1)].
Definition: libcris.h:123
cris_l1_t::rad_mw
float rad_mw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_MW]
Midwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:108
cris_l1_t::nedn_mw
float nedn_mw[L1_NFOV][L1_NCHAN_MW]
Midwave radiance noise [W/(m^2 sr cm^-1)].
Definition: libcris.h:111
cris_l1_t::nu_mw
double nu_mw[L1_NCHAN_MW]
Midwave channel frequencies [cm^-1].
Definition: libcris.h:105
cris_l1_t::rad_sw
float rad_sw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_SW]
Shortwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:120
cris_l1_t::nedn_lw
float nedn_lw[L1_NFOV][L1_NCHAN_LW]
Longwave radiance noise [W/(m^2 sr cm^-1)].
Definition: libcris.h:99
cris_l1_t::qual_mw
short qual_mw[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Midwave quality flag (0=best, 1=good, 2=don't use).
Definition: libcris.h:114
cris_l1_t::lon
double lon[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Footprint longitude [deg].
Definition: libcris.h:78
cris_l1_t::rad_lw
float rad_lw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_LW]
Longwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:96
cris_l1_t::sat_lat
double sat_lat[L1_NTRACK]
Satellite latitude [deg].
Definition: libcris.h:90
cris_l1_t::time
double time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:75
cris_l1_t::qual_sw
short qual_sw[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Shortwave quality flag (0=best, 1=good, 2=don't use).
Definition: libcris.h:126
cris_l1_t::sat_lon
double sat_lon[L1_NTRACK]
Satellite longitude [deg].
Definition: libcris.h:87
cris_l1_t::nu_lw
double nu_lw[L1_NCHAN_LW]
Longwave channel frequencies [cm^-1].
Definition: libcris.h:93
cris_l1_t::lat
double lat[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Footprint latitude [deg].
Definition: libcris.h:81
cris_l1_t::qual_lw
short qual_lw[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Longwave quality flag (0=best, 1=good, 2=don't use).
Definition: libcris.h:102
pert_t
Perturbation data.
Definition: libcris.h:131
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature (cloud channel) [K].
Definition: libcris.h:152
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Latitude [deg].
Definition: libcris.h:149
pert_t::nfov
int nfov
Number of field of views.
Definition: libcris.h:140
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libcris.h:158
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libcris.h:134
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libcris.h:161
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:143
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature (4.3 or 15 micron) [K].
Definition: libcris.h:155
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libcris.h:137
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Longitude [deg].
Definition: libcris.h:146
ret_t
Retrieval results.
Definition: libcris.h:166
ret_t::t_apr
double t_apr[NDS][NPG]
Temperature (a priori data) [K].
Definition: libcris.h:193
ret_t::chisq
double chisq[NDS]
Chi^2.
Definition: libcris.h:211
ret_t::t_noise
double t_noise[NDS][NPG]
Temperature (noise error) [K].
Definition: libcris.h:199
ret_t::lat
double lat[NDS][NPG]
Latitude [deg].
Definition: libcris.h:184
ret_t::t_cont
double t_cont[NDS][NPG]
Temperature (measurement content).
Definition: libcris.h:205
ret_t::t
double t[NDS][NPG]
Temperature [K].
Definition: libcris.h:190
ret_t::t_fm
double t_fm[NDS][NPG]
Temperature (forward model error) [K].
Definition: libcris.h:202
ret_t::p
double p[NDS][NPG]
Pressure [hPa].
Definition: libcris.h:187
ret_t::z
double z[NDS][NPG]
Altitude [km].
Definition: libcris.h:178
ret_t::t_res
double t_res[NDS][NPG]
Temperature (resolution).
Definition: libcris.h:208
ret_t::nds
int nds
Number of data sets.
Definition: libcris.h:169
ret_t::t_tot
double t_tot[NDS][NPG]
Temperature (total error) [K].
Definition: libcris.h:196
ret_t::np
int np
Number of data points.
Definition: libcris.h:172
ret_t::time
double time[NDS][NPG]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:175
ret_t::lon
double lon[NDS][NPG]
Longitude [deg].
Definition: libcris.h:181
wave_t
Wave analysis data.
Definition: libcris.h:216
wave_t::nx
int nx
Number of across-track values.
Definition: libcris.h:219
wave_t::ny
int ny
Number of along-track values.
Definition: libcris.h:222
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libcris.h:255
wave_t::temp
double temp[WX][WY]
Temperature [K].
Definition: libcris.h:243
wave_t::x
double x[WX]
Across-track distance [km].
Definition: libcris.h:237
wave_t::lon
double lon[WX][WY]
Longitude [deg].
Definition: libcris.h:231
wave_t::y
double y[WY]
Along-track distance [km].
Definition: libcris.h:240
wave_t::lat
double lat[WX][WY]
Latitude [deg].
Definition: libcris.h:234
wave_t::z
double z
Altitude [km].
Definition: libcris.h:228
wave_t::bg
double bg[WX][WY]
Background [K].
Definition: libcris.h:246
wave_t::fit
double fit[WX][WY]
Fit [K].
Definition: libcris.h:252
wave_t::time
double time
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:225
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libcris.h:249
Generated on Wed Feb 5 2025 20:00:06 for CrIS Code Collection by doxygen 1.9.4