diff --git a/cmd.c b/cmd.c deleted file mode 100644 index 8d99a68..0000000 --- a/cmd.c +++ /dev/null @@ -1,76 +0,0 @@ - -int cmd(int kstar, double rstar, double lstar, double Rgal, double *abvmag, double *vmag, double *BV, double *Teff, double *dvmag, double *dBV) - { - - double lTeff, BC, kb; - double bvc[8], bcc[8]; - double dbmag; - double BCsun, abvmagsun; - double rand1, rand2, prand; - - kb = 5.6704E-08*0.5*1.3914E9*0.5*1.3914E9/3.846E26; //Stefan-Boltzmann constant in Lsun Rsun^-2 K^-4 - - bvc[0] = -654597.405559323; - bvc[1] = 1099118.61158915; - bvc[2] = -789665.995692672; - bvc[3] = 314714.220932623; - bvc[4] = -75148.4728506455; - bvc[5] = 10751.803394526; - bvc[6] = -853.487897283685; - bvc[7] = 28.9988730655392; - - bcc[0] = -4222907.80590972; - bcc[1] = 7209333.13326442; - bcc[2] = -5267167.04593882; - bcc[3] = 2134724.55938336; - bcc[4] = -518317.954642773; - bcc[5] = 75392.2372207101; - bcc[6] = -6082.7301194776; - bcc[7] = 209.990478646363; - - BCsun = 0.11; //sun's bolometric correction - abvmagsun = 4.83; //sun's absolute V magnitude - - if( rstar && (kstar<14) ) - { - *Teff = pow(lstar/(4.0*Pi*rstar*rstar*kb),0.25); - - if( (*Teff>3000.0) && (*Teff<55000.0) ) - { - lTeff = log10(*Teff); - *BV = bvc[0] + bvc[1]*lTeff + bvc[2]*pow(lTeff,2) + bvc[3]*pow(lTeff,3) + bvc[4]*pow(lTeff,4) + bvc[5]*pow(lTeff,5) + bvc[6]*pow(lTeff,6) + bvc[7]*pow(lTeff,7); - BC = bcc[0] + bcc[1]*lTeff + bcc[2]*pow(lTeff,2) + bcc[3]*pow(lTeff,3) + bcc[4]*pow(lTeff,4) + bcc[5]*pow(lTeff,5) + bcc[6]*pow(lTeff,6) + bcc[7]*pow(lTeff,7); - if(lstar) *abvmag = -2.5*log10(lstar)-BC+BCsun+abvmagsun; - *vmag = *abvmag + 5.0*log10(Rgal) - 5.0; - - do{ - rand1 = 2.0*drand48()-1.0; - rand2 = 2.0*drand48()-1.0; - } while (rand1*rand1+rand2*rand2 > 1.0); - - prand = sqrt(-2.0*log(rand1*rand1+rand2*rand2)/(rand1*rand1+rand2*rand2)); - *dvmag = rand1*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2)); - dbmag = rand2*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2)); - *dBV = *dvmag + dbmag; - } - else - { - *vmag = 9999.9; - *abvmag = 9999.9; - *BV = 9999.9; - *dvmag = 0.0; - *dBV = 0.0; - } - } - else - { - *Teff = 0.0; - *vmag = 9999.9; - *abvmag = 9999.9; - *BV = 9999.9; - *dvmag = 0.0; - *dBV = 0.0; - } - - return(0); -} diff --git a/def_DEN.c b/def_DEN.c deleted file mode 100644 index 7cb217b..0000000 --- a/def_DEN.c +++ /dev/null @@ -1,38 +0,0 @@ - -/*************************************************************************/ -/*************************************************************************/ -/*************************************************************************/ - -void DEF_DEN(double x_gas[][3], double m_gas[], double h_gas[], int sosed_gas[][NB], double DEN_gas[]) -{ - -double Xij, Yij, Zij, Rij, tmp; - -/* SPH style */ - -for(i=0;i= 0.0) && (v <= 1.0) ) - { - tmp = 1.0 - 1.5*v2 + 0.75*v3; - } - else - { - if( (v > 1.0) && (v < 2.0) ) - { - tmp = 0.25*(2.0-v)*(2.0-v)*(2.0-v); - } - else - { - tmp = 0.0; /* if ( v >= 2.0 ) */ - } - } - - tmp = 1.0/(Pi*h*h*h)*tmp; - - return(tmp); - -} - -/*************************************************************************/ - -void DEF_H(double x_gas[][3], double h_gas[], int sosed_gas[][NB]) -{ - -int sosed_tmp[N_GAS_MAX]; -double Xij, Yij, Zij; - -for(i=0;i N_GAS_MAX vsegda ! -double x_star, y_star, z_star, DEN_star; -double h_min = 1.0E-04, h_max = 1.0; - -#include "def_H.c" -#include "def_DEN.c" -#include "def_DEN_STAR.c" - -#endif // SPH_GAS - -#include "drag_force.c" - -// calculate the SG for the set of active particles which is deepley INSIDE the EXT BH infl. rad. -// EXPERIMENTAL feature !!! - -//int SG_CALC = 0; -//double summ_tmp_r = 0.0, aver_tmp_r = 0.0, R_INT_CUT = 1.0E-03; - -#endif // STARDISK - - - -/* GMC data... */ - -#ifdef GMC - -double dt_gmc, E_pot_gmc, E_pot_gmc_gmc, E_pot_ext_gmc, E_kin_gmc; - -#define N_gmc_MAX (500) - -int N_gmc, ind_gmc[N_gmc_MAX], name_gmc[N_gmc_MAX]; - -double m_gmc[N_gmc_MAX], pot_gmc_gmc[N_gmc_MAX], pot_ext_gmc[N_gmc_MAX], - x_gmc[N_gmc_MAX][3], v_gmc[N_gmc_MAX][3], a_gmc[N_gmc_MAX][3], - eps_gmc[N_gmc_MAX]; - -double pot_gmc[N_MAX]; - -#endif // GMC - - /****************************************************************************/ /****************************************************************************/ @@ -788,349 +603,18 @@ gettimeofday(&tv, NULL); } /****************************************************************************/ -#ifdef STEVOL -/****************************************************************************/ -double t_dead(double m, double Z) -/* -m - mass of the star (in Msol) -Z - metallicity (absolut values - i.e. Zsol = 0.0122) -t - return value of star lifetime (in Year) - -Raiteri, Villata & Navarro, 1996, A&A, 315, 105 -*/ -{ -double a0, a1, a2, lZ, lm, tmp=0.0; - -m *= (m_norm/Msol); - -lZ = log10(Z); -lm = log10(m); - -a0 = 10.130 + 0.07547*lZ - 0.008084*lZ*lZ; -a1 = -4.424 - 0.79390*lZ - 0.118700*lZ*lZ; -a2 = 1.262 + 0.33850*lZ + 0.054170*lZ*lZ; - -tmp = a0 + a1*lm + a2*lm*lm; - -tmp = pow(10.0,tmp); - -tmp *= (Year/t_norm); - -return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_loss_old(double m, double Z) -/* -m - initial mass of the star (in Msol) -Z - metallicity (e.g. Zsol = 0.02) -m_ret - returned maass to ISM from the star (in Msol) - -approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305 -*/ -{ -double tmp=0.0; - -m *= (m_norm/Msol); - -// tmp = (-0.46168 + 0.912274 * m); // first approx. - -tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m; - -tmp *= (Msol/m_norm); - -return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_curr_old(double m0, double Z, double t0, double t) -{ -double td, ml, tmp=0.0; - -td = t_dead(m0, Z); -ml = m_loss_old(m0, Z); - -if ( (t-t0) < td ) - { - tmp = m0; // instant mass loss -// tmp = m0 - (t-t0)/(td-t0) * ml; // linear mass loss - } -else - { - tmp = m0 - ml; - - if (event[i] == 0) - { - num_dead++; - event[i] = 1; - } - - } /* if ( (t-t0) < td ) */ - -return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_loss(double m, double Z) -/* -m - initial mass of the star (in norm. units) -Z - metallicity (absolut value - i.e. Zsol = 0.0122) -m_loss - returned maass to ISM from the star (in norm. units) -*/ -{ -double tmp=0.0; - -m *= (m_norm/Msol); - -if ( m < 8.0) - { -// approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305 - -// tmp = (-0.46168 + 0.912274 * m); // first approx. - tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m; // second approx. - } -else - { -// approx. from data of Woosley & Weaver, 1995, ApJS, 110, 181 - -// tmp = 0.813956*pow(m, 1.04214); - tmp = 0.813956*pow(m, 1.04); - } - -tmp *= (Msol/m_norm); - -return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_curr(double m0, double Z, double t0, double t) -{ -double td, ml, tmp=0.0; - -tmp = m0; - -if ( event[i] == 0 ) - { - - td = t_dead(m0, Z); - - if ( (t-t0) > td ) - { - - ml = m_loss(m0, Z); - - tmp = m0 - ml; - - num_dead++; - -// wd http://en.wikipedia.org/wiki/Chandrasekhar_limit 1.38 - - if ( tmp * (m_norm/Msol) < 1.38 ) - { - event[i] = 1; - } - else - { - -// ns http://en.wikipedia.org/wiki/Neutron_stars 1.38 - ~2.0 ??? - - if ( tmp * (m_norm/Msol) > 1.38 ) event[i] = 2; - -// bh http://en.wikipedia.org/wiki/Tolman-Oppenheimer-Volkoff_limit ~1.5 - ~3.0 ??? - - if ( tmp * (m_norm/Msol) > 2.00 ) event[i] = 3; - - } - - } /* if ( (t-t0) > td ) */ - - } /* if ( event[i] == 0 ) */ - -return(tmp); -} -/****************************************************************************/ -#endif - - -#ifdef STEVOL_SSE -/****************************************************************************/ -double m_curr(int ipart, double m0, double Z, double t0, double t) -{ -//set up parameters and variables for SSE (Hurley, Pols & Tout 2002) -int kw=0; //stellar type -double mass; //initial mass -double mt=0.0; //actual mass -double r=0.0; //radius -double lum=0.0; //luminosity -double mc=0.0; //core mass -double rc=0.0; //core radius -double menv=0.0; //envelope mass -double renv=0.0; //envelope radius -double ospin=0.0; //spin -double epoch=0.0; //time spent in current evolutionary state -double tms=0.0; //main-sequence lifetime -double tphys; //initial age -double tphysf; //final age -double dtp=0.0; //data store value, if dtp>tphys no data will be stored -double z; //metallicity -double zpars[20]; //metallicity parameters -double vkick=0.0; //kick velocity for compact remnants -double vs[3]; //kick velocity componets - -// M_V_[mag] V_[mag] B-V_[mag] T_eff_[K] dV_[mag] d(B-V) -double abvmag, vmag, BV, Teff, dvmag, dBV; - -mass = m0*(m_norm/Msol); // m0[i] initial mass of the stars Msol -mt = mass; // current mass for output in Msol -tphys = t0*(t_norm/Myr); // t0[i] time of the star formation Myr -tphysf = t*(t_norm/Myr); // t[i] current time for the star Myr -dtp = 0.0; // output parameter, just set to 0.0 !!! -z = Z; // ZZZ[i] of the star - -/* -if (ipart > 9995) - { -printf("Initial parameters: \n"); -printf("M = %g [Mo] \n", mass); -printf("Z = %g \n", z); -printf("kw = %d \n", kw); -printf("spin = %g \n", ospin); -printf("epoch = %g [Myr] \n", epoch); -printf("t_beg = %g [Myr] \n", tphys); -printf("t_end = %g [Myr] \n", tphysf); -printf("R = %g [R_o] \n", r); -printf("L = %g [L_o] \n", lum); -printf("V_kick = %g [km/s] \n", vkick); - } -*/ - -for (k=0; k<20; k++) zpars[k] = 0; -zcnsts_(&z,zpars); //get metallicity parameters - -evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, &epoch, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick, vs); -cmd(kw, r, lum, Rgal, &abvmag, &vmag, &BV, &Teff, &dvmag, &dBV); - -/* -if (ipart > 9995) - { -printf("Final parameters: \n"); -printf("M = %g [Mo] \n", mt); -printf("Z = %g \n", z); -printf("kw = %d \n", kw); -printf("spin = %g \n", ospin); -printf("epoch = %g [Myr] \n", epoch); -printf("t_beg = %g [Myr] \n", tphys); -printf("t_end = %g [Myr] \n", tphysf); -printf("R = %g [R_o] \n", r); -printf("L = %g [L_o] \n", lum); -printf("V_kick = %g [km/s] \n", vkick); -printf("M_V_[mag]\tV_[mag]\t\tB-V_[mag]\tT_eff_[K]\tdV_[mag]\td(B-V)\n"); -printf("%.3E\t%.3E\t%.3E\t%.3E\t%.3E\t%.3E\n", abvmag, vmag, BV, Teff, dvmag, dBV); -// printf("%g\t%g\t%g\t%g\t%g\t%g\n", abvmag, vmag, BV, Teff, dvmag, dBV); - } -*/ - -event[ipart] = kw; // Evolution type (0 - 15). -SSEMass[ipart] = mt; // Mass in Msol -SSESpin[ipart] = ospin; // Spin in SI ??? [kg*m*m/s] -SSERad[ipart] = r; // Radius in Rsol -SSELum[ipart] = lum; // Luminosity in Lsol -SSETemp[ipart] = Teff; // Temperature - -SSEMV[ipart] = abvmag; // M_V absolute V magnitude -SSEBV[ipart] = BV; // B-V color index - -/* -c STELLAR TYPES - KW -c -c 0 - deeply or fully convective low mass MS star -c 1 - Main Sequence star -c 2 - Hertzsprung Gap -c 3 - First Giant Branch -c 4 - Core Helium Burning -c 5 - First Asymptotic Giant Branch -c 6 - Second Asymptotic Giant Branch -c 7 - Main Sequence Naked Helium star -c 8 - Hertzsprung Gap Naked Helium star -c 9 - Giant Branch Naked Helium star -c 10 - Helium White Dwarf -c 11 - Carbon/Oxygen White Dwarf -c 12 - Oxygen/Neon White Dwarf -c 13 - Neutron Star -c 14 - Black Hole -c 15 - Massless Supernova -*/ - -if ( van_kick == 1 ) // we have a kick ??? - { - -if ( (event_old[ipart] < 10) && ( (event[ipart] == 13) || (event[ipart] == 14) ) ) // NS or BH - { - - if (myRank == rootRank) - { -// printf("KICK: %06d %.6E [Myr] %02d (%02d) %.6E [Mo] %.6E [km/s] [% .3E, % .3E, % .3E] OLD: [% .3E, % .3E, % .3E] %.6E [Mo] \n", ipart, tphysf, kw, event_old[ipart], mt, vkick, vs[0], vs[1], vs[2], v[ipart][0]*(v_norm/km), v[ipart][1]*(v_norm/km), v[ipart][2]*(v_norm/km), mass); - printf("KICK: %06d %.6E %02d %02d %.6E %.6E % .3E % .3E % .3E % .3E % .3E % .3E %.6E \n", ipart, tphysf, kw, event_old[ipart], mt, vkick, vs[0], vs[1], vs[2], v[ipart][0]*(v_norm/km), v[ipart][1]*(v_norm/km), v[ipart][2]*(v_norm/km), mass); - fflush(stdout); - } /* if (myRank == rootRank) */ - - v[ipart][0] += vs[0]*(km/v_norm); - v[ipart][1] += vs[1]*(km/v_norm); - v[ipart][2] += vs[2]*(km/v_norm); - } - - } // we have a kick ??? - - -event_old[ipart] = event[ipart]; - - -tmp = mt*(Msol/m_norm); // return the current mass(t-t0) of the star in NB units - -return(tmp); -} -/****************************************************************************/ -#endif - - /****************************************************************************/ void read_data() { inp = fopen(inp_fname,"r"); fscanf(inp,"%d \n", &diskstep); -//#ifdef STEVOL -// fscanf(inp,"%d \n", &N); -//#endif -// fscanf(inp,"%d %d %d %d \n", &N, &N_bh, &N_gmc, &N_star); - fscanf(inp,"%d \n", &N); fscanf(inp,"%lE \n", &time_cur); for (i=0; i= (mass_frac[k]*m_tot) ) - { - lagr_rad[k] = var_sort[j]; - k++; - } -#endif - } - -fclose(out); - -// write out lagr-rad.dat - -#ifdef LAGR_RAD -out = fopen("lagr-rad.dat","a"); -fprintf(out,"%.6E ", time_cur); -for (k=0; k 0 allows velocity kick at BH formation - - value3_.idum = 19640916; // random number seed !!! - - //BSE internal parameters (see Hurley, Pols & Tout 2002) - - flags_.ceflag = 0; //ceflag > 0 activates spin-energy correction in common-envelope (0) #defunct# - flags_.tflag = 1; //tflag > 0 activates tidal circularisation (1) - flags_.ifflag = 0; //ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0) - flags_.nsflag = 1; //nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1) - flags_.wdflag = 1; //wdflag > 0 uses modified-Mestel cooling for WDs (0) - - value5_.beta = 0.125; //beta is wind velocity factor: proportional to vwind**2 (1/8) - value5_.xi = 1.0; //xi is the wind accretion efficiency factor (1.0) - value5_.acc2 = 1.5; //acc2 is the Bondi-Hoyle wind accretion factor (3/2) - value5_.epsnov = 0.001; //epsnov is the fraction of accreted matter retained in nova eruption (0.001) - value5_.eddfac = 1.0; //eddfac is Eddington limit factor for mass transfer (1.0) - value5_.gamma = -1.0; //gamma is the angular momentum factor for mass lost during Roche (-1.0) - - value2_.alpha1 = 1.0; //alpha1 is the common-envelope efficiency parameter (1.0) - value2_.lambda = 0.5; //lambda is the binding energy factor for common envelope evolution (0.5) - - #endif // STEVOL_SSE - - /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); @@ -4655,31 +2778,6 @@ int main(int argc, char *argv[]) fflush(stdout); } /* if (myRank == rootRank) */ - #ifdef STEVOL - dt_stevol = dt_max; - t_stevol = dt_stevol*(1.0 + floorl(time_cur/dt_stevol)); - - if (myRank == rootRank) - { - printf("t_stevol = %.6E dt_stevol = %.6E \n", t_stevol, dt_stevol); - printf("\n"); - fflush(stdout); - } /* if (myRank == rootRank) */ - #endif - - #ifdef STEVOL_SSE - dt_stevol = dt_max; - t_stevol = dt_stevol*(1.0 + floorl(time_cur/dt_stevol)); - - if (myRank == rootRank) - { - printf("t_stevol = %.6E dt_stevol = %.6E \n", t_stevol, dt_stevol); - printf("\n"); - fflush(stdout); - } /* if (myRank == rootRank) */ - #endif - - for (i=0; i 0.1) dt[i] = min_dt; - } /* i */ - - #endif // GMC2 - - - #ifdef GMC2222 - - /* define the max. dt for the GMC... */ - - for (i=1; i 0.1) dt[i] = dt_max; - } /* i */ - - #endif // GMC2 - - - #ifdef ACT_DEF_LL CreateLinkList(); - #ifdef DEBUG111 - if ( myRank == rootRank ) check_linklist(0); - #endif #endif @@ -5311,45 +3236,6 @@ int main(int argc, char *argv[]) g6_flush_jp_buffer(clusterid); #endif - - - #ifdef STEVOL_SSE - /* Define the current mass of all star particles... */ - - for (i=0; i (R_CRIT*R_CRIT)) - { - hz = HZ; - } - else - { - hz = HZ*(sqrt(r_new_drag2)/R_CRIT); - } - - if (fabs(x_act_new[i][2]) > hz*2.2) - { - - if (z_new_drag * x_act_new[i][2] < 0.0) - { - dt_tmp *= 0.125; - } - else - { - if (fabs(z_new_drag) < hz*1.8) - { - dt_tmp *= 0.5; - } - } - - } - - #endif - - if ( (dt_new > 2.0*dt_tmp) && (fmod(min_t, 2.0*dt_tmp) == 0.0) && (2.0*dt_tmp <= dt_max) ) @@ -5989,55 +3778,6 @@ int main(int argc, char *argv[]) - #ifdef GMC - - min_dt = dt_act[0]; - - for (i=1; i 0.1) dt_act[i] = min_dt; - } /* i */ - - - #endif // GMC2 - - - #ifdef GMC2222 - - /* define the max. dt for the GMC... */ - - for (i=1; i 0.1) dt_act[i] = dt_max; - } /* i */ - - #endif // GMC2 - - - - - - #ifdef BBH_INF if (myRank == rootRank) { @@ -6354,113 +4094,6 @@ int main(int argc, char *argv[]) n_act_sum += n_act; n_act_distr[n_act-1]++; - - /* STEVOL routine on all the nodes */ - - #ifdef STEVOL - - #ifdef TIMING - get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); - #endif - - if (time_cur >= t_stevol) - { - - /* Define the current mass of all star particles... */ - - for (i=0; i= t_stevol) */ - - #ifdef TIMING - get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); - DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0); - #endif - - #endif - - - /* STEVOL_SSE routine on all the nodes */ - - #ifdef STEVOL_SSE - - - #ifdef TIMING - get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); - #endif - - if (time_cur >= t_stevol) - { - - /* Define the current mass of all star particles... */ - - for (i=0; i= t_stevol) */ - - #ifdef TIMING - get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); - DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0); - #endif - - - - - /* Update ALL the new m,r,v "j" part. to the GRAPE/GPU memory */ - - - for (j=0; j= t_bh) { @@ -6547,15 +4180,6 @@ int main(int argc, char *argv[]) #endif // STARDESTR_EXT - - - #ifdef ACT_DEF_LL - #ifdef DEBUG111 - if (myRank == rootRank) check_linklist((int)Timesteps); - #endif - #endif - - energy_contr(); #ifdef CMCORR111 @@ -6576,12 +4200,6 @@ int main(int argc, char *argv[]) write_cont_data(); - #ifdef GMC - write_cont_GMC(); - #endif - - - /* possible OUT for timing !!! */ #ifdef TIMING @@ -6608,24 +4226,6 @@ int main(int argc, char *argv[]) fclose(out); #endif - - - #ifdef DEBUG111 - tmp_file = fopen("n_act_distr.dat","w"); - - fprintf(tmp_file,"\n"); - fprintf(tmp_file,"Total Timesteps = %.0f Total sum of integrated part. = %.0f g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls); - fprintf(tmp_file,"\n"); - fprintf(tmp_file,"Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09); - - for (i=1;i