Kidger triplet

rd = [ 0 18.609 460.967 -29.934 18.605 0 67.594 -22.162 0];
th = [ 0 2.5 5.157 1 1.207 4.007 2.5 41.60525 0];
rn = [ 1  1.620408 1 1.636355 1 1 1.620408 1 1];
sa = [ 0  7  7 5.1 5.1 5 6.2 6.2 0];
gn = char('SK16','F6','SK16');
ast = 6;

pya = [50/8 0];
pyc = [-5.12787 tan(22*pi/180)];

This is the paraxial raytrace from Kidger

cv = curvature(rd);
pya = parax(pya,cv,th,rn);
pyc = parax(pyc,cv,th,rn);

[pya pyc]

ans =

    6.2500         0   -5.1279    0.4040
    6.2500   -0.1286   -5.1279    0.3548
    5.9285   -0.2004   -4.2408    0.5693
    4.8951   -0.0589   -1.3050    0.3309
    4.8362    0.0691   -0.9741    0.5082
    4.9196    0.0691   -0.3606    0.5082
    5.1965    0.0132    1.6758    0.3041
    5.2295   -0.1250    2.4361    0.4246
    0.0287   -0.1250   20.1034    0.4246

This paraxial raytrace sets the stop per the lens drawing/specifications

k = -pyc(ast,1)/pya(ast,1)
pyc = pyc + k*pya;
[pya pyc]

k =

    0.0733


ans =

    6.2500         0   -4.6697    0.4040
    6.2500   -0.1286   -4.6697    0.3454
    5.9285   -0.2004   -3.8062    0.5546
    4.8951   -0.0589   -0.9462    0.3266
    4.8362    0.0691   -0.6195    0.5133
    4.9196    0.0691         0    0.5133
    5.1965    0.0132    2.0567    0.3051
    5.2295   -0.1250    2.8195    0.4155
    0.0287   -0.1250   20.1055    0.4155

Seidel wavefront coefficients

ford(pya, pyc, cv, th, rn, 0.58756e-3);

L =

    2.5252


W =

         0         0         0         0         0
   11.9025   21.7012    9.8917   60.7667   55.3965
    5.3684  -62.5597  182.2553   88.8742 -517.8351
  -22.6766  146.1077 -235.3472 -152.9207  492.6427
  -11.7029  -68.2875  -99.6162 -106.5180 -310.7718
         0         0         0         0         0
    1.4354   21.3874   79.6660   55.2008  411.2352
   19.3016  -61.6536   49.2338   71.4887 -114.1754
    3.6285   -3.3044  -13.9165   16.8916   16.4921


      W040       W131       W222       W220       W311
      3.63      -3.30     -13.92      16.89      16.49

load oldschott;

% look-up position of each glass in Schott glass table
n = size(gn);
for k=1:n(1)
    idx(k) = strmatch(gn(k,:),schott96.name);
end

wvln = linspace(0.4,0.7,21);
% calculate refractive index as function of wavelength
for k=1:length(idx)
    n = idx(k);
    rg(k,:) = sellmeier(schott96.disperse(n,:),wvln);
end

Find focus shift vs. wavelength

dz = zeros(size(wvln));
rnsave = rn;
n = length(rn);
lp = find(rn>1);
% use paraxial trace to find focus shift vs. wavelength
for k=1:length(wvln)
    for m = 1:length(lp)
        rn(lp(m)) = rg(m,k);
    end
    py = parax(pya,cv,th,rn);
    dz(k) = -py(n,1)/py(n,2);
end
rn = rnsave;

plot(dz,wvln,'k','LineWidth',2);
grid
xlabel('focus shift');
ylabel('wavelength ( \mum)');