src/gen/genglaze.c

/*
        Glazing system multi-layer optics calculation.
        Source equations: LBNL WINDOW technical documention.
        System calcuation: Section 7.6.1
        Anuglar calculation Section 7.7
        T. Wang
*/

#include "rtmath.h"
#include "color.h"
#include "data.h"

enum {
        NTHETA = 59,
        MAX_NAME = 256,
};

const double DBL_EPSILON = 1E-9;
const double FLT_MAX = 3.40282347E+38;
const double RAD_PER_DEG = 0.0174532925199;
const double COEFFS_TAU_CLEAR[] = {-0.0015, 3.355, -3.84, 1.46, 0.0288};
const double COEFFS_TAU_BRONZE[] = {0.002, 2.813, -2.341, -0.05725, 0.599};
const double COEFFS_RHO_CLEAR[] = {0.999, -0.563, 2.043, -2.532, 1.054};
const double COEFFS_RHO_BRONZE[] = {0.997, -1.868, 6.513, -7.862, 3.225};

const double THETAS[NTHETA] = {
        0., 5., 10., 15., 20., 25., 30., 35., 40.,
        41., 42., 43., 44., 45., 46., 47., 48., 49, 50.,
        51., 52., 53., 54., 55., 56., 57., 58., 59, 60.,
        61., 62., 63., 64., 65., 66., 67., 68., 69, 70.,
        71., 72., 73., 74., 75., 76., 77., 78., 79, 80.,
        81., 82., 83., 84., 85., 86., 87., 88., 89, 90.,
};


char *progname;

typedef struct {
        char *filename;
        int is_mono;
        double thickness_m;

        /* Interpolated data at standard wavelengths */
        double *t_0;
        double *rf_0;
        double *rb_0;

        /* Calculated angular data */
        double *t_lambda_theta;
        double *rf_lambda_theta;
        double *rb_lambda_theta;

} GlazingLayer;


inline double polynomial_5(double x, const double coeffs[5]) {
        return x * ( x * ( x * (coeffs[4] * x + coeffs[3]) + coeffs[2]) + coeffs[1]) + coeffs[0];
}


/* Calculate rho0 (unpolarized reflectance of a single surface) */
static void rho0_calc(const double *t_0, const double *r_0, const int nwvl, double *rho0) {
        for (int i = 0; i < nwvl; ++i) {
                const double beta = t_0[i] * t_0[i] - r_0[i] * r_0[i] + 2.0 * r_0[i] + 1.0;
                double denominator = 2.0 * (2.0 - r_0[i]);
                if (fabs(denominator) < DBL_EPSILON) 
                        denominator = DBL_EPSILON;
                double discriminant = beta * beta - 2.0 * denominator * r_0[i];
                if (discriminant < 0.0) 
                        discriminant = 0.0;
                rho0[i] = (beta - sqrt(discriminant)) / denominator;
                if (rho0[i] < 0.0) 
                        rho0[i] = 0.0;
                if (rho0[i] > 1.0) 
                        rho0[i] = 1.0;
        }
}

/* refraction index calc */
static void n_calc(const double *rho_0, const int nwvl, double *n)
{
        for (int i=0; i < nwvl; ++i) {
                const double sqrt_rho = sqrt(rho_0[i]);
                n[i] = ((1.0 - sqrt_rho) < DBL_EPSILON) ? 1.0 : (1.0 + sqrt_rho) / (1.0 - sqrt_rho);
        }
}

/* Calculate absoprtion coefficient */
static void alpha_calc(const double *t0, const double *r0, const double *rho0,
                                        const double thickness_m, const int nwvl, double *abs_coeff)
{
        for (int i = 0; i < nwvl; ++i) {
                const double numerator = r0[i] - rho0[i];
                const double denominator = rho0[i] * t0[i];
                if (denominator > DBL_EPSILON && numerator > DBL_EPSILON) {
                        const double log_val = log(numerator / denominator);
                        abs_coeff[i] = -log_val / thickness_m;
                } else {
                        abs_coeff[i] = 1.0;
                }
                if (abs_coeff[i] < 0.0 || isnan(abs_coeff[i]) || isinf(abs_coeff[i])) {
                        abs_coeff[i] = 0.0;
                }
        }
}


void angular_monolithic(GlazingLayer *layer, const int nwvl) {
        double *rho_0 = malloc(nwvl * sizeof(double));
        double *refrac = malloc(nwvl * sizeof(double));
        double *abs_coeff = malloc(nwvl * sizeof(double));

        rho0_calc(layer->t_0, layer->rf_0, nwvl, rho_0);
        n_calc(rho_0, nwvl, refrac);
        alpha_calc(layer->t_0, layer->rf_0, rho_0, layer->thickness_m, nwvl, abs_coeff);

        for (int itheta = 0; itheta < NTHETA; ++itheta) {
                double phi = THETAS[itheta] * RAD_PER_DEG;
                double cos_phi = cos(phi);
                double sin_phi = sin(phi);

                for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
                        double sin_phi_prime_sq = sin_phi * sin_phi / (refrac[iwvl] * refrac[iwvl]);

                        double cos_phi_prime = (sin_phi_prime_sq >= 1.0) ? 0.0 : sqrt(1.0 - sin_phi_prime_sq);

                        double rho_lambda;
                        if (cos_phi_prime <= DBL_EPSILON) {
                                rho_lambda = 1.0;
                        } else {
                                double nf_cos_phi = refrac[iwvl] * cos_phi;
                                double nf_cos_phi_p = refrac[iwvl] * cos_phi_prime;
                                double par_num = nf_cos_phi - cos_phi_prime;
                                double par_den = nf_cos_phi + cos_phi_prime;
                                double per_num = nf_cos_phi_p - cos_phi;
                                double per_den = nf_cos_phi_p + cos_phi;
                                if (fabs(par_den) < DBL_EPSILON || fabs(per_den) < DBL_EPSILON) {
                                        rho_lambda = 1.0;
                                } else {
                                        double par_comp = pow(par_num / par_den, 2.0);
                                        double per_comp = pow(per_num / per_den, 2.0);
                                        rho_lambda = 0.5 * (par_comp + per_comp);
                                }
                        }
                        if (rho_lambda < 0.0)
                                rho_lambda = 0.0;
                        if (rho_lambda > 1.0)
                                rho_lambda = 1.0;

                        double exp_arg = (cos_phi_prime < DBL_EPSILON) ? -FLT_MAX : (-abs_coeff[iwvl] * layer->thickness_m / cos_phi_prime);
                        /* Clamp exponent argument to prevent overflow in exp() */
                        if (exp_arg < -700.0) 
                                exp_arg = -700.0;

                        double T_internal = exp(exp_arg);

                        double denominator = 1.0 - rho_lambda * rho_lambda * T_internal * T_internal;

                        double t_lambda, r_lambda;

                        if (fabs(denominator) < DBL_EPSILON || rho_lambda >= 1.0) {
                                t_lambda = 0.0;
                                r_lambda = 1.0;
                        } else {
                                double tau_lambda = 1.0 - rho_lambda;
                                t_lambda = tau_lambda * tau_lambda * T_internal / denominator;
                                r_lambda = rho_lambda * (1.0 + t_lambda * T_internal);

                        }
                        if (t_lambda < 0.0)
                                t_lambda = 0.0;
                        if (t_lambda > 1.0)
                                t_lambda = 1.0;

                        const int flat_idx = itheta * nwvl + iwvl;
                        layer->t_lambda_theta[flat_idx] = t_lambda;
                        layer->rf_lambda_theta[flat_idx] = r_lambda;
                        layer->rb_lambda_theta[flat_idx] = r_lambda;
                }
        }
}


void angular_coated(GlazingLayer *layer, const int nwvl) {
        for (int itheta = 0; itheta < NTHETA; ++itheta) {
                const double phi = THETAS[itheta] * RAD_PER_DEG;
                const double cos_phi = cos(phi);

                for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
                        const double t_0_lambda = layer->t_0[iwvl];
                        const double *coeffs_tau = (t_0_lambda > 0.645) ? COEFFS_TAU_CLEAR : COEFFS_TAU_BRONZE;
                        const double *coeffs_rho = (t_0_lambda > 0.645) ? COEFFS_RHO_CLEAR : COEFFS_RHO_BRONZE;

                        const double tau_bar = polynomial_5(cos_phi, coeffs_tau);
                        const double rho_bar_term = polynomial_5(cos_phi, coeffs_rho);
                        const double rho_bar = rho_bar_term - tau_bar;

                        const int flat_idx = itheta * nwvl + iwvl;
                        layer->t_lambda_theta[flat_idx] = t_0_lambda * tau_bar;
                        layer->rf_lambda_theta[flat_idx] = layer->rf_0[iwvl] * (1.0 - rho_bar) + rho_bar;
                        layer->rb_lambda_theta[flat_idx] = layer->rb_0[iwvl] * (1.0 - rho_bar) + rho_bar;
                }
        }
}


void multi_layer_calc(
        GlazingLayer *layers,
        int nlayers,
        int nwvl,
        double *total_t,
        double *total_rf,
        double *total_rb)
{
        if (nlayers <= 0) 
                return;

        size_t total_size = (size_t)NTHETA * nwvl * sizeof(double);

        memcpy(total_t, layers[0].t_lambda_theta, total_size);
        memcpy(total_rf, layers[0].rf_lambda_theta, total_size);
        memcpy(total_rb, layers[0].rb_lambda_theta, total_size);

        if (nlayers == 1) 
                return;

        double *prev_t = malloc(total_size);
        double *prev_rf = malloc(total_size);
        double *prev_rb = malloc(total_size);

        if (!prev_t || !prev_rf || !prev_rb)
                perror("Failed to allocate temporary storage in multi_layer");

        for (int j = 1; j < nlayers; ++j) {
                memcpy(prev_t, total_t, total_size);
                memcpy(prev_rf, total_rf, total_size);
                memcpy(prev_rb, total_rb, total_size);

                double *t_j = layers[j].t_lambda_theta;
                double *rf_j = layers[j].rf_lambda_theta;
                double *rb_j = layers[j].rb_lambda_theta;

                for (int itheta = 0; itheta < NTHETA; ++itheta) {
                        for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
                                const int idx = itheta * nwvl + iwvl;
                                double denominator = 1.0 - rf_j[idx] * prev_rb[idx];
                                if (fabs(denominator) < DBL_EPSILON) {
                                        denominator = DBL_EPSILON;
                                }
                                total_t[idx] = prev_t[idx] * t_j[idx] / denominator;
                                total_rf[idx] = prev_rf[idx] + prev_t[idx] * prev_t[idx] * rf_j[idx] / denominator;
                                total_rb[idx] = rb_j[idx] + t_j[idx] * t_j[idx] * prev_rb[idx] / denominator;
                        }
                }
        }
        free(prev_t);
        free(prev_rf);
        free(prev_rb);
}


int add_layer(GlazingLayer **layers_array, int *count, int *capacity,
                          const char *filename, double thickness_m, int is_monolithic) {
        if (*count >= *capacity) {
                *capacity = (*capacity == 0) ? 4 : *capacity * 2;
                GlazingLayer *new_layers = realloc(*layers_array, *capacity * sizeof(GlazingLayer));
                if (!new_layers) {
                        perror("Failed to reallocate memory for layers");
                        return 0;
                }
                *layers_array = new_layers;
        }

        memset(&((*layers_array)[*count]), 0, sizeof(GlazingLayer));

        #ifdef _WIN32
                (*layers_array)[*count].filename = _strdup(filename);
        #else
                (*layers_array)[*count].filename = strdup(filename);
        #endif

        if (!(*layers_array)[*count].filename) {
                 perror("Failed to duplicate filename");
                 return 0;
        }
        (*layers_array)[*count].is_mono = is_monolithic;
        (*layers_array)[*count].thickness_m = thickness_m;
        (*layers_array)[*count].t_0 = NULL;
        (*layers_array)[*count].rf_0 = NULL;
        (*layers_array)[*count].rb_0 = NULL;
        (*layers_array)[*count].t_lambda_theta = NULL;
        (*layers_array)[*count].rf_lambda_theta = NULL;
        (*layers_array)[*count].rb_lambda_theta = NULL;

        (*count)++;
        return 1;
}


int interpolate(GlazingLayer *layer, const double wvl_start, const double wvl_end, const double wvl_interval, const int nwvl) {
        DATARRAY *dp = getdata(layer->filename);
        if (!dp) {
                fprintf(stderr, "Error: Cannot open file '%s'\n", layer->filename);
                return 0;
        }

        layer->t_0 = malloc(nwvl * sizeof(double));
        layer->rf_0 = malloc(nwvl * sizeof(double));
        layer->rb_0 = malloc(nwvl * sizeof(double));

        double wvl = wvl_start;
        int i = 0;
        while (wvl <= wvl_end) {
                double t_pt[2] = {2., wvl};
                layer->t_0[i] = datavalue(dp, t_pt);
                double rf_pt[2] = {0., wvl};
                layer->rf_0[i] = datavalue(dp, rf_pt);
                double rb_pt[2] = {1., wvl};
                layer->rb_0[i] = datavalue(dp, rb_pt);
                wvl = wvl + wvl_interval;
                i = i + 1;

        }
        freedata(dp);
        return 1;
}


int write_output_file(
                const char *tfilename, 
                const char *rfilename, 
                const double *tdata, 
                const double *rfdata, 
                const double *rbdata, 
                const double wvl_start, 
                const double wvl_end, 
                const int nwvl,
                int argc, char *argv[])
{
        FILE *tfp = fopen(tfilename, "w");
        FILE *rfp = fopen(rfilename, "w");
        if (!tfp || !rfp) {
                fprintf(stderr, "Error: Cannot open output files\n");
                return(0);
        }

        fprintf(tfp, "# ");
        for (int i = 0; i< argc; ++i) {
                fprintf(tfp, "%s ", argv[i]);
        }
        fprintf(tfp, "\n2\n0 0 %d", NTHETA);
        for (int i = 0; i < NTHETA; i++) {
                fprintf(tfp, " %d", (int)THETAS[i]);
        }
        fprintf(tfp, "\n%.0f %.0f %d\n", wvl_start, wvl_end, nwvl);

        for (int itheta = 0; itheta < NTHETA; ++itheta) {
                for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
                        fprintf(tfp, "%.6f\n", tdata[itheta*nwvl+iwvl]);
                }
        }
        fclose(tfp);

        fprintf(rfp, "# ");
        for (int i = 0; i< argc; ++i) {
                fprintf(rfp, "%s ", argv[i]);
        }
        fprintf(rfp, "\n2\n0 0 %d", NTHETA * 2 - 1);
        for (int i = 0; i < NTHETA; i++) {
                fprintf(rfp, " %d", (int)THETAS[i]);
        }
        for (int i = NTHETA-2; i >= 0; --i) {
                fprintf(rfp, " %d", (int)(180 - THETAS[i]));
        }
        fprintf(rfp, "\n%.0f %.0f %d\n", wvl_start, wvl_end, nwvl);

        for (int itheta = 0; itheta < NTHETA; ++itheta) {
                for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
                        fprintf(rfp, "%.6f\n", rfdata[itheta * nwvl + iwvl]);
                }
        }
        for (int itheta = NTHETA - 2; itheta >= 0 ; --itheta) {
                for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
                        fprintf(rfp, "%.6f\n", rbdata[itheta * nwvl + iwvl]);
                }
        }

        fclose(rfp);
        return(1);
}


void print_usage() {
        fprintf(stderr, "Usage: %s [-m monolithic_layer.dat thickness | -c coated_layer.dat] ...\n", progname);
        fprintf(stderr, "Calculate multi-layer glazing optics from spectral data files.\n");
        fprintf(stderr, "Options:\n");
        fprintf(stderr, "  -m <filename> thickness    Specify an uncoated (monolithic) glazing layer .dat file and its thickness (meter).\n");
        fprintf(stderr, "  -c <filename>    Specify a coated or laminate glazing layer .dat file.\n");
        fprintf(stderr, "  -s start_wvl end_wvl interval    Specify wavelength range and interval.\n");
        fprintf(stderr, "  -p prefix Specify prefix name to the output files.\n");
        fprintf(stderr, "  -h, --help        Show this help message.\n");
        fprintf(stderr, "Layer order is determined by the sequence of options on the command line.\n");
        fprintf(stderr, "Output Files: prefix_t.dat, prefix_r.dat\n");
}

void cleanup_layers(GlazingLayer *layers, int num_layers) {
         if (!layers) return;
         for (int i = 0; i < num_layers; ++i) {
                if (layers[i].filename)
                        free(layers[i].filename);
        }
        free(layers);
}


int main(int argc, char *argv[])
{
        GlazingLayer *layers = NULL;
        int num_layers = 0;
        int layer_capacity = 0;
        int success = 1;
        double wvl_start_nm = 380.;
        double wvl_end_nm = 780.;
        double wvl_interval_nm = 5.;
        int nwvl = 81;
        double thickness_m = 0.003;
        char *filename;
        char *prefix = "unnamed";
        char file_t[MAX_NAME];
        char file_r[MAX_NAME];
        progname = argv[0];
        if (argc <= 1) {
          print_usage();
          return EXIT_FAILURE;
        }

        for (int i=1; i < argc; ++i) {
                switch (argv[i][1]) {
                case 'm':
                        filename = argv[++i];
                        thickness_m = atof(argv[++i]);
                        if (!add_layer(&layers, &num_layers, &layer_capacity, filename, thickness_m, 1)) {
                                 cleanup_layers(layers, num_layers);
                                 return EXIT_FAILURE;
                        }
                        break;
                case 'c':
                        filename = argv[++i];
                        if (!add_layer(&layers, &num_layers, &layer_capacity, filename, thickness_m, 0)) {
                                 cleanup_layers(layers, num_layers);
                                 return EXIT_FAILURE;
                        }
                        break;
                case 'p':
                        prefix = argv[++i];
                        break;
                case 's':
                        wvl_start_nm = atof(argv[++i]);
                        wvl_end_nm = atof(argv[++i]);
                        wvl_interval_nm = atof(argv[++i]);
                        if (wvl_start_nm > wvl_end_nm) {
                                fprintf(stderr, "Starting wavelength > End wavelength\n");
                                return EXIT_FAILURE;
                        }
                        if (((int)(wvl_end_nm - wvl_start_nm) % (int)wvl_interval_nm) > 0) {
                                fprintf(stderr, 
                                                "Error: Wavelength range (%f to %f nm) must be evenly divisible by the interval (%f nm).\n", 
                                                wvl_start_nm, wvl_end_nm, wvl_interval_nm);
                                return EXIT_FAILURE;
                        }
                        break;
                case 'h':
                        print_usage();
                        return 1;
                default:
                        break;
                }
        }
        snprintf(file_t, MAX_NAME, "%s_t.dat", prefix);
        snprintf(file_r, MAX_NAME, "%s_r.dat", prefix);

        if (fabs(wvl_start_nm - wvl_end_nm) < DBL_EPSILON && wvl_interval_nm > DBL_EPSILON) {
                 nwvl = 1;
        } else {
                 nwvl = (int)((wvl_end_nm - wvl_start_nm) / wvl_interval_nm + 1.5);
        }

        if (num_layers == 0) {
                fprintf(stderr, "Error: No layers specified.\n");
                print_usage();
                return EXIT_FAILURE;
        }

        for (int i = 0; i < num_layers; ++i) {
                if (!interpolate(&layers[i], wvl_start_nm, wvl_end_nm, wvl_interval_nm, nwvl)) {
                        fprintf(stderr, "Error: Failed to parse or interpolate data for layer %d.\n", i + 1);
                        success = 0;
                        break;
                }
                layers[i].t_lambda_theta = malloc(NTHETA * nwvl * sizeof(double));
                layers[i].rf_lambda_theta = malloc(NTHETA * nwvl * sizeof(double));
                layers[i].rb_lambda_theta = malloc(NTHETA * nwvl * sizeof(double));

                if (layers[i].is_mono) {
                        angular_monolithic(&layers[i], nwvl);
                } else {
                        angular_coated(&layers[i], nwvl);
                }
        }

        if (!success) {
                cleanup_layers(layers, num_layers);
                return EXIT_FAILURE;
        }


        size_t total_flat_size = (size_t)NTHETA * nwvl * sizeof(double);
        double *total_t = malloc(total_flat_size);
        double *total_rf = malloc(total_flat_size);
        double *total_rb = malloc(total_flat_size);

        if (!total_t || !total_rf || !total_rb) {
                perror("Failed to allocate memory for final results");
                cleanup_layers(layers, num_layers);
                free(total_t);
                free(total_rf);
                free(total_rb);
                return EXIT_FAILURE;
        }

        multi_layer_calc(layers, num_layers, nwvl, total_t, total_rf, total_rb);

        success &= write_output_file(file_t, file_r, total_t, total_rf, total_rb, wvl_start_nm, wvl_end_nm, nwvl, argc, argv);

        if (success) {
                printf("# ");
                for (int i = 0;i < argc; ++i) {
                        printf("%s ", argv[i]);
                }
                printf("\nvoid specdata refl_spec\n");
                printf("4 noop %s . 'Acos(Rdot)/DEGREE'\n0\n0\n\n", file_r);
                printf("void specdata trans_spec\n");
                printf("4 noop %s . 'Acos(abs(Rdot))/DEGREE'\n0\n0\n\n", file_t);
                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");
        } else {
                fprintf(stderr, "Error: Failed to write one or more output files.\n");
        }

        cleanup_layers(layers, num_layers);
        free(total_t);
        free(total_rf);
        free(total_rb);

        return success ? EXIT_SUCCESS : EXIT_FAILURE;
}
Revision:	2.1
Committed:	Mon May 5 19:20:16 2025 UTC (5 weeks ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	feat(genglaze): New tool created by Taoning Wang added to distro
#	User	Rev	Content
1	greg	2.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			}