!     Last change:  JE    4 Jan 2001   10:12 am
      PROGRAM MAPSL
!-------------------------------------------------------------------------------
!     This program uses the Karman-Trefftz mapping function to compute the flow
!     field and pressure distribution around a foil section. The geometry of the
!     foil is determined by specifying the position of the center of the circle,
!     (XC,YC) and the trailing edge angle (TAU). The circulation is obtained by
!     satisfying the Kutta condition, but the user may override this value to
!     explore the effect on the local trailing edge flow. The output consists
!     of two TECPLOT formatted files. MAPSL.PLT is used to generate a field plot
!     showing the streamlines and pressure field. PLOTCP.PLT is used to generate
!     an x-y plot of the presure distribution on the foil surface.
!     WRITTEN BY: Justin E. Kerwin    December, 2000
!-------------------------------------------------------------------------------

!     Declare the variables-----------------------------------------------------
      IMPLICIT NONE
      CHARACTER*1 YESNO
      REAL, PARAMETER :: PI=3.1415927E00, RAD=1.7453293E-02
      REAL :: X,Y,THETA,RSQ,U,V,R,XC,YC,ALPHA_DEG,ALPHA,RC,BETA,RCSQ,G,GAMMA,   &
              LAMBDA,TAU_DEG,XI,ETA,X_CIRCLE,Y_CIRCLE,FL2,U_FOIL,V_FOIL,        &
              THETA1,THETA2,DTHETA,RMAX,DR,CP,XLE,XTE,CHORD,CL
      REAL, ALLOCATABLE, DIMENSION(:) :: XI_FOIL,CP_FOIL
      COMPLEX :: Z,ZP1,ZM1,ZETA,W_CIRCLE,DZETA_DZ,W_FOIL,ZPM1
      INTEGER :: N,M,NTHETA,NRADIAL
 !    Input the variables defining foil geometry and angle of attack -----------
      WRITE(*,'(A)') ' ENTER XC,YC   '   ! Position of center of circle---------
      READ(*,*) XC,YC                    ! Note: XC must be <= 0   -------------
      WRITE(*,'(A)') ' ENTER TRAILING EDGE ANGLE (DEG)   '
      READ(*,*) TAU_DEG
      LAMBDA=2.0-TAU_DEG/180.0     ! The exponent in the K-T mapping function --
      FL2=4.0*LAMBDA**2            ! Pre-comute a mapping function constant-----
      WRITE(*,'(A)') ' ENTER ALPHA (DEG)  ' ! Note: Enter 180 to get a circle---
      READ(*,*) ALPHA_DEG          ! User selected angle of attack.-------------
      ALPHA=ALPHA_DEG*RAD
      RCSQ=(1.0-XC)**2+YC**2       ! Radius of circle passing through (1,0) ----
      RC=SQRT(RCSQ)
      BETA=ATAN(YC/(1.0-XC))
      GAMMA=-4.0*PI*RC*SIN(BETA+ALPHA) ! Circulation required to satisfy Kutta--
      WRITE(*,'('' GAMMA FOR KUTTA CONDITION='',F12.7)') GAMMA
      WRITE(*,'(A)') ' ENTER "Y" TO ACCEPT, "N" TO OVERRIDE  '
      READ(*,'(A)') YESNO
      IF(YESNO/='Y'.AND.YESNO/='y') THEN
         WRITE(*,'(A)') ' ENTER DESIRED VALUE  '
         READ(*,*) GAMMA
      END IF
      G=GAMMA/(2.0*PI)

!----- Define grid parameters --------------------------------------------------
      WRITE(*,'(A)') ' ENTER THETA1,THETA2(DEG)  '   ! Angular extent of grid---
      READ(*,*) THETA1,THETA2
      WRITE(*,'(A)') ' ENTER DTHETA(DEG)   '         ! Angular resolution-------
      READ(*,*) DTHETA
      NTHETA=NINT((THETA2+THETA1)/DTHETA)+1        ! Number of radial grid lines
      WRITE(*,'(A)') ' ENTER MAXIMUM GRID RADIUS, RMAX   ' ! Radial extent of
      READ(*,*) RMAX                                       ! the grid-----------
      WRITE(*,'(A)') ' ENTER RDIAL SPACING, DR   '
      READ(*,*) DR                                   ! Radial resolution -------
      NRADIAL=NINT((RMAX-RC)/DR)+1
      ALLOCATE(XI_FOIL(NTHETA),CP_FOIL(NTHETA))      ! Arrays for foil Cp ------

      OPEN(1,FILE='MAPSL.PLT',STATUS='UNKNOWN',FORM='FORMATTED')  ! Output files
      WRITE(1,'(A)') ' VARIABLES = X,Y,U,V,CP '
      OPEN(2,FILE='PLOTCP.PLT',STATUS='UNKNOWN',FORM='FORMATTED')
      WRITE(2,'(A)') ' VARIABLES = X,CP '

!-----Generate velocity and presure field --------------------------------------
      WRITE(1,'('' ZONE   I='',I4,'' J='',I4)') NTHETA-1,NRADIAL
      DO M=1,NRADIAL

         R=RC+REAL(M-1)*DR   ! Grid radial coordinate in the Z (circle) plane---
         RSQ=R**2

         DO N=2,NTHETA

            THETA=RAD*(THETA1-REAL(N-1)*DTHETA)     ! Grid angular coordinate---
            X=XC+R*COS(THETA)                       ! Convert to cartesian------
            Y=YC+R*SIN(THETA)

!-----------Compute the velocity field around the circle------------------------
            U=COS(ALPHA)-RCSQ*COS(2.0*THETA-ALPHA)/RSQ-G*SIN(THETA)/R
            V=SIN(ALPHA)-RCSQ*SIN(2.0*THETA-ALPHA)/RSQ+G*COS(THETA)/R

!-----------Express the field point position and velocity in complex form-------
            Z=CMPLX(X,Y)                    ! Form complex number (X+iY)--------
            W_CIRCLE=CMPLX(U,-V)            ! Complex velocity in Z plane-------

!-----------Use the Karman-Trefftz transformation to map points to ZETA plane---
            ZP1=(Z+1.0)**LAMBDA
            ZM1=(Z-1.0)**LAMBDA
            ZPM1=ZP1-ZM1
            ZETA=LAMBDA*(ZP1+ZM1)/(ZP1-ZM1)
            XI=REAL(ZETA)                   ! Transformed X coordinate----------
            ETA=AIMAG(ZETA)                 ! Transformed Y coordinate ---------

!-----------Compute the derivative of the mapping function; Transform velocities
            DZETA_DZ=FL2*((Z-1.0)**(LAMBDA-1.0)*(Z+1.0)**(LAMBDA-1.0))/         &
                     (ZP1-ZM1)**2
            IF(ABS(DZETA_DZ)>0.0001) THEN
               W_FOIL=W_CIRCLE/DZETA_DZ  ! Complex velocity in foil plane-------
            ELSE
               W_FOIL=0.0               ! Avoids divide by zero at trailing edge
            END IF
            U_FOIL=REAL(W_FOIL)
            V_FOIL=-AIMAG(W_FOIL)

!-----------Compute the presure coefficient, CP, and output to plotting file----
            CP=ABS(W_FOIL)**2-1.0
            WRITE(1,'(5F10.5)') XI,ETA,U_FOIL,V_FOIL,CP
            IF(M==1) THEN        ! Save presure on foil surface for later use---
               XI_FOIL(N)=XI
               CP_FOIL(N)=CP
            END IF
         END DO
      END DO

!-----Scale the chordwise coordinate to (0,1) and output foil pressure dist.----
      XLE=10.0
      XTE=-10.0
      DO N=2,NTHETA
         XLE=MIN(XLE,XI_FOIL(N))
         XTE=MAX(XTE,XI_FOIL(N))
      END DO
      CHORD=XTE-XLE
      DO N=2,NTHETA
         XI_FOIL(N)=(XI_FOIL(N)-XLE)/CHORD
      END DO

      WRITE(2,'(2F10.5)') (XI_FOIL(N),CP_FOIL(N),N=2,NTHETA)
      CLOSE(2)

!-----
      CL=-2.0*GAMMA/CHORD
      WRITE(*,'('' LIFT COEFFICIENT='',F8.4)') CL

!-----Generate coordinates of circle -------------------------------------------
      WRITE(1,'(A)') ' ZONE   I= 360, J=1 '
      DO M=2,361
         THETA=-BETA+REAL(M-1)*RAD
         X_CIRCLE=XC+RC*COS(THETA)
         Y_CIRCLE=YC+RC*SIN(THETA)
         Z=CMPLX(X_CIRCLE,Y_CIRCLE)
         ZP1=(Z+1.0)**LAMBDA
         ZM1=(Z-1.0)**LAMBDA
         ZETA=LAMBDA*(ZP1+ZM1)/(ZP1-ZM1)
         XI=REAL(ZETA)
         ETA=AIMAG(ZETA)
         WRITE(1,'(5F10.5)') XI,ETA,0.0,0.0,0.0
      END DO


      CLOSE(1)
      STOP
      END PROGRAM MAPSL







