% This program is used to optimize multilayer thin film structure in order 
% to get high transmittance or reflectance in specified regions. 

% The equations of front reflection, back reflection, front transmittance,
% back transmittance are all provided in function Caltrans. Users can
% modify the CalError code to ultilize these optical properties and design
% other applications, such as DBR.

% The program is adapted from TransferMatrix matlab code and previous
% DBR_transp matlab code on January 24, 2016. 
% Reference imformation: 
% Griffiths "Intro to Electrodynamics 3rd Ed. Eq. 9.3
% J. Appl. Phys Vol 86, No. 1 (1999) p.487
% JAP 93 No. 7 p. 3693

function DBR_transp1
% --------  User parameters input ---------- %
% Set high-transparency spectra region
global lambda;
lambda = 300:10:750;
% Set Layer Structure
layers = {};
% Set initial thickness 
iniThickness = [];
% Set low bound & up bound to thickness
LB = [];
UB = [];
% --------  End of input section ----------- %

% Read refractive index from excel file
global n;
n = zeros(size(layers,2),size(lambda,2));
for index = 1:size(layers,2)
    n(index,:) = LoadRefrIndex(layers{index},lambda);
end

% Find optimal thickness
% X is the local optimal thickness
% RESNORM is the total deviation from desired transmittance curve
[X RESNORM] = lsqnonlin(@CalError,iniThickness,LB,UB)

% define wavelength region for ploting curve
wavelength =300:2:850;
% get new n value
n = zeros(size(layers,2),size(wavelength,2));
for index = 1:size(layers,2)
    n(index,:) = LoadRefrIndex(layers{index},wavelength);
end
% Calculate transmittace spetra
newTrans = CalTrans(X, wavelength);
oldTrans = CalTrans(iniThickness,wavelength);
% plot initial curve and optimized curve
plot(wavelength,oldTrans,wavelength,newTrans);
% hold on;
fid = fopen('wavelength.txt','wt');
fprintf(fid,'%g\n',wavelength);
fclose(fid);
fid = fopen('new.txt','wt');
fprintf(fid,'%g\n',newTrans);
fclose(fid);
legend('old','new');
% End of the main function



%  --------------- Help Functions --                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    -------------------  %

function Error = CalError(iniThickness)
global lambda;
calRegion = lambda;
Transmittance = CalTrans(iniThickness,calRegion);
dT = 0.002;
stepT = 1:length(Transmittance);
Error = abs((Transmittance(stepT)-1)/dT);


function Trans = CalTrans(thickness, Region)
lambda1 = Region;
tempT=thickness;
% interfacial glass layer
t(1) = 0;
% glass layer in nm
t(2) = 1e6;
for i=1:length(tempT)
t(i+2)=tempT(i);
end
% air thickness equals 100 nm
t(length(tempT)+3)=100;
global n;
% Calculate Incoherent power transmission through substrate
T_glass=abs(4*1*n(1,:)./(1+n(1,:)).^2); 
R_glass=abs((1-n(1,:))./(1+n(1,:))).^2;

% initialize Reflectance and transmittance data
R=lambda1*0;
T=lambda1*0;
T2=lambda1*0;
bR=lambda1*0;
bT=lambda1*0;
% Calculate transfer matrix
for l = 1:length(lambda1)
    % Calculate transfer matrices for incoherent reflection and transmission at the first interface
    S=I_mat(n(1,l),n(2,l));
    for matindex=2:(length(t)-1)
        S=S*L_mat(n(matindex,l),t(matindex),lambda1(l))*I_mat(n(matindex,l),n(matindex+1,l));
    end
    R(l)=abs(S(2,1)/S(1,1))^2; %JAP Vol 86 p.487 Eq 9 Power Reflection from layers other than substrate
    T(l)=abs(2/(1+n(1,l)))/sqrt(1-R_glass(l)*R(l)); %Transmission of field through glass substrate Griffiths 9.85 + multiple reflection geometric series
    T2(l) = abs(1/S(1,1))^2/n(1,l);
    bR(l)=abs(-S(1,2)/S(1,1))^2;
    bT(l)=abs(det(S)/S(1,1))^2*n(1,1);
end
% Overall Reflection from device with incoherent reflections at first
% interface (typically air-glass)
Reflection = R_glass+T_glass.^2.*R./(1-R_glass.*R);
% Overall front transmittance with incoherent reflections at first 
% interface
Transmission = T_glass.*T2./(1-R_glass.*R);
% Overall back reflection with incoherent reflections at last
% interface
bReflection = bR+bT.*R_glass.*T2./(1-R_glass.*R);
% Overall back transmittance with incoherent reflections at last
% interface
bTransmittance = bT.*T_glass./(1-R_glass.*R);

% Return optical property
Trans= Transmission;

function ntotal = LoadRefrIndex(name,wavelengths)

%Data in IndRefr, Column names in IndRefr_names
[IndRefr,IndRefr_names]=xlsread('Index_of_Refraction_library.xls');

% Load index of refraction data in spread sheet, will crash if misspelled
file_wavelengths=IndRefr(:,strmatch('Wavelength',IndRefr_names));
n=IndRefr(:,strmatch(strcat(name,'_n'),IndRefr_names));
k=IndRefr(:,strmatch(strcat(name,'_k'),IndRefr_names));

% Interpolate/Extrapolate data linearly to desired wavelengths
n_interp=interp1(file_wavelengths, n, wavelengths, 'nearest', 'extrap');
k_interp=interp1(file_wavelengths, k, wavelengths, 'nearest', 'extrap');

%Return interpolated complex index of refraction data
ntotal = n_interp+1i*k_interp;

function I = I_mat(n1,n2)
r=(n1-n2)/(n1+n2);
t=2*n1/(n1+n2);
I=[1 r; r 1]/t;
function L = L_mat(n,d,lambda)
xi=2*pi*n/lambda;
L=[exp(-1i*xi*d) 0; 0 exp(1i*xi*d)];