C-----------------------------------------------------
C
C  Program name:  PPLOT
C
C  Date:          March 1988
C
C  Revised:       August 1988
C
C-----------------------------------------------------
C
      REAL*8      DATA(100), SDATA(100)
      INTEGER     N, SEL
      CHARACTER*8 XLABEL
C
      WRITE(*,100)
  100 FORMAT(////,32X,'PPLOT',
     +///,20X,'A Probability Plotting Program')
C
      CALL INPUT(N,DATA,SDATA,XLABEL)
C
  101 WRITE(*,103)
  103 FORMAT(8(/),1X,'CHOOSE OPTION FOR EDITING OR ',
     +'PLOTTING DATA:  ',
     +////,5X,' (1) ENTER DATA FROM KEYBOARD OR FILE',
     +/,5X,' (2) EDIT DATA',
     +/,5X,' (3) SAVE DATA',
     +/,5X,' (4) NORMAL PLOT ',
     +/,5X,' (5) TWO-PARAMETER LOGNORMAL PLOT ',
     +/,5X,' (6) THREE-PARAMETER LOGNORMAL PLOT',
     +/,5X,' (7) RIGHT-TAIL HALF-NORMAL PLOT ',
     +/,5X,' (8) LEFT-TAIL HALF-NORMAL PLOT ',
     +/,5X,' (9) EXPONENTIAL PLOT ',
     +/,5X,'(10) TWO-PARAMETER WEIBULL PLOT ',
     +/,5X,'(11) THREE-PARAMETER WEIBULL PLOT ',
     +/,5X,'(12) EXIT FROM THE PROGRAM',
     +////,1X,'ENTER NUMBER OF SELECTION:  ',/)
C
      READ(*,*) SEL
      IF (SEL.GE.1.AND.SEL.LE.12) GOTO 107
C
      WRITE(*,105)
  105 FORMAT(//,1X,'INVALID SELECTION.  PLEASE TRY ',
     +'AGAIN')
      GOTO 101
C
  107 GOTO (110,120,130,140,150,160,
     +      170,180,190,200,210,220), SEL
C
  110 CALL INPUT(N,DATA,SDATA,XLABEL)
      GOTO 101
C
  120 CALL EDIT(N,DATA,SDATA)
      GOTO 101
C
  130 CALL SAVE(N,SDATA,XLABEL)
      GOTO 101
C
  140 CALL NPLOT(N,SDATA,XLABEL)
      GOTO 101
C
  150 CALL LN2PLT(N,SDATA,XLABEL)
      GOTO 101
C
  160 CALL LN3PLT(N,SDATA,XLABEL)
      GOTO 101
C
  170 CALL HRPLOT(N,SDATA,XLABEL)
      GOTO 101
C
  180 CALL HLPLOT(N,SDATA,XLABEL)
      GOTO 101
C
  190 CALL EPLOT(N,SDATA,XLABEL)
      GOTO 101
C
  200 CALL W2PLOT(N,SDATA,XLABEL)
      GOTO 101
C
  210 CALL W3PLOT(N,SDATA,XLABEL)
      GOTO 101
C
  220 STOP
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE INPUT(N,DATA,SDATA,XLABEL)
C
      INTEGER     N, SEL
      REAL*8      DATA(100), SDATA(100)
      CHARACTER*8 XLABEL
C
  201 WRITE(*,203)
  203 FORMAT(8(/),1X,'CHOOSE OPTION FOR ENTERING DATA:  ',
     +////,5X,' (1) ENTER DATA FROM KEYBOARD',
     +/,5X,' (2) FETCH STORED DATA FILE ',
     +////,1X,'ENTER NUMBER OF SELECTION:  ',/)
C
      READ(*,*) SEL
      IF (SEL.GE.1.AND.SEL.LE.2) GOTO 207
C
      WRITE(*,205)
  205 FORMAT(//,1X,'INVALID SELECTION.  PLEASE TRY ',
     +'AGAIN')
      GOTO 201
C
  207 GOTO (210,220), SEL
C
  210 CALL KEYBD(N,DATA,XLABEL)
      GOTO 230
C
  220 CALL FETCH(N,DATA,XLABEL)
C
  230 CALL PRDATA(N,DATA)
      CALL PAUSE
C
      CALL SORT(N,DATA,SDATA)
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE KEYBD(N,DATA,XLABEL)
C
      REAL*8      DATA(100)
      INTEGER     N
      CHARACTER*8 XLABEL
C
  300 WRITE(*,305)
  305 FORMAT(////,' ENTER NUMBER OF DATA VALUES:  ',/)
      READ(*,*) N
      IF (N.GE.1.AND.N.LE.100) GOTO 315
C
      WRITE(*,310)
  310 FORMAT(//,1X,'THE NUMBER OF DATA VALUES MUST ',
     +'BE BETWEEN 1 AND 100.',
     +/,1X,'PLEASE TRY AGAIN.')
      GOTO 300
C
  315 WRITE(*,320)
  320 FORMAT(//,' ENTER NAME FOR DATA VARIABLE:  ',/)
      READ(*,325) XLABEL
  325 FORMAT(A8)
C
      WRITE(*,330)
  330 FORMAT(//,1X,'ENTER DATA VALUES, ',
     +'LEAVING A BLANK BETWEEN ENTRIES:')
      READ(*,*) (DATA(I),I=1,N)
C
      RETURN
      END
C
C
C----------------------------------------------------
C
      SUBROUTINE FETCH(N,DATA,XLABEL)
C
      INTEGER     N
      REAL*8      DATA(100)
      CHARACTER*8 XLABEL, FILENAME
C
      WRITE(*,61)
   61 FORMAT(//,' PLEASE ENTER DATA FILE NAME:  ',/)
      READ(*,63) FILENAME
   63 FORMAT(A8)
      IF (FILENAME.EQ.' ') GO TO 201
      OPEN (1,FILE=FILENAME,STATUS='OLD')
         READ(1,191) XLABEL
  191    FORMAT(A8)
         READ(1,192) N
  192    FORMAT(I3)
         READ(1,*) DATA
         CLOSE(1)
C
  201 RETURN
      END
C
C
C----------------------------------------------------
C
      SUBROUTINE PRDATA(N,DATA)
C
      INTEGER N
      REAL*8 DATA(100)
C
      WRITE(*,400)
  400 FORMAT(////,' OBSERVATION NUMBER',10X,'VALUE',/)
C
      DO 410 I=1,N
        WRITE(*,405) I,DATA(I)
  405   FORMAT(8X,I3,15X,F10.4)
  410 CONTINUE
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE EDIT(N,DATA,SDATA)
C
      INTEGER N,CH1
      REAL*8 DATA(100),SDATA(100)
C
  500 WRITE(*,505)
  505 FORMAT(////,' SELECT OPTION FOR EDITING DATA:',
     +//,5X,'(1) DELETE OBSERVATIONS',
     +/,5X,'(2) ADD OBSERVATIONS',
     +/,5X,'(3) CHANGE OBSERVATIONS',
     +/,5X,'(4) STOP EDITING',
     +//,1X,'ENTER NUMBER OF SELECTION:  ',/)
      READ(*,*) CH1
      IF (CH1.GE.1.AND.CH1.LE.4) GOTO 515
      WRITE(*,510)
  510 FORMAT(/,'INVALID RESPONSE.  PLEASE TRY AGAIN')
      GOTO 500
C
  515 GOTO (520,525,530,540),CH1
C
  520 CALL DELETE(N,DATA)
      GOTO 500
C
  525 CALL ADD(N,DATA)
      GOTO 500
C
  530 CALL CHANGE(N,DATA)
      GOTO 500
C
  540 CALL SORT(N,DATA,SDATA)
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE DELETE(N,DATA)
C
      INTEGER N,NUM,FLAG,INDEX(100),SINDEX(100)
      REAL*8 DATA(100)
C
      CALL PRDATA(N,DATA)
  600 WRITE(*,605)
  605 FORMAT(//,1X,'HOW MANY OF THESE OBSERVATIONS ',
     +'DO YOU WISH TO DELETE?   (0=NONE)  '/)
      READ(*,*) NUM
      IF (NUM.GE.1.AND.NUM.LE.N) GOTO 615
C
      IF (NUM .EQ. 0) GOTO 650
C
      WRITE(*,610)
  610 FORMAT(//,1X,'INVALID RESPONSE.  ',
     +'PLEASE TRY AGAIN')
      GOTO 600
C
  615 WRITE(*,620)
  620 FORMAT(//,1X,'ENTER THE OBSERVATION NUMBERS ',
     +'OF THOSE YOU WISH TO DELETE.',
     +/,1X,'LEAVE A BLANK BETWEEN ENTRIES.')
      READ(*,*) (INDEX(I),I=1,NUM)
      CALL CHECK(N,NUM,INDEX,FLAG)
      IF (FLAG.EQ.0) GOTO 630
C
      WRITE(*,625) N
  625 FORMAT(/,1X,'THESE NUMBERS MUST BE BETWEEN ',
     +'1 AND ',I3,'.',/,1X,'PLEASE TRY AGAIN.')
      GOTO 615
C
  630 CALL ISORT(NUM,INDEX,SINDEX)
      DO 640 I=1,NUM
        DO 635 J=SINDEX(I),N-1
          DATA(J-I+1)=DATA(J-I+2)
  635   CONTINUE
  640 CONTINUE
C
      N=N-NUM
C
      CALL PRDATA(N,DATA)
C
  650 RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE ADD(N,DATA)
C
      INTEGER N,NMORE
      REAL*8 DATA(100)
C
      CALL PRDATA(N,DATA)
  650 WRITE(*,655)
  655 FORMAT(//,1X,'HOW MANY OBSERVATIONS ',
     +'DO YOU WISH TO ADD?  (0=NONE)  '/)
      READ(*,*) NMORE
      IF ((NMORE+N).LE.100.AND.NMORE.GE.1) GOTO 665
C
      IF (NMORE .EQ. 0) GOTO 680
C
      WRITE(*,660)
  660 FORMAT(//,1X,'INVALID RESPONSE.  ',
     +'PLEASE TRY AGAIN')
      GOTO 650
C
  665 WRITE(*,670) NMORE
  670 FORMAT(//,1X,'ENTER THE ',I2,' NEW DATA '
     +'VALUES, LEAVING A BLANK BETWEEN ENTRIES:')
      READ(*,*) (DATA(I),I=N+1,N+NMORE)
C
      N=N+NMORE
C
      CALL PRDATA(N,DATA)
C
  680 RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE CHANGE(N,DATA)
C
      INTEGER N,NUM,FLAG,INDEX(100),SINDEX(100)
      REAL*8 DATA(100)
C
      CALL PRDATA(N,DATA)
  700 WRITE(*,705)
  705 FORMAT(//,1X,'HOW MANY OF THESE OBSERVATIONS ',
     +'DO YOU WISH TO CHANGE?  (0=NONE)  '/)
      READ(*,*) NUM
      IF (NUM.GE.1.AND.NUM.LE.N) GOTO 715
C
      IF (NUM .EQ. 0) GOTO 750
C
      WRITE(*,710)
  710 FORMAT(//,1X,'INVALID RESPONSE.  ',
     +'PLEASE TRY AGAIN')
      GOTO 700
C
  715 WRITE(*,720)
  720 FORMAT(//,1X,'ENTER THE OBSERVATION NUMBERS ',
     +'OF THOSE YOU WISH TO CHANGE.',
     +/,1X,'LEAVE A BLANK BETWEEN ENTRIES.')
      READ(*,*) (INDEX(I),I=1,NUM)
      CALL CHECK(N,NUM,INDEX,FLAG)
      IF (FLAG.EQ.0) GOTO 730
C
      WRITE(*,725) N
  725 FORMAT(/,1X,'THESE NUMBERS MUST BE BETWEEN ',
     +'1 AND ',I3,'.',/,1X,'PLEASE TRY AGAIN.')
      GOTO 715
C
  730 CALL ISORT(NUM,INDEX,SINDEX)
      WRITE(*,735)
  735 FORMAT(/,1X,'ENTER THE NEW DATA VALUES:',/)
      DO 745 I=1,NUM
        WRITE(*,740) SINDEX(I),DATA(SINDEX(I))
  740   FORMAT(1X,'OBSERVATION NUMBER ',I3,
     +':  OLD VALUE = ',F10.4,4X,'NEW VALUE = ',/)
        READ(*,*) DATA(SINDEX(I))
  745 CONTINUE
C
      CALL PRDATA(N,DATA)
C
  750 RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE SAVE(N,SDATA,XLABEL)
C
      REAL*8      SDATA(100)
      INTEGER     N
      CHARACTER*8 XLABEL, FILENAME
C
   70 WRITE(*,80)
   80 FORMAT(//,' PLEASE ENTER DATA FILE NAME (UP TO 8 CHARS):  ',/)

      READ(*,81) FILENAME
   81 FORMAT(A8)
      IF (FILENAME.EQ.' ') GO TO 70

      OPEN (1,FILE=FILENAME,STATUS='NEW')
        WRITE(1,91) XLABEL
   91   FORMAT(A8)
        WRITE(1,92) N
   92   FORMAT(I3)
        WRITE(1,*)  SDATA
        CLOSE(1)
C
      RETURN
      END
C
C
C----------------------------------------------------
C
      SUBROUTINE CHECK(N,NUM,INDEX,FLAG)
C
      INTEGER N,NUM,INDEX(100),FLAG
C
      FLAG=0
      DO 810 I=1,NUM
        IF (INDEX(I).GE.1.AND.INDEX(I).LE.N) GOTO 810
        FLAG=1
  810 CONTINUE
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE PAUSE
C
      CHARACTER*1 GO
C
C  Pause to let user view the display
C
      WRITE(*,70)
   70 FORMAT(/,1X,'STRIKE CARRIAGE RETURN WHEN YOU ',
     +'WISH TO CONTINUE',/)
      READ(*,75) GO
   75 FORMAT(A1)
C
      CALL WRLINE(2)
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE NPLOT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),ARG,PN,Y(100)
      CHARACTER*8 XLABEL
C
      WRITE(*,900)
  900 FORMAT(////,1X,'NORMAL PROBABILITY PLOT:',17X,
     +'Y',14X,'X',/)
C
      DO 910 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=PN(ARG)
        WRITE(*,905) Y(I),SDATA(I)
  905   FORMAT(35X,F10.4,5X,F10.4)
  910 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,920) XLABEL
  920 FORMAT(//,20X,'PLOT OF Y=NORMAL QUANTILE VS. ',
     +'X=',A8)
      CALL PLTTR(N,SDATA,Y)
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE LN2PLT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),Y(100),LDATA(100),ARG,PN
      CHARACTER*8 XLABEL
C
      IF (SDATA(1).GT.0) GOTO 2005
      WRITE(*,2000)
 2000 FORMAT(//,1X,'A TWO-PARAMETER LOGNORMAL PLOT ',
     +'IS NOT POSSIBLE',/
     +,1X,'BECAUSE OF NEGATIVE DATA VALUES.',//
     +,1X,'PLEASE CHOOSE A DIFFERENT TYPE OF PLOT.')
      CALL WRLINE(4)
      CALL PAUSE
      GOTO 2090
C
 2005 WRITE(*,2007)
 2007 FORMAT(////,1X,'TWO-PARAMETER LOGNORMAL ',
     +'PROBABILITY PLOT:',
     +13X,'Y',14X,'X',/)
C
      DO 2010 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=PN(ARG)
        LDATA(I)=DLOG(SDATA(I))
        WRITE(*,2008) Y(I),LDATA(I)
 2008   FORMAT(48X,F10.4,5X,F10.4)
 2010 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,2020) XLABEL
 2020 FORMAT(//,20X,'PLOT OF Y=NORMAL QUANTILE VS. ',
     +'X=LOG ',A8)
      CALL PLTTR(N,LDATA,Y)
C
 2090 RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE LN3PLT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),Y(100),ARG,PN,SHAPE
      CHARACTER*8 XLABEL
C
      WRITE(*,2100)
 2100 FORMAT(//,1X,'ENTER THE VALUE OF THE SHAPE '
     +'PARAMETER (SIGMA)',
     +/,1X,'FOR THE LOGNORMAL DISTRIBUTION:   '/)
      READ(*,*) SHAPE
C
 2105 WRITE(*,2107)
 2107 FORMAT(////,1X,'THREE-PARAMETER LOGNORMAL ',
     +'PROBABILITY PLOT:',
     +13X,'Y',14X,'X',/)
C
      DO 2110 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=DEXP(SHAPE*PN(ARG))
        WRITE(*,2108) Y(I),SDATA(I)
 2108   FORMAT(50X,F10.4,5X,F10.4)
 2110 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,2120) XLABEL,SHAPE
 2120 FORMAT(//,11X,'PLOT OF Y=LOGNORMAL QUANTILE VS.',
     +' X=',A8,5X,'** SIGMA = ',F6.2,' **')
      CALL PLTTR(N,SDATA,Y)
C
 2190 RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE HRPLOT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),Y(100),ARG,PN
      CHARACTER*8 XLABEL
C
      WRITE(*,1105)
 1105 FORMAT(////,1X,'RIGHT-TAIL HALF-NORMAL ',
     +'PROBABILITY PLOT:',12X,'Y',14X,'X',/)
C
      DO 1115 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=PN((ARG/2.0D0)+.50D0)
        WRITE(*,1110) Y(I),SDATA(I)
 1110   FORMAT(46X,F10.4,5X,F10.4)
 1115 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,1120) XLABEL
 1120 FORMAT(//,15X,'PLOT OF Y=RIGHT-TAIL NORMAL ',
     +'QUANTILE VS. X = ',A8)
      CALL PLTTR(N,SDATA,Y)
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE HLPLOT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),Y(100),ARG,PN
      CHARACTER*8 XLABEL
C
      WRITE(*,1205)
 1205 FORMAT(////,1X,'LEFT-TAIL HALF-NORMAL ',
     +'PROBABILITY PLOT:',12X,'Y',14X,'X',/)
C
      DO 1215 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=PN(ARG/2.0D0)
        WRITE(*,1210) Y(I),SDATA(I)
 1210   FORMAT(45X,F10.4,5X,F10.4)
 1215 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,1220) XLABEL
 1220 FORMAT(//,15X,'PLOT OF Y=LEFT-TAIL NORMAL ',
     +'QUANTILE VS. X = ',A8)
      CALL PLTTR(N,SDATA,Y)
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE EPLOT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),ARG,Y(100)
      CHARACTER*8 XLABEL
C
      WRITE(*,1300)
 1300 FORMAT(////,1X,'EXPONENTIAL PROBABILITY PLOT:',
     +17X,'Y',14X,'X',/)
C
      DO 1310 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=-DLOG(1.0D0-ARG)
        WRITE(*,1305) Y(I),SDATA(I)
 1305   FORMAT(40X,F10.4,5X,F10.4)
 1310 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,1320) XLABEL
 1320 FORMAT(//,20X,'PLOT OF Y=EXPONENTIAL ',
     +'QUANTILE VS. X=',A8)
      CALL PLTTR(N,SDATA,Y)
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE W2PLOT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),LDATA(100),ARG,Y(100)
      CHARACTER*8 XLABEL
C
      IF (SDATA(1).GT.0) GOTO 1402
      WRITE(*,1400)
 1400 FORMAT(//,1X,'A TWO-PARAMETER WEIBULL PLOT ',
     +'IS NOT POSSIBLE',/
     +,1X,'BECAUSE OF NEGATIVE DATA VALUES.',//
     +,1X,'PLEASE CHOOSE A DIFFERENT TYPE OF PLOT.')
      CALL WRLINE(4)
      CALL PAUSE
      GOTO 1490
C
 1402 WRITE(*,1403)
 1403 FORMAT(////,1X,'TWO-PARAMETER WEIBULL ',
     +'PROBABILITY PLOT:',
     +12X,'Y',14X,'X',/)
C
      DO 1410 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=DLOG(-DLOG(1.0D0-ARG))
        LDATA(I)=DLOG(SDATA(I))
        WRITE(*,1405) Y(I),LDATA(I)
 1405   FORMAT(45X,F10.4,5X,F10.4)
 1410 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,1420) XLABEL
 1420 FORMAT(//,13X,'PLOT OF Y = TYPE I EXTREME ',
     +'VALUE QUANTILE VS. X = LOG ',A8)
      CALL PLTTR(N,LDATA,Y)
C
 1490 RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE W3PLOT(N,SDATA,XLABEL)
C
      INTEGER N
      REAL*8 SDATA(100),C,Y(100),ARG
      CHARACTER*8 XLABEL
C
      WRITE(*,1500)
 1500 FORMAT(//,1X,'ENTER THE VALUE OF THE SHAPE '
     +'PARAMETER (C)',
     +/,1X,'FOR THE WEIBULL DISTRIBUTION:   '/)
      READ(*,*) C
C
      WRITE(*,1501)
 1501 FORMAT(////,1X,'THREE-PARAMETER WEIBULL ',
     +'PROBABILITY PLOT:',
     +12X,'Y',14X,'X',/)
C
      DO 1510 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=(-DLOG(1.0D0-ARG))**(1/C)
        WRITE(*,1505) Y(I),SDATA(I)
 1505   FORMAT(47X,F10.4,5X,F10.4)
 1510 CONTINUE
      CALL WRLINE(3)
      CALL PAUSE
C
      WRITE(*,1520) XLABEL,C
 1520 FORMAT(//,13X,'PLOT OF Y=WEIBULL QUANTILE ',
     +'VS. X=',A8,5X,'** C = ',F6.2,' **')
      CALL PLTTR(N,SDATA,Y)
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE SORT(N,OBS,SOBS)
C
C  Sort double precision numbers into ascending order
C
C  Source:  Loeser, Communications of the ACM,
C           Vol. 17, No. 3, page 143.
C
      REAL*8 OBS(N),SOBS(N),S
C
      DO 190 J=1,N
        SOBS(J)=OBS(J)
  190 CONTINUE
C
      I=1
  101 IF (I-N) 102,102,103
  102 I=I+I
      GO TO 101
  103 M=I-1
  104 M=M/2
      IF (M) 110,110,105
  105 K=N-M
      DO 109 J=1,K
        I=J+M
  106   I=I-M
        IF (I) 109,109,107
  107   L=I+M
        IF (SOBS(L)-SOBS(I)) 108,108,109
  108   S=SOBS(I)
        SOBS(I)=SOBS(L)
        SOBS(L)=S
        GO TO 106
  109 CONTINUE
      GO TO 104
  110 RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE ISORT(N,OBS,SOBS)
C
C  Sort integers into ascending order
C
C  Source:  Loeser, Communications of the ACM,
C           Vol. 17, No. 3, page 143.
C
      INTEGER OBS(N),SOBS(N),S
C
      DO 190 J=1,N
        SOBS(J)=OBS(J)
  190 CONTINUE
C
      I=1
  101 IF (I-N) 102,102,103
  102 I=I+I
      GO TO 101
  103 M=I-1
  104 M=M/2
      IF (M) 110,110,105
  105 K=N-M
      DO 109 J=1,K
        I=J+M
  106   I=I-M
        IF (I) 109,109,107
  107   L=I+M
        IF (SOBS(L)-SOBS(I)) 108,108,109
  108   S=SOBS(I)
        SOBS(I)=SOBS(L)
        SOBS(L)=S
        GO TO 106
  109 CONTINUE
      GO TO 104
  110 RETURN
      END
C
C-----------------------------------------------------
C
      DOUBLE PRECISION FUNCTION PN(P)
C
C   Algorithm AS 111, Applied Statistics (1977),
C   Vol. 26, pp 118-121.
C
      REAL*8 A0,A1,A2,A3,B1,B2,B3,B4,C0,C1,C2,C3,D1,
     +D2,P,Q,R
      A0=2.50662823884D0
      A1=-18.61500062529D0
      A2=41.39119773534D0
      A3=-25.44106049637D0
      B1=-8.47351093090D0
      B2=23.08336743743D0
      B3=-21.06224101826D0
      B4=3.13082909833D0
      C0=-2.78718931138D0
      C1=-2.29796479138D0
      C2=4.85014127135D0
      C3=2.32121276858D0
      D1=3.54388924762D0
      D2=1.63706781897D0
      Q=P-.50D0
      IF (DABS(Q) .GT. .42D0) GO TO 102
      R=Q*Q
      PN=Q*(((A3*R + A2)*R + A1)*R + A0)/
     +((((B4*R + B3)*R + B2)*R + B1)*R + 1.0D0)
      RETURN
C
  102 R=P
      IF (Q .GT. 0.0D0) R=1.0D0-P
      R=DSQRT(-DLOG(R))
      PN=(((C3*R + C2)*R +C1)*R + C0)/
     +((D2*R + D1)*R + 1.0D0)
      IF (Q .LT. 0.0D0) PN=-PN
      RETURN
C
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE PLTTR(N,HORIZ,VERT)
C
      INTEGER I,J,K,COUNT,PRNTFLG(75)
      REAL*8 VERT(128),HORIZ(128),UP,LO,MX,MN,X(128),
     +DX,DY,XMIN,XMAX
      CHARACTER*1 PRNTLN(68),PRNTCHR(6),GO
      DATA PRNTCHR/' ','*','*','*','*','*'/
C
C  Find largest and smallest horizontal coordinates
C  and determine the horizontal scaling factor, DX.
C
      MN=9999999999.9D0
      MX=-9999999999.9D0
      DO 900 I=1,N
       IF (HORIZ(I) .GT. MX) MX=HORIZ(I)
       IF (HORIZ(I) .LT. MN) MN=HORIZ(I)
  900 CONTINUE
      IF (MX.NE.MN) GOTO 902
      MX=MX+1.0D0
      MN=MN-1.0D0
  902 XMIN=MN
      XMAX=MX
      DX=(MX-MN)/64.0D0
C
C  Do the same for the vertical scale.
C
      MN=9999999999.9D0
      MX=-9999999999.9D0
      DO 905 I=1,N
       IF (VERT(I) .GT. MX) MX=VERT(I)
       IF (VERT(I) .LT. MN) MN=VERT(I)
  905 CONTINUE
      IF (MX.NE.MN) GOTO 908
      MX=MX+1.0D0
      MN=MN-1.0D0
  908 DY=(MX-MN)/15.0D0
C
C  Draw the plot a line at a time.
C
      WRITE (*,910)
  910 FORMAT(/,'   Y',5X,'|',/,9X,'|')
      LO=MX
      COUNT=3
C
      DO 940 I=1,17
        DO 915 J=1,65
         PRNTFLG(J)=1
  915   CONTINUE
        UP=LO
        LO=UP-DY
        DO 917 K=1,N
          IF ((VERT(K) .GT. LO) .AND.
     +        (VERT(K) .LE. UP)) THEN
            J=INT(SNGL((HORIZ(K)-XMIN)/DX))+1
            PRNTFLG(J)=PRNTFLG(J)+1
            IF (PRNTFLG(J) .GT. 6) PRNTFLG(J)=6
          END IF
  917   CONTINUE
      IF (COUNT.EQ.3) THEN
          WRITE (*,920) UP,
     +                  (PRNTCHR(PRNTFLG(J)),J=1,65)
  920     FORMAT (1X,F8.2,'+ ',65A1)
          COUNT=0
        ELSE
          WRITE (*,925) (PRNTCHR(PRNTFLG(J)),J=1,65)
  925     FORMAT (9X,'| ',65A1)
          COUNT=COUNT+1
        ENDIF
  940 CONTINUE
C
C  Put the labels under the horizontal axis.
C
      WRITE (*,945)
  945 FORMAT (10X,'-+-------+-------+-------+----',
     +'---+-------+-------+-------+-------+----')
      DX=(XMAX-XMIN)/8.0D0
      DO 950 I=1,9
        X(I)=XMIN+DBLE(I-1)*DX
  950 CONTINUE
C
      WRITE (*,960) (X(I),I=1,9)
  960 FORMAT(7X,9(F6.1,2X),'X')
C
      CALL PAUSE
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE WRLINE(LINES)
C
      INTEGER LINES
C
      DO 3010 I=1,LINES
        WRITE(*,3000)
 3000   FORMAT(' ')
 3010 CONTINUE
C
      RETURN
      END
C
C
C-----------------------------------------------------