%Copyright (C) 2004  Mihai Pruna

%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.

%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%http://www.gnu.org/licenses/gpl.html

%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.











%initializing global variables, freestream conditions
global difference dummyro dummyu dummyv dummyEt dummye dummyT dummyp
dummyro=0;
dummyT=0;
dummyu=0;
global oldro
oldro=zeros(72,72)
global deltatcfl
deltatcfl=zeros(72,72);
global tauxx
tauxx=zeros(72,72);
global tauxy
tauxy=zeros(72,72);
global tauyy
tauyy=zeros(72,72);
global qx
qx=zeros(72,72);
global qy
qy=zeros(72,72);
global p
p=zeros(72,72);
global ro
ro=zeros(72,72);
global u
u=zeros(72,72);
global v
v=zeros(72,72);
global T
T=zeros(72,72);
global e
e=zeros(72,72);
global Et
Et=zeros(72,72);
global k
k=zeros(72,72);
global U1
U1=zeros(72,72);
global U2
U2=zeros(72,72);
global U3
U3=zeros(72,72);
global U5
U5=zeros(72,72);
global E1
E1=zeros(72,72);
global E2
E2=zeros(72,72);
global E3
E3=zeros(72,72);
global E5
E5=zeros(72,72);
global F1
F1=zeros(72,72);
global F2
F2=zeros(72,72);
global F3
F3=zeros(72,72);
global F5
F5=zeros(72,72);
global mew
mew=zeros(72,72);

global deltat
global IMIN
IMIN=1;
global JMIN
JMIN=1;
global IMAX
IMAX=70;
global JMAX
JMAX=70;
global MAXIT
MAXIT=3000;
global Minf
Minf=4;
global ainf
ainf=340.28;
global pinf
pinf=101325;
global Tinf
Tinf=288.16;
timesave=zeros(100,1);
global Tw
Tw=288.16;
global gamma
gamma=1.4;
global Pr
Pr=0.71;
global Vinf
Vinf=Minf*ainf;
global uinf
uinf=Vinf;
global vinf
vinf=0;
global meaw0
meaw0=0.000017894;
global T0
T0=288.16;
global R
R=287;
global LHORI
LHORI=0.00001;
global cv
cv=R./(gamma-1);
global cp
cp=gamma.*cv;
global roinf
roinf=pinf./(R.*Tinf);
%here we call dynvis to get meawinf
global meawinf
meawinf=dynvis(Tinf);

global Rel
Rel=roinf.*Vinf.*LHORI./meawinf;
global einf
einf=cv.*Tinf;
global delta
delta=5*LHORI./(Rel)^(0.5);
global LVERT
LVERT=5*delta;
global deltax
deltax=LHORI./(IMAX-1);
global deltay
deltay=LVERT./(JMAX-1);
global kinf
kinf=thermc(meawinf);
global lambda
%initializing deltat with a small value. deltat will be recalculated each
%iteration
deltat=1*10^(-11)
%iterating with respect to time
for ITER=1:MAXIT
        %for the first iteration, the initial flowfield conditions are put
        %in the grid
        
    if (ITER==1)
    
        for I=IMIN:IMAX
        for J=JMIN:JMAX
      u(I,J)=uinf;
      v(I,J)=vinf;
      p(I,J)=pinf;
      T(I,J)=Tinf;
        end
    
        end
    
        for I=IMIN:IMAX
        for J=JMIN:JMAX
            if (I==IMIN) & (J==JMIN)
            u(I,J)=0;
            v(I,J)=0;
            elseif (J==JMIN)
            u(I,J)=0;
            v(I,J)=0;
            p(I,J)=2*p(I,2)-p(I,3);
            T(I,J)=Tw;
            
            elseif (I==IMAX)
            u(I,J)=2*u(I-1,J)-u(I-2,J);        
            v(I,J)=2*v(I-1,J)-v(I-2,J);
            p(I,J)=2*p(I-1,J)-p(I-2,J);
            T(I,J)=2*T(I-1,J)-T(I-2,J);
            end
        
        end
    end    
    end 
    
    
    %calling the maccormak subroutine
   somanydummyvars=mac(1);
%recalculating the time step each iteration
   deltat=tstep(1);
%checking for convergence
   checkforconvergence=conver(1);
ITER;
if (checkforconvergence==1) break
end
    oldro=ro;
    
    

end
%checking for mass continuity upon convergence
verifx=MDOT(somanydummyvars);
if (verifx==0) fprintf('\nMass flow discontinuity!!!\n')
end
%plots, and generating additional variables used in plots
I=1:70;
figure(1)
plot((I-1)*deltax,p(I,1)/pinf)
for I=2:70
for J=1:70
ybar(I,J)=(J-1)*deltay./((I-1)*deltax)*sqrt(ro(I,J)*sqrt(u(I,J)^2+v(I,J)^2)*(I-1)*deltax./mew(I,J));
end
end
J=1:70;
figure(2)
plot(p(35,J)/pinf,ybar(35,J))
figure(3)
plot(p(9,J)/pinf,ybar(9,J))
figure(4)
plot(p(70,J)/pinf,ybar(70,J))
figure(5)
plot(T(9,J)/Tinf,ybar(9,J))
figure(6)
plot(T(9,J)/Tinf,ybar(9,J))
figure(7)
plot(T(35,J)/Tinf,ybar(35,J))
figure(8)
plot(u(35,J)/uinf,ybar(35,J))
figure(9)
plot(u(70,J)/uinf,ybar(70,J))
figure(10)
plot(u(9,J)/uinf,ybar(9,J))
for I=1:70
for J=1:70
Vmag(I,J)=sqrt(u(I,J)^2+v(I,J)^2);
end
end
J=1:70
figure(11)
plot(Vmag(9,J)/Vinf,ybar(9,J))
figure(12)
plot(Vmag(35,J)/Vinf,ybar(35,J))
figure(13)
plot(Vmag(70,J)/Vinf,ybar(70,J))
figure(14)
plot(v(70,J),ybar(70,J))
figure(15)
plot(v(35,J),ybar(35,J))
figure(16)
plot(v(9,J),ybar(9,J))
for I=1:70
for J=1:70
Mach(I,J)=Vmag(I,J)/sqrt(gamma.*R.*T(I,J));
end
end
J=1:70
figure(17)
plot(Mach(70,J),ybar(70,J))
figure(18)
plot(Mach(35,J),ybar(35,J))
figure(19)
plot(Mach(9,J),ybar(9,J))
figure(20)
plot(ybar(9,J),p(9,J)/pinf,ybar(35,J),p(35,J)/pinf,ybar(70,J),p(70,J)/pinf) 
figure(21)
plot(ybar(9,J),T(9,J)/Tinf,ybar(35,J),T(35,J)/Tinf,ybar(70,J),T(70,J)/Tinf) 
figure(22)
plot(ybar(9,J),Mach(9,J),ybar(35,J),Mach(35,J),ybar(70,J),Mach(70,J)) 
figure(23)
plot(ybar(9,J),u(9,J)/uinf,ybar(35,J),u(35,J)/uinf,ybar(70,J),u(70,J)/uinf)
figure(24)
plot(ybar(9,J),v(9,J),ybar(35,J),v(35,J),ybar(70,J),v(70,J))


%functions, put them in separate files before compiling!

%bc.m

function boundcond=bc(anotherdummyvar)
%computes boundary conditions(the 4 cases described in the book)
global ro u v T Tw e Et meaw0 mew k p IMIN IMAX JMIN JMAX roinf pinf uinf vinf Tinf meawinf kinf cv R
for I=IMIN:IMAX
    ro(I,JMAX)=roinf;
    u(I,JMAX)=uinf;
    v(I,JMAX)=0;
    p(I,JMAX)=pinf;
    T(I,JMAX)=Tinf;
    mew(I,JMAX)=dynvis(T(I,JMAX));
    k(I,JMAX)=thermc(mew(I,JMAX));
    e(I,JMAX)=cv.*T(I,JMAX);
    Et(I,JMAX)=ro(I,JMAX)*(e(I,JMAX)+(u(I,JMAX)^2+v(I,JMAX)^2)/2);
end
 for J=JMIN+1:JMAX-1
    ro(1,J)=roinf;
    u(1,J)=uinf;
    v(1,J)=0;
    p(1,J)=pinf;
    T(1,J)=Tinf;
    mew(1,J)=dynvis(T(1,J));
    k(1,J)=thermc(mew(1,J));
    e(1,J)=cv.*T(1,J);
    Et(1,J)=ro(1,J)*(e(1,J)+(u(1,J)^2+v(1,J)^2)/2);
 end   


    ro(1,1)=roinf;
 u(1,1)=0;
    v(1,1)=0;
    p(1,1)=pinf;
    T(1,1)=Tinf;
    mew(1,1)=meaw0;
    k(1,1)=thermc(meaw0);
    e(1,1)=cv.*T(1,1);
    Et(1,1)=ro(1,1)*(e(1,1)+(u(1,1)^2+v(1,1)^2)/2);
    
    for I=IMIN:IMAX
   
    u(I,1)=0;
    v(I,1)=0;
    p(I,1)=2*p(I,2)-p(I,3);
    T(I,1)=Tw;
     ro(I,1)=p(I,1)/(R.*T(I,1));
    
     mew(I,1)=meaw0;
    k(I,1)=thermc(mew(I,1));
    e(I,1)=cv.*T(I,1);
    Et(I,1)=ro(I,1)*(e(I,1)+(u(I,1)^2+v(I,1)^2)/2);
    end
    
    for J=JMIN+1:JMAX-1
        u(IMAX,J)=2*u(IMAX-1,J)-u(IMAX-2,J);
        v(IMAX,J)=2*v(IMAX-1,J)-v(IMAX-2,J);
        p(IMAX,J)=2*p(IMAX-1,J)-p(IMAX-2,J);
        T(IMAX,J)=2*T(IMAX-1,J)-T(IMAX-2,J);
        ro(IMAX,J)=p(IMAX,J)/(R.*T(IMAX,J));
        mew(IMAX,J)=dynvis(T(IMAX,J));
        k(IMAX,J)=thermc(mew(IMAX,J));
    e(IMAX,J)=cv.*T(IMAX,J);
    Et(IMAX,J)=ro(IMAX,J)*(e(IMAX,J)+(u(IMAX,J)^2+v(IMAX,J)^2)/2);
    end
    
    boundcond=0;

%------EOF-------------------
%calcqx.m

function theqx=calcqx(i,j,CASE)
%calculates qx(2 cases)
global IMIN IMAX JMIN JMAX T deltax deltay k
%predictor E
if (CASE==1)
    if (i==IMIN)
    dTdx=(T(i+1,j)-T(i,j))/deltax;
else 
    dTdx=(T(i,j)-T(i-1,j))/deltax;
end
%corrector E
elseif (CASE==2)
    if (i==IMAX)
        dTdx=(T(i,j)-T(i-1,j))/deltax;
    else
        dTdx=(T(i+1,j)-T(i,j))/deltax;
    end
end
theqx=-k(i,j)*dTdx;

%------EOF-------------------    

%calcqy.m

function theqy=calcqy(i,j,CASE)
%calculates qy (2 cases)
global IMIN IMAX JMIN JMAX T deltax deltay k
%predictor F
if (CASE==1)
    if (j==JMIN)
    dTdy=(T(i,j+1)-T(i,j))/deltay;
else 
    dTdy=(T(i,j)-T(i,j-1))/deltay;
end
%corrector F
elseif (CASE==2)
    if (j==JMAX)
        dTdy=(T(i,j)-T(i,j-1))/deltay;
    else
        dTdy=(T(i,j+1)-T(i,j))/deltay;
    end
end
theqy=-k(i,j)*dTdy;

%------EOF-------------------
%calctauxx.m

function thetauxx=calctauxx(i,j,CASE)
%calculates tauxx (2 cases)
global u v mew deltax deltay IMIN IMAX JMIN JMAX lambda
%predictor E
if (CASE==1)
    if (i==IMIN)
        dudx=(u(i+1,j)-u(i,j))/deltax;
    else dudx=(u(i,j)-u(i-1,j))/deltax;
    end
    
    %corrector E
elseif (CASE==2)
    if(i==IMAX)
        dudx=(u(i,j)-u(i-1,j))/deltax;
    else dudx=(u(i+1,j)-u(i,j))/deltax;
    end
end
    if (j==JMIN)
        dvdy=(v(i,j+1)-v(i,j))/deltay;
    elseif  (j==JMAX)
        dvdy=(v(i,j)-v(i,j-1))/deltay;
    else
        dvdy=(v(i,j+1)-v(i,j-1))/(2*deltay);
    end
    lambda=-2/3*mew(i,j);
    thetauxx=lambda.*(dudx+dvdy)+2*mew(i,j)*dudx;

%------EOF-------------------
%calctauxy.m

function thetauxy=calctauxy(i,j,CASE)
%calculates tauxy (4 cases)
global IMAX IMIN JMAX JMIN u v deltax deltay deltat mew
%predictor E
if (CASE==1) 
    if (j==JMIN) dudy=(u(i,j+1)-u(i,j))/deltay; 
    elseif (j==JMAX) dudy=(u(i,j)-u(i,j-1))/deltay;
    else dudy=(u(i,j+1)-u(i,j-1))/(2*deltay);
    end
    if (i==IMIN) dvdx=(v(i+1,j)-v(i,j))/deltax;
    else dvdx=(v(i,j)-v(i-1,j))/deltax; 
    end
    %corrector E
elseif (CASE==2)
    if (j==JMIN) dudy=(u(i,j+1)-u(i,j))/deltay; 
    elseif (j==JMAX) dudy=(u(i,j)-u(i,j-1))/deltay;
    else dudy=(u(i,j+1)-u(i,j-1))/(2*deltay);
    end
    if (i==IMAX) dvdx=(v(i,j)-v(i-1,j))/deltax;
    else dvdx=(v(i+1,j)-v(i,j))/deltax; 
    end
%predictor F    
elseif (CASE==3)
    if (i==IMIN) dvdx=(v(i+1,j)-v(i,j))/deltax; 
    elseif (i==IMAX) dvdx=(v(i,j)-v(i-1,j))/deltax;
    else dvdx=(v(i+1,j)-v(i-1,j))/(2*deltax);
    end
    if (j==JMIN) dudy=(u(i,j+1)-u(i,j))/deltay;
    else dudy=(u(i,j)-u(i,j-1))/deltay; 
    end
 %corrector F
 elseif (CASE==4)
    if (i==IMIN) dvdx=(v(i+1,j)-v(i,j))/deltax; 
    elseif (i==IMAX) dvdx=(v(i,j)-v(i-1,j))/deltax;
    else dvdx=(v(i+1,j)-v(i-1,j))/(2*deltax);
    end
    if (j==JMAX) dudy=(u(i,j)-u(i,j-1))/deltay;
    else dudy=(u(i,j+1)-u(i,j))/deltay; 
    end
end
thetauxy=mew(i,j)*(dudy+dvdx);

%------EOF-------------------  
   
%calctauyy.m

function thetauyy=calctauyy(i,j,CASE)
%calculates tauyy (2 cases)
global u v mew deltax deltay IMIN IMAX JMIN JMAX lambda
%predictor F
if (CASE==1)
    if (j==JMIN)
        dvdy=(v(i,j+1)-v(i,j))/deltay;
    else dvdy=(v(i,j)-v(i,j-1))/deltay;
    end
   
    %corrector F
elseif (CASE==2)
    if(j==JMAX)
        dvdy=(v(i,j)-v(i,j-1))/deltay;
    else dvdy=(v(i,j+1)-v(i,j))/deltay;
    end
end
    if (i==IMIN)
        dudx=(u(i+1,j)-u(i,j))/deltax;
    elseif  (i==IMAX)
        dudx=(u(i,j)-u(i-1,j))/deltax;
    else
        dudx=(u(i+1,j)-u(i-1,j))/(2*deltax);
    end
    lambda=-2/3*mew(i,j);
    thetauyy=lambda.*(dudx+dvdy)+2*mew(i,j)*dvdy; 

%------EOF-------------------
%conver.m

function convergence=conver(doomvar)
%checks for convergence by checking that the change in density between the last
%two time steps at each grid point does not exceed 10^(-8)
global IMIN IMAX JMIN JMAX ro oldro
convergence=1;
for I=IMIN:IMAX
    for J=JMIN:JMAX
        if (abs(ro(I,J)-oldro(I,J))>10^(-8))
            convergence=0;
        end
    end
end

%------EOF-------------------   
%decode.m

function dec=decode(UU1,UU2,UU3,UU5)
%decodes the flux variables U into primitive variables
global cv R dummyro dummyu dummyv dummyEt dummye dummyT dummyp
dummyro=UU1;
dummyu=UU2./UU1;
dummyv=UU3./UU1;
dummyEt=UU5;
dummye=UU5./UU1-(dummyu.^2+dummyv.^2)/2;
dummyT=dummye./cv;
dummyp=dummyro.*R.*dummyT;
dec=0;

%------EOF-------------------
%dynvis.m

function dynviscosity= dynvis(temp)          % Subfunction dynvis
% Calculates dynamic viscosity
global meaw0
global T0
dynviscosity=meaw0.*(temp./T0)^(1.5)*(T0+110)/(temp+110);

%------EOF-------------------
%mac.m

function maccormack=mac(dummyvar)
%the maccormack subroutine
global Pr U1 U2 U3 U5 E1 E2 E3 E5 F1 F2 F3 F5 e u v p T ro Et tauxx tauxy tauyy R cv qx qy mew k IMIN IMAX JMIN JMAX deltat deltax deltay dummyro dummyu dummyv dummyEt dummye dummyT dummyp
U1bar=zeros(72,72);
U2bar=zeros(72,72);
U3bar=zeros(72,72);
U5bar=zeros(72,72);
dummmy=0;
%Us,Es and Fs are calculated at the current time at each grid point
for I=IMIN:IMAX
    for J=JMIN:JMAX
              ro(I,J)=p(I,J)/(R.*T(I,J));
        mew(I,J)=dynvis(T(I,J));
         k(I,J)=thermc(mew(I,J));
    e(I,J)=cv.*T(I,J);
    Et(I,J)=ro(I,J)*(e(I,J)+(u(I,J)^2+v(I,J)^2)/2);
       
        U1(I,J)=p(I,J)/(R.*T(I,J));
        U2(I,J)=U1(I,J)*u(I,J);
        U3(I,J)=U1(I,J)*v(I,J);
        U5(I,J)=U1(I,J)*(cv.*T(I,J)+(u(I,J)^2+v(I,J)^2)/2);
        
        tauxy(I,J)=calctauxy(I,J,1);
        
        tauxx(I,J)=calctauxx(I,J,1);
        qx(I,J)=calcqx(I,J,1);
        
        E1(I,J)=U2(I,J);
        E2(I,J)=E1(I,J)*u(I,J)+p(I,J)-tauxx(I,J);
        E3(I,J)=E1(I,J)*v(I,J)-tauxy(I,J);
        E5(I,J)=(U5(I,J)+p(I,J))*u(I,J)-u(I,J)*tauxx(I,J)-v(I,J)*tauxy(I,J)+qx(I,J);
        
        tauxy(I,J)=calctauxy(I,J,3);
        tauyy(I,J)=calctauyy(I,J,1);
        qy(I,J)=calcqy(I,J,1);
        
        F1(I,J)=U3(I,J);
        F2(I,J)=E1(I,J)*v(I,J)-tauxy(I,J);
        F3(I,J)=U3(I,J)*v(I,J)+p(I,J)-tauyy(I,J);
        F5(I,J)=(U5(I,J)+p(I,J))*v(I,J)-u(I,J)*tauxy(I,J)-v(I,J)*tauyy(I,J)+qy(I,J);
    end
end
%Us are predicted for interior points
for I=IMIN+1:IMAX-1
    for J=JMIN+1:JMAX-1

        U1bar(I,J)=U1(I,J)-deltat./deltax.*(E1(I+1,J)-E1(I,J))-deltat./deltay.*(F1(I,J+1)-F1(I,J));
        U2bar(I,J)=U2(I,J)-deltat./deltax.*(E2(I+1,J)-E2(I,J))-deltat./deltay.*(F2(I,J+1)-F2(I,J));
        U3bar(I,J)=U3(I,J)-deltat./deltax.*(E3(I+1,J)-E3(I,J))-deltat./deltay.*(F3(I,J+1)-F3(I,J));
        U5bar(I,J)=U5(I,J)-deltat./deltax.*(E5(I+1,J)-E5(I,J))-deltat./deltay.*(F5(I,J+1)-F5(I,J));

        deeee=decode(U1bar(I,J),U2bar(I,J),U3bar(I,J),U5bar(I,J));
        %primitive variables are recalculated
        ro(I,J)=dummyro;
        u(I,J)=dummyu;
        v(I,J)=dummyv;
        Et(I,J)=dummyEt;
        e(I,J)=dummye;
        T(I,J)=dummyT;
        p(I,J)=dummyp;
        mew(I,J)=dynvis(T(I,J));
        k(I,J)=thermc(mew(I,J));
    end
end
%boundary conditions are recomputed
dummybc=bc(1);
tauxx=zeros(72,72);
qx=zeros(72,72);
qy=zeros(72,72);
tauxy=zeros(72,72);
%Es and Fs are calculated with predicted values
for I=IMIN:IMAX
    for J=JMIN:JMAX
        
        tauxx(I,J)=calctauxx(I,J,2);
        
        tauxy(I,J)=calctauxy(I,J,2);
        qx(I,J)=calcqx(I,J,2);
        E1(I,J)=ro(I,J)*u(I,J);
        E2(I,J)=ro(I,J)*u(I,J)^2+p(I,J)-tauxx(I,J);
        E3(I,J)=ro(I,J)*u(I,J)*v(I,J)-tauxy(I,J);
        E5(I,J)=(Et(I,J)+p(I,J))*u(I,J)-u(I,J)*tauxx(I,J)-v(I,J)*tauxy(I,J)+qx(I,J);
        tauxy(I,J)=calctauxy(I,J,4);
        tauyy(I,J)=calctauyy(I,J,2);
        qy(I,J)=calcqy(I,J,2);
        F1(I,J)=ro(I,J)*v(I,J);
        F2(I,J)=ro(I,J)*u(I,J)*v(I,J)-tauxy(I,J);
        F3(I,J)=ro(I,J)*v(I,J)^2+p(I,J)-tauyy(I,J);
        F5(I,J)=(Et(I,J)+p(I,J))*v(I,J)-u(I,J)*tauxy(I,J)-v(I,J)*tauyy(I,J)+qy(I,J);
    end 
end
%Us are corrected
for I=IMIN+1:IMAX-1
    for J=JMIN+1:JMAX-1
        Ubuffer=U1(I,J);
        U1(I,J)=1/2*(Ubuffer+U1bar(I,J)-deltat./deltax.*(E1(I,J)-E1(I-1,J))-deltat./deltay.*(F1(I,J)-F1(I,J-1)));
        Ubuffer=U2(I,J);
        U2(I,J)=1/2*(Ubuffer+U2bar(I,J)-deltat./deltax.*(E2(I,J)-E2(I-1,J))-deltat./deltay.*(F2(I,J)-F2(I,J-1)));
        Ubuffer=U3(I,J);
        U3(I,J)=1/2*(Ubuffer+U3bar(I,J)-deltat./deltax.*(E3(I,J)-E3(I-1,J))-deltat./deltay.*(F3(I,J)-F3(I,J-1)));
        Ubuffer=U5(I,J);
        U5(I,J)=1/2*(Ubuffer+U5bar(I,J)-deltat./deltax.*(E5(I,J)-E5(I-1,J))-deltat./deltay.*(F5(I,J)-F5(I,J-1)));
    end 
end
%primitive variables are recomputed
for I=IMIN+1:IMAX-1
    for J=JMIN+1:JMAX-1
morefreavar=decode(U1(I,J),U2(I,J),U3(I,J),U5(I,J));
        
        ro(I,J)=dummyro;
        u(I,J)=dummyu;
        v(I,J)=dummyv;
        Et(I,J)=dummyEt;
        e(I,J)=dummye;
        T(I,J)=dummyT;
        p(I,J)=dummyp;
        
        mew(I,J)=dynvis(T(I,J));
        
        k(I,J)=thermc(mew(I,J));
        
    end 
end
   dummybc=bc(1);
   maccormack=0;
        

%------EOF-------------------
% MDOT.m
function verify=MDOT(somanydummyvarsss)
%verifies mass continuity
global difference IMIN IMAX JMIN JMAX ro u
I=IMIN
MDOTL=0;
MDOTR=0;
for J=JMIN:JMAX
    MDOTL=MDOTL+ro(I,J)*u(I,J);
end

I=IMAX
for J=JMIN:JMAX
    MDOTR=MDOTR+ro(I,J)*u(I,J);
end

verify=1;
difference=abs(MDOTL-MDOTR)/(MDOTL)*100;
if (difference>1) verify=0;
end

%------EOF-------------------
%thermc.m

function kay=thermc(meaw)
%calculates thermal coefficient
global Pr cp
kay=meaw.*cp./Pr;

%------EOF-------------------

%tstep.m
function timestep=tstep(dummyvar)
%calculates time step
global Pr u v R gamma ro mew T IMIN IMAX JMIN JMAX deltax deltay 




deltatcfl=zeros(IMIN-2,IMAX-2);
for I=IMIN+1:IMAX-1
    for J=JMIN+1:JMAX-1
        deltatcfl(I-1,J-1)=0.5*(abs(u(I,J))/deltax + abs(v(I,J))/deltay + sqrt(gamma.*R.*T(I,J))*sqrt(1/deltax.^2 + 1/deltay.^2))^(-1);
    end
end

timestep=min(min(deltatcfl));



%------EOF-------------------