script2

Contents

fit ray errors

clear;
load rayfan_f2doublet

yfit =  [ y y.^3 y.^5 y.^7 y.^9 ];
sz=size(yfit);
n = sz(2);
rtable = zeros(n,n);
close all
for k=1:n
    c = yfit(:,1:k)\dy;
    rtable(k,1:k)=c';
    yf = yfit(:,1:k)*c;
    subplot(3,2,k);
    %plot(y,dy,y,yf);
    fanplot(y,[dy yf]);
    title(sprintf('order %d',2*k-1));
end
rtable

rtable =

    0.0023         0         0         0         0
   -0.1436    0.2318         0         0         0
    0.0049   -0.4311    0.5702         0         0
   -0.0286   -0.1406   -0.0440    0.3650         0
   -0.0234   -0.2150    0.2371   -0.0231    0.1770

quality of fit

close all;
order = 1:2:9;
qual = zeros(size(order));
ref = max(abs(dy));
for k=1:length(order)
    c = rtable(k,1:k)';
    yf = yfit(:,1:k)*c;
    diff = dy - yf;
    qual(k) = max(abs(diff))/ref;
end
stem(order,qual,'k','LineWidth',2);
grid;
xlabel('order');
ylabel('maximum relative difference');
axis([0 10 0 1]);

fit opd

yfit =  [ ones(size(y)) y.^2 y.^4 y.^6 y.^8 y.^10];
sz=size(yfit);
n = sz(2);
table = zeros(n,n);
for k=1:n
    c = yfit(:,1:k)\opd;
    table(k,1:k)=c';
    yf = yfit(:,1:k)*c;
    subplot(3,2,k);
    %plot(y,opd,y,yf);
    fanplot(y,[opd yf],10);
    title(sprintf('order %d',2*k-2));
end
table

table =

    3.2729         0         0         0         0         0
    0.3335    8.6031         0         0         0         0
   -0.7104   19.0291  -11.7089         0         0         0
    0.0978    1.6588   38.8640  -35.7055         0         0
   -0.0093    5.7093   17.0337    1.0004  -19.0008         0
    0.0019    5.0442   22.7466  -15.7789    0.8558   -8.1496

quality of fit

close all;
order = 0:2:10;
qual = zeros(size(order));
ref = max(abs(opd));
for k=1:length(order)
    c = table(k,1:k)';
    yf = yfit(:,1:k)*c;
    diff = opd - yf;
    qual(k) = max(abs(diff))/ref;
end
stem(order,qual,'k','LineWidth',2);
grid;
xlabel('order');
ylabel('maximum relative difference');
axis([0 12 0 1]);

get OPD from ray error

fact = [2 4 6 8 10 ];
fact = ones(size(fact'))*fact;

% ua for f/2 system
ua = 0.25;
lambda = 0.5876e-3;
epsilon = lambda/ua;
ifact = (1./fact);
otable = -rtable.*ifact/epsilon;
otable

otable =

   -0.4940         0         0         0         0
   30.5405  -24.6501         0         0         0
   -1.0492   45.8578  -40.4347         0         0
    6.0936   14.9500    3.1198  -19.4128         0
    4.9777   22.8732  -16.8106    1.2287   -7.5312

% show differences
for k=1:4
    subplot(2,2,k);
    n = k+1;
    c = table(n,1:n)';
    y1 = yfit(:,1:n)*c;
    c = [0; otable(k,1:k)'];
    y2 = yfit(:,1:k+1)*c;
    fanplot(y,[opd y1 y2],10);
    title(sprintf('order %d',2*k+2));
end

Integrating the ray error fit to get an OPD fit does not work well unless the quality of fit is good.

Published with MATLAB® 7.6