yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
phigrape/star_destr_ext.c
Yohai aaf40315c8 first commit
2019-08-25 16:22:03 +08:00

286 lines
6.7 KiB
C
Raw Blame History

/*****************************************************************************
File Name : "Star Destr EXT.c"
Contents : star "destruction" by tidal field of EXT fixed BH
Coded by : Peter Berczik
Last redaction : 2010.IX.14 1:43PM
*****************************************************************************/
void star_destr_ext(double time,
int n,
int ind[],
double m[],
double x[][3],
double v[][3],
double pot[],
double a[][3],
double adot[][3],
double t[],
double dt[],
int N,
double m_N[],
double x_N[][3],
double v_N[][3],
double a_N[][3],
double adot_N[][3],
double t_N[],
double *m_bh,
int *num_bh)
{
int n_end, k_act, N_end;
double R_t, e_kin=0.0, e_pot_BH=0.0, e_pot=0.0, e_corr=0.0;
double eps_bh, eps_bh2, eps2, rsb, rkb2, rks2, xp[3];
//double vir = 0.6, gamma = 1000.0;
double x_max = 1.0E+03, v_max = 1.0E-08,
a_max = 1.0E-08, adot_max = 1.0E-08;
double a_ext_i[3], adot_ext_i[3];
if(time < t_diss_on) return;
eps2 = SQR(eps);
eps_bh = b_bh;
eps_bh2 = SQR(eps_bh);
/*
m_s = 1.0/N;
R_s = 2.52E-08/2.507328103;
R_t = gamma * R_s * pow( 2.0*(*m_bh/m_s), over3 );
*/
R_t = R_TIDAL;
/*
if(myRank == rootRank)
{
printf("%.6E \t %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, i_bh, ind[i_bh], R_t, *m_bh, *num_bh);
fflush(stdout);
}
*/
n_end = n;
N_end = N;
/*
if(myRank == rootRank)
{
printf("%.6E \t %06d %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, n_end, i_bh, ind[i_bh], R_t, *m_bh, *num_bh);
fflush(stdout);
}
*/
for(i=0; i<n_end; i++)
{
rsb = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
// R_t = (R_s/5.848035) * pow( (2.0*(*m_bh)/m_s), over3 );
// if( (r < R_t) && (e_kin < ABS(vir*e_pot)) )
if( (rsb < R_t) )
{
if(myRank == rootRank)
{
x_ij = x[i][0];
y_ij = x[i][1];
z_ij = x[i][2];
vx_ij = v[i][0];
vy_ij = v[i][1];
vz_ij = v[i][2];
r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij) + eps_bh2;
r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = *m_bh / (r * r2);
a_ext_i[0] = -tmp * x_ij;
a_ext_i[1] = -tmp * y_ij;
a_ext_i[2] = -tmp * z_ij;
adot_ext_i[0] = -tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot_ext_i[1] = -tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot_ext_i[2] = -tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
out = fopen("accr-data.dat","a");
ii = ind[i];
#ifndef STARDISK
a_drag[ii][0] = 0.0;
a_drag[ii][1] = 0.0;
a_drag[ii][2] = 0.0;
adot_drag[ii][0] = 0.0;
adot_drag[ii][1] = 0.0;
adot_drag[ii][2] = 0.0;
#endif // no STARDISK
fprintf(out, "%04d %04d \t %.6E %.6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \n",
i, ii,
time_cur, *m_bh,
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
a_ext_i[0], a_ext_i[1], a_ext_i[2],
adot_ext_i[0], adot_ext_i[1], adot_ext_i[2],
a_drag[ii][0], a_drag[ii][1], a_drag[ii][2],
adot_drag[ii][0], adot_drag[ii][1], adot_drag[ii][2]);
fclose(out);
} /* if(myRank == rootRank) */
e_kin = 0.5*m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
e_pot_BH = - *m_bh * m[i] / sqrt( SQR(rsb) + eps_bh2 );
tmp = 0.0;
for(k=0; k<N_end; k++)
{
if( k != ind[i] )
{
/* define if k prinadlezhit ACT ili net */
k_act = 0;
for(ii=0; ii<n_end; ii++)
{
if(k==ind[ii])
{
k_act=1; break;
}
} /* ii */
/*
if(myRank == rootRank)
{
printf("0: \t %04d %04d \t %04d %04d \n", i, ind[i], k, k_act);
fflush(stdout);
}
*/
if(k_act == 1) /* k prinadlezhit ACT */
{
rkb2 = SQR(x_N[k][0]) + SQR(x_N[k][1]) + SQR(x_N[k][2]);
tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
rks2 = SQR(x_N[k][0]-x[i][0]) + SQR(x_N[k][1]-x[i][1]) + SQR(x_N[k][2]-x[i][2]);
tmp += m_N[k]/sqrt( rks2 + eps2 );
}
else /* k ne prinadlezhit ACT */
{
dt_tmp = time - t_N[k];
dt2half = over2*dt_tmp*dt_tmp;
dt3over6 = over3*dt_tmp*dt2half;
xp[0] = x_N[k][0] + v_N[k][0]*dt_tmp + a_N[k][0]*dt2half + adot_N[k][0]*dt3over6;
xp[1] = x_N[k][1] + v_N[k][1]*dt_tmp + a_N[k][1]*dt2half + adot_N[k][1]*dt3over6;
xp[2] = x_N[k][2] + v_N[k][2]*dt_tmp + a_N[k][2]*dt2half + adot_N[k][2]*dt3over6;
rkb2 = SQR(xp[0]) + SQR(xp[1]) + SQR(xp[2]);
tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
rks2 = SQR(xp[0]-x[i][0]) + SQR(xp[1]-x[i][1]) + SQR(xp[2]-x[i][2]);
tmp += m_N[k]/sqrt( rks2 + eps2 );
} /* if(k_act == 1) */
/*
if(myRank == rootRank)
{
printf("\t %04d %04d \t %04d \t % .8E \n", k, ind[i], k_act, tmp);
fflush(stdout);
}
*/
} /* if( k != ind[i] ) */
} /* k */
e_pot = - m[i] * tmp;
e_corr = e_kin + e_pot_BH + e_pot;
E_corr += e_corr;
*m_bh += m[i]; // add the star mass to the BH mass
*num_bh = *num_bh + 1;
if(myRank == rootRank)
{
printf("ACCR: %.4E %.4E %.4E %04d %06d %.4E % .4E % .4E % .4E % .4E %.4E %04d \n",
time, rsb, R_t, i, ind[i], e_kin, e_pot_BH, e_pot, e_corr, E_corr, *m_bh, *num_bh);
fflush(stdout);
}
// tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
// tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
// tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
/* for part "i" */
x[i][0] = x_max + 1.0*my_rand2();
x[i][1] = x_max + 1.0*my_rand2();
x[i][2] = x_max + 1.0*my_rand2();
v[i][0] = v_max*my_rand2();
v[i][1] = v_max*my_rand2();
v[i][2] = v_max*my_rand2();
pot[i] = 0.0;
a[i][0] = a_max*my_rand2();
a[i][1] = a_max*my_rand2();
a[i][2] = a_max*my_rand2();
adot[i][0] = adot_max*my_rand2();
adot[i][1] = adot_max*my_rand2();
adot[i][2] = adot_max*my_rand2();
m[i] = 0.0;
t[i] = t[i] + 0.125;
dt[i] = 0.125;
} /* if( (r < R_t) ) */
} /* i */
}
/***************************************************************/
/***************************************************************/
/***************************************************************/
Reference in a new issue View git blame Copy permalink