/* lens_isw.c - plots for the iSW-review in Physics Reports */
/* by Bjoern Malte Schaefer, spirou@ita.uni-heidelberg.de */


// gcc-mp-4.9 -o kappa kappa.c -L/sw/lib -I/sw/include -lgsl -lgslcblas -lm

// add lensing bits


/* --- includes --- */
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_const_num.h>
#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_blas.h>


/* --- definitions --- */
#define pi M_PI
#define pi2 (pi * pi)
#define twopi (2.0 * M_PI)
#define fourpi (4.0 * M_PI)
#define fpisq (twopi * twopi)
#define eps 1.0e-4
#define epsabs eps
#define epsrel eps
#define NEVAL 1000

#define spirou_clight GSL_CONST_MKSA_SPEED_OF_LIGHT	
#define spirou_parsec GSL_CONST_MKSA_PARSEC	
#define spirou_kparsec (1.0e3 * spirou_parsec)
#define spirou_mparsec (1.0e6 * spirou_parsec)
#define spirou_gparsec (1.0e9 * spirou_parsec)
#define spirou_hubble (1.0e5 / spirou_mparsec)	
#define spirou_dhubble (spirou_clight / spirou_hubble / spirou_mparsec)	
#define spirou_thubble (1.0 / spirou_hubble / 3600.0 / 24.0 / 365.25 / 1e9)	
#define spirou_tcmb 2.726	
#define spirou_gnewton GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT	
#define spirou_hplanck GSL_CONST_MKSA_PLANCKS_CONSTANT_H	
#define spirou_kboltzmann GSL_CONST_MKSA_BOLTZMANN		
#define spirou_msolar GSL_CONST_MKSA_SOLAR_MASS	
#define spirou_qelectron GSL_CONST_MKSA_ELECTRON_CHARGE
#define spirou_melectron GSL_CONST_MKSA_MASS_ELECTRON

#define ASTEP 200
#define AMIN 0.001
#define AMAX 1.0
#define LMIN 10
#define LMAX 10000
#define CHISTEP 100
#define CHIMIN 10
#define CHIMAX 10000
#define ZMIN 0.01
#define ZMAX 10.0
#define ZSTEP 100

#define eps_derivative 0.05

#define spirou_rad2deg (360.0 / 2.0 / pi)
#define spirou_deg2rad (1.0 / spirou_rad2deg)
#define spirou_rad2min (spirou_rad2deg * 60.0)
#define spirou_min2rad (1.0 / spirou_rad2min)
#define spirou_rad2sec (spirou_rad2deg * 3600.0)
#define spirou_sec2rad (1.0 / spirou_rad2sec)
#define spirou_hubble_volume (spirou_dhubble * spirou_dhubble * spirou_dhubble)
#define spirou_rhocrit (3.0e11)
#define spirou_sigma2fwhm (2.0 * sqrt(2.0 * log(2.0)))
#define spirou_fwhm2sigma (1.0 / spirou_sigma2fwhm)

#define mode_s2n 0
#define mode_fisher 1
#define mode_chi2 2
#define mode_redshift 3


/* --- structures --- */
struct COSMOLOGY
{
	double omega_m;			// cosmological parameters
	double omega_k;
	double omega_q;
	double omega_b;
	double h;
	double w,wprime;
	double n_s;
	double sigma8;
	double rsmooth;
	double zmean,beta;		// observation
	
	gsl_interp_accel *acc_dplus;	// accelerators for growth
	gsl_spline *spline_dplus;
	gsl_interp_accel *acc_ddplus;
	gsl_spline *spline_ddplus;
	
	gsl_interp_accel *acc_a;	// a(chi) inversion
	gsl_spline *spline_a;
	
	double anorm;
	double fsky;
};

struct DATA
{
	struct COSMOLOGY *cosmology;
	
	double dp[ASTEP];					// growth function etc.
	
	double cpp[LMAX];					// auto-spectra
	
	double a[ASTEP],chi[ASTEP];				// a(chi) inversion

	double kappa_cmb[ZSTEP];
};


/* --- prototypes --- */
// cosmology
int INIT_COSMOLOGY(struct DATA *data);
int FIX_CDM_SPECTRUM();

// acceleration
int ALLOCATE_INTERPOLATION_MEMORY(struct DATA *data);
int FREE_INTERPOLATION_MEMORY(struct DATA *data);

// stepping
double chistep(int i);
double astep(int i);

// main functions
int COMPUTE_COMOVING_DISTANCE(struct DATA *data);
int COMPUTE_SPECTRUM(struct DATA *data);
int WRITE_SPECTRUM2DISK(struct DATA *data);

/* linear cdm spectrum */
double cdm_spectrum(double k);
double cdm_transfer(double k);
double sigma8();
double aux_dsigma8(double k,void *params);

/* cosmometry */
double a2z(double a);
double z2a(double z);
double a2com(double a);
double aux_dcom(double a,void *params);
double hubble(double a);
double dhubble(double a);
double w_de(double a);

/* cdm growth */
int d_plus_function(double t, const double y[], double f[],void *params);
double aux_d_plus(double a,double *result_d_plus,double *result_d_plus_prime);
double d_plus(double a);
double d_plus_prime(double a);
double x_plus(double a);
double aux_r(double a);
double aux_q(double a);

/* lensing spectrum */
double cpp_spectrum(double l);
double aux_cpp(double chi,void *params);
double p_weighting(double a,double chi);


/* --- global variables --- */;
struct COSMOLOGY *gcosmo;


/* --- main --- */
int main()
{
	struct DATA *data;
	
	printf("kappa: lensing for a fixed source redshift\n");
	printf("by Bjoern Malte Schaefer, spirou@ita.uni-heidelberg.de\n");
	
	data = (struct DATA *)calloc(1,sizeof(struct DATA));
	data->cosmology = (struct COSMOLOGY *)calloc(1,sizeof(struct COSMOLOGY));
	gcosmo = data->cosmology;
	
	INIT_COSMOLOGY(data);
	
	FIX_CDM_SPECTRUM();
	
	ALLOCATE_INTERPOLATION_MEMORY(data);
	
	COMPUTE_GROWTH(data);
	
	COMPUTE_COMOVING_DISTANCE(data);
	
	COMPUTE_SPECTRUM(data);
	
	WRITE_SPECTRUM2DISK(data);
	
	FREE_INTERPOLATION_MEMORY(data);
	
	free(data);
	
	exit(0);
}
/* --- end of main --- */


/* --- function allocate_interpolation_memory --- */
int ALLOCATE_INTERPOLATION_MEMORY(struct DATA *data)
{
	printf("lens_isw: allocating memory for interpolation functions....\n");

	gcosmo->acc_dplus = gsl_interp_accel_alloc();
	gcosmo->spline_dplus = gsl_spline_alloc(gsl_interp_cspline,ASTEP);
	gcosmo->acc_ddplus = gsl_interp_accel_alloc();
	gcosmo->spline_ddplus = gsl_spline_alloc(gsl_interp_cspline,ASTEP);
	gcosmo->acc_a = gsl_interp_accel_alloc();
	gcosmo->spline_a = gsl_spline_alloc(gsl_interp_cspline,ASTEP);
	
	return(0);
}
/* --- end of function allocate_interpolation_memory --- */


/* --- function free_interpolation_memory --- */
int FREE_INTERPOLATION_MEMORY(struct DATA *data)
{
	printf("lens_isw: deallocating memory for interpolation functions....\n");
	
	gsl_interp_accel_free(gcosmo->acc_dplus);
	gsl_spline_free(gcosmo->spline_dplus);
	gsl_interp_accel_free(gcosmo->acc_ddplus);
	gsl_spline_free(gcosmo->spline_ddplus);
	gsl_interp_accel_free(gcosmo->acc_a);
	gsl_spline_free(gcosmo->spline_a);
	
	return(0);
}
/* --- end of function free_interpolation memory --- */


/* --- function compute_comoving_distance [a(chi)-relation]--- */
int COMPUTE_COMOVING_DISTANCE(struct DATA *data)
{
	int i;
	double a,temp[ASTEP];
	
	printf("isw_bias: swapping (chi,a)-pairs....\n");
	
	// compute chi(a)
	for(i=0;i<ASTEP;i++)
	{
		data->a[i] = a = astep(i);
		data->chi[i] = a2com(a);
	}
	
	// swap chi
	for(i=0;i<ASTEP;i++)
	{
		temp[i] = data->chi[i];
	}
	
	for(i=0;i<ASTEP;i++)
	{
		data->chi[i] = temp[ASTEP-i-1];
	}
	
	// swap a
	for(i=0;i<ASTEP;i++)
	{
		temp[i] = data->a[i];
	}
	
	for(i=0;i<ASTEP;i++)
	{
		data->a[i] = temp[ASTEP-i-1];
	}
	
	gsl_spline_init(gcosmo->spline_a,data->chi,data->a,ASTEP);
	
	return(0);
}
/* --- end of function compute_comoving distance --- */


/* --- function compute_growth --- */
int COMPUTE_GROWTH(struct DATA *data)
{
	int i;
	double aux,a;
	
	printf("lens_isw: computing growth function and derivatives....\n");
	
	for(i=0;i<ASTEP;i++)
	{
		data->a[i] = a = astep(i);
		
		data->dp[i] = d_plus(a);
	}

	gsl_spline_init(gcosmo->spline_dplus,data->a,data->dp,ASTEP);
	
	return(0);
}
/* --- end of function compute_growth --- */


/* ---  --- */
int COMPUTE_SPECTRUM(struct DATA *data)
{
	int l;
	
	printf("lens_isw: computing all auto- and cross-spectra....\n");
	
	// compute spectra
	
#pragma omp parallel for
	for(l=LMIN;l<LMAX;l++)
	{
		data->cpp[l] = cpp_spectrum(l);
	}
	
	return(0);
}
/* ---  --- */


// ---  ---
int WRITE_SPECTRUM2DISK(struct DATA *data)
{
	int l;
	FILE *file;
	
	printf("lens_isw: saving all auto- and cross-spectra to disk....\n");
	
	file = fopen("/Users/spirou/Seafile/data/projects/flrw_cone/spectra/data/spectrum.data","w+");
	for(l=LMIN;l<LMAX;l++)
	{
		fprintf(file,"%i  %e\n",l,data->cpp[l]);
	}
	fclose(file);
	
	return(0);
}
// ---  ---


/* --- function astep [linear equidistant stepping in scale factor a] --- */
double astep(int i)
{
	double result,delta;
	
	delta = (AMAX - AMIN) / ((double)(ASTEP-1));
	result = AMIN + delta * i;
	
	return(result);
}
/* --- end of function astep --- */


/* --- function init_cosmology [fiducial values for cosmological parameters] --- */
int INIT_COSMOLOGY(struct DATA *data)
{
	printf("lens_isw: initialising cosmological model....\n");
	
	data->cosmology->omega_m = 0.25;		// matter density
	data->cosmology->omega_k = 0.0;			// curvature - spatial flatness!
	data->cosmology->omega_q = 1.0 - data->cosmology->omega_m;
	data->cosmology->omega_b = 0.04;		// baryonic density, for shape parameter
	data->cosmology->sigma8 = 0.8;			// fluctuation amplitude
	data->cosmology->n_s = 1.0;			// spectral index
	data->cosmology->h = 0.72;			// hubble parameter, for shape parameter
	
	data->cosmology->w = -1.0;			// de eos parameter
	data->cosmology->wprime = 0.0;			// time evolution of de eos
	data->cosmology->rsmooth = 8.0;			// sigma8-smoothing
	
	data->cosmology->zmean = 0.64;			
	data->cosmology->beta = 1.5;
	data->cosmology->fsky = 0.5;
	
	return(0);
}
/* --- end of function init_cosmology --- */


/* --- function fix_cdm_spectrum [sigma8-normalisation] --- */
int FIX_CDM_SPECTRUM()
{
	printf("lens_isw: fixing cdm spectrum normalisation....\n");
	
	gcosmo->anorm = 1.0;
	gcosmo->anorm = gcosmo->sigma8 / sigma8();
	
	return(0);
}
/* --- end of function fix_cdm_spectrum --- */


/* --- function sigma8 --- */
double sigma8()
{
	double result,aux,error;
	gsl_integration_workspace *w = gsl_integration_workspace_alloc(NEVAL);
	gsl_function F;
	
	F.function = &aux_dsigma8;
	F.params = NULL;
	
	gsl_integration_qagiu(&F,0.0,eps,eps,NEVAL,w,&aux,&error);
	result = sqrt(aux / 2.0 / pi2);
	
	gsl_integration_workspace_free(w);
	
	return(result);
}
/* --- end of function sigma8 --- */


/* --- function aux_dsigma8 --- */
double aux_dsigma8(double k,void *params)
{
	double result,window,aux;
	
	aux = gcosmo->rsmooth * k;
	window = 3.0 * gsl_sf_bessel_j1(aux) / aux;
	result = cdm_spectrum(k) * gsl_pow_int(k * window,2);
	
	return(result);
}
/* --- end of function aux_dsigma8 --- */


/* --- function cdm_spectrum --- */
double cdm_spectrum(double k)
{
	double result;
	
	result = pow(k,gcosmo->n_s) * gsl_pow_int(gcosmo->anorm * cdm_transfer(k),2);
	
	return(result);
}
/* --- end of function cdm_spectrum --- */


/* --- function cdm_transfer [Bardeen - CDM transfer function] --- */
double cdm_transfer(double k)
{
	double result,q,shape,fb,aux;
	double alpha = 2.34, beta = 3.89, gamma = 16.1, delta = 5.46, epsilon = 6.71;

	fb = gcosmo->omega_b / gcosmo->omega_m;
	shape = gcosmo->omega_m * gcosmo->h * exp(-gcosmo->omega_b - sqrt(2.0 * gcosmo->h) * fb);
	q = k / shape;
	
	aux = 1.0 + beta * q + gsl_pow_int(gamma * q,2) + gsl_pow_int(delta * q,3) + gsl_pow_int(epsilon * q,4);
	result = log(1.0 + alpha * q) / (alpha * q) * pow(aux,-1.0/4.0);
	
	return(result);
}
/* --- end of function cdm_transfer --- */


/* --- function a2z --- */
double a2z(double a)
{
	double result;
	
	result = 1.0 / a - 1.0;
	
	return(result);
}
/* --- end of function a2z --- */


/* --- function z2a --- */
double z2a(double z)
{
	double result;
	
	result = 1.0 / (1.0 + z);
	
	return(result);
}
/* --- end of function z2a --- */


/* --- function a2com [comoving distance in Mpc/h] --- */
double a2com(double a)
{
	double result,error;
	gsl_integration_workspace *w = gsl_integration_workspace_alloc(NEVAL);
	gsl_function F;
	
	F.function = &aux_dcom;
	F.params = NULL;
	
	gsl_integration_qag(&F,1.0,a,epsabs,epsrel,NEVAL,GSL_INTEG_GAUSS61,w,&result,&error);
	
	gsl_integration_workspace_free(w);
	
	return(result);
}
/* --- end of function a2com --- */


/* --- function aux_dcom [auxiliary function for a2com] --- */
double aux_dcom(double a,void *params)
{
	double result;
	
	result = - spirou_dhubble / gsl_pow_int(a,2) / hubble(a);
	
	return(result);
}
/* --- end of function aux_dcom --- */


/* --- function hubble [parameterisation of dark energy eos: w(a) = w0 + (1-a)w'] --- */
double hubble(double a)
{
	double result;
	double aux;
	
	/* matter and curvature */
	result = gcosmo->omega_m / gsl_pow_int(a,3);
	result += gcosmo->omega_k / gsl_pow_int(a,2);
	
	/* dark energy , parameterised with eos w(a) = w0 + (1-a) * w' */
	aux = -(1.0 + gcosmo->w + gcosmo->wprime) * log(a) + gcosmo->wprime * (a - 1.0);
	result += gcosmo->omega_q * exp(3.0 * aux);
	
	result = sqrt(result);
	
	return(result);
}
/* --- end of function hubble --- */


/* --- function dhubble [derivative dlnH/dlna of the Hubble function, only constant w] --- */
double dhubble(double a)
{
	double h,aux,result;
	
	h = hubble(a);
	aux = gcosmo->omega_m / gsl_pow_int(a,3) + 
		2.0 / 3.0 * gcosmo->omega_k / gsl_pow_int(a,2) +
		(1.0 + gcosmo->w) * gcosmo->omega_q / gsl_pow_int(a,3.0 * (1.0 + gcosmo->w));
	result = -3.0 / 2.0 * aux / gsl_pow_int(h,2);
	
	return(result);
}
/* --- end of function dhubble --- */


/* --- function w [dark energy eos parameterisation] --- */
double w_de(double a)
{
	double result;
	
	result = gcosmo->w + (1.0 - a) * gcosmo->wprime;
	
	return(result);
}
/* --- end of function w --- */


/* --- function p_galaxy [redshift distribution] --- */
double p_galaxy(double z)
{
	double norm,aux,result;
	
	aux = z / gcosmo->zmean;
	norm = gcosmo->zmean / gcosmo->beta * gsl_sf_gamma(3.0 / gcosmo->beta);
	
	// Bartelmann + Schneider
	result = gsl_pow_int(aux,2) * exp(-pow(aux,gcosmo->beta)) / norm;
		
	return(result);
}
/* --- end of function p_galaxy --- */


/* --- function omega_m [matter density parameter as a function of a] --- */
double omega_m(double a)
{
	double h,result;
	
	h = hubble(a);
	result = gcosmo->omega_m / gsl_pow_int(h,2) / gsl_pow_int(a,3);
	
	return(result);
}
/* --- end of function omega_m --- */


/* --- d_plus_function [defines f0 = dy1/da and f1 = dy2/da] --- */
int d_plus_function(double t, const double y[], double f[],void *params)
{
	/* derivatives f_i = dy_i/dt */
	f[0] = y[1];
	f[1] = aux_r(t) * y[0] - aux_q(t) * y[1];
	
	return(GSL_SUCCESS);
}
/* --- end of function dplus_function --- */


/* --- function dplus --- */
double aux_d_plus(double a,double *result_d_plus,double *result_d_plus_prime)
{
	double result;
	int status;
	const gsl_odeiv_step_type *T = gsl_odeiv_step_rkf45;
	gsl_odeiv_step *s = gsl_odeiv_step_alloc(T,2);
	gsl_odeiv_control *c = gsl_odeiv_control_y_new(0.0,epsrel);
	gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(2);
	double t=1e-6,h = 1e-4;			/* t = initial scale factor, h = absolute accuracy */
	double y[2] = {1.0,0.0};		/* initial conditions, dy1(0)/da = 1, dy2(0)/da=0 */
  
	/* result from solution of a 2nd order differential equation, transformed to a system of 2 1st order deqs */
	gsl_odeiv_system sys = {d_plus_function,NULL,2,NULL};
	
	while(t < a)
	{
		status = gsl_odeiv_evolve_apply(e,c,s,&sys,&t,a,&h,y);
		if(status != GSL_SUCCESS) break;
	}
	
	gsl_odeiv_evolve_free(e);
	gsl_odeiv_control_free(c);
	gsl_odeiv_step_free(s);
	
	result = *result_d_plus = y[0];		/* d_plus */
	*result_d_plus_prime = y[1];		/* d(d_plus)/da */
	
	return(result);
}
/* --- end of function dplus --- */


/* --- function d_plus --- */
double d_plus(double a)
{
	double dummy,norm,result;
	
	aux_d_plus(a,&result,&dummy);
	aux_d_plus(1.0,&norm,&dummy);
	
	result /= norm;
	
	return(result);
}
/* --- end of function d_plus --- */


/* --- function d_plus_prime = d(d_plus)/da--- */
double d_plus_prime(double a)
{
	double dummy,norm,result;
	
	aux_d_plus(a,&dummy,&result);
	aux_d_plus(1.0,&norm,&dummy);
	
	result /= norm;
	
	return(result);
}
/* --- end of function d_plus_prime --- */


/* --- function x_plus [auxiliary function, see Linder+Jenkins MNRAS 346, 573-583 (2003) for definition] --- */
double x_plus(double a)
{
	double result,aux;

	aux = 3.0 * gcosmo->wprime * (1.0 - a);
	result = gcosmo->omega_m / (1.0 - gcosmo->omega_m) * pow(a,3.0 * (gcosmo->w + gcosmo->wprime)) * exp(aux);
	
	return(result);
}
/* --- end of function x_plus --- */


/* --- function aux_r --- */
double aux_r(double a)
{
	double aux,result;
	
	aux = x_plus(a);
	result = 3.0 / 2.0 * aux / (1.0 + aux) / gsl_pow_int(a,2);
	
	return(result);
}
/* --- end of function aux_r --- */


/* --- function aux_q --- */
double aux_q(double a)
{
	double result;
	
	result = 3.0 / 2.0 * (1.0 - w_de(a) / (1.0 + x_plus(a))) / a;
	
	return(result);
}
/* --- end of function aux_q --- */


/* --- function cpp_spectrum --- */
double cpp_spectrum(double l)
{
	double result,aux,error;
	gsl_integration_workspace *w = gsl_integration_workspace_alloc(NEVAL);
	gsl_function F;
	
	F.function = &aux_cpp;
	F.params = &l;
	
	gsl_integration_qag(&F,CHIMIN,CHIMAX,eps,eps,NEVAL,GSL_INTEG_GAUSS61,w,&aux,&error);
	result = aux;
	
	gsl_integration_workspace_free(w);
	
	return(result);
}
/* --- end of function cpp_spectrum --- */


/* --- function aux_cpp [integrand for cpp_spectrum] --- */
double aux_cpp(double chi,void *params)
{
	double a,aux,k,result;
	double l = *(double *)params;
	
	a = gsl_spline_eval(gcosmo->spline_a,chi,gcosmo->acc_a);
	k = l / chi;
	
	aux = cdm_spectrum(k);
	result = gsl_pow_int(p_weighting(a,chi),2) * aux / gsl_pow_int(chi,2);
	
	return(result);
}
/* --- end of function aux_cgg --- */


/* --- function p_weighting [weighting function for the CMB-lensing potential] --- */
double p_weighting(double a,double chi)
{
	double dp,aux,result;

	dp = gsl_spline_eval(gcosmo->spline_dplus,a,gcosmo->acc_dplus);
	aux = (CHIMAX - chi) / CHIMAX;
	
	result = 1.5 * gcosmo->omega_m / gsl_pow_int(spirou_dhubble,2) * aux * dp / a * chi;
	
	return(result);
}
/* --- end of function p_weighting --- */