vert  = load('te_vertical_pol.txt');
horiz = load('tm_horiz_pol.txt');


hold off;
plot(vert(:,1),vert(:,2),'o','markersize',10);
hold on;
%plot(vert(:,1),vert(:,2));
plot(horiz(:,1),horiz(:,2),'or','markersize',10);
%plot(horiz(:,1),horiz(:,2),'r');

%grid on;
xlabel('Angle of incidence, {\theta}_i');
ylabel('Absolute intensity {\mu}W');
title('Reflected intensity as a function of angle of both polarizations');

% Fitting of data
theta_vert = pi/180*vert(:,1);
function_Rvert = @(a,x) a(1)*((cos(x)-sqrt(a(2)^2-sin(x).^2))./(cos(x)+sqrt(a(2)^2-sin(x).^2))).^2;
av0 = [vert(1,2)/0.04 1.54]; % Estimated n=1.54 from measurement of Brewster's angle
av = lsqcurvefit(function_Rvert, av0, theta_vert, vert(:,2));
plot(vert(:,1),function_Rvert(av,theta_vert));

theta_horiz = pi/180*horiz(:,1);
function_Rhoriz = @(a,x) a(1)*((a(2)^2*cos(x)-sqrt(a(2)^2-sin(x).^2))./(a(2)^2*cos(x)+sqrt(a(2)^2-sin(x).^2))).^2;
ah0 = [horiz(1,2)/0.04 1.54];
ah = lsqcurvefit(function_Rhoriz, ah0, theta_horiz, horiz(:,2));
plot(horiz(:,1),function_Rhoriz(av,theta_horiz),'r');
text(15,2100,'Fits are based on the Fresnel equations');
text(15,1950,'Blue: Vertical polarization (TE)');
text(15,1800,'Red: Horizontal polarization (TM)');