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, LERR
      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...Try to read default file "pplot.def" for stored plot setup
      LU = 10
      OPEN(LU,FILE='pplot.def',STATUS='OLD',ERR=2)
      CALL RDDEF(LU,LERR)
      CLOSE(LU)
      IF(LERR) THEN
       WRITE(*,*)
       WRITE(*,*) 'Read error on file  pplot.def'
       WRITE(*,*) 'Using default settings'
       WRITE(*,*)
       CALL GETDEF
      ELSE
       WRITE(*,*)
       WRITE(*,*) 'Settings read from file  pplot.def'
       WRITE(*,*)
      ENDIF
      GO TO 3
C
 2    CONTINUE
      WRITE(*,*)
      WRITE(*,*) 'No  pplot.def  file found'
      WRITE(*,*) 'Using default settings'
      WRITE(*,*)
C
 3    CONTINUE
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 5
          ENDIF
        ENDIF
      END DO
C
    5 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=5) IOPTS
      IOPT = ABS(IOPTS)
C
      GO TO (900, 10, 20, 30, 40, 50, 60, 70, 80, 90, 900,
     &            10, 20,130,  5,  5,  5,  5,  5,190,   5 ), IOPT+1
      GO TO 5
C
C=============================================
C---- read polars
 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(1,IP),XTRIP(1,IP),
     &    PTRAT(IP),ETAP(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
        IF(ACRIT(1,IP).EQ.ACRIT(2,IP)) THEN
         WRITE(*,8030) ACRIT(1,IP)
        ELSE
         WRITE(*,8030) ACRIT(1,IP), ACRIT(2,IP)
        ENDIF
        IF(PTRAT(IP).NE.0.0) THEN
         WRITE(*,8040) PTRAT(IP)
         WRITE(*,8041) ETAP(IP)
        ENDIF
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 =', 10F6.2       )
 8040   FORMAT('          pi_p  =', F8.4,       $)
 8041   FORMAT('          eta_p =', F8.4         )
C
 105  CONTINUE
      IP = IP2+1
C
 108  CONTINUE
      NPOL = IP-1
      IP2 = MIN(IP2,NPOL)
C
      DO IP = IP1, IP2
        CALL STRIP(NAME(IP),NNAME)
cccc        CALL GETTYP(NAX,NA(IP),CPOL(1,1,IP),IMATYP(IP),IRETYP(IP))
cc        ICOL(IP) = 2 + IP
cc        ILIN(IP) = IP
      ENDDO
CCC   CALL MINMAX(NAX,NPOL,NA,CPOL,CPOLPLF)
Co
C---- are these dimensional polars?
      DO IP = IP1, IP2
        CALL GETCLEN(NAME(IP),CCLEN,NCLEN)
        IF(NCLEN.GT.0) THEN
         LCLEN = .TRUE.
         GO TO 5
        ENDIF
      ENDDO
      IF(.NOT.LPLOT) GO TO 30
      GO TO 5
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
        IF(LGETFN) THEN
         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
        ENDIF
C
ccc     IFCOL(ID) = NCOLOR - ID + 1
c       IFCOL(ID) = 2 + ID
c       IFSYM(ID) = MOD(ID,10)
 25   CONTINUE
 27   CONTINUE
      GO TO 5
C
C=============================================
C---- Make the CD(CL) Plot
 30   IF (NPOL.EQ.0 .AND. NDAT.EQ.0) GO TO 5
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,ICM,IP) = CPOL(IA,ICM,IP)
     &                     + DXMREF(IP)*CPOL(IA,ICL,IP)
C
          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 (LAUTO) 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,PTRAT,ETAP,
     &            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,LCDH,LCMDOT,
     &            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 5
C
C=============================================
C---- hardcopy output
 40   IF(LPLOT) CALL PLEND
      LPLOT = .FALSE.
      CALL REPLOT(IDEVRP)
      GO TO 5
C
C=============================================
C---- change settings
 50   CALL GETSET
      GO TO 5
C
C=============================================
C---- zoom
 60   CALL USETZOOM(.FALSE.,.TRUE.)
      CALL REPLOT(IDEV)
      GO TO 5
C
C=============================================
C---- unzoom
 70   CALL CLRZOOM
      CALL REPLOT(IDEV)
      GO TO 5
C
C=============================================
C---- annotate plot
 80   IF(.NOT.LPLOT) THEN
       WRITE(*,*) 'No active plot to annotate'
       GO TO 5
      ENDIF
      CALL ANNOT(CH)
      GO TO 5
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 5
C
C=============================================
C---- Make the Vz(V) Plot
 130  IF (NPOL.EQ.0 .AND. NDAT.EQ.0) GO TO 5
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 (LAUTO) 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, ISX,NBL,
     &            REYN,MACH,ACRIT,PTRAT,ETAP,
     &            NAME ,ICOL,ILIN,
     &            IMATYP,IRETYP,
     &            TITLE,CODE,VERSION,
     &            PLOTAR, CH,CH2, 
     &            LGRID,LLIST,LEGND,
     &            VPOLPLF)
      GO TO 5
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 5
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'
      LOGICAL LERR
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
C---- default polar line types and colors
C
C     1  *****************************  SOLID
C     2  **** **** **** **** **** ****  LONG DASHED
C     3  ** ** ** ** ** ** ** ** ** **  SHORT DASHED
C     4  * * * * * * * * * * * * * * *  DOTTED
C     5  ***** * ***** * ***** * *****  DASH-DOT
C     6  ***** * * ***** * * ***** * *  DASH-DOT-DOT
C     7  ***** * * * ***** * * * *****  DASH-DOT-DOT-DOT
C     8  **** **** * * **** **** * *    DASH-DASH-DOT-DOT
C
C     3  red
C     4  orange
C     5  yellow
C     6  green
C     7  cyan
C     8  blue
C     9  violet
C    10  magenta
C
      DO IP=1, NPX
ccc       ILIN(IP) = 1 + MOD(IP-1,8)
ccc       ICOL(IP) = 3 + MOD(IP-1,8)
C
C------ normally solid, going to dashed after IP=7
        ILIN(IP) = 1 + (IP-1)/7
C
C------ skip yellow (hard to see on white background)
        ICOL(IP) = 3 + MOD(IP-1,7)
        IF(ICOL(IP) .GE. 5) ICOL(IP) = ICOL(IP) + 1
      ENDDO

      DO ID=1, NDX
ccc        IFCOL(ID) = 3 + MOD(ID-1,8)
        IFSYM(ID) = MOD(ID-1,14)

C------ skip yellow (hard to see on white background)
        IFCOL(ID) = 3 + MOD(ID-1,7)
        IF(IFCOL(ID) .GE. 5) IFCOL(ID) = IFCOL(ID) + 1
      ENDDO
C
      LGRID = .TRUE.
      LCDW  = .FALSE.
      LLIST = .TRUE.
      LEGND = .TRUE.
      LCLIP = .FALSE.
      LCLEN = .FALSE.
      LAECEN = .FALSE.
      LCDH = .FALSE.
      LCMDOT = .FALSE.
C
C---- automatic scaling for axes
      LAUTO  = .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
C---- no CM location shift by default
      DO IP = 1, NPX
        DXMREF(IP) = 0.0
      ENDDO
C
      TITLE = '                                '
CCC            12345678901234567890123456789012
C
      RETURN
      END ! GETDEF


 
 
      SUBROUTINE RDDEF(LU,LERR)
C--- Read PPLOT plot parameters from save file
      INCLUDE 'PPLOT.INC'
      LOGICAL LERR
C
      CHARACTER*256 LINE
      LOGICAL LCOLH
C
      LCOLH = IDEVRP .EQ. 4
      SIZE0 = SIZE
C
 1000 FORMAT(A)
C
      READ(LU,*,ERR=90,END=80)
     &     CPOLPLF(1,ICL), CPOLPLF(2,ICL), CPOLPLF(3,ICL)
      READ(LU,*,ERR=90,END=80)
     &     CPOLPLF(1,ICD), CPOLPLF(2,ICD), CPOLPLF(3,ICD)
      READ(LU,*,ERR=90,END=80)
     &     CPOLPLF(1,ICM), CPOLPLF(2,ICM), CPOLPLF(3,ICM)
      READ(LU,*,ERR=90,END=80)
     &     CPOLPLF(1,IAL), CPOLPLF(2,IAL), CPOLPLF(3,IAL)
      READ(LU,*,ERR=90,END=80) XCD, XAL, XOC
      READ(LU,*,ERR=90,END=80) SIZE, PLOTAR
      READ(LU,*,ERR=90,END=80) CH, CH2
C
      READ(LU,*,ERR=90,END=80) LAUTO , LCDW 
      READ(LU,*,ERR=90,END=80) LLIST , LEGND
      READ(LU,*,ERR=90,END=80) LAECEN       
      READ(LU,*,ERR=90,END=80) LCMDOT, LCDH 
      READ(LU,*,ERR=90,END=80) LGRID , LCOLH
C
      READ(LU,1000,ERR=90,END=80) LINE
      KBAR = INDEX(LINE,'|') - 1
      IF(KBAR.LE.0) KBAR = LEN(LINE)
      READ(LINE(1:KBAR),*,ERR=90,END=80) (DXMREF(IP), IP=1, NPX)
C
      READ(LU,1000,ERR=90,END=80) LINE
      KBAR = INDEX(LINE,'|') - 1
      IF(KBAR.LE.0) KBAR = LEN(LINE)
      READ(LINE(1:KBAR),*,ERR=90,END=80) (ICOL(IP), IP=1, NPX)
C
      READ(LU,1000,ERR=90,END=80) LINE
      KBAR = INDEX(LINE,'|') - 1
      IF(KBAR.LE.0) KBAR = LEN(LINE)
      READ(LINE(1:KBAR),*,ERR=90,END=80) (ILIN(IP), IP=1, NPX)
C
C
      READ(LU,1000,ERR=90,END=80) LINE
      KBAR = INDEX(LINE,'|') - 1
      IF(KBAR.LE.0) KBAR = LEN(LINE)
      READ(LINE(1:KBAR),*,ERR=90,END=80) (IFCOL(ID), ID=1, NDX)
C
      READ(LU,1000,ERR=90,END=80) LINE
      KBAR = INDEX(LINE,'|') - 1
      IF(KBAR.LE.0) KBAR = LEN(LINE)
      READ(LINE(1:KBAR),*,ERR=90,END=80) (IFSYM(ID), ID=1, NDX)
C
C
      READ(LU,*,ERR=90,END=80) (VPOLPLF(K,1), K=1, 3)
      READ(LU,*,ERR=90,END=80) (VPOLPLF(K,2), K=1, 3)
C
C
      IF(LCOLH) THEN
       IDEVRP = 4
      ELSE
       IDEVRP = 2
      ENDIF
      IF(SIZE.LE.0.0) SIZE = SIZE0
C
      LERR = .FALSE.
      RETURN
C
 80   CONTINUE
 90   CONTINUE
      LERR = .TRUE.
      RETURN
      END ! RDDEF

 
      SUBROUTINE WRTDEF(LU)
C--- Write PPLOT plot parameters to save file
      INCLUDE 'PPLOT.INC'
      CHARACTER*256 LINE
      LOGICAL LCOLH
C
      LCOLH = IDEVRP .EQ. 4
C
 1000 FORMAT(A)
C
      WRITE(LU,1030) CPOLPLF(1,ICL), CPOLPLF(2,ICL), CPOLPLF(3,ICL),
     &               'CLmin   CLmax    dCL'
      WRITE(LU,1030) CPOLPLF(1,ICD), CPOLPLF(2,ICD), CPOLPLF(3,ICD),
     &               'CDmin   CDmax    dCD'
      WRITE(LU,1030) CPOLPLF(1,ICM), CPOLPLF(2,ICM), CPOLPLF(3,ICM),
     &               'CMmin   CMmax    dCM'
      WRITE(LU,1030) CPOLPLF(1,IAL), CPOLPLF(2,IAL), CPOLPLF(3,IAL),
     &               'ALmin   ALmax    dAL'
      WRITE(LU,1030) XCD, XAL, XOC,
     &               'CL-CD   CL-alpha  CL-Xtr  (widths)'
      WRITE(LU,1020) SIZE, PLOTAR,
     &               'width    height/width'
      WRITE(LU,1020) CH, CH2,
     &               'char_ht1 char_ht2'
 1010 FORMAT(1X, F9.4,9X,9X,' | ', A)
 1020 FORMAT(1X,2F9.4,9X   ,' | ', A)
 1030 FORMAT(1X,3F9.4      ,' | ', A)
C
      WRITE(LU,1120) LAUTO , LCDW  , 'auto_scale?  CDp_plot?'
      WRITE(LU,1120) LLIST , LEGND , 'airf_list?   legend_box?'
      WRITE(LU,1110) LAECEN,         'x_AC_plot?'
      WRITE(LU,1120) LCMDOT, LCDH  , 'HX_mass?     HX_CD?'
      WRITE(LU,1120) LGRID , LCOLH , 'grid_plot?   color_PS?'
 1110 FORMAT(1X,  L4,1X ,5X,5X,' | ', A)
 1120 FORMAT(1X,2(L4,1X),5X,   ' | ', A)
 1130 FORMAT(1X,3(L4,1X),      ' | ', A)
C
      WRITE(LINE,1300) (DXMREF(IP), IP=1, NPX)
 1300 FORMAT(1X,80(F7.3))
      CALL STRIP(LINE,NLINE)
      LINE = ' ' // LINE(1:NLINE) // ' | dXmom_ref'
      CALL STRIP(LINE,NLINE)
      WRITE(LU,1000) LINE(1:NLINE)
C
      WRITE(LINE,1400) (ICOL(IP), IP=1, NPX)
      CALL STRIP(LINE,NLINE)
      LINE = ' ' // LINE(1:NLINE) // ' | line_color'
      CALL STRIP(LINE,NLINE)
      WRITE(LU,1000) LINE(1:NLINE)
C
      WRITE(LINE,1400) (ILIN(IP), IP=1, NPX)
      CALL STRIP(LINE,NLINE)
      LINE = ' ' // LINE(1:NLINE) // ' | line_type'
      CALL STRIP(LINE,NLINE)
      WRITE(LU,1000) LINE(1:NLINE)
C
C
      WRITE(LINE,1400) (IFCOL(ID), ID=1, NDX)
      CALL STRIP(LINE,NLINE)
      LINE = ' ' // LINE(1:NLINE) // ' | data_point_color'
      CALL STRIP(LINE,NLINE)
      WRITE(LU,1000) LINE(1:NLINE)
C
      WRITE(LINE,1400) (IFSYM(ID), ID=1, NDX)
      CALL STRIP(LINE,NLINE)
      LINE = ' ' // LINE(1:NLINE) // ' | data_point_type'
      CALL STRIP(LINE,NLINE)
      WRITE(LU,1000) LINE(1:NLINE)
C
C
 1400 FORMAT(1X,80I4)
C
      WRITE(LU,1030) (VPOLPLF(K,1), K=1, 3), 'V_scale'
      WRITE(LU,1030) (VPOLPLF(K,2), K=1, 3), 'Vz_scale'
C
      RETURN
      END ! WRTDEF
 


      SUBROUTINE GETSET
      INCLUDE 'PPLOT.INC'
      LOGICAL OK, ERROR
      LOGICAL LCOLH
      CHARACTER*1 ANS
      CHARACTER*2 OPTION
      CHARACTER*80 LINE
      REAL RINP(10)
C
C---- Change plotting parameters
C
    1 CONTINUE
      LCOLH = IDEVRP .EQ. 4
C
      WRITE(*,1000)
     &    LAUTO,
     &    LCDW,
     &    LLIST,
     &    LEGND,
     &    LAECEN,
     &    LCMDOT,
     &    LCDH,
     &    LGRID,
     &    LCLIP,
     &    LCOLH
C
 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'
     &      // ' 10',L3,' autoscaling?'
     &       / ' 11',L3,' plot pressure-CD?'
     &       / ' 12',L3,' plot airfoil list?'
     &       / ' 13',L3,' plot CL-CD legend box?'
     &       / ' 14',L3,' plot aero. center?'
     &       / ' 15',L3,' plot streamtube mass coeff.?'
     &       / ' 16',L3,' plot streamtube thrust?'
     &       / ' 17',L3,' plot grid overlay?'
     &       / ' 18',L3,' clip to plot boundaries?'
     &       / ' 19',L3,' color hardcopy?'
     &      // ' 20   Rescale forces by chord factor'
     &       / ' 21   Change reference-length unit'
     &       / ' 22   Change moment-reference x/c'
     &       / ' 23   Change polar line colors'
     &       / ' 24   Change polar line styles'
     &       / ' 25   Change data point colors'
     &       / ' 26   Change data point shapes'
     &      // ' 27   Change V  scaling'
     &       / ' 28   Change Vz scaling'
     &      // ' 30   Read  settings from defaults file'
     &       / ' 31   Write settings to   defaults file'
     &      // '    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
        LAUTO = .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
        LAUTO = .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
        LAUTO = .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
        LAUTO = .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.'10') THEN
        LAUTO = .NOT. LAUTO
C
      ELSE IF(OPTION.EQ.'11') THEN
        LCDW = .NOT. LCDW
C
      ELSE IF(OPTION.EQ.'12') THEN
        LLIST = .NOT. LLIST
C
      ELSE IF(OPTION.EQ.'13') THEN
        LEGND = .NOT. LEGND
C
      ELSE IF(OPTION.EQ.'14') THEN
        LAECEN = .NOT. LAECEN
C
      ELSE IF(OPTION.EQ.'15') THEN
        LCMDOT = .NOT. LCMDOT
C
      ELSE IF(OPTION.EQ.'16') THEN
        LCDH = .NOT. LCDH
C
      ELSE IF(OPTION.EQ.'17') THEN
        LGRID = .NOT. LGRID
C
      ELSE IF(OPTION.EQ.'18') THEN
        LCLIP = .NOT. LCLIP
C
      ELSE IF(OPTION.EQ.'19') THEN
C--- Color hardcopy toggle
        IF(IDEVRP.EQ.2) THEN
          IDEVRP = 4
         ELSE
          IDEVRP = 2
        ENDIF
C
      ELSE IF(OPTION.EQ.'20') 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.'21') 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.'22') THEN
C--- change moment reference locations
        IF(NPOL.EQ.0) THEN
         WRITE(*,*) 'No current polars'
         GO TO 1
        ELSE
         WRITE(*,*)
         WRITE(*,*) 'Enter new moment-reference location shifts...'
         DO IP = 1, NPOL
           WRITE(*,8010) IP, DXMREF(IP)
 8010      FORMAT('  New dXref for polar', I3,' [',F9.4,' ] : ', $)
           CALL READR(1,DXMREF(IP),ERROR)
         ENDDO
        ENDIF
C
      ELSE IF(OPTION.EQ.'23') THEN
C------ change polar colors
        IF(NPOL.EQ.0) THEN
         WRITE(*,*) 'No current polars to change'
         GO TO 1
        ELSE
         WRITE(*,5020)
 5020    FORMAT(
     &   / '   1  black (white in revVideo)'
     &   / '   2  white (invisible)'
     &   / '   3  red'
     &   / '   4  orange'
     &   / '   5  yellow'
     &   / '   6  green'
     &   / '   7  cyan'
     &   / '   8  blue'
     &   / '   9  violet'
     &   / '  10  magenta' )
C
 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.'24') THEN
C------ change polar line styles
        IF(NPOL.EQ.0) THEN
         WRITE(*,*) 'No current polars to change'
         GO TO 1
        ELSE
         WRITE(*,5030)
 5030    FORMAT(
     &   / '   1  -----------------------------  solid'
     &   / '   2  ---- ---- ---- ---- ---- ----  long dashed'
     &   / '   3  -- -- -- -- -- -- -- -- -- --  short dashed'
     &   / '   4  - - - - - - - - - - - - - - -  dotted'
     &   / '   5  ----- - ----- - ----- - -----  dash-dot'
     &   / '   6  ----- - - ----- - - ----- - -  dash-dot-dot'
     &   / '   7  ----- - - - ----- - - - -----  dash-dot-dot-dot'
     &   / '   8  ---- ---- - - ---- ---- - -    dash-dash-dot-dot')
C
 830     WRITE(LINE,3100) 'current polar line styles', 
     &                     (ILIN(IP), IP=1, NPOL)
         WRITE(*,1005) LINE
         WRITE(*,3105)    ' select 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.'25') THEN
C------ change data colors
        IF(NDAT.EQ.0) THEN
         WRITE(*,*) 'No current data to change'
         GO TO 1
        ELSE
         WRITE(*,5040)
 5040    FORMAT(
     &   / '   1  black (white in revVideo)'
     &   / '   2  white (invisible)'
     &   / '   3  red'
     &   / '   4  orange'
     &   / '   5  yellow'
     &   / '   6  green'
     &   / '   7  cyan'
     &   / '   8  blue'
     &   / '   9  violet'
     &   / '  10  magenta' )
C
 840     WRITE(LINE,3100) 'data colors',
     &                    (IFCOL(IP), IP=1, NDAT)
         WRITE(*,1005) LINE
         WRITE(*,3105)    'data colors'
         READ(*,1005) LINE
         NINP = NDAT
         CALL GETINT(LINE,IFCOL,NINP,ERROR)
         IF(ERROR) GO TO 840
        ENDIF
C
      ELSE IF(OPTION.EQ.'26') THEN
C------ change data symbols
        IF(NDAT.EQ.0) THEN
         WRITE(*,*) 'No current data to change'
         GO TO 1
        ELSE
         WRITE(*,5050)
 5050    FORMAT(
     &   / '  0  square         7  Y         '
     &   / '  1  circle         8  flipped Y '
     &   / '  2  triangle       9  *         '
     &   / '  3  +             10  flipped * '
     &   / '  4  x             11  hourglass '
     &   / '  5  diamond       12  bowtie    '
     &   / '  6  yield sign    13  star      ' )
C
 850     WRITE(LINE,3100) 'data point symbols',
     &                    (IFSYM(IP), IP=1, NDAT)
         WRITE(*,1005) LINE
         WRITE(*,3105)    'data point symbols'
         READ(*,1005) LINE
         NINP = NDAT
         CALL GETINT(LINE,IFSYM,NINP,ERROR)
         IF(ERROR) GO TO 850
        ENDIF
C
      ELSE IF(OPTION.EQ.'27') 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
        LAUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'28') 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
        LAUTO = .FALSE.
C
      ELSE IF(OPTION.EQ.'30') THEN
C--- Read defaults from pplot.def file
        LINE = 'Enter settings filename  [pplot.def] ^'
        CALL ASKS(LINE,FNAME)
        IF(FNAME.EQ.' ') FNAME = 'pplot.def'
        LU = 10
        OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=703)
        CALL RDDEF(LU,ERROR)
        CLOSE(LU)
        GO TO 1
 703    WRITE(*,*) 
        WRITE(*,*) 'Open error on pplot defaults file'
        GO TO 1
C
      ELSE IF(OPTION.EQ.'31') THEN
C--- Save defaults to parameter file
        LU = 10
        LINE = 'Enter settings filename  [ pplot.def ] ^'
        CALL ASKS(LINE,FNAME)
        IF(FNAME.EQ.' ') FNAME = 'pplot.def'
        OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=803)
        WRITE(*,*)
        WRITE(*,*) 'File exists.  Overwrite?  Y'
        READ(*,1001) ANS
        IF(INDEX('Nn',ANS) .EQ. 0) GO TO 806
        WRITE(*,*)
        WRITE(*,*) 'No action taken'
        CLOSE(LU)
        GO TO 1
C
 803    OPEN(LU,FILE=FNAME,STATUS='UNKNOWN')
 806    REWIND(LU)
        CALL WRTDEF(LU)
        WRITE(*,*)
        WRITE(*,*) 'PPLOT plot settings written to file'
        CLOSE(LU)
C
      ENDIF
      GO TO 1
C
 1001 FORMAT(A)
 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: ',$)
 1910 FORMAT(/'Enter moment reference x/c [', F8.3, ' ]: ',$)
 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
 3000 FORMAT(' Currently ',A,' =', 20F8.4)
 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
        DXMREF(IP) = DXMREF(IP)*SCAL
      END DO
C
      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

