// these routines call numerical recipes memory allocation routines.
// svd_invert(revcov,NB,C05_icov,&logdet) inverts the covariance
// matrix and calculates log of the determinant.
// Pk(k) calculates the power spectrum at k / hMpc^-1

// read in 2dFGRS data and set up covariance matrix
void setup_C05()
{
  // allocate file pointers for data files
  FILE *fp_data_2dF, *fp_win_2dF, *fp_cov_2dF;

  // memory to read in a line at a time
  const int bsz=120; char buf[bsz];

  // memory for model power spectrum
  double *pk_mod = dvector(1,NK);
  double *pk_mod_conv = dvector(1,NB);
  double sigsq = 0.0;
  int i,j;

  // allocate memory for 2dFGRS data
  C05_k_conv = dvector(1,NB);
  C05_k_orig = dvector(1,NK);
  C05_pk     = dvector(1,NB);
  C05_win    = dmatrix(1,NB,1,NK);
  C05_cov    = dmatrix(1,NB,1,NB);
  
  // read in power
  if((fp_data_2dF=fopen("2dF_C05_data.txt","r"))==NULL) 
    err_handler("cannot open data");
  i=1;
  while(fgets(buf,bsz,fp_data_2dF)) if(strncmp(buf,"#",1)!=0) {
    double tmp_err=0.0;
    if(sscanf(buf,"%*f %lf %lf\n",&C05_pk[i], &tmp_err)!=2)
      err_handler("data input 1");
    if(tmp_err>0.0) i++;
  }
  fclose(fp_data_2dF);

  // read in window
  if((fp_win_2dF=fopen("2dF_C05_win.txt","r"))==NULL) 
    err_handler("cannot open win");
  for(i=1;i<=NB;i++) fscanf(fp_win_2dF,"%lf",&C05_k_conv[i]);
  for(j=1;j<=NK;j++) {
    fscanf(fp_win_2dF,"%lf",&C05_k_orig[j]);
    for(i=1;i<=NB;i++) fscanf(fp_win_2dF,"%lf",&C05_win[i][j]);
  }
  fclose(fp_win_2dF);

  // read in covariance matrix
  if((fp_cov_2dF=fopen("2dF_C05_cov.txt","r"))==NULL) 
    err_handler("cannot open cov");
  for(i=1;i<=NB;i++) for(j=1;j<=NB;j++) fscanf(fp_cov_2dF,"%lf",&C05_cov[i][j]);
  fclose(fp_cov_2dF);

  // divide covariance matrix by LN power squared to get approximate
  // correlation matrix. We then multiply by the theoretical power
  // at each node to get back to an approximate covariance matrix
  // for that theory. We will need to include logdet in the -2lnL
  // calculation...

  // these are the parameters used for the LN catalogues
  h        = 0.7;             // hubble constant/100 kms^-1Mpc^-1
  n_s      = 1.0;             // scalar spectral index
  om_b     = 0.0408;          // baryon density
  om_c     = 0.1992;          // CDM density
  sig8     = 1.125;           // normalisation

  // normalise P(k) to sigma_8 value
  pow_norm = 1.;
  sigsq    = calc_sigmasq();
  pow_norm = (sig8*sig8)/sigsq;

  // set up model power to match LN catalogues
  for(j=1;j<=NK;j++) pk_mod[j] = Pk(C05_k_orig[j]);
  for(i=1;i<=NB;i++) {
    pk_mod_conv[i] = 0.0;
    for(j=1;j<=NK;j++) pk_mod_conv[i] += C05_win[i][j]*pk_mod[j];
  }

  // divide covariance matrix by convolved LN power
  for(i=1;i<=NB;i++) for(j=1;j<=NB;j++) 
    C05_cov[i][j]/=(pk_mod_conv[i]*pk_mod_conv[j]);

  free_dvector(pk_mod,1,NK);
  free_dvector(pk_mod_conv,1,NB);
}

// calculate likelihood of 2dFGRS data
double calc_chisq_C05()
{
  int i,j,k;

  // memory for arrays we need to calculate on the fly
  double *pk_orig   = dvector(1,NK);
  double *pk_conv   = dvector(1,NB); 
  double *diff      = dvector(1,NB);
  double **revcov   = dmatrix(1,NB,1,NB);
  double **C05_icov = dmatrix(1,NB,1,NB);
  double logdet=0.0, loglike=0.0; 

  // We assume an redshift space galaxy power spectrum of the form
  // (1.+Qk^2)/(1+Ak)*P(k)_{lin,mass}. A is fixed by the halo model +
  // semi-analytic simulations for redshift space to be A=1.4. Q is
  // also fixed to be Q=4.6. A & Q could instead be marginalised over
  // as nuisance parameters.

  // calculate model power
  for(j=1;j<=NK;j++) pk_orig[j] = Pk(C05_k_orig[j])*( (1.+4.6*C05_k_orig[j]*C05_k_orig[j])
						       / (1.+1.4*C05_k_orig[j]));

  // convolve model power with window
  for(i=1;i<=NB;i++) {
    pk_conv[i]=0.0;
    for(j=1;j<=NK;j++) pk_conv[i]+=C05_win[i][j]*pk_orig[j];
  }

  // set up model differences
  for(j=1;j<=NB;j++) diff[j] = C05_pk[j] - pk_conv[j];

  // renormalise covariance matrix
  for(i=1;i<=NB;i++) for(j=1;j<=NB;j++) 
    revcov[i][j]=C05_cov[i][j]*pk_conv[i]*pk_conv[j];

  // invert covariance matrix
  svd_invert(revcov,NB,C05_icov,&logdet);

  // log(likelihood) - chi^2 component
  loglike = 0.0;
  for(i=1;i<=NB;i++) for(j=i+1;j<=NB;j++) loglike -= diff[i]*C05_icov[i][j]*diff[j];
  for(i=1;i<=NB;i++) loglike -= 0.5*diff[i]*C05_icov[i][i]*diff[i];
    
  // log(likelihood) - determinant component 
  loglike -= 0.5*logdet;
    
  free_dvector(pk_orig,1,NK);
  free_dvector(pk_conv,1,NB);
  free_dvector(diff,1,NB);
  free_dmatrix(revcov,1,NB,1,NB);
  free_dmatrix(C05_icov,1,NB,1,NB);

  return -2.*loglike; // marginalised log(likelihood)
}