CrIS Code Collection
volcano.c
Go to the documentation of this file.
1#include "libcris.h"
2
3/* ------------------------------------------------------------
4 Functions...
5 ------------------------------------------------------------ */
6
7/* Estimate noise. */
8double get_noise(
9 double rad,
10 double nesr,
11 double nu);
12
13/* ------------------------------------------------------------
14 Main...
15 ------------------------------------------------------------ */
16
17int main(
18 int argc,
19 char *argv[]) {
20
21 FILE *out;
22
23 static cris_l1_t l1;
24
25 const int ai_low_nu1 = 341, ai_low_nu2 = 499, ai_high_nu1 =
26 40, ai_high_nu2 = 698, ai_old_nu1 = 295, ai_old_nu2 = 499, ci_nu =
27 36, si_low_nu1 = 326, si_low_nu2 = 260, si_high_nu1 = 328, si_high_nu2 =
28 282, si_old_nu1 = 318, si_old_nu2 = 260, si_oper_nu1 = 359, si_oper_nu2 =
29 244;
30
31 /* Check arguments... */
32 if (argc < 4)
33 ERRMSG("Give parameters: <ctl> <out.tab> <l1b_file1> [<l1b_file2> ...]");
34
35 /* Get control parameters... */
36 const int apo = (int) scan_ctl(argc, argv, "APO", -1, "0", NULL);
37
38 /* Create file... */
39 printf("Write volcanic emission data: %s\n", argv[2]);
40 if (!(out = fopen(argv[2], "w")))
41 ERRMSG("Cannot create file!");
42
43 /* Loop over CrIS files... */
44 for (int iarg = 3; iarg < argc; iarg++) {
45
46 /* Read CrIS data... */
47 if (!read_cris_l1(argv[iarg], &l1, apo))
48 continue;
49
50 /* Write header... */
51 if (iarg == 3) {
52 fprintf(out,
53 "# $1 = time [s]\n"
54 "# $2 = footprint longitude [deg]\n"
55 "# $3 = footprint latitude [deg]\n"
56 "# $4 = satellite altitude [km]\n"
57 "# $5 = satellite longitude [deg]\n"
58 "# $6 = satellite latitude [deg]\n");
59 fprintf(out,
60 "# $7 = cloud index, BT(%.4f/cm) [K]\n"
61 "# $8 = cloud index error [K]\n"
62 "# $9 = ash index (low wavenumbers),"
63 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
64 "# $10 = ash index (low wavenumbers) error [K]\n"
65 "# $11 = ash index (high wavenumbers),"
66 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
67 "# $12 = ash index (high wavenumbers) error [K]\n"
68 "# $13 = ash index (Hoffmann et al., 2014),"
69 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
70 "# $14 = ash index (Hoffmann et al., 2014) error [K]\n",
71 l1.nu_mw[ci_nu],
72 l1.nu_lw[ai_low_nu1], l1.nu_lw[ai_low_nu2],
73 l1.nu_mw[ai_high_nu1], l1.nu_lw[ai_high_nu2],
74 l1.nu_lw[ai_old_nu1], l1.nu_lw[ai_old_nu2]);
75 fprintf(out,
76 "# $15 = SO2 index (low concentrations),"
77 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
78 "# $16 = SO2 index (low concentrations) error [K]\n"
79 "# $17 = SO2 index (high concentrations),"
80 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
81 "# $18 = SO2 index (high concentrations) error [K]\n"
82 "# $19 = SO2 index (operational),"
83 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
84 "# $20 = SO2 index (operational) error [K]\n"
85 "# $21 = SO2 index (Hoffmann et al., 2014),"
86 " BT(%.4f/cm) - BT(%.4f/cm) [K]\n"
87 "# $22 = SO2 index (Hoffmann et al., 2014) error [K]\n\n",
88 l1.nu_mw[si_low_nu1], l1.nu_mw[si_low_nu2],
89 l1.nu_mw[si_high_nu1], l1.nu_mw[si_high_nu2],
90 l1.nu_mw[si_oper_nu1], l1.nu_mw[si_oper_nu2],
91 l1.nu_mw[si_old_nu1], l1.nu_mw[si_old_nu2]);
92 }
93
94 /* Loop over footprints... */
95 for (int track = 0; track < L1_NTRACK; track++)
96 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
97 for (int ifov = 0; ifov < L1_NFOV; ifov++) {
98
99 /* cloud index... */
100 double ci = BRIGHT(l1.rad_mw[track][xtrack][ifov][ci_nu] * 0.001,
101 l1.nu_mw[ci_nu]);
102 double ci_err =
103 get_noise(l1.rad_mw[track][xtrack][ifov][ci_nu] * 0.001,
104 l1.nedn_mw[ifov][ci_nu] * 0.001,
105 l1.nu_mw[ci_nu]);
106
107 /* ash index (low wavenumbers)... */
108 double ai_low_bt1 =
109 BRIGHT(l1.rad_lw[track][xtrack][ifov][ai_low_nu1] *
110 0.001, l1.nu_lw[ai_low_nu1]);
111 double ai_low_bt2 =
112 BRIGHT(l1.rad_lw[track][xtrack][ifov][ai_low_nu2] *
113 0.001, l1.nu_lw[ai_low_nu2]);
114 double ai_low = ai_low_bt1 - ai_low_bt2;
115 double ai_low_err =
116 sqrt(gsl_pow_2
117 (get_noise
118 (l1.rad_lw[track][xtrack][ifov][ai_low_nu1] * 0.001,
119 l1.nedn_lw[ifov][ai_low_nu1] * 0.001,
120 l1.nu_lw[ai_low_nu1]))
121 +
122 gsl_pow_2(get_noise
123 (l1.rad_lw[track][xtrack][ifov][ai_low_nu2] *
124 0.001, l1.nedn_lw[ifov][ai_low_nu2] * 0.001,
125 l1.nu_lw[ai_low_nu2])));
126
127 /* ash index (high wavenumbers)... */
128 double ai_high_bt1 =
129 BRIGHT(l1.rad_mw[track][xtrack][ifov][ai_high_nu1] *
130 0.001, l1.nu_mw[ai_high_nu1]);
131 double ai_high_bt2 =
132 BRIGHT(l1.rad_lw[track][xtrack][ifov][ai_high_nu2] *
133 0.001, l1.nu_lw[ai_high_nu2]);
134 double ai_high = ai_high_bt1 - ai_high_bt2;
135 double ai_high_err =
136 sqrt(gsl_pow_2
137 (get_noise
138 (l1.rad_mw[track][xtrack][ifov][ai_high_nu1] * 0.001,
139 l1.nedn_mw[ifov][ai_high_nu1] * 0.001,
140 l1.nu_mw[ai_high_nu1]))
141 +
142 gsl_pow_2(get_noise
143 (l1.rad_lw[track][xtrack][ifov][ai_high_nu2] *
144 0.001, l1.nedn_lw[ifov][ai_high_nu2] * 0.001,
145 l1.nu_lw[ai_high_nu1])));
146
147 /* ash index (old)... */
148 double ai_old_bt1 =
149 BRIGHT(l1.rad_lw[track][xtrack][ifov][ai_old_nu1] *
150 0.001, l1.nu_lw[ai_old_nu1]);
151 double ai_old_bt2 =
152 BRIGHT(l1.rad_lw[track][xtrack][ifov][ai_old_nu2] *
153 0.001, l1.nu_lw[ai_old_nu2]);
154 double ai_old = ai_old_bt1 - ai_old_bt2;
155 double ai_old_err =
156 sqrt(gsl_pow_2
157 (get_noise
158 (l1.rad_lw[track][xtrack][ifov][ai_old_nu1] * 0.001,
159 l1.nedn_lw[ifov][ai_old_nu1] * 0.001,
160 l1.nu_lw[ai_old_nu1]))
161 +
162 gsl_pow_2(get_noise
163 (l1.rad_lw[track][xtrack][ifov][ai_old_nu2] *
164 0.001, l1.nedn_lw[ifov][ai_old_nu2] * 0.001,
165 l1.nu_lw[ai_old_nu2])));
166
167 /* SO2 index (low concentrations)... */
168 double si_low_bt1 =
169 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_low_nu1] *
170 0.001, l1.nu_mw[si_low_nu1]);
171 double si_low_bt2 =
172 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_low_nu2] *
173 0.001, l1.nu_mw[si_low_nu2]);
174 double si_low = si_low_bt1 - si_low_bt2;
175 double si_low_err =
176 sqrt(gsl_pow_2
177 (get_noise
178 (l1.rad_mw[track][xtrack][ifov][si_low_nu1] * 0.001,
179 l1.nedn_mw[ifov][si_low_nu1] * 0.001,
180 l1.nu_mw[si_low_nu1]))
181 +
182 gsl_pow_2(get_noise
183 (l1.rad_mw[track][xtrack][ifov][si_low_nu2] *
184 0.001, l1.nedn_mw[ifov][si_low_nu2] * 0.001,
185 l1.nu_mw[si_low_nu2])));
186
187 /* SO2 index (high concentrations)... */
188 double si_high_bt1 =
189 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_high_nu1] *
190 0.001, l1.nu_mw[si_high_nu1]);
191 double si_high_bt2 =
192 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_high_nu2] *
193 0.001, l1.nu_mw[si_high_nu2]);
194 double si_high = si_high_bt1 - si_high_bt2;
195 double si_high_err =
196 sqrt(gsl_pow_2
197 (get_noise
198 (l1.rad_mw[track][xtrack][ifov][si_high_nu1] * 0.001,
199 l1.nedn_mw[ifov][si_high_nu1] * 0.001,
200 l1.nu_mw[si_high_nu1]))
201 +
202 gsl_pow_2(get_noise
203 (l1.rad_mw[track][xtrack][ifov][si_high_nu2] *
204 0.001, l1.nedn_mw[ifov][si_high_nu2] * 0.001,
205 l1.nu_mw[si_high_nu2])));
206
207 /* SO2 index (operational)... */
208 double si_oper_bt1 =
209 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_oper_nu1] *
210 0.001, l1.nu_mw[si_oper_nu1]);
211 double si_oper_bt2 =
212 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_oper_nu2] *
213 0.001, l1.nu_mw[si_oper_nu2]);
214 double si_oper = si_oper_bt1 - si_oper_bt2;
215 double si_oper_err =
216 sqrt(gsl_pow_2
217 (get_noise
218 (l1.rad_mw[track][xtrack][ifov][si_oper_nu1] * 0.001,
219 l1.nedn_mw[ifov][si_oper_nu1] * 0.001,
220 l1.nu_mw[si_oper_nu1]))
221 +
222 gsl_pow_2(get_noise
223 (l1.rad_mw[track][xtrack][ifov][si_oper_nu2] *
224 0.001, l1.nedn_mw[ifov][si_oper_nu2] * 0.001,
225 l1.nu_mw[si_oper_nu2])));
226
227 /* SO2 index (old)... */
228 double si_old_bt1 =
229 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_old_nu1] *
230 0.001, l1.nu_mw[si_old_nu1]);
231 double si_old_bt2 =
232 BRIGHT(l1.rad_mw[track][xtrack][ifov][si_old_nu2] *
233 0.001, l1.nu_mw[si_old_nu2]);
234 double si_old = si_old_bt1 - si_old_bt2;
235 double si_old_err =
236 sqrt(gsl_pow_2
237 (get_noise
238 (l1.rad_mw[track][xtrack][ifov][si_old_nu1] * 0.001,
239 l1.nedn_mw[ifov][si_old_nu1] * 0.001,
240 l1.nu_mw[si_old_nu1]))
241 +
242 gsl_pow_2(get_noise
243 (l1.rad_mw[track][xtrack][ifov][si_old_nu2] *
244 0.001, l1.nedn_mw[ifov][si_old_nu2] * 0.001,
245 l1.nu_mw[si_old_nu2])));
246
247 /* Write output... */
248 fprintf(out,
249 "%.2f %.4f %.4f %.3f %.4f %.4f %.2f %.2f %.2f %.2f %.2f %.2f "
250 "%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",
251 l1.time[track][xtrack] - 220838400,
252 l1.lon[track][xtrack][ifov],
253 l1.lat[track][xtrack][ifov],
254 l1.sat_z[track],
255 l1.sat_lon[track],
256 l1.sat_lat[track],
257 ci, ci_err, ai_low, ai_low_err, ai_high, ai_high_err,
258 ai_old, ai_old_err, si_low, si_low_err, si_high,
259 si_high_err, si_oper, si_oper_err, si_old, si_old_err);
260 }
261 }
262
263 /* Close file... */
264 fclose(out);
265
266 return EXIT_SUCCESS;
267}
268
269/************************************************************************/
270
271double get_noise(
272 double rad,
273 double nesr,
274 double nu) {
275
276 return BRIGHT(rad + nesr, nu) - BRIGHT(rad, nu);
277}
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5278
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:127
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:238
read_cris_l1
int read_cris_l1(char *filename, cris_l1_t *l1, int apo)
Read CrIS Level-1 data.
Definition: libcris.c:1100
L1_NXTRACK
#define L1_NXTRACK
Across-track size of CrIS radiance granule.
Definition: libcris.h:24
L1_NTRACK
#define L1_NTRACK
Along-track size of CrIS radiance granule.
Definition: libcris.h:21
L1_NFOV
#define L1_NFOV
Number of field of views of CrIS radiance granule.
Definition: libcris.h:27
cris_l1_t
CrIS Level-1 data.
Definition: libcris.h:72
cris_l1_t::sat_z
double sat_z[L1_NTRACK]
Satellite altitude [km].
Definition: libcris.h:84
cris_l1_t::rad_mw
float rad_mw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_MW]
Midwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:108
cris_l1_t::nedn_mw
float nedn_mw[L1_NFOV][L1_NCHAN_MW]
Midwave radiance noise [W/(m^2 sr cm^-1)].
Definition: libcris.h:111
cris_l1_t::nu_mw
double nu_mw[L1_NCHAN_MW]
Midwave channel frequencies [cm^-1].
Definition: libcris.h:105
cris_l1_t::nedn_lw
float nedn_lw[L1_NFOV][L1_NCHAN_LW]
Longwave radiance noise [W/(m^2 sr cm^-1)].
Definition: libcris.h:99
cris_l1_t::lon
double lon[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Footprint longitude [deg].
Definition: libcris.h:78
cris_l1_t::rad_lw
float rad_lw[L1_NTRACK][L1_NXTRACK][L1_NFOV][L1_NCHAN_LW]
Longwave radiance [W/(m^2 sr cm^-1)].
Definition: libcris.h:96
cris_l1_t::sat_lat
double sat_lat[L1_NTRACK]
Satellite latitude [deg].
Definition: libcris.h:90
cris_l1_t::time
double time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libcris.h:75
cris_l1_t::sat_lon
double sat_lon[L1_NTRACK]
Satellite longitude [deg].
Definition: libcris.h:87
cris_l1_t::nu_lw
double nu_lw[L1_NCHAN_LW]
Longwave channel frequencies [cm^-1].
Definition: libcris.h:93
cris_l1_t::lat
double lat[L1_NTRACK][L1_NXTRACK][L1_NFOV]
Footprint latitude [deg].
Definition: libcris.h:81
main
int main(int argc, char *argv[])
Definition: volcano.c:17
get_noise
double get_noise(double rad, double nesr, double nu)
Definition: volcano.c:271
Generated on Mon Aug 18 2025 19:49:24 for CrIS Code Collection by doxygen 1.9.4