AIRS Code Collection
Data Structures | Macros | Functions
libairs.h File Reference

AIRS Code Collection library declarations. More...

#include <netcdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_spline.h>
#include <airs_rad_typ.h>
#include <airs_rad_struct.h>
#include <airs_ret_typ.h>
#include <airs_ret_struct.h>
#include "jurassic.h"

Go to the source code of this file.

Data Structures

struct  airs_l1_t
 AIRS Level-1 data. More...
 
struct  airs_l2_t
 AIRS Level-2 data. More...
 
struct  pert_t
 Perturbation data. More...
 
struct  ret_t
 Retrieval results. More...
 
struct  wave_t
 Wave analysis data. More...
 

Macros

#define NDS   200000
 Maximum number of data sets per granule. More...
 
#define NPG   30
 Maximum number of data points per granule. More...
 
#define L1_NCHAN   34
 Number of AIRS radiance channels (don't change). More...
 
#define L1_NTRACK   135
 Along-track size of AIRS radiance granule (don't change). More...
 
#define L1_NXTRACK   90
 Across-track size of AIRS radiance granule (don't change). More...
 
#define L2_NLAY   27
 Number of AIRS pressure layers (don't change). More...
 
#define L2_NTRACK   45
 Along-track size of AIRS retrieval granule (don't change). More...
 
#define L2_NXTRACK   30
 Across-track size of AIRS retrieval granule (don't change). More...
 
#define PERT_NTRACK   132000
 Along-track size of perturbation data. More...
 
#define PERT_NXTRACK   360
 Across-track size of perturbation data. More...
 
#define WX   300
 Across-track size of wave analysis data. More...
 
#define WY   33000
 Along-track size of wave analysis data. More...
 
#define PMAX   512
 Maximum number of data points for spectral analysis. More...
 
#define NC(cmd)
 Execute netCDF library command and check result. More...
 

Functions

void add_att (int ncid, int varid, const char *unit, const char *long_name)
 Add variable attributes to netCDF file. More...
 
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. More...
 
void background_poly (wave_t *wave, int dim_x, int dim_y)
 Get background based on polynomial fits. More...
 
void background_poly_help (double *xx, double *yy, int n, int dim)
 Get background based on polynomial fits. More...
 
void background_smooth (wave_t *wave, int npts_x, int npts_y)
 Smooth background. More...
 
void create_background (wave_t *wave)
 Set background... More...
 
void create_noise (wave_t *wave, double nedt)
 Add noise to perturbations and temperatures... More...
 
void create_wave (wave_t *wave, double amp, double lx, double ly, double phi, double fwhm)
 Add linear wave pattern... More...
 
void day2doy (int year, int mon, int day, int *doy)
 Get day of year from date. More...
 
void doy2day (int year, int doy, int *mon, int *day)
 Get date from day of year. More...
 
void fit_wave (wave_t *wave, double amp, double phi, double kx, double ky, double *chisq)
 Evaluate wave fit... More...
 
void fft_help (double *fcReal, double *fcImag, int n)
 Calculate 1-D FFT... More...
 
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... More...
 
void gauss (wave_t *wave, double fwhm)
 Apply Gaussian filter to perturbations... More...
 
void hamming (wave_t *wave, int nit)
 Apply Hamming filter to perturbations... More...
 
void intpol_x (wave_t *wave, int n)
 Interpolate to regular grid in x-direction. More...
 
void median (wave_t *wave, int dx)
 Apply median filter to perturbations... More...
 
void merge_y (wave_t *wave1, wave_t *wave2)
 Merge wave structs in y-direction. More...
 
void noise (wave_t *wave, double *mu, double *sig)
 Estimate noise. More...
 
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. More...
 
void pert2wave (pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
 Convert radiance perturbation data to wave analysis struct. More...
 
void read_l1 (char *filename, airs_l1_t *l1)
 Read AIRS Level-1 data. More...
 
void read_l2 (char *filename, airs_l2_t *l2)
 Read AIRS Level-2 data. More...
 
void read_pert (char *filename, char *pertname, pert_t *pert)
 Read radiance perturbation data. More...
 
void read_retr (char *filename, ret_t *ret)
 Read AIRS retrieval data. More...
 
void read_retr_help (double *help, int nds, int np, double mat[NDS][NPG])
 Convert array. More...
 
void read_wave (char *filename, wave_t *wave)
 Read wave analysis data. More...
 
void rad2wave (airs_rad_gran_t *airs_rad_gran, double *nu, int nd, wave_t *wave)
 Convert AIRS radiance data to wave analysis struct. More...
 
void ret2wave (ret_t *ret, wave_t *wave, int dataset, int ip)
 Convert AIRS retrieval results to wave analysis struct. More...
 
void variance (wave_t *wave, double dh)
 Compute local variance. More...
 
void write_l1 (char *filename, airs_l1_t *l1)
 Write AIRS Level-1 data. More...
 
void write_l2 (char *filename, airs_l2_t *l2)
 Write AIRS Level-2 data. More...
 
void write_wave (char *filename, wave_t *wave)
 Write wave analysis data. More...
 

Detailed Description

AIRS Code Collection library declarations.

Definition in file libairs.h.

Macro Definition Documentation

◆ NDS

#define NDS   200000

Maximum number of data sets per granule.

Definition at line 99 of file libairs.h.

◆ NPG

#define NPG   30

Maximum number of data points per granule.

Definition at line 102 of file libairs.h.

◆ L1_NCHAN

#define L1_NCHAN   34

Number of AIRS radiance channels (don't change).

Definition at line 105 of file libairs.h.

◆ L1_NTRACK

#define L1_NTRACK   135

Along-track size of AIRS radiance granule (don't change).

Definition at line 108 of file libairs.h.

◆ L1_NXTRACK

#define L1_NXTRACK   90

Across-track size of AIRS radiance granule (don't change).

Definition at line 111 of file libairs.h.

◆ L2_NLAY

#define L2_NLAY   27

Number of AIRS pressure layers (don't change).

Definition at line 114 of file libairs.h.

◆ L2_NTRACK

#define L2_NTRACK   45

Along-track size of AIRS retrieval granule (don't change).

Definition at line 117 of file libairs.h.

◆ L2_NXTRACK

#define L2_NXTRACK   30

Across-track size of AIRS retrieval granule (don't change).

Definition at line 120 of file libairs.h.

◆ PERT_NTRACK

#define PERT_NTRACK   132000

Along-track size of perturbation data.

Definition at line 123 of file libairs.h.

◆ PERT_NXTRACK

#define PERT_NXTRACK   360

Across-track size of perturbation data.

Definition at line 126 of file libairs.h.

◆ WX

#define WX   300

Across-track size of wave analysis data.

Definition at line 129 of file libairs.h.

◆ WY

#define WY   33000

Along-track size of wave analysis data.

Definition at line 132 of file libairs.h.

◆ PMAX

#define PMAX   512

Maximum number of data points for spectral analysis.

Definition at line 135 of file libairs.h.

◆ NC

#define NC (   cmd)
Value:
{ \
int nc_result=(cmd); \
if(nc_result!=NC_NOERR) \
ERRMSG("%s", nc_strerror(nc_result)); \
}

Execute netCDF library command and check result.

Definition at line 142 of file libairs.h.

Function Documentation

◆ add_att()

void add_att ( int  ncid,
int  varid,
const char *  unit,
const char *  long_name 
)

Add variable attributes to netCDF file.

Definition at line 30 of file libairs.c.

34 {
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}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35

◆ 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.

Add variable to netCDF file.

Definition at line 45 of file libairs.c.

53 {
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}

◆ background_poly()

void background_poly ( wave_t *  wave,
int  dim_x,
int  dim_y 
)

Get background based on polynomial fits.

Definition at line 126 of file libairs.c.

129 {
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}
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
WX
#define WX
Across-track size of wave analysis data.
Definition: libairs.h:129
WY
#define WY
Along-track size of wave analysis data.
Definition: libairs.h:132
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::temp
double temp[WX][WY]
Temperature [K].
Definition: libairs.h:314
wave_t::bg
double bg[WX][WY]
Background [K].
Definition: libairs.h:317
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libairs.h:320
Here is the call graph for this function:

◆ background_poly_help()

void background_poly_help ( double *  xx,
double *  yy,
int  n,
int  dim 
)

Get background based on polynomial fits.

Definition at line 72 of file libairs.c.

76 {
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}

◆ background_smooth()

void background_smooth ( wave_t *  wave,
int  npts_x,
int  npts_y 
)

Smooth background.

Definition at line 176 of file libairs.c.

179 {
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}
wave_t::x
double x[WX]
Across-track distance [km].
Definition: libairs.h:308
wave_t::y
double y[WY]
Along-track distance [km].
Definition: libairs.h:311

◆ create_background()

void create_background ( wave_t *  wave)

Set background...

Definition at line 227 of file libairs.c.

228 {
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}

◆ create_noise()

void create_noise ( wave_t *  wave,
double  nedt 
)

Add noise to perturbations and temperatures...

Definition at line 254 of file libairs.c.

256 {
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}

◆ create_wave()

void create_wave ( wave_t *  wave,
double  amp,
double  lx,
double  ly,
double  phi,
double  fwhm 
)

Add linear wave pattern...

Definition at line 275 of file libairs.c.

281 {
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}

◆ day2doy()

void day2doy ( int  year,
int  mon,
int  day,
int *  doy 
)

Get day of year from date.

Definition at line 304 of file libairs.c.

308 {
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}

◆ doy2day()

void doy2day ( int  year,
int  doy,
int *  mon,
int *  day 
)

Get date from day of year.

Definition at line 324 of file libairs.c.

328 {
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}

◆ fit_wave()

void fit_wave ( wave_t *  wave,
double  amp,
double  phi,
double  kx,
double  ky,
double *  chisq 
)

Evaluate wave fit...

Definition at line 354 of file libairs.c.

360 {
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}
POW2
#define POW2(x)
Compute x^2.
Definition: jurassic.h:187
wave_t::fit
double fit[WX][WY]
Fit [K].
Definition: libairs.h:326

◆ fft_help()

void fft_help ( double *  fcReal,
double *  fcImag,
int  n 
)

Calculate 1-D FFT...

Definition at line 380 of file libairs.c.

383 {
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}
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
PMAX
#define PMAX
Maximum number of data points for spectral analysis.
Definition: libairs.h:135

◆ 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 at line 419 of file libairs.c.

428 {
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}
fft_help
void fft_help(double *fcReal, double *fcImag, int n)
Calculate 1-D FFT...
Definition: libairs.c:380
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libairs.h:323
wave_t::lon
double lon[WX][WY]
Longitude [deg].
Definition: libairs.h:302
wave_t::lat
double lat[WX][WY]
Latitude [deg].
Definition: libairs.h:305
wave_t::z
double z
Altitude [km].
Definition: libairs.h:299
Here is the call graph for this function:

◆ gauss()

void gauss ( wave_t *  wave,
double  fwhm 
)

Apply Gaussian filter to perturbations...

Definition at line 579 of file libairs.c.

581 {
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}

◆ hamming()

void hamming ( wave_t *  wave,
int  nit 
)

Apply Hamming filter to perturbations...

Definition at line 619 of file libairs.c.

621 {
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}

◆ intpol_x()

void intpol_x ( wave_t *  wave,
int  n 
)

Interpolate to regular grid in x-direction.

Definition at line 648 of file libairs.c.

650 {
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}
cart2geo
void cart2geo(const double *x, double *z, double *lon, double *lat)
Convert Cartesian coordinates to geolocation.
Definition: jurassic.c:108
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Convert geolocation to Cartesian coordinates.
Definition: jurassic.c:3500
LIN
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
Definition: jurassic.h:164
Here is the call graph for this function:

◆ median()

void median ( wave_t *  wave,
int  dx 
)

Apply median filter to perturbations...

Definition at line 726 of file libairs.c.

728 {
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}

◆ merge_y()

void merge_y ( wave_t *  wave1,
wave_t *  wave2 
)

Merge wave structs in y-direction.

Definition at line 765 of file libairs.c.

767 {
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}

◆ noise()

void noise ( wave_t *  wave,
double *  mu,
double *  sig 
)

Estimate noise.

Definition at line 798 of file libairs.c.

801 {
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}

◆ 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 at line 842 of file libairs.c.

853 {
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}

◆ 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 at line 999 of file libairs.c.

1005 {
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}
DIST
#define DIST(a, b)
Compute Cartesian distance between two vectors.
Definition: jurassic.h:134
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::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
wave_t::time
double time
Time (seconds since 2000-01-01T00:00Z).
Definition: libairs.h:296
Here is the call graph for this function:

◆ read_l1()

void read_l1 ( char *  filename,
airs_l1_t *  l1 
)

Read AIRS Level-1 data.

Definition at line 1073 of file libairs.c.

1075 {
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}
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

◆ read_l2()

void read_l2 ( char *  filename,
airs_l2_t *  l2 
)

Read AIRS Level-2 data.

Definition at line 1107 of file libairs.c.

1109 {
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}
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

◆ read_pert()

void read_pert ( char *  filename,
char *  pertname,
pert_t *  pert 
)

Read radiance perturbation data.

Definition at line 1137 of file libairs.c.

1140 {
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}
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libairs.h:126
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libairs.h:123
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libairs.h:223

◆ read_retr()

void read_retr ( char *  filename,
ret_t *  ret 
)

Read AIRS retrieval data.

Definition at line 1219 of file libairs.c.

1221 {
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}
read_retr_help
void read_retr_help(double *help, int nds, int np, double mat[NDS][NPG])
Convert array.
Definition: libairs.c:1385
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
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
Here is the call graph for this function:

◆ read_retr_help()

void read_retr_help ( double *  help,
int  nds,
int  np,
double  mat[NDS][NPG] 
)

Convert array.

Definition at line 1385 of file libairs.c.

1389 {
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}

◆ read_wave()

void read_wave ( char *  filename,
wave_t *  wave 
)

Read wave analysis data.

Definition at line 1400 of file libairs.c.

1402 {
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}

◆ rad2wave()

void rad2wave ( airs_rad_gran_t *  airs_rad_gran,
double *  nu,
int  nd,
wave_t *  wave 
)

Convert AIRS radiance data to wave analysis struct.

Definition at line 1461 of file libairs.c.

1465 {
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}
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:126
Here is the call graph for this function:

◆ ret2wave()

void ret2wave ( ret_t *  ret,
wave_t *  wave,
int  dataset,
int  ip 
)

Convert AIRS retrieval results to wave analysis struct.

Definition at line 1531 of file libairs.c.

1535 {
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}
Here is the call graph for this function:

◆ variance()

void variance ( wave_t *  wave,
double  dh 
)

Compute local variance.

Definition at line 1584 of file libairs.c.

1586 {
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}

◆ write_l1()

void write_l1 ( char *  filename,
airs_l1_t *  l1 
)

Write AIRS Level-1 data.

Definition at line 1634 of file libairs.c.

1636 {
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}
L1_NXTRACK
#define L1_NXTRACK
Definition: diff_apr.c:52
L1_NTRACK
#define L1_NTRACK
Definition: diff_apr.c:49
L1_NCHAN
#define L1_NCHAN
Definition: diff_apr.c:46
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
Here is the call graph for this function:

◆ write_l2()

void write_l2 ( char *  filename,
airs_l2_t *  l2 
)

Write AIRS Level-2 data.

Definition at line 1692 of file libairs.c.

1694 {
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}
L2_NXTRACK
#define L2_NXTRACK
Definition: diff_apr.c:61
L2_NLAY
#define L2_NLAY
Definition: diff_apr.c:55
L2_NTRACK
#define L2_NTRACK
Definition: diff_apr.c:58
Here is the call graph for this function:

◆ write_wave()

void write_wave ( char *  filename,
wave_t *  wave 
)

Write wave analysis data.

Definition at line 1741 of file libairs.c.

1743 {
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}
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4