%Heterogeneous Fuel Grain Burn Simulator v1.5
%Simple version with standard algorithm
%By Graham Orr
%Revised on 25 Feb 2009
%www.thirdcritical.com

%This outputs an animation, "thrust curve", and "topo map" for an input geometry.
%It uses the finite difference method.
%The input matrix is just 1s and 0s arranged for fuel and bore,
%   respectively.
%All lengths are nondimensionalized to R. that way this is applicable to
%   larger and smaller motors.
%Surface area is nondimensionalized as SA=A/(GrainLength*2*pi*R).
%Time is nondimensionalized as tau=n/R (relative to the time it takes a
%   pinpoint to burn outward to radius R).
%GrainLength input is also nondimensionized and relative to R. Example:
%   grainLength=2 means that the grain length is twice the radius.
%MakeBore.m can make input arrays for simple shapes. Just run the script in
%   the command window in form: [B SA]=burnArea(makeBore(0.5,3),4);.
%Have fun!
%///////////////////////////E80 ROCK(et)S!!!!!!!!!!!!!!!


function [SA D Graham_Score] = burnArea(A,grainLength)
runAnimation=1;
close(figure(1));
close(figure(2));
[w h]=size(A);
if (runAnimation)
    f= figure('Name','BurnSim for E80, HMC. (c) Graham Orr, 2009.','Position', [100, 200, 3*w+100, h+100],'Resize','off');
end

N=150;  %maximum length of simulation
B=zeros(w,h,N); %burn matrix (contains geometry data indexed at time N)
C=zeros(w,h);   %will become the initial matrix to determine corners
SA=zeros(N,2);  %the compiled surface area array with time index
count=zeros(w,h);   %an array that comes from C...like the numbers in minesweeper
R=w/2;  %sets the discrete radius to half the width of the input array A
outerlimit=zeros(w,h);  %this matrix will just be ones in circle of radius R
B(:,:,1)=A(:,:);    %makes the first entry in the burn matrix equal to initial conditions
for i=1:w;
    for j=1:h
        if (w/2-i)^2+(h/2-j)^2 <= R^2
            outerlimit(i,j)=1;  %here we define this array as 1s and 0s
        end
    end
end
    for i=2:w-1
        for j=2:h-1
            if (A(i,j)==1) 
                count(i,j)=0;
                if A(i+1,j,1)==1
                    count(i,j)=count(i,j)+1;    %this is basically the minesweeper count routine
                else end
                if A(i,j,1)==1
                    count(i,j)=count(i,j)+1;    %how many 1s are around at point i,j
                else end                
                if A(i-1,j,1)==1
                    count(i,j)=count(i,j)+1;    %which can be used to find corners
                else end
                if (A(i,j+1,1)==1) 
                    count(i,j)=count(i,j)+1;    %because corners expand as circles
                else end
                if (A(i,j-1,1)==1) 
                    count(i,j)=count(i,j)+1;    %will all other things, we'll just assume linear expansion
                else end

            else end
        if (count(i,j)==3 || count(i,j)==1)     %corner is defined by there being 1 or 3 in proximity
            C(i,j)=1;
        else end
       
        end
    end
activation=0;
edgeIndex=0;
for n=2:N
    if runAnimation
    wb=waitbar(2*n/N);
    end
    for i=2:w-1
        for j=2:h-1
            if (B(i,j,n-1)==1)
                B(i+1,j,n)=1;   %if it is not a corner, we advance B in all directions
                B(i-1,j,n)=1;   %B just contains 1s and 0s
                B(i,j+1,n)=1;   %if point i,j =1, then the spaces immediately around it become 1
                B(i,j-1,n)=1;   %it's like the game of life.
                B(i,j,n)=1;
            end
            if (C(i,j)==1 || C(i,j)==3)
                for x=2:w-1
                    for y=2:h-1
                        if ((x-i)^2+(y-j)^2 < (n-1)^2)
                            B(x,y,n)=1; %in the case of an edge, it expands as a circle
                        else end
                    end
                end
            else end
        end
    end
    D(:,:,n-1)=B(:,:,n)-B(:,:,n-1); %we can find the difference between each iteration to find the "surface area"
    for i=1:w
        for j=1:h
            burnarea=outerlimit(i,j)*D(i,j,n-1)*(R*grainLength-2*(n-2))+2*outerlimit(i,j)-2*B(i,j,n);
            SA(n,2)=SA(n,2)+burnarea; %this ^^ big long equation takes into consideration the
        end%                          end burning effects of the grain. to make the ends inhibited change to burnarea=outerlimit(i,j)*D(i,j,n-1)*R
                                        
    end
    B_frame=B(:,:,n);
    if (1-isequal(B_frame.*outerlimit,B_frame) & (activation==0))
        edgeIndex=n-1;
        activation=1;
    end
        
    if SA(n,2)<=0   %this is basically the break condition
        k=n;    %when the surface area drops below 0 the sim stops
        break
        break
    else
        k=N;
    end
end

if runAnimation
    close(wb)
    wb=waitbar(1,'Complete');
end
%clf;
SA(:,2)=smooth(SA(:,2),6)/(2*R^2*pi()*grainLength); %SA is smoothed and nondimensionalized in terms of R
for n=1:k+3
    SA(n,1)=(n-1)/R;
end
SA=SA(1:k+3,:);
if runAnimation==1
    contour=zeros(w,h);
    subplot(1,2,1)
    title('burn animation','FontSize',14);
    
    for n=1:k-1
        subplot(1,2,1)
        image([0 w],[0 h],100*D(:,:,n)+15*outerlimit+30*B(:,:,n));  %this just outputs the burn animation in different colors
        %f(n)=getframe; %uncomment this to export animation movie
        subplot(1,2,2)
        plot(SA(n,1),SA(n,2),'.','MarkerFaceColor',[0 0 0],'MarkerSize', 10);
        xlabel('non-dimensional time [t/(r/R)]','FontSize',14);
        ylabel('non-dimentional flux [A_s/(2*pi*R)]','FontSize',14);
        pause(0.033);
        %     g(n)=getframe; %uncomment this to export animation movie
        hold on
        if (mod(n,3)==1)
            contour=contour+D(:,:,n);   %contour is just D but it skips every 3 so it is like a topo map
        end
    end
title('Thrust Curve','FontSize',14);
subplot(1,2,1)
pause(1);
image([0 w],[0 h],30*contour+15*outerlimit); 
title('Contour Plot','FontSize',14);
pause(.01)
%f(k)=getframe; %uncomment this to export animation movie
subplot(1,2,2)
plot(SA(:,1),SA(:,2));
%movie2avi(f,'BurnAnimation') %uncomment this to export animation movie
%movie2avi(g,'BurnCurve') %uncomment this to export animation movie
end
if runAnimation
close(wb);
end
packing_Efficiency=2*sum(SA(:,2))*(SA(2,1)-SA(1,1));
sliver_Time=(k-edgeIndex)/R
if (edgeIndex==1)
    sliver_Time=0;
end
F_Nom_stD=std(SA(5:k-10,2));
Graham_Score = packing_Efficiency*(1-F_Nom_stD/mean(SA(3:k-10,2)))^2
if runAnimation
    packing_Efficiency
end