#include "tag.h"
#define dA (0.5 * sq(Delta))

trace
void vorticity3d(vector u, vector omega){
  vector du[], dv[], dw[]; // Would be nice to use a tensor here.
  foreach(){
    du.x[] = (u.x[1] - u.x[-1])/(2.*Delta);
    dv.x[] = (u.y[1] - u.y[-1])/(2.*Delta);
    dw.x[] = (u.z[1] - u.z[-1])/(2.*Delta);

    du.y[] = (u.x[0,1] - u.x[0,-1])/(2.*Delta);
    dv.y[] = (u.y[0,1] - u.y[0,-1])/(2.*Delta);
    dw.y[] = (u.z[0,1] - u.z[0,-1])/(2.*Delta);

    du.z[] = (u.x[0,0,1] - u.x[0,0,-1])/(2.*Delta);
    dv.z[] = (u.y[0,0,1] - u.y[0,0,-1])/(2.*Delta);
    dw.z[] = (u.z[0,0,1] - u.z[0,0,-1])/(2.*Delta);
  }

  foreach(){
    omega.x[] = dw.y[] - dv.z[];
    omega.y[] = du.z[] - dw.x[];
    omega.z[] = dv.x[] - du.y[];
  }
  boundary ((scalar*){omega});
}

trace
void vorticity_moments_plane(scalar omega, scalar m, int nm, char * name, coord n, double _alpha){

  double circulation, M20, M02, M11, va, vb, vc, ve, omax;
  coord centroids={0.,0.,0.};

  for (int nj = 0; nj < nm; nj++){
    // Compute the circulation $\Gamma = \iint \omega dA$ and position of the
    // vortex centroids $\vec{\mu} = \iint \vec{x}\omega dA$
    double var1 = 0., var2 = 0., var3 = 0., var4 = 0.;
    foreach(reduction(+:var1), reduction(+:var2), reduction(+:var3), reduction(+:var4)){
      double alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;
      if (fabs(alpha) > 0.87)
        continue;
      if (is_active(cell) && is_local(cell)) {
        coord p = {x,y,z};

        double sval;
        if (n.x == 1)
          sval = 0.5*(val(omega) + val(omega,1,0,0));
        else if (n.y == 1)
          sval = 0.5*(val(omega) + val(omega,0,1,0));
        else
          sval = 0.5*(val(omega) + val(omega,0,0,1));

        if (m[] == nj+1){
          var1 += dA * sval;
          var2 += dA * sval * p.x;
          var3 += dA * sval * p.y;
          var4 += dA * sval * p.z;
        }
      }
    }
    circulation = var1;


    if (circulation != 0) {
      centroids = (coord) {var2/var1, var3/var1, var4/var1};

      // For the ellipticity we require the centered vorticity moments
      // $ M_{p,q} = \iint (x - \mu_x)^p (y - \mu_y)^q \omega dA $

      var1 = 0., var2 = 0., var3 = 0., var4 = 0.;
      double var5 = 0., var6 = 0.;

      foreach(reduction(+:var1), reduction(+:var2), reduction(+:var3), reduction(+:var4), reduction(+:var5), reduction(+:var6)){
        double alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;
        if (fabs(alpha) > 0.87)
          continue;
        if (is_active(cell) && is_local(cell)) {
          coord p = {x,y,z};

          double sval;
          if (n.x == 1)
            sval = 0.5*(val(omega) + val(omega,1,0,0));
          else if (n.y == 1)
            sval = 0.5*(val(omega) + val(omega,0,1,0));
          else
            sval = 0.5*(val(omega) + val(omega,0,0,1));

          if (m[] == nj+1){
            var1 += dA * sval * sq(p.x - centroids.x);
            var2 += dA * sval * sq(p.y - centroids.y);
            var3 += dA * sval * sq(p.z - centroids.z);
            var4 += dA * sval * (p.x - centroids.x) * (p.y - centroids.y);
            var5 += dA * sval * (p.y - centroids.y) * (p.z - centroids.z);
            var6 += dA * sval * (p.z - centroids.z) * (p.x - centroids.x);
          }
        }
      }

      if (n.x == 1){
        M20 = var2/circulation;
        M02 = var3/circulation;
        M11 = var5/circulation;
      }
      else if (n.y == 1){
        M20 = var1/circulation;
        M02 = var3/circulation;
        M11 = var6/circulation;
      }
      else{
        M20 = var1/circulation;
        M02 = var2/circulation;
        M11 = var4/circulation;
      }

      /** The mean radius is obtained from $a=\sqrt{M_{2,0} + M_{0,2}}$, whereas
      the major and minos axes of the equivalent ellipse are given by
      $b = \sqrt{M_{2,0} + M_{0,2} + \sqrt{4M_{1,1}^2 + (M_{2,0} - M_{0,2})^2}}$
      and
      $c = \sqrt{M_{2,0} + M_{0,2} - \sqrt{4M_{1,1}^2 + (M_{2,0} - M_{0,2})^2}}$,
      respectively, while the ellipticity reads as $e = \sqrt{1 - c^2/b^2}$.
      */

      va = sqrt(M20 + M02);
      vb = sqrt((M20 + M02 + sqrt(4*sq(M11) + sq(M20 - M02))));
      vc = sqrt((M20 + M02 - sqrt(4*sq(M11) + sq(M20 - M02))));
      ve = sqrt(1 - sq(vc)/sq(vb));
    }
    //fprintf (stdout, "%.15g %d %d %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g \n", t, nj, pid(), circulation, centroids.x, centroids.y, centroids.z, M20, M02, M11, va, vb, vc, ve);

    omax = interpolate (omega, centroids.x, centroids.y, centroids.z);
    if (pid() == 0){
      @if _MPI
        MPI_Reduce (MPI_IN_PLACE, &omax, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
      @endif
    }
    @if _MPI
    else
      MPI_Reduce (&omax, NULL, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
    @endif

    FILE * fp;
    if ((pid()==0) & (circulation != 0)) {
      fp = fopen(name, "a");
      fprintf (fp, "%.4f %d %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g \n", t, nj, circulation, centroids.x, centroids.y, centroids.z, M20, M02, M11, va, vb, vc, ve, omax);
      fclose(fp);
    }
  }
}























trace
void vorticity_moments2_plane(vector omega, scalar m, int nm, char * name, coord n, double _alpha){

  double circulation, M20, M02, M11, va, vb, vc, ve, omax;
  coord centroids={0.,0.,0.}, tvec={0.,0.,0.};

  for (int nj = 0; nj < nm; nj++){
    // Compute the circulation $\Gamma = \iint \omega dA$ and position of the
    // vortex centroids $\vec{\mu} = \iint \vec{x}\omega dA$
    double var1 = 0., var2 = 0., var3 = 0., var4 = 0.;
    foreach(reduction(+:var1), reduction(+:var2), reduction(+:var3), reduction(+:var4)){
      double alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;
      if (fabs(alpha) > 0.87)
        continue;
      if (is_active(cell) && is_local(cell)) {
        coord sval;
        if (n.x == 1){
          sval.x = val(omega.x) + val(omega.x,1,0,0);
          sval.y = val(omega.y) + val(omega.y,1,0,0);
          sval.z = val(omega.z) + val(omega.z,1,0,0);
        }
        else if (n.y == 1){
          sval.x = val(omega.x) + val(omega.x,0,1,0);
          sval.y = val(omega.y) + val(omega.y,0,1,0);
          sval.z = val(omega.z) + val(omega.z,0,1,0);
        }
        else{
          sval.x = val(omega.x) + val(omega.x,0,0,1);
          sval.y = val(omega.y) + val(omega.y,0,0,1);
          sval.z = val(omega.z) + val(omega.z,0,0,1);
        }

        if (m[] == nj+1){
          var1 += dA * 0.5 * sval.x;
          var2 += dA * 0.5 * sval.y;
          var3 += dA * 0.5 * sval.z;
          var4 += dA * 0.25 * sqrt(sq(sval.x) + sq(sval.y) + sq(sval.z));
        }
      }
    }

    circulation = var4;

    double tol = 1e-14;
    if (circulation != 0) {
      tvec=(coord){var1, var2, var3};
      normalize(&tvec);

      double var1 = 0., var2 = 0., var3 = 0., var4 = 0.;
      foreach(reduction(+:var1), reduction(+:var2), reduction(+:var3), reduction(+:var4)){
        double alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;
        if (fabs(alpha) > 0.87)
          continue;
        if (is_active(cell) && is_local(cell)) {
          coord p = {x,y,z};

          double sval;
          if (n.x == 1)
            sval = 0.5*(val(omega.x) + val(omega.x,1,0,0)) * tvec.x + 0.5*(val(omega.y) + val(omega.y,1,0,0)) * tvec.y + 0.5*(val(omega.z) + val(omega.z,1,0,0)) * tvec.z;
          else if (n.y == 1)
            sval = 0.5*(val(omega.x) + val(omega.x,0,1,0)) * tvec.x + 0.5*(val(omega.y) + val(omega.y,0,1,0)) * tvec.y + 0.5*(val(omega.z) + val(omega.z,0,1,0)) * tvec.z;
          else
            sval = 0.5*(val(omega.x) + val(omega.x,0,0,1)) * tvec.x + 0.5*(val(omega.y) + val(omega.y,0,0,1)) * tvec.y + 0.5*(val(omega.z) + val(omega.z,0,0,1)) * tvec.z;

          if (m[] == nj+1){
            var1 += dA * sval;
            var2 += dA * sval * p.x;
            var3 += dA * sval * p.y;
            var4 += dA * sval * p.z;
          }
        }
      }

      centroids = (coord) {var2/var1, var3/var1, var4/var1};

      // For the ellipticity we require the centered vorticity moments
      // $ M_{p,q} = \iint (x - \mu_x)^p (y - \mu_y)^q \omega dA $

      var1 = 0., var2 = 0., var3 = 0., var4 = 0.;
      double var5 = 0., var6 = 0.;

      foreach(reduction(+:var1), reduction(+:var2), reduction(+:var3), reduction(+:var4), reduction(+:var5), reduction(+:var6)){
        double alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;
        if (fabs(alpha) > 0.87)
          continue;
        if (is_active(cell) && is_local(cell)) {
          coord p = {x,y,z};

          double sval;
          if (n.x == 1)
            sval = 0.5*(val(omega.x) + val(omega.x,1,0,0)) * tvec.x + 0.5*(val(omega.y) + val(omega.y,1,0,0)) * tvec.y + 0.5*(val(omega.z) + val(omega.z,1,0,0)) * tvec.z;
          else if (n.y == 1)
            sval = 0.5*(val(omega.x) + val(omega.x,0,1,0)) * tvec.x + 0.5*(val(omega.y) + val(omega.y,0,1,0)) * tvec.y + 0.5*(val(omega.z) + val(omega.z,0,1,0)) * tvec.z;
          else
            sval = 0.5*(val(omega.x) + val(omega.x,0,0,1)) * tvec.x + 0.5*(val(omega.y) + val(omega.y,0,0,1)) * tvec.y + 0.5*(val(omega.z) + val(omega.z,0,0,1)) * tvec.z;

          if (m[] == nj+1){
            var1 += dA * sval * sq(p.x - centroids.x);
            var2 += dA * sval * sq(p.y - centroids.y);
            var3 += dA * sval * sq(p.z - centroids.z);
            var4 += dA * sval * (p.x - centroids.x) * (p.y - centroids.y);
            var5 += dA * sval * (p.y - centroids.y) * (p.z - centroids.z);
            var6 += dA * sval * (p.z - centroids.z) * (p.x - centroids.x);
          }
        }
      }

      if (n.x == 1){
        M20 = var2/circulation;
        M02 = var3/circulation;
        M11 = var5/circulation;
      }
      else if (n.y == 1){
        M20 = var1/circulation;
        M02 = var3/circulation;
        M11 = var6/circulation;
      }
      else{
        M20 = var1/circulation;
        M02 = var2/circulation;
        M11 = var4/circulation;
      }

      /** The mean radius is obtained from $a=\sqrt{M_{2,0} + M_{0,2}}$, whereas
      the major and minos axes of the equivalent ellipse are given by
      $b = \sqrt{M_{2,0} + M_{0,2} + \sqrt{4M_{1,1}^2 + (M_{2,0} - M_{0,2})^2}}$
      and
      $c = \sqrt{M_{2,0} + M_{0,2} - \sqrt{4M_{1,1}^2 + (M_{2,0} - M_{0,2})^2}}$,
      respectively, while the ellipticity reads as $e = \sqrt{1 - c^2/b^2}$.
      */

      va = sqrt(M20 + M02);
      vb = sqrt((M20 + M02 + sqrt(4*sq(M11) + sq(M20 - M02))));
      vc = sqrt((M20 + M02 - sqrt(4*sq(M11) + sq(M20 - M02))));
      ve = sqrt(1 - sq(vc)/sq(vb));
      //fprintf (stdout, "%.15g %d %d %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g \n", t, nj, pid(), circulation, centroids.x, centroids.y, centroids.z, M20, M02, M11, va, vb, vc, ve);
    }

    var1 = interpolate (omega.x, centroids.x, centroids.y, centroids.z);
    var2 = interpolate (omega.y, centroids.x, centroids.y, centroids.z);
    var3 = interpolate (omega.z, centroids.x, centroids.y, centroids.z);
    omax = (var1 * tvec.x) + (var2 * tvec.y) + (var3 * tvec.z);
    if (pid() == 0){
      @if _MPI
        MPI_Reduce (MPI_IN_PLACE, &omax, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
      @endif
    }
    @if _MPI
    else
      MPI_Reduce (&omax, NULL, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
    @endif

    FILE * fp;
    if ((pid()==0) & (circulation != 0)) {
      fp = fopen(name, "a");
      fprintf (fp, "%.4f %d %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g \n", t, nj, circulation, centroids.x, centroids.y, centroids.z, M20, M02, M11, va, vb, vc, ve, omax);
      fclose(fp);
    }
  }
}