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


}

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