NMAX=101;
fid = fopen('ogrid.inp.dat','r')
%l=fscanf(fid,'%s')
%l=fscanf(fid','%s',1)
isym=fscanf(fid,'%d',1)
nu=fscanf(fid,'%d',1)
nl=fscanf(fid,'%d',1)
imax=fscanf(fid,'%d',1)
jmax=fscanf(fid,'%d',1)
if(imax >81) end
if(isym>0)
  nl = nu;
end
n = nu + nl - 1;
%
% Allocate storage for x and y
%
x = zeros(1,100);
y = zeros(1,100);
r = zeros(1,n);
theta=zeros(1,n);
r1 = zeros(1,imax);
theta1=zeros(1,imax);
%
% Read the points on the upper surface
%
for l = nl:nl+nu-1
   a=fscanf(fid,'%f',1);
   b = fscanf(fid,'%f',1);
   x(l) = a;
   y(l) = b;
end
if isym == 0 
%
% If the airfoil is not symmetric, read lower side ordinates too..
%
    for l = 1:nl
      a=fscanf(fid, '%f',1);
      b = fscanf(fid, '%f', 1);
      x(nl+1-l) = a;
      y(nl+1-l) = b;
    end
else

    for l =1:nl 
       x(nl+1-l) = x(nl-1+l);
       y(nl+1-l) = - y(nl-1+l);
    end
 end
 xn = fscanf(fid,'%f',1)
 yn = fscanf(fid,'%f',1)
 fclose(fid);
 slop1 = atan2((y(1)-y(2)),(x(1)-x(2)));
      slop2 = - atan2((y(n)-y(n-1)),(x(n)-x(n-1)));             
      eps   = abs(slop1 + slop2);
      pi    = 3.14158;
      power = 2.0 - eps / pi;
      xt    = 0.5 * (x(1) + x(n));
      yt    = 0.5 * (y(1) + y(n));
      z1    = xn + i * yn
      z2    = xt + i * yt
      z7    = 0.5 * (z1+z2);
      u     = 1.0;
      v     = 0.0;
      angl  = 8.0 * atan(1.0);
      power = 1.0 / power;
      for m=1:n
         
         z3 = x(m) + j * y(m);
         z4 = (z3 - z2)/(z3-z1);
         z4 = z4 ^ power;
         z4=(1+j*0++z4)*z7/(1.-j*0.-z4);
         x4=real(z4);
         y4=imag(z4);
         r(m)= abs(z4);
%         fprintf('r(i) %g. \n',r(i));
         a = x4;
         b = y4;
         angl = angl + atan2((u*b-v*a),(u*a+v*b));
         u = a;
         v = b;
         theta(m) = angl;
   %      fprintf('r(m), theta(m) %g %g. \n',r(m),theta(m));
      end
%      r(0) = 0.5 * (r(1)+r(m));
%      theta(0) = 8. * atan(1.);
%      r(m+1) = r(1);
%      theta(m+1) = 0.;
      dtheta = 8. * atan(1.0)/(imax-3);
           
      
      for m = 3: imax - 2
         theta11 = (m-2) * dtheta;
         theta1(m) = theta11;
      end
%      theta1
r1 = interp1(theta,r,theta1);
r1(2) = 0.5 * (r1(3)+r1(imax-2));
r1(imax-1)= r1(2);
r1(imax) = r1(3);
r1(1) = r1(imax-2);
theta1(2) = 0.0;
theta1(imax -1) = 8. * atan(1.0);
theta1(1) = - dtheta;
theta1(imax) = theta1(imax-1) + dtheta;
xxx = zeros(1,imax);
zzz = zeros(1,imax);
dr = 0.99 /(jmax-2);
fid = fopen('GRID.DAT','w');
fprintf(fid,'%d %d \n',imax,jmax); 
xg=zeros(81,31);
zg=zeros(81,31);

for n = 1:jmax
   eta = (n-2) * dr;
   eta1 = 1. / (1.-eta);
for m = 1: imax
   x4 = r1(m)*eta1 * cos(theta1(m));
   y4 = r1(m) * eta1 * sin(theta1(m));
   z4 = x4 + j * y4;
   z4 = ((z4-z7)/(z4+z7))^ (1./power+j*0);
   z4 = (z2 - z4*z1)/(1.+0.*j-z4);
   x4 = real(z4);
   y4 = imag(z4);
   fprintf(fid,'%d %d %f %f\n',m,n,x4,y4); 
   xg(m,n) = x4;
   zg(m,n) = y4;
   if(n ==2)
    xxx(m) = x4;
    zzz(m) = y4;
   end
end
end
fclose(fid);
% Plot the airfoil for world to see..
xxx3=zeros(1,imax-2);
zzz3=zeros(1,imax-2);
for n = 2 : 20
   for m = 2:imax-1
      xxx3(m-1)=xg(m,n);
      zzz3(m-1)=zg(m,n);
   end
   hold on
   plot(xxx3,zzz3);
end
xxx2=zeros(1,19);
zzz2=zeros(1,19);
for m = 2 : imax-1
   for n = 2:20
      xxx2(n-1)=xg(m,n);
      zzz2(n-1)=zg(m,n);
   end
   hold on
   plot(xxx2,zzz2);
   end