51 {
52
53 FILE *in;
54
55 static char *line = NULL;
56
58 nu, nu0, nu1, nuf[
MAXNF], press, temp, u;
59
60 static int i, idx, nf, npts;
61
62 size_t line_buf_size = 0;
63
64
66
67
68 if (argc != 5)
69 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
70 "Usage: tblgen <press> <temp> <spec> <filter>\n\n"
71 "Use -h for full help.");
72 sscanf(argv[1], "%lg", &press);
73 sscanf(argv[2], "%lg", &temp);
74
75
76 u = 1e-6 * press * 100 / (1.380658e-23 * temp) * 1000 / 1e4;
77
78
79 if (!(in = fopen(argv[4], "r")))
80 ERRMSG(
"Cannot open filter file!");
81 while (getline(&line, &line_buf_size, in) != -1)
82 if (sscanf(line, "%lg %lg", &nuf[nf], &filt[nf]) == 2)
84 ERRMSG(
"Too many points in filter function");
85 fclose(in);
86
87
88 if (!(in = fopen(argv[3], "r")))
89 ERRMSG(
"Cannot open spectrum!");
90 for (i = 0; i < 4; ++i)
91 if (getline(&line, &line_buf_size, in) == -1)
92 ERRMSG(
"Error while reading spectrum!");
93 sscanf(line, "%d %lg %lg %lg", &npts, &nu0, &dnu, &nu1);
95 ERRMSG(
"Too many points in optical depth spectrum!");
96 i = 0;
97 while (getline(&line, &line_buf_size, in) != -1) {
98 char *tok = strtok(line, " \t\n");
99 while (tok != NULL) {
101 ERRMSG(
"Too many data points in spectrum file!");
102 sscanf(tok, "%lg", &abs[i]);
103 abs[i] /= u;
104 ++i;
105 tok = strtok(NULL, " \t\n");
106 }
107 }
108 fclose(in);
109
110
111 dnu = (nu1 - nu0) / ((double) npts - 1.0);
112 const double r0 = (nuf[0] - nu0) / (nu1 - nu0) * (double) npts;
113 const int i0 = (int) r0;
114
115
116 for (u = 1.0; u <= 1e30; u *= 1.122) {
117
118
119 double epssum = 0, fsum = 0;
120 for (i = i0; i < npts; i++) {
121 nu = nu0 + dnu * (double) i;
122 if (nu < nuf[0])
123 continue;
124 else if (nu > nuf[nf - 1])
125 break;
126 else {
127 if (nu < nuf[idx] || nu > nuf[idx + 1])
129 f =
LIN(nuf[idx], filt[idx], nuf[idx + 1], filt[idx + 1], nu);
130 fsum += f;
131 epssum += f * exp(-abs[i] * u);
132 }
133 }
134 epssum = 1 - epssum / fsum;
135
136
137 if (epssum >= 1e-6 && epssum <= 0.999999 && epssum > epsold)
138 printf("%g %g %g %g\n", press, temp, u, epssum);
139 epsold = epssum;
140
141
142 if (epssum > 0.999999)
143 return EXIT_SUCCESS;
144 }
145
146 return EXIT_SUCCESS;
147}
int locate_irr(const double *xx, const int n, const double x)
Locate index for interpolation on an irregular grid.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define USAGE
Print usage information on -h or --help.
#define LIN(x0, y0, x1, y1, x)
Compute linear interpolation.
1.9.4