---

markers.c

Go to the documentation of this file.

#include <malloc.h>
#include "utilf.h"
#include "markers.h"

/* Periodic conditions for the index i on a nx grid */
#define PERIODIC(i, nx, ip) {ip=(i-2)%(nx-2)+2; if(ip<2) ip+=nx-2;}
#define PERIODIC_DL(dl, nx) {if(fabs(dl)>0.5*(nx-2))\
                              dl=dl>0.0?dl-nx+2.:dl+nx-2.;}

int interface_arc(interface in, real t)
{
  int imin = 0, imax = in.n - 1, i;
  
  if (t < 0.0 || t > in.t[in.n - 1]) return -1;
  
  while(imax - imin > 1) {
    i = (imin + imax)/2;
    if (t == in.t[i]) return i;
    if (t > in.t[i]) imin = i;
    else imax = i;
  }
  return imin;
}


void interface_write_markers(interface in, FILE *fptr)
{
  int i;

  fprintf(fptr, "\n");
  for (i = 0; i <= in.n-1; i++)
    fprintf(fptr, "%g %g\n", in.x[i], in.y[i]); 
  fprintf(fptr, "&\n"); 
}

void interface_write_markers_speed(real time, interface in, real2D u, 
                real2D v, real tau, real du, 
                int nx, int ny, FILE *fptr)
{
/* Physical velocity*/
  int i;
  real h=1./(nx-2);

  fprintf(fptr, "\n");
  for (i = 0; i <= in.n-1; i++)
    fprintf(fptr, "%g %g %g %g %g\n", time, in.x[i], in.y[i], 
                    bilinu(u, in.x[i], in.y[i])*du*h/tau, 
                    bilinu(u, in.x[i], in.y[i])*du*h/tau); 
  fprintf(fptr, "&\n"); 
}

void interface_write_splines(interface in, FILE *fptr, int nb)
{
  int i;
  real t;

  fprintf(fptr, "\n");
  for (i = 0; i < nb; i++) {
    t = (real) i * in.t[in.n-1] / (real) (nb - 1);
    fprintf(fptr, "%g %g\n", 
            splint_find(in.t, in.spx, in.n, t),
            splint_find(in.t, in.spy, in.n, t));
  }
  fprintf(fptr, "&\n");
}


void curvature_write(interface in, FILE *fptr)
{
  int i;
  real x1, y1, dl, kappa;

  fprintf(fptr, "\n");
  for (i = 0; i < in.n - 1; i++) {
    x1 = splint1(in.spx[i], in.t[i]);
    y1 = splint1(in.spy[i], in.t[i]);    
    dl = sqrt(sq(x1) + sq(y1));
    kappa = (x1*splint2(in.spy[i], in.t[i]) - y1*splint2(in.spx[i], in.t[i]))
            / (dl*dl*dl);
    fprintf(fptr, "%g %g\n", in.t[i]/in.t[in.n-1], kappa);
  }
}


real interfaces_area(interface *in, int n)
{
  int i, j;
  real ar = 0.0, a1, a2, a3, a4, a5, a6, t2, t1;
  interface in1;

  for (j = 0; j < n; j++) {
    in1 = in[j];
    for (i = 0; i < in1.n - 1; i++) 
      {
        a1 = in1.spx[i].a*in1.spy[i].b; 
        a2 = 0.5*(in1.spx[i].b*in1.spy[i].b + 2.*in1.spx[i].a*in1.spy[i].c);
        a3 = (in1.spx[i].c*in1.spy[i].b + 2.*in1.spx[i].b*in1.spy[i].c 
              + 3.*in1.spx[i].a*in1.spy[i].d)/3.;
        a4 = 0.25*(in1.spx[i].d*in1.spy[i].b + 2.*in1.spx[i].c*in1.spy[i].c 
                   + 3.*in1.spx[i].b*in1.spy[i].d);
        a5 = (3.*in1.spx[i].c*in1.spy[i].d + 2.*in1.spx[i].d*in1.spy[i].c)/5.;
        a6 = 0.5*in1.spx[i].d*in1.spy[i].d;
        t1 = in1.t[i]; t2 = in1.t[i+1];
        ar += t2*(a1 + t2*(a2 + t2*(a3 + t2*(a4 + t2*(a5 + t2*a6))))) -
          t1*(a1 + t1*(a2 + t1*(a3 + t1*(a4 + t1*(a5 + t1*a6)))));
      }
  }
  return ar;
}


real interfaces_axivolume(interface *in, int n)
{
  int i, j;
  real ar = 0.0, a1, a2, a3, a4, a5, a6, a7, a8, a9, t2, t1;
  real ax, bx, cx, dx, ay, by, cy, dy;
  real x1, y1, x2, y2;
  interface in1;

  for (j = 0; j < n; j++) {
    in1 = in[j];
    for (i = 0; i < in1.n - 1; i++) {
#ifdef POLYVOLUME
      /* integration of xy dy along the polynomial segment */
      ax = in1.spx[i].a; bx = in1.spx[i].b;
      cx = in1.spx[i].c; dx = in1.spx[i].d;
      ay = in1.spy[i].a - 2.0; by = in1.spy[i].b;
      cy = in1.spy[i].c; dy = in1.spy[i].d;
      
      a1 = ax*ay*by;
      a2 = ax*ay*cy+(ax*by+bx*ay)*by/2.;
      a3 = ax*ay*dy+2.0/3.0*(ax*by+bx*ay)*cy+(ax*cy+bx*by+cx*ay)*by/3.;
      a4 = 3.0/4.0*(ax*by+bx*ay)*dy+(ax*cy+bx*by+cx*ay)*cy/2.
        +(ax*dy+bx*cy+cx*by+dx*ay)*by/4.;
      a5 = 3.0/5.0*(ax*cy+bx*by+cx*ay)*dy
        +2.0/5.0*(ax*dy+bx*cy+cx*by+dx*ay)*cy+(bx*dy+cx*cy+dx*by)*by/5.;
      a6 = (ax*dy+bx*cy+cx*by+dx*ay)*dy/2+(bx*dy+cx*cy+dx*by)*cy/3
        +(cx*dy+dx*cy)*by/6.;
      a7 = 3.0/7.0*(bx*dy+cx*cy+dx*by)*dy+2.0/7.0*(cx*dy+dx*cy)*cy
        +dx*dy*by/7.;
      a8 = 3.0/8.0*(cx*dy+dx*cy)*dy+dx*dy*cy/4.;
      a9 = dx*dy*dy/3.;
      
      t1 = in1.t[i]; t2 = in1.t[i+1];
      ar += t2*(a1 + t2*(a2 + t2*(a3 + t2*(a4 + t2*(a5 + t2*(a6 + t2*(a7 + t2*(a8 + t2*a9)))))))) -
        t1*(a1 + t1*(a2 + t1*(a3 + t1*(a4 + t1*(a5 + t1*(a6 + t1*(a7 + t1*(a8 + t1*a9))))))));
#else
      /* integration of xy dy along the segment */
      x1 = in1.x[i]; y1 = in1.y[i] - 2.;
      x2 = in1.x[i+1]; y2 = in1.y[i+1] - 2.;
      ar += ((x2 - x1)*(y2 - y1)/3. + 
             (y1*(x2 - x1) + x1*(y2 - y1))/2. + x1*y1)*(y2 - y1);
#endif
    }
  }
  return 2.*PI*ar;
}


real interfaces_xmoment(interface *in, int n)
{
  /* compute Sum{x^2.dx.dy} */
  /* coefficients of Newton-Cotes formula */
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i, j, k;
  real x, yp, h, t1, t2, t, sum, total = 0.0;
  interface in1;
  polynom3 px, py;

  for (k = 0; k < n; k++) {
    in1 = in[k];
    for (j = 0; j < in1.n - 1; j++) {
      t1 = in1.t[j]; t2 = in1.t[j+1];
      px = in1.spx[j]; py = in1.spy[j];
      sum = 0.0; h = (t2 - t1)/6.;
      for (i = 0; i < 7; i++) {
        t = t1 + h*i;
        x = splint(px, t) - 2.;
        yp = splint1(py, t);
        sum += coef[i]*x*x*x*yp;
      }
      total += 3.*h*sum/8.;
    }
  }
  return total/3.;
}


real interfaces_ymoment(interface *in, int n)
{
  /* compute Sum{y^2.dx.dy} */
  /* coefficients of Newton-Cotes formula */
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i, j, k;
  real y, xp, h, t1, t2, t, sum, total = 0.0;
  interface in1;
  polynom3 px, py;

  for (k = 0; k < n; k++) {
    in1 = in[k];
    for (j = 0; j < in1.n - 1; j++) {
      t1 = in1.t[j]; t2 = in1.t[j+1];
      px = in1.spx[j]; py = in1.spy[j];
      sum = 0.0; h = (t2 - t1)/6.;
      for (i = 0; i < 7; i++) {
        t = t1 + h*i;
        y = splint(py, t) - 2.;
        xp = splint1(px, t);
        sum += coef[i]*y*y*y*xp;
      }
      total += 3.*h*sum/8.;
    }
  }
  return total/3.;
}


real interfaces_x_amp(interface *in, int n)
{
  int i, j;
  real max = -1e6, min = 1e6;

  for (j = 0; j < n; j++)
    for (i = 0; i < in[j].n; i++) {
      max = in[j].x[i] > max ? in[j].x[i] : max;
      min = in[j].x[i] < min ? in[j].x[i] : min;
    }
  return fabs(max - min);
}


real interfaces_y_amp(interface *in, int n)
{
  int i, j;
  real max = -1e6, min = 1e6;

  for (j = 0; j < n; j++)
    for (i = 0; i < in[j].n; i++) {
      max = in[j].y[i] > max ? in[j].y[i] : max;
      min = in[j].y[i] < min ? in[j].y[i] : min;
    }
  return fabs(max - min);
}


void interface_splines(interface in, int nx, int ny)
{
  real *d2y, l = 0.0, d0, d1, det, dy0, dy1, yp;
  real frac;
  int i;

  d2y = (real *)malloc(in.n * sizeof(real));

  in.t[0] = 0.0;
  for (i = 1; i < in.n; i++) {
    /* Move slightly the points when they are near the grid lines */
    frac = fabs(modf(in.x[i], &d0));
    if (frac < INTERPREC || 1. - frac < INTERPREC || 
        fabs(0.5 - frac) < INTERPREC) {
      printf("WARNING: interfaces_splines(): moving x[%d]: %.10f to %.10f\n",
             i, in.x[i], in.x[i] + 10.*INTERPREC);
      in.x[i] += 10.*INTERPREC;
    }
    frac = fabs(modf(in.y[i], &d0));
    if (frac < INTERPREC || 1. - frac < INTERPREC || 
        fabs(0.5 - frac) < INTERPREC) {
      printf("WARNING: interfaces_splines(): moving y[%d]: %.10f to %.10f\n", 
             i, in.y[i], in.y[i] + 10.*INTERPREC);
      in.y[i] += 10.*INTERPREC;
    }

    d0 = sqrt(sq(in.x[i] - in.x[i-1]) + sq(in.y[i] - in.y[i-1]));
/*
    if (d0 >= 1.0)
      printf("WARNING: interfaces_splines(): d0=%f >= 1.0\n", d0);
*/
    l += d0;
    in.t[i] = l;
  }

  switch(in.bc) {
  case TRUE_PERIODIC:
    Pspline(in.t, in.x, d2y, in.t[in.n - 1], in.n - 1, in.spx);
    Pspline(in.t, in.y, d2y, in.t[in.n - 1], in.n - 1, in.spy);
    break;
  case APPROX_PERIODIC:
    /* We need an estimate of the derivative for the periodic ends 
       We suppose a 2nd order interpolation for the points 1, 0, n-2
       will do */
    d0 = in.t[in.n-1] - in.t[in.n-2];  
    d1 = d0 + in.t[1] - in.t[0];
    det = d1*d0*(d0 - d1);
    
    dy0 = in.x[0] - in.x[in.n-2];
    dy1 = in.x[1] - in.x[in.n-2];
    PERIODIC_DL(dy0, nx);
    PERIODIC_DL(dy1, nx);
    yp = (2.*(dy0*d1 - dy1*d0)*d0 + dy1*d0*d0 - dy0*d1*d1) / det;
    spline(in.t, in.x, in.n, yp, yp, d2y, in.spx, SP_DERIVATIVE);
    
    dy0 = in.y[0] - in.y[in.n-2];
    dy1 = in.y[1] - in.y[in.n-2];
    yp = (2.*(dy0*d1 - dy1*d0)*d0 + dy1*d0*d0 - dy0*d1*d1) / det;
    spline(in.t, in.y, in.n, yp, yp, d2y, in.spy, SP_DERIVATIVE);
    break;
  case AXITRUE:
    spline(in.t, in.x, in.n, 0.0, 0.0, d2y, in.spx, SP_DERIVATIVE);
    spline(in.t, in.y, in.n, 1.0, -1.0, d2y, in.spy, SP_DERIVATIVE);
    break;
  case AXINATURAL:
    spline(in.t, in.x, in.n, 0.0, 0.0, d2y, in.spx, SP_NATURAL);
    spline(in.t, in.y, in.n, 0.0, 0.0, d2y, in.spy, SP_NATURAL);
    break;
  case AXIBOTH:
    spline(in.t, in.x, in.n, 0.0, 0.0, d2y, in.spx, SP_DERIVATIVE);
    spline(in.t, in.y, in.n, 1.0, 1.0, d2y, in.spy, SP_DERIVATIVE);
    break;
  default:
    printf("interface_splines(): bad value for the bc parameter: bc = %d\n", in.bc);
    exit(1);
  }
    
  free(d2y);
}


void interface_redis(interface *in, real space)
{
  int i, j, nb, k = 0;
  real l, dl, xe, ye;

  nb = (int)(in->t[in->n-1] / space) + 2;
  dl = in->t[in->n-1] / (real)(nb - 1);
  if (dl >= 1.0)
    printf("WARNING interface_redis(): dl=%g\n", dl);

  xe = in->x[in->n-1];  /* saving the coordinates of the last point */
  ye = in->y[in->n-1];  

  in->x = (real *)realloc(in->x, nb * sizeof(real));
  in->y = (real *)realloc(in->y, nb * sizeof(real));

  l = dl;               /* Starting at i=1 the first point is unchanged */
  for (i = 1; i < nb - 1; i++, l += dl) {
    /* find the next interval [j;j+1] such as t[j] <= l <= t[j+1] */
    while(in->t[k+1] < l)
      k++;

#ifdef DEBUG_REDIS
    if (k >= in->n-1) {
      printf("interface_redis: Error i=%d/%d k too large: %d/%d\n",
              i, nb - 1, k, in->n - 1);
      exit(1);
    }
    if (in->t[k] > l || in->t[k+1] < l) {
      printf("interface_redis: Error i=%d/%d k=%d, %g <= %g <= %g\n",
              i, nb - 1, k, in->t[k], l, in->t[k+1]);
      exit(1);
    }
#endif    
    
    in->x[i] = splint(in->spx[k], l);
    in->y[i] = splint(in->spy[k], l);
  }

  in->x[nb - 1] = xe;  /* The last point is unchanged */
  in->y[nb - 1] = ye;  /* Copying the saved coordinates */
  in->t = (real *)realloc(in->t, nb * sizeof(real));
  in->spx = (polynom3 *)realloc(in->spx, (nb - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)realloc(in->spy, (nb - 1) * sizeof(polynom3));
#ifdef FREE
  in->p = (real *)realloc(in->p, nb * sizeof(real));
#endif
  in->n = nb;
}


void interface_redis_scale(interface *in, interface src, real space, 
                           real scale)
{
  int i, j, nb, k = 0;
  real l, dl, xe, ye;

  nb = (int)(src.t[src.n-1] / space) + 2;
  dl = src.t[src.n-1] / (real)(nb - 1);

  in->x = (real *)realloc(in->x, nb * sizeof(real));
  in->y = (real *)realloc(in->y, nb * sizeof(real));
  in->t = (real *)realloc(in->t, nb * sizeof(real));
  in->spx = (polynom3 *)realloc(in->spx, (nb - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)realloc(in->spy, (nb - 1) * sizeof(polynom3));
#ifdef FREE
  in->p = (real *)realloc(in->p, nb * sizeof(real));
#endif
  in->n = nb;

  l = dl;               /* Starting at i=1 the first point is unchanged */
  for (i = 1; i < nb - 1; i++, l += dl) {
    /* find the next interval [j;j+1] such as t[j] <= l <= t[j+1] */
    while(src.t[k+1] < l)
      k++;
    in->x[i] = scale*(splint(src.spx[k], l) - 2.) + 2.;
    in->y[i] = scale*(splint(src.spy[k], l) - 2.) + 2.;
  }
  in->x[0] = scale*(src.x[0] - 2.) + 2.;
  in->y[0] = scale*(src.y[0] - 2.) + 2.;
  in->x[nb - 1] = scale*(src.x[src.n - 1] - 2.) + 2.;
  in->y[nb - 1] = scale*(src.y[src.n - 1] - 2.) + 2.;
}


void interface_filter(interface *in, int nx, int ny)
{
  int i;
  real t;

  for (i = 0; i < in->n - 1; i++) {
    t = (in->t[i] + in->t[i+1])/2.;
    in->x[i] = splint(in->spx[i], t);
    in->y[i] = splint(in->spy[i], t);
  }
  in->x[in->n-1] = in->x[0];
  in->y[in->n-1] = in->y[0];
  
  interface_splines(*in, nx, ny);
}

/*cld
void interface_axisymmetry(interface *in, real space, int nx, int ny)
{
  int i, k, nb = 2 * (in->n - 1) + 1, nold = in->n, n, nnew;
  real t, ts, te, ti[4], xe, xs, l, dl;
  char old_bc = in->bc;

  if (in->bc == TRUE_PERIODIC || in->bc == APPROX_PERIODIC) {
    interface_filter(in, nx, ny);
    interface_redis(in, space);
    return;
  }

  in->x = (real *)realloc(in->x, nb * sizeof(real));
  in->y = (real *)realloc(in->y, nb * sizeof(real));
  in->t = (real *)realloc(in->t, nb * sizeof(real));
  in->spx = (polynom3 *)realloc(in->spx, (nb - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)realloc(in->spy, (nb - 1) * sizeof(polynom3));

  for (i = 0; i < in->n - 1; i++) {
    t = (in->t[i] + in->t[i+1])/2.;
    in->x[i] = splint(in->spx[i], t);
    in->y[i] = splint(in->spy[i], t);
  }
  for (i = in->n - 1; i < nb - 1; i++) {
    in->x[i] = in->x[nb - 2 - i];
    in->y[i] = 4. - in->y[nb - 2 - i];
  }
  in->x[nb-1] = in->x[0]; in->y[nb-1] = in->y[0];
  in->n = nb;
  in->bc = TRUE_PERIODIC;
  
  interface_splines(*in, nx, ny);
  n = polyroots(in->spy[nb-2], in->t[nb-2], in->t[nb-1], 
                ti, 2.0, INTERPREC);
  if (n != 1) {
    FILE *fptr = fopen("aximarkers", "wt");
    interface_write_markers(*in, fptr);
    fclose(fptr);
    printf("ERROR: interface_axisymmetry(): n = %d\n", n);    
    exit(1);
  }
  xs = splint(in->spx[nb-2], ti[0]);
  ts = ti[0] - in->t[nb-1];
  n = polyroots(in->spy[nold-2], in->t[nold-2], in->t[nold-1], 
                ti, 2.0, INTERPREC);
  if (n != 1) {
    FILE *fptr = fopen("aximarkers", "wt");
    interface_write_markers(*in, fptr);
    fclose(fptr);
    printf("ERROR: interface_axisymmetry(): n = %d\n", n);
    exit(1);
  }
  xe = splint(in->spx[nold-2], ti[0]);
  te = ti[0];
  nnew = (int)((te - ts)/space) + 2;
  if (nnew % 2) nnew++;
  dl = (te - ts)/(real)(nnew - 1);
  if (dl >= 1.0)
    printf("WARNING interface_redis(): dl=%g\n", dl);

  in->x = (real *)realloc(in->x, nnew * sizeof(real));
  in->y = (real *)realloc(in->y, nnew * sizeof(real));

  l = ts + dl; k = 0;
  if (l < 0.0) printf("interface_axisymmetry(): l = %g\n", l);
  for (i = 1; i < nnew - 1; i++, l += dl) {
    /* find the next interval [j;j+1] such as t[j] <= l <= t[j+1] 
    while(in->t[k+1] < l)
      k++;
    in->x[i] = splint(in->spx[k], l);
    in->y[i] = splint(in->spy[k], l);
  }
  in->x[0] = xs; in->y[0] = 1.9999;
  in->x[nnew-1] = xe; in->y[nnew-1] = 1.9999;
  in->t = (real *)realloc(in->t, nnew * sizeof(real));
  in->spx = (polynom3 *)realloc(in->spx, (nnew - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)realloc(in->spy, (nnew - 1) * sizeof(polynom3));
#ifdef FREE
  in->p = (real *)realloc(in->p, nnew * sizeof(real));
#endif
  in->n = nnew;
  in->bc = old_bc;
}
cld*/

void interface_copy(interface *dst, interface in)
{
  int l;
  dst->x = (real *)malloc(in.n * sizeof(real));
  dst->y = (real *)malloc(in.n * sizeof(real));
  dst->t = (real *)malloc(in.n * sizeof(real));
  dst->spx = (polynom3 *)malloc((in.n - 1) * sizeof(polynom3));
  dst->spy = (polynom3 *)malloc((in.n - 1) * sizeof(polynom3));
#ifdef FREE
  dst->p = (real *)calloc(in.n, sizeof(real));
#endif
  dst->n = in.n;
  for (l = 0; l < in.n - 1; l++) {
    dst->x[l] = in.x[l];
    dst->y[l] = in.y[l];
    dst->t[l] = in.t[l];
    dst->spx[l] = in.spx[l];
    dst->spy[l] = in.spy[l];
  }
  dst->x[in.n-1] = in.x[in.n-1];
  dst->y[in.n-1] = in.y[in.n-1];
  dst->t[in.n-1] = in.t[in.n-1];
}


void interface_init(interface *dst)
{
  dst->x = (real *)malloc(sizeof(real));
  dst->y = (real *)malloc(sizeof(real));
  dst->t = (real *)malloc(sizeof(real));
  dst->spx = (polynom3 *)malloc(sizeof(polynom3));
  dst->spy = (polynom3 *)malloc(sizeof(polynom3));
#ifdef FREE
  dst->p = (real *)calloc(1, sizeof(real));
#endif
  dst->n = 0;
}


void interface_free(interface in)
{
  free(in.x); free(in.y); free(in.t);
  free(in.spx); free(in.spy);
#ifdef FREE
  free(in.p);
#endif
}


void interface_advect(interface in, real2D u, real2D v, 
#ifdef FREE
                      real2D a, real2D c, real2D cc,
#endif
                      int nx, int ny)
{
#ifdef BICUB
  int i, j, l, ip;
  int iuold = -1, juold = -1, ivold = -1, jvold = -1;
  real cu[16], cv[16], dx, dy;

  for (l = 0; l < in.n; l++)
    {
      i = in.x[l]; j = (int) ((real) in.y[l] - 0.5);
      PERIODIC(i, nx, ip);
      if (ip != iuold || j != juold) {
        bcucof(u, ip, j, cu, nx, ny);
        iuold = ip; juold = j;
      }
      dx = bcuint(in.x[l] - (real) i, in.y[l] - 0.5 - (real) j, cu);

      i = (int) ((real) in.x[l] - 0.5); j = in.y[l];
      PERIODIC(i, nx, ip);
      if (ip != ivold || j != jvold) {
        bcucof(v, ip, j, cv, nx, ny);
        ivold = ip; jvold = j;
      }
      dy = bcuint(in.x[l] - 0.5 - (real) i, in.y[l] - (real) j, cv);

      in.x[l] += dx; in.y[l] += dy;
    }
#endif
#ifdef BILIN
  int i, j, ip, l;
  real dx, dy, x, y;

  for (l = 0; l < in.n; l++)
    {
      i = in.x[l]; j = (int) ((real) in.y[l] - 0.5);
      PERIODIC(i, nx, ip);
      x = in.x[l] - (real) i;
      y = in.y[l] - 0.5 - (real) j;
      dx = u[ip][j]*(1. - x - y + x*y) + u[ip+1][j]*x*(1. - y) +
        u[ip][j+1]*y*(1. - x) + u[ip+1][j+1]*x*y;

      i = (int) ((real) in.x[l] - 0.5); j = in.y[l];
      PERIODIC(i, nx, ip);
      x = in.x[l] - 0.5 - (real) i;
      y = in.y[l] - (real) j;
      dy = v[ip][j]*(1. - x - y + x*y) + v[ip+1][j]*x*(1. - y) +
        v[ip][j+1]*y*(1. - x) + v[ip+1][j+1]*x*y;
      
#ifdef AXI
      /* Don't move the first and last markers vertically */
      if (l > 0 && l < in.n - 1)
        in.y[l] += dy;
      in.x[l] += dx;
#else
      in.x[l] += dx; in.y[l] += dy;
#endif
    }
#endif
#ifdef FREE
  int l;
  real dx, dy;
#ifdef AXI
  int i, j, ip;
  real x, y;

  /* Don't move the first marker vertically */
  i = in.x[0]; j = (int) ((real) in.y[0] - 0.5);
  PERIODIC(i, nx, ip);
  x = in.x[0] - (real) i;
  y = in.y[0] - 0.5 - (real) j;
  in.x[0] += u[ip][j]*(1. - x - y + x*y) + u[ip+1][j]*x*(1. - y) +
    u[ip][j+1]*y*(1. - x) + u[ip+1][j+1]*x*y;

  for (l = 1; l < in.n - 1; l++)
    {
      i = in.x[l]; j = (int) ((real) in.y[l] - 0.5);
      PERIODIC(i, nx, ip);
      x = in.x[l] - (real) i;
      y = in.y[l] - 0.5 - (real) j;
      dx = u[ip][j]*(1. - x - y + x*y) + u[ip+1][j]*x*(1. - y) +
        u[ip][j+1]*y*(1. - x) + u[ip+1][j+1]*x*y;

      i = (int) ((real) in.x[l] - 0.5); j = in.y[l];
      PERIODIC(i, nx, ip);
      x = in.x[l] - 0.5 - (real) i;
      y = in.y[l] - (real) j;
      dy = v[ip][j]*(1. - x - y + x*y) + v[ip+1][j]*x*(1. - y) +
        v[ip][j+1]*y*(1. - x) + v[ip+1][j+1]*x*y;

      in.x[l] += dx; in.y[l] += dy;
    }
  /* Don't move the last marker vertically */
  i = in.x[in.n-1]; j = (int) ((real) in.y[in.n-1] - 0.5);
  PERIODIC(i, nx, ip);
  x = in.x[in.n-1] - (real) i;
  y = in.y[in.n-1] - 0.5 - (real) j;
  in.x[in.n-1] += u[ip][j]*(1. - x - y + x*y) + u[ip+1][j]*x*(1. - y) +
    u[ip][j+1]*y*(1. - x) + u[ip+1][j+1]*x*y;

#else
  real xn, yn, dl;

  xn = in.y[in.n-2] - in.y[1];
  yn = in.x[1] - in.x[in.n-2];
  dl = sqrt(sq(xn) + sq(yn));
  xn /= dl; yn /= dl;
  dx = linear(u, a, in.x[0], in.y[0], xn, yn, bilin_u);
  dy = linear(v, c, in.x[0], in.y[0], xn, yn, bilin_v);
  /*
    dx = interplate_u(u, a, cc, in.x[0], in.y[0], nx, ny, INTERPLATE_INTER);
    dy = interplate_v(v, c, cc, in.x[0], in.y[0], nx, ny, INTERPLATE_INTER);
  */
  in.x[0] += dx; in.y[0] += dy;

  for (l = 1; l < in.n - 1; l++) {
    
    /* xn = -splint1(in.spy[l], in.t[l]); */
      xn = in.y[l-1] - in.y[l+1];
      /*      yn = splint1(in.spx[l], in.t[l]); */
      yn = in.x[l+1] - in.x[l-1];
      dl = sqrt(sq(xn) + sq(yn));
      xn /= dl; yn /= dl;
      
      dx = linear(u, a, in.x[l], in.y[l], xn, yn, bilin_u);
      dy = linear(v, c, in.x[l], in.y[l], xn, yn, bilin_v);
    /*
    dx = interplate_u(u, a, cc, in.x[l], in.y[l], nx, ny, INTERPLATE_INTER);
    dy = interplate_v(v, c, cc, in.x[l], in.y[l], nx, ny, INTERPLATE_INTER);
    */
    in.x[l] += dx; in.y[l] += dy;
  }  
  in.x[in.n-1] = in.x[0];
  in.y[in.n-1] = in.y[0];
#endif
#endif
}


static void intersec_update(intersec *in1, intersec *in2, int o, int j)
{
  if (in2->n == -1) {
    in2->j2[0] = j; in2->in[0] = o;
  }
  else {
    in2->j2[in2->n - 1] = j; in2->in[in2->n - 1] = o;
  }

  if (in1->n == -1) {
    in1->j1[0] = j; in1->in[0] = o;
    in1->n = 1;
  }
  else {
    in1->in[in1->n] = o;
    in1->j1[in1->n++] = j;
  }
}


#ifdef AXI
real stensionaxix(polynom3 px, polynom3 py, real t1, real t2)
{
  /* coefficients of Newton-Cotes formula */
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i;
  real xp, yp, y, h = (t2 - t1)/6., t, sum = 0.0;

  for (i = 0; i < 7; i++) {
    t = t1 + h*i;
    xp = splint1(px, t);
    yp = splint1(py, t);
    y = splint(py, t) - 2.; /* Position of the curve relative to the axis of 
                               symetry: y = splint(py, t) + yo - 2 */
    if (fabs(y) <= 0.01) /* second order approximation near the axis */
      sum += coef[i]*splint2(px, t)/sqrt(xp*xp + yp*yp);
    else
      sum += coef[i]*xp*yp/(y*sqrt(xp*xp + yp*yp));
  }
  return 3.*h*sum/8.;
}


real stensionaxiy(polynom3 px, polynom3 py, real t1, real t2)
{
  /* coefficients of Newton-Cotes formula */
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i;
  real xp, yp, y, h = (t2 - t1)/6., t, sum = 0.0;

  for (i = 0; i < 7; i++) {
    t = t1 + h*i;
    xp = splint1(px, t);
    yp = splint1(py, t);
    y = splint(py, t) - 1.5; /* Position of the curve relative to the axis of 
                                 symetry: y = splint(py, t) + yo - 2 */
    if (fabs(y) <= 0.01) /* second order approximation near the axis */
      sum -= coef[i]*splint2(px, t)*xp/(yp*sqrt(xp*xp + yp*yp));
    else
      sum -= coef[i]*xp*xp/(y*sqrt(xp*xp + yp*yp));
  }
  return 3.*h*sum/8.;
}
#endif


void interface_to_c(interface *in, int nint, 
#ifdef RECON
                    intersec **inter, 
#endif
                    real2D work, real2D c, real xo, real yo,
                    int nx, int ny, int sflag,
                    real2D rhox, real2D p, real sigma, real tau)
{
  int i, ip, ip1, ip2, j, k, l, m, o;
  int n, nin;
  real a, a1, a2, a3, a4, a5, a6, ti[14], t, t1, t2;
  real x, y, y1, dx, dy, h;
  polynom3 px, py;
  typedef struct {
    real t;
    int type, coord;
  } listinter;
  listinter listi[14], xl;

  h = 1.0 / (nx - 2);
  sigma *= tau * tau / (h * h * h);
  
  /* No intersections, no volume fraction */
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {      
#ifdef RECON
      inter[i][j].n = 0;
#endif
      work[i][j] = 0.0;
    }
      
  /* For each interface */
  for (o = 0; o < nint; o++)
    /* For each segment */  
    for (l = 0; l < in[o].n - 1; l++)
      {
         /* the positions are shifted (-xo, -yo) */
        i = (real)(in[o].x[l] - xo); 
        j = (real)(in[o].y[l] - yo);
        t1 = in[o].t[l]; t2 = in[o].t[l+1];
        px = in[o].spx[l]; py = in[o].spy[l];
        px.a -= xo; py.a -= yo;
        
        /* first point */
        listi[0].t = t1; listi[0].type = MARK_TYPE;
        nin = 1;
        
        /* intersections with the horizontals */
        n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
        for (k = 0; k < n; k++) {
          listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
          listi[nin++].coord = j;
        }
        n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
        for (k = 0; k < n; k++) {
          listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
          listi[nin++].coord = j + 1;
        }
        /* intersections with the verticals */
        n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
        for (k = 0; k < n; k++) {
          listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
          listi[nin++].coord = i;
        }
        n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
        for (k = 0; k < n; k++) {
          listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
          listi[nin++].coord = i + 1;
        }
        
        /* last point */
        listi[nin].t = t2; listi[nin++].type = MARK_TYPE;
        
        /* sort the intersections in ascending order for t */
        for (k = 1; k < nin; k++) {
          xl = listi[k]; i = k - 1;
          while (i >= 0 && listi[i].t > xl.t) {
            listi[i + 1] = listi[i]; i--;
          }
          listi[i + 1] = xl;
        }
        
        /* coefficients for the computation of the circulation */
        a1 = px.a*py.b; 
        a2 = 0.5*(px.b*py.b + 2.*px.a*py.c);
        a3 = (px.c*py.b + 2.*px.b*py.c + 3.*px.a*py.d)/3.;
        a4 = 0.25*(px.d*py.b + 2.*px.c*py.c + 3.*px.b*py.d);
        a5 = (3.*px.c*py.d + 2.*px.d*py.c)/5.;
        a6 = 0.5*px.d*py.d;      
        
        for (k = 0; k < nin - 1; k++) 
          {
            t1 = listi[k].t; t2 = listi[k+1].t;
            t = 0.5*(t1 + t2);
            i = splint(px, t);
            j = splint(py, t);
            
            x = splint(px, t1);
            y = splint(py, t1);
            
            /* Intersection with the verticals */
            if (listi[k].type == VER_TYPE) {
              /* add the circulation of the vertical segment */     
              if (listi[k].coord == i) {
                t = y - j - 1.0; m = -1;
                PERIODIC(i, nx, ip);
                work[ip][j] += t*ip;
                PERIODIC(i - 1, nx, ip1);
                work[ip1][j] -= t*(ip1 + 1);
              } 
              else {
                t = y - j; m = 1;
                PERIODIC(i, nx, ip);
                work[ip][j] += t*(ip + 1);
                PERIODIC(i + 1, nx, ip1);
                work[ip1][j] -= t*ip1;
              } 
              
              /* update the intersections structure */
#ifdef RECON
              if (inter[ip][j].n >= NSEC)
                printf("WARNING intersec_update: inter[%d][%d].n = %d\n",
                       ip, j, inter[ip][j].n);
              else
                intersec_update(&inter[ip][j], &inter[ip1][j], o, l);
#endif
              /* Add the surface tension term and the pressure
                 gradient correction term */
              switch(sflag) {
              case TENSION_U:
                dx = splint1(px, t1); dy = splint1(py, t1);
                a = sigma*dx / sqrt(dx*dx + dy*dy);
                PERIODIC(i + 1, nx, ip);
                rhox[ip][j] -= a;
                PERIODIC(i + m + 1, nx, ip1);
                rhox[ip1][j] += a;
                /* pressure gradient correction */
                a = y - (real) j;
                if (m == 1) {
                  if (a > 0.5) 
                    a = (1.0 - a)*(p[ip][j+1] - p[ip][j]);
                  else
                    a = a*(p[ip][j-1] - p[ip][j]);
                  rhox[ip][j] -= a;
                  rhox[ip1][j] += a;
                } 
                else {
                  PERIODIC(i, nx, ip2);
                  if (a > 0.5)
                    a = (1.0 - a)*(p[ip2][j+1] - p[ip2][j]);
                  else
                    a = a*(p[ip2][j-1] - p[ip2][j]);
                  rhox[ip][j] += a;
                  rhox[ip1][j] -= a;
                }
                break;
              case TENSION_V:
                dx = splint1(px, t1); dy = splint1(py, t1);
                a = sigma*dy / sqrt(dx*dx + dy*dy);
                rhox[ip][j + 1] -= a;
                rhox[ip1][j + 1] += a;
                break;
              }
            }

            /* Intersection with the horizontals */
            if (listi[k].type == HOR_TYPE) { 
              if (listi[k].coord == j)
                m = -1;
              else
                m = 1;
              
              /* update the intersections structure */
              PERIODIC(i, nx, ip);
#ifdef RECON
              if (inter[ip][j].n >= NSEC)
                printf("WARNING intersec_update: inter[%d][%d].n = %d\n",
                       ip, j, inter[ip][j].n);
              else
                intersec_update(&inter[ip][j], &inter[ip][j + m], o, l);
#endif
              /* Add the surface tension term and the pressure gradient
                 correction term */
              switch(sflag) {
              case TENSION_U:
                dx = splint1(px, t1); dy = splint1(py, t1);
                a = sigma*dx / sqrt(dx*dx + dy*dy);
                PERIODIC(i + 1, nx, ip);
                rhox[ip][j] -= a;
                rhox[ip][j + m] += a;
                break;
              case TENSION_V:
                dx = splint1(px, t1); dy = splint1(py, t1);
                a = sigma*dy / sqrt(dx*dx + dy*dy);
                rhox[ip][j + 1] -= a;
                rhox[ip][j + m + 1] += a;
                /* pressure gradient correction */
                a = x - (real) i;
                if (m == 1) {
                  if (a > 0.5) {
                    PERIODIC(i + 1, nx, ip1);             
                    a = (1.0 - a)*(p[ip1][j+1] - p[ip][j+1]);
                  }
                  else {
                    PERIODIC(i - 1, nx, ip1);
                    a = a*(p[ip1][j+1] - p[ip][j+1]);
                  }
                  rhox[ip][j + 1] -= a;
                  rhox[ip][j + 2] += a;
                } 
                else {
                  if (a > 0.5) {
                    PERIODIC(i + 1, nx, ip1);
                    a = (1.0 - a)*(p[ip1][j] - p[ip][j]);
                  }
                  else {
                    PERIODIC(i - 1, nx, ip1);
                    a = a*(p[ip1][j] - p[ip][j]);
                  }
                  rhox[ip][j + 1] += a;
                  rhox[ip][j] -= a;
                }
                break;
              }
            }

#ifdef AXI
            /* add the axisymetric contribution to the surface tension */
            switch(sflag) {
            case TENSION_U:
              PERIODIC(i + 1, nx, ip);
              rhox[ip][j] += sigma*stensionaxix(px, py, t1, t2);
              break;
            case TENSION_V:
              PERIODIC(i, nx, ip);
              rhox[ip][j+1] += sigma*stensionaxiy(px, py, t1, t2);
              break;
            }
#endif
            PERIODIC(i, nx, ip);            
            /* add the circulation along the polynomial segment */
            work[ip][j] += 
              t2*(a1 + t2*(a2 + t2*(a3 + t2*(a4 + t2*(a5 + t2*a6))))) -
              t1*(a1 + t1*(a2 + t1*(a3 + t1*(a4 + t1*(a5 + t1*a6)))));
            a = ip - i; /* we shift the polynom by a to insure the
                           periodicity of the circulation */
            if (a != 0.0) /* correction if periodic */
              work[ip][j] += a*(t2*(py.b + t2*(py.c + t2*py.d)) -
                                t1*(py.b + t1*(py.c + t1*py.d)));
#ifdef RECON        
            if (inter[ip][j].n == 0)
              inter[ip][j].n = -1;
#endif
          }
      }

  /* magic formula */
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) 
      {
        work[i][j] = modf(work[i][j], &y);
        if (work[i][j] >= c[i][j] + 0.5)
          work[i][j] -= 1.0;
        if (work[i][j] <= c[i][j] - 0.5)
          work[i][j] += 1.0;
        if (work[i][j] <= c[i][j] - 0.5)
          work[i][j] += 1.0;
        
        /* limiters used to avoid round-off error problems */
#ifdef LIMITERS
        c[i][j] = work[i][j] > 1.0 ? 1.0 : work[i][j] < 0.0 ? 0.0 : work[i][j];
#else
        c[i][j] = work[i][j];
#endif
      }
}

void interface_reconnect(interface *in, int nint, intersec **inter, 
                         int nx, int ny)
{
  int i, j, i1, i2;
  int j11, j12, j21, j22;
  real t11, t12, t21, t22;

  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      if (inter[i][j].n > 1 && inter[i][j].in[0] != inter[i][j].in[1]) {
        i1 = inter[i][j].in[0];
        i2 = inter[i][j].in[1];
        j11 = inter[i][j].j1[0];
        j12 = inter[i][j].j2[0];
        j21 = inter[i][j].j1[1];
        j22 = inter[i][j].j2[1];
        printf("inter[%d][%d].n: %d interfaces: %d %d (%d,%d)(%d,%d)\n", 
               i, j, inter[i][j].n, i1, i2, j11, j12, j21, j22);
      }
    }
}


real interface_inside(interface in, real x, real y)
{
  int i;
  real a, d, dmin = 1e30;
  for (i = 0; i < in.n - 1; i++) {
    d = sqrt(sq(x - 0.5*(in.x[i] + in.x[i+1])) + 
             sq(y - 0.5*(in.y[i] + in.y[i+1])));
    a = (in.x[i+1] - in.x[i])*(y - in.y[i]) - 
      (x - in.x[i])*(in.y[i+1] - in.y[i]);
    if (d < fabs(dmin)) dmin = a >= 0.0 ? d : -d;
  }
  if (dmin >= 0.0) return dmin;
  return 0.0;
}


void interface_fill(interface *in, int nin, real2D ap, int nx, int ny)
{
  int i, j, k;

  fill0(ap, nx, ny);
  for (i = 2; i <= nx; i++)
    for (j = 2; j <= ny; j++)
      for (k = 0; k < nin; k++)
        if (interface_inside(in[k], (real)i, (real)j)) {
          ap[i][j] += 0.25; ap[i-1][j] += 0.25;
          ap[i][j-1] += 0.25; ap[i-1][j-1] += 0.25;
        }
  for (j = 2; j <= ny; j++)
    for (k = 0; k < nin; k++)
      if (interface_inside(in[k], 1.0, (real)j)) {
        ap[1][j] += 0.25; ap[1][j-1] += 0.25;
      }
  for (i = 2; i <= nx; i++)
    for (k = 0; k < nin; k++)
      if (interface_inside(in[k], (real)i, 1.0)) {
        ap[i][1] += 0.25; ap[i-1][1] += 0.25;
      }
}


int interface_split_axis(interface **inlist, real length)
{
  interface *in = *inlist;
  int i;
  for (i = 2; i < in->n - 2; i++)
    if (in->y[i] - 2.0 < length) {
      /* new interface */
      interface *newin;
      int j;
      real x0, y0, x1, y1, x2, y2;
      *inlist = (interface *)realloc(*inlist, sizeof(interface)*2);
      in = &((*inlist)[0]);
      newin = &((*inlist)[1]);
      newin->n = in->n - i;
      newin->x = (real *)malloc(newin->n * sizeof(real));
      newin->y = (real *)malloc(newin->n * sizeof(real));
      newin->t = (real *)malloc(newin->n * sizeof(real));
      newin->spx = (polynom3 *)malloc((newin->n - 1) * sizeof(polynom3));
      newin->spy = (polynom3 *)malloc((newin->n - 1) * sizeof(polynom3));
      newin->x[0] = in->x[i+1];
      newin->y[0] = 1.999;
      for (j = i + 1; j < in->n; j++) {
        newin->x[j - i] = in->x[j];
        newin->y[j - i] = in->y[j];
      }
      newin->bc = AXIBOTH;
      printf("newin->n: %d i: %d\n", newin->n, i);
      /* trim old interface */
      in->x[i] = in->x[i-1];
      in->y[i] = 1.999;
      in->n = i + 1;
      in->bc = AXITRUE;
      {
        FILE *fptr = fopen("old", "wt");
        for (i = 0; i < in->n; i++)
          fprintf(fptr, "%g %g\n", in->x[i], in->y[i]);
        fclose(fptr);
        fptr = fopen("new", "wt");
        for (i = 0; i < newin->n; i++)
          fprintf(fptr, "%g %g\n", newin->x[i], newin->y[i]);
        fclose(fptr);
      }
      return 2;
    }
  return 1;
}


int interface_merge_axis(interface **inlist, real length)
{
  interface *right = *inlist, *left = &((*inlist)[1]);
  if (right->x[right->n - 1] - left->x[0] < length) {
    int i, n = right->n + left->n - 2;
    right->x = (real *)realloc(right->x, n * sizeof(real));
    right->y = (real *)realloc(right->y, n * sizeof(real));
    right->t = (real *)realloc(right->t, n * sizeof(real));
    right->spx = (polynom3 *)realloc(right->spx, (n - 1) * sizeof(polynom3));
    right->spy = (polynom3 *)realloc(right->spy, (n - 1) * sizeof(polynom3));
    for (i = 1; i < left->n; i++) {
      right->x[right->n + i - 2] = left->x[i];
      right->y[right->n + i - 2] = left->y[i];
    }
    right->n = n;
    right->bc = AXIBOTH;
    interface_free(*left);
    {
      FILE *fptr = fopen("merged", "wt");
      for (i = 0; i < right->n; i++)
        fprintf(fptr, "%g %g\n", right->x[i], right->y[i]);
      fclose(fptr);
    }
    return 1;
  }
  return 2;
}

/*
int interface_merge_axis2(interface **inlist, real length)
{
  interface *temp
  interface *right = *inlist, *left = &((*inlist)[1]);
  if (right->x[0] - left->x[0] < length) {
    int i, n = right->n + left->n - 2;
    temp->x = (real *)malloc(temp->x, n * sizeof(real));
    temp->y = (real *)malloc(temp->y, n * sizeof(real));
    temp->t = (real *)malloc(temp->t, n * sizeof(real));
    temp->spx = (polynom3 *)malloc(temp->spx, (n - 1) * sizeof(polynom3));
    temp->spy = (polynom3 *)malloc(temp->spy, (n - 1) * sizeof(polynom3));
    for (i = 0; i < left->n ; i++) {
      temp->x[i] = left->x[i];
      temp->y[i] = left->y[i];
    }
    right->n = n;
    right->bc = AXIBOTH;
    interface_free(*left);
    {
      FILE *fptr = fopen("merged", "wt");
      for (i = 0; i < right->n; i++)
      fprintf(fptr, "%g %g\n", right->x[i], right->y[i]);
      fclose(fptr);
      }
    return 1;
    }
    return 2;
}
*/




void interface_axisymmetry(interface *in, real space, int nx, int ny)
{
  int i, k, nb = 2 * (in->n - 1) , nold = in->n, n, nnew;
  real t, ts, te, ti[4], xe, xs, l, dl, temp;
  char old_bc = in->bc;

  if (in->x[2] != in->x[1]){
          in->x[0] = (sq(in->y[2]-2.)-sq(in->y[1]-2.)+
                          sq(in->x[2])-sq(in->x[1])) / 
                  (2. * (in->x[2]-in->x[1]));
          temp = in->x[0];
          if (in->x[1] > in->x[2])
          in->x[0] += sqrt(sq(in->y[1]-2.) + sq(in->x[1] - temp));
          else 
          in->x[0] -= sqrt(sq(in->y[1]-2.) + sq(in->x[1] - temp));
  }
  else in->x[0] = in->x[1];

  interface_splines(*in, nx, ny);

  interface_redis(in, space);

}



/*cld
void interface_axisymmetry(interface *in, real space, int nx, int ny)
{
  int i, k, nb = 2 * (in->n - 1) - 2, nold = in->n, n, nnew;
  real t, ts, te, ti[4], xe, xs, l, dl, xold, yold;
  real1D tempx, tempy;
  char old_bc = in->bc;

  if (in->bc == TRUE_PERIODIC || in->bc == APPROX_PERIODIC) {
    interface_filter(in, nx, ny);
    interface_redis(in, space);
    return;
  }

  xold = in->x[in->n-1];
  yold = in->y[in->n-1];

  in->x = (real *)realloc(in->x, nb * sizeof(real));
  in->y = (real *)realloc(in->y, nb * sizeof(real));
  in->t = (real *)realloc(in->t, nb * sizeof(real));
  in->spx = (polynom3 *)realloc(in->spx, (nb - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)realloc(in->spy, (nb - 1) * sizeof(polynom3));

  for (i = 1; i <= in->n - 2; i++) {
          t = (in->t[i] + in->t[i+1])/2.;
          in->x[i] = splint(in->spx[i], t);
          in->y[i] = splint(in->spy[i], t);
  }

  for (i = in->n - 2; i <= 2*in->n - 5; i++){
          in->x[i] = in->x[i - in->n + 3];
          in->y[i] = in->y[i - in->n + 3];
  }

  for (i = 0; i <= in->n - 3; i++){
          in->x[i] = in->x[nb-1-i];
          in->y[i] = 4. - in->y[nb-1-i];
  }

//  in->x[nb-1] = xold;
//  in->y[nb-1] = yold + 2. * (real)ny - 4.;

  in->n = nb;
  in->bc = AXIBOTH;

  interface_splines(*in, nx, ny);

#if 0
  
  n = polyroots(in->spy[nb-2], in->t[nb-2], in->t[nb-1], 
                ti, 2.0, INTERPREC);
  if (n != 1) {
    FILE *fptr = fopen("aximarkers", "wt");
    interface_write_markers(*in, fptr);
    fclose(fptr);
    printf("ERROR: interface_axisymmetry(): n = %d\n", n);    
    exit(1);
  }
  xs = splint(in->spx[nb-2], ti[0]);
  ts = ti[0] - in->t[nb-1]; 

#endif
  
  n = polyroots(in->spy[nold-3], in->t[nold-3], in->t[nold-2], 
                ti, (real)ny, INTERPREC);
  if (n != 1) {
    FILE *fptr = fopen("aximarkers", "wt");
    interface_write_markers(*in, fptr);
    fclose(fptr);
    printf("ERROR: interface_axisymmetry(): n = %d\n", n);
    exit(1);
  }
  xe = splint(in->spx[nold-3], ti[0]);
  te = ti[0];

/*  n = polyroots(in->spy[2*nold-2], in->t[2*nold-2], in->t[2*nold-1],
                  ti, 2.*(real)ny-2., INTERPREC);
  if (n != 1) {
          FILE *fptr = fopen("aximarkers", "wt");
          interface_write_markers(*in, fptr);
          fclose(fptr);
          printf("ERROR: interface_axisymmetry(): n = %d\n", n);
          exit(1);
  }
  xs = splint(in->spx[2*nold-2], ti[0]);
  ts = ti[0];
*/

/*  dl = te / (real)(nold-1); 
  if (dl >= 1.0)
    printf("WARNING interface_redis(): dl=%g\n", dl);

  in->x = (real *)realloc(in->x, nold * sizeof(real));
  in->y = (real *)realloc(in->y, nold * sizeof(real));

  l = te + dl;
  k = 0;

  for (i = 1; i <= nold - 2; i++, l += dl) { 
/* find the next interval [j;j+1] such as t[j] <= l <= t[j+1] 
    while(in->t[k+1] < l)
      k++;
    in->x[i] = splint(in->spx[k+nold-2], l);
    in->y[i] = splint(in->spy[k+nold-2], l);
  }
  
  in->x[0] = xs; in->y[0] = 1.9999;
  in->x[nold-1] = xe; in->y[nold-1] = (real)ny + 0.001;
  in->t = (real *)realloc(in->t, nold * sizeof(real));
  in->spx = (polynom3 *)realloc(in->spx, (nold - 1) * sizeof(polynom3));
  in->spy = (polynom3 *)realloc(in->spy, (nold - 1) * sizeof(polynom3));
#ifdef FREE
  in->p = (real *)realloc(in->p, nold * sizeof(real));
#endif
  in->n = nold;
  in->bc = old_bc;

}
cld*/

---