function [ output_args ] = run(file,moviefile)
%run(file,moviefile)
%  The main file
%
%     file: The input file to use
%moviefile: The name of the file to save the AVI to.
clc;
close all;
clear all;

disp('Starting program')
tic;
starttime=now;
disp('Program started at')
disp(datestr(starttime))

%Read in the data from the data file given
disp('Reading in input file')
[node,elem,mg,mp,M,gamma,ki,time,dt,solvertype,acc,mindt,maxdt,minel,maxel] = readdata('inputfile.dat');


% for kk=1:1:37
% scatter(node(kk,1),node(kk,2));
% hold on;
% end

% 
% for i=1:length(elem(:,1))
% %     for k=1:i
%         n1=elem(i,1); n2=elem(i,2);
%         x11=[node(n1,1),node(n2,1)];
%         y11=[node(n1,2),node(n2,2)];
%         if elem(i,3)==1
%         plot(x11,y11,'k.-')
%         else
%           plot(x11,y11,'r.-')
%         end
%         hold on
% end



%Find the number of nodes
nnode=size(node,1);

%Initialize the vectors used to plot the shape
x={};
y={};

nel=size(elem,1);

%Set the counter for the movie frames
counter=1;

%Save the initial nodal positions
x{counter}=node(1:nnode,1);
%x{counter}(nnode+1)=node(1,1);
y{counter}=node(1:nnode,2);
%y{counter}(nnode+1)=node(1,2);


%Find the inital min and max for printing
xmin=min(x{counter});
xmax=max(x{counter});
ymin=min(y{counter});
ymax=max(y{counter});

%The initial time
t=0;

%Flag used later
flag=0;

%Set the first delta time to infinity to ensure it does not affect later
%calculation.
savedt(counter)=inf;

%Save the time for printing later
savet(counter)=0;

while and(t<time,flag~=1)
    
    %Increment the counter
    counter=counter+1;
    disp('********************************************************************')
    disp('Starting step number');
    disp(counter);
    
    if counter==1
        a=1;
    end
    
    %Check the elements to make sure they are not either too short or too
    %long
    disp('Checking the elements and adjusting sizes');
    [node,elem,flag]=check(node,elem,minel,maxel);
    %% add a loop to detect the fake element???????? markers
    
    
    % 
for i=1:length(elem(:,1))
%     for k=1:i
        n1=elem(i,1); n2=elem(i,2);
        x11=[node(n1,1),node(n2,1)];
        y11=[node(n1,2),node(n2,2)];
        if elem(i,3)==1
        plot(x11,y11,'k.-')
        else
          plot(x11,y11,'r.-')
        end
        hold on
end

    hold off;
    
    
    %If there are less than three elements, it is not possible to have an
    %enclosed volume
    if flag==1
        disp('There are less than three elements. This is no longer an enclosed surface.')
        disp('Results up to this point will be displayed')
        disp('The final time is:')
        disp(t)
        input('Press enter to continue')
        counter=counter-1;
    end
	
    if flag~=1;
        
        %Apply the specified solver to find the new position
        disp('Finding the new nodal positions');
        if solvertype==2 %Euler's method
            [vnodal]=fea(node,elem,mg,mp,M,gamma,ki);
            node=node+vnodal*dt;
            delt=dt;
        elseif solvertype==1 %Fixed Runge-Kutta
            [vnodal]=fea(node,elem,mg,mp,M,gamma,ki);
            k1=vnodal*dt;
            [vnodal]=fea(node+k1/2,elem,mg,mp,M,gamma,ki);
            k2=vnodal*dt;    
            [vnodal]=fea(node+k2/2,elem,mg,mp,M,gamma,ki);
            k3=vnodal*dt;    
            [vnodal]=fea(node+k3,elem,mg,mp,M,gamma,ki);
            k4=vnodal*dt;
            node=node+(1/6)*(k1+2*k2+2*k3+k4);
            delt=dt;
        else %Use variable Runge-Kutta by default
            [node,delt]=varRK(node,elem,mg,mp,M,gamma,ki,dt,acc,maxdt);
        end
        
        %Check to see is the presribed minimum time step is meet.
        if delt<mindt
            disp('The minimum time step was not achieved.')
            disp('Results up to this point will be displayed')
            disp('The final time is:')
            disp(t)
            input('Press enter to continue')
            counter=counter-1;
            flag=1;
        else
            savedt(counter)=delt;
            t=t+delt;
        end
        
    end
    
    %Only do the following if the current time, t, is under the final time
    %and the minimum step time was achieved
    if and(t<=time,flag~=1)

        disp('Finished calculating the following time:');
        disp(t);
        disp('Step size was:')
        disp(delt);
	
        nnode=size(node,1);
        
        %Save the nodal positions for printing later
        x{counter}=node(1:nnode,1);
        %x{counter}(nnode+1)=node(1,1);
        y{counter}=node(1:nnode,2);
        %y{counter}(nnode+1)=node(1,2);
         element{counter}=elem;      
        %Find the maximum to set the printout later
        if max(x{counter})>xmax
            xmax=max(x{counter});
        end
        if min(x{counter})<xmin
            xmin=min(x{counter});
        end
        if max(y{counter})>ymax
            ymax=max(y{counter});
        end
        if min(y{counter})<ymin
            ymin=min(y{counter});
        end
        
        savet(counter)=t;
        
        savefile=['moviefile','.mat'];
        save savefile xmin xmax ymin ymax x y savedt savet element;
        
    elseif and(t>time,flag~=1)
        disp('This time goes over the prescribed limiting time. This step will be ignored.');
    end
end

disp('Simulation complete. Data will now be displayed and saved')
disp('Simulation started at')
disp(datestr(starttime))
disp('Simulation ended at')
disp(datestr(now))
disp('The total running time was (in seconds):')
disp(toc)

%Determine the maximum and minimum of the x and y values. This is used for
%printing.
xmax=max(xmax,ymax);
xmin=min(xmin,ymin);
ymax=xmax;
ymin=xmin;

%Pad the sides of the plot
xsize=xmax-xmin;
ysize=ymax-ymin;
xmin=xmin-xsize/10;
xmax=xmax+xsize/10;
ymin=ymin-ysize/10;
ymax=ymax+ysize/10;

%Save the file in case something goes wrong with the AVI creation, the data
%is not lost.

disp('Saving the data')
save savefile xmin xmax ymin ymax x y savedt savet;

%Plot out the nodal positions from the simulation
close all;
figure;
format short;
% plot(x{1},y{1},'k.-')
% axis([xmin xmax ymin ymax]);
% printtitle=['Time is: ',num2str(savet(1),'%0.7g'),' sec'];
% text(xmin+xsize/20,ymax-ysize/20,printtitle);
% frames(1)=getframe;
% for i=2:size(x,2)
%     plot(x{i},y{i},'k.-')
%     axis([xmin xmax ymin ymax]);
%     printtitle=['Time is: ',num2str(savet(i),'%0.7g'),' sec'];
%     text(xmin+xsize/20,ymax-ysize/20,printtitle);
%     frames(i)=getframe;
% end


for i=2:size(x,2)
    for k=1:length(element{1,i}(:,1))
        n1=element{1,i}(k,1); n2=element{1,i}(k,2);
        x11=[x{i}(n1),x{i}(n2)];
        y11=[y{i}(n1),y{i}(n2)];
        
         if element{1,i}(k,3)==1
        plot(x11,y11,'k.-')
        else
          plot(x11,y11,'r.-')
        end
        
        hold on
    end
    xlim([-6 6]);
ylim([4 16]);
    hold off
%     axis([xmin xmax ymin ymax]);
    printtitle=['Time is: ',num2str(savet(i),'%0.7g'),' sec'];
    text(xmin+xsize/20,ymax-ysize/20,printtitle);
    frames(i-1)=getframe;
end

%Write the AVI file
movie2avi(frames,'moviefile.avi','quality',90);
%aviobj = avifile('moviefile.avi','compression','None');

