function run()
% Read in expt data and compare to expected shape

% RI versus wavelength data
lp=[2325.4,1970.1,1529.6,1060.0,1014.0,852.1,706.5,656.3,643.8,632.8,589.3,587.6,546.1,486.1,480.0,435.8,404.7];
np=[1.68218,1.68750,1.69378,1.70227,1.70345,1.70887,1.71661,1.72085,1.72200,1.72309,1.72803,1.72825,1.73430,1.74648,1.74805,1.76198,1.77579];

% Large angle data messes up
%ind=n<1.78;
%t=t(ind); n=n(ind);

% Resample data to regular spacing
l=410:10:2300;

n=1:10;

for i=n
  [v{i}]=getrun(i,lp,np,l);
end

if(1)
figure(5)
clf
hold on
for i=n
  plot(l,v{i});
  ind=l>600&l<1700;
  [p(i,:),pe,chi2(i)]=matmin('chisq',[2100,0.36e-11],[1,0],'planck2_fit',v{i}(ind),ones(size(v{i}(ind))),l(ind)*1e-9);
  plot(l(ind),planck2_fit(p(i,:),l(ind)*1e-9),'r');
  plot(l,planck2_fit(p(i,:),l*1e-9),'r--');
end
hold off
xlim([400,2300]);
grid; box on
xlabel('Wavelength (nm)'); ylabel('Light Intensity (arb units)');
end

sum(chi2)

figure(6)
plot(n,p(:,1),'*');
grid
xlabel('Lamp voltage (V)'); ylabel('Inferred Filament Temperature (K)');

% Draw what Planck curves should look like:
figure(7)
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}');

return

function [v]=getrun(i,lp,np,l)

% Read data
[t,v]=textread(sprintf('run1/%dv.txt',i));

% Throw out redundant obs
[t,i]=unique(t);
v=v(i);

% Translate angle to RI
n=sqrt(((2/sqrt(3)).*sin(t*pi/180)+1/2).^2+3/4);

% Disguard data out of range
ind=n>np(1)&n<np(end);
t=t(ind); n=n(ind); v=v(ind);

% Translate RI to wavelength
ldat=interp1(np,lp,n,'spline');

% Fit from pasco guy
if(0)
A=6015783.2978286;
B=-13717512.2572366;
C=11725406.1052594;
D=-4452989.0254548;
E=633985.0106661;
ldat=3000./sqrt(A+B*n+C*n.^2+D*n.^3+E*n.^4);
end

figure(1)
plot(lp,np,'*');
hold on;
plot(ldat,n,'r')
hold off;
grid
xlabel('Wavelength (nm)'); ylabel('Refractive index');

figure(2)
plot(t,n)
grid
xlabel('Sensor angle (deg)'); ylabel('Refractive Index');

figure(3)
plot(t,ldat)
grid
xlabel('Sensor angle (deg)'); ylabel('Wavelength (nm)');

% Resample to regular spacing
v=interp1(ldat,v,l);

% Slit width is constant in angle - estimate width
% in units of wavelength
tup=t+0.1;
ldatup=interp1(t,ldat,tup);
wid=ldat-ldatup;

figure(4)
plot(ldat,wid);
grid
xlabel('Wavelength'); ylabel('Wavelength range corresponding to 0.1 deg (nm)');

% Correct for slit width
wid=interp1(ldat,wid,l);
wid=wid/mean(wid);
v=v./wid;

return

