/*                   
FRESNEL EQUATIONS
This file contains Fresnel equations and supporting calculations.
Function summary:
	Fresnel amplitude equations (s and p polarized)
	Fresnel Irradiance equations (s and p polarized)
	Snell's law
	retardance and diattenuation calculations
	example code to test calculation speed
	test cases against known results and published literature

Decreasing phase convention assumed (n + i k)
Incident material assumed real-valued refractive index.
Transmitted irradiance defined immediately after interface.
Angles have units of radians.

"Inc" refers to incident material (assumed real-valued)
"Sub" refers to the transmitted (substrate) material (may be complex-valued)
"thetaInc" is the incident angle (assumed real-valued in range 0 <= thetaInc < pi/2) [units = radians]

timing and automated tests evaluated with build on Linux:
gcc -Wall Fresnel.c -o Fresnel.exe -lm
./Fresnel.exe




© 2022 Arizona Board of Regents on behalf of the University Of Arizona 
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Arizona (USA) Required Clauses:
1.1. Arbitration. The Parties hereby acknowledge that this Agreement may be subject to arbitration in accordance with applicable law and court rules.
1.2. Applicable Law and Venue. This Agreement shall be interpreted pursuant to the laws of the State of Arizona. Any arbitration or litigation between the Parties shall be conducted in Pima County, Arizona, and LICENSEE hereby submits to venue and jurisdiction in Pima County, Arizona.
1.3. Non-Discrimination. The Parties agree to be bound by state and federal laws and regulations governing equal opportunity and non-discrimination and immigration.
1.4. Appropriation of Funds. The Parties recognize that performance by ARIZONA may depend upon appropriation of funds by the State Legislature of Arizona. If the Legislature fails to appropriate the necessary funds, or if ARIZONA’S appropriation is reduced during the fiscal year, ARIZONA may cancel this Agreement without further duty or obligation. ARIZONA will notify LICENSEE as soon as reasonably possible after it knows of the loss of funds.
1.5. Conflict of Interest. This Agreement is subject to the provisions of A.R.S. 38-511 and other conflict of interest regulations.  Within three years of the EFFECTIVE DATE, ARIZONA may cancel this Agreement if any person significantly involved in initiating, negotiating, drafting, securing, or creating this Agreement for or on behalf of ARIZONA becomes an employee or consultant in any capacity of LICENSEE with respect to the subject matter of this Agreement.
*/

// 2022-10-12 - fix: conservation of energy for transmission into metal - Greg Smith
// 2022-Sept - initial implementation - Greg Smith

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double complex Snell(double n1, double complex n2, double thetaInc){
	// Snell's law with complex refractive index and complex angle of refraction
	// incident material (n1) is assumed to be non-absorbing (i.e. real-valued)
	//		if incident material has small absorption, use only real-part
	// because incident material is non-absorbing, the incident angle must also be real-valued

	// small bias needed to force the desired branch cut (casin) when n2 is real-valued
	if (abs(cimag(n2))<0.00000000000001){n2 = creal(n2)+0.00000000000001I;}

	return casin(n1*sin(thetaInc)/n2);
}

double complex rs(double n1, double complex n2, double thetaInc) {
	// s-polarized amplitude reflection
	complex double thetaRef = Snell(n1, n2, thetaInc);
	double term1 = n1*cos(thetaInc);
	complex double term2 = n2*ccos(thetaRef);
	return (term1-term2)/(term1+term2);
}

double complex rp(double n1, double complex n2, double thetaInc) {
	// p-polarized amplitude reflection
	complex double thetaRef = Snell(n1, n2, thetaInc);
	complex double term1 = n2*cos(thetaInc);
	complex double term2 = n1*ccos(thetaRef);
	return (term1-term2)/(term1+term2);
}

double complex ts(double n1, double complex n2, double thetaInc) {
	// s-polarized amplitude transmission
	complex double cosThetaRef = ccos(Snell(n1, n2, thetaInc));
	double cosThetaInc = cos(thetaInc);
	return 2.0*n1*cosThetaInc/(n1*cosThetaInc + n2*cosThetaRef);
}

double complex tp(double n1, double complex n2, double thetaInc) {
	// p-polarized amplitude transmission
	complex double cosThetaRef = ccos(Snell(n1, n2, thetaInc));
	double cosThetaInc = cos(thetaInc);
	return 2.0*n1*cosThetaInc/(n2*cosThetaInc + n1*cosThetaRef);
}

double Rs(double n1, double complex n2, double thetaInc) {
	// s-polarized irradiance reflection
	double rsVal = cabs(rs(n1, n2, thetaInc));
	return rsVal*rsVal;
}

double Rp(double n1, double complex n2, double thetaInc) {
	// p-polarized irradiance reflection
	double rpVal = cabs(rp(n1, n2, thetaInc));
	return rpVal*rpVal;
}

double Ts(double n1, double complex n2, double thetaInc) {
	// s-polarized irradiance transmission
	// absorbing materials require careful consideration of scale factor for energy conservation
	complex double cosThetaRef = ccos(Snell(n1, n2, thetaInc));
	double numerator = creal(n2 * cosThetaRef);
	double tsVal = cabs(ts(n1, n2, thetaInc));
	return tsVal*tsVal*numerator/n1/cos(thetaInc);
}

double Tp(double n1, double complex n2, double thetaInc) {
	// p-polarized irradiance transmission
	// absorbing materials require careful consideration of scale factor for energy conservation
	complex double cosThetaRef = ccos(Snell(n1, n2, thetaInc));
	double numerator = creal(n2)*creal(cosThetaRef) + cimag(n2)*cimag(cosThetaRef);
	double tpVal = cabs(tp(n1, n2, thetaInc));
	return tpVal*tpVal*numerator/n1/cos(thetaInc);
}


double Diattenuation(double irrad1, double irrad2) {
	// irrad1 and irrad2 are assumed to be irradiances of two orthogonal polarizations
	// diattenuation = 1 is a perfect polarizer
	return (fabs(irrad1-irrad2)/(irrad1+irrad2));
}

double Retardance (double complex amp1, double complex amp2){
	// amp1 and amp2 are assumed to be amplitudes of two orthogonal electric fields
	// retardance magnitude output is expressed in units of radians in range [-pi, +pi]
	double retardanceWrapped = carg(amp1) - carg(amp2);
	if (retardanceWrapped < -M_PI){retardanceWrapped += (2.0*M_PI);}
	if (retardanceWrapped >  M_PI){retardanceWrapped -= (2.0*M_PI);}
	return retardanceWrapped;
}



///////////////////// TEST CASES AND EXMAPLE ///////////////////////////
// threshold for test result accuracy
const double testPrecision = 0.001;	// some references have limited precision

// helper function to compare two complex numbers
int CompareComplex(double complex var1, double complex var2){
	return(
		fabs(creal(var1)-creal(var2)) < testPrecision && 
		fabs(cimag(var1)-cimag(var2)) < testPrecision
	);
};

// ideal glass n = 1.5 at normal incidence
// expected values for amplitude and irradiance known from theory
// see "Principles of Optics", Born and Wolf, 6th ed. sec. 1.5.2., eqn. 22-23, p.41
int testIdealGlassNormal(){
	double n1 = 1.0;
	complex double n2 = 1.5;
	double thetaInc = 0.0;
	return(
		CompareComplex(ts(n1,n2,thetaInc), 0.8) &&
		CompareComplex(tp(n1,n2,thetaInc), 0.8) &&
		CompareComplex(rs(n1,n2,thetaInc), -0.2) &&
		CompareComplex(rp(n1,n2,thetaInc), 0.2) &&
		fabs(Rs(n1,n2,thetaInc) - 0.04)<testPrecision &&
		fabs(Rp(n1,n2,thetaInc) - 0.04)<testPrecision &&
		fabs(Ts(n1,n2,thetaInc) - 0.96)<testPrecision &&
		fabs(Tp(n1,n2,thetaInc) - 0.96)<testPrecision &&
		fabs(Diattenuation(Rs(n1,n2,thetaInc), Rp(n1,n2,thetaInc)))<testPrecision &&
		fabs(Diattenuation(Ts(n1,n2,thetaInc), Tp(n1,n2,thetaInc)))<testPrecision &&
		fabs(Retardance(ts(n1,n2,thetaInc), tp(n1,n2,thetaInc)))<testPrecision &&
		fabs(Retardance(rs(n1,n2,thetaInc), rp(n1,n2,thetaInc))-M_PI)<testPrecision
	);
};

// Brewster angle with semi-arbitrary (real-valued) refractive indices (high-to-low)
// p-polarized reflection should be zero and transmitted irradiance 100%
int testBrewster(){
	double n1 = 1.72;
	complex double n2 = 1.15;
	double thetaB = atan(creal(n2)/n1);
	return(
		CompareComplex(rp(n1,n2,thetaB), 0.0) &&
		fabs(Rp(n1,n2,thetaB))<testPrecision &&
		fabs(Tp(n1,n2,thetaB) - 1.0)<testPrecision &&
		fabs(Diattenuation(Rs(n1,n2,thetaB), Rp(n1,n2,thetaB))-1.0)<testPrecision 
	);
};

// Critical angle with semi-arbitrary (real-valued) refractive indices (high-to-low)
// both polarization irradiance should be 100% reflected and 0% transmitted
int testCritical(){
	double n1 = 1.72;
	complex double n2 = 1.15;
	double thetaC = asin(creal(n2)/n1);
	return(
		fabs(Rs(n1,n2,thetaC) - 1.0)<testPrecision &&
		fabs(Rp(n1,n2,thetaC) - 1.0)<testPrecision &&
		fabs(Ts(n1,n2,thetaC))<testPrecision &&
		fabs(Tp(n1,n2,thetaC))<testPrecision &&
		fabs(Diattenuation(Rs(n1,n2,thetaC), Rp(n1,n2,thetaC)))<testPrecision
	);
};

// angles where Fresnel rhomb gives 1/8-wave retardance
// from text book: Born and Wolf, "Principles of Optics", 6th ed. section 1.5, below equation 63, p.50.
int testFresnelRhomb(){
	double n1 = 1.51;
	complex double n2 = 1.0;
	double theta1 = 0.84852090;	// 48 degrees, 37 arc-min
	double theta2 = 0.95324066;	// 54 degrees, 37 arc-min
	return(
		fabs(Retardance(rs(n1,n2,theta1), rp(n1,n2,theta1)) - (M_PI/4.0))<testPrecision &&
		fabs(Rs(n1,n2,theta1) + Ts(n1,n2,theta1) - 1.0)<testPrecision &&
		fabs(Rp(n1,n2,theta1) + Tp(n1,n2,theta1) - 1.0)<testPrecision &&
		fabs(Retardance(rs(n1,n2,theta2), rp(n1,n2,theta2)) - (M_PI/4.0))<testPrecision &&
		fabs(Rs(n1,n2,theta2) + Ts(n1,n2,theta2) - 1.0)<testPrecision &&
		fabs(Rp(n1,n2,theta2) + Tp(n1,n2,theta2) - 1.0)<testPrecision
	);
};

// metal reflection: gold at 582 nm using data interpolated from https://RefractiveIndex.info
// theory results from P. B. Johnson and R. W. Christy. "Optical constants of the noble metals", Phys. Rev. B 6, 4370-4379 (1972). https://doi.org/10.1103/PhysRevB.6.4370
int testFresnelMetal1(){
	double n1 = 1.0;
	complex double n2 = 0.29006+2.8628*I;
	double thetaInc = 0.87266463;	// 50 degrees
	return(
		fabs(Rs(n1,n2,thetaInc) - 0.92516)<testPrecision &&
		fabs(Rp(n1,n2,thetaInc) - 0.83233)<testPrecision &&
		fabs(Rs(n1,n2,thetaInc) + Ts(n1,n2,thetaInc) - 1.0)<testPrecision &&
		fabs(Rp(n1,n2,thetaInc) + Tp(n1,n2,thetaInc) - 1.0)<testPrecision
	);
};

// metal reflection: nickel at 632.8 nm
// theory results from M. Chiu, J. Lee, and D. Su. "Complex refractive-index measurement based on Fresnel's equations and the uses of heterodyne interferometry," Appl.Opt. 38, 4047-4052 (1999). https://doi.org/10.1364/AO.38.004047
int testFresnelMetal2(){
	double n1 = 1.0;
	complex double n2 = 2.007+3.781*I;
	double thetaInc = 1.0471976;	// 60 degrees
	return(
		CompareComplex(rs(n1,n2,thetaInc), -0.882373-0.183978*I) &&
		CompareComplex(rp(n1,n2,thetaInc), 0.462133+0.492465*I) &&
		fabs(Retardance(rp(n1,n2,thetaInc), rs(n1,n2,thetaInc)) + 2.53003)<testPrecision &&
		fabs(Rs(n1,n2,thetaInc) + Ts(n1,n2,thetaInc) - 1.0)<testPrecision &&
		fabs(Rp(n1,n2,thetaInc) + Tp(n1,n2,thetaInc) - 1.0)<testPrecision
	);
};


// high-level test functions
int testFresnel(){
	if (!testIdealGlassNormal()){printf("testIdealGlassNormal: failed\n");}
	if (!testBrewster()){printf("testBrewster: failed\n");}
	if (!testCritical()){printf("testCritical: failed\n");}
	if (!testFresnelRhomb()){printf("testFresnelRhomb: failed\n");}
	if (!testFresnelMetal1()){printf("testFresnelMetal1: failed\n");}
	if (!testFresnelMetal2()){printf("testFresnelMetal2: failed\n");}
	int combinedTest = 
		testIdealGlassNormal() &&
		testBrewster() &&
		testCritical() &&
		testFresnelRhomb() &&
		testFresnelMetal1() &&
		testFresnelMetal2();
	printf("all FRESNEL TESTS passed? (true = 1): %d\n",combinedTest);
	return combinedTest;
};


double randfrom(double min, double max) {
	// helper function for random number creation during speed testing
	double range = (max - min); 
	double div = RAND_MAX / range;
	return min + (rand() / div);
}


int main(int argc, char *argv[]) {
	// run test cases
	testFresnel();


	// TIMING TEST
	// note: to keep code simple, manual stopwatch must be used
	// note: Readability and calculation speed are both priorities for these
	// Fresnel functions. The boundary between the two priorities was chosen 
	// arbitrarily and I am not an expert in C, so the speed is not optimal.
	long nArray = 1000;		// 1D array of data to test possible speed

	// initialize values and arrays
	double *n1 = (double*) malloc(nArray * sizeof(double));
	double complex *n2 = (double complex*) malloc(nArray * sizeof(double complex));
	double *thetaInc = (double*) malloc(nArray * sizeof(double));
	double *result = (double*) malloc(nArray * sizeof(double));
	for (long ii=0; ii<nArray; ii++){
		n1[ii] = randfrom(0.0,2.0);
		n2[ii] = randfrom(0.0,2.0) + randfrom(0.0,10.0)*I;
		thetaInc[ii] = randfrom(0,M_PI/2.0);
	}

	// actual Fresnel calculation
	// result is written to array, but otherwise no analysis performed
	printf("random data qty: %ld\n",nArray);
	for (long ii=0; ii<nArray; ii++){
		result[ii] = Ts(n1[ii],n2[ii],thetaInc[ii]);
	}
	printf("Fresnel results complete!\n\n");

	// cleanup
	free(n1);
	free(n2);
	free(thetaInc);
	free(result);

	return EXIT_SUCCESS;
}