c--------------------------------------------------c
c     pmcc.f computes power                        c
c       for 2-sided MCC inference                  c
c     last revision: 02-06-97 by jch               c
c--------------------------------------------------c
      program main
      implicit doubleprecision (a-h,o-z)
c6    implicit real*8 (a-h,o-z)
c     real fact1
      real conf,sdbs
      external glv,ghc1
c     dimension fact1(100)
      dimension a(100),b(100),c(100)
      dimension ph(48),wh(48),whln(48),zh(48)
      dimension pl(48),wl(48)
      dimension ah(7),bh(12)
      common/list2/a,b,c
      common/list5/dinfnu
      common/list6/ph,wh,zh,whln
      common/list7/pl,wl
      common/list8/zero,zsr,zlog,zln
      common/list9/sr2,pi,srpi,pib2
      common/list10/iin,iout,iouti
      common/list12/ah,bh
c     ----------------------------------
c     initialize commonly used constants
c     ----------------------------------
      call commgq
c     -----------------------------------------
c     re-initialize infinite degrees of freedom
c     -----------------------------------------
      dinfnu=750d0
c     ----------------------------------------
c     initialize default tolerence (eps) on cc
c     ----------------------------------------
      eps=0.20d-04
c     --------------------------------------------------------------
c     initialize default tolerence (dstop) on half width of interval
c     --------------------------------------------------------------
      dstop=0.5d-03
c     -------------------------------------------------------
c     initialize default bound on number of iteration (itmax)
c     -------------------------------------------------------
      itmax=9
      imc=32
c     
c     -----------
c     batch input
c     -----------
c     read(5,1005) k,conf,n,sdbs
c1005 format(i3,2x,f3.2,2x,i3,2x,f5.2)
c
c     ---------------
c     prompted inputs
c     ---------------
      write(6,1010)
 1010 format(' ','input number of treatments')
      read(5,*) k
      write(6,1020)
 1020 format(' ','input confidence coefficient')
      read(5,*) conf
c     write(6,1060)
c1060 format(' ','input array of factor pattern = sqrt(communality)')
c     read(5,*) (fact1(j),j=1,km1)
      write(6,1070)
 1070 format(' ','input sample size n')
      read(5,*) n
c
      write(6,1080)
 1080 format(' ','input dbs')
      read(5,*) sdbs
c
      km1=k-1
      cc=dble(conf)
      nu=k*(n-1)
      dnu=dble(float(nu))
      dbs=dble(sdbs)
c     
c     -----------------------
c     initialize correlations
c     -----------------------
      do 10 j=1,km1
c       ----------------------
c       a = factor pattern
c       c = 1/sqrt(uniqueness)
c       ----------------------
        a(j)=1.0d0/dsqrt(2.0d0)
        c(j)=1.0d0/dsqrt(1.0d0-a(j)*a(j))
 10   continue
c
c     -------------------
c     find critical value
c     -------------------
c
c         -----------------------------------------------
c         obtain upper and lower bounds on critical value
c         -----------------------------------------------
          call bdb12(km1,cc,nu,dnu,dhigh,dlow)
c         call bdmcc2(k,cc,nu,dnu,dhigh,dlow)
c
c         ------------------------------------------------------------
c         compute critical value d by the modified regula falsi method
c         ------------------------------------------------------------
          call mrgfls(ghc1,imc,k,glv,cc,dhigh,dlow,w,dnu,-1.0d0,
     *                dstop,eps,itmax,iflag,iter)
          if(iflag.lt.0) stop
          d=w
c
      write(6,2000) k,conf,sdbs,n,nu
 2000 format('number of treatments =',i2,1x, /,
     *'confidence coefficient =',f4.3,1x, /,
     *'delta divided by sigma =',f10.5,1x, /, 
     *'individual sample size =',i3,1x, /,
     *'error degrees of freedom =',i3)
c
c
c     -------------------------------------------------------
c     compute power
c     d in Hsu(1988) corresponds to unstandardised statistics
c     -------------------------------------------------------
      su=dsqrt(dble(float(n)))*dbs/(dsqrt(2.0d0)*d)
c     --------------------------------------------------
c     glv integrates with respect to chi/sqrt(2) density
c     --------------------------------------------------
      u=dsqrt(dnu*0.5d0)*su
      power=glv(ghc1,imc,k,d,dnu,u)
      write(6,3000) power
 3000 format('prob(correct and useful inference)=',f5.3)
c
      stop
      end