program pocock

C      PROGRAM TO COMPUTE CRITICAL VALUES FOR POCOCK'S METHOD

       implicit double precision (a-h, o-z)

       character *59 input_query(6)
       character *55 out_note(2)

       data input_query/
     $   'Number of looks:                                    ',
     $   'Number of grid points in numerical integration:     ',
     $   'Significance level:                                 ',
     $   'One-sided (0) or two-sided (1):                     ',
     $   'Lower and upper search bounds for critical level:   ',
     $   'Bisection stopping criterion:                       '/

       data out_note/
     $   'Critical level:                                 ',
     $   'Computed Type I error probability:              '/

C      READ IN INPUT PARAMETERS
       write(6,*)
       write(6,*) 
     $   'PROGRAM TO COMPUTE CRITICAL LEVELS FOR THE POCOCK METHOD'
       write(6,*)
       write(6,881) input_query(1)
       read(5,*) k
       write(6,881) input_query(2)
       read(5,*) l
       write(6,881) input_query(3)
       read(5,*) alpha
       write(6,881) input_query(4)
       read(5,*) iside
       write(6,881) input_query(5)
       read(5,*) cl, cu
       write(6,881) input_query(6)
       read(5,*) tol
 881   format(x,a,x,$)
       write(6,*)     
       gam = 1.0d0 - alpha
            
C      SIMPLE BISECTION SEARCH FOR CRITICAL LEVEL
       iflg = 0
       do while (iflg .eq. 0)
         c = (cl+cu)/2.0d0
         call prcmp(c,k,l,iside,pr)
         if (pr .le. gam)
     $     then 
             cl = c
           else
             cu = c
          endif
         diff = cu-cl
         if (diff .le. tol) iflg=1
        enddo

C      PRINTOUT
       pr = 1.0d0 - pr
       write(6,882) out_note(1), c
       write(6,882) out_note(2), pr
       write(6,*)
 882   format(x,a,f8.4)
       write(6,*)
     $   'If the computed probability is not close to the desired ',
     $   'probability,'
       write(6,*)
     $   'try changing the initial bisection search limits or the ',
     $   'bisection '
       write(6,*)
     $   'stopping criterion.'
       write(6,*)

       end
       

C      SUBROUTINE TO CALCULATE THE PROBABILITY 
C      Pr(Z_k < c for all k) (one-sided case)
C      OR Pr(|Z_k| < c for all k) (two-sided case)
C      RECURSIVE ALGORITHM AS DESCRIBED IN CLASS
C      NUMERICAL INTEGRATION USING SIMPSON'S RULE
C      ARRAY z CONTAINS THE GRID POINTS
C      ARRAY w CONTAINS THE SIMPSON WEIGHTS
C      WE TAKE vmininf AS A LOWER CUTOFF POINT FOR THE INTEGRAL

       subroutine prcmp(c,k,l,iside,prob)

       implicit double precision (a-h, o-z)       
       parameter (maxgrid=5001, maxk=20, vmininf=-6.0d0)       

       dimension z(maxgrid), w(maxgrid), f(maxk,maxgrid)

       a = vmininf
       if (iside .eq. 1) a=-c
       b = c
       dl = l
       h = (b-a)/dl

C     SET UP THE GRID POINTS AND WEIGHTS
       do ll = 1, l+1
         dll = ll
         z(ll) = a + (dll-1.0d0)*h
         mm = mod(ll,2)
         dm = mm
         w(ll) = 2.0d0*(1.0d0+dm)
         f(1,ll) = phi(z(ll))
        enddo
       w(1) = 1.0d0
       w(l+1) = 1.0d0

C     CARRY OUT THE RECURSIVE CALCULATION OF THE F FUNCTIONS
       do kk = 2, k
         dkk = kk
         rootkk = dsqrt(dkk)
         rkkprv = dsqrt(dkk-1.0d0)
         do ll1 = 1, l+1
           f(kk,ll1) = 0.0d0
           do ll2 = 1, l+1
             pharg = rootkk*z(ll1)-rkkprv*z(ll2)
             addtrm = phi(pharg)*f(kk-1,ll2)*rootkk
             f(kk,ll1) = f(kk,ll1) + w(ll2)*addtrm
            enddo
           f(kk,ll1) = (h/3.0d0)*f(kk,ll1)
          enddo
        enddo

C      CALCULATE THE FINAL PROBABILITY
        prob = 0.0d0
        do ll = 1, l
          prob = prob + w(ll)*f(k,ll)
         enddo
        prob = (h/3.0)*prob

        end


C       STANDARD NORMAL DENSITY FUNCTION
        double precision function phi(u)
        double precision u
        phi = dexp(-0.5d0*(u**2))/2.506628275d0
        return
        end