26#include <gsl/gsl_fit.h>
48 static atm_t atm, atm2;
52 static FILE *in, *out;
54 static char line[
LEN];
57 obs_sim, scl = 1.0, scl_err, c0, c1, cov00, cov01, cov11, sumsq,
64 static int i, ig, n, nl, ndata[
NMAX], nprof;
68 ERRMSG(
"Give parameters: <ctl> <prof> <inv> <atm> <kernel>");
72 const double dt =
scan_ctl(argc, argv,
"INVERT_DT", -1,
"86400", NULL);
73 const double obs_err =
74 scan_ctl(argc, argv,
"INVERT_OBS_ERR", -1,
"1.0", NULL);
75 const int data = (int)
scan_ctl(argc, argv,
"INVERT_DATA", -1,
"2", NULL);
76 const int fit = (int)
scan_ctl(argc, argv,
"INVERT_FIT", -1,
"3", NULL);
78 (int)
scan_ctl(argc, argv,
"INVERT_ITMAX", -1,
"10", NULL);
79 const double tol =
scan_ctl(argc, argv,
"INVERT_TOL", -1,
"1e-4", NULL);
84 if (strcmp(ctl.
emitter[0],
"SO2") != 0)
85 ERRMSG(
"Set EMITTER[0] = SO2!");
86 if (strcmp(ctl.
emitter[1],
"H2O") != 0)
87 ERRMSG(
"Set EMITTER[1] = H2O!");
88 if (strcmp(ctl.
emitter[2],
"O3") != 0)
89 ERRMSG(
"Set EMITTER[2] = O3!");
90 if (strcmp(ctl.
emitter[3],
"CO2") != 0)
91 ERRMSG(
"Set EMITTER[3] = CO2!");
108 LOG(1,
"Read profile data: %s", argv[2]);
111 if (!(in = fopen(argv[2],
"r")))
112 ERRMSG(
"Cannot open file!");
115 while (fgets(line,
LEN, in)) {
118 if (sscanf(line,
"%lg %lg %lg %lg %g %g %g %g %g %g", &rtime[nl],
119 &rz[nl], &rlon[nl], &rlat[nl], &rp[nl], &rt[nl], &rso2[nl],
120 &rh2o[nl], &ro3[nl], &robs[nl]) != 10)
123 ERRMSG(
"Too many profile data points!");
134 for (
int it = 0; it < itmax; it++) {
138 for (i = 0; i <
NMAX; i++) {
144 for (
int il = 0; il < nl; il++) {
147 if ((rtime[il] != atm.
time[0]
148 || rlon[il] != atm.
lon[0]
149 || rlat[il] != atm.
lat[0])
154 obs_sim = obs.
rad[0][0] - obs.
rad[1][0];
157 i = (int) ((atm.
time[0] - rtime[0]) / dt);
158 if (i < 0 && i >=
NMAX)
159 ERRMSG(
"Time index out of range!");
163 x[i] = (isfinite(x[i]) ?
MAX(x[i], obs_sim) : obs_sim);
164 y[i] = (isfinite(y[i]) ?
MAX(y[i], obs_meas) : obs_meas);
166 if (isfinite(x[i]) && isfinite(y[i]))
171 else if (data == 2) {
175 y_err[i] =
POW2(obs_meas);
179 y_err[i] +=
POW2(obs_meas);
187 for (
int ip = 0; ip < atm.
np; ip++) {
189 atm2.
z[ip] += atm.
z[ip];
190 atm2.
lon[ip] += atm.
lon[ip];
191 atm2.
lat[ip] += atm.
lat[ip];
192 atm2.
p[ip] += atm.
p[ip];
193 atm2.
t[ip] += atm.
t[ip];
194 for (ig = 1; ig < ctl.
ng; ig++)
195 atm2.
q[ig][ip] += atm.
q[ig][ip];
204 atm.
time[atm.
np] = rtime[il];
205 atm.
z[atm.
np] = rz[il];
206 atm.
lon[atm.
np] = rlon[il];
207 atm.
lat[atm.
np] = rlat[il];
208 atm.
p[atm.
np] = rp[il];
209 atm.
t[atm.
np] = rt[il];
210 atm.
q[0][atm.
np] = rso2[il] * scl;
211 atm.
q[1][atm.
np] = rh2o[il];
212 atm.
q[2][atm.
np] = ro3[il];
213 atm.
q[3][atm.
np] = 371.789948e-6 + 2.026214e-6
214 * (atm.
time[atm.
np] - 63158400.) / 31557600.;
216 ERRMSG(
"Too many data points!");
221 for (i = 0; i <
NMAX; i++)
225 y_err[i] = sqrt(
MAX(y_err[i] / ndata[i] -
POW2(y[i]), 0.0))
231 for (i = 0; i <
NMAX; i++)
232 if (ndata[i] > 0 && isfinite(x[i]) && isfinite(y[i])
233 && isfinite(y_err[i])) {
236 y2_err[n] = y_err[i];
237 w2[n] = 1. /
POW2(y_err[i]);
243 gsl_fit_mul(x2, 1, y2, 1, (
size_t) n, &c1, &cov11, &sumsq);
245 gsl_fit_wmul(x2, 1, w2, 1, y2, 1, (
size_t) n, &c1, &cov11, &sumsq);
247 gsl_fit_linear(x2, 1, y2, 1, (
size_t) n, &c0, &c1, &cov00, &cov01,
250 gsl_fit_wlinear(x2, 1, w2, 1, y2, 1, (
size_t) n, &c0, &c1, &cov00,
251 &cov01, &cov11, &sumsq);
253 ERRMSG(
"Check INVERT_FIT!");
256 double scl_old = scl;
257 scl_err = scl * sqrt(cov11);
261 LOG(1,
" it= %d | scl= %g +/- %g | RMSE= %g",
262 it, scl, scl_err, sqrt(sumsq / n));
265 if (fabs(2.0 * (scl - scl_old) / (scl + scl_old)) < tol)
274 LOG(1,
"Write inversion data: %s", argv[3]);
277 if (!(out = fopen(argv[3],
"w")))
278 ERRMSG(
"Cannot create file!");
282 "# $1 = time (seconds since 2000-01-01T00:00Z)\n"
283 "# $2 = simulated SO2 index [K]\n"
284 "# $3 = scaled simulated SO2 index [K]\n"
285 "# $4 = error of scaled simulated SO2 index [K]\n"
286 "# $5 = observed SO2 index [K]\n"
287 "# $6 = error of observed SO2 index [K]\n\n");
290 for (i = 0; i < n; i++) {
293 if (fit == 1 || fit == 2)
294 gsl_fit_mul_est(x2[i], c1, cov11, &y2_sim[i], &y2_sim_err[i]);
295 else if (fit == 3 || fit == 4)
296 gsl_fit_linear_est(x2[i], c0, c1, cov00, cov01, cov11, &y2_sim[i],
300 fprintf(out,
"%.2f %g %g %g %g %g\n", rtime[0] + (i + 0.5) * dt,
301 x2[i], y2_sim[i], y2_sim_err[i], y2[i], y2_err[i]);
306 fprintf(out,
"# scl= %g +/- %g\n", scl, scl_err);
307 fprintf(out,
"# c1= %g +/- %g\n", c1, sqrt(cov11));
308 if (fit == 3 || fit == 4) {
309 fprintf(out,
"# c0= %g +/- %g\n", c0, sqrt(cov00));
310 fprintf(out,
"# corr= %g\n", cov01 / (sqrt(cov00) * sqrt(cov11)));
312 fprintf(out,
"# RMSE= %g\n", sqrt(sumsq / n));
313 fprintf(out,
"# n= %d\n", n);
323 for (
int ip = 0; ip < atm2.
np; ip++) {
324 atm2.
time[ip] /= nprof;
326 atm2.
lon[ip] /= nprof;
327 atm2.
lat[ip] /= nprof;
330 for (ig = 0; ig < ctl.
ng; ig++)
331 atm2.
q[ig][ip] /= nprof;
335 const size_t nk =
atm2x(&ctl, &atm2, NULL, NULL, NULL);
336 const size_t mk =
obs2y(&ctl, &obs, NULL, NULL, NULL);
339 gsl_matrix *k = gsl_matrix_alloc(mk, nk);
342 kernel(&ctl, &atm2, &obs, k);
348 write_matrix(NULL, argv[5], &ctl, k, &atm2, &obs,
"y",
"x",
"r");
int main(int argc, char *argv[])
#define NMAX
Maximum number of data points...
#define NLMAX
Maximum number of data lines...
void write_atm(const char *dirname, const char *filename, const ctl_t *ctl, const atm_t *atm)
Write atmospheric data.
void read_ctl(int argc, char *argv[], ctl_t *ctl)
Read forward model control parameters.
void formod(const ctl_t *ctl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Compose measurement vector.
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Compose state vector or parameter vector.
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)
Write matrix.
void kernel(ctl_t *ctl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
JURASSIC library declarations.
#define LEN
Maximum length of ASCII data lines.
#define POW2(x)
Compute x^2.
#define ERRMSG(...)
Print error message and quit program.
#define NP
Maximum number of atmospheric data points.
#define LOG(level,...)
Print log message.
#define MAX(a, b)
Macro to determine the maximum of two values.
double time[NP]
Time (seconds since 2000-01-01T00:00Z).
double lat[NP]
Latitude [deg].
double lon[NP]
Longitude [deg].
double t[NP]
Temperature [K].
int np
Number of data points.
double z[NP]
Altitude [km].
double q[NG][NP]
Volume mixing ratio [ppv].
double p[NP]
Pressure [hPa].
Forward model control parameters.
int write_matrix
Write matrix file (0=no, 1=yes).
int ng
Number of emitters.
int nd
Number of radiance channels.
char emitter[NG][LEN]
Name of each emitter.
int write_bbt
Use brightness temperature instead of radiance (0=no, 1=yes).
Observation geometry and radiance data.
double rad[ND][NR]
Radiance [W/(m^2 sr cm^-1)].
double obsz[NR]
Observer altitude [km].
int nr
Number of ray paths.