/* UDFs for the deposition of sewage sludge particles and generation of slag flow */
/* !!!UMDIs and DPM-UDS must be set before loading!!! */

#define T_crit 1473		/* Temperature of critical viscosity in Kelvin - particle is starting to melt */
#define mu_crit	10.0		/* Viscosity at Tcrit */
#define T_melt 1673		/* Temperature in Kelvin at which particle is assumed to be fully molten */
#define mu_melt	0.1		/* Viscosity at T_melt */
#define Wecrit 1000000		/* Critical Weber number above which a particle will be reflected instead of trapped */
#define sigmasp 0.42		/* Surface tension for We number definition */
#define CUTOFF 0.5		/* If a particle encounters a slag volume fraction above this value, it will be deposited in the respective cell */
#define FACTOR 1		/* Factor to speed up the build-up of the slag flow */

/* Execute command to be called at every timestep: /define/user-defined execute-on-demand "reset_udmi::libudf" */
/* Execute command to be called at every timestep: /define/user-defined execute-on-demand "kill_phase_outlet::libudf" */
/* Execute command to be called at every timestep: /define/user-defined execute-on-demand "kill_phase_internal::libudf" */
/* Execute command to be called at every timestep: /define/user-defined execute-on-demand "kill_phase_burner::libudf" */

#include "udf.h"
#include "math.h"
#include "dpm_types.h"
#include "dpm_laws.h"
#include "mem.h"
#include "dpm.h"
#include "materials.h"

void deposit_func(Tracked_Particle *tp, cell_t c, Thread *t)
{
	/* calculate strenght of particle */
 	real strenght;
	if (dpm_par.unsteady_tracking)
 		strenght = TP_N(tp) / CURRENT_TIMESTEP;
 	else
 	{
		strenght = TP_FLOW_RATE(tp) / TP_INIT_MASS(tp);
		if (TP_STOCHASTIC(tp))
			strenght /= (real)TP_STOCHASTIC_NTRIES(tp);
	}

	/* REMOVAL OF PARTICLE */
	/*---------------------------------- */
	
	/* set UDMI with mass source */
		C_UDMI(c, t, 0) +=  strenght * TP_MASS(tp);

 	/* set UDMI with energy source */
 		C_UDMI(c, t, 1) +=  strenght * TP_MASS(tp) * 1561 * (TP_T(tp)-298.0);

	/* Mark particle for removal */
		MARK_PARTICLE(tp , P_FL_REMOVED);

	/* add particle momentum to continuous phase */
		dpms_t *cs = &(tp ->source);					
			cs ->momentum_s [0] += TP_MASS(tp)*TP_VEL(tp)[0];
			cs ->momentum_s [1] += TP_MASS(tp)*TP_VEL(tp)[1];	
			cs ->momentum_s [2] += TP_MASS(tp)*TP_VEL(tp)[2];	
	
	/* End tracking */
		TrackStepEnd(tp);

	/* Adds the particles to "trapped" counter of DPM iteration */
		tp ->gvtp.n_trapped += 1;

	/*Message("\nTrapped part. | cell temp.: %.1f K, part. temp.: %.1f K strenght: %.2f mass = %.6f energysource: %.6f", C_T(c,t), TP_T(tp), strenght, TP_MASS(tp), strenghth * TP_MASS(tp) * 1680 * (TP_T(tp) - C_T(c,t)));*/

}

DEFINE_DPM_BC(bc_deposit, tp, t, f, f_normal, dim)
{

	cell_t c;
	Thread *t1 , *tmix;
	
	c = F_C0(f, t); 								/* cell index of adjacent cell */
	t1 = THREAD_T0(t); 								/* cell thread of adjacent cell */
	
	tmix = (NNULLP(THREAD_SUPER_THREAD(t1))) ? THREAD_SUPER_THREAD(t1) : t1; 	/* mixture thread */
	Thread *phase_t = THREAD_SUB_THREAD(tmix , 1); 					/* secondary phase thread (0 for primary , 1 for secondary) */
	

	real alpha; 									/* angle of particle path with face normal */
	real vn=0.;
	real nor_coeff = 0.5;
	real tan_coeff = 0.5;
	real normal[3];
	int i, idim = dim;
	real NV_VEC(x);
	
	real traptime = 0.;
	real velmag = 0.;
	real distance = 0.;

#if RP_2D
	/* dim is always 2 in 2D compilation. Need special treatment for 2d
	axisymmetric and swirl flows */
	if (rp_axi_swirl)
	{
		real R = sqrt(TP_POS(tp)[1]*TP_POS(tp)[1] +
		TP_POS(tp)[2]*TP_POS(tp)[2]);

		if (R > 1.e-20)
		{
			idim = 3;
			normal[0] = f_normal[0];
			normal[1] = (f_normal[1]*TP_POS(tp)[1])/R;
			normal[2] = (f_normal[1]*TP_POS(tp)[2])/R;
		}
		else
		{
			for (i=0; i<idim; i++)
			normal[i] = f_normal[i];
		}
	}
	else
#endif
	for (i=0; i<idim; i++)
		normal[i] = f_normal[i];

		alpha = M_PI/2. - acos(MAX(-1.,MIN(1.,NV_DOT(normal,TP_VEL(tp))/
			MAX(NV_MAG(TP_VEL(tp)),DPM_SMALL))));
		if ((NNULLP(t)) && (THREAD_TYPE(t) == THREAD_F_WALL))
			F_CENTROID(x,f,t);

			real T_wall = C_T(f,t);
			real Weber = TP_USER_REAL(tp,2);
			/* Message("\nPoint of impact: x-coord = %f y-coord = %f z-coord = %f", x[0], x[1], x[2]); */		
			/* Message("\nThe temperature of particle %f at the time of impact is %f K", TP_ID(tp), TP_T(tp));*/
			/* Message("\nThe velocity at the time of impact is %f m/s", sqrt(TP_VEL(tp)[0]*TP_VEL(tp)[0]+TP_VEL(tp)[1]*TP_VEL(tp)[1]+TP_VEL(tp)[2]*TP_VEL(tp)[2]));*/
		
		
		/* Non-sticky particle - non sticky wall --> Reflect */
		if( (TP_T(tp) < T_crit) && (T_wall < T_crit) )
		{	

			/* calculate the normal component, rescale its magnitude by the coefficient of restitution and subtract the change */
			
			for(i=0; i<idim; i++)
				vn += TP_VEL(tp)[i]*normal[i];
			
			for(i=0; i<idim; i++)				
				TP_VEL(tp)[i] -= vn*normal[i];				/* Subtract off normal velocity. */
			
			for(i=0; i<idim; i++)
				TP_VEL(tp)[i] *= tan_coeff;				/* Apply tangential coefficient of restitution. */
			
			for(i=0; i<idim; i++)
				TP_VEL(tp)[i] -= nor_coeff*vn*normal[i];		/* Add reflected normal velocity. */
				
			for(i=0; i<idim; i++)
				TP_VEL0(tp)[i] = TP_VEL(tp)[i];				/* Store new velocity in TP_VEL0 of particle */
			
			Message("\nReflected (NSP - NSW) - particle temp. %.1f K, wall temp. %.1f K, Weber Nr. %.2f", TP_T(tp), T_wall, Weber);
			return PATH_ACTIVE;
		}

		/* Non-sticky particle - sticky wall */
		else if( (TP_T(tp) < T_crit) && (T_wall >= T_crit) )
		{					
			/* We > Wecrit --> Reflect */
			if ( Weber > Wecrit)
			{ 
				/* calculate the normal component, rescale its magnitude by the coefficient of restitution and subtract the change */
			
				for(i=0; i<idim; i++)
					vn += TP_VEL(tp)[i]*normal[i];
			
				for(i=0; i<idim; i++)				
					TP_VEL(tp)[i] -= vn*normal[i];				/* Subtract off normal velocity. */
			
				for(i=0; i<idim; i++)
					TP_VEL(tp)[i] *= tan_coeff;				/* Apply tangential coefficient of restitution. */
			
				for(i=0; i<idim; i++)
					TP_VEL(tp)[i] -= nor_coeff*vn*normal[i];		/* Add reflected normal velocity. */
				
				for(i=0; i<idim; i++)
					TP_VEL0(tp)[i] = TP_VEL(tp)[i];				/* Store new velocity in TP_VEL0 of particle */
			
			Message("\nReflected (NSP - SW) - particle temp. %.1f K, wall temp. %.1f K, Weber Nr. %.2f", TP_T(tp), T_wall, Weber);
			return PATH_ACTIVE;
			}
			/* We < Wecrit --> Deposit after full burnout*/
			else
			{
				deposit_func(tp, c, tmix);
			} 
		}

		/* Sticky particle - non-sticky wall */
		else if( (TP_T(tp) > T_crit) && (T_wall <= T_crit) )
		{
			/* We > Wecrit --> Reflect */
			if ( Weber > Wecrit)
			{ 
				/* calculate the normal component, rescale its magnitude by the coefficient of restitution and subtract the change */
			
				for(i=0; i<idim; i++)
					vn += TP_VEL(tp)[i]*normal[i];
			
				for(i=0; i<idim; i++)				
					TP_VEL(tp)[i] -= vn*normal[i];				/* Subtract off normal velocity. */
			
				for(i=0; i<idim; i++)
					TP_VEL(tp)[i] *= tan_coeff;				/* Apply tangential coefficient of restitution. */
			
				for(i=0; i<idim; i++)
					TP_VEL(tp)[i] -= nor_coeff*vn*normal[i];		/* Add reflected normal velocity. */
				
				for(i=0; i<idim; i++)
					TP_VEL0(tp)[i] = TP_VEL(tp)[i];				/* Store new velocity in TP_VEL0 of particle */
			
			Message("\nReflected (SP - NSW) - particle temp. %.1f K, wall temp. %.1f K, Weber Nr. %.2f", TP_T(tp), T_wall, Weber);
			return PATH_ACTIVE;
			}
			/* We < Wecrit --> Deposit after full burnout */
			else
			{
				deposit_func(tp, c, tmix);

			} 
		}
		/* Sticky Particle - sticky wall --> Deposit after full burnout */
		else
		{	
				deposit_func(tp, c, tmix);
		}		
	
}

DEFINE_ON_DEMAND(reset_udmi)
  {
	Domain *domain = Get_Domain(ROOT_DOMAIN_ID);
	Thread *t;
	cell_t c;

	Alloc_Storage_Vars(domain, SV_RTMP_0, SV_NULL);

	thread_loop_c(t,domain) {
		if (THREAD_ID(t)!=0) {
			begin_c_loop(c,t) {
				C_UDMI(c,t,0) = 0;
				C_UDMI(c,t,1) = 0;
			}
		}
	end_c_loop(c,t)

	}
  }

DEFINE_ON_DEMAND(kill_phase_surface)
  {
	Domain *domain = Get_Domain(ROOT_DOMAIN_ID);
	Thread *t;
	cell_t c;

	int fluid_ID = 70;
	Thread *tf = Lookup_Thread(domain,fluid_ID);
	Thread *mphase_tf = THREAD_SUB_THREAD(tf , 1);

	thread_loop_c(t,domain) 
	{
			begin_c_loop(c,tf) 
			{
				C_VOF(c,mphase_tf) = 0;
			}
	}
	end_c_loop(c,tf)
	Message0("\nSlag phase fraction on surface zone reset.");
	
  }

DEFINE_ON_DEMAND(kill_phase_internal)
  {
	Domain *domain = Get_Domain(ROOT_DOMAIN_ID);
	Thread *t;
	cell_t c;

	int fluid_ID = 19;
	Thread *tf = Lookup_Thread(domain,fluid_ID);
	Thread *mphase_tf = THREAD_SUB_THREAD(tf , 1);

	thread_loop_c(t,domain) 
	{
			begin_c_loop(c,tf) 
			{
				C_VOF(c,mphase_tf) = 0;
			}
	}
	end_c_loop(c,tf)
	Message0("\nSlag phase fraction in reactor zone away from wall reset.");
	
  }

DEFINE_ON_DEMAND(set_udmi)
  {
	Domain *domain = Get_Domain(ROOT_DOMAIN_ID);
	Thread *t;
	cell_t c;

	Alloc_Storage_Vars(domain, SV_RTMP_0, SV_NULL);

	thread_loop_c(t,domain) {
		if (THREAD_ID(t)!=0) {
			begin_c_loop(c,t) {
				C_UDMI(c,t,0) = 1;
			}
		}
	end_c_loop(c,t)
	}
  }

DEFINE_SOURCE(slag_mass , c, t, dS , eqn)
{
	real source;
	/*real dS; */
	if (C_UDMI(c, t, 0))
	{
		source = FACTOR * C_UDMI(c, t, 0) / MAX(C_VOLUME(c, t), 1e-14);
	}
	else
		source = 0.;

	dS[eqn] = 0;
	return source;
}


DEFINE_SOURCE(slag_energy , c, t, dS , eqn)
{
	real source;
	/*real dS; */
	if (C_UDMI(c, t, 0))
	{
		source = FACTOR * C_UDMI(c, t, 1) / MAX(C_VOLUME(c, t), 1e-14);
	}
	else
		source = 0.;

	dS[eqn] = 0;
	return source;
}

DEFINE_SOURCE(gas_mass , c, t, dS , eqn)
{
	real source;
	/*real dS; */
	if (C_UDMI(c, t, 0))
	{
		Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) :t; /* mixture thread */
		Thread *mphase_t = THREAD_SUB_THREAD(tmix , 1);
		Thread *gphase_t = THREAD_SUB_THREAD(tmix , 0);

		source = FACTOR * (-C_UDMI(c, t, 0)) / MAX(C_VOLUME(c, t), 1e-14)*C_R(c, gphase_t) / C_R(c, mphase_t);
	}
	else
		source = 0.;

	dS[eqn] = 0;
	return source;
}

DEFINE_DPM_DRAG(modifieddraglaw_spherical, Re, tp)
{
	real Cd, drag_force = 0.0;
	real a1, a2, a3;
   	real corr = 1.0;  
 	real Rec;

    	/* Cell Thread check */
    	Thread *t = TP_CELL_THREAD(tp);
    	cell_t c = TP_CELL(tp);
    	Thread *mphase_t = THREAD_SUB_THREAD(t, 1); 
    	Thread *gphase_t = THREAD_SUB_THREAD(t, 0); 


	if(C_VOF(c,mphase_t) < CUTOFF)
	{
		/* Modify Re number */
		/* Re = rho *d_p* abs(u_p -u)/mu */
 		Rec = Re * C_R(c, gphase_t) / C_R(c, t) * C_MU_L(c, t) / C_MU_L(c, gphase_t);
 		Re = Rec;
 		corr = C_MU_L(c, gphase_t) / C_MU_L(c, t);
	}	

    		if (Re < 0.01)
    		{
        		a1 = 0.0;
        		a2 = 24.0;
        		a3 = 0.0;
    		}
    		else if (Re < 1.0)
   		{
    	    		a1 = 3.69;
        		a2 = 22.73;
        		a3 = 0.0903;
    		}
    		else if (Re < 10.0)
    		{
        		a1 = 1.222;
        		a2 = 29.1667;
        		a3 = -3.8889;
    		}
    		else if (Re < 100.0)
    		{
        		a1 = 0.6167;
        		a2 = 46.5;
        		a3 = -116.67;
    		}
    		else if (Re < 1000.0)
    		{
        		a1 = 0.3644;
        		a2 = 98.33;
        		a3 = -2778.0;
    		}
    		else if (Re < 5000.0)
    		{
        		a1 = 0.357;
        		a2 = 148.62;
        		a3 = -47500.0;
    		}
    		else if (Re < 10000.0)
    		{
        		a1 = 0.46;
        		a2 = -490.546;
        		a3 = 578700.0;
    		}
    		else
    		{
        		a1 = 0.5191;
        		a2 = -1662.5;
        		a3 = 5416700.0;
    		}

    		Cd = a1 + a2 / Re + a3 / SQR(Re);

    		/* Corrected drag force calculation */

		if (TP_USER_REAL(tp,9) == 0)
		{
    			drag_force = 18.0 * Cd * Re / 24.0 * corr;
		}
		else
		{
			drag_force = 0;
		}
		return drag_force;
}


DEFINE_INIT(udm_setup, domain)
{
  /* if memory for the particle variable titles and UDMs has not been
   * allocated yet, do it now */
  	if (NULLP(user_particle_vars)) Init_User_Particle_Vars();
 		/* now set the name and label */
  		strcpy(user_particle_vars[0].name,"biot-number");
  		strcpy(user_particle_vars[0].label,"Biot number");
  		strcpy(user_particle_vars[1].name,"particle-charfrac");
  		strcpy(user_particle_vars[1].label,"Particle char fraction");
		strcpy(user_particle_vars[2].name,"weber-number");
  		strcpy(user_particle_vars[2].label,"Weber number");
		strcpy(user_particle_vars[3].name,"impact-x");
  		strcpy(user_particle_vars[3].label,"impact-x");
		strcpy(user_particle_vars[4].name,"impact-y");
  		strcpy(user_particle_vars[4].label,"impact-y");
		strcpy(user_particle_vars[5].name,"impact-z");
  		strcpy(user_particle_vars[5].label,"impact-z");
		strcpy(user_particle_vars[6].name,"impact-vel-x");
  		strcpy(user_particle_vars[6].label,"impact-vel-x");
		strcpy(user_particle_vars[7].name,"impact-vel-y");
  		strcpy(user_particle_vars[7].label,"impact-vel-y");
		strcpy(user_particle_vars[8].name,"impact-vel-z");
  		strcpy(user_particle_vars[8].label,"impact-vel-z");
		strcpy(user_particle_vars[9].name,"deposition-flag");
  		strcpy(user_particle_vars[9].label,"Deposition flag");
		strcpy(user_particle_vars[10].name,"dist-from-impact");
  		strcpy(user_particle_vars[10].label,"Dist. fr. init. wall cont.");
		strcpy(user_particle_vars[11].name,"part-res-time");
  		strcpy(user_particle_vars[11].label,"Particle residence time");
		Set_User_Memory_Name(0,"DPM mass source");
		Set_User_Memory_Name(1,"DPM enthalpy source");
		Message0("\nDPM-UDV and UDM names have been assigned.");
}

DEFINE_DPM_PROPERTY(diffusivity_wall_burning, c, t, tp, T)
{	

    	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
    	Thread *mphase_t = THREAD_SUB_THREAD(tmix, 1);
    	Thread *gphase_t = THREAD_SUB_THREAD(tmix, 0);

	/* Diffusivity */
	real D;
	
	if ( C_VOF(c,mphase_t) > 0)
	{
		D = 0.0;
	}
	else
	{
		if(TP_USER_REAL(tp,9) == 0)
		{
			/* D = 4e-5; */
				D = 4e-5;
		}
		else
		{
			/* D = 2e-5; */
				D = 2e-5;
		}
	}

    	return D;
}

/* Used to update the user scalar variables and to handle in-cell deposition (not at wall)*/
DEFINE_DPM_SCALAR_UPDATE(scalars_cell_deposit, cell, thread, initialize, tp)
{

	/*Cell Thread check ?*/
	Thread *t = TP_CELL_THREAD(tp);
	cell_t c = TP_CELL(tp);
	Thread *mphase_t = THREAD_SUB_THREAD(t, 1); 					/* melt phase */
	Thread *gphase_t = THREAD_SUB_THREAD(t, 0); 					/* gas phase */
				
  	/* cphase_state_t *c = &(tp->cphase); */

	
	real lambdaP = 1.0;						/* Thermal conductivity of particle material (order of magnitude, INSPYRO analysis) */	
	real psi = 0.806;						/* Sphericity of the particles (TUG measurement)*/				
										
		real v_rel = sqrt(pow((P_VEL(tp)[0] - C_U(c,gphase_t)),2) + pow((P_VEL(tp)[1] - C_V(c,gphase_t)),2) + pow((P_VEL(tp)[2] - C_W(c,gphase_t)),2)); 
									/* Determine relative velocity between particle and fluid for Reynolds number determination */							
		real rhoF = C_R(c,gphase_t); 				/* Determine gas density at particle location */	
		real mueF = C_MU_L(c,gphase_t);				/* Determine gas viscosity at particle location */	
		real lambdaF = C_K_L(c,gphase_t);			/* Determine gas thermal conductivity */	
		real cpF = C_CP(c,gphase_t);				/* Determine gas specific heat capacity */
		real ReP = P_DIAM(tp)*v_rel*rhoF/mueF;  		/* Particle reynolds number */	
		real Pr =  mueF*cpF/lambdaF;				/* Prandtl number */
		real Nu = (2+0.6*pow(ReP,0.5)*pow(Pr,0.3333));		/* Nusselt number (Ranz-Marshall correlation) */	
		real VP = P_MASS(tp)/P_RHO(tp);				/* Particle volume */	
		real AP = P_DIAM(tp)*P_DIAM(tp)*3.1415926/psi;		/* Approximate particle surface area, accounting for non-sphericity */	
		real L = VP/AP;						/* Characteristic lenght for Biot number */	
      		TP_USER_REAL(tp,0) = (Nu*lambdaF/P_DIAM(tp))*L/lambdaP;	/* assign biot number to the first user value */
		
		real particle_charfrac = TP_CF(tp);
		TP_USER_REAL(tp,1) = particle_charfrac;			/* assign charfrac time to the second user value */
		
		real v_part = sqrt(P_VEL(tp)[0]*P_VEL(tp)[0]+P_VEL(tp)[1]*P_VEL(tp)[1]+P_VEL(tp)[2]*P_VEL(tp)[2]);	
		real weber_number = P_RHO(tp)*P_DIAM(tp)*v_part*v_part/sigmasp;		
		TP_USER_REAL(tp,2) = weber_number;			/* assign weber number to the second user value */
		TP_USER_REAL(tp,11) = TP_TIME(tp)-TP_TIME_OF_BIRTH(tp);			/* assign particle residence time to user real 10 */


	real traptime;
	real velmag;
	real distance = 0;

	if (C_VOF(c, mphase_t) > CUTOFF && TP_USER_REAL(tp,9) == 0) 		
	{
		/* Store impact positions x, y, z at time of impact in user-reals - if particle is not deposited yet (Flag user real 6) */
			TP_USER_REAL(tp,3) = TP_POS(tp)[0];
			TP_USER_REAL(tp,4) = TP_POS(tp)[1];
			TP_USER_REAL(tp,5) = TP_POS(tp)[2];
		/* Store velocity components u, v, w at time of impact in user-reals - if particle is not deposited yet (Flag user real 6) */
			TP_USER_REAL(tp,6) = TP_VEL(tp)[0];
			TP_USER_REAL(tp,7) = TP_VEL(tp)[1];
			TP_USER_REAL(tp,8) = TP_VEL(tp)[2];
		/* Time at which particle was trapped */
			traptime = TP_TIME(tp) - CURRENT_TIME;
			/* Message("\nTraptime: %.5f", traptime);*/
		/* Velocity magnitude at the time when the particle was trapped */
			velmag = sqrt(TP_USER_REAL(tp,3)*TP_USER_REAL(tp,3)+TP_USER_REAL(tp,4)*TP_USER_REAL(tp,4)+TP_USER_REAL(tp,5)*TP_USER_REAL(tp,5));
			/* Message("\nvelmag: %.5f", velmag);*/
		/* Set flag = contact with slag detected */
			TP_USER_REAL(tp,9) = 1;		
	}

	/* Calculate distance the particle is travelling from its impact point with the slag surface */
			distance = sqrt((TP_POS(tp)[0]-TP_USER_REAL(tp,3))*(TP_POS(tp)[0]-TP_USER_REAL(tp,3)) + (TP_POS(tp)[1]-TP_USER_REAL(tp,4))*(TP_POS(tp)[1]-TP_USER_REAL(tp,4)) + (TP_POS(tp)[2]-TP_USER_REAL(tp,5))*(TP_POS(tp)[2]-TP_USER_REAL(tp,5)));
			TP_USER_REAL(tp,10) = distance;
	
	/* If contact with slag surface has been detected: */
	if(TP_USER_REAL(tp,9) == 1)
	{
		/* If particle is not finished, reverse back a little to avoid VOF stability issue. */
		if(particle_charfrac != 0)
		{						
			if( C_VOF(c, mphase_t) > 0 || distance < 0.001)
			{	
				/*Message("\nSending back from slag zone.");*/
				TP_VEL(tp)[0] = - TP_USER_REAL(tp,6);
				TP_VEL(tp)[1] = - TP_USER_REAL(tp,7);
				TP_VEL(tp)[2] = - TP_USER_REAL(tp,8);					
			}
			else
			{	
				/*Message("\nBack from slag zone and sticking.");*/
				TP_VEL(tp)[0] = 1e-5;
				TP_VEL(tp)[1] = 1e-5;
				TP_VEL(tp)[2] = 1e-5;			
			}
		}
		else
		{	
				/*Message("\nSending back to slag zone.");*/
				TP_VEL(tp)[0] = TP_USER_REAL(tp,6);
				TP_VEL(tp)[1] = TP_USER_REAL(tp,7);
				TP_VEL(tp)[2] = TP_USER_REAL(tp,8);		
		}
		
		if((C_VOF(c, mphase_t) > CUTOFF) && (TP_USER_REAL(tp,9) == 1) && particle_charfrac == 0)
		{	
			/*Message("\nA particle was deposited after full burnout.");*/
			deposit_func(tp, c, t);
		}

		if((TP_USER_REAL(tp,11) > 60) && (TP_USER_REAL(tp,9) == 1))
		{	
			/*Message("\nA particle was deleted after exceeding the maximum tracking time.");*/
			TP_VF(tp) = 0.;
			TP_CF(tp) = 0.;
			deposit_func(tp, c, t);
		}

		
	}
}

DEFINE_SOURCE(gravity_x,c,t,dS,eqn)
{

	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
	Thread *mphase_t = THREAD_SUB_THREAD(tmix, 1);
    	Thread *gphase_t = THREAD_SUB_THREAD(tmix, 0);
	
	real gravity = 0.;
	real source =0.;

	real x[ND_ND];
	C_CENTROID(x,c,t);
	
 	if( C_VOF(c, mphase_t) > 0.01 )
		gravity = 0.5*cos(atan(x[1]/x[0]));
	else
		gravity = 0.0;

	/* gravity*cr(ct) is gravity*density = N/m³ = specific weight = F(g) on a unit of volume of a fluid */
	source = gravity*C_R(c,t);
		/* C_UDMI(c, t, 2) = source/C_R(c,t); */

	dS[eqn] = 0;
	return source;
}

DEFINE_SOURCE(gravity_y,c,t,dS,eqn)
{

	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
    	Thread *mphase_t = THREAD_SUB_THREAD(tmix, 1);
    	Thread *gphase_t = THREAD_SUB_THREAD(tmix, 0);
	
	real gravity = 0.;
	real source =0.;

	real x[ND_ND];
	C_CENTROID(x,c,t);

 	if( C_VOF(c, mphase_t) > 0.01 )
		gravity = 0.5*sin(atan(x[1]/x[0]));
	else
		gravity = 0.0;

	/* gravity*cr(ct) is gravity*density = N/m³ = specific weight = F(g) on a unit of volume of a fluid */
	source = gravity*C_R(c,t);
		/* C_UDMI(c, t, 3) = source/C_R(c,t); */	
		/* C_UDMI(c, t, 4) = pow( C_UDMI(c, t, 2)*C_UDMI(c, t, 2)+C_UDMI(c, t, 3)*C_UDMI(c, t, 3),2); */


	dS[eqn] = 0;
	return source;
}


DEFINE_SOURCE(gravity_z,c,t,dS,eqn)
{

	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
	
	real gravity = 0.;
	real source = 0.;

	real x[ND_ND];
	C_CENTROID(x,c,tmix);

 	gravity = -9.81;

	/* gravity*cr(ct) is gravity*density = N/m³ = specific weight = F(g) on a unit of volume of a fluid */
	source = gravity*C_R(c,tmix);
		

	dS[eqn] = 0;
	return source;
}

DEFINE_DPM_BODY_FORCE(particle_gravity,tp,i)

{
	real bforce=0;
	/* i == 1 --> 1 direction */
	if(i==1)
	{
		if(TP_USER_REAL(tp,9) == 1)
		{	
			/* Set zero if wall burning is active*/
			bforce = 0;
		}
		else
			bforce = -9.81*TP_MASS(tp);
	}

	/* an acceleration should be returned, this one will counter gravity to hold particle in place */
	return bforce/TP_MASS(tp);
}

DEFINE_PROPERTY(viscosity_geom_avg, c, t)
{	

    	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
    	Thread *mphase_t = THREAD_SUB_THREAD(tmix, 1);
    	Thread *gphase_t = THREAD_SUB_THREAD(tmix, 0);

	real mu_melt_visc = mu_melt;
	real mu_crit_visc = mu_crit;

	real T_melt_visc = T_melt;
	real T_crit_visc = T_crit;

	real mu_gas = 1e-5;
    	real mu_slag = mu_melt_visc;

    	real k;
	real d;

    	real mu_mix;
	real alpha;

    	real temp = C_T(c, tmix);
	alpha = C_VOF(c, mphase_t);

	/* Temperature dependent viscosity of the slag */ 

    		k = (mu_melt_visc - mu_crit_visc) / (T_melt_visc - T_crit_visc);
    		d = mu_crit_visc - k * T_crit_visc;

		if (temp > 0 && temp < T_crit_visc)
   			mu_slag = mu_crit_visc;
		else if (temp >= T_crit_visc && temp < T_melt_visc)
   		 	mu_slag = k * temp + d;
		else
    			mu_slag = mu_melt_visc;
      
    	/* Temperature dependent viscosity of the gas */ 
		if (temp <= 303.15) 
		{         		
       			mu_gas = 17.15e-6;
    		} 
		else if (temp <= 373.15) 
		{
        		mu_gas = 17.15e-6 + (temp - 303.15) * (21.74e-6 - 17.15e-6) / (373.15 - 303.15);
    		} 
		else if (temp <= 673.15) 
		{
        		mu_gas = 21.74e-6 + (temp - 373.15) * (32.5e-6 - 21.74e-6) / (673.15 - 373.15);
    		} 
		else if (temp <= 973.15) 
		{
       		 	mu_gas = 32.5e-6 + (temp - 673.15) * (40.9e-6 - 32.5e-6) / (973.15 - 673.15);
   		} 
		else if (temp <= 1273.15) 
		{
        		mu_gas = 40.9e-6 + (temp - 973.15) * (47.9e-6 - 40.9e-6) / (1273.15 - 973.15);
    		} 
		else if (temp <= 1573.15) 
		{
        		mu_gas = 47.9e-6 + (temp - 1273.15) * (54.0e-6 - 47.9e-6) / (1573.15 - 1273.15);
    		} 
		else if (temp <= 1873.15) 
		{
        		mu_gas = 54.0e-6 + (temp - 1573.15) * (59.6e-6 - 54.0e-6) / (1873.15 - 1573.15);
    		} 
		else 
		{
        		// Extend linearly above 1873.15 K
        		mu_gas = 59.6e-6 + (temp - 1873.15) * (59.6e-6 - 54.0e-6) / (1873.15 - 1573.15);
    		}

    	/* Geometric averaging of the properties, */
		if(alpha <= 1 && alpha >= 0)
			mu_mix = pow(mu_slag, alpha) * pow(mu_gas, 1.0 - alpha);
		else
			mu_mix = 17.15e-6;

    	return mu_mix;
}

DEFINE_PROPERTY(viscosity_harm_avg, c, t)
{	

    	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
    	Thread *mphase_t = THREAD_SUB_THREAD(tmix, 1);
    	Thread *gphase_t = THREAD_SUB_THREAD(tmix, 0);

	real mu_melt_visc = mu_melt;
	real mu_crit_visc = mu_crit;

	real T_melt_visc = T_melt;
	real T_crit_visc = T_crit;

	real mu_gas = 1e-5;
    	real mu_slag = mu_melt_visc;

    	real k;
	real d;

    	real mu_mix;
	real alpha;

    	real temp = C_T(c, tmix);
	alpha = C_VOF(c, mphase_t);

	/* Temperature dependent viscosity of the slag */ 

    		k = (mu_melt_visc - mu_crit_visc) / (T_melt_visc - T_crit_visc);
    		d = mu_crit_visc - k * T_crit_visc;

		if (temp > 0 && temp < T_crit_visc)
   			mu_slag = mu_crit_visc;
		else if (temp >= T_crit_visc && temp < T_melt_visc)
   		 	mu_slag = k * temp + d;
		else
    			mu_slag = mu_melt_visc;
      
    	/* Temperature dependent viscosity of the gas */ 
		if (temp <= 303.15) 
		{         		
       			mu_gas = 17.15e-6;
    		} 
		else if (temp <= 373.15) 
		{
        		mu_gas = 17.15e-6 + (temp - 303.15) * (21.74e-6 - 17.15e-6) / (373.15 - 303.15);
    		} 
		else if (temp <= 673.15) 
		{
        		mu_gas = 21.74e-6 + (temp - 373.15) * (32.5e-6 - 21.74e-6) / (673.15 - 373.15);
    		} 
		else if (temp <= 973.15) 
		{
       		 	mu_gas = 32.5e-6 + (temp - 673.15) * (40.9e-6 - 32.5e-6) / (973.15 - 673.15);
   		} 
		else if (temp <= 1273.15) 
		{
        		mu_gas = 40.9e-6 + (temp - 973.15) * (47.9e-6 - 40.9e-6) / (1273.15 - 973.15);
    		} 
		else if (temp <= 1573.15) 
		{
        		mu_gas = 47.9e-6 + (temp - 1273.15) * (54.0e-6 - 47.9e-6) / (1573.15 - 1273.15);
    		} 
		else if (temp <= 1873.15) 
		{
        		mu_gas = 54.0e-6 + (temp - 1573.15) * (59.6e-6 - 54.0e-6) / (1873.15 - 1573.15);
    		} 
		else 
		{
        		// Extend linearly above 1873.15 K
        		mu_gas = 59.6e-6 + (temp - 1873.15) * (59.6e-6 - 54.0e-6) / (1873.15 - 1573.15);
    		}

    	/* Harmonic averaging of the properties, */
			mu_mix = mu_slag*mu_gas/(mu_gas*alpha+(1-alpha)*mu_slag);

    	return mu_mix;
}

DEFINE_PROPERTY(viscosity_arith_avg, c, t)
{	

    	Thread *tmix = (NNULLP(THREAD_SUPER_THREAD(t))) ? THREAD_SUPER_THREAD(t) : t; /* mixture thread */
    	Thread *mphase_t = THREAD_SUB_THREAD(tmix, 1);
    	Thread *gphase_t = THREAD_SUB_THREAD(tmix, 0);

	real mu_melt_visc = mu_melt;
	real mu_crit_visc = mu_crit;

	real T_melt_visc = T_melt;
	real T_crit_visc = T_crit;

	real mu_gas = 1e-5;
    	real mu_slag = mu_melt_visc;

    	real k;
	real d;

    	real mu_mix;
	real alpha;

    	real temp = C_T(c, tmix);
	alpha = C_VOF(c, mphase_t);

	/* Temperature dependent viscosity of the slag */ 

    		k = (mu_melt_visc - mu_crit_visc) / (T_melt_visc - T_crit_visc);
    		d = mu_crit_visc - k * T_crit_visc;

		if (temp > 0 && temp < T_crit_visc)
   			mu_slag = mu_crit_visc;
		else if (temp >= T_crit_visc && temp < T_melt_visc)
   		 	mu_slag = k * temp + d;
		else
    			mu_slag = mu_melt_visc;
      
    	/* Temperature dependent viscosity of the gas */ 
		if (temp <= 303.15) 
		{         		
       			mu_gas = 17.15e-6;
    		} 
		else if (temp <= 373.15) 
		{
        		mu_gas = 17.15e-6 + (temp - 303.15) * (21.74e-6 - 17.15e-6) / (373.15 - 303.15);
    		} 
		else if (temp <= 673.15) 
		{
        		mu_gas = 21.74e-6 + (temp - 373.15) * (32.5e-6 - 21.74e-6) / (673.15 - 373.15);
    		} 
		else if (temp <= 973.15) 
		{
       		 	mu_gas = 32.5e-6 + (temp - 673.15) * (40.9e-6 - 32.5e-6) / (973.15 - 673.15);
   		} 
		else if (temp <= 1273.15) 
		{
        		mu_gas = 40.9e-6 + (temp - 973.15) * (47.9e-6 - 40.9e-6) / (1273.15 - 973.15);
    		} 
		else if (temp <= 1573.15) 
		{
        		mu_gas = 47.9e-6 + (temp - 1273.15) * (54.0e-6 - 47.9e-6) / (1573.15 - 1273.15);
    		} 
		else if (temp <= 1873.15) 
		{
        		mu_gas = 54.0e-6 + (temp - 1573.15) * (59.6e-6 - 54.0e-6) / (1873.15 - 1573.15);
    		} 
		else 
		{
        		// Extend linearly above 1873.15 K
        		mu_gas = 59.6e-6 + (temp - 1873.15) * (59.6e-6 - 54.0e-6) / (1873.15 - 1573.15);
    		}

    	/* Harmonic averaging of the properties, */
			mu_mix = mu_slag*alpha+mu_gas*(1-alpha);

    	return mu_mix;
}

/* Sludge Data */
#define Ns 21
/* Temperature in K */
real sludge_temp_data[Ns] = 		{293.15, 323.15, 423.15, 523.15, 623.15, 723.15, 823.15, 923.15, 1023.15, 1123.15, 1223.15, 1323.15, 1373.15, 1423.15, 1523.15, 1573.15, 1623.15, 1723.15, 1823.15, 1873.15, 2773.15};

/* Specific heat including melting enthalpy in J/(kg*K), data from INSPYRO */
real sludge_specific_heat_data[Ns] = 	{750, 770, 866, 959, 1027, 1118, 1140, 1170, 1226, 1260, 1290, 1362, 3308, 3303, 3295, 3291, 1353, 1331, 1327, 1325, 1325};

/* Sewage sludge particle enthalpy */
real sludge_enthalpy_data[Ns] = 	{0, 22.8e3, 104.6e3, 195.9e3, 295.2e3, 402.4e3, 515.3e3, 630.8e3, 750.6e3, 874.9e3, 1002.4e3, 1135.0e3, 1251.8e3, 1417.0e3, 1746.9e3, 1991.6e3, 2027.7e3, 2161.9e3, 2294.8e3, 2361.1e3, 3553.6e3};

DEFINE_DPM_PROPERTY(sludge_cp, c, t, tp, T)
{
	real cp;
    	real enthalpy_sludge = 0;
    	int i;


    	if (NULLP(tp->pp) || NULLP(TP_CELL_THREAD(tp)))
    	{
        	cp = 1680;
        	return cp;
    	}

    	/* Handle cases where temperature is outside the data range */
    	if (T <= sludge_temp_data[0]) 
    	{
        	cp = sludge_specific_heat_data[0];
        	enthalpy_sludge = sludge_enthalpy_data[0];
    	} 
    	else if (T >= sludge_temp_data[Ns-1]) 
    	{
        	cp = sludge_specific_heat_data[Ns-1];
        	enthalpy_sludge = sludge_enthalpy_data[Ns-1];
    	} 
    	else 
    	{
        /* Perform piecewise linear interpolation */
        	for (i = 0; i < Ns-1; i++) 
        	{
            		if (T >= sludge_temp_data[i] && T < sludge_temp_data[i+1]) 
            		{
                		cp = sludge_specific_heat_data[i] + (sludge_specific_heat_data[i+1] - sludge_specific_heat_data[i]) * 
                		(T - sludge_temp_data[i]) / (sludge_temp_data[i+1] - sludge_temp_data[i]);
                
                		enthalpy_sludge = sludge_enthalpy_data[i] + (sludge_enthalpy_data[i+1] - sludge_enthalpy_data[i]) * 
                		(T - sludge_temp_data[i]) / (sludge_temp_data[i+1] - sludge_temp_data[i]);
            		}
        	}
    	}

    	tp->enthalpy = enthalpy_sludge;
    	return cp;
}

/* Limestone Data */
#define Nls 7

/* Temperature in K */
real limestone_temp_data[Nls] = 		{293.15, 1123.15, 1223.15, 1473.15, 1673.15, 1873.15, 2473.15};

/* Specific heat including calcination enthalpy in J/(kg*K) */
real limestone_specific_heat_data[Nls] = 	{856, 1200, 3650, 3650, 800, 750, 750};

/* Limestone particle enthalpy */
real limestone_enthalpy_data[Nls] = 		{0, 853.2e3, 1095.7e3, 2008.2e3, 2453.2e3, 2608.2e3, 3058.2e3};

DEFINE_DPM_PROPERTY(limestone_cp, c, t, tp, T)
{
	real cp;
	real rho0 = 3000;			/* Density of limestone particle at room temperature */	 
	real rho;				/* Density of limestone particle that is assumed to decrease due to calcination */
    	real enthalpy_limestone = 0;
    	int i;

   	if (NULLP(tp->pp) || NULLP(TP_CELL_THREAD(tp)))
    	{
        	cp = 1680;
        	return cp;
    	}


    	/* Handle cases where temperature is outside the data range */
    	if (T <= limestone_temp_data[0]) 
    	{
        	cp = limestone_specific_heat_data[0];
        	enthalpy_limestone = limestone_enthalpy_data[0];
   	} 
    	else if (T >= limestone_temp_data[Nls-1]) 
    	{
        	cp = limestone_specific_heat_data[Nls-1];
        	enthalpy_limestone = limestone_enthalpy_data[Nls-1];
    	} 
    	else 
    	{
        	/* Perform piecewise linear interpolation */
        	for (i = 0; i < Nls-1; i++) 
        	{
            		if (T >= limestone_temp_data[i] && T < limestone_temp_data[i+1]) 
            		{
                		cp = limestone_specific_heat_data[i] + (limestone_specific_heat_data[i+1] - limestone_specific_heat_data[i]) * 
                		(T - limestone_temp_data[i]) / (limestone_temp_data[i+1] - limestone_temp_data[i]);

                		enthalpy_limestone = limestone_enthalpy_data[i] + (limestone_enthalpy_data[i+1] - limestone_enthalpy_data[i]) * 
                		(T - limestone_temp_data[i]) / (limestone_temp_data[i+1] - limestone_temp_data[i]);
            		}
        	}
    	}

	rho = TP_RHO(tp);
    	tp->enthalpy = enthalpy_limestone*rho0/TP_RHO(tp);
    	return cp*rho0/TP_RHO(tp);
}


DEFINE_ON_DEMAND(kill_phase_burner)
  {
	Domain *domain = Get_Domain(ROOT_DOMAIN_ID);
	Thread *t;
	cell_t c;

	int fluid_ID = 146;
	Thread *tf = Lookup_Thread(domain,fluid_ID);
	Thread *mphase_tf = THREAD_SUB_THREAD(tf , 1);

	thread_loop_c(t,domain) 
	{
			begin_c_loop(c,tf) 
			{
				C_VOF(c,mphase_tf) = 0;
			}
	}
	end_c_loop(c,tf)
	Message0("\nSlag phase fraction in burner zone reset.");
	
  }

DEFINE_DPM_DRAG(modifieddraglaw_nonspherical, Re, tp)
{
	real Cd, drag_force = 0.0;
	real b1, b2, b3, b4;
   	real corr = 1.0;  
 	real Rec;

	/* phi = sphericity (wadell): 0.806 ~ cubic */
	real phi = 0.806;

    	/* Cell Thread check */
    	Thread *t = TP_CELL_THREAD(tp);
    	cell_t c = TP_CELL(tp);
    	Thread *mphase_t = THREAD_SUB_THREAD(t, 1); 
    	Thread *gphase_t = THREAD_SUB_THREAD(t, 0); 


	if(C_VOF(c,mphase_t) < CUTOFF)
	{
		/* Modify Re number */
		/* Re = rho *d_p* abs(u_p -u)/mu */
 		Rec = Re * C_R(c, gphase_t) / C_R(c, t) * C_MU_L(c, t) / C_MU_L(c, gphase_t);
 		Re = Rec;
 		corr = C_MU_L(c, gphase_t) / C_MU_L(c, t);
	}	
	
	/*Drag corrleation Haider-Levenspiel*/

		b1 = exp(2.3288 - 6.4581 * phi + 2.4486 * pow(phi, 2));
		b2 = 0.0964 + 0.5565 * phi;
		b3 = exp(4.905 - 13.8944 * phi + 18.4222 * pow(phi, 2) - 10.2599 * pow(phi, 3));
		b4 = exp(1.4681 + 12.254 * phi - 20.7322 * pow(phi, 2) + 15.8855 * pow(phi, 3));

		Cd = (24 / Re) * (1 + b1 * pow(Re, b2)) + (b3 * Re / (b4 + Re));

    		/* Corrected drag force calculation */

		if (TP_USER_REAL(tp,9) == 0)
		{
    			drag_force = 18.0 * Cd * Re / 24.0 * corr;
		}
		else
		{
			drag_force = 0;
		}
		return drag_force;
}