%%fiber dispersion calculation 24/11/2014 made by Masaki Tokurakawa UEC
%%less approxomation
clear all
c=3*10^8; %% velocity of light
a=4.0*10^(-6); %%fiber core 
syms x; %%symbol for a function of refractive index n(x) 
syms vx; %%symbol for a function of b parameter b(V)
%%%core sellmeir's equation SiO2 93.7% GeO2 6.3% Trans IECE Japan65-E 166
A0=0.7083952;
A1=0.4203993;
A2=0.8663412;
A3=sqrt(0.007290464);
A4=sqrt(0.01050294);
A5=sqrt(97.93428);
n1=symfun(sqrt(1+A0/(1-(A3/x)^2)+A1/(1-(A4/x)^2)+A2/(1-(A5/x)^2)),[x])-0.004;%%to reduce NA, -0.00x
GV1=diff(n1);
GVD1=diff(n1,2);
Nm=0;%%mode number
%%%clad sellmeir's equation SiO2 100% JOSA 55 1205
A0=0.6961663;
A1=0.4079426;
A2=0.8974794;
A3=0.0684043;
A4=0.1162414;
A5=9.896161;
n2=symfun(sqrt(1+A0/(1-(A3/x)^2)+A1/(1-(A4/x)^2)+A2/(1-(A5/x)^2)),[x]);
GV2=diff(n2);
GVD2=diff(n2,2);
%%%V para peter & b
x = 1:0.01:2.5; %%wavelength range
V=double(subs(2.*pi.*a.*sqrt((n1.^2-n2.^2))));
V=V./x*10^6;

b = 0.1:0.0001:0.80; %% range for b
V1=1./sqrt(1-b).*(atan(sqrt(b./(1-b)))+(pi+2*Nm)/4); %% V as a functon of b
VC = polyfit(V1,b,8); %% fitting by 8 order, be careful if there is huge error!! 
bV=symfun(VC(9)+VC(8).*vx+VC(7).*vx.^2+VC(6).*vx.^3+VC(5).*...
    vx.^4+VC(4).*vx.^5+VC(3).*vx.^6+VC(2).*vx.^7+VC(1).*vx.^8,[vx]);
%%% b is expressed as a function of V
vx=V; %% insert V to vx at corrresponding wavelength
difV=diff(bV);              %%db/dV
difVV=diff(bV,2);           %%d2b/dV2
X1=double(subs(n1,x));      %%table of core refractive index
X2=double(subs(n2,x));      %%table of clad refractive index
X3=double(subs(bV,vx));     %%table of b at X(x)
X4=double(subs(GVD1,x));    %%d2n/dx2 core
X5=double(subs(GVD2,x));    %%d2n/dx2 clad
X6=double(subs(difV,V));    %% table of db/dV
X7=double(subs(difVV,V));   %% table of d2b/dV2
X8=double(subs(GV1,x));     %%dn/dx core
X9=double(subs(GV2,x));     %%dn/dx clad
X10 = double(subs((-1./c./x.*(x.^2.*X4+(x.^2.*X4-x.^2.*X5).*(X3+1/2.*vx.*X6))))*10^12);
NA=sqrt(X1.^2-X2.^2);
%%%stracture dispersion 
beta=2*pi./x.*sqrt(X1.^2.*X3+X2.^2.*(1-X3));%% propagation constant
neq=beta.*x/2/pi;
N1=X1-x.*X8;
N2=X2-x.*X9;
X11=-1./c./x./2.*(X1.*N1-X2.*N2).^2./neq./(X1.^2-X2.^2).*(vx.^2.*X7+2*vx.*X6)*10^12;

figure(1)
subplot(2,2,1);
plot(x,X1,x,X2)
axis ([1 2.5 1.44 1.47 ])
legend('n1','n2',1)
xlabel('wavelength','Color','b','fontsize',18);
ylabel('refractive index','fontsize',18);

subplot(2,2,2);
plot(x,V,'r',x,X3,'b')
axis ([1 2.5 -inf inf ])
legend('V','b',1)
xlabel('wavelength','Color','b','fontsize',18);
ylabel('V & b','fontsize',18);

subplot(2,2,3);
plot(V1,b,V,X3)
axis ([0 5 0 1])
legend('V(b)','b(V) fitted',2)
set(legend,'FontSize',18);
set(gca,'FontSize',14);
xlabel('V','Color','b','fontsize',18);
ylabel('b','Color','b','fontsize',18);

subplot(2,2,4);
plot(x,NA)
axis ([1 2.2 0.11 inf])
legend('NA',2)
set(legend,'FontSize',18);
set(gca,'FontSize',14);
xlabel('wavelength(um)','Color','b','fontsize',18);
ylabel('NA','Color','b','fontsize',18);

figure(2)
plot(x,X10,'r',x,X11,'b',x,X10+X11,'g',x,0);
axis ([1 2.2 -30 inf ])
legend('material','structure','total',2)
set(legend,'FontSize',18);
set(gca,'FontSize',14);
xlabel('wavelength(um)','Color','b','fontsize',18);
ylabel('dispersion ps/nm/km','Color','b','fontsize',18);

%%soliton pulse duration
dis_beta=-(X10+X11).*x.^2./2/pi/c*10^(-18);
Pout=20*10^(-3);
rr=20*10^(6);
Loss=0.1;
Ep=Pout/rr/Loss;
w0 = a.*(0.65+1.619.*vx-1.5+2.879.*vx.^(-6));
Area=w0.^2*pi;
Nnonre=2.1*10^(-20);

ganma=Nnonre*2*pi./x/Area*10^6;
t=3.11.*dis_beta./ganma/Ep*10^12;

figure(3)
plot(x,dis_beta,'r');
axis ([1 2.2 -inf inf ])
legend('second orderdispersion',2)
set(legend,'FontSize',18);
set(gca,'FontSize',14);
xlabel('wavelength(um)','Color','b','fontsize',18);
ylabel('dispersion s2/m','Color','b','fontsize',18);

figure(4)
plot(x,abs(t),'r');
axis ([1 2.2 -inf inf ])
legend('soliton',2)
set(legend,'FontSize',18);
set(gca,'FontSize',14);
xlabel('wavelength(um)','Color','b','fontsize',18);
ylabel('pulse duration (ps)','Color','b','fontsize',18);


% figure(5)
% plot(x,vx.^2.*X7+vx.*X6,'r',x,0,x,N2,'b');
% axis ([1 1.7 -inf inf ])