global rou T V A;
time=MC_TM_1D_C(61,0.015,1600,0.6784,0.1)

function time=MC_TM_1D_C(N,dt,t,pN,Cx)
% MACCORMACK`S TECHNIQUE
% Conservation Form
% Time Marching

tic;
initial_data(t,N);
for h=1:t
    preCal(h,N,Cx);
    predictor_step(N,h,dt);
    BoundaryConditions(h+1,N,pN);
    preCal(h+1,N,Cx);
    corrector_step(N,h,dt);
    BoundaryConditions(h+1,N,pN);
    Cal_rVT(h);
end
time=toc;

function preCal(h,N,Cx)
global U1 U2 U3 F1 F2 F3 J2 S1 S2 S3 A pA_px gama;
 F1(h,:)=U2(h,:);
 F2(h,:)=U2(h,:).^2./U1(h,:)+(gama-1)./(gama).*(U3(h,:)-(gama./2).*(U2(h,:).^2./U1(h,:)));
 F3(h,:)=gama.*(U2(h,:).*U3(h,:)./U1(h,:))-(gama.*(gama-1)./2).*(U2(h,:).^3./U1(h,:).^2);
 p(h,:)=(gama-1).*(U3(h,:)-(gama./2).*(U2(h,:).^2./U1(h,:)))./A;
 J2(h,:)=p(h,:).*pA_px./gama;
 for i=2:(N-1)
     temp=Cx*abs(p(h,i+1)-2*p(h,i)+p(h,i-1))/(p(h,i+1)+2*p(h,i)+p(h,i-1));
     S1(h,i)=temp*(U1(h,i+1)-2*U1(h,i)+U1(h,i-1));
     S2(h,i)=temp*(U2(h,i+1)-2*U2(h,i)+U2(h,i-1));
     S3(h,i)=temp*(U3(h,i+1)-2*U3(h,i)+U3(h,i-1));
 end

function predictor_step(N,h,dt);
global U1 U2 U3 F1 F2 F3 J2 S1 S2 S3 dx pU1_pt pU2_pt pU3_pt;
 for i=2:(N-1)
     pU1_pt(h,i)=-((F1(h,i+1)-F1(h,i))/dx);
     pU2_pt(h,i)=-((F2(h,i+1)-F2(h,i))/dx)+J2(h,i);
     pU3_pt(h,i)=-((F3(h,i+1)-F3(h,i))/dx);
 end
 U1(h+1,:)=U1(h,:)+pU1_pt(h,:).*dt+S1(h,:);
 U2(h+1,:)=U2(h,:)+pU2_pt(h,:).*dt+S2(h,:);
 U3(h+1,:)=U3(h,:)+pU3_pt(h,:).*dt+S3(h,:);

function BoundaryConditions(h,N,pN)
global U1 U2 U3 A gama;
 U1(h,1)=A(1);
 U2(h,1)=2*U2(h,2)-U2(h,3);
 U3(h,1)=U1(h,1)*(1/(gama-1)+(gama/2)*((U2(h,1)/U1(h,1))^2));

 U1(h,N)=2*U1(h,N-1)-U1(h,N-2);
 U2(h,N)=2*U2(h,N-1)-U2(h,N-2);
 VN=U2(h,N)/U1(h,N);
 U3(h,N)=pN*A(N)/(gama-1)+(gama/2)*U2(h,N)*VN;

function corrector_step(N,h,dt);
global U1 U2 U3 F1 F2 F3 J2 S1 S2 S3 dx pU1_pt pU2_pt pU3_pt;
 for i=2:(N-1)
     pU1_pt(h+1,i)=-((F1(h+1,i)-F1(h+1,i-1))/dx);
     pU2_pt(h+1,i)=-((F2(h+1,i)-F2(h+1,i-1))/dx)+J2(h+1,i);
     pU3_pt(h+1,i)=-((F3(h+1,i)-F3(h+1,i-1))/dx);
 end
 
 pU1_pt(h,:)=0.5.*(pU1_pt(h+1,:)+pU1_pt(h,:));
 pU2_pt(h,:)=0.5.*(pU2_pt(h+1,:)+pU2_pt(h,:));
 pU3_pt(h,:)=0.5.*(pU3_pt(h+1,:)+pU3_pt(h,:));
   
 U1(h+1,:)=U1(h,:)+pU1_pt(h,:).*dt+S1(h+1,:);
 U2(h+1,:)=U2(h,:)+pU2_pt(h,:).*dt+S2(h+1,:);
 U3(h+1,:)=U3(h,:)+pU3_pt(h,:).*dt+S3(h+1,:);

function initial_data(t,N)
global rou T V U1 U2 U3 F1 F2 F3 J2 S1 S2 S3 x dx A pA_px pU1_pt pU2_pt pU3_pt gama;
gama=1.4;
tU2=0.59;
dx=3/(N-1);
x=0:dx:3;

A=NF_Shape(N);
pA_px=zeros(1,N);
for i=2:(N-1)
    pA_px(i)=(A(i+1)-A(i))/dx;
end

rou=zeros(t+1,N);
T=zeros(t+1,N);
V=zeros(t+1,N);
init_condition(N,tU2);
U1=zeros(t+1,N);
U2=zeros(t+1,N);
U3=zeros(t+1,N);
F1=zeros(t+1,N);
F2=zeros(t+1,N);
F3=zeros(t+1,N);
J2=zeros(t+1,N);
S1=zeros(t+1,N);
S2=zeros(t+1,N);
S3=zeros(t+1,N);
U1(1,:)=rou(1,:).*A;
U2(1,:)=U1(1,:).*V(1,:);
U3(1,:)=rou(1,:).*(T(1,:)./(gama-1)+(gama./2).*V(1,:).^2).*A;

pU1_pt=zeros(t,N);
pU2_pt=zeros(t,N);
pU3_pt=zeros(t,N);

function y=NF_Shape(N)
global x;
  y=1+2.2*(x-1.5).^2;
 
function init_condition(N,tU2)
global rou T V x A;
  N1=floor(N/6);
  N2=floor(N/2);
  N3=floor(N*2/3);
 
  rou(1,1:N1)=1;
  rou(1,N1+1:N2)=1-0.366.*(x(N1+1:N2)-0.5);
  rou(1,N2+1:N3)=0.634-0.702.*(x(N2+1:N3)-1.5);
  rou(1,N3+1:N)=0.5892-0.10228.*(x(N3+1:N)-2.1);
 
  T(1,1:N1)=1;
  T(1,N1+1:N2)=1-0.167.*(x(N1+1:N2)-0.5);
  T(1,N2+1:N3)=0.833-0.4908.*(x(N2+1:N3)-1.5);
  T(1,N3+1:N)=0.93968+0.0622.*(x(N3+1:N)-2.1);
 
  V(1,:)=tU2./(rou(1,:).*A);

function Cal_rVT(h);
global rou V T U1 U2 U3 A gama;
 rou(h+1,:)=U1(h+1,:)./A;
 V(h+1,:)=U2(h+1,:)./U1(h+1,:);
 T(h+1,:)=(gama-1).*(U3(h+1,:)./U1(h+1,:)-(gama/2).*V(h+1,:).^2);