#include "filaments.h"


double test_vorticity_filament0(coord pcar, int n_seg, double a, double* t0, coord* c, coord* tvec, coord* nvec, coord* bvec, int period){
  /* Compute the local coordinates required for the vorticity field.
  Each position P(x,y,z) is projected into the local Frenet-Serret frame to
  obtain a set of local coordinates, such that:
  i)  (P - X . T = 0
  ii) (P - X) . N = x_n
  ii) (P - X) . B = x_b
  This requires finding the value of X(t0) along each space curve that verifies
  i) through a minization process.
  Then, we use the local coordinates (x_n, x_b) to define a radial coordinate
  rho required to compute the vorticity of a Lamb-Oseen vortex as
  omega = Gamma/(pi a²) exp(-rho²/a²) . T
  where Gamma is the circulation and a the core size.
  */

  /* First, we approximate the minimal distance between the point P and each
  segment of the curve. */
  double dmin = 1e30, tmin=0;
  double rho_loc;
  for (int i = 0; i < n_seg; i++){
    if (vecdist2(pcar, c[i]) < dmin){
      dmin = vecdist2(pcar, c[i]);
      tmin = t0[i];
    }
  }

  if (period != 0)
    tmin = fmod(tmin + period*2*pi, period*2*pi);

  // If P is close to the vortex, we refine the initial guess
  coord ploc, ccar, frenet[3];
  double tq = frenet_projection_min(n_seg, a, t0, c, tvec, nvec, bvec, pcar, tmin);

  ccar = gsl_interp1d( n_seg, t0, c, tq);

  // Then, compute the local coordinates for the vortex
  frenet[0] = gsl_interp1d( n_seg, t0, tvec, tq);
  frenet[1] = gsl_interp1d( n_seg, t0, nvec, tq);
  frenet[2] = gsl_interp1d( n_seg, t0, bvec, tq);

  ploc.x = vecdot(vecdiff(pcar, ccar), frenet[0]);
  ploc.y = vecdot(vecdiff(pcar, ccar), frenet[1]);
  ploc.z = vecdot(vecdiff(pcar, ccar), frenet[2]);
  rho_loc = sqrt(vecdot(ploc, ploc));

  return rho_loc;
}


coord test_vorticity_filament1(coord pcar, int n_seg, double a, double* t0, coord* c, coord* tvec, coord* nvec, coord* bvec, int period){
  double dmin = 1e30, tmin=0;
  for (int i = 0; i < n_seg; i++){
    if (vecdist2(pcar, c[i]) < dmin){
      dmin = vecdist2(pcar, c[i]);
      tmin = t0[i];
    }
  }

  if (period != 0)
    tmin = fmod(tmin + period*2*pi, period*2*pi);

  // If P is close to the vortex, we refine the initial guess
  coord ccar, frenet[1];
  double tq = frenet_projection_min(n_seg, a, t0, c, tvec, nvec, bvec, pcar, tmin);

  ccar = gsl_interp1d( n_seg, t0, c, tq);

  // Then, compute the local coordinates for the vortex
  frenet[0] = gsl_interp1d( n_seg, t0, tvec, tq);

  return (coord) {frenet[0].x, frenet[0].y, frenet[0].z};
}































coord get_vorticity_filament(coord pcar, int n_seg, double a, double* t0, coord* c, coord* tvec, coord* nvec, coord* bvec, int period){
  /* Compute the local coordinates required for the vorticity field.
  Each position P(x,y,z) is projected into the local Frenet-Serret frame to
  obtain a set of local coordinates, such that:
  i)  (P - X . T = 0
  ii) (P - X) . N = x_n
  ii) (P - X) . B = x_b
  This requires finding the value of X(t0) along each space curve that verifies
  i) through a minization process.
  Then, we use the local coordinates (x_n, x_b) to define a radial coordinate
  rho required to compute the vorticity of a Lamb-Oseen vortex as
  omega = Gamma/(pi a²) exp(-rho²/a²) . T
  where Gamma is the circulation and a the core size.
  */

  /* First, we approximate the minimal distance between the point P and each
  segment of the curve. */
  coord omega;
  double dmin = 1e30, tmin=0;
  double rho_loc, omega_mag;
  for (int i = 0; i < n_seg; i++){
    if (vecdist2(pcar, c[i]) < dmin){
      dmin = vecdist2(pcar, c[i]);
      tmin = t0[i];
    }
  }

  if (period != 0)
    tmin = fmod(tmin + period*2*pi, period*2*pi);

  if (dmin < L0){
    // If P is close to the vortex, we refine the initial guess
    coord ploc, ccar, frenet[3];
    double tq = frenet_projection_min(n_seg, a, t0, c, tvec, nvec, bvec, pcar, tmin);

    ccar = gsl_interp1d( n_seg, t0, c, tq);

    // Then, compute the local coordinates for the vortex
    frenet[0] = gsl_interp1d( n_seg, t0, tvec, tq);
    frenet[1] = gsl_interp1d( n_seg, t0, nvec, tq);
    frenet[2] = gsl_interp1d( n_seg, t0, bvec, tq);

    ploc.x = vecdot(vecdiff(pcar, ccar), frenet[0]);
    ploc.y = vecdot(vecdiff(pcar, ccar), frenet[1]);
    ploc.z = vecdot(vecdiff(pcar, ccar), frenet[2]);
    rho_loc = sqrt(vecdot(ploc, ploc));

    // Last, we compute the vorticity for the vortex 1
    omega_mag = exp(-sq(rho_loc)/sq(a))/(pi*sq(a));
    omega = (coord) {omega_mag * frenet[0].x, omega_mag * frenet[0].y, omega_mag * frenet[0].z};
  }
  else {
    // Otherwise, if the point is too far, we set the vorticity to zero.
    foreach_dimension()
      omega.x = 0;
  }
  return omega;
}





coord get_vorticity_filament2(double Gamma, double Uc, coord pcar, int n_seg, double a, double* t0, coord* c, coord* tvec, coord* nvec, coord* bvec, int period){
  /* Compute the local coordinates required for the vorticity field.
  Each position P(x,y,z) is projected into the local Frenet-Serret frame to
  obtain a set of local coordinates, such that:
  i)  (P - X . T = 0
  ii) (P - X) . N = x_n
  ii) (P - X) . B = x_b
  This requires finding the value of X(t0) along each space curve that verifies
  i) through a minization process.
  Then, we use the local coordinates (x_n, x_b) to define a radial coordinate
  rho required to compute the vorticity of a Lamb-Oseen vortex as
  omega = Gamma/(pi a²) exp(-rho²/a²) . T + 2 Uc/a (rho/a) exp(-rho²/a²)
  where Gamma is the circulation and a the core size and Uc is the
  axial centerline velocity.
  */

  /* First, we approximate the minimal distance between the point P and each
  segment of the curve. */
  coord omega;
  double dmin = 1e30, tmin=0;
  double rho_loc, omega_mag;
  for (int i = 0; i < n_seg; i++){
    if (vecdist2(pcar, c[i]) < dmin){
      dmin = vecdist2(pcar, c[i]);
      tmin = t0[i];
    }
  }

  if (period != 0)
    tmin = fmod(tmin + period*2*pi, period*2*pi);

  if (dmin < L0){
    // If P is close to the vortex, we refine the initial guess
    coord ploc, ccar, frenet[3];
    double tq = frenet_projection_min(n_seg, a, t0, c, tvec, nvec, bvec, pcar, tmin);

    ccar = gsl_interp1d( n_seg, t0, c, tq);

    // Then, compute the local coordinates for the vortex
    frenet[0] = gsl_interp1d( n_seg, t0, tvec, tq);
    frenet[1] = gsl_interp1d( n_seg, t0, nvec, tq);
    frenet[2] = gsl_interp1d( n_seg, t0, bvec, tq);

    ploc.x = vecdot(vecdiff(pcar, ccar), frenet[0]);
    ploc.y = vecdot(vecdiff(pcar, ccar), frenet[1]);
    ploc.z = vecdot(vecdiff(pcar, ccar), frenet[2]);
    rho_loc = sqrt(vecdot(ploc, ploc));

    // Last, we compute the vorticity for the vortex 1
    omega_mag = exp(-sq(rho_loc)/sq(a))/(sq(a));
    omega.x = omega_mag * (Gamma/pi * frenet[0].x + (2.0 * Uc) * (-ploc.z * frenet[1].x + ploc.y * frenet[2].x));
    omega.y = omega_mag * (Gamma/pi * frenet[0].y + (2.0 * Uc) * (-ploc.z * frenet[1].y + ploc.y * frenet[2].y));
    omega.z = omega_mag * (Gamma/pi * frenet[0].z + (2.0 * Uc) * (-ploc.z * frenet[1].z + ploc.y * frenet[2].z));
  }
  else {
    // Otherwise, if the point is too far, we set the vorticity to zero.
    foreach_dimension()
      omega.x = 0;
  }
  return omega;
}