<commit_before>#include <math.h>
#include "allvars.h"
#include "proto.h"
#include <map>

static double R8;
static double r_tophat;

static double AA, BB, CC;
static double nu;
static double Norm;

/*This stores the conversion between the tables and the internal units. By default 1e-3*/
static double kctog;
/*This prints the CAMB transfer function*/
double tk_CAMB(double k, int Type);
/*Power spectra*/
double PowerSpec_CAMB(double k, int Type);
static int NPowerTable;
#define APRIM 2.43e-9
#define PIVOT_SCALE (0.05*kctog/HubbleParam)

/*Structure for matter power table*/
static struct pow_table
{
  double logk, logD,logDb;
  double logD2nd;
  double logDtot;
}
 *PowerTable;

static struct pow_matter
{ 
  double kmat,pmat;
 }
*PowerMatter;

 /*Structure for transfer table*/
struct trans_row{
	double T_CDM;
	double T_b;
	double T_g;
	double T_r;
	double T_n;
	double T_t;
};

static std::map<double, struct trans_row> transfer_tables;

double PowerSpec(double k, int Type)
{
  double power, alpha, Tf;
  switch (WhichSpectrum)
    {
    case 1:
      power = PowerSpec_EH(k);
      break;
    case 2:
      if (Type == 2)
#ifdef NEUTRINOS
        power = PowerSpec_Tabulated2nd(k);
#else
        power = PowerSpec_DM_2ndSpecies(k);
#endif
      else if(Type == 0)
          power = PowerSpec_Tabulated_b(k);
      else
          power = PowerSpec_Tabulated(k);
      break;
#ifdef DIFFERENT_TRANSFER_FUNC
    case 3:
      /*This is my reimplementation of Spectrum 2. */
  /*ADD THE FACTOR OF (2π)^3 to convert from CAMB conventions to GADGET conventions!!*/
      power = PowerSpec_CAMB(k,Type)/pow(2*M_PI,3);
      break;
#endif
    default:
      power = PowerSpec_Efstathiou(k);
      break;
    }

  if(WDM_On == 1)
    {
      /* Eqn. (A9) in Bode, Ostriker & Turok (2001), assuming gX=1.5  */
      alpha =
	0.048 * pow((Omega - OmegaBaryon) / 0.4, 0.15) * pow(HubbleParam / 0.65,
							     1.3) * pow(1.0 / WDM_PartMass_in_kev, 1.15);
      Tf = pow(1 + pow(alpha * k * (3.085678e24 / UnitLength_in_cm), 2 * 1.2), -5.0 / 1.2);
      power *= Tf * Tf;
    }

#if defined(DIFFERENT_TRANSFER_FUNC)
/* Type 100 recomputes normalization over the total matter p(k) */
  if(Type == 100 && WhichSpectrum ==2) 
	power = PowerSpec_TOTAL(k);
#endif 
  if(WhichSpectrum < 2)
    power *= pow(k, PrimordialIndex - 1.0);

  return power/ (Dplus*Dplus);
}


double PowerSpec_DM_2ndSpecies(double k)
{
  /* at the moment, we simply call the Eistenstein & Hu spectrum
   * for the second DM species, but this could be replaced with
   * something more physical, say for neutrinos
   */

  double power;

  power = Norm * k * pow(tk_eh(k), 2);

  return power;
}


/*Read the transfer tables from CAMB*/
void read_transfer_table(void)
{
  /* Transfer file format:  
   * k/h Delta_CDM/k2 Delta_b/k2 Delta_g/k2 Delata_r/k2 Delta_nu/k2 Delta_tot/k2
   *     CDM, baryon,photon,massless neutrino, massive neutrinos, and total (massive)*/
	FILE *trans;
	char buf[200];
	struct trans_row tmp_row;
        double tmp_k;
        kctog=UnitLength_in_cm/InputSpectrum_UnitLength_in_cm;
	/*Open file*/
	sprintf(buf,FileWithTransfer);
	if(!(trans=fopen(buf,"r"))) {
		fprintf(stderr,"Can't open transfer file %s! You probably forgot to create it!\n",buf);
		exit(17);
	} 
	/*Work out how many lines in file*/
	if(!transfer_tables.empty()) {
		fprintf(stderr,"read_CAMB_tables has been called more than once! Failing.");
		exit(19);
	}
	/*Read line by line.*/
	while(fscanf(trans," %lg %lg %lg %lg %lg %lg %lg\n",&tmp_k,&tmp_row.T_CDM,&tmp_row.T_b, &tmp_row.T_g,&tmp_row.T_r,&tmp_row.T_n,&tmp_row.T_t)==7)
	{
		/* k needs to go from (h/Mpc) units to internal Gadget units (h/kpc) by default.
                 * kctog is by default 1e-3 */
		tmp_k *= kctog;
		/*Append line to table.*/
		transfer_tables[tmp_k]=tmp_row;
		if(feof(trans))
			break;
	}
	if(ferror(trans)) {
		fprintf(stderr,"Error reading transfer file: %d",errno);
		exit(21);
	}
        printf("Found %d rows in input CAMB transfer file\n",(int)transfer_tables.size());
	fclose(trans);
	/*The CAMB T_f/k_c is in units of Mpc^2! NOTE NO h!*/
	double tctog=(HubbleParam*HubbleParam)/(kctog*kctog);
        std::map<double, struct trans_row>::iterator it;
	for(it=transfer_tables.begin(); it != transfer_tables.end(); ++it) {
	/*The transfer function should be normalized to about 1 on large scales.*/
		((*it).second).T_CDM *= tctog;
		((*it).second).T_g *= tctog;
		((*it).second).T_b *= tctog;
		((*it).second).T_r *= tctog;
		((*it).second).T_n *= tctog;
		((*it).second).T_t *= tctog;
	}
	return;
}

void read_power_table(void)
{
  FILE *fd;
  char buf[500];
  double k;

  sprintf(buf, FileWithTransfer);
  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input TRANSFER in file '%s'\n", buf);
      exit(17);
    }

  NPowerTable = 0;
  do
    {
      double T_cdm, T_b, dummy, T_nu, T_tot;

      /* read transfer function file from CAMB */
      if(fscanf(fd, " %lg %lg %lg %lg %lg %lg %lg", &k, &T_cdm, &T_b, &dummy, &dummy, &T_nu, &T_tot) == 7)
	NPowerTable++;
      else
	break;
    }
  while(1);
  fclose(fd);
      printf("found %d rows in input TRANSFER table\n", NPowerTable);
      fflush(stdout);



sprintf(buf, FileWithInputSpectrum);
  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input SPECTRUM in file '%s'\n", buf);
      exit(17);
    }
  NPowerTable = 0;
  do
    {
      double ktab, Pktab;
      /* read TOTAL matter power spectrum from CAMB*/
      if(fscanf(fd, " %lg %lg ", &ktab, &Pktab) == 2)
	NPowerTable++;
      else
	break;
    }
  while(1);
  fclose(fd);
  printf("found %d rows in input SPECTRUM table\n", NPowerTable);
  fflush(stdout);
  PowerTable = (pow_table *) malloc(NPowerTable * sizeof(struct pow_table));  
  PowerMatter = (pow_matter *) malloc(NPowerTable * sizeof(struct pow_matter)); 
  /* define matter array */
  sprintf(buf, FileWithInputSpectrum);
  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input SPECTRUM in file '%s'\n", buf);
      exit(18);
    }

  NPowerTable = 0;
  do
    {
      double kmat, pmat;
      if(fscanf(fd, " %lg %lg", &kmat, &pmat) == 2)
	{
	  PowerMatter[NPowerTable].kmat = kmat;
	  PowerMatter[NPowerTable].pmat = pmat;
	  NPowerTable++;
	}
      else
	break;
    }
  while(1);

  fclose(fd);

  qsort(PowerMatter, NPowerTable, sizeof(struct pow_matter), compare_logk);

   sprintf(buf, FileWithTransfer);
   if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input spectrum in file '%s'\n", buf);
      exit(18);
    }

  NPowerTable = 0;
  do
    {
      double T_cdmnew, T_b, dummy, T_nu, T_tot, T_cdm;
      double delta_cdm, delta_nu, delta_tot, delta_b;

      if(fscanf(fd, " %lg %lg %lg %lg %lg %lg %lg", &k, &T_cdm, &T_b, &dummy, &dummy, &T_nu, &T_tot) == 7)
	{
	  PowerTable[NPowerTable].logk = log10(k);

          /* The dark matter may incorporate other particle types as well, eg,
           * fake neutrinos, or baryons.
           * NOTE that CAMB defines T_tot = (T_CDM M_CDM + T_b M_b +T_nu M_nu) / (M_CDM + M_b +M_nu)
           * HOWEVER, when relativistic (ie, in the early universe), neutrinos will redshift
           * like radiation. The CAMB transfer function takes this into account when calculating T_tot,
           * so instead of summing the transfer functions, use the total transfer function and
           * optionally subtract the baryons.*/
          if(neutrinos_ks){
                T_cdm = T_tot;
                /*If we have separate gas particles, subtract
                 * the baryon transfer function */
                if(!no_gas){
                    T_cdm -=T_b*OmegaBaryon/Omega;
                }
          }
          /*Add baryons to the CDM if there are no gas particles*/
          else if(no_gas){
                T_cdm = (T_cdm*(Omega-OmegaDM_2ndSpecies - OmegaBaryon) + T_b*OmegaBaryon)/(Omega-OmegaDM_2ndSpecies);
          }
          /*This should be equivalent to the above in almost all cases,
           * but perhaps we want to see the effect of changing only the ICs in CAMB for the neutrinos.*/
          if(no_gas && OmegaDM_2ndSpecies == 0){
                  T_cdm = T_tot;
          }
	  /* obtain P(k) from transfer function ratios like suggested by JL*/
	  /*NOTE for this to work CAMB's transfer_k_interp_matterpower must be off!!*/
	  delta_b   = k * k * k * pow(T_b/T_tot,2)* PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI);
	  delta_cdm = k * k * k * pow(T_cdm/T_tot,2)* PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI);
	  delta_nu = k * k * k * pow(T_nu/T_tot,2) * PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI); 
	  delta_tot = k * k * k * PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI);

	  PowerTable[NPowerTable].logD = log10(delta_cdm);
	  // printf("NT,d_cdm,log10(d_cdm),k %d %g %g %g \n",NPowerTable,delta_cdm,log10(delta_cdm),k);
	  PowerTable[NPowerTable].logD2nd = log10(delta_nu);
          PowerTable[NPowerTable].logDtot = log10(delta_tot);
	  PowerTable[NPowerTable].logDb = log10(delta_b);

	  NPowerTable++;
	}
      else
	break;
    }
  while(1);<commit_msg>Print an error if the k from the matterpower file differs from that in the transfer file<commit_after>#include <math.h>
#include "allvars.h"
#include "proto.h"
#include <map>

static double R8;
static double r_tophat;

static double AA, BB, CC;
static double nu;
static double Norm;

/*This stores the conversion between the tables and the internal units. By default 1e-3*/
static double kctog;
/*This prints the CAMB transfer function*/
double tk_CAMB(double k, int Type);
/*Power spectra*/
double PowerSpec_CAMB(double k, int Type);
static int NPowerTable;
#define APRIM 2.43e-9
#define PIVOT_SCALE (0.05*kctog/HubbleParam)

/*Structure for matter power table*/
static struct pow_table
{
  double logk, logD,logDb;
  double logD2nd;
  double logDtot;
}
 *PowerTable;

static struct pow_matter
{ 
  double kmat,pmat;
 }
*PowerMatter;

 /*Structure for transfer table*/
struct trans_row{
	double T_CDM;
	double T_b;
	double T_g;
	double T_r;
	double T_n;
	double T_t;
};

static std::map<double, struct trans_row> transfer_tables;

double PowerSpec(double k, int Type)
{
  double power, alpha, Tf;
  switch (WhichSpectrum)
    {
    case 1:
      power = PowerSpec_EH(k);
      break;
    case 2:
      if (Type == 2)
#ifdef NEUTRINOS
        power = PowerSpec_Tabulated2nd(k);
#else
        power = PowerSpec_DM_2ndSpecies(k);
#endif
      else if(Type == 0)
          power = PowerSpec_Tabulated_b(k);
      else
          power = PowerSpec_Tabulated(k);
      break;
#ifdef DIFFERENT_TRANSFER_FUNC
    case 3:
      /*This is my reimplementation of Spectrum 2. */
  /*ADD THE FACTOR OF (2π)^3 to convert from CAMB conventions to GADGET conventions!!*/
      power = PowerSpec_CAMB(k,Type)/pow(2*M_PI,3);
      break;
#endif
    default:
      power = PowerSpec_Efstathiou(k);
      break;
    }

  if(WDM_On == 1)
    {
      /* Eqn. (A9) in Bode, Ostriker & Turok (2001), assuming gX=1.5  */
      alpha =
	0.048 * pow((Omega - OmegaBaryon) / 0.4, 0.15) * pow(HubbleParam / 0.65,
							     1.3) * pow(1.0 / WDM_PartMass_in_kev, 1.15);
      Tf = pow(1 + pow(alpha * k * (3.085678e24 / UnitLength_in_cm), 2 * 1.2), -5.0 / 1.2);
      power *= Tf * Tf;
    }

#if defined(DIFFERENT_TRANSFER_FUNC)
/* Type 100 recomputes normalization over the total matter p(k) */
  if(Type == 100 && WhichSpectrum ==2) 
	power = PowerSpec_TOTAL(k);
#endif 
  if(WhichSpectrum < 2)
    power *= pow(k, PrimordialIndex - 1.0);

  return power/ (Dplus*Dplus);
}


double PowerSpec_DM_2ndSpecies(double k)
{
  /* at the moment, we simply call the Eistenstein & Hu spectrum
   * for the second DM species, but this could be replaced with
   * something more physical, say for neutrinos
   */

  double power;

  power = Norm * k * pow(tk_eh(k), 2);

  return power;
}


/*Read the transfer tables from CAMB*/
void read_transfer_table(void)
{
  /* Transfer file format:  
   * k/h Delta_CDM/k2 Delta_b/k2 Delta_g/k2 Delata_r/k2 Delta_nu/k2 Delta_tot/k2
   *     CDM, baryon,photon,massless neutrino, massive neutrinos, and total (massive)*/
	FILE *trans;
	char buf[200];
	struct trans_row tmp_row;
        double tmp_k;
        kctog=UnitLength_in_cm/InputSpectrum_UnitLength_in_cm;
	/*Open file*/
	sprintf(buf,FileWithTransfer);
	if(!(trans=fopen(buf,"r"))) {
		fprintf(stderr,"Can't open transfer file %s! You probably forgot to create it!\n",buf);
		exit(17);
	} 
	/*Work out how many lines in file*/
	if(!transfer_tables.empty()) {
		fprintf(stderr,"read_CAMB_tables has been called more than once! Failing.");
		exit(19);
	}
	/*Read line by line.*/
	while(fscanf(trans," %lg %lg %lg %lg %lg %lg %lg\n",&tmp_k,&tmp_row.T_CDM,&tmp_row.T_b, &tmp_row.T_g,&tmp_row.T_r,&tmp_row.T_n,&tmp_row.T_t)==7)
	{
		/* k needs to go from (h/Mpc) units to internal Gadget units (h/kpc) by default.
                 * kctog is by default 1e-3 */
		tmp_k *= kctog;
		/*Append line to table.*/
		transfer_tables[tmp_k]=tmp_row;
		if(feof(trans))
			break;
	}
	if(ferror(trans)) {
		fprintf(stderr,"Error reading transfer file: %d",errno);
		exit(21);
	}
        printf("Found %d rows in input CAMB transfer file\n",(int)transfer_tables.size());
	fclose(trans);
	/*The CAMB T_f/k_c is in units of Mpc^2! NOTE NO h!*/
	double tctog=(HubbleParam*HubbleParam)/(kctog*kctog);
        std::map<double, struct trans_row>::iterator it;
	for(it=transfer_tables.begin(); it != transfer_tables.end(); ++it) {
	/*The transfer function should be normalized to about 1 on large scales.*/
		((*it).second).T_CDM *= tctog;
		((*it).second).T_g *= tctog;
		((*it).second).T_b *= tctog;
		((*it).second).T_r *= tctog;
		((*it).second).T_n *= tctog;
		((*it).second).T_t *= tctog;
	}
	return;
}

void read_power_table(void)
{
  FILE *fd;
  char buf[500];
  double k;

  sprintf(buf, FileWithTransfer);
  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input TRANSFER in file '%s'\n", buf);
      exit(17);
    }

  NPowerTable = 0;
  do
    {
      double T_cdm, T_b, dummy, T_nu, T_tot;

      /* read transfer function file from CAMB */
      if(fscanf(fd, " %lg %lg %lg %lg %lg %lg %lg", &k, &T_cdm, &T_b, &dummy, &dummy, &T_nu, &T_tot) == 7)
	NPowerTable++;
      else
	break;
    }
  while(1);
  fclose(fd);
      printf("found %d rows in input TRANSFER table\n", NPowerTable);
      fflush(stdout);



sprintf(buf, FileWithInputSpectrum);
  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input SPECTRUM in file '%s'\n", buf);
      exit(17);
    }
  NPowerTable = 0;
  do
    {
      double ktab, Pktab;
      /* read TOTAL matter power spectrum from CAMB*/
      if(fscanf(fd, " %lg %lg ", &ktab, &Pktab) == 2)
	NPowerTable++;
      else
	break;
    }
  while(1);
  fclose(fd);
  printf("found %d rows in input SPECTRUM table\n", NPowerTable);
  fflush(stdout);
  PowerTable = (pow_table *) malloc(NPowerTable * sizeof(struct pow_table));  
  PowerMatter = (pow_matter *) malloc(NPowerTable * sizeof(struct pow_matter)); 
  /* define matter array */
  sprintf(buf, FileWithInputSpectrum);
  if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input SPECTRUM in file '%s'\n", buf);
      exit(18);
    }

  NPowerTable = 0;
  do
    {
      double kmat, pmat;
      if(fscanf(fd, " %lg %lg", &kmat, &pmat) == 2)
	{
	  PowerMatter[NPowerTable].kmat = kmat;
	  PowerMatter[NPowerTable].pmat = pmat;
	  NPowerTable++;
	}
      else
	break;
    }
  while(1);

  fclose(fd);

  qsort(PowerMatter, NPowerTable, sizeof(struct pow_matter), compare_logk);

   sprintf(buf, FileWithTransfer);
   if(!(fd = fopen(buf, "r")))
    {
      printf("can't read input spectrum in file '%s'\n", buf);
      exit(18);
    }

  NPowerTable = 0;
  do
    {
      double T_cdmnew, T_b, dummy, T_nu, T_tot, T_cdm;
      double delta_cdm, delta_nu, delta_tot, delta_b;

      if(fscanf(fd, " %lg %lg %lg %lg %lg %lg %lg", &k, &T_cdm, &T_b, &dummy, &dummy, &T_nu, &T_tot) == 7)
	{
          if(fabs(PowerMatter[NPowerTable].kmat - k) > 0.01 *k){
                  fprintf(stderr, "Error: Input spectrum row %d has k=%g, transfer has k=%g.\n",NPowerTable,PowerMatter[NPowerTable].kmat,k);
                  fprintf(stderr, "Remember you need transfer_k_interp_matterpower = F in CAMB\n");
                  exit(47);
          }
	  PowerTable[NPowerTable].logk = log10(k);

          /* The dark matter may incorporate other particle types as well, eg,
           * fake neutrinos, or baryons.
           * NOTE that CAMB defines T_tot = (T_CDM M_CDM + T_b M_b +T_nu M_nu) / (M_CDM + M_b +M_nu)
           * HOWEVER, when relativistic (ie, in the early universe), neutrinos will redshift
           * like radiation. The CAMB transfer function takes this into account when calculating T_tot,
           * so instead of summing the transfer functions, use the total transfer function and
           * optionally subtract the baryons.*/
          if(neutrinos_ks){
                T_cdm = T_tot;
                /*If we have separate gas particles, subtract
                 * the baryon transfer function */
                if(!no_gas){
                    T_cdm -=T_b*OmegaBaryon/Omega;
                }
          }
          /*Add baryons to the CDM if there are no gas particles*/
          else if(no_gas){
                T_cdm = (T_cdm*(Omega-OmegaDM_2ndSpecies - OmegaBaryon) + T_b*OmegaBaryon)/(Omega-OmegaDM_2ndSpecies);
          }
          /*This should be equivalent to the above in almost all cases,
           * but perhaps we want to see the effect of changing only the ICs in CAMB for the neutrinos.*/
          if(no_gas && OmegaDM_2ndSpecies == 0){
                  T_cdm = T_tot;
          }
	  /* obtain P(k) from transfer function ratios like suggested by JL*/
	  /*NOTE for this to work CAMB's transfer_k_interp_matterpower must be off!!*/
	  delta_b   = k * k * k * pow(T_b/T_tot,2)* PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI);
	  delta_cdm = k * k * k * pow(T_cdm/T_tot,2)* PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI);
	  delta_nu = k * k * k * pow(T_nu/T_tot,2) * PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI); 
	  delta_tot = k * k * k * PowerMatter[NPowerTable].pmat/(2*M_PI*M_PI);

          PowerTable[NPowerTable].logD = log10(delta_cdm);
	  // printf("NT,d_cdm,log10(d_cdm),k %d %g %g %g \n",NPowerTable,delta_cdm,log10(delta_cdm),k);
	  PowerTable[NPowerTable].logD2nd = log10(delta_nu);
          PowerTable[NPowerTable].logDtot = log10(delta_tot);
	  PowerTable[NPowerTable].logDb = log10(delta_b);

	  NPowerTable++;
	}
      else
	break;
    }
  while(1);
{'fraction_non_alphanumeric': 0.1356179957771642, 'fraction_numerical': 0.016295815061447675, 'mean_word_length': 2.963098047629264, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5055100', 'n_tokens_mistral': 7252, 'n_tokens_neox': 6824, 'n_words': 1781}