theta0 = pi/2;
theta = linspace(0,pi);

w = cos(theta-theta0).*cos(theta);
wsq = w.^2;
subplot(211);
plot(180/pi*theta,wsq);
grid on;
xlabel('Relative angle between input and central polarizers, \theta');
ylabel('Relative transmission');
title('Expected intensity at detector, ideal case of {\theta}_0 = 90^{\circ}');

subplot(212);
data = load('central_polarizer.txt');
data(:,1) = data(:,1)+45;
plot(data(:,1),data(:,2),'o','markersize',10);
xlabel('Relative angle between input and central polarizers, \theta');
ylabel('Absolute transmission ({\mu}W)');
grid on;
hold on;

% curve fitting
function_polarizer = @(a,x) a(1)*cos(x-a(2)).^2.*cos(x).^2;
x = pi/180*data(:,1); % Turn into radians
y = data(:,2);
a0 = [1600 pi/2]; % Initial guess;
options = optimset('TolFun',1e-8);
a = lsqcurvefit(function_polarizer, a0, x, y, [], [], options);
theta = linspace(pi/4,3*pi/4);
fit = function_polarizer(a,theta);
plot(180/pi*theta,fit,'r');
headline = sprintf('Measured data and fitted curve with {\theta}_0 = %f0.5^{\circ}', 180/pi*a(2));
title('Measured data and fitted curve with {\theta}_0 = 92.4^{\circ}');