*** Fisher matrix calculation ***
August 30, 2008: E.Komatsu

Here we provide a simple program for computing the Fisher matrix
of the angular diameter distance, Da, and the Hubble expansion rate, H,
for an idealized cosmic-variance limited galaxy survey (i.e., the number
density of galaxies is infinite) with a volume of 1 h^-3 Gpc^3.
For other volumes the errors may be scaled simply by sqrt(volue/[1 h^-3 Gpc^3]).

The assumption that the survey is limited only by cosmic variance simplifies
the program greatly, as the result becomes independent of the absolute value
of power spectrum, P(k), the number density, or redshifts.

The only necessary inputs are:

(1) The logarithmic derivative of P(k), i.e., n_eff(k)=dln[P(k)]/dlnk,
(2) The linear redshift distortion parameter, beta, and
(3) The maximum wavenumber included in the analysis, k_max.

We provide the data for (1), "wmap5baosn_max_likelihood_dlnpkdlnk.dat,"
computed from "wmap5baosn_max_likelihood_matterpower.dat," which was
generated using CAMB code for the maximum likelihood parameters
given in Table I of Komatsu et al.(2008) [WMAP 5-year interpretation paper] 
with "WMAP5+BAO+SN". The input file for CAMB is also provided 
(wmap5baosn_max_likelihood_params.ini). You may compute dlnPk/dlnk by yourself
using another program supplied in this directory, "compute_dlnpkdlnk.f90".

The program reports on the errors in ln(Da) and ln(H) in percent
as well as the correlation coefficient between them, with some other parameters 
marginalized over: ln(amplitude), beta, and/or ns.

- To compile and use the sample programs (given below), edit Makefile
and simply "./make"
- It will generate executables called "compute_dlnpkdlnk" and "fisher_distance_cv"

=========================================================================
A simple program to generate dlnPk/dlnk is given as "compute_dlnpkdlnk.f90":

PROGRAM Compute_dlnPkdlnk
 USE linearpk
 ! A program for computing the logarithmic derivative of
 ! the linear power spectrum, dln(P(k))/dlnk
 ! August 30, 2008: E.Komatsu
 double precision :: dlnpkdlnk,k_ov_h,dummy
 character(len=128) :: infile,outfile
 integer :: i,n
 external dlnpkdlnk
 ! input power spectrum data
 infile='wmap5baosn_max_likelihood_matterpower.dat'
 n=896 ! # of lines in the file
 CALL open_linearpk(infile,n)
 ! output 
 outfile='dlnpkdlnk.txt'
 open(1,file=infile,status='old')
 open(2,file=outfile,status='unknown')
 read(1,*)k_ov_h,dummy
 do i=1,n-2
    read(1,*)k_ov_h,dummy
    write(2,'(2E16.5)')k_ov_h,dlnpkdlnk(k_ov_h)
 enddo 
 close(1)
 close(2)
 CALL close_linearpk
END PROGRAM Compute_dlnPkdlnk

The Fisher matrix code ("fisher_distance_cv.f90") is

PROGRAM Fisher_Distance_CV
  ! A sample program for computing the errors on ln(Da) and ln(H) for
  ! a cosmic-variance limited survey with 1 h^-3 Gpc^3 volume, with
  ! ln(amplitude), beta, and/or ns marginalized over.
  ! August 30, 2008: E.Komatsu
  IMPLICIT none
  double precision :: kmax_ov_h,beta
  integer :: npara=5 ! # of parameters
  double precision, dimension(:), allocatable :: der
  double precision, dimension(:,:), allocatable :: cov,fis
  integer, dimension(:), allocatable :: work,ind
  double precision :: k_ov_h,neff,k0_ov_h,dlnk,mu,mu2,dmu,factor,wkmu
  character(len=128) :: filename
  integer :: i,j
  print*,'k_max in units of h Mpc^-1?'
  read*,kmax_ov_h
  print*,'beta?'
  read*,beta
  filename='wmap5baosn_max_likelihood_dlnpkdlnk.dat'
  open(1,file=filename,status='old')
  read(1,*)k0_ov_h,neff
  ALLOCATE(der(npara),fis(npara,npara))
  fis=0d0
  do while (k_ov_h<=kmax_ov_h)
     read(1,*)k_ov_h,neff
     dlnk=dlog(k_ov_h)-dlog(k0_ov_h)
     factor=(k_ov_h)**3d0*dlnk/(2d0*3.1415926535d0**2d0) ! h^3 Mpc^-3
     mu=0d0
     dmu=1d-3
     do while (mu<=1d0)
        wkmu=0.5d0*1d9 ! for a survey volume of 1 h^-3 Gpc^3
        mu2=mu*mu
        der(1)=neff*(1d0-mu2)-4d0*beta*mu2*(1d0-mu2)/(1d0+beta*mu2) ! dlnPkdlnDa
        der(2)=neff*(   -mu2)-4d0*beta*mu2*(1d0-mu2)/(1d0+beta*mu2) ! dlnPkdlnH
        der(3)=1d0                    ! dlnPkdln(amplitude)
        der(4)=2d0*mu2/(1d0+beta*mu2) ! dlnPkd(beta)
        der(5)=dlog(k_ov_h)           ! dlnPkd(ns)
        do i=1,npara
           fis(i,:)=fis(i,:)+factor*wkmu*der(i)*der(:)*dmu
        enddo
        mu=mu+dmu
     enddo
     k0_ov_h=k_ov_h
  enddo
  close(1)
  DEALLOCATE(der)
  print*,'k_max=',kmax_ov_h,' h Mpc^-1'
  print*,'beta=',beta
  print*,'*** Errors are cosmic variance only, and reported for a volume of 1 h^-3 Gpc^3. ***'
  print*,'*** Scale the result by sqrt(volume/[1 h^-3 Gpc^3]) for other volumes. ***'
  ! lnDa, lnH
  print*,'no marginalization over the other parameters'
  ALLOCATE(cov(2,2),work(2),ind(2))
  ind=(/1,2/)
  cov=fis(ind,ind)
  CALL DVERT(cov,2,2,work)
  CALL output(cov)
  DEALLOCATE(cov,work,ind)
  ! lnDa, lnH, ln(amplitude)
  print*,'marginalized over ln(amplitude)'
  ALLOCATE(cov(3,3),work(3),ind(3))
  ind=(/1,2,3/)
  cov=fis(ind,ind)
  CALL DVERT(cov,3,3,work)
  CALL output(cov(1:2,1:2))
  DEALLOCATE(cov,work,ind)
  ! lnDa, lnH, beta
  print*,'marginalized over beta'
  ALLOCATE(cov(3,3),work(3),ind(3))
  ind=(/1,2,4/)
  cov=fis(ind,ind)
  CALL DVERT(cov,3,3,work)
  CALL output(cov(1:2,1:2))
  DEALLOCATE(cov,work,ind)
  ! lnDa, lnH, ln(amplitude), beta
  print*,'marginalized over ln(amplitude) and beta'
  ALLOCATE(cov(4,4),work(4),ind(4))
  ind=(/1,2,3,4/)
  cov=fis(ind,ind)
  CALL DVERT(cov,4,4,work)
  CALL output(cov(1:2,1:2))
  DEALLOCATE(cov,work,ind)
  ! lnDa, lnH, ns
  print*,'marginalized over ns'
  ALLOCATE(cov(3,3),work(3),ind(3))
  ind=(/1,2,5/)
  cov=fis(ind,ind)
  CALL DVERT(cov,3,3,work)
  CALL output(cov(1:2,1:2))
  DEALLOCATE(cov,work,ind)
  ! lnDa, lnH, ln(amplitude), ns
  print*,'marginalized over ln(amplitude) and ns'
  ALLOCATE(cov(4,4),work(4),ind(4))
  ind=(/1,2,3,5/)
  cov=fis(ind,ind)
  CALL DVERT(cov,4,4,work)
  CALL output(cov(1:2,1:2))
  DEALLOCATE(cov,work,ind)
  ! lnDa, lnH, ln(amplitude), beta, ns
  print*,'marginalized over ln(amplitude), beta and ns'
  ALLOCATE(cov(5,5),work(5),ind(5))
  ind=(/1,2,3,4,5/)
  cov=fis(ind,ind)
  CALL DVERT(cov,5,5,work)
  CALL output(cov(1:2,1:2))
  DEALLOCATE(cov,work,ind)
!
  DEALLOCATE(fis)
END PROGRAM Fisher_Distance_CV
SUBROUTINE output(c)
  IMPLICIT none
  double precision, intent(IN) :: c(2,2)
  double precision :: r12,err_lnda,err_lnh
  r12=c(1,2)/sqrt(c(1,1)*c(2,2))
  err_lnda=sqrt(c(1,1))
  err_lnh=sqrt(c(2,2))
  print'(1A15,1F9.5)','Err[lnDa](%) =',err_lnda*1d2
  print'(1A15,1F9.5)','Err[lnH](%)  =',err_lnh*1d2
  print'(1A15,1F9.5)','r(lnDa,lnH)  =',r12
  return
END SUBROUTINE output