function burn(N,fps)
%
%  Burning heart (ver 1.01)
%
% --------------------------------------------------------------------------------
% Copyright (c) 2013, Alessandro Foi
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%

if ~exist('N','var')
    N=384;   %% image has size NxN
end
if ~exist('fps','var')
    fps=20;  %% maximum frames per second
end
maxIterPerFrame=3;

disp(' Close figure to terminate execution ' );

%% assemble heart DY DX
L=2; C=sqrt(2); T=3/4; O=(T*sqrt(2*C+1)-0.5)/2;
[X,Y]=meshgrid(L*linspace(-1,1,N),O+L*linspace(-1,1,N));
DY=min(abs(Y-T*sqrt(-C^2+(-abs(X)+1+C).^2).*(abs(X)<1)),abs(-Y-0.5*sqrt(abs(1-abs(1-2*abs(X)).^2)).*(abs(X)<1)));
DX=abs(-abs(X)-sqrt((Y/T).^2+C^2)+C+1).*((Y)>=0)+min(abs(-abs(X)+0.5-0.5*sqrt((Y>-0.5).*abs(1-abs(1-2*abs(Y-0.5)).^2))),abs(abs(X)-0.5-0.5*sqrt((Y>-0.5).*abs(1-abs(1-2*abs(Y-0.5)).^2)))).*((Y)<0);

%% heart
matlabOctaveVer=ver;
a=2; b=1/3; c=-1/a*2;
theta=atan2(X,Y-O);
D=(abs(DY.*cos(theta)).^a+abs(DX.*sin(theta)).^a).^b;
clear DX DY X Y
D=D.^c*4;

%% burning
P0=D.*randn(size(D)).^2;
KS=[1 0.5];
KN=ceil(2*KS);
[KX,KY]=meshgrid(-KN:KN,-KN:KN);
KC=[-KS(1)*1.5,0];
kernel2=exp(-(((KY-KC(1))/(2*KS(1))).^2+((KX-KC(2))/(2*KS(2))).^2))/(2*pi*KS(1)*KS(2));
kernel2=kernel2/sum(kernel2(:));
figEveryNsecs=1/fps;
oldTime=now-figEveryNsecs/(60*60*24);
jj=0; ii=1;
dampedNCycles=1;
waitExtra=0;
figH=figure;
axisH=subplot(1,1,1);
probImpulses=0.1;
sigmaImpulses=5;
Nimpulses=ceil(numel(D)*probImpulses);
lambda=.9;
while ishandle(figH)
    jj=jj+1;
    P1=D.*randn(size(D)).^2;
    impp=randi(numel(P1),Nimpulses,1);
    P1(impp)=P1(impp)+rand(Nimpulses,1).*D(impp).*randn(Nimpulses,1)*sigmaImpulses;
    P0=P0*lambda+P1*(1-lambda)/2;
    ii=ii+1;
    currentTime=now;
    timeElapsed=(currentTime-oldTime)*60*60*24;
    if timeElapsed > figEveryNsecs
        oldTime=currentTime;
        if jj==1;
            imageH=image(P0,'parent',axisH);
        elseif ishandle(imageH)
            set(imageH,'CData',P0);
        end
        if jj==1
            colormap(hot); axis equal tight off;
            if ~strcmp(matlabOctaveVer(1).Name,'Octave')
                truesize;
            end
        end
        dampedNCycles=dampedNCycles*0.3+0.7*ii;
        drawnow
        pause(eps)
        waitExtra=max(0,waitExtra+0.3*(timeElapsed-(currentTime-oldTime)*60*60*24-waitExtra)/dampedNCycles*(ii-maxIterPerFrame));
        if waitExtra
            pause(waitExtra);
        end
        ii=0;
    end
    P0=conv2(P0,kernel2,'same');
end

end

