function [Q rect] = m1drift()

d = imread('face.png');
d = double(rgb2gray(d));
d = d/max(d(:));
d = abs(1-d); d = flipud(d);

%h = figure(1);
%imagesc(d); axis xy; caxis([0 1]); colormap(bluewhitered(256));
%truesize(h);

% inverse mapping function
%m = [1 -0.5; 0 1];
%m = [1 0; 1 1];
%m = [1 0; -1 1];

% setup beamline
c.L = [1 1 1 1 1 1 1 1 1];
c.k = [0 1 0 1 0 1 0 1 0];

z = linspace(0, sum(c.L), 700);
d1 = padarray(d, [40 160]);
dn = phasespacetrans(d1, eye(2));
dns = size(dn);
pprof = zeros(dns(2),length(z));
profscale = 165;
profheight = 130;
for ii = 1:length(z)

    dn = phasespacetrans(d1, whichelem(c,z(ii)));
    
    dn = fliplr(dn);
    prof = sum(dn,1).';
    pprof(:,ii) = -prof/profscale;
    pp = resizem(pprof, [profheight, dns(2)]);
    dn = padarray(dn, [100 0], 'post');
    dn = padarray(dn, [50 0], 'pre');
    pp = padarray(pp, [size(dn,1)-size(pp,1) 0], 'pre');

    f3 = figure(3); set(f3, 'Color', [1 1 1]);
    imagesc(dn+pp); axis xy; caxis([-1 1]); colormap(graywhitered(256));
    truesize(f3);
    hold on;
    plot(1:length(prof), prof*0.5,'k','LineWidth',2);
    %plot(sum(dn,2),1:size(dn,1),'k');
    %imagesc(d); colormap(gray);
    drawbeamline(c, size(dn), profheight);
    hold off; axis off; drawnow; 
    rect = get(f3, 'Position');
    rect(1:2) = [1 1];
    Q(:,ii) = getframe(f3, rect); clf;
    %Q(:,ii) = getframe(f3);
end

end

function drawbeamline(c, sd, ph)
Ltot = cumsum(c.L);
ind = c.k ~= 0;
lbeg = Ltot(find(ind)-1)/Ltot(end)*sd(2);
lend = Ltot(find(ind))/Ltot(end)*sd(2);
plot([1 sd(2)],[sd(1)-ph/2 sd(1)-ph/2],'b--');
pp = 0:ph; xx = (pp/(ph/2)-1).^2;
for ii = 1:length(lbeg)
    plot(-xx*(lend(ii)-lbeg(ii))/2+lend(ii), sd(1)-pp, 'b', 'Linewidth',2);
    plot(xx*(lend(ii)-lbeg(ii))/2+lbeg(ii), sd(1)-pp, 'b', 'Linewidth',2);    
end
    
end

function m = whichelem(c, z)
Ltot = cumsum(c.L);

aftel = z > Ltot;
Ltotplus = [0 Ltot];

m = eye(2);
for ii = 1:length(aftel)
    if aftel(ii)
        m = m*genmat2(c.L(ii), c.k(ii));
    else
        m = m*genmat2(z-Ltotplus(ii), c.k(ii));
        break;
    end
end

%det(m)

end

function m = genmat2(L, k)
m = [cos(k*L) sinc(k*L/pi)*L; -k*sin(k*L) cos(k*L)];
end

function out = phasespacetrans(in, m)

f = @(x,unused) [m(1,1)*x(:,1)+m(1,2)*x(:,2), m(2,1)*x(:,1)+m(2,2)*x(:,2)];
% maketform arguments
ndims_in = 2;
ndims_out = 2;
forward_mapping = [];
inverse_mapping = f;
tdata = [];
tform = maketform('custom', ndims_in, ndims_out, ...
    forward_mapping, inverse_mapping, tdata);
out = imtransform(in, tform, 'UData', [-1 1], 'VData', [-1 1], ...
    'XData', [-1.4 1.4], 'YData', [-1 1]);

out = resizem(out,[size(out,1) fix(size(out,2)/1.5)-1]);

end