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