B1=1.5039759;                     
B2=0.55069141;
B3=6.5927379;
C1=0.0054804113;
C2=0.014799428;
C3=402.89514;
c=3*10^8;
x=0.4:0.01:2.0;
n=(1+(B1.*x.^2)./(x.^2-C1)+(B2.*x.^2)./(x.^2-C2)+(B3.*x.^2)./(x.^2-C3)).^0.5;
T1=gradient(n);                                                                    %dn/dx
T2=gradient(T1);                                                                   %d2n/dx2
T3=gradient(T2);                                                                   %d3n/dx3
D1=-x.^2/(2*pi.*c).*T1;                                                            %dn/dw
D2=2*x.^3/(4*pi.*pi.*c.*c).*T1+x.^4/(4*pi.*pi.*c.*c).*T2;                          %d2n/dw2
D3=-x.^2.*(4.*x.^2.*T1+5*x.^3.*T2+x.^4.*T3)./(8*pi.^3.*c.^3);                      %d3n/dw3
b3=(3*D2+(2*pi.*c)./x.*D3)./c;
plot(x,b3.*10^30,'b')
hold on;

B1=1.4313493;
B2=0.65054713;
B3=5.3414021;
C1=0.0052799261;
C2=0.014238265;
C3=325.01783;
n=(1+(B1.*x.^2)./(x.^2-C1)+(B2.*x.^2)./(x.^2-C2)+(B3.*x.^2)./(x.^2-C3)).^0.5;
T1=gradient(n);                                                                    %dn/dx
T2=gradient(T1);                                                                   %d2n/dx2
T3=gradient(T2);                                                                   %d3n/dx3
D1=-x.^2/(2*pi.*c).*T1;                                                            %dn/dw
D2=2*x.^3/(4*pi.*pi.*c.*c).*T1+x.^4/(4*pi.*pi.*c.*c).*T2;                          %d2n/dw2
D3=-x.^2.*(4.*x.^2.*T1+5*x.^3.*T2+x.^4.*T3)./(8*pi.^3.*c.^3);                      %d3n/dw3
b3=(3*D2+(2*pi.*c)./x.*D3)./c;
plot(x,b3.*10^30,'r')
hold on;

B1=0.6961663;
B2=0.4079426;
B3=0.8974794;
C1=0.0046791483;
C2=0.013512063;
C3=97.934003;
n=(1+(B1.*x.^2)./(x.^2-C1)+(B2.*x.^2)./(x.^2-C2)+(B3.*x.^2)./(x.^2-C3)).^0.5;
T1=gradient(n);                                                                    %dn/dx
T2=gradient(T1);                                                                   %d2n/dx2
T3=gradient(T2);                                                                   %d3n/dx3
D1=-x.^2/(2*pi.*c).*T1;                                                            %dn/dw
D2=2*x.^3/(4*pi.*pi.*c.*c).*T1+x.^4/(4*pi.*pi.*c.*c).*T2;                          %d2n/dw2
D3=-x.^2.*(4.*x.^2.*T1+5*x.^3.*T2+x.^4.*T3)./(8*pi.^3.*c.^3);                      %d3n/dw3
b3=(3*D2+(2*pi.*c)./x.*D3)./c;
plot(x,b3.*10^30,'g')
hold on;

B1=1.0396121;
B2=0.23179234;
B3=1.0104695;
C1=0.0060006987;
C2=0.020017914;
C3=103.56065;
n=(1+(B1.*x.^2)./(x.^2-C1)+(B2.*x.^2)./(x.^2-C2)+(B3.*x.^2)./(x.^2-C3)).^0.5;
T1=gradient(n);                                                                    %dn/dx
T2=gradient(T1);                                                                   %d2n/dx2
T3=gradient(T2);                                                                   %d3n/dx3
D1=-x.^2/(2*pi.*c).*T1;                                                            %dn/dw
D2=2*x.^3/(4*pi.*pi.*c.*c).*T1+x.^4/(4*pi.*pi.*c.*c).*T2;                          %d2n/dw2
D3=-x.^2.*(4.*x.^2.*T1+5*x.^3.*T2+x.^4.*T3)./(8*pi.^3.*c.^3);                      %d3n/dw3
b3=(3*D2+(2*pi.*c)./x.*D3)./c;
plot(x,b3.*10^30,'m')
hold on;

B1=1.6162598;
B2=0.25922933;
B3=1.1749087;
C1=0.013606860;
C2=0.061596046;
C3=121.92271;
n=(1+(B1.*x.^2)./(x.^2-C1)+(B2.*x.^2)./(x.^2-C2)+(B3.*x.^2)./(x.^2-C3)).^0.5;
T1=gradient(n);                                                                    %dn/dx
T2=gradient(T1);                                                                   %d2n/dx2
T3=gradient(T2);                                                                   %d3n/dx3
D1=-x.^2/(2*pi.*c).*T1;                                                            %dn/dw
D2=2*x.^3/(4*pi.*pi.*c.*c).*T1+x.^4/(4*pi.*pi.*c.*c).*T2;                          %d2n/dw2
D3=-x.^2.*(4.*x.^2.*T1+5*x.^3.*T2+x.^4.*T3)./(8*pi.^3.*c.^3);                      %d3n/dw3
b3=(3*D2+(2*pi.*c)./x.*D3)./c;
plot(x,b3.*10^30,'c')
title('三階色散曲線');
xlabel('波長');
ylabel('三階色散');
legend('Sapphire(n_o)','Sapphire(n_e)','Fuse silica','BK7','SF10')
hold on;