PROGRAM RINOTT
C
C     PROGRAM TO CALCULATE G(T,PSTAR,NU) OF RINOTT PROCEDURE
C
      DOUBLE PRECISION ZCDF,CHIPDF
      DOUBLE PRECISION PSTAR,H
      DOUBLE PRECISION X(32),W(32),WEX(32),DUMMY,LNGAM(50),
     *LOWERH,UPPERH,ANS,TMP
      INTEGER I,T,NU,LOOPH
      EXTERNAL ZCDF
C
      INCLUDE 'GLQUAD'
C
C     CALCULATE LNGAM(K) = LN(GAMMA(K/2))
C
      LNGAM(1) = 0.5723649429
      LNGAM(2) = 0.0
      DO 100 I = 2,25
      LNGAM(2*I - 1) = DLOG(DBLE(I)-1.50D0) + LNGAM(2*I - 3)
      LNGAM(2*I) = DLOG(DBLE(I) - 1.0D0) + LNGAM(2*I - 2)
100   CONTINUE
C
C    READ PARAMETERS FOR WHICH TO DETERMINE G(T,PSTAR,\NU)
C
110   WRITE(6,*) 'ENTER T,P-STAR AND NU'
      READ(5,*) T,PSTAR,NU
      IF (NU .LT. 5) THEN
        WRITE(6,520)
        STOP
      ELSE
           IF (T .LT. 2) THEN
             WRITE(6,530)
             STOP
           ENDIF
      ENDIF
      DUMMY = 1.0D0
C
C
      H = 4.0D0
      LOWERH = 0.0D0
      UPPERH = 20.00
C
      DO 120 LOOPH = 1,50
C
        ANS = 0.0D0
        DO 130 J = 1,32
          TMP = 0.0D0
          DO 140 I = 1,32
            TMP = TMP + WEX(I)
     *      * ZCDF(H/DSQRT(DBLE(NU)*(1.0D0/X(I) + 1.0D0/X(J))/DUMMY))
     *      * CHIPDF(NU,DUMMY*X(I),LNGAM) * DUMMY
140       CONTINUE
          TMP = TMP ** (T-1)
130       ANS = ANS + WEX(J) * TMP * CHIPDF(NU,DUMMY*X(J),LNGAM)*DUMMY
C
          IF(DABS(ANS-PSTAR) .GT. 0.000001D0) GOTO 150
C
C   IF NOT, THEN DONE SO WRITE G(T,PSTAR,\NU)
C
          WRITE(6,500) PSTAR,T,NU,H
          WRITE(6,*) 'DO YOU WANT TO CONTINUE? PLEASE ANSWER Y/N'
          READ(5,510) CHOICE2
          IF (CHOICE2 .EQ. 'Y' .OR. CHOICE2 .EQ. 'y') GOTO 110
        STOP
150   IF(ANS .GT. PSTAR) GOTO 160
      LOWERH = H
      H = (LOWERH + UPPERH)/2.0D0
      GOTO 120
160   UPPERH = H
      H = (LOWERH + UPPERH)/2.0D0
120   CONTINUE
C
500   FORMAT(1X,' G^{',F5.3,'}_{',I3,',',I3,'} = ',F6.3)
510   FORMAT(A1)
520   FORMAT(1X,'THIS PROGRAM WILL NOT COMPUTE  G  FOR LESS',/,
     *' THAN 5 DEGREES OF FREEDOM')
530   FORMAT(1X,' T MUST BE AN INTEGER GREATER THAN OR EQUAL TO 2')
C
      END
C
C
      DOUBLE PRECISION FUNCTION CHIPDF(N,C,LNGAM)
C
C     THIS GIVES THE PDF OF THE CHI^2 DISTN WITH N DF FOR N .LE. 50.
C     LNGAM(N) IS LN(GAMMA(N/2))
C
      DOUBLE PRECISION C,TMP,LNGAM(50),FLN2
      INTEGER N
      FLN2 = DBLE(N)/2.
      TMP = -FLN2*DLOG(2.0D0) - LNGAM(N) +
     *    (FLN2 - 1.0D0)*DLOG(C) - C/2.0D0
      CHIPDF = DEXP(TMP)
      RETURN
      END