/* simple program to demonstrate how to convolve P(k) with the   */
/* Feb '01 2dFGRS window function both numerically, and using a  */
/* window matrix in k. Note that using this matrix intrinsically */
/* that P(k) is smooth on scales of 0.0045 hMpc^-1.              */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int main() {
  FILE *fp;
  float PkconvW(float), Pk(float);
  float PkmcW, k_conv[33], k_orig[101], win[33][101], pk_orig[101];
  int i,j;

  /* read in window matrix and relevant k-values */
  if((fp = fopen("2dFGRS_pk_win.txt","r"))==NULL) 
    { printf("cannot find window matrix"); exit(0); }
  for(i=1;i<=32;i++) fscanf(fp,"%f",&k_conv[i]);
  for(j=1;j<=100;j++) {
    fscanf(fp,"%f",&k_orig[j]);
    for(i=1;i<=32;i++) fscanf(fp,"%f",&win[i][j]);
  }

  /* output numerically calculated convolved P(k) */
  /* and convolved P(k) calculated using the window matrix */
  for(j=1;j<=100;j++) pk_orig[j] = Pk(k_orig[j]);
  for(i=1;i<=32;i++) {
    PkmcW = 0.0;
    for(j=1;j<=100;j++) PkmcW += win[i][j]*pk_orig[j];
    printf("k=%g, P(k) conv W(k) = %g, %g\n",k_conv[i],PkconvW(k_conv[i]),PkmcW);
  }

  return 1;
}

/* return normalised Power spectrum */
float Pk(float k)  
{  
  float gam,pow_n;
  float q, pn, tf, tfsq;
  float tk_eh(float);

  if(k==0.) return 0.;

  /* put in normalised power spectrum here */
  /* example - unnormalised BBKS */
  gam   = 0.2;
  pow_n = 1.0;
  q     = k/gam;
  tf    = log(1.+2.34*q)/(2.34*q);
  tfsq  = tf*tf/sqrt(1.0+q*(3.89+q*(259.21+q*(162.77+q*2027.17))));
  pn    = pow(k,pow_n)*tfsq;

  return pn;
}

/* perform 3D convolution in k */
float PkconvW(float rk)
{
  double s, sp, oa, a;
  int ik, imu, nint1=1, nint2=18;
  float dk, dmu, rkmax=0.3;
  float rmu, xk, xksq, rksq, yk;

  float Pk(float);
  float Wksq(float);

  dmu  = 2.0/(float)nint2;
  rksq = rk*rk; 

  s = 0.0;
 
  /* first data point in trapezoidal rule */
  dk=0.5*rkmax;
  for(ik=1;ik<=2;ik++) {
    if(ik==1) { xk=xksq=0.0; } else { xk=rkmax; xksq=xk*xk; }
    for(sp=0.0,imu=1;imu<=nint2;imu++) {
      rmu = -1.0+((float)imu-0.5)*dmu;
      yk  = sqrt(rksq+xksq-2.*xk*rk*rmu);
      sp += Pk(yk);
    }
    s += sp*xksq*Wksq(xk);
  }
  a = s*dk;

  /* loop through next iterations until convergence */
  do {
    oa = a; s = 0.0;
    dk = rkmax/(float)nint1;
    for(ik=1;ik<=nint1;ik++) {
      xk   = ((float)ik-0.5)*dk;
      xksq = xk*xk;
      for(sp=0.0,imu=1;imu<=nint2;imu++) {
	rmu = -1.0+((float)imu-0.5)*dmu;
	yk  = sqrt(rksq+xksq-2.*xk*rk*rmu);
	sp  += Pk(yk);
      }
      s += sp*xksq*Wksq(xk);
    }
    a = 0.5*(oa+s*dk);
    nint1 *= 2;
  } while( fabs(oa-a)>1.0e-4*a );

  return a*dmu;
}

/* return normalised window function (W_k)^2 */
float Wksq(float k)
{
  float ksq, kp1, kp2, norm;
  ksq  = k*k;
  kp1  = 1.17011e-05; 
  kp2  = 9.33106e-09; 
  norm = 1150292.2;

  return norm/(1.+(ksq/kp1)+(ksq*ksq/kp2));
}