!------------------------------------------------------------
!                
!               Get correlation function from P(k)
!                      
!------------------------------------------------------------
Module Pars
  IMPLICIT REAL*8 (A-H,O-Z)
    Integer*4, parameter :: NtabMax = 10000
    Real*8, parameter :: &
                         pi      = 3.141592653d0, &
                         h       = 0.70
Real*8,  dimension(NtabMax) :: Ppk,kt,Pp2                
Real*8                :: dlogK,dR,Rmax,Rzero,alog0,Rf,kmax,Pkmax
Real*8                :: OmBar,OmC
Integer*4             :: Ntab
character*80          :: FileName,Line
character*10           :: Txt1,Txt2      
Contains
!-----------------------------------------------------------
!
Subroutine ReadTable                                 ! read header                          
!
  Real*8 :: d0, d1, d2
         !Do i = 1, 5 
         !   READ (50, '(a)', end = 100, err = 100) Line 
         !EndDo 
          Read(50,'(a)')Line
          i0 =scan(Line,'=') ;  Line = adjustl(Line(i0+1:)) ; read(Line,*) OmBar
          Read(50,'(a)')Line
          i0 =scan(Line,'=') ;  Line = adjustl(Line(i0+1:)) ; read(Line,*) OmC
          Do i=1,3
           Read(50,'(a)')Line
          EndDo
         Ntab = 0 
         kmax = 0. ; Pkmax =0.
   12 READ (50, *, end = 32, err = 32) xx, pp 
        Ntab = Ntab + 1 
       kt (Ntab) = xx
       d1 = 2.*pi/250.
       d0 = xx/d1
       d2 = (2.+d0**4)/(1.+d0 **4)  !!!!! remove it
       !write(*,'(1p,12g12.4)') xx,d1,d0,d2,pp
       pp = d2*pp                   !!!! remove it
       Ppk (Ntab) = pp                !  Pk
       Pp2 (Ntab) = pp*xx**2          ! store k**2P(k)
       If(pp>Pkmax)Then
          kmax = xx ; Pkmax = pp
       EndIf
      GOTO 12 
   32   continue 
          Write (*,*) ' Read ', Ntab, ' lines from P(k) table '
                                                                        
      If (Ntab.le.1) stop 'wrong table for p(k)' 
      dlogK = log10 (kt (Ntab) / kt (1) ) /(Ntab-1)
      alog0 = log10(kt(1))
      !write(*,'(a,g12.4)') ' dLogK = ',dlogK
                            ! test that spacing is the same             
      Do k = 2, Ntab - 1 
        ss = log10 (kt (k + 1) / kt (k) ) 
        If (abs (ss /dlogK - 1.) .gt.2.e-2) Then 
        Write (*,*) ' error in K spacing. k=', k, kt (k + 1),kt (k)
        STOP 
        EndIf 
      EndDo 

      Return 
  100 Write (*,*) ' Error reading table of P(k)' 
      STOP 

End Subroutine ReadTable
!---------------------------------------
      REAL*8 FUNCTION Pk(x)
!                     interpolate table with p(k)
!                      
!---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)

      If(x.ge.kt(Ntab))Then  ! slope is ns =-3
         Pk =Ppk(Ntab)/(x/kt(Ntab))**3
         Return
      EndIf
      If(x.lt.kt(1))Then                ! slope is ns=1
         Pk =Ppk(1)*(x/kt(1))
         Return
      EndIf
      ind = INT((log10(x)-alog0)/dlogK)+1
      dk  = kt(ind+1)-kt(ind)
        Pk   = (Ppk(ind)*(kt(ind+1)-x)+Ppk(ind+1)*(x-kt(ind)))/dk
    !     write(*,'(i5,1p,10g12.4)')ind,x,log10(x),log10(x)-alog0, &
    !     kt(ind),x,kt(ind+1), Ppk(ind),Pk,Ppk(ind+1)

    End FUNCTION Pk
!---------------------------------------
      REAL*8 FUNCTION Pk2(x)
!                     interpolate table with p(k)
!                      
!---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)

      If(x.ge.kt(Ntab))Then  ! slope is ns =-3
         Pk2 =Pp2(Ntab)/(x/kt(Ntab))
         Return
      EndIf
      If(x.lt.kt(1))Then                ! slope is ns=1
         Pk2 =Pp2(1)*(x/kt(1))**3
         Return
      EndIf
      ind = INT((log10(x)-alog0)/dlogK)+1
      dk  = kt(ind+1)-kt(ind)
        Pk2   = (Pp2(ind)*(kt(ind+1)-x)+Pp2(ind+1)*(x-kt(ind)))/dk
    !     write(*,'(i5,1p,10g12.4)')ind,x,log10(x),log10(x)-alog0, &
    !     kt(ind),x,kt(ind+1), Ppk(ind),Pk,Ppk(ind+1)

    End FUNCTION Pk2
!----------------------------------------------------------
Real*8 Function Faux(x)
!              integrant for correlation function
  Real*8 x

    Faux = sin(x*Rf)*Pk2(x)/x
  End Function Faux   
!----------------------------------------------------------
Real*8 Function Faux2(x)
!              integrant for correlation function 2
  Real*8 x

    Faux2 = Pk2(x)/x
  End Function Faux2   

!----------------------------------------------------------
Subroutine FindPkPeak
!              find zero crossing and max of Xi(r)
  Real*8 a,b,c,P_max,xm,x0,x1,Pm,P0,P1,dd
  Integer*4 :: k_mm(1)

  k_mm = MAXLOC(Ppk)            ! position of maximum
  k_max = k_mm(1)
  x0    =  kt(k_max)
  P0    = Ppk(k_max)    
  If(k_max ==1 .or. k_max == Ntab)Then
     kmax = x0 ; Pmax =P0 ; return
  EndIf
  xm = kt(k_max-1) ; Pm =Ppk(k_max-1)
  x1 = kt(k_max+1) ; P1 =Ppk(k_max+1)

  c = (P1*(x0-xm)-P0*(x1-xm)+Pm*(x1-x0))/(x0-xm)/(x1-x0)/(x1-xm)
  b = (P1-P0)/(x1-x0) - c*(x1+x0)
  a = P1 -b*x1-c*x1**2

  Pmax = a -b**2/4./c
  kmax = -b/2./c

  print *,k_max-1,k_max,k_max+1
  print *,xm,x0,x1
  print *,Pm,P0,P1
  print *,' max =',kmax,Pmax
  end Subroutine FindPkPeak
!----------------------------------------------------------
Subroutine FindZeroPeak(R0,Rmx,Ximax)
!              find zero crossing and max of Xi(r)
  Real*8 R,dR,yy,Rp,xi0,xi1
    R = 30.     ! initial radius
    dR = 0.5   ! initial step
    Rp =0.
    xi0 = Xi2(R)
    xi1 = Xi2(R+dR)
    If(xi0>xi1)Then
       iFlag  = 0  ! branch with declining xi
    Else
       iFlag  = 1  ! 
    EndIf
    !write(*,'(3g12.5,3x,i3)') R,xi0,xi1,iFlag

    do while(R<300.)      ! find first approximation BAO peak
       R = R + dR
       yy = Xi2(R)
       If(yy>xi0.and.iFlag ==0)iFlag = 1 ! flip iFlag for rising brunch
       If(iFlag ==1)Then
          If(yy>xi0)Then
             xi0 =yy
             Rp  =R
          Else
             exit
          EndIf
       EndIf
       !write(*,'(1p12g14.5)') R,Rp,yy,xi0
       xi0 = yy
    end do

    R = 250.     ! initial radius. Find zero crossing
    R0 = 0. 
    xi0 = Xi2(R)
    xi1 = Xi2(R-dR)
    write(*,'(a,1p,2g12.4,a,2g12.4)') ' R=',R,R-dR,' xi =', xi0,xi1
    If(xi0<0.)Then
       iFlag  = 0  ! branch with R> R_peak
    Else
       iFlag  = 1  ! R< R_peak
    EndIf
    if(iFlag == 1)Stop ' STOP: Initial radius too small: less than R_peak '

    xi0 = Xi2(R)
    xi1 = Xi2(R-dR)

    xi1 = Xi2(R)
    do while(R>30.)      ! find first approximation zero point
       R = R - dR
       yy = Xi2(R)
       If(yy>0.d0.and.xi0.le.0.d0)Then
            R0 = R
            exit
       EndIf
       !write(*,'(1p12g14.5)') R,R0,Rp,yy,xi1,xi0
       xi1 = yy
       !if(R>138.)exit
    end do

    Rm = R0+2.*dR     ! Improve: Find zero point
    R = R0
    dR0 = dR/100.
    xi0 = Xi2(R)
    do while(R<=Rm)
       R = R + dR0
       yy = Xi2(R)
       If(yy<=0.)exit
    end do
    R0 = R-dR0/2.

    Rm = Rp+dR  
    R = Rp-dR      ! Improve find peak Xi(r)
    dRp = dR/100.
    xi0 = Xi2(R)
    do while(R<=Rm)
       R = R + dRp
       yy = Xi2(R)
       If(yy>xi0)Then
          xi0 =yy
          Rp  =R
       EndIf
    end do
    Rmx   = Rp
    Ximax = xi0
End Subroutine FindZeroPeak

End Module Pars
!-----------------------------------------------------------
Program Correlation
use Pars
  IMPLICIT REAL*8 (A-H,O-Z)
  logical :: ex
    !write(*,'(a,$)') ' Enter file name = '
    !read(*,*) FileName

    Call Getarg(1,FileName)                ! open files -------------
    Inquire(file=TRIM(FileName),exist=ex)
    If(.not.ex)Stop ' --- Input File does not exist ----'
    Open(50,file=TRIM(FileName))
    Call ReadTable                        ! read power spectrum
    ii = INT(OmBar/h**2*10000.+0.5)
    write(Line,'(a,i4.4,a)')'TableXiZero.OmBar.0.',ii,'.dat'
    !Line =Trim('TableXiZero.n1.OmBar.0.0500.dat')
    Inquire(file=Line,exist=ex)
    If(.not.ex)Then
       Open(1,file=Line)
       write(1,'(2(3x,a))') 'OmBarh**2     OmDMh**2     OmLambda', &
                      '     Ro(Mpch)     R(xi=peak)    Xi(peak)      kmax       P(kmax)'
    Else
       Open(1,file=Line,position='append')
    EndIf

    close(50)
    Call FindZeroPeak(R0,Rmx,Ximax)       ! find zero crossing and BAO peak
    Call FindPkPeak
    
    write(*,'(2(3x,a))') 'OmBarh**2     OmDMh**2     OmLambda', &
                      '     Ro(Mpch)     R(xi=peak)    Xi(peak)'
    write(*,'(1p,8g14.6)') OmBar,OmC+OmBar,1.-(OmC+OmBar)/h**2,R0,Rmx,Ximax,kmax,Pkmax
    write(1,'(1p,8g14.6)') OmBar,OmC+OmBar,1.-(OmC+OmBar)/h**2,R0,Rmx,Ximax,kmax,Pkmax
    close (1)

    Line='xi'//FileName(6:)               ! write correlation function
    open(10,file=TRIM(Line))
    write(10,'(a,3x,a)') '#   OmBarh**2     OmDMh**2     OmLambda', &
                      '     Ro(Mpch)     R(xi=peak)    Xi(peak)'
    write(10,'("# ",1p,10g14.6)') OmBar,OmC+OmBar,1.-(OmC+OmBar)/h**2,R0, &
                                 Rmx,Ximax,kmax,Pkmax
    write(10,'(a)') '# R(Mpch)     Xi(R)' 
    r = 1. ; dr =1.
    Do while(r<150.)
       xx =Xi2(r)
       !yy =XiTr(r)
       !zz =Xi(r)
       write(*,'(1p,8g14.6)') r,xx,xx*r**2
       r = r+dr
    End Do
end Program Correlation
!----------------------------------------------------------
Real*8 Function Xi(R)
          !    correlation function
  use Pars
  Real*8, parameter :: twoPisq= 5.06606d-2, twopi =6.283185d0
  Real*8 x,r,ANS,ANSa,ANSb,ANSc,ANSd,er,A,B,XupA,XupB,XupC,XupD
 !External Faux

  Rf = R
  Er = 1.d-8
  B   = twopi/Rf
  A   = B
     Call DGAUS8 (Faux, 6.28d-3, A, ER, ANSa, IERR)
  XupA = A
  Do i=1,1000
     XupB = XupA + 10.d0*B
     Call DGAUS8 (Faux, XupA, XupB,ER, ANSb, IERR)
     ANSa = ANSa + ANSb
     If(abs(ANSb)<1.d-6*abs(ANSa))exit
     !write(*,'(1p,10g15.5)') XupA,XupB,ANSb,ANSa,ER
     XupA = XupB
  EndDo
     Xi = twoPisq*(ANSa)/Rf

    !write(*,'(1p,16g12.4)')R,ANSa,ANSb,ANSc,ANSd,ANS,ANSa+ANSb+ANSc+ANSd+ANS
   End Function Xi
!----------------------------------------------------------
Real*8 Function Xi2(R)
!              correlation function
  use Pars
  Real*8, parameter :: twoPisq= 5.06606d-2, &
                       twopi =6.283185d0,   &
                       epsabs =1.d-6, epsrel=2.d-5
  Real*8 R,ANS,er,XupA,XupB,XupC,XupD,abserr
  !External Faux

  Rf = R
  XupD = 1000.*twopi/Rf
  !XupD =100.
  Call qawo ( Faux2, 6.28d-3, XupD, Rf, 2, epsabs, epsrel, ANS, abserr, &
             neval, ier )
     Xi2 = twoPisq*(ANS)/Rf
!     Xi = twoPisq*(ANSa)

    ! write(*,'(/1p,4g12.4,0p,3i7)')R,Xi2,ANS,abserr,neval,ier
   End Function Xi2
!----------------------------------------------------------
Real*8 Function XiTr(R)
!              correlation function
  use Pars
  Real*8, parameter :: twoPisq= 5.06606d-2, &
                       twopi =6.283185d0,   &
                       epsabs =1.d-6, epsrel=1.d-5
  Real*8 R,ANS,er,XupA,XupB,XupC,XupD,abserr
  !External Faux

  Rf = R
  XupD = 100.*twopi/Rf
  XupA = 10.*twopi/Rf
  !XupD =100.
  Call qawo ( Faux2, XupA, XupD, Rf, 2, epsabs, epsrel, ANS, abserr, &
             neval, ier )
     XiTr = twoPisq*(ANS)/Rf
!     Xi = twoPisq*(ANSa)

    ! write(*,'(/1p,4g12.4,0p,3i7)')R,Xi2,ANS,abserr,neval,ier
   End Function XiTr