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