/*ddrift compute grains drift velocity */
#include "cddefines.h"
#include "physconst.h"
#include "sqrne.h"
#include "grainvar.h"
#include "rfield.h"
#include "phycon.h"
#include "tepowers.h"
#include "h.h"
#include "he.h"
#include "negdrg.h"
#include "trace.h"
#include "ddrift.h"

void ddrift(void)
{
	long int i, 
	  loop, 
	  nd;

	double alam, 
	  corr, 
	  dmomen, 
	  fac, 
	  fdrag, 
	  g0, 
	  g2, 
	  phi2lm, 
	  psi, 
	  rdust, 
	  si, 
	  vdold, 
	  volmom;

#	ifdef DEBUG_FUN
	fputs( "<+>ddrift()\n", debug_fp );
#	endif

	for( nd=0; nd < NDUST; nd++ )
	{
		if( GrainVar.lgDustOn1[nd] )
		{
			/* find momentum absorbed by grain */
			dmomen = 0.;
			for( i=0; i < rfield.nflux; i++ )
			{
				dmomen += (rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i])*
				  rfield.anu[i]*(GrainVar.dqabs[nd][i] + GrainVar.dqscat[nd][i]);
			}
			dmomen *= EN1RYD;

			/* now find force on grain, and drift velocity
			 * factor is 2*k */
			fac = 2.7612e-16*phycon.te;

			/* now PSI defined by Draine and Salpeter 79 Ap.J. 231, 77 (1979) */
			psi = (GrainVar.dstpot[nd]*EVRYD)*1.1606e4/phycon.te;
			if( psi > 0. )
			{
				rdust = 1e-6;
				alam = log(20.702/rdust/psi*TePowers.sqrte/sqrne.SqrtEden);
			}
			else
			{
				alam = 0.;
			}

			phi2lm = POW2(psi)*alam;
			corr = 2.;
			loop = 0;
			while( loop < 10 && fabs(corr-1.) > 0.03 )
			{
				loop += 1;
				vdold = GrainVar.DustDftVel[nd];

				/* interactions with protons */
				si = GrainVar.DustDftVel[nd]/TePowers.sqrte*7.755e-5;
				g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
				g2 = si/(1.329 + POW3(si));

				/* drag force due to protons, both linear and square in velocity
				 * equation 4 from D+S Ap.J. 231, p77. */
				fdrag = fac*h.hii*(g0 + phi2lm*g2);

				/* drag force due to interactions with electrons */
				si = GrainVar.DustDftVel[nd]/TePowers.sqrte*1.816e-6;
				g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
				g2 = si/(1.329 + POW3(si));
				fdrag += fac*phycon.eden*(g0 + phi2lm*g2);

				/* drag force due to collisions with hydrogen and helium atoms */
				si = GrainVar.DustDftVel[nd]/TePowers.sqrte*7.755e-5;
				g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
				fdrag += fac*(h.hi + 1.1*he.hei)*g0;

				/* drag force due to interactions with helium ions */
				si = GrainVar.DustDftVel[nd]/TePowers.sqrte*1.551e-4;
				g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
				g2 = si/(1.329 + POW3(si));
				fdrag += fac*he.heii*(g0 + phi2lm*g2);

				/* this term does not work
				 *  2      HEIII*(G0+4.*PSI**2*(ALAM-0.693)*G2) )
				 * this is total momentum absorbed by dust per unit vol */
				volmom = dmomen/2.998e10;

				if( fdrag > 0. )
				{
					corr = sqrt(volmom/fdrag);
					GrainVar.DustDftVel[nd] = (float)(vdold*corr);
				}
				else
				{
					corr = 1.;
					negdrg.lgNegGrnDrg = TRUE;
				}

				if( trace.lgTrace && trace.lgDustBug )
				{
					fprintf( ioQQQ, "     %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n", 
					  loop, GrainVar.DustDftVel[nd], volmom );
				}
			}
		}
	}


#	ifdef DEBUG_FUN
	fputs( " <->ddrift()\n", debug_fp );
#	endif
	return;
}