text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#include "ccl.h"

#ifdef HAVE_ANGPOW
#include "Angpow/angpow_ccl.h"
#endif

#define CCL_FRAC_RELEVANT 5E-4
//#define CCL_FRAC_RELEVANT 1E-3
//Gets the x-interval where the values of y are relevant
//(meaning, that the values of y for those x are at least above a fraction frac of its maximum)
static void get_support_interval(int n,double *x,double *y,double frac,
				 double *xmin_out,double *xmax_out)
{
  int ix;
  double ythr=-1000;

  //Initialize as the original edges in case we don't find an interval
  *xmin_out=x[0];
  *xmax_out=x[n-1];

  //Find threshold
  for(ix=0;ix<n;ix++) {
    if(y[ix]>ythr) ythr=y[ix];
  }
  ythr*=frac;

  //Find minimum
  for(ix=0;ix<n;ix++) {
    if(y[ix]>=ythr) {
      *xmin_out=x[ix];
      break;
    }
  }

  //Find maximum
  for(ix=n-1;ix>=0;ix--) {
    if(y[ix]>=ythr) {
      *xmax_out=x[ix];
      break;
    }
  }
}

//Wrapper around spline_eval with GSL function syntax
static double speval_bis(double x,void *params)
{
  return ccl_spline_eval(x,(SplPar *)params);
}

void ccl_cl_workspace_free(CCL_ClWorkspace *w)
{
  free(w->l_arr);
  free(w);
}

CCL_ClWorkspace *ccl_cl_workspace_default(int lmax,int l_limber,int non_limber_method,
					  double l_logstep,int l_linstep,
					  double dchi,double dlk,double zmin,int *status)
{
  CCL_ClWorkspace *w=(CCL_ClWorkspace *)malloc(sizeof(CCL_ClWorkspace));
  if(w==NULL) {
    *status=CCL_ERROR_MEMORY;
    //Can't access cosmology object
    //    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_workspace_default(); memory allocation\n");
    return NULL;
  }

  //Set params
  w->dchi=dchi;
  w->dlk =dlk ;
  if(zmin<=0) { //We should make sure that zmin is always strictly positive
    free(w);
    *status=CCL_ERROR_INCONSISTENT;
    return NULL;
  }
  w->zmin=zmin;
  w->lmax=lmax;
  if((non_limber_method!=CCL_NONLIMBER_METHOD_NATIVE) && (non_limber_method!=CCL_NONLIMBER_METHOD_ANGPOW)) {
    free(w);
    *status=CCL_ERROR_INCONSISTENT;
    //Can't access cosmology object
    //    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_workspace_default(); unknown non-limber method\n");
    return NULL;
  }
  w->nlimb_method=non_limber_method;
  w->l_limber=l_limber;
  w->l_logstep=l_logstep;
  w->l_linstep=l_linstep;

  //Compute number of multipoles
  int i_l=0,l0=0;
  int increment=CCL_MAX(((int)(l0*(w->l_logstep-1.))),1);
  while((l0 < w->lmax) && (increment < w->l_linstep)) {
    i_l++;
    l0+=increment;
    increment=CCL_MAX(((int)(l0*(w->l_logstep-1))),1);
  }
  increment=w->l_linstep;
  while(l0 < w->lmax) {
    i_l++;
    l0+=increment;
  }

  //Allocate array of multipoles
  w->n_ls=i_l+1;
  w->l_arr=(int *)malloc(w->n_ls*sizeof(int));
  if(w->l_arr==NULL) {
    free(w);
    *status=CCL_ERROR_MEMORY;
    //Can't access cosmology object
    //    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_workspace_default(); memory allocation\n");
    return NULL;
  }

  //Redo the computation above and store values of ell
  i_l=0; l0=0;
  increment=CCL_MAX(((int)(l0*(w->l_logstep-1.))),1);
  while((l0 < w->lmax) && (increment < w->l_linstep)) {
    w->l_arr[i_l]=l0;
    i_l++;
    l0+=increment;
    increment=CCL_MAX(((int)(l0*(w->l_logstep-1))),1);
  }
  increment=w->l_linstep;
  while(l0 < w->lmax) {
    w->l_arr[i_l]=l0;
    i_l++;
    l0+=increment;
  }
  //Don't go further than lmaw
  w->l_arr[w->n_ls-1]=w->lmax;

  return w;
}

CCL_ClWorkspace *ccl_cl_workspace_default_limber(int lmax,double l_logstep,int l_linstep,
						 double dlk,int *status)
{
  return ccl_cl_workspace_default(lmax,-1,CCL_NONLIMBER_METHOD_NATIVE,l_logstep,l_linstep,3.,dlk,0.05,status);
}

//Params for lensing kernel integrand
typedef struct {
  double chi;
  SplPar *spl_pz;
  ccl_cosmology *cosmo;
  int *status;
} IntLensPar;

//Integrand for lensing kernel
static double integrand_wl(double chip,void *params)
{
  IntLensPar *p=(IntLensPar *)params;
  double chi=p->chi;
  double a=ccl_scale_factor_of_chi(p->cosmo,chip, p->status);
  double z=1./a-1;
  double pz=ccl_spline_eval(z,p->spl_pz);
  double h=p->cosmo->params.h*ccl_h_over_h0(p->cosmo,a, p->status)/CLIGHT_HMPC;

  if(chi==0)
    return h*pz;
  else
    return h*pz*ccl_sinn(p->cosmo,chip-chi,p->status)/ccl_sinn(p->cosmo,chip,p->status);
}

//Integral to compute lensing window function
//chi     -> comoving distance
//cosmo   -> ccl_cosmology object
//spl_pz  -> normalized N(z) spline
//chi_max -> maximum comoving distance to which the integral is computed
//win     -> result is stored here
static int window_lensing(double chi,ccl_cosmology *cosmo,SplPar *spl_pz,double chi_max,double *win)
{
  int gslstatus =0, status =0;
  double result,eresult;
  IntLensPar ip;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);

  ip.chi=chi;
  ip.cosmo=cosmo;
  ip.spl_pz=spl_pz;
  ip.status = &status;
  F.function=&integrand_wl;
  F.params=&ip;
  gslstatus=gsl_integration_qag(&F, chi, chi_max, 0,
                                ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                w, &result, &eresult);
  *win=result;
  gsl_integration_workspace_free(w);
  if(gslstatus!=GSL_SUCCESS || *ip.status) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: window_lensing():");
    return 1;
  }
  //TODO: chi_max should be changed to chi_horizon
  //we should precompute this quantity and store it in cosmo by default

  return 0;
}

//Params for lensing kernel integrand
typedef struct {
  double chi;
  SplPar *spl_pz;
  SplPar *spl_sz;
  ccl_cosmology *cosmo;
  int *status;
} IntMagPar;

//Integrand for magnification kernel
static double integrand_mag(double chip,void *params)
{
  IntMagPar *p=(IntMagPar *)params;
//EK: added local status here as the status testing is done in routines called from this function
  double chi=p->chi;
  double a=ccl_scale_factor_of_chi(p->cosmo,chip, p->status);
  double z=1./a-1;
  double pz=ccl_spline_eval(z,p->spl_pz);
  double sz=ccl_spline_eval(z,p->spl_sz);
  double h=p->cosmo->params.h*ccl_h_over_h0(p->cosmo,a, p->status)/CLIGHT_HMPC;

  if(chi==0)
    return h*pz*(1-2.5*sz);
  else
    return h*pz*(1-2.5*sz)*ccl_sinn(p->cosmo,chip-chi,p->status)/ccl_sinn(p->cosmo,chip,p->status);
}

//Integral to compute magnification window function
//chi     -> comoving distance
//cosmo   -> ccl_cosmology object
//spl_pz  -> normalized N(z) spline
//spl_pz  -> magnification bias s(z)
//chi_max -> maximum comoving distance to which the integral is computed
//win     -> result is stored here
static int window_magnification(double chi,ccl_cosmology *cosmo,SplPar *spl_pz,SplPar *spl_sz,
				double chi_max,double *win)
{
  int gslstatus =0, status =0;
  double result,eresult;
  IntMagPar ip;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);

  ip.chi=chi;
  ip.cosmo=cosmo;
  ip.spl_pz=spl_pz;
  ip.spl_sz=spl_sz;
  ip.status = &status;
  F.function=&integrand_mag;
  F.params=&ip;
  gslstatus=gsl_integration_qag(&F, chi, chi_max, 0,
                                ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                w, &result, &eresult);
  *win=result;
  gsl_integration_workspace_free(w);
  if(gslstatus!=GSL_SUCCESS || *ip.status) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: window_magnification():");
    return 1;
  }
  //TODO: chi_max should be changed to chi_horizon
  //we should precompute this quantity and store it in cosmo by default

  return 0;
}

//CCL_ClTracer creator
//cosmo   -> ccl_cosmology object
//tracer_type -> type of tracer. Supported: CL_TRACER_NC, CL_TRACER_WL
//nz_n -> number of points for N(z)
//z_n  -> array of z-values for N(z)
//n    -> corresponding N(z)-values. Normalization is irrelevant
//        N(z) will be set to zero outside the range covered by z_n
//nz_b -> number of points for b(z)
//z_b  -> array of z-values for b(z)
//b    -> corresponding b(z)-values.
//        b(z) will be assumed constant outside the range covered by z_n
static CCL_ClTracer *cl_tracer(ccl_cosmology *cosmo,int tracer_type,
				   int has_rsd,int has_magnification,int has_intrinsic_alignment,
				   int nz_n,double *z_n,double *n,
				   int nz_b,double *z_b,double *b,
				   int nz_s,double *z_s,double *s,
				   int nz_ba,double *z_ba,double *ba,
				   int nz_rf,double *z_rf,double *rf,
				   double z_source, int * status)
{
  int clstatus=0, gslstatus;
  CCL_ClTracer *clt=(CCL_ClTracer *)malloc(sizeof(CCL_ClTracer));
  if(clt==NULL) {

    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): memory allocation\n");
    return NULL;
  }

  if ( ((cosmo->params.N_nu_mass)>0) && tracer_type==CL_TRACER_NC && has_rsd){
	  free(clt);
	  *status=CCL_ERROR_NOT_IMPLEMENTED;
	  ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer_new(): Number counts tracers with rsd not yet implemented in cosmologies with massive neutrinos.");
	  return NULL;
  }

  clt->tracer_type=tracer_type;
  clt->computed_transfer=0;

  double hub=cosmo->params.h*ccl_h_over_h0(cosmo,1.,status)/CLIGHT_HMPC;
  clt->prefac_lensing=1.5*hub*hub*cosmo->params.Omega_m;

  if((tracer_type==CL_TRACER_NC)||(tracer_type==CL_TRACER_WL)) {
    get_support_interval(nz_n,z_n,n,CCL_FRAC_RELEVANT,&(clt->zmin),&(clt->zmax));
    clt->chimax=ccl_comoving_radial_distance(cosmo,1./(1+clt->zmax),status);
    clt->chimin=ccl_comoving_radial_distance(cosmo,1./(1+clt->zmin),status);
    clt->spl_nz=ccl_spline_init(nz_n,z_n,n,0,0);
    if(clt->spl_nz==NULL) {
      free(clt);
      *status=CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): error initializing spline for N(z)\n");
      return NULL;
    }

    //Normalize n(z)
    gsl_function F;
    double nz_norm,nz_enorm;
    double *nz_normalized=(double *)malloc(nz_n*sizeof(double));
    if(nz_normalized==NULL) {
      ccl_spline_free(clt->spl_nz);
      free(clt);
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): memory allocation\n");
      return NULL;
    }

    gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);
    F.function=&speval_bis;
    F.params=clt->spl_nz;
    gslstatus=gsl_integration_qag(&F, z_n[0], z_n[nz_n-1], 0,
                                  ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                  ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                  w, &nz_norm, &nz_enorm);
    gsl_integration_workspace_free(w);
    if(gslstatus!=GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: cl_tracer():");
      ccl_spline_free(clt->spl_nz);
      free(clt);
      *status=CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): integration error when normalizing N(z)\n");
      return NULL;
    }
    for(int ii=0;ii<nz_n;ii++)
      nz_normalized[ii]=n[ii]/nz_norm;
    ccl_spline_free(clt->spl_nz);
    clt->spl_nz=ccl_spline_init(nz_n,z_n,nz_normalized,0,0);
    free(nz_normalized);
    if(clt->spl_nz==NULL) {
      free(clt);
      *status=CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): error initializing normalized spline for N(z)\n");
      return NULL;
    }

    if(tracer_type==CL_TRACER_NC) {
      //Initialize bias spline
      clt->spl_bz=ccl_spline_init(nz_b,z_b,b,b[0],b[nz_b-1]);
      if(clt->spl_bz==NULL) {
	ccl_spline_free(clt->spl_nz);
	free(clt);
	*status=CCL_ERROR_SPLINE;
	ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): error initializing spline for b(z)\n");
	return NULL;
      }
      clt->has_rsd=has_rsd;
      clt->has_magnification=has_magnification;
      if(clt->has_magnification) {
	//Compute weak lensing kernel
	int nchi;
	double *x,*y;
	double dchi_here=5.;
	double zmax=clt->spl_nz->xf;
	double chimax=ccl_comoving_radial_distance(cosmo,1./(1+zmax),status);
	//TODO: The interval in chi (5. Mpc) should be made a macro

	//In this case we need to integrate all the way to z=0. Reset zmin and chimin
	clt->zmin=0;
	clt->chimin=0;
	clt->spl_sz=ccl_spline_init(nz_s,z_s,s,s[0],s[nz_s-1]);
	if(clt->spl_sz==NULL) {
	  ccl_spline_free(clt->spl_nz);
	  ccl_spline_free(clt->spl_bz);
	  free(clt);
	  *status=CCL_ERROR_SPLINE;
	  ccl_cosmology_set_status_message(cosmo,
		 "ccl_cls.c: ccl_cl_tracer(): error initializing spline for s(z)\n");
	  return NULL;
	}

	nchi=(int)(chimax/dchi_here)+1;
	x=ccl_linear_spacing(0.,chimax,nchi);
	dchi_here=chimax/nchi;
	if(x==NULL || (fabs(x[0]-0)>1E-5) || (fabs(x[nchi-1]-chimax)>1e-5)) {
	  ccl_spline_free(clt->spl_nz);
	  ccl_spline_free(clt->spl_bz);
	  ccl_spline_free(clt->spl_sz);
	  free(clt);
	  *status=CCL_ERROR_LINSPACE;
	  ccl_cosmology_set_status_message(cosmo,
		 "ccl_cls.c: ccl_cl_tracer(): Error creating linear spacing in chi\n");
	  return NULL;
	}
	y=(double *)malloc(nchi*sizeof(double));
	if(y==NULL) {
	  free(x);
	  ccl_spline_free(clt->spl_nz);
	  ccl_spline_free(clt->spl_bz);
	  ccl_spline_free(clt->spl_sz);
	  free(clt);
	  *status=CCL_ERROR_MEMORY;
	  ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): memory allocation\n");
	  return NULL;
	}

	for(int j=0;j<nchi;j++)
	  clstatus|=window_magnification(x[j],cosmo,clt->spl_nz,clt->spl_sz,chimax,&(y[j]));
	if(clstatus) {
	  free(y);
	  free(x);
	  ccl_spline_free(clt->spl_nz);
	  ccl_spline_free(clt->spl_bz);
	  ccl_spline_free(clt->spl_sz);
	  free(clt);
	  *status=CCL_ERROR_INTEG;
	  ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): error computing lensing window\n");
	  return NULL;
	}

	clt->spl_wM=ccl_spline_init(nchi,x,y,y[0],0);
	if(clt->spl_wM==NULL) {
	  free(y);
	  free(x);
	  ccl_spline_free(clt->spl_nz);
	  ccl_spline_free(clt->spl_bz);
	  ccl_spline_free(clt->spl_sz);
	  free(clt);
	  *status=CCL_ERROR_SPLINE;
	  ccl_cosmology_set_status_message(cosmo,
		 "ccl_cls.c: ccl_cl_tracer(): error initializing spline for lensing window\n");
	  return NULL;
	}
	free(x); free(y);
      }
    }
    else if(tracer_type==CL_TRACER_WL) {
      //Compute weak lensing kernel
      int nchi;
      double *x,*y;
      double dchi_here=5.;
      double zmax=clt->spl_nz->xf;
      double chimax=ccl_comoving_radial_distance(cosmo,1./(1+zmax),status);
      //TODO: The interval in chi (5. Mpc) should be made a macro

      //In this case we need to integrate all the way to z=0. Reset zmin and chimin
      clt->zmin=0;
      clt->chimin=0;
      nchi=(int)(chimax/dchi_here)+1;
      x=ccl_linear_spacing(0.,chimax,nchi);
      dchi_here=chimax/nchi;
      if(x==NULL || (fabs(x[0]-0)>1E-5) || (fabs(x[nchi-1]-chimax)>1e-5)) {
	ccl_spline_free(clt->spl_nz);
	free(clt);
	*status=CCL_ERROR_LINSPACE;
	ccl_cosmology_set_status_message(cosmo,
	       "ccl_cls.c: ccl_cl_tracer(): Error creating linear spacing in chi\n");
	return NULL;
      }
      y=(double *)malloc(nchi*sizeof(double));
      if(y==NULL) {
	free(x);
	ccl_spline_free(clt->spl_nz);
	free(clt);
	*status=CCL_ERROR_MEMORY;
	ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): memory allocation\n");
	return NULL;
      }

      for(int j=0;j<nchi;j++)
	clstatus|=window_lensing(x[j],cosmo,clt->spl_nz,chimax,&(y[j]));
      if(clstatus) {
	free(y);
	free(x);
	ccl_spline_free(clt->spl_nz);
	free(clt);
	*status=CCL_ERROR_INTEG;
	ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): error computing lensing window\n");
	return NULL;
      }

      clt->spl_wL=ccl_spline_init(nchi,x,y,y[0],0);
      if(clt->spl_wL==NULL) {
	free(y);
	free(x);
	ccl_spline_free(clt->spl_nz);
	free(clt);
	*status=CCL_ERROR_SPLINE;
	ccl_cosmology_set_status_message(cosmo,
	       "ccl_cls.c: ccl_cl_tracer(): error initializing spline for lensing window\n");
	return NULL;
      }
      free(x); free(y);

      clt->has_intrinsic_alignment=has_intrinsic_alignment;
      if(clt->has_intrinsic_alignment) {
	clt->spl_rf=ccl_spline_init(nz_rf,z_rf,rf,rf[0],rf[nz_rf-1]);
	if(clt->spl_rf==NULL) {
	  ccl_spline_free(clt->spl_nz);
	  free(clt);
	  *status=CCL_ERROR_SPLINE;
	  ccl_cosmology_set_status_message(cosmo,
		 "ccl_cls.c: ccl_cl_tracer(): error initializing spline for rf(z)\n");
	  return NULL;
	}
	clt->spl_ba=ccl_spline_init(nz_ba,z_ba,ba,ba[0],ba[nz_ba-1]);
	if(clt->spl_ba==NULL) {
	  ccl_spline_free(clt->spl_rf);
	  ccl_spline_free(clt->spl_nz);
	  free(clt);
	  *status=CCL_ERROR_SPLINE;
	  ccl_cosmology_set_status_message(cosmo,
		 "ccl_cls.c: ccl_cl_tracer(): error initializing spline for ba(z)\n");
	  return NULL;
	}
      }
    }
  }
  else if(tracer_type==CL_TRACER_CL) {
    clt->chi_source=ccl_comoving_radial_distance(cosmo,1./(1+z_source),status);
    clt->chimax=clt->chi_source;
    clt->chimin=0;
  }
  else {
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_tracer(): unknown tracer type\n");
    return NULL;
  }

  return clt;
}

//CCL_ClTracer constructor with error checking
//cosmo   -> ccl_cosmology object
//tracer_type -> type of tracer. Supported: CL_TRACER_NC, CL_TRACER_WL
//nz_n -> number of points for N(z)
//z_n  -> array of z-values for N(z)
//n    -> corresponding N(z)-values. Normalization is irrelevant
//        N(z) will be set to zero outside the range covered by z_n
//nz_b -> number of points for b(z)
//z_b  -> array of z-values for b(z)
//b    -> corresponding b(z)-values.
//        b(z) will be assumed constant outside the range covered by z_n
CCL_ClTracer *ccl_cl_tracer(ccl_cosmology *cosmo,int tracer_type,
				int has_rsd,int has_magnification,int has_intrinsic_alignment,
				int nz_n,double *z_n,double *n,
				int nz_b,double *z_b,double *b,
				int nz_s,double *z_s,double *s,
				int nz_ba,double *z_ba,double *ba,
				int nz_rf,double *z_rf,double *rf,
				double z_source, int * status)
{
  CCL_ClTracer *clt=cl_tracer(cosmo,tracer_type,has_rsd,has_magnification,has_intrinsic_alignment,
			      nz_n,z_n,n,nz_b,z_b,b,nz_s,z_s,s,
			      nz_ba,z_ba,ba,nz_rf,z_rf,rf,z_source,status);
  ccl_check_status(cosmo,status);
  return clt;
}

//CCL_ClTracer destructor
void ccl_cl_tracer_free(CCL_ClTracer *clt)
{
  if((clt->tracer_type==CL_TRACER_NC) || (clt->tracer_type==CL_TRACER_WL))
    ccl_spline_free(clt->spl_nz);

  if(clt->tracer_type==CL_TRACER_NC) {
    ccl_spline_free(clt->spl_bz);
    if(clt->has_magnification) {
      ccl_spline_free(clt->spl_sz);
      ccl_spline_free(clt->spl_wM);
    }
  }
  else if(clt->tracer_type==CL_TRACER_WL) {
    ccl_spline_free(clt->spl_wL);
    if(clt->has_intrinsic_alignment) {
      ccl_spline_free(clt->spl_ba);
      ccl_spline_free(clt->spl_rf);
    }
  }
  if(clt->computed_transfer) {
    int il;
    free(clt->n_k);
    for(il=0;il<clt->n_ls;il++)
      ccl_spline_free(clt->spl_transfer[il]);
    free(clt->spl_transfer);
  }
  free(clt);
}

CCL_ClTracer *ccl_cl_tracer_cmblens(ccl_cosmology *cosmo,double z_source,int *status)
{
  return ccl_cl_tracer(cosmo,CL_TRACER_CL,
			   0,0,0,
			   0,NULL,NULL,0,NULL,NULL,0,NULL,NULL,
			   0,NULL,NULL,0,NULL,NULL,z_source,status);
}

CCL_ClTracer *ccl_cl_tracer_number_counts(ccl_cosmology *cosmo,
					      int has_rsd,int has_magnification,
					      int nz_n,double *z_n,double *n,
					      int nz_b,double *z_b,double *b,
					      int nz_s,double *z_s,double *s, int * status)
{
  return ccl_cl_tracer(cosmo,CL_TRACER_NC,has_rsd,has_magnification,0,
			   nz_n,z_n,n,nz_b,z_b,b,nz_s,z_s,s,
			   -1,NULL,NULL,-1,NULL,NULL,0, status);
}

CCL_ClTracer *ccl_cl_tracer_number_counts_simple(ccl_cosmology *cosmo,
						     int nz_n,double *z_n,double *n,
						     int nz_b,double *z_b,double *b, int * status)
{
  return ccl_cl_tracer(cosmo,CL_TRACER_NC,0,0,0,
			   nz_n,z_n,n,nz_b,z_b,b,-1,NULL,NULL,
			   -1,NULL,NULL,-1,NULL,NULL,0, status);
}

CCL_ClTracer *ccl_cl_tracer_lensing(ccl_cosmology *cosmo,
					int has_alignment,
					int nz_n,double *z_n,double *n,
					int nz_ba,double *z_ba,double *ba,
					int nz_rf,double *z_rf,double *rf, int * status)
{
  return ccl_cl_tracer(cosmo,CL_TRACER_WL,0,0,has_alignment,
			   nz_n,z_n,n,-1,NULL,NULL,-1,NULL,NULL,
			   nz_ba,z_ba,ba,nz_rf,z_rf,rf,0, status);
}

CCL_ClTracer *ccl_cl_tracer_lensing_simple(ccl_cosmology *cosmo,
					       int nz_n,double *z_n,double *n, int * status)
{
  return ccl_cl_tracer(cosmo,CL_TRACER_WL,0,0,0,
			   nz_n,z_n,n,-1,NULL,NULL,-1,NULL,NULL,
			   -1,NULL,NULL,-1,NULL,NULL,0, status);
}

static void limits_bessel(double l,double thr,double *xmin,double *xmax)
{
  double thrb=thr*0.5635/pow(l+0.53,0.834);
  *xmax=1./thrb;
  if(l<=0)
    *xmin=0;
  else {
    double logxmin=((l+1.)*(log(l+1.)+M_LN2-1.)+0.5*M_LN2+log(thrb))/l;
    *xmin=exp(logxmin);
  }
}

static double j_bessel_limber(int l,double k)
{
  return sqrt(M_PI/(2*l+1.))/k;
}

static double f_dens(double a,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double z=1./a-1;
  double pz=ccl_spline_eval(z,clt->spl_nz);
  double bz=ccl_spline_eval(z,clt->spl_bz);
  double h=cosmo->params.h*ccl_h_over_h0(cosmo,a,status)/CLIGHT_HMPC;

  return pz*bz*h;
}

static double f_rsd(double a,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double z=1./a-1;
  double pz=ccl_spline_eval(z,clt->spl_nz);
  double fg=ccl_growth_rate(cosmo,a,status);
  double h=cosmo->params.h*ccl_h_over_h0(cosmo,a,status)/CLIGHT_HMPC;

  return pz*fg*h;
}

static double f_mag(double a,double chi,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double wM=ccl_spline_eval(chi,clt->spl_wM);

  if(wM<=0)
    return 0;
  else
    return wM/(a*chi);
}

//Transfer function for number counts
//l -> angular multipole
//k -> wavenumber modulus
//cosmo -> ccl_cosmology object
//w -> CCL_ClWorskpace object
//clt -> CCL_ClTracer object (must be of the CL_TRACER_NC type)
static double transfer_nc(int l,double k,
			  ccl_cosmology *cosmo,CCL_ClWorkspace *w,CCL_ClTracer *clt, int * status)
{
  double ret=0;
  if(l>w->l_limber) {
    double x0=(l+0.5);
    double chi0=x0/k;
    if(chi0<=clt->chimax) {
      double a0=ccl_scale_factor_of_chi(cosmo,chi0,status);
      double pk0=ccl_nonlin_matter_power(cosmo,k,a0,status);
      double jl0=j_bessel_limber(l,k);
      double f_all=f_dens(a0,cosmo,clt,status)*jl0;
      if(clt->has_rsd) {
	double x1=(l+1.5);
	double chi1=x1/k;
	if(chi1<=clt->chimax) {
	  double a1=ccl_scale_factor_of_chi(cosmo,chi1,status);
	  double pk1=ccl_nonlin_matter_power(cosmo,k,a1,status);
	  double fg0=f_rsd(a0,cosmo,clt,status);
	  double fg1=f_rsd(a1,cosmo,clt,status);
	  double jl1=j_bessel_limber(l+1,k);
	  f_all+=fg0*(1.-l*(l-1.)/(x0*x0))*jl0-fg1*2.*jl1*sqrt(pk1/pk0)/x1;
	}
      }
      if(clt->has_magnification)
      	f_all+=-2*clt->prefac_lensing*l*(l+1)*f_mag(a0,chi0,cosmo,clt,status)*jl0/(k*k);
      ret=f_all*sqrt(pk0);
    }
  }
  else {
    int i,nchi=(int)((clt->chimax-clt->chimin)/w->dchi)+1;
    for(i=0;i<nchi;i++) {
      double chi=clt->chimin+w->dchi*(i+0.5);
      if(chi<=clt->chimax) {
	double a=ccl_scale_factor_of_chi(cosmo,chi,status);
	double pk=ccl_nonlin_matter_power(cosmo,k,a,status);
	double jl=ccl_j_bessel(l,k*chi);
	double f_all=f_dens(a,cosmo,clt,status)*jl;
	if(clt->has_rsd) {
	  double ddjl,x=k*chi;
	  if(x<1E-10) {
	    if(l==0) ddjl=0.3333-0.1*x*x;
	    else if(l==2) ddjl=-0.13333333333+0.05714285714285714*x*x;
	    else ddjl=0;
          }
	  else {
	    double jlp1=ccl_j_bessel(l+1,x);
	    ddjl=((x*x-l*(l-1))*jl-2*x*jlp1)/(x*x);
	  }
	  f_all+=f_rsd(a,cosmo,clt,status)*ddjl;
	}
	if(clt->has_magnification)
	  f_all+=-2*clt->prefac_lensing*l*(l+1)*f_mag(a,chi,cosmo,clt,status)*jl/(k*k);

	ret+=f_all*sqrt(pk); //TODO: is it worth splining this sqrt?
      }
    }
    ret*=w->dchi;
  }

  return ret;
}

static double f_lensing(double a,double chi,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  double wL=ccl_spline_eval(chi,clt->spl_wL);

  if(wL<=0)
    return 0;
  else
    return clt->prefac_lensing*wL/(a*chi);
}

static double f_IA_NLA(double a,double chi,ccl_cosmology *cosmo,CCL_ClTracer *clt, int * status)
{
  if(chi<=1E-10)
    return 0;
  else {
    double a=ccl_scale_factor_of_chi(cosmo,chi, status);
    double z=1./a-1;
    double pz=ccl_spline_eval(z,clt->spl_nz);
    double ba=ccl_spline_eval(z,clt->spl_ba);
    double rf=ccl_spline_eval(z,clt->spl_rf);
    double h=cosmo->params.h*ccl_h_over_h0(cosmo,a,status)/CLIGHT_HMPC;

    return pz*ba*rf*h/(chi*chi);
  }
}

//Transfer function for shear
//l -> angular multipole
//k -> wavenumber modulus
//cosmo -> ccl_cosmology object
//w -> CCL_ClWorskpace object
//clt -> CCL_ClTracer object (must be of the CL_TRACER_WL type)
static double transfer_wl(int l,double k,
			  ccl_cosmology *cosmo,CCL_ClWorkspace *w,CCL_ClTracer *clt, int * status)
{
  double ret=0;
  if(l>w->l_limber) {
    double chi=(l+0.5)/k;
    if(chi<=clt->chimax) {
      double a=ccl_scale_factor_of_chi(cosmo,chi,status);
      double pk=ccl_nonlin_matter_power(cosmo,k,a,status);
      double jl=j_bessel_limber(l,k);
      double f_all=f_lensing(a,chi,cosmo,clt,status)*jl;
      if(clt->has_intrinsic_alignment)
	f_all+=f_IA_NLA(a,chi,cosmo,clt,status)*jl;

      ret=f_all*sqrt(pk);
    }
  }
  else {
    int i,nchi=(int)((clt->chimax-clt->chimin)/w->dchi)+1;
    for(i=0;i<nchi;i++) {
      double chi=clt->chimin+w->dchi*(i+0.5);
      if(chi<=clt->chimax) {
	double a=ccl_scale_factor_of_chi(cosmo,chi,status);
	double pk=ccl_nonlin_matter_power(cosmo,k,a,status);
	double jl=ccl_j_bessel(l,k*chi);
	double f_all=f_lensing(a,chi,cosmo,clt,status)*jl;
	if(clt->has_intrinsic_alignment)
	  f_all+=f_IA_NLA(a,chi,cosmo,clt,status)*jl;

	ret+=f_all*sqrt(pk); //TODO: is it worth splining this sqrt?
      }
    }
    ret*=w->dchi;
  }

  return sqrt((l+2.)*(l+1.)*l*(l-1.))*ret/(k*k);
  //return (l+1.)*l*ret/(k*k);
}

static double transfer_cmblens(int l,double k,ccl_cosmology *cosmo,CCL_ClTracer *clt,int *status)
{
  double chi=(l+0.5)/k;
  if(chi>=clt->chi_source)
    return 0;

  if(chi<=clt->chimax) {
    double a=ccl_scale_factor_of_chi(cosmo,chi,status);
    double w=1-chi/clt->chi_source;
    double jl=j_bessel_limber(l,k);
    double pk=ccl_nonlin_matter_power(cosmo,k,a,status);
    return clt->prefac_lensing*l*(l+1.)*w*sqrt(pk)*jl/(a*chi*k*k);
  }
  return 0;
}

//Wrapper for transfer function
//l -> angular multipole
//k -> wavenumber modulus
//cosmo -> ccl_cosmology object
//clt -> CCL_ClTracer object
static double transfer_wrap(int il,double lk,ccl_cosmology *cosmo,
			    CCL_ClWorkspace *w,CCL_ClTracer *clt, int * status)
{
  double transfer_out=0;
  double k=pow(10.,lk);

  if(clt->tracer_type==CL_TRACER_NC)
    transfer_out=transfer_nc(w->l_arr[il],k,cosmo,w,clt,status);
  else if(clt->tracer_type==CL_TRACER_WL)
    transfer_out=transfer_wl(w->l_arr[il],k,cosmo,w,clt,status);
  else if(clt->tracer_type==CL_TRACER_CL)
    transfer_out=transfer_cmblens(w->l_arr[il],k,cosmo,clt,status);
  else
    transfer_out=-1;
  return transfer_out;
}

static double *get_lkarr(ccl_cosmology *cosmo,CCL_ClWorkspace *w,
			 double l,double chimin,double chimax,int *nk,
			 int *status)
{
  int ik;

  //First compute relevant k-range for this ell
  double kmin,kmax,lkmin,lkmax;
  if(l>w->l_limber) {
    kmin=CCL_MAX(ccl_splines->K_MIN,0.8*(l+0.5)/chimax);
    kmax=CCL_MIN(ccl_splines->K_MAX,1.2*(l+0.5)/chimin);
  }
  else {
    double xmin,xmax;
    limits_bessel(l,CCL_FRAC_RELEVANT,&xmin,&xmax);
    kmin=CCL_MAX(ccl_splines->K_MIN,xmin/chimax);
    kmax=CCL_MIN(ccl_splines->K_MAX,xmax/chimin);
    //Cap by maximum meaningful argument of the Bessel function
    kmax=CCL_MIN(kmax,2*(w->l_arr[w->n_ls-1]+0.5)/chimin); //Cap by 2 x inverse scale corresponding to l_max
  }
  lkmin=log10(kmin);
  lkmax=log10(kmax);

  //Allocate memory for transfer function
  double *lkarr;
  double lknew=lkmin;
  double k_period=2*M_PI/chimax;
  double dklin=0.45;
  ik=0;
  while(lknew<=lkmax) {
    double kk=pow(10.,lknew);
    double dk1=k_period*dklin;
    double dk2=kk*(pow(10.,w->dlk)-1.);
    double dk=CCL_MIN(dk1,dk2);
    lknew=log10(kk+dk);
    ik++;
  }

  *nk=CCL_MAX(10,ik+1);
  lkarr=(double *)malloc((*nk)*sizeof(double));
  if(lkarr==NULL) {
    return NULL;
  }

  if((*nk)==10) {
    for(ik=0;ik<(*nk);ik++)
      lkarr[ik]=lkmin+(lkmax-lkmin)*ik/((*nk)-1.);
  }
  else {
    ik=0;
    lkarr[ik]=lkmin;
    while(lkarr[ik]<=lkmax) {
      double kk=pow(10.,lkarr[ik]);
      double dk1=k_period*dklin;
      double dk2=kk*(pow(10.,w->dlk)-1.);
      double dk=CCL_MIN(dk1,dk2);
      ik++;
      lkarr[ik]=log10(kk+dk);
    }
  }

  return lkarr;
}

static void compute_transfer(CCL_ClTracer *clt,ccl_cosmology *cosmo,CCL_ClWorkspace *w,int *status)
{
  int il;
  double zmin=CCL_MAX(w->zmin,clt->zmin);
  double chimin=ccl_comoving_radial_distance(cosmo,1./(1+zmin),status);
  double chimax=clt->chimax;

  //Get how many multipoles and allocate info for each of them
  clt->n_ls=w->n_ls;
  clt->n_k=(int *)malloc(clt->n_ls*sizeof(int));
  if(clt->n_k==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: compute_transfer(): memory allocation\n");
    return;
  }
  clt->spl_transfer=(SplPar **)malloc(clt->n_ls*sizeof(SplPar *));
  if(clt->spl_transfer==NULL) {
    free(clt->n_k);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: compute_transfer(): memory allocation\n");
    return;
  }

  //Loop over multipoles and compute transfer function for each
  for(il=0;il<clt->n_ls;il++) {
    int ik,nk;
    double l=(double)(w->l_arr[il]);
    double *lkarr=get_lkarr(cosmo,w,l,chimin,chimax,&nk,status);
    if(lkarr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: compute_transfer(): memory allocation\n");
      break;
    }

    double *tkarr=(double *)malloc(nk*sizeof(double));
    if(tkarr==NULL) {
      free(lkarr);
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: compute_transfer(): memory allocation\n");
      break;
    }
    clt->n_k[il]=nk;

    //Loop over k and compute transfer function
    for(ik=0;ik<nk;ik++)
      tkarr[ik]=transfer_wrap(il,lkarr[ik],cosmo,w,clt,status);
    if(*status) {
      free(clt->n_k);
      free(lkarr);
      free(tkarr);
      break;
    }

    //Initialize spline for this ell
    clt->spl_transfer[il]=ccl_spline_init(nk,lkarr,tkarr,0,0);
    if(clt->spl_transfer[il]==NULL) {
      free(lkarr);
      free(tkarr);
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: compute_transfer(): memory allocation\n");
      break;
    }
    free(lkarr);
    free(tkarr);
  }
  if(*status) {
    free(clt->n_k);
    int ill;
    for(ill=0;ill<il;ill++)
      ccl_spline_free(clt->spl_transfer[il]);
    free(clt->spl_transfer);
  }

  clt->computed_transfer=1;
}

static double transfer(int il,double lk,ccl_cosmology *cosmo,
		       CCL_ClWorkspace *w,CCL_ClTracer *clt,int *status)
{
  if(il<w->l_limber) {
    if(!(clt->computed_transfer))
      compute_transfer(clt,cosmo,w,status);

    return ccl_spline_eval(lk,clt->spl_transfer[il]);
  } else {
    return transfer_wrap(il,lk,cosmo,w,clt,status);
  }
}

//Params for power spectrum integrand
typedef struct {
  int il;
  ccl_cosmology *cosmo;
  CCL_ClWorkspace *w;
  CCL_ClTracer *clt1;
  CCL_ClTracer *clt2;
  int *status;
} IntClPar;

//Integrand for integral power spectrum
static double cl_integrand(double lk,void *params)
{
  double d1,d2;
  IntClPar *p=(IntClPar *)params;
  d1=transfer(p->il,lk,p->cosmo,p->w,p->clt1,p->status);
  d2=transfer(p->il,lk,p->cosmo,p->w,p->clt2,p->status);

  return pow(10.,3*lk)*d1*d2;
}

//Figure out k intervals where the Limber kernel has support
//clt1 -> tracer #1
//clt2 -> tracer #2
//l    -> angular multipole
//lkmin, lkmax -> log10 of the range of scales where the transfer functions have support
static void get_k_interval(ccl_cosmology *cosmo,CCL_ClWorkspace *w,
			   CCL_ClTracer *clt1,CCL_ClTracer *clt2,int l,
			   double *lkmin,double *lkmax)
{
  double chimin,chimax;
  int cut_low_1=0,cut_low_2=0;

  //Define a minimum distance only if no lensing is needed
  if((clt1->tracer_type==CL_TRACER_NC) && (clt1->has_magnification==0)) cut_low_1=1;
  if((clt2->tracer_type==CL_TRACER_NC) && (clt2->has_magnification==0)) cut_low_2=1;

  if(l<w->l_limber) {
    chimin=2*(l+0.5)/ccl_splines->K_MAX;
    chimax=0.5*(l+0.5)/ccl_splines->K_MIN;
  }
  else {
    if(cut_low_1) {
      if(cut_low_2) {
	chimin=fmax(clt1->chimin,clt2->chimin);
	chimax=fmin(clt1->chimax,clt2->chimax);
      }
      else {
	chimin=clt1->chimin;
	chimax=clt1->chimax;
      }
    }
    else if(cut_low_2) {
      chimin=clt2->chimin;
      chimax=clt2->chimax;
    }
    else {
      chimin=0.5*(l+0.5)/ccl_splines->K_MAX;
      chimax=2*(l+0.5)/ccl_splines->K_MIN;
    }
  }

  if(chimin<=0)
    chimin=0.5*(l+0.5)/ccl_splines->K_MAX;

  *lkmax=fmin( 2,log10(2  *(l+0.5)/chimin));
  *lkmin=fmax(-4,log10(0.5*(l+0.5)/chimax));
}

//Compute angular power spectrum between two bins
//cosmo -> ccl_cosmology object
//il -> index in angular multipole array
//clt1 -> tracer #1
//clt2 -> tracer #2
static double ccl_angular_cl_native(ccl_cosmology *cosmo,CCL_ClWorkspace *cw,int il,
				    CCL_ClTracer *clt1,CCL_ClTracer *clt2,int * status)
{
  int clastatus=0, gslstatus;
  IntClPar ipar;
  double result=0,eresult;
  double lkmin,lkmax;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);

  ipar.il=il;
  ipar.cosmo=cosmo;
  ipar.w=cw;
  ipar.clt1=clt1;
  ipar.clt2=clt2;
  ipar.status = &clastatus;
  F.function=&cl_integrand;
  F.params=&ipar;
  get_k_interval(cosmo,cw,clt1,clt2,cw->l_arr[il],&lkmin,&lkmax);
  gslstatus=gsl_integration_qag(&F, lkmin, lkmax, 0,
                                ccl_gsl->INTEGRATION_LIMBER_EPSREL, ccl_gsl->N_ITERATION,
                                ccl_gsl->INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS,
                                w, &result, &eresult);
  gsl_integration_workspace_free(w);

  // Test if a round-off error occured in the evaluation of the integral
  // If so, try another integration function, more robust but potentially slower
  if(gslstatus == GSL_EROUND) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: ccl_angular_cl_native(): Default GSL integration failure, attempting backup method.");
    gsl_integration_cquad_workspace *w_cquad= gsl_integration_cquad_workspace_alloc (ccl_gsl->N_ITERATION);
    size_t nevals=0;
    gslstatus=gsl_integration_cquad(&F, lkmin, lkmax, 0,
				    ccl_gsl->INTEGRATION_LIMBER_EPSREL,
				    w_cquad, &result, &eresult, &nevals);
    gsl_integration_cquad_workspace_free(w_cquad);
  }
  if(gslstatus!=GSL_SUCCESS || *ipar.status) {
    ccl_raise_gsl_warning(gslstatus, "ccl_cls.c: ccl_angular_cl_native():");
    // If an error status was already set, don't overwrite it.
    if(*status == 0){
        *status=CCL_ERROR_INTEG;
        ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_angular_cl_native(): error integrating over k\n");
    }
    return -1;
  }
  ccl_check_status(cosmo,status);

  return result*M_LN10*2./M_PI;
}

void ccl_angular_cls(ccl_cosmology *cosmo,CCL_ClWorkspace *w,
		     CCL_ClTracer *clt1,CCL_ClTracer *clt2,
		     int nl_out,int *l_out,double *cl_out,int *status)
{
  int ii;
  //First check if ell range is within workspace
  for(ii=0;ii<nl_out;ii++) {
    if(l_out[ii]>w->lmax) {
      *status=CCL_ERROR_SPLINE_EV;
      ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_angular_cls(); "
	     "requested l beyond range allowed by workspace\n");
      return;
    }
  }

  //Allocate array for power spectrum at interpolation nodes
  double *l_nodes=(double *)malloc(w->n_ls*sizeof(double));
  if(l_nodes==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_angular_cls(); memory allocation\n");
    return;
  }
  double *cl_nodes=(double *)malloc(w->n_ls*sizeof(double));
  if(cl_nodes==NULL) {
    free(l_nodes);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_angular_cls(); memory allocation\n");
    return;
  }
  for(ii=0;ii<w->n_ls;ii++)
    l_nodes[ii]=(double)(w->l_arr[ii]);

  //Now check if angpow is needed at all
  int method_use=w->nlimb_method;
  if(method_use==CCL_NONLIMBER_METHOD_ANGPOW) {
    int do_angpow=0;
    for(ii=0;ii<w->n_ls;ii++) {
      if(w->l_arr[ii]<=w->l_limber)
	do_angpow=1;
    }
    //Resort to native method if we have lensing (this will hopefully only be temporary)
    if(clt1->tracer_type==CL_TRACER_WL || clt2->tracer_type==CL_TRACER_WL ||
       clt1->has_magnification || clt2->has_magnification) {
      do_angpow=0;
      method_use=CCL_NONLIMBER_METHOD_NATIVE;
    }
    //Use angpow if non-limber is needed
#ifdef HAVE_ANGPOW
    if(do_angpow)
      ccl_angular_cls_angpow(cosmo,w,clt1,clt2,cl_nodes,status);
    ccl_check_status(cosmo,status);
#else
    do_angpow=0;
    method_use=CCL_NONLIMBER_METHOD_NATIVE;
#endif
  }

  //Compute limber nodes
  for(ii=0;ii<w->n_ls;ii++) {
    if((method_use==CCL_NONLIMBER_METHOD_NATIVE) || (w->l_arr[ii]>w->l_limber))
      cl_nodes[ii]=ccl_angular_cl_native(cosmo,w,ii,clt1,clt2,status);
  }

  //Interpolate into ells requested by user
  SplPar *spcl_nodes=ccl_spline_init(w->n_ls,l_nodes,cl_nodes,0,0);
  if(spcl_nodes==NULL) {
    free(cl_nodes);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: ccl_cl_angular_cls(); memory allocation\n");
    return;
  }
  for(ii=0;ii<nl_out;ii++)
    cl_out[ii]=ccl_spline_eval((double)(l_out[ii]),spcl_nodes);

  //Cleanup
  ccl_spline_free(spcl_nodes);
  free(cl_nodes);
  free(l_nodes);
}

static int check_clt_fa_inconsistency(CCL_ClTracer *clt,int func_code)
{
  if(((func_code==CCL_CLT_NZ) && (clt->tracer_type==CL_TRACER_CL)) || //Lensing has no N(z)
     (((func_code==CCL_CLT_BZ) || (func_code==CCL_CLT_SZ) || (func_code==CCL_CLT_WM)) &&
      (clt->tracer_type!=CL_TRACER_NC)) || //bias and magnification only for clustering
     (((func_code==CCL_CLT_RF) || (func_code==CCL_CLT_BA) || (func_code==CCL_CLT_WL)) &&
      (clt->tracer_type!=CL_TRACER_WL))) //IAs only for weak lensing
    return 1;
  if((((func_code==CCL_CLT_SZ) || (func_code==CCL_CLT_WM)) &&
      (clt->has_magnification==0)) || //Correct combination, but no magnification
     (((func_code==CCL_CLT_RF) || (func_code==CCL_CLT_BA)) &&
      (clt->has_intrinsic_alignment==0))) //Correct combination, but no IAs
    return 1;
  return 0;
}

double ccl_get_tracer_fa(ccl_cosmology *cosmo,CCL_ClTracer *clt,double a,int func_code,int *status)
{
  SplPar *spl;
  double x=1./a-1; //x-variable is redshift by default

  if(check_clt_fa_inconsistency(clt,func_code)) {
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: inconsistent combination of tracer and internal function to be evaluated");
    return -1;
  }

  if(func_code==CCL_CLT_NZ)
    spl=clt->spl_nz;
  if(func_code==CCL_CLT_BZ)
    spl=clt->spl_bz;
  if(func_code==CCL_CLT_SZ)
    spl=clt->spl_sz;
  if(func_code==CCL_CLT_RF)
    spl=clt->spl_rf;
  if(func_code==CCL_CLT_BA)
    spl=clt->spl_ba;
  if((func_code==CCL_CLT_WL) || (func_code==CCL_CLT_WM)) {
    x=ccl_comoving_radial_distance(cosmo,a,status);
    if(func_code==CCL_CLT_WL)
      spl=clt->spl_wL;
    if(func_code==CCL_CLT_WM)
      spl=clt->spl_wM;
  }

  return ccl_spline_eval(x,spl);
}

int ccl_get_tracer_fas(ccl_cosmology *cosmo,CCL_ClTracer *clt,int na,double *a,double *fa,
		       int func_code,int *status)
{
  SplPar *spl;

  if(check_clt_fa_inconsistency(clt,func_code)) {
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_cls.c: inconsistent combination of tracer and internal function to be evaluated");
    return -1;
  }

  if(func_code==CCL_CLT_NZ)
    spl=clt->spl_nz;
  if(func_code==CCL_CLT_BZ)
    spl=clt->spl_bz;
  if(func_code==CCL_CLT_SZ)
    spl=clt->spl_sz;
  if(func_code==CCL_CLT_RF)
    spl=clt->spl_rf;
  if(func_code==CCL_CLT_BA)
    spl=clt->spl_ba;
  if(func_code==CCL_CLT_WL)
    spl=clt->spl_wL;
  if(func_code==CCL_CLT_WM)
    spl=clt->spl_wM;

  int compchi=0;
  if((func_code==CCL_CLT_WL) || (func_code==CCL_CLT_WM))
    compchi=1;

  int ia;
  for(ia=0;ia<na;ia++) {
    double x;
    if(compchi)
      x=ccl_comoving_radial_distance(cosmo,a[ia],status);
    else
      x=1./a[ia]-1;
    fa[ia]=ccl_spline_eval(x,spl);
  }

  return 0;
}

ext: c

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
/*****************************************************************************
*
* Copyright (c) 2003-2018 by The University of Queensland
* http://www.uq.edu.au
*
* Primary Business: Queensland, Australia
* Licensed under the Apache License, version 2.0
* http://www.apache.org/licenses/LICENSE-2.0
*
* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
* Development 2012-2013 by School of Earth Sciences
* Development from 2014 by Centre for Geoscience Computing (GeoComp)
*
*****************************************************************************/

#ifndef __PASO_BLOCKOPS_H__
#define __PASO_BLOCKOPS_H__

#include "Paso.h"
#include "PasoException.h"

#include <cstring> // memcpy

#ifdef ESYS_HAVE_LAPACK
   #ifdef ESYS_MKL_LAPACK
      #include <mkl_lapack.h>
      #include <mkl_cblas.h>
   #else
      extern "C" {
      #include <clapack.h>
      #include <cblas.h>
      }
   #endif
#endif

namespace paso {

inline void BlockOps_Cpy_N(dim_t N, double* R, const double* V)
{
    memcpy((void*)R, (void*)V, N*sizeof(double));
}

/// performs operation R=R-mat*V (V and R are not overlapping) - 2x2
inline void BlockOps_SMV_2(double* R, const double* mat, const double* V)
{
    const double S1 = V[0];
    const double S2 = V[1];
    const double A11 = mat[0];
    const double A12 = mat[2];
    const double A21 = mat[1];
    const double A22 = mat[3];
    R[0] -= A11 * S1 + A12 * S2;
    R[1] -= A21 * S1 + A22 * S2;
}

/// performs operation R=R-mat*V (V and R are not overlapping) - 3x3
inline void BlockOps_SMV_3(double* R, const double* mat, const double* V)
{
    const double S1 = V[0];
    const double S2 = V[1];
    const double S3 = V[2];
    const double A11 = mat[0];
    const double A21 = mat[1];
    const double A31 = mat[2];
    const double A12 = mat[3];
    const double A22 = mat[4];
    const double A32 = mat[5];
    const double A13 = mat[6];
    const double A23 = mat[7];
    const double A33 = mat[8];
    R[0] -= A11 * S1 + A12 * S2 + A13 * S3;
    R[1] -= A21 * S1 + A22 * S2 + A23 * S3;
    R[2] -= A31 * S1 + A32 * S2 + A33 * S3;
}

#define PASO_MISSING_CLAPACK throw PasoException("You need to install a LAPACK version to enable operations on block sizes > 3.")

/// performs operation R=R-mat*V (V and R are not overlapping) - NxN
inline void BlockOps_SMV_N(dim_t N, double* R, const double* mat, const double* V)
{
#ifdef ESYS_HAVE_LAPACK
    cblas_dgemv(CblasColMajor,CblasNoTrans, N, N, -1., mat, N, V, 1, 1., R, 1);
#else
    PASO_MISSING_CLAPACK;
#endif
}

inline void BlockOps_MV_N(dim_t N, double* R, const double* mat, const double* V)
{
#ifdef ESYS_HAVE_LAPACK
    cblas_dgemv(CblasColMajor,CblasNoTrans, N, N, 1., mat, N, V, 1, 0., R, 1);
#else
    PASO_MISSING_CLAPACK;
#endif
}

inline void BlockOps_invM_2(double* invA, const double* A, int* failed)
{
    const double A11 = A[0];
    const double A12 = A[2];
    const double A21 = A[1];
    const double A22 = A[3];
    double D = A11*A22-A12*A21;
    if (std::abs(D) > 0) {
        D = 1./D;
        invA[0] =  A22*D;
        invA[1] = -A21*D;
        invA[2] = -A12*D;
        invA[3] =  A11*D;
    } else {
        *failed = 1;
    }
}

inline void BlockOps_invM_3(double* invA, const double* A, int* failed)
{
    const double A11 = A[0];
    const double A21 = A[1];
    const double A31 = A[2];
    const double A12 = A[3];
    const double A22 = A[4];
    const double A32 = A[5];
    const double A13 = A[6];
    const double A23 = A[7];
    const double A33 = A[8];
    double D = A11*(A22*A33-A23*A32) +
               A12*(A31*A23-A21*A33) +
               A13*(A21*A32-A31*A22);
    if (std::abs(D) > 0) {
        D = 1./D;
        invA[0] = (A22*A33-A23*A32)*D;
        invA[1] = (A31*A23-A21*A33)*D;
        invA[2] = (A21*A32-A31*A22)*D;
        invA[3] = (A13*A32-A12*A33)*D;
        invA[4] = (A11*A33-A31*A13)*D;
        invA[5] = (A12*A31-A11*A32)*D;
        invA[6] = (A12*A23-A13*A22)*D;
        invA[7] = (A13*A21-A11*A23)*D;
        invA[8] = (A11*A22-A12*A21)*D;
    } else {
        *failed = 1;
    }
}

/// LU factorization of NxN matrix mat with partial pivoting
inline void BlockOps_invM_N(dim_t N, double* mat, index_t* pivot, int* failed)
{
#ifdef ESYS_HAVE_LAPACK
#ifdef ESYS_MKL_LAPACK
    int res = 0;
    dgetrf(&N, &N, mat, &N, pivot, &res);
    if (res != 0)
        *failed = 1;
#else
    int res = clapack_dgetrf(CblasColMajor, N, N, mat, N, pivot);
    if (res != 0)
        *failed = 1;
#endif // ESYS_MKL_LAPACK
#else
    PASO_MISSING_CLAPACK;
#endif
}

/// solves system of linear equations A*X=B
inline void BlockOps_solve_N(dim_t N, double* X, double* mat, index_t* pivot, int* failed)
{
#ifdef ESYS_HAVE_LAPACK
#ifdef ESYS_MKL_LAPACK
    int res = 0;
    int ONE = 1;
    dgetrs("N", &N, &ONE, mat, &N, pivot, X, &N, &res);
    if (res != 0)
        *failed = 1;
#else
    int res = clapack_dgetrs(CblasColMajor, CblasNoTrans, N, 1, mat, N, pivot, X, N);
    if (res != 0)
        *failed = 1;
#endif // ESYS_MKL_LAPACK
#else
    PASO_MISSING_CLAPACK;
#endif
}

/// inplace matrix vector product - order 2
inline void BlockOps_MViP_2(const double* mat, double* V)
{
    const double S1 = V[0];
    const double S2 = V[1];
    const double A11 = mat[0];
    const double A12 = mat[2];
    const double A21 = mat[1];
    const double A22 = mat[3];
    V[0] = A11 * S1 + A12 * S2;
    V[1] = A21 * S1 + A22 * S2;
}

/// inplace matrix vector product - order 3
inline void BlockOps_MViP_3(const double* mat, double* V)
{
    const double S1 = V[0];
    const double S2 = V[1];
    const double S3 = V[2];
    const double A11 = mat[0];
    const double A21 = mat[1];
    const double A31 = mat[2];
    const double A12 = mat[3];
    const double A22 = mat[4];
    const double A32 = mat[5];
    const double A13 = mat[6];
    const double A23 = mat[7];
    const double A33 = mat[8];
    V[0] = A11 * S1 + A12 * S2 + A13 * S3;
    V[1] = A21 * S1 + A22 * S2 + A23 * S3;
    V[2] = A31 * S1 + A32 * S2 + A33 * S3;
}

inline void BlockOps_solveAll(dim_t n_block, dim_t n, double* D,
                              index_t* pivot, double* x)
{
    if (n_block == 1) {
#pragma omp parallel for
        for (dim_t i=0; i<n; ++i)
            x[i] *= D[i];
    } else if (n_block == 2) {
#pragma omp parallel for
        for (dim_t i=0; i<n; ++i)
            BlockOps_MViP_2(&D[4*i], &x[2*i]);
    } else if (n_block == 3) {
#pragma omp parallel for
        for (dim_t i=0; i<n; ++i)
            BlockOps_MViP_3(&D[9*i], &x[3*i]);
    } else {
        int failed = 0;
#pragma omp parallel for
        for (dim_t i=0; i<n; ++i) {
            const dim_t block_size = n_block*n_block;
            BlockOps_solve_N(n_block, &x[n_block*i], &D[block_size*i], &pivot[n_block*i], &failed);
        }
        if (failed > 0) {
            throw PasoException("BlockOps_solveAll: solution failed.");
        }
    }
}

} // namespace paso

#endif // __PASO_BLOCKOPS_H__

ext: h


text:
#include <stdio.h>
#include <math.h>
#include <ascii.h>
#include <cdisort.h>
#include <c_tzs.h>
#include "locate.h"
#if HAVE_LIBGSL
#include <gsl/gsl_sf_expint.h>
#endif

// include DISORT.MXD

/*  nlyr   = number of atmospheric layers                     */
/*  nlev_c = number of atmospheric levels on common z grid    */
/*  zd_c   = atmospheric levels on common z grid              */
/*  nzout  = number of output levels = number of black clouds */
/*  zout   = vertical level heights in km above surface       */

int c_tzs (int      nlyr,
           float*   dtauc,
           int      nlev_c,
           float*   zd_c,
           int      nzout,
           float*   zout,
           float*   ssalb,
           float*   temper,
           float    wvnmlo,
           float    wvnmhi,
           int      usrtau,
           int      ntau,
           float*   utau,
           int      usrang,
           int      numu,
           float*   umu,
           int      nphi,
           float*   phi,
           float    albedo,
           double   btemp,
           float    ttemp,
           float    temis,
           int      planck,
           int*     prnt,
           char*    header,  //  int prndis[7]; char header[127]
           float*   rfldir,
           float*   rfldn,
           float*   flup,
           float*   dfdt,
           float*   uavg,
           float*** uu,
           int      quiet) {
  /*    .. Internal variables  */
  int     lc, j, iu, lu;
  double* tauc = NULL;
  double  tau1 = 0, tau2 = 0, etau1 = 0, etau2 = 0, etaumu = 0, ei1 = 0, ei2 = 0, *rad = NULL;
  double *radem = NULL, *pkag = NULL, *etau = NULL, rademtot = 0.0;
  double  radrefl, radrefl1, radrefl2, radrefl11, radrefl22;
  double *b0 = NULL, *b1 = NULL;
  double *coeff1 = NULL, *coeff2 = NULL;
  double  sqtau1 = 0.0, sqtau2 = 0.0, mu = 0.0, cutau1 = 0.0, tauetau1;
  int *   nlayers = NULL, maxnlyr = 0;
  float * alb = NULL, *bplanck = NULL;

  /* nlyr     = number of layers */
  /* nlyr + 1 = number of levels */

  /* initialise */

  if ((rad = (double*)calloc (ntau, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((nlayers = (int*)calloc (ntau, sizeof (int))) == NULL)
    return ASCII_NO_MEMORY;

  if ((alb = (float*)calloc (ntau, sizeof (float))) == NULL)
    return ASCII_NO_MEMORY;
  if ((bplanck = (float*)calloc (ntau, sizeof (float))) == NULL)
    return ASCII_NO_MEMORY;

  if (prnt[0] && strlen (header) != 0) {
    fprintf (stderr, "**************\n TZS: %s\n*****************\n", header);
  }

  /* Initialise everything to 0. */
  for (iu = 0; iu < numu; iu++) {
    for (lu = 0; lu < nzout; lu++) {
      for (j = 0; j < nphi; j++) {
        uu[j][lu][iu] = 0.;
      }
    }
  }
  for (lu = 0; lu < nzout; lu++) {
    dfdt[lu]   = 0.;
    flup[lu]   = 0.;
    rfldir[lu] = 0.;
    rfldn[lu]  = 0.;
    uavg[lu]   = 0.;
  }

  /*    Layers: 0...NLYR-1 (# layers = NLYR)                                                           */
  /*    Levels: 0...NLYR   (# levels = NLYR+1)                                                         */
  /*    LC = 0...NYR-1                                                                                */
  /*    LC = 0 : TOA                                                                                */
  /*    -------------------------------------------------- Level LC, TAUC(LC), TEMPER(LC)              */
  /*                                                                                */
  /*                  Layer LC                                                                         */
  /*                                                                                */
  /*    -------------------------------------------------- Level LC+1  , TAUC(LC+1)  , TEMPER(LC+1)    */

  /*    Planck radiation from surface with emissivity=1-ALBEDO and temperature BTEMP */
  /*     BPLANCK = (1-ALBEDO)*PLKAVG2( WVNMLO, WVNMHI, BTEMP )  */

  /*    Max number of layers counted from TOA used for all UTAU  */
  maxnlyr = -1;

  /*    Loop over user defined tau = black cloud top heights:  */
  /*    determine layer number corresponding to user tau       */
  /*    Rem.: NTAU = NZOUT checked in CHEKIN_TZS               */

  if (ntau != nzout) {
    errmsg_tzs ("TZS--nzout != ntau", TZS_ERROR);
  }

  for (lu = 0; lu < ntau; lu++) {

    /* zd_c[0] = TOA                                                         */
    /* zd_c[nlyr+1] = SUR                                                    */
    /* locate gives back lc such that zd_c(lc) >= zout(lu) > zd_lc[lc+1]     */
    lc = flocate (zd_c, nlev_c, zout[lu]);

    //    fprintf(stderr,"lc = %d  zd_c[lc] = %f  zd_c[lc+1] = %f  zout[lu] = %f\n",lc,zd_c[lc],zd_c[lc+1],zout[lu]);

    if (lc == -1 || lc == nlev_c) {
      errmsg_tzs ("TZS--Zout outside atmospheric levels", TZS_ERROR);
    }

    /* Check real presence of this level in tauc ( max diff = 1mm )                         */
    /* Determine layer above given user level and store layer index w.r.t. zd_c in nlayers  */
    /* nlayers contains the layer index whose bottom is at zout[lu]                         */

    if (fabs (zd_c[lc + 1] - zout[lu]) <= 1.e-6) {

      /*          check lower boundary           */

      /*          nlevel = nlayers[lu] = LC+1-1  */
      /*          nlayer = nlevel - 1            */
      nlayers[lu] = lc;

    } else if (fabs (zd_c[lc] - zout[lu]) <= 1.e-6) {

      /*          check upper boundary           */

      /*          nlevel = nlayers[lu] = LC      */
      /*          nlayer = nlevel - 1            */
      nlayers[lu] = (int)fmax (lc - 1, 0);
      /* nlayers could be -1 when black cloud is at TOA */

    } else {

      /* zout_interpolate */

      //      fprintf(stderr,"%d4 %4d %9.3f %9.3f %9.3f %9.3f %9.3f\n",
      //	      lc,nlev_c,zd_c[lc+1],zd_c[lc],zout[lu],fabs(zd_c[lc+1]-zout[lu]),fabs(zd_c[lc]-zout[lu]));
      errmsg_tzs ("TZS--User height is not contained in atmospheric levels: use zout_interpolate!", TZS_ERROR);
    }

    //    fprintf(stderr,"nlev_c = %4d   zout[lu] = %9.3f  zd_c[nlayers[lu]] = %9.3f   zd_c[nlayers[lu]+1] = %9.3f   lu = %4d\n",
    //	    nlev_c,zout[lu],zd_c[nlayers[lu]],zd_c[nlayers[lu]+1],lu);
    //
    //    fprintf(stderr,"lu = %4d   maxnlyr = %4d   nlayers[lu] = %4d\n",lu,maxnlyr,nlayers[lu]);

    maxnlyr = (int)fmax (maxnlyr, nlayers[lu]);
  }

  /* Above we defined layer indices, now we add +1 because we want the number of layers */
  maxnlyr += 1;

  //  fprintf(stderr,"%s %4d\n","Max # layers = ",maxnlyr);

  if (maxnlyr < 0) {
    errmsg_tzs ("TZS--wrong zout", TZS_ERROR);
  }

  if ((radem = (double*)calloc (maxnlyr, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((tauc = (double*)calloc (maxnlyr + 1, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((pkag = (double*)calloc (maxnlyr + 1, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((etau = (double*)calloc (maxnlyr + 1, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;

  if ((b0 = (double*)calloc (maxnlyr, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((b1 = (double*)calloc (maxnlyr, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((coeff1 = (double*)calloc (maxnlyr, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;
  if ((coeff2 = (double*)calloc (maxnlyr, sizeof (double))) == NULL)
    return ASCII_NO_MEMORY;

  /*    calculate cumulative optical depth tauc[lc] up to level lc and      */
  /*    Planck radiation for every level (not layer!) from 0 to MAXNLYR     */
  /*    Remark: layer 1 is the first layer from above=TOA??                   */
  /*            level 0 is the first level from above=TOA, level 0 = TOA    */

  /* lc = level count */
  lc       = 0;
  tauc[lc] = 0.0;
  pkag[lc] = (double)c_planck_func1 (wvnmlo, wvnmhi, temper[lc]);

  //  fprintf(stderr,"lc = %d  maxnlyr = %d  pkag[lc] = %f  tauc[lc] = %13.6e dtauc[lc-1] = %13.6e temper[lc] = %9.3f\n",lc,maxnlyr,pkag[lc],tauc[lc],dtauc[lc-1],temper[lc]);

  for (lc = 1; lc <= maxnlyr; lc++) {
    /* We work here with atmospheric optical thickness as it is provided by the user. */
    /* If you want to use aborption optical thickness the you have to define it in    */
    /* the uvspec input file.                                                         */
    tauc[lc] = tauc[lc - 1] + ((double)dtauc[lc - 1]);

    pkag[lc] = (double)c_planck_func1 (wvnmlo, wvnmhi, temper[lc]);

    //    fprintf(stderr,"lc = %d  maxnlyr = %d  pkag[lc] = %f  tauc[lc] = %13.6e dtauc[lc-1] = %13.6e temper[lc] = %9.3f\n",lc,maxnlyr,pkag[lc],tauc[lc],dtauc[lc-1],temper[lc]);
  }

  /*    Print input information  */
  if (prnt[0]) {

    prtinp_tzs (nlyr,
                maxnlyr,
                dtauc,
                ssalb,
                temper,
                wvnmlo,
                wvnmhi,
                ntau,
                utau,
                numu,
                umu,
                nphi,
                phi,
                albedo,
                btemp,
                ttemp,
                temis,
                tauc,
                nlev_c,
                zd_c,
                nzout,
                zout);
  }

  /*    check input dimensions and variables */
  chekin_tzs (nlyr,
              maxnlyr,
              dtauc,
              nlev_c,
              zd_c,
              nzout,
              zout,
              ssalb,
              temper,
              wvnmlo,
              wvnmhi,
              usrtau,
              ntau,
              utau,
              usrang,
              numu,
              umu,
              nphi,
              phi,
              albedo,
              btemp,
              ttemp,
              temis,
              planck,
              tauc,
              quiet);

  for (lu = 0; lu < ntau; lu++) {
    //    fprintf(stderr,"lu = %4d  nlayers[lu] = %4d   nlyr = %4d\n",lu,nlayers[lu],nlyr);
    /* check layer bottom level */
    if (nlayers[lu] + 1 == nlyr) {
      /*          real surface albedo = 1 - emissivity                                         */
      alb[lu] = albedo;
      /*          Planck radiation from surface with emissivity=1-ALBEDO and temperature BTEMP */
      bplanck[lu] = (1 - alb[lu]) * c_planck_func1 (wvnmlo, wvnmhi, btemp);
      //      fprintf(stderr,"%4d  %4d  %12.10f  %12.9f  %12.8f\n",lu, nlayers[lu]+1,alb[lu],bplanck[lu],btemp);
    } else {
      /*          Real surface albedo = 1 - emissivity ==> emissivity=1 (black cloud)          */
      alb[lu] = 0.0;
      /*          Planck radiation of bottom level with emissivity=1 (black cloud)             */
      /*          and temperature temper( nlayers[lu]+1 ) of bottom level                      */
      bplanck[lu] = (float)pkag[nlayers[lu] + 1];
    }
  }

  if (prnt[0]) {
    fprintf (stderr, "\n%4d %s", ntau, " User optical depths :");
    for (lu = 0; lu < ntau; lu++)
      fprintf (stderr, "%14.8f", utau[lu]);
    fprintf (stderr, "\n");

    fprintf (stderr, "\n%4d %s", ntau, " User level numbers :");
    for (lu = 0; lu < ntau; lu++)
      fprintf (stderr, "%4d", nlayers[lu] + 1);
    fprintf (stderr, "\n");

    fprintf (stderr, "\n%4d %s", ntau, " User levels / km :");
    for (lu = 0; lu < ntau; lu++)
      fprintf (stderr, "%14.8f", zout[lu]);
    fprintf (stderr, "\n");

    fprintf (stderr, "\n%4d %s", ntau, " User levels / km :");
    for (lu = 0; lu < ntau; lu++)
      fprintf (stderr, "%14.8f", zd_c[nlayers[lu] + 1]);
    fprintf (stderr, "\n");
  }

  /* loop from TOA to lowest layer/level needed for calculations */
  for (lc = 0; lc < maxnlyr; lc++) {
    /*          assume linear variation of Planck radiation from one level                     */
    /*          to the next one: Planck=b0+tauc*b1                                             */
    b1[lc] = 0.;

    /* CE: this threshold was commented out, but it seems that it is needed because otherwise we get
     negative emissions in some layers */
    if (dtauc[lc] > 1.e-4)
      b1[lc] = (pkag[lc + 1] - pkag[lc]) / dtauc[lc];
    /* In case of thin atmosphere but high temperature variation from one level to the next, dtauc is very small  */
    /* but Planck varies sensibly, so we have to take care. The threshold 1.e10 is not crucial because thin       */
    /* atmospheric layers emit few radiation.                                                                     */
    if (fabs (b1[lc]) > 1.e10) {
      /* set b0 to average Planck emission, also not crucial */
      b0[lc] = 0.5 * (pkag[lc + 1] + pkag[lc]);
      b1[lc] = 0.0;
    } else {
      b0[lc] = pkag[lc] - tauc[lc] * b1[lc];
    }

    /* Compute coefficients like this even if coeff1[lc]=pkag[lc+1] and coeff2[lc]=pkag[lc]    */
    coeff1[lc] = b0[lc] + tauc[lc + 1] * b1[lc];
    coeff2[lc] = b0[lc] + tauc[lc] * b1[lc];

    //    fprintf(stderr," LC %d Deltatau %e DIFFCOEFF1 %e COEFF1 %e pkag[lc+1] = %e\n",lc,dtauc[lc],coeff1[lc]-pkag[lc+1],coeff1[lc],pkag[lc+1]);
    //    fprintf(stderr," LC %d Deltatau %e DIFFCOEFF2 %e COEFF2 %e pkag[lc]   = %e\n",lc,dtauc[lc],coeff2[lc]-pkag[lc],coeff2[lc],pkag[lc]);
  }

  /*    Loop over UMUs    */
  for (iu = 0; iu < numu; iu++) {

    /*       Double precision UMU  */
    mu = (double)umu[iu];

    //    /* loop from TOA to lowest layer/level needed for calculations */
    //    for (lc=0;lc<=maxnlyr;lc++){
    //      /*          Initialise ETAU = Attenuation factors from TOA to levels LC in direction MU  */
    //      etau[lc]=exp(-tauc[lc]/mu);
    //    }

    /*       Emission of atmosphere:   */
    /*       ***********************   */

    /*          Initialise ETAU = Attenuation factors from TOA to level LC in direction MU   */
    lc       = 0;
    etau[lc] = exp (-tauc[lc] / mu);

    /* loop from TOA to lowest layer/level needed for calculations */

    for (lc = 0; lc < maxnlyr; lc++) {

      /*          Initialise ETAU = Attenuation factors from TOA to level LC+1 in direction MU */
      etau[lc + 1] = exp (-tauc[lc + 1] / mu);

      /*          Attenuation factors from TOA to levels LC+1 resp. LC:   */
      /*           ETAU[lc+1]                                             */
      /*           ETAU[lc]                                               */

      /*          Radiation emitted from layer LC (between levels LC and LC+1) towards TOA   */
      radem[lc] = etau[lc] * coeff2[lc] - etau[lc + 1] * coeff1[lc] + b1[lc] * mu * (etau[lc] - etau[lc + 1]);
    }

    /*       Loop over user defined tau = black cloud top heights   */
    for (lu = 0; lu < ntau; lu++) {

      //      fprintf(stderr,"lu = %4d  %4d\n",lu,nlayers[lu]);

      /*          Initialise RAD to 0.              */
      rad[lu] = 0.;
      /*          Initialise TAU1                   */
      /*          to optical thickness from TOA up to bottom of layer nlayers[lu] */
      /*           TAU1  = TAUC(NLYR)               */
      tau1 = tauc[nlayers[lu] + 1];
      /*          Initialise ETAU1                  */
      etau1    = exp (-tau1);
      tauetau1 = tau1 * etau1;
      /*          Initialise EI1                    */

#if HAVE_LIBGSL
      /* exponential integral does not converge for values larger about 600*/
      /* the integral for this case becomes 0*/
      if (tau1 > 500.0)
        ei1 = 0.0;
      else
        ei1 = gsl_sf_expint_Ei (-tau1);
        /* fprintf(stderr, "lu %d  tau1 %g ei1 %g\n",lu,  tau1, ei1);  */
#endif
#if !HAVE_LIBGSL
      fprintf (stderr, "libRadtran was built without the gsl-library. Thus \n");
      fprintf (stderr, "tzs may not be used. Please install gsl on your system\n");
      fprintf (stderr, "and rebuild libRadtran.\n");
      exit (0);
#endif

      /*          Initialise SQTAU1                 */
      sqtau1 = tau1 * tau1;
      /*          Initialise CUTAU1                 */
      cutau1 = tau1 * sqtau1;
      /*          initialise RADREFL11              */
      radrefl11 = -2.0 * etau1 + tauetau1 - tau1 * tauetau1 - cutau1 * ei1;

      /*          Attenuation factor due to full atmospheric absorption in direction MU  */
      etaumu = etau[nlayers[lu] + 1];

      /* total atmospheric emission i.e. intergral from TOA to bottom */
      rademtot = 0.0;

      /*          Contribution from surface + atmosphere              */
      for (lc = 0; lc <= nlayers[lu]; lc++) {

        if (alb[lu] > 0.0) {

          /*             Reflection from surface:   */
          /*             ************************   */
          tau2      = tau1;
          etau2     = etau1;
          ei2       = ei1;
          sqtau2    = sqtau1;
          radrefl22 = radrefl11;
          radrefl2  = etau2 * coeff2[lc] * (1. - tau2) - coeff2[lc] * sqtau2 * ei2;

          if (lc == nlayers[lu]) {
            radrefl1  = coeff1[lc];
            radrefl11 = -2.0;
          } else {
            tau1     = tauc[nlayers[lu] + 1] - tauc[lc + 1];
            etau1    = exp (-tau1);
            tauetau1 = tau1 * etau1;

#if HAVE_LIBGSL
            /* exponential integral does not converge for values larger about 600*/
            /* the integral for this case becomes 0*/
            if (tau1 > 500.0)
              ei1 = 0.0;
            else
              ei1 = gsl_sf_expint_Ei (-tau1);
#endif

            sqtau1 = tau1 * tau1;
            cutau1 = tau1 * sqtau1;

            radrefl11 = -2.0 * etau1 + tauetau1 - tau1 * tauetau1 - cutau1 * ei1;

            radrefl1 = etau1 * coeff1[lc] * (1. - tau1) - coeff1[lc] * sqtau1 * ei1;
          }

          radrefl = radrefl1 - radrefl2;

          /* fprintf(stderr,"%s %4d %4d %13.6e\n","LC NLAYERS RADREFL[lc]",lc,nlayers[lu],radrefl);  */
          /* fprintf(stderr,"%s %4d %13.6e\n","LC RADREFL1",lc,radrefl1);                          */
          /* fprintf(stderr,"%s %4d %13.6e\n","LC RADREFL2",lc,radrefl2);                          */
          /* fprintf(stderr,"%s %4d %13.6e %13.6e\n","LC COEFF1 COEFF2",lc,coeff1[lc],coeff2[lc]); */
          /* fprintf(stderr,"%s %4d %13.6e %13.6e %13.6e %13.6e\n","LC TAU1 TAU2 SQTAU1 SQTAU2",lc,tau1,tau2,sqtau1,sqtau2); */
          /* fprintf(stderr,"%s %4d %13.6e %13.6e\n","LC EI1 EI2",lc,ei1,ei2); */

          radrefl += b1[lc] / 3. * (radrefl11 - radrefl22);

          //	  fprintf(stderr,"%s %4d %13.6e\n","LC RADREFL[lc]",lc,radrefl);
          //	  fprintf(stderr,"%s %4d %13.6e\n","LC B1[lc]",lc,b1[lc]);
          //	  fprintf(stderr,"%s %4d %13.6e\n","LC RADREFL11",lc,radrefl11);
          //	  fprintf(stderr,"%s %4d %13.6e\n","LC RADREFL22",lc,radrefl22);

          /*              RAD[lu]   = RAD[lu]+RADEM[lc]+ALBEDO*ETAUMU*RADREFL      */
          rad[lu] += radrefl;
        }

        /* emission of atmosphere */
        rademtot += radem[lc];
      }

      /* fprintf(stderr,"lu = %4d   rad[lu] refl = %13.6e\n",lu,rad[lu]); */

      /* after computation of radiance reflected at the surface we multiply by the albedo and the extinction from bottom to TOA */
      /*              RAD[lu]   = RAD[lu]+RADEM[lc]+ALBEDO*ETAUMU*RADREFL      */
      rad[lu] *= alb[lu] * etaumu;

      /* fprintf(stderr,"lu = %4d   rad[lu] refl attenuated = %13.6e\n",lu,rad[lu]); */

      /*          Contribution from atmospheric emission of all levels           */
      rad[lu] += rademtot;

      /* fprintf(stderr,"lu = %4d   rad[lu] refl attenuated + emission = %13.6e\n",lu,rad[lu]); */

      /*          contribution from surface emission (bottom):     */
      /*          ********************************************     */
      rad[lu] += bplanck[lu] * etaumu;

      /* fprintf(stderr,"lu = %4d   rad[lu] refl attenuated + emission + surface = %13.6e\n",lu,rad[lu]); */

      /* isotropic radiance                               */
      for (j = 0; j < nphi; j++) {
        uu[j][lu][iu] = (float)rad[lu];
      }
      uavg[lu] = (float)rad[lu];
      /*       End loop over user tau          */
    }

    /*    End loop over user umu    */
  }

  free (radem);
  free (tauc);
  free (pkag);
  free (etau);
  free (b0);
  free (b1);
  free (coeff1);
  free (coeff2);

  return 0;
}

/*============================= errmsg_tzs() ===============================*/

/*
 * Print out a warning or error message;  abort if type == TZS_ERROR
 */

#define MAX_WARNINGS 100

void errmsg_tzs (char* messag, int type) {
  static int warning_limit = FALSE, num_warnings = 0;

  if (type == TZS_ERROR) {
    fprintf (stderr, "\n ******* ERROR >>>>>>  %s\n", messag);
    exit (1);
  }

  if (warning_limit)
    return;

  if (++num_warnings <= MAX_WARNINGS) {
    fprintf (stderr, "\n ******* WARNING >>>>>>  %s\n", messag);
  } else {
    fprintf (stderr, "\n\n >>>>>>  TOO MANY WARNING MESSAGES --  They will no longer be printed  <<<<<<<\n\n");
    warning_limit = TRUE;
  }

  return;
}

#undef MAX_WARNINGS

void prtinp_tzs (int     nlyr,
                 int     maxnlyr,
                 float*  dtauc,
                 float*  ssalb,
                 float*  temper,
                 float   wvnmlo,
                 float   wvnmhi,
                 int     ntau,
                 float*  utau,
                 int     numu,
                 float*  umu,
                 int     nphi,
                 float*  phi,
                 float   albedo,
                 double  btemp,
                 float   ttemp,
                 float   temis,
                 double* tauc,
                 int     nlev_c,
                 float*  zd_c,
                 int     nzout,
                 float*  zout) {
  /*     print values of input variables                                 */
  /*                                                                     */
  /*   Called by- TZS                                                    */
  /* --------------------------------------------------------------------*/

  /*     .. Local Scalars ..  */
  int lc, lu;

  fprintf (stderr, "\n%s%4d\n", " No. computational layers        = ", nlyr);
  fprintf (stderr, "\n%s%4d\n", " No. computational layers needed = ", maxnlyr);

  fprintf (stderr, "\n%4d %s", ntau, " User optical depths :");
  for (lu = 0; lu < ntau; lu++)
    fprintf (stderr, "%13.6e", utau[lu]);
  fprintf (stderr, "\n");

  fprintf (stderr, "\n%4d %s", numu, " User polar angle cosines :");
  for (lu = 0; lu < numu; lu++)
    fprintf (stderr, "%9.5f", umu[lu]);
  fprintf (stderr, "\n");

  fprintf (stderr, "\n%4d %s", nphi, " User azimuthal angles :");
  for (lu = 0; lu < nphi; lu++)
    fprintf (stderr, "%9.2f", phi[lu]);
  fprintf (stderr, "\n");

  fprintf (stderr,
           "%s%14.4f%14.4f\n%s%10.2f%s%8.4f\n%s%8.4f\n",
           "    Thermal emission in wavenumber interval :",
           wvnmlo,
           wvnmhi,
           "    Bottom temperature = ",
           btemp,
           "    Top temperature = ",
           ttemp,
           "    Top emissivity = ",
           temis);

  fprintf (stderr, "%s%8.4f\n", "    Bottom albedo (Lambertian) =", albedo);

  /* Print layer variables */

  fprintf (stderr,
           "\n%s\n%s\n%s\n",
           "            Layer        Total        Single",
           "          Optical      Optical    Scattering",
           "            Depth        Depth        Albedo   Temperature");

  for (lc = 0; lc < maxnlyr; lc++) {
    fprintf (stderr, "%4d%13.6e%13.6e%14.3f%14.3f\n", lc, dtauc[lc], tauc[lc], ssalb[lc], temper[lc]);
  }

  fprintf (stderr, "%4d%13.6e%13.6e%14.3f%14.3f\n", maxnlyr, 0.0, tauc[maxnlyr], 0.0, temper[maxnlyr]);

  fprintf (stderr,
           "\n%s\n%s\n%s\n",
           "            Layer        Single",
           "          Optical    Scattering",
           "            Depth        Albedo   Temperature");

  for (lc = 0; lc < nlyr; lc++) {
    fprintf (stderr, "%4d%13.6e%14.3f%14.3f\n", lc, dtauc[lc], ssalb[lc], temper[lc]);
  }

  fprintf (stderr, "%4d%13.6e%14.3f%14.3f\n", nlyr, -1.0, -1.0, temper[nlyr]);

  fprintf (stderr, "\n%s\n%s\n%s\n", "      Atmospheric", "            Level", "             [km]");

  for (lc = 0; lc < nlev_c; lc++) {
    fprintf (stderr, "%4d     %8.3f\n", lc, zd_c[lc]);
  }

  fprintf (stderr, "\n%s\n%s\n%s\n", "           Output", "            Level", "             [km]");

  for (lc = 0; lc < nzout; lc++) {
    fprintf (stderr, "%4d     %8.3f\n", lc, zout[lc]);
  }
}

void chekin_tzs (int     nlyr,
                 int     maxnlyr,
                 float*  dtauc,
                 int     nlev_c,
                 float*  zd_c,
                 int     nzout,
                 float*  zout,
                 float*  ssalb,
                 float*  temper,
                 float   wvnmlo,
                 float   wvnmhi,
                 int     usrtau,
                 int     ntau,
                 float*  utau,
                 int     usrang,
                 int     numu,
                 float*  umu,
                 int     nphi,
                 float*  phi,
                 float   albedo,
                 double  btemp,
                 float   ttemp,
                 float   temis,
                 int     planck,
                 double* tauc,
                 int     quiet) {
  /*     Checks the input dimensions and variables  */
  int inperr = FALSE;
  int lc, j, iu, lu;

  if (planck == FALSE)
    inperr = c_write_bad_var (VERBOSE, "planck");

  if (nlyr < 1)
    inperr = c_write_bad_var (VERBOSE, "NLYR");

  if (maxnlyr < 1)
    inperr = c_write_bad_var (VERBOSE, "MAXNLYR");

  if (nlyr + 1 != nlev_c)
    inperr = c_write_bad_var (VERBOSE, "NLYR/NLEV_C");

  for (lc = 0; lc < nlyr; lc++) {
    if (dtauc[lc] < 0.) {
      inperr = c_write_bad_var (VERBOSE, "dtauc");
    }
    if (zd_c[lc] < 0.) {
      inperr = c_write_bad_var (VERBOSE, "zd");
    }
    if (ssalb[lc] < 0.0 || ssalb[lc] > 1.0) {
      inperr = c_write_bad_var (VERBOSE, "ssalb");
    }
    if (temper[lc] < 0.) {
      inperr = c_write_bad_var (VERBOSE, "temper");
    }
  }

  if (usrtau) {
    if (ntau < 1) {
      inperr = c_write_bad_var (VERBOSE, "ntau");
    }
    if (nzout < ntau)
      inperr = c_write_bad_var (VERBOSE, "ntau/nzout");
    for (lu = 0; lu < ntau; lu++) {
      /* Do a relative check to see if we are just beyond the bottom boundary */
      /* This might happen due to numerical rounding off problems.  ak20110224*/

      //      fprintf(stderr,"lu = %d  utau = %13.6e  tauc[maxnlyr] = %13.6e  zout = %f\n",lu,utau[lu],tauc[maxnlyr],zout[lu]);
      if (fabs (utau[lu] - tauc[maxnlyr]) <= 1.e-6 * tauc[maxnlyr]) {
        utau[lu] = tauc[maxnlyr];
      }
      //      fprintf(stderr,"lu = %d  utau = %13.6e  tauc[maxnlyr] = %13.6e  zout = %f\n",lu,utau[lu],tauc[maxnlyr],zout[lu]);

      //      if(utau[lu] < 0. || utau[lu] > tauc[maxnlyr]) {
      if (utau[lu] < 0. || (utau[lu] - tauc[maxnlyr]) > 1.e-6 * tauc[maxnlyr]) {
        fprintf (stderr,
                 "WARNING: lu = %d  utau = %13.6e  tauc[maxnlyr] = %13.6e  zout = %f\n",
                 lu,
                 utau[lu],
                 tauc[maxnlyr],
                 zout[lu]);
        inperr = c_write_bad_var (VERBOSE, "utau");
      }
      if (zout[lu] < 0.0)
        inperr = c_write_bad_var (VERBOSE, "zout");
    }
  } else {
    ntau    = 1;
    utau    = (float*)calloc (1, sizeof (float));
    utau[0] = 0.;
    nzout   = 1;
    zout    = (float*)calloc (1, sizeof (float));
    zout[0] = zd_c[0];
  }

  if (usrang) {
    if (numu < 0) {
      inperr = c_write_bad_var (VERBOSE, "numu");
    }
    for (iu = 0; iu < numu; iu++) {
      if (umu[iu] < -1. || umu[iu] > 1. || umu[iu] == 0.) {
        inperr = c_write_bad_var (VERBOSE, "umu");
      }
      if (iu >= 1) {
        if (umu[iu] < umu[iu - 1]) {
          inperr = c_write_bad_var (VERBOSE, "umu");
        }
      }
    }
    if (nphi <= 0)
      inperr = c_write_bad_var (VERBOSE, "nphi");

    for (j = 0; j < nphi; j++) {
      if (phi[j] < 0. || phi[j] > 360.) {
        inperr = c_write_bad_var (VERBOSE, "phi");
      }
    }
  } else {
    numu   = 1;
    umu    = (float*)calloc (1, sizeof (float));
    umu[0] = 1.;
    nphi   = 1;
    phi    = (float*)calloc (1, sizeof (float));
    phi[0] = 0.;
  }

  if (albedo < 0. || albedo > 1.)
    inperr = c_write_bad_var (VERBOSE, "albedo");

  if (wvnmlo < 0.0 || wvnmhi <= wvnmlo)
    inperr = c_write_bad_var (VERBOSE, "wvnmlo,hi");

  if (temis < 0. || temis > 1.)
    inperr = c_write_bad_var (VERBOSE, "temis");

  if (btemp < 0.0)
    inperr = c_write_bad_var (VERBOSE, "btemp");

  if (ttemp < 0.0)
    inperr = c_write_bad_var (VERBOSE, "ttemp");

  if (inperr)
    c_errmsg ("TZS--input and/or dimension errors", TZS_ERROR);

  if (!quiet) {
    for (lc = 1; lc <= nlyr; lc++) {
      if (fabs (temper[lc] - temper[lc - 1]) > 10.) {
        c_errmsg ("check_inputs--vertical temperature step may be too large for good accuracy", TZS_WARNING);
      }
    }
  }
}


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "Data\Types.h"
#include "MageSettings.h"
#include "Image\ORBDescriptor.h"
#include "BaseFeatureMatcher.h"

#include <vector>
#include <gsl\gsl>
#include <memory>

namespace mage
{
    class Keyframe;
    class AnalyzedImage;

    class BaseBow
    {
    public:
        struct QueryMatch
        {
            Id<Keyframe> Id;
            float Score;

            operator mage::Id<Keyframe>() const
            {
                return Id;
            }
        };

        BaseBow(const BagOfWordsSettings& settings)
            :m_settings{ settings }
        {}

        virtual ~BaseBow() {}

        virtual void AddTrainingDescriptors(gsl::span<const ORBDescriptor>) {}

        virtual void AddImage(const Id<Keyframe>, const AnalyzedImage& image) = 0;

        virtual void RemoveImage(const Id<Keyframe>) = 0;

        virtual size_t QueryFeatures(const ORBDescriptor& descriptor, const Id<Keyframe>& keyframe, std::vector<ptrdiff_t>& features) const = 0;

        virtual std::unique_ptr<BaseFeatureMatcher> CreateFeatureMatcher(const Id<Keyframe>& id, gsl::span<const ORBDescriptor> features) const = 0;

        virtual std::unique_ptr<BaseBow> CreateTemporaryBow() const = 0;

        virtual std::vector<QueryMatch> QueryUnknownImage(gsl::span<const ORBDescriptor> descriptors, size_t maxResults) const = 0;

        virtual void Clear() = 0;

        virtual bool IsTrainingDone() const { return true;}

    protected:
        const BagOfWordsSettings& m_settings;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam