jurassic.h

Go to the documentation of this file.

/*
  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-2021 Forschungszentrum Juelich GmbH
*/
 
#ifndef JURASSIC_H
#define JURASSIC_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>
 
/* ------------------------------------------------------------
   Macros...
   ------------------------------------------------------------ */
 
#define ALLOC(ptr, type, n)                             \
  if((ptr=malloc((size_t)(n)*sizeof(type)))==NULL)      \
    ERRMSG("Out of memory!");
 
#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 POW2(x) ((x)*(x))
 
#define POW3(x) ((x)*(x)*(x))
 
#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);
 
/* ------------------------------------------------------------
   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 RE
#define RE 6367.421
#endif
 
#ifndef RI
#define RI 8.3144598
#endif
 
#ifndef P0
#define P0 1013.25
#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+5)
#endif
 
#ifndef NQ
#define NQ (7+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
 
#ifndef RFMLINE
#define RFMLINE 100000
#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 IDXSFZ (4+ctl->ng+ctl->nw+ctl->ncl)
 
#define IDXSFP (5+ctl->ng+ctl->nw+ctl->ncl)
 
#define IDXSFT (6+ctl->ng+ctl->nw+ctl->ncl)
 
#define IDXSFEPS(isf) (7+ctl->ng+ctl->nw+ctl->ncl+isf)
 
/* ------------------------------------------------------------
   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 sfz;
 
  double sfp;
 
  double sft;
 
  double sfeps[NSF];
 
} atm_t;
 
typedef struct {
 
  int ng;
 
  char emitter[NG][LEN];
 
  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 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_sfz;
 
  int ret_sfp;
 
  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 {
 
  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...
   ------------------------------------------------------------ */
 
size_t atm2x(
  ctl_t * ctl,
  atm_t * atm,
  gsl_vector * x,
  int *iqa,
  int *ipa);
 
void atm2x_help(
  double value,
  int value_iqa,
  int value_ip,
  gsl_vector * x,
  int *iqa,
  int *ipa,
  size_t *n);
 
double brightness(
  double rad,
  double nu);
 
void cart2geo(
  double *x,
  double *z,
  double *lon,
  double *lat);
 
void climatology(
  ctl_t * ctl,
  atm_t * atm_mean);
 
double ctmco2(
  double nu,
  double p,
  double t,
  double u);
 
double ctmh2o(
  double nu,
  double p,
  double t,
  double q,
  double u);
 
double ctmn2(
  double nu,
  double p,
  double t);
 
double ctmo2(
  double nu,
  double p,
  double t);
 
void copy_atm(
  ctl_t * ctl,
  atm_t * atm_dest,
  atm_t * atm_src,
  int init);
 
void copy_obs(
  ctl_t * ctl,
  obs_t * obs_dest,
  obs_t * obs_src,
  int init);
 
int find_emitter(
  ctl_t * ctl,
  const char *emitter);
 
void formod(
  ctl_t * ctl,
  atm_t * atm,
  obs_t * obs);
 
void formod_continua(
  ctl_t * ctl,
  los_t * los,
  int ip,
  double *beta);
 
void formod_fov(
  ctl_t * ctl,
  obs_t * obs);
 
void formod_pencil(
  ctl_t * ctl,
  atm_t * atm,
  obs_t * obs,
  int ir);
 
void formod_rfm(
  ctl_t * ctl,
  atm_t * atm,
  obs_t * obs);
 
void formod_srcfunc(
  ctl_t * ctl,
  tbl_t * tbl,
  double t,
  double *src);
 
void geo2cart(
  double z,
  double lon,
  double lat,
  double *x);
 
void hydrostatic(
  ctl_t * ctl,
  atm_t * atm);
 
void idx2name(
  ctl_t * ctl,
  int idx,
  char *quantity);
 
void init_srcfunc(
  ctl_t * ctl,
  tbl_t * tbl);
 
void intpol_atm(
  ctl_t * ctl,
  atm_t * atm,
  double z,
  double *p,
  double *t,
  double *q,
  double *k);
 
void intpol_tbl_cga(
  ctl_t * ctl,
  tbl_t * tbl,
  los_t * los,
  int ip,
  double tau_path[ND][NG],
  double tau_seg[ND]);
 
void intpol_tbl_ega(
  ctl_t * ctl,
  tbl_t * tbl,
  los_t * los,
  int ip,
  double tau_path[ND][NG],
  double tau_seg[ND]);
 
double intpol_tbl_eps(
  tbl_t * tbl,
  int ig,
  int id,
  int ip,
  int it,
  double u);
 
double intpol_tbl_u(
  tbl_t * tbl,
  int ig,
  int id,
  int ip,
  int it,
  double eps);
 
void jsec2time(
  double jsec,
  int *year,
  int *mon,
  int *day,
  int *hour,
  int *min,
  int *sec,
  double *remain);
 
void kernel(
  ctl_t * ctl,
  atm_t * atm,
  obs_t * obs,
  gsl_matrix * k);
 
int locate_irr(
  double *xx,
  int n,
  double x);
 
int locate_reg(
  double *xx,
  int n,
  double x);
 
int locate_tbl(
  float *xx,
  int n,
  double x);
 
size_t obs2y(
  ctl_t * ctl,
  obs_t * obs,
  gsl_vector * y,
  int *ida,
  int *ira);
 
double planck(
  double t,
  double nu);
 
void raytrace(
  ctl_t * ctl,
  atm_t * atm,
  obs_t * obs,
  los_t * los,
  int ir);
 
void read_atm(
  const char *dirname,
  const char *filename,
  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,
  ctl_t * ctl,
  obs_t * obs);
 
double read_obs_rfm(
  const char *basename,
  double z,
  double *nu,
  double *f,
  int n);
 
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);
 
void read_tbl(
  ctl_t * ctl,
  tbl_t * tbl);
 
double refractivity(
  double p,
  double t);
 
double scan_ctl(
  int argc,
  char *argv[],
  const char *varname,
  int arridx,
  const char *defvalue,
  char *value);
 
double sza(
  double sec,
  double lon,
  double lat);
 
void tangent_point(
  los_t * los,
  double *tpz,
  double *tplon,
  double *tplat);
 
void time2jsec(
  int year,
  int mon,
  int day,
  int hour,
  int min,
  int sec,
  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,
  ctl_t * ctl,
  atm_t * atm);
 
void write_atm_rfm(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
void write_matrix(
  const char *dirname,
  const char *filename,
  ctl_t * ctl,
  gsl_matrix * matrix,
  atm_t * atm,
  obs_t * obs,
  const char *rowspace,
  const char *colspace,
  const char *sort);
 
void write_obs(
  const char *dirname,
  const char *filename,
  ctl_t * ctl,
  obs_t * obs);
 
void write_shape(
  const char *filename,
  double *x,
  double *y,
  int n);
 
void write_tbl(
  ctl_t * ctl,
  tbl_t * tbl);
 
void x2atm(
  ctl_t * ctl,
  gsl_vector * x,
  atm_t * atm);
 
void x2atm_help(
  double *value,
  gsl_vector * x,
  size_t *n);
 
void y2obs(
  ctl_t * ctl,
  gsl_vector * y,
  obs_t * obs);
 
#endif