Untitled

/* singlay.c
   Calculates the reflectance of a single dielectric layer deposited
   on a glass substrate using the transfer matrix method.

   The program computes reflectance for wavelengths between
   300 nm and 800 nm in steps of 10 nm and writes the results
   to spectrum.csv for plotting in Excel.
*/



#include <math.h>
#include <stdio.h>
#include "clib.h"



#define PI 3.14159265 // Defining pi as a constant
#define ALPHA sqrt(8.854e-12 / (4*PI*1e-7)) // sqrt(epsilon_0/mu_0)



// Named indices for readability
#define INCIDENT 0
#define LAYER    1
#define GLASS 2
#define NUM_MATERIALS 3
#define NUM_WAVELENGTHS 51



// Defining a 2x2 matrix type
typedef fcomplex Matrix2x2[2][2];



// Procedure to calculate the wave vectors (k = 2pi/wavelenght) and storing them as an array.
void getWaveVector(double wavelengths[], double k[], int size) {

	for (int i = 0; i < size; i++) {
		k[i] = (2 * PI) / wavelengths[i];
	}
}



// Procedure to calculate the y values and storing them as an array.
void getAdmittances(double n[], double y[], int size) {

	for (int i = 0; i < size; i++) {
		y[i] = ALPHA * n[i];
	}
}



// Procedure to fill up the matrix elements 
void getM_l(Matrix2x2 M, double k, double y_l, double h_l) {
	double phase = k * h_l; // Defining the phase inside the cos / sin terms

	M[0][0] = complex(cos(phase), 0.0);        // m11
	M[0][1] = complex(0.0, sin(phase) / y_l);  // m12
	M[1][0] = complex(0.0, y_l * sin(phase));  // m21
	M[1][1] = complex(cos(phase), 0.0);         // m22
}



// Function to calculate reflectance from a 2x2 Matrix
double getReflectance(Matrix2x2 M, double y_0, double y_s) {
	// Calculating terms from Eq. 7 in coursework document.

	fcomplex term1 = rcmult(y_0, M[0][0]);           // Y_0 * m11
	fcomplex term2 = rcmult(y_0 * y_s, M[0][1]);     // Y_0 * Y_s * m12
	fcomplex term3 = M[1][0];                          // m21
	fcomplex term4 = rcmult(y_s, M[1][1]);            // Y_s * m22

	// Numerator = term1 + term2 - term3 - term4
	fcomplex numerator = cadd(term1, term2);
	numerator = csub(numerator, term3);
	numerator = csub(numerator, term4);

	// Denominator = term1 + term2 + term3 + term4
	fcomplex denominator = cadd(term1, term2);
	denominator = cadd(denominator, term3);
	denominator = cadd(denominator, term4);

	// r = numerator / denominator
	fcomplex r = cdiv(numerator, denominator);

	// Reflectance R = |r|²
	double R = ccabs(r) * ccabs(r);

	return R;
}



int main() {


	// Defining wavelengths [m] in an array
	double wavelengths[NUM_WAVELENGTHS];

	int j = 0; // Parameter to see what position in the array we're at

	for (double w = 300e-9; w < 810e-9 && j < NUM_WAVELENGTHS; w += 10e-9){  // Iterating from 300nm to 800nm in steps of 10nm
		wavelengths[j] = w;
		j++;
	}
	

	// Defining refractive indices in an array
	double n[NUM_MATERIALS] = {
		[INCIDENT] = 1.0,
		[LAYER] = 1.38,
		[GLASS] = 1.5
	};

	double d_l = 100e-9; // Thickness of Mgf_2 [m]


	// Declaring empty variables 
	double h_l;

	double k[NUM_WAVELENGTHS]; // Empty array for wave vectors
	double y[NUM_MATERIALS]; // Empty array for admittances 


	// Calculating the wavevectors
	getWaveVector(wavelengths, k, NUM_WAVELENGTHS);


	// Calculating admittances
	getAdmittances(n, y, NUM_MATERIALS);


	// Calculating the optical thickness
	h_l = n[LAYER] * d_l;


	// Generating array of 2x2 M_l for each wavelength
	Matrix2x2 m_l[NUM_WAVELENGTHS]; // Array of 2x2 matrices

	for (int i = 0; i < NUM_WAVELENGTHS; i++) {
		getM_l(m_l[i], k[i], y[LAYER], h_l);
	}


	// Defining a csv file with headers
	FILE* file = fopen("spectrum.csv", "w");
	fprintf(file, "Wavelength_nm,Reflectance\n");


	// Getting the reflectances for each wavelength.
	double reflectance[NUM_WAVELENGTHS];

	for (int i = 0; i < NUM_WAVELENGTHS; i++) {

		reflectance[i] = getReflectance(m_l[i], y[INCIDENT], y[GLASS]);

		// Write wavelength and reflectance to CSV file for plotting	
		fprintf(file, "%lf,%lf\n", wavelengths[i] * 1e9, reflectance[i]);
	}


	fclose(file); // Close output file


	return 0;
}