%Phase Matching

a=load('SHGPhseMatchingBetterData2.txt');
% for i=1:length(a)
%     if a(i+1, 3)>a(i, 3)
%         

b=a(:, 3:6);
time=[1:length(b)]/60;
figure(1)
[AX,H1,H2]=plotyy(time, b(:,2), time, b(:,4)*570/409);
set(get(AX(1),'Ylabel'),'String','Temperature [C] rough calibration') 
set(get(AX(2),'Ylabel'),'String','Green Power (uW) rough calibration') 
xlabel('time(min)')
title('SHG phase matching')
%stability
c=load('SHGstability2.txt');
stab=c(:,5:6);

for i=1:(length(stab)-1)
    if stab(i+1,1)>0.5
        if stab(i,1)<0.5
            onRow=i+1;
        end
    end
    if stab(i+1, 1)<0.5
        if stab(i,1)>0.5
            offRow=i+1;
        end
    end
end
%cut out rows while it is transitioning
buffer=1;
on1=onRow+buffer;
off1=offRow-buffer;
on2=onRow-buffer;
off2=offRow+buffer;
lengthStab=length(stab);
on=stab(on1:off1, :);%vector of powers measured while laser was on
off=[stab(1:on2, :); stab(off2:lengthStab, :)]; %vecotr of powers measured while laser was off
convEff=on(:,1)./on(:,2);
convEffP=convEff./mean(convEff);
EffNL=on(:,1)./(on(:, 2).^2);
EffNLP=EffNL/mean(EffNL);

onP=[on(:,1)./mean(on(:,1)),on(:,2)./mean(on(:,2))]; 
offP=[off(:,1)./mean(off(:,1)),off(:,2)./mean(off(:,2))]; 

stdonP=std(onP)
stdon=std(on)
stdoff=std(off)
stdCEff=std(convEffP)

timestep=6.463057-3.464479;
time=([1:length(on)])/(60^2);
figure(2)
plot(time, EffNLP, 'k')
hold on
plot(time, convEffP, 'b')
plot(time, onP(:,1), 'g')
plot(time, onP(:,2), 'r')
legend('Effective nonlinearity [532 power/(1064 power)^2] normalized', 'Conversion efficiceny normalized', 'Green power normalized', 'IR power normalized')
hold off
xlabel('Time (hours)')
ylabel('Normalized')

bins=0.8:0.0005:1.1;
[cEHist, bout]=hist(convEffP, bins);
[onHist, bout2]=hist(onP, bins);
GreenHist=onHist(:,1);
IRHist=onHist(:,2);
cEfit=fit(bout', cEHist', 'gauss1')
GreenFit=fit(bout', GreenHist, 'gauss1')
IRFit=fit(bout', IRHist, 'gauss1')

figure(3)
plot(bout, cEHist, 'b.', bout2, onHist(:,1), 'g.', bout2, onHist(:,2), 'r.')
xlim([0.95, 1.05])
hold on
plot(cEfit, 'b')
plot(GreenFit, 'g')
plot(IRFit, 'r')
legend('Conversion Efficiency', '532 Power', '1064 Power', ...
    strcat(['Conversion efficency fit \sigma = ', num2str(100*cEfit.c1, 3), '%']),...
    strcat(['532 power fit \sigma = ', num2str(100*GreenFit.c1, 3), '%']),...
    strcat([' 1064 power fit \sigma = ', num2str(100*IRFit.c1, 3), '%']))
% 

%caclulate exptected conversion efficiency, boyd pg 121
K=1.068; %gaussian units
clight=299792458*100;
omega1=2*pi*clight/1064e-7;
deff=(13.8e-12/299.79)/100;
L=1; %crystal is 10 mm, I hope this is true it is what label said on box
n1=1.8;%1.74;
n2=1.8%1.78;
Pomega=0.24e7;
% in mW
P532=K*Pomega^2*128*pi^2*omega1^3*deff^2*L/(clight^4*n1*n2)*10^7*10^3

%somehow, this must be wrong
