/standalone/matlab/fuel_burnt_fd.m

function fuel_fraction=fuel_burnt_fd(lfn,tign,tnow,fd,fuel_time)

% lfn       - level set function, size (2,2)

% tign      - time of ignition, size (2,2)

% tnow      - time now

% fd        - mesh cell size (2)

% fuel_time - time constant of fuel

%

% output: approximation of

%                              /\

%                     1        |                   T(x,y)

%  fuel_burnt  =  ----------   |  ( 1  -  exp( - ----------- ) ) dxdy

%                 fd(1)*fd(2)  |                  fuel_time

%                             \/

%                        (x,y) in R;

%                         L(x,y)<=0

%

% where R is the rectangle [0,fd(1)] * [0,fd(2)]

%       T is given by tign and L is given by lfn at the 4 vertices of R

%

% note that fuel_left = 1 - fuel_burnt

%

% Requirements:

%       exact in the case when T and L are linear

%       varies continuously with input

%       value of tign-tnow where lfn>0 ignored

%       assume T=0 when lfn=0

figs=1:4;

for i=figs,figure(i),clf,hold off,end

% tign

% lfn

if all(lfn(:)>=0),

    % nothing burning, nothing to do - most common case, put it first

    fuel_fraction=0;

elseif all(lfn(:)<=0),

    % all burning

    % T=u(1)*x+u(2)*y+u(3)

    % t(0,0)=tign(1,1)

% t(0,msh_sz)=tign(1,2)

% t(msh_sz,0)=tign(2,1)

% t(msh_sz,msh_sz)=tign(2,2)

% msh_sz=h

% t(g/2,h/2)=sum(tign(i,i))/4

% dt/dx=(1/2h)*(t10-t00+t11-t01)

% dt/dy=(1/2h)*(t01-t00+t11-t10)

% approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences

% t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy

% u(1)=dt/dx

% u(2)=dt/dy

% u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy)



 tign_middle=(tign(1,1)+tign(1,2)+tign(2,1)+tign(2,2))/4;

 dt_dx=(tign(1,1)-tign(2,1)+tign(1,2)-tign(2,2))/(2*fd(1));

 dt_dy=(tign(1,1)-tign(1,2)+tign(2,1)-tign(2,2))/(2*fd(2));

%  dt_dx=(tign(2,1)-tign(1,1)+tign(2,2)-tign(1,2))/(2*fd(1));

%  dt_dy=(tign(1,2)-tign(1,1)+tign(2,2)-tign(2,1))/(2*fd(2));



 u(1)=dt_dx;

 u(2)=dt_dy;

 u(3)=tign_middle-(dt_dx*fd(1)+dt_dy*fd(2))/2;

 u1=u'

 display('fuel_burnt_fd coefficients of u')

    % integrate

    uu = -u / fuel_time;

    fuel_fraction = 1 - exp(uu(3)) * intexp(uu(1)*fd(1)) * intexp(uu(2)*fd(2));

    if(fuel_fraction<0 | fuel_fraction>1),

        warning('fuel_fraction should be between 0 and 1')

        stop

    end

else

    % part of cell is burning - the interesting case

    % set up a list of points with the corresponding values of T and L

    xylist=zeros(8,2);    % allocate list of points

    Tlist=zeros(8,1);

    Llist=zeros(8,1);

    xx=[-fd(1) fd(1) fd(1) -fd(1) -fd(1)]/2; % make center (0,0)

    yy=[-fd(2) -fd(2) fd(2) fd(2) -fd(2)]/2; % cyclic, counterclockwise

    ii=[1 2 2 1 1]; % indices of corners, cyclic, counterclockwise

    jj=[1 1 2 2 1];

    points=0;

    for k=1:4

        lfn0=lfn(ii(k),jj(k));

        lfn1=lfn(ii(k+1),jj(k+1));

        if(lfn0<=0),

            points=points+1;

            xylist(points,:)=[xx(k),yy(k)]; % add corner to list

            Tlist(points)=tnow-tign(ii(k),jj(k)); % time since ignition

            Llist(points)=lfn0;

        end

        if(lfn0*lfn1<0),

            points=points+1;

            % coordinates of intersection of fire line with segment k k+1

            %lfn(t)=lfn0 + t*(lfn1-lfn0)=0

            t=lfn0/(lfn0-lfn1);

            x0=xx(k)+(xx(k+1)-xx(k))*t;

            y0=yy(k)+(yy(k+1)-yy(k))*t;

            xylist(points,:)=[x0,y0];

            Tlist(points)=0; % now at ignition

            Llist(points)=0; % at fireline

     end

    end

    % make the lists circular and trim to size

    Tlist(points+1)=Tlist(1);

    Tlist=Tlist(1:points+1);

    Llist(points+1)=Llist(1);

    Llist=Llist(1:points+1);

    xylist(points+1,:)=xylist(1,:);

    xylist=xylist(1:points+1,:);

    for k=1:5,lfnk(k)=lfn(ii(k),jj(k));end

    figure(1)

    plot3(xx,yy,lfnk,'k')

    hold on

    plot3(xylist(:,1),xylist(:,2),Tlist,'m--o')

    plot3(xylist(:,1),xylist(:,2),Llist,'g--o')

    plot(xx,yy,'b')

    plot(xylist(:,1),xylist(:,2),'-.r+')

    legend('lfn on whole cell','tign-tnow on burned area',...

        'lfn on burned area', 'cell boundary','burned area boundary')

    patch(xylist(:,1),xylist(:,2),zeros(points+1,1),250,'FaceAlpha',0.3)

    patch(xylist(:,1),xylist(:,2),Tlist,100,'FaceAlpha',0.3)

    hold off,grid on,drawnow,pause(0.5)

 %%%%% Approximation of the plane for lfn using finite differences

% approximate L(x,y)=u(1)*x+u(2)*y+u(3)

 lfn_middle=(lfn(1,1)+lfn(1,2)+lfn(2,1)+lfn(2,2))/4;

 dt_dx=(lfn(1,1)-lfn(2,1)+lfn(1,2)-lfn(2,2))/(2*fd(1));

 dt_dy=(lfn(1,1)-lfn(1,2)+lfn(2,1)-lfn(2,2))/(2*fd(2));

 % I feel that it is right but not least squares tnow-tign

 u(1)=dt_dx;

 u(2)=dt_dy;

 u(3)=lfn_middle-(dt_dx*fd(1)+dt_dy*fd(2))/2;

% finding the coefficient c, reminder we work over one subcell only

% T(x,y)=c*L(x,y) in the sence of least squares

    a=0; b=0;

    for i=1:2

        for j=1:2

    if (lfn(i,j)<=0) 

                        a=a+lfn(i,j)*lfn(i,j);

                        if (tign(i,j)>tnow)

                        disp('fuel_burnt_fd: tign(i1) should be less then time_now');

                        else

                        b=b+(tign(i,j)-tnow)*lfn(i,j);

                        end

    end

        end

        end

                     if (a==0)

                        display('fuel_burnt_fd: if c is on fire then one of cells should be on fire');

                     end

                        c=b/a;

                        u(1)=u(1)*c;

                        u(2)=u(2)*c;

                        u(3)=u(3)*c;

    u=u'

    display('fuel_burnt_fd coefficients of u')

    nr=sqrt(u(1)*u(1)+u(2)*u(2));

    c=u(1)/nr;

    s=u(2)/nr;

    Q=[c,  s; -s, c];  % rotation such that Q*u(1:2)=[something;0]

     ut=Q*u(1:2);

    errQ=ut(2);  % should be zero

    ae=-ut(1)/fuel_time;

    ce=-u(3)/fuel_time;     %  -T(xt,yt)/fuel_time=ae*xt+ce

    cet=ce;  %keep ce for later

    xytlist=xylist*Q'; % rotate the points in the list

    xt=sort(xytlist(1:points,1)); % sort ascending in x

    fuel_fraction=0;

    for k=1:points-1

        % integrate the vertical slice from xt1 to xt2

        figure(2)

        plot(xytlist(:,1),xytlist(:,2),'-o'),grid on,hold on

        xt1=xt(k);

        xt2=xt(k+1);

        plot([xt1,xt1],sum(fd)*[-1,1]/3,'y')

        plot([xt2,xt2],sum(fd)*[-1,1]/3,'y')

        if(xt2-xt1>100*eps*max(fd)), % slice of nonzero width

            % find slice height as h=ah*x+ch

            upper=0;lower=0;

            ah=0;ch=0;

            for s=1:points % pass counterclockwise

                xts=xytlist(s,1); % start point of the line

                yts=xytlist(s,2);

                xte=xytlist(s+1,1); % end point of the line

                yte=xytlist(s+1,2);

                if (xts>xt1 & xte > xt1) | (xts<xt2 & xte < xt2)

                    plot([xts,xte],[yts,yte],'--k')

                    % that is not the one

                else % line y=a*x+c through (xts,yts), (xte,yte)

                    a=(yts-yte)/(xts-xte);

                    c=(xts*yte-xte*yts)/(xts-xte);

                    if xte<xts % upper boundary

                        aupp=a;

                        cupp=c;

                        plot([xts,xte],[yts,yte],'g')

                        ah=ah+a;

                        ch=ch+c;

                        upper=upper+1;

                    else % lower boundary

                        alow=a;

                        clow=c;

      plot([xts,xte],[yts,yte],'m')

                        lower=lower+1;

                    end

                end

            end

            if(lower~=1|upper~=1),

                error('slice does not have one upper and one lower line')

            end

            ce=cet;% use kept ce

            % shift samll amounts to avoid negative fuel burnt

            if ae*xt1+ce > 0;%disp('ae*xt1+ce');disp(ae*xt1+ce);pause;

               ce=ce-(ae*xt1+ce);end;

            if ae*xt2+ce > 0;%disp('ae*xt2+ce');disp(ae*xt2+ce);pause;

               ce=ce-(ae*xt2+ce);end;



            ah=aupp-alow;

            ch=cupp-clow;

            % debug only

            patch([xt1,xt2,xt2,xt1],...

                [alow*[xt1,xt2]+clow,aupp*[xt2,xt1]+cupp,],k*10)

            axis equal,drawnow, pause(0.5)

            figure(3)

            x=[xt1:(xt2-xt1)/100:xt2];

            plot(x,1-exp(ae*x+ce),x,ah*x+ch,x,...

                (ah*x+ch).*(1-exp(ae*x+ce)),[xt1,xt2],[0,0],'+k')

            xlabel x,legend('burned frac','slice height','height*burned','dividers')

            hold on,grid on,drawnow,pause(0.5)

            % integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2

            % numerically sound for ae->0, ae -> infty

            % this can be important for different model scales

            % esp. if someone runs the model in single precision!!

            % s1=int((ah*x+ch),x,xt1,xt2)

            s1=(xt2-xt1)*((1/2)*ah*(xt2+xt1)+ch);

            % s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)

            ceae=ce/ae;

            s2=-ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1));

            % s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)

            % s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)

            % expand in Taylor series around ae=0

            % collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)

            % =(1/8*xt1^4+1/3*xt1^3*ceae+1/4*xt1^2*ceae^2-1/8*xt2^4-1/3*xt2^3*ceae-1/4*xt2^2*ceae^2)*ae^2

 %     + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae

            %     + 1/2*xt1^2-1/2*xt2^2

            %

            % coefficient at ae^2 in the expansion, after some algebra

            a2=(xt1-xt2)*((1/4)*(xt1+xt2)*ceae^2+(1/3)*(xt1^2+xt1*xt2+xt2^2)*ceae+(1/8)*(xt1^3+xt1*xt2^2+xt1^2*xt2+xt2^3));

            d=ae^4*a2;

            if abs(d)>eps

                % since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae

                % for ae, ce -> 0 rounding error approx eps/ae^2

                s3=(exp(ae*(xt1+ceae))*(ae*xt1-1)-exp(ae*(xt2+ceae))*(ae*xt2-1))/ae^2;

                % we do not worry about rounding as xt1 -> xt2, then s3 -> 0

            else

                % coefficient at ae^1 in the expansion

                a1=(xt1-xt2)*((1/2)*ceae*(xt1+xt2)+(1/3)*(xt1^2+xt1*xt2+xt2^2));

                % coefficient at ae^0 in the expansion for ae->0

                a0=(1/2)*(xt1-xt2)*(xt1+xt2);

                s3=a0+a1*ae+a2*ae^2; % approximate the integral

            end

            s3=ah*s3;

            fuel_fraction=fuel_fraction+s1+s2+s3;

            if(fuel_fraction<0 | fuel_fraction>(fd(1)*fd(2))),

                fuel_fraction,(fd(1)*fd(2)),s1,s2,s3

                warning('fuel_fraction should be between 0 and 1')

                stop

            end

        end

    end

    fuel_fraction=fuel_fraction/(fd(1)*fd(2));

    % display('fuel_burnt_fd')

    %fuel_fraction=1-fuel_fraction;

end

for i=figs,figure(i),hold off,end

end % function



function s=intexp(ab)

% function s=intexp(ab)

% s = (1/b)*int(exp(a*x),x,0,b)  ab = a*b

%    = 1 if a==0

if eps < abs(ab)^3/6,

    s = (exp(ab)-1)/(ab);  % rounding error approx eps/(a*b) for small a*b

else

    s = 1 + ab/2+ab^2/6;   % last term (a*b)^2/6 for small a*b

end

end