/*****************************************************************************
  File Name      : "Star Destr.c"
  Contents       : star "destruction" by tidal field of "live" BH (1 or 2)
  Coded by       : Peter Berczik
  Last redaction : 2010.IX.14 1:43PM
*****************************************************************************/

void star_destr(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 i_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, eps2, eps_bh2, rsb, rkb2, rks2, xp[3], v_bh[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;


if( time < t_diss_on ) return;


eps_bh = eps;

eps2    = SQR(eps);
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);
  }
*/

#ifdef ADD_BH2
  n_end = n-2;
  N_end = N-2;
#else
#ifdef ADD_BH1
  n_end = n-1;
  N_end = N-1;
#endif		// ADD_BH1
#endif		// ADD_BH2

/*
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]-x[i_bh][0]) + SQR(x[i][1]-x[i_bh][1]) + SQR(x[i][2]-x[i_bh][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) && (m[i] != 0.0) )
    {


    if(myRank == rootRank)
      {

      out = fopen("accr-data.dat","a");

      fprintf(out,"%04d %04d %04d %04d \t %.6E %.6E %.6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \n",
      i, ind[i], i_bh, ind[i_bh], 
      time,        
      m[i],        *m_bh, 
      x[i][0],     x[i][1],     x[i][2], 
      x[i_bh][0],  x[i_bh][1],  x[i_bh][2], 
      v[i][0],     v[i][1],     v[i][2],
      v[i_bh][0],  v[i_bh][1],  v[i_bh][2],
      pot[i], 
      a[i][0],     a[i][1],     a[i][2],
      adot[i][0],  adot[i][1],  adot[i][2]);

      fclose(out);

      } /* if(myRank == rootRank) */


    v_bh[0] = (m[i] * v[i][0] + *m_bh * v[i_bh][0]) / (m[i] + *m_bh);
    v_bh[1] = (m[i] * v[i][1] + *m_bh * v[i_bh][1]) / (m[i] + *m_bh);
    v_bh[2] = (m[i] * v[i][2] + *m_bh * v[i_bh][2]) / (m[i] + *m_bh);

    e_kin  =  m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
    e_kin += *m_bh * ( SQR(v[i_bh][0]) + SQR(v[i_bh][1]) + SQR(v[i_bh][2]) );
    e_kin -= (m[i] + *m_bh) * ( SQR(v_bh[0]) + SQR(v_bh[1]) + SQR(v_bh[2]) );
    e_kin *= 0.5;

//    e_kin = 0.5*m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) - SQR(v[i_bh][0]) - SQR(v[i_bh][1]) - SQR(v[i_bh][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]-x[i_bh][0]) + SQR(x_N[k][1]-x[i_bh][1]) + SQR(x_N[k][2]-x[i_bh][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]-x[i_bh][0]) + SQR(xp[1]-x[i_bh][1]) + SQR(xp[2]-x[i_bh][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: %.6E %04d %06d   %.4E %.4E %04d %06d   %.4E % .4E % .4E % .4E % .4E   %.4E %04d \n", 
            time, i_bh, ind[i_bh], rsb, R_t, i, ind[i], 
	    e_kin, e_pot_BH, e_pot, e_corr, E_corr, *m_bh, *num_bh);
      fflush(stdout);
      }


    v[i_bh][0] = v_bh[0];
    v[i_bh][1] = v_bh[1];
    v[i_bh][2] = v_bh[2];


//    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 */


  m[i_bh] = *m_bh;

/*
  if(myRank == rootRank)
    {
    printf("%06d \t %.6E \n", i_bh, m[i_bh]);
    fflush(stdout);
    }
*/

}

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