C        
C                                  PC-SIZE
C                  A Program for Sample Size Determinations
C                                Version 2.0
C                              August 8, 1985
C
C                              Gerard E. Dallal
C        
C               USDA Human Nutrition Research Center on Aging
C                            at Tufts University
C                           711 Washington Street
C                             Boston, MA  02111
C
C                                    and
C        
C                    Tufts University School of Nutrition
C                             132 Curtis Street
C                             Medford, MA  02155
C        
C
C
C                                  NOTICE
C
C       Documentation and original code copyright 1985 by  Gerard  E.  
C       Dallal.  Reproduction of material for non-commercial purposes 
C       is   permitted,   without  charge,   provided  that  suitable 
C       reference is made to PC-SIZE and its author.  
C
C       Neither PC-SIZE nor its documentation should be  modified  in 
C       any way without permission from the author,  except for those 
C       changes  that  are  essential  to  move  PC-SIZE  to  another 
C       computer.  
C
C       Please acknowledge PC-SIZE in any manuscript  that  uses  its 
C       calculations.  
C
C
C        IIN   -- INPUT UNIT NUMBER (SCREEN)
C        IOUT  -- OUTPUT UNIT NUMBER (SCREEN)
C        IWOUT -- SAVED OUTPUT UNIT NUMBER 
C        NMAX0 -- LARGEST PERMISSIBLE INTEGER
C        NEMAX -- LARGEST NUMBER OF EFFECTS THAT CAN BE ENTERED
C                    INDIVIDUALLY  (DIMENSION OF 'EFFECT')
C
C        OPEN STATEMENT LOCATED JUST BEFORE STATEMENT 10
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL QFILE, QAB
      CHARACTER*1 QUERY, BLANK
      CHARACTER FNAME*40, CF2C0*18, CEFECT*18
      DIMENSION EFFECT(50)
C
      DATA IIN /0/, IOUT /0/, IWOUT /11/, NMAX0 /2147483647/
      DATA QFILE /.FALSE./, BLANK /' '/, NEMAX /50/
C
C        INITIALIZE CHOICES
C
      MODE = 1
      ALPHA = 0.05
      RPOWER = 0.80
      NG = 2
      NROW = 2
      NCOL = 2
      EDIFF = 0.0
      WSD = 0.0
      F1 = 0.0
      F2C1 = 0.0
      F2C0 = 0.0
      XLAM0 = 0.0
      XLAM01 = 0.0
      IMETH = 3
      RANGE = 0.0
      RHO = 0.0
      P1 = 0.0
      P2 = 0.0
C
      WRITE(IOUT,1)
    1 FORMAT (//36X,'PC-SIZE'/
     *   20X,'A Program for Sample Size Determinations'/
     *   34X,'Version 2.0'/32X,'August 8, 1985'//32X,'Gerard E. Dallal'/
     *   17X,'USDA Human Nutrition Research Center on Aging'/
     *   30X,'at Tufts University'/29X,'711 Washington Street'/
     *   31X,'Boston, MA  02111'///37X,'NOTICE'//9X,'Please '
     *   'acknowledge  PC-SIZE  in any manuscript  that uses its',/
     *   9X,'calculations.')

      WRITE(IOUT,2)
    2 FORMAT(//' Do you wish to save the results of this session?  ',
     *  '[N]:  '$)
      READ (IIN,3) QUERY
    3 FORMAT (A1)
      IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QFILE = .TRUE.
      IF (.NOT. QFILE) GOTO 10
      WRITE (IOUT,4)
    4 FORMAT (' Enter output filename:  '$)
      READ (IIN,'(A)') FNAME
      OPEN (IWOUT, FILE = FNAME, STATUS = 'NEW')
      WRITE (IWOUT,1)

   10 WRITE(IOUT,20) 
      IF (QFILE) WRITE (IWOUT,20) 
   20 FORMAT(///' This program operates in seven modes:'//
     *   '     1  --  single factor design'/
     *   '     2  --  two factor design'/
     *   '     3  --  randomized blocks design'/
     *   '     4  --  paired t-test'/
     *   '     5  --  generic F'/
     *   '     6  --  correlation coefficient'/
     *   '     7  --  proportions'/)
      IF (QFILE) WRITE (IWOUT,22) MODE
   21 WRITE(IOUT,22) MODE
   22 FORMAT(' Choose a mode [',I1,']:  '$)
      READ (IIN,50,ERR=21) XMODE
      IF (XMODE.NE.DINT(XMODE) .OR. XMODE.LT.0.0D0
     *   .OR. XMODE.GT.7.0D0) GOTO 21
      IF (XMODE.NE.0.0D0) MODE = XMODE
      IF (QFILE) WRITE(IWOUT,29) MODE
   29 FORMAT (1X,I6)
      IF (QFILE) WRITE (IWOUT,40) ALPHA
   30 WRITE (IOUT,40) ALPHA
   40 FORMAT (/' Enter the level of the test [',F5.3,']:  '$)
      READ (IIN,50,ERR=30) DALPHA
      IF (DALPHA.LT.0.0D0 .OR. DALPHA.GE.1.0D0) GOTO 30
      IF (DALPHA.NE.0.0D0) ALPHA = DALPHA
      IF (QFILE) WRITE (IWOUT,49) ALPHA
   49 FORMAT (1X,G12.5)
   50 FORMAT (BN,F18.0)

      IF (RPOWER.LT.1.0D0 .AND. QFILE) WRITE (IWOUT,61) RPOWER
      IF (RPOWER.GE.1.0D0 .AND. QFILE) WRITE (IWOUT,62) RPOWER
   60 IF (RPOWER.LT.1.0D0) WRITE (IOUT,61) RPOWER
   61 FORMAT (' Enter the required power [',F5.3,']:  '$)
      IF (RPOWER.GE.1.0D0) WRITE (IOUT,62) RPOWER
   62 FORMAT (' Enter the required power [',F8.0,']:  '$)
      READ (IIN,50,ERR=60) DPOWER
      IF (DPOWER.LT.0.0D0 .OR. 
     *   (DPOWER.GT.1.D0 .AND. DPOWER.NE.DINT(DPOWER))) GOTO 60
      IF (MODE.EQ.6 .AND. DPOWER.GE.1.0D0 .AND. DPOWER.LT.3.0D0)
     *   GOTO 60
      IF (DPOWER.NE.0.0D0) RPOWER = DPOWER
      IF (QFILE) WRITE (IWOUT,49) RPOWER
      IF (RPOWER.LT.1.0D0) GOTO 79

      NSTART = RPOWER
      IF (QFILE) WRITE (IWOUT,69) NSTART
      WRITE (IOUT,69) NSTART
   69 FORMAT (' Power calculations will start with a sample size of:',
     *   I6)
      IF (QFILE) WRITE (IWOUT,71) 
   70 WRITE (IOUT,71) 
   71 FORMAT(' Enter the largest sample size to be considered:  '$)
      READ (IIN,50,ERR=70) XNMAX
      IF (XNMAX.NE.DINT(XNMAX) .OR. XNMAX.LT.RPOWER) GOTO 70
      NMAX = XNMAX
      IF (QFILE) WRITE (IWOUT,29) NMAX
      NSTEP = 1
      IF (NMAX.EQ.NSTART .OR. NMAX.EQ.NSTART+1) GOTO 79
C
      IF (QFILE) WRITE (IWOUT,73) 
   72 WRITE (IOUT,73) 
   73 FORMAT (' Enter increment:  '$)
      READ (IIN,50,ERR=72) XNSTEP
      IF(XNSTEP.LT.1.0D0 .OR. XNSTEP.NE.DINT(XNSTEP)) GOTO 72
      NSTEP = XNSTEP
      IF (QFILE) WRITE (IWOUT,29) NSTEP
C
C        GET NUMERATOR DEGREES OF FREEDOM
C
   79 GOTO (80, 110, 210, 230, 300, 2000, 3000), MODE
C
C        SINGLE FACTOR DESIGN
C
   80 IF (QFILE) WRITE (IWOUT,90) NG
   81 WRITE (IOUT,90) NG
   90 FORMAT (' Enter the number of groups [',I2,']:  '$)
      READ (IIN,50,ERR=81) XNG
      IF (XNG.LT.0.0D0 .OR. XNG.NE.DINT(XNG)) GOTO 81
      IF (XNG.NE.0.0D0) NG = XNG
      XNG = NG
      IF (NG.LT.2) GOTO 81
      IF (QFILE) WRITE (IWOUT,29) NG
      F1 = NG - 1
      GOTO 400
C
C        TWO FACTOR DESIGN
C
  110 IF (QFILE) WRITE (IWOUT,120) NROW
  111 WRITE (IOUT,120) NROW
  120 FORMAT (/' Enter the number of levels for Factor A [',I2,']:  '$)
      READ (IIN,50,ERR=111) XNROW
      IF (XNROW.LT.0.0D0 .OR. XNROW.NE.DINT(XNROW)) GOTO 111
      IF (XNROW.NE.0.0D0) NROW = XNROW
      IF (NROW.LT.2) GOTO 111
      IF (QFILE) WRITE (IWOUT,29) NROW
      XNROW = NROW

      IF (QFILE) WRITE (IWOUT,140) NCOL
  130 WRITE (IOUT,140) NCOL
  140 FORMAT (' Enter the number of levels for Factor B [',I2,']:  '$)
      READ (IIN,50,ERR=130) XNCOL
      IF (XNCOL.LT.0.0D0 .OR. XNCOL.NE.DINT(XNCOL)) GOTO 130
      IF (XNCOL.NE.0.0D0) NCOL = XNCOL
      IF (NCOL.LT.2) GOTO 130
      IF (QFILE) WRITE (IWOUT,29) NCOL
      XNCOL = NCOL
      NG = NROW
      XNG = NG
      F1 = NG - 1

      IF (QFILE) WRITE (IWOUT,160)
  150 WRITE (IOUT,160)
  160 FORMAT (' Do interactions appear in the model and the ANOVA ',
     *   'table?  (Y or N):  '$)
      READ (IIN,3) QUERY
      IF (QUERY.NE.'Y' .AND. QUERY.NE.'y' .AND. QUERY.NE.'N'
     *   .AND. QUERY.NE.'n') GOTO 150
      QAB = .FALSE.
      IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QAB = .TRUE.
      IF (QFILE) WRITE (IWOUT,161) QUERY
  161 FORMAT(1X,A1)
      GOTO 400
C
C        RANDOMIZED BLOCKS DESIGN
C      
  210 IF (QFILE) WRITE (IWOUT,220) NG
  211 WRITE (IOUT,220) NG
  220 FORMAT(' Enter the number of levels for the treatment factor [',
     *   I2,']:  '$)
      READ (IIN,50,ERR=211) XNG
      IF (XNG.LT.0.0D0 .OR. XNG.NE.DINT(XNG)) GOTO 211
      IF (XNG.NE.0.0) NG = XNG
      XNG = NG
      IF (NG.LT.2) GOTO 211
      IF (QFILE) WRITE (IWOUT,29) NG
      F1 = NG - 1
      GOTO 400
C
C        PAIRED T-TEST
C
  230 IF (EDIFF.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,240) EDIFF
  240 FORMAT (' Enter expected difference [',G12.5,']:  '$)
      IF (EDIFF.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,241) 
  241 FORMAT (' Enter expected difference [0.00]:  '$)
  242 IF (EDIFF.NE.0.0D0) WRITE (IOUT,240) EDIFF
      IF (EDIFF.EQ.0.0D0) WRITE (IOUT,241) 
      READ (IIN,50,ERR=242) XDIFF
      IF (XDIFF.NE.0.0D0) EDIFF = XDIFF
      IF (EDIFF.EQ.0.0D0) GOTO 242
      IF (QFILE) WRITE (IWOUT,49) EDIFF
      F1 = 1
      NG = 2
      XNG = 2.0

      IF (WSD.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,250) WSD
  250 FORMAT (' Enter estimate of standard deviation of difference [',
     *   G12.5,']:  '$)
      IF (WSD.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,255) 
  255 FORMAT (' Enter estimate of standard deviation of difference [',
     *   '0.00]:  '$)
  256 IF (WSD.NE.0.0D0) WRITE (IOUT,250) WSD
      IF (WSD.EQ.0.0D0) WRITE (IOUT,255) 
      READ (IIN,50,ERR=256)  XWSD
      IF (XWSD.LT.0.0D0) GOTO 256
      IF (XWSD.GT.0.0D0) WSD = XWSD
      IF (QFILE) WRITE (IWOUT,49) WSD
      IF (WSD.EQ.0.0D0) GOTO 256
      XLAM0 =  (EDIFF / WSD) ** 2
      GOTO 900
C
C        GENERIC
C

  300 NF1 = F1
      IF (QFILE) WRITE (IWOUT,301) NF1
  301 FORMAT (' Enter the numerator degrees of freedom [',I2,']:  '$)
  310 NF1 = F1
  311 WRITE (IOUT,301) NF1
      READ (IIN,50,ERR=311) XF1
      IF (XF1.LT.0.0D0 .OR. XF1.NE.DINT(XF1)) GOTO 311
      IF (XF1.GT.0.0D0) F1 = XF1
      IF (F1.LE.0.0D0) GOTO 310
      IF (QFILE) WRITE (IWOUT,49) F1

      WRITE (IOUT,313) 
      IF (QFILE) WRITE (IWOUT,313) 
  313 FORMAT (' The denominator degrees of freedom are calculated as'/
     *   '   the linear function  C1 * N - C0  of the sample size, N.')
      NF2C1 = F2C1
      IF (QFILE) WRITE (IWOUT,317) NF2C1
  314 NF2C1 = F2C1
  316 WRITE (IOUT,317) NF2C1
  317 FORMAT(' Enter C1 [',I3,']:  '$)
      READ (IIN,50,ERR=316) XF2C1
      IF (XF2C1.LT.0.0D0 .OR. XF2C1.NE.DINT(XF2C1)) GOTO 316
      IF (XF2C1.NE.0.0D0) F2C1 = XF2C1
      IF (F2C1.LE.0.0D0) GOTO 314
      IF (QFILE) WRITE (IWOUT,49) F2C1

      NF2C0 = F2C0
      IF (QFILE) WRITE (IWOUT,322) NF2C0
  320 NF2C0 = F2C0
  321 WRITE (IOUT,322) NF2C0
  322 FORMAT (' Enter C0 [',I3,']:  '$)
      READ (IIN,'(A)') CF2C0

C
C        CHECK FOR EMPTY
C
      DO 326 I = 1,18
      IF (CF2C0(I:I).NE.BLANK) GOTO 327
  326 CONTINUE
      GOTO 330

  327 READ (CF2C0,50,ERR=321) XF2C0
      IF (XF2C0.NE.DINT(XF2C0)) GOTO 320
      F2C0 = XF2C0
  330 IF (QFILE) WRITE (IWOUT,49) F2C0

      WRITE (IOUT,340) 
      IF (QFILE) WRITE (IWOUT,340) 
  340 FORMAT (' The noncentrality parameter '
     *   'is calculated as  C * N.')
      IF (XLAM0.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,350) 
  350 FORMAT (' Enter  C [0.00]:  '$) 
      IF (XLAM0.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,360) XLAM0
  360 FORMAT (' Enter  C [',G12.5,']:  '$) 
  365 IF (XLAM0.EQ.0.0D0) WRITE (IOUT,350) 
      IF (XLAM0.NE.0.0D0) WRITE (IOUT,360) XLAM0
      READ (IIN,50,ERR=365) XXLAM0
      IF (XXLAM0.LT.0.0D0) GOTO 365
      IF (XXLAM0.NE.0.0D0) XLAM0 = XXLAM0
      IF (XLAM0.LE.0.0D0) GOTO 365
      IF (QFILE) WRITE (IWOUT,49) XLAM0

      GOTO 900

  400 WRITE (IOUT,410) 
      IF (QFILE) WRITE (IWOUT,410) 
  410 FORMAT (/' There are three ways to specify the alternative:'//
     *   '    1 -- Specify individual effects.'/
     *   '    2 -- Specify a range throughout which group'/
     *   '            means are uniformly distrbuted.'/
     *   '    3 -- Specify average squared effect divided'/
     *   '            by the error variance.'/)
C    
      IF (QFILE) WRITE (IWOUT,420) IMETH
  415 WRITE (IOUT,420) IMETH
  420 FORMAT (' Choose a method [',I1,']:  '$)  
C
      READ (IIN,50,ERR=415) XIMETH
      IF (XIMETH.LT.0.0D0 .OR. 
     *   XIMETH.GT.3.0D0 .OR. XIMETH.NE.DINT(XIMETH)) GOTO 415
      IF (XIMETH.NE.0.0D0) IMETH = XIMETH
      IF (QFILE) WRITE (IWOUT,29) IMETH

      GOTO (500,600,800),IMETH

  500 EAVG = 0.
      IF (NG.LE.NEMAX) GOTO 502
      WRITE(IOUT,501) NEMAX
  501 FORMAT (/' NO MORE THAN',I3,' EFFECTS MAY BE SPECIFIED')
      GOTO 400

  502 DO 520  I = 1, NG
      IF (MODE.EQ.1 .AND. QFILE) WRITE (IWOUT,510) I,EFFECT(I)
  510 FORMAT (' Enter estimate of effect for group ',I2,' [',
     *   G12.5,']:  '$)
      IF (MODE.EQ.2 .AND. QFILE) WRITE (IWOUT,511) I,EFFECT(I)
  511 FORMAT (' Enter estimate of effect for Factor A, level',I2,' [',
     *   G12.5,']:  '$)
      IF (MODE.EQ.3 .AND. QFILE) WRITE (IWOUT,512) I,EFFECT(I)
  512 FORMAT (' Enter estimate of effect for treatment ',I2,' [',
     *   G12.5,']:  '$)
  513 IF (MODE.EQ.1) WRITE (IOUT,510) I,EFFECT(I)
      IF (MODE.EQ.2) WRITE (IOUT,511) I,EFFECT(I)
      IF (MODE.EQ.3) WRITE (IOUT,512) I,EFFECT(I)

      READ (IIN,'(A)')  CEFECT
C
C        CHECK FOR EMPTY
C
      DO 514 II = 1, 18
      IF (CEFECT(II:II).NE.BLANK) GOTO 515
  514 CONTINUE
      GOTO 519

  515 READ (CEFECT,50,ERR=513) EFFECT(I)

  519 IF (QFILE) WRITE (IWOUT,49) EFFECT(I)
      EAVG = EAVG + EFFECT(I)
  520 CONTINUE
      EAVG = EAVG / XNG
      XLAM0 = 0.
      DO 530 I = 1, NG
  530 XLAM0 = XLAM0 + (EFFECT(I) - EAVG) ** 2
      IF (XLAM0.GT.0.0D0) GOTO 700
      WRITE (IOUT,540)
      IF (QFILE) WRITE (IWOUT,540)
  540 FORMAT (/' EFFECTS ARE ALL EQUAL.  PLEASE RESPECIFY.'/)
      GOTO 500

  600 IF (RANGE.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,610) 
  610 FORMAT (' Enter range [0.00]:  '$)
      IF (RANGE.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,620) RANGE
  620 FORMAT (' Enter range [',G12.5,']:  '$)
  630 IF (RANGE.EQ.0.0D0) WRITE (IOUT,610) 
      IF (RANGE.NE.0.0D0) WRITE (IOUT,620) RANGE
      READ (IIN,50,ERR=630) XRANGE
      IF (XRANGE.LT.0.0) GOTO 630
      IF (XRANGE.GT.0.0D0) RANGE = XRANGE
      IF (RANGE.LE.0.0D0) GOTO 630
      IF (QFILE) WRITE (IWOUT,49) RANGE
      DELTA = RANGE / (XNG - 1.0D0)
      XLAM0 = DELTA ** 2 * (XNG ** 3 - XNG) / 12.0D0 

  700 IF (WSD.NE.0.0D0 .AND. QFILE) WRITE (IWOUT,720) WSD
  720 FORMAT (' Enter estimate of within cell standard deviation [',
     *   G12.5,']:  '$)
      IF (WSD.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,730) 
  730 FORMAT (' Enter estimate of within cell standard deviation [',
     *   '0.00]:  '$)
  740 IF (WSD.NE.0.0D0) WRITE (IOUT,720) WSD
      IF (WSD.EQ.0.0D0) WRITE (IOUT,730) 
      READ (IIN,50,ERR=740) XWSD
      IF (XWSD.LT.0.0D0) GOTO 740
      IF (XWSD.GT.0.0D0) WSD = XWSD
      IF (WSD.LE.0.0D0) GOTO 740
      IF (QFILE) WRITE (IWOUT,49) WSD
      XLAM0 =  XLAM0 / WSD ** 2
      GOTO 900
     
  800 IF (XLAM01.GT.0.0D0 .AND. QFILE) WRITE (IWOUT,810) XLAM01
  810 FORMAT (' Specify average squared effect divided'/
     *   '    by the error variance [',G12.5,']:  '$)
      IF (XLAM01.EQ.0.0D0 .AND. QFILE) WRITE (IWOUT,820) 
  820 FORMAT (' Specify average squared effect divided'/
     *   '    by the error variance [0.00]:  '$)
  830 IF (XLAM01.GT.0.0D0) WRITE (IOUT,810) XLAM01
      IF (XLAM01.EQ.0.0D0) WRITE (IOUT,820) 
      READ (IIN,50,ERR=830) XXLAM0
      IF (XXLAM0.LT.0.0D0) GOTO 830
      IF (XXLAM0.GT.0.0D0) XLAM01 = XXLAM0
      IF (XLAM01.EQ.0.0D0) GOTO 830
      IF (QFILE) WRITE (IWOUT,49) XLAM01
      XLAM0 = XLAM01 * XNG
C
  900 IF (MODE.EQ.2) XLAM0 =  XNCOL * XLAM0
      WRITE (IOUT,910) ALPHA,XLAM0
      IF (QFILE) WRITE (IWOUT,910) ALPHA,XLAM0
  910 FORMAT(//' Level of the test:  ',F6.4/
     *   ' Non-centrality parameter: ',G12.5,' *  sample size')
C
C        LARGE SAMPLE APPROXIMATION
C
      IF (RPOWER .GE. 1.0D0) GOTO 979
      ALPHAC = 1.0D0 - ALPHA
      CHISQ = PPCHI2(ALPHAC, F1)
      CALL CHISQL(CHISQ, F1, 1.0D0 - RPOWER, FLLRGE)
      XNLRGE = FLLRGE / XLAM0 
      NLRGE = XNLRGE
      NLRGE1 = 0.05D0 * XNLRGE
C
      NSTART = MAX0(NLRGE - 10, 1)
      IF (NLRGE1.GT.10) NSTART = NLRGE - NLRGE1
      NMAX = NMAX0
      NSTEP = DLOG10(XNLRGE) - DLOG10(50.D0)
      IF (NSTEP.LT.0) NSTEP = 0
      NSTEP = 10 ** NSTEP
      NSTART = NSTART / NSTEP * NSTEP
C
  979 WRITE (IOUT,980)
      IF (QFILE) WRITE (IWOUT,980)
  980 FORMAT (//'     SAMPLE'/'      SIZE     POWER'/)

  989 DO 1000 N = NSTART, NMAX, NSTEP
      XN = N
      XLAM = XN * XLAM0
      IF (MODE.EQ.1) F2 = XNG * (XN - 1.0D0)
      IF (MODE.EQ.2 .AND. QAB) F2 = XNROW * XNCOL * (XN - 1.0D0)
      IF (MODE.EQ.2 .AND. .NOT.QAB)
     *    F2 = XNROW * XNCOL * XN - XNROW - XNCOL + 1.0D0
      IF (MODE.EQ.3 .OR. MODE.EQ.4) F2 = (XNG - 1.0D0) * (XN - 1.0D0)
      IF (MODE.EQ.5) F2 = F2C1 * XN - F2C0
      IF (F2 .LT. 1.0D0) GOTO 1000
      AA = F2 / 2.0D0
      BB = F1 / 2.0D0
      XALPHA = XINBTA(AA, BB, ALPHA)
      FALPHA = (F2 / F1) *((1.0D0 - XALPHA) / XALPHA)
      POWER = FNCDFC(FALPHA, F1, F2, XLAM)
      WRITE (IOUT,990) N, POWER
      IF (QFILE) WRITE (IWOUT,990) N, POWER
  990 FORMAT (1X, I8, 3X, F8.5)
      IF (POWER .GT. RPOWER) GOTO 1010
 1000 CONTINUE
 1010 WRITE (IOUT,1020)
 1020 FORMAT (///' Do you wish to continue?  [Y]:  '$)
      READ (IIN,3) QUERY
      IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
      GOTO 10
C
C        CORRELATION
C
 2000 IF (QFILE) WRITE (IWOUT,2080) RHO
 2080 FORMAT (/' Enter the correlation coefficient'/
     *   '    at the alternative [',F5.3,']:  '$)
 2081 WRITE (IOUT,2080) RHO
      READ (IIN,50,ERR=2081) DRHO
      IF (DABS(DRHO).GE.1.0D0) GOTO 2081
      IF (DRHO.NE.0.0D0) RHO = DRHO
      IF (RHO.EQ.0.0D0) GOTO 2081
      IF (QFILE) WRITE (IWOUT,49) RHO
      IF (RPOWER.LT.1.0D0) GOTO 2200
C
      WRITE (IOUT,980)
      IF (QFILE) WRITE (IWOUT,980)

      DO 2100 N = NSTART, NMAX, NSTEP
      XN = N
      R = RCRIT(ALPHA,XN - 2.0D0)
      POWER = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
     *   RDF(N, R, RHO, IFAULT) 
      WRITE (IOUT,990) N, POWER
      IF (QFILE) WRITE (IWOUT,990) N, POWER
 2100 CONTINUE
      WRITE (IOUT,1020)
      READ (IIN,3) QUERY
      IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
      GOTO 10


C
C        N = 3
C
 2200 WRITE (IOUT,980)
      IF (QFILE) WRITE (IWOUT,980)
      N = 3
      XN = N
      R = RCRIT(ALPHA,XN - 2.0D0)
      P3 = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
     *   RDF(N, R, RHO, IFAULT)
      WRITE (IOUT,990) N, P3
      IF (QFILE) WRITE (IWOUT,990) N, P3
      IF (P3.GT.RPOWER) GOTO 230
C
C        N = 4
C
      N = 4
      XN = N
      R = RCRIT(ALPHA,XN - 2.0D0)
      P4 = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
     *   RDF(N, R, RHO, IFAULT) 
C
      WRITE (IOUT,990) N, P4
      IF (QFILE) WRITE (IWOUT,990) N, P4
      IF (P4.GT.RPOWER) GOTO 2230
C
 2210 N = 2 * N
      XN = N
      R = RCRIT(ALPHA,XN - 2.0D0)
      POWER = RDF(N, -R, RHO, IFAULT) + 1.0D0 -
     *   RDF(N, R, RHO, IFAULT) 
      WRITE (IOUT,990) N, POWER
      IF (QFILE) WRITE (IWOUT,990) N, POWER
      IF (POWER.LT.RPOWER) GOTO 2210
C
      NHI = N
      HIPOWR = POWER
      NLO = N / 2
 2220 NM = (NLO + NHI) / 2
      XN = NM
      R = RCRIT(ALPHA,XN - 2.0D0)
      POWER = RDF(NM, -R, RHO, IFAULT) + 1.0D0 -
     *   RDF(NM, R, RHO, IFAULT) 
      IF (POWER.LT.RPOWER) NLO = NM
      IF (POWER.GE.RPOWER) NHI = NM
      IF (POWER.GE.RPOWER) HIPOWR = POWER
      IF (NHI - NLO.GT.1) GOTO 2220
      WRITE (IOUT,990) NHI, HIPOWR
      IF (QFILE) WRITE (IWOUT,990) NHI, HIPOWR

 2230 WRITE (IOUT,1020)
      READ (IIN,3) QUERY
      IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
      GOTO 10
C
C        PROPORTIONS -- TWO BY TWO TABLES
C
 3000 IF (QFILE) WRITE (IWOUT,3080) P1
 3080 FORMAT (/' Enter the first proportion under'/
     *   '    the alternative to equality [',F5.3,']:  '$)
 3081 WRITE (IOUT,3080) P1
      READ (IIN,50,ERR=3081) DP1
      IF (DP1.LT.0.0D0 .OR. DP1.GE.1.0D0) GOTO 3081
      IF (DP1.NE.0.0D0) P1 = DP1
      IF (P1.EQ.0.0D0) GOTO 3081
      IF (QFILE) WRITE (IWOUT,49) P1
C
      IF (QFILE) WRITE (IWOUT,3090) P2
 3090 FORMAT (/' Enter the second proportion under'/
     *   '    the alternative to equality [',F5.3,']:  '$)
 3091 WRITE (IOUT,3090) P2
      READ (IIN,50,ERR=3091) DP2
      IF (DP2.LT.0.0D0 .OR. DP2.GE.1.0D0) GOTO 3091
      IF (DP2.NE.0.0D0) P2 = DP2
      IF (P2.EQ.0.0D0) GOTO 3091
      IF (QFILE) WRITE (IWOUT,49) P2
C
      IF (P1.NE.P2) GOTO 3095
      WRITE(IOUT,3094)
 3094 FORMAT(/' The two proportions are equal.  Please respecify.'/)
      GOTO 3000
C
 3095 PBAR = (P1 + P2) / 2.0D0
      DELTA = DABS(P2 - P1)
C
      IF (RPOWER.LT.1.0D0) GOTO 3200
C
      WRITE (IOUT,980)
      IF (QFILE) WRITE (IWOUT,980)

      DO 3100 N = NSTART, NMAX, NSTEP
      XN = N
      XM = XN - 2.0D0 / DELTA + 1.0D0 / (XN * DELTA**2)
      PC1 = GAUINV(1.0D0 - ALPHA / 2.0D0, IFAULT)
      PC2 = (PC1 * DSQRT(2.0D0 * PBAR * (1.0D0 - PBAR)) - DELTA *
     *   DSQRT(XM)) / DSQRT(P1 * (1.0D0 - P1) + P2 * (1.0D0 - P2))
      POWER = ALNORM(PC2,.TRUE.)
      WRITE (IOUT,990) N, POWER
      IF (QFILE) WRITE (IWOUT,990) N, POWER
 3100 CONTINUE
      WRITE (IOUT,1020)
      READ (IIN,3) QUERY
      IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
      GOTO 10
C
 3200 PC1 = GAUINV(1.0D0 - ALPHA / 2.0D0, IFAULT)
      PC2 = GAUINV(1.0D0 - RPOWER,IFAULT) 
      XN = ((PC1 * DSQRT(2.0D0 * 
     *   PBAR * (1.0D0 - PBAR)) - PC2 * DSQRT(P1 *
     *   (1.0D0 - P1) + P2 * (1.0D0 - P2))) / DELTA)**2
      XN = XN / 4.0D0 * (1.0D0 + DSQRT(1.0D0 + 4.0D0 /
     *   (XN * DELTA)))**2 + 0.9999999999D0
C
      N = XN
      XN = N
      XM = XN - 2.0D0 / DELTA + 1.0D0 / (XN * DELTA**2)
      PC1 = GAUINV(1.0D0 - ALPHA / 2.0D0, IFAULT)
      PC2 = (PC1 * DSQRT(2.0D0 * PBAR * (1.0D0 - PBAR)) - DELTA *
     *   DSQRT(XM)) / DSQRT(P1 * (1.0D0 - P1) + P2 * (1.0D0 - P2))
      POWER = ALNORM(PC2,.TRUE.)

      WRITE (IOUT,980)
      IF (QFILE) WRITE (IWOUT,980)
      WRITE (IOUT,990) N, POWER
      IF (QFILE) WRITE (IWOUT,990) N, POWER

      WRITE (IOUT,1020)
      READ (IIN,3) QUERY
      IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
      GOTO 10

      END
C
      DOUBLE PRECISION FUNCTION ALGAMA(S)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        CACM ALGORITHM 291
C        LOGARITHM OF GAMMA FUNCTION
C        BY PIKE, M.C. AND HILL, I.D.
C
C        TRANSLATED INTO FORTRAN BY
C                   Gerard E. Dallal
C                   USDA-Human Nutrition Research Center on Aging
C                              at Tufts University
C                   711 Washington Street
C                   Boston, MA  02111
C
C
      X = S
      ALGAMA = 0.0
      IF (X .LE. 0.0D0) RETURN
      F = 0.0
      IF (X .GE. 7.0D0) GOTO 30
C
      F = 1.0
      Z = X - 1.0D0
C
   10 Z = Z + 1.0D0
      IF (Z .GE. 7.0D0) GOTO 20
      X = Z
      F = F * Z
      GOTO 10
C
   20 X = X + 1.0D0
      F = -DLOG(F)
C
   30 Z = 1.0D0 / X ** 2
      ALGAMA = F + ( X - 0.5D0) * DLOG(X) - X + 0.918938533204673D0 +
     1   (((-0.000595238095238D0 * Z + 0.000793650793651D0) * Z
     2   - 0.002777777777778D0) * Z + 0.083333333333333D0) / X
      RETURN
      END

      DOUBLE PRECISION FUNCTION ALNORM(X,UPPER)
C
C        ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3
C
C        EVALUATES THE TAIL AREA OF THE STANDARD NORMAL CURVE
C        FROM X TO INFINITY IF UPPER IS .TRUE. OR
C        FROM MINUS INFINITY TO X IF UPPER IS .FALSE.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LTONE
      LOGICAL UPPER, UP
C
C        LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR
C        COMPUTER (SEE INTRODUCTORY TEXT)
C
      DATA LTONE, UTZERO/7.0D0, 18.66D0/
      DATA ZERO, HALF, ONE, CON/0.0D0, 0.5D0, 1.0D0, 1.28D0/
      UP = UPPER
      Z = X
      IF (Z .GE. ZERO) GOTO 10
      UP = .NOT. UP
      Z = -Z
   10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20
      ALNORM = ZERO
      GOTO 40
   20 Y = HALF * Z * Z
      IF (Z .GT. CON) GOTO 30
C
      ALNORM = HALF - Z * (0.398942280444D0 - 0.399903438504D0 * Y /
     1   (Y + 5.75885480458D0 - 29.8213557808D0 /
     2   (Y + 2.62433121679D0 + 48.6959930692D0 /
     3   (Y + 5.92885724438D0))))
      GOTO 40
C
   30 ALNORM = .398942280385D0 * EXP(-Y) /
     1   (Z - 3.8052D-8 + 1.00000615302D0 /
     2   (Z + 3.98064794D-4 + 1.98615381364D0 /
     3   (Z - 0.151679116635D0 + 5.29330324926D0 /
     4   (Z + 4.8385912808D0 - 15.1508972451D0 /
     5   (Z + 0.742380924027D0 + 30.789933034D0 / 
     6   (Z + 3.99019417011D0))))))
C
   40 IF (.NOT. UP) ALNORM = ONE - ALNORM
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION BETAIN(X,P,Q)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        ALGORITHM AS 63  APPL.STATIST. (1973), VOL.22, NO.3
C        MODIFIED AS PER REMARK ASR 19  (1977), VOL.26, NO.1
C        ***[FAULT INDICATOR REMOVED BY G.E.D.]***
C
C        COMPUTES INCOMPLETE BETA FUNCTION RATIO FOR ARGUMENTS
C        X BETWEEN ZERO AND ONE, P AND Q POSITIVE.
C        LOG OF COMPLETE BETA FUNCTION, BETA, ASSUMED TO BE KNOWN.
C        ***[CALCULATION OF BETA ADDED BY G.E.D]***
C
      LOGICAL INDEX
C
C        DEFINE ACCURACY AND INITIALIZE
C
      DATA ACU /0.1D-7/
      BETAIN = X 
C
C        TEST FOR ADMISIBILITY OF AGRUMENTS
C
      IF (P .LE. 0.0D0 .OR. Q .LE. 0.0D0) STOP 40
      IF (X .LT. 0.0D0 .OR .X .GT. 1.0D0) STOP 41
      IF (X .EQ .0.0D0 .OR .X .EQ. 1.0D0) RETURN
C
C     CHANGE TAIL IF NECESSARY AND DETERMINE S
C
      PSQ = P + Q
      CX = 1.0D0 - X
      IF (P .GE .PSQ * X) GOTO 1
      XX = CX
      CX = X
      PP = Q
      QQ = P
      INDEX = .TRUE.
      GOTO 2
    1 XX = X
      PP = P
      QQ = Q
      INDEX = .FALSE.
    2 TERM = 1.0
      AI = 1.0
      BETAIN = 1.0
      NS = QQ + CX * PSQ
C
C        USE SOPER'S REDUCTION FORMULAE
C
      RX = XX / CX
    3 TEMP = QQ - AI
      IF (NS .EQ. 0) RX = XX
    4 TERM = TERM * TEMP * RX / (PP + AI)
      BETAIN = BETAIN + TERM
      TEMP = ABS(TERM)
      IF (TEMP .LE. ACU .AND. TEMP .LE. ACU * BETAIN) GOTO 5
      AI = AI + 1.0D0
      NS = NS - 1
      IF (NS .GE. 0) GOTO 3
      TEMP = PSQ
      PSQ = PSQ + 1.0D0
      GOTO 4
C
C        CALCULATE RESULT
C
    5 BETA = ALGAMA(P) + ALGAMA(Q) - ALGAMA(P+Q)
      BETAIN = BETAIN * 
     *     DEXP(PP * DLOG(XX) + (QQ - 1.0D0) * DLOG(CX) - BETA) / PP
      IF (INDEX) BETAIN = 1.0D0 - BETAIN
      RETURN
      END
C
      SUBROUTINE CHISQL(X, DF, FX, FL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C
C        DEFINE ACCURACY AND INITIALIZE
C
C        N SHOULD BE SECIFIED SUCH THAT ACC IS GREATER THAN 
C        OR EQUAL TO (AU - AL) / 2 ** N
C
C
      DATA ACC, N /1.0D-6, 30/
      AL = 0.0
      AINC = 80.0
      AU = 80.0
C
      IF (DF .LE. 0.0D0) STOP 71
      IF (X .LT. 0.0D0) STOP 72
      IF (FX .LE. 0.0D0) STOP 73
C
    1 APROX = CHI2NC(X, DF, AU)
      IF (FX .GT. APROX) GOTO 2
      IF (FX .LT. APROX) AL = AU
      AU = AU + AINC
      GOTO 1
    2 DO 3 J = 1, N
      FL = 0.5D0 * (AL + AU)
      APROX = CHI2NC(X, DF, FL)
      IF ((AU - AL)/AU .LT. ACC) GOTO 4
      IF (FX .LT. APROX) AL = FL
      IF (FX .GE. APROX) AU = FL
    3 CONTINUE
    4 RETURN
      END
C
      DOUBLE PRECISION FUNCTION CHI2NC(FQUAN, DFN, XLAM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        THE CDF OF THE NON-CENTRAL CHI-SQUARE DISTRIBUTION
C        WITH DFN DEGREES OF FREEDOM
C        AND NON-CENTRALITY PARAMETER XLAM
C        EVALUATED AT FQUAN.
C
C        CALCULATED BY SUMMING AN 'INFINITE SERIES' 
C
C        THE SERIES IS A WEIGHTED SUM OF CENTRAL CHI-SQUARE
C        PROBABILITIES 
C        WITH WEIGHTS GIVEN BY THE POISSON DISTRIBUTION WITH 
C        MEAN XLAM / 2.  SUMMATION STARTS AT XLAM / 2 AND CONTINUES
C        UP AN DOWN UNTIL THE CONTRIBUTION IN BOTH DIRECTIONS IS LESS 
C        THAN ACU.
C
C
      DATA ACU / 1.0D-10/
C
      XLAMH = XLAM / 2.0D0
      NSPLIT = XLAMH + 0.499999D0
      IF (NSPLIT .LT. 1) NSPLIT = 1
   
      CHI2NC = 0.0
      ESQ = FQUAN / 2.0D0

      XI = NSPLIT
      DO 100 I = 1, NSPLIT
      XI = XI - 1.0D0
      TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
      IF (TERM .LT. -80.0D0) GOTO 200
      TERM = DEXP(TERM)
      TERM = TERM * GAMAIN (ESQ, DFN/2.0D0 + XI)
      CHI2NC = CHI2NC + TERM
      IF (TERM .LT. ACU) GOTO 200
  100 CONTINUE

  200 XI = NSPLIT - 1

  300 XI = XI + 1.0D0
      TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
      IF (TERM .LT. -80.0D0) RETURN
      TERM = DEXP(TERM)
      TERM = TERM * GAMAIN (ESQ, DFN / 2.0D0 + XI)
      CHI2NC = CHI2NC + TERM
      IF (TERM .LT. ACU) RETURN
      GOTO 300
      END
C
      DOUBLE PRECISION FUNCTION FNCDFC(FQUAN, DFN, DFD, XLAM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        THE COMPLEMENT OF CDF FOR THE NON-CENTRAL F DISTRIBUTION
C        WITH DFN NUMERATOR DEGREES OF FREEDOM, DFD DENOMINATOR 
C        DEGREES OF FREEDOM, AND NON-CENTRALITY PARAMETER XLAM
C        EVALUATED AT FQUAN.
C
C        CALCULATED BY SUMMING AN 'INFINITE SERIES' 
C        (GRAYBILL (1961), AN INTRODUCTION TO LINEAR STATISTICAL
C        MODELS, VOLUME 1, P.79, EQU.4.5.).  
C
C        THE SERIES IS A WEIGHTED SUM OF CENTRAL F PROBABILITIES 
C        WITH WEIGHTS GIVEN BY THE POISSON DISTRIBUTION WITH 
C        MEAN XLAM / 2.  SUMMATION STARTS AT XLAM / 2 AND CONTINUES
C        UP AN DOWN UNTIL THE CONTRIBUTION IN BOTH DIRECTIONS IS LESS 
C        THAN ACU.
C
C        NMAX -- MAXIMUM NUMBER OF TERMS SUMMED
C
      DATA ACU / 1.0D-10/
C
      XLAMH = XLAM / 2.0D0
      NSPLIT = XLAMH + 0.499999D0
      IF (NSPLIT .LT. 1) NSPLIT = 1
   
      FNCDFC = 1.0
      ESQ = (DFN * FQUAN) / (DFD + DFN * FQUAN)

      XI = NSPLIT
      DO 100 I = 1, NSPLIT
      XI = XI - 1.0
      TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
      IF (TERM .LT. -80.0D0) GOTO 200
      TERM = DEXP(TERM)
      TERM = TERM * BETAIN (ESQ, DFN / 2.0D0 + XI, DFD / 2.0D0)
      FNCDFC = FNCDFC - TERM
      IF (TERM .LT. ACU) GOTO 200
  100 CONTINUE

  200 XI = NSPLIT - 1

  300 XI = XI + 1.0D0
      TERM = -XLAMH + XI * DLOG(XLAMH) - ALGAMA (XI + 1.0D0)
      IF (TERM .LT. -80.0D0) RETURN
      TERM = DEXP(TERM)
      TERM = TERM * BETAIN (ESQ, DFN / 2.0D0 + XI, DFD / 2.0D0)
      FNCDFC = FNCDFC - TERM
      IF (TERM .LT. ACU) RETURN
      GOTO 300
      END
C
      DOUBLE PRECISION FUNCTION GAMAIN(X,P)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        ALGORITHM AS 32 J.R.STATIST.SOC. C. (1970) VOL.19 NO.3
C
C        COMPUTES INCOMPLETE GAMMA RATIO FOR POSITIVE VALUES OF
C        ARGUMENTS X AND P.  G MUST BE SUPPLIED AND SHOULD BE EQUAL TO
C        LN(GAMMA(P)).
C        IFAULT = 1 IF P.LE.0 ELSE 2 IF X.LT.0 ELSE 0
C        USES SERIES EXPANSION IF P.GT.X OR X.LE.1, OTHERWISE A 
C        CONTINUED FRACTION APPROXIMATION.C
C
C        [FAULT INDICATOR REMOVED BY G.E.D
C        CALCULATION OF G INSERTED BY G.E.D.]
C
      DIMENSION PN(6)
C
C        DEFINE ACCURACY AND INITIALIZE
C
      ACU = 1.0D-8
      OFLO = 1.0D30
      GIN = 0.0
C
C        TEST FOR ADMISSIBILITY OF ARGUMENTS
C
      IF (P .LE. 0.0D0) STOP 21
      IF (X .LT. 0.0D0) STOP 22
      IF (X .EQ. 0.0D0) GOTO 50
      G = ALGAMA(P)
      FACTOR = DEXP(P * DLOG(X) - X - G)
      IF (X .GT. 1.0D0 .AND. X .GE. P) GOTO 30
C
C        CALCULATION BY SERIES EXPANSION
C
      GIN = 1.0
      TERM = 1.0
      RN = P
   20 RN = RN + 1.0D0
      TERM = TERM * X / RN
      GIN = GIN + TERM
      IF (TERM .GT. ACU) GOTO 20
      GIN = GIN * FACTOR / P
      GOTO 50
C
C        CALCULATION BY CONTINUED FRACTION
C
   30 A = 1.0D0 - P
      B = A + X + 1.0D0
      TERM = 0.0
      PN(1) = 1.0D0
      PN(2) = X
      PN(3) = X + 1.0D0
      PN(4) = X * B
      GIN = PN(3) / PN(4)
   32 A = A + 1.0D0
      B = B + 2.0D0
      TERM = TERM + 1.0D0
      AN = A * TERM
      DO 33 I = 1, 2
   33 PN(I+4) = B * PN(I+2) - AN * PN(I)
      IF (PN(6) .EQ. 0.0D0) GOTO 35
      RN = PN(5) / PN(6)
      DIF = ABS(GIN - RN)
      IF (DIF .GT. ACU) GOTO 34
      IF (DIF .LE. ACU * RN) GOTO 42
   34 GIN = RN
   35 DO 36 I = 1, 4
   36 PN(I) = PN(I+2)
      IF (ABS(PN(5)).LT.OFLO) GOTO 32
      DO 41 I = 1, 4
   41 PN(I) = PN(I) / OFLO
      GOTO 32
   42 GIN = 1.0D0 - FACTOR * GIN
C
   50 GAMAIN = GIN
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION GAUINV(P,IFAULT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C              ALGORITHM AS 70 APPL. STATIST. (1974) VOL. 23, NO.1
C              GAUINV FINDS PERCENTAGE POINTS OF THE NORMAL DISTRIBUTION
C
      DATA ZERO,ONE,HALF,ALIMIT/0.0D0, 1.0D0, 0.5D0, 1.0D-20/
C
      DATA               P0,    P1,               P2,                P3
     1   / -.322232431088D0, -1.D0, -.342242088547D0, -.204231210245D-1/
C
      DATA                P4,               Q0,              Q1
     1   / -.453642210148D-4, .993484626060D-1, .588581570495D0/
C
      DATA              Q2,              Q3,              Q4
     1   / .531103462366D0, .103537752850D0, .38560700634D-2/
C
      IFAULT=1
      GAUINV=ZERO
      PS=P
      IF (PS.GT.HALF) PS=ONE-PS
      IF (PS.LT.ALIMIT) RETURN
      IFAULT=0
      IF (PS.EQ.HALF) RETURN
      YI=DSQRT(DLOG(ONE/(PS*PS)))
      GAUINV=YI+((((YI*P4+P3)*YI+P2)*YI+P1)*YI+P0)
     1         /((((YI*Q4+Q3)*YI+Q2)*YI+Q1)*YI+Q0)
      IF (P.LT.HALF) GAUINV=-GAUINV
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION PPCHI2(P,V)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        ALGORITHM AS 91  APPL. STATIST. (1975) VOL.24, NO.3
C
C        TO EVALUATE THE PECENTAGE POINTS OF THE CHI-SQUARED
C        PROBABILITY DISTRIBUTION FUNCTION.
C        P MUST LIE IN THE RANGE 0.000002 TO 0.999998, V MUST BE POSITIVE,
C        G MUST BE SUPPLIED AND SHOULD BE EQUAL TO LN(GAMMA(V/2.0))
C
C        [FAULT INDICATOR REMOVED BY G.E.D.  
C        CALCULATION OF G ADDED BY G.E.D.]
C
      DATA E, AA /0.5D-6, 0.6931471805D0/
C
C        AFTER DEFINING THE ACCURACY AND LN(2), TEST ARGUMENTS AND INITIALIZE
C
      PPCHI2 = -1.0
      IF (P .LT. 0.000002D0 .OR. P .GT. 0.999998D0) STOP 31
      IF (V .LE. 0.0D0) STOP 32
      G = ALGAMA(V / 2.0D0)
      XX = 0.5D0 * V
      C = XX - 1.0D0
C
C        STARTING APPROXIMATION FOR SMALL CHI-SQUARED
C
      IF (V .GE. -1.24D0 * DLOG(P)) GOTO 1
      CH = (P * XX * DEXP(G + XX * AA)) ** (1.0D0 / XX)
      IF (CH - E) 6, 4, 4
C
C        STARTING APPROXIMATION FOR V LESS THAN OR EQUAL TO 0.32
C
    1 IF (V .GT. 0.32D0) GOTO 3
      CH = 0.4
      A = DLOG(1.0D0 - P)
    2 Q = CH
      P1 = 1.0D0 + CH * (4.67D0 + CH)
      P2 = CH * (6.73D0 + CH * (6.66D0 + CH))
      T = -0.5D0 + (4.67D0 + 2.0D0 * CH) / P1 -
     *   (6.73D0 + CH * (13.32D0 + 3.0D0 * CH)) / P2
      CH = CH - (1.0D0 - DEXP(A + G + 0.5D0 * CH + C * AA) * P2 / P1)/T
      IF (ABS(Q / CH - 1.0D0) - 0.01D0) 4, 4, 2
C
C        CALL TO ALGORITHM AS 70 - NOTE THAT P HAS BEEN TESTED ABOVE
C
    3 X = GAUINV(P,IF1)
C
C        STARTING APPROXIMATION USING WILSON AND HILFERTY ESTIMATE
C
      P1 = (2.0D0 / 9.0D0) / V
      CH = V * (X * DSQRT(P1) + 1.0D0 - P1) ** 3
C
C        STARTING APPROXIMATION FOR P TENDING TO 1
C
      IF (CH .GT. 2.2D0 * V + 6.0D0)
     *   CH = -2.0D0 * (DLOG(1.0D0 - P) - C * DLOG(0.5D0 * CH) + G)
C
C        CALL TO ALGORITHM AS 32 AND CALCULATION OF SEVEN 
C        TERM TAYLOR SERIES
C
    4 Q = CH
      P1 = 0.5D0 * CH
      P2 = P - GAMAIN(P1, XX)
      T = P2 * DEXP(XX * AA + G + P1 - C * DLOG(CH))
      B = T / CH
      A = 0.5D0 * T - B * C
      S1 = (210.0D0+A*(140.0D0+A*(105.0D0+A*
     *   (84.0D0+A*(70.0D0+60.0D0*A))))) / 420.0D0
      S2 = (420.0D0+A*(735.0D0+A*(966.0D0+A*(1141.0D0+1278.0D0
     *   *A)))) / 2520.0D0
      S3 = (210.0D0 + A * (426.0D0 + A * (707.0D0 + 932.0D0 
     *   * A))) / 2520.0D0
      S4 =(252.0D0+A*(672.0D0+1182.0D0*A)+C*(294.0D0+A*
     *   (889.0D0+1740.0D0*A)))/5040.0D0
      S5 = (84.0D0 + 264.0D0 * A + C * (175.0D0 
     *   + 606.0D0 * A)) / 2520.0D0
      S6 = (120.0D0 + C * (346.0D0 + 127.0D0 * C)) / 5040.0D0
      CH = CH+T*(1.0+0.5D0*T*S1-B*C*(S1-B*(S2-B
     *   *(S3-B*(S4-B*(S5-B*S6))))))
      IF (ABS(Q / CH - 1.0D0) .GT. E) GOTO 4
C
    6 PPCHI2 = CH
      RETURN
      END

      DOUBLE PRECISION FUNCTION RCRIT(ALPHA,F2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      F1 = 1.0
      AA = F2 / 2.0D0
      BB = F1 / 2.0D0
      XALPHA = XINBTA(AA, BB, ALPHA)
      FALPHA = (F2 / F1) *((1.0D0 - XALPHA) / XALPHA)
      RCRIT = DSQRT(FALPHA / (FALPHA + F2))
      RETURN
      END

      DOUBLE PRECISION FUNCTION RDF(NOBS, R, RHO, IFAULT)
C
C        THE INTEGRAL OF THE DISTRIBUTION FUNCTION OF THE CORRELATION 
C        COEFFICIENT IN A SAMPLE OF SIZE 'N' AND POPULATION VALUE 'RHO'
C        FROM -1 TO 'R'
C
C        INCORPORATES A BINARY SEARCH TO SHORTEN THE LIMITS OF NUMERICAL
C        INTEGRATION
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA PI /3.141592653589793238D0/, NBIN /4/, SMALL /1.0D-10/
C
C        CHECK ARGUMENTS
C
      IFAULT = 0
      IF (NOBS.LE.2) IFAULT = 1
      IF (DABS(R).GT.1.0D0) IFAULT = 2
      IF (DABS(RHO).GE.1.0D0) IFAULT = 3
      IF (IFAULT.GT.0) RETURN
C
C        PATHOLOGICAL CASES
C
      RDF = (R + 1.0D0) / 2.0D0
      IF (DABS(R).EQ.1.0D0) RETURN
C
      IF (NOBS - 4) 10, 20, 30
C
C         THREE OBSERVATIONS
C
   10 RDF = (DACOS(-R) - RHO * DSQRT((1.0D0 - R**2) / 
     *   (1.0D0 - RHO**2 * R**2)) * DACOS(-RHO * R)) / PI
      RETURN
C  
C         FOUR OBSERVATIONS
C
   20 RDF = (DACOS(RHO) + R * DACOS(-RHO * R) *
     *   ((1.0D0 - RHO**2) / (1.0D0 - RHO**2 * R**2))**1.5) / PI
      IF (RHO.NE.0.0D0) RDF = RDF + ((1.0D0 - RHO**2)**1.5 /
     *   (1.0D0 - RHO**2 * R**2) - DSQRT(1.0D0 - RHO**2)) /
     *   (PI * RHO)
      RETURN
C
C         FIVE OR MORE OBSERVATIONS
C
   30 IF (R.GT.RHO) GOTO 100
      RDF = 0.0
      IF (RORD(NOBS, R, RHO).LE.SMALL) RETURN
C
C         LOWER TAIL
C
C         BINARY SEARCH
C
      XLO = -1.0
      XHI = R
C
      DO 40 I = 1, NBIN
      XM = (XLO + XHI) / 2.0D0
      YM = RORD(NOBS, XM, RHO)
      IF (YM.GT.SMALL) XHI = XM
      IF (YM.LE.SMALL) XLO = XM
   40 CONTINUE
C
      RDF = SIMPSN(XLO,R,NOBS,RHO)
      RETURN
C
C         UPPER TAIL
C
C         BINARY SEARCH
C
  100 RDF = 1.0
      IF (RORD(NOBS, R, RHO).LE.SMALL) RETURN
      XLO = R
      XHI = 1.0
C
      DO 140 I = 1, NBIN
      XM = (XLO + XHI) / 2.0D0
      YM = RORD(NOBS, XM, RHO)
      IF (YM.GT.SMALL) XLO = XM
      IF (YM.LE.SMALL) XHI = XM
  140 CONTINUE
C
      RDF = 1.0D0 - SIMPSN(R,XHI,NOBS,RHO)
      RETURN
      END

      DOUBLE PRECISION FUNCTION RORD(NOBS, R, RHO)

C
C        ORDINATE OF THE DISTRIBUTION OF THE CORRELATION COEFFICIENT
C        IN A SAMPLE OF SIZE  'N' AND POPULATION VALUE 'RHO' AT 'R'
C
C        F.N. DAVID RECURRENCE FORMULA
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA PI /3.141592653589793238D0/
C
      RORD = 0.0
      IF (R.EQ.-1.0D0 .OR. R.EQ.1.0D0) RETURN
C
      R2 = R ** 2
      RHO2 = RHO ** 2
      RRHO = R * RHO
      R2RHO2 = R2 * RHO2
      X1 = RRHO * DSQRT((1.0D0 - RHO2) * (1.0D0 - R2)) /
     *   (1.0D0 - R2RHO2)
      X2 = (1.0D0 - RHO2) * (1.0D0 - R2) / (1.0D0 - R2RHO2)

C
      YNM2 = (1.0D0 - RHO2) / (PI * DSQRT(1.0D0 - R2)) *
     *   (1.0D0 / (1.0D0 - R2RHO2) + RRHO * DACOS(-RRHO) /
     *   (1.0D0 - R2RHO2)**1.5)
      YNM1 = (1.0D0 - RHO2)**1.5 / PI * (3.0D0 * RRHO /
     *   (1.0D0 - R2RHO2)**2 + (1.0D0 + 2.0D0 * R2RHO2) * 
     *   DACOS(-RRHO) / (1.0D0 - R2RHO2)**2.5)
C
      DO 100 I = 5, NOBS
      XI = I
      YN = (2.0D0 * XI - 5.0D0) / (XI - 3.0D0) * X1 * YNM1
     *   + (XI - 3.0D0) / (XI - 4.0D0) * X2 * YNM2
      YNM2 = YNM1
      YNM1 = YN
  100 CONTINUE
C
      RORD = YN
      RETURN
      END

      DOUBLE PRECISION FUNCTION SIMPSN(XLO,XHI,NOBS,RHO)
C
C        SIMPSON'S RULE... THE INTEGRAL OF F(X) FROM
C           X = XLO TO X = XHI
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA ERR /1.0D-6/, RERR /1.0D-4/
C
C        AT EACH STEP ADD NEW POINTS BETWEEN THE OLD POINTS,
C        WHEN A SET OF POINTS FIRST ENTERS IT COMES IN WITH
C           A COEFFICIENT OF 4 .  AT THE NEXT AND ALL 
C           SUBSEQUENT STEPS IT HAS A COEFFICIENT OF 2.
C
      LOOP = -1
      RANGE = XHI - XLO
      H = RANGE
      SUM = RORD(NOBS,XLO,RHO) + RORD(NOBS,XHI,RHO)
      SIMPSN = 0.0
      IF (SUM.LT.ERR) RETURN
      SIMP0 = H / 3.0D0 * SUM
C
   10 LOOP = LOOP + 1
      NSTOP = 2 ** LOOP
      X2STOP = 2 * NSTOP
      H = H / 2.0D0
      TERM = 0.0
C
      DO 100 I = 1, NSTOP
      XORD = XLO + DBLE(2 * I - 1) / X2STOP * RANGE
      TERM = TERM + RORD(NOBS,XORD,RHO)
  100 CONTINUE
C
      SUM = SUM + 4.0D0 * TERM
      SIMPSN = H / 3.0D0 * SUM
      DIFF = DABS(SIMPSN - SIMP0)
      IF (DIFF.LE.ERR .AND. DIFF/SIMPSN.LE.RERR) RETURN
C
      SIMP0 = SIMPSN
      SUM = SUM - 2.0D0 * TERM
      GOTO 10
      END

C
      DOUBLE PRECISION FUNCTION XINBTA(P, Q, ALPHA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C        ALGORITHM AS 109 APPL.STATIST. (1977), VOL.26, NO.1
C        (REPLACING ALGORITHM AS 64  APPL. STATIST. (1973),
C        VOL.22, NO.3)
C        ***[FAULT INDICATOR REMOVED BY G.E.D.]***
C
C        COMPUTES INVERSE OF INCOMPLETE BETA FUNCTION
C        RATIO FOR GIVEN POSITIVE VALUES OF THE ARGUMENTS
C        P AND Q, ALPHA BETWEEN ZERO AND ONE.
C        LOG OF COMPLETE BETA FUNCTION, BETA, IS ASSUMED TO BE KNOWN.
C        ***[CALCULATION OF BETA ADDED BY G.E.D]***
C
      LOGICAL INDEX
C
C        DEFINE ACCURACY AND INITIALIZE
C
      DATA ACU /1.0D-14/
      XINBTA = ALPHA
C
C        TEST FOR ADMISIBILITY OF PARAMETERS
C
      IF (P.LE.0.0D0.OR.Q.LE.0.0D0) STOP 50
      IF (ALPHA.LT.0.0D0.OR.ALPHA.GT.1.0D0) STOP 51
      IF (ALPHA.EQ.0.0D0.OR.ALPHA.EQ.1.0D0) RETURN
C
      BETA=ALGAMA(P)+ALGAMA(Q)-ALGAMA(P+Q)
C
C        CHANGE TAIL IF NECESSARY
C
      IF (ALPHA .LE. 0.5D0) GOTO 1
      A = 1.0D0 - ALPHA
      PP = Q
      QQ = P
      INDEX = .TRUE.
      GOTO 2
    1 A = ALPHA
      PP = P 
      QQ = Q
      INDEX = .FALSE.
C
C        CALCULATE INITIAL APPROXIMATION
C
    2 R = DSQRT(-DLOG(A * A))
      Y = R - (2.30753D0 + 0.27061D0 * R) / 
     *             (1.0D0 + (0.99229D0 + 0.04481D0 * R) * R)
      IF (PP .GT. 1.0D0 .AND. QQ .GT. 1.0D0) GOTO 5
      R = QQ + QQ
      T = 1.0D0 / (9.0D0 * QQ)
      T = R * (1.0D0 - T + Y * DSQRT(T)) ** 3
      IF (T .LE. 0.0D0) GOTO 3
      T = (4.0D0 * PP + R - 2.0D0) / T
      IF (T .LE. 1.0D0) GOTO 4
      XINBTA = 1.0D0 - 2.0D0 / (T + 1.0D0)
      GOTO 6
    3 XINBTA = 1.0D0 - DEXP((DLOG((1.0D0 - A) * QQ) + BETA) / QQ)
      GOTO 6
    4 XINBTA = DEXP((DLOG(A * PP) + BETA) / PP)
      GOTO 6
    5 R = (Y * Y - 3.0D0) / 6.0D0
      S = 1.0D0 / (PP + PP - 1.0D0)
      T = 1.0D0 / (QQ + QQ - 1.0D0)
      H = 2.0D0 / (S + T)
      W = Y * DSQRT(H + R) / H - (T - S) * 
     *   (R + 5.0D0 / 6.0D0 - 2.0D0 / (3.0D0 * H))
      XINBTA = PP / (PP + QQ * DEXP(W + W))
C
C        SOLVE FOR X BY A MODIFIED NEWTON-RAPHSON METHOD,
C        USING THE FUNCTION BETAIN.
C
    6 R = 1.0D0 - PP
      T = 1.0D0 - QQ
      YPREV = 0.0
      SQ = 1.0
      PREV = 1.0
      IF (XINBTA .LT. 0.0001D0) XINBTA = 0.0001D0
      IF (XINBTA .GT. 0.9999D0) XINBTA = 0.9999D0
    7 Y = BETAIN(XINBTA, PP, QQ)
      Y = (Y - A) * DEXP(BETA +
     *             R * DLOG(XINBTA) + T * DLOG(1.0D0 - XINBTA))
      IF (Y * YPREV .LE. 0.0D0) PREV = SQ
      G = 1.0
    9 ADJ = G * Y
      SQ = ADJ * ADJ
      IF (SQ .GE. PREV) GOTO 10
      TX = XINBTA - ADJ
      IF (TX .GE. 0.0D0 .AND. TX .LE. 1.0D0) GOTO 11
   10 G = G / 3.0D0
      GOTO 9
   11 IF (PREV .LE. ACU) GOTO 12
      IF (Y * Y .LE. ACU) GOTO 12
      IF (TX .EQ. 0.0D0 .OR. TX .EQ. 1.0D0) GOTO 10
      IF (TX .EQ. XINBTA) GOTO 12
      XINBTA = TX
      YPREV = Y
      GOTO 7
   12 IF (INDEX) XINBTA = 1.0D0 - XINBTA
      RETURN
      END
