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 LARGE_DOUBE = 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;


static 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) ? -LARGE_DOUBE : (-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;
                }
        }
        free(rho_0);
        free(refrac);
        free(abs_coeff);
}


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");
}


static void free_layer_data(GlazingLayer *layer) {
    if (!layer) return;
    
    free(layer->filename);
    free(layer->t_0);
    free(layer->rf_0);
    free(layer->rb_0);
    free(layer->t_lambda_theta);
    free(layer->rf_lambda_theta);
    free(layer->rb_lambda_theta);
    
    // Set pointers to NULL to prevent double-free errors
    layer->filename = NULL;
    layer->t_0 = NULL;
    layer->rf_0 = NULL;
    layer->rb_0 = NULL;
    layer->t_lambda_theta = NULL;
    layer->rf_lambda_theta = NULL;
    layer->rb_lambda_theta = NULL;
}


void cleanup_layers(GlazingLayer *layers, int num_layers) {
         if (!layers) return;
         for (int i = 0; i < num_layers; ++i) {
                if (layers[i].filename)
                        free_layer_data(&layers[i]);
        }
        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");
                                cleanup_layers(layers, num_layers);
                                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);
                                cleanup_layers(layers, num_layers);
                                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_%s\n", prefix);
                printf("4 noop %s . 'Acos(Rdot)/DEGREE'\n0\n0\n\n", file_r);
                printf("void specdata trans_spec_%s\n", prefix);
                printf("4 noop %s . 'Acos(abs(Rdot))/DEGREE'\n0\n0\n\n", file_t);
                printf("void WGMDfunc glaze_mat_%s\n13\n\trefl_spec_%s 1 0 0\n\ttrans_spec_%s 1 0 0\n\tvoid\n\t0 0 1 .\n0\n9 0 0 0 0 0 0 0 0 0\n", prefix, prefix, prefix);
        } 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.2
Committed:	Tue Jun 10 23:13:13 2025 UTC (5 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad6R0, HEAD
Changes since 2.1:	+37 -9 lines
Log Message:	fix(genglaze): TW fixed a memory leak and some minor improvements
#	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	greg	2.2	const double LARGE_DOUBE = 3.40282347E+38;
20	greg	2.1	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	greg	2.2	static inline double polynomial_5(double x, const double coeffs[5]) {
57	greg	2.1	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	greg	2.2	double exp_arg = (cos_phi_prime < DBL_EPSILON) ? -LARGE_DOUBE : (-abs_coeff[iwvl] * layer->thickness_m / cos_phi_prime);
151	greg	2.1	/* 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	greg	2.2	free(rho_0);
182			free(refrac);
183			free(abs_coeff);
184	greg	2.1	}
185
186
187			void angular_coated(GlazingLayer *layer, const int nwvl) {
188			for (int itheta = 0; itheta < NTHETA; ++itheta) {
189			const double phi = THETAS[itheta] * RAD_PER_DEG;
190			const double cos_phi = cos(phi);
191
192			for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
193			const double t_0_lambda = layer->t_0[iwvl];
194			const double *coeffs_tau = (t_0_lambda > 0.645) ? COEFFS_TAU_CLEAR : COEFFS_TAU_BRONZE;
195			const double *coeffs_rho = (t_0_lambda > 0.645) ? COEFFS_RHO_CLEAR : COEFFS_RHO_BRONZE;
196
197			const double tau_bar = polynomial_5(cos_phi, coeffs_tau);
198			const double rho_bar_term = polynomial_5(cos_phi, coeffs_rho);
199			const double rho_bar = rho_bar_term - tau_bar;
200
201			const int flat_idx = itheta * nwvl + iwvl;
202			layer->t_lambda_theta[flat_idx] = t_0_lambda * tau_bar;
203			layer->rf_lambda_theta[flat_idx] = layer->rf_0[iwvl] * (1.0 - rho_bar) + rho_bar;
204			layer->rb_lambda_theta[flat_idx] = layer->rb_0[iwvl] * (1.0 - rho_bar) + rho_bar;
205			}
206			}
207			}
208
209
210			void multi_layer_calc(
211			GlazingLayer *layers,
212			int nlayers,
213			int nwvl,
214			double *total_t,
215			double *total_rf,
216			double *total_rb)
217			{
218			if (nlayers <= 0)
219			return;
220
221			size_t total_size = (size_t)NTHETA * nwvl * sizeof(double);
222
223			memcpy(total_t, layers[0].t_lambda_theta, total_size);
224			memcpy(total_rf, layers[0].rf_lambda_theta, total_size);
225			memcpy(total_rb, layers[0].rb_lambda_theta, total_size);
226
227			if (nlayers == 1)
228			return;
229
230			double *prev_t = malloc(total_size);
231			double *prev_rf = malloc(total_size);
232			double *prev_rb = malloc(total_size);
233
234			if (!prev_t \|\| !prev_rf \|\| !prev_rb)
235			perror("Failed to allocate temporary storage in multi_layer");
236
237			for (int j = 1; j < nlayers; ++j) {
238			memcpy(prev_t, total_t, total_size);
239			memcpy(prev_rf, total_rf, total_size);
240			memcpy(prev_rb, total_rb, total_size);
241
242			double *t_j = layers[j].t_lambda_theta;
243			double *rf_j = layers[j].rf_lambda_theta;
244			double *rb_j = layers[j].rb_lambda_theta;
245
246			for (int itheta = 0; itheta < NTHETA; ++itheta) {
247			for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
248			const int idx = itheta * nwvl + iwvl;
249			double denominator = 1.0 - rf_j[idx] * prev_rb[idx];
250			if (fabs(denominator) < DBL_EPSILON) {
251			denominator = DBL_EPSILON;
252			}
253			total_t[idx] = prev_t[idx] * t_j[idx] / denominator;
254			total_rf[idx] = prev_rf[idx] + prev_t[idx] * prev_t[idx] * rf_j[idx] / denominator;
255			total_rb[idx] = rb_j[idx] + t_j[idx] * t_j[idx] * prev_rb[idx] / denominator;
256			}
257			}
258			}
259			free(prev_t);
260			free(prev_rf);
261			free(prev_rb);
262			}
263
264
265			int add_layer(GlazingLayer *layers_array, int count, int *capacity,
266			const char *filename, double thickness_m, int is_monolithic) {
267			if (count >= capacity) {
268			capacity = (capacity == 0) ? 4 : capacity 2;
269			GlazingLayer new_layers = realloc(layers_array, capacity sizeof(GlazingLayer));
270			if (!new_layers) {
271			perror("Failed to reallocate memory for layers");
272			return 0;
273			}
274			*layers_array = new_layers;
275			}
276
277			memset(&((layers_array)[count]), 0, sizeof(GlazingLayer));
278
279			#ifdef _WIN32
280			(layers_array)[count].filename = _strdup(filename);
281			#else
282			(layers_array)[count].filename = strdup(filename);
283			#endif
284
285			if (!(layers_array)[count].filename) {
286			perror("Failed to duplicate filename");
287			return 0;
288			}
289			(layers_array)[count].is_mono = is_monolithic;
290			(layers_array)[count].thickness_m = thickness_m;
291			(layers_array)[count].t_0 = NULL;
292			(layers_array)[count].rf_0 = NULL;
293			(layers_array)[count].rb_0 = NULL;
294			(layers_array)[count].t_lambda_theta = NULL;
295			(layers_array)[count].rf_lambda_theta = NULL;
296			(layers_array)[count].rb_lambda_theta = NULL;
297
298			(*count)++;
299			return 1;
300			}
301
302
303			int interpolate(GlazingLayer *layer, const double wvl_start, const double wvl_end, const double wvl_interval, const int nwvl) {
304			DATARRAY *dp = getdata(layer->filename);
305			if (!dp) {
306			fprintf(stderr, "Error: Cannot open file '%s'\n", layer->filename);
307			return 0;
308			}
309
310			layer->t_0 = malloc(nwvl * sizeof(double));
311			layer->rf_0 = malloc(nwvl * sizeof(double));
312			layer->rb_0 = malloc(nwvl * sizeof(double));
313
314			double wvl = wvl_start;
315			int i = 0;
316			while (wvl <= wvl_end) {
317			double t_pt[2] = {2., wvl};
318			layer->t_0[i] = datavalue(dp, t_pt);
319			double rf_pt[2] = {0., wvl};
320			layer->rf_0[i] = datavalue(dp, rf_pt);
321			double rb_pt[2] = {1., wvl};
322			layer->rb_0[i] = datavalue(dp, rb_pt);
323			wvl = wvl + wvl_interval;
324			i = i + 1;
325
326			}
327			freedata(dp);
328			return 1;
329			}
330
331
332			int write_output_file(
333			const char *tfilename,
334			const char *rfilename,
335			const double *tdata,
336			const double *rfdata,
337			const double *rbdata,
338			const double wvl_start,
339			const double wvl_end,
340			const int nwvl,
341			int argc, char *argv[])
342			{
343			FILE *tfp = fopen(tfilename, "w");
344			FILE *rfp = fopen(rfilename, "w");
345			if (!tfp \|\| !rfp) {
346			fprintf(stderr, "Error: Cannot open output files\n");
347			return(0);
348			}
349
350			fprintf(tfp, "# ");
351			for (int i = 0; i< argc; ++i) {
352			fprintf(tfp, "%s ", argv[i]);
353			}
354			fprintf(tfp, "\n2\n0 0 %d", NTHETA);
355			for (int i = 0; i < NTHETA; i++) {
356			fprintf(tfp, " %d", (int)THETAS[i]);
357			}
358			fprintf(tfp, "\n%.0f %.0f %d\n", wvl_start, wvl_end, nwvl);
359
360			for (int itheta = 0; itheta < NTHETA; ++itheta) {
361			for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
362			fprintf(tfp, "%.6f\n", tdata[itheta*nwvl+iwvl]);
363			}
364			}
365			fclose(tfp);
366
367			fprintf(rfp, "# ");
368			for (int i = 0; i< argc; ++i) {
369			fprintf(rfp, "%s ", argv[i]);
370			}
371			fprintf(rfp, "\n2\n0 0 %d", NTHETA * 2 - 1);
372			for (int i = 0; i < NTHETA; i++) {
373			fprintf(rfp, " %d", (int)THETAS[i]);
374			}
375			for (int i = NTHETA-2; i >= 0; --i) {
376			fprintf(rfp, " %d", (int)(180 - THETAS[i]));
377			}
378			fprintf(rfp, "\n%.0f %.0f %d\n", wvl_start, wvl_end, nwvl);
379
380			for (int itheta = 0; itheta < NTHETA; ++itheta) {
381			for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
382			fprintf(rfp, "%.6f\n", rfdata[itheta * nwvl + iwvl]);
383			}
384			}
385			for (int itheta = NTHETA - 2; itheta >= 0 ; --itheta) {
386			for (int iwvl = 0; iwvl < nwvl; ++iwvl) {
387			fprintf(rfp, "%.6f\n", rbdata[itheta * nwvl + iwvl]);
388			}
389			}
390
391			fclose(rfp);
392			return(1);
393			}
394
395
396			void print_usage() {
397			fprintf(stderr, "Usage: %s [-m monolithic_layer.dat thickness \| -c coated_layer.dat] ...\n", progname);
398			fprintf(stderr, "Calculate multi-layer glazing optics from spectral data files.\n");
399			fprintf(stderr, "Options:\n");
400			fprintf(stderr, " -m <filename> thickness Specify an uncoated (monolithic) glazing layer .dat file and its thickness (meter).\n");
401			fprintf(stderr, " -c <filename> Specify a coated or laminate glazing layer .dat file.\n");
402			fprintf(stderr, " -s start_wvl end_wvl interval Specify wavelength range and interval.\n");
403			fprintf(stderr, " -p prefix Specify prefix name to the output files.\n");
404			fprintf(stderr, " -h, --help Show this help message.\n");
405			fprintf(stderr, "Layer order is determined by the sequence of options on the command line.\n");
406			fprintf(stderr, "Output Files: prefix_t.dat, prefix_r.dat\n");
407			}
408
409	greg	2.2
410			static void free_layer_data(GlazingLayer *layer) {
411			if (!layer) return;
412
413			free(layer->filename);
414			free(layer->t_0);
415			free(layer->rf_0);
416			free(layer->rb_0);
417			free(layer->t_lambda_theta);
418			free(layer->rf_lambda_theta);
419			free(layer->rb_lambda_theta);
420
421			// Set pointers to NULL to prevent double-free errors
422			layer->filename = NULL;
423			layer->t_0 = NULL;
424			layer->rf_0 = NULL;
425			layer->rb_0 = NULL;
426			layer->t_lambda_theta = NULL;
427			layer->rf_lambda_theta = NULL;
428			layer->rb_lambda_theta = NULL;
429			}
430
431
432	greg	2.1	void cleanup_layers(GlazingLayer *layers, int num_layers) {
433			if (!layers) return;
434			for (int i = 0; i < num_layers; ++i) {
435			if (layers[i].filename)
436	greg	2.2	free_layer_data(&layers[i]);
437	greg	2.1	}
438			free(layers);
439			}
440
441
442			int main(int argc, char *argv[])
443			{
444			GlazingLayer *layers = NULL;
445			int num_layers = 0;
446			int layer_capacity = 0;
447			int success = 1;
448			double wvl_start_nm = 380.;
449			double wvl_end_nm = 780.;
450			double wvl_interval_nm = 5.;
451			int nwvl = 81;
452			double thickness_m = 0.003;
453			char *filename;
454			char *prefix = "unnamed";
455			char file_t[MAX_NAME];
456			char file_r[MAX_NAME];
457			progname = argv[0];
458			if (argc <= 1) {
459			print_usage();
460			return EXIT_FAILURE;
461			}
462
463			for (int i=1; i < argc; ++i) {
464			switch (argv[i][1]) {
465			case 'm':
466			filename = argv[++i];
467			thickness_m = atof(argv[++i]);
468			if (!add_layer(&layers, &num_layers, &layer_capacity, filename, thickness_m, 1)) {
469	greg	2.2	cleanup_layers(layers, num_layers);
470			return EXIT_FAILURE;
471	greg	2.1	}
472			break;
473			case 'c':
474			filename = argv[++i];
475			if (!add_layer(&layers, &num_layers, &layer_capacity, filename, thickness_m, 0)) {
476			cleanup_layers(layers, num_layers);
477			return EXIT_FAILURE;
478			}
479			break;
480			case 'p':
481			prefix = argv[++i];
482			break;
483			case 's':
484			wvl_start_nm = atof(argv[++i]);
485			wvl_end_nm = atof(argv[++i]);
486			wvl_interval_nm = atof(argv[++i]);
487			if (wvl_start_nm > wvl_end_nm) {
488			fprintf(stderr, "Starting wavelength > End wavelength\n");
489	greg	2.2	cleanup_layers(layers, num_layers);
490	greg	2.1	return EXIT_FAILURE;
491			}
492			if (((int)(wvl_end_nm - wvl_start_nm) % (int)wvl_interval_nm) > 0) {
493			fprintf(stderr,
494			"Error: Wavelength range (%f to %f nm) must be evenly divisible by the interval (%f nm).\n",
495			wvl_start_nm, wvl_end_nm, wvl_interval_nm);
496	greg	2.2	cleanup_layers(layers, num_layers);
497	greg	2.1	return EXIT_FAILURE;
498			}
499			break;
500			case 'h':
501			print_usage();
502			return 1;
503			default:
504			break;
505			}
506			}
507			snprintf(file_t, MAX_NAME, "%s_t.dat", prefix);
508			snprintf(file_r, MAX_NAME, "%s_r.dat", prefix);
509
510			if (fabs(wvl_start_nm - wvl_end_nm) < DBL_EPSILON && wvl_interval_nm > DBL_EPSILON) {
511			nwvl = 1;
512			} else {
513			nwvl = (int)((wvl_end_nm - wvl_start_nm) / wvl_interval_nm + 1.5);
514			}
515
516			if (num_layers == 0) {
517			fprintf(stderr, "Error: No layers specified.\n");
518			print_usage();
519			return EXIT_FAILURE;
520			}
521
522			for (int i = 0; i < num_layers; ++i) {
523			if (!interpolate(&layers[i], wvl_start_nm, wvl_end_nm, wvl_interval_nm, nwvl)) {
524			fprintf(stderr, "Error: Failed to parse or interpolate data for layer %d.\n", i + 1);
525			success = 0;
526			break;
527			}
528			layers[i].t_lambda_theta = malloc(NTHETA * nwvl * sizeof(double));
529			layers[i].rf_lambda_theta = malloc(NTHETA * nwvl * sizeof(double));
530			layers[i].rb_lambda_theta = malloc(NTHETA * nwvl * sizeof(double));
531
532			if (layers[i].is_mono) {
533			angular_monolithic(&layers[i], nwvl);
534			} else {
535			angular_coated(&layers[i], nwvl);
536			}
537			}
538
539			if (!success) {
540			cleanup_layers(layers, num_layers);
541			return EXIT_FAILURE;
542			}
543
544
545			size_t total_flat_size = (size_t)NTHETA * nwvl * sizeof(double);
546			double *total_t = malloc(total_flat_size);
547			double *total_rf = malloc(total_flat_size);
548			double *total_rb = malloc(total_flat_size);
549
550			if (!total_t \|\| !total_rf \|\| !total_rb) {
551			perror("Failed to allocate memory for final results");
552			cleanup_layers(layers, num_layers);
553			free(total_t);
554			free(total_rf);
555			free(total_rb);
556			return EXIT_FAILURE;
557			}
558
559			multi_layer_calc(layers, num_layers, nwvl, total_t, total_rf, total_rb);
560
561			success &= write_output_file(file_t, file_r, total_t, total_rf, total_rb, wvl_start_nm, wvl_end_nm, nwvl, argc, argv);
562
563			if (success) {
564			printf("# ");
565			for (int i = 0;i < argc; ++i) {
566			printf("%s ", argv[i]);
567			}
568	greg	2.2	printf("\nvoid specdata refl_spec_%s\n", prefix);
569	greg	2.1	printf("4 noop %s . 'Acos(Rdot)/DEGREE'\n0\n0\n\n", file_r);
570	greg	2.2	printf("void specdata trans_spec_%s\n", prefix);
571	greg	2.1	printf("4 noop %s . 'Acos(abs(Rdot))/DEGREE'\n0\n0\n\n", file_t);
572	greg	2.2	printf("void WGMDfunc glaze_mat_%s\n13\n\trefl_spec_%s 1 0 0\n\ttrans_spec_%s 1 0 0\n\tvoid\n\t0 0 1 .\n0\n9 0 0 0 0 0 0 0 0 0\n", prefix, prefix, prefix);
573	greg	2.1	} else {
574			fprintf(stderr, "Error: Failed to write one or more output files.\n");
575			}
576
577			cleanup_layers(layers, num_layers);
578			free(total_t);
579			free(total_rf);
580			free(total_rb);
581
582			return success ? EXIT_SUCCESS : EXIT_FAILURE;
583			}