program pocock

C      PROGRAM TO PERFORM COMPUTATIONS FOR POCOCK'S METHOD

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

       character *59 input_query(9)
       character *55 out_note(5)

       data input_query/
     $   'Number of looks:                                    ',
     $   'Number of grid points in numerical integration:     ',
     $   'Significance level:                                 ',
     $   'Lower and upper search bounds for critical level:   ',
     $   'Bisection stopping criterion for Type I error:      ',
     $   'Bisection stopping criterion for power:             ',
     $   'Per-Group Sample Size:                              ',
     $   'Variance:                                           ',
     $   'Desired Power:                                      '/

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

C      READ IN INPUT PARAMETERS
       write(6,*)
       write(6,*)
     $   'PROGRAM TO PERFORM COMPUTATIONS 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,*) cl, cu
       write(6,881) input_query(5)
       read(5,*) tol1
       write(6,881) input_query(6)
       read(5,*) tol2
       write(6,881) input_query(7)
       read(5,*) nsam
       write(6,881) input_query(8)
       read(5,*) sig2
       write(6,881) input_query(9)
       read(5,*) power
  881   format(1x,a,1x,$)
       write(6,*)

C      SIMPLE BISECTION SEARCH FOR CRITICAL LEVEL
       gam = 1.0d0 - alpha
       theta = 0
       iflg = 0
       do while (iflg .eq. 0)
         c = (cl+cu)/2.0d0
         call prcmp(c,k,l,theta,pr)
         if (pr .le. gam)
     $     then
             cl = c
           else
             cu = c
          endif
         diff = cu-cl
         if (diff .le. tol1) iflg=1
        enddo
        alf = 1.0d0 - pr

C      SIMPLE BISECTION SEARCH FOR THETA VALUE THAT GIVES DESIRED POWER
       gam = 1.0d0 - power
       th_l = 0.0d0
       th_u = 5.0d0
       iflg = 0
       do while (iflg .eq. 0)
         theta = (th_l+th_u)/2.0d0
         call prcmp(c,k,l,theta,pr)
         if (pr .ge. gam)
     $     then
             th_l = theta
           else
             th_u = theta
          endif
         diff = th_u-th_l
         if (diff .le. tol2) iflg=1
        enddo
        pow = 1.0d0 - pr
        delta = theta/dsqrt(nsam/(2.0d0*sig2))

C      PRINTOUT
       pr = 1.0d0 - pr
       write(6,882) out_note(1), c
       write(6,882) out_note(2), alf
       write(6,882) out_note(3), pow
       write(6,882) out_note(4), theta
       write(6,882) out_note(5), delta
       write(6,*)
 882   format(1x,a,f8.4)
       write(6,*)
     $   'If the computed Type I error or power is not close ',
     $   'to the desired levels,'
       write(6,*)
     $   'try changing the initial bisection search limits or the ',
     $   'bisection stopping criterion.'
       write(6,*)

       end


C      SUBROUTINE TO CALCULATE THE PROBABILITY
C      Pr(Z_k < c - theta*sqrt(k/K) for all k)
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,theta,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
       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
             zz = z(ll2)
             upper = c - dsqrt((dkk-1.0d0)/dble(k))*theta
             if (zz .le. upper) then
               pharg = rootkk*z(ll1)-rkkprv*z(ll2)
               addtrm = phi(pharg)*f(kk-1,ll2)*rootkk
               f(kk,ll1) = f(kk,ll1) + w(ll2)*addtrm
              endif
            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