
lightspeed=3e8;
h=6.626e-34;
k=1.38e-23;

v0=lightspeed/cenwavelenth; 
wlmax=489.45;  % maximal wavelength in nm
wlmin=450.41;  % 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;  % wavelength axis from 450nm to 489nm

A21=1.25e8;
% W, D electron impact parameters for stark width and stark shift
% A ion correction parameter

SiI=degecyAndE('E:\new code\degenercy\Si I.txt');
gI=1+2*SiI(:,1);
EnergyI=SiI(:,2);

SiII=degecyAndE('E:\new code\degenercy\Si II.txt');
gII=1+2*SiII(:,1);
EnergyII=SiII(:,2);     % Si II ionization energy is 8.152eV

SiIII=degecyAndE('E:\new code\degenercy\Si III.txt');
gIII=1+2*SiIII(:,1);
EnergyIII=SiIII(:,2);

SiIV=degecyAndE('E:\new code\degenercy\Si IV.txt');
gIV=1+2*SiIV(:,1);
EnergyIV=SiIV(:,2);

SiV=degecyAndE('E:\new code\degenercy\Si V.txt');
gV=1+2*SiV(:,1);
EnergyV=SiV(:,2);

SiVI=degecyAndE('E:\new code\degenercy\Si VI.txt');
gVI=1+2*SiVI(:,1);
EnergyVI=SiVI(:,2);



Experimentdata0=fitsread('E:\work\phd\Silicon_PLA\experiment data\hao_Si_Starkbroadening_1200lines\1200lines_40ns_240us.fits');
expdata=Experimentdata0(1:233,:);     % offset

zzz=0.4;  %spectrum position at distance axis zzz
pixesposition=round((1-zzz/5.94)*233);

wlLeftPix=round((cenwavelenth-wlmin)/wlstep)-65;
wlRightPix=round((cenwavelenth-wlmin)/wlstep)+65;
wlLeft=wlmin+wlLeftPix*wlstep;
wlRight=wlmin+wlRightPix*wlstep;
wlll=wlLeft:wlstep:wlRight;       %wavelength axis of Si III centraled at cenwavelength

Experimentdata00=Experimentdata0(pixesposition,:);
Experimentdata=Experimentdata00(wlLeftPix:wlRightPix);
Nor1=Experimentdata(:)-min(Experimentdata);
NormExper=Nor1./max(Nor1);

%Ne=1.25e22;
%Te=0.79552*11600;   d=1.8 mm


Ne=0.69e22;
Ti=0.5*28*1.67e-27*(2.6e6/40)^2/1.6e-19;
Te=0.46;
% d=2.6mm

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;
v_width=3e8/(456.78e-9)^2*2*w*Ne/1e22*(1+1.75*a*(Ne/1e22)^0.25*(1-0.75*Nd^(-1/3)));
v_shift=d*3e8/(456.78e-9)^2*Ne/1e22*(1+2*a*w/d*(Ne/1e22)^0.25*(1-0.75*(Nd)^(-1/3)));
v1=3e8/455.2622e-9;
v2=3e8/456.784e-9;
v3=3e8/457.4757e-9;


B21=3.14*(456.78e-9)^3/6.63e-34/3e8*A21;
g2_4552=5;
g1_4552=3;
g2_4568=3;
g1_4568=3;
g2_4575=1;
g1_4575=3;
B12_4552=B21*5/3;
B12_4568=B21*1;
B12_4575=B21*1/3;
Eup=21.739+8.15+16.35;
Elow=19.016+8.15+16.35;

N2_4552=5*exp(-Kb*(((Eup)^b-(8.15+16.35)^b)/Ti^b))*NsiI0*ratioIII0_I0;
N2_4568=3*exp(-Kb*(((Eup)^b-(8.15+16.35)^b)/Ti^b))*NsiI0*ratioIII0_I0;
N2_4575=1*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;

v_left=3e17/wlLeft;
v_right=3e17/wlRight;
v_stepnum=wlRightPix-wlLeftPix+1;
v_step=(-v_right+v_left)/v_stepnum;
for i=1:v_stepnum
    v=v_left-(i+1)*v_step;
    Y1=1/3.14/2*v_width./(v_width^2/4+(v+v_shift-v1)^2);
    Y2=1/3.14/2*v_width./(v_width^2/4+(v+v_shift-v2)^2);
    Y3=1/3.14/2*v_width./(v_width^2/4+(v+v_shift-v3)^2);
    Kk=1/4/3.14*6.63e-34*v;
    Iemission(i)=60*3e8*Kk*A21*(N2_4552*Y1+N2_4568*Y2+N2_4575*Y3)./(1-Kk*B21*(N2_4552*Y1+N2_4568*Y2+N2_4575*Y3)+Kk*N1*(B12_4552*Y1+B12_4568*Y2+B12_4575*Y3));
end
NorIem=Iemission./max(Iemission);
figure (1)
plot(wlll,NorIem,'r')
hold on
plot(wlll,NormExper)

