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


   //Boundary Conditions
c1[top] = dirichlet(0.);
c1[bottom] = dirichlet(0.);
c1[left] = dirichlet(0.);
c1[right] = dirichlet(0.);

   //Boundary Conditions
pfield[top] = dirichlet(0.);
pfield[bottom] = dirichlet(0.);
pfield[left] = dirichlet(0.);
pfield[right] = dirichlet(0.);

int counter2 = 0;
double db = 0.4 [*];
double Pe = 160. [*];
double omega = 10. [*];

int main (int argc, char * argv[])
{

   //Redistance at each comp step
  reinit_skip_steps = 10; 
  //If redistancing not done every step then use transport equation for phase-field
  advect_diff_phase_field = 1; 

  //Peclet number
  
  D_s = sq(db)*omega/Pe;

  //Do not advect velocity (since it is reinitialised at each timestep)
  stokes = true;
  
  //Use strict minmod limiting when advecting c1
  theta = 1.;
  c1.gradient = minmod2;
 
  // Iterate over 3 resolutions
  for (int LEVEL = 5; LEVEL <= 7; LEVEL++) {

    N = 1 << LEVEL;
    origin (-L0/2. + L0/N/2., -L0/2 + L0/N/2.);
    counter2 = 0;
    counter = 0;
    run();
  }

}

event init (i = 0)
{
    fraction(f, -(sq(db/2.)- pow(x,2) - pow(y,2)) );
    event ("properties2");

  
 // Thickness of the phase-field over which surfactant is distributed
  // Initialize surfactant 
  foreach(){ 
    d2[] = (db/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 = 2. + sin(atan2(y,x));
    c1[] = gamma0*deltas;
    gamma2[] = c1[]*4.*EPSILON;
  }
  event("stability");
  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", min = 0, max = 3, spread = -1., linear = true);
  
    
    save ("rotating_fields_initial.png");
  }

}

// Radial Velocity 
event stability (i++) {
  foreach_face (x)
    uf.x[] = omega*y;
  foreach_face (y)
    uf.y[] = - omega*x;
  foreach(){
    u.x[] =   omega*y;
    u.y[] = - omega*x;
  }
}

// event logfile(i += 1){
event logfile(t = 0; t += 0.01; t <= 2.){
  double cnet = 0;
  scalar deltas[];
  foreach() {
    deltas[] = pfield[]*(1. - pfield[])/EPSILON;
    cnet += c1[]*dv();
  }
  
  char filename[200];
  sprintf(filename,"rotating-circle-test-%d",(1 << grid->maxdepth));
  static FILE * fp2 = fopen(filename,"w");
  double theta = M_PI/2.;
  double c_piby2 = interpolate(c1, db/2.*cos(theta), db/2.*sin(theta), 0.);
  double deltas_piby2 = interpolate(deltas, db/2.*cos(theta), db/2.*sin(theta), 0.);
  fprintf (fp2, "%.17g %.17g %.17g\n", t, c_piby2/deltas_piby2,cnet);
  //  fprintf (fp2, "%.17g %.17g %.17g\n", t, gamma2_theta,cnet);
  fflush (fp2);



  counter2++;
//  exit(0);

}

event end( t = 2.){

  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", min = 0, max = 3, spread = -1., linear = true);
    // vectors (u = "u", scale = 1);
    save ("rotating_fields_end.png");
  }
  // save();
}