Homec4science
Diffusion Gear (master) Browse

Gear/src/chimie.c
b9a4fa236396master

Last Change
Raw File

chimie.c
View Options

This document is not UTF8. It was detected as ISO-8859-1 (Latin 1) and converted to UTF8 for display.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"
#ifdef CHIMIE
#ifdef CHIMIE_FROM_HDF5
#include <hdf5.h>
#include "hdf5io.h"
#endif
#ifdef PYCHEM
#include <Python.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <numpy/arrayobject.h>
/*
****************************************************
these variables are already defined in Gadget (or not needed)
****************************************************
*/
#define TO_DOUBLE(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_DOUBLE) ,0) )
#endif /* PYCHEM */
const char *get_filename_ext(const char *filename) {
const char *dot = strrchr(filename, '.');
if(!dot || dot == filename) return "";
return dot + 1;
}
/****************************************************************************************/
/*
/*
/*
/* COMMON CHIMIE PART
/*
/*
/*
/****************************************************************************************/
#define MAXPTS 10
#define MAXDATASIZE 200
#define MAXDATAZSIZE 110
/********************************************************
hdf5 reading routines
*********************************************************/
#ifdef CHIMIE_FROM_HDF5
#define HEADER_GRP "/Header"
#define DATA_GRP "/Data"
/********************************************************
Chemistry functions
*********************************************************/
struct ChemistryHeaderStruct
{
char * version;
char * author;
char * date;
};
struct ChemistryDataAttributesStruct
{
int nelts;
char **elts;
double *SolarMassAbundances;
double MeanWDMass;
char *SNIIYieldsFile;
char *SNIIHeliumCoreFile;
char *DYINYieldsFile;
char *DYINHeliumCoreFile;
char *SNIaFile;
char *SolarAbundancesFile;
};
struct ChemistryYieldsTableStruct
{
int nbins;
double min;
double step;
char *label;
double *data;
};
struct ChemistryYieldsTable2DStruct
{
int nx;
int ny;
double x0;
double y0;
double dx;
double dy;
char *label;
double **data;
};
struct ChemistryLiveTimesStruct
{
int nx;
int ny;
double **coeff_z;
};
struct ChemistryIMFStruct
{
int n;
double *ms;
double *as;
char *bs;
double Mmin;
double Mmax;
};
struct ChemistrySNIaStruct
{
double Mpl;
double Mpu;
double a;
double Mdl1;
double Mdu1;
double bb1;
double Mdl2;
double Mdu2;
double bb2;
struct ElementsDataStruct Metals;
};
struct ChemistrySNIIStruct
{
double Mmin;
double Mmax;
int npts;
int nelts;
char **elts;
struct ChemistryYieldsTableStruct *table;
};
struct ChemistryDYINStruct
{
double Mmin;
double Mmax;
int npts;
int nelts;
char **elts;
struct ChemistryYieldsTableStruct *table;
struct ChemistryYieldsTable2DStruct *tableZ;
int Zflag;
int nptZs;
};
/*! Read a table containing yields
* stored in the dataset named "name"
*/
struct ChemistryYieldsTableStruct readYieldsTable(hid_t group,char *name)
{
struct ChemistryYieldsTableStruct table;
hid_t dset;
herr_t status;
table.data = readDatasetAsArrayDouble(group,name);
/* read attributes */
dset = H5Dopen( group , name, H5P_DEFAULT);
table.nbins = readAttributeAsInt(dset,"nbins");
table.min = readAttributeAsDouble(dset,"min");
table.step = readAttributeAsDouble(dset,"step");
table.label = readAttributeAsString(dset,"label");
//printYieldsTable(table);
status = H5Dclose(dset);
return table;
}
/*! Read a table containing 2d yields
* stored in the dataset named "name"
*/
struct ChemistryYieldsTable2DStruct readYieldsTable2D(hid_t group,char *name)
{
struct ChemistryYieldsTable2DStruct table;
hid_t dset;
herr_t status;
table.data = readDatasetAsArray2DDouble_v0(group,name);
/* read attributes */
dset = H5Dopen( group , name, H5P_DEFAULT);
table.nx = readAttributeAsInt(dset,"nx");
table.ny = readAttributeAsInt(dset,"ny");
table.x0 = readAttributeAsDouble(dset,"x0");
table.y0 = readAttributeAsDouble(dset,"y0");
table.dx = readAttributeAsDouble(dset,"dx");
table.dy = readAttributeAsDouble(dset,"dy");
table.label = readAttributeAsString(dset,"label");
//printYieldsTable(table);
status = H5Dclose(dset);
return table;
}
/*! Read the attributes linked to the data group
*/
int readDataAttributes(hid_t table,struct ChemistryDataAttributesStruct *Param)
{
hid_t group;
herr_t status;
group = H5Gopen(table, DATA_GRP,H5P_DEFAULT);
Param->nelts = readAttributeAsInt(group,"nelts");
Param->elts = readAttributeAsArrayString(group,"elts");
Param->SolarMassAbundances = readAttributeAsArrayDouble(group,"SolarMassAbundances");
Param->MeanWDMass = readAttributeAsDouble(group,"MeanWDMass");
Param->elts = readAttributeAsArrayString(group,"elts");
Param->SNIIYieldsFile = readAttributeAsString(group,"SNIIYieldsFile");
Param->SNIIHeliumCoreFile = readAttributeAsString(group,"SNIIHeliumCoreFile");
Param->DYINYieldsFile = readAttributeAsString(group,"DYINYieldsFile");
Param->DYINHeliumCoreFile = readAttributeAsString(group,"DYINHeliumCoreFile");
Param->SNIaFile = readAttributeAsString(group,"SNIaFile");
Param->SolarAbundancesFile = readAttributeAsString(group,"SolarAbundancesFile");
status = H5Gclose (group);
return 0;
}
/*! Print the attributes linked to the data group
*/
int printDataAttributes(struct ChemistryDataAttributesStruct Parameters)
{
int i;
printf("\n");
printf("Data attribute content:\n\n");
printf("\t nelts = %d\n",Parameters.nelts);
printf("\t elts = ");
for (i=0; i<Parameters.nelts; i++)
printf ("%s ",Parameters.elts[i]);
printf ("\n");
printf("\t SolarMassAbundances = ");
for (i=0; i<Parameters.nelts; i++)
printf ("%g ",Parameters.SolarMassAbundances[i]);
printf ("\n");
printf("\t MeanWDMass = %g\n",Parameters.MeanWDMass);
printf("\n");
printf("\t SNIIYieldsFile = %s\n",Parameters.SNIIYieldsFile);
printf("\t SNIIHeliumCoreFile = %s\n",Parameters.SNIIHeliumCoreFile);
printf("\n");
printf("\t DYINYieldsFile = %s\n",Parameters.DYINYieldsFile);
printf("\t DYINHeliumCoreFile = %s\n",Parameters.DYINHeliumCoreFile);
printf("\n");
printf("\t SNIaFile = %s\n",Parameters.SNIaFile);
printf("\t SolarAbundancesFile = %s\n",Parameters.SolarAbundancesFile);
printf("\n");
return 0;
}
/*! Read the attributes linked to the header group
*/
int readHeader(hid_t table,struct ChemistryHeaderStruct *Header)
{
hid_t group;
herr_t status;
group = H5Gopen(table, HEADER_GRP,H5P_DEFAULT);
/* read attribute */
Header->version = readAttributeAsString(group,"version");
Header->author = readAttributeAsString(group,"author");
Header->date = readAttributeAsString(group,"date");
status = H5Gclose (group);
return 0;
}
/*! Print the attributes linked to the header group
*/
int printHeader(struct ChemistryHeaderStruct Header)
{
printf("\n");
printf("Header content:\n\n");
printf("\t version = %s\n",Header.version);
printf("\t author = %s\n",Header.author);
printf("\t date = %s\n",Header.date);
printf("\n");
return 0;
}
/*! Print the content of a yield table
*/
int printYieldsTable(struct ChemistryYieldsTableStruct table,char * space)
{
int i;
printf("%s%s:\n\n",space,table.label);
printf("%s\t nbins = %d\n",space,table.nbins);
printf("%s\t min = %g\n",space,table.min);
printf("%s\t step = %g\n",space,table.step);
printf("%s\t label = %s\n",space,table.label);
printf("\n");
for(i=0;i<3;i++)
printf("%s\t [%d] %g\n",space,i,table.data[i]);
printf("%s\t ...\n",space);
for(i=table.nbins-3;i<table.nbins;i++)
printf("%s\t [%d] %g\n",space,i,table.data[i]);
printf("\n");
return 0;
}
/*! Print the content of a yield table
*/
int printYieldsTable2D(struct ChemistryYieldsTable2DStruct table,char * space)
{
int i,j;
printf("%s%s:\n\n",space,table.label);
printf("%s\t nx = %d\n",space,table.nx);
printf("%s\t x0 = %g\n",space,table.x0);
printf("%s\t dx = %g\n",space,table.dx);
printf("%s\t ny = %d\n",space,table.ny);
printf("%s\t y0 = %g\n",space,table.y0);
printf("%s\t dy = %g\n",space,table.dy);
printf("%s\t label = %s\n",space,table.label);
printf("\n");
for(i=0;i<3;i++)
{
printf("%s\t [%d] ",space,i);
for (j=0;j<table.ny;j++)
printf("%g ",table.data[i][j]);
printf("\n");
}
printf("%s\t ...\n",space);
for(i=table.nx-3;i<table.nx;i++)
{
printf("%s\t [%d] ",space,i);
for (j=0;j<table.ny;j++)
printf("%g ",table.data[i][j]);
printf("\n");
}
printf("\n");
return 0;
}
int readLiveTimes(hid_t table,struct ChemistryLiveTimesStruct *livetimes)
{
hid_t group;
hid_t subgroup;
hid_t dset;
herr_t status;
group = H5Gopen(table, DATA_GRP,H5P_DEFAULT);
subgroup = H5Gopen(group,"LiveTimes",H5P_DEFAULT);
livetimes->coeff_z = readDatasetAsArray2DDouble_v0(subgroup,"coeff_z");
dset = H5Dopen( subgroup , "coeff_z", H5P_DEFAULT);
livetimes->nx = readAttributeAsInt(dset,"nx");
livetimes->ny = readAttributeAsInt(dset,"ny");
status = H5Dclose (dset);
status = H5Gclose (subgroup);
status = H5Gclose (group);
return 0;
}
/*! Print the content of an LiveTimes struct
*/
int printLiveTimes(struct ChemistryLiveTimesStruct livetimes)
{
int i,j;
printf("Output from LiveTimes:\n\n");
printf("\t coeff_z = \n");
for (i=0; i<livetimes.nx; i++) {
printf ("\t\t [");
for (j=0; j<livetimes.ny; j++)
printf (" %6.4f", livetimes.coeff_z[i][j]);
printf ("]\n");
}
printf("\n");
}
int readIMF(hid_t table,struct ChemistryIMFStruct *imf)
{
hid_t group;
hid_t subgroup;
herr_t status;
group = H5Gopen(table, DATA_GRP,H5P_DEFAULT);
subgroup = H5Gopen(group,"IMF",H5P_DEFAULT);
imf->Mmin = readAttributeAsDouble(subgroup,"Mmin");
imf->Mmax = readAttributeAsDouble(subgroup,"Mmax");
imf->n = readAttributeAsInt(subgroup,"n");
imf->ms = readAttributeAsArrayDouble(subgroup,"ms");
imf->as = readAttributeAsArrayDouble(subgroup,"as");
status = H5Gclose (subgroup);
status = H5Gclose (group);
return 0;
}
/*! Print the content of an IMF
*/
int printIMF(struct ChemistryIMFStruct imf)
{
int i;
printf("Output from IMF:\n\n");
printf("\t n = %d\n",imf.n);
printf("\t Mmin = %g\n",imf.Mmin);
printf("\t Mmax = %g\n",imf.Mmax);
printf("\t as = ");
for (i=0; i<imf.n+1; i++)
printf ("%g ",imf.as[i]);
printf ("\n");
//printf("\t bs = ");
//for (i=0; i<imf.n+1; i++)
// printf ("%g ",imf.bs[i]);
//printf ("\n");
printf("\t ms = ");
for (i=0; i<imf.n; i++)
printf ("%g ",imf.ms[i]);
printf ("\n");
printf("\n");
}
int readSNIa(hid_t table,struct ChemistrySNIaStruct *snia)
{
hid_t group;
hid_t subgroup;
herr_t status;
group = H5Gopen(table, DATA_GRP,H5P_DEFAULT);
subgroup = H5Gopen(group,"SNIa",H5P_DEFAULT);
snia->Mpl = readAttributeAsDouble(subgroup,"Mpl");
snia->Mpu = readAttributeAsDouble(subgroup,"Mpu");
snia->a = readAttributeAsDouble(subgroup,"a");
snia->Mdl1 = readAttributeAsDouble(subgroup,"Mdl1");
snia->Mdu1 = readAttributeAsDouble(subgroup,"Mdu1");
snia->bb1 = readAttributeAsDouble(subgroup,"bb1");
snia->Mdl2 = readAttributeAsDouble(subgroup,"Mdl2");
snia->Mdu2 = readAttributeAsDouble(subgroup,"Mdu2");
snia->bb2 = readAttributeAsDouble(subgroup,"bb2");
snia->Metals = readGroupAsElementsData(subgroup,"Metals");
status = H5Gclose (subgroup);
status = H5Gclose (group);
return 0;
}
/*! Print the content of an SNIa struct
*/
int printSNIa(struct ChemistrySNIaStruct snia)
{
int i;
printf("Output from SNIa:\n\n");
printf("\t Mpl = %g\n",snia.Mpl );
printf("\t Mpu = %g\n",snia.Mpu );
printf("\t a = %g\n",snia.a );
printf("\t Mdl1 = %g\n",snia.Mdl1);
printf("\t Mdu1 = %g\n",snia.Mdu1);
printf("\t bb1 = %g\n",snia.bb1 );
printf("\t Mdl2 = %g\n",snia.Mdl2);
printf("\t Mdu2 = %g\n",snia.Mdu2);
printf("\t bb2 = %g\n",snia.bb2 );
printf("\n");
printf("\t Metals:\n\n");
for (i=0; i<snia.Metals.nelts; i++)
printf ("\t\t %s\t= %g\n",snia.Metals.elts[i],snia.Metals.data[i]);
printf("\n");
}
int readSNII(hid_t table,struct ChemistrySNIIStruct *snii)
{
hid_t group;
hid_t subgroup;
herr_t status;
int i;
group = H5Gopen(table, DATA_GRP,H5P_DEFAULT);
subgroup = H5Gopen(group,"SNII",H5P_DEFAULT);
snii->Mmin = readAttributeAsDouble(subgroup,"Mmin");
snii->Mmax = readAttributeAsDouble(subgroup,"Mmax");
snii->npts = readAttributeAsInt(subgroup,"npts");
snii->nelts = readAttributeAsInt(subgroup,"nelts");
snii->elts = readAttributeAsArrayString(subgroup,"elts");
snii->table = malloc( snii->nelts * sizeof(struct ChemistryYieldsTableStruct) );
/* read the yields */
for(i=0;i<snii->nelts;i++)
snii->table[i]=readYieldsTable(subgroup,snii->elts[i]);
status = H5Gclose (subgroup);
status = H5Gclose (group);
return 0;
}
/*! Print the content of an SNII struct
*/
int printSNII(struct ChemistrySNIIStruct snii)
{
int i;
printf("Output from SNII:\n\n");
printf("\t Mmin = %g\n",snii.Mmin);
printf("\t Mmax = %g\n",snii.Mmax);
printf("\t npts = %d\n",snii.npts);
printf("\t nelts = %d\n",snii.nelts);
printf("\t elts = ");
for (i=0; i<snii.nelts; i++)
printf ("%s ",snii.elts[i]);
printf ("\n");
printf("\n");
printf("\t Metals:\n\n");
for(i=0;i<snii.nelts;i++)
printYieldsTable(snii.table[i],"\t\t");
}
int readDYIN(hid_t table,struct ChemistryDYINStruct *dyin)
{
hid_t group;
hid_t subgroup;
hid_t subsubgroup;
herr_t status;
int i;
group = H5Gopen(table, DATA_GRP,H5P_DEFAULT);
subgroup = H5Gopen(group,"DYIN",H5P_DEFAULT);
dyin->Mmin = readAttributeAsDouble(subgroup,"Mmin");
dyin->Mmax = readAttributeAsDouble(subgroup,"Mmax");
dyin->npts = readAttributeAsInt(subgroup,"npts");
dyin->nelts = readAttributeAsInt(subgroup,"nelts");
dyin->elts = readAttributeAsArrayString(subgroup,"elts");
dyin->table = malloc( dyin->nelts * sizeof(struct ChemistryYieldsTableStruct) );
/* read the yields */
for(i=0;i<dyin->nelts;i++)
dyin->table[i]=readYieldsTable(subgroup,dyin->elts[i]);
dyin->Zflag = readAttributeAsInt(subgroup,"Zflag");
if (dyin->Zflag==1)
{
/* read data with Z dependences */
subsubgroup = H5Gopen(subgroup,"MetallicityDependent",H5P_DEFAULT);
dyin->nptZs = readAttributeAsInt(subsubgroup,"nptZs");
dyin->tableZ = malloc( dyin->nelts * sizeof(struct ChemistryYieldsTable2DStruct) );
/* read the yields */
for(i=0;i<dyin->nelts;i++)
dyin->tableZ[i]=readYieldsTable2D(subsubgroup,dyin->elts[i]);
status = H5Gclose (subsubgroup);
/* end read data with Z dependences */
}
status = H5Gclose (subgroup);
status = H5Gclose (group);
return 0;
}
/*! Print the content of an DYIN struct
*/
int printDYIN(struct ChemistryDYINStruct dyin)
{
int i;
printf("Output from DYIN:\n\n");
printf("\t Mmin = %g\n",dyin.Mmin);
printf("\t Mmax = %g\n",dyin.Mmax);
printf("\t npts = %d\n",dyin.npts);
printf("\t nelts = %d\n",dyin.nelts);
if (dyin.Zflag==1)
printf("\t nptZs = %d\n",dyin.nptZs);
printf("\t elts = ");
for (i=0; i<dyin.nelts; i++)
printf ("%s ",dyin.elts[i]);
printf ("\n");
printf("\n");
printf("\t Metals:\n\n");
for(i=0;i<dyin.nelts;i++)
printYieldsTable(dyin.table[i],"\t\t");
if (dyin.Zflag==1)
{
printf("\t Metals Z:\n\n");
for(i=0;i<dyin.nelts;i++)
printYieldsTable2D(dyin.tableZ[i],"\t\t");
}
}
#endif // CHIMIE_FROM_HDF5
/********************************************************
endof hdf5 reading routines
*********************************************************/
static int verbose=0;
static double *MassFracSNII;
static double *MassFracSNIa;
static double *MassFracDYIN;
static double *SingleMassFracSNII;
static double *SingleMassFracSNIa;
static double *SingleMassFracDYIN;
static double *EjectedMass;
static double *SingleEjectedMass;
static double **MassFracSNIIs;
static double **MassFracSNIas;
static double **MassFracDYINs;
static double **SingleMassFracSNIIs;
static double **SingleMassFracSNIas;
static double **SingleMassFracDYINs;
static double **EjectedMasss;
static double **SingleEjectedMasss;
/* intern global variables */
static struct local_params_chimie
{
float coeff_z[3][3];
float Mmin,Mmax;
int n;
float ms[MAXPTS];
float as[MAXPTS+1];
float bs[MAXPTS+1];
float fs[MAXPTS];
double imf_Ntot;
#ifdef CHIMIE_OPTIMAL_SAMPLING
float k;
float m_max;
#endif
float SNII_Mmin;
float SNII_Mmax;
//float SNII_cte;
//float SNII_a;
float SNIa_Mpl;
float SNIa_Mpu;
float SNIa_a;
float SNIa_cte;
float SNIa_Mdl1;
float SNIa_Mdu1;
float SNIa_a1;
float SNIa_b1;
float SNIa_cte1;
float SNIa_bb1;
float SNIa_Mdl2;
float SNIa_Mdu2;
float SNIa_a2;
float SNIa_b2;
float SNIa_cte2;
float SNIa_bb2;
float SNIa_NWDN;
float DYIN_Mmin;
float DYIN_Mmax;
//float DYIN_cte;
//float DYIN_a;
float Mco;
int nptsSNII; /* number of division in mass SNII */
int nptsDYIN; /* number of division in mass DYIN */
int nptZsSNII; /* number of division in Z SNII */
int nptZsDYIN; /* number of division in Z DYIN */
int nelts;
int DYIN_Zflag;
float MassMinForEjecta;
}
*Cps,*Cp;
static struct local_elts_chimie
{
/* SNII */
float MminSNII; /* minimal mass */
float StepSNII; /* log of mass step */
float ArraySNII[MAXDATASIZE]; /* data */
float MetalSNII[MAXDATASIZE]; /* data */
/* DYIN */
float MminDYIN; /* minimal mass */
float StepDYIN; /* log of mass step */
float ArrayDYIN[MAXDATASIZE]; /* data */
float MetalDYIN[MAXDATASIZE]; /* data */
float MminZDYIN; /* minimal Z */
float StepZDYIN; /* step in Z */
float ArrayZDYIN[MAXDATASIZE][MAXDATAZSIZE]; /* data Z */
float MetalZDYIN[MAXDATASIZE][MAXDATAZSIZE]; /* data Z */
float MSNIa;
float SolarMassAbundance;
char label[72];
}
**Elts,*Elt;
/*! This function allocate all varaiables related to the chemistry
*/
void allocate_chimie()
{
int j;
/* allocate Cp */
Cps = malloc((All.ChimieNumberOfParameterFiles) * sizeof(struct local_params_chimie));
/* allocate elts */
Elts = malloc((All.ChimieNumberOfParameterFiles) * sizeof(struct local_elts_chimie));
//for (j=0;j<All.ChimieNumberOfParameterFiles;j++)
// Elt[j] = malloc((nelts) * sizeof(struct local_elts_chimie));
MassFracSNIIs = malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
MassFracSNIas = malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
MassFracDYINs = malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
EjectedMasss = malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
SingleMassFracSNIIs= malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
SingleMassFracSNIas= malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
SingleMassFracDYINs= malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
SingleEjectedMasss = malloc((All.ChimieNumberOfParameterFiles) * sizeof(double));
}
/*! Set the chemistry table to use
*/
void set_table(int i)
{
if (i>=All.ChimieNumberOfParameterFiles)
{
printf("\n set_table : i>= %d !!!\n\n",All.ChimieNumberOfParameterFiles);
endrun(88809);
}
else
{
Cp = &Cps[i];
Elt = Elts[i];
MassFracSNII = MassFracSNIIs[i]; /* all this is useless, no ?*/
MassFracSNIa = MassFracSNIas[i];
MassFracDYIN = MassFracDYINs[i];
SingleMassFracSNII = SingleMassFracSNIIs[i];
SingleMassFracSNIa = SingleMassFracSNIas[i];
SingleMassFracDYIN = SingleMassFracDYINs[i];
EjectedMass = EjectedMasss[i];
SingleEjectedMass = SingleEjectedMasss[i];
}
}
/*! Read the chemistry table (hdf5 format)
*/
#ifdef CHIMIE_FROM_HDF5
void read_chimie_h5(char * filename,int it)
{
hid_t table;
herr_t status;
struct ChemistryHeaderStruct ChemistryHeader;
struct ChemistryDataAttributesStruct ChemistryBasicParameters;
struct ChemistryLiveTimesStruct ChemistryLiveTimes;
struct ChemistryIMFStruct ChemistryIMF;
struct ChemistrySNIaStruct ChemistrySNIa;
struct ChemistrySNIIStruct ChemistrySNII;
struct ChemistryDYINStruct ChemistryDYIN;
table = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
readHeader(table,&ChemistryHeader);
if (verbose && ThisTask==0)
printHeader(ChemistryHeader);
readDataAttributes(table,&ChemistryBasicParameters);
if (verbose && ThisTask==0)
printDataAttributes(ChemistryBasicParameters);
readLiveTimes(table,&ChemistryLiveTimes);
if (verbose && ThisTask==0)
printLiveTimes(ChemistryLiveTimes);
readIMF(table,&ChemistryIMF);
if (verbose && ThisTask==0)
printIMF(ChemistryIMF);
readSNIa(table,&ChemistrySNIa);
if (verbose && ThisTask==0)
printSNIa(ChemistrySNIa);
readSNII(table,&ChemistrySNII);
if (verbose && ThisTask==0)
printSNII(ChemistrySNII);
readDYIN(table,&ChemistryDYIN);
if (verbose && ThisTask==0)
printDYIN(ChemistryDYIN);
status = H5Fclose(table);
Cps[it].MassMinForEjecta = MAX_REAL_NUMBER;
/* convert to chimie struct */
int i,j,k;
/* Livetimes */
for(i=0;i<3;i++)
for(j=0;j<3;j++)
Cps[it].coeff_z[i][j] = ChemistryLiveTimes.coeff_z[i][j];
/* IMF Parameters */
Cps[it].Mmin = ChemistryIMF.Mmin;
Cps[it].Mmax = ChemistryIMF.Mmax;
Cps[it].n = ChemistryIMF.n;
if (Cps[it].n>0)
for (i=0;i<Cps[it].n;i++)
Cps[it].ms[i] = ChemistryIMF.ms[i];
for (i=0;i<Cps[it].n+1;i++)
Cps[it].as[i]= ChemistryIMF.as[i];
/* Parameters for SNII Rates */
Cps[it].SNII_Mmin = ChemistrySNII.Mmin;
Cps[it].SNII_Mmax = ChemistrySNII.Mmax;
Cps[it].MassMinForEjecta = dmin(Cps[it].MassMinForEjecta,Cps[it].SNII_Mmin);
/* Parameters for DYIN Rates */
Cps[it].DYIN_Mmin = ChemistryDYIN.Mmin;
Cps[it].DYIN_Mmax = ChemistryDYIN.Mmax;
Cps[it].DYIN_Zflag= ChemistryDYIN.Zflag;
Cps[it].MassMinForEjecta = dmin(Cps[it].MassMinForEjecta,Cps[it].DYIN_Mmin);
/* Parameters for SNIa Rates */
Cps[it].SNIa_Mpl = ChemistrySNIa.Mpl;
Cps[it].SNIa_Mpu = ChemistrySNIa.Mpu;
Cps[it].SNIa_a = ChemistrySNIa.a;
Cps[it].SNIa_Mdl1 = ChemistrySNIa.Mdl1;
Cps[it].SNIa_Mdu1 = ChemistrySNIa.Mdu1;
Cps[it].SNIa_bb1 = ChemistrySNIa.bb1;
Cps[it].SNIa_Mdl2 = ChemistrySNIa.Mdl2;
Cps[it].SNIa_Mdu2 = ChemistrySNIa.Mdu2;
Cps[it].SNIa_bb2 = ChemistrySNIa.bb2;
Cps[it].MassMinForEjecta = dmin(Cps[it].MassMinForEjecta,Cps[it].SNIa_Mdl1);
Cps[it].MassMinForEjecta = dmin(Cps[it].MassMinForEjecta,Cps[it].SNIa_Mdl2);
/* Metal injection SNII */
Cps[it].nptsSNII = ChemistrySNII.npts;
//if (Cps[it].SNII_Zflag==1)
// Cps[it].nptZsSNII = ChemistrySNII.nptZs;
/* Metal injection DYIN */
Cps[it].nptsDYIN = ChemistryDYIN.npts;
if (Cps[it].DYIN_Zflag==1)
Cps[it].nptZsDYIN = ChemistryDYIN.nptZs;
Cps[it].nelts = ChemistryBasicParameters.nelts;
/* allocate memory for elts */
if ((Cps[it].nptsSNII<=MAXDATASIZE)&&(Cps[it].nptsDYIN<=MAXDATASIZE))
{
Elts[it] = malloc((Cps[it].nelts+2) * sizeof(struct local_elts_chimie));
}
else
{
printf("\n Cps[it].nptsSNII = %d > MAXDATASIZE = %d !!!\n\n",Cps[it].nptsSNII,MAXDATASIZE);
printf("\n Cps[it].nptsDYIN = %d > MAXDATASIZE = %d !!!\n\n",Cps[it].nptsDYIN,MAXDATASIZE);
endrun(88800);
}
/* allocate memory */
MassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); /* really needed ? */
MassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double));
MassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
EjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleEjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double));
/**/
/* injected metals SNII */
/**/
for (i=0;i<Cps[it].nelts+2;i++)
{
//printf("%s\n",ChemistrySNII.table[i].label);
strcpy(Elts[it][i].label,ChemistrySNII.table[i].label);
Elts[it][i].MminSNII = ChemistrySNII.table[i].min;
Elts[it][i].StepSNII = ChemistrySNII.table[i].step;
for (j=0;j<Cps[it].nptsSNII;j++)
Elts[it][i].MetalSNII[j] = ChemistrySNII.table[i].data[j];
}
/* integral of injected metals SNII */
for (i=0;i<Cps[it].nelts+2;i++)
{
//printf("%s\n",ChemistrySNII.table[i+Cps[it].nelts+2].label);
for (j=0;j<Cps[it].nptsSNII;j++)
Elts[it][i].ArraySNII[j] = ChemistrySNII.table[i+Cps[it].nelts+2].data[j];
}
/**/
/* injected metals DYIN */
/**/
for (i=0;i<Cps[it].nelts+2;i++)
{
//printf("%s\n",ChemistryDYIN.table[i].label);
strcpy(Elts[it][i].label,ChemistryDYIN.table[i].label);
Elts[it][i].MminDYIN = ChemistryDYIN.table[i].min;
Elts[it][i].StepDYIN = ChemistryDYIN.table[i].step;
for (j=0;j<Cps[it].nptsDYIN;j++)
Elts[it][i].MetalDYIN[j] = ChemistryDYIN.table[i].data[j];
}
/* integral of injected metals DYIN */
for (i=0;i<Cps[it].nelts+2;i++)
{
//printf("%s\n",ChemistryDYIN.table[i+Cps[it].nelts+2].label);
for (j=0;j<Cps[it].nptsDYIN;j++)
Elts[it][i].ArrayDYIN[j] = ChemistryDYIN.table[i+Cps[it].nelts+2].data[j];
}
/* metallicity dependent table */
if (Cps[it].DYIN_Zflag==1)
{
for (i=0;i<Cps[it].nelts+2;i++)
{
Elts[it][i].MminZDYIN = ChemistryDYIN.tableZ[i].y0;
Elts[it][i].StepZDYIN = ChemistryDYIN.tableZ[i].dy;
for (j=0;j<Cps[it].nptsDYIN;j++)
for (k=0;k<Cps[it].nptZsDYIN;k++)
{
Elts[it][i].MetalZDYIN[j][k] = ChemistryDYIN.tableZ[i].data[j][k];
Elts[it][i].ArrayZDYIN[j][k] = ChemistryDYIN.tableZ[i+Cps[it].nelts+2].data[j][k];
}
}
}
/* Metal injection SNIa */
Cps[it].Mco = ChemistryBasicParameters.MeanWDMass;
/* for each elt in */
for (i=0;i<Cps[it].nelts+2;i++)
{
if (strcmp(Elts[it][i].label,ChemistrySNIa.Metals.elts[i])!=0)
{
printf("%s != %s\n",Elts[it][i].label,ChemistrySNIa.Metals.elts[i]);
endrun(888012);
}
Elts[it][i].MSNIa = ChemistrySNIa.Metals.data[i];
}
/* Solar Mass Abundances */
for (i=0;i<Cps[it].nelts;i++)
{
if (strcmp(Elts[it][i+2].label,ChemistryBasicParameters.elts[i])!=0)
{
printf("%s != %s\n",Elts[it][i+2].label,ChemistryBasicParameters.elts[i]);
endrun(888013);
}
Elts[it][i+2].SolarMassAbundance = ChemistryBasicParameters.SolarMassAbundances[i];
}
//if (verbose && ThisTask==0)
// info(it);
}
#endif // CHIMIE_FROM_HDF5
/*! Read the chemistry table
*/
void read_chimie(char * filename,int it)
{
char line[72],buffer[72];
FILE *fd;
int i,j;
if (verbose && ThisTask==0)
printf("reading %s ...\n",filename);
fd = fopen(filename,"r");
/* read Lifetime */
/* #### Livetime #### */
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%g %g %g\n", &Cps[it].coeff_z[0][0],&Cps[it].coeff_z[0][1],&Cps[it].coeff_z[0][2]);
fscanf(fd, "%g %g %g\n", &Cps[it].coeff_z[1][0],&Cps[it].coeff_z[1][1],&Cps[it].coeff_z[1][2]);
fscanf(fd, "%g %g %g\n", &Cps[it].coeff_z[2][0],&Cps[it].coeff_z[2][1],&Cps[it].coeff_z[2][2]);
fgets(line, sizeof(line), fd);
/* IMF Parameters */
/* #### IMF Parameters #### */
fgets(line, sizeof(line), fd);
fscanf(fd, "%g %g\n",&Cps[it].Mmin,&Cps[it].Mmax);
fscanf(fd, "%d\n",&Cps[it].n);
if (Cps[it].n>0)
for (i=0;i<Cps[it].n;i++)
fscanf(fd,"%g",&Cps[it].ms[i]);
else
fgets(line, sizeof(line), fd);
for (i=0;i<Cps[it].n+1;i++)
fscanf(fd,"%g",&Cps[it].as[i]);
fgets(line, sizeof(line), fd);
/* Parameters for SNII Rates */
/* #### SNII Parameters #### */
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%g \n",&Cps[it].SNII_Mmin);
fgets(line, sizeof(line), fd);
/* Parameters for SNIa Rates */
/* #### SNIa Parameters #### */
fgets(line, sizeof(line), fd);
fscanf(fd, "%g %g\n",&Cps[it].SNIa_Mpl,&Cps[it].SNIa_Mpu);
fscanf(fd, "%g \n",&Cps[it].SNIa_a);
fscanf(fd, "%g %g %g\n",&Cps[it].SNIa_Mdl1,&Cps[it].SNIa_Mdu1,&Cps[it].SNIa_bb1);
fscanf(fd, "%g %g %g\n",&Cps[it].SNIa_Mdl2,&Cps[it].SNIa_Mdu2,&Cps[it].SNIa_bb2);
fgets(line, sizeof(line), fd);
/* Metal injection SNII */
/* #### Metal Parameters ####*/
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%d %d\n",&Cps[it].nptsSNII,&Cps[it].nelts);
/* allocate memory for elts */
if (Cps[it].nptsSNII<=MAXDATASIZE)
{
Elts[it] = malloc((Cps[it].nelts+2) * sizeof(struct local_elts_chimie));
}
else
{
printf("\n Cps[it].nptsSNII = %d > MAXDATASIZE = %d !!!\n\n",Cps[it].nptsSNII,MAXDATASIZE);
endrun(88800);
}
/* allocate memory */
MassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); /* really needed ? */
MassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double));
MassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
EjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleEjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double));
/* injected metals */
for (i=0;i<Cps[it].nelts+2;i++)
{
fgets(line, sizeof(line), fd);
/* strip trailing line */
for (j = 0; j < strlen(line); j++)
if ( line[j] == '\n' || line[j] == '\r' )
line[j] = '\0';
/* copy labels */
strcpy(Elts[it][i].label,line);
/* probleme */
strcpy(buffer,&Elts[it][i].label[2]);
strcpy(Elts[it][i].label,buffer);
fgets(line, sizeof(line), fd);
fscanf(fd, "%g %g\n",&Elts[it][i].MminSNII,&Elts[it][i].StepSNII);
for (j=0;j<Cps[it].nptsSNII;j++)
{
fscanf(fd, "%g\n",&Elts[it][i].MetalSNII[j]);
}
}
/* integrals of injected metals */
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%d %d\n",&Cps[it].nptsSNII,&Cps[it].nelts);
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
/* integrals of injected metals */
for (i=0;i<Cps[it].nelts+2;i++)
{
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%g %g\n",&Elts[it][i].MminSNII,&Elts[it][i].StepSNII);
for (j=0;j<Cps[it].nptsSNII;j++)
{
fscanf(fd, "%g\n",&Elts[it][i].ArraySNII[j]);
}
}
/* Metal injection SNIa */
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%g\n",&Cps[it].Mco);
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
int nelts;
char label[72];
fscanf(fd, "%d\n",&nelts);
/* check */
if (nelts != Cps[it].nelts)
{
printf("\nThe number of elements in SNII (=%d) is not identical to the on of SNIa (=%d) !!!\n\n",Cps[it].nelts,nelts);
printf("This is not supported by the current implementation !!!\n");
endrun(88805);
}
for (i=0;i<Cps[it].nelts+2;i++)
{
fgets(line, sizeof(line), fd); /* label */
/* check label */
/* strip trailing line */
for (j = 0; j < strlen(line); j++)
if ( line[j] == '\n' || line[j] == '\r' )
line[j] = '\0';
strcpy(label,line);
strcpy(buffer,&label[2]);
strcpy(label,buffer);
if (strcmp(label,Elts[it][i].label)!=0)
{
printf("\nLabel of SNII element %d (=%s) is different from the SNIa one (=%s) !!!\n\n",i,Elts[it][i].label,label);
endrun(88806);
}
//fgets(line, sizeof(line), fd);
fscanf(fd, "%g\n",&Elts[it][i].MSNIa);
}
/* Solar Mass Abundances */
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fgets(line, sizeof(line), fd);
fscanf(fd, "%d\n",&nelts);
/* check */
if (nelts != Cps[it].nelts)
{
printf("\nThe number of elements in SolarMassAbundances (=%d) is not identical to the on of SNIa (=%d) !!!\n\n",Cps[it].nelts,nelts);
printf("This is not supported by the current implementation !!!\n");
endrun(88805);
}
for (i=0;i<Cps[it].nelts;i++)
{
fgets(line, sizeof(line), fd); /* label */
/* check label */
/* strip trailing line */
for (j = 0; j < strlen(line); j++)
if ( line[j] == '\n' || line[j] == '\r' )
line[j] = '\0';
strcpy(label,line);
strcpy(buffer,&label[2]);
strcpy(label,buffer);
if (strcmp(label,Elts[it][i+2].label)!=0)
{
printf("\nLabel of SNII element %d (=%s) is different from the SNIa one (=%s) !!!\n\n",i,Elts[it][i+2].label,label);
endrun(88806);
}
//fgets(line, sizeof(line), fd);
fscanf(fd, "%g\n",&Elts[it][i+2].SolarMassAbundance);
}
fclose(fd);
if (verbose && ThisTask==0)
info(it);
}
/*! This function returns the mass fraction of a star of mass m
* using the current IMF
*/
static double get_imf(double m)
{
int i;
int n;
n = Cp->n;
/* convert m in msol */
m = m*All.CMUtoMsol;
if (n==0)
return Cp->bs[0]* pow(m,Cp->as[0]);
else
{
for (i=0;i<n;i++)
if (m < Cp->ms[i])
return Cp->bs[i]* pow(m,Cp->as[i]);
return Cp->bs[n]* pow(m,Cp->as[n]);
}
}
/*! This function returns the mass fraction between m1 and m2
* per mass unit, using the current IMF
*/
static double get_imf_M(double m1, double m2)
{
int i;
int n;
double p;
double integral=0;
double mmin,mmax;
n = Cp->n;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
if (n==0)
{
p = Cp->as[0]+1;
integral = (Cp->bs[0]/p) * ( pow(m2,p) - pow(m1,p) );
//printf("--> %g %g %g %g int=%g\n",m1,m2,pow(m2,p), pow(m1,p),integral);
}
else
{
integral = 0;
/* first */
if (m1<Cp->ms[0])
{
mmin = m1;
mmax = dmin(Cp->ms[0],m2);
p = Cp->as[0] + 1;
integral += (Cp->bs[0]/p) * ( pow(mmax,p) - pow(mmin,p) );
}
/* last */
if (m2>Cp->ms[n-1])
{
mmin = dmax(Cp->ms[n-1],m1);
mmax = m2;
p = Cp->as[n] + 1;
integral += (Cp->bs[n]/p) * ( pow(mmax,p) - pow(mmin,p) );
}
/* loop over other segments */
for (i=0;i<n-1;i++)
{
mmin = dmax(Cp->ms[i ],m1);
mmax = dmin(Cp->ms[i+1],m2);
if (mmin<mmax)
{
p = Cp->as[i+1] + 1;
integral += (Cp->bs[i+1]/p) * ( pow(mmax,p) - pow(mmin,p) );
}
}
}
return integral;
}
/*! This function returns the number fraction between m1 and m2
* per mass unit, using the current IMF
*/
static double get_imf_N(double m1, double m2)
{
int i;
int n;
double p;
double integral=0;
double mmin,mmax;
n = Cp->n;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
if (n==0)
{
p = Cp->as[0];
integral = (Cp->bs[0]/p) * ( pow(m2,p) - pow(m1,p) );
}
else
{
integral = 0;
/* first */
if (m1<Cp->ms[0])
{
mmin = m1;
mmax = dmin(Cp->ms[0],m2);
p = Cp->as[0];
integral += (Cp->bs[0]/p) * ( pow(mmax,p) - pow(mmin,p) );
}
/* last */
if (m2>Cp->ms[n-1])
{
mmin = dmax(Cp->ms[n-1],m1);
mmax = m2;
p = Cp->as[n];
integral += (Cp->bs[n]/p) * ( pow(mmax,p) - pow(mmin,p) );
}
/* loop over other segments */
for (i=0;i<n-1;i++)
{
mmin = dmax(Cp->ms[i ],m1);
mmax = dmin(Cp->ms[i+1],m2);
if (mmin<mmax)
{
p = Cp->as[i+1];
integral += (Cp->bs[i+1]/p) * ( pow(mmax,p) - pow(mmin,p) );
}
}
}
/* convert into mass unit mass unit */
integral = integral *All.CMUtoMsol;
return integral;
}
/*! Sample the imf using monte carlo approach
*/
static double imf_sampling()
{
int i;
int n;
double m;
double f;
double pmin,pmax;
n = Cp->n;
/* init random */
//srandom(irand);
f = (double)random()/(double)RAND_MAX;
if (n==0)
{
pmin = pow(Cp->Mmin,Cp->as[0]);
pmax = pow(Cp->Mmax,Cp->as[0]);
m = pow(f*(pmax - pmin) + pmin ,1./Cp->as[0]);
return m* All.MsoltoCMU;
}
else
{
if (f<Cp->fs[0])
{
pmin = pow(Cp->Mmin ,Cp->as[0]);
m = pow(Cp->imf_Ntot*Cp->as[0]/Cp->bs[0]* (f-0) + pmin ,1./Cp->as[0]);
return m* All.MsoltoCMU;
}
for (i=0;i<n-1;i++)
{
if (f<Cp->fs[i+1])
{
pmin = pow(Cp->ms[i] ,Cp->as[i+1]);
m = pow(Cp->imf_Ntot*Cp->as[i+1]/Cp->bs[i+1]* (f-Cp->fs[i]) + pmin ,1./Cp->as[i+1]);
return m* All.MsoltoCMU;
}
}
/* last portion */
pmin = pow(Cp->ms[n-1] ,Cp->as[n]);
m = pow(Cp->imf_Ntot*Cp->as[n]/Cp->bs[n]* (f-Cp->fs[n-1]) + pmin ,1./Cp->as[n]);
return m* All.MsoltoCMU;
}
}
static double imf_sampling_from_random(double f)
{
int i;
int n;
double m;
double pmin,pmax;
n = Cp->n;
if (n==0)
{
pmin = pow(Cp->Mmin,Cp->as[0]);
pmax = pow(Cp->Mmax,Cp->as[0]);
m = pow(f*(pmax - pmin) + pmin ,1./Cp->as[0]);
return m* All.MsoltoCMU;
}
else
{
if (f<Cp->fs[0])
{
pmin = pow(Cp->Mmin ,Cp->as[0]);
m = pow(Cp->imf_Ntot*Cp->as[0]/Cp->bs[0]* (f-0) + pmin ,1./Cp->as[0]);
return m* All.MsoltoCMU;
}
for (i=0;i<n-1;i++)
{
if (f<Cp->fs[i+1])
{
pmin = pow(Cp->ms[i] ,Cp->as[i+1]);
m = pow(Cp->imf_Ntot*Cp->as[i+1]/Cp->bs[i+1]* (f-Cp->fs[i]) + pmin ,1./Cp->as[i+1]);
return m* All.MsoltoCMU;
}
}
/* last portion */
pmin = pow(Cp->ms[n-1] ,Cp->as[n]);
m = pow(Cp->imf_Ntot*Cp->as[n]/Cp->bs[n]* (f-Cp->fs[n-1]) + pmin ,1./Cp->as[n]);
return m* All.MsoltoCMU;
}
}
/*! This function initializes the imf parameters
defined in the chemistry file
*/
void init_imf(void)
{
float integral = 0;
float p;
float cte;
int i,n;
double mmin,mmax;
n = Cp->n;
if (n==0)
{
p = Cp->as[0]+1;
integral = integral + ( pow(Cp->Mmax,p)-pow(Cp->Mmin,p))/(p) ;
Cp->bs[0] = 1./integral ;
}
else
{
cte = 1.0;
if (Cp->Mmin < Cp->ms[0])
{
p = Cp->as[0]+1;
integral = integral + (pow(Cp->ms[0],p) - pow(Cp->Mmin,p))/p;
}
for (i=0;i<n-1;i++)
{
cte = cte* pow( Cp->ms[i],( Cp->as[i] - Cp->as[i+1] ));
p = Cp->as[i+1]+1;
integral = integral + cte*(pow(Cp->ms[i+1],p) - pow(Cp->ms[i],p))/p;
}
if (Cp->Mmax > Cp->ms[-1])
{
cte = cte* pow( Cp->ms[n-1] , ( Cp->as[n-1] - Cp->as[n] ) );
p = Cp->as[n]+1;
integral = integral + cte*(pow(Cp->Mmax,p) - pow(Cp->ms[n-1],p))/p;
}
/* compute all b */
Cp->bs[0] = 1./integral;
for (i=0;i<n;i++)
{
Cp->bs[i+1] = Cp->bs[i] * pow( Cp->ms[i],( Cp->as[i] - Cp->as[i+1] ));
}
}
//if (verbose && ThisTask==0)
// {
// printf("-- bs -- \n");
// for (i=0;i<n+1;i++)
// printf("%g ",Cp->bs[i]);
// printf("\n");
// }
mmin = Cp->Mmin *All.MsoltoCMU; /* in CMU */
mmax = Cp->Mmax *All.MsoltoCMU; /* in CMU */
Cp->imf_Ntot = get_imf_N(mmin,mmax) *All.MsoltoCMU; /* in CMU */
/* init fs : mass fraction at ms */
if (n>0)
{
for (i=0;i<n+1;i++)
{
mmax = Cp->ms[i] *All.MsoltoCMU; /* in CMU */
Cp->fs[i] = All.MsoltoCMU*get_imf_N(mmin,mmax)/Cp->imf_Ntot;
}
}
}
#ifdef CHIMIE_OPTIMAL_SAMPLING
/*
integral of m\xi with the kroupa normalisation ( in k)
input : mass in Msol
*/
double imf_int_mxi(double left,double right, double k)
{
return get_imf_M(left*All.MsoltoCMU,right*All.MsoltoCMU)/(Cp->bs[0])*k;
}
/*
integral of \xi with the kroupa normalisation ( in k)
input : mass in Msol
*/
double imf_int_xi(double left,double right, double k)
{
return get_imf_N(left*All.MsoltoCMU,right*All.MsoltoCMU)/ (Cp->bs[0]*All.CMUtoMsol) * k;
}
/*
compute the normalisation and the maximum stellar mass
i.e, Cp->k and Cp->m_max
input : mass in Msol
output : mass in Msol
*/
double optimal_init_norm(double Mecl)
{
double a,b,c;
double mb;
double Mmax;
Mmax = Cp->Mmax;
Cp->k = 1.; /* kroupa normalisation */
//we use 150Msun for the maximum possible stellar mass (TODO: put it as a parameter)
a = Cp->Mmin;
c = Mmax;
b = (c+a)/2.;
//solve the two equations to find m_max and k iteratively
while(((c/b)-(a/b)) > 0.00001)
{
mb = imf_int_mxi(Cp->Mmin,b,Cp->k)/imf_int_xi(b,Mmax,Cp->k)+b;
if(mb < Mecl)
a = b;
else
c = b;
b = (c+a)/2;
}
//store the results
Cp->m_max = b;
Cp->k = (Mecl-Cp->m_max) / imf_int_mxi(Cp->Mmin,Cp->m_max,Cp->k);
return Cp->m_max;
}
/*
compute the normalisation and the maximum stellar mass
for one single stellar particle
output in StP[m].OptIMF_CurrentMass, in Msol
*/
double init_optimal_imf_sampling(int i)
{
double a,b,c;
double mb;
double Mmax;
double Mecl;
int m;
m = P[i].StPIdx;
Mmax = Cp->Mmax;
Mecl = StP[m].InitialMass*All.CMUtoMsol; /* to Msol */
StP[m].OptIMF_k = 1.; /* kroupa normalisation */
StP[m].OptIMF_N_WD = 0;
//we use 150Msun for the maximum possible stellar mass (TODO: put it as a parameter)
a = Cp->Mmin;
c = Mmax;
b = (c+a)/2.;
//solve the two equations to find m_max and k iteratively
while(((c/b)-(a/b)) > 0.00001)
{
mb = imf_int_mxi(Cp->Mmin,b,StP[m].OptIMF_k)/imf_int_xi(b,Mmax,StP[m].OptIMF_k)+b;
if(mb < Mecl)
a = b;
else
c = b;
b = (c+a)/2;
}
//store the results
StP[m].OptIMF_m_max = b;
StP[m].OptIMF_k = (Mecl-StP[m].OptIMF_m_max) / imf_int_mxi(Cp->Mmin,StP[m].OptIMF_m_max,StP[m].OptIMF_k);
StP[m].OptIMF_CurrentMass = StP[m].OptIMF_m_max;
return StP[m].OptIMF_m_max;
}
/*
compute the next mass given the previous
input : mass in Msol
output : mass in Msol
*/
double optimal_get_next_mass(double m){
double a,b,c,mb;
const double err = 1.e-8;
a = Cp->Mmin;
c = m;
b = (c+a)/2;
mb = imf_int_mxi(b,m,Cp->k);
while(fabs((mb-b)/mb) > err)
{
mb = imf_int_mxi(b,m,Cp->k);
if(mb < b)
c = b;
else
a = b;
b = (c+a)/2;
}
return b;
}
/*
comute the next mass given the previous,
according to the imf of a given particle
output : mass in Msol
*/
double optimal_get_next_mass_for_one_particle(int j){
double a,b,c,mb;
const double err = 1.e-8;
double m;
m = StP[j].OptIMF_CurrentMass;
a = Cp->Mmin;
c = m;
b = (c+a)/2;
mb = imf_int_mxi(b,m,StP[j].OptIMF_k);
while(fabs((mb-b)/mb) > err)
{
mb = imf_int_mxi(b,m,StP[j].OptIMF_k);
if(mb < b)
c = b;
else
a = b;
b = (c+a)/2;
}
return b;
}
/*
stop the imf loop
input : mass in Msol
*/
int optimal_stop_loop(double m)
{
int bool;
if(imf_int_mxi(Cp->Mmin,m,Cp->k) < Cp->Mmin)
bool = 1;
else
bool = 0;
return bool;
}
/*
find m1 such as the stellar mass, according to the IMF, between m1 and m2 is mp.
input : masses in Msol
*/
double optimal_get_m1_from_m2(double m2, double mp)
{
double a,b,c;
double m1,mb;
a = Cp->Mmin;
c = m2;
b = (c+a)/2.;
//solve the two equations to find m_max and k iteratively
while(((c/b)-(a/b)) > 0.00001)
{
mb = imf_int_mxi(b,m2,Cp->k);
if(mb < mp)
c = b;
else
a = b;
b = (c+a)/2;
}
//store the results
m1 = b;
return m1;
}
float optimal_get_k_normalisation_factor()
{
return Cp->k;
}
#endif /* CHIMIE_OPTIMAL_SAMPLING */
/*! This function initializes the chemistry parameters
*/
void init_chimie(void)
{
int i,nf;
double u_lt;
double UnitLength_in_kpc;
double UnitMass_in_Msol;
char filename[500];
char ext[100];
/* check some flags */
#ifndef COSMICTIME
if (All.ComovingIntegrationOn)
{
if(ThisTask == 0)
printf("Code wasn't compiled with COSMICTIME support enabled!\n");
endrun(-88800);
}
#endif
#ifdef CHIMIE_AGB
if(ThisTask == 0)
printf("AGB is enabled\n");
#else
if(ThisTask == 0)
printf("AGB is disabled\n");
#endif
UnitLength_in_kpc = All.UnitLength_in_cm / KPC_IN_CM;
UnitMass_in_Msol = All.UnitMass_in_g / SOLAR_MASS;
//u_lt = -log10( 4.7287e11*sqrt(pow(UnitLength_in_kpc,3)/UnitMass_in_Msol));
/*Sat Dec 25 23:27:10 CET 2010 */
u_lt = -log10(All.UnitTime_in_Megayears*1e6);
allocate_chimie();
for (nf=0;nf<All.ChimieNumberOfParameterFiles;nf++)
{
if (All.ChimieNumberOfParameterFiles==1)
sprintf(filename,"%s",All.ChimieParameterFile);
else
sprintf(filename,"%s.%d",All.ChimieParameterFile,nf);
/* check the file extension and read */
#ifdef CHIMIE_FROM_HDF5
if ( (strcmp("h5",get_filename_ext(filename))==0) || (strcmp("hdf5",get_filename_ext(filename))==0))
read_chimie_h5(filename,nf);
else
#endif //CHIMIE_FROM_HDF5
read_chimie(filename,nf);
/* set the table */
set_table(nf);
/* Conversion into program time unit */
Cp->coeff_z[2][2] = Cp->coeff_z[2][2] + u_lt;
for (i=0;i<3;i++)
Cp->coeff_z[1][i] = Cp->coeff_z[1][i]/2.0;
/* init imf parameters */
init_imf();
/* init SNII parameters */
//if (Cp->n==0)
// {
// //Cp->SNII_cte[0] = Cp->bs[0]/Cp->as[0];
// Cp->SNII_cte = Cp->bs[0]/Cp->as[0];
// Cp->SNII_a = Cp->as[0];
// }
//else
// {
// //for (i=0;i<Cp->n+1;i++) /* if multiple power law in the SNII mass range */
// // Cp->SNII_cte[i] = Cp->bs[i]/Cp->as[i];
// Cp->SNII_cte = Cp->bs[Cp->n]/Cp->as[Cp->n];
// Cp->SNII_a = Cp->as[Cp->n];
// }
/* init DYIN parameters */
//if (Cp->n==0)
// {
// //Cp->DYIN_cte[0] = Cp->bs[0]/Cp->as[0];
// Cp->DYIN_cte = Cp->bs[0]/Cp->as[0];
// Cp->DYIN_a = Cp->as[0];
// }
//else
// {
// //for (i=0;i<Cp->n+1;i++) /* if multiple power law in the DYIN mass range */
// // Cp->DYIN_cte[i] = Cp->bs[i]/Cp->as[i];
// Cp->DYIN_cte = Cp->bs[Cp->n]/Cp->as[Cp->n];
// Cp->DYIN_a = Cp->as[Cp->n];
// }
/* init SNIa parameters */
Cp->SNIa_a1 = Cp->SNIa_a;
Cp->SNIa_b1 = (Cp->SNIa_a1+1)/(pow(Cp->SNIa_Mdu1,Cp->SNIa_a1+1)-pow(Cp->SNIa_Mdl1,Cp->SNIa_a1+1));
Cp->SNIa_cte1 = Cp->SNIa_b1/Cp->SNIa_a1;
Cp->SNIa_a2 = Cp->SNIa_a;
Cp->SNIa_b2 = (Cp->SNIa_a2+1)/(pow(Cp->SNIa_Mdu2,Cp->SNIa_a2+1)-pow(Cp->SNIa_Mdl2,Cp->SNIa_a2+1));
Cp->SNIa_cte2 = Cp->SNIa_b2/Cp->SNIa_a2;
/* init SNIa parameters */
if (Cp->n==0)
{
Cp->SNIa_cte = Cp->bs[0]/Cp->as[0];
Cp->SNIa_a = Cp->as[0];
}
else
{
Cp->SNIa_cte = Cp->bs[Cp->n]/Cp->as[Cp->n];
Cp->SNIa_a = Cp->as[Cp->n];
}
//for (i=0;i<Cp->nelts+2;i++)
// Elt[i].MminSNII = log10(Elt[i].MminSNII);
//
//for (i=0;i<Cp->nelts+2;i++)
// Elt[i].MminDYIN = log10(Elt[i].MminDYIN);
/* output info */
//if (verbose && ThisTask==0)
// {
// //printf("-- SNII_cte -- \n");
// //for (i=0;i<Cp->n+1;i++)
// // printf("%g ",Cp->SNII_cte[i]);
// //printf("%g ",Cp->SNII_cte);
//
// printf("\n");
//}
/* output info */
//if (verbose && ThisTask==0)
// {
// printf("-- DYIN_cte -- \n");
// //for (i=0;i<Cp->n+1;i++)
// // printf("%g ",Cp->DYIN_cte[i]);
// //printf("%g ",Cp->DYIN_cte);
//
// printf("\n");
// }
/* check that the masses are higher than the last IMF elbow */
if (Cp->n>0)
{
if (Cp->SNIa_Mpl < Cp->ms[Cp->n-1])
{
printf("\nSNIa_Mpl = %g < ms[n-1] = %g !!!\n\n",Cp->SNIa_Mpl,Cp->ms[Cp->n-1]);
//printf("This is not supported by the current implementation !!!\n");
//endrun(88801);
printf("We assume that there will be no SNIa !!!\n");
Cp->SNIa_Mpl = -1;
Cp->SNIa_Mpu = -1;
}
if (Cp->SNIa_Mpu < Cp->ms[Cp->n-1])
{
printf("\nSNIa_Mpu = %g < ms[n-1] = %g !!!\n\n",Cp->SNIa_Mpu,Cp->ms[Cp->n-1]);
//printf("This is not supported by the current implementation !!!\n");
//endrun(88802);
printf("We assume that there will be no SNIa !!!\n");
Cp->SNIa_Mpl = -1;
Cp->SNIa_Mpu = -1;
}
if (Cp->SNII_Mmin < Cp->ms[Cp->n-1])
{
printf("\nSNII_Mmin = %g < ms[n-1] = %g !!!\n\n",Cp->SNII_Mmin,Cp->ms[Cp->n-1]);
printf("This is not supported by the current implementation !!!\n");
endrun(88803);
}
if (Cp->SNII_Mmax < Cp->ms[Cp->n-1])
{
printf("\nSNII_Mmax = %g < ms[n-1] = %g !!!\n\n",Cp->SNII_Mmax,Cp->ms[Cp->n-1]);
printf("This is not supported by the current implementation !!!\n");
endrun(88804);
}
if (Cp->DYIN_Mmin < Cp->ms[Cp->n-1])
{
printf("\nDYIN_Mmin = %g < ms[n-1] = %g !!!\n\n",Cp->DYIN_Mmin,Cp->ms[Cp->n-1]);
printf("This is not supported by the current implementation !!!\n");
endrun(88805);
}
if (Cp->DYIN_Mmax < Cp->ms[Cp->n-1])
{
printf("\nDYIN_Mmax = %g < ms[n-1] = %g !!!\n\n",Cp->DYIN_Mmax,Cp->ms[Cp->n-1]);
printf("This is not supported by the current implementation !!!\n");
endrun(88806);
}
}
/* number of WD per solar mass */
Cp->SNIa_NWDN = Cp->SNIa_cte * ( pow(Cp->SNIa_Mpu,Cp->SNIa_a)-pow(Cp->SNIa_Mpl,Cp->SNIa_a) );
}
/*
init short cuts to some elements
*/
FE = get_ElementIndex("Fe");
METALS = get_ElementIndex("Metals");
/* write chime stat header */
#ifdef CHIMIE_STATS
int k;
fprintf(FdChimieStatsSNs,"# ThisTask TimeStep Time ID Mass NSNII NSNIa MassSN x y z ");
fprintf(FdChimieStatsSNs,"EjectedGasMass ");
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsSNs,"EjectedMass_%s ",&Elt[k+2].label);
fprintf(FdChimieStatsSNs,"\n");
fflush(FdChimieStatsSNs);
fprintf(FdChimieStatsGas,"# ThisTask TimeStep Time ID IDStar aij x y z ");
fprintf(FdChimieStatsGas,"Mass_0 ");
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsGas,"%s_0 ",&Elt[k+2].label);
fprintf(FdChimieStatsGas,"Mass_1 ");
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsGas,"%s_1 ",&Elt[k+2].label);
fprintf(FdChimieStatsGas,"\n");
fflush(FdChimieStatsGas);
#endif
}
/*! This function performe simple checks
* to validate the chemistry initialization
*/
void check_chimie(void)
{
int i;
if (ThisTask==0)
{
printf("Number of elts : %d\n",Cp->nelts);
for(i=2;i<Cp->nelts+2;i++)
printf("%s ",&Elt[i].label);
printf("\n");
}
/* check number of elements */
if (NELEMENTS != Cp->nelts)
{
printf("(Taks=%d) NELEMENTS (=%d) != Cp->nelts (=%d) : please check !!!\n\n",ThisTask,NELEMENTS,Cp->nelts);
endrun(88807);
}
/* check that iron is the first element */
//if ((strcmp("Fe",Elt[2].label))!=0)
// {
// printf("(Taks=%d) first element (=%s) is not %s !!!\n\n",ThisTask,Elt[2].label,FIRST_ELEMENT);
// endrun(88808);
// }
}
/*! This function print some info on the chimie parameters
*/
void info(int it)
{
int i,j;
printf("\nTable %d\n\n",it);
printf("%g %g %g\n", Cps[it].coeff_z[0][0],Cps[it].coeff_z[0][1],Cps[it].coeff_z[0][2]);
printf("%g %g %g\n", Cps[it].coeff_z[1][0],Cps[it].coeff_z[1][1],Cps[it].coeff_z[1][2]);
printf("%g %g %g\n", Cps[it].coeff_z[2][0],Cps[it].coeff_z[2][1],Cps[it].coeff_z[2][2]);
printf("\n");
printf("\nIMF\n");
printf("%g %g\n",Cps[it].Mmin,Cps[it].Mmax);
printf("%d\n",Cps[it].n);
for (i=0;i<Cps[it].n;i++)
printf( "ms : %g ",Cps[it].ms[i]);
printf("\n");
for (i=0;i<Cps[it].n+1;i++)
printf( "as : %g ",Cps[it].as[i]);
printf("\n");
printf("\nRate SNII\n");
printf("%g ",Cps[it].SNII_Mmin);
printf("%g ",Cps[it].SNII_Mmax);
printf("\n");
printf("\nRate DYIN\n");
printf("%g ",Cps[it].DYIN_Mmin);
printf("%g ",Cps[it].DYIN_Mmax);
printf("\n");
printf("\nRate SNIa\n");
printf("%g %g\n",Cps[it].SNIa_Mpl,Cps[it].SNIa_Mpu);
printf("%g \n",Cps[it].SNIa_a);
printf("%g %g %g\n",Cps[it].SNIa_Mdl1,Cps[it].SNIa_Mdu1,Cps[it].SNIa_b1);
printf("%g %g %g\n",Cps[it].SNIa_Mdl2,Cps[it].SNIa_Mdu2,Cps[it].SNIa_b2);
printf("\n");
for (i=0;i<Cps[it].nelts+2;i++)
{
printf("> %g %g\n",Elts[it][i].MminSNII,Elts[it][i].StepSNII);
for (j=0;j<Cps[it].nptsSNII;j++)
{
printf(" %g\n",Elts[it][i].ArraySNII[j]);
}
}
for (i=0;i<Cps[it].nelts+2;i++)
{
printf("> %g %g\n",Elts[it][i].MminDYIN,Elts[it][i].StepDYIN);
for (j=0;j<Cps[it].nptsDYIN;j++)
{
printf(" %g\n",Elts[it][i].ArrayDYIN[j]);
}
}
printf("\n");
printf("%g\n",Cps[it].Mco);
for (i=0;i<Cps[it].nelts+2;i++)
printf("%g\n",Elts[it][i].MSNIa);
printf("\n");
}
/*! Return the number of elements considered
*/
int get_nelts()
{
return Cp->nelts;
}
/*! Return the solar mass abundance of elt i
*/
float get_SolarMassAbundance(i)
{
return Elt[i+2].SolarMassAbundance;
}
/*! Return the label of element i
*/
char* get_Element(i)
{
return Elt[i+2].label;
}
int get_ElementIndex(char *elt)
{
/*
example:
get_SolarMassAbundance(get_ElementIndex("Fe"));
*/
int i=-1;
for(i=0;i<Cp->nelts;i++)
if( strcmp(Elt[i+2].label, elt) == 0 )
return i;
printf("element %s is unknown\n",elt);
endrun(777123);
return i;
}
/*! Return the lifetime of a star of mass m and metallicity z
*/
double star_lifetime(double z,double m)
{
/* z is the mass fraction of metals, ie, the metallicity */
/* m is the stellar mass in code unit */
/* Return t in code time unit */
int i;
double a,b,c;
double coeff[3];
double logm,twologm,logm2,time;
/* convert m in msol */
m = m*All.CMUtoMsol;
for (i=0;i<3;i++)
coeff[i] = ( Cp->coeff_z[i][0]*z+Cp->coeff_z[i][1] )*z+Cp->coeff_z[i][2];
a = coeff[0];
b = coeff[1];
c = coeff[2];
logm = log10(m);
twologm = 2.0 * logm;
logm2 = logm*logm;
time = pow(10.,(a*logm2+b*twologm+c));
time = time * All.HubbleParam; /* correct from hubble as not done in coeff_z */
return time;
}
/*! Return the mass of a star having a livetime t and a metallicity z
*/
double star_mass_from_age(double z,double t)
{
/* z is the mass fraction of metals, ie, the metallicity */
/* t is the star life time (in code unit) */
/* return the stellar mass (in code unit) that has a lifetime equal to t */
/* this is the inverse of star_lifetime */
int i;
double a,b,c;
double coeff[3];
double m;
t = t/All.HubbleParam; /* correct from hubble as not done in coeff_z */
for (i=0;i<3;i++)
coeff[i] = ( Cp->coeff_z[i][0]*z+Cp->coeff_z[i][1] )*z+Cp->coeff_z[i][2];
a = coeff[0];
b = coeff[1];
c = coeff[2];
m = -(b+sqrt(b*b-a*(c-log10(t))))/a;
m = pow(10,m); /* here, m is in solar mass */
m = m*All.MsoltoCMU; /* Msol to mass unit */
return m;
}
/****************************************************************************************/
/*
/* Supernova rate : number of supernova per mass unit
/*
/****************************************************************************************/
double DYIN_rate(double m1,double m2)
{
/*
compute the number of stars between m1 and m2
masses in code unit
*/
double RDYIN;
double md,mu;
RDYIN = 0.0;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
/* find md, mu */
md = dmax(m1,Cp->DYIN_Mmin);
mu = dmin(m2,Cp->DYIN_Mmax);
if (mu<=md) /* no dying stars in that mass range */
return 0.0;
/* to code units */
md = md * All.MsoltoCMU;
mu = mu * All.MsoltoCMU;
/* compute number in the mass range */
RDYIN = get_imf_N(md,mu);
return RDYIN;
}
double SNII_rate(double m1,double m2)
{
/*
compute the number of SNII between m1 and m2
masses in code unit
*/
double RSNII;
double md,mu;
RSNII = 0.0;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
/* (1) find md, mu */
md = dmax(m1,Cp->SNII_Mmin);
mu = dmin(m2,Cp->SNII_Mmax);
if (mu<=md) /* no SNII in that mass range */
return 0.0;
/* !!!!! here we should use get_imf_N !!!! */
/* to ensure the full imf */
//RSNII = Cp->SNII_cte * (pow(mu,Cp->SNII_a)-pow(md,Cp->SNII_a)); /* number per solar mass */
///* convert in number per solar mass to number per mass unit */
//RSNII = RSNII *All.CMUtoMsol;
/* to code units */
md = md * All.MsoltoCMU;
mu = mu * All.MsoltoCMU;
/* compute number in the mass range */
RSNII = get_imf_N(md,mu);
return RSNII;
}
double SNIa_rate(double m1,double m2)
{
/*
compute the number of SNIa between m1 and m2
masses in code unit
*/
double RSNIa;
double md,mu;
RSNIa = 0.0;
if (Cp->SNIa_Mpl<0)
return RSNIa;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
/* RG contribution */
md = dmax(m1,Cp->SNIa_Mdl1);
mu = dmin(m2,Cp->SNIa_Mdu1);
if (md<mu)
RSNIa = RSNIa + Cp->SNIa_bb1 * Cp->SNIa_cte1 * (pow(mu,Cp->SNIa_a1)-pow(md,Cp->SNIa_a1));
/* MS contribution */
md = dmax(m1,Cp->SNIa_Mdl2);
mu = dmin(m2,Cp->SNIa_Mdu2);
if (md<mu)
RSNIa = RSNIa + Cp->SNIa_bb2 * Cp->SNIa_cte2 * (pow(mu,Cp->SNIa_a2)-pow(md,Cp->SNIa_a2));
/* WD contribution */
md = dmax(m1,Cp->SNIa_Mpl); /* select stars that have finished their life -> WD */
mu = Cp->SNIa_Mpu; /* no upper bond */
if (mu<=md) /* no SNIa in that mass range */
return 0.0;
RSNIa = RSNIa * Cp->SNIa_cte * (pow(mu,Cp->SNIa_a)-pow(md,Cp->SNIa_a)); /* number per solar mass */
/* convert in number per solar mass to number per mass unit */
RSNIa = RSNIa *All.CMUtoMsol;
return RSNIa;
}
double SNIa_single_rate(double m)
{
/*
compute the number of SNIa per mass unit (in CMU) if the mass of the star is m
m is given in code units.
NOTE : as single stars are follow a normal imf distribution and here
the RG and MS follow another IMF, we need to introduce the normalisation factor
*/
double RSNIa;
RSNIa = 0.0;
if (Cp->SNIa_Mpl<0)
return RSNIa;
/* convert m in msol */
m = m*All.CMUtoMsol;
/* the star mass is greater than the WD interval */
if (m <= Cp->SNIa_Mpl )
{
/* according to its mass the star is a RG */
if ( ( m > Cp->SNIa_Mdl1 ) && ( m <= Cp->SNIa_Mdu1 ) )
RSNIa = Cp->SNIa_bb1 * Cp->SNIa_NWDN * Cp->SNIa_b1 * pow(m,Cp->SNIa_a1 ) / get_imf(m*All.MsoltoCMU);
/* according to its mass the star is a MS */
if ( ( m > Cp->SNIa_Mdl2 ) && ( m <= Cp->SNIa_Mdu2 ) )
RSNIa = Cp->SNIa_bb2 * Cp->SNIa_NWDN * Cp->SNIa_b2 * pow(m,Cp->SNIa_a2 ) / get_imf(m*All.MsoltoCMU);
}
return RSNIa;
}
void DYIN_mass_ejection(double m1,double m2)
{
/*
Compute the mass fraction and yields of dying stars with masses between m1 and m2.
Store the result in the global variable`` MassFracDYIN``::
MassFracDYIN[0] = total gas
MassFracDYIN[1] = helium core (i.e. alpha(m))
MassFracDYIN[i] = frac mass elt i.
*/
double l1,l2;
int i1,i2,i1p,i2p,j;
double f1,f2;
double v1,v2;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
/* this was not in Poirier... */
m1 = dmax(m1,Cp->DYIN_Mmin);
m2 = dmin(m2,Cp->DYIN_Mmax);
if ( m2<=m1 )
{
for (j=0;j<Cp->nelts+2;j++)
MassFracDYIN[j] = 0;
return;
}
j = 0;
l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ;
l2 = ( log10(m2) - Elt[j].MminDYIN) / Elt[j].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<0) i1=0;
if (i2<0) i2=0;
/* --------- TOTAL GAS ---------- */
j = 0;
v1 = f1 * ( Elt[j].ArrayDYIN[i1p] - Elt[j].ArrayDYIN[i1] ) + Elt[j].ArrayDYIN[i1];
v2 = f2 * ( Elt[j].ArrayDYIN[i2p] - Elt[j].ArrayDYIN[i2] ) + Elt[j].ArrayDYIN[i2];
MassFracDYIN[j] = v2-v1;
/* --------- He core therm ---------- */
j = 1;
v1 = f1 * ( Elt[j].ArrayDYIN[i1p] - Elt[j].ArrayDYIN[i1] ) + Elt[j].ArrayDYIN[i1];
v2 = f2 * ( Elt[j].ArrayDYIN[i2p] - Elt[j].ArrayDYIN[i2] ) + Elt[j].ArrayDYIN[i2];
MassFracDYIN[j] = v2-v1;
/* ---------------------------- */
/* --------- Metals ---------- */
/* ---------------------------- */
j = 2;
l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ;
l2 = ( log10(m2) - Elt[j].MminDYIN) / Elt[j].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<0) i1=0;
if (i2<0) i2=0;
for (j=2;j<Cp->nelts+2;j++)
{
v1 = f1 * ( Elt[j].ArrayDYIN[i1p] - Elt[j].ArrayDYIN[i1] ) + Elt[j].ArrayDYIN[i1];
v2 = f2 * ( Elt[j].ArrayDYIN[i2p] - Elt[j].ArrayDYIN[i2] ) + Elt[j].ArrayDYIN[i2];
MassFracDYIN[j] = v2-v1;
}
}
void DYIN_mass_ejection_Z(double m1,double m2, double Z)
{
/*
Compute the mass fraction and yields of dying stars with masses between m1 and m2.
Store the result in the global variable`` MassFracDYIN``::
MassFracDYIN[0] = total gas
MassFracDYIN[1] = helium core (i.e. alpha(m))
MassFracDYIN[i] = frac mass elt i.
*/
double l1,l2;
int i1,i2,i1p,i2p,k;
int j1,j1p;
double f1,f2,g1;
double v1,v2,v11,v12;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
/* this was not in Poirier... */
m1 = dmax(m1,Cp->DYIN_Mmin);
m2 = dmin(m2,Cp->DYIN_Mmax);
if ( m2<=m1 )
{
for (k=0;k<Cp->nelts+2;k++)
MassFracDYIN[k] = 0;
return;
}
/* find i1,i1p, index in m1 and m2 */
k = 0;
l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ;
l2 = ( log10(m2) - Elt[k].MminDYIN) / Elt[k].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<0) i1=0;
if (i2<0) i2=0;
/* find j1,j1p, index in Z */
l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ;
if (l1 < 0.0) l1 = 0.0;
j1 = (int)l1;
j1p = j1 + 1;
g1 = l1 - j1;
/* check (yr) */
if (j1<0) j1=0;
/* --------- TOTAL GAS ---------- */
k = 0;
/* mass 1 */
v11 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1] - Elt[k].ArrayZDYIN[i1][j1] ) + Elt[k].ArrayZDYIN[i1][j1];
v12 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1p] - Elt[k].ArrayZDYIN[i1][j1p] ) + Elt[k].ArrayZDYIN[i1][j1p];
v1 = g1 * ( v12 - v11 ) + v11;
/* mass 2 */
v11 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1] - Elt[k].ArrayZDYIN[i2][j1] ) + Elt[k].ArrayZDYIN[i2][j1];
v12 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1p] - Elt[k].ArrayZDYIN[i2][j1p] ) + Elt[k].ArrayZDYIN[i2][j1p];
v2 = g1 * ( v12 - v11 ) + v11;
MassFracDYIN[k] = v2-v1;
/* --------- He core therm ---------- */
k = 1;
/* mass 1 */
v11 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1] - Elt[k].ArrayZDYIN[i1][j1] ) + Elt[k].ArrayZDYIN[i1][j1];
v12 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1p] - Elt[k].ArrayZDYIN[i1][j1p] ) + Elt[k].ArrayZDYIN[i1][j1p];
v1 = g1 * ( v12 - v11 ) + v11;
/* mass 2 */
v11 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1] - Elt[k].ArrayZDYIN[i2][j1] ) + Elt[k].ArrayZDYIN[i2][j1];
v12 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1p] - Elt[k].ArrayZDYIN[i2][j1p] ) + Elt[k].ArrayZDYIN[i2][j1p];
v2 = g1 * ( v12 - v11 ) + v11;
MassFracDYIN[k] = v2-v1;
/* ---------------------------- */
/* --------- Metals ---------- */
/* ---------------------------- */
k = 2;
l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ;
l2 = ( log10(m2) - Elt[k].MminDYIN) / Elt[k].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<0) i1=0;
if (i2<0) i2=0;
/* find j1,j1p, index in Z */
l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ;
if (l1 < 0.0) l1 = 0.0;
j1 = (int)l1;
j1p = j1 + 1;
g1 = l1 - j1;
/* check (yr) */
if (j1<0) j1=0;
for (k=2;k<Cp->nelts+2;k++)
{
/* mass 1 */
v11 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1] - Elt[k].ArrayZDYIN[i1][j1] ) + Elt[k].ArrayZDYIN[i1][j1];
v12 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1p] - Elt[k].ArrayZDYIN[i1][j1p] ) + Elt[k].ArrayZDYIN[i1][j1p];
v1 = g1 * ( v12 - v11 ) + v11;
/* mass 2 */
v11 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1] - Elt[k].ArrayZDYIN[i2][j1] ) + Elt[k].ArrayZDYIN[i2][j1];
v12 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1p] - Elt[k].ArrayZDYIN[i2][j1p] ) + Elt[k].ArrayZDYIN[i2][j1p];
v2 = g1 * ( v12 - v11 ) + v11;
MassFracDYIN[k] = v2-v1;
}
}
void DYIN_single_mass_ejection(double m1)
{
/*
Compute the mass fraction and yields of a dying stars of masse m1.
Store the result in the global variable ``SingleMassFracDYIN``::
SingleMassFracDYIN[0] = total gas
SingleMassFracDYIN[1] = helium core (i.e. alpha(m))
SingleMassFracDYIN[i] = frac mass elt i.
*/
double l1;
int i1,i1p,j;
double f1;
double v1;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
/* this was not in Poirier... */
if ( (m1<=Cp->DYIN_Mmin) || (m1>=Cp->DYIN_Mmax) )
{
for (j=0;j<Cp->nelts+2;j++)
SingleMassFracDYIN[j] = 0;
return;
}
j = 0;
l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
i1 = (int)l1;
i1p = i1 + 1;
f1 = l1 - i1;
/* check (yr) */
if (i1<0) i1=0;
/* --------- TOTAL GAS ---------- */
j = 0;
v1 = f1 * ( Elt[j].MetalDYIN[i1p] - Elt[j].MetalDYIN[i1] ) + Elt[j].MetalDYIN[i1];
SingleMassFracDYIN[j] = v1;
/* --------- He core therm ---------- */
j = 1;
v1 = f1 * ( Elt[j].MetalDYIN[i1p] - Elt[j].MetalDYIN[i1] ) + Elt[j].MetalDYIN[i1];
SingleMassFracDYIN[j] = v1;
/* ---------------------------- */
/* --------- Metals ---------- */
/* ---------------------------- */
j = 2;
l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
i1 = (int)l1;
i1p = i1 + 1;
f1 = l1 - i1;
/* check (yr) */
if (i1<0) i1=0;
for (j=2;j<Cp->nelts+2;j++)
{
v1 = f1 * ( Elt[j].MetalDYIN[i1p] - Elt[j].MetalDYIN[i1] ) + Elt[j].MetalDYIN[i1];
SingleMassFracDYIN[j] = v1;
}
}
void DYIN_single_mass_ejection_Z(double m1,double Z)
{
/*
Compute the mass fraction and yields of a dying stars of masse m1.
Store the result in the global variable ``SingleMassFracDYIN``::
SingleMassFracDYIN[0] = total gas
SingleMassFracDYIN[1] = helium core (i.e. alpha(m))
SingleMassFracDYIN[i] = frac mass elt i.
*/
double l1;
int i1,i1p,k;
int j1,j1p;
double f1,g1;
double v1,v11,v12;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
/* this was not in Poirier... */
if ( (m1<=Cp->DYIN_Mmin) || (m1>=Cp->DYIN_Mmax) )
{
for (k=0;k<Cp->nelts+2;k++)
SingleMassFracDYIN[k] = 0;
return;
}
k = 0;
/* find i1,i1p, index in m */
l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
i1 = (int)l1;
i1p = i1 + 1;
f1 = l1 - i1;
/* check (yr) */
if (i1<0) i1=0;
/* find j1,j1p, index in Z */
l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ;
if (l1 < 0.0) l1 = 0.0;
j1 = (int)l1;
j1p = j1 + 1;
g1 = l1 - j1;
/* check (yr) */
if (j1<0) j1=0;
/* --------- TOTAL GAS ---------- */
k = 0;
v11 = f1 * ( Elt[k].MetalZDYIN[i1p][j1] - Elt[k].MetalZDYIN[i1][j1] ) + Elt[k].MetalZDYIN[i1][j1];
v12 = f1 * ( Elt[k].MetalZDYIN[i1p][j1p] - Elt[k].MetalZDYIN[i1][j1p] ) + Elt[k].MetalZDYIN[i1p][j1];
v1 = g1 * ( v12 - v11 ) + v11;
SingleMassFracDYIN[k] = v1;
/* --------- He core therm ---------- */
k = 1;
v11 = f1 * ( Elt[k].MetalZDYIN[i1p][j1] - Elt[k].MetalZDYIN[i1][j1] ) + Elt[k].MetalZDYIN[i1][j1];
v12 = f1 * ( Elt[k].MetalZDYIN[i1p][j1p] - Elt[k].MetalZDYIN[i1][j1p] ) + Elt[k].MetalZDYIN[i1p][j1];
v1 = g1 * ( v12 - v11 ) + v11;
SingleMassFracDYIN[k] = v1;
/* ---------------------------- */
/* --------- Metals ---------- */
/* ---------------------------- */
k = 2;
/* find i1,i1p, index in m */
l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ;
if (l1 < 0.0) l1 = 0.0;
i1 = (int)l1;
i1p = i1 + 1;
f1 = l1 - i1;
/* check (yr) */
if (i1<0) i1=0;
/* find j1,j1p, index in Z */
l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ;
if (l1 < 0.0) l1 = 0.0;
j1 = (int)l1;
j1p = j1 + 1;
g1 = l1 - j1;
/* check (yr) */
if (j1<0) j1=0;
for (k=2;k<Cp->nelts+2;k++)
{
v11 = f1 * ( Elt[k].MetalZDYIN[i1p][j1] - Elt[k].MetalZDYIN[i1][j1] ) + Elt[k].MetalZDYIN[i1][j1];
v12 = f1 * ( Elt[k].MetalZDYIN[i1p][j1p] - Elt[k].MetalZDYIN[i1][j1p] ) + Elt[k].MetalZDYIN[i1p][j1];
v1 = g1 * ( v12 - v11 ) + v11;
SingleMassFracDYIN[k] = v1;
}
}
void SNII_mass_ejection(double m1,double m2)
{
/*
.. warning:: here, we we do not limit the computation to SNII !!!
Compute the mass fraction and yields of SNII stars with masses between m1 and m2.
Store the result in the global variable ``MassFracSNII``::
MassFracSNII[0] = total gas
MassFracSNII[1] = 1-helium core (i.e. non processed elts)
MassFracSNII[i] = frac mass elt i.
*/
double l1,l2;
int i1,i2,i1p,i2p,j;
double f1,f2;
double v1,v2;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
m2 = m2*All.CMUtoMsol;
/* this was not in Poirier... */
m1 = dmax(m1,Cp->SNII_Mmin);
m2 = dmin(m2,Cp->SNII_Mmax);
if ( m2<=m1 )
{
for (j=0;j<Cp->nelts+2;j++)
MassFracSNII[j] = 0;
return;
}
j = 0;
l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ;
l2 = ( log10(m2) - Elt[j].MminSNII) / Elt[j].StepSNII ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<0) i1=0;
if (i2<0) i2=0;
/* --------- TOTAL GAS ---------- */
j = 0;
v1 = f1 * ( Elt[j].ArraySNII[i1p] - Elt[j].ArraySNII[i1] ) + Elt[j].ArraySNII[i1];
v2 = f2 * ( Elt[j].ArraySNII[i2p] - Elt[j].ArraySNII[i2] ) + Elt[j].ArraySNII[i2];
MassFracSNII[j] = v2-v1;
/* --------- He core therm ---------- */
j = 1;
v1 = f1 * ( Elt[j].ArraySNII[i1p] - Elt[j].ArraySNII[i1] ) + Elt[j].ArraySNII[i1];
v2 = f2 * ( Elt[j].ArraySNII[i2p] - Elt[j].ArraySNII[i2] ) + Elt[j].ArraySNII[i2];
MassFracSNII[j] = v2-v1;
/* ---------------------------- */
/* --------- Metals ---------- */
/* ---------------------------- */
j = 2;
l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ;
l2 = ( log10(m2) - Elt[j].MminSNII) / Elt[j].StepSNII ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<0) i1=0;
if (i2<0) i2=0;
for (j=2;j<Cp->nelts+2;j++)
{
v1 = f1 * ( Elt[j].ArraySNII[i1p] - Elt[j].ArraySNII[i1] ) + Elt[j].ArraySNII[i1];
v2 = f2 * ( Elt[j].ArraySNII[i2p] - Elt[j].ArraySNII[i2] ) + Elt[j].ArraySNII[i2];
MassFracSNII[j] = v2-v1;
}
}
// void SNII_mass_ejection_Z(double m1,double m2, double Z)
// {
//
// /*
// Compute the mass fraction and yields of SNII with masses between m1 and m2.
// Store the result in the global variable`` MassFracSNII``::
//
// MassFracSNII[0] = total gas
// MassFracSNII[1] = helium core (i.e. alpha(m))
// MassFracSNII[i] = frac mass elt i.
//
// */
//
// double l1,l2;
// int i1,i2,i1p,i2p,k;
// int j1,j1p;
// double f1,f2,g1;
// double v1,v2,v11,v12;
//
// /* convert m in msol */
// m1 = m1*All.CMUtoMsol;
// m2 = m2*All.CMUtoMsol;
//
// /* this was not in Poirier... */
// m1 = dmax(m1,Cp->SNII_Mmin);
// m2 = dmin(m2,Cp->SNII_Mmax);
//
//
// if ( m2<=m1 )
// {
// for (k=0;k<Cp->nelts+2;k++)
// MassFracSNII[k] = 0;
// return;
// }
//
//
//
//
// /* find i1,i1p, index in m1 and m2 */
//
// k = 0;
// l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ;
// l2 = ( log10(m2) - Elt[k].MminSNII) / Elt[k].StepSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
// if (l2 < 0.0) l2 = 0.0;
//
// i1 = (int)l1;
// i2 = (int)l2;
//
// i1p = i1 + 1;
// i2p = i2 + 1;
//
// f1 = l1 - i1;
// f2 = l2 - i2;
//
// /* check (yr) */
// if (i1<0) i1=0;
// if (i2<0) i2=0;
//
//
// /* find j1,j1p, index in Z */
//
// l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
//
// j1 = (int)l1;
// j1p = j1 + 1;
// g1 = l1 - j1;
//
// /* check (yr) */
// if (j1<0) j1=0;
//
//
//
// /* --------- TOTAL GAS ---------- */
// k = 0;
//
// /* mass 1 */
// v11 = f1 * ( Elt[k].ArrayZSNII[i1p][j1] - Elt[k].ArrayZSNII[i1][j1] ) + Elt[k].ArrayZSNII[i1][j1];
// v12 = f1 * ( Elt[k].ArrayZSNII[i1p][j1p] - Elt[k].ArrayZSNII[i1][j1p] ) + Elt[k].ArrayZSNII[i1p][j1];
// v1 = g1 * ( v12 - v11 ) + v11;
//
// /* mass 2 */
// v11 = f2 * ( Elt[k].ArrayZSNII[i2p][j1] - Elt[k].ArrayZSNII[i2][j1] ) + Elt[k].ArrayZSNII[i2][j1];
// v12 = f2 * ( Elt[k].ArrayZSNII[i2p][j1p] - Elt[k].ArrayZSNII[i2][j1p] ) + Elt[k].ArrayZSNII[i2p][j1];
// v2 = g1 * ( v12 - v11 ) + v11;
//
// MassFracSNII[k] = v2-v1;
//
// /* --------- He core therm ---------- */
// k = 1;
//
// /* mass 1 */
// v11 = f1 * ( Elt[k].ArrayZSNII[i1p][j1] - Elt[k].ArrayZSNII[i1][j1] ) + Elt[k].ArrayZSNII[i1][j1];
// v12 = f1 * ( Elt[k].ArrayZSNII[i1p][j1p] - Elt[k].ArrayZSNII[i1][j1p] ) + Elt[k].ArrayZSNII[i1p][j1];
// v1 = g1 * ( v12 - v11 ) + v11;
//
// /* mass 2 */
// v11 = f2 * ( Elt[k].ArrayZSNII[i2p][j1] - Elt[k].ArrayZSNII[i2][j1] ) + Elt[k].ArrayZSNII[i2][j1];
// v12 = f2 * ( Elt[k].ArrayZSNII[i2p][j1p] - Elt[k].ArrayZSNII[i2][j1p] ) + Elt[k].ArrayZSNII[i2p][j1];
// v2 = g1 * ( v12 - v11 ) + v11;
//
// MassFracSNII[k] = v2-v1;
//
// /* ---------------------------- */
// /* --------- Metals ---------- */
// /* ---------------------------- */
//
// k = 2;
//
// l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ;
// l2 = ( log10(m2) - Elt[k].MminSNII) / Elt[k].StepSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
// if (l2 < 0.0) l2 = 0.0;
//
// i1 = (int)l1;
// i2 = (int)l2;
//
// i1p = i1 + 1;
// i2p = i2 + 1;
//
// f1 = l1 - i1;
// f2 = l2 - i2;
//
// /* check (yr) */
// if (i1<0) i1=0;
// if (i2<0) i2=0;
//
// /* find j1,j1p, index in Z */
//
// l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
//
// j1 = (int)l1;
// j1p = j1 + 1;
// g1 = l1 - j1;
//
// /* check (yr) */
// if (j1<0) j1=0;
//
// for (k=2;k<Cp->nelts+2;k++)
// {
//
// /* mass 1 */
// v11 = f1 * ( Elt[k].ArrayZSNII[i1p][j1] - Elt[k].ArrayZSNII[i1][j1] ) + Elt[k].ArrayZSNII[i1][j1];
// v12 = f1 * ( Elt[k].ArrayZSNII[i1p][j1p] - Elt[k].ArrayZSNII[i1][j1p] ) + Elt[k].ArrayZSNII[i1p][j1];
// v1 = g1 * ( v12 - v11 ) + v11;
//
// /* mass 2 */
// v11 = f2 * ( Elt[k].ArrayZSNII[i2p][j1] - Elt[k].ArrayZSNII[i2][j1] ) + Elt[k].ArrayZSNII[i2][j1];
// v12 = f2 * ( Elt[k].ArrayZSNII[i2p][j1p] - Elt[k].ArrayZSNII[i2][j1p] ) + Elt[k].ArrayZSNII[i2p][j1];
// v2 = g1 * ( v12 - v11 ) + v11;
//
// MassFracSNII[k] = v2-v1;
//
// }
//
// }
void SNII_single_mass_ejection(double m1)
{
/*
.. warning:: here, we we do not limit the computation to SNII !!!
Compute the mass fraction and yields of a SNII stars of masse m1.
Store the result in the global variable ``SingleMassFracSNII``::
SingleMassFracSNII[0] = total gas
SingleMassFracSNII[1] = 1-helium core (i.e. non processed elts)
SingleMassFracSNII[i] = frac mass elt i.
*/
double l1;
int i1,i1p,j;
double f1;
double v1;
/* convert m in msol */
m1 = m1*All.CMUtoMsol;
/* this was not in Poirier... */
if ( (m1<=Cp->SNII_Mmin) || (m1>=Cp->SNII_Mmax) )
{
for (j=0;j<Cp->nelts+2;j++)
SingleMassFracSNII[j] = 0;
return;
}
j = 0;
l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ;
if (l1 < 0.0) l1 = 0.0;
i1 = (int)l1;
i1p = i1 + 1;
f1 = l1 - i1;
/* check (yr) */
if (i1<0) i1=0;
/* --------- TOTAL GAS ---------- */
j = 0;
v1 = f1 * ( Elt[j].MetalSNII[i1p] - Elt[j].MetalSNII[i1] ) + Elt[j].MetalSNII[i1];
SingleMassFracSNII[j] = v1;
/* --------- He core therm ---------- */
j = 1;
v1 = f1 * ( Elt[j].MetalSNII[i1p] - Elt[j].MetalSNII[i1] ) + Elt[j].MetalSNII[i1];
SingleMassFracSNII[j] = v1;
/* ---------------------------- */
/* --------- Metals ---------- */
/* ---------------------------- */
j = 2;
l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ;
if (l1 < 0.0) l1 = 0.0;
i1 = (int)l1;
i1p = i1 + 1;
f1 = l1 - i1;
/* check (yr) */
if (i1<0) i1=0;
for (j=2;j<Cp->nelts+2;j++)
{
v1 = f1 * ( Elt[j].MetalSNII[i1p] - Elt[j].MetalSNII[i1] ) + Elt[j].MetalSNII[i1];
SingleMassFracSNII[j] = v1;
}
}
// void SNII_single_mass_ejection_Z(double m1,double Z)
// {
//
//
// /*
// Compute the mass fraction and yields of a SNII stars of masse m1.
// Store the result in the global variable ``SingleMassFracSNII``::
//
// SingleMassFracSNII[0] = total gas
// SingleMassFracSNII[1] = helium core (i.e. alpha(m))
// SingleMassFracSNII[i] = frac mass elt i.
// */
//
// double l1;
// int i1,i1p,k;
// int j1,j1p;
// double f1,g1;
// double v1,v11,v12;
//
// /* convert m in msol */
// m1 = m1*All.CMUtoMsol;
//
// /* this was not in Poirier... */
// if ( (m1<=Cp->SNII_Mmin) || (m1>=Cp->SNII_Mmax) )
// {
// for (k=0;k<Cp->nelts+2;k++)
// SingleMassFracSNII[k] = 0;
// return;
// }
//
//
// k = 0;
//
//
// /* find i1,i1p, index in m */
//
// l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
//
// i1 = (int)l1;
// i1p = i1 + 1;
// f1 = l1 - i1;
//
// /* check (yr) */
// if (i1<0) i1=0;
//
//
// /* find j1,j1p, index in Z */
//
// l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
//
// j1 = (int)l1;
// j1p = j1 + 1;
// g1 = l1 - j1;
//
// /* check (yr) */
// if (j1<0) j1=0;
//
//
//
// /* --------- TOTAL GAS ---------- */
// k = 0;
// v11 = f1 * ( Elt[k].MetalZSNII[i1p][j1] - Elt[k].MetalZSNII[i1][j1] ) + Elt[k].MetalZSNII[i1][j1];
// v12 = f1 * ( Elt[k].MetalZSNII[i1p][j1p] - Elt[k].MetalZSNII[i1][j1p] ) + Elt[k].MetalZSNII[i1p][j1];
// v1 = g1 * ( v12 - v11 ) + v11;
//
// SingleMassFracSNII[k] = v1;
//
// /* --------- He core therm ---------- */
// k = 1;
// v11 = f1 * ( Elt[k].MetalZSNII[i1p][j1] - Elt[k].MetalZSNII[i1][j1] ) + Elt[k].MetalZSNII[i1][j1];
// v12 = f1 * ( Elt[k].MetalZSNII[i1p][j1p] - Elt[k].MetalZSNII[i1][j1p] ) + Elt[k].MetalZSNII[i1p][j1];
// v1 = g1 * ( v12 - v11 ) + v11;
//
// SingleMassFracSNII[k] = v1;
//
// /* ---------------------------- */
// /* --------- Metals ---------- */
// /* ---------------------------- */
//
// k = 2;
//
// /* find i1,i1p, index in m */
//
// l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
//
// i1 = (int)l1;
// i1p = i1 + 1;
// f1 = l1 - i1;
//
// /* check (yr) */
// if (i1<0) i1=0;
//
//
// /* find j1,j1p, index in Z */
//
// l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ;
//
// if (l1 < 0.0) l1 = 0.0;
//
// j1 = (int)l1;
// j1p = j1 + 1;
// g1 = l1 - j1;
//
// /* check (yr) */
// if (j1<0) j1=0;
//
// for (k=2;k<Cp->nelts+2;k++)
// {
//
// v11 = f1 * ( Elt[k].MetalZSNII[i1p][j1] - Elt[k].MetalZSNII[i1][j1] ) + Elt[k].MetalZSNII[i1][j1];
// v12 = f1 * ( Elt[k].MetalZSNII[i1p][j1p] - Elt[k].MetalZSNII[i1][j1p] ) + Elt[k].MetalZSNII[i1p][j1];
// v1 = g1 * ( v12 - v11 ) + v11;
//
// SingleMassFracSNII[k] = v1;
// }
//
// }
void SNIa_mass_ejection(double m1,double m2)
{
/*
Compute the total mass and element mass per mass unit of SNIa stars with masses between m1 and m2.
Store the result in the global variable ``MassFracSNIa``::
MassFracSNIa[0] = total gas
MassFracSNIa[1] = unused
MassFracSNIa[i] = frac mass elt i.
*/
int j;
double NSNIa;
/* number of SNIa per mass unit between time and time+dt */
NSNIa = SNIa_rate(m1,m2);
/* ejected mass in gas per mass unit */
MassFracSNIa[0] = Cp->Mco * All.MsoltoCMU * NSNIa;
/* ejected elements in gas per mass unit */
for (j=2;j<Cp->nelts+2;j++)
MassFracSNIa[j] = NSNIa* Elt[j].MSNIa * All.MsoltoCMU;
/* unused */
MassFracSNIa[1]=-1;
}
void SNIa_single_mass_ejection(double m1)
{
/*
Compute the total mass mass of element of a SNIa stars of masse m1.
Store the result in the global variable ``SingleMassFracSNIa``::
SingleMassFracSNIa[0] = total gas
SingleMassFracSNIa[1] = unused
SingleMassFracSNIa[i] = frac mass elt i.
*/
int j;
if (Cp->SNIa_Mpl<0) // no SNIa
{
/* total ejected gas mass */
SingleMassFracSNIa[0] = 0;
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
SingleMassFracSNIa[j] = 0;
/* unused */
SingleMassFracSNIa[1] = 0;
}
else
{
/* total ejected gas mass */
SingleMassFracSNIa[0] = Cp->Mco * All.MsoltoCMU;
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
SingleMassFracSNIa[j] = Elt[j].MSNIa * All.MsoltoCMU;
/* unused */
SingleMassFracSNIa[1] = -1;
}
}
void Total_mass_ejection(double m1,double m2,double M0,double *z)
{
/*
Sum the contribution in mass and yields of all stars in the mass range m1,m2.
Store the result in the global variable EjectedMass::
EjectedMass[0] = total gas
EjectedMass[1] = UNUSED
EjectedMass[i+2] = frac mass elt i.
FOR THE MOMENT::
- contrib of SNII (= all stars)
- contrib of SNIa
EjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa
EjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts)
*/
int j;
/* compute SNII mass ejection -> MassFracSNII */
SNII_mass_ejection(m1,m2);
/* compute SNIa mass ejection -> MassFracSNIa */ /* not really a mass fraction */
SNIa_mass_ejection(m1,m2);
/* compute DYIN mass ejection -> MassFracDYIN */ /* not really a mass fraction */
#ifdef CHIMIE_AGB
DYIN_mass_ejection(m1,m2);
#endif
/* total ejected gas mass */
#ifdef CHIMIE_AGB
EjectedMass[0] = M0 * ( MassFracDYIN[0] + MassFracSNII[0] + MassFracSNIa[0] );
#else
EjectedMass[0] = M0 * ( MassFracSNII[0] + MassFracSNIa[0] );
#endif
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
#ifdef CHIMIE_AGB
EjectedMass[j] = M0*( MassFracDYIN[j] +z[j-2]*MassFracDYIN[1] + MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] );
#else
EjectedMass[j] = M0*( MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] );
#endif
/* not used */
EjectedMass[1] = -1;
}
void Total_mass_ejection_Z(double m1,double m2,double M0,double *z)
{
/*
Sum the contribution in mass and yields of all stars in the mass range m1,m2.
Store the result in the global variable EjectedMass::
EjectedMass[0] = total gas
EjectedMass[1] = UNUSED
EjectedMass[i+2] = frac mass elt i.
FOR THE MOMENT::
- contrib of SNII (= all stars)
- contrib of SNIa
EjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa
EjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts)
*/
int j;
float Z; /* metalicity */
Z = z[METALS];
/* compute SNII mass ejection -> MassFracSNII */
SNII_mass_ejection(m1,m2);
/* compute SNIa mass ejection -> MassFracSNIa */ /* not really a mass fraction */
SNIa_mass_ejection(m1,m2);
/* compute DYIN mass ejection -> MassFracDYIN */ /* not really a mass fraction */
#ifdef CHIMIE_AGB
DYIN_mass_ejection_Z(m1,m2,Z);
#endif
/* total ejected gas mass */
#ifdef CHIMIE_AGB
EjectedMass[0] = M0 * ( MassFracDYIN[0] + MassFracSNII[0] + MassFracSNIa[0] );
#else
EjectedMass[0] = M0 * ( MassFracSNII[0] + MassFracSNIa[0] );
#endif
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
#ifdef CHIMIE_AGB
EjectedMass[j] = M0*( MassFracDYIN[j] +z[j-2]*MassFracDYIN[1] + MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] );
#else
EjectedMass[j] = M0*( MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] );
#endif
/* not used */
EjectedMass[1] = -1;
}
void DYIN_Total_single_mass_ejection(double m1,double *z)
{
/*
Mass and element ejected by a single dying stars of mass m1.
This takes into account processed and non processed gas
The results are stored in::
SingleEjectedMass[0] = gas mass
SingleEjectedMass[1] = unsued
SingleEjectedMass[i+2] = frac mass elt i
*/
int j;
float M0;
M0 = m1;
/* compute dying stars mass ejection -> SingleMassFracDYIN */
DYIN_single_mass_ejection(m1);
/* total ejected gas mass */
SingleEjectedMass[0] = M0 * SingleMassFracDYIN[0];
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
SingleEjectedMass[j] = M0*(SingleMassFracDYIN[j] +z[j-2]*SingleMassFracDYIN[1]);
/* not used */
SingleEjectedMass[1] = -1;
}
void SNII_Total_single_mass_ejection(double m1,double *z)
{
/*
Mass and element ejected by a single SNII of mass m1.
This takes into account processed and non processed gas
The results are stored in::
SingleEjectedMass[0] = gas mass
SingleEjectedMass[1] = unsued
SingleEjectedMass[i+2] = frac mass elt i
*/
int j;
float M0;
M0 = m1;
/* compute SNII mass ejection -> SingleMassFracSNII */
SNII_single_mass_ejection(m1);
/* total ejected gas mass */
SingleEjectedMass[0] = M0 * SingleMassFracSNII[0];
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
SingleEjectedMass[j] = M0*(SingleMassFracSNII[j] +z[j-2]*SingleMassFracSNII[1]);
/* not used */
SingleEjectedMass[1] = -1;
}
void SNIa_Total_single_mass_ejection(double m1, double *z)
{
int j;
/*
Mass and element ejected by a single SNIa of mass m1.
The results are stored in::
SingleEjectedMass[0] = gas mass
SingleEjectedMass[1] = unsued
SingleEjectedMass[i+2] = frac mass elt i
*/
/* compute SNIa mass ejection -> SingleMassFracSNIa */
SNIa_single_mass_ejection(m1);
/* total ejected gas mass */
SingleEjectedMass[0] = SingleMassFracSNIa[0];
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
SingleEjectedMass[j] = SingleMassFracSNIa[j];
}
void Total_single_mass_ejection(double m1,double *z,double NSNII,double NSNIa,double NDYIN)
{
/*
Sum the contribution in mass and yields of one star for mass m1.
Store the result in the global variable EjectedMass::
SingleEjectedMass[0] = total gas
SingleEjectedMass[1] = UNUSED
SingleEjectedMass[i+2] = frac mass elt i.
FOR THE MOMENT::
- contrib of SNII (= all stars)
- contrib of SNIa
SingleEjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa
SingleEjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts)
*/
int j;
float M0;
M0 = m1;
/* compute SNII mass ejection -> SingleMassFracSNII */
SNII_single_mass_ejection(m1);
/* compute SNII mass ejection -> SingleMassFracSNIa */
SNIa_single_mass_ejection(m1);
/* compute DYIN mass ejection -> SingleMassFracDYIN */
#ifdef CHIMIE_AGB
DYIN_single_mass_ejection(m1);
#endif
/* total ejected gas mass */
#ifdef CHIMIE_AGB
SingleEjectedMass[0] = M0 * ( SingleMassFracDYIN[0]*NDYIN + SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa;
#else
SingleEjectedMass[0] = M0 * ( SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa;
#endif
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
#ifdef CHIMIE_AGB
SingleEjectedMass[j] = M0*( SingleMassFracDYIN[j]*NDYIN +z[j-2]*SingleMassFracDYIN[1]*NDYIN + SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa;
#else
SingleEjectedMass[j] = M0*( SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa;
#endif
/* not used */
SingleEjectedMass[1] = -1;
}
void Total_single_mass_ejection_Z(double m1,double *z,double NSNII,double NSNIa,double NDYIN)
{
/*
Sum the contribution in mass and yields of one star for mass m1.
Store the result in the global variable EjectedMass::
SingleEjectedMass[0] = total gas
SingleEjectedMass[1] = UNUSED
SingleEjectedMass[i+2] = frac mass elt i.
FOR THE MOMENT::
- contrib of SNII (= all stars)
- contrib of SNIa
SingleEjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa
SingleEjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts)
*/
int j;
float M0;
float Z; /* metalicity */
M0 = m1;
Z = z[METALS];
//z[METALS]=0;
/* compute SNII mass ejection -> SingleMassFracSNII */
SNII_single_mass_ejection(m1);
/* compute SNII mass ejection -> SingleMassFracSNIa */
SNIa_single_mass_ejection(m1);
/* compute DYIN mass ejection -> SingleMassFracDYIN */
#ifdef CHIMIE_AGB
DYIN_single_mass_ejection_Z(m1,Z);
#endif
/* total ejected gas mass */
#ifdef CHIMIE_AGB
SingleEjectedMass[0] = M0 * ( SingleMassFracDYIN[0]*NDYIN + SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa;
#else
SingleEjectedMass[0] = M0 * ( SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa;
#endif
/* ejected mass per element */
for (j=2;j<Cp->nelts+2;j++)
#ifdef CHIMIE_AGB
SingleEjectedMass[j] = M0*( SingleMassFracDYIN[j]*NDYIN +z[j-2]*SingleMassFracDYIN[1]*NDYIN + SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa;
#else
SingleEjectedMass[j] = M0*( SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa;
#endif
/* not used */
SingleEjectedMass[1] = -1;
}
/****************************************************************************************/
/*
/*
/*
/* GADGET ONLY PART
/*
/*
/*
/****************************************************************************************/
static double hubble_a, atime, hubble_a2, fac_mu, fac_vsic_fix, a3inv, fac_egy;
#ifdef FEEDBACK
static double fac_pow;
#endif
#ifdef PERIODIC
static double boxSize, boxHalf;
#ifdef LONG_X
static double boxSize_X, boxHalf_X;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static double boxSize_Y, boxHalf_Y;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static double boxSize_Z, boxHalf_Z;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif
#if defined(CHIMIE_THERMAL_FEEDBACK) && defined(CHIMIE_COMPUTE_THERMAL_FEEDBACK_ENERGY)
void chimie_compute_energy_int(int mode)
{
int i;
double DeltaEgyInt;
double Tot_DeltaEgyInt;
DeltaEgyInt = 0;
Tot_DeltaEgyInt = 0;
if (mode==1)
{
LocalSysState.EnergyInt1 = 0;
LocalSysState.EnergyInt2 = 0;
}
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (mode==1)
#ifdef DENSITY_INDEPENDENT_SPH
LocalSysState.EnergyInt1 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1);
#else
LocalSysState.EnergyInt1 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].Density*a3inv, GAMMA_MINUS1);
#endif
else
#ifdef DENSITY_INDEPENDENT_SPH
LocalSysState.EnergyInt2 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1);
#else
LocalSysState.EnergyInt2 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].Density*a3inv, GAMMA_MINUS1);
#endif
}
}
if (mode==2)
{
DeltaEgyInt = LocalSysState.EnergyInt2 - LocalSysState.EnergyInt1;
MPI_Reduce(&DeltaEgyInt, &Tot_DeltaEgyInt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
LocalSysState.EnergyThermalFeedback -= DeltaEgyInt;
}
}
#endif
#if defined(CHIMIE_KINETIC_FEEDBACK) && defined(CHIMIE_COMPUTE_KINETIC_FEEDBACK_ENERGY)
void chimie_compute_energy_kin(int mode)
{
int i;
double DeltaEgyKin;
double Tot_DeltaEgyKin;
DeltaEgyKin = 0;
Tot_DeltaEgyKin = 0;
if (mode==1)
{
LocalSysState.EnergyKin1 = 0;
LocalSysState.EnergyKin2 = 0;
}
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (mode==1)
LocalSysState.EnergyKin1 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
else
LocalSysState.EnergyKin2 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
}
}
if (mode==2)
{
DeltaEgyKin = LocalSysState.EnergyKin2 - LocalSysState.EnergyKin1;
MPI_Reduce(&DeltaEgyKin, &Tot_DeltaEgyKin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
LocalSysState.EnergyKineticFeedback -= DeltaEgyKin;
}
}
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
void chimie_apply_thermal_feedback(void)
{
int i;
double EgySpec,NewEgySpec,DeltaEntropy;
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (SphP[i].DeltaEgySpec > 0)
{
//printf("(%d) Step=%d i=%08d particle receive feedback\n",ThisTask,All.NumCurrentTiStep,i);
/* spec energy at current step (allways compute energy budget based on predicted entropy) */
#ifdef DENSITY_INDEPENDENT_SPH
EgySpec = SphP[i].EntropyPred / GAMMA_MINUS1 * pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1);
#else
EgySpec = SphP[i].EntropyPred / GAMMA_MINUS1 * pow(SphP[i].Density *a3inv, GAMMA_MINUS1);
#endif
/* new egyspec */
NewEgySpec = EgySpec + SphP[i].DeltaEgySpec;
LocalSysState.EnergyThermalFeedback -= SphP[i].DeltaEgySpec*P[i].Mass;
/* new entropy */
#ifdef DENSITY_INDEPENDENT_SPH
DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1) - SphP[i].EntropyPred;
#else
DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[i].Density *a3inv, GAMMA_MINUS1) - SphP[i].EntropyPred;
#endif
SphP[i].EntropyPred += DeltaEntropy;
SphP[i].Entropy += DeltaEntropy;
#ifdef DENSITY_INDEPENDENT_SPH
SphP[i].EntVarPred = pow(SphP[i].EntropyPred, 1/GAMMA);
#endif
/* set the adiabatic period for SNIa */
if (SphP[i].NumberOfSNIa>0)
SphP[i].SNIaThermalTime = All.Time;
/* set the adiabatic period for SNII */
if (SphP[i].NumberOfSNII>0)
SphP[i].SNIIThermalTime = All.Time;
/* reset variables */
SphP[i].DeltaEgySpec = 0;
SphP[i].NumberOfSNIa = 0;
SphP[i].NumberOfSNII = 0;
}
}
}
}
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
void chimie_apply_wind(void)
{
/* apply wind */
int i;
double e1,e2;
double phi,costh,sinth,vx,vy,vz;
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (SphP[i].WindFlag)
{
phi = get_ChimieKineticFeedback_random_number(P[i].ID)*PI*2.;
costh = 1.-2.*get_ChimieKineticFeedback_random_number(P[i].ID+1);
sinth = sqrt(1.-pow(costh,2));
vx = All.ChimieWindSpeed*sinth*cos(phi);
vy = All.ChimieWindSpeed*sinth*sin(phi);
vz = All.ChimieWindSpeed*costh;
e1 = 0.5*P[i].Mass * ( SphP[i].VelPred[0]*SphP[i].VelPred[0] + SphP[i].VelPred[1]*SphP[i].VelPred[1] + SphP[i].VelPred[2]*SphP[i].VelPred[2]);
P[i].Vel[0] += vx;
P[i].Vel[1] += vy;
P[i].Vel[2] += vz;
SphP[i].VelPred[0] += vx;
SphP[i].VelPred[1] += vy;
SphP[i].VelPred[2] += vz;
e2 = 0.5*P[i].Mass * ( SphP[i].VelPred[0]*SphP[i].VelPred[0] + SphP[i].VelPred[1]*SphP[i].VelPred[1] + SphP[i].VelPred[2]*SphP[i].VelPred[2]);
LocalSysState.EnergyKineticFeedback -= e2-e1;
SphP[i].WindFlag = 0;
}
}
}
}
#endif
/*! This function is the driver routine for the calculation of chemical evolution
*/
void chimie(void)
{
double t0, t1;
t0 = second(); /* measure the time for the full chimie computation */
#ifdef CHIMIE_TIMEBET
if((All.Time - All.ChimieTimeLastChimie) >= All.ChimieTimeBetChimie)
{
#endif
if (ThisTask==0)
printf("Start Chimie computation.\n");
if(All.ComovingIntegrationOn)
{
/* Factors for comoving integration of hydro */
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);
hubble_a2 = All.Time * All.Time * hubble_a;
fac_mu = pow(All.Time, 3 * (GAMMA - 1) / 2) / All.Time;
fac_egy = pow(All.Time, 3 * (GAMMA - 1));
fac_vsic_fix = hubble_a * pow(All.Time, 3 * GAMMA_MINUS1);
a3inv = 1 / (All.Time * All.Time * All.Time);
atime = All.Time;
#ifdef FEEDBACK
fac_pow = fac_egy*atime*atime;
#endif
}
else
{
hubble_a = hubble_a2 = atime = fac_mu = fac_vsic_fix = a3inv = fac_egy = 1.0;
#ifdef FEEDBACK
fac_pow = 1.0;
#endif
}
stars_density(); /* compute density */
#ifdef CHIMIE_ONE_SN_ONLY
if(All.ChimieOneSN==0) /* explode only if not one sn only*/
#endif
do_chimie(); /* chimie */
#ifdef CHIMIE_TIMEBET
All.ChimieTimeLastChimie += All.ChimieTimeBetChimie;
All.ChimieTi_begstep = All.Ti_Current;
if (ThisTask==0)
printf("Chimie: next Chimie=%g \n",All.ChimieTimeLastChimie);
#endif
if (ThisTask==0)
printf("Chimie computation done.\n");
t1 = second();
All.CPU_Chimie += timediff(t0, t1);
#ifdef CHIMIE_TIMEBET
}
#endif
}
/*! This function is the driver routine for the calculation of chemical evolution
*/
void do_chimie(void)
{
long long ntot, ntotleft;
int i, j, k, n, m, ngrp, maxfill, source, ndone;
int *nbuffer, *noffset, *nsend_local, *nsend, *numlist, *ndonelist;
int level, sendTask, recvTask, nexport, place;
double tstart, tend, sumt, sumcomm;
double timecomp = 0, timecommsumm = 0, timeimbalance = 0, sumimbalance;
int flag_chimie;
MPI_Status status;
int do_it;
int Ti0,Ti1,Ti2;
double t1,t2,t01,t02;
double tmin,tmax;
double minlivetime,maxlivetime;
double m1,m2,M0;
double NSNIa,NSNII,NDYIN;
double NSNIa_tot,NSNII_tot,NDYIN_tot,NSNIa_totlocal,NSNII_totlocal,NDYIN_totlocal;
double EgySN,EgySNlocal;
double EgySNThermal,EgySNKinetic;
int Nchim,Nchimlocal;
int Nwind,Nwindlocal;
int Nflag,Nflaglocal;
int Noldwind,Noldwindlocal;
double metals[NELEMENTS];
double FeH;
float MinRelMass=1e-3;
#ifdef DETAILED_CPU_OUTPUT_IN_CHIMIE
double *timecomplist;
double *timecommsummlist;
double *timeimbalancelist;
#endif
#ifdef PERIODIC
boxSize = All.BoxSize;
boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
boxHalf_X = boxHalf * LONG_X;
boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
boxHalf_Y = boxHalf * LONG_Y;
boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
boxHalf_Z = boxHalf * LONG_Z;
boxSize_Z = boxSize * LONG_Z;
#endif
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
double v1m,v2m;
#endif
/* `NumStUpdate' gives the number of particles on this processor that want a chimie computation */
for(n = 0, NumStUpdate = 0; n < N_gas+N_stars; n++)
{
#ifndef CHIMIE_TIMEBET
if(P[n].Ti_endstep == All.Ti_Current)
#else
/*nothing*/ /* as this is performed every ChimieTimeBetChimie, do it over all particles */
#endif
if(P[n].Type == ST)
{
m = P[n].StPIdx;
#ifdef FOF_SFR
if (StP[m].NStars != 0)
#endif
if ( (P[n].Mass/StP[m].InitialMass) > MinRelMass)
NumStUpdate++;
}
if(P[n].Type == 0)
SphP[n].dMass = 0.;
}
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumStUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
free(numlist);
noffset = malloc(sizeof(int) * NTask); /* offsets of bunches in common list */
nbuffer = malloc(sizeof(int) * NTask);
nsend_local = malloc(sizeof(int) * NTask);
nsend = malloc(sizeof(int) * NTask * NTask);
ndonelist = malloc(sizeof(int) * NTask);
i = 0; /* first gas particle, because stars may be hidden among gas particles */
ntotleft = ntot; /* particles left for all tasks together */
NSNIa_tot = 0;
NSNII_tot = 0;
NDYIN_tot = 0;
NSNIa_totlocal = 0;
NSNII_totlocal = 0;
NDYIN_totlocal = 0;
EgySN = 0;
EgySNlocal =0;
Nchimlocal = 0;
Nchim = 0;
Nwindlocal = 0;
Nwind = 0;
Noldwindlocal = 0;
Noldwind = 0;
Nflaglocal = 0;
Nflag = 0;
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
tstart = second();
for(nexport = 0, ndone = 0; i < N_gas+N_stars && nexport < All.BunchSizeChimie - NTask; i++)
{
/* only active particles and stars */
#ifndef CHIMIE_TIMEBET
if((P[i].Ti_endstep == All.Ti_Current)&&(P[i].Type == ST))
#else
if(P[i].Type == ST) /* as this is performed every ChimieTimeBetChimie, do it over all particles */
#endif
{
if(P[i].Type != ST)
{
printf("P[i].Type != ST, we better stop.\n");
printf("N_gas=%d (type=%d) i=%d (type=%d)\n",N_gas,P[N_gas].Type,i,P[i].Type);
printf("Please, check that you do not use PEANOHILBERT\n");
endrun(777001);
}
m = P[i].StPIdx;
#ifdef FOF_SFR
if (StP[m].NStars != 0)
#endif
if ( (P[i].Mass/StP[m].InitialMass) > MinRelMass)
{
flag_chimie = 0;
/******************************************/
/* do chimie */
/******************************************/
/*****************************************************/
/* look if a SN may have explode during the last step
/*****************************************************/
/***********************************************/
/***********************************************/
/* set the right table base of the metallicity */
set_table(0);
/***********************************************/
/***********************************************/
#ifdef FOF_SFR
/* minimum live time for a given metallicity */
minlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],StP[m].MassMax);
/* maximum live time for a given metallicity */
maxlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],StP[m].MassMin);
#else
/* minimum live time for a given metallicity */
minlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmax*All.MsoltoCMU);
/* maximum live time for a given metallicity */
maxlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmin*All.MsoltoCMU);
#endif
if (All.ComovingIntegrationOn)
{
/* FormationTime on the time line */
Ti0 = log(StP[m].FormationTime/All.TimeBegin) / All.Timebase_interval;
/* Beginning of time step on the time line */
#ifndef CHIMIE_TIMEBET
Ti1 = P[i].Ti_begstep;
#else
Ti1 = All.ChimieTi_begstep;
#endif
/* End of time step on the time line */
Ti2 = All.Ti_Current;
#ifdef COSMICTIME
t01 = get_cosmictime_difference(Ti0,Ti1);
t02 = get_cosmictime_difference(Ti0,Ti2);
#endif
}
else
{
#ifndef CHIMIE_TIMEBET
t1 = All.TimeBegin + (P[i].Ti_begstep * All.Timebase_interval);
#else
t1 = All.TimeBegin + (All.ChimieTi_begstep * All.Timebase_interval);
#endif
t2 = All.TimeBegin + (All.Ti_Current * All.Timebase_interval);
t01 = t1-StP[m].FormationTime;
t02 = t2-StP[m].FormationTime;
}
/* now treat all cases */
do_it=1;
#ifdef CHIMIE_ONE_SN_ONLY
if (All.Time<0.1)
do_it=0;
else
m1=m2=1; /* fix to 1 in order numerical problems with lifetime */
#else
/* beginning of interval */
if (t01>=minlivetime)
if (t01>=maxlivetime)
do_it=0; /* nothing to do */
else
m2 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t01);
else
m2 = Cp->Mmax*All.MsoltoCMU;
/* end of interval */
if (t02<=maxlivetime)
if (t02<=minlivetime)
do_it=0; /* nothing to do */
else
m1 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t02);
else
m1 = Cp->Mmin*All.MsoltoCMU;
#endif
//printf("Time=%g t01=%g t02=%g id=%d minlivetime=%g maxlivetime=%g \n",All.Time,t01,t02,P[i].ID,minlivetime,maxlivetime);
/* if some of the stars in the SSP explode between t1 and t2 */
if (do_it)
{
#ifdef CHIMIE_SN_ON_SINGLE_PART
printf("(%d) CHIMIE_SN_ON_SINGLE_PART : %d -> %d dist=%g",ThisTask,P[i].ID,StP[m].IDMinDist,StP[m].MinDist);
#endif
Nchimlocal++;
StP[m].Flag = 1; /* mark it as active */
if (m1>m2)
{
printf("m1=%g (%g Msol) > m2=%g (%g Msol) !!!\n\n",m1,m1*All.CMUtoMsol,m2,m2*All.CMUtoMsol);
endrun(777002);
}
M0 = StP[m].InitialMass;
for (k=0;k<NELEMENTS;k++)
metals[k] = StP[m].Metal[k];
#ifdef FOF_SFR
float current_mass;
NSNII = 0;
NSNIa = 0;
NDYIN = 0;
/* initialize */
StP[m].TotalEjectedGasMass = 0;
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] = 0;
/**************************************************************************
first case: we are dealing with a particle containing individual stars
**************************************************************************/
if (StP[m].NStars>0)
{
current_mass = StP[m].MassMax; /* start by the most massive star represented by the stellar particle */
printf("(%d) FOF_SFR : INDIVIDUAL STARS id=%8d %8.4g %8.4g %8d current_mass=%g MassSSP=%g\n",ThisTask,P[i].ID,StP[m].MassMin*All.CMUtoMsol,StP[m].MassMax*All.CMUtoMsol,StP[m].NStars,current_mass*All.CMUtoMsol,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol);
fprintf(FdFOF_Chimie,"%15g : id=%8d INDIVIDUAL STARS %8.4g %8.4g %8d current_mass=%g MassSSP=%g %g\n",All.Time,P[i].ID,StP[m].MassMin*All.CMUtoMsol,StP[m].MassMax*All.CMUtoMsol,StP[m].NStars,current_mass*All.CMUtoMsol,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol);
if(current_mass > m2)
{
printf(" current_mass = %g > m2 = %g, this seems odd (id=%d)!!!\n ",current_mass*All.CMUtoMsol,m2*All.CMUtoMsol,P[i].ID);
fprintf(FdFOF_Chimie,"%15g : id=%8d !!! current_mass = %g > m2 = %g, this seems odd !!!\n ",All.Time,P[i].ID,current_mass*All.CMUtoMsol,m2*All.CMUtoMsol);
endrun(888029);
}
while ( (current_mass > m1) && (StP[m].NStars>0))
{
/* stop the loop if the mass goes beyond the minimum stellar mass responsible of ejetas */
if (current_mass < Cp->MassMinForEjecta*All.MsoltoCMU)
{
StP[m].NStars=0;
printf("(%d) FOF_SFR : id=%8d current_mass=%g ---> mass below MassMinForEjecta=%g\n",ThisTask,P[i].ID,current_mass*All.CMUtoMsol,Cp->MassMinForEjecta);
fprintf(FdFOF_Chimie,"%15g : id=%8d current_mass=%g ---> mass below MassMinForEjecta=%g\n",All.Time,P[i].ID,current_mass*All.CMUtoMsol,Cp->MassMinForEjecta);
break;
}
/****************
SNII
****************/
printf("(%d) FOF_SFR : id=%8d current_mass=%g m1=%g\n",ThisTask,P[i].ID,current_mass*All.CMUtoMsol,m1*All.CMUtoMsol);
fprintf(FdFOF_Chimie,"%15g : id=%8d current_mass=%g m1=%g\n",All.Time,P[i].ID,current_mass*All.CMUtoMsol,m1*All.CMUtoMsol);
if ( (Cp->SNII_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->SNII_Mmax*All.MsoltoCMU) )
{
NSNII++;
//StP[m].NStars--;
printf("(%d) FOF_SFR : id=%8d one SNII\n",ThisTask,P[i].ID);
fprintf(FdFOF_Chimie,"%15g : id=%8d one SNII\n",All.Time,P[i].ID);
SNII_Total_single_mass_ejection(current_mass,metals);
StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] += SingleEjectedMass[k+2]; /* metal mass */
}
#ifdef CHIMIE_AGB
else
{
/****************
DYIN
****************/
if ( (Cp->DYIN_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->DYIN_Mmax*All.MsoltoCMU) )
{
NDYIN++;
//StP[m].NStars--;
printf("(%d) FOF_SFR : id=%8d one DYIN\n",ThisTask,P[i].ID);
fprintf(FdFOF_Chimie,"%15g : id=%8d one DYIN\n",All.Time,P[i].ID);
DYIN_Total_single_mass_ejection(current_mass,metals);
StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] += SingleEjectedMass[k+2]; /* metal mass */
}
}
#endif /* CHIMIE_AGB */
/****************
SNIa
****************/
NSNIa = NSNIa + SNIa_single_rate(current_mass);
/* find the next mass */
current_mass = optimal_get_next_mass_for_one_particle(m)*All.MsoltoCMU;
StP[m].OptIMF_CurrentMass = current_mass*All.CMUtoMsol;
StP[m].MassMax = current_mass;
/* remove the star */
StP[m].NStars--;
}
printf("(%d) FOF_SFR : id=%8d remaining stars= %8d \n",ThisTask,P[i].ID,StP[m].NStars);
fprintf(FdFOF_Chimie,"%15g : id=%8d remaining stars= %8d \n",All.Time,P[i].ID,StP[m].NStars);
// this is the standard way of doing
//NSNIa = SNIa_rate(m1,m2)*M0;
#ifdef CHIMIE_MC_SUPERNOVAE
double fNSNIa = NSNIa-floor(NSNIa);
NSNIa = floor(NSNIa);
if (get_Chimie_random_number(P[i].ID) < fNSNIa)
NSNIa = NSNIa+1;
#endif
/* compute ejectas for SNIa */
SNIa_Total_single_mass_ejection(0.5*(m1+m2),metals);
StP[m].TotalEjectedGasMass += SingleEjectedMass[0]*NSNIa; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] += SingleEjectedMass[k+2]*NSNIa; /* metal mass */
if (NSNIa>0)
fprintf(FdFOF_Chimie,"%15g : id=%8d NSNIa=%d, m1=%g m2=%g\n",All.Time,P[i].ID,NSNIa,m1/All.MsoltoCMU,m2/All.MsoltoCMU);
}
/**************************************************************************
second case: we are dealing with a particle containing a portion of IMF
**************************************************************************/
else
{
m2 = dmin(m2,StP[m].MassMax);
m1 = dmax(m1,StP[m].MassMin);
/* number of SNIa */
NSNIa = SNIa_rate(m1,m2)*StP[m].MassSSP;
/* number of SNII */
NSNII = SNII_rate(m1,m2)*StP[m].MassSSP;
/* number of DYIN */
#ifdef CHIMIE_AGB
NDYIN = DYIN_rate(m1,m2)*StP[m].MassSSP;
#endif
printf("(%d) FOF_SFR : IMF PORTION id=%8d Mass=%g MassSSP=%g NSNII=%g NDYIN=%g NSNIa=%g\n",ThisTask,P[i].ID,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol,NSNII,NDYIN,NSNIa);
fprintf(FdFOF_Chimie,"%15g : id=%8d IMF PORTION Mass=%g MassSSP=%g NSNII=%g NDYIN=%g NSNIa=%g\n",All.Time,P[i].ID,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol,NSNII,NDYIN,NSNIa);
#ifdef CHIMIE_MC_SUPERNOVAE
double fNSNIa,fNSNII,fNDYIN;
/* discretize SNIa */
fNSNIa = NSNIa-floor(NSNIa);
NSNIa = floor(NSNIa);
if (get_Chimie_random_number(P[i].ID) < fNSNIa)
NSNIa = NSNIa+1;
/* discretize SNII */
fNSNII = NSNII-floor(NSNII);
NSNII = floor(NSNII);
//if (get_Chimie_random_number(P[i].ID) < fNSNII)
// NSNII = NSNII+1;
/* discretize DYIN */
#ifdef CHIMIE_AGB
fNDYIN = NDYIN-floor(NDYIN);
NDYIN = floor(NDYIN);
//if (get_Chimie_random_number(P[i].ID) < fNDYIN)
// NDYIN = NDYIN+1;
#endif /* CHIMIE_AGB */
#endif
/* compute ejectas */
Total_single_mass_ejection(0.5*(m1+m2),metals,NSNII,NSNIa,NDYIN);
StP[m].TotalEjectedGasMass = SingleEjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] = SingleEjectedMass[k+2]; /* metal mass */
printf("(%d) FOF_SFR : IMF PORTION id=%8d NSNII=%g NDYIN=%g NSNIa=%g\n",ThisTask,P[i].ID,NSNII,NDYIN,NSNIa);
fprintf(FdFOF_Chimie,"%15g : id=%8d IMF PORTION NSNII=%g NDYIN=%g NSNIa=%g\n",All.Time,P[i].ID,NSNII,NDYIN,NSNIa);
}
#else /* FOF_SFR */
#ifdef CHIMIE_OPTIMAL_SAMPLING /* no FOF_SFR */
float current_mass;
NSNII = 0;
NSNIa = 0;
NDYIN = 0;
/* initialize */
StP[m].TotalEjectedGasMass = 0;
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] = 0;
current_mass = StP[m].OptIMF_CurrentMass*All.MsoltoCMU;
if(current_mass > m2)
{
printf(" current_mass = %g > m2 = %g, this seems odd (id=%d)!!!\n ",current_mass*All.CMUtoMsol,m2*All.CMUtoMsol,P[i].ID);
endrun(888027);
}
while (current_mass > m1)
{
/****************
SNII
****************/
if ( (Cp->SNII_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->SNII_Mmax*All.MsoltoCMU) )
{
NSNII++;
SNII_Total_single_mass_ejection(current_mass,metals);
StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] += SingleEjectedMass[k+2]; /* metal mass */
}
#ifdef CHIMIE_AGB
else
{
/****************
DYIN
****************/
if ( (Cp->DYIN_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->DYIN_Mmax*All.MsoltoCMU) )
{
NDYIN++;
DYIN_Total_single_mass_ejection(current_mass,metals);
StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] += SingleEjectedMass[k+2]; /* metal mass */
}
}
#endif
/* find the next mass */
current_mass = optimal_get_next_mass_for_one_particle(m)*All.MsoltoCMU;
StP[m].OptIMF_CurrentMass = current_mass*All.CMUtoMsol;
}
/****************
SNIa (similar than whith CHIMIE_MC_SUPERNOVAE)
****************/
/* compute the number of SNIa (this remains the same) */
NSNIa = SNIa_rate(m1,m2)*M0;
double fNSNIa = NSNIa-floor(NSNIa);
NSNIa = floor(NSNIa);
if (get_Chimie_random_number(P[i].ID) < fNSNIa)
NSNIa = NSNIa+1;
/* compute ejectas for SNIa */
SNIa_Total_single_mass_ejection(0.5*(m1+m2),metals);
StP[m].TotalEjectedGasMass += SingleEjectedMass[0]*NSNIa; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] += SingleEjectedMass[k+2]*NSNIa; /* metal mass */
#else /* no CHIMIE_OPTIMAL_SAMPLING */
/* number of SNIa */
NSNIa = SNIa_rate(m1,m2)*M0;
/* number of SNII */
NSNII = SNII_rate(m1,m2)*M0;
/* number of DYIN */
#ifdef CHIMIE_AGB
NDYIN = DYIN_rate(m1,m2)*M0;
#endif
#ifdef CHIMIE_MC_SUPERNOVAE
double fNSNIa,fNSNII,fNDYIN;
/* discretize SNIa */
fNSNIa = NSNIa-floor(NSNIa);
NSNIa = floor(NSNIa);
if (get_Chimie_random_number(P[i].ID) < fNSNIa)
NSNIa = NSNIa+1;
/* discretize SNII */
fNSNII = NSNII-floor(NSNII);
NSNII = floor(NSNII);
if (get_Chimie_random_number(P[i].ID) < fNSNII)
NSNII = NSNII+1;
/* discretize DYIN */
#ifdef CHIMIE_AGB
fNDYIN = NDYIN-floor(NDYIN);
NDYIN = floor(NDYIN);
if (get_Chimie_random_number(P[i].ID) < fNDYIN)
NDYIN = NDYIN+1;
#endif
#ifdef CHIMIE_ONE_SN_ONLY
/* here, we force to explode one and only one */
NSNIa=0;
NSNII=1;
NDYIN=0;
m1 = m2 = All.ChimieOneSNMass;
printf("CHIMIE_ONE_SN_ONLY : explode %g Msol SN.\n",m1*All.CMUtoMsol);
#endif
/* compute ejectas */
Total_single_mass_ejection(0.5*(m1+m2),metals,NSNII,NSNIa,NDYIN);
StP[m].TotalEjectedGasMass = SingleEjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] = SingleEjectedMass[k+2]; /* metal mass */
#else /* no CHIMIE_MC_SUPERNOVAE */
/* compute ejectas */
Total_mass_ejection(m1,m2,M0,metals);
StP[m].TotalEjectedGasMass = EjectedMass[0]; /* gas mass */
for (k=0;k<NELEMENTS;k++)
StP[m].TotalEjectedEltMass[k] = EjectedMass[k+2]; /* metal mass */
#endif /* CHIMIE_MC_SUPERNOVAE */
#endif /* CHIMIE_OPTIMAL_SAMPLING */
#endif /* FOF_SFR */
if (P[i].Mass-StP[m].TotalEjectedGasMass>0)
{
#ifdef CHIMIE_STATS
if (StP[m].TotalEjectedGasMass>0)
{
fprintf(FdChimieStatsSNs,"%4d %10d %8.4f %8d %10.3g %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f ",ThisTask,All.NumCurrentTiStep,All.Time, P[i].ID,P[i].Mass,NSNII,NSNIa,0.5*(m1+m2)*All.CMUtoMsol,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2]);
fprintf(FdChimieStatsSNs,"%10.3g ",StP[m].TotalEjectedGasMass);
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsSNs,"%10.3g ",StP[m].TotalEjectedEltMass[k]);
fprintf(FdChimieStatsSNs,"\n");
}
#endif
if (StP[m].TotalEjectedGasMass>0)
flag_chimie=1;
/* compute injected energy */
StP[m].TotalEjectedEgySpec = All.ChimieSupernovaEnergy* (NSNIa + NSNII) /StP[m].TotalEjectedGasMass;
StP[m].NumberOfSNIa = NSNIa;
StP[m].NumberOfSNII = NSNII;
EgySNlocal += All.ChimieSupernovaEnergy* (NSNIa + NSNII);
NSNIa_totlocal += NSNIa;
NSNII_totlocal += NSNII;
NDYIN_totlocal += NDYIN;
/* correct mass particle */
P[i].Mass = P[i].Mass-StP[m].TotalEjectedGasMass;
if(P[i].Mass<0)
endrun(777023);
}
else
{
printf("CHIMIE : mass wants to be less than zero...\n");
printf("CHIMIE : ID=%d Mass=%g StP[m].TotalEjectedGasMass=%g\n",P[i].ID,P[i].Mass*All.CMUtoMsol,StP[m].TotalEjectedGasMass*All.CMUtoMsol);
}
}
/******************************************/
/* end do chimie */
/******************************************/
ndone++;
if (flag_chimie)
{
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
chimie_evaluate(i, 0);
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
for(k = 0; k < 3; k++)
{
ChimieDataIn[nexport].Pos[k] = P[i].Pos[k];
ChimieDataIn[nexport].Vel[k] = P[i].Vel[k];
}
ChimieDataIn[nexport].ID = P[i].ID;
ChimieDataIn[nexport].Timestep = P[i].Ti_endstep - P[i].Ti_begstep;
ChimieDataIn[nexport].Hsml = StP[m].Hsml;
ChimieDataIn[nexport].Density = StP[m].Density;
ChimieDataIn[nexport].Y = StP[m].Y;
#ifdef CHIMIE_KINETIC_FEEDBACK
ChimieDataIn[nexport].NgbMass = StP[m].NgbMass;
#endif
ChimieDataIn[nexport].TotalEjectedGasMass = StP[m].TotalEjectedGasMass;
for(k = 0; k < NELEMENTS; k++)
ChimieDataIn[nexport].TotalEjectedEltMass[k] = StP[m].TotalEjectedEltMass[k];
ChimieDataIn[nexport].TotalEjectedEgySpec = StP[m].TotalEjectedEgySpec;
ChimieDataIn[nexport].NumberOfSNIa = StP[m].NumberOfSNIa;
ChimieDataIn[nexport].NumberOfSNII = StP[m].NumberOfSNII;
#ifdef WITH_ID_IN_HYDRA
ChimieDataIn[nexport].ID = P[i].ID;
#endif
#if CHIMIE_EJECTA_RADIUS == 2
ChimieDataIn[nexport].Pressure = StP[m].Pressure;
#endif
ChimieDataIn[nexport].Index = i;
ChimieDataIn[nexport].Task = j;
nexport++;
nsend_local[j]++;
}
}
}
}
}
}
tend = second();
timecomp += timediff(tstart, tend);
qsort(ChimieDataIn, nexport, sizeof(struct chimiedata_in), chimie_compare_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
tstart = second();
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeChimie)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
/* get the particles */
MPI_Sendrecv(&ChimieDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct chimiedata_in), MPI_BYTE,
recvTask, TAG_CHIMIE_A,
&ChimieDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct chimiedata_in), MPI_BYTE,
recvTask, TAG_CHIMIE_A, MPI_COMM_WORLD, &status);
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
/* now do the imported particles */
tstart = second();
for(j = 0; j < nbuffer[ThisTask]; j++)
chimie_evaluate(j, 1);
tend = second();
timecomp += timediff(tstart, tend);
/* do a block to measure imbalance */
tstart = second();
MPI_Barrier(MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* get the result */
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeChimie)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
/* send the results */
MPI_Sendrecv(&ChimieDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct chimiedata_out),
MPI_BYTE, recvTask, TAG_CHIMIE_B,
&ChimieDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct chimiedata_out),
MPI_BYTE, recvTask, TAG_CHIMIE_B, MPI_COMM_WORLD, &status);
/* add the result to the particles */
for(j = 0; j < nsend_local[recvTask]; j++)
{
source = j + noffset[recvTask];
place = ChimieDataIn[source].Index;
// for(k = 0; k < 3; k++)
// SphP[place].HydroAccel[k] += HydroDataPartialResult[source].Acc[k];
//
// SphP[place].DtEntropy += HydroDataPartialResult[source].DtEntropy;
//#ifdef FEEDBACK
// SphP[place].DtEgySpecFeedback += HydroDataPartialResult[source].DtEgySpecFeedback;
//#endif
// if(SphP[place].MaxSignalVel < HydroDataPartialResult[source].MaxSignalVel)
// SphP[place].MaxSignalVel = HydroDataPartialResult[source].MaxSignalVel;
//#ifdef COMPUTE_VELOCITY_DISPERSION
// for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
// SphP[place].VelocityDispersion[k] += HydroDataPartialResult[source].VelocityDispersion[k];
//#endif
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
level = ngrp - 1;
}
MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
for(j = 0; j < NTask; j++)
ntotleft -= ndonelist[j];
}
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
/* do final operations on results */
tstart = second();
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
P[i].Mass += SphP[i].dMass;
SphP[i].dMass = 0.;
}
}
tend = second();
timecomp += timediff(tstart, tend);
/* collect some timing information */
MPI_Reduce(&timecomp, &sumt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timecommsumm, &sumcomm, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timeimbalance, &sumimbalance, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if(ThisTask == 0)
{
All.CPU_ChimieCompWalk += sumt / NTask;
All.CPU_ChimieCommSumm += sumcomm / NTask;
All.CPU_ChimieImbalance += sumimbalance / NTask;
}
#ifdef DETAILED_CPU_OUTPUT_IN_CHIMIE
numlist = malloc(sizeof(int) * NTask);
timecomplist = malloc(sizeof(double) * NTask);
timecommsummlist = malloc(sizeof(double) * NTask);
timeimbalancelist = malloc(sizeof(double) * NTask);
MPI_Gather(&NumStUpdate, 1, MPI_INT, numlist, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Gather(&timecomp, 1, MPI_DOUBLE, timecomplist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(&timecommsumm, 1, MPI_DOUBLE, timecommsummlist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(&timeimbalance, 1, MPI_DOUBLE, timeimbalancelist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(ThisTask == 0)
{
fprintf(FdTimings, "\n chimie\n\n");
fprintf(FdTimings, "Nupdate ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12d ",numlist[i]); /* nombre de part par proc */
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timecomp ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timecomplist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timecommsumm ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timecommsummlist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timeimbalance ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timeimbalancelist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "\n");
}
free(timeimbalancelist);
free(timecommsummlist);
free(timecomplist);
free(numlist);
#endif
/* collect some chimie informations */
MPI_Reduce(&NSNIa_totlocal, &NSNIa_tot, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&NSNII_totlocal, &NSNII_tot, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&NDYIN_totlocal, &NDYIN_tot, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&EgySNlocal, &EgySN, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&Nchimlocal, &Nchim, 1, MPI_INT , MPI_SUM, 0, MPI_COMM_WORLD);
#ifdef CHIMIE_THERMAL_FEEDBACK
EgySNThermal = EgySN*(1-All.ChimieKineticFeedbackFraction);
#else
EgySNThermal = 0;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
EgySNKinetic = EgySN*All.ChimieKineticFeedbackFraction;
/* count number of wind particles */
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (SphP[i].WindTime >= (All.Time-All.ChimieWindTime))
Nwindlocal++;
//else
// if (SphP[i].WindTime > All.TimeBegin-2*All.ChimieWindTime)
// Noldwindlocal++;
if (SphP[i].WindFlag)
Nflaglocal++;
}
}
MPI_Reduce(&Nwindlocal, &Nwind, 1, MPI_INT , MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&Noldwindlocal, &Noldwind, 1, MPI_INT , MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Allreduce(&Nflaglocal, &Nflag, 1, MPI_INT , MPI_SUM, MPI_COMM_WORLD);
#else
EgySNKinetic = 0;
#endif
/* write some info */
if (ThisTask==0)
{
fprintf(FdChimie, "%15g %10d %15g %15g %15g %15g %15g %10d %10d %10d\n",All.Time,Nchim,NSNIa_tot,NSNII_tot,EgySN,EgySNThermal,EgySNKinetic,Nwind,Noldwind,Nflag);
fflush(FdChimie);
}
/* this is no longer used */
// if (Nflag>0)
// {
// SetMinTimeStepForActives=1;
// if (ThisTask==0)
// fprintf(FdLog,"%g : !!! set min timestep for active particles !!!\n",All.Time);
// }
#ifdef CHIMIE_ONE_SN_ONLY
if (EgySN>0)
All.ChimieOneSN=1;
MPI_Bcast(&All.ChimieOneSN, 1, MPI_INT, 0, MPI_COMM_WORLD);
#endif
}
/*! This function is the 'core' of the Chemie computation. A target
* particle is specified which may either be local, or reside in the
* communication buffer.
*/
void chimie_evaluate(int target, int mode)
{
int j, n, startnode, numngb,numngb_inbox,k;
FLOAT *pos,*vel;
//FLOAT *vel;
//FLOAT mass;
double h, h2;
double acc[3];
double dx, dy, dz;
double wk, r, r2, u=0;
double hinv=1, hinv3;
int target_stp;
double density;
double Y;
#ifdef CHIMIE_KINETIC_FEEDBACK
double ngbmass;
double p;
#endif
double aij;
double ejectedGasMass;
double ejectedEltMass[NELEMENTS];
double ejectedEgySpec;
double NumberOfSNIa;
double NumberOfSNII;
double mass_k;
double NewMass;
double fv,vi2,vj2;
double EgySpec,NewEgySpec;
double DeltaEntropy;
double DeltaVel[3];
#ifndef LONGIDS
unsigned int id; /*!< particle identifier */
#else
unsigned long long id; /*!< particle identifier */
#endif
#if CHIMIE_EJECTA_RADIUS == 2
double rho,E51,n0,P04,pressure,P0;
#endif
if(mode == 0)
{
pos = P[target].Pos;
vel = P[target].Vel;
id = P[target].ID;
target_stp = P[target].StPIdx;
h = StP[target_stp].Hsml;
density = StP[target_stp].Density;
Y = StP[target_stp].Y;
#ifdef CHIMIE_KINETIC_FEEDBACK
ngbmass = StP[target_stp].NgbMass;
#endif
ejectedGasMass = StP[target_stp].TotalEjectedGasMass;
for(k=0;k<NELEMENTS;k++)
ejectedEltMass[k] = StP[target_stp].TotalEjectedEltMass[k];
ejectedEgySpec = StP[target_stp].TotalEjectedEgySpec;
NumberOfSNIa = StP[target_stp].NumberOfSNIa;
NumberOfSNII = StP[target_stp].NumberOfSNII;
#if CHIMIE_EJECTA_RADIUS == 2
pressure = StP[target_stp].Pressure;
#endif
}
else
{
pos = ChimieDataGet[target].Pos;
vel = ChimieDataGet[target].Vel;
id = ChimieDataGet[target].ID;
h = ChimieDataGet[target].Hsml;
density = ChimieDataGet[target].Density;
Y = ChimieDataGet[target].Y;
#ifdef CHIMIE_KINETIC_FEEDBACK
ngbmass = ChimieDataGet[target].NgbMass;
#endif
ejectedGasMass = ChimieDataGet[target].TotalEjectedGasMass;
for(k=0;k<NELEMENTS;k++)
ejectedEltMass[k] = ChimieDataGet[target].TotalEjectedEltMass[k];
ejectedEgySpec = ChimieDataGet[target].TotalEjectedEgySpec;
NumberOfSNIa = ChimieDataGet[target].NumberOfSNIa;
NumberOfSNII = ChimieDataGet[target].NumberOfSNII;
#if CHIMIE_EJECTA_RADIUS == 2
pressure = ChimieDataGet[target].Pressure;
#endif
}
#if CHIMIE_EJECTA_RADIUS == 1
h = All.ChimieEjectaRadius;
#elif CHIMIE_EJECTA_RADIUS == 2 /* Maximum Blast radius computed according to Chevalier 1974 */
/* density in cgs */
rho = density *All.UnitDensity_in_cgs*(All.HubbleParam*All.HubbleParam);
/* pressure in cgs */
P0 = pressure * All.UnitDensity_in_cgs*(All.HubbleParam*All.HubbleParam) * All.UnitEnergy_in_cgs/All.UnitMass_in_g;
E51 = 1.0;
n0 = rho/PROTONMASS;
P04 = 1e-4 * P0/BOLTZMANN;
h = pow(10,1.74) * pow(E51,0.32) * pow(n0,-0.16) * pow(P04,-0.2); /* [in pc] Chevalier 1974 */
h = h/1e3; /* [in kpc] */
h = dmin(h,1);
h = h* KPC_IN_CM/All.UnitLength_in_cm * All.HubbleParam; /* [in code unit length] */
#endif
/* initialize variables before SPH loop is started */
acc[0] = acc[1] = acc[2] = 0;
vi2 = 0;
for(k=0;k<3;k++)
vi2 += vel[k]*vel[k];
h2 = h * h;
hinv = 1.0 / h;
#ifndef TWODIMS
hinv3 = hinv * hinv * hinv;
#else
hinv3 = hinv * hinv / boxSize_Z;
#endif
/* Now start the actual SPH computation for this particle */
startnode = All.MaxPart;
numngb = 0;
do
{
numngb_inbox = ngb_treefind_variable_for_chimie(&pos[0], h, &startnode);
for(n = 0; n < numngb_inbox; n++)
{
j = Ngblist[n];
dx = pos[0] - P[j].Pos[0];
dy = pos[1] - P[j].Pos[1];
dz = pos[2] - P[j].Pos[2];
#ifdef PERIODIC /* now find the closest image in the given box size */
if(dx > boxHalf_X)
dx -= boxSize_X;
if(dx < -boxHalf_X)
dx += boxSize_X;
if(dy > boxHalf_Y)
dy -= boxSize_Y;
if(dy < -boxHalf_Y)
dy += boxSize_Y;
if(dz > boxHalf_Z)
dz -= boxSize_Z;
if(dz < -boxHalf_Z)
dz += boxSize_Z;
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(r2 < h2)
{
numngb++;
r = sqrt(r2);
u = r * hinv;
if(u < 0.5)
{
wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
}
else
{
wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
}
/* normalisation using mass */
//aij = P[j].Mass*wk/density;
/* normalisation using volume */
/* !!! si on utilise, il faut stoquer une nouvelle variable : OldDensity, car density est modifié plus bas... */
//aij = P[j].Mass/SphP[j].Density*wk/volume;
/* !!! we shoud better use an OldDensity here, as .Density may slighly change after */
#if CHIMIE_ELTS_WEIGHTING == 0
aij = P[j].Mass * wk; // 0: mass + kernel weighting (default, Revaz & Jablonka 2012)
#elif CHIMIE_ELTS_WEIGHTING == 1
aij = P[j].Mass/SphP[j].Density * wk; // 1: volume + kernel weighting
#elif CHIMIE_ELTS_WEIGHTING == 2
aij = wk; // 2: number + kernel weighting
#elif CHIMIE_ELTS_WEIGHTING == 3
aij = P[j].Mass; // 3: mass weighting
#elif CHIMIE_ELTS_WEIGHTING == 4
aij = P[j].Mass/SphP[j].Density; // 4: volume weighting
#elif CHIMIE_ELTS_WEIGHTING == 5
aij = 1; // 5: number weighting
#endif
aij = aij/Y; // normalisation
#ifdef CHIMIE_STATS /* before the interaction */
fprintf(FdChimieStatsGas,"%4d %10d %8.4f %8d %8d %10.3g %10.3g %10.3g %10.3g ",ThisTask,All.NumCurrentTiStep,All.Time, P[j].ID,id,aij,P[j].Pos[0],P[j].Pos[1],P[j].Pos[2]);
fprintf(FdChimieStatsGas,"%10.3g ",P[j].Mass);
#ifdef CHIMIE_SMOOTH_METALS
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsGas,"%10.3g ",SphP[j].MassMetal[k]*P[j].Mass);
#else
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsGas,"%10.3g ",SphP[j].Metal[k]*P[j].Mass);
#endif
fprintf(FdChimieStatsGas," ");
#endif
/* metal injection */
for(k=0;k<NELEMENTS;k++)
{
#ifdef CHIMIE_SMOOTH_METALS
mass_k = SphP[j].MassMetal[k]*P[j].Mass; /* mass of elt k */
SphP[j].MassMetal[k] = ( mass_k + aij*ejectedEltMass[k] )/( P[j].Mass + aij*ejectedGasMass );
#else
mass_k = SphP[j].Metal[k]*P[j].Mass; /* mass of elt k */
SphP[j].Metal[k] = ( mass_k + aij*ejectedEltMass[k] )/( P[j].Mass + aij*ejectedGasMass );
#endif
}
/* new mass */
NewMass = P[j].Mass + aij*ejectedGasMass;
#ifdef CHIMIE_STATS /* after the interaction !!! slighly wrong, as the mass is not the final ones (could receive mass from other particles) */
fprintf(FdChimieStatsGas,"%10.3g ",NewMass);
#ifdef CHIMIE_SMOOTH_METALS
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsGas,"%10.3g ",SphP[j].MassMetal[k]*NewMass);
#else
for (k=0;k<NELEMENTS;k++)
fprintf(FdChimieStatsGas,"%10.3g ",SphP[j].Metal[k]*NewMass);
#endif
fprintf(FdChimieStatsGas,"\n");
#endif
/* new velocity */
vj2 = 0;
for(k=0;k<3;k++)
vj2 += SphP[j].VelPred[k]*SphP[j].VelPred[k];
fv = sqrt( (P[j].Mass/NewMass) + aij*(ejectedGasMass/NewMass) * (vi2/vj2) );
for(k=0;k<3;k++)
{
DeltaVel[k] = fv*SphP[j].VelPred[k] - SphP[j].VelPred[k];
SphP[j].VelPred[k] += DeltaVel[k];
P[j].Vel [k] += DeltaVel[k];
}
/* spec energy at current step */
#ifdef DENSITY_INDEPENDENT_SPH
EgySpec = SphP[j].EntropyPred / GAMMA_MINUS1 * pow(SphP[j].EgyWtDensity*a3inv, GAMMA_MINUS1);
#else
EgySpec = SphP[j].EntropyPred / GAMMA_MINUS1 * pow(SphP[j].Density*a3inv, GAMMA_MINUS1);
#endif
/* new egyspec */
NewEgySpec = (EgySpec )*(P[j].Mass/NewMass);
#ifdef CHIMIE_TIMEBET
/* avoid problematic values */ /* added with the new chimie timing */
if(NewEgySpec<All.MinEgySpec)
NewEgySpec = All.MinEgySpec;
#endif
/* new density */
#ifdef DENSITY_INDEPENDENT_SPH
SphP[j].Density = SphP[j].Density*NewMass/P[j].Mass;
SphP[j].EgyWtDensity = SphP[j].EgyWtDensity*NewMass/P[j].Mass;
#else
SphP[j].Density = SphP[j].Density*NewMass/P[j].Mass;
#endif
/* new entropy */
#ifdef DENSITY_INDEPENDENT_SPH
DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[j].EgyWtDensity*a3inv, GAMMA_MINUS1) - SphP[j].EntropyPred;
#else
DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[j].Density*a3inv, GAMMA_MINUS1) - SphP[j].EntropyPred;
#endif
SphP[j].EntropyPred += DeltaEntropy;
SphP[j].Entropy += DeltaEntropy;
#ifdef CHIMIE_THERMAL_FEEDBACK
SphP[j].DeltaEgySpec += (1.-All.ChimieKineticFeedbackFraction)*(ejectedGasMass*ejectedEgySpec)* aij/NewMass;
SphP[j].NumberOfSNII += NumberOfSNII*aij;
SphP[j].NumberOfSNIa += NumberOfSNIa*aij;
#ifdef TIMESTEP_UPDATE_FOR_FEEDBACK
if(P[j].Ti_endstep != All.Ti_Current)
make_particle_active(j);
#endif
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
p = (All.ChimieKineticFeedbackFraction*ejectedEgySpec*ejectedGasMass)/(0.5*ngbmass*All.ChimieWindSpeed*All.ChimieWindSpeed);
double r;
r = get_Chimie_random_number(P[j].ID+id);
if ( r < p) /* we should maybe have a 2d table here... */
{
if (SphP[j].WindTime < (All.Time-All.ChimieWindTime)) /* not a wind particle */
{
SphP[j].WindFlag = 1;
SphP[j].WindTime = All.Time;
}
}
#endif
#ifdef CHECK_ENTROPY_SIGN
if ((SphP[j].EntropyPred < 0)||(SphP[j].Entropy < 0))
{
printf("\ntask=%d: entropy less than zero in chimie_evaluate !\n", ThisTask);
printf("ID=%d Entropy=%g EntropyPred=%g DeltaEntropy=%g\n",P[j].ID,SphP[j].Entropy,SphP[j].EntropyPred,DeltaEntropy);
fflush(stdout);
endrun(777003);
}
#endif
/* store mass diff. */
SphP[j].dMass += NewMass-P[j].Mass;
}
}
}
while(startnode >= 0);
/* Now collect the result at the right place */
if(mode == 0)
{
// for(k = 0; k < 3; k++)
// SphP[target].HydroAccel[k] = acc[k];
// SphP[target].DtEntropy = dtEntropy;
//#ifdef FEEDBACK
// SphP[target].DtEgySpecFeedback = dtEgySpecFeedback;
//#endif
// SphP[target].MaxSignalVel = maxSignalVel;
//#ifdef COMPUTE_VELOCITY_DISPERSION
// for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
// SphP[target].VelocityDispersion[k] = VelocityDispersion[k];
//#endif
}
else
{
// for(k = 0; k < 3; k++)
// HydroDataResult[target].Acc[k] = acc[k];
// HydroDataResult[target].DtEntropy = dtEntropy;
//#ifdef FEEDBACK
// HydroDataResult[target].DtEgySpecFeedback = dtEgySpecFeedback;
//#endif
// HydroDataResult[target].MaxSignalVel = maxSignalVel;
//#ifdef COMPUTE_VELOCITY_DISPERSION
// for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
// HydroDataResult[target].VelocityDispersion[k] = VelocityDispersion[k];
//#endif
}
}
/*! This is a comparison kernel for a sort routine, which is used to group
* particles that are going to be exported to the same CPU.
*/
int chimie_compare_key(const void *a, const void *b)
{
if(((struct chimiedata_in *) a)->Task < (((struct chimiedata_in *) b)->Task))
return -1;
if(((struct chimiedata_in *) a)->Task > (((struct chimiedata_in *) b)->Task))
return +1;
return 0;
}
/****************************************************************************************/
/*
/*
/*
/* PYTHON INTERFACE
/*
/*
/*
/****************************************************************************************/
#ifdef PYCHEM
static PyObject *
chemistry_CodeUnits_to_SolarMass_Factor(PyObject *self, PyObject *args)
{
return Py_BuildValue("d",All.CMUtoMsol);
}
static PyObject *
chemistry_SolarMass_to_CodeUnits_Factor(PyObject *self, PyObject *args)
{
return Py_BuildValue("d",All.MsoltoCMU);
}
static PyObject *
chemistry_SetVerbosityOn(PyObject *self, PyObject *args)
{
verbose=1;
return Py_BuildValue("i",0);
}
static PyObject *
chemistry_SetVerbosityOff(PyObject *self, PyObject *args)
{
verbose=0;
return Py_BuildValue("i",0);
}
static PyObject * chemistry_InitDefaultParameters(void)
{
/* list of Gadget parameters */
/* System of units */
All.UnitLength_in_cm = 3.085e+21; /* 1.0 kpc */
All.UnitMass_in_g = 1.989e+43; /* 1.0e10 solar masses */
All.UnitVelocity_in_cm_per_s = 20725573.785998672; /* 207 km/sec */
All.GravityConstantInternal = 0;
All.HubbleParam = 1;
/* other usefull constants */
All.UnitTime_in_s = All.UnitLength_in_cm / All.UnitVelocity_in_cm_per_s;
All.UnitTime_in_Megayears=All.UnitTime_in_s / SEC_PER_MEGAYEAR;
All.CMUtoMsol = All.UnitMass_in_g/SOLAR_MASS/All.HubbleParam; /* convertion factor from Code Mass Unit to Solar Mass */
All.MsoltoCMU = 1/All.CMUtoMsol; /* convertion factor from Solar Mass to Code Mass Unit */
return Py_BuildValue("i",1);
}
static PyObject * SetParameters(PyObject *dict)
{
PyObject *key;
PyObject *value;
int ivalue;
float fvalue;
double dvalue;
/* check that it is a PyDictObject */
if(!PyDict_Check(dict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
if (PyDict_Size(dict)==0)
return Py_BuildValue("i",0);
Py_ssize_t pos=0;
while(PyDict_Next(dict,&pos,&key,&value))
{
if(PyString_Check(key))
{
/* System of units */
if(strcmp(PyString_AsString(key), "UnitLength_in_cm")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitLength_in_cm = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitMass_in_g")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitMass_in_g = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitVelocity_in_cm_per_s")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitVelocity_in_cm_per_s = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "HubbleParam")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.HubbleParam = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "GravityConstantInternal")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.GravityConstantInternal = PyFloat_AsDouble(value);
}
}
}
/* other usefull constants */
All.UnitTime_in_s = All.UnitLength_in_cm / All.UnitVelocity_in_cm_per_s;
All.UnitTime_in_Megayears=All.UnitTime_in_s / SEC_PER_MEGAYEAR;
All.CMUtoMsol = All.UnitMass_in_g/SOLAR_MASS/All.HubbleParam; /* convertion factor from Code Mass Unit to Solar Mass */
All.MsoltoCMU = 1/All.CMUtoMsol; /* convertion factor from Solar Mass to Code Mass Unit */
return Py_BuildValue("i",1);
}
static PyObject * chemistry_SetParameters(PyObject *self, PyObject *args)
{
PyObject *dict;
/* here, we can have either arguments or dict directly */
if(PyDict_Check(args))
{
dict = args;
}
else
{
if (! PyArg_ParseTuple(args, "O",&dict))
return NULL;
}
SetParameters(dict);
return Py_BuildValue("i",1);
}
static PyObject * chemistry_GetParameters(void)
{
PyObject *dict;
PyObject *key;
PyObject *value;
dict = PyDict_New();
/* System of units */
key = PyString_FromString("UnitLength_in_cm");
value = PyFloat_FromDouble(All.UnitLength_in_cm);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitMass_in_g");
value = PyFloat_FromDouble(All.UnitMass_in_g);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitVelocity_in_cm_per_s");
value = PyFloat_FromDouble(All.UnitVelocity_in_cm_per_s);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("HubbleParam");
value = PyFloat_FromDouble(All.HubbleParam);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("GravityConstantInternal");
value = PyFloat_FromDouble(All.GravityConstantInternal);
PyDict_SetItem(dict,key,value);
return Py_BuildValue("O",dict);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_init_chimie(PyObject *self, PyObject *args, PyObject *kwds)
{
int NumberOfTables=1;
int DefaultTable=0;
PyObject *paramsDict=NULL;
paramsDict= PyDict_New();
//PyObject *filename;
//if (! PyArg_ParseTuple(args, "Oii",&filename,&NumberOfTables,&DefaultTable))
// {
// PyErr_SetString(PyExc_ValueError,"init_chimie, error in parsing.");
// return NULL;
// }
static char *kwlist[] = {"filename","NumberOfTables","DefaultTable","params", NULL};
PyObject *filename=PyString_FromString("chimie.yr.dat");
/* this fails with python2.6, I do not know why ??? */
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OiiO",kwlist,&filename,&NumberOfTables,&DefaultTable,&paramsDict))
{
PyErr_SetString(PyExc_ValueError,"init_chimie, error in parsing arguments.");
return NULL;
}
if (!PyString_Check(filename))
{
PyErr_SetString(PyExc_ValueError,"Argument must be a string.");
return NULL;
}
/* copy filename */
All.ChimieParameterFile = PyString_AsString(filename);
/* set number of tables */
All.ChimieNumberOfParameterFiles = NumberOfTables;
/* check if the file exists */
if(!(fopen(All.ChimieParameterFile, "r")))
{
PyErr_SetString(PyExc_ValueError,"The parameter file does not exists.");
return NULL;
}
/* use default parameters */
chemistry_InitDefaultParameters();
/* check if units are given */
/* check that it is a PyDictObject */
if(!PyDict_Check(paramsDict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
else
{
SetParameters(paramsDict);
}
init_chimie();
/* by default, set the first one */
set_table(DefaultTable);
return Py_BuildValue("O",Py_None);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_info(PyObject *self, PyObject *args, PyObject *kwds)
{
int DefaultTable=0;
static char *kwlist[] = {"DefaultTable", NULL};
/* this fails with python2.6, I do not know why ??? */
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|i",kwlist,&DefaultTable))
{
PyErr_SetString(PyExc_ValueError,"init_chimie, error in parsing arguments.");
return NULL;
}
info(0);
return Py_BuildValue("O",Py_None);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_set_table(PyObject *self, PyObject *args, PyObject *kwds)
{
int i;
if (! PyArg_ParseTuple(args, "i",&i))
return PyString_FromString("error");
/* set the table */
set_table(i);
return Py_BuildValue("d",0);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_get_imf(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *m,*imf;
int i;
if (! PyArg_ParseTuple(args, "O",&m))
return PyString_FromString("error");
m = TO_DOUBLE(m);
/* create an output */
imf = (PyArrayObject *) PyArray_SimpleNew(m->nd,m->dimensions,PyArray_DOUBLE);
//printf("--> %g\n",Cp->bs[0]);
//for (i=0;i<Cp->n;i++)
// printf("%g %g\n",Cp->ms[i],Cp->as[i]);
for(i = 0; i < m->dimensions[0]; i++)
{
*(double *)(imf->data + i*(imf->strides[0])) = get_imf(*(double *)(m->data + i*(m->strides[0])));
}
return PyArray_Return(imf);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_get_imf_M(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *m1,*m2,*imf;
int i;
if (! PyArg_ParseTuple(args, "OO",&m1,&m2))
return PyString_FromString("error");
m1 = TO_DOUBLE(m1);
m2 = TO_DOUBLE(m2);
/* create an output */
imf = (PyArrayObject *) PyArray_SimpleNew(m1->nd,m1->dimensions,PyArray_DOUBLE);
for(i = 0; i < imf->dimensions[0]; i++)
{
*(double *)(imf->data + i*(imf->strides[0])) = get_imf_M( *(double *)(m1->data + i*(m1->strides[0])), *(double *)(m2->data + i*(m2->strides[0])) );
}
return PyArray_Return(imf);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_get_imf_N(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *m1,*m2,*imf;
int i;
if (! PyArg_ParseTuple(args, "OO",&m1,&m2))
return PyString_FromString("error");
m1 = TO_DOUBLE(m1);
m2 = TO_DOUBLE(m2);
/* create an output */
imf = (PyArrayObject *) PyArray_SimpleNew(m1->nd,m1->dimensions,PyArray_DOUBLE);
for(i = 0; i < imf->dimensions[0]; i++)
{
*(double *)(imf->data + i*(imf->strides[0])) = get_imf_N( *(double *)(m1->data + i*(m1->strides[0])), *(double *)(m2->data + i*(m2->strides[0])) );
}
return PyArray_Return(imf);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_star_lifetime(self, args)
PyObject *self;
PyObject *args;
{
/* z is the mass fraction of metals, ie, the metallicity */
/* m is the star mass in code unit */
/* Return t in time unit */
double time,z,m;
if (!PyArg_ParseTuple(args, "dd", &z, &m))
return NULL;
time = star_lifetime(z,m);
return Py_BuildValue("d",time);
}
static PyObject *
chemistry_star_mass_from_age(self, args)
PyObject *self;
PyObject *args;
{
/* t : life time (in code unit) */
/* return the stellar mass (in code unit) that has a lifetime equal to t */
double time,z,m;
if (!PyArg_ParseTuple(args, "dd", &z, &time))
return NULL;
m = star_mass_from_age(z,time);
return Py_BuildValue("d",m);
}
static PyObject *
chemistry_DYIN_rate(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2;
double RDYIN;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1,&m2))
return NULL;
RDYIN = DYIN_rate(m1,m2);
return Py_BuildValue("d",RDYIN);
}
static PyObject *
chemistry_SNII_rate(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2;
double RSNII;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1,&m2))
return NULL;
RSNII = SNII_rate(m1,m2);
return Py_BuildValue("d",RSNII);
}
static PyObject *
chemistry_SNIa_rate(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2;
double RSNIa;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1,&m2))
return NULL;
RSNIa = SNIa_rate(m1,m2);
return Py_BuildValue("d",RSNIa);
}
static PyObject *
chemistry_SNIa_single_rate(self, args)
PyObject *self;
PyObject *args;
{
double m;
double RSNIa;
/* parse arguments */
if (!PyArg_ParseTuple(args, "d", &m))
return NULL;
RSNIa = SNIa_single_rate(m);
return Py_BuildValue("d",RSNIa);
}
static PyObject *
chemistry_DYIN_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2;
PyArrayObject *ArrMassDYIN;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1,&m2))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute dying stars ejection */
DYIN_mass_ejection(m1,m2);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = MassFracDYIN[i];
/* convert in array */
return PyArray_Return(ArrMassDYIN);
}
static PyObject *
chemistry_DYIN_mass_ejection_Z(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2,Z;
PyArrayObject *ArrMassDYIN;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "ddd", &m1,&m2,&Z))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute dying stars ejection */
DYIN_mass_ejection_Z(m1,m2, Z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = MassFracDYIN[i];
/* convert in array */
return PyArray_Return(ArrMassDYIN);
}
static PyObject *
chemistry_DYIN_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
PyArrayObject *ArrMassDYIN;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "d", &m1))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
DYIN_single_mass_ejection(m1);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = SingleMassFracDYIN[i];
/* convert in array */
return PyArray_Return(ArrMassDYIN);
}
static PyObject *
chemistry_DYIN_single_mass_ejection_Z(self, args)
PyObject *self;
PyObject *args;
{
double m1,Z;
PyArrayObject *ArrMassDYIN;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1, &Z))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
DYIN_single_mass_ejection_Z(m1,Z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = SingleMassFracDYIN[i];
/* convert in array */
return PyArray_Return(ArrMassDYIN);
}
static PyObject *
chemistry_SNII_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2;
PyArrayObject *ArrMassSNII;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1,&m2))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
SNII_mass_ejection(m1,m2);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = MassFracSNII[i];
/* convert in array */
return PyArray_Return(ArrMassSNII);
}
static PyObject *
chemistry_SNII_mass_ejection_Z(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2,Z;
PyArrayObject *ArrMassSNII;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "ddd", &m1,&m2,&Z))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SNII stars ejection */
SNII_mass_ejection_Z(m1,m2, Z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = MassFracSNII[i];
/* convert in array */
return PyArray_Return(ArrMassSNII);
}
static PyObject *
chemistry_SNII_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
PyArrayObject *ArrMassSNII;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "d", &m1))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
SNII_single_mass_ejection(m1);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = SingleMassFracSNII[i];
/* convert in array */
return PyArray_Return(ArrMassSNII);
}
static PyObject *
chemistry_SNII_single_mass_ejection_Z(self, args)
PyObject *self;
PyObject *args;
{
double m1,Z;
PyArrayObject *ArrMassSNII;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1, &Z))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
SNII_single_mass_ejection_Z(m1,Z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = SingleMassFracSNII[i];
/* convert in array */
return PyArray_Return(ArrMassSNII);
}
static PyObject *
chemistry_SNIa_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2;
PyArrayObject *ArrMassSNIa;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m1,&m2))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassSNIa = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
SNIa_mass_ejection(m1,m2);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassSNIa->data + (i)*(ArrMassSNIa->strides[0])) = MassFracSNIa[i];
/* convert in array */
return PyArray_Return(ArrMassSNIa);
}
static PyObject *
chemistry_SNIa_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
PyArrayObject *ArrMassSNIa;
npy_intp ld[1];
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "d", &m1))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
ArrMassSNIa = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* compute SN ejection */
SNIa_single_mass_ejection(m1);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(ArrMassSNIa->data + (i)*(ArrMassSNIa->strides[0])) = SingleMassFracSNIa[i];
/* convert in array */
return PyArray_Return(ArrMassSNIa);
}
static PyObject *
chemistry_Total_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2,M;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dddO", &m1,&m2,&M,&zs))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute SN ejection */
Total_mass_ejection(m1,m2,M,z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = EjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
static PyObject *
chemistry_Total_mass_ejection_Z(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2,M;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dddO", &m1,&m2,&M,&zs))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute SN ejection */
Total_mass_ejection_Z(m1,m2,M,z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = EjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
static PyObject *
chemistry_DYIN_Total_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dO", &m1,&zs))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute dying stars ejection */
DYIN_Total_single_mass_ejection(m1,z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
static PyObject *
chemistry_SNII_Total_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dO", &m1,&zs))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute SN ejection */
SNII_Total_single_mass_ejection(m1,z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
static PyObject *
chemistry_SNIa_Total_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dO", &m1,&zs))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute SN ejection */
SNIa_Total_single_mass_ejection(m1,z);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
static PyObject *
chemistry_Total_single_mass_ejection(self, args)
PyObject *self;
PyObject *args;
{
double m1;
double NSNII,NSNIa,NDYIN;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dOddd", &m1,&zs,&NSNII,&NSNIa,&NDYIN))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute SN ejection */
Total_single_mass_ejection(m1,z,NSNII,NSNIa,NDYIN);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
static PyObject *
chemistry_Total_single_mass_ejection_Z(self, args)
PyObject *self;
PyObject *args;
{
double m1;
double NSNII,NSNIa,NDYIN;
PyArrayObject *zs;
PyArrayObject *EMass;
npy_intp ld[1];
int i;
double *z;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dOddd", &m1,&zs,&NSNII,&NSNIa,&NDYIN))
return NULL;
/* create output array */
ld[0]= Cp->nelts+2;
EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* allocate memory for the metallicity array */
z = malloc((Cp->nelts) * sizeof(double));
/* export values */
for (i=0;i<Cp->nelts;i++)
z[i]= *(double *)(zs->data + (i)*(zs->strides[0]));
/* compute SN ejection */
Total_single_mass_ejection_Z(m1,z,NSNII,NSNIa,NDYIN);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i];
free(z);
/* convert in array */
return PyArray_Return(EMass);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_cooling_function(self, args)
PyObject *self;
PyObject *args;
{
/*
on gives :
u_energy
metal = metal(i,2)
parameters
t_const,zmin,zmax,slz,tmin,tmax,slt,FeHSolar,cooling_data_max
*/
PyArrayObject *cooling_data;
double u_energy,metal;
double t_const,zmin,zmax,slz,tmin,tmax,slt,FeHSolar,cooling_data_max;
double cooling,u_cutoff,T,Z;
double rt, rz, ft, fz, v1, v2, v;
int it,iz,itp,izp;
/* parse arguments */
if (!PyArg_ParseTuple(args, "ddOddddddddd", &u_energy, &metal, &cooling_data,&t_const,&zmin,&zmax,&slz,&tmin,&tmax,&slt,&FeHSolar,&cooling_data_max))
return NULL;
u_cutoff=(100)/t_const;
cooling = 0.0;
if (u_energy > u_cutoff)
{
T = log10( t_const*u_energy );
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*(cooling_data( 1, itp)-cooling_data( 1,it) ) + cooling_data( 1,it );
v1 = ft * (*(double *) (cooling_data->data + 1*(cooling_data->strides[0]) + itp*cooling_data->strides[1])
- *(double *) (cooling_data->data + 1*(cooling_data->strides[0]) + it *cooling_data->strides[1]))
+ *(double *) (cooling_data->data + 1*(cooling_data->strides[0]) + it *cooling_data->strides[1]);
//v2 = ft*(cooling_data( 0,itp )-cooling_data( 0, it ) ) + cooling_data( 0, it );
v2 = ft * (*(double *) (cooling_data->data + 0*(cooling_data->strides[0]) + itp*cooling_data->strides[1])
- *(double *) (cooling_data->data + 0*(cooling_data->strides[0]) + it *cooling_data->strides[1]))
+ *(double *) (cooling_data->data + 0*(cooling_data->strides[0]) + it *cooling_data->strides[1]);
v = v2 + fz*(v1-v2);
}
else
{
rt = (T-tmin)/slt;
rz = (Z-zmin)/slz+1.0;
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*(cooling_data( izp, itp)-cooling_data(izp,it)) + cooling_data( izp, it );
v1 = ft * (*(double *) (cooling_data->data + izp*(cooling_data->strides[0]) + itp*cooling_data->strides[1])
- *(double *) (cooling_data->data + izp*(cooling_data->strides[0]) + it *cooling_data->strides[1]))
+ *(double *) (cooling_data->data + izp*(cooling_data->strides[0]) + it *cooling_data->strides[1]);
//v2 = ft*(cooling_data( iz, itp )-cooling_data(iz,it )) + cooling_data( iz, it );
v2 = ft * (*(double *) (cooling_data->data + iz *(cooling_data->strides[0]) + itp*cooling_data->strides[1])
- *(double *) (cooling_data->data + iz *(cooling_data->strides[0]) + it *cooling_data->strides[1]))
+ *(double *) (cooling_data->data + iz *(cooling_data->strides[0]) + it *cooling_data->strides[1]);
v = v2 + fz*(v1-v2);
}
cooling = pow(10,v);
}
return Py_BuildValue("d",cooling);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_get_Mmax(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->Mmax * All.MsoltoCMU);
}
static PyObject *
chemistry_get_Mmin(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->Mmin * All.MsoltoCMU);
}
static PyObject *
chemistry_get_Mco(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->Mco * All.MsoltoCMU);
}
static PyObject *
chemistry_get_SNIa_Mpl(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->SNIa_Mpl * All.MsoltoCMU);
}
static PyObject *
chemistry_get_SNIa_Mpu(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->SNIa_Mpu * All.MsoltoCMU);
}
static PyObject *
chemistry_get_SNII_Mmin(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->SNII_Mmin * All.MsoltoCMU);
}
static PyObject *
chemistry_get_SNII_Mmax(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->SNII_Mmax * All.MsoltoCMU);
}
static PyObject *
chemistry_get_DYIN_Mmin(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->DYIN_Mmin * All.MsoltoCMU);
}
static PyObject *
chemistry_get_DYIN_Mmax(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->DYIN_Mmax * All.MsoltoCMU);
}
static PyObject *
chemistry_get_imf_Ntot(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("d",(double)Cp->imf_Ntot*All.CMUtoMsol); /* in code mass unit */
}
static PyObject *
chemistry_get_as(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *as;
npy_intp ld[1];
int i;
/* create output array */
ld[0]= Cp->n+1;
as = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* import values */
for (i=0;i<Cp->n+1;i++)
*(double *)(as->data + (i)*(as->strides[0])) = Cp->as[i];
return PyArray_Return(as);
}
static PyObject *
chemistry_get_bs(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *bs;
npy_intp ld[1];
int i;
/* create output array */
ld[0]= Cp->n+1;
bs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* import values */
for (i=0;i<Cp->n+1;i++)
*(double *)(bs->data + (i)*(bs->strides[0])) = Cp->bs[i];
return PyArray_Return(bs);
}
static PyObject *
chemistry_get_fs(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *fs;
npy_intp ld[1];
int i;
/* create output array */
ld[0]= Cp->n;
fs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* import values */
for (i=0;i<Cp->n;i++)
*(double *)(fs->data + (i)*(fs->strides[0])) = Cp->fs[i];
return PyArray_Return(fs);
}
static PyObject *
chemistry_get_allnelts(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("i",(int)Cp->nelts+2);
}
static PyObject *
chemistry_get_nelts(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("i",(int)Cp->nelts);
}
static PyObject *
chemistry_get_allelts_labels(self, args)
PyObject *self;
PyObject *args;
{
int i;
PyObject *LabelList,*LabelString;
LabelList = PyList_New((Py_ssize_t)Cp->nelts+2);
for(i=0;i<Cp->nelts+2;i++)
{
LabelString = PyString_FromString(Elt[i].label);
PyList_SetItem(LabelList, (Py_ssize_t)i,LabelString);
}
return Py_BuildValue("O",LabelList);
}
static PyObject *
chemistry_get_elts_labels(self, args)
PyObject *self;
PyObject *args;
{
int i;
PyObject *LabelList,*LabelString;
LabelList = PyList_New((Py_ssize_t)Cp->nelts);
for(i=2;i<Cp->nelts+2;i++)
{
LabelString = PyString_FromString(Elt[i].label);
PyList_SetItem(LabelList, (Py_ssize_t)i-2,LabelString);
}
return Py_BuildValue("O",LabelList);
}
/* static PyObject *
chemistry_get_elts_SolarMassAbundances(self, args)
PyObject *self;
PyObject *args;
{
int i;
npy_intp ld[1];
PyArrayObject *AbList;
ld[0] = Cp->nelts;
AbList = (PyArrayObject *) PyArray_SimpleNew(1,ld,NPY_FLOAT);
for(i=2;i<Cp->nelts+2;i++)
*(float*)(AbList->data + (i-2)*(AbList->strides[0])) = (float) Elt[i].SolarMassAbundance;
return PyArray_Return(AbList);
} */
static PyObject *
chemistry_get_elts_SolarMassAbundances(self, args)
PyObject *self;
PyObject *args;
{
int i;
PyObject *AbDict,*LabelString,*AbVal;
AbDict = PyDict_New();
for(i=2;i<Cp->nelts+2;i++)
{
AbVal = PyFloat_FromDouble(Elt[i].SolarMassAbundance);
LabelString = PyString_FromString(Elt[i].label);
PyDict_SetItem(AbDict,LabelString, AbVal);
}
return Py_BuildValue("O",AbDict);
}
static PyObject *
chemistry_get_MSNIa(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *MSNIa;
npy_intp ld[1];
int i;
/* create output array */
ld[0]= Cp->nelts;
MSNIa = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* import values */
for (i=0;i<Cp->nelts;i++)
*(double *)(MSNIa->data + (i)*(MSNIa->strides[0])) = Elt[i+2].MSNIa*All.MsoltoCMU;
return PyArray_Return(MSNIa);
}
static PyObject *
chemistry_get_MassFracSNII(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *MassFrac;
npy_intp ld[1];
int i;
/* create output array */
ld[0]= Cp->nelts+2;
MassFrac = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(MassFrac->data + (i)*(MassFrac->strides[0])) = MassFracSNII[i];
return PyArray_Return(MassFrac);
}
static PyObject *
chemistry_get_SingleMassFracSNII(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *MassFrac;
npy_intp ld[1];
int i;
/* create output array */
ld[0]= Cp->nelts+2;
MassFrac = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
/* import values */
for (i=0;i<Cp->nelts+2;i++)
*(double *)(MassFrac->data + (i)*(MassFrac->strides[0])) = SingleMassFracSNII[i];
return PyArray_Return(MassFrac);
}
static PyObject *
chemistry_imf_sampling(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *ms;
npy_intp ld[1];
int i;
int n,seed;
/* parse arguments */
if (!PyArg_ParseTuple(args, "ii", &n,&seed))
return NULL;
/* create output array */
ld[0]= n;
ms = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
srandom(seed);
/* import values */
for (i=0;i<n;i++)
*(double *)(ms->data + (i)*(ms->strides[0])) = imf_sampling();
return PyArray_Return(ms);
}
static PyObject *
chemistry_imf_init_seed(self, args)
PyObject *self;
PyObject *args;
{
int n,seed;
/* parse arguments */
if (!PyArg_ParseTuple(args, "i",&seed))
return NULL;
srandom(seed);
return Py_BuildValue("i",1);
}
static PyObject *
chemistry_imf_sampling_single(self, args)
PyObject *self;
PyObject *args;
{
return Py_BuildValue("f",imf_sampling());
}
static PyObject *
chemistry_imf_sampling_single_from_random(self, args)
PyObject *self;
PyObject *args;
{
double f;
if (!PyArg_ParseTuple(args, "f",&f))
return NULL;
return Py_BuildValue("f",imf_sampling_from_random(f));
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_SNII_rate_P(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *ConstSN,*Msn;
double m1,m2,md;
double powSN1,powSN2;
double RSNII;
RSNII = 0.0;
/* parse arguments */
if (!PyArg_ParseTuple(args, "ddddOO", &m1,&m2,&powSN1,&powSN2,&ConstSN,&Msn))
return NULL;
if ( m1 < *(double *) (Msn->data + 2*(Msn->strides[0]) + 1*(Msn->strides[1])) )
md = *(double *) (Msn->data + 2*(Msn->strides[0]) + 1*(Msn->strides[1]));
else
md = m1;
if (md >= m2)
RSNII = 0;
else
RSNII = *(double *) (ConstSN->data + 2*ConstSN->strides[0]) *(pow(m2,powSN1)-pow(md,powSN1));
return Py_BuildValue("d",RSNII);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_SNIa_rate_P(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *ConstSN,*Msn;
double m1,m2,md,mu;
double powSN1,powSN2;
double RSNIa;
double rate;
int i;
/* parse arguments */
if (!PyArg_ParseTuple(args, "ddddOO", &m1,&m2,&powSN1,&powSN2,&ConstSN,&Msn))
return NULL;
RSNIa = 0.0;
for (i=0;i<2;i++)
{
if ( m1 < *(double *) (Msn->data + i*(Msn->strides[0]) + 0*(Msn->strides[1])) )
md = *(double *) (Msn->data + i*(Msn->strides[0]) + 0*(Msn->strides[1]));
else
md = m1;
if ( m2 > *(double *) (Msn->data + i*(Msn->strides[0]) + 1*(Msn->strides[1])) )
mu = *(double *) (Msn->data + i*(Msn->strides[0]) + 1*(Msn->strides[1]));
else
mu = m2;
if (md<mu)
RSNIa = RSNIa+ *(double *) (ConstSN->data + i*ConstSN->strides[0])*(pow(mu,powSN2)-pow(md,powSN2));
}
if ( m1 < *(double *) (Msn->data + 2*(Msn->strides[0]) + 0*(Msn->strides[1])) )
md = *(double *) (Msn->data + 2*(Msn->strides[0]) + 0*(Msn->strides[1]));
else
md = m1;
mu = *(double *) (Msn->data + 2*(Msn->strides[0]) + 1*(Msn->strides[1]));
if (md >= mu)
RSNIa = 0.0;
else
RSNIa = RSNIa*(pow(mu,powSN1)-pow(md,powSN1));
if (RSNIa<0)
RSNIa = 0;
return Py_BuildValue("d",RSNIa);
}
/*********************************/
/* */
/*********************************/
static PyObject *
chemistry_SNII_mass_ejection_P(self, args)
PyObject *self;
PyObject *args;
{
PyArrayObject *ArrayOrigin,*ArrayStep,*ChemArray;
double m1,m2;
int NbElement;
double l1,l2;
int i1,i2,i1p,i2p,j;
double f1,f2;
double v1,v2;
PyArrayObject *ArrMassSNII;
npy_intp ld[1];
/* parse arguments */
if (!PyArg_ParseTuple(args, "ddOOOi", &m1,&m2,&ArrayOrigin,&ArrayStep,&ChemArray,&NbElement))
return NULL;
/* create output array */
ld[0]= NbElement+2;
ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
l1 = ( log10(m1) - *(double *)(ArrayOrigin->data + 0*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 0*(ArrayStep->strides[0])) ;
l2 = ( log10(m2) - *(double *)(ArrayOrigin->data + 0*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 0*(ArrayStep->strides[0])) ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<1) i1=1;
if (i2<1) i2=1;
/* --------- TOTAL GAS ---------- */
j = NbElement;
v1=f1* (*(double *)(ChemArray->data + (i1p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))
- *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])))
+ *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]));
v2=f2* (*(double *)(ChemArray->data + (i2p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))
- *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])))
+ *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]));
*(double *)(ArrMassSNII->data + (j)*(ArrMassSNII->strides[0])) = v2-v1;
/* --------- He core therm ---------- */
j = NbElement+1;
v1=f1* (*(double *)(ChemArray->data + (i1p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))
- *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])))
+ *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]));
v2=f2* (*(double *)(ChemArray->data + (i2p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))
- *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])))
+ *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]));
*(double *)(ArrMassSNII->data + (j)*(ArrMassSNII->strides[0])) = v2-v1;
/* --------- Metals ---------- */
l1 = ( log10(m1) - *(double *)(ArrayOrigin->data + 1*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 1*(ArrayStep->strides[0])) ;
l2 = ( log10(m2) - *(double *)(ArrayOrigin->data + 1*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 1*(ArrayStep->strides[0])) ;
if (l1 < 0.0) l1 = 0.0;
if (l2 < 0.0) l2 = 0.0;
i1 = (int)l1;
i2 = (int)l2;
i1p = i1 + 1;
i2p = i2 + 1;
f1 = l1 - i1;
f2 = l2 - i2;
/* check (yr) */
if (i1<1) i1=1;
if (i2<1) i2=1;
for (j=0;j<NbElement;j++)
{
v1=f1* (*(double *)(ChemArray->data + (i1p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))
- *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])))
+ *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]));
v2=f2* (*(double *)(ChemArray->data + (i2p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))
- *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])))
+ *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]));
*(double *)(ArrMassSNII->data + (j)*(ArrMassSNII->strides[0])) = v2-v1;
}
return PyArray_Return(ArrMassSNII);
}
#ifdef CHIMIE_OPTIMAL_SAMPLING
static PyObject *
chemistry_optimal_init_norm(self, args)
PyObject *self;
PyObject *args;
{
double m_max,Mecl;
if (!PyArg_ParseTuple(args, "d", &Mecl))
return NULL;
m_max = optimal_init_norm(Mecl);
return Py_BuildValue("d",m_max);
}
static PyObject *
chemistry_optimal_get_next_mass(self, args)
PyObject *self;
PyObject *args;
{
double m;
double m_next;
/* parse arguments */
if (!PyArg_ParseTuple(args, "d", &m))
return NULL;
m_next = optimal_get_next_mass(m);
return Py_BuildValue("d",m_next);
}
static PyObject *
chemistry_stop_loop(self, args)
PyObject *self;
PyObject *args;
{
double m;
int bool;
/* parse arguments */
if (!PyArg_ParseTuple(args, "d", &m))
return NULL;
bool = optimal_stop_loop(m);
return Py_BuildValue("i",bool);
}
static PyObject *
chemistry_optimal_get_m1_from_m2(self, args)
PyObject *self;
PyObject *args;
{
double m1,m2,mp;
/* parse arguments */
if (!PyArg_ParseTuple(args, "dd", &m2,&mp))
return NULL;
m1 = optimal_get_m1_from_m2(m2,mp);
return Py_BuildValue("d",m1);
}
#endif
/* definition of the method table */
static PyMethodDef chemistryMethods[] = {
{"CodeUnits_to_SolarMass_Factor", chemistry_CodeUnits_to_SolarMass_Factor, METH_VARARGS,
"convertion factor : CodeUnits -> SolarMass"},
{"SolarMass_to_CodeUnits_Factor", chemistry_SolarMass_to_CodeUnits_Factor, METH_VARARGS,
"convertion factor : SolarMass -> CodeUnits"},
{"SetVerbosityOn", (PyCFunction)chemistry_SetVerbosityOn, METH_VARARGS,
"Set verbosity to on"},
{"SetVerbosityOff", (PyCFunction)chemistry_SetVerbosityOff, METH_VARARGS,
"Set verbosity to off"},
{"InitDefaultParameters", (PyCFunction)chemistry_InitDefaultParameters, METH_VARARGS,
"Init default parameters"},
{"SetParameters", (PyCFunction)chemistry_SetParameters, METH_VARARGS,
"Set gadget parameters"},
{"GetParameters", (PyCFunction)chemistry_GetParameters, METH_VARARGS,
"get some gadget parameters"},
{"init_chimie", chemistry_init_chimie, METH_VARARGS| METH_KEYWORDS,
"Init chimie."},
{"info", chemistry_info, METH_VARARGS| METH_KEYWORDS,
"Get info on tables."},
{"set_table", chemistry_set_table, METH_VARARGS,
"Set the chimie table."},
{"get_imf", chemistry_get_imf, METH_VARARGS,
"Compute corresponding imf value."},
{"get_imf_M", chemistry_get_imf_M, METH_VARARGS,
"Compute the mass fraction between m1 and m2."},
{"get_imf_N", chemistry_get_imf_N, METH_VARARGS,
"Compute the fraction number between m1 and m2."},
{"star_lifetime", chemistry_star_lifetime, METH_VARARGS,
"Compute star life time."},
{"star_mass_from_age", chemistry_star_mass_from_age, METH_VARARGS,
"Return the stellar mass that has a lifetime equal to t."},
{"DYIN_rate", chemistry_DYIN_rate, METH_VARARGS,
"Return the number of dying stars per unit mass with masses between m1 and m2."},
{"SNII_rate", chemistry_SNII_rate, METH_VARARGS,
"Return the number of SNII per unit mass with masses between m1 and m2."},
{"SNIa_rate", chemistry_SNIa_rate, METH_VARARGS,
"Return the number of SNIa per unit mass with masses between m1 and m2."},
{"SNIa_single_rate", chemistry_SNIa_single_rate, METH_VARARGS,
"Return the number of SNIa per unit mass for a star of mass m."},
{"DYIN_mass_ejection", chemistry_DYIN_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of dying stars with masses between m1 and m2."},
{"DYIN_mass_ejection_Z", chemistry_DYIN_mass_ejection_Z, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of dying stars with masses between m1 and m2 (metallicity dependency)."},
{"DYIN_single_mass_ejection", chemistry_DYIN_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements due to the explotion of one dying star of mass m."},
{"DYIN_single_mass_ejection_Z", chemistry_DYIN_single_mass_ejection_Z, METH_VARARGS,
"Mass fraction of ejected elements due to the explotion of one dying star of mass m (metallicity dependency)."},
{"SNII_mass_ejection", chemistry_SNII_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of SNII with masses between m1 and m2."},
//{"SNII_mass_ejection_Z", chemistry_SNII_mass_ejection_Z, METH_VARARGS,
// "Mass fraction of ejected elements, per unit mass due to the explotion of SNII with masses between m1 and m2 (metallicity dependency)."},
{"SNII_single_mass_ejection", chemistry_SNII_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements due to the explotion of one SNII of mass m."},
//{"SNII_single_mass_ejection_Z", chemistry_SNII_single_mass_ejection_Z, METH_VARARGS,
// "Mass fraction of ejected elements due to the explotion of one SNII star of mass m (metallicity dependency)."},
{"SNIa_mass_ejection", chemistry_SNIa_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of SNIa with masses between m1 and m2."},
{"SNIa_single_mass_ejection", chemistry_SNIa_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements due to the explotion of one SNIa of mass m."},
{"Total_mass_ejection", chemistry_Total_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of SNIa and SNII with masses between m1 and m2."},
{"Total_mass_ejection_Z", chemistry_Total_mass_ejection_Z, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of SNIa and SNII with masses between m1 and m2 (metallicity dependency)."},
{"DYIN_Total_single_mass_ejection", chemistry_DYIN_Total_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements (including processed and non processed elements) due to the explotion of one dying star of mass m."},
{"SNII_Total_single_mass_ejection", chemistry_SNII_Total_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements (including processed and non processed elements) due to the explotion of one SNII of mass m."},
{"SNIa_Total_single_mass_ejection", chemistry_SNIa_Total_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements (including processed and non processed elements) due to the explotion of one SNIa of mass m."},
{"Total_single_mass_ejection", chemistry_Total_single_mass_ejection, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of one dying star of mass m1."},
{"Total_single_mass_ejection_Z", chemistry_Total_single_mass_ejection_Z, METH_VARARGS,
"Mass fraction of ejected elements, per unit mass due to the explotion of one dying star of mass m1 (metallicity dependency)."},
{"get_Mmax", chemistry_get_Mmax, METH_VARARGS,
"Get max star mass of the IMF, in code unit."},
{"get_Mmin", chemistry_get_Mmin, METH_VARARGS,
"Get min star mass of the IMF, in code unit."},
{"get_Mco", chemistry_get_Mco, METH_VARARGS,
"Get mean WD mass, in code unit."},
{"get_SNIa_Mpl", chemistry_get_SNIa_Mpl, METH_VARARGS,
"Get min mass of SNIa, in code unit."},
{"get_SNIa_Mpu", chemistry_get_SNIa_Mpu, METH_VARARGS,
"Get max mass of SNIa, in code unit."},
{"get_SNII_Mmin", chemistry_get_SNII_Mmin, METH_VARARGS,
"Get min mass of SNII, in code unit."},
{"get_SNII_Mmax", chemistry_get_SNII_Mmax, METH_VARARGS,
"Get max mass of SNII, in code unit."},
{"get_DYIN_Mmin", chemistry_get_DYIN_Mmin, METH_VARARGS,
"Get min mass of DYIN, in code unit."},
{"get_DYIN_Mmax", chemistry_get_DYIN_Mmax, METH_VARARGS,
"Get max mass of DYIN, in code unit."},
{"get_imf_Ntot", chemistry_get_imf_Ntot, METH_VARARGS,
"Get number of stars in the imf, per unit mass."},
{"get_as", chemistry_get_as, METH_VARARGS,
"Get power coefficients."},
{"get_bs", chemistry_get_bs, METH_VARARGS,
"Get normalisation coefficients."},
{"get_fs", chemistry_get_fs, METH_VARARGS,
"Get fs, mass fraction at ms."},
{"get_allnelts", chemistry_get_allnelts, METH_VARARGS,
"Get the number of element considered, including ejected mass (Ej) and non processed ejected mass (Ejnp).."},
{"get_nelts", chemistry_get_nelts, METH_VARARGS,
"Get the number of element considered."},
{"get_allelts_labels", chemistry_get_allelts_labels, METH_VARARGS,
"Get the labels of elements, including ejected mass (Ej) and non processed ejected mass (Ejnp)."},
{"get_elts_labels", chemistry_get_elts_labels, METH_VARARGS,
"Get the labels of elements."},
{"get_elts_SolarMassAbundances", chemistry_get_elts_SolarMassAbundances, METH_VARARGS,
"Get the solar mass abundance of elements."},
{"get_MassFracSNII", chemistry_get_MassFracSNII, METH_VARARGS,
"Get the mass fraction per element ejected by a set of SNII."},
{"get_SingleMassFracSNII", chemistry_get_SingleMassFracSNII, METH_VARARGS,
"Get the mass fraction per element ejected by a SNII."},
{"get_MSNIa", chemistry_get_MSNIa, METH_VARARGS,
"Get the mass per element ejected by a SNIa."},
{"cooling_function", chemistry_cooling_function, METH_VARARGS,
"Compute cooling."},
{"imf_sampling", chemistry_imf_sampling, METH_VARARGS,
"Sample imf with n points."},
{"imf_sampling_single", chemistry_imf_sampling_single, METH_VARARGS,
"Sample imf with a single point."},
{"imf_init_seed", chemistry_imf_init_seed, METH_VARARGS,
"Init the random seed."},
{"imf_sampling_single_from_random", chemistry_imf_sampling_single_from_random, METH_VARARGS,
"Sample imf with a single point from a given random number."},
/* old poirier */
{"SNIa_rate_P", chemistry_SNIa_rate_P, METH_VARARGS,
"Return the number of SNIa per unit mass and time. (Poirier version)"},
{"SNII_rate_P", chemistry_SNII_rate_P, METH_VARARGS,
"Return the number of SNII per unit mass and time. (Poirier version)"},
{"SNII_mass_ejection_P", chemistry_SNII_mass_ejection_P, METH_VARARGS,
"Mass ejection due to SNII per unit mass and time. (Poirier version)"},
#ifdef CHIMIE_OPTIMAL_SAMPLING
{"optimal_sampling_init_norm", chemistry_optimal_init_norm, METH_VARARGS,
"init normalistation"},
{"optimal_sampling_get_next_mass", chemistry_optimal_get_next_mass, METH_VARARGS,
"next star mass"},
{"optimal_sampling_stop_loop", chemistry_stop_loop, METH_VARARGS,
"return 1 if the loop for optimal sampling must be stopped "},
{"optimal_sampling_get_m1_from_m2", chemistry_optimal_get_m1_from_m2, METH_VARARGS,
"for a givent mass m2, return m1, such that the mass in the IMF between m1 and m2 is mp."},
#endif
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initchemistry(void)
{
(void) Py_InitModule("chemistry", chemistryMethods);
import_array();
}
#endif /* PYCHEM */
#endif /* CHIMIE */
c4science · Help