#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "surfactant-transport.h"
#include "view.h"


double R0 = 0.2;
double K = 1.;
double Pe = 16.;
int LEVEL = 6;
int main (int argc, char * argv[])
{

  origin (-L0/2., -L0/2.);

   //Redistance at each comp step
  reinit_skip_steps = 1; 
  //If redistancing not done every step then use transport equation for phase-field
  advect_diff_phase_field = 1; 
  D_s = sq(2.*R0)*K/Pe;
  // Iterate over 3 resolutions
  for (LEVEL = 5; LEVEL <= 7; LEVEL++) {
    N = 1 << LEVEL;
    counter = 0;
    run();
  }

}

event init (i = 0)
{
     refine( (sq(0.3)- pow(x,2) - pow(y,2)) > 0 && level < LEVEL);

    fraction(f, -(sq(0.2)- pow(x,2) - pow(y,2)) );
    event ("properties2");
  
 // Initialize surfactant 
  foreach(){ 
  //  d2[] = (0.2) - sqrt(pow(x,2) + pow(y,2)) ;
  //  pfield[] = 0.5*(1. - tanh((d2[])/2./EPSILON));
  //  pfield[] = clamp(pfield[], 0., 1.);
    double deltas = (pfield[]*(1. - pfield[]))/EPSILON;
    double gamma0 = 1.;
    c1[] = gamma0*deltas;
  }

  if (grid->maxdepth == 7){
    //  view (fov = 30, near = 0.01, far = 1000,
	  // tx = 0.009, ty = -0.076, tz = -0.291,
	  // width = 1239, height = 575);
    draw_vof (c = "f", lw = 2);
    squares (color = "c1", spread = -1, linear = true);
    vectors (u = "u", scale = 1);
    
    save ("fields-initial.png");
  }
 
}

// Radial Velocity 
event stability (i++) {

	trash ({u,uf});

	foreach_face(y) uf.x[] = K*(sqrt(x*x + y*y))*cos(atan2(y,x));
	foreach_face(x) uf.y[] = K*(sqrt(x*x + y*y))*sin(atan2(y,x));

	foreach(){
		u.x[] = K*(sqrt(x*x + y*y))*cos(atan2(y,x));
		u.y[] = K*(sqrt(x*x + y*y))*sin(atan2(y,x));
	}

}

/*
event adapt (i++) {

	scalar pfield2[];
	foreach(){
     pfield2[] = 0.5*(1. - tanh((d2[])/4./EPSILON));
     //c2[] = clamp(pfield[], 0., 1.);	
 }

 double uxemax = 0.1*statsf(u.x).max;
 double uyemax = 0.1*statsf(u.y).max;
// double uzemax = 0.1*statsf(u.z).max;
 double c1emax = 0.001*statsf(c1).max;
 double pfemax = 0.001*statsf(pfield).max;
 double pf2emax = 0.001*statsf(pfield2).max;
 double d2emax = 0.1*statsf(d2).max;

  adapt_wavelet ({c1,pfield2}, (double[]){c1emax,pf2emax}, LEVEL); 
//  printf("\n%d %g",i,t); fflush(stdout);
}
*/

int count = 0;
// event logfile(i += 1){
event logfile(t = 0; t += 0.01; t <= 0.7){
double cnet = 0.; //net surfactant
scalar deltas[];
  foreach(){
   // pfield[] = 0.5*(1. - tanh((d2[])/2./EPSILON));
   // pfield[] = clamp(pfield[], 0., 1.);
    deltas[] = pfield[]*(1. - pfield[])/EPSILON;
    cnet += c1[]*dv();
    gamma2[] = c1[]*4.*EPSILON;
  }
  char filename[200];
  sprintf(filename,"expanding-circle-test-%d",(1 << grid->maxdepth));
  static FILE * fp2 = fopen(filename,"w");
  //double theta = M_PI/2.;
  fprintf (fp2, "%.17g %.17g %.17g\n", t, statsf(c1).max/statsf(deltas).max, cnet);
  fflush (fp2);

	
  count++;

}

event end( t = 0.7){

  if (grid->maxdepth == 7){
    //  view (fov = 30, near = 0.01, far = 1000,
	  // tx = 0.009, ty = -0.076, tz = -0.291,
	  // width = 1239, height = 575);
    draw_vof (c = "f", lw = 2);
    squares (color = "gamma2", spread = -1., min = 0, max = 1, linear = true);
    // vectors (u = "u", scale = 1);
    
    save ("fields-end.png");
  }
  // save();
}