% This program plots the variation of refractive index against lamda
% Taken from Malitson, JOSA 1965. 
clear;clc;close all;
inc = 5e-3;lamarr = 0.5:inc:2;
n = sqrt(1 + lamarr.^2.*(0.6961663./(lamarr.^2 - 0.0684043^2) + 0.4079426./(lamarr.^2 - 0.1162414^2) + 0.8974794./(lamarr.^2 - 9.896161^2)));
ng = n(1:end-1) - lamarr(1:end-1).*diff(n)/inc;ng = [ng ng(end)];
ngstr = '\itn\rm\bf_2\it_g';n2str = '\itn\rm\bf_2';
figure(1)
plot(lamarr,n,'-k ','LineWidth',2);hold on
plot(lamarr,ng,'--r ','LineWidth',2);hold off
set(gcf,'Color','w');set(gca,'FontSize',16);
ylabel('refractive index','FontSize',18,'FontWeight','bold');
xlabel('wavelength in \mum','FontSize',18,'FontWeight','bold')
text(min(lamarr)*1.2,max(n)*1,ngstr,'FontSize',22,'FontWeight','bold','Edgecolor','white','BackgroundColor','white');
text(min(lamarr)*1.2,max(n)*1,n2str,'FontSize',22,'FontWeight','bold','Edgecolor','white','BackgroundColor','white');
dng = diff(ng(1:end-1))/inc;
figure(2)
plot(lamarr(1:end-2),dng,'-k ','LineWidth',2);
set(gcf,'Color','w');set(gca,'FontSize',16);
ylabel('\itdn\rm\bf_2\it_g / d\lambda','FontSize',18,'FontWeight','bold');
xlabel('wavelength in \mum','FontSize',18,'FontWeight','bold')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
incV = 0.01;V = 0:incV:3;
%%%% Numeric differentiation
bnV = 1 - (1 + sqrt(2))^2./(1 + (4 + V.^4).^0.25).^2;
dVbn = diff(V.*bnV)/incV;
Vd2Vbn = V(1:end - 2).*diff(diff(V.*bnV))/incV^2;
figure(3)
plot(V(1:end - 2),bnV(1:end - 2),'-k ','LineWidth',2);hold on
plot(V(1:end - 2),dVbn(1:end - 1),'--r ','LineWidth',2);hold on
plot(V(1:end - 2),Vd2Vbn,'--g ','LineWidth',2);hold off
set(gcf,'Color','w');set(gca,'FontSize',16);
ylabel('\itb_n, d\rm( \itVb_n\rm )/\itdV','FontSize',18,'FontWeight','bold');
xlabel('\itV','FontSize',18,'FontWeight','bold')
axis([min(V) max(V) 0 max(Vd2Vbn)*1.05])
%%%% Analytic differentiation
syms V
bnV = 1 - (1 + sqrt(2))^2./(1 + (4 + V.^4).^0.25).^2;
dVbn = diff(V.*bnV, V);
Vd2Vbn = V.*diff(dVbn, V);
V = 0:incV:3;bnV = eval(bnV);dVbn = eval(dVbn);Vd2Vbn = eval(Vd2Vbn);
figure(4)
plot(V,bnV,'-k ','LineWidth',2);hold on
plot(V,dVbn,'--r ','LineWidth',2);hold on
plot(V,Vd2Vbn,'--g ','LineWidth',2);hold off
set(gcf,'Color','w');set(gca,'FontSize',16);
ylabel('\itb_n, d \rm\bf( \itV b_n\rm\bf ) / \itdV','FontSize',18,'FontWeight','bold');
xlabel('\itV','FontSize',18,'FontWeight','bold')
axis([min(V) max(V) 0 max(Vd2Vbn)*1.05])