% http://www.uccs.edu/~pbhallam/documents/PES315/PES315_LB_Blackbody.pdf
% http://hyperphysics.phy-astr.gsu.edu/hbase/mod6.html#c3

l=[404.7,435.8,480.0,486.1,546.1,587.6,643.8,656.3,706.5,852.1,1014.0];
n=[1.776,1.762,1.748,1.746,1.734,1.728,1.722,1.721,1.717,1.709,1.703];

%ind=l>600;
%l=l(ind); n=n(ind);

% Fit RI data to Cauchy form
%p=matmin('chisq',[1,1],[],'cauchy',n,ones(size(n)),l);
p=matmin('chisq',[1.7,0],[],'poly',n,ones(size(n)),l);

plot(l,n,'.-')
hold on
l=400:20:4000;
%plot(l,cauchy(p,l),'r');
plot(l,poly(p,l),'r');
hold off

l=400:20:40000;

% Come in at fixed angle
t1=ones(size(l))*pi/3; % 60deg
n=cauchy(p,l);

% Propogate through prism
t2=asin(sin(t1)./n);
t3=pi/3-t2;

% Angle out of prism
t=asin(n.*sin(t3));

% Wavelength as a fn of angle
plot(t*180/pi,l); grid

% Convert to linear position along screen in cm at 1m distance
dt=t-t(1);
x=-100*tan(dt);

% Wavelength as a function of linear position
plot(x,l);

% Intensity in W m^-2 nm^-1
I=planck2(3000,l/1e9)/1e9;

% Planck function versus wavelength
plot(l,I);

% This is what it measured curve should look like!
plot(x,I);


% 11/06/02
% Got the fancy Pasco setup but it's shit!

% Draw what Planck curves should look like:
l=400:20:2500;
lmax=[661,741,803,862,933,1046]*1e-9;
lmax=[1700,1400,1300,1200]*1e-9;
T=0.0029./lmax;
for i=1:length(T)
  I(i,:)=pi*planck2(T(i),l/1e9)/1e9;
end
plot(l,I'); grid
xlabel('\lambda (nm)'); ylabel('Emission from surface in W m^{-2} nm^{-1}');

l=200:20:1200;
I=pi*planck2(5800,l/1e9);
