/**********************************************************************
   UDF to Look at the distribution of the energy dissipation rate             
***********************************************************************/
#include "udf.h"

DEFINE_ON_DEMAND(on_demand_calc)
{
   Domain *d;

   real velocity;
   Thread *t;
   cell_t c;
   Node *nod;
   int n;
   double pos_x;
   double pos_y;
   double xvel;
   double yvel;
   double xprev=0;
   double yprev=0;
   double MeshArea=0;
   double cell_area=0;
      
   //bin analysis
   double checkmax = 0; // check max energy dissipation
   int binlength = 20;
   int binsize = 20;
   double bin[20]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
   int i = 0;
   int j = 0;
   int nondimensionE = 1; // should be 100 if normalize
   double area[20]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
   double start = 0.001;
   double end = 0.01;
  

   d = Get_Domain(1);     /* Get the domain using Fluent utility */

   /* Loop over all cell threads in the domain */
   thread_loop_c(t,d)
     {

     begin_c_loop(c,t)
       {
		double delta_x = 0;
		double delta_y = 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);
					
					if(delta_y < fabs(pos_y-yprev))
						delta_y = fabs(pos_y-yprev);
					
				}
				run_number = run_number +1;
				xprev = pos_x;
				yprev = pos_y;

				}
				MeshArea = MeshArea+delta_x*delta_y;	
        }
     end_c_loop(c,t)
     printf("\nMeshArea = %g\n", MeshArea);
     /* Loop over all cells  */
     begin_c_loop(c,t)
       {
		    double delta_x = 0;
			double delta_y = 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);
					
					if(delta_y < fabs(pos_y-yprev))
						delta_y = fabs(pos_y-yprev);
					
				}
				run_number = run_number +1;
				xprev = pos_x;
				yprev = pos_y;

				}
			


       }
    
	 end_c_loop(c,t)	


		
			checkmax = end - start;
		printf("\nmax e = %g\n\n", checkmax);

		for (i = 0; i<=binsize-1; i++)
		{
			bin[i] = start + checkmax/binlength*i;
		}

		//chaterizing bins by looping through bins
		
		for (i=0; i<=binsize-1; i++)
		{
			//printf("\n bin = %g\n", bin[i]);
		
		begin_c_loop(c, t)
  		{ 		
			if(bin[i]<C_D(c,t) && C_D(c,t)<bin[i+1] )
			{
				double delta_x = 0;
				double delta_y = 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);
					
					if(delta_y < fabs(pos_y-yprev))
						delta_y = fabs(pos_y-yprev);
					
				}
				run_number = run_number +1;
				xprev = pos_x;
				yprev = pos_y;
				
				area[i] =  delta_x*delta_y+area[i];


				}


			}

			


		}
		end_c_loop(c,t)	
	
			area[i]=area[i]/2/MeshArea;
			//if (i<=19)
			//{printf("\n area = %g\n", area[i]);}
		}

printf("bin\n");
for (j=0; j<=binsize-1; j++)
{printf(" %g\n", bin[j]);}
printf("\n area\n");
for (j=0; j<=binsize-1; j++)
{printf(" %g\n", area[j]);}

		
	}

}