clc;
clear;
% Vsb FWHM of the Stark broadened line
% Vss Stark shift
%v frequency 
%v0 central frequency
wlmin=627.06;  % maximal wavelength in nm
wlmax=663.21;  % minimum wavelength in nm
wnmin=1e7/wlmax;  % minimum wavenumber in cm-1
wnmax=1e7/wlmin; % maximal wavenumber in cm-1
wlstep=(wlmax-wlmin)/625;  % CCD pixel 626*255
wl=wlmin:wlstep:wlmax;
spastep=6.5/255;
offset_zero=0.45;
spadistance=(6.5-offset_zero):-spastep:0;      % offset target place
lowwavelength=633;
highwavelength=635.8;
target1Low=round((lowwavelength-wlmin)/wlstep);
target1High=round((highwavelength-wlmin)/wlstep);
wlLeft=wlmin+wlstep*(target1Low-1);
wlRight=wlmin+wlstep*(target1High-1);
wlll=wlLeft:wlstep:wlRight;
[wlllcol wlllraw]=size(wlll);
lightspeed=3e8;
h=6.626e-34;
k=1.38e-23;
jjj=40;
centrawave1=634.711;
v1=3e8/centrawave1/1e-9;
SiI=degecyAndE('Si I.txt');
gI=1+2*SiI(:,1);
EnergyI=SiI(:,2);
SiII=degecyAndE('Si II.txt');
gII=1+2*SiII(:,1);
EnergyII=SiII(:,2);     % Si II ionization energy is 8.152eV
SiIII=degecyAndE('Si III.txt');
gIII=1+2*SiIII(:,1);
EnergyIII=SiIII(:,2);
SiIV=degecyAndE('Si IV.txt');
gIV=1+2*SiIV(:,1);
EnergyIV=SiIV(:,2);
SiV=degecyAndE('Si V.txt');
gV=1+2*SiV(:,1);
EnergyV=SiV(:,2);
SiVI=degecyAndE('Si VI.txt');
gVI=1+2*SiVI(:,1);
EnergyVI=SiVI(:,2);
ExperimentdataHV0=fitsread('SiC in vacuo, 532nm, 70ns, 645nm centre wavelength.fits');
offset_pixel=(6.5-offset_zero)/6.5*255;
expdataHV01=ExperimentdataHV0(1:offset_pixel,target1Low:target1High);     % offset
ypix=255*(1-(jjj*0.025+offset_zero)/6.5);
ypixel=round(ypix);
Experimentdata=expdataHV01(ypixel,:);     % High Voltage data
Nor1=(Experimentdata(:)-min(Experimentdata))';
NormExper=Nor1./max(Nor1);
a1=0.11111;
A21=0.584e8;
B21=3.14*(centrawave1*1e-9)^3/6.63e-34/3e8*A21;
g2_peak1=4;
g1_peak1=2;
g2_peak2=2;
g1_peak2=2;
B12=B21*4/2;
B12=B21*2/2;
Ne=9.167e22;
Te=0.5892;
Eup=10.073879+8.151683;
Eup2=10.066443+8.151683;
Elow=8.121023+8.151683;
Elow2=8.121023+8.151683;
Ti=0.5*28*1.67e-27*(80000)^2/1.6e-19;
v_left=3e17/wlLeft;
v_right=3e17/wlRight;
v_stepnum=target1High-target1Low+1;
v_step=(-v_right+v_left)/v_stepnum;

b=2;
% EEDF=Ka*sqrt(Energy)./(Te).^1.5.*exp(-Kb*(Energy/Te).^b)
Kb=(2*gamma(2.5/b)/3/gamma(1.5/b))^b;
Ka=b*Kb^(1.5/b)/gamma(1.5/b);
        % Si II
        J_IIintegral=@(v) v.^2.*1e12.*exp(-Kb.*((8.15+2.8.*v.^2).^b)./(Ti)^b);
        [J_II,err1II]=quadgk(J_IIintegral,0,100);
        ratioII0_I0=8*3.14*(9.11e-31/6.63e-34)^3/Ne*2/1*J_II;   % ground level population ratio of sigle charged species and neutral species   
        for i=2:size(gII)
            ratioII(i)=ratioII0_I0*gII(i)/gII(1)*exp(-Kb*(((EnergyII(i))^b-(EnergyII(1))^b)/Ti^b));
            % i-th level population ratio of sigle charged species and ground level
            % population neutral species
        end
        ratioII(1)=ratioII0_I0;
        totalSiII_I0=sum(ratioII);
        % Si III
        J_IIIintegral=@(v) v.^2.*1e12.*exp(-Kb.*((8.15+16.34+2*2.8.*v.^2).^b)./(Ti)^b);
        [J_III,err1III]=quadgk(J_IIIintegral,0,100);
        ratioIII0_I0=8*3.14*(9.11e-31/6.63e-34)^3/Ne*2/1*J_III;   % ground level population ratio of sigle charged species and neutral species   
        for i=2:size(gIII)
            ratioIII(i)=ratioIII0_I0*gIII(i)/gIII(1)*exp(-Kb*(((EnergyIII(i))^b-(EnergyIII(1))^b)/Ti^b));
            % i-th level population ratio of twice charged species and ground level
            % population neutral species
        end
        ratioIII(1)=ratioIII0_I0;
        totalSiIII_I0=sum(ratioIII);
        % si IV
        J_IVintegral=@(v) v.^2.*1e12.*exp(-Kb.*((8.15+16.34+33.49+3*2.8.*v.^2).^b)./(Ti)^b);
        [J_IV,err1IV]=quadgk(J_IVintegral,0,100);
        ratioIV0_I0=8*3.14*(9.11e-31/6.63e-34)^3/Ne*2/1*J_IV;   % ground level population ratio of sigle charged species and neutral species   
        for i=2:size(gIV)
            ratioIV(i)=ratioIV0_I0*gIV(i)/gIV(1)*exp(-Kb*(((EnergyIV(i))^b-(EnergyIV(1))^b)/Ti^b));
            % i-th level population ratio of triple charged species and ground level
            % population neutral species
        end
        ratioIV(1)=ratioIV0_I0;
        totalSiIV_I0=sum(ratioIV);
        % si V
        J_Vintegral=@(v) v.^2.*1e12.*exp(-Kb.*((8.15+16.34+33.49+45.14+3*2.8.*v.^2).^b)./(Ti)^b);
        [J_V,err1V]=quadgk(J_Vintegral,0,100);
        ratioV0_I0=8*3.14*(9.11e-31/6.63e-34)^3/Ne*2/1*J_V;   % ground level population ratio of sigle charged species and neutral species   
        for i=2:size(gV)
            ratioV(i)=ratioV0_I0*gV(i)/gV(1)*exp(-Kb*(((EnergyV(i))^b-(EnergyV(1))^b)/Ti^b));
            % i-th level population ratio of four charged species and ground level
            % population neutral species
        end
        ratioV(1)=ratioV0_I0;
        totalSiV_I0=sum(ratioV);
        NsiI0=Ne/(ratioII0_I0+2*totalSiIII_I0+3*totalSiIV_I0+4*totalSiV_I0);
        NsiII0=NsiI0*ratioII0_I0;
        NsiIII0=NsiI0*totalSiIII_I0;
        NsiIV0=NsiI0*totalSiIV_I0;
        NsiV0=NsiI0*totalSiV_I0;
        Nd=(8.9e-12*Te*1.6e-19)^1.5/(1.6e-19)^3/(Ne)^0.5;      
        for iwd=1:501
            for jwd=1:501
                w1=0.1e-9+(iwd-1)*0.001e-9;
                d1=40e-12+(jwd-1)*0.1e-12;
                v_width1=3e8/(centrawave1*1e-9)^2*sqrt((2*w1*Ne/1e23*(1+1.75*a1*(Ne/1e23)^0.25*(1-0.75*Nd^(-1/3)))/0.6)^2+(0.062e-9)^2);
                v_shift1=-2*d1*3e8/(centrawave1*1e-9)^2*Ne/1e23*(1-2*a1*w1/d1*(Ne/1e23)^0.25*(1-0.75*(Nd)^(-1/3))/0.6);
                error1=0.0045e-9; 
                N2_peak1=5*exp(-Kb*(((Eup)^b-(8.15+16.35)^b)/Ti^b))*NsiI0*ratioIII0_I0;
                N1=3*exp(-Kb*(((Elow)^b-(8.15+16.35)^b)/Ti^b))*NsiI0*ratioIII0_I0;
                for i=1:v_stepnum
                    v=v_left-(i-1)*v_step;
                    Y1(i)=1/3.14/2*v_width1./(v_width1^2/4+(v-v_shift1-v1)^2);            
                    Kk=0.28*6.63e-34*v*(3e9);
                    Iemission(i)=Kk*A21*(N2_peak1*Y1(i))-(6.63e-34*v/3/1e8)*(B21*(N2_peak1*Y1(i))-N1*(B12*Y1(i)));
                end
                NormIemi=Iemission./max(Iemission);
                [mmmm nnnn]=rsquare(NormExper,NormIemi);
                r2insta(iwd,jwd)=mmmm;
            end
        end
        
        [rw cd]=find(r2insta==max(max(r2insta(:,:))));
        w1=0.1e-9+(rw-1)*0.001e-9
        d1=40e-12+(cd-1)*0.1e-12
