aminf=input('Please enter the Mach number= ');
alpha=input('Please enter the angle of attack in degrees= ');
itmax = input('Please enter the number of iterations to solve continuity= ');
% Initialize all arrays to zero
x=zeros(81,31);
y=zeros(81,31);
phi=zeros(81,31);
res=zeros(85,31);
delphi=zeros(81,31);
rhoa1=zeros(81,31);
rhoa3=zeros(81,31);
rhoout=zeros(1,91);
rho=zeros(1,91);
ucon=zeros(1,91);
vcont=zeros(1,91);
msqloc=zeros(1,91);
rhobar=zeros(1,91);
rhohat=zeros(1,91);
rhouqj=zeros(1,91);
rhs=zeros(1,91);
aaa=zeros(1,91);
bbb = zeros(1,91);
ccc = zeros(1,91);
xxx = zeros(1,91);
cp=zeros(1,78);
gamma=zeros(1,31);

% Initialize some constants
pi = 4. * atan(1.);
ga = 1.4;
time = 0.0;

% Open the file and read the grid

fid=fopen('GRID.DAT','r');
imax = fscanf(fid,'%d',1)
jmax = fscanf(fid,'%d',1)
if(imax > 81)
   fprintnf('imax too large');
end
if(jmax > 31)
   fprintnf('jmax too large');
end
for j=1:jmax
   for n=1:imax
      i=imax+1-n;
      iii=fscanf(fid,'%d',1);
      iii=fscanf(fid,'%d',1);
      a = fscanf(fid,'%f',1);
      b = fscanf(fid,'%f',1);
      x(i,j)=a;
      y(i,j) = b;
   end
end
% Flow and calculation parameters
alpha = alpha * pi/180;
cosang=cos(alpha);
sinang=sin(alpha);
numsup=0;
% initialize everything to zero..
for  j=1:jmax
   gamma(j)= 0.0;
end
for j=1:jmax
   for i=1:imax
      phi(i,j)=0.0;
      delphi(i,j) = 0.0;
   end
end
lcount = 0;

% Start iterating..

for itn=1:itmax
   lcount = lcount + 1;
   if (lcount > 5) 
      lcount = 1;
   end
   maxdel=0.;
   maxmach=0.;
   maxres=0.;
   delt = 0.1* (1.0 / 0.1) ^ (0.25* (lcount-1));
          
   numsup=0;
   ga    = 1.4;
   
   % Compute residuals
   for j=2:jmax-1
      %
      %     compute metrics, flow properties and ucon at (i+1/2,j)
      %
      for i=2:imax-1
      xx     = x(i+1,j) - x(i,j);
      xy     = 0.25 * (x(i+1,j+1)-x(i+1,j-1)+x(i,j+1)-x(i,j-1));
      yx     = y(i+1,j) - y(i,j);
      yy     = 0.25 * (y(i+1,j+1)-y(i+1,j-1)+y(i,j+1)-y(i,j-1));
      yac1   = 1. / (xx * yy - xy * yx);
      ax     =   yy * yac1;
      ay     = - xy * yac1;
      bx     = - yx * yac1;
      by     =   xx * yac1;
      px       = phi(i+1,j)-phi(i,j);
      py       =.25*(phi(i+1,j+1)-phi(i+1,j-1)+phi(i,j+1)-phi(i,j-1));      
      if (j == 2)
        py  =  -(px*(ax*bx + ay*by)+ cosang*bx + sinang*by)/(bx*bx+by*by);
     end
     %   calculate the u's and v's
     u = px * ax + py * bx+cosang;
     v = px * ay + py * by+sinang;
     %   calculate density
     b=(1./(ga-1.));
     a=1. - u * u -v * v;
     c=(1+((ga-1.)/2)*(aminf*aminf)*a);
     rho(i)=c^b;
     %   calculate contravariant velocity over jacobian
     ucon(i)    = (u * ax + v * ay) / yac1;
     rhoa1(i,j) = (ax * ax + ay * ay) / yac1;
     %   calculate local mach number squared
     msqloc(i)=(u*u+v*v)/c;
     if (msqloc(i) < 0.) 
        fprintf('uh oh... msq<0');
        break;
     else
        if (sqrt(msqloc(i))> maxmach)
           maxmach=sqrt(msqloc(i));
           mmxi=i;
           mmxj=j;
        end
     end
     %   supersonic counter
     if (msqloc(i) > 1.)
        numsup=numsup+1;
     end
  end
  %   set ghost boundaries
  msqloc(1)=msqloc(imax-2);
  msqloc(imax)=msqloc(3);
  rho(1)=rho(imax-2);
  rho(imax)=rho(3);
  %   loop to calculate rhobar
  for i=2:imax-1
     if (ucon(i) > 0.) 
        mu=max(0.,1.-2./(msqloc(i-1)+msqloc(i)));
        rhobar(i)=rho(i)-mu*(rho(i)-rho(i-1));
     else
        mu=max(0.,1.-2./(msqloc(i)+msqloc(i+1)));
        rhobar(i)=rho(i)-mu*(rho(i)-rho(i+1));
     end
  end
  %   calculate quantities (rhobar*u/j) and (rhobar*a1)
  for  i=2:imax-1
     rhouqj(i)=rhobar(i)*ucon(i);
     rhoa1(i,j)=rhobar(i) * rhoa1(i,j);
  end
  %   set ghost boundaries
  rhouqj(1)=rhouqj(imax-2);
  rhoa1(1,j)=rhoa1(imax-2,j);
  for i=2:imax-1
     res(i,j)=rhouqj(i)-rhouqj(i-1);
  end
      
  if (j == 2) 
     for i=2:imax-1
        rhoout(i) = rho(i);
     end
        
     rhoout(imax) = rho(3);
     rhoout(1) = rho(imax-2);
  end
end
%   now repeat process going from bottom to top, calculating
%   the j+1/2 half-point quantities, then move over a line and repeat.
for i=2:imax-1
   for j=2:jmax-1
      %   calculation of the phi-derivatives wrt xsi and eta:
      px = .25*(phi(i+1,j+1)-phi(i-1,j+1)+phi(i+1,j)-phi(i-1,j));
      py = phi(i,j+1) - phi(i,j);
      xx = 0.25 * (x(i+1,j+1)-x(i-1,j+1)+x(i+1,j)-x(i-1,j));
      yx = 0.25 * (y(i+1,j+1)-y(i-1,j+1)+y(i+1,j)-y(i-1,j));
      xy = x(i,j+1) - x(i,j);
      yy = y(i,j+1) - y(i,j);
      yac1 = 1. / (xx * yy - xy * yx);
      ax   =   yy * yac1;
      ay   = - xy * yac1;
      bx   = - yx * yac1;
      by   =   xx * yac1;
      %   calculate the u's and v's
      u = px * ax+py * bx +cosang;
      v = px * ay+py * by+sinang;
      %   calculate density, called rhohat here
      b=(1./(ga-1.));
      a=(1.-u*u-v*v );
      c=(1+((ga-1.)/2)*(aminf^2.)*a);
      rhohat(j)=c^b;
      %   calculate contravariant velocities
      vcont(j)= rhohat(j) * ( u * bx + v * by) / yac1;
      rhoa3(i,j) = rhohat(j) * (bx * bx + by * by) / yac1;
   end
   %---------------------------------------------------------------------
   %   add the contribution to residual(i,j)
   %   (treat j=2 differently here)
   
   res(i,2)=res(i,2)+(2.*vcont(2)) ;
   
   if (abs(res(i,2)) > maxres)
      maxres=abs(res(i,2));
      rmxi=i;
      rmxj=2;
   end
   %   loop from j=3 to j=jmax-1
   for j=3:jmax-1
      res(i,j)=res(i,j)+ ( vcont(j) - vcont(j-1));
      if (abs(res(i,j)) > maxres)
         maxres=abs(res(i,j));
         rmxi=i;
         rmxj=j;
      end
   end
end
%

% begin solve
      
%     af2 scheme
%
maxdel = 0.0;
ga     = 1.4;
omega  = 1.5;
%
%     xi - sweep
%
%for i = 2:imax -1
 %  for j = 2: jmax
  %    delphi(i,j) = 0.0;
   %end
%end

for j1 = 2: jmax-1
   j=jmax+1-j1;
   %
   %     assemble tri-diagonal matrix coefficients
   %
   for i = 2 : imax - 1
      bbb(i)    = -1./delt - (rhoa1(i,j)+rhoa1(i-1,j));
      aaa(i)    = rhoa1(i-1,j);
      ccc(i)    = rhoa1(i,j);
      rhs(i)  = omega * res(i,j)- 1./delt * delphi(i,j+1) ;
   end
   
   %
   %     perform thomas algorithm
   %      call thomas(a,b,c,2,imax-1,rhs)
   
   istart = 2;
   iend = imax - 1;
   rhs(istart-1) = 0.0;
   rhs(iend+1)   = 0.0;
   xxx(istart-1)   = 0.0;
   aaa(istart)     = 0.0;
   ccc(iend)       = 0.0;
   
   for i = istart : iend
      dd      = bbb(i) - aaa(i) * xxx(i-1);
      xxx(i)    = ccc(i) / dd;
      rhs(i)  = (rhs(i) - aaa(i) * rhs(i-1)) / dd;
   end
   
   for  m = istart : iend
      i = istart + iend - m;
      rhs(i)  = rhs(i) - xxx(i) * rhs(i+1);
   end
   
   for i     = istart : iend
      delphi(i,j) =  rhs(i);
   end
end

%
%     eta sweep
%
for i = 2 : imax - 1
   for   j = 2 : jmax - 1
      bbb(j)      = -1./delt -(rhoa3(i,j)+rhoa3(i,j-1));              
      aaa(j)      = rhoa3(i,j-1) ;
      ccc(j)      = rhoa3(i,j) +1./delt;
      rhs(j)    = -(delphi(i,j+1) - delphi(i,j))/delt;
   end
   %  call thomas(a,b,c,2,jmax-1,rhs)
   istart = 2;
   iend = jmax - 1;
   rhs(istart-1) = 0.0;
   rhs(iend+1)   = 0.0;
   xxx(istart-1)   = 0.0;
   aaa(istart)     = 0.0;
   ccc(iend)       = 0.0;
   for m = istart : iend
      dd      = bbb(m) - aaa(m) * xxx(m-1);
      xxx(m)    = ccc(m) / dd;
      rhs(m)  = (rhs(m) - aaa(m) * rhs(m-1)) / dd;
   end
   for m = istart : iend 
      n = istart + iend - m;
      rhs(n)  = rhs(n) - xxx(n) * rhs(n+1);
   end
   for j=2:jmax-1
      delphi(i,j) = rhs(j);
      phi(i,j) = phi(i,j) + delphi(i,j);
      maxdel   =max(maxdel,abs(delphi(i,j)/delt));
   end
end

%
%fprintf ('maxdel= %g \n', maxdel);
gamma(1) = phi(imax-2,2) - phi(3,2);
gamma(2) = gamma(1);
%     update top boundary
for i=2:(imax-1)/2
   a=-atan2(((1.-aminf^2.)^.5)*y(i,jmax),x(i,jmax));
   phi(i,jmax)=(gamma(1)/(2.*pi))*a;
end

for i=(imax-1)/2+1:imax-1
   b=2*pi-atan2(((1.-aminf^2.)^.5)*y(i,jmax),x(i,jmax));
   phi(i,jmax)=(gamma(1)/(2.*pi))*b;
end

% update left and right boundaries
for j=2:jmax
   phi(1,j)=phi(imax-2,j)-gamma(1);
   phi(imax,j)=phi(3,j)+gamma(1);
end

%  do bottom boundary
for i=2:imax-1
   %  metrics at i,j:
   xxs=(x(i+1,2)-x(i-1,2))/2.;
   yxs=(y(i+1,2)-y(i-1,2))/2.;
   xet=(x(i,3)-x(i,1))/2.;
   yet=(y(i,3)-y(i,1))/2.;
   yacl=1./(xxs*yet-xet*yxs);
   xsx=yet*yacl;
   xsy=-xet*yacl;
   etx=-yxs*yacl;
   ety=xxs*yacl;
   
   %---------------------------------------------------------------------
   %  compute phi(i,1):
   %
   %--------------------------------------------------------------------
   phixs=(phi(i+1,2)-phi(i-1,2))/2.;
   at=xsx*etx+xsy*ety;
   bt=etx*cosang+ety*(sinang);
   ct=etx^2.+ety^2.;
   %-------ft=fcorr*yacl-------------------------------------------------
   phiet=-(phixs*at+bt)/ct;
   phi(i,1)=phi(i,3)-2.*phiet;
end
%--------------------------------------------------------------------
% final update of ghost boundaries:
phi(1,1)=phi(imax-2,1)-gamma(1);
phi(imax,1)=phi(3,1)+gamma(1);
time = time + delt;

% Print Cp at mid-points
% Also compute Cl, Cd
cy = 0.0;
cx = 0.0;
cm = 0.0;
xmid=zeros(1,78);
for i=2:imax-2
   xmid(i-1)=0.5*(x(i,2)+x(i+1,2));
   cp(i-1) = (rhoout(i) ^ 1.4 - 1.0)/(0.7*aminf^2);
   cy = cy - cp(i-1) * (x(i+1,2)-x(i,2));
   cx = cx + cp(i-1) * (y(i+1,2)-y(i,2));
   cm = cm - cp(i-1) * (xmid(i-1)-0.25) * (x(i+1,2)-x(i,2));
end   
cl = cy * cos(alpha) - cx * sin(alpha);
cd = cx *cos(alpha) + cy * sin(alpha);
if(lcount==1)
   fprintf('Iteration=%4d Error in continuity=%14.6g Cl=%14.6g Cd=%14.6g Cm=%14.6g \n', itn,maxres,cl,cd,cm);
end
end

% Write Cp values to a file

fid=fopen('CP.DAT','w');
fprintf(fid,'           X       Y    CP \n');
for i=3:imax-2
      cpnode = 0.5 *(cp(i-1)+cp(i-2));
      fprintf(fid,'%14.7g %14.7g %14.7g \n',x(i,2),y(i,2),cpnode);
end
fprintf(fid,'Cl= %g \n', cl);
fprintf(fid,'Cd= %g \n', cd);
fprintf(fid,'Cm c/4= %g \n', cm);
fclose(fid);
plot(xmid,cp);