*     FPL file

      SUBROUTINE FPL
**********************************************************************
*
*     SUBROUTINE FPL
*     --------------
*     PRINTER PLOTS OF DATA Y,X,... ON ONE PAGE (60 * 120  UNITS)
*                                                50 * 120  UNITS)
*
*     ADDITION OF DATA TO PLOT BUFFER (IN ANY ORDER)
*     ----------------------------------------------
*
*     CH = CHARACTER, E.G. '*'
*
*     CALL FPXY(CH,X,Y)               ADD POINT X,Y
*
*     CALL FPNXY(CH,N,X,Y)            ADD POLYGON X(1),Y(1)  X(2),Y(2)
*                                     ... (X(N),Y(N))
*                                     A SEQUENCE OF CALLS WITH N=1 AND
*                                     IDENTICAL CHARACTER CH WILL ALSO
*                                     BE PLOTTED AS A POLYGON. IF CH=' '
*                                     THE X Y VALUES ARE TAKEN INTO
*                                     ACCOUNT IN LIMIT CALCULATION,
*                                     BUT THE POINT IS NOT PLOTTED.
*
*     CALL FPXYD(CH,X,Y,D)            ADD POINT X, Y+-DY
*
*     CALL FPXDYD(CH,X,DX,Y,DY)       ADD POINT X+-DX, Y+-DY
*
*     CALL FPC(TEXT)                  TEXT=CHARACTER STRING, E.G. 'MASS'
*                                     TO BE PRINTED ABOVE PLOT.
*                                     SEVERAL CALLS POSSIBLE, 132 CHARS
*                                     IN TOTAL ACCEPTED.
*     CALL FLX(X1,X2)                 LIMITS X1 AND X2
*
*     CALL FLY(Y1,Y2)                 LIMITS Y1 AND Y2
*
*
*     PRINTOUT
*     --------
*
*     CALL FPL                        PLOT DATA AND RESET BUFFER
*
*     CALL FPS                        RESET DATA (WITHOUT PLOTTING)
*
*     CALL FPLUN(LUN)                 CHANGE PRINT-UNIT TO LUN
*                                     LUN = 0 MEANS SUPPRESS PRINTOUT
*
*     ALL  ARGUMENTS  ARE  INPUT  ARGUMENTS.   SCALES   ARE   DETERMINED
*     AUTOMATICALLY, TO INCLUDE ALL DATA OR THE LIMITS CAN BE DEFINED
*     BY USER CALLS (FLX, FLY).           FPL USES A  BUFFER  WITH  1602
*     WORDS, WHICH IS SUFFICIENT FOR 800 COORDINATES OF ONE  POLYGON  OR
*     400 COORDINATE PAIRS FOR ONE CALL OF FPNXY.
*
*     EXAMPLE: PLOT OF FUNCTION EXP(-X) FOR X BETWEEN 0 AND 3
*
*         CALL FPC('PLOT OF FUNCTION EXP(-X)')
*         DO 10 I=1,31
*         X=0.1*FLOAT(I-1)
*         Y=EXP(-X)
*      10 CALL FPNYX('*',1,Y,X)
*         CALL FPL
**********************************************************************

*     PARAMETER
      PARAMETER (LINENR = 50)
*     ARGUMENTS
         REAL X(*),Y(*)
         CHARACTER*1 CH
         CHARACTER*(*) TEXT
*     STORAGE FOR CHARACTERS
         CHARACTER*120 LINE,SLINE
         CHARACTER*1   C
         COMMON/CHACOM/LINE(LINENR),C
         EQUIVALENCE (SLINE,LINE(1))
*     ARRAYS FOR SCALE NUMBERS
         REAL RNUM(13)
         CHARACTER*8 NUMB(12)
         CHARACTER*4 CHXP
*     STORAGE FOR COMMENT
         CHARACTER*132 COMENT/' '/
*     STORAGE FOR COORDINATES
         COMMON/FICOM/IST(1602)
         REAL         RST(1602)
         EQUIVALENCE (IST(1),RST(1))
         COMMON/CFPL1/XA,XB,YA,YB,XLA,XLB,YLA,YLB
         DATA    ILIM/1602/,I/0/,ILAST/0/,ILX/0/,ILY/0/,ILZ/0/
         DATA IC/0/,ICLAST/0/,LUP/6/,IM/0/

*     ...

      IF(IM.EQ.0.OR.LUP.EQ.0) GOTO 60
      CALL FPL3(I,ILX,ILY,ILZ)
*     COMMENT LINE
      IF(IC.EQ.0) COMENT=' '
      WRITE(LUP,101) COMENT
      IC=0
*
      DO  2 J=1,13
    2 RNUM(J)=(XLA*FLOAT(13-J)+XLB*FLOAT(J-1))/12.0
      CALL CFMT(RNUM,12,NUMB,CHXP)
      WRITE(LUP,105) (NUMB(J),J=1,12),CHXP
      DO 4 J=1,LINENR/10+1
    4 RNUM(J)=(YLB*FLOAT(LINENR/10+1-J)+YLA*FLOAT(J-1))/FLOAT(LINENR/10)
      CALL CFMT(RNUM,LINENR/10+1,NUMB,CHXP)
      WRITE(LUP,102) CHXP
      CHXP=' '
*     PRINT THE PLOT ARRAY
      DO 40 J=1,LINENR
      IF(MOD(J,10).EQ.1) THEN
         WRITE(LUP,103) NUMB(1+J/10),LINE(J)
      ELSE
         C='I'
         IF(MOD(J,10).EQ.6) C='T'
         WRITE(LUP,104) C,LINE(J),C
      END IF
   40 CONTINUE
      WRITE(LUP,102) CHXP
      DO 55 J=1,13
   55 RNUM(J)=(XLA*FLOAT(13-J)+XLB*FLOAT(J-1))/12.0
      CALL CFMT(RNUM,12,NUMB,CHXP)
      WRITE(LUP,105) (NUMB(J),J=1,12),CHXP
*     RESET BUFFER
      ENTRY FPS
   60 ILX=0
      ILY=0
      ILZ=0
      IM=0
      IC=0
      I =0
      GOTO 100

*     STORE COMMENT

      ENTRY FPC(TEXT)
      IF(IC.EQ.0) COMENT=' '
      M=LEN(TEXT)
      IA=IC+1
      IC=MIN0(132,IC+M)
      IF(IA.GT.IC) GOTO 100
      COMENT(IA:IC)=TEXT
      IA=IC+1
      IC=MIN0(132,IC+3)
      IF(IA.GT.IC) GOTO 100
      COMENT(IA:IC)='   '
      GOTO 100

*     STORE SINGLE PAIR X,Y

      ENTRY FPXY(CH,X,Y)
      N=1

*     STORE USER Y LIMITS

      LCN=0
      IF(I+6.GT.ILIM) CALL FPL3(I,ILX,ILY,ILZ)
      GOTO 85

*     STORE NN PAIRS X,Y

      ENTRY FPNXY(CH,NN,X,Y)
      IF(NN.LE.0) GOTO 100
      N=NN
      IF(I+4+2*N.GT.ILIM) CALL FPL3(I,ILX,ILY,ILZ)
      IF(N.NE.1.OR.LCN.NE.1)     GOTO 80
      IF(II.GT.I) GOTO 80
      IF(IST(II-1).NE.ICHAR(CH)) GOTO 80
      IF(I+2.GT.ILIM) GOTO 100
*     APPEND TO PREVIOUS DATA
      RST(I+1)=X(1)
      RST(I+2)=Y(1)
      I=I+2
      IST(II)=IST(II)+1
      XA=MIN(XA,X(1))
      XB=MAX(XB,X(1))
      YA=MIN(YA,Y(1))
      YB=MAX(YB,Y(1))
      RETURN

   80 LCN=NN
   85 IF(IM.EQ.0) THEN
         YA=Y(1)
         XA=X(1)
         XB=X(1)
         YB=Y(1)
         IM=1
      END IF
      IF(I+4.GT.ILIM) GOTO 100
      IF(CH.EQ.' ') GOTO 100
*     STORAGE OF X,Y PAIRS
*     ICHAR / NR OF PAIRS / X / Y / X / Y / ...

      IST(I+1)=ICHAR(CH)
      IST(I+2)=0
      I=I+2
      II=I
      DO 90 J=1,N
      IF(I+2.GT.ILIM) GOTO 100
      RST(I+1)=X(J)
      RST(I+2)=Y(J)
      I=I+2
      IST(II)=IST(II)+1
      XA=MIN(XA,X(J))
      XB=MAX(XB,X(J))
      YA=MIN(YA,Y(J))
   90 YB=MAX(YB,Y(J))
      IF(LCN.EQ.1) RETURN
      GOTO 100

*     STORE PAIR X,Y WITH Y ERROR BAR

      ENTRY FPXYD(CH,X,Y,DY)
      IF(IM.EQ.0) THEN
         XA=X(1)
         XB=X(1)
         YA=Y(1)
         YB=Y(1)+DY
         IM=1
      END IF
      IF(I+5.GT.ILIM) CALL FPL3(I,ILX,ILY,ILZ)
*     STORAGE
*     ICHAR / -3 / X / Y / DY /

      IST(I+1)=ICHAR(CH)
      IST(I+2)=-3
      RST(I+3)=X(1)
      RST(I+4)=Y(1)
      RST(I+5)=DY
      I=I+5
      XA=MIN(XA,X(1))
      XB=MAX(XB,X(1))
      YA=MIN(YA,Y(1)   )
      YB=MAX(YB,Y(1)+DY)
      GOTO 100

*     STORE PAIR X,Y WITH X HALF BIN WIDTH AND Y ERROR BAR

      ENTRY FPXDYD(CH,X,DX,Y,DY)
      IF(IM.EQ.0) THEN
         XA=X(1)-DX
         XB=X(1)+DX
         YA=Y(1)
         YB=Y(1)+DY
         IM=1
      END IF
      IF(I+6.GT.ILIM) CALL FPL3(I,ILX,ILY,ILZ)

*     STORAGE
*     ICHAR / -4 / X / DX / Y / DY /

      IST(I+1)=ICHAR(CH)
      IST(I+2)=-4
      RST(I+3)=X(1)
      RST(I+4)=DX
      RST(I+5)=Y(1)
      RST(I+6)=DY
      I=I+6
      XA=MIN(XA,X(1)-DX)
      XB=MAX(XB,X(1)+DX)
      YA=MIN(YA,Y(1)   )
      YB=MAX(YB,Y(1)+DY)
      GOTO 100

*     STORE USER X LIMITS

      ENTRY FLX(X1,X2)
      XLA=MIN(X1,X2)
      XLB=MAX(X1,X2)
      IF(XLA.EQ.XLB) XLB=XLA+1.0
      ILX=1
      RETURN
      ENTRY FLY(Y1,Y2)
      YLA=MIN(Y1,Y2)
      YLB=MAX(Y1,Y2)
      IF(YLA.EQ.YLB) YLB=YLA+1.0
      ILY=1
      RETURN
C
  100 LCN=0
      RETURN
      ENTRY FPLUN(LUN)
      LUP=LUN
      RETURN
  101 FORMAT('1',A132/)
  102 FORMAT(5X,A4,1X,'+',12('L----L----'),'+')
  103 FORMAT(1X,A8,1X,'T',A120,'T')
  104 FORMAT(10X,A1,A120,A1)
  105 FORMAT(4X,12(2X,A8),A4)
      END
      SUBROUTINE FPL1(X1,Y1)

*     PLOT POINT X,Y

      PARAMETER (LINENR = 50)
      COMMON/CFPL1/XA,XB,YA,YB,XLA,XLB,YLA,YLB
      CHARACTER*120 LINE
      CHARACTER*1 C
      COMMON/CHACOM/LINE(LINENR),C
C     ...
      RJ=1.0+ FLOAT(LINENR)*(Y1-YLA)/(YLB-YLA)
      IF(RJ.LT.0.99.OR.RJ.GT.FLOAT(LINENR)+1.01) GOTO 100
      LJ=RJ
      LJ=MIN0(LINENR,MAX0(1,LJ))
      RI=1.0+120.0*(X1-XLA)/(XLB-XLA)
      IF(RI.LT.0.99.OR.RI.GT.121.01) GOTO 100
      LI=RI
      LI=MIN0(120,MAX0(1,LI))
      LINE(LINENR+1-LJ)(LI:LI)=C
      GOTO 100

*     PLOT LINE X0,Y0 TO X2,Y2

      ENTRY FPL2(X0,Y0,X2,Y2)
      SI=2.0+120.0*(X2-X0)/(XLB-XLA)
      SJ=2.0+FLOAT(LINENR)*(Y2-Y0)/(YLB-YLA)
      KP=SQRT(SI*SI+SJ*SJ)+1.5
      IF(KP.GT.1000) KP=2
      DO 10 K=1,KP
      XX=(X0*FLOAT(KP-K)+X2*FLOAT(K-1))/FLOAT(KP-1)
      YY=(Y0*FLOAT(KP-K)+Y2*FLOAT(K-1))/FLOAT(KP-1)
      LJ=1.0+FLOAT(LINENR)*(YY-YLA)/(YLB-YLA)
      IF(LJ.LT.1.OR.LJ.GT.LINENR) GOTO 10
      LI=1.0+120.0*(XX-XLA)/(XLB-XLA)
      IF(LI.LT.1.OR.LI.GT.120) GOTO 10
      LINE(LINENR+1-LJ)(LI:LI)=C
   10 CONTINUE
C
  100 RETURN
      END
      SUBROUTINE FPL3(ILAST,ILX,ILY,ILZ)

*     CONVERT STORED NUMERICAL DATA TO PLOT BUFFER

      PARAMETER (LINENR = 50)
*     STORAGE FOR CHARACTERS
         CHARACTER*120 LINE,SLINE
         CHARACTER*1   C
         COMMON/CHACOM/LINE(LINENR),C
         EQUIVALENCE (SLINE,LINE(1))
*     STORAGE FOR COORDINATES
         COMMON/FICOM/IST(1602)
         REAL         RST(1602)
         EQUIVALENCE (IST(1),RST(1))
      COMMON/CFPL1/XA,XB,YA,YB,XLA,XLB,YLA,YLB
*     ...
*     SCALES IN X
      IF(ILX.EQ.0) THEN
         CALL ROUNDB(XA,XB,120,XLA,XLB)
         ILX=1
      END IF
*     SCALES IN Y
      IF(ILY.EQ.0) THEN
         CALL ROUNDB(YA,YB,LINENR,YLA,YLB)
         ILY=1
      END IF
*     BLANK PLOT BUFFER
      IF(ILZ.EQ.0) THEN
         ILZ=1
         DO 8 J=1,LINENR
    8    LINE(J)=' '
      END IF
*     LOOP OVER NUMERICAL DATA
      I=0
   10 IF(I.GE.ILAST) GOTO 30
      C=CHAR(IST(I+1))
      M=IST(I+2)
      I=I+2
      IF(M.EQ.(-3)) THEN
*        POINT WITH Y ERROR BAR
         C='I'
         CALL FPL2(RST(I+1),RST(I+2)-RST(I+3),
     1             RST(I+1),RST(I+2)+RST(I+3))
         C=CHAR(IST(I-1))
         CALL FPL1(RST(I+1),RST(I+2))
         I=I+3
      ELSE IF(M.EQ.(-4)) THEN
*        POINT WITH X BIN WIDTH AND Y ERROR BAR
         C='-'
         CALL FPL2(RST(I+1)-RST(I+2),RST(I+3),
     1             RST(I+1)+RST(I+2),RST(I+3))
         C='I'
         CALL FPL2(RST(I+1),RST(I+3)-RST(I+4),
     1             RST(I+1),RST(I+3)+RST(I+4))
         C=CHAR(IST(I-1))
         CALL FPL1(RST(I+1),RST(I+3))
         I=I+4
      ELSE IF(M.EQ.1) THEN
*        POINT
         CALL FPL1(RST(I+1),RST(I+2))
         I=I+2
      ELSE
*        POLYGON
         DO 20 J=1,M-1
         CALL FPL2(RST(I+1),RST(I+2),RST(I+3),RST(I+4))
   20    I=I+2
         I=I+2
      END IF
      GOTO 10
*     CLEAR NUMERICAL BUFFER
   30 ILAST=0
  100 RETURN
      END
      SUBROUTINE CFMT(X,N,XCHAR,ECHAR)
C
C     SUBROUTINE CFMT
C     ---------------
C     PREPARE PRINTOUT OF ARRAY OF REAL NUMBERS AS CHARACTER STRINGS
C
C                  - -
C        CALL CFMT(X,N,XCHAR,ECHAR)
C                      ----- -----
C     WHERE X( )     = ARRAY OF N REAL VALUES
C           XCHAR( ) = ARRAY OF N CHARACTER*8 VARIABLES
C           ECHAR    = CHARACTER*4 VARIABLE
C
C     CFMT CONVERTS AN ARRAY OF  N  REAL  VALUES  INTO  N  CHARACTER*  8
C     VARIABLES (CONTAINING THE VALUES AS TEXT) AND  A  COMMON  EXPONENT
C     FOR PRINTING. UNNECCESSARY ZEROS ARE SUPPRESSED.
C
C
C     EXAMPLE - X(1)=1200.0, X(2)=1700.0 WITH N=2 ARE CONVERTED TO
C               XCHAR(1)='  1.2   ', XCHAR(2)='  1.7   ', ECHAR='E 03'
C
C
      REAL X(*)
      CHARACTER*1 CH(10)/'0','1','2','3','4','5','6','7','8','9'/
      CHARACTER*4 ECHAR
      CHARACTER*8 XCHAR(N),FXF,SXF,NULL/'00000000'/
C     ...
C
C     DETERMINE FACTOR FC, SO THAT FC*XMAX IS 5 DIGIT NUMBER
C
      IP=0
      XM=0.0
      DO 10 I=1,N
   10 XM=AMAX1(XM,ABS(X(I)))
      IF(XM.NE.0.0) THEN
         JP=104-IFIX(ALOG10(ABS(XM))+100.04)
         FC=10.0**JP
      ELSE
         JP=5
         FC=1.0
      END IF
C
C     STORE DIGITS AS CHARACTERS AND FIND JM = FIRST NONZERO DIGIT
C
      IM=6
      DO 30 I=1,N
      FXF=NULL
      IJ=FC*ABS(X(I))+0.5
      JM=6
      IF(IJ.NE.0) THEN
         DO 20 J=1,5
         JN=MOD(IJ,10)
         IJ=IJ/10
         IF(JN.NE.0.AND.JM.EQ.6) JM=J
   20    FXF(J:J)=CH(JN+1)
         IM=MIN(IM,JM)
      END IF
   30 XCHAR(I)=FXF
      JM=IM
C
C     DETERMINE EXPONENT AS A MULTIPLE OF 3

   32 IF(JP.LT.1) THEN
         JP=JP+3
         IP=IP+3
         GOTO 32
      END IF
   34 IF(JP.GT.JM+4.OR.JP.GE.8) THEN
         JP=JP-3
         IP=IP-3
         GOTO 34
      END IF
C
C     LOOP TO CONVERT TO PRINT FORMAT
C
      JA=MIN(JM,JP)
      JB=MAX(6,JP+1)
      DO 90 I=1,N
      FXF=XCHAR(I)
      SXF=' '
      IB=7+(JB-JA)/2
      DO 80 J=JA,JB
      IF(FXF(J:J).NE.CH(1)) GOTO 70
      IF(J.GT.JP+1) GOTO 50
      IF(FXF.NE.NULL.OR.J.GE.JP) GOTO 70
      IB=IB-1
      GOTO 80
   50 DO 60 K=J,JB
      IF(FXF(K:K).NE.CH(1)) GOTO 70
   60 CONTINUE
      GOTO 80
C     INSERT DIGIT
   70 IB=IB-1
      SXF(IB:IB)=FXF(J:J)
      IF(J.EQ.JP) THEN
C     INSERT DECIMAL DOT
         IB=IB-1
         SXF(IB:IB)='.'
      END IF
   80 CONTINUE
C     INSERT - SIGN
      IF(X(I).LT.0.0) SXF(IB-1:IB-1)='-'
   90 XCHAR(I)=SXF
C
C     PREPARE PRINT FORMAT FOR EXPONENT
C
      ECHAR=' '
      IF(IP.NE.0) THEN
         ECHAR='E 0'
         IF(IP.LE.0) THEN
            ECHAR(2:2)='-'
            IP=IABS(IP)
         END IF
         J=MOD(IP,10)
         ECHAR(4:4)=CH(J+1)
         IP=IP-10*J
         IF(IP.NE.0) THEN
            J=MOD(IP,10)
            ECHAR(3:3)=CH(J+1)
         END IF
      END IF
  100 RETURN
      END
      SUBROUTINE HPRB(H,N)

*     determine bins

      REAL H(N),H2(N),B(121),C(128)
      INTEGER NP/0/
      CHARACTER*1 CH
      CALL HCOMPX(H,N,B,NP)
      GOTO 100

      ENTRY HPRH(H,H2,N,XA,XB)

*     add hsitogram

      CH='O'
      IF(NP.LE.0) GOTO 100
      BX=0.5*(XB-XA)/FLOAT(N)
      M=N
      IF(3*N.LE.120) M=3*N
      IF(M.LT.120) THEN
         XH=(XA*FLOAT(M-120)+XB*120.0)/FLOAT(M)
      ELSE
         XH=XB
      END IF
      XL=XA
      XU=XH
C     XL=(119.0*XA+XH)/120.0
C     XU=(XA+119.0*XH)/120.0
      CALL FLX(XL,XU)
      DO 20 I=1,NP
      X =XA+BX*(B(I)+B(I+1))
      DX=BX*(B(I+1)-B(I))
      Y =0.0
      DY=0.0
      JA=B(I)+1.5
      JB=B(I+1)+0.5
      DO 10 J=JA,JB
      Y =Y +H(J)
   10 DY=DY+H2(J)
      Y=Y/(B(I+1)-B(I))
      DY=SQRT(DY/(B(I+1)-B(I)))
   20 CALL FPXDYD(CH ,X,DX,Y,DY)
      GOTO 100

      ENTRY HPRC(H,N,XA,XB)

      DO 25 I=1,N
   25 C(I)=H(I)

*     add curve

      DO 30 I=1,N
      IF(C(I).NE.0.0) GOTO 32
   30 CONTINUE
      GOTO 100
   32 IA=I
      DO 34 I=N,1,-1
      IF(C(I).NE.0.0) GOTO 36
   34 CONTINUE
   36 IB=I
      CALL SPLFT(C,N)
      DO 38 I=IA,IB
      Y=C(I)
      X=XA+(FLOAT(I)-0.5)*(XB-XA)/FLOAT(N)
   38 CALL FPNXY('.',1,X,Y)
  100 RETURN
      END