Physics307L:People/Osinski/Photoelectric/Matlab Code

%% Charge Time
% raw data
intens=[20 40 60 80 100];
tYellow20=[8.25 8.75 9.20 9.06 9.28];
tYellow40=[8.00 8.42 6.72 8.75 7.59];
tYellow60=[6.66 5.25 6.00 5.06 5.75];
tYellow80=[4.37 5.19 4.13 4.25 4.03];
tYellow100=[3.69 3.50 4.28 4.50 4.19];
tViolet20=[16.50 19.06 19.30 17.12 16.80];
tViolet40=[8.72 9.68 9.68 9.75 8.84];
tViolet60=[7.06 7.66 7.31 6.97 7.19];
tViolet80=[4.41 5.44 5.22 5.25 5.46];
tViolet100=[4.88 4.56 4.44 4.62 4.65];

% mean for each intensity
MeanY20=mean(tYellow20);
MeanY40=mean(tYellow40);
MeanY60=mean(tYellow60);
MeanY80=mean(tYellow80);
MeanY100=mean(tYellow100);
MeanV20=mean(tViolet20);
MeanV40=mean(tViolet40);
MeanV60=mean(tViolet60);
MeanV80=mean(tViolet80);
MeanV100=mean(tViolet100);
MeanY=[MeanY20 MeanY40 MeanY60 MeanY80 MeanY100];
MeanV=[MeanV20 MeanV40 MeanV60 MeanV80 MeanV100];

% standard deviation of the mean for each intensity
stdmY20=sqrt(sum((tYellow20-MeanY20).^2)/20);
stdmY40=sqrt(sum((tYellow40-MeanY40).^2)/20);
stdmY60=sqrt(sum((tYellow60-MeanY60).^2)/20);
stdmY80=sqrt(sum((tYellow80-MeanY80).^2)/20);
stdmY100=sqrt(sum((tYellow100-MeanY100).^2)/20);
stdmV20=sqrt(sum((tViolet20-MeanV20).^2)/20);
stdmV40=sqrt(sum((tViolet40-MeanV40).^2)/20);
stdmV60=sqrt(sum((tViolet60-MeanV60).^2)/20);
stdmV80=sqrt(sum((tViolet80-MeanV80).^2)/20);
stdmV100=sqrt(sum((tViolet100-MeanV100).^2)/20);

% total stdm
for f=20:20:100
    stringy=sprintf('sum(stdmY%d + stdmV%d)',f,f)
    eval(stringy)
end
stdmfull=(0.7785+0.5803+0.4045+0.3997+0.2599)/10
stdmbetter=(0.1877+0.3529+0.2843+0.2070+0.1878+0.2273+ 0.1202+0.1927+0.0721)/9

% error bar plots
figure(1)
title('charge time vs. intensity for yellow and violet spectra')
xlabel('intensity (percentage)')
ylabel('charge time (s)')
for f=20:20:100
    hold on
    plot(intens,MeanY,'Color','yellow')
    plot(intens,MeanV,'Color','magenta')
    Ystring=sprintf('errorbar(%d,MeanY%d,stdmY%d,''*'')',f,f,f)
    eval(Ystring)
    Vstring=sprintf('errorbar(%d,MeanV%d,stdmV%d,''*'',''Color'',''red'')',f,f,f)
    eval(Vstring)
end    
legend('yellow','violet')
hold off

%% Stop Voltage vs. Frequency
MaxV=[2.044 1.71 1.49 .847 .716];
Color=[8.20264E14 7.40858E14 6.87858E14 5.48896E14 5.18672E14];
polV=polyfit(Color, MaxV, 1);
funV=polyval(polV,Color);
figure(2)
hold on
plot(Color,MaxV,'o'); plot(Color, funV,'color','green')
title('Stop Voltage vs. Frequency')
xlabel('frequency(Hz)')
ylabel('stopping voltage(V)')

%% Determination of h
% raw data
data1 = [2.064, 1.718, 1.493, 0.849, 0.717];
data2 = [2.062, 1.731, 1.520, 0.842, 0.735];
data3 = [2.069, 1.725, 1.495, 0.850, 0.714];
data4 = [2.061, 1.725, 1.515, 0.849, 0.735];
freq = [8.20264E14, 7.40858E14, 6.87858E14, 5.48896E14, 5.18672E14];

ev = 1.602E-19; %electron charge.
data1 = data1 .* ev;
data2 = data2 .* ev;
data3 = data3 .* ev;
data4 = data4 .* ev;

% least squares fitting
poly1=polyfit(freq, data1, 1)
fun1=polyval(poly1,freq);
poly2=polyfit(freq, data2, 1)
fun2=polyval(poly2,freq);
poly3=polyfit(freq, data3, 1)
fun3=polyval(poly3,freq);
poly4=polyfit(freq, data4, 1)
fun4=polyval(poly4,freq);

% plots
figure(3);
subplot(2,2,1); scatter(freq, data1); hold on, plot(freq, fun1,'color','green');...
    title('Run1, 1st order'),xlabel('fq(Hz)'),ylabel('V'), hold off
subplot(2,2,2); scatter(freq, data2); hold on, plot(freq, fun2,'color','green');...
    title('Run1, 2nd order'),xlabel('fq(Hz)'),ylabel('V'), hold off
subplot(2,2,3); scatter(freq, data3); hold on, plot(freq, fun3,'color','green');...
    title('Run2, 1st order'),xlabel('fq(Hz)'),ylabel('V'), hold off
subplot(2,2,4); scatter(freq, data4); hold on, plot(freq, fun4,'color','green');...
    title('Run2, 2nd order'),xlabel('fq(Hz)'),ylabel('V'), hold off

%calculating h
h1=(fun1(end)-fun1(1))/(freq(end)-freq(1))
h2=(fun2(end)-fun2(1))/(freq(end)-freq(1))
h3=(fun3(end)-fun3(1))/(freq(end)-freq(1))
h4=(fun4(end)-fun4(1))/(freq(end)-freq(1))

%calculating omega in J
omeg1J=-interp1(freq,fun1,0,'pchip')
omeg2J=-interp1(freq,fun2,0,'pchip')
omeg3J=-interp1(freq,fun3,0,'pchip')
omeg4J=-interp1(freq,fun4,0,'pchip')

%calculating omega in eV
omeg1eV=-interp1(freq,fun1,0,'pchip')/ev
omeg2eV=-interp1(freq,fun2,0,'pchip')/ev
omeg3eV=-interp1(freq,fun3,0,'pchip')/ev
omeg4eV=-interp1(freq,fun4,0,'pchip')/ev

%% Error Analysis

% error for eV from least squares fits coefficients h and omeg
sigmay1=sqrt((1/15)*sum((data1-(h1*freq)+omeg1J).^2))
sigmay2=sqrt((1/15)*sum((data2-(h2*freq)+omeg2J).^2))
sigmay3=sqrt((1/15)*sum((data3-(h3*freq)+omeg3J).^2))
sigmay4=sqrt((1/15)*sum((data4-(h4*freq)+omeg4J).^2))

Delta=(5*sum(freq.^2))-(sum(freq)^2)

% calculating stdm for h and omega
sigmah1=sigmay1*sqrt(5/Delta)
sigmah2=sigmay2*sqrt(5/Delta)
sigmah3=sigmay3*sqrt(5/Delta)
sigmah4=sigmay4*sqrt(5/Delta)
sigmaoeV1=sigmay1*sqrt(sum(freq.^2)/Delta)/ev
sigmaoeV2=sigmay2*sqrt(sum(freq.^2)/Delta)/ev
sigmaoeV3=sigmay3*sqrt(sum(freq.^2)/Delta)/ev
sigmaoeV4=sigmay4*sqrt(sum(freq.^2)/Delta)/ev

% average h and omega with their average stdms
Havg=(h1+h2+h3+h4)/4
STDMH=(sigmah1+sigmah2+sigmah3+sigmah4)/4
OMEGavg=(omeg1eV+omeg2eV+omeg3eV+omeg4eV)/4
STDMOMEG=(sigmaoeV1+sigmaoeV2+sigmaoeV3+sigmaoeV4)/4

%% Extra Formal Report Measurements
freq = [5.18672E14, 5.48896E14, 6.87858E14, 7.40858E14, 8.20264E14];
newdata1=[.713 .845 1.488 1.669 2.009]
newdata2=[.710 .841 1.480 1.705 2.045]
newdata20=[.707 .840 1.480 1.696 2.035]
newdata60=[.707 .842 1.481 1.704 2.035]

ev = 1.602E-19; %electron charge.
newdata1 = newdata1 .* ev;
newdata2 = newdata2 .* ev;
newdata20 = newdata20 .* ev;
newdata60 = newdata60 .* ev;

% least squares fitting
newpoly1=polyfit(freq, newdata1, 1)
newfun1=polyval(newpoly1,freq);
newpoly2=polyfit(freq, newdata2, 1)
newfun2=polyval(newpoly2,freq);
newpoly20=polyfit(freq, newdata20, 1)
newfun20=polyval(newpoly20,freq);
newpoly60=polyfit(freq, newdata60, 1)
newfun60=polyval(newpoly60,freq);

% plots
figure(4);
subplot(2,2,1); scatter(freq, newdata1); hold on, plot(freq, newfun1,'color','red');...
    title('Run1, 100% intensity'),xlabel('fq(Hz)'),ylabel('V'), hold off
subplot(2,2,2); scatter(freq, newdata2); hold on, plot(freq, newfun2,'color','red');...
    title('Run2, 100% intensity'),xlabel('fq(Hz)'),ylabel('V'), hold off
subplot(2,2,3); scatter(freq, newdata20); hold on, plot(freq, newfun20,'color','red');...
    title('Run3, 20% intensity'),xlabel('fq(Hz)'),ylabel('V'), hold off
subplot(2,2,4); scatter(freq, newdata60); hold on, plot(freq, newfun60,'color','red');...
    title('Run4, 60% intensity'),xlabel('fq(Hz)'),ylabel('V'), hold off

% calculating h
h1f=(newfun1(end)-newfun1(1))/(freq(end)-freq(1))
h2f=(newfun2(end)-newfun2(1))/(freq(end)-freq(1))
h20f=(newfun20(end)-newfun20(1))/(freq(end)-freq(1))
h60f=(newfun60(end)-newfun60(1))/(freq(end)-freq(1))

% calculating w
omeg1Jf=-interp1(freq,newfun1,0,'pchip');
omeg2Jf=-interp1(freq,newfun2,0,'pchip');
omeg20Jf=-interp1(freq,newfun20,0,'pchip');
omeg60Jf=-interp1(freq,newfun60,0,'pchip');

% error for eV from least squares fits coefficients h and omeg
sigmay1f=sqrt((1/15)*sum((newdata1-(h1f*freq)+omeg1Jf).^2));
sigmay2f=sqrt((1/15)*sum((newdata2-(h2f*freq)+omeg2Jf).^2));
sigmay20f=sqrt((1/15)*sum((newdata20-(h20f*freq)+omeg20Jf).^2));
sigmay60f=sqrt((1/15)*sum((newdata60-(h60f*freq)+omeg60Jf).^2));

Delta=(5*sum(freq.^2))-(sum(freq)^2);

% calculating stdm for h and omega
sigmah1f=sigmay1f*sqrt(5/Delta)
sigmah2f=sigmay2f*sqrt(5/Delta)
sigmah20f=sigmay20f*sqrt(5/Delta)
sigmah60f=sigmay60f*sqrt(5/Delta)

STDMHf=(sigmah1f+sigmah2f+sigmah20f+sigmah60f)/4