DEFINE_CG_MOTION(foilright,dt,vel,omega,time,dtime)
{
   real x_A=0.0,y_A=0.0,x_N=0.0,y_N=0.0;
   real theta_A,theta_N;
   real omegap = 0.0;
   real presure;
   real x[ND_ND];
/////////////////////////////////////////////////////////////////////////////
    face_t f ;
    cell_t c ;
        real f_glob[ND_ND],m_glob[ND_ND],x_cg[ND_ND];
    real f1,f3,f5;
    Domain *domain= Get_Domain (1);
    Thread *tf1 = Lookup_Thread (domain, 33);
/////////////////////////////////////////////////////////////////////////////
   theta_A = -theta0*cos(w*CURRENT_TIME);
   theta_N = -theta0*cos(w*(CURRENT_TIME+CURRENT_TIMESTEP));
   omegap = (theta_N-theta_A)/CURRENT_TIMESTEP;

   x_A= Radius*cos(w*CURRENT_TIME);
   x_N= Radius*cos(w*(CURRENT_TIME+CURRENT_TIMESTEP));
   y_A= Radius*sin(w*CURRENT_TIME);
   y_N= Radius*sin(w*(CURRENT_TIME+CURRENT_TIMESTEP));
  
   vel[0] =(x_N - x_A) / CURRENT_TIMESTEP;
   vel[1] =(y_N - y_A) / CURRENT_TIMESTEP;
   omega[0] =omega[1] =0.0;
   omega[2] = omegap;
   // f1  f3  f5 output /////////////////////////////////////////////////////
   if(time>=dtime){
   x_cg[0]=DT_CG(dt)[0];
    x_cg[1]=DT_CG(dt)[1];
    Compute_Force_And_Moment (domain, tf1, x_cg, f_glob, m_glob,1);
    f1 = f_glob[0]/(0.5*rou*cord*U*U);
    f3 = f_glob[1]/(0.5*rou*cord*U*U);
    f5 = m_glob[2]/(0.25*rou*cord*cord*U*U);
    fpright=fopen("foil_right.dat", "a+");
    fprintf(fpright,"%.6f %.5f %.5f %.5f\n",CURRENT_TIME*w/(2*PI),f1,f3,f5);
        fclose(fpright);}
   //////////////////////////////////////////////////////////////////////////
  return;
}