PROGRAM MAIN
      DIMENSION X(101) , Y(101) , R(101) , THETA(101),
     1          R1(101), THETA1(101)
      DIMENSION XG(81,31),YG(81,31)
      COMPLEX Z1,Z2,Z3,Z4,Z5,Z6,Z7
C
C     ALGEBRAIC O- GRID GENERATOR BASED ON 
C     KARMAN-TREFFZ MAPPING
C
C     ******************************************************
C
C     EXPALNATION OF INPUT QUANTITIES
C
C     FSYM : FLAG THAT INDICATES IF AIRFOIL IS SYMMETRIC 
C     FSYM > 0 MEANS AIRFOIL IS SYMMETRIC
C     FSYM = 0 MEANS CAMBERED AIRFOIL
C     NU     NO. OF DATA POINTS WHERE UPPER SURFACE
C            ORDINATES ARE SPECIFIED
C     NL     NO. OF DATA POINTS WHERE LOWER SURFACE
C            ORDINATES ARE SPECIFIED
C     X,Y    INPUT X,Y COORDINATES OF THE AIRFOIL
C
C     XN,YN  COORDINATES OF A SINGULAR POINT INSIDE THE
C            AIRFOIL NEEDED BY THE KARMAN-TREFFZ MAPPING
C
C    IMAX    TOTAL NUMBER OF I- NODES ON THE BODY-FITTED
C             GRID. MUST BE LESS THAN 81
C            I = IMAX-1 IS UPPER-SURFACE TRAILING EDGE
C            I = 2 IS LOWER-SURFACE TRAILING EDGE
C
C    JMAX    TOTAL NUMBER OF J- NODES ON THE BODY FITTTED
C            GRID MUST BE LESS THAN 31
C            NOTE; J = 2 IS AIRFOIL SURFACE
C
C    ************************************************************
      READ (5,1)
    1 FORMAT(1X)
      READ (5,1)
      READ (5,*) FSYM ,NU,NL , IMAX , JMAX
      IF(IMAX.GT.81) STOP
      IF(JMAX.GT.31) STOP
      IF(FSYM.GT.0) NL = NU
      N = NU + NL - 1
C
C     READ UPPER SURFACE ORDINATES
C
      READ (5,1)
      DO 1000 I =  NL , N
 1000 READ (5,*) X(I),Y(I)
C
      IF(FSYM.GT.0.0) THEN
C
C     SYMMETRIC AIRFOIL... FLIP SIGN TO GET
C     LOWER SURFACE COORDINATES
C
      DO 2000 I = NL , 1 , -1
      X(I) = X(N+1-I)
      Y(I) = - Y(N+1-I)
 2000 CONTINUE
      ELSE
C
C     CAMBERED AIRFOIL... READ LOWER SURFACE DATA TOO
C
      READ (5,1)
      DO 2010 I = NL , 1 , -1
 2010 READ (5,*) X(I) , Y(I)      
      ENDIF
      READ (5,1)
C
C     READ LEADING ENDGE SINGULAR POINT
C
      READ (5,*) XN,YN
      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    = CMPLX(XN,YN)
      Z2    = CMPLX(XT,YT)
      Z7    = 0.5 * (Z1+Z2)
      U     = 1.0
      V     = 0.0
      ANGL  = 8.0 * ATAN(1.0)
      POWER = 1.0 / POWER 
      DO 3000 I = 1 , N
      Z3        = CMPLX(X(I),Y(I))
      Z4        = (Z3 - Z2) / (Z3-Z1)
      Z4        = CLOG(Z4) * CMPLX(POWER,0.0)
      Z4        = CEXP(Z4)
      Z5        = CMPLX(1.0,0.0)
      Z4        = (Z5+Z4) * Z7 / (Z5-Z4)
      X4        = REAL(Z4)
      Y4        = AIMAG(Z4)
      R(I)      = SQRT(X4*X4+Y4*Y4)
      A         = X4
      B         = Y4
      ANGL      = ANGL + ATAN2((U*B-V*A),(U*A+V*B))
      U         = A
      V         = B 
      WRITE (6,*) 'MAG=',R(I),'THETA=',ANGL
      THETA(I)  = ANGL
 3000 CONTINUE
C
C     LINEARLY INTERPOLATE R VS. THETA TO GET R AT UNIFORMLY
C     SPACED THETA LOCATIONS
C
      DTHETA = 8. * ATAN(1.0) / (IMAX-3)
      DO 3001 I = 3 , IMAX - 2
      THETA1(I) = (I-2) * DTHETA
      DO 3010 L = 1 , N-1
      IF((THETA(L) - THETA1(I)) * 
     1   (THETA(L+1) - THETA1(I)).LE.0.0) THEN
      DRDT      = (R(L+1) - R(L)) / (THETA(L+1)-THETA(L))
      R1(I)     = R(L) + DRDT * (THETA1(I) - THETA(L))
      WRITE (6,*) I ,THETA1(I) , R1(I)
      GOTO 3011
      END IF
 3010 CONTINUE
      WRITE (6,*) 'INTERPOLATION FAILED....'
      STOP
 3011 CONTINUE
 3001 CONTINUE
      R1(2)      = 0.5 * (R(1) + R(N))
      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

C
C     GENERATE GRID
C
      DR        = 0.99 / (JMAX-2)
      DO 4000 J = 2 , JMAX
      ETA       = (J-2) * DR
      ETA1      = 1. / (1. - ETA)
      DO 4000 I = 1 , IMAX
      X4        = R1(I) * ETA1 * COS(THETA1(I))
      Y4        = R1(I) * ETA1 * SIN(THETA1(I))
      Z4        = CMPLX(X4,Y4)
      Z4        = (Z4-Z7) / (Z4+Z7)
      Z4        = CMPLX(1.0/POWER,0.0) * CLOG(Z4)
      Z4        = CEXP(Z4)
      Z4        = (Z2 - Z4 * Z1) / (CMPLX(1.0,0.0) - Z4)
      X4        = REAL(Z4)
      Y4        = AIMAG(Z4)
      IF(J.EQ.2)  WRITE (6,*) I , X4 , Y4
      XG(I,J)   = X4
      YG(I,J)   = Y4
 4000 CONTINUE
C
C     EXTRAPOLATE GRID TO THE INTERIOR OF THE AIRFOIL
C     TO OBTAIN COORDINATES AT J = 1
C
      DO 4001 I = 1 , IMAX
      XG(I,1)   = 2. * XG(I,2) - XG(I,3)
 4001 YG(I,1)   = 2. * YG(I,2) - YG(I,3)
      OPEN(UNIT=9,FILE='GRID.DAT',STATUS='NEW')
      WRITE (9,*) IMAX,JMAX
      WRITE (9,4002) ((I,J,XG(I,J),YG(I,J),I=1,IMAX),
     1                 J=1,JMAX)
      STOP
 4002 FORMAT(2I5,2F10.4)
      END