# sandbox/ariviere/tools/gparabola.h

I slighly modified parabola.h of the main branch to add a fit for the gaussian curvature. *

#include "utils.h"

// Define this to use a x^iy^j polynomial with i = 0...NP-1, j = 0...NP-1
// #define NP 3

static double parabola_fit_gcurvature (ParabolaFit * p,
double gkappamax)
{
double kappa;
#if dimension == 2//just the mean curvature
double dnm = 1. + sq(p->a[1]);
kappa = - 2.*p->a[0]/pow(dnm, 3/2.);
//  if (kmax)
//    *kmax = fabs (kappa);
#else /* 3D */
# ifdef NP
double hxx = 2.*p->a[2*NP], hyy = 2.*p->a[2], hxy = p->a[NP + 1];
double hx = p->a[NP], hy = p->a[1];
# else
double hxx = 2.*p->a[0], hyy = 2.*p->a[1], hxy = p->a[2];
double hx = p->a[3], hy = p->a[4];
# endif
double dnm = 1. + sq(hx) + sq(hy);
kappa = (hxx*hyy - hxy*hxy)/(dnm*dnm);
//  kappa = - (hxx*(1. + sq(hy)) + hyy*(1. + sq(hx)) - 2.*hxy*hx*hy)
//    /sqrt (dnm*dnm*dnm); //where does the - come from?
//  if (kmin) {
//      double km = - (hxx*(1. + sq(hy)) + hyy*(1. + sq(hx)) - 2.*hxy*hx*hy)
//    /sqrt (dnm*dnm*dnm);
//      double a = km*km/4. - kappa;
//      *kmin = fabs (km/2.);
//    if (a >= 0.)
//      *kmax += sqrt (a);
//  }
#endif /* 3D */
if (fabs (kappa) > gkappamax) {
//    if (kmax)
//      *kmax = kappamax;
return kappa > 0. ? gkappamax : - gkappamax;
}
return kappa;
}

#if AXI
static void parabola_fit_axi_gcurvature (const ParabolaFit * p,
double r, double h,
double * kappa)
{
double nr = (p->m.x*p->a[1] + p->m.y)/sqrt (1. + sq(p->a[1]));
/* limit the minimum radius to half the grid size */
double kaxi = nr/max(r, h/2.);
*kappa *= kaxi;//check that
//  if (kmax)
//    *kmax = max (*kmax, fabs (kaxi));
}
#endif /* 2D */