| 1 |
/* |
| 2 |
Glazing system multi-layer optics calculation. |
| 3 |
Source equations: LBNL WINDOW technical documention. |
| 4 |
System calcuation: Section 7.6.1 |
| 5 |
Anuglar calculation Section 7.7 |
| 6 |
T. Wang |
| 7 |
*/ |
| 8 |
|
| 9 |
#include "rtmath.h" |
| 10 |
#include "color.h" |
| 11 |
#include "data.h" |
| 12 |
|
| 13 |
enum { |
| 14 |
NTHETA = 59, |
| 15 |
MAX_NAME = 256, |
| 16 |
}; |
| 17 |
|
| 18 |
const double DBL_EPSILON = 1E-9; |
| 19 |
const double FLT_MAX = 3.40282347E+38; |
| 20 |
const double RAD_PER_DEG = 0.0174532925199; |
| 21 |
const double COEFFS_TAU_CLEAR[] = {-0.0015, 3.355, -3.84, 1.46, 0.0288}; |
| 22 |
const double COEFFS_TAU_BRONZE[] = {0.002, 2.813, -2.341, -0.05725, 0.599}; |
| 23 |
const double COEFFS_RHO_CLEAR[] = {0.999, -0.563, 2.043, -2.532, 1.054}; |
| 24 |
const double COEFFS_RHO_BRONZE[] = {0.997, -1.868, 6.513, -7.862, 3.225}; |
| 25 |
|
| 26 |
const double THETAS[NTHETA] = { |
| 27 |
0., 5., 10., 15., 20., 25., 30., 35., 40., |
| 28 |
41., 42., 43., 44., 45., 46., 47., 48., 49, 50., |
| 29 |
51., 52., 53., 54., 55., 56., 57., 58., 59, 60., |
| 30 |
61., 62., 63., 64., 65., 66., 67., 68., 69, 70., |
| 31 |
71., 72., 73., 74., 75., 76., 77., 78., 79, 80., |
| 32 |
81., 82., 83., 84., 85., 86., 87., 88., 89, 90., |
| 33 |
}; |
| 34 |
|
| 35 |
|
| 36 |
char *progname; |
| 37 |
|
| 38 |
typedef struct { |
| 39 |
char *filename; |
| 40 |
int is_mono; |
| 41 |
double thickness_m; |
| 42 |
|
| 43 |
/* Interpolated data at standard wavelengths */ |
| 44 |
double *t_0; |
| 45 |
double *rf_0; |
| 46 |
double *rb_0; |
| 47 |
|
| 48 |
/* Calculated angular data */ |
| 49 |
double *t_lambda_theta; |
| 50 |
double *rf_lambda_theta; |
| 51 |
double *rb_lambda_theta; |
| 52 |
|
| 53 |
} GlazingLayer; |
| 54 |
|
| 55 |
|
| 56 |
inline double polynomial_5(double x, const double coeffs[5]) { |
| 57 |
return x * ( x * ( x * (coeffs[4] * x + coeffs[3]) + coeffs[2]) + coeffs[1]) + coeffs[0]; |
| 58 |
} |
| 59 |
|
| 60 |
|
| 61 |
/* Calculate rho0 (unpolarized reflectance of a single surface) */ |
| 62 |
static void rho0_calc(const double *t_0, const double *r_0, const int nwvl, double *rho0) { |
| 63 |
for (int i = 0; i < nwvl; ++i) { |
| 64 |
const double beta = t_0[i] * t_0[i] - r_0[i] * r_0[i] + 2.0 * r_0[i] + 1.0; |
| 65 |
double denominator = 2.0 * (2.0 - r_0[i]); |
| 66 |
if (fabs(denominator) < DBL_EPSILON) |
| 67 |
denominator = DBL_EPSILON; |
| 68 |
double discriminant = beta * beta - 2.0 * denominator * r_0[i]; |
| 69 |
if (discriminant < 0.0) |
| 70 |
discriminant = 0.0; |
| 71 |
rho0[i] = (beta - sqrt(discriminant)) / denominator; |
| 72 |
if (rho0[i] < 0.0) |
| 73 |
rho0[i] = 0.0; |
| 74 |
if (rho0[i] > 1.0) |
| 75 |
rho0[i] = 1.0; |
| 76 |
} |
| 77 |
} |
| 78 |
|
| 79 |
/* refraction index calc */ |
| 80 |
static void n_calc(const double *rho_0, const int nwvl, double *n) |
| 81 |
{ |
| 82 |
for (int i=0; i < nwvl; ++i) { |
| 83 |
const double sqrt_rho = sqrt(rho_0[i]); |
| 84 |
n[i] = ((1.0 - sqrt_rho) < DBL_EPSILON) ? 1.0 : (1.0 + sqrt_rho) / (1.0 - sqrt_rho); |
| 85 |
} |
| 86 |
} |
| 87 |
|
| 88 |
/* Calculate absoprtion coefficient */ |
| 89 |
static void alpha_calc(const double *t0, const double *r0, const double *rho0, |
| 90 |
const double thickness_m, const int nwvl, double *abs_coeff) |
| 91 |
{ |
| 92 |
for (int i = 0; i < nwvl; ++i) { |
| 93 |
const double numerator = r0[i] - rho0[i]; |
| 94 |
const double denominator = rho0[i] * t0[i]; |
| 95 |
if (denominator > DBL_EPSILON && numerator > DBL_EPSILON) { |
| 96 |
const double log_val = log(numerator / denominator); |
| 97 |
abs_coeff[i] = -log_val / thickness_m; |
| 98 |
} else { |
| 99 |
abs_coeff[i] = 1.0; |
| 100 |
} |
| 101 |
if (abs_coeff[i] < 0.0 || isnan(abs_coeff[i]) || isinf(abs_coeff[i])) { |
| 102 |
abs_coeff[i] = 0.0; |
| 103 |
} |
| 104 |
} |
| 105 |
} |
| 106 |
|
| 107 |
|
| 108 |
void angular_monolithic(GlazingLayer *layer, const int nwvl) { |
| 109 |
double *rho_0 = malloc(nwvl * sizeof(double)); |
| 110 |
double *refrac = malloc(nwvl * sizeof(double)); |
| 111 |
double *abs_coeff = malloc(nwvl * sizeof(double)); |
| 112 |
|
| 113 |
rho0_calc(layer->t_0, layer->rf_0, nwvl, rho_0); |
| 114 |
n_calc(rho_0, nwvl, refrac); |
| 115 |
alpha_calc(layer->t_0, layer->rf_0, rho_0, layer->thickness_m, nwvl, abs_coeff); |
| 116 |
|
| 117 |
for (int itheta = 0; itheta < NTHETA; ++itheta) { |
| 118 |
double phi = THETAS[itheta] * RAD_PER_DEG; |
| 119 |
double cos_phi = cos(phi); |
| 120 |
double sin_phi = sin(phi); |
| 121 |
|
| 122 |
for (int iwvl = 0; iwvl < nwvl; ++iwvl) { |
| 123 |
double sin_phi_prime_sq = sin_phi * sin_phi / (refrac[iwvl] * refrac[iwvl]); |
| 124 |
|
| 125 |
double cos_phi_prime = (sin_phi_prime_sq >= 1.0) ? 0.0 : sqrt(1.0 - sin_phi_prime_sq); |
| 126 |
|
| 127 |
double rho_lambda; |
| 128 |
if (cos_phi_prime <= DBL_EPSILON) { |
| 129 |
rho_lambda = 1.0; |
| 130 |
} else { |
| 131 |
double nf_cos_phi = refrac[iwvl] * cos_phi; |
| 132 |
double nf_cos_phi_p = refrac[iwvl] * cos_phi_prime; |
| 133 |
double par_num = nf_cos_phi - cos_phi_prime; |
| 134 |
double par_den = nf_cos_phi + cos_phi_prime; |
| 135 |
double per_num = nf_cos_phi_p - cos_phi; |
| 136 |
double per_den = nf_cos_phi_p + cos_phi; |
| 137 |
if (fabs(par_den) < DBL_EPSILON || fabs(per_den) < DBL_EPSILON) { |
| 138 |
rho_lambda = 1.0; |
| 139 |
} else { |
| 140 |
double par_comp = pow(par_num / par_den, 2.0); |
| 141 |
double per_comp = pow(per_num / per_den, 2.0); |
| 142 |
rho_lambda = 0.5 * (par_comp + per_comp); |
| 143 |
} |
| 144 |
} |
| 145 |
if (rho_lambda < 0.0) |
| 146 |
rho_lambda = 0.0; |
| 147 |
if (rho_lambda > 1.0) |
| 148 |
rho_lambda = 1.0; |
| 149 |
|
| 150 |
double exp_arg = (cos_phi_prime < DBL_EPSILON) ? -FLT_MAX : (-abs_coeff[iwvl] * layer->thickness_m / cos_phi_prime); |
| 151 |
/* Clamp exponent argument to prevent overflow in exp() */ |
| 152 |
if (exp_arg < -700.0) |
| 153 |
exp_arg = -700.0; |
| 154 |
|
| 155 |
double T_internal = exp(exp_arg); |
| 156 |
|
| 157 |
double denominator = 1.0 - rho_lambda * rho_lambda * T_internal * T_internal; |
| 158 |
|
| 159 |
double t_lambda, r_lambda; |
| 160 |
|
| 161 |
if (fabs(denominator) < DBL_EPSILON || rho_lambda >= 1.0) { |
| 162 |
t_lambda = 0.0; |
| 163 |
r_lambda = 1.0; |
| 164 |
} else { |
| 165 |
double tau_lambda = 1.0 - rho_lambda; |
| 166 |
t_lambda = tau_lambda * tau_lambda * T_internal / denominator; |
| 167 |
r_lambda = rho_lambda * (1.0 + t_lambda * T_internal); |
| 168 |
|
| 169 |
} |
| 170 |
if (t_lambda < 0.0) |
| 171 |
t_lambda = 0.0; |
| 172 |
if (t_lambda > 1.0) |
| 173 |
t_lambda = 1.0; |
| 174 |
|
| 175 |
const int flat_idx = itheta * nwvl + iwvl; |
| 176 |
layer->t_lambda_theta[flat_idx] = t_lambda; |
| 177 |
layer->rf_lambda_theta[flat_idx] = r_lambda; |
| 178 |
layer->rb_lambda_theta[flat_idx] = r_lambda; |
| 179 |
} |
| 180 |
} |
| 181 |
} |
| 182 |
|
| 183 |
|
| 184 |
void angular_coated(GlazingLayer *layer, const int nwvl) { |
| 185 |
for (int itheta = 0; itheta < NTHETA; ++itheta) { |
| 186 |
const double phi = THETAS[itheta] * RAD_PER_DEG; |
| 187 |
const double cos_phi = cos(phi); |
| 188 |
|
| 189 |
for (int iwvl = 0; iwvl < nwvl; ++iwvl) { |
| 190 |
const double t_0_lambda = layer->t_0[iwvl]; |
| 191 |
const double *coeffs_tau = (t_0_lambda > 0.645) ? COEFFS_TAU_CLEAR : COEFFS_TAU_BRONZE; |
| 192 |
const double *coeffs_rho = (t_0_lambda > 0.645) ? COEFFS_RHO_CLEAR : COEFFS_RHO_BRONZE; |
| 193 |
|
| 194 |
const double tau_bar = polynomial_5(cos_phi, coeffs_tau); |
| 195 |
const double rho_bar_term = polynomial_5(cos_phi, coeffs_rho); |
| 196 |
const double rho_bar = rho_bar_term - tau_bar; |
| 197 |
|
| 198 |
const int flat_idx = itheta * nwvl + iwvl; |
| 199 |
layer->t_lambda_theta[flat_idx] = t_0_lambda * tau_bar; |
| 200 |
layer->rf_lambda_theta[flat_idx] = layer->rf_0[iwvl] * (1.0 - rho_bar) + rho_bar; |
| 201 |
layer->rb_lambda_theta[flat_idx] = layer->rb_0[iwvl] * (1.0 - rho_bar) + rho_bar; |
| 202 |
} |
| 203 |
} |
| 204 |
} |
| 205 |
|
| 206 |
|
| 207 |
void multi_layer_calc( |
| 208 |
GlazingLayer *layers, |
| 209 |
int nlayers, |
| 210 |
int nwvl, |
| 211 |
double *total_t, |
| 212 |
double *total_rf, |
| 213 |
double *total_rb) |
| 214 |
{ |
| 215 |
if (nlayers <= 0) |
| 216 |
return; |
| 217 |
|
| 218 |
size_t total_size = (size_t)NTHETA * nwvl * sizeof(double); |
| 219 |
|
| 220 |
memcpy(total_t, layers[0].t_lambda_theta, total_size); |
| 221 |
memcpy(total_rf, layers[0].rf_lambda_theta, total_size); |
| 222 |
memcpy(total_rb, layers[0].rb_lambda_theta, total_size); |
| 223 |
|
| 224 |
if (nlayers == 1) |
| 225 |
return; |
| 226 |
|
| 227 |
double *prev_t = malloc(total_size); |
| 228 |
double *prev_rf = malloc(total_size); |
| 229 |
double *prev_rb = malloc(total_size); |
| 230 |
|
| 231 |
if (!prev_t || !prev_rf || !prev_rb) |
| 232 |
perror("Failed to allocate temporary storage in multi_layer"); |
| 233 |
|
| 234 |
for (int j = 1; j < nlayers; ++j) { |
| 235 |
memcpy(prev_t, total_t, total_size); |
| 236 |
memcpy(prev_rf, total_rf, total_size); |
| 237 |
memcpy(prev_rb, total_rb, total_size); |
| 238 |
|
| 239 |
double *t_j = layers[j].t_lambda_theta; |
| 240 |
double *rf_j = layers[j].rf_lambda_theta; |
| 241 |
double *rb_j = layers[j].rb_lambda_theta; |
| 242 |
|
| 243 |
for (int itheta = 0; itheta < NTHETA; ++itheta) { |
| 244 |
for (int iwvl = 0; iwvl < nwvl; ++iwvl) { |
| 245 |
const int idx = itheta * nwvl + iwvl; |
| 246 |
double denominator = 1.0 - rf_j[idx] * prev_rb[idx]; |
| 247 |
if (fabs(denominator) < DBL_EPSILON) { |
| 248 |
denominator = DBL_EPSILON; |
| 249 |
} |
| 250 |
total_t[idx] = prev_t[idx] * t_j[idx] / denominator; |
| 251 |
total_rf[idx] = prev_rf[idx] + prev_t[idx] * prev_t[idx] * rf_j[idx] / denominator; |
| 252 |
total_rb[idx] = rb_j[idx] + t_j[idx] * t_j[idx] * prev_rb[idx] / denominator; |
| 253 |
} |
| 254 |
} |
| 255 |
} |
| 256 |
free(prev_t); |
| 257 |
free(prev_rf); |
| 258 |
free(prev_rb); |
| 259 |
} |
| 260 |
|
| 261 |
|
| 262 |
int add_layer(GlazingLayer **layers_array, int *count, int *capacity, |
| 263 |
const char *filename, double thickness_m, int is_monolithic) { |
| 264 |
if (*count >= *capacity) { |
| 265 |
*capacity = (*capacity == 0) ? 4 : *capacity * 2; |
| 266 |
GlazingLayer *new_layers = realloc(*layers_array, *capacity * sizeof(GlazingLayer)); |
| 267 |
if (!new_layers) { |
| 268 |
perror("Failed to reallocate memory for layers"); |
| 269 |
return 0; |
| 270 |
} |
| 271 |
*layers_array = new_layers; |
| 272 |
} |
| 273 |
|
| 274 |
memset(&((*layers_array)[*count]), 0, sizeof(GlazingLayer)); |
| 275 |
|
| 276 |
#ifdef _WIN32 |
| 277 |
(*layers_array)[*count].filename = _strdup(filename); |
| 278 |
#else |
| 279 |
(*layers_array)[*count].filename = strdup(filename); |
| 280 |
#endif |
| 281 |
|
| 282 |
if (!(*layers_array)[*count].filename) { |
| 283 |
perror("Failed to duplicate filename"); |
| 284 |
return 0; |
| 285 |
} |
| 286 |
(*layers_array)[*count].is_mono = is_monolithic; |
| 287 |
(*layers_array)[*count].thickness_m = thickness_m; |
| 288 |
(*layers_array)[*count].t_0 = NULL; |
| 289 |
(*layers_array)[*count].rf_0 = NULL; |
| 290 |
(*layers_array)[*count].rb_0 = NULL; |
| 291 |
(*layers_array)[*count].t_lambda_theta = NULL; |
| 292 |
(*layers_array)[*count].rf_lambda_theta = NULL; |
| 293 |
(*layers_array)[*count].rb_lambda_theta = NULL; |
| 294 |
|
| 295 |
(*count)++; |
| 296 |
return 1; |
| 297 |
} |
| 298 |
|
| 299 |
|
| 300 |
int interpolate(GlazingLayer *layer, const double wvl_start, const double wvl_end, const double wvl_interval, const int nwvl) { |
| 301 |
DATARRAY *dp = getdata(layer->filename); |
| 302 |
if (!dp) { |
| 303 |
fprintf(stderr, "Error: Cannot open file '%s'\n", layer->filename); |
| 304 |
return 0; |
| 305 |
} |
| 306 |
|
| 307 |
layer->t_0 = malloc(nwvl * sizeof(double)); |
| 308 |
layer->rf_0 = malloc(nwvl * sizeof(double)); |
| 309 |
layer->rb_0 = malloc(nwvl * sizeof(double)); |
| 310 |
|
| 311 |
double wvl = wvl_start; |
| 312 |
int i = 0; |
| 313 |
while (wvl <= wvl_end) { |
| 314 |
double t_pt[2] = {2., wvl}; |
| 315 |
layer->t_0[i] = datavalue(dp, t_pt); |
| 316 |
double rf_pt[2] = {0., wvl}; |
| 317 |
layer->rf_0[i] = datavalue(dp, rf_pt); |
| 318 |
double rb_pt[2] = {1., wvl}; |
| 319 |
layer->rb_0[i] = datavalue(dp, rb_pt); |
| 320 |
wvl = wvl + wvl_interval; |
| 321 |
i = i + 1; |
| 322 |
|
| 323 |
} |
| 324 |
freedata(dp); |
| 325 |
return 1; |
| 326 |
} |
| 327 |
|
| 328 |
|
| 329 |
int write_output_file( |
| 330 |
const char *tfilename, |
| 331 |
const char *rfilename, |
| 332 |
const double *tdata, |
| 333 |
const double *rfdata, |
| 334 |
const double *rbdata, |
| 335 |
const double wvl_start, |
| 336 |
const double wvl_end, |
| 337 |
const int nwvl, |
| 338 |
int argc, char *argv[]) |
| 339 |
{ |
| 340 |
FILE *tfp = fopen(tfilename, "w"); |
| 341 |
FILE *rfp = fopen(rfilename, "w"); |
| 342 |
if (!tfp || !rfp) { |
| 343 |
fprintf(stderr, "Error: Cannot open output files\n"); |
| 344 |
return(0); |
| 345 |
} |
| 346 |
|
| 347 |
fprintf(tfp, "# "); |
| 348 |
for (int i = 0; i< argc; ++i) { |
| 349 |
fprintf(tfp, "%s ", argv[i]); |
| 350 |
} |
| 351 |
fprintf(tfp, "\n2\n0 0 %d", NTHETA); |
| 352 |
for (int i = 0; i < NTHETA; i++) { |
| 353 |
fprintf(tfp, " %d", (int)THETAS[i]); |
| 354 |
} |
| 355 |
fprintf(tfp, "\n%.0f %.0f %d\n", wvl_start, wvl_end, nwvl); |
| 356 |
|
| 357 |
for (int itheta = 0; itheta < NTHETA; ++itheta) { |
| 358 |
for (int iwvl = 0; iwvl < nwvl; ++iwvl) { |
| 359 |
fprintf(tfp, "%.6f\n", tdata[itheta*nwvl+iwvl]); |
| 360 |
} |
| 361 |
} |
| 362 |
fclose(tfp); |
| 363 |
|
| 364 |
fprintf(rfp, "# "); |
| 365 |
for (int i = 0; i< argc; ++i) { |
| 366 |
fprintf(rfp, "%s ", argv[i]); |
| 367 |
} |
| 368 |
fprintf(rfp, "\n2\n0 0 %d", NTHETA * 2 - 1); |
| 369 |
for (int i = 0; i < NTHETA; i++) { |
| 370 |
fprintf(rfp, " %d", (int)THETAS[i]); |
| 371 |
} |
| 372 |
for (int i = NTHETA-2; i >= 0; --i) { |
| 373 |
fprintf(rfp, " %d", (int)(180 - THETAS[i])); |
| 374 |
} |
| 375 |
fprintf(rfp, "\n%.0f %.0f %d\n", wvl_start, wvl_end, nwvl); |
| 376 |
|
| 377 |
for (int itheta = 0; itheta < NTHETA; ++itheta) { |
| 378 |
for (int iwvl = 0; iwvl < nwvl; ++iwvl) { |
| 379 |
fprintf(rfp, "%.6f\n", rfdata[itheta * nwvl + iwvl]); |
| 380 |
} |
| 381 |
} |
| 382 |
for (int itheta = NTHETA - 2; itheta >= 0 ; --itheta) { |
| 383 |
for (int iwvl = 0; iwvl < nwvl; ++iwvl) { |
| 384 |
fprintf(rfp, "%.6f\n", rbdata[itheta * nwvl + iwvl]); |
| 385 |
} |
| 386 |
} |
| 387 |
|
| 388 |
fclose(rfp); |
| 389 |
return(1); |
| 390 |
} |
| 391 |
|
| 392 |
|
| 393 |
void print_usage() { |
| 394 |
fprintf(stderr, "Usage: %s [-m monolithic_layer.dat thickness | -c coated_layer.dat] ...\n", progname); |
| 395 |
fprintf(stderr, "Calculate multi-layer glazing optics from spectral data files.\n"); |
| 396 |
fprintf(stderr, "Options:\n"); |
| 397 |
fprintf(stderr, " -m <filename> thickness Specify an uncoated (monolithic) glazing layer .dat file and its thickness (meter).\n"); |
| 398 |
fprintf(stderr, " -c <filename> Specify a coated or laminate glazing layer .dat file.\n"); |
| 399 |
fprintf(stderr, " -s start_wvl end_wvl interval Specify wavelength range and interval.\n"); |
| 400 |
fprintf(stderr, " -p prefix Specify prefix name to the output files.\n"); |
| 401 |
fprintf(stderr, " -h, --help Show this help message.\n"); |
| 402 |
fprintf(stderr, "Layer order is determined by the sequence of options on the command line.\n"); |
| 403 |
fprintf(stderr, "Output Files: prefix_t.dat, prefix_r.dat\n"); |
| 404 |
} |
| 405 |
|
| 406 |
void cleanup_layers(GlazingLayer *layers, int num_layers) { |
| 407 |
if (!layers) return; |
| 408 |
for (int i = 0; i < num_layers; ++i) { |
| 409 |
if (layers[i].filename) |
| 410 |
free(layers[i].filename); |
| 411 |
} |
| 412 |
free(layers); |
| 413 |
} |
| 414 |
|
| 415 |
|
| 416 |
int main(int argc, char *argv[]) |
| 417 |
{ |
| 418 |
GlazingLayer *layers = NULL; |
| 419 |
int num_layers = 0; |
| 420 |
int layer_capacity = 0; |
| 421 |
int success = 1; |
| 422 |
double wvl_start_nm = 380.; |
| 423 |
double wvl_end_nm = 780.; |
| 424 |
double wvl_interval_nm = 5.; |
| 425 |
int nwvl = 81; |
| 426 |
double thickness_m = 0.003; |
| 427 |
char *filename; |
| 428 |
char *prefix = "unnamed"; |
| 429 |
char file_t[MAX_NAME]; |
| 430 |
char file_r[MAX_NAME]; |
| 431 |
progname = argv[0]; |
| 432 |
if (argc <= 1) { |
| 433 |
print_usage(); |
| 434 |
return EXIT_FAILURE; |
| 435 |
} |
| 436 |
|
| 437 |
for (int i=1; i < argc; ++i) { |
| 438 |
switch (argv[i][1]) { |
| 439 |
case 'm': |
| 440 |
filename = argv[++i]; |
| 441 |
thickness_m = atof(argv[++i]); |
| 442 |
if (!add_layer(&layers, &num_layers, &layer_capacity, filename, thickness_m, 1)) { |
| 443 |
cleanup_layers(layers, num_layers); |
| 444 |
return EXIT_FAILURE; |
| 445 |
} |
| 446 |
break; |
| 447 |
case 'c': |
| 448 |
filename = argv[++i]; |
| 449 |
if (!add_layer(&layers, &num_layers, &layer_capacity, filename, thickness_m, 0)) { |
| 450 |
cleanup_layers(layers, num_layers); |
| 451 |
return EXIT_FAILURE; |
| 452 |
} |
| 453 |
break; |
| 454 |
case 'p': |
| 455 |
prefix = argv[++i]; |
| 456 |
break; |
| 457 |
case 's': |
| 458 |
wvl_start_nm = atof(argv[++i]); |
| 459 |
wvl_end_nm = atof(argv[++i]); |
| 460 |
wvl_interval_nm = atof(argv[++i]); |
| 461 |
if (wvl_start_nm > wvl_end_nm) { |
| 462 |
fprintf(stderr, "Starting wavelength > End wavelength\n"); |
| 463 |
return EXIT_FAILURE; |
| 464 |
} |
| 465 |
if (((int)(wvl_end_nm - wvl_start_nm) % (int)wvl_interval_nm) > 0) { |
| 466 |
fprintf(stderr, |
| 467 |
"Error: Wavelength range (%f to %f nm) must be evenly divisible by the interval (%f nm).\n", |
| 468 |
wvl_start_nm, wvl_end_nm, wvl_interval_nm); |
| 469 |
return EXIT_FAILURE; |
| 470 |
} |
| 471 |
break; |
| 472 |
case 'h': |
| 473 |
print_usage(); |
| 474 |
return 1; |
| 475 |
default: |
| 476 |
break; |
| 477 |
} |
| 478 |
} |
| 479 |
snprintf(file_t, MAX_NAME, "%s_t.dat", prefix); |
| 480 |
snprintf(file_r, MAX_NAME, "%s_r.dat", prefix); |
| 481 |
|
| 482 |
if (fabs(wvl_start_nm - wvl_end_nm) < DBL_EPSILON && wvl_interval_nm > DBL_EPSILON) { |
| 483 |
nwvl = 1; |
| 484 |
} else { |
| 485 |
nwvl = (int)((wvl_end_nm - wvl_start_nm) / wvl_interval_nm + 1.5); |
| 486 |
} |
| 487 |
|
| 488 |
if (num_layers == 0) { |
| 489 |
fprintf(stderr, "Error: No layers specified.\n"); |
| 490 |
print_usage(); |
| 491 |
return EXIT_FAILURE; |
| 492 |
} |
| 493 |
|
| 494 |
for (int i = 0; i < num_layers; ++i) { |
| 495 |
if (!interpolate(&layers[i], wvl_start_nm, wvl_end_nm, wvl_interval_nm, nwvl)) { |
| 496 |
fprintf(stderr, "Error: Failed to parse or interpolate data for layer %d.\n", i + 1); |
| 497 |
success = 0; |
| 498 |
break; |
| 499 |
} |
| 500 |
layers[i].t_lambda_theta = malloc(NTHETA * nwvl * sizeof(double)); |
| 501 |
layers[i].rf_lambda_theta = malloc(NTHETA * nwvl * sizeof(double)); |
| 502 |
layers[i].rb_lambda_theta = malloc(NTHETA * nwvl * sizeof(double)); |
| 503 |
|
| 504 |
if (layers[i].is_mono) { |
| 505 |
angular_monolithic(&layers[i], nwvl); |
| 506 |
} else { |
| 507 |
angular_coated(&layers[i], nwvl); |
| 508 |
} |
| 509 |
} |
| 510 |
|
| 511 |
if (!success) { |
| 512 |
cleanup_layers(layers, num_layers); |
| 513 |
return EXIT_FAILURE; |
| 514 |
} |
| 515 |
|
| 516 |
|
| 517 |
size_t total_flat_size = (size_t)NTHETA * nwvl * sizeof(double); |
| 518 |
double *total_t = malloc(total_flat_size); |
| 519 |
double *total_rf = malloc(total_flat_size); |
| 520 |
double *total_rb = malloc(total_flat_size); |
| 521 |
|
| 522 |
if (!total_t || !total_rf || !total_rb) { |
| 523 |
perror("Failed to allocate memory for final results"); |
| 524 |
cleanup_layers(layers, num_layers); |
| 525 |
free(total_t); |
| 526 |
free(total_rf); |
| 527 |
free(total_rb); |
| 528 |
return EXIT_FAILURE; |
| 529 |
} |
| 530 |
|
| 531 |
multi_layer_calc(layers, num_layers, nwvl, total_t, total_rf, total_rb); |
| 532 |
|
| 533 |
success &= write_output_file(file_t, file_r, total_t, total_rf, total_rb, wvl_start_nm, wvl_end_nm, nwvl, argc, argv); |
| 534 |
|
| 535 |
if (success) { |
| 536 |
printf("# "); |
| 537 |
for (int i = 0;i < argc; ++i) { |
| 538 |
printf("%s ", argv[i]); |
| 539 |
} |
| 540 |
printf("\nvoid specdata refl_spec\n"); |
| 541 |
printf("4 noop %s . 'Acos(Rdot)/DEGREE'\n0\n0\n\n", file_r); |
| 542 |
printf("void specdata trans_spec\n"); |
| 543 |
printf("4 noop %s . 'Acos(abs(Rdot))/DEGREE'\n0\n0\n\n", file_t); |
| 544 |
printf("void WGMDfunc glaze_mat\n13\n\trefl_spec 1 0 0\n\ttrans_spec 1 0 0\n\tvoid\n\t0 0 1 .\n0\n9 0 0 0 0 0 0 0 0 0\n"); |
| 545 |
} else { |
| 546 |
fprintf(stderr, "Error: Failed to write one or more output files.\n"); |
| 547 |
} |
| 548 |
|
| 549 |
cleanup_layers(layers, num_layers); |
| 550 |
free(total_t); |
| 551 |
free(total_rf); |
| 552 |
free(total_rb); |
| 553 |
|
| 554 |
return success ? EXIT_SUCCESS : EXIT_FAILURE; |
| 555 |
} |
