% Fit index of refraction  versus wavelength to the Cauchy
% formula a+b/lambda^2. Find the index of refraction n and the
% dispersion dn/dlambda at 550 nm.

% Enter the data (these are dummy results). y are the indices of refraction
% and x are the wavelengths measured with the grating.
% Make sure you have the same number of points in the x and y vectors.
y  = [1.6201 1.6231 1.6287 1.6383 1.6393 1.6408 1.6433 1.6483 1.6507]
x  = [707.9  669.3  588.9  506.2  502.6  493.1  472.4  447.6  439.4]
dtheta = sqrt(power(0.1,2)+power(0.1,2))/2*pi()/180
% Since n=2*sin(theta/2+pi/6),
% dn=cos(theta/2+pi/6)= sqrt(1-sin²(theta/2+pi/6))*dtheta=sqrt(1-(n/2)²)*dtheta.
dy = dtheta*sqrt(1-power(y/2,2))
% Fit the data to a + b/lambda^2
% The resulting coefficients are stored in coefs.
% The starting values are just a very rough guess to help
% speed up the fitting procedure. The procedure assumes equal weights for
% the moment - have to include the errors in a more advanced version.
fun = @(p,lambda)p(1) + p(2)*power(lambda,-2);
startingValues = [1.5,10];
coeffs = nlinfit(x,y,fun,startingValues)
a = coeffs(1)
b = coeffs(2)

% Plot the data as o's and the fitted function as a curve
plot(x,y,'o')
errorbar(x,y,dy,'Marker','o','LineStyle','none')
hold on
fittedfun = @(lambda)a + b*power(lambda,-2)
fplot(fittedfun,xlim)

% Find the index of refraction at 589 nm.
n589=fittedfun(589)
ylim = get(gca,'Ylim')
line([589,589],[ylim(1),n589])
xlim = get(gca,'Xlim')
line([xlim(1),589],[n589,n589])
text(589+(xlim(2)-xlim(1))*0.040,n589 + 0.120*(ylim(2)-ylim(1)),strcat('n_p(\lambda=589 nm)= ',num2str(n589)))

% Find the disperson at 589 nm.
dndl589 = -2*b*power(589,-3)
text(589+(xlim(2)-xlim(1))*0.040,n589 + 0.040*(ylim(2)-ylim(1)),strcat('dn_p/d\lambda(\lambda=589 nm)= ',num2str(dndl589)))

% Plot the dispersion as a short tangent
tangent = @(lambda) n589 + (lambda-589)*dndl589
line([500,700],[tangent(500),tangent(700)])

% Choose your own histogram axis labels and title
xlabel('\lambda (nm)')
ylabel('n_p')
title('Index of refraction n_p versus wavelength \lambda for prism')
