clear all; clf; clc;

sd_id=hdfsd('start','./out/q00000.hdf','read');
sds_index=0;
sds_id=hdfsd('select',sd_id,sds_index);
start=[0 0]; stride=[1 1]; edge=[32 32];
[q0,status]=hdfsd('readdata',sds_id,start,stride,edge);
status = hdfsd('end',sd_id);
sd_id=hdfsd('start','./out/q00000.hdf','read');
sds_index=1;
sds_id=hdfsd('select',sd_id,sds_index);
start=[0 0]; stride=[1 1]; edge=[32 32];
[q1,status]=hdfsd('readdata',sds_id,start,stride,edge);
status = hdfsd('end',sd_id);

sd_id=hdfsd('start','./out/q00000.hdf','read');
sds_index=2;
sds_id=hdfsd('select',sd_id,sds_index);
start=[0 0]; stride=[1 1]; edge=[32 32];
[q2,status]=hdfsd('readdata',sds_id,start,stride,edge);
status = hdfsd('end',sd_id);

sd_id=hdfsd('start','./out/q00000.hdf','read');
sds_index=3;
sds_id=hdfsd('select',sd_id,sds_index);
start=[0 0]; stride=[1 1]; edge=[32 32];
[q3,status]=hdfsd('readdata',sds_id,start,stride,edge);
status = hdfsd('end',sd_id);

gamma=1.4; gamma1=gamma-1;

mx=32; my=32;
dt=0.001; dx=1./mx; dy=1./32; dtdx=dt/dx; dtdy=dt/dy;
qLF=zeros([mx my 4]); f=zeros([mx my 4]); g=zeros([mx my 4]);
qnew=zeros([mx my 4]); 
rho=zeros([mx my]); u=zeros([mx my]); v=zeros([mx my]);
p=zeros([mx my]); u2v2=zeros([mx my]);
qLF(:,:,1)=q0; qLF(:,:,2)=q1; qLF(:,:,3)=q2; qLF(:,:,4)=q3;

for n=2:6
  rho=qLF(:,:,1); invrho=1./rho;
  u=qLF(:,:,2).*invrho; v=qLF(:,:,3).*invrho;
  u2v2=u.^2+v.^2;
  p=gamma1*(qLF(:,:,4)-0.5*u2v2);  
  f(:,:,1)=qLF(:,:,2);
  f(:,:,2)=qLF(:,:,2).*u+p;
  f(:,:,3)=qLF(:,:,2).*v;
  f(:,:,4)=(qLF(:,:,4)+p).*u;
  g(:,:,1)=qLF(:,:,3);
  g(:,:,2)=qLF(:,:,3).*u;
  g(:,:,3)=qLF(:,:,3).*v+p;
  g(:,:,4)=(qLF(:,:,4)+p).*v;
  disp('Iter ')
  n-2
  contourf(qLF(n-1:mx-n+2,n-1:my-n+2,1),20);
  pause(0.5)
  for i=2:mx-1
      for j=2:my-1
          %qLF(i,j,n)=0.25*(qLF(i+1,j,n-1)+qLF(i-1,j,n-1)+qLF(i,j+1,n-1)+qLF(i,j-1,n-1));
          f1=0.5*(f(i,j,:)+f(i+1,j,:))-0.25*dx/dt*(qLF(i+1,j,:)-qLF(i,j,:));
          f0=0.5*(f(i-1,j,:)+f(i,j,:))-0.25*dx/dt*(qLF(i,j,:)-qLF(i-1,j,:));
          g1=0.5*(g(i,j,:)+g(i,j+1,:))-0.25*dy/dt*(qLF(i,j+1,:)-qLF(i,j,:));
          g0=0.5*(g(i,j-1,:)+g(i,j,:))-0.25*dy/dt*(qLF(i,j,:)-qLF(i,j-1,:));          
          qnew(i,j,:)=qLF(i,j,:)-dtdx*(f1-f0)-dtdy*(g1-g0);
      end
  end
  qLF=qnew;
end