/github/workspace/src/jurassic.h

Analyze averaging kernel (AVK) matrix for retrieval diagnostics.

Analyze averaging kernel (AVK) matrix for retrieval diagnostics.Decomposes and evaluates the averaging kernel matrix \(\mathbf{A}\) to quantify the retrieval sensitivity and vertical resolution for each retrieved quantity (pressure, temperature, trace gases, extinction, etc.). The results are written to diagnostic atmospheric files for visualization.

Parameters

[in]	ret	Retrieval control and configuration structure (ret_t).
[in]	ctl	Global control structure (ctl_t) defining retrieval setup and quantities.
[in]	atm	Retrieved atmospheric state structure (atm_t).
[in]	iqa	Array mapping state vector indices to physical quantities (e.g., IDXP, IDXT, IDXQ(...)).
[in]	ipa	Array mapping state vector indices to atmospheric grid points.
[in]	avk	Averaging kernel matrix \(\mathbf{A}\) (gsl_matrix, n×n).

The averaging kernel matrix \(\mathbf{A}\) describes the sensitivity of the retrieved state \(\hat{x}\) to the true state \(x\):

\[ \hat{x} - x_a = \mathbf{A} (x - x_a) \]

where \(x_a\) is the a priori state.
This function separates and evaluates:

• Contribution profiles: diagonal dominance or total response (sensitivity).
• Resolution profiles: width or spread of the averaging kernels.

Internally, this function:

1. Identifies submatrices of \(\mathbf{A}\) corresponding to each physical quantity (pressure, temperature, VMRs, extinction, clouds, surface).
2. Calls analyze_avk_quantity() for each retrieved parameter type to compute quantitative measures of contribution and resolution.
3. Writes the results as atmospheric profiles:
   • atm_cont.tab — contribution functions (sensitivity).
   • atm_res.tab — resolution functions (vertical response width).

See also

analyze_avk_quantity, write_atm, set_cov_apr, set_cov_meas

Note

• Submatrices are identified by matching the quantity index iqa[i] to the quantity constants (e.g., IDXT, IDXQ(ig), etc.).
• Cloud and surface quantities are treated as scalar elements.
• Output files are written to the retrieval working directory (ret->dir).

gsl_matrix *A = gsl_matrix_alloc(n, n);
// ... compute averaging kernel ...
analyze_avk(&ret, &ctl, &atm, iqa, ipa, A);
// Creates atm_cont.tab and atm_res.tab for diagnostics

Warning

• The averaging kernel must be fully computed and dimensionally consistent with the state vector before calling this function.
• File writing will overwrite existing diagnostics in ret->dir.

@references Rodgers, C. D. (2000). Inverse Methods for Atmospheric Sounding: Theory and Practice. World Scientific.

Author

Lars Hoffmann

/*
  This file is part of JURASSIC.
  
  JURASSIC is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
  
  JURASSIC is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with JURASSIC. If not, see <http://www.gnu.org/licenses/>.
  
  Copyright (C) 2003-2025 Forschungszentrum Juelich GmbH
*/
 
#ifndef JURASSIC_H
#define JURASSIC_H
 
/* ------------------------------------------------------------
   Includes...
   ------------------------------------------------------------ */
 
#include <errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>
#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
 
/* ------------------------------------------------------------
   Constants...
   ------------------------------------------------------------ */
 
#ifndef C1
#define C1 1.19104259e-8
#endif
 
#ifndef C2
#define C2 1.43877506
#endif
 
#ifndef EPSMIN
#define EPSMIN 0
#endif
 
#ifndef EPSMAX
#define EPSMAX 1
#endif
 
#ifndef G0
#define G0 9.80665
#endif
 
#ifndef H0
#define H0 7.0
#endif
 
#ifndef KB
#define KB 1.3806504e-23
#endif
 
#ifndef ME
#define ME 5.976e24
#endif
 
#ifndef NA
#define NA 6.02214199e23
#endif
 
#ifndef N2
#define N2 0.78084
#endif
 
#ifndef O2
#define O2 0.20946
#endif
 
#ifndef P0
#define P0 1013.25
#endif
 
#ifndef RE
#define RE 6367.421
#endif
 
#ifndef RI
#define RI 8.3144598
#endif
 
#ifndef T0
#define T0 273.15
#endif
 
#ifndef TMIN
#define TMIN 100.
#endif
 
#ifndef TMAX
#define TMAX 400.
#endif
 
#ifndef TSUN
#define TSUN 5780.
#endif
 
#ifndef UMIN
#define UMIN 0
#endif
 
#ifndef UMAX
#define UMAX 1e30
#endif
 
/* ------------------------------------------------------------
   Dimensions...
   ------------------------------------------------------------ */
 
#ifndef NCL
#define NCL 8
#endif
 
#ifndef ND
#define ND 128
#endif
 
#ifndef NG
#define NG 8
#endif
 
#ifndef NP
#define NP 256
#endif
 
#ifndef NR
#define NR 256
#endif
 
#ifndef NSF
#define NSF 8
#endif
 
#ifndef NW
#define NW 4
#endif
 
#ifndef LEN
#define LEN 10000
#endif
 
#ifndef M
#define M (NR*ND)
#endif
 
#ifndef N
#define N ((2 + NG + NW) * NP + NCL + NSF + 3)
#endif
 
#ifndef NQ
#define NQ (5 + NG + NW + NCL + NSF)
#endif
 
#ifndef NLOS
#define NLOS 4096
#endif
 
#ifndef NSHAPE
#define NSHAPE 20000
#endif
 
#ifndef NFOV
#define NFOV 5
#endif
 
#ifndef TBLNP
#define TBLNP 41
#endif
 
#ifndef TBLNT
#define TBLNT 30
#endif
 
#ifndef TBLNU
#define TBLNU 320
#endif
 
#ifndef TBLNS
#define TBLNS 1200
#endif
 
#ifndef RFMNPTS
#define RFMNPTS 10000000
#endif
 
/* ------------------------------------------------------------
   Quantity indices...
   ------------------------------------------------------------ */
 
#define IDXP 0
 
#define IDXT 1
 
#define IDXQ(ig) (2 + (ig))
 
#define IDXK(iw) (2 + (ctl->ng) + (iw))
 
#define IDXCLZ (2 + (ctl->ng) + (ctl->nw))
 
#define IDXCLDZ (3 + (ctl->ng) + (ctl->nw))
 
#define IDXCLK(icl) (4 + (ctl->ng) + (ctl->nw) + (icl))
 
#define IDXSFT (4 + (ctl->ng) + (ctl->nw) + (ctl->ncl))
 
#define IDXSFEPS(isf) (5 + (ctl->ng) + (ctl->nw) + (ctl->ncl) + (isf))
 
/* ------------------------------------------------------------
   Macros...
   ------------------------------------------------------------ */
 
#define ALLOC(ptr, type, n) \
  if((ptr=malloc((size_t)(n)*sizeof(type)))==NULL) \
    ERRMSG("Out of memory!");
 
#define BRIGHT(rad, nu) \
  (C2 * (nu) / gsl_log1p(C1 * POW3(nu) / (rad)))
 
#define DEG2RAD(deg) ((deg) * (M_PI / 180.0))
 
#define DIST(a, b) sqrt(DIST2(a, b))
 
#define DIST2(a, b) \
  ((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])+(a[2]-b[2])*(a[2]-b[2]))
 
#define DOTP(a, b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
 
#define FREAD(ptr, type, size, out) { \
    if(fread(ptr, sizeof(type), size, out)!=size) \
      ERRMSG("Error while reading!"); \
  }
 
#define FWRITE(ptr, type, size, out) { \
    if(fwrite(ptr, sizeof(type), size, out)!=size) \
      ERRMSG("Error while writing!"); \
  }
 
#define MAX(a,b) (((a)>(b))?(a):(b))
 
#define MIN(a,b) (((a)<(b))?(a):(b))
 
#define LIN(x0, y0, x1, y1, x) \
  ((y0)+((y1)-(y0))/((x1)-(x0))*((x)-(x0)))
 
#define LOGX(x0, y0, x1, y1, x) \
  (((x)/(x0)>0 && (x1)/(x0)>0) \
   ? ((y0)+((y1)-(y0))*log((x)/(x0))/log((x1)/(x0))) \
   : LIN(x0, y0, x1, y1, x))
 
#define LOGY(x0, y0, x1, y1, x) \
  (((y1)/(y0)>0) \
   ? ((y0)*exp(log((y1)/(y0))/((x1)-(x0))*((x)-(x0)))) \
   : LIN(x0, y0, x1, y1, x))
 
#define NORM(a) sqrt(DOTP(a, a))
 
#define PLANCK(T, nu) \
  (C1 * POW3(nu) / gsl_expm1(C2 * (nu) / (T)))
 
#define POW2(x) ((x)*(x))
 
#define POW3(x) ((x)*(x)*(x))
 
#define RAD2DEG(rad) ((rad) * (180.0 / M_PI))
 
#define REFRAC(p, T) (7.753e-05 * (p) / (T))
 
#define TIMER(name, mode) \
  {timer(name, __FILE__, __func__, __LINE__, mode);}
 
#define TOK(line, tok, format, var) { \
    if(((tok)=strtok((line), " \t"))) { \
      if(sscanf(tok, format, &(var))!=1) continue; \
    } else ERRMSG("Error while reading!"); \
  }
 
/* ------------------------------------------------------------
   Log messages...
   ------------------------------------------------------------ */
 
#ifndef LOGLEV
#define LOGLEV 2
#endif
 
#define LOG(level, ...) {                                               \
    if(level >= 2)                                                      \
      printf("  ");                                                     \
    if(level <= LOGLEV) {                                               \
      printf(__VA_ARGS__);                                              \
      printf("\n");                                                     \
    }                                                                   \
  }
 
#define WARN(...) {                                                     \
    printf("\nWarning (%s, %s, l%d): ", __FILE__, __func__, __LINE__);  \
    LOG(0, __VA_ARGS__);                                                \
  }
 
#define ERRMSG(...) {                                                   \
    printf("\nError (%s, %s, l%d): ", __FILE__, __func__, __LINE__);    \
    LOG(0, __VA_ARGS__);                                                \
    exit(EXIT_FAILURE);                                                 \
  }
 
#define PRINT(format, var)                                              \
  printf("Print (%s, %s, l%d): %s= "format"\n",                         \
         __FILE__, __func__, __LINE__, #var, var);
 
/* ------------------------------------------------------------
   Structs...
   ------------------------------------------------------------ */
 
typedef struct {
 
  int np;
 
  double time[NP];
 
  double z[NP];
 
  double lon[NP];
 
  double lat[NP];
 
  double p[NP];
 
  double t[NP];
 
  double q[NG][NP];
 
  double k[NW][NP];
 
  double clz;
 
  double cldz;
 
  double clk[NCL];
 
  double sft;
 
  double sfeps[NSF];
 
} atm_t;
 
typedef struct {
 
  int ng;
 
  char emitter[NG][LEN];
 
  int ig_co2;
 
  int ig_h2o;
 
  int ig_n2;
 
  int ig_o2;
 
  int nd;
 
  double nu[ND];
 
  int nw;
 
  int window[ND];
 
  int ncl;
 
  double clnu[NCL];
 
  int nsf;
 
  double sfnu[NSF];
 
  int sftype;
 
  double sfsza;
 
  char tblbase[LEN];
 
  int tblfmt;
 
  double hydz;
 
  int ctm_co2;
 
  int ctm_h2o;
 
  int ctm_n2;
 
  int ctm_o2;
 
  int refrac;
 
  double rayds;
 
  double raydz;
 
  char fov[LEN];
 
  double fov_dz[NSHAPE];
 
  double fov_w[NSHAPE];
 
  int fov_n;
 
  double retp_zmin;
 
  double retp_zmax;
 
  double rett_zmin;
 
  double rett_zmax;
 
  double retq_zmin[NG];
 
  double retq_zmax[NG];
 
  double retk_zmin[NW];
 
  double retk_zmax[NW];
 
  int ret_clz;
 
  int ret_cldz;
 
  int ret_clk;
 
  int ret_sft;
 
  int ret_sfeps;
 
  int write_bbt;
 
  int write_matrix;
 
  int formod;
 
  char rfmbin[LEN];
 
  char rfmhit[LEN];
 
  char rfmxsc[NG][LEN];
 
} ctl_t;
 
typedef struct {
 
  int np;
 
  double z[NLOS];
 
  double lon[NLOS];
 
  double lat[NLOS];
 
  double p[NLOS];
 
  double t[NLOS];
 
  double q[NLOS][NG];
 
  double k[NLOS][ND];
 
  double sft;
 
  double sfeps[ND];
 
  double ds[NLOS];
 
  double u[NLOS][NG];
 
  double cgp[NLOS][NG];
 
  double cgt[NLOS][NG];
 
  double cgu[NLOS][NG];
 
  double eps[NLOS][ND];
 
  double src[NLOS][ND];
 
} los_t;
 
typedef struct {
 
  int nr;
 
  double time[NR];
 
  double obsz[NR];
 
  double obslon[NR];
 
  double obslat[NR];
 
  double vpz[NR];
 
  double vplon[NR];
 
  double vplat[NR];
 
  double tpz[NR];
 
  double tplon[NR];
 
  double tplat[NR];
 
  double tau[ND][NR];
 
  double rad[ND][NR];
 
} obs_t;
 
typedef struct {
 
  char dir[LEN];
 
  int kernel_recomp;
 
  int conv_itmax;
 
  double conv_dmin;
 
  int err_ana;
 
  double err_formod[ND];
 
  double err_noise[ND];
 
  double err_press;
 
  double err_press_cz;
 
  double err_press_ch;
 
  double err_temp;
 
  double err_temp_cz;
 
  double err_temp_ch;
 
  double err_q[NG];
 
  double err_q_cz[NG];
 
  double err_q_ch[NG];
 
  double err_k[NW];
 
  double err_k_cz[NW];
 
  double err_k_ch[NW];
 
  double err_clz;
 
  double err_cldz;
 
  double err_clk[NCL];
 
  double err_sft;
 
  double err_sfeps[NSF];
 
} ret_t;
 
typedef struct {
 
  int np[ND][NG];
 
  int nt[ND][NG][TBLNP];
 
  int nu[ND][NG][TBLNP][TBLNT];
 
  double p[ND][NG][TBLNP];
 
  double t[ND][NG][TBLNP][TBLNT];
 
  float u[ND][NG][TBLNP][TBLNT][TBLNU];
 
  float eps[ND][NG][TBLNP][TBLNT][TBLNU];
 
  double st[TBLNS];
 
  double sr[TBLNS][ND];
 
} tbl_t;
 
/* ------------------------------------------------------------
   Functions...
   ------------------------------------------------------------ */
 
void analyze_avk(
  ret_t * ret,
  ctl_t * ctl,
  atm_t * atm,
  int *iqa,
  int *ipa,
  gsl_matrix * avk);
 
void analyze_avk_quantity(
  gsl_matrix * avk,
  int iq,
  int *ipa,
  size_t *n0,
  size_t *n1,
  double *cont,
  double *res);
 
size_t atm2x(
  const ctl_t * ctl,
  const atm_t * atm,
  gsl_vector * x,
  int *iqa,
  int *ipa);
 
void atm2x_help(
  const double value,
  const int value_iqa,
  const int value_ip,
  gsl_vector * x,
  int *iqa,
  int *ipa,
  size_t *n);
 
void cart2geo(
  const double *x,
  double *z,
  double *lon,
  double *lat);
 
void climatology(
  const ctl_t * ctl,
  atm_t * atm_mean);
 
double cost_function(
  gsl_vector * dx,
  gsl_vector * dy,
  gsl_matrix * s_a_inv,
  gsl_vector * sig_eps_inv);
 
double ctmco2(
  const double nu,
  const double p,
  const double t,
  const double u);
 
double ctmh2o(
  const double nu,
  const double p,
  const double t,
  const double q,
  const double u);
 
double ctmn2(
  const double nu,
  const double p,
  const double t);
 
double ctmo2(
  const double nu,
  const double p,
  const double t);
 
void copy_atm(
  const ctl_t * ctl,
  atm_t * atm_dest,
  const atm_t * atm_src,
  const int init);
 
void copy_obs(
  const ctl_t * ctl,
  obs_t * obs_dest,
  const obs_t * obs_src,
  const int init);
 
int find_emitter(
  const ctl_t * ctl,
  const char *emitter);
 
void formod(
  const ctl_t * ctl,
  const tbl_t * tbl,
  atm_t * atm,
  obs_t * obs);
 
void formod_continua(
  const ctl_t * ctl,
  const los_t * los,
  const int ip,
  double *beta);
 
void formod_fov(
  const ctl_t * ctl,
  obs_t * obs);
 
void formod_pencil(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const atm_t * atm,
  obs_t * obs,
  const int ir);
 
void formod_rfm(
  const ctl_t * ctl,
  const atm_t * atm,
  obs_t * obs);
 
void formod_srcfunc(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const double t,
  double *src);
 
void geo2cart(
  const double z,
  const double lon,
  const double lat,
  double *x);
 
void hydrostatic(
  const ctl_t * ctl,
  atm_t * atm);
 
void idx2name(
  const ctl_t * ctl,
  const int idx,
  char *quantity);
 
void init_srcfunc(
  const ctl_t * ctl,
  tbl_t * tbl);
 
void intpol_atm(
  const ctl_t * ctl,
  const atm_t * atm,
  const double z,
  double *p,
  double *t,
  double *q,
  double *k);
 
void intpol_tbl_cga(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const los_t * los,
  const int ip,
  double tau_path[ND][NG],
  double tau_seg[ND]);
 
void intpol_tbl_ega(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const los_t * los,
  const int ip,
  double tau_path[ND][NG],
  double tau_seg[ND]);
 
double intpol_tbl_eps(
  const tbl_t * tbl,
  const int ig,
  const int id,
  const int ip,
  const int it,
  const double u);
 
double intpol_tbl_u(
  const tbl_t * tbl,
  const int ig,
  const int id,
  const int ip,
  const int it,
  const double eps);
 
void jsec2time(
  const double jsec,
  int *year,
  int *mon,
  int *day,
  int *hour,
  int *min,
  int *sec,
  double *remain);
 
void kernel(
  const ctl_t * ctl,
  const tbl_t * tbl,
  atm_t * atm,
  obs_t * obs,
  gsl_matrix * k);
 
int locate_irr(
  const double *xx,
  const int n,
  const double x);
 
int locate_reg(
  const double *xx,
  const int n,
  const double x);
 
int locate_tbl(
  const float *xx,
  const int n,
  const double x);
 
void matrix_invert(
  gsl_matrix * a);
 
void matrix_product(
  gsl_matrix * a,
  gsl_vector * b,
  int transpose,
  gsl_matrix * c);
 
size_t obs2y(
  const ctl_t * ctl,
  const obs_t * obs,
  gsl_vector * y,
  int *ida,
  int *ira);
 
void raytrace(
  const ctl_t * ctl,
  const atm_t * atm,
  obs_t * obs,
  los_t * los,
  const int ir);
 
void read_atm(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
void read_ctl(
  int argc,
  char *argv[],
  ctl_t * ctl);
 
void read_matrix(
  const char *dirname,
  const char *filename,
  gsl_matrix * matrix);
 
void read_obs(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  obs_t * obs);
 
double read_obs_rfm(
  const char *basename,
  const double z,
  double *nu,
  double *f,
  int n);
 
void read_ret(
  int argc,
  char *argv[],
  ctl_t * ctl,
  ret_t * ret);
 
void read_rfm_spec(
  const char *filename,
  double *nu,
  double *rad,
  int *npts);
 
void read_shape(
  const char *filename,
  double *x,
  double *y,
  int *n);
 
tbl_t *read_tbl(
  const ctl_t * ctl);
 
double scan_ctl(
  int argc,
  char *argv[],
  const char *varname,
  const int arridx,
  const char *defvalue,
  char *value);
 
void set_cov_apr(
  ret_t * ret,
  ctl_t * ctl,
  atm_t * atm,
  int *iqa,
  int *ipa,
  gsl_matrix * s_a);
 
void set_cov_meas(
  ret_t * ret,
  ctl_t * ctl,
  obs_t * obs,
  gsl_vector * sig_noise,
  gsl_vector * sig_formod,
  gsl_vector * sig_eps_inv);
 
double sza(
  double sec,
  double lon,
  double lat);
 
void tangent_point(
  const los_t * los,
  double *tpz,
  double *tplon,
  double *tplat);
 
void time2jsec(
  const int year,
  const int mon,
  const int day,
  const int hour,
  const int min,
  const int sec,
  const double remain,
  double *jsec);
 
void timer(
  const char *name,
  const char *file,
  const char *func,
  int line,
  int mode);
 
void write_atm(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_atm_rfm(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_matrix(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const gsl_matrix * matrix,
  const atm_t * atm,
  const obs_t * obs,
  const char *rowspace,
  const char *colspace,
  const char *sort);
 
void write_obs(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const obs_t * obs);
 
void write_shape(
  const char *filename,
  const double *x,
  const double *y,
  const int n);
 
void write_stddev(
  const char *quantity,
  ret_t * ret,
  ctl_t * ctl,
  atm_t * atm,
  gsl_matrix * s);
 
void write_tbl(
  const ctl_t * ctl,
  const tbl_t * tbl);
 
void x2atm(
  const ctl_t * ctl,
  const gsl_vector * x,
  atm_t * atm);
 
void x2atm_help(
  double *value,
  const gsl_vector * x,
  size_t *n);
 
void y2obs(
  const ctl_t * ctl,
  const gsl_vector * y,
  obs_t * obs);
 
#endif