diff --git a/src/chimie.c b/src/chimie.c index 4e28b86..c3bf1d6 100644 --- a/src/chimie.c +++ b/src/chimie.c @@ -1,8339 +1,8358 @@ #include #include #include #include #include #include #include "allvars.h" #include "proto.h" #ifdef CHIMIE #ifdef CHIMIE_FROM_HDF5 #include #include "hdf5io.h" #endif #ifdef PYCHEM #include #include #include #include #include /* **************************************************** 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 #define KPC_IN_CM 3.085e+21 /******************************************************** 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; iversion = 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;icoeff_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; iMmin = 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; iMpl = 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; iMmin = 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;inelts;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; iMmin = 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;inelts;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;inelts;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=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 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;i0) for (i=0;i 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;in; /* convert m in msol */ m = m*All.CMUtoMsol; if (n==0) return Cp->bs[0]* pow(m,Cp->as[0]); else { for (i=0;ims[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 (m1ms[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;ims[i ],m1); mmax = dmin(Cp->ms[i+1],m2); if (mminas[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 (m1ms[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;ims[i ],m1); mmax = dmin(Cp->ms[i+1],m2); if (mminas[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 (ffs[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;ifs[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 (ffs[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;ifs[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;ims[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;ibs[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;ibs[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;ims[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 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;nfcoeff_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;in+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;in+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;inelts+2;i++) // Elt[i].MminSNII = log10(Elt[i].MminSNII); // //for (i=0;inelts+2;i++) // Elt[i].MminDYIN = log10(Elt[i].MminDYIN); /* output info */ //if (verbose && ThisTask==0) // { // //printf("-- SNII_cte -- \n"); // //for (i=0;in+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;in+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); } 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); } 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"); } /*! This function performe simple checks * to validate the chemistry initialization */ void check_chimie(void) { int i; if (ThisTask==0) { printf("(Taks=%d) Number of elts : %d\n",ThisTask,Cp->nelts); for(i=2;inelts+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 %g %g\n",Elts[it][i].MminSNII,Elts[it][i].StepSNII); for (j=0;j %g %g\n",Elts[it][i].MminDYIN,Elts[it][i].StepDYIN); for (j=0;jnelts; } /*! 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;inelts;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; /* 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 (mdSNIa_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 (mdSNIa_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; /* 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;jnelts+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;jnelts+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;knelts+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;knelts+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;jnelts+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;jnelts+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;knelts+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;knelts+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;jnelts+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;jnelts+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;knelts+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;knelts+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;jnelts+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;jnelts+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;knelts+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;knelts+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;jnelts+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; /* total ejected gas mass */ SingleMassFracSNIa[0] = Cp->Mco * All.MsoltoCMU; /* ejected mass per element */ for (j=2;jnelts+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 */ DYIN_mass_ejection(m1,m2); /* total ejected gas mass */ EjectedMass[0] = M0 * ( MassFracDYIN[0] + MassFracSNII[0] + MassFracSNIa[0] ); /* ejected mass per element */ for (j=2;jnelts+2;j++) EjectedMass[j] = M0*( MassFracDYIN[j] +z[j-2]*MassFracDYIN[1] + MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] ); /* 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 */ DYIN_mass_ejection_Z(m1,m2,Z); /* total ejected gas mass */ EjectedMass[0] = M0 * ( MassFracDYIN[0] + MassFracSNII[0] + MassFracSNIa[0] ); /* ejected mass per element */ for (j=2;jnelts+2;j++) EjectedMass[j] = M0*( MassFracDYIN[j] +z[j-2]*MassFracDYIN[1] + MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] ); /* 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;jnelts+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;jnelts+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;jnelts+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 */ DYIN_single_mass_ejection(m1); /* total ejected gas mass */ SingleEjectedMass[0] = M0 * ( SingleMassFracDYIN[0]*NDYIN + SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleEjectedMass[j] = M0*( SingleMassFracDYIN[j]*NDYIN +z[j-2]*SingleMassFracDYIN[1]*NDYIN + SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa; /* 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 */ DYIN_single_mass_ejection_Z(m1,Z); /* total ejected gas mass */ SingleEjectedMass[0] = M0 * ( SingleMassFracDYIN[0]*NDYIN + SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleEjectedMass[j] = M0*( SingleMassFracDYIN[j]*NDYIN +z[j-2]*SingleMassFracDYIN[1]*NDYIN + SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa; /* 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].Density*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].Density*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 */ 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 */ if (ThisTask==0) printf("Chimie computation done.\n"); t1 = second(); All.CPU_Chimie += timediff(t0, t1); } /*! 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++) { if(P[n].Ti_endstep == All.Ti_Current) 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 */ if((P[i].Ti_endstep == All.Ti_Current)&&(P[i].Type == ST)) { 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 */ Ti1 = P[i].Ti_begstep; /* 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 { t1 = All.TimeBegin + (P[i].Ti_begstep * All.Timebase_interval); 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) { 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;k0) { current_mass = StP[m].MassMax; printf("(%d) FOF_SFR : INDIVIDUAL STARS id=%8d %8.4g %8.4g %8d current_mass=%g Mass=%g MassSSP\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); 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(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); + 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); if ( (Cp->SNII_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->SNII_Mmax*All.MsoltoCMU) ) { NSNII++; - StP[m].NStars--; + //StP[m].NStars--; printf("(%d) FOF_SFR : id=%8d one SNII\n",ThisTask,P[i].ID); SNII_Total_single_mass_ejection(current_mass,metals); StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */ for (k=0;kDYIN_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->DYIN_Mmax*All.MsoltoCMU) ) { NDYIN++; - StP[m].NStars--; + //StP[m].NStars--; printf("(%d) FOF_SFR : id=%8d one DYIN\n",ThisTask,P[i].ID); DYIN_Total_single_mass_ejection(current_mass,metals); StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */ for (k=0;k 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;kDYIN_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;k0) 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 */ if (P[i].Mass-StP[m].TotalEjectedGasMass<0) { printf("mass wants to be less than zero...\n"); printf("ID=%d Mass=%g StP[m].TotalEjectedGasMass=%g\n",P[i].ID,P[i].Mass,StP[m].TotalEjectedGasMass); endrun(777100); } P[i].Mass = P[i].Mass-StP[m].TotalEjectedGasMass; if(P[i].Mass<0) endrun(777023); } /******************************************/ /* 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].Volume = StP[m].Volume; #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 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= (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 volume; #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(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; volume = StP[target_stp].Volume; #ifdef CHIMIE_KINETIC_FEEDBACK ngbmass = StP[target_stp].NgbMass; #endif ejectedGasMass = StP[target_stp].TotalEjectedGasMass; for(k=0;k 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; /* metal injection */ for(k=0;k= 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,¶msDict)) { 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;in;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;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = MassFracDYIN[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = MassFracDYIN[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = SingleMassFracDYIN[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = SingleMassFracDYIN[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = MassFracSNII[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = MassFracSNII[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = SingleMassFracSNII[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = SingleMassFracSNII[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassSNIa->data + (i)*(ArrMassSNIa->strides[0])) = MassFracSNIa[i]; /* convert in array */ return Py_BuildValue("O",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;inelts+2;i++) *(double *)(ArrMassSNIa->data + (i)*(ArrMassSNIa->strides[0])) = SingleMassFracSNIa[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = EjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = EjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;inelts;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;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; /* convert in array */ return Py_BuildValue("O",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;in+1;i++) *(double *)(as->data + (i)*(as->strides[0])) = Cp->as[i]; return Py_BuildValue("O",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;in+1;i++) *(double *)(bs->data + (i)*(bs->strides[0])) = Cp->bs[i]; return Py_BuildValue("O",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;in;i++) *(double *)(fs->data + (i)*(fs->strides[0])) = Cp->fs[i]; return Py_BuildValue("O",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;inelts+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;inelts+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;inelts+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;inelts+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;inelts;i++) *(double *)(MSNIa->data + (i)*(MSNIa->strides[0])) = Elt[i+2].MSNIa*All.MsoltoCMU; return Py_BuildValue("O",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;inelts+2;i++) *(double *)(MassFrac->data + (i)*(MassFrac->strides[0])) = MassFracSNII[i]; return Py_BuildValue("O",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;inelts+2;i++) *(double *)(MassFrac->data + (i)*(MassFrac->strides[0])) = SingleMassFracSNII[i]; return Py_BuildValue("O",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;idata + (i)*(ms->strides[0])) = imf_sampling(); return Py_BuildValue("O",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 (mddata + 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;jdata + (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 Py_BuildValue("O",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 */