/standalone/matlab/prop_ls_cir.m

function [t,phi]=prop_ls(phi0,time0,time1,vx,vy,r,dx,dy,spread_rate)

% a jm/mkim version of propagate

% follows Osher-Fedkiw and Mitchell toolbox



% phi       level function out

% phi0      level function in

% time0     starting time

% time1     end time

% vx        component x of velocity field, passed to spread_rate

% vy        component y of velocity field, passed to spread_rate

% r         propagation speed in normal direction

% dx        mesh step in x direction

% dy        mesh step in y direction

% 



% initialize

phi=phi0; % allocate level function

[m,n]=size(phi);

t=time0;  % current time

tol=300*eps;

istep=0;

msteps=1000;

diffLx=zeros(m,n);

diffLy=zeros(m,n);

diffRx=zeros(m,n);

diffRy=zeros(m,n);



method='godunov'

split='spread';

% split='wind';

% split = 'orig';



% time loop

while t<time1-tol & istep < msteps

    istep=istep+1;

    % one-sided differences

    [diffLx,diffRx,diffLy,diffRy,diffCx,diffCy]=get_diff(phi,dx,dy);         

    % gradient of level function is the normal direction

    scale=sqrt(diffCx.*diffCx+diffCy.*diffCy+eps); 

    nvx=diffCx./scale;

    nvy=diffCy./scale;

    % propagation speed 

    %speed=spread_rate(r,vx,vy,nvx,nvy,scale);

    %

    speed=vx.*nvx + vy.*nvy;

    speed=max(speed,0);

    % to recover advection vv and spread r, transition between:

    % r = 0 & (normal,v) > const*speed => vv=speed*v/(normal,v) & rr=0

    % speed >> (normal,v) => vv=0 & rr=speed    

    nvv = nvx .* vx + nvy .* vy;

    % the CIR method

    switch method

        case 'cir'

            dt=min(min(dx./max(abs(speed.*nvx),eps)))+min(min(dy./max(abs(speed.*nvy),eps)));

            dt=min(time1-t,0.5*dt)

            % locations back on characteristics, in mesh index coordinates

            [xii,yii]=meshgrid(1:m,1:n);

            xinc=dt.*speed.*nvx./dx;

            xi = xii - xinc;

            yinc=dt.*speed.*nvy./dy;

            yi = yii - yinc; 

            phi_new=interp2(xii,yii,phi,xi,yi);

            phi=min(phi_new,phi);

        case 'godunov'

            switch split

            case 'wind'

                a=2*nvv>speed;

                nvv(a==0)=1;

                ! a=zeros(size(a));

                %no_wind_cases=sum(sum(1-a))

                rr=speed.*(1-a);

                corr = a .* speed ./ nvv

                vvx = vx .* corr;

                vvy = vy .* corr;

            case 'orig'

                vvx = vx;vvy= vy;rr=r;  % all original, if r=0 pure upwinding

            case 'spread'

                rr = speed; vvx=0*vx; vvy=0*vy; % all in normal speed, Godunov meth.

                % vvx=vx;vvy=vy;

            otherwise

                error('unknown split')

            end

            % Godunov scheme for normal motion

            flowLx=(diffLx>=0 & diffCx>=0);

            flowRx=(diffRx<=0 & diffCx<0);

            diff2x=diffLx.*flowLx + diffRx.*flowRx; %

            flowLy=(diffLy>=0 & diffCy>=0);

            flowRy=(diffRy<=0 & diffCy<0);

            diff2y=diffLy.*flowLy + diffRy.*flowRy; %

            grad=sqrt(diff2x.*diff2x + diff2y.*diff2y);

            nz=find(grad);

            %tbound_n(1)=max(abs(r.*sqrt(diff2x(nz)))./grad(nz))/dx; % time step bnd

            %tbound_n(2)=max(abs(r.*sqrt(diff2y(nz)))./grad(nz))/dy;

            tbound_np=rr.*(abs(diff2x)./dx+abs(diff2y)./dy); % pointwise time step bnd

            tbound_n=max(tbound_np(nz)./abs(grad(nz))); % worst case

            tend_n =-rr.*grad;   % the result, tendency

            % standard upwinding for advection

            tend_a=-(diffLx.*max(vvx,0)+diffRx.*min(vvx,0)+...

            diffLy.*max(vvy,0)+diffRy.*min(vvy,0));

            tbound_ax=max(abs(vvx(:)))/dx;

            tbound_ay=max(abs(vvy(:)))/dy;

            % complete right-hand side

            tend=tend_n + tend_a;

            % complete time step bound

            tbound = 1/(tbound_n+tbound_ax+tbound_ay+eps); 

            % decide on timestep

            dt=min(time1-t,0.5*tbound);

            % trailing edge correction - do not allow fireline to go backwards

            tend=min(tend,0);

            % advance

            phi=phi+dt*tend;

    end

    t=t+dt;

end

fprintf('prop_ls: %g steps from %g to %g, split %s\n',istep,time0,t,split)

end % of the function prop_ls