#ifndef WindingEnergy_H
#define WindingEnergy_H

/* Last edited on DATE TIME by USER */

#include <Energy.h>

typedef struct WindingEnergy_data_t 
{
} WindingEnergy_data_t;

Energy_t VarWindingEnergy_new(void);
  /* The {Winding} energy penalizes @{edge->?}s whose star, when projected to
    {R^3}, is too irregular, or does not make a single complete turn
    around the edge's line. The projection to R^3 is performed by
    dropping the last coordinate.

    Let {e} be an @{edge->?} with endpoints {u,v}, and let {w[i]} (for {i}
    in {0..N-1}) be the third corner of the {i}th triangle incident to
    {e}, in topological order.  Let {L} be the line supporting {e},
    and {P} be a plane orthogonal to {L}.
    
    Define {h == v - u}, {r[i] == h × (w[i] - u)}. Note that {|r[i]|} is
    the distance from {w[i]} to the line {L} (projected onto the plane
    {P}), times {|h|}. Define also {R == sum{ |r[i]|^2 }/RNorm}, where
    {RNorm} is a normalization factor that depends only on the 
    local topology (see below).  Note that, for a fixed topology,
    {R} is proportional to the mean squared distance of the points
    {w[i]} to the line {L}, times {|h|^2}

    Define furthermore {s[i] == h × (w[i] - w[i-1])}, and {S == sum
    {
    |s[i]|^2 }/SNorm}, where {SNorm} is a normalization factor. Note
    that {|s[i]|} is the length of the edge {w[i] - w[i-1]}, projected
    onto {P} and scaled by {|h|}. Thus, for a fixed topology, {S} is
    proportional to the mean squared length of the sides of the
    polygon {w[0..N-1]}, projected onto {P}, times {|h|^2}.
    
    The energy contributed by @{edge->?} {e} is then defined as {S/R + R/S -
    2}. This formula has minimum value (zero) when {R == S}, and tends
    to infinity when {R/S -> 0} or {S/R -> 0}. Note, furthermore, that
    this energy is invariant with respect to the length of the edge
    {e}, to arbitrary displacement of the {w[i]} parallel to {L}, and
    to uniform scaling of the configuration.

    The normalization factors {RNorm} and {SNorm} are chosen so that
    {R} and {S} have the same value when the {w[i]}, projected onto {P},
    have an ideal configuration around the line {L}. If {e} is an
    interior @{edge->?}, the ideal configuration is a regular {N}-gon of
    arbitrary radius centered on the edge. The normalization 
    factors are then
    
    | RNorm == N
    | SNorm == 4·N·(sin(Pi/N))^2
    
    If {e} lies on the manifold's boundary, and there are {G} gaps
    (missing cells) in the ring of elements surrounding {e}, then the
    ideal configuration is a regular {2(N-G)}-gon, consisting of {N-G}
    `filled' sectors (corresponding to the non-missing cells) and
    {N-G} `empty' ones (corresponding to missing cells, or fractions
    thereof). The normalization factors are then
    
    | RNorm == N
    | SNorm == 4·(N-G)·(sin(Pi/2/(N-G)))^2
     */
  
#endif