c
C     MAIN PROGRAM TO READ AND POST PROCESS PGM PANELv2 DATA

C     W.H. MASON, FEBRUARY 1988 (mod for PANELv2 - March 1992)

c     this version, Mar. 8, 1994

      character*1  ifile
      character*24 filenm, fileln
      character*79 filet1, filet2

      DIMENSION xc(140),yc(140),S(140),CP(140),ue(140)

      iwrit = 6
      iread = 5

      WRITE(IWRIT,900)

c     read a file containing output from PANELv2

      write (iwrit,910)
      read  (iread,*) filenm

      open(unit = 3, file = filenm, status = 'old')

      read  (3,920)     filet1
      write (iwrit,925) filet1

      read  (3,920)     filet2
      write (iwrit,935) filet2

      read(3,990)       alpha,cl,cmc4,cd
      write(iwrit,995)  alpha,cl,cmc4,cd

      read(3,930) xmc
      mc = xmc - 2.0
      write(iwrit,930)  xmc

      read  (3,920)     filet2
      write (iwrit,935) filet2


      do  10 i = 1,mc
      read(3,930)         xc(i),yc(i),cp(i),ue(i)
   10 write(iwrit,932) i, xc(i),yc(i),cp(i),ue(i)

      KWRIT = 1
C
      CALL STAGPT(xc,yc,ue,cp,s,XSP,YSP,SSP,MC,JS,JLE,KWRIT,IWRIT)
C
      ISTAGP = 1
      JMC    = MC
C
      E1     = SSP - S(JS)
      E2     = S(JS+1) - SSP
      WRITE(IWRIT,125) E1,E2
  125 FORMAT(8X,'E1     = ',F11.7,5X,'E2    = ',F11.7)
C
      IF( E1.LT.0.00001 .OR. E2.LT.0.00001 ) ISTAGP = 0
      WRITE(IWRIT,127) ISTAGP
  127 FORMAT(8X,'ISTAGP =  ',I2/)
C
      JMX    = JMC + ISTAGP
      JX1    = JS + 1

      jlower = js  + istagp
      jupper = jmc - jx1 + istagp + 1

C     WRITE OUT RESULTS TO A NEW FILE

      WRITE(IWRIT,130)
  130 FORMAT(/8X,'OUTPUT OF POSTp RESULTS')

C     SEND OUTPUT TO A FILE?

      WRITE(IWRIT,1020)
      READ(IREAD,*) IFILE
      IF(IFILE .EQ. 'N' .or. IFILE .EQ. 'n') stop

      write(iwrit,1040) 
      read(iread,*) filenm
      open(unit = 2, file = filenm, status = 'new')
      ipun  = 2

      fmach   = 0.1
      CPSTAG  = 1.0

      WRITE(IPUN,920) filet1
      WRITE(IPUN,112) alpha,fmach,cl,cmc4,cd
C
      WRITE(IWRIT,161) alpha,fmach,cl,cmc4,cd,jupper,jlower
  161 FORMAT(/8X,'ALPHA    = ',F10.5/8X,'MACH NO. = ',F10.5/
     1 8X,       'CL       = ',F10.5/8X,'CMC4     = ',F10.5/
     2 8X,       'CD       = ',F10.5/
     3 8X,   'No. of upper surface values in x/c, Cp table = ',i4/
     4 8X,   'No. of lower surface values in x/c, Cp table = ',i4//)

      STAG   = 0.0
      ww     = 0.0
      sspref = 0.0

c     output lower surface

      write(iwrit,950)
      write(ipun,955)
      WRITE(IWRIT,162)

      jout = 0

      IF(ISTAGP.EQ.1) then
                      jout = 1
                      WRITE(IPUN,240)  SSPref,CPSTAG,ww
                      WRITE(IWRIT,262) jout,XSP,YSP,SSPref,STAG,CPSTAG
                      end if

      DO 200 J = 1, JS
      jw       = js -j + 1
      jout     = jout + 1
      sarc = ssp - s(jw)
      WRITE(IPUN,240)  sarc,CP(Jw),ww
      WRITE(IWRIT,262) jout,xc(Jw),yc(Jw),sarc,ue(Jw),CP(Jw)
  200 CONTINUE

c     output upper surface

      write(iwrit,980)
      write(ipun,985)
      WRITE(IWRIT,162)

      jout = 0

      IF(ISTAGP.EQ.1) then
                      jout = 1
                      WRITE(IPUN,240)  SSPref,CPSTAG,ww
                      WRITE(IWRIT,262) Jout,XSP,YSP,SSPref,STAG,CPSTAG
                      end if

      DO 210 JJ= JX1, JMC
      jout      = jout + 1
      J         = JJ   + ISTAGP
      sarc      = s(jj) - ssp
      WRITE(IPUN,240) sarc,CP(JJ),ww
      WRITE(IWRIT,262) jout,xc(JJ),yc(JJ),sarc,ue(JJ),CP(JJ)
  210 CONTINUE
 
      stop

  112 format(6f11.5)
  162 FORMAT(/8X,'J',6X,'X/C',8X,'Y/C',8X,'S/C',7X,'U/UINF',
     1       7X,'CP')
  240 FORMAT(3F10.6)
  262 FORMAT(6X,I3,6F11.6)
  900 FORMAT(/5X,'PGM POSTP  - POST PROCESS DATA FROM PGM. PANELv2'//
     1        5X,'ECHO OF INPUT DATA:'/)
  910 format(/2x,'Enter name of file to be read:')
  920 format(a79)
  925 format(/3x,'Input data:'/1x,a79)
  930 format(4f20.7)
  932 format(5x,i3,f12.7,3f20.7)
  935 format(1x,a79)
  950 format(/5x,'lower surface'/)
  955 format(/5x,'lower surface'/
     1        5x,'s/c',7x,'Cp',7x,'dT/dy')
  980 format(/5x,'upper surface'/)
  985 format(/5x,'upper surface'/
     1        5x,'s/c',7x,'Cp',7x,'dT/dy')
  990 format(5f10.4)
  995 format(5f10.5)
 1020 format(/2x,'send output to a file? (Y/N):')
 1040 format(/2x,'enter file name:')

      END

      SUBROUTINE STAGPT(X, Y, SV,CP,S,XSP,YSP,SSP,IC,JS,JLE,KWRIT,IWRIT)

C    ROUTINE TO DEFINE LOCATION OF STAGNATION POINT
C    W.H. MASON, SEPT. 18, 1984

      DIMENSION S(1),X(1),Y(1),SV(1),CP(1)

C      DEFINE ARC LENGTH

      CALL ARCLNG(IC,X,Y,S)

C      LOCATE SIGN CHANGE in velocity

      DO 10 J=2,IC
      JS = J
      IF(SV(J-1).LT.0.0.AND.SV(J).GE.0.0) GO TO 30
   10 CONTINUE

   30 CONTINUE

      A   = (SV(JS-1)-SV(JS))/(S(JS-1)-S(JS))
      B   =  SV(JS-1)-A*S(JS-1)
      SSP = -B/A
      A   = (X(JS-1)-X(JS))/(S(JS-1)-S(JS))
      B   =  X(JS-1)-A*S(JS-1)
      XSP =  A*SSP+B
      A   = (Y(JS-1)-Y(JS))/(S(JS-1)-S(JS))
      B   =  Y(JS-1)-A*S(JS-1)
      YSP =  A*SSP+B

      JS = JS - 1
C
C      LOCATE LEADING EDGE POINT
C
      DO 50 J=2,IC
      JLE = J - 1
      IF(X(J).GE.X(J-1)) GO TO 60
   50 CONTINUE

   60 CONTINUE
C
      IF ( KWRIT .EQ. 0 ) RETURN
C
      WRITE(IWRIT,180)
      J1=JS-5
      J2=JS+5
      DO 20 J=J1,J2
   20 WRITE(IWRIT,190) J,X(J),Y(J),S(J),SV(J),CP(J)
      WRITE(IWRIT,200) XSP,YSP,SSP,JS,JLE
 
      RETURN

  190 FORMAT(6X,I4,6F11.6)
  180 FORMAT(/8X,'STAGNATION PT. SEARCH'//9X,'J',6X,'X/C',
     1 8X,'Y/C',8X,'SARC',5X,'UE/UINF',7X,'CP'/)
  200 FORMAT(//8X,'STAG PT: XSP=',F9.6,2X,'YSP=',F9.6,
     1 3X,'SSP=',F8.6,2X,'JS=',I2,2X,'JLE=',I2/)

      END

      SUBROUTINE ARCLNG(N,X,Y,SARC)
C
C----- CALCULATES ARC-LENGTH FOR TWO-DIMENSIONAL CURVE
C-----                                      REVISION DATE 4-7-78
C
      DIMENSION X(1),Y(1),SARC(1)
      DOUBLE PRECISION SSUM
      ATRI(XA,YA,XB,YB,XC,YC)=.5*((XB-XA)*(YC-YB)+(XB-XC)*(YB-YA))
      XLNGTH(XA,YA,XB,YB)=SQRT((XB-XA)*(XB-XA)+(YB-YA)*(YB-YA))
      IF(N.LE.0) GO TO 999
      SSUM=0.D0
      NM=N-1
      DO 110 I=1,N
  110 SARC(I)=0.
      BLNG=XLNGTH(X(1),Y(1),X(2),Y(2))
      IF(N.GT.2) GO TO 120
      SARC(2)=BLNG
      GO TO 999
C----- CALCULATE LEFT AND RIGHT ESTIMATES OF ARC LENGTH
  120 DO 130 J=2,NM
      I=J-1
      K=J+1
      ALNG=BLNG
      BLNG=XLNGTH(X(J),Y(J),X(K),Y(K))
      CLNG=XLNGTH(X(I),Y(I),X(K),Y(K))
      AARC=ALNG
      BARC=BLNG
      AREA=ABS(ATRI(X(I),Y(I),X(J),Y(J),X(K),Y(K)))
C-----  TEST FOR THREE POINTS CO-LINEAR
      IF(AREA.EQ.0.) GO TO 140
      RADIUS=ALNG*BLNG*CLNG/(4.*AREA)
      BARC=2.*RADIUS*ASIN(.5*BLNG/RADIUS)
      AARC=2.*RADIUS*ASIN(.5*ALNG/RADIUS)
  140 CONTINUE
      SARC(J)=SARC(J)+AARC
      SARC(K)=SARC(K)+BARC
  130 CONTINUE
      IF(N.GT.3) GO TO 150
      SARC(3)=SARC(2)+0.E0+SARC(3)
      GO TO 999
C----- AVERAGE LEFT AND RIGHT ESTIMATES
  150 DO 160 I=3,NM
  160 SARC(I)=.5*SARC(I)
C-----  SUM UP RUNNING ARC LENGTH
      SSUM=SSUM+SARC(2)
      DO 170 I=3,NM
      SSUM=SSUM+SARC(I)
  170 SARC(I)=SSUM
      SARC(N)=SSUM+SARC(N)
  999 RETURN
      end
