function [Q rect] = m3spacecharge()
d = imread('face.png');
d = double(rgb2gray(d));
d = d/max(d(:));
d = abs(1-d); d = flipud(d);
d = d(1:2:end-2,1:2:end);
d = padarray(d, [50 90]);


size(d)

v = linspace(-1,1,size(d,1)+1); v = mean([v(1:end-1); v(2:end)]);
x = linspace(-1,1,size(d,2)+1); x = mean([x(1:end-1); x(2:end)]);
[v0,x0] = ndgrid(v,x);
xscale = 1.5; x0 = x0*xscale;
vscale = 1; v0 = v0*vscale;

%surf(x,v,d);

t = linspace(0, 10, 801);
dt = t(2)-t(1);

x = x0; v = v0;  a = getacc(x,v,d,x0,v0);
e(1) = smc(d);

pprof = zeros(534,length(t));
profscale = 77.2118;
profheight = 200;

for ii = 1:length(t)
    xn = x - (v + 0.5*a*dt)*dt;
    an = getacc(xn,v,d,x0,v0);
    vn = v + 0.5*(a + an)*dt;

    dn = interp2(x0, v0, d, xn, vn);
    dn(isnan(dn)) = 0;
    e(ii) = smc(dn);
    
    dn = padarray(dn, [0 100]);
    dn = padarray(dn, [50 0], 'pre');
    dn = padarray(dn, [150 0], 'post');

    prof = sum(dn,1).';
    pprof(:,ii) = -prof/profscale;

    pp = resizem(pprof, [profheight, size(dn,2)]);
    pp = padarray(pp, [size(dn,1)-size(pp,1) 0], 'pre');
    
    f1 = figure(1); set(f1, 'Color', [1 1 1]);
    set(f1,'DoubleBuffer','on');
    imagesc(dn+pp);
    %title(sprintf('t=%.2f',t(ii)));
    %s(ii) = sum(dn(:));
    axis xy; caxis([-1 1]); colormap(graywhitered(256));
    hold on;
    plot(1:size(dn,2), prof,'k','LineWidth', 2);
    plot([1 size(dn,2)], [0 0],'k');
    plot([1 size(dn,2)],size(dn,1)-profheight/2*[1 1],'b--');
    hold off;
    axis off;
    truesize(f1);
    
    set(gcf,'PaperPositionMode','auto');
    drawnow;    
    rect = get(f1, 'Position');
    rect(1:2) = [1 1];
    rect = [50 50 rect(3)-100 rect(4)-60];
    Q(:,ii) = getframe(f1, rect);

    x = xn;
    v = vn;
    a = an;
%    figure(2)
    

    %disp('press a key');
    %pause();
end

%figure; plot(t, e/e(1));

end

function a = getacc(x,v,d,x0,v0)
% focusing term
foc = 0.75;
% space charge term
sch = 0.01;

dn = interp2(x0, v0, d, x, v);
dn(isnan(dn)) = 0;

x0 = x0(1,:); dx = x0(2)-x0(1);
c = sum(dn,1); cs = cumsum(c); cs = [0 cs cs(end)]*dx;

%figure(2)
%plot(x0, c);

f = zeros(size(x0));
for ii = 1:length(x0)
    f(ii) = cs(ii) + cs(ii+1) - cs(end);
end

sp = interp1(x0, f, x);
p = isnan(sp); p2 = x > 0;
sp(p & p2) = f(end);
sp(p & ~p2) = f(1);

%figure(3)
%plot(x0, f);

a = foc*x - sch*sp;
end