cooling.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, Sep 17, 15:04

cooling.c
View Options

	#include <stdio.h>
	#include <stdlib.h>
	#include <string.h>
	#include <math.h>
	#include <mpi.h>
	#include <gsl/gsl_errno.h>
	#include <gsl/gsl_spline.h>
	#include <gsl/gsl_math.h>
	#include <gsl/gsl_roots.h>

	#include "allvars.h"
	#include "proto.h"


	#ifdef COOLING




	/*! initialize cooling function (the metallicity is fixed)
	*
	* T = temperature
	* L0 = m000 primordial metallicity
	* L1 = m-30
	* L2 = m-20
	* L3 = m-10
	* L4 = m-05
	* L5 = m-00 solar metallicity
	* L6 = m+05
	*/
	int init_cooling(FLOAT metallicity)
	{

	FILE *fd;
	int n,i;
	char line[72];
	float T,L0,L1,L2,L3,L4,L5,L6;
	int MetallicityIndex=4;


	/* find the right index */
	if (All.InitGasMetallicity<-3)
	MetallicityIndex = 0;
	else
	{
	if (All.InitGasMetallicity<-2)
	MetallicityIndex = 1;
	else
	{
	if (All.InitGasMetallicity<-1)
	MetallicityIndex = 2;
	else
	{
	if (All.InitGasMetallicity<-0.5)
	MetallicityIndex = 3;
	else
	{
	if (All.InitGasMetallicity<0)
	MetallicityIndex = 4;
	else
	{
	MetallicityIndex = 5;
	}
	}
	}
	}
	}



	fd = fopen(All.CoolingFile,"r");
	fscanf(fd, "# %6d\n", &n);

	/* allocate memory */
	All.logT = malloc(n * sizeof(double));
	All.logL = malloc(n * sizeof(double));

	/* read empty line */
	fgets(line, sizeof(line), fd);


	/* read file */

	for (i=0;i<n;i++){

	fscanf(fd, "%f %f %f %f %f %f %f %f\n",&T,&L0,&L1,&L2,&L3,&L4,&L5,&L6);

	//printf("%8.3f %8.3f\n",T,L0);

	/* keep only solar values */
	All.logT[i] = (double)T;

	switch (MetallicityIndex)
	{
	case 0:
	All.logL[i] = (double)L0;
	break;
	case 1:
	All.logL[i] = (double)L1;
	break;
	case 2:
	All.logL[i] = (double)L2;
	break;
	case 3:
	All.logL[i] = (double)L3;
	break;
	case 4:
	All.logL[i] = (double)L4;
	break;
	case 5:
	All.logL[i] = (double)L5;
	break;
	case 6:
	All.logL[i] = (double)L6;
	break;
	}
	}

	fclose(fd);

	/* init interpolation */
	All.acc_cooling_spline = gsl_interp_accel_alloc ();
	All.cooling_spline = gsl_spline_alloc (gsl_interp_cspline, n);
	gsl_spline_init (All.cooling_spline, All.logT, All.logL, n);


	#ifdef OUTPUT_COOLING_FUNCTION
	/* test cooling */
	double logT;
	double l;
	logT = 1.;
	while(logT<8)
	{

	T = pow(10,logT);
	l = log10(cooling_function(T));

	if(ThisTask == 0)
	printf("%8.3f %8.3f\n",logT,l);

	logT = logT + 0.05;
	}
	#endif

	return 0;
	}




	/*! This function return the normalized cooling function (no metallicity dependency)
	*/
	double cooling_function(double temperature)
	{

	double logT;


	if (temperature >= All.CutofCoolingTemperature)
	{

	logT = log10(temperature);

	if (logT>8.5)
	logT = 8.5;


	return pow(10,gsl_spline_eval (All.cooling_spline, logT, All.acc_cooling_spline));
	}
	else
	return 1e-100;


	}







	/***************************************************************************

	METALLICITY DEPENDENT COOLING

	**************************************************************************/




	int init_cooling_with_metals()
	{

	/*

	zmin zmax slz
	tmin tmax slt
	FeHSolar
	p k

	*/


	FILE *fd;
	int p,k,i,j;
	float zmin,zmax,slz,tmin,tmax,slt,FeHSolar;
	float lbd;


	if (ThisTask==0)
	{

	fd = fopen(All.CoolingFile,"r");
	fscanf(fd, "%f %f %f\n", &zmin,&zmax,&slz);
	fscanf(fd, "%f %f %f\n", &tmin,&tmax,&slt);
	fscanf(fd, "%f\n" , &FeHSolar);
	fscanf(fd, "%d %d\n" , &p,&k);


	All.CoolingParameters_zmin = zmin;
	All.CoolingParameters_zmax = zmax;
	All.CoolingParameters_slz = slz;
	All.CoolingParameters_tmin = tmin;
	All.CoolingParameters_tmax = tmax;
	All.CoolingParameters_slt = slt;
	//All.CoolingParameters_FeHSolar = FEH_SOLAR; /* instead of FeHSolar*/
	All.CoolingParameters_cooling_data_max = k-2;


	for (i=0;i<p;i++)
	for (j=0;j<k;j++)
	{
	fscanf(fd, "%f\n" ,&lbd);
	All.CoolingParameters_cooling_data[i][j]=lbd;
	}

	fclose(fd);

	}


	/* now broadcast */
	/*
	This is quite bad to do it like this, however, there is no other solution, as
	All has already be sent.
	The other solution will be to create a structure devoted to the cooling (like Cps in chimie.c)
	avoiding to link the parameters to All.
	*/
	MPI_Bcast(&All.CoolingParameters_zmin, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_zmax, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_slz, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_tmin, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_tmax, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_slt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_cooling_data_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&All.CoolingParameters_cooling_data , (COOLING_NMETALICITIES*COOLING_NTEMPERATURES), MPI_DOUBLE, 0, MPI_COMM_WORLD);



	#ifdef OUTPUT_COOLING_FUNCTION
	/* test cooling */
	double logT,T;
	double l;
	double metal;
	logT = 1.;



	metal = (pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar;
	while(logT<8)
	{

	T = pow(10,logT);
	l = log10(cooling_function_with_metals(T,metal));

	if(ThisTask == 0)
	printf("%8.3f %8.3f\n",logT,l);

	logT = logT + 0.05;
	}
	#endif







	return 0;


	}

	/*! This function return the normalized cooling function, that depends on metallicity
	*/
	double cooling_function_with_metals(double temperature,double metal)
	{


	double cooling;
	double T,Z;
	double rt, rz, ft, fz, v1, v2, v;
	int it,iz,itp,izp;

	double zmin,zmax,slz,tmin,tmax,slt,FeHSolar,cooling_data_max;

	zmin = All.CoolingParameters_zmin;
	zmax = All.CoolingParameters_zmax;
	slz = All.CoolingParameters_slz;
	tmin = All.CoolingParameters_tmin;
	tmax = All.CoolingParameters_tmax;
	slt = All.CoolingParameters_slt;
	FeHSolar = All.CoolingParameters_FeHSolar;
	cooling_data_max = All.CoolingParameters_cooling_data_max;


	cooling = 0.0;


	T = log10( temperature );
	Z = log10( metal/FeHSolar + 1.e-10 );


	if (Z>zmax)
	{
	/print ,'Warning: Z>Zmax for',i*/
	Z=zmax;
	}

	if (Z < zmin)
	{

	rt = (T-tmin)/slt;

	it = (int)rt;

	if (it < cooling_data_max )
	it = (int)rt;
	else
	it = cooling_data_max;

	itp = it+1;

	ft = rt - it;

	fz = ( 10. + Z )/( 10. + zmin);


	v1 = ft*( All.CoolingParameters_cooling_data[1][itp] -All.CoolingParameters_cooling_data[1][it] ) + All.CoolingParameters_cooling_data[1][it];
	v2 = ft*( All.CoolingParameters_cooling_data[0][itp] -All.CoolingParameters_cooling_data[0][it] ) + All.CoolingParameters_cooling_data[0][it];
	v = v2 + fz*(v1-v2);

	}
	else
	{

	rt = (T-tmin)/slt;
	rz = (Z-zmin)/slz+1.0;

	it = (int)rt;

	if (it < cooling_data_max )
	it = (int)rt;
	else
	it = cooling_data_max;


	iz = (int)rz;

	itp = it+1;
	izp = iz+1;

	ft = rt - it;
	fz = rz - iz;

	v1 = ft*(All.CoolingParameters_cooling_data[izp][itp] - All.CoolingParameters_cooling_data[izp][it]) + All.CoolingParameters_cooling_data[izp][it];
	v2 = ft*(All.CoolingParameters_cooling_data[iz][itp] - All.CoolingParameters_cooling_data[iz][it]) + All.CoolingParameters_cooling_data[iz][it];
	v = v2 + fz*(v1-v2);


	}

	cooling = pow(10,v);



	return cooling;

	}








	/***************************************************************************

	END OF METALLICITY DEPENDENT COOLING

	**************************************************************************/















	/*! \file cooling.c
	* \brief Compute gas cooling
	*
	*/

	static double hubble_a, a3inv;
	static double eV = 1.6022000e-12;
	static double normfacJ0 = 0.74627;
	static double J0min = 1.e-29;
	static double alpha = 1.0;

	static int Norderweinberg = 7; /* polynom order+1 */
	static double coefweinberg[7][6];
	static double z;
	static double J0;

	static double Cte_G_gHI;
	static double Cte_G_gHeI;
	static double Cte_G_gHeII;
	static double Cte_heating_radiative_HI;
	static double Cte_heating_radiative_HeI;
	static double Cte_heating_radiative_HeII;


	/*
	* init some variables that depends only on redshift
	*/
	void init_from_new_redshift(double Redshift)
	{

	/* init weinberg coeff */

	coefweinberg[0][0] = -0.31086729929951613e+002;
	coefweinberg[1][0] = 0.34803667059463761e+001;
	coefweinberg[2][0] = -0.15145716066316397e+001;
	coefweinberg[3][0] = 0.54649951450632972e+000;
	coefweinberg[4][0] = -0.16395924120387340e+000;
	coefweinberg[5][0] = 0.25197466148524143e-001;
	coefweinberg[6][0] = -0.15352763785487806e-002;

	coefweinberg[0][1] = -0.31887274113252204e+002;
	coefweinberg[1][1] = 0.44178493140927095e+001;
	coefweinberg[2][1] = -0.20158132553082293e+001;
	coefweinberg[3][1] = 0.64080497292269134e+000;
	coefweinberg[4][1] = -0.15981267091909040e+000;
	coefweinberg[5][1] = 0.22056900050237707e-001;
	coefweinberg[6][1] = -0.12837570029562849e-002;

	coefweinberg[0][2] = -0.35693331167978656e+002;
	coefweinberg[1][2] = 0.20207245722165794e+001;
	coefweinberg[2][2] = -0.76856976101363744e-001;
	coefweinberg[3][2] = -0.75691470654320359e-001;
	coefweinberg[4][2] = -0.54502220282734729e-001;
	coefweinberg[5][2] = 0.20633345104660583e-001;
	coefweinberg[6][2] = -0.18410307456285177e-002;

	coefweinberg[0][3] = -0.56967559787460921e+002;
	coefweinberg[1][3] = 0.38601174525546353e+001;
	coefweinberg[2][3] = -0.18318926655684415e+001;
	coefweinberg[3][3] = 0.67360594266440688e+000;
	coefweinberg[4][3] = -0.18983466813215341e+000;
	coefweinberg[5][3] = 0.27768907786915147e-001;
	coefweinberg[6][3] = -0.16330066969315893e-002;

	coefweinberg[0][4] = -0.56977907250821026e+002;
	coefweinberg[1][4] = 0.38686249565302266e+001;
	coefweinberg[2][4] = -0.13330942368518774e+001;
	coefweinberg[3][4] = 0.33988839029092172e+000;
	coefweinberg[4][4] = -0.98997915675929332e-001;
	coefweinberg[5][4] = 0.16781612113050747e-001;
	coefweinberg[6][4] = -0.11514328893746039e-002;

	coefweinberg[0][5] = -0.59825233828609278e+002;
	coefweinberg[1][5] = 0.21898162706563347e+001;
	coefweinberg[2][5] = -0.42982055888598525e+000;
	coefweinberg[3][5] = 0.50312144291614215e-001;
	coefweinberg[4][5] = -0.61550639239553132e-001;
	coefweinberg[5][5] = 0.18017109270959387e-001;
	coefweinberg[6][5] = -0.15438891584271634e-002;



	z = Redshift;
	J0 = J_0();

	/* here, we initialize the ctes that uses J_nu(z) */

	/* Tessier */
	/*
	Cte_G_gHI = G_gHI();
	Cte_G_gHeI = G_gHeI();
	Cte_G_gHeII = G_gHeII();
	Cte_heating_radiative_HI = heating_radiative_HI();
	Cte_heating_radiative_HeI = heating_radiative_HeI();
	Cte_heating_radiative_HeII = heating_radiative_HeII();
	*/

	/* Theuns */
	/*
	Cte_G_gHI = G_gHI_t(J0);
	Cte_G_gHeI = G_gHeI_t(J0);
	Cte_G_gHeII = G_gHeII_t(J0);
	Cte_heating_radiative_HI = heating_radiative_HI_t(J0);
	Cte_heating_radiative_HeI = heating_radiative_HeI_t(J0);
	Cte_heating_radiative_HeII = heating_radiative_HeII_t(J0);
	*/

	/* Weinberg */
	Cte_G_gHI = G_gHI_w();
	Cte_G_gHeI = G_gHeI_w();
	Cte_G_gHeII = G_gHeII_w();
	Cte_heating_radiative_HI = heating_radiative_HI_w();
	Cte_heating_radiative_HeI = heating_radiative_HeI_w();
	Cte_heating_radiative_HeII = heating_radiative_HeII_w();

	}

	/*
	* J0
	*/

	double J_0()
	{
	double Fz;

	if (z > 6)
	Fz = 0;
	else
	{
	if (z > 3)
	Fz = 4/(z+1);
	else
	{
	if (z > 2)
	Fz = 1;
	else
	Fz = pow(((1+z)/3.),3);
	}
	}

	return 1.0e-22*Fz;
	}

	/*
	* UV background intensity
	*/

	double J_nu(double e)
	{
	double e_L;
	e_L = 13.598*eV;
	return (e_L/e)*J_0();
	}

	/*
	* sigma_rad
	*/

	double sigma_rad_HI(double e)
	{
	double xxx,alph,e_i;

	e_i = 13.598 *eV;
	xxx = e/e_i;
	alph = sqrt(xxx-1.0);
	return 6.30e-18/pow(xxx,4)exp(4.0-4.0atan(alph)/alph) /(1.0-exp(-TWOPI/alph));
	}


	double sigma_rad_HeI(double e)
	{
	double xxx,alph,e_i;

	e_i = 24.587 *eV;
	xxx = e/e_i;
	alph = sqrt(xxx-1.0);
	return 7.42e-18*(1.660/pow(xxx,2.050)-0.660/pow(xxx,3.050));
	}

	double sigma_rad_HeII(double e)
	{
	double xxx,alph,e_i;

	e_i = 54.416 *eV;
	xxx = e/e_i;
	alph = sqrt(xxx-1.0);
	return 1.58e-18/pow(xxx,4)exp(4.0-4.0atan(alph)/alph)/(1.0-exp(-TWOPI/alph));
	}


	/*
	* cooling rates
	*/

	/* Bremstrahlung */
	double cooling_bremstrahlung_HI(double T)
	{
	return 1.42e-27sqrt(T)(1.10+0.340*exp(-pow((5.50-log10(T)),2) /3.0));
	}

	double cooling_bremstrahlung_HeI(double T)
	{
	return 1.42e-27sqrt(T)(1.10+0.340*exp(-pow((5.50-log10(T)),2) /3.0));
	}

	double cooling_bremstrahlung_HeII(double T)
	{
	return 5.68e-27sqrt(T)(1.10+0.340*exp(-pow((5.50-log10(T)),2) /3.0));
	}

	/* Ionization */
	double cooling_ionization_HI(double T)
	{
	double T5;
	T5 = T/1e5;
	return 2.54e-21sqrt(T)exp(-157809.1/T)/(1+sqrt(T5));
	}

	double cooling_ionization_HeI(double T)
	{
	double T5;
	T5 = T/1e5;
	return 1.88e-21sqrt(T)exp(-285335.4/T)/(1+sqrt(T5));
	}

	double cooling_ionization_HeII(double T)
	{
	double T5;
	T5 = T/1e5;
	return 9.90e-22sqrt(T)exp(-631515.0/T)/(1+sqrt(T5));
	}

	/* Recombination */
	double cooling_recombination_HI(double T)
	{
	double T3,T6;
	T3 = T/1e3;
	T6 = T/1e6;
	return 8.70e-27*sqrt(T)/pow(T3,0.2)/(1.0+pow(T6,0.7));
	}

	double cooling_recombination_HeI(double T)
	{
	return 1.55e-26*pow(T,0.3647);
	}

	double cooling_recombination_HeII(double T)
	{
	double T3,T6;
	T3 = T/1e3;
	T6 = T/1e6;
	return 3.48e-26*sqrt(T)/pow(T3,0.2)/(1.0+pow(T6,0.7));
	}

	/* Dielectric Recombination */
	double cooling_dielectric_recombination(double T)
	{
	return 1.24e-13pow(T,-1.5)exp(-470000.0/T)(1.0+0.3exp(-94000.0/T));
	}

	/* Ecitation cooling (line cooling) */
	double cooling_excitation_HI(double T)
	{
	double T5;
	T5 = T/1e5;
	return 7.50e-19*exp(-118348.0/T)/(1+sqrt(T5));
	}

	double cooling_excitation_HII(double T)
	{
	double T5;
	T5 = T/1e5;
	return 5.54e-17/pow(T,0.397)*exp(-473638.0/T)/(1+sqrt(T5));
	}

	/* Compton cooling */
	double cooling_compton(double T)
	{
	return 5.406e-36(T-2.7(1+z))*pow((1+z),4);
	}

	/*
	* recombination rates (taux_rec)
	*/

	double A_HII(double T)
	{
	double T3,T6;
	T3 = T/1e3;
	T6 = T/1e6;
	return 6.30e-11/sqrt(T)/pow(T3,0.2)/(1+pow(T6,0.7));
	}

	double A_HeIId(double T)
	{
	return 1.9e-3/pow(T,1.50)exp(-470000.0/T)(1.0+0.30*exp(-94000.0/T));
	}

	double A_HeII(double T)
	{
	return 1.5e-10/pow(T,0.6353) + A_HeIId(T);
	}

	double A_HeIII(double T)
	{
	double T3,T6;
	T3 = T/1e3;
	T6 = T/1e6;
	return 3.36e-10/sqrt(T)/pow(T3,0.2)/(1.0+pow(T6,0.7));
	}

	/*
	* collisional rates (taux_ion)
	*/

	double G_HI(double T)
	{
	double T5;
	T5 = T/1e5;
	return 1.17e-10sqrt(T)exp(-157809.1/T)/(1.0+sqrt(T5));
	}

	double G_HeI(double T)
	{
	double T5;
	T5 = T/1e5;
	return 2.38e-11sqrt(T)exp(-285335.4/T)/(1.0+sqrt(T5));
	}

	double G_HeII(double T)
	{
	double T5;
	T5 = T/1e5;
	return 5.68e-12sqrt(T)exp(-631515.0/T)/(1.0+sqrt(T5));
	}

	/*
	* photoionisation rates (depend only on z)
	*/

	double G_gHI()
	{
	double e_i,integ,e,de,error;

	e_i = 13.598*eV;
	integ = 0.0;
	e = e_i;
	de = e/100.0;
	error = 1.0;
	while (error>1.e-6)
	{
	e = e + de;
	de = e/100.0;
	error = 2TWOPIJ_nu(e)sigma_rad_HI(e)de/e;
	integ = integ + error;
	error = error/fabs(integ);
	}

	return integ/PLANCK;
	}


	double G_gHeI()
	{
	double e_i,integ,e,de,error;

	e_i = 24.587*eV;
	integ = 0.0;
	e = e_i;
	de = e/100.0;
	error = 1.0;
	while (error>1.e-6)
	{
	e = e + de;
	de = e/100.0;
	error = 2TWOPIJ_nu(e)sigma_rad_HeI(e)de/e;
	integ = integ + error;
	error = error/fabs(integ);
	}

	return integ/PLANCK;
	}

	double G_gHeII()
	{
	double e_i,integ,e,de,error;

	e_i = 54.416*eV;
	integ = 0.0;
	e = e_i;
	de = e/100.0;
	error = 1.0;
	while (error>1.e-6)
	{
	e = e + de;
	de = e/100.0;
	error = 2TWOPIJ_nu(e)sigma_rad_HeII(e)de/e;
	integ = integ + error;
	error = error/fabs(integ);
	}

	return integ/PLANCK;
	}


	double G_gHI_t(double J0)
	{
	return 1.26e10*J0/(3.0+alpha);

	}

	double G_gHeI_t(double J0)
	{
	return 1.48e10J0pow(0.5530,alpha) *(1.660/(alpha+2.050)-0.660/(alpha+3.050));
	}

	double G_gHeII_t(double J0)
	{
	return 3.34e9J0pow(0.2490,alpha)/(3.0+alpha);
	}




	double G_gHI_w()
	{
	double taux_rad_weinbergint;
	double hh,tt,zz;
	int i;

	if (z < 8.50)
	{
	hh=0.0;
	zz=dmax(z,1.0e-15);
	for (i=0;i<Norderweinberg;i++)
	hh=hh+coefweinberg[i][0]*pow(zz,i);
	taux_rad_weinbergint=normfacJ0*exp(hh);
	}
	else
	taux_rad_weinbergint=0.0;

	tt=G_gHI_t(J0min);
	if (taux_rad_weinbergint < tt)
	taux_rad_weinbergint=tt;

	return taux_rad_weinbergint;
	}

	double G_gHeI_w()
	{
	double taux_rad_weinbergint;
	double hh,tt,zz;
	int i;

	if (z < 8.50)
	{
	hh=0.0;
	zz=dmax(z,1.0e-15);
	for (i=0;i<Norderweinberg;i++)
	hh=hh+coefweinberg[i][1]*pow(zz,i);
	taux_rad_weinbergint=normfacJ0*exp(hh);
	}
	else
	taux_rad_weinbergint=0.0;

	tt=G_gHeI_t(J0min);
	if (taux_rad_weinbergint < tt)
	taux_rad_weinbergint=tt;

	return taux_rad_weinbergint;
	}


	double G_gHeII_w()
	{
	double taux_rad_weinbergint;
	double hh,tt,zz;
	int i;

	if (z < 8.50)
	{
	hh=0.0;
	zz=dmax(z,1.0e-15);
	for (i=0;i<Norderweinberg;i++)
	hh=hh+coefweinberg[i][2]*pow(zz,i);
	taux_rad_weinbergint=normfacJ0*exp(hh);
	}
	else
	taux_rad_weinbergint=0.0;

	tt=G_gHeII_t(J0min);
	if (taux_rad_weinbergint < tt)
	taux_rad_weinbergint=tt;

	return taux_rad_weinbergint;
	}





	/*
	* heating rates (depend only on z)
	*/

	double heating_radiative_HI() /* use J_nu */
	{
	double e_i,integ,e,de,error;

	e_i = 13.598*eV;
	integ = 0.0;
	e = e_i;
	de = e/100.0;
	error = 1.0;

	while(error>1.e-6)
	{
	e = e + de;
	de = e/100.0;
	error = 2.0TWOPIJ_nu(e)sigma_rad_HI(e)(e/e_i-1.0)*de/e;
	integ = integ + error;
	error=error/fabs(integ);
	}
	return integ/PLANCK*e_i;
	}

	double heating_radiative_HeI() /* use J_nu */
	{
	double e_i,integ,e,de,error;

	e_i = 24.587*eV;
	integ = 0.0;
	e = e_i;
	de = e/100.0;
	error = 1.0;

	while(error>1.e-6)
	{
	e = e + de;
	de = e/100.0;
	error = 2.0TWOPIJ_nu(e)sigma_rad_HeI(e)(e/e_i-1.0)*de/e;
	integ = integ + error;
	error=error/fabs(integ);
	}
	return integ/PLANCK*e_i;
	}


	double heating_radiative_HeII() /* use J_nu */
	{
	double e_i,integ,e,de,error;

	e_i = 54.416*eV;
	integ = 0.0;
	e = e_i;
	de = e/100.0;
	error = 1.0;

	while(error>1.e-6)
	{
	e = e + de;
	de = e/100.0;
	error = 2.0TWOPIJ_nu(e)sigma_rad_HeII(e)(e/e_i-1.0)*de/e;
	integ = integ + error;
	error=error/fabs(integ);
	}
	return integ/PLANCK*e_i;
	}


	double heating_radiative_HI_t(double J0) /* use Theuns */
	{
	return (2.91e-1*J0/(2.0+alpha))/(3.0+alpha);
	}

	double heating_radiative_HeI_t(double J0) /* use Theuns */
	{
	return 5.84e-1J0pow(0.5530,alpha)*(1.660/(alpha+1.050)-2.320/(alpha+2.050)+0.660/(alpha+3.050));
	}

	double heating_radiative_HeII_t(double J0) /* use Theuns */
	{
	return (2.92e-1J0pow(0.2490,alpha)/(2.0+alpha))/(3.0+alpha);
	}




	double heating_radiative_HI_w() /* use weinberg coeff */
	{
	double heat_rad_weinbergint;
	double hh,tt,zz;
	int i;

	if (z < 8.50)
	{
	hh=0.0;
	zz=dmax(z,1.0e-15);
	for (i=0;i<Norderweinberg;i++)
	hh=hh+coefweinberg[i][3]*pow(zz,i);
	heat_rad_weinbergint=normfacJ0*exp(hh);
	}
	else
	heat_rad_weinbergint=0.0;

	tt=heating_radiative_HI_t(J0min);
	if (heat_rad_weinbergint < tt)
	heat_rad_weinbergint=tt;

	return heat_rad_weinbergint;
	}

	double heating_radiative_HeI_w() /* use weinberg coeff */
	{
	double heat_rad_weinbergint;
	double hh,tt,zz;
	int i;

	if (z < 8.50)
	{
	hh=0.0;
	zz=dmax(z,1.0e-15);
	for (i=0;i<Norderweinberg;i++)
	hh=hh+coefweinberg[i][4]*pow(zz,i);
	heat_rad_weinbergint=normfacJ0*exp(hh);
	}
	else
	heat_rad_weinbergint=0.0;

	tt=heating_radiative_HeI_t(J0min);
	if (heat_rad_weinbergint < tt)
	heat_rad_weinbergint=tt;

	return heat_rad_weinbergint;
	}

	double heating_radiative_HeII_w() /* use weinberg coeff */
	{
	double heat_rad_weinbergint;
	double hh,tt,zz;
	int i;

	if (z < 8.50)
	{
	hh=0.0;
	zz=dmax(z,1.0e-15);
	for (i=0;i<Norderweinberg;i++)
	hh=hh+coefweinberg[i][5]*pow(zz,i);
	heat_rad_weinbergint=normfacJ0*exp(hh);
	}
	else
	heat_rad_weinbergint=0.0;

	tt=heating_radiative_HeII_t(J0min);
	if (heat_rad_weinbergint < tt)
	heat_rad_weinbergint=tt;

	return heat_rad_weinbergint;
	}







	double heating_compton()
	{
	/* Abel, Tom; Haehnelt, Martin G.Apj 520 */
	//return 5.406e-362.726pow((1+z),5); /* from Ramses */

	//if (z>6)
	// return 0;
	//else
	// return 1.25e-31*pow((1+z),13/3.);

	return 0;
	}




	void compute_densities(double T,double X,double pn_H, double pn_HI,double pn_HII,double pn_HEI,double pn_HEII,double pn_HEIII,double pn_E,double pmu)
	{

	double Y,yy,x1;
	double t_rad_HI,t_rec_HI,t_ion_HI;
	double t_rad_HEI,t_rec_HEI,t_ion_HEI;
	double t_rad_HEII,t_rec_HEII,t_ion_HEII;
	double t_ion2_HI,t_ion2_HEI,t_ion2_HEII;
	double err_nE;
	double n_T;
	double n_H,n_HI,n_HII,n_HEI,n_HEII,n_HEIII,n_E,mu;

	Y = 1-X;
	yy = Y/(4-4*Y);

	t_rad_HI = Cte_G_gHI;
	t_rec_HI = A_HII(T);
	t_ion_HI = G_HI(T);

	t_rad_HEI = Cte_G_gHeI;
	t_rec_HEI = A_HeII(T);
	t_ion_HEI = G_HeI(T);

	t_rad_HEII = Cte_G_gHeII;
	t_rec_HEII = A_HeIII(T);
	t_ion_HEII = G_HeII(T);

	n_H = *pn_H;

	n_E = n_H;
	err_nE = 1.;


	while(err_nE > 1.e-8)
	{

	/* compute densities (Ramses implementation) */
	t_ion2_HI = t_ion_HI + t_rad_HI /dmax(n_E,1e-15*n_H);
	t_ion2_HEI = t_ion_HEI + t_rad_HEI /dmax(n_E,1e-15*n_H);
	t_ion2_HEII = t_ion_HEII + t_rad_HEII/dmax(n_E,1e-15*n_H);

	n_HI = t_rec_HI/(t_ion2_HI+t_rec_HI)*n_H;
	n_HII = t_ion2_HI/(t_ion2_HI+t_rec_HI)*n_H;

	x1 = (t_rec_HEIIt_rec_HEI+t_ion2_HEIt_rec_HEII+t_ion2_HEII*t_ion2_HEI);

	n_HEIII = yyt_ion2_HEIIt_ion2_HEI/x1*n_H;
	n_HEII = yyt_ion2_HEI t_rec_HEII/x1*n_H;
	n_HEI = yyt_rec_HEII t_rec_HEI /x1*n_H;

	err_nE = fabs((n_E - (n_HII + n_HEII + 2.*n_HEIII))/n_H);
	n_E = 0.5n_E+0.5(n_HII + n_HEII + 2.*n_HEIII);

	}

	n_T = n_HI + n_HII+ n_HEI+ n_HEII+ n_HEIII+ n_E;
	mu = n_H/X/n_T;

	*pn_H = n_H;
	*pn_HI = n_HI;
	*pn_HII = n_HII;
	*pn_HEI = n_HEI;
	*pn_HEII = n_HEII;
	*pn_HEIII = n_HEIII;
	*pn_E = n_E;
	*pmu = mu;


	}

	void print_cooling(double T,
	double c1,double c2,double c3,double c4,double c5,double c6,double c7,double c8,double c9,
	double c10,double c11,double c12,double c13,double h1, double h2, double h3, double h4)
	{

	double ctot,htot,chtot;

	ctot = c1+c2+c3+c4+c5+c6+c7+c8+c9+c10+c11+c12+c13;
	htot = h1+h2+h3+h4;
	chtot= ctot - htot;

	printf("%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",T,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,h1,h2,h3,h4,ctot,htot,chtot);

	}


	void compute_cooling_from_T_and_Nh(double T,double X,double n_H,
	double c1,double c2,double c3,double c4,double c5,double c6,double c7,double c8,double *c9,
	double c10,double c11,double c12,double c13,double h1, double h2, double h3, double h4)
	{

	double n_HI,n_HII,n_HEI,n_HEII,n_HEIII,n_E,mu;
	double nH2;
	//double c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13;
	//double h1,h2,h3,h4;

	compute_densities(T,X,&n_H,&n_HI,&n_HII,&n_HEI,&n_HEII,&n_HEIII,&n_E,&mu);

	nH2 = n_H*n_H;

	/*
	* compute cooling
	*/

	/* Bremstrahlung (cool_bre) */
	c1 = cooling_bremstrahlung_HI(T) n_E*n_HII /nH2;
	c2 = cooling_bremstrahlung_HeI(T) n_E*n_HEII /nH2;
	c3 = cooling_bremstrahlung_HeII(T) n_E*n_HEIII/nH2;


	/* Ionization cooling (cool_ion) */
	c4 = cooling_ionization_HI(T) n_E*n_HI /nH2;
	c5 = cooling_ionization_HeI(T) n_E*n_HEI /nH2;
	c6 = cooling_ionization_HeII(T) n_E*n_HEII /nH2;


	/* Recombination cooling (cool_rec) */
	c7 = cooling_recombination_HI(T) n_E*n_HII /nH2;
	c8 = cooling_recombination_HeI(T) n_E*n_HEII /nH2;
	c9 = cooling_recombination_HeII(T) n_E*n_HEIII/nH2;

	/* Dielectric recombination cooling (cool_die) */
	c10 = cooling_dielectric_recombination(T) n_E*n_HEII /nH2;

	/* Line cooling (cool_exc) */
	c11 = cooling_excitation_HI(T) n_E*n_HI /nH2;
	c12 = cooling_excitation_HII(T) n_E*n_HEII /nH2;

	/* Compton cooling (cool_com) */
	c13 = cooling_compton(T) n_E /nH2; /* !! dep on z */


	/*
	* compute heating
	*/

	/* Radiative heating (h_rad_spec) */
	h1 = Cte_heating_radiative_HI n_HI /nH2;
	h2 = Cte_heating_radiative_HeI n_HEI /nH2;
	h3 = Cte_heating_radiative_HeII n_HEII /nH2;


	/* Compton heating (heat_com) */
	h4 = heating_compton() n_E /nH2; /* !! dep on z */

	}



	double compute_cooling_from_Egyspec_and_Density(double Egyspec,double Density,double *MeanWeight)
	{

	double T,mu,n_H;
	double n_HI,n_HII,n_HEI,n_HEII,n_HEIII,n_E;

	double err_mu,mu_left,mu_right,mu_old;
	int niter;

	double c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13;
	double h1,h2,h3,h4;

	double nH2;

	/* Hydrogen density (cgs) */
	n_H = HYDROGEN_MASSFRAC * Density/ PROTONMASS;


	/* itterate to find the right mu and T */

	err_mu=1.;
	mu_left=0.5;
	mu_right=1.3;
	niter=0;

	while ( (err_mu > 1.e-4) && (niter <= 50) )
	{

	mu_old=0.5*(mu_left+mu_right);

	/* compute temperature */
	T = GAMMA_MINUS1 mu_oldPROTONMASS/BOLTZMANN *Egyspec;

	/* compute all */
	compute_densities(T,HYDROGEN_MASSFRAC,&n_H,&n_HI,&n_HII,&n_HEI,&n_HEII,&n_HEIII,&n_E,&mu);

	err_mu = (mu-mu_old)/mu_old;


	if(err_mu>0.)
	{
	mu_left =0.5*(mu_left+mu_right);
	mu_right=mu_right;
	}
	else
	{
	mu_left =mu_left;
	mu_right=0.5*(mu_left+mu_right);
	}

	err_mu=fabs(err_mu);
	niter=niter+1;

	}

	if (niter > 50)
	printf("ERROR : too many iterations.");

	MeanWeight = 0.5(mu_left+mu_right);

	/* now, compute cooling */

	nH2 = n_H*n_H;

	/*
	* compute cooling
	*/

	/* Bremstrahlung (cool_bre) */
	c1 = cooling_bremstrahlung_HI(T) n_En_HII /nH2;
	c2 = cooling_bremstrahlung_HeI(T) n_En_HEII /nH2;
	c3 = cooling_bremstrahlung_HeII(T) n_En_HEIII/nH2;


	/* Ionization cooling (cool_ion) */
	c4 = cooling_ionization_HI(T) n_En_HI /nH2;
	c5 = cooling_ionization_HeI(T) n_En_HEI /nH2;
	c6 = cooling_ionization_HeII(T) n_En_HEII /nH2;


	/* Recombination cooling (cool_rec) */
	c7 = cooling_recombination_HI(T) n_En_HII /nH2;
	c8 = cooling_recombination_HeI(T) n_En_HEII /nH2;
	c9 = cooling_recombination_HeII(T) n_En_HEIII/nH2;

	/* Dielectric recombination cooling (cool_die) */
	c10 = cooling_dielectric_recombination(T) n_En_HEII /nH2;

	/* Line cooling (cool_exc) */
	c11 = cooling_excitation_HI(T) n_En_HI /nH2;
	c12 = cooling_excitation_HII(T) n_En_HEII /nH2;

	/* Compton cooling (cool_com) */
	c13 = cooling_compton(T) n_E /nH2; / !! dep on z */


	/*
	* compute heating
	*/

	/* Radiative heating (h_rad_spec) */
	h1 = Cte_heating_radiative_HI *n_HI /nH2;
	h2 = Cte_heating_radiative_HeI *n_HEI /nH2;
	h3 = Cte_heating_radiative_HeII *n_HEII /nH2;


	/* Compton heating (heat_com) */
	h4 = heating_compton() n_E /nH2; / !! dep on z */

	/* output info */
	//print_cooling(T,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,h1,h2,h3,h4);

	c1 = dmax(c1,0);
	c2 = dmax(c2,0);
	c3 = dmax(c3,0);
	c4 = dmax(c4,0);
	c5 = dmax(c5,0);
	c6 = dmax(c6,0);
	c7 = dmax(c7,0);
	c8 = dmax(c8,0);
	c9 = dmax(c9,0);
	c10 = dmax(c10,0);
	c11 = dmax(c11,0);
	c12 = dmax(c12,0);
	c13 = dmax(c13,0);
	h1 = dmax(h1,0);
	h2 = dmax(h2,0);
	h3 = dmax(h3,0);
	h4 = dmax(h4,0);

	return (c1+c2+c3+c4+c5+c6+c7+c8+c9+c10+c11+c12+c13) - (h1+h2+h3+h4);

	}



	struct cooling_solver_params
	{
	double Entropy;
	double Density;
	int Phase;
	int i;
	double DtEntropyVisc;
	double dt;
	double hubble_a;
	};

	double cooling_solver_function(double EntropyVar, void *params)
	{
	struct cooling_solver_params p = (struct cooling_solver_params ) params;
	double Entropy = p->Entropy;
	double Density = p->Density;
	int Phase = p->Phase;
	int i = p->i;
	double DtEntropyVisc = p->DtEntropyVisc;
	double dt = p->dt;
	double hubble_a = p->hubble_a;

	double DtEntropyRadSph=0;

	#ifdef MULTIPHASE
	switch (Phase)
	{
	case GAS_SPH:
	DtEntropyRadSph = -GAMMA_MINUS1pow(Density,-GAMMA)lambda(Density,EntropyVar,Phase,i)/hubble_a;
	break;
	case GAS_STICKY:
	case GAS_DARK:
	DtEntropyRadSph = -1/(Density * a3inv) *lambda(Density,EntropyVar,Phase,i)/hubble_a;
	break;
	}
	#else
	DtEntropyRadSph = -GAMMA_MINUS1pow(Density,-GAMMA)lambda(Density,EntropyVar,Phase,i)/hubble_a;
	#endif


	return Entropy + (DtEntropyVisc + DtEntropyRadSph)*dt - EntropyVar;
	};





	/*! This function compute the new Entropy due to isochoric cooling
	* using an implicit iteration scheme
	*
	* !!! here Density is already expressed in comobile coord
	*
	*/
	double DoCooling(FLOAT Density,FLOAT Entropy,int Phase,int i,FLOAT DtEntropyVisc, double dt, double hubble_a)
	{

	double EntropyNew;
	double Entropy_lo=0, Entropy_hi=0;
	double lo,hi;

	int status;
	int iter = 0;
	int max_iter = 100;

	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	gsl_function F;
	struct cooling_solver_params params = {(double)Entropy,(double)Density,(int)Phase,(int)i,(double)DtEntropyVisc,(double)dt,(double)hubble_a};


	F.function = &cooling_solver_function;
	F.params = &params;

	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc (T);


	Entropy_lo = 0.5*Entropy;
	Entropy_hi = 1.1*Entropy;


	lo = cooling_solver_function(Entropy_lo,&params);
	hi = cooling_solver_function(Entropy_hi,&params);


	if (lo*hi>0)
	{

	do
	{
	Entropy_hi = 2* Entropy_hi;
	Entropy_lo = 0.5*Entropy_lo;

	lo = cooling_solver_function(Entropy_lo,&params);
	hi = cooling_solver_function(Entropy_hi,&params);

	//printf("here, we need to iterate...\n");


	}
	while (lo*hi>0);
	}





	gsl_root_fsolver_set (s, &F, Entropy_lo, Entropy_hi);


	do
	{
	iter++;
	status = gsl_root_fsolver_iterate (s);
	EntropyNew = gsl_root_fsolver_root (s);
	Entropy_lo = gsl_root_fsolver_x_lower (s);
	Entropy_hi = gsl_root_fsolver_x_upper (s);
	status = gsl_root_test_interval (Entropy_lo, Entropy_hi,0, 0.001);


	}
	while (status == GSL_CONTINUE && iter < max_iter);

	gsl_root_fsolver_free (s);


	if (status!=GSL_SUCCESS)
	{
	printf("WARNING, HERE WE DO NOT CONVERGE...%g %g\n",Entropy_lo,Entropy_hi);
	endrun(3737);
	}



	return EntropyNew;

	}



	/*! This function computes the entropy variation due to the cooling.
	* Cooling is computed only for sph active particles.
	*/
	void cooling()
	{

	int i;
	double dt=0;
	double EntropyNew;

	/* set the right Redshift and compute value indep of Temperature */
	if (All.CoolingType==1)
	init_from_new_redshift(1.0 / (All.Time) - 1);

	if(All.ComovingIntegrationOn)
	{

	hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
	+ (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;

	hubble_a = All.Hubble * sqrt(hubble_a);

	a3inv = 1 / (All.Time * All.Time * All.Time);
	}
	else
	a3inv = hubble_a = 1;


	for(i = 0; i < NumPart; i++)
	{
	if(P[i].Ti_endstep == All.Ti_Current) /* active particles */
	{


	if(P[i].Type == 0) /* SPH stuff */
	{


	//DtEntropyRadSph=0.;
	//SphP[i].DtEntropyRadSph = 0;


	#ifdef MULTIPHASE
	if (SphP[i].Phase == GAS_SPH)
	{
	#endif
	/* note : SphP[i].DtEntropyRadSph should not be necessary */


	dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval;
	//SphP[i].DtEntropyRadSph = -GAMMA_MINUS1pow(SphP[i].Density a3inv,-GAMMA)lambda(SphP[i].Density a3inv,SphP[i].Entropy,0,i)/hubble_a;

	//if (fabs((SphP[i].DtEntropyRadSph+SphP[i].DtEntropy)dt) > 0.1fabs(SphP[i].Entropy))
	{
	/* do implicit isochoric cooling */
	EntropyNew = DoCooling(SphP[i].Density*a3inv,SphP[i].Entropy,0,i,SphP[i].DtEntropy,dt,hubble_a);

	if(dt > 0)
	SphP[i].DtEntropy = (EntropyNew - SphP[i].Entropy)/dt;
	}
	//else
	{
	// SphP[i].DtEntropy += SphP[i].DtEntropyRadSph;
	}

	//SphP[i].DtEgySpecRadSph = - 1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1) * (SphP[i].DtEntropyRadSph);


	#ifdef MULTIPHASE
	}
	else /* STICKY OR DARK */
	{

	//SphP[i].DtEntropyRadSph = -1/(SphP[i].Density * a3inv)lambda(SphP[i].Density a3inv,SphP[i].Entropy,SphP[i].Phase,i)/hubble_a;
	//SphP[i].DtEntropy += SphP[i].DtEntropyRadSph;

	/* do implicit isochoric cooling */
	dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval;
	EntropyNew = DoCooling(SphP[i].Density*a3inv,SphP[i].Entropy,SphP[i].Phase,i,SphP[i].DtEntropy,dt,hubble_a);

	if(dt > 0)
	SphP[i].DtEntropy = (EntropyNew - SphP[i].Entropy)/dt;

	/* !!! here, we do not take into account the energy variation !!! */
	/* SphP[i].DtEgySpecRadSph = SphP[i].DtEntropy, no ?*/

	}
	#endif

	/* finally sum to the entropy variation */

	/* WARNING : we do not compute DtEntropy here, it is updated in timestep.c */
	/* no, because, it is updated juste above, no ? */
	//SphP[i].DtEntropy += SphP[i].DtEntropyRadSph;
	//SphP[i].DtEntropy += SphP[i].DtEgySpecRadSph / (-1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1));

	}
	}
	}

	}





	/*! This function computes the new entropy due to the cooling,
	* between step t0 and t1.
	*/
	void CoolingForOne(int i,int tstart,int tend,double a3inv,double hubble_a)
	{

	double dt,dadt,tcool,dt_entr;
	double MinSizeTimestep,ErrTolIntAccuracy;
	int ti_current,istep;
	int ti_step;
	double minentropy;

	FLOAT Entropy,DEntropy,DtEntropy,DtEgySpec;



	if(All.MinEgySpec)
	minentropy = All.MinEgySpec * GAMMA_MINUS1 / pow(SphP[i].Density * a3inv, GAMMA_MINUS1);



	/* compute dt */
	/* here we use the convention of Gadget */
	/* this assume that DtEntropy = dA/dt/hubble_a */
	/* and not only dA/dt */
	dt_entr = (tend - tstart) * All.Timebase_interval;


	ErrTolIntAccuracy = 0.02;
	MinSizeTimestep = 0.01*dt_entr;


	/* compute da/dt */
	//dadt = fabs( -GAMMA_MINUS1pow(SphP[i].Density a3inv,-GAMMA)lambda(SphP[i].Density a3inv,SphP[i].Entropy,0,i)/hubble_a );

	/* compute cooling time */
	//tcool = SphP[i].Entropy / dadt;

	//if (ErrTolIntAccuracy*tcool/dt_entr < 1)
	// printf("** %g %g\n",ErrTolIntAccuracytcool,dt_entr); / --> verifier le cooling time */



	/***************************************/
	/* integrate with adaptative timesteps */
	/***************************************/


	Entropy = SphP[i].Entropy;
	ti_current = tstart;
	istep = 0;

	#ifdef CHIMIE_THERMAL_FEEDBACK
	int no_cooling_SNII,no_cooling_SNIa,no_cooling;
	int Tis,Tic;
	double td;

	no_cooling=0;
	no_cooling_SNIa=0;
	no_cooling_SNII=0;

	/* check if we are in an adiabatic phase or not */

	if (All.ComovingIntegrationOn)
	{
	Tic = All.Ti_Current;


	if (SphP[i].SNIaThermalTime>0) /* only if the time has been set at least once, it is negative instead (see init.c) */
	{
	Tis = log(SphP[i].SNIaThermalTime/All.TimeBegin) / All.Timebase_interval;
	td = get_cosmictime_difference(Tis,Tic);

	if(td<All.ChimieSNIaThermalTime)
	no_cooling_SNIa=1;
	}

	if (SphP[i].SNIIThermalTime>0) /* only if the time has been set at least once, it is negative instead (see init.c) */
	{
	Tis = log(SphP[i].SNIIThermalTime/All.TimeBegin) / All.Timebase_interval;
	td = get_cosmictime_difference(Tis,Tic);

	if(td<All.ChimieSNIIThermalTime)
	no_cooling_SNII=1;
	}

	}
	else
	{
	if (SphP[i].SNIaThermalTime>0) /* only if the time has been set at least once, it is negative instead (see init.c) */
	if ((All.Time-SphP[i].SNIaThermalTime)<All.ChimieSNIaThermalTime)
	no_cooling_SNIa=1;

	if (SphP[i].SNIIThermalTime>0) /* only if the time has been set at least once, it is negative instead (see init.c) */
	if ((All.Time-SphP[i].SNIIThermalTime)<All.ChimieSNIIThermalTime)
	no_cooling_SNII=1;
	}


	no_cooling=no_cooling_SNIa+no_cooling_SNII;


	if(no_cooling)
	{
	Entropy = Entropy + SphP[i].DtEntropy* dt_entr;

	/* avoid Entropy to be less than minentropy */
	if(All.MinEgySpec)
	if(Entropy < minentropy)
	Entropy = minentropy;

	/* update particle */
	SphP[i].DtEntropy = (Entropy-SphP[i].Entropy)/dt_entr;
	SphP[i].Entropy = Entropy;

	return ;
	}
	#endif



	while (ti_current<tend)
	{

	/* compute da/dt */
	dadt = fabs( -GAMMA_MINUS1pow(SphP[i].Density a3inv,-GAMMA)lambda(SphP[i].Density a3inv,Entropy,0,i)/hubble_a );

	/* compute cooling time */
	/* this is similar in comobile integraction */
	tcool = Entropy / dadt;


	/* find dt */
	dt = dmax(MinSizeTimestep, tcool*ErrTolIntAccuracy);
	dt = dmin(dt,dt_entr);


	ti_step = dt / All.Timebase_interval;

	ti_step = imax(1,ti_step);
	ti_step = imin(ti_step,tend-ti_current);

	dt = ti_step* All.Timebase_interval;


	#ifndef IMPLICIT_COOLING_INTEGRATION
	/* normal integration of Entropy */
	Entropy += SphP[i].DtEntropy* dt; /* viscosity */
	Entropy += -GAMMA_MINUS1pow(SphP[i].Density a3inv,-GAMMA)lambda(SphP[i].Density a3inv,Entropy,0,i)/hubble_a dt; / cooling */
	#else
	/* or use implicit integration of Entropy */
	/* need this if there is also heating like UV */

	if(All.ComovingIntegrationOn)
	{
	printf("CoolingForOne : this must be checked !\n");
	endrun(123321);
	}
	Entropy = DoCooling(SphP[i].Density*a3inv,Entropy,0,i,SphP[i].DtEntropy,dt,hubble_a);
	#endif



	/* avoid Entropy to be less than minentropy */
	if(All.MinEgySpec)
	if(Entropy < minentropy)
	{
	Entropy = minentropy;
	break;
	}

	ti_current += ti_step;
	istep = istep+1;

	}

	/* entropy only due to cooling */
	DEntropy = Entropy-SphP[i].Entropy - SphP[i].DtEntropy* dt_entr;
	DtEntropy = DEntropy/dt_entr;


	//if (istep>1)
	// printf("* %d %d %d %g %g\n",istep,ti_current,tend,SphP[i].DtEntropy,DtEntropy);



	/* update particle */

	SphP[i].Entropy = Entropy;
	SphP[i].DtEntropy = SphP[i].DtEntropy + DtEntropy;

	DtEgySpec = - 1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1) * (DtEntropy);


	LocalSysState.RadiatedEnergy += DtEgySpec * dt_entr * P[i].Mass;



	#ifdef CHECK_ENTROPY_SIGN
	if (SphP[i].Entropy < 0)
	{
	printf("\ntask=%d: entropy less than zero in CoolingForOne !\n", ThisTask);
	printf("ID=%d Entropy=%g EntropyPred=%g DtEntropy=%g\n",P[i].ID,SphP[i].Entropy,SphP[i].EntropyPred,SphP[i].DtEntropy);
	fflush(stdout);
	endrun(444001);

	}
	#endif




	}













	/*! cooling function
	*
	*/
	double lambda(FLOAT Density,FLOAT Entropy,int phase,int i)
	{
	/*
	* These function returns the Lambda (not the Lambda_n)
	* Here, we assume that Lambda may also contain the heating term.
	*
	* Here, the Entropy and Density are physical, but in h units
	*
	* The function is used only in cooling.c
	*
	*/

	double EgySpec;
	double MeanWeight;
	double T=0,nH=0,nH2=0,l=0;
	double nHcgs=0,nH2cgs=0;

	#ifdef HEATING
	double Gpe=0;
	double X,XTne,eps,ne,flux_in_cgs;
	#endif



	/* number of Hydrogen atoms per unit volume (user units, not corrected from h : [nH] = h2/cm^3 ) */
	#ifndef DO_NO_USE_HYDROGEN_MASSFRAC_IN_COOLING
	nH = HYDROGEN_MASSFRAC*Density/All.ProtonMass;
	#else
	nH = 1 *Density/All.ProtonMass;
	#endif
	nH2 = nH*nH;


	/* in cgs, corrected from h */
	nHcgs = nH/pow(All.UnitLength_in_cm, 3)(All.HubbleParamAll.HubbleParam);
	nH2cgs = nHcgs*nHcgs;


	/* compute temperature */
	#ifdef MULTIPHASE
	switch(phase)
	{
	case GAS_SPH:
	T = All.mumh/All.Boltzmann * Entropy * pow(Density,GAMMA_MINUS1);
	break;
	case GAS_STICKY:
	case GAS_DARK:
	T = All.mumh/All.Boltzmann * GAMMA_MINUS1 * Entropy;
	break;
	}
	#else
	T = All.mumh/All.Boltzmann * Entropy * pow(Density,GAMMA_MINUS1);
	#endif


	/*******************

	* * * COOLING * * *

	*******************/

	if (All.CoolingType==0 \|\| All.CoolingType==2)
	{

	/**************/
	/* Sutherland */
	/**************/

	#ifdef MULTIPHASE
	switch(phase)
	{
	case GAS_SPH:
	if (T > All.CutofCoolingTemperature)
	if (All.CoolingType==0)
	l = cooling_function(T);
	else
	#ifdef CHIMIE
	l = cooling_function_with_metals(T,SphP[i].Metal[FE]);
	#else
	l = cooling_function_with_metals(T,(pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar);
	#endif
	else
	l = 0;
	break;
	case GAS_STICKY:
	case GAS_DARK:
	if (T > All.CutofCoolingTemperature)
	if (All.CoolingType==0)
	l = cooling_function(T);
	else
	#ifdef CHIMIE
	l = cooling_function_with_metals(T,SphP[i].Metal[FE]);
	#else
	l = cooling_function_with_metals(T,(pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar);
	#endif
	else
	l = 0;
	break;
	}

	#else
	/* here, lambda' is in ergcm^3/s = kgm^5/s^3 */
	if (T > All.CutofCoolingTemperature)
	if (All.CoolingType==0)
	l = cooling_function(T);
	else
	#ifdef CHIMIE
	l = cooling_function_with_metals(T,SphP[i].Metal[FE]);
	#else
	l = cooling_function_with_metals(T,(pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar);
	#endif
	else
	l = 0;
	#endif
	}
	else
	{

	/******************************/
	/* cooling with UV background */
	/******************************/

	/* get the right density and egyspec in cgs */
	/* entropy and density are already physical */

	#ifdef MULTIPHASE

	/* WARNING, HERE, WE MUST DIFERENCIATE ACORDING TO THE PHASE... */
	printf("WARNING, HERE, WE MUST DIFERENCIATE ACORDING TO THE PHASE...\n");
	exit(0);

	// if (phase == GAS_SPH)
	// EgySpec = Entropy / GAMMA_MINUS1 * pow(Density, GAMMA_MINUS1);
	// else
	// EgySpec = Entropy;
	#else
	EgySpec = Entropy / GAMMA_MINUS1 * pow(Density, GAMMA_MINUS1);
	#endif

	/* into cgs, corrected from h */
	EgySpec *= All.UnitEnergy_in_cgs/All.UnitMass_in_g;
	Density *= All.UnitDensity_in_cgs;
	//if(All.ComovingIntegrationOn)
	// Density = (All.HubbleParamAll.HubbleParam);

	/* compute cooling from EnergySpec and Density */
	l = compute_cooling_from_Egyspec_and_Density(EgySpec,Density,&MeanWeight);

	/* compute temperature */
	/*
	Temperature = GAMMA_MINUS1 MeanWeightPROTONMASS/BOLTZMANN *EgySpec;
	//printf("%g %g %g\n",Temperature,MeanWeight,Lambda);
	logT = log10(Temperature);
	*/

	}



	/*******************

	* * * HEATING * * *

	*******************/
	#ifdef HEATING


	#ifdef HEATING_PE
	/**************************/
	/* Photo-electric heating / / all must be in cgs */
	/**************************/
	X = 0;

	#ifdef STELLAR_FLUX
	flux_in_cgs = SphP[i].EnergyFlux* All.UnitEnergy_in_cgs/All.UnitTime_in_s/pow(All.UnitLength_in_cm, 2);
	X = X + flux_in_cgs/C / All.HeatingPeSolarEnergyDensity;
	#endif

	#ifdef EXTERNAL_FLUX
	X = X + All.HeatingExternalFLuxEnergyDensity/All.HeatingPeSolarEnergyDensity ;
	#endif

	ne = nHcgs*All.HeatingPeElectronFraction;
	XTne = X*sqrt(T)/ne;
	eps = 4.87e-2/(1+4e-3pow(XTne,0.73)) + 3.65e-2(T/1e4)/(1+2e-4*XTne);

	Gpe = (1e-24 * eps * X * nHcgs)/ nH2cgs ;

	l = l - Gpe;

	#endif /HEATING_PE/








	#endif



	/**********************************

	* * * final unit conversions * * *

	***********************************/

	/* convert lambda' in user units */
	l = l / All.UnitEnergy_in_cgs /pow(All.UnitLength_in_cm,3) * All.UnitTime_in_s;

	/* in unit with h */
	l = l*All.HubbleParam;


	/* correct from h */
	/*
	* [ Lambda / H / rho_p ] = [u] = cm^2/s^2
	*
	* [H] = h/s
	* [rho_p] = g/cm^3 * h^2
	* [Lambda_n] = g * m^5 / s^3
	* [n] = h^2/m^5
	*
	* => Lambda_n must be multiplied by h (in order to remove one h !! not a unit !!)
	*
	*/
	//if(All.ComovingIntegrationOn)
	// l = l * All.HubbleParam;

	/* get the final lambda by multiplying lambda' by nH2 (all in user units) */
	l = l*nH2;

	return l;
	}


	#endif

cooling.cNo OneTemporaryActions

File Metadata

cooling.cView Options

Event Timeline

cooling.c
No OneTemporary
Actions

cooling.c
View Options