clear all;
accepted = 6.626068E-34
% First order data
data1 = [2.064 1.718 1.493 0.849 0.717];
data2 = [2.069 1.725 1.495 0.850 0.714];

% Second order data, with UV filter correction on green
data3 = [2.062 1.731 1.520 0.849 0.735];
data4 = [2.061 1.725 1.515 0.847 0.735];

% Accepted value for frequency of spectrum
freq = [8.20264E14, 7.40858E14, 6.87858E14, 5.48896E14, 5.18672E14];

% Charge of an electron. Also one electron volt
ev = 1.602E-19;

% Convert measured stopping voltage to energy in the 
data1 = data1' .* ev;
data2 = data2' .* ev;
data3 = data3' .* ev;
data4 = data4' .* ev;

freq = freq';
%data1 = data1';

% It's always good to take a look at the data
figure(1);
scatter(freq, data1);
hold on;
scatter(freq, data2);
scatter(freq, data3);
scatter(freq, data4);
hold off;
xlabel('freq');
ylabel('eV');

% Assuming that Plank's proportionality is accurate, the constant of
% proportionality for the line fitting the energy to frequency curve is
% Plank's constant, the y intercept is the work function of the material in
% the photo-diode. Use standard fit for first order polynomial and print to
% screen.
[p1] = fit(freq, data1, 'poly1');
[p2] = fit(freq, data2, 'poly1');
[p3] = fit(freq, data3, 'poly1');
[p4] = fit(freq, data4, 'poly1');


disp(sprintf('Planks: %.3E, %.3E', p1.p1, p2.p1));
disp(sprintf('Planks: %.3E, %.3E', p3.p1, p4.p1));

planks = mean([p1.p1 p2.p1 p3.p1 p4.p1]);
disp(sprintf('Mean value for planks constant: %.3E', planks));
disp(sprintf('Error: planks/accepted = %.3f', planks/accepted));

disp('');

data5 = (data1 + data2 + data3 + data4)/4;
p5 = fit(freq,data5,'poly1');
disp(sprintf('Method 2: Planks: %.3E', p5.p1));
disp(sprintf('Work Function: %.3E', abs(p5.p2)));
disp(sprintf('Error: planks/accepted = %.3f', p5.p1/accepted));

for i=1:5
    S(i) =  std([data1(i) data2(i) data3(i) data4(i)]);
end

data6 = data5 + S';
p6 = fit(freq,data6,'poly1');
error = abs(p5.p1-p6.p1);
disp(sprintf('Error : Planks: %.3E', error));