/*****************************************************************************
  File Name      : "Drag Force.c"
  Contents       : Drag force calculus.
                 : Converted & optimized from Chingis Omarov Fortran code.
  Coded by       : Taras Panamarev, 
                 : Denis Yurin, Maxim Makukov & Peter Berczik :-)
  Last redaction : 2012-11-15 15:33
*****************************************************************************/

// main parameters

const double Md = 0.01;
const double s=4, beta_s=0.7, alpha=0.75, h=0.001;
const double r_lim2=(0.5*0.5), h_lim2=(5e-3*5e-3); 
const double gradP_rho = 0.90;	// account for -grad(P)/rho in gas disk... 

const double Qtot = 0.01;	//  16e3

/*
const double Qtot = 0.00547000;	//   8e3
const double Qtot = 0.00296975;	//  16e3
const double Qtot = 0.00193372;	//  32e3
const double Qtot = 0.00054215;	// 128KB
*/

int    K; 

double RDISC, //projection of radius vector to the disc plane
       h_z, // disc profile
       RDISC2, Z2, RDISC_3over2, R, RR0, RR02, RR0_s, RDISC_5over2, R_zeta, R0_zeta, RRcrit, RRcrit_zeta,Rcrit_zeta, R_zetaR, 
       inner_hole, temp, sqrt_m_bh, pot_drag, 
       sigma0, R02, hR02, Rcrit, R0, zeta, 
       R_DOT, //the first derivative of absolute value of projection
       COEFDENS, 
       RHO, RHODOT, //disc density and it's firs derivative with respect to time 
       VXDISC, VYDISC,//disc particle velocity components 
       VXRE, VYRE, VZRE, VRE,//the absolute value of relative (star-disc) velocity and it's components  
       coef, 
       FDRAGX, FDRAGY, FDRAGZ,//drag force components 
       DVXRE, DVYRE, DVZRE, DVRE,//first derivative of relative velocity and it's absolute value 
       KFDOT;

void init_drag_force() 
  {
  R0 = R_0;
  Rcrit = R_CRIT;

  sigma0 = (2-alpha)*Md/(TWOPi*sqrt_TWOPi*h*R0);
  R02 = SQR(R0);
  hR02 = SQR(h*R0);
  }

void calc_drag_force() 
  {

  if (min_t < t_diss_on) return;
//  if (min_t == t_diss_on) init_drag_force();
  init_drag_force();

  for (i=0; i<n_act; i++) 
    {
  
    K = ind_act[i];

    RDISC2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]);
    Z2 = SQR(x_act_new[i][2]);

    if( (RDISC2 < r_lim2) && (Z2 < h_lim2) )
      {

      RDISC2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]);
      R = sqrt(RDISC2 + SQR(x_act_new[i][2]));
        
      RDISC = sqrt(RDISC2);
      R_DOT = ( x_act_new[i][0]*v_act_new[i][0] + x_act_new[i][1]*v_act_new[i][1] )/RDISC;

      RR0 = RDISC/R0;
      RR02 = SQR(RR0);
//    RR0_s = pow(RR0, s);
      RR0_s = SQR(RR02);
      RRcrit = RDISC/Rcrit;

/*
      if (RDISC < Rcrit) zeta = 1; else zeta = 0;

      RRcrit_zeta = pow(RRcrit, zeta);
      R_zeta = pow(RDISC, -2*zeta);
      R0_zeta = pow(R0, -2*zeta);
      Rcrit_zeta = pow(Rcrit, -2*zeta);
*/

      if (RDISC < Rcrit)
        {
        zeta = 1;
        RRcrit_zeta = RRcrit;
        R_zeta = 1.0/RDISC2;
        R0_zeta = 1.0/R02;
        Rcrit_zeta = 1.0/SQR(Rcrit);
        }
      else 
        {
        zeta = 0; 
        RRcrit_zeta = 1.0;
        R_zeta = 1.0;
        R0_zeta = 1.0;
        Rcrit_zeta = 1.0;
        }

      R_zetaR = R_zeta/R;

      COEFDENS = pow(RR0, -alpha) * exp(-beta_s*RR0_s)/RRcrit_zeta;
      RHO = sigma0*COEFDENS*exp(-Z2/(2.0*hR02*SQR(RRcrit_zeta)));

      RHODOT = - R_DOT*(zeta+alpha+beta_s*s*RR0_s )/RDISC 
               - (v_act_new[i][2]*x_act_new[i][2]*R_zeta-zeta*R_zetaR*Z2*R_DOT)/(hR02*Rcrit_zeta);

      RDISC_3over2 = sqrt(RDISC2*RDISC);
      RDISC_5over2 = RDISC_3over2*RDISC;

      sqrt_m_bh = sqrt(m_bh);

      sqrt_m_bh *= gradP_rho;			// account for -grad(P)/rho in gas disk... 

      VXDISC = -sqrt_m_bh * x_act_new[i][1]/RDISC_3over2;
      VYDISC = +sqrt_m_bh * x_act_new[i][0]/RDISC_3over2;

      VXRE = VXDISC - v_act_new[i][0];
      VYRE = VYDISC - v_act_new[i][1];
      VZRE =        - v_act_new[i][2];

      VRE = sqrt( SQR(VXRE) + SQR(VYRE) + SQR(VZRE) );

      coef = Qtot*RHO*VRE;

      FDRAGX = coef*VXRE;
      FDRAGY = coef*VYRE;
      FDRAGZ = coef*VZRE;

      DVXRE = -sqrt_m_bh * v_act_new[i][1]/RDISC_3over2
      	      +sqrt_m_bh * 1.5*x_act_new[i][1]*R_DOT/RDISC_5over2
      	      -a_act_new[i][0];

      DVYRE = +sqrt_m_bh * v_act_new[i][0]/RDISC_3over2
      	      -sqrt_m_bh * 1.5*x_act_new[i][0]*R_DOT/RDISC_5over2
      	      -a_act_new[i][1];

      DVZRE = -a_act_new[i][2];

      DVRE  = (VXRE*DVXRE + VYRE*DVYRE + VZRE*DVZRE)/VRE;

      KFDOT = Qtot*RHO;

      adot_drag[K][0] = KFDOT*(RHODOT*VXRE*VRE + DVRE*VXRE + DVXRE*VRE);
      adot_drag[K][1] = KFDOT*(RHODOT*VYRE*VRE + DVRE*VYRE + DVYRE*VRE);
      adot_drag[K][2] = KFDOT*(RHODOT*VZRE*VRE + DVRE*VZRE + DVZRE*VRE);

      pot_drag = -0.5*( (a_drag[K][0]+FDRAGX)*(x_act_new[i][0]-x_act[i][0]) +
                (a_drag[K][1]+FDRAGY)*(x_act_new[i][1]-x_act[i][1]) +
                (a_drag[K][2]+FDRAGZ)*(x_act_new[i][2]-x_act[i][2]) );

      E_sd += m[K]*pot_drag;

      a_drag[K][0] = FDRAGX;
      a_drag[K][1] = FDRAGY;
      a_drag[K][2] = FDRAGZ;

      }
    else
      {
      adot_drag[K][0] = 0.0;
      adot_drag[K][1] = 0.0;
      adot_drag[K][2] = 0.0;

      a_drag[K][0] = 0.0;
      a_drag[K][1] = 0.0;
      a_drag[K][2] = 0.0;
      } /* if( (RDISC2 < r_lim2) && (Z2 < h_lim2) ) */

  } /* i */


}
/***************************************************************/
/***************************************************************/
/***************************************************************/