#include "udf.h"

#define WALL_ID 4  // the wall id to be computed collision angle
#define PARTICLE_PHASE_ID 1 // the phase id of particle phase

DEFINE_ON_DEMAND(report_angle)
{
        face_t f;
        cell_t c0;
        real cos_angle,angle;
        real NV_VEC(f_normal),NV_VEC(particle_velocity);
        Domain*mixture_domain,*particle_domain;
        Thread*wall_thread,*particle_t0;

        mixture_domain=Get_Domain(1);
        particle_domain=DOMAIN_SUB_DOMAIN(mixture_domain,PARTICLE_PHASE_ID);
        wall_thread=Lookup_Thread(particle_domain,WALL_ID);
        begin_f_loop(f,wall_thread)
        {
                c0=F_C0(f,wall_thread);
                particle_t0=THREAD_T0(wall_thread);
#if RP_3D
                NV_D(particle_velocity,=,C_U(c0,particle_t0),C_V(c0,particle_t0),C_W(c0,particle_t0));
#else
                NV_D(particle_velocity,=,C_U(c0,particle_t0),C_V(c0,particle_t0),0);
#endif
                F_AREA(f_normal,f,wall_thread);
                cos_angle=NV_DOT(particle_velocity,f_normal)/MAX(NV_MAG(particle_velocity)*NV_MAG(f_normal),1E-10);
                angle=90.0-acos(cos_angle)/M_PI*180.0;
                Message0("face id=%d,angle between particle velocity and wall face=%f degree\n",f,angle);
        }
        end_f_loop(f,wall_thread);
}