/**********************************************************************
   UDF to calculate standard deviation and mean of the Monroe Parameter.             
***********************************************************************/
#include "udf.h"

DEFINE_ON_DEMAND(on_demand2_calc)
{
   Domain *d;

   real velocity;
   Thread *t;
   cell_t c;
   Node *nod;
   int n;
   int j;
   int density=1000;
   double flow=.01;
   double fh=1;
   double bl_thick=.01;
   double visc=1e-6;
   double pos_x;
   double pos_y;
   double xvel;
   double yvel;
   double xprev=0;
   double yprev=0;
   double MeshArea=0;
   double cell_area=0;
   double magvel=0;
   double flowtot=0;
   double baffle[7] = {0, .1, .2, .3, .4, .5, .6};
   double Monroe_sum1[6]={0, 0, 0, 0, 0, 0};
   double Monroe_sum2[6]={0, 0, 0, 0, 0, 0};
   double Monroe_sum3[6]={0, 0, 0, 0, 0, 0};
   double Manroe[6] = {0, 0, 0, 0, 0, 0};
   double display2 =0;
   double display1 =0;
   int counter=0;
      int q=0;
   
   d = Get_Domain(1);     /* Get the domain using Fluent utility */

   /* Loop over all cell threads in the domain */
   for (q = 0; q <=5; q++)
   {
	thread_loop_c(t,d)
		{
		begin_c_loop(c,t)
			{
		    double delta_x = 0;
			double xpos_ave = 0;
			double delta_y = 0;
			double ypos_ave =0;
			int run_number = 0;
		    for (n=0; n < cell_type_nnodes[(int)C_TYPE(c,t)]; n++) 	 
				{
					nod = C_NODE(c, t, n); 
					pos_x=NODE_X(nod);
					pos_y=NODE_Y(nod);
					
				if(run_number>0)
				{
					if(delta_x<fabs(pos_x-xprev))
						delta_x= fabs(pos_x-xprev);
						xpos_ave=pos_x*.5+xprev*.5;
					
					if(delta_y < fabs(pos_y-yprev))
						delta_y = fabs(pos_y-yprev);
						ypos_ave=pos_y*.5+yprev*.5;
					
				}
				run_number = run_number +1;
				xprev = pos_x;
				yprev = pos_y;

			}
			if(xpos_ave>(baffle[q]+bl_thick) & xpos_ave<(baffle[q+1]-bl_thick) & ypos_ave>bl_thick & ypos_ave<(fh-bl_thick))
			{	Manroe[q]=Manroe[q]+delta_x*delta_y;
				Monroe_sum1[q] = Monroe_sum1[q] + pow(C_D(c,t),.333)*delta_x*delta_y;
				Monroe_sum2[q] = Monroe_sum2[q] + pow(C_D(c,t),0.5)*delta_x*delta_y;
				Monroe_sum3[q] = Monroe_sum3[q] + C_D(c,t)*delta_x*delta_y;
			}
		}
     		end_c_loop(c,t);
	}
	
		
   }
	printf("Eps_onethird_pheta %g\n");	
   for (q=0; q<=5; q++)
	{printf("%g\n",Monroe_sum1[q]/flow);}
   printf("Eps_onehalf_pheta = %g\n");
   for (q=0; q<=5; q++)
   {printf("%g\n",Monroe_sum2[q]/flow);}
   printf("Eps_one_pheta = %g\n");
   for (q=0; q<=5; q++)
	printf("%g\n",Monroe_sum3[q]/flow);
   for (q=0; q<=5; q++)
   printf("Area = %g\n", Manroe[q]);
}