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 fit_wave(
280 wave_t *wave,
281 double amp,
282 double phi,
283 double kx,
284 double ky,
285 double *chisq) {
286
287 /* Init... */
288 *chisq = 0;
289
290 /* Calculate fit... */
291 for (int ix = 0; ix < wave->nx; ix++)
292 for (int iy = 0; iy < wave->ny; iy++) {
293 wave->fit[ix][iy]
294 = amp * cos(kx * wave->x[ix] + ky * wave->y[iy]
295 - phi * M_PI / 180.);
296 *chisq += POW2(wave->fit[ix][iy] - wave->pt[ix][iy]);
297 }
298
299 /* Calculate chisq... */
300 *chisq /= (double) (wave->nx * wave->ny);
301}
302
303/*****************************************************************************/
304
305void fft_help(
306 double *fcReal,
307 double *fcImag,
308 int n) {
309
310 double data[2 * PMAX];
311
312 /* Check size... */
313 if (n > PMAX)
314 ERRMSG("Too many data points!");
315
316 /* Allocate... */
317 gsl_fft_complex_wavetable *wavetable =
318 gsl_fft_complex_wavetable_alloc((size_t) n);
319 gsl_fft_complex_workspace *workspace =
320 gsl_fft_complex_workspace_alloc((size_t) n);
321
322 /* Set data (real, complex)... */
323 for (int i = 0; i < n; i++) {
324 data[2 * i] = fcReal[i];
325 data[2 * i + 1] = fcImag[i];
326 }
327
328 /* Calculate FFT... */
329 gsl_fft_complex_forward(data, 1, (size_t) n, wavetable, workspace);
330
331 /* Copy data... */
332 for (int i = 0; i < n; i++) {
333 fcReal[i] = data[2 * i];
334 fcImag[i] = data[2 * i + 1];
335 }
336
337 /* Free... */
338 gsl_fft_complex_wavetable_free(wavetable);
339 gsl_fft_complex_workspace_free(workspace);
340}
341
342/*****************************************************************************/
343
344void fft(
345 wave_t *wave,
346 double *Amax,
347 double *phimax,
348 double *lhmax,
349 double *kxmax,
350 double *kymax,
351 double *alphamax,
352 double *betamax,
353 char *filename) {
354
355 static double A[PMAX][PMAX], phi[PMAX][PMAX], kx[PMAX], ky[PMAX],
356 cutReal[PMAX], cutImag[PMAX], boxImag[PMAX][PMAX], boxReal[PMAX][PMAX];
357
358 int imin, imax, jmin, jmax;
359
360 /* Find box... */
361 imin = jmin = 9999;
362 imax = jmax = -9999;
363 for (int i = 0; i < wave->nx; i++)
364 for (int j = 0; j < wave->ny; j++)
365 if (gsl_finite(wave->var[i][j])) {
366 imin = GSL_MIN(imin, i);
367 imax = GSL_MAX(imax, i);
368 jmin = GSL_MIN(jmin, j);
369 jmax = GSL_MAX(jmax, j);
370 }
371 const int nx = imax - imin + 1;
372 const int ny = jmax - jmin + 1;
373
374 /* Copy data... */
375 for (int i = imin; i <= imax; i++)
376 for (int j = jmin; j <= jmax; j++) {
377 if (gsl_finite(wave->pt[i][j]))
378 boxReal[i - imin][j - jmin] = wave->pt[i][j];
379 else
380 boxReal[i - imin][j - jmin] = 0.0;
381 boxImag[i - imin][j - jmin] = 0.0;
382 }
383
384 /* FFT of the rows... */
385 for (int i = 0; i < nx; i++) {
386 for (int j = 0; j < ny; j++) {
387 cutReal[j] = boxReal[i][j];
388 cutImag[j] = boxImag[i][j];
389 }
390 fft_help(cutReal, cutImag, ny);
391 for (int j = 0; j < ny; j++) {
392 boxReal[i][j] = cutReal[j];
393 boxImag[i][j] = cutImag[j];
394 }
395 }
396
397 /* FFT of the columns... */
398 for (int j = 0; j < ny; j++) {
399 for (int i = 0; i < nx; i++) {
400 cutReal[i] = boxReal[i][j];
401 cutImag[i] = boxImag[i][j];
402 }
403 fft_help(cutReal, cutImag, nx);
404 for (int i = 0; i < nx; i++) {
405 boxReal[i][j] = cutReal[i];
406 boxImag[i][j] = cutImag[i];
407 }
408 }
409
410 /* Get frequencies, amplitude, and phase... */
411 for (int i = 0; i < nx; i++)
412 kx[i] = 2. * M_PI * ((i < nx / 2) ? (double) i : -(double) (nx - i))
413 / (nx * fabs(wave->x[imax] - wave->x[imin]) / (nx - 1.0));
414 for (int j = 0; j < ny; j++)
415 ky[j] = 2. * M_PI * ((j < ny / 2) ? (double) j : -(double) (ny - j))
416 / (ny * fabs(wave->y[jmax] - wave->y[jmin]) / (ny - 1.0));
417 for (int i = 0; i < nx; i++)
418 for (int j = 0; j < ny; j++) {
419 A[i][j]
420 = (i == 0 && j == 0 ? 1.0 : 2.0) / (nx * ny)
421 * sqrt(gsl_pow_2(boxReal[i][j]) + gsl_pow_2(boxImag[i][j]));
422 phi[i][j]
423 = 180. / M_PI * atan2(boxImag[i][j], boxReal[i][j]);
424 }
425
426 /* Check frequencies... */
427 for (int i = 0; i < nx; i++)
428 for (int j = 0; j < ny; j++)
429 if (kx[i] == 0 || ky[j] == 0) {
430 A[i][j] = GSL_NAN;
431 phi[i][j] = GSL_NAN;
432 }
433
434 /* Find maximum... */
435 *Amax = 0;
436 for (int i = 0; i < nx; i++)
437 for (int j = 0; j < ny / 2; j++)
438 if (gsl_finite(A[i][j]) && A[i][j] > *Amax) {
439 *Amax = A[i][j];
440 *phimax = phi[i][j];
441 *kxmax = kx[i];
442 *kymax = ky[j];
443 }
444
445 /* Get horizontal wavelength... */
446 *lhmax = 2 * M_PI / sqrt(gsl_pow_2(*kxmax) + gsl_pow_2(*kymax));
447
448 /* Get propagation direction in xy-plane... */
449 *alphamax = 90. - 180. / M_PI * atan2(*kxmax, *kymax);
450
451 /* Get propagation direction in lon,lat-plane... */
452 *betamax = *alphamax
453 +
454 180. / M_PI *
455 atan2(wave->lat[wave->nx / 2 >
456 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
457 - wave->lat[wave->nx / 2 <
458 wave->nx - 1 ? wave->nx / 2 +
459 1 : wave->nx / 2][wave->ny / 2],
460 wave->lon[wave->nx / 2 >
461 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
462 - wave->lon[wave->nx / 2 <
463 wave->nx - 1 ? wave->nx / 2 +
464 1 : wave->nx / 2][wave->ny / 2]);
465
466 /* Save FFT data... */
467 if (filename != NULL) {
468
469 /* Write info... */
470 LOG(1, "Write FFT data: %s", filename);
471
472 /* Create file... */
473 FILE *out;
474 if (!(out = fopen(filename, "w")))
475 ERRMSG("Cannot create file!");
476
477 /* Write header... */
478 fprintf(out,
479 "# $1 = altitude [km]\n"
480 "# $2 = wavelength in x-direction [km]\n"
481 "# $3 = wavelength in y-direction [km]\n"
482 "# $4 = wavenumber in x-direction [1/km]\n"
483 "# $5 = wavenumber in y-direction [1/km]\n"
484 "# $6 = amplitude [K]\n" "# $7 = phase [rad]\n");
485
486 /* Write data... */
487 for (int i = nx - 1; i > 0; i--) {
488 fprintf(out, "\n");
489 for (int j = ny / 2; j > 0; j--) {
490 int i2 = (i == nx / 2 ? 0 : i);
491 int j2 = (j == ny / 2 ? 0 : j);
492 fprintf(out, "%g %g %g %g %g %g %g\n", wave->z,
493 (kx[i2] != 0 ? 2 * M_PI / kx[i2] : 0),
494 (ky[j2] != 0 ? 2 * M_PI / ky[j2] : 0),
495 kx[i2], ky[j2], A[i2][j2], phi[i2][j2]);
496 }
497 }
498
499 /* Close file... */
500 fclose(out);
501 }
502}
503
504/*****************************************************************************/
505
506void gauss(
507 wave_t *wave,
508 double fwhm) {
509
510 static double help[WX][WY];
511
512 /* Check parameters... */
513 if (fwhm <= 0)
514 return;
515
516 /* Compute sigma^2... */
517 const double sigma2 = gsl_pow_2(fwhm / 2.3548);
518
519 /* Loop over data points... */
520 for (int ix = 0; ix < wave->nx; ix++)
521 for (int iy = 0; iy < wave->ny; iy++) {
522
523 /* Init... */
524 double wsum = 0;
525 help[ix][iy] = 0;
526
527 /* Average... */
528 for (int ix2 = 0; ix2 < wave->nx; ix2++)
529 for (int iy2 = 0; iy2 < wave->ny; iy2++) {
530 double d2 = gsl_pow_2(wave->x[ix] - wave->x[ix2])
531 + gsl_pow_2(wave->y[iy] - wave->y[iy2]);
532 if (d2 <= 9 * sigma2) {
533 double w = exp(-d2 / (2 * sigma2));
534 wsum += w;
535 help[ix][iy] += w * wave->pt[ix2][iy2];
536 }
537 }
538
539 /* Normalize... */
540 wave->pt[ix][iy] = help[ix][iy] / wsum;
541 }
542}
543
544/*****************************************************************************/
545
546void hamming(
547 wave_t *wave,
548 int niter) {
549
550 static double help[WX][WY];
551
552 /* Iterations... */
553 for (int iter = 0; iter < niter; iter++) {
554
555 /* Filter in x direction... */
556 for (int ix = 0; ix < wave->nx; ix++)
557 for (int iy = 0; iy < wave->ny; iy++)
558 help[ix][iy]
559 = 0.23 * wave->pt[ix > 0 ? ix - 1 : ix][iy]
560 + 0.54 * wave->pt[ix][iy]
561 + 0.23 * wave->pt[ix < wave->nx - 1 ? ix + 1 : ix][iy];
562
563 /* Filter in y direction... */
564 for (int ix = 0; ix < wave->nx; ix++)
565 for (int iy = 0; iy < wave->ny; iy++)
566 wave->pt[ix][iy]
567 = 0.23 * help[ix][iy > 0 ? iy - 1 : iy]
568 + 0.54 * help[ix][iy]
569 + 0.23 * help[ix][iy < wave->ny - 1 ? iy + 1 : iy];
570 }
571}
572
573/*****************************************************************************/
574
575void intpol_x(
576 wave_t *wave,
577 int n) {
578
579 double dummy, x[WX], xc[WX][3], xc2[WX][3], y[WX];
580
581 /* Check parameters... */
582 if (n <= 0)
583 return;
584 if (n > WX)
585 ERRMSG("Too many data points!");
586
587 /* Set new x-coordinates... */
588 for (int i = 0; i < n; i++)
589 x[i] = LIN(0.0, wave->x[0], n - 1.0, wave->x[wave->nx - 1], i);
590
591 /* Allocate... */
592 gsl_interp_accel *acc = gsl_interp_accel_alloc();
593 gsl_spline *spline =
594 gsl_spline_alloc(gsl_interp_cspline, (size_t) wave->nx);
595
596 /* Loop over scans... */
597 for (int iy = 0; iy < wave->ny; iy++) {
598
599 /* Interpolate Cartesian coordinates... */
600 for (int ix = 0; ix < wave->nx; ix++)
601 geo2cart(0, wave->lon[ix][iy], wave->lat[ix][iy], xc[ix]);
602 for (int ic = 0; ic < 3; ic++) {
603 for (int ix = 0; ix < wave->nx; ix++)
604 y[ix] = xc[ix][ic];
605 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
606 for (int i = 0; i < n; i++)
607 xc2[i][ic] = gsl_spline_eval(spline, x[i], acc);
608 }
609 for (int i = 0; i < n; i++)
610 cart2geo(xc2[i], &dummy, &wave->lon[i][iy], &wave->lat[i][iy]);
611
612 /* Interpolate temperature... */
613 for (int ix = 0; ix < wave->nx; ix++)
614 y[ix] = wave->temp[ix][iy];
615 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
616 for (int i = 0; i < n; i++)
617 wave->temp[i][iy] = gsl_spline_eval(spline, x[i], acc);
618
619 /* Interpolate background... */
620 for (int ix = 0; ix < wave->nx; ix++)
621 y[ix] = wave->bg[ix][iy];
622 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
623 for (int i = 0; i < n; i++)
624 wave->bg[i][iy] = gsl_spline_eval(spline, x[i], acc);
625
626 /* Interpolate perturbations... */
627 for (int ix = 0; ix < wave->nx; ix++)
628 y[ix] = wave->pt[ix][iy];
629 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
630 for (int i = 0; i < n; i++)
631 wave->pt[i][iy] = gsl_spline_eval(spline, x[i], acc);
632
633 /* Interpolate variance... */
634 for (int ix = 0; ix < wave->nx; ix++)
635 y[ix] = wave->var[ix][iy];
636 gsl_spline_init(spline, wave->x, y, (size_t) wave->nx);
637 for (int i = 0; i < n; i++)
638 wave->var[i][iy] = gsl_spline_eval(spline, x[i], acc);
639 }
640
641 /* Free... */
642 gsl_spline_free(spline);
643 gsl_interp_accel_free(acc);
644
645 /* Set new x-coordinates... */
646 for (int i = 0; i < n; i++)
647 wave->x[i] = x[i];
648 wave->nx = n;
649}
650
651/*****************************************************************************/
652
653void median(
654 wave_t *wave,
655 int dx) {
656
657 static double data[WX * WY], help[WX][WY];
658
659 size_t n;
660
661 /* Check parameters... */
662 if (dx <= 0)
663 return;
664
665 /* Loop over data points... */
666 for (int ix = 0; ix < wave->nx; ix++)
667 for (int iy = 0; iy < wave->ny; iy++) {
668
669 /* Init... */
670 n = 0;
671
672 /* Get data... */
673 for (int ix2 = GSL_MAX(ix - dx, 0);
674 ix2 < GSL_MIN(ix + dx, wave->nx - 1); ix2++)
675 for (int iy2 = GSL_MAX(iy - dx, 0);
676 iy2 < GSL_MIN(iy + dx, wave->ny - 1); iy2++) {
677 data[n] = wave->pt[ix2][iy2];
678 n++;
679 }
680
681 /* Normalize... */
682 gsl_sort(data, 1, n);
683 help[ix][iy] = gsl_stats_median_from_sorted_data(data, 1, n);
684 }
685
686 /* Loop over data points... */
687 for (int ix = 0; ix < wave->nx; ix++)
688 for (int iy = 0; iy < wave->ny; iy++)
689 wave->pt[ix][iy] = help[ix][iy];
690}
691
692/*****************************************************************************/
693
694void merge_y(
695 wave_t *wave1,
696 wave_t *wave2) {
697
698 /* Check data... */
699 if (wave1->nx != wave2->nx)
700 ERRMSG("Across-track sizes do not match!");
701 if (wave1->ny + wave2->ny > WY)
702 ERRMSG("Too many data points!");
703
704 /* Get offset in y direction... */
705 const double y =
706 wave1->y[wave1->ny - 1] + (wave1->y[wave1->ny - 1] -
707 wave1->y[0]) / (wave1->ny - 1);
708
709 /* Merge data... */
710 for (int ix = 0; ix < wave2->nx; ix++)
711 for (int iy = 0; iy < wave2->ny; iy++) {
712 wave1->y[wave1->ny + iy] = y + wave2->y[iy];
713 wave1->lon[ix][wave1->ny + iy] = wave2->lon[ix][iy];
714 wave1->lat[ix][wave1->ny + iy] = wave2->lat[ix][iy];
715 wave1->temp[ix][wave1->ny + iy] = wave2->temp[ix][iy];
716 wave1->bg[ix][wave1->ny + iy] = wave2->bg[ix][iy];
717 wave1->pt[ix][wave1->ny + iy] = wave2->pt[ix][iy];
718 wave1->var[ix][wave1->ny + iy] = wave2->var[ix][iy];
719 }
720
721 /* Increment counter... */
722 wave1->ny += wave2->ny;
723}
724
725/*****************************************************************************/
726
727void noise(
728 wave_t *wave,
729 double *mu,
730 double *sig) {
731
732 int n = 0;
733
734 /* Init... */
735 *mu = 0;
736 *sig = 0;
737
738 /* Estimate noise (Immerkaer, 1996)... */
739 for (int ix = 1; ix < wave->nx - 1; ix++)
740 for (int iy = 1; iy < wave->ny - 1; iy++) {
741
742 /* Check data... */
743 int okay = 1;
744 for (int ix2 = ix - 1; ix2 <= ix + 1; ix2++)
745 for (int iy2 = iy - 1; iy2 <= iy + 1; iy2++)
746 if (!gsl_finite(wave->temp[ix2][iy2]))
747 okay = 0;
748 if (!okay)
749 continue;
750
751 /* Get mean noise... */
752 n++;
753 *mu += wave->temp[ix][iy];
754 *sig += gsl_pow_2(+4. / 6. * wave->temp[ix][iy]
755 - 2. / 6. * (wave->temp[ix - 1][iy]
756 + wave->temp[ix + 1][iy]
757 + wave->temp[ix][iy - 1]
758 + wave->temp[ix][iy + 1])
759 + 1. / 6. * (wave->temp[ix - 1][iy - 1]
760 + wave->temp[ix + 1][iy - 1]
761 + wave->temp[ix - 1][iy + 1]
762 + wave->temp[ix + 1][iy + 1]));
763 }
764
765 /* Normalize... */
766 *mu /= (double) n;
767 *sig = sqrt(*sig / (double) n);
768}
769
770/*****************************************************************************/
771
772void noise_pert(
773 pert_t *pert,
774 int track0,
775 int track1,
776 double *mu,
777 double *sig) {
778
779 int n = 0;
780
781 /* Init... */
782 *mu = 0;
783 *sig = 0;
784
785 /* Estimate noise (Immerkaer, 1996)... */
786 for (int track = track0; track < track1; track++)
787 for (int xtrack = 0; xtrack < pert->nxtrack; xtrack++) {
788
789 /* Check data... */
790 int okay = 1;
791 for (int ifov = 0; ifov < pert->nfov; ifov++)
792 if (!(gsl_finite(pert->bt[track][xtrack][ifov])
793 && gsl_finite(pert->pt[track][xtrack][ifov])))
794 okay = 0;
795 if (!okay)
796 continue;
797
798 /* Get mean noise... */
799 n++;
800 *mu += pert->bt[track][xtrack][4];
801 *sig += gsl_pow_2(+4. / 6. * pert->pt[track][xtrack][4]
802 - 2. / 6. * (pert->pt[track][xtrack][1]
803 + pert->pt[track][xtrack][3]
804 + pert->pt[track][xtrack][5]
805 + pert->pt[track][xtrack][7])
806 + 1. / 6. * (pert->pt[track][xtrack][0]
807 + pert->pt[track][xtrack][2]
808 + pert->pt[track][xtrack][6]
809 + pert->pt[track][xtrack][8]));
810 }
811
812 /* Normalize... */
813 *mu /= (double) n;
814 *sig = sqrt(*sig / (double) n);
815}
816
817/*****************************************************************************/
818
819void period(
820 wave_t *wave,
821 double lxymax,
822 double dlxy,
823 double *Amax,
824 double *phimax,
825 double *lhmax,
826 double *kxmax,
827 double *kymax,
828 double *alphamax,
829 double *betamax,
830 char *filename) {
831
832 static double kx[PMAX], ky[PMAX], A[PMAX][PMAX], phi[PMAX][PMAX],
833 cx[PMAX][WX], cy[PMAX][WY], sx[PMAX][WX], sy[PMAX][WY], a, b, c;
834
835 int imin, imax, jmin, jmax, lmax = 0, mmax = 0;
836
837 /* Compute wavenumbers and periodogram coefficients... */
838 for (double lx = -lxymax; lx <= lxymax; lx += dlxy) {
839 kx[lmax] = (lx != 0 ? 2 * M_PI / lx : 0);
840 for (int i = 0; i < wave->nx; i++) {
841 cx[lmax][i] = cos(kx[lmax] * wave->x[i]);
842 sx[lmax][i] = sin(kx[lmax] * wave->x[i]);
843 }
844 if ((++lmax) > PMAX)
845 ERRMSG("Too many wavenumbers for periodogram!");
846 }
847 for (double ly = 0; ly <= lxymax; ly += dlxy) {
848 ky[mmax] = (ly != 0 ? 2 * M_PI / ly : 0);
849 for (int j = 0; j < wave->ny; j++) {
850 cy[mmax][j] = cos(ky[mmax] * wave->y[j]);
851 sy[mmax][j] = sin(ky[mmax] * wave->y[j]);
852 }
853 if ((++mmax) > PMAX)
854 ERRMSG("Too many wavenumbers for periodogram!");
855 }
856
857 /* Find area... */
858 imin = jmin = 9999;
859 imax = jmax = -9999;
860 for (int i = 0; i < wave->nx; i++)
861 for (int j = 0; j < wave->ny; j++)
862 if (gsl_finite(wave->var[i][j])) {
863 imin = GSL_MIN(imin, i);
864 imax = GSL_MAX(imax, i);
865 jmin = GSL_MIN(jmin, j);
866 jmax = GSL_MAX(jmax, j);
867 }
868
869 /* Get Nyquist frequencies... */
870 const double kx_ny =
871 M_PI / fabs((wave->x[imax] - wave->x[imin]) /
872 ((double) imax - (double) imin));
873 const double ky_ny =
874 M_PI / fabs((wave->y[jmax] - wave->y[jmin]) /
875 ((double) jmax - (double) jmin));
876
877 /* Loop over wavelengths... */
878 for (int l = 0; l < lmax; l++)
879 for (int m = 0; m < mmax; m++) {
880
881 /* Check frequencies... */
882 if (kx[l] == 0 || fabs(kx[l]) > kx_ny ||
883 ky[m] == 0 || fabs(ky[m]) > ky_ny) {
884 A[l][m] = GSL_NAN;
885 phi[l][m] = GSL_NAN;
886 continue;
887 }
888
889 /* Compute periodogram... */
890 a = b = c = 0;
891 for (int i = imin; i <= imax; i++)
892 for (int j = jmin; j <= jmax; j++)
893 if (gsl_finite(wave->var[i][j])) {
894 a += wave->pt[i][j] * (cx[l][i] * cy[m][j] - sx[l][i] * sy[m][j]);
895 b += wave->pt[i][j] * (sx[l][i] * cy[m][j] + cx[l][i] * sy[m][j]);
896 c++;
897 }
898 a *= 2. / c;
899 b *= 2. / c;
900
901 /* Get amplitude and phase... */
902 A[l][m] = sqrt(gsl_pow_2(a) + gsl_pow_2(b));
903 phi[l][m] = atan2(b, a) * 180. / M_PI;
904 }
905
906 /* Find maximum... */
907 *Amax = 0;
908 for (int l = 0; l < lmax; l++)
909 for (int m = 0; m < mmax; m++)
910 if (gsl_finite(A[l][m]) && A[l][m] > *Amax) {
911 *Amax = A[l][m];
912 *phimax = phi[l][m];
913 *kxmax = kx[l];
914 *kymax = ky[m];
915 }
916
917 /* Get horizontal wavelength... */
918 *lhmax = 2 * M_PI / sqrt(gsl_pow_2(*kxmax) + gsl_pow_2(*kymax));
919
920 /* Get propagation direction in xy-plane... */
921 *alphamax = 90. - 180. / M_PI * atan2(*kxmax, *kymax);
922
923 /* Get propagation direction in lon,lat-plane... */
924 *betamax = *alphamax
925 +
926 180. / M_PI *
927 atan2(wave->lat[wave->nx / 2 >
928 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
929 - wave->lat[wave->nx / 2 <
930 wave->nx - 1 ? wave->nx / 2 +
931 1 : wave->nx / 2][wave->ny / 2],
932 wave->lon[wave->nx / 2 >
933 0 ? wave->nx / 2 - 1 : wave->nx / 2][wave->ny / 2]
934 - wave->lon[wave->nx / 2 <
935 wave->nx - 1 ? wave->nx / 2 +
936 1 : wave->nx / 2][wave->ny / 2]);
937
938 /* Save periodogram data... */
939 if (filename != NULL) {
940
941 /* Write info... */
942 LOG(1, "Write periodogram data: %s", filename);
943
944 /* Create file... */
945 FILE *out;
946 if (!(out = fopen(filename, "w")))
947 ERRMSG("Cannot create file!");
948
949 /* Write header... */
950 fprintf(out,
951 "# $1 = altitude [km]\n"
952 "# $2 = wavelength in x-direction [km]\n"
953 "# $3 = wavelength in y-direction [km]\n"
954 "# $4 = wavenumber in x-direction [1/km]\n"
955 "# $5 = wavenumber in y-direction [1/km]\n"
956 "# $6 = amplitude [K]\n" "# $7 = phase [rad]\n");
957
958 /* Write data... */
959 for (int l = 0; l < lmax; l++) {
960 fprintf(out, "\n");
961 for (int m = 0; m < mmax; m++)
962 fprintf(out, "%g %g %g %g %g %g %g\n", wave->z,
963 (kx[l] != 0 ? 2 * M_PI / kx[l] : 0),
964 (ky[m] != 0 ? 2 * M_PI / ky[m] : 0),
965 kx[l], ky[m], A[l][m], phi[l][m]);
966 }
967
968 /* Close file... */
969 fclose(out);
970 }
971}
972
973/*****************************************************************************/
974
975#if 0
976void pert2wave(
977 pert_t *pert,
978 wave_t *wave,
979 int track0,
980 int track1,
981 int xtrack0,
982 int xtrack1) {
983
984 double x0[3], x1[3];
985
986 /* Check ranges... */
987 track0 = GSL_MIN(GSL_MAX(track0, 0), pert->ntrack - 1);
988 track1 = GSL_MIN(GSL_MAX(track1, 0), pert->ntrack - 1);
989 xtrack0 = GSL_MIN(GSL_MAX(xtrack0, 0), pert->nxtrack - 1);
990 xtrack1 = GSL_MIN(GSL_MAX(xtrack1, 0), pert->nxtrack - 1);
991
992 /* Set size... */
993 wave->nx = xtrack1 - xtrack0 + 1;
994 if (wave->nx > WX)
995 ERRMSG("Too many across-track values!");
996 wave->ny = track1 - track0 + 1;
997 if (wave->ny > WY)
998 ERRMSG("Too many along-track values!");
999
1000 /* Loop over footprints... */
1001 for (int itrack = track0; itrack <= track1; itrack++)
1002 for (int ixtrack = xtrack0; ixtrack <= xtrack1; ixtrack++) {
1003
1004 /* Get distances... */
1005 if (itrack == track0) {
1006 wave->x[0] = 0;
1007 if (ixtrack > xtrack0) {
1008 geo2cart(0, pert->lon[itrack][ixtrack - 1],
1009 pert->lat[itrack][ixtrack - 1], x0);
1010 geo2cart(0, pert->lon[itrack][ixtrack],
1011 pert->lat[itrack][ixtrack], x1);
1012 wave->x[ixtrack - xtrack0] =
1013 wave->x[ixtrack - xtrack0 - 1] + DIST(x0, x1);
1014 }
1015 }
1016 if (ixtrack == xtrack0) {
1017 wave->y[0] = 0;
1018 if (itrack > track0) {
1019 geo2cart(0, pert->lon[itrack - 1][ixtrack],
1020 pert->lat[itrack - 1][ixtrack], x0);
1021 geo2cart(0, pert->lon[itrack][ixtrack],
1022 pert->lat[itrack][ixtrack], x1);
1023 wave->y[itrack - track0] =
1024 wave->y[itrack - track0 - 1] + DIST(x0, x1);
1025 }
1026 }
1027
1028 /* Save geolocation... */
1029 wave->time = pert->time[(track0 + track1) / 2][(xtrack0 + xtrack1) / 2];
1030 wave->z = 0;
1031 wave->lon[ixtrack - xtrack0][itrack - track0] =
1032 pert->lon[itrack][ixtrack];
1033 wave->lat[ixtrack - xtrack0][itrack - track0] =
1034 pert->lat[itrack][ixtrack];
1035
1036 /* Save temperature data... */
1037 wave->temp[ixtrack - xtrack0][itrack - track0]
1038 = pert->bt[itrack][ixtrack];
1039 wave->bg[ixtrack - xtrack0][itrack - track0]
1040 = pert->bt[itrack][ixtrack] - pert->pt[itrack][ixtrack];
1041 wave->pt[ixtrack - xtrack0][itrack - track0]
1042 = pert->pt[itrack][ixtrack];
1043 wave->var[ixtrack - xtrack0][itrack - track0]
1044 = pert->var[itrack][ixtrack];
1045 }
1046}
1047#endif
1048
1049/*****************************************************************************/
1050
1051int read_cris_l1(
1052 char *filename,
1053 cris_l1_t *l1,
1054 int apo) {
1055
1056 int ncid, varid, dimid;
1057
1058 size_t n;
1059
1060 /* Write info... */
1061 LOG(1, "Read CrIS Level-1B file: %s", filename);
1062
1063 /* Open netCDF file... */
1064 if (nc_open(filename, NC_NOWRITE, &ncid) != NC_NOERR) {
1065 WARN("Cannot open file!");
1066 return 0;
1067 }
1068
1069 /* Check dimensions... */
1070 NC(nc_inq_dimid(ncid, "atrack", &dimid));
1071 NC(nc_inq_dimlen(ncid, dimid, &n));
1072 if (n != L1_NTRACK)
1073 ERRMSG("Level-1B file dimension mismatch!");
1074
1075 NC(nc_inq_dimid(ncid, "xtrack", &dimid));
1076 NC(nc_inq_dimlen(ncid, dimid, &n));
1077 if (n != L1_NXTRACK)
1078 ERRMSG("Level-1B file dimension mismatch!");
1079
1080 NC(nc_inq_dimid(ncid, "fov", &dimid));
1081 NC(nc_inq_dimlen(ncid, dimid, &n));
1082 if (n != L1_NFOV)
1083 ERRMSG("Level-1B file dimension mismatch!");
1084
1085 NC(nc_inq_dimid(ncid, "wnum_lw", &dimid));
1086 NC(nc_inq_dimlen(ncid, dimid, &n));
1087 if (n != L1_NCHAN_LW)
1088 ERRMSG("Level-1B file dimension mismatch!");
1089
1090 NC(nc_inq_dimid(ncid, "wnum_mw", &dimid));
1091 NC(nc_inq_dimlen(ncid, dimid, &n));
1092 if (n != L1_NCHAN_MW)
1093 ERRMSG("Level-1B file dimension mismatch!");
1094
1095 NC(nc_inq_dimid(ncid, "wnum_sw", &dimid));
1096 NC(nc_inq_dimlen(ncid, dimid, &n));
1097 if (n != L1_NCHAN_SW)
1098 ERRMSG("Level-1B file dimension mismatch!");
1099
1100 /* Read satellite position... */
1101 NC(nc_inq_varid(ncid, "obs_time_tai93", &varid));
1102 NC(nc_get_var_double(ncid, varid, l1->time[0]));
1103 NC(nc_inq_varid(ncid, "lon", &varid));
1104 NC(nc_get_var_double(ncid, varid, l1->lon[0][0]));
1105 NC(nc_inq_varid(ncid, "lat", &varid));
1106 NC(nc_get_var_double(ncid, varid, l1->lat[0][0]));
1107 NC(nc_inq_varid(ncid, "sat_alt", &varid));
1108 NC(nc_get_var_double(ncid, varid, l1->sat_z));
1109 NC(nc_inq_varid(ncid, "subsat_lon", &varid));
1110 NC(nc_get_var_double(ncid, varid, l1->sat_lon));
1111 NC(nc_inq_varid(ncid, "subsat_lat", &varid));
1112 NC(nc_get_var_double(ncid, varid, l1->sat_lat));
1113
1114 /* Convert units... */
1115 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1116 l1->sat_z[itrack] /= 1e3;
1117
1118 /* Read spectra... */
1119 NC(nc_inq_varid(ncid, "wnum_lw", &varid));
1120 NC(nc_get_var_double(ncid, varid, l1->nu_lw));
1121 NC(nc_inq_varid(ncid, "rad_lw", &varid));
1122 NC(nc_get_var_float(ncid, varid, l1->rad_lw[0][0][0]));
1123 NC(nc_inq_varid(ncid, "wnum_mw", &varid));
1124 NC(nc_get_var_double(ncid, varid, l1->nu_mw));
1125 NC(nc_inq_varid(ncid, "rad_mw", &varid));
1126 NC(nc_get_var_float(ncid, varid, l1->rad_mw[0][0][0]));
1127 NC(nc_inq_varid(ncid, "wnum_sw", &varid));
1128 NC(nc_get_var_double(ncid, varid, l1->nu_sw));
1129 NC(nc_inq_varid(ncid, "rad_sw", &varid));
1130 NC(nc_get_var_float(ncid, varid, l1->rad_sw[0][0][0]));
1131
1132 /* Read noise... */
1133 NC(nc_inq_varid(ncid, "nedn_lw", &varid));
1134 NC(nc_get_var_float(ncid, varid, l1->nedn_lw[0]));
1135 NC(nc_inq_varid(ncid, "nedn_mw", &varid));
1136 NC(nc_get_var_float(ncid, varid, l1->nedn_mw[0]));
1137 NC(nc_inq_varid(ncid, "nedn_sw", &varid));
1138 NC(nc_get_var_float(ncid, varid, l1->nedn_sw[0]));
1139
1140 /* Check quality flags... */
1141 NC(nc_inq_varid(ncid, "rad_lw_qc", &varid));
1142 NC(nc_get_var_short(ncid, varid, l1->qual_lw[0][0]));
1143 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1144 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1145 for (int ifov = 0; ifov < L1_NFOV; ifov++)
1146 if (l1->qual_lw[itrack][ixtrack][ifov] > 1)
1147 for (int ichan = 0; ichan < L1_NCHAN_LW; ichan++)
1148 l1->rad_lw[itrack][ixtrack][ifov][ichan] = GSL_NAN;
1149
1150 NC(nc_inq_varid(ncid, "rad_mw_qc", &varid));
1151 NC(nc_get_var_short(ncid, varid, l1->qual_mw[0][0]));
1152 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1153 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1154 for (int ifov = 0; ifov < L1_NFOV; ifov++)
1155 if (l1->qual_mw[itrack][ixtrack][ifov] > 1)
1156 for (int ichan = 0; ichan < L1_NCHAN_MW; ichan++)
1157 l1->rad_mw[itrack][ixtrack][ifov][ichan] = GSL_NAN;
1158
1159 NC(nc_inq_varid(ncid, "rad_sw_qc", &varid));
1160 NC(nc_get_var_short(ncid, varid, l1->qual_sw[0][0]));
1161 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1162 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1163 for (int ifov = 0; ifov < L1_NFOV; ifov++)
1164 if (l1->qual_sw[itrack][ixtrack][ifov] > 1)
1165 for (int ichan = 0; ichan < L1_NCHAN_SW; ichan++)
1166 l1->rad_sw[itrack][ixtrack][ifov][ichan] = GSL_NAN;
1167
1168 /* Apodization... */
1169 if (apo)
1170 for (int itrack = 0; itrack < L1_NTRACK; itrack++)
1171 for (int ixtrack = 0; ixtrack < L1_NXTRACK; ixtrack++)
1172 for (int ifov = 0; ifov < L1_NFOV; ifov++) {
1173
1174 double help[L1_NCHAN_MW];
1175 for (int ichan = 0; ichan < L1_NCHAN_LW; ichan++)
1176 help[ichan] = 0;
1177 for (int ichan = 1; ichan < L1_NCHAN_LW - 1; ichan++)
1178 help[ichan]
1179 = 0.23 * l1->rad_lw[itrack][ixtrack][ifov][ichan - 1]
1180 + 0.54 * l1->rad_lw[itrack][ixtrack][ifov][ichan]
1181 + 0.23 * l1->rad_lw[itrack][ixtrack][ifov][ichan + 1];
1182 for (int ichan = 1; ichan < L1_NCHAN_LW - 1; ichan++)
1183 l1->rad_lw[itrack][ixtrack][ifov][ichan] = (float) help[ichan];
1184 l1->rad_lw[itrack][ixtrack][ifov][0] =
1185 l1->rad_lw[itrack][ixtrack][ifov][L1_NCHAN_LW - 1] = GSL_NAN;
1186
1187 for (int ichan = 0; ichan < L1_NCHAN_MW; ichan++)
1188 help[ichan] = 0;
1189 for (int ichan = 1; ichan < L1_NCHAN_MW - 1; ichan++)
1190 help[ichan]
1191 = 0.23 * l1->rad_mw[itrack][ixtrack][ifov][ichan - 1]
1192 + 0.54 * l1->rad_mw[itrack][ixtrack][ifov][ichan]
1193 + 0.23 * l1->rad_mw[itrack][ixtrack][ifov][ichan + 1];
1194 for (int ichan = 1; ichan < L1_NCHAN_MW - 1; ichan++)
1195 l1->rad_mw[itrack][ixtrack][ifov][ichan] = (float) help[ichan];
1196 l1->rad_mw[itrack][ixtrack][ifov][0] =
1197 l1->rad_mw[itrack][ixtrack][ifov][L1_NCHAN_MW - 1] = GSL_NAN;
1198
1199 for (int ichan = 0; ichan < L1_NCHAN_SW; ichan++)
1200 help[ichan] = 0;
1201 for (int ichan = 1; ichan < L1_NCHAN_SW - 1; ichan++)
1202 help[ichan]
1203 = 0.23 * l1->rad_sw[itrack][ixtrack][ifov][ichan - 1]
1204 + 0.54 * l1->rad_sw[itrack][ixtrack][ifov][ichan]
1205 + 0.23 * l1->rad_sw[itrack][ixtrack][ifov][ichan + 1];
1206 for (int ichan = 1; ichan < L1_NCHAN_SW - 1; ichan++)
1207 l1->rad_sw[itrack][ixtrack][ifov][ichan] = (float) help[ichan];
1208 l1->rad_sw[itrack][ixtrack][ifov][0] =
1209 l1->rad_sw[itrack][ixtrack][ifov][L1_NCHAN_SW - 1] = GSL_NAN;
1210 }
1211
1212 /* Close file... */
1213 NC(nc_close(ncid));
1214
1215 /* Return success... */
1216 return 1;
1217}
1218
1219/*****************************************************************************/
1220
1221void read_pert(
1222 char *filename,
1223 char *pertname,
1224 int dc,
1225 pert_t *pert) {
1226
1227 static char varname[LEN];
1228
1229 static double help[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV];
1230
1231 static int dimid[3], ncid, varid;
1232
1233 static size_t ntrack, nxtrack, nfov, start[3] = { 0, 0, 0 }, count[3] =
1234 { 1, 1, 1 };
1235
1236 /* Write info... */
1237 LOG(1, "Read perturbation data: %s", filename);
1238
1239 /* Open netCDF file... */
1240 NC(nc_open(filename, NC_NOWRITE, &ncid));
1241
1242 /* Get dimensions... */
1243 NC(nc_inq_dimid(ncid, "NTRACK", &dimid[0]));
1244 NC(nc_inq_dimlen(ncid, dimid[0], &ntrack));
1245 if (ntrack > PERT_NTRACK)
1246 ERRMSG("Too many scans!");
1247
1248 NC(nc_inq_dimid(ncid, "NXTRACK", &dimid[1]));
1249 NC(nc_inq_dimlen(ncid, dimid[1], &nxtrack));
1250 if (nxtrack > PERT_NXTRACK)
1251 ERRMSG("Too many tracks!");
1252
1253 NC(nc_inq_dimid(ncid, "NFOV", &dimid[2]));
1254 NC(nc_inq_dimlen(ncid, dimid[2], &nfov));
1255 if (nfov > PERT_NFOV)
1256 ERRMSG("Too many field of views!");
1257
1258 pert->ntrack = (int) ntrack;
1259 pert->nxtrack = (int) nxtrack;
1260 pert->nfov = (int) nfov;
1261 count[2] = nfov;
1262
1263 /* Read data... */
1264 NC(nc_inq_varid(ncid, "time", &varid));
1265 for (size_t itrack = 0; itrack < ntrack; itrack++)
1266 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1267 start[0] = itrack;
1268 start[1] = ixtrack;
1269 NC(nc_get_vara_double
1270 (ncid, varid, start, count, pert->time[itrack][ixtrack]));
1271 }
1272
1273 NC(nc_inq_varid(ncid, "lon", &varid));
1274 for (size_t itrack = 0; itrack < ntrack; itrack++)
1275 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1276 start[0] = itrack;
1277 start[1] = ixtrack;
1278 NC(nc_get_vara_double
1279 (ncid, varid, start, count, pert->lon[itrack][ixtrack]));
1280 }
1281
1282 NC(nc_inq_varid(ncid, "lat", &varid));
1283 for (size_t itrack = 0; itrack < ntrack; itrack++)
1284 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1285 start[0] = itrack;
1286 start[1] = ixtrack;
1287 NC(nc_get_vara_double
1288 (ncid, varid, start, count, pert->lat[itrack][ixtrack]));
1289 }
1290
1291 NC(nc_inq_varid(ncid, "bt_8mu", &varid));
1292 for (size_t itrack = 0; itrack < ntrack; itrack++)
1293 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1294 start[0] = itrack;
1295 start[1] = ixtrack;
1296 NC(nc_get_vara_double
1297 (ncid, varid, start, count, pert->dc[itrack][ixtrack]));
1298 }
1299
1300 NC(nc_inq_varid(ncid, "bt_10mu", &varid));
1301 for (size_t itrack = 0; itrack < ntrack; itrack++)
1302 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1303 start[0] = itrack;
1304 start[1] = ixtrack;
1305 NC(nc_get_vara_double
1306 (ncid, varid, start, count, help[itrack][ixtrack]));
1307 }
1308
1309 for (size_t itrack = 0; itrack < ntrack; itrack++)
1310 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++)
1311 for (size_t ifov = 0; ifov < nfov; ifov++)
1312 if (dc == 2
1313 || (dc == 0 && !gsl_finite(pert->dc[itrack][ixtrack][ifov])))
1314 pert->dc[itrack][ixtrack][ifov] = help[itrack][ixtrack][ifov];
1315
1316 sprintf(varname, "bt_%s", pertname);
1317 NC(nc_inq_varid(ncid, varname, &varid));
1318 for (size_t itrack = 0; itrack < ntrack; itrack++)
1319 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1320 start[0] = itrack;
1321 start[1] = ixtrack;
1322 NC(nc_get_vara_double
1323 (ncid, varid, start, count, pert->bt[itrack][ixtrack]));
1324 }
1325
1326 sprintf(varname, "bt_%s_pt", pertname);
1327 NC(nc_inq_varid(ncid, varname, &varid));
1328 for (size_t itrack = 0; itrack < ntrack; itrack++)
1329 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1330 start[0] = itrack;
1331 start[1] = ixtrack;
1332 NC(nc_get_vara_double
1333 (ncid, varid, start, count, pert->pt[itrack][ixtrack]));
1334 }
1335
1336 sprintf(varname, "bt_%s_var", pertname);
1337 NC(nc_inq_varid(ncid, varname, &varid));
1338 for (size_t itrack = 0; itrack < ntrack; itrack++)
1339 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1340 start[0] = itrack;
1341 start[1] = ixtrack;
1342 NC(nc_get_vara_double
1343 (ncid, varid, start, count, pert->var[itrack][ixtrack]));
1344 }
1345
1346 /* Close file... */
1347 NC(nc_close(ncid));
1348}
1349
1350/*****************************************************************************/
1351
1352void read_retr(
1353 char *filename,
1354 retr_t *ret) {
1355
1356 static double help[NDS * NPG];
1357
1358 int dimid, ids, ncid, varid;
1359
1360 size_t nds, np, ntrack, nxtrack;
1361
1362 /* Write info... */
1363 LOG(1, "Read retrieval data: %s", filename);
1364
1365 /* Open netCDF file... */
1366 NC(nc_open(filename, NC_NOWRITE, &ncid));
1367
1368 /* Read new retrieval file format... */
1369 if (nc_inq_dimid(ncid, "L1_NTRACK", &dimid) == NC_NOERR) {
1370
1371 /* Get dimensions... */
1372 NC(nc_inq_dimid(ncid, "RET_NP", &dimid));
1373 NC(nc_inq_dimlen(ncid, dimid, &np));
1374 ret->np = (int) np;
1375 if (ret->np > NPG)
1376 ERRMSG("Too many data points!");
1377
1378 NC(nc_inq_dimid(ncid, "L1_NTRACK", &dimid));
1379 NC(nc_inq_dimlen(ncid, dimid, &ntrack));
1380 NC(nc_inq_dimid(ncid, "L1_NXTRACK", &dimid));
1381 NC(nc_inq_dimlen(ncid, dimid, &nxtrack));
1382 ret->nds = (int) (ntrack * nxtrack);
1383 if (ret->nds > NDS)
1384 ERRMSG("Too many data sets!");
1385
1386 /* Read time... */
1387 NC(nc_inq_varid(ncid, "l1_time", &varid));
1388 NC(nc_get_var_double(ncid, varid, help));
1389 ids = 0;
1390 for (size_t itrack = 0; itrack < ntrack; itrack++)
1391 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1392 for (int ip = 0; ip < ret->np; ip++)
1393 ret->time[ids][ip] = help[ids];
1394 ids++;
1395 }
1396
1397 /* Read altitudes... */
1398 NC(nc_inq_varid(ncid, "ret_z", &varid));
1399 NC(nc_get_var_double(ncid, varid, help));
1400 ids = 0;
1401 for (size_t itrack = 0; itrack < ntrack; itrack++)
1402 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1403 for (int ip = 0; ip < ret->np; ip++)
1404 ret->z[ids][ip] = help[ip];
1405 ids++;
1406 }
1407
1408 /* Read longitudes... */
1409 NC(nc_inq_varid(ncid, "l1_lon", &varid));
1410 NC(nc_get_var_double(ncid, varid, help));
1411 ids = 0;
1412 for (size_t itrack = 0; itrack < ntrack; itrack++)
1413 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1414 for (int ip = 0; ip < ret->np; ip++)
1415 ret->lon[ids][ip] = help[ids];
1416 ids++;
1417 }
1418
1419 /* Read latitudes... */
1420 NC(nc_inq_varid(ncid, "l1_lat", &varid));
1421 NC(nc_get_var_double(ncid, varid, help));
1422 ids = 0;
1423 for (size_t itrack = 0; itrack < ntrack; itrack++)
1424 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1425 for (int ip = 0; ip < ret->np; ip++)
1426 ret->lat[ids][ip] = help[ids];
1427 ids++;
1428 }
1429
1430 /* Read temperatures... */
1431 NC(nc_inq_varid(ncid, "ret_temp", &varid));
1432 NC(nc_get_var_double(ncid, varid, help));
1433 ids = 0;
1434 for (size_t itrack = 0; itrack < ntrack; itrack++)
1435 for (size_t ixtrack = 0; ixtrack < nxtrack; ixtrack++) {
1436 for (int ip = 0; ip < ret->np; ip++)
1437 ret->t[ids][ip] =
1438 help[(itrack * nxtrack + ixtrack) * (size_t) np + (size_t) ip];
1439 ids++;
1440 }
1441 }
1442
1443 /* Read old retrieval file format... */
1444 if (nc_inq_dimid(ncid, "np", &dimid) == NC_NOERR) {
1445
1446 /* Get dimensions... */
1447 NC(nc_inq_dimid(ncid, "np", &dimid));
1448 NC(nc_inq_dimlen(ncid, dimid, &np));
1449 ret->np = (int) np;
1450 if (ret->np > NPG)
1451 ERRMSG("Too many data points!");
1452
1453 NC(nc_inq_dimid(ncid, "nds", &dimid));
1454 NC(nc_inq_dimlen(ncid, dimid, &nds));
1455 ret->nds = (int) nds;
1456 if (ret->nds > NDS)
1457 ERRMSG("Too many data sets!");
1458
1459 /* Read data... */
1460 NC(nc_inq_varid(ncid, "time", &varid));
1461 NC(nc_get_var_double(ncid, varid, help));
1462 read_retr_help(help, ret->nds, ret->np, ret->time);
1463
1464 NC(nc_inq_varid(ncid, "z", &varid));
1465 NC(nc_get_var_double(ncid, varid, help));
1466 read_retr_help(help, ret->nds, ret->np, ret->z);
1467
1468 NC(nc_inq_varid(ncid, "lon", &varid));
1469 NC(nc_get_var_double(ncid, varid, help));
1470 read_retr_help(help, ret->nds, ret->np, ret->lon);
1471
1472 NC(nc_inq_varid(ncid, "lat", &varid));
1473 NC(nc_get_var_double(ncid, varid, help));
1474 read_retr_help(help, ret->nds, ret->np, ret->lat);
1475
1476 NC(nc_inq_varid(ncid, "press", &varid));
1477 NC(nc_get_var_double(ncid, varid, help));
1478 read_retr_help(help, ret->nds, ret->np, ret->p);
1479
1480 NC(nc_inq_varid(ncid, "temp", &varid));
1481 NC(nc_get_var_double(ncid, varid, help));
1482 read_retr_help(help, ret->nds, ret->np, ret->t);
1483
1484 NC(nc_inq_varid(ncid, "temp_apr", &varid));
1485 NC(nc_get_var_double(ncid, varid, help));
1486 read_retr_help(help, ret->nds, ret->np, ret->t_apr);
1487
1488 NC(nc_inq_varid(ncid, "temp_total", &varid));
1489 NC(nc_get_var_double(ncid, varid, help));
1490 read_retr_help(help, ret->nds, ret->np, ret->t_tot);
1491
1492 NC(nc_inq_varid(ncid, "temp_noise", &varid));
1493 NC(nc_get_var_double(ncid, varid, help));
1494 read_retr_help(help, ret->nds, ret->np, ret->t_noise);
1495
1496 NC(nc_inq_varid(ncid, "temp_formod", &varid));
1497 NC(nc_get_var_double(ncid, varid, help));
1498 read_retr_help(help, ret->nds, ret->np, ret->t_fm);
1499
1500 NC(nc_inq_varid(ncid, "temp_cont", &varid));
1501 NC(nc_get_var_double(ncid, varid, help));
1502 read_retr_help(help, ret->nds, ret->np, ret->t_cont);
1503
1504 NC(nc_inq_varid(ncid, "temp_res", &varid));
1505 NC(nc_get_var_double(ncid, varid, help));
1506 read_retr_help(help, ret->nds, ret->np, ret->t_res);
1507
1508 NC(nc_inq_varid(ncid, "chisq", &varid));
1509 NC(nc_get_var_double(ncid, varid, ret->chisq));
1510 }
1511
1512 /* Close file... */
1513 NC(nc_close(ncid));
1514}
1515
1516/*****************************************************************************/
1517
1518void read_retr_help(
1519 double *help,
1520 int nds,
1521 int np,
1522 double mat[NDS][NPG]) {
1523
1524 int n = 0;
1525
1526 for (int ip = 0; ip < np; ip++)
1527 for (int ids = 0; ids < nds; ids++)
1528 mat[ids][ip] = help[n++];
1529}
1530
1531/*****************************************************************************/
1532
1533void read_wave(
1534 char *filename,
1535 wave_t *wave) {
1536
1537 FILE *in;
1538
1539 char line[LEN];
1540
1541 double rtime, rz, rlon, rlat, rx, ry, ryold = -1e10,
1542 rtemp, rbg, rpt, rvar, rfit;
1543
1544 /* Init... */
1545 wave->nx = 0;
1546 wave->ny = 0;
1547
1548 /* Write info... */
1549 LOG(1, "Read wave data: %s", filename);
1550
1551 /* Open file... */
1552 if (!(in = fopen(filename, "r")))
1553 ERRMSG("Cannot open file!");
1554
1555 /* Read data... */
1556 while (fgets(line, LEN, in))
1557 if (sscanf(line, "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", &rtime,
1558 &rz, &rlon, &rlat, &rx, &ry, &rtemp, &rbg, &rpt,
1559 &rvar, &rfit) == 11) {
1560
1561 /* Set index... */
1562 if (ry != ryold) {
1563 if ((++wave->ny >= WY))
1564 ERRMSG("Too many y-values!");
1565 wave->nx = 0;
1566 } else if ((++wave->nx) >= WX)
1567 ERRMSG("Too many x-values!");
1568 ryold = ry;
1569
1570 /* Save data... */
1571 wave->time = rtime;
1572 wave->z = rz;
1573 wave->lon[wave->nx][wave->ny] = rlon;
1574 wave->lat[wave->nx][wave->ny] = rlat;
1575 wave->x[wave->nx] = rx;
1576 wave->y[wave->ny] = ry;
1577 wave->temp[wave->nx][wave->ny] = rtemp;
1578 wave->bg[wave->nx][wave->ny] = rbg;
1579 wave->pt[wave->nx][wave->ny] = rpt;
1580 wave->var[wave->nx][wave->ny] = rvar;
1581 wave->var[wave->nx][wave->ny] = rfit;
1582 }
1583
1584 /* Increment counters... */
1585 wave->nx++;
1586 wave->ny++;
1587
1588 /* Close file... */
1589 fclose(in);
1590}
1591
1592/*****************************************************************************/
1593
1594void ret2wave(
1595 retr_t *ret,
1596 wave_t *wave,
1597 int dataset,
1598 int ip) {
1599
1600 double x0[3], x1[3];
1601
1602 /* Initialize... */
1603 wave->nx = 90;
1604 if (wave->nx > WX)
1605 ERRMSG("Too many across-track values!");
1606 wave->ny = 135;
1607 if (wave->ny > WY)
1608 ERRMSG("Too many along-track values!");
1609 if (ip < 0 || ip >= ret->np)
1610 ERRMSG("Altitude index out of range!");
1611
1612 /* Loop over data sets and data points... */
1613 for (int ids = 0; ids < ret->nds; ids++) {
1614
1615 /* Get horizontal indices... */
1616 int ix = ids % 90;
1617 int iy = ids / 90;
1618
1619 /* Get distances... */
1620 if (iy == 0) {
1621 geo2cart(0.0, ret->lon[0][0], ret->lat[0][0], x0);
1622 geo2cart(0.0, ret->lon[ids][ip], ret->lat[ids][ip], x1);
1623 wave->x[ix] = DIST(x0, x1);
1624 }
1625 if (ix == 0) {
1626 geo2cart(0.0, ret->lon[0][0], ret->lat[0][0], x0);
1627 geo2cart(0.0, ret->lon[ids][ip], ret->lat[ids][ip], x1);
1628 wave->y[iy] = DIST(x0, x1);
1629 }
1630
1631 /* Save geolocation... */
1632 wave->time = ret->time[0][0];
1633 if (ix == 0 && iy == 0)
1634 wave->z = ret->z[ids][ip];
1635 wave->lon[ix][iy] = ret->lon[ids][ip];
1636 wave->lat[ix][iy] = ret->lat[ids][ip];
1637
1638 /* Save temperature... */
1639 if (dataset == 1)
1640 wave->temp[ix][iy] = ret->t[ids][ip];
1641 else if (dataset == 2)
1642 wave->temp[ix][iy] = ret->t_apr[ids][ip];
1643 }
1644}
1645
1646/*****************************************************************************/
1647
1648void variance(
1649 wave_t *wave,
1650 double dh) {
1651
1652 /* Check parameters... */
1653 if (dh <= 0)
1654 return;
1655
1656 /* Compute squared radius... */
1657 const double dh2 = gsl_pow_2(dh);
1658
1659 /* Get sampling distances... */
1660 const int dx =
1661 (int) (dh / fabs(wave->x[wave->nx - 1] - wave->x[0]) * (wave->nx - 1.0) +
1662 1);
1663 const int dy =
1664 (int) (dh / fabs(wave->y[wave->ny - 1] - wave->y[0]) * (wave->ny - 1.0) +
1665 1);
1666
1667 /* Loop over data points... */
1668 for (int ix = 0; ix < wave->nx; ix++)
1669 for (int iy = 0; iy < wave->ny; iy++) {
1670
1671 /* Init... */
1672 double mu = 0, help = 0;
1673 int n = 0;
1674
1675 /* Get data... */
1676 for (int ix2 = GSL_MAX(ix - dx, 0);
1677 ix2 <= GSL_MIN(ix + dx, wave->nx - 1); ix2++)
1678 for (int iy2 = GSL_MAX(iy - dy, 0);
1679 iy2 <= GSL_MIN(iy + dy, wave->ny - 1); iy2++)
1680 if ((gsl_pow_2(wave->x[ix] - wave->x[ix2])
1681 + gsl_pow_2(wave->y[iy] - wave->y[iy2])) <= dh2)
1682 if (gsl_finite(wave->pt[ix2][iy2])) {
1683 mu += wave->pt[ix2][iy2];
1684 help += gsl_pow_2(wave->pt[ix2][iy2]);
1685 n++;
1686 }
1687
1688 /* Compute local variance... */
1689 if (n > 1)
1690 wave->var[ix][iy] = help / n - gsl_pow_2(mu / n);
1691 else
1692 wave->var[ix][iy] = GSL_NAN;
1693 }
1694}
1695
1696/*****************************************************************************/
1697
1698void write_wave(
1699 char *filename,
1700 wave_t *wave) {
1701
1702 FILE *out;
1703
1704 /* Write info... */
1705 LOG(1, "Write wave data: %s", filename);
1706
1707 /* Create file... */
1708 if (!(out = fopen(filename, "w")))
1709 ERRMSG("Cannot create file!");
1710
1711 /* Write header... */
1712 fprintf(out,
1713 "# $1 = time (seconds since 2000-01-01T00:00Z)\n"
1714 "# $2 = altitude [km]\n"
1715 "# $3 = longitude [deg]\n"
1716 "# $4 = latitude [deg]\n"
1717 "# $5 = across-track distance [km]\n"
1718 "# $6 = along-track distance [km]\n"
1719 "# $7 = temperature [K]\n"
1720 "# $8 = background [K]\n"
1721 "# $9 = perturbation [K]\n"
1722 "# $10 = variance [K^2]\n" "# $11 = fitting model [K]\n");
1723
1724 /* Write data... */
1725 for (int j = 0; j < wave->ny; j++) {
1726 fprintf(out, "\n");
1727 for (int i = 0; i < wave->nx; i++)
1728 fprintf(out, "%.2f %g %g %g %g %g %g %g %g %g %g\n",
1729 wave->time, wave->z, wave->lon[i][j], wave->lat[i][j],
1730 wave->x[i], wave->y[j], wave->temp[i][j], wave->bg[i][j],
1731 wave->pt[i][j], wave->var[i][j], wave->fit[i][j]);
1732 }
1733
1734 /* Close file... */
1735 fclose(out);
1736}
cart2geo
void cart2geo(const double *x, double *z, double *lon, double *lat)
Converts Cartesian coordinates to geographic coordinates.
Definition: jurassic.c:185
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
Definition: jurassic.c:3853
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:267
POW2
#define POW2(x)
Compute the square of a value.
Definition: jurassic.h:753
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: jurassic.h:948
WARN
#define WARN(...)
Print a warning message with contextual information.
Definition: jurassic.h:915
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: jurassic.h:878
DIST
#define DIST(a, b)
Compute Cartesian distance between two 3D vectors.
Definition: jurassic.h:447
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:548
merge_y
void merge_y(wave_t *wave1, wave_t *wave2)
Merge wave structs in y-direction.
Definition: libcris.c:694
fit_wave
void fit_wave(wave_t *wave, double amp, double phi, double kx, double ky, double *chisq)
Evaluate wave fit...
Definition: libcris.c:279
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:772
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:305
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:575
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:344
noise
void noise(wave_t *wave, double *mu, double *sig)
Estimate noise.
Definition: libcris.c:727
read_retr_help
void read_retr_help(double *help, int nds, int np, double mat[NDS][NPG])
Convert array.
Definition: libcris.c:1518
ret2wave
void ret2wave(retr_t *ret, wave_t *wave, int dataset, int ip)
Convert CrIS retrieval results to wave analysis struct.
Definition: libcris.c:1594
hamming
void hamming(wave_t *wave, int niter)
Apply Hamming filter to perturbations...
Definition: libcris.c:546
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:1051
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libcris.c:1648
read_retr
void read_retr(char *filename, retr_t *ret)
Read CrIS retrieval data.
Definition: libcris.c:1352
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
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:1221
read_wave
void read_wave(char *filename, wave_t *wave)
Read wave analysis data.
Definition: libcris.c:1533
gauss
void gauss(wave_t *wave, double fwhm)
Apply Gaussian filter to perturbations...
Definition: libcris.c:506
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:653
write_wave
void write_wave(char *filename, wave_t *wave)
Write wave analysis data.
Definition: libcris.c:1698
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:819
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libcris.h:61
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
retr_t
Retrieval results.
Definition: libcris.h:166
retr_t::lat
double lat[NDS][NPG]
Latitude [deg].
Definition: libcris.h:184
retr_t::t_noise
double t_noise[NDS][NPG]
Temperature (noise error) [K].
Definition: libcris.h:199
retr_t::p
double p[NDS][NPG]
Pressure [hPa].
Definition: libcris.h:187
retr_t::time
double time[NDS][NPG]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:175
retr_t::lon
double lon[NDS][NPG]
Longitude [deg].
Definition: libcris.h:181
retr_t::t_fm
double t_fm[NDS][NPG]
Temperature (forward model error) [K].
Definition: libcris.h:202
retr_t::np
int np
Number of data points.
Definition: libcris.h:172
retr_t::t_res
double t_res[NDS][NPG]
Temperature (resolution).
Definition: libcris.h:208
retr_t::z
double z[NDS][NPG]
Altitude [km].
Definition: libcris.h:178
retr_t::t_cont
double t_cont[NDS][NPG]
Temperature (measurement content).
Definition: libcris.h:205
retr_t::nds
int nds
Number of data sets.
Definition: libcris.h:169
retr_t::t_tot
double t_tot[NDS][NPG]
Temperature (total error) [K].
Definition: libcris.h:196
retr_t::t_apr
double t_apr[NDS][NPG]
Temperature (a priori data) [K].
Definition: libcris.h:193
retr_t::chisq
double chisq[NDS]
Chi^2.
Definition: libcris.h:211
retr_t::t
double t[NDS][NPG]
Temperature [K].
Definition: libcris.h:190
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 by doxygen 1.9.4