/**********************************************************************/
/* unsteady.c                                                         */
/* UDF for specifying a Mechanical Ventilation waveform */
/*Raghu Arambakam*/
/**********************************************************************/
 
#include "udf.h"
#define tidalVol 420e-6 /*tidal volume of lung m^3*/
/*#define flowRate 0.001 /*tidal volume of lung m^3/s i.e. 0.0005 m^3/s = 30 l/min*/
real counter=0.0;
 
DEFINE_PROFILE(unsteady_velocity, thread, position)
{
  face_t f;
  real crossSect=0.0;
  real inhaleTime = 0.4;       /*time constant s*/
  real exhaleTime = 1.6;       /*time constant s*/
  real flowRate=tidalVol/inhaleTime; /*flow rate*/
  real taut3=0.0;
  real tau = tidalVol/flowRate;       /*time constant s*/
  real generation=15;
  real NV_VEC(farea);

  begin_f_loop(f, thread) /*compute area of inlet*/
   {
    F_AREA(farea,f,thread);   /* computes area of each face  */
    crossSect +=NV_MAG(farea);  /* computes total face area by accumulating  magnitude of each face's area  Units defined in FLUENT*/
   }
  end_f_loop(f, thread)
   
  Message ("inlet area %g \n", crossSect);

  inhaleTime=0.4;
  exhaleTime=1.6; /*2 for 30 l/m and 4 for 60 l/m according to Koombua MCB 2009*/
  taut3= (inhaleTime+exhaleTime);
  /*tau= inhaleTime;*/
  begin_f_loop(f, thread) /*assign velocity*/
    {
      real t = RP_Get_Real("flow-time");
      	if (t<=inhaleTime) /*inhalation i.e. positive velocity*/
         {
          F_PROFILE(f, thread, position) = flowRate/(crossSect*pow(2,generation-1));
         }
        else /*exhalatioin negative velocity*/
         {
          taut3=t-inhaleTime;
          F_PROFILE(f, thread, position) = -(tidalVol*exp(-(taut3)/tau))/(0.3926737*crossSect*tau*pow(2,generation-1));
         }
     }
  end_f_loop(f, thread)
}