#ifndef VarWindingEnergy_H
#define VarWindingEnergy_H

/* Last edited on DATE TIME by USER */

#include <Energy.h>

typedef enum { VWE_Perimeter, VWE_Area } VarWindingEnergyKind_t;

typedef struct VarWindingEnergy_data_t 
{
  VarWindingEnergyKind_t kind; /* Kind ??? */
} VarWindingEnergy_data_t;

Energy_t VarWindingEnergy_new(void);
  /* The {VarWinding} 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).

    For the {Perimeter} formula, we define also {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|}.
    
    For the {Area} formula, we define instead {s[i] == ((r[i] ×
    r[i-1])}, and {S == sum{ s[i] }/SNorm}, where {SNorm} is a
    (different) normalization factor. Note that {s[i]} here is the
    signed area of the triangle {u,w[i-1],w[i]}, projected onto {P},
    scaled by {|h|^2}, and multiplied by the direction of {h}; so that
    {S} is the area of the polygon {w[0],.. w[N-1]}, also projected
    onto {P} and scaled by {|h|^2}.
    
    In either case, the energy contributed by @{edge->?} {e} is {(S/R +
    R/S)^2 - 4}. This formula has minimum value (zero) when {|R| ==
    |S|}, and tends to infinity when {|R/S| -> 0} or {|S/R| -> 0}.

    The normalization factors {RNorm} and {SNorm} are chosen so that
    {R,S,T} 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. 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, therefore
    
    | RNorm == 
    In either case, the quantities {R} and {S} will have the same value 
    
     {b^2|h|^2}
    
    In the {Perimeter} formula, .
    In the {Area} formula, the energy contributed to {e} is {T/R + R/T - 2}.
    Note that both formulas are not affected by changes of scale in the 
    {r[i]} or {h}.  

    The above defintions assume that {e} is an interior @{edge->?}.
    If {e} lies on the manifold's border, the {S} and {T} sums ignore the 
    `gaps' (@places {w[i], w[i-1]} that are separated by a null cell).
    In that case, the factors {SNorm,TNorm} are chosen to give 
    {S == T == b^2|h|^2} when the {w[i]} are {N} corners of a regular 
    {2(N-G)}-gon, where {G} is the number of gaps. */
  
#endif