C***********************************************************************
C  POLAR PLOTTING FACILITY FOR MSES AND XFOIL
C
C    INPUT:
C     * Polar file(s) generated by MSES or XFOIL
C     * Reference data files in the format:
C
C         CD(1)  CL(1)
C         CD(2)  CL(2)
C          .      .
C          .      .
C         999.0  999.0
C         alpha(1)  CL(1)
C         alpha(2)  CL(2)
C           .       .
C           .       .
C         999.0   999.0
C         alpha(1)  Cm(1)
C         alpha(2)  Cm(2)
C           .       .
C           .       .
C         999.0   999.0
C         Xtr/c(1)  CL(1)
C         Xtr/c(2)  CL(2)
C           .        .
C           .        .
C         999.0   999.0
C         
C         The number of points in each set (CD-CL, alpha-CL, etc.) 
C         is arbitrary, and can be zero.
C
C     * pplot.def  plot parameter file (optional)
C
C***********************************************************************
C
      PROGRAM PPLOT
      INCLUDE 'PPLOT.INC'
C
      LOGICAL ERROR, LGETFN
      REAL RINP(10)
      REAL CPOLO(NAX,IPTOT,NPX), VPOLO(NAX,2,NPX)
C
      LPLOT = .FALSE.
C
      PI = 4.0*ATAN(1.0)
C
      CALL PLINITIALIZE
C
C...Get default settings
      CALL GETDEF
C
C---- Check for command line args (load file names)
      NPOL = 0
      DO II=1, NPX
        FNAME = ' '
        CALL GETARG0(II,FNAME)
        IF(FNAME.NE.' ') THEN
          NPOL = NPOL + 1
          FNPOL(NPOL) = FNAME
         ELSE
          IF(NPOL.GT.0) THEN
            IOPTS = 11
            GO TO 10
           ELSE
            GO TO 1
          ENDIF
        ENDIF
      END DO
C
    1 WRITE(*,1000)
      IF(NPOL.GT.0) WRITE(*,1010)
      WRITE(*,1020)
      WRITE(*,1050)
C
 1000 FORMAT(/'  1  Read polars          (-1 for new set)'
     &       /'  2  Read reference data  (-2 for new set)'
     &       /'  3  Plot CD(CL)'
     &       /'  4  Hardcopy current plot'
     &       /'  5  Change plot settings'
     &       /'  6  Zoom'
     &       /'  7  Unzoom'
     &       /'  8  Annotation menu'
     &       /'  9  Set CD(CL) modifiers')
 1010 FORMAT( ' 11  Re-read current polars'
     &       /' 12  Re-read current reference data')
 1020 FORMAT( ' 13  Plot Vz(V)'
     &       /' 19  Set aicraft parameters')
 1050 FORMAT(/'   Select option (0=quit): ', $)
C
      READ(*,*,ERR=1) IOPTS
      IOPT = ABS(IOPTS)
C
      GO TO (900, 10, 20, 30, 40, 50, 60, 70, 80, 90, 900,
     &            10, 20,130,900,900,900,900,900,190, 900 ), IOPT+1
      GO TO 1
C
C=============================================
C---- read polars and assign colors
 10   CONTINUE
      IF    (IOPTS.EQ.-1) THEN
C----- read new polars
       IP1 = 1
       IP2 = NPX
      ELSEIF(IOPTS.EQ. 1) THEN
C----- read additional polars
       IP1 = NPOL+1
       IP2 = NPX
      ELSE
C----- re-read old polars
       IP1 = 1
       IP2 = NPOL
      ENDIF
C
      DO 105 IP = IP1, IP2
        IF(IOPTS.EQ.1 .OR. IOPTS.EQ.-1) THEN
          CALL ASKS('Enter polar data filename or <return>^',FNPOL(IP))
        ENDIF
        IF(FNPOL(IP)(1:1) .EQ. ' ') GO TO 108
C
        LU = 9
        CALL POLREAD(LU,FNPOL(IP),ERROR,
     &    NAX,NA(IP),CPOL(1,1,IP), 
     &    REYN(IP),MACH(IP),ACRIT(IP),XTRIP(1,IP),
     &    NAME(IP),IRETYP(IP),IMATYP(IP),
     &    ISX,NBL(IP),CPOLSD(1,1,1,IP),
     &    CODE,VERSION )
        IF(ERROR) THEN
         WRITE(*,*) 'Polar file READ error'
         GO TO 108
        ENDIF
C
        WRITE(*,8000) NAME(IP)
        IF(IMATYP(IP).EQ.1) WRITE(*,8011) MACH(IP)
        IF(IMATYP(IP).EQ.2) WRITE(*,8012) MACH(IP)
        IF(IMATYP(IP).EQ.3) WRITE(*,8013) MACH(IP)
        IF(IRETYP(IP).EQ.1) WRITE(*,8021) REYN(IP)/1.0E6 
        IF(IRETYP(IP).EQ.2) WRITE(*,8022) REYN(IP)/1.0E6
        IF(IRETYP(IP).EQ.3) WRITE(*,8023) REYN(IP)/1.0E6
        WRITE(*,8030) ACRIT(IP)
C
 8000   FORMAT(1X,A)
 8011   FORMAT('             Ma =', F7.3,       $)
 8012   FORMAT('    sqrt(CL)*Ma =', F7.3,       $)
 8013   FORMAT('          CL*Ma =', F7.3,       $)
 8021   FORMAT('             Re =', F7.3,' e 6',$)
 8022   FORMAT('    sqrt(CL)*Re =', F7.3,' e 6',$)
 8023   FORMAT('          CL*Re =', F7.3,' e 6',$)
 8030   FORMAT('          Ncrit =', F6.2         )
C
 105  CONTINUE
      IP = IP2+1
C
 108  CONTINUE
      NPOL = IP-1
      IP2 = MIN(IP2,NPOL)
C
      DO IP = IP1, IP2
ccc        CALL GETTYP(NAX,NA(IP),CPOL(1,1,IP),IMATYP(IP),IRETYP(IP))
        CALL STRIP(NAME(IP),NNAME)
        ICOL(IP) = 2 + IP
        ILIN(IP) = IP
      ENDDO
CCC   CALL MINMAX(NAX,NPOL,NA,CPOL,CPOLPLF)
C
C---- are these dimensional polars?
      DO IP = IP1, IP2
        CALL GETCLEN(NAME(IP),CCLEN,NCLEN)
        IF(NCLEN.GT.0) THEN
         LCLEN = .TRUE.
         GO TO 1
        ENDIF
      ENDDO
      IF(.NOT.LPLOT) GO TO 30
      GO TO 1
C
C=============================================
C---- read reference data
 20   CONTINUE
      IF(IOPTS.EQ.12 ) THEN
C------ re-read old data sets
        ID1 = 1
        ID2 = NDAT
        LGETFN = .FALSE.
      ELSEIF(IOPTS.GT.0) THEN
C------ read additional data sets
        ID1 = NDAT+1
        ID2 = NDX
        LGETFN = .TRUE.
      ELSE
C------ read new data sets
        ID1 = 1
        ID2 = NDX
        LGETFN = .TRUE.
      ENDIF
C
      DO 25 ID = ID1, ID2
        IF(LGETFN) THEN
         CALL ASKS('Enter reference data filename or <return>^',
     &              FNREF(ID))
         IF(FNREF(ID)(1:1) .EQ. ' ') GO TO 27
        ENDIF
C
        LU = 9
        OPEN(LU,FILE=FNREF(ID),STATUS='OLD',ERR=27)
        CALL POLREF(LU, FNREF(ID), ERROR,
     &              NFX, NF(1,ID), XYREF(1,1,1,ID), LABREF(ID) )
        CLOSE(LU)
        IF(ERROR) GO TO 27
C
        NDAT = ID
C
        CALL STRIP(LABREF(ID),NLAB)
        IF(NLAB.EQ.0) THEN
          CALL ASKS('Enter label for reference data^',LABREF(ID))
          CALL STRIP(LABREF(ID),NLAB)
        ENDIF
C
ccc     IFCOL(ID) = NCOLOR - ID + 1
        IFCOL(ID) = 2 + ID
        IFSYM(ID) = MOD(ID,10)
 25   CONTINUE
 27   CONTINUE
      GO TO 1
C
C=============================================
C---- Make the CD(CL) Plot
 30   IF (NPOL.EQ.0 .AND. NDAT.EQ.0) GO TO 1
C
C---- sort each polar by increasing alpha
      DO IP=1, NPOL
        CALL PLRSRT(IP,IAL)
      ENDDO
C
C---- set modified polars
      DO IP = 1, NPOL
        DO IA = 1, NA(IP)
          DO I = 1, IPTOT
            CPOLO(IA,I,IP) = CPOL(IA,I,IP)
          ENDDO
          CPOLO(IA,ICD,IP) = CPOL(IA,ICD,IP)
     &                     + CDLMOD(1,IP)
     &                     + CDLMOD(2,IP)*CPOL(IA,ICL,IP)
     &                     + CDLMOD(3,IP)*CPOL(IA,ICL,IP)**2
          IF(CDLMOD(4,IP) .NE. 1.0) THEN
           CPOLO(IA,ICL,IP) = ABS(CPOL(IA,ICL,IP))**CDLMOD(4,IP)
          ENDIF
          IF(CDLMOD(5,IP) .NE. 0.0) THEN
           CPOLO(IA,ICD,IP) = CPOLO(IA,ICD,IP)
     &                      * ABS(CPOL(IA,ICL,IP))**CDLMOD(5,IP)
          ENDIF
        ENDDO
      ENDDO

      IF (AUTO) THEN 
        CALL MINMAX(NAX,NPOL,NA,CPOLO,CPOLPLF)
        CALL SETINC
      ENDIF
C
      IF (LPLOT) CALL PLEND
      CALL PLOPEN(SCRNFR,IPSLU,IDEV)
      LPLOT = .TRUE.
C
C---- set 0.3" left,bottom margins
      CALL PLOTABS(0.3,0.3,-3)
      CALL NEWFACTOR(SIZE)
      CALL PLOT(6.0*CH,6.0*CH,-3)
C

c      WRITE(*,*) CPOLPLF(1,ICL),CPOLPLF(2,ICL),CPOLPLF(3,ICL)
c      write(*,*)

      CALL POLPLT(NAX,NPOL,NA,CPOLO,
     &            REYN,MACH,ACRIT, NAME ,ICOL,ILIN,
     &            NFX,NDAT,NF,XYREF,LABREF,IFCOL,IFSYM,
     &            ISX,NBL,CPOLSD, IMATYP,IRETYP,
     &            TITLE,CODE,VERSION,
     &            PLOTAR, XCD,XAL,XOC, CH,CH2, CDLMOD(4,1),
     &            LGRID,LCDW,LLIST,LEGND,LAECEN,
     &            CPOLPLF, CCLEN,NCLEN )
C
c      CALL POLFIT(NAX,NPOL,NA,CPOLO,
c     &            REYN,MACH,ACRIT, NAME ,ICOL,ILIN,
c     &            IMATYP,IRETYP,
c     &            PLOTAR, XCD,XAL,XOC, CH,CH2, CDLMOD(4,1),
c     &            CPOLPLF, CCLEN,NCLEN )
C
      GO TO 1
C
C=============================================
C---- hardcopy output
 40   IF(LPLOT) CALL PLEND
      LPLOT = .FALSE.
      CALL REPLOT(IDEVRP)
      GO TO 1
C
C=============================================
C---- change settings
 50   CALL GETSET
      GO TO 1
C
C=============================================
C---- zoom
 60   CALL USETZOOM(.FALSE.,.TRUE.)
      CALL REPLOT(IDEV)
      GO TO 1
C
C=============================================
C---- unzoom
 70   CALL CLRZOOM
      CALL REPLOT(IDEV)
      GO TO 1
C
C=============================================
C---- annotate plot
 80   IF(.NOT.LPLOT) THEN
       WRITE(*,*) 'No active plot to annotate'
       GO TO 1
      ENDIF
      CALL ANNOT(CH)
      GO TO 1
C=============================================
C---- get modifiers
 90   CONTINUE
      WRITE(*,4900)
 4900 FORMAT(/' CD_plotted  = (CD  +  CD0 + CD1*CL + CD2*CL^2)*CL^expD'
     &       /' CL_plotted  =  CL^exp')
      DO IP = 1, NPOL
 91     WRITE(*,4910) IP, (CDLMOD(K,IP), K=1, 5)
 4910   FORMAT(/' Polar', I3,'...'
     &         /' Currently CD0,CD1,CD2,exp,expD = ', 3F10.6, 2F10.4,
     &         /' Input new CD0,CD1,CD2,exp,expD:    ', $)
        CALL READR(5,CDLMOD(1,IP),ERROR)
        IF(ERROR) GO TO 91
      ENDDO
      GO TO 1
C
C=============================================
C---- Make the Vz(V) Plot
 130  IF (NPOL.EQ.0 .AND. NDAT.EQ.0) GO TO 1
C
C---- sort each polar by increasing alpha
      DO IP=1, NPOL
        CALL PLRSRT(IP,IAL)
      ENDDO
C
C---- set V and Vz for plotting
      DO IP = 1, NPOL
        WOS = VPPARS(1,IP)
        RHO = VPPARS(2,IP)
        AR  = VPPARS(3,IP)
        CD0 = VPPARS(4,IP)
        REF = VPPARS(5,IP)
        REX = VPPARS(6,IP)
C
        IF(WOS .EQ. 0.0) THEN
         WRITE(*,*) 'Wing loading W/S not defined.  Using 1.0'
         WOS = 1.0
        ENDIF
        IF(RHO .EQ. 0.0) THEN
         WRITE(*,*) 'Air density RHO not defined.  Using 1.0'
         RHO = 1.0
        ENDIF
        IF(AR .EQ. 0.0) THEN
         WRITE(*,*) 'Aspect ratio AR not defined.  Using 1.0'
         AR = 1.0
        ENDIF
        IF(REF .EQ. 0.0) THEN
         WRITE(*,*) 'Reference REref not defined.  Using 10^6'
         REF = 1.0E6
        ENDIF
C
        DO IA = 1, NA(IP)
          CDP = CPOL(IA,ICD,IP)
          CL  = CPOL(IA,ICL,IP)
          RE  = CPOL(IA,IRE,IP)
C
          CLM = MAX( CL , 0.001 )
          VEL = SQRT( 2.0*WOS/(RHO*CLM) )
C
          CD = CDP
     &       + CL*CL/(PI*AR)
     &       + CD0*(RE/REF)**REX
C
          VZ = -VEL * CD/CL
C
          VPOLO(IA,1,IP) = VEL
          VPOLO(IA,2,IP) = VZ
        ENDDO
      ENDDO

      IF (AUTO) THEN 
        CALL MINMAX(NAX,NPOL,NA,VPOLO,VPOLPLF)
        CALL SETINCV
      ENDIF
C
      IF (LPLOT) CALL PLEND
      CALL PLOPEN(SCRNFR,IPSLU,IDEV)
      LPLOT = .TRUE.
C
C---- set 0.3" left,bottom margins
      CALL PLOTABS(0.3,0.3,-3)
      CALL NEWFACTOR(SIZE)
      CALL PLOT(6.0*CH,6.0*CH,-3)
C
      CALL VEPPLT(NAX,NPOL,NA,VPOLO,
     &            REYN,MACH,ACRIT, NAME ,ICOL,ILIN,
     &            IMATYP,IRETYP,
     &            TITLE,CODE,VERSION,
     &            PLOTAR, CH,CH2, 
     &            LGRID,LLIST,LEGND,
     &            VPOLPLF)
      GO TO 1
C
C=============================================
C---- get velocity-polar parameters
 190  CONTINUE
      DO IP = 1, NPOL
 191    WRITE(*,5910) IP, (VPPARS(K,IP), K=1, 6)
 5910   FORMAT(
     &  /' Polar', I3,'...'
     &  /' Currently W/S,rho,AR,CDo,REref,REexp = ', 
     &              G12.4,G12.4,F7.2,F10.6,G12.4,F6.2
     &  /' Input new W/S,rho,AR,CDo,REref,REexp:    ', $)
        CALL READR(6,VPPARS(1,IP),ERROR)
        IF(ERROR) GO TO 191
      ENDDO
      GO TO 1
C
C=============================================
  900 CALL PLCLOSE
      STOP
      END ! PPLOT


      SUBROUTINE GETCLEN(NAME,CLEN,NCLEN)
      CHARACTER*(*) NAME, CLEN
C--------------------------------------------------
C     Looks for substring  "(c=01234***)"
C     in the NAME string.  If found, then
C     the "***" string is returned in CLEN.
C     If not found, then CLEN is returned blank.
C--------------------------------------------------
C
      CLEN = ' '
C
      K1 = INDEX( NAME , '(c=' )
      IF(K1.EQ.0) RETURN
C
      NNAME = LEN(NAME)
      K2 = INDEX( NAME(K1:NNAME) , ')' ) + K1 - 2
      IF(K2-K1.LT.3) RETURN
C
      DO K = K1+3, K2
        IF(INDEX( '0123456789.,)' , NAME(K:K) ) .EQ. 0) THEN
         CLEN = NAME(K:K2)
         NCLEN = K2-K+1
         RETURN
        ENDIF
      ENDDO
C
      RETURN
      END

 

      SUBROUTINE MINMAX(NAX,NPOL,NA,CPOL,CPOLPLF)
      INCLUDE 'PINDEX.INC'
      DIMENSION NA(NPOL)
      DIMENSION CPOL(NAX,IPTOT,NPOL), CPOLPLF(3,*)
C--------------------------------------------
C     Determines max and min limits of polar
C     quantities among all polars passed in.
C--------------------------------------------
C
      IF(NPOL.LT.1) RETURN
C
      DO K = 1, 4
        CPOLPLF(1,K) = CPOL(1,K,1)
        CPOLPLF(2,K) = CPOL(1,K,1)
      END DO
C
      DO IP=1, NPOL
        DO K=1, 4
          DO I=1, NA(IP)
            CPOLPLF(1,K) = MIN( CPOL(I,K,IP) , CPOLPLF(1,K) )
            CPOLPLF(2,K) = MAX( CPOL(I,K,IP) , CPOLPLF(2,K) )
          END DO
        END DO
      END DO
C
      RETURN
      END ! MINMAX


 
      SUBROUTINE GETDEF
      INCLUDE 'PPLOT.INC'
C
C---- Plotting flag
      IDEV = 1   ! X11 window only
c     IDEV = 2   ! B&W PostScript output file only (no color)
c     IDEV = 3   ! both X11 and B&W PostScript file
c     IDEV = 4   ! Color PostScript output file only 
c     IDEV = 5   ! both X11 and Color PostScript file 
C
C---- Re-plotting flag (for hardcopy)
c      IDEVRP = 2   ! B&W PostScript
      IDEVRP = 4   ! Color PostScript
C
C---- PostScript output logical unit and file specification
      IPSLU = 0  ! output to file  plot.ps   on LU 4    (default case)
c     IPSLU = ?  ! output to file  plot?.ps  on LU 80+?
C
C---- screen fraction taken up by plot window upon opening
      SCRNFR = 0.70
C
C---- Default plot size in inches
C-    (Default plot window is 11.0 x 8.5)
      SIZE = 10.0
C
C---- plot aspect ratio V/H
      PLOTAR = 0.60
C
C---- character height
      CH  = 0.014
      CH2 = 0.012
C
C---- set default color table and get number of colors
      CALL COLORMAPDEFAULT
      CALL GETNUMCOLOR(NCOLOR)
C
      LGRID = .TRUE.
      LCDW  = .FALSE.
      LLIST = .TRUE.
      LEGND = .TRUE.
      LCLEN = .FALSE.
      LAECEN = .FALSE.
C
C---- automatic scaling for axes
      AUTO  = .TRUE.
C
      CPOLPLF(1,ICL) = 0.0   ! CLmax
      CPOLPLF(2,ICL) = 1.5   ! CLmin
      CPOLPLF(3,ICL) = 0.5   ! Axis CL increment
C         
      CPOLPLF(1,ICD) = 0.0   ! CDmax
      CPOLPLF(2,ICD) = 0.02  ! CDmin
      CPOLPLF(3,ICD) = 0.01  ! Axis CD increment
C         
      CPOLPLF(1,ICM) =  0.0  ! CMmax
      CPOLPLF(2,ICM) = -0.25 ! CMmin
      CPOLPLF(3,ICM) =  0.05 ! Axis CM increment
C         
      CPOLPLF(1,IAL) = -4.0  ! ALmax
      CPOLPLF(2,IAL) = 10.0  ! ALmin
      CPOLPLF(3,IAL) =  2.0  ! Axis AL increment
C
C---- Plot layout (relative X size to CL-CD, CL-alfa, transition plots)
      XCD = 0.45
      XAL = 0.25
      XOC = 0.20
C
C---- Set CL,CD modifiers
      DO IP = 1, NPX
        CDLMOD(1,IP) = 0.
        CDLMOD(2,IP) = 0.
        CDLMOD(3,IP) = 0.
        CDLMOD(4,IP) = 1.0
        CDLMOD(5,IP) = 0.
      ENDDO
C
cC---- Set CL,CD quadratic-fit polar parameters
c      DO IP = 1, NPX
c        CDLFIT(1,IP) = 0.
c        CDLFIT(2,IP) = 0.
c        CDLFIT(3,IP) = 0.
c        CDLFIT(4,IP) = 1.0
c        CDLFIT(5,IP) = 0.
c      ENDDO
C
C---- velocity polar plot axis parameters
      VPOLPLF(1,1) =  0.0  ! Vmin
      VPOLPLF(2,1) = 20.0  ! Vmax
      VPOLPLF(3,1) =  2.0  ! Vdel
C
      VPOLPLF(1,2) = -5.0  ! Vzmin
      VPOLPLF(2,2) =  1.0  ! Vzmax
      VPOLPLF(3,2) =  0.5  ! Vzdel
C
C---- Set Vz(V) parameters
      DO IP = 1, NPX
        VPPARS(1,IP) = 0.
        VPPARS(2,IP) = 0.
        VPPARS(3,IP) = 0.
        VPPARS(4,IP) = 0.
        VPPARS(5,IP) = 0.
        VPPARS(6,IP) = 0.
      ENDDO
C
      TITLE = '                                '
CCC            12345678901234567890123456789012
C
C...Try to read default file "pplot.def" for stored plot setup
      OPEN(UNIT=10,FILE='pplot.def',STATUS='OLD',ERR=900)
      CALL RDDEF(10)
      RETURN
C
  900 WRITE(*,*)
      WRITE(*,*) 'No  pplot.def  file found'
      WRITE(*,*) 'Hard-wired defaults used'
      WRITE(*,*)
      RETURN
      END ! GETDEF


 
 
      SUBROUTINE RDDEF(LU)
C--- Read PPLOT plot parameters from save file
      INCLUDE 'PPLOT.INC'
C
      READ(LU,*,ERR=90) CPOLPLF(1,ICL), CPOLPLF(2,ICL), CPOLPLF(3,ICL)
      READ(LU,*,ERR=90) CPOLPLF(1,ICD), CPOLPLF(2,ICD), CPOLPLF(3,ICD)
      READ(LU,*,ERR=90) CPOLPLF(1,ICM), CPOLPLF(2,ICM), CPOLPLF(3,ICM)
      READ(LU,*,ERR=90) CPOLPLF(1,IAL), CPOLPLF(2,IAL), CPOLPLF(3,IAL)
      READ(LU,*,ERR=90) XCD, XAL, XOC
      READ(LU,*,ERR=90) SIZ, PLOTAR, LGRID
      READ(LU,*,ERR=90) CH, CH2
      IF(SIZ.GT.0.0) SIZE = SIZ
      AUTO = .FALSE.
      RETURN
C
 90   WRITE(*,*)
      WRITE(*,*) '*** Error reading PPLOT parameter file'
      WRITE(*,*)
      RETURN
      END

 
      SUBROUTINE WRTDEF(LU)
C--- Write PPLOT plot parameters to save file
      INCLUDE 'PPLOT.INC'
      CHARACTER*1 CGRID
C
      CGRID = 'F'
      IF(LGRID) CGRID = 'T'
C
      WRITE(LU,1010) CPOLPLF(1,ICL), CPOLPLF(2,ICL), CPOLPLF(3,ICL)
      WRITE(LU,1020) CPOLPLF(1,ICD), CPOLPLF(2,ICD), CPOLPLF(3,ICD)
      WRITE(LU,1030) CPOLPLF(1,ICM), CPOLPLF(2,ICM), CPOLPLF(3,ICM)
      WRITE(LU,1040) CPOLPLF(1,IAL), CPOLPLF(2,IAL), CPOLPLF(3,IAL)
      WRITE(LU,1050) XCD, XAL, XOC
      WRITE(LU,1070) SIZE, PLOTAR, CGRID
      WRITE(LU,1080) CH, CH2
      RETURN
C
C...............................................
 1010 FORMAT(1X,3F9.4          ,' | CLmin   CLmax    dCL')
 1020 FORMAT(1X,3F9.4          ,' | CDmin   CDmax    dCD')
 1030 FORMAT(1X,3F9.4          ,' | CMmin   CMmax    dCM')
 1040 FORMAT(1X,3F9.4          ,' | ALmin   ALmax    dAL')
 1050 FORMAT(1X,3F9.4          ,' | CL-CD   CL-alpha  CL-Xtr  (widths)')
 1070 FORMAT(1X,2F9.4,1X,A1,7X ,' | width height/width  grid_plot_flag')
 1080 FORMAT(1X,F9.4,F9.4  ,9X ,' | char.height1  char.height2')
      END ! WRTDEF
 


      SUBROUTINE GETSET
      INCLUDE 'PPLOT.INC'
      LOGICAL OK, ERROR
      CHARACTER OPTION*2, LINE*80
      REAL RINP(10)
C
C---- Change plotting parameters
C
    1 WRITE(*,1000)
 1000 FORMAT(/ '  1   Change CL scaling'
     &       / '  2   Change CD scaling'
     &       / '  3   Change CM scaling'
     &       / '  4   Change ALPHA scaling'
     &       / '  5   Plot Layout'
     &       / '  6   Plot Title'
     &       / '  7   Plot Size'
     &      // '  8   Toggle autoscaling'
     &       / '  9   Toggle Pressure-CD plot'
     &       / ' 10   Toggle airfoil list'
     &       / ' 11   Toggle plot CL-CD legend box'
     &       / ' 12   Toggle aerodynamic center plotting'
     &       / ' 13   Toggle color Hardcopy'
     &      // ' 14   Read  settings from defaults file'
     &       / ' 15   Write settings to   defaults file'
     &       / ' 16   Rescale forces by chord factor'
     &       / ' 17   Change reference-length unit'
     &       / ' 18   Polar colors'
     &       / ' 19   Polar line styles'
     &      // ' 21   Change V  scaling'
     &       / ' 22   Change Vz scaling'
     &      // '    Select option:  ',$)
C
      READ(*,1005) OPTION
 1005 FORMAT(A)
C
      IF(OPTION .EQ. ' ' .OR. OPTION.EQ.'0') THEN
C
        RETURN
C
      ELSE IF(OPTION.EQ.'1') THEN
C--- Get CL min,max,delta
        WRITE(*,1100) (CPOLPLF(K,ICL), K=1, 3)
 20     READ(*,1005)  LINE
        NINP = 3
        CALL GETFLT(LINE,CPOLPLF(1,ICL),NINP,ERROR)
        IF(ERROR) GO TO 20
        IF(NINP.EQ.0) GO TO 1
        AUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'2') THEN
C--- Get CD min,max,delta
        WRITE(*,1200) (CPOLPLF(K,ICD), K=1, 3)
 30     READ(*,1005)  LINE
        NINP = 3
        CALL GETFLT(LINE,CPOLPLF(1,ICD),NINP,ERROR)
        IF(ERROR) GO TO 30
        IF(NINP.EQ.0) GO TO 1
        AUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'3') THEN
C--- Get CM min,max,delta
        WRITE(*,1300) (CPOLPLF(K,ICM), K=1, 3)
 40     READ(*,1005)  LINE
        NINP = 3
        CALL GETFLT(LINE,CPOLPLF(1,ICM),NINP,ERROR)
        IF(ERROR) GO TO 40
        IF(NINP.EQ.0) GO TO 1
        AUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'4') THEN
C--- Get ALFA min,max,delta
        WRITE(*,1400) (CPOLPLF(K,IAL), K=1, 3)
 50     READ(*,1005)  LINE
        NINP = 3
        CALL GETFLT(LINE,CPOLPLF(1,IAL),NINP,ERROR)
        IF(ERROR) GO TO 50
        IF(NINP.EQ.0) GO TO 1
        AUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'5') THEN
C--- Get Layout offsets for CL-CD,CL-alfa,transition plot sections
 80     WRITE(*,1700) XCD,XAL,XOC
        READ(*,1005)  LINE
        RINP(1) = XCD
        RINP(2) = XAL
        RINP(3) = XOC
        NINP = 3
        CALL GETFLT(LINE,RINP,NINP,ERROR)
        IF(ERROR) GO TO 80
        IF(NINP.EQ.0) GO TO 1
        XCD = RINP(1)
        XAL = RINP(2)
        XOC = RINP(3)
C
      ELSE IF(OPTION.EQ.'6') THEN
C--- Get plot title
        TITLE = ' '
        CALL ASKS('Enter plot title (80 chars)^',TITLE)
        CALL STRIP(TITLE,NTITLE)
C
      ELSE IF(OPTION.EQ.'7') THEN
C--- Get plot size
 60     WRITE(*,1500) SIZE
        READ(*,1005)  LINE
        IF(LINE.EQ.' ') GO TO 1
        READ(LINE,*,ERR=60) SIZE
C
      ELSE IF(OPTION.EQ.'8') THEN
C--- Toggle autoscaling
        AUTO = .NOT. AUTO
        IF(AUTO) THEN 
          WRITE(*,*) 'Axes auto-scaling enabled'
         ELSE
          WRITE(*,*) 'Axes auto-scaling disabled'
        ENDIF
C
      ELSE IF(OPTION.EQ.'9') THEN
C--- CDp plot toggle
        LCDW = .NOT. LCDW
        IF(     LCDW) WRITE(*,*) 'CDp will be plotted'
        IF(.NOT.LCDW) WRITE(*,*) 'CDp will not be plotted'
C
      ELSE IF(OPTION.EQ.'10') THEN
C--- List of polars plotted toggle
        LLIST = .NOT. LLIST
        IF(     LLIST) WRITE(*,*) 'List of polars will be plotted'
        IF(.NOT.LLIST) WRITE(*,*) 'List of polars will not be plotted'
C
      ELSE IF(OPTION.EQ.'11') THEN
C--- Legend plotted toggle
        LEGND = .NOT. LEGND
        IF(     LEGND) WRITE(*,*) 'Legend will be plotted'
        IF(.NOT.LEGND) WRITE(*,*) 'Legend will not be plotted'
C
      ELSE IF(OPTION.EQ.'12') THEN
C--- Legend plotted toggle
        LAECEN = .NOT. LAECEN
        IF(     LAECEN) WRITE(*,*) 'Aero.center will be plotted'
        IF(.NOT.LAECEN) WRITE(*,*) 'Aero.center will not be plotted'
C
      ELSE IF(OPTION.EQ.'13') THEN
C--- Color hardcopy toggle
        IF(IDEVRP.EQ.2) THEN
          IDEVRP = 4
          WRITE(*,*) 'Switching to color hardcopy'
         ELSE
          IDEVRP = 2
          WRITE(*,*) 'Switching to B&W hardcopy'
        ENDIF
C
      ELSE IF(OPTION.EQ.'14') THEN
C--- Read defaults from pplot.def file
        LINE = 'Enter settings filename  [pplot.def] ^'
        CALL ASKS(LINE,FNAME)
        IF(FNAME.EQ.' ') FNAME = 'pplot.def'
        OPEN(10,FILE=FNAME,STATUS='OLD',ERR=703)
        CALL RDDEF(10)
        CLOSE(10)
        GO TO 1
 703    WRITE(*,*) 
        WRITE(*,*) 'Open error on pplot defaults file'
        GO TO 1
C
      ELSE IF(OPTION.EQ.'15') THEN
C--- Save defaults to parameter file
        LINE = 'Enter settings filename  [pplot.def] ^'
        CALL ASKS(LINE,FNAME)
        IF(FNAME.EQ.' ') FNAME = 'pplot.def'
        OPEN(10,FILE=FNAME,STATUS='OLD',ERR=803)
        CALL ASKL('File exists.  Overwrite ?^',OK)
        IF(OK) GO TO 806
        WRITE(*,*)
        WRITE(*,*) 'No action taken'
        CLOSE(10)
        GO TO 1
 803    OPEN(10,FILE=FNAME,STATUS='UNKNOWN')
 806    CALL WRTDEF(10)
        WRITE(*,*)
        WRITE(*,*) 'PPLOT plot settings written to file'
        CLOSE(10)
C
      ELSE IF(OPTION.EQ.'16') THEN
C--- rescale forces and moments
      WRITE(*,1900)
      SCAL = 1.0
      READ(*,1005) LINE
      IF(LINE.EQ.' ') GO TO 1
      READ(LINE,*,ERR=1,END=1) CSCAL
      IF(SCAL.NE.0.0) CALL RESCAL(1.0/SCAL)
C
      ELSE IF(OPTION.EQ.'17') THEN
C--- change reference length unit
        WRITE(*,2000)
        CALL ASKS(
     &   'Enter new reference length unit (<return> if none)^',CCLEN)
        CALL STRIP(CCLEN,NCLEN)
C
      ELSE IF(OPTION.EQ.'18') THEN
C------ change polar colors
        IF(NPOL.EQ.0) THEN
         WRITE(*,*) 'No current polars to change'
         GO TO 1
        ELSE
 820     WRITE(LINE,3100) 'polar colors', (ICOL(IP), IP=1, NPOL)
         WRITE(*,1005) LINE
         WRITE(*,3105)    'polar colors'
         READ(*,1005) LINE
         NINP = NPOL
         CALL GETINT(LINE,ICOL,NINP,ERROR)
         IF(ERROR) GO TO 820
        ENDIF
C
      ELSE IF(OPTION.EQ.'19') THEN
C------ change polar line styles
        IF(NPOL.EQ.0) THEN
         WRITE(*,*) 'No current polars to change'
         GO TO 1
        ELSE
 830     WRITE(LINE,3100) 'polar line styles', (ILIN(IP), IP=1, NPOL)
         WRITE(*,1005) LINE
         WRITE(*,3105)    'polar line styles'
         READ(*,1005) LINE
         NINP = NPOL
         CALL GETINT(LINE,ILIN,NINP,ERROR)
         IF(ERROR) GO TO 830
        ENDIF
C
      ELSE IF(OPTION.EQ.'21') THEN
C--- Get V min,max,delta
        WRITE(*,2100) (VPOLPLF(K,1), K=1, 3)
 210    READ(*,1005)  LINE
        NINP = 3
        CALL GETFLT(LINE,VPOLPLF(1,1),NINP,ERROR)
        IF(ERROR) GO TO 210
        IF(NINP.EQ.0) GO TO 1
        AUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'22') THEN
C--- Get Vz min,max,delta
        WRITE(*,2200) (VPOLPLF(K,2), K=1, 3)
 220    READ(*,1005)  LINE
        NINP = 3
        CALL GETFLT(LINE,VPOLPLF(1,2),NINP,ERROR)
        IF(ERROR) GO TO 220
        IF(NINP.EQ.0) GO TO 1
        AUTO = .FALSE.
C
      ENDIF
      GO TO 1
C
 1100 FORMAT(/' Current   CLmin, CLmax, dCL = ',3F10.5
     &       /' Enter new CLmin, CLmax, dCL:  ',$)
 1200 FORMAT(/' Current   CDmin, CDmax, dCD = ',3F10.5
     &       /' Enter new CDmin, CDmax, dCD:  ',$)
 1300 FORMAT(/' Current   CMmin, CMmax, dCM = ',3F10.5
     &       /' Enter new CMmin, CMmax, dCM:  ',$)
 1400 FORMAT(/' Current   ALmin, ALmax, dAL = ',3F10.5
     &       /' Enter new ALmin, ALmax, dAL:  ',$)
 1500 FORMAT(/' Current   plot size = ', F10.5
     &       /' Enter new plot size:  ',$)
 1700 FORMAT(/'  Current layout offsets  xCD =',F8.4,
     &        '  xALPHA = ',F8.4,'  xTR = ',F8.4/
     &        ' Enter new xCD, xALPHA, xTR:  ',$)
 1800 FORMAT(/' Default settings file: ',A)
 1900 FORMAT(/'Enter chord scale factor for forces: ',$)
 2000 FORMAT(/'Current reference length unit: ', A)

 2100 FORMAT(/' Current   Vmin, Vmax, dV = ',3F10.5
     &       /' Enter new Vmin, Vmax, dV:  ',$)
 2200 FORMAT(/' Current   Vzmin, Vzmax, dVz = ',3F10.5
     &       /' Enter new Vzmin, Vzmax, dVz:  ',$)
C
 3100 FORMAT(' Currently ',A,' =', 20I3)
 3105 FORMAT(' Enter new ',A,':  ',$)
C
      END ! GETSET



      SUBROUTINE PLRSRT(IP,IDSORT)
      INCLUDE 'PPLOT.INC'
      DIMENSION INDX(NAX), ATMP(NAX)
C
C---- sort polar in increasing variable IDSORT
      CALL HSORT(NA(IP),CPOL(1,IDSORT,IP),INDX)
C
C---- do the actual reordering
      DO ID = 1, IPTOT
        CALL ASORT(NA(IP),CPOL(1,ID,IP),INDX,ATMP)
      ENDDO
      DO ID = 1, JPTOT
        DO IS = 1, 2
          CALL ASORT(NA(IP),CPOLSD(1,IS,ID,IP),INDX,ATMP)
        ENDDO
      ENDDO
C
      RETURN
      END ! PLRSRT
 


      SUBROUTINE GETTYP(NAX,NA,CPOL, IMATYP,IRETYP )
C
C---- Determines type of Ma(CL) and Re(CL) dependence
C
      INCLUDE 'PINDEX.INC'
C
      DIMENSION CPOL(NAX,IPTOT)
C
      IF(CPOL(NA,ICL)*CPOL(1,ICL) .LE. 0.0) THEN
        IMATYP = 1
        IRETYP = 1
        RETURN
      ENDIF
C
      IF(CPOL(NA,IMA)*CPOL(1,IMA) .LE. 0.0) THEN
        IMATYP = 1
      ELSE
        EX = LOG( CPOL(NA,IMA)/CPOL(1,IMA) )
     &     / LOG( CPOL(NA,ICL)/CPOL(1,ICL) )
        IF     (ABS(EX) .LT. 0.25) THEN
          IMATYP = 1
        ELSEIF (ABS(EX) .LT. 0.75) THEN
          IMATYP = 2
        ELSE
          IMATYP = 3
        ENDIF
      ENDIF
C
      IF(CPOL(NA,IRE)*CPOL(1,IRE) .LE. 0.0) THEN
        IRETYP = 1
      ELSE
        EX = LOG( CPOL(NA,IRE)/CPOL(1,IRE) )
     &     / LOG( CPOL(NA,ICL)/CPOL(1,ICL) )
        IF     (ABS(EX) .LT. 0.25) THEN
          IRETYP = 1
        ELSEIF (ABS(EX) .LT. 0.75) THEN
          IRETYP = 2
        ELSE
          IRETYP = 3
        ENDIF
      ENDIF
C
      RETURN
      END ! GETTYP


      SUBROUTINE RESCAL(SCAL)
      INCLUDE 'PPLOT.INC'
C--------------------------------------------
C     Rescales forces and moments
C--------------------------------------------
C---- rescale polar forces by SCAL, moments by SCAL**2
      DO IP=1, NPOL
        DO I=1, NA(IP)
          CPOL(I,ICL,IP) = CPOL(I,ICL,IP)*SCAL
          CPOL(I,ICD,IP) = CPOL(I,ICD,IP)*SCAL
          CPOL(I,ICW,IP) = CPOL(I,ICW,IP)*SCAL
          CPOL(I,ICM,IP) = CPOL(I,ICM,IP)*SCAL*SCAL
        END DO
      END DO
      RETURN
      END

 
      SUBROUTINE SETINC
      INCLUDE 'PPLOT.INC'
C--------------------------------------------
C     Determines axes increments for polars
C     from quantities for all polars read in.
C--------------------------------------------
C
      CLMAX = CPOLPLF(2,ICL)
      CLMIN = CPOLPLF(1,ICL)
      CDMAX = CPOLPLF(2,ICD)
      CDMIN = CPOLPLF(1,ICD)
      CMMAX = CPOLPLF(2,ICM)
      CMMIN = CPOLPLF(1,ICM)
      ALMAX = CPOLPLF(2,IAL)
      ALMIN = CPOLPLF(1,IAL)
C
C--- CL axes
      CALL AXISADJ2(CLMIN,CLMAX,CLSPAN,DCL,NCLTICS)
C--- CD axes
      CDMIN = 0.0
      CALL AXISADJ2(CDMIN,CDMAX,CDSPAN,DCD,NCDTICS)
C--- CM axes
      IF(ABS(CMMAX).GT.ABS(CMMIN)) THEN
        CMMIN = 0.0
       ELSE
        CMMAX = 0.0
      ENDIF
      CALL AXISADJ2(CMMIN,CMMAX,CMSPAN,DCM,NCMTICS)
c      write(*,*) 'cmmin,cmmax ',cmmin,cmmax
c      write(*,*) 'dcm,ncmtics ',dcm,ncmtics
C--- ALFA axes
      ALMIN = MIN(0.0,ALMIN)
      CALL AXISADJ2(ALmin,ALmax,ALspan,dAL,nALtics)
      IF(ALMIN.EQ.0.0) ALMIN = -DAL
C
      CPOLPLF(2,ICL) = CLMAX
      CPOLPLF(1,ICL) = CLMIN
      CPOLPLF(3,ICL) = DCL
      CPOLPLF(2,ICD) = CDMAX
      CPOLPLF(1,ICD) = CDMIN
      CPOLPLF(3,ICD) = DCD
      CPOLPLF(2,ICM) = CMMAX
      CPOLPLF(1,ICM) = CMMIN
      CPOLPLF(3,ICM) = DCM
      CPOLPLF(2,IAL) = ALMAX
      CPOLPLF(1,IAL) = ALMIN
      CPOLPLF(3,IAL) = DAL
C
      RETURN
      END ! SETINC


      SUBROUTINE SETINCV
      INCLUDE 'PPLOT.INC'
C--------------------------------------------
C     Determines axes increments for polars
C     from quantities for all polars read in.
C--------------------------------------------
C
      VHMAX = VPOLPLF(2,1)
      VHMIN = VPOLPLF(1,1)
      VZMAX = VPOLPLF(2,2)
      VZMIN = VPOLPLF(1,2)
C
C---- V axes
      CALL AXISADJ2(VHMIN,VHMAX,VHSPAN,DVH,NVHTICS)
C
C---- Vz axes
      VZMIN = 0.0
      CALL AXISADJ2(VZMIN,VZMAX,VZSPAN,DVZ,NVZTICS)
C
      VPOLPLF(2,1) = VHMAX
      VPOLPLF(1,1) = VHMIN
      VPOLPLF(3,1) = DVH
      VPOLPLF(2,2) = VZMAX
      VPOLPLF(1,2) = VZMIN
      VPOLPLF(3,2) = DVZ
C
      RETURN
      END ! SETINCV


      subroutine AXISADJ2(xmin,xmax,xspan,deltax,ntics)
C...Make scaled axes with engineering increments between tics
C
C   Input:    xmin, xmax - input range for which scaled axis is desired
C
C   Output:   xmin, xmax - adjusted range for scaled axis
C             xspan      - adjusted span of scaled axis
C             deltax     - increment to be used for scaled axis
C             nincr      - number of tics to be used on axis
C                          note that ntics=1+(xspan/deltax)
C
      real    xmin,xmax,xspan,deltax,xinc,xinctbl(4)
      integer ntics,i
      data    xinctbl / 0.1, 0.2, 0.5, 1. /
c
      xspan1 = xmax-xmin
      if (xspan1.eq.0.) xspan1 = 1.
c
      xpon = ifix(log10(xspan1))
      xspan = xspan1 / 10.**xpon
c
      do i = 1, 4
        xinc = xinctbl(i)
        ntics = 1 + ifix(xspan/xinc + 0.1)
        if (ntics.LE.6) go to 1
      end do
c
   1  deltax = xinc*10.**xpon
      xmin = deltax*  ifloor2(xmin/deltax)
      xmax = deltax*iceiling2(xmax/deltax)
      xspan = xmax - xmin
      ntics = 1 + ifix(xspan/xinc + 0.1)
      return
      end

      function iceiling2(x)
c--- returns next highest integer value if fraction is non-zero 
      integer iceiling2
      real x
      i = ifix(x)
      if(x-i.GT.0.) i = i+1
      iceiling2 = i
      return
      end

      function ifloor2(x)
c--- returns next lowest integer value if fraction is negative, non-zero
      integer ifloor2
      real x
      i = ifix(x)
      if(x-i.LT.0.) i = i-1
      ifloor2 = i
      return
      end

