clear all;
c = 3 * 10^8; %Speed of light in millimeters
h = (6.62606885*10^-34)/(2*pi); % h bar technically.
%h = (4.13566733 * 10^-15)/(2*pi); % h bar in eV-s
e = 1.60217646 * 10^-19; % Charge of an electron in coulombs
em = 9.10938188 * 10^-31 ; % Electron mass in killograms 
%emc2 = e * em * c^2; %electron mass time speed of light squared 


%%%%%%%%%%%%%%%%%%%%%%
% Calculating projected diameter via geometry of the bulb
% Using known circular radius of curvature of the bulb at 66.7mm (to outside of bulb)
% See Ediffraction drawing1.tif for references

L = .130; % Length of bulb in meters. +/- 2mm
D1 = sym('2*L*tan(theta)');
theta = sym('atan(R/(L-x))');
x = sym('0.0675 - a');
a = sym('0.0675 * cos( asin(R/0.0675) )');
R = 33; % Need a temporary value for R to build final sym
Dsym = eval(eval(eval(D1))); % Build the final equation which is a function of R

eval(Dsym);

V5000.Inner = [24.59 25.10 22.45 24.97 22.90 23.48 23.88 22.15 22.75];
V5000.Outer = [38.9 38.45 38.27 38.65 38.22 38.43 38.44 38.35 38.10 38.17];
% Adjust to meters from mm
V5000.Inner = V5000.Inner ./ 1000; 
V5000.Outer = V5000.Outer ./ 1000;

V4500.Inner = [25.37 24.49 23.98 24.22 24.06 24.68 24.42 24.15 24.13 24.51];
V4500.Outer = [41.88 41.33 42.15 40.19 41.13 40.42 40.33 40.25 40.94 40.63];
V4500.Inner = V4500.Inner ./ 1000; 
V4500.Outer = V4500.Outer ./ 1000;

V4000.Inner = [25.58 25.36 26.69 25.20 25.36 25.26 24.92 25.78 25.69 24.88];
V4000.Outer = [42.92 42.86 43.40 42.45 42.98 43.26 42.51 42.88 42.79 42.60];
V4000.Inner = V4000.Inner ./ 1000; 
V4000.Outer = V4000.Outer ./ 1000;

V3500.Inner = [28.72 29.03 29.21 27.18 26.57 27.52 26.60 25.94 26.28 26.62 ];
V3500.Outer = [48.32 46.65 49.77 48.90 47.20 48.42 44.84 46.28 45.47 45.20];
V3500.Inner = V3500.Inner ./ 1000; 
V3500.Outer = V3500.Outer ./ 1000;

V3000.Inner = [32.49 34.82 33.65 32.83 35.45 31.54 32.56 33.40 32.86 32.70];
V3000.Outer = [54.01 54.02 53.45 54.71 54.25 53.74 54.01 53.87 53.93 54.65 ];
V3000.Inner = V3000.Inner ./ 1000; 
V3000.Outer = V3000.Outer ./ 1000;

f(1) = figure(1);
clf(f(1));
colordef(f(1), 'black');
set(f(1), 'WindowStyle', 'docked');

hold on;
plot(V5000.Inner, 'r');
plot(V5000.Outer, 'r');

plot(V4500.Inner, 'y');
plot(V4500.Outer, 'y');

plot(V4000.Inner, 'g');
plot(V4000.Outer, 'g');

plot(V3500.Inner, 'b');
plot(V3500.Outer, 'b');

plot(V3000.Inner, 'w');
plot(V3000.Outer, 'w');
hold off;
title('e Diffraction Lab Raw Data');
xlabel('sample number');
ylabel('ring diameter in mm');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the mean values, standard deviation and 
display('Mean values by accelerating voltage');

V5000.Imean = mean(V5000.Inner);
V5000.Omean = mean(V5000.Outer);
V5000.ISE = std(V5000.Inner)./sqrt(5);
V5000.OSE = std(V5000.Outer)./sqrt(5);
display(sprintf('5.0KV - Inner Ring - mean: %.4Em, SE: %.4E',V5000.Imean, V5000.ISE));
display(sprintf('5.0KV - Outer Ring - mean: %.4Em, SE: %.4E',V5000.Imean, V5000.OSE));
D = V5000.Imean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('5.0KV - inner ring calculations for d = %.3E', d));
d_inner(1) = d;

D = V5000.Omean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('5.0KV - outer ring calculations for d = %.3E', d));
d_outer(1) =d ;

%% Propagate error estimate
D = V5000.Imean + V5000.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_hi(1) = d_inner(1) - d;

D = V5000.Imean - V5000.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_lo(1) = d_inner(1) - d;

D = V5000.Omean + V5000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_hi(1) = d_outer(1) - d;

D = V5000.Omean - V5000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_lo(1) = d_outer(1) - d;

display(' ')

V4500.Imean = mean(V4500.Inner);
V4500.Omean = mean(V4500.Outer);
V4500.ISE = std(V4500.Inner)./sqrt(5);
V4500.OSE = std(V4500.Outer)./sqrt(5);
display(sprintf('4.5KV - Inner Ring - mean: %.4Em, SE: %.4E',V4500.Imean, V4500.ISE));
display(sprintf('4.5KV - Outer Ring - mean: %.4Em, SE: %.4E',V4500.Imean, V4500.OSE));

D = V4500.Imean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4500*em*c^2));
display(sprintf('4.5KV - inner ring calculations for d = %.3E', d));
d_inner(2) = d;

D = V4500.Omean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4500*em*c^2));
display(sprintf('4.5KV - outer ring calculations for d = %.3E', d));
d_outer(2) =d ;

%% Propagate error estimate 4500
D = V4500.Imean + V4500.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_hi(2) = d_inner(2) - d;

D = V4500.Imean - V4500.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_lo(2) = d_inner(2) - d;

D = V4500.Omean + V4500.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_hi(2) = d_outer(2) - d;

D = V5000.Omean - V5000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*5000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_lo(2) = d_outer(2) - d;

display(' ')

display(sprintf('4.5KV - inner ring - d = %.3E', d));
display(' ')

%% V4000
V4000.Imean = mean(V4000.Inner);
V4000.Omean = mean(V4000.Outer);
V4000.ISE = std(V4000.Inner)./sqrt(5);
V4000.OSE = std(V4000.Outer)./sqrt(5);
display(sprintf('4.0KV - Inner Ring - mean: %.4Em, SE: %.4E',V4000.Imean, V4000.ISE));
display(sprintf('4.0KV - Outer Ring - mean: %.4Em, SE: %.4E',V4000.Imean, V4000.OSE));

D = V4000.Imean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4000*em*c^2));
display(sprintf('4.0KV - inner ring calculations for d = %.3E', d));
d_inner(3) = d;

D = V4000.Omean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4000*em*c^2));
display(sprintf('4.0KV - outer ring calculations for d = %.3E', d));
d_outer(3) =d ;
display(' ')

%% Propagate error estimate
D = V4000.Imean + V4000.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_hi(3) = d_inner(3) - d;

D = V4000.Imean - V4000.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_lo(3) = d_inner(3) - d;

D = V4000.Omean + V4000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_hi(3) = d_outer(3) - d;

D = V4000.Omean - V4000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*4000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_lo(3) = d_outer(3) - d

%% V3500
V3500.Imean = mean(V3500.Inner);
V3500.Omean = mean(V3500.Outer);
V3500.ISE = std(V3500.Inner)./sqrt(5);
V3500.OSE = std(V3500.Outer)./sqrt(5);
display(sprintf('3.5KV - Inner Ring - mean: %.4Em, SE: %.4E',V3500.Imean, V3500.ISE));
display(sprintf('3.5KV - Outer Ring - mean: %.4Em, SE: %.4E',V3500.Imean, V3500.OSE));

D = V3500.Imean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3500*em*c^2));
display(sprintf('3.5KV - inner ring calculations for d = %.3E', d));
d_inner(4) = d;

D = V3500.Omean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3500*em*c^2));
display(sprintf('3.5KV - outer ring calculations for d = %.3E', d));
d_outer(4) =d ;

%% Propagate error estimate
D = V3500.Imean + V3500.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_hi(4) = d_inner(4) - d;

D = V3500.Imean - V3500.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_lo(4) = d_inner(4) - d;

D = V3500.Omean + V3500.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_hi(4) = d_outer(4) - d;

D = V3500.Omean - V3500.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3500*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_lo(4) = d_outer(4) - d;

display(' ')

%% V3000
V3000.Imean = mean(V3000.Inner);
V3000.Omean = mean(V3000.Outer);
V3000.ISE = std(V3000.Inner)./sqrt(5);
V3000.OSE = std(V3000.Outer)./sqrt(5);
display(sprintf('3.0KV - Inner Ring - mean: %.4Em, SE: %.4E',V3000.Imean, V3000.ISE));
display(sprintf('3.0KV - Outer Ring - mean: %.4Em, SE: %.4E',V3000.Imean, V3000.OSE));

D = V3000.Imean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3000*em*c^2));
display(sprintf('3.5KV - inner ring calculations for d = %.3E', d));
d_inner(5) = d;

D = V3000.Omean;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3000*em*c^2));
display(sprintf('3.0KV - outer ring calculations for d = %.3E', d));
d_outer(5) =d;

%% Propagate error estimate
D = V3000.Imean + V3000.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_hi(5) = d_inner(5) - d;

D = V3000.Imean - V3000.ISE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_inner_lo(5) = d_inner(5) - d;

D = V3000.Omean + V3000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_hi(5) = d_outer(5) - d;

D = V3000.Omean - V3000.OSE;
d = (4*pi*L*h*c)/(D*sqrt(2*e*3000*em*c^2));
display(sprintf('+error d = %.3E', d));
d_outer_lo(5) = d_outer(5) - d;
display(' ')


f(2) = figure(2);
clf(f(2));
colordef(f(2), 'black');
set(f(2), 'WindowStyle', 'docked');
V = [5000 4500 4000 3500 3000];
errorbar(V, d_inner.*10^-9, d_inner_hi.*10^-9, d_inner_lo.*10^-9);
hold on;
errorbar(V,d_outer.*10^-9, d_outer_hi.*10^-9, d_outer_lo.*10^-9,'r');
hold off;
title('Values for lattice');
xlabel('voltage');
ylabel('lattice spacing (nm)');
l =legend('inner ring', 'outer ring', -1);
set(l, 'TextColor', 'white');
grid on;

display(sprintf('Inner ring lattice value: %.3E', mean(d_inner)));
display(sprintf('Outer ring lattice value: %.3E', mean(d_outer)));

inner_accepted = .213 * 10^-9;
outer_accepted = .123 * 10^-9;

display(sprintf('Real Errors, Inner: %.2f, Outer: %.2f', mean(d_inner)/inner_accepted, mean(d_outer)/outer_accepted));


f(3) = figure(3);
clf(f(3));
colordef(f(3), 'black');
set(f(3), 'WindowStyle', 'docked');

DI = [V5000.Imean V4500.Imean V4000.Imean V3500.Imean V3000.Imean];
DO = [V5000.Omean V4500.Omean V4000.Omean V3500.Omean V3000.Omean];

plot(V, DI)
hold on;
plot(V, DO, 'r');
hold off;
grid on;
title('Voltage vs Diameter');
xlabel('Voltage');
ylabel('Ring Diameter');
l = legend('Inner', 'Outer', -1);
set(l, 'TextColor', 'white');