/**
 * BasicCSTCurve -- An implementation of a Class-Shape-Transform planar curve.
 *
 * Copyright (C) 2014-2015, by Joseph A. Huwaldt. All rights reserved.
 *
 * This library is free software; you can redistribute it and/or modify it under the terms
 * of the GNU Lesser General Public License as published by the Free Software Foundation;
 * either version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
 * PARTICULAR PURPOSE. See the GNU Library General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place -
 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html
 */
package geomss.geom.cst;

import geomss.geom.*;
import geomss.geom.nurbs.*;
import jahuwaldt.js.math.BinomialCoef;
import jahuwaldt.js.param.Parameter;
import jahuwaldt.tools.math.MathTools;
import java.text.MessageFormat;
import java.util.List;
import java.util.Objects;
import static java.util.Objects.requireNonNull;
import javax.measure.converter.ConversionException;
import javax.measure.quantity.Dimensionless;
import javax.measure.quantity.Length;
import javax.measure.unit.SI;
import javax.measure.unit.Unit;
import javolution.context.ArrayFactory;
import javolution.context.ObjectFactory;
import javolution.context.StackContext;
import javolution.lang.ValueType;
import javolution.util.FastTable;
import javolution.xml.XMLFormat;
import javolution.xml.stream.XMLStreamException;

/**
 * An implementation of the Class-Shape-Transform (CST) planar curve. This curve type is
 * useful for defining smooth aircraft component shapes with a minimum number of defining
 * parameters.
 *
 * <p> Modified by: Joseph A. Huwaldt </p>
 *
 * @author Joseph A. Huwaldt, Date: March 14, 2014
 * @version November 28, 2015
 */
@SuppressWarnings({"serial", "CloneableImplementsClone"})
public final class BasicCSTCurve extends CSTCurve implements ValueType {
    // Reference: Kulfan, B.M., Bussoletti, J.E., "'Fundamental' Parametric Geometry
    //          Representations for Aircraft Component Shapes", AIAA-2006-6948.

    /**
     * The class function for this curve.
     */
    private CSTClassFunction Cf;

    /**
     * The shape function for this curve.
     */
    private CSTShapeFunction Sf;

    /**
     * The trailing edge offset for this curve.
     */
    private double offsetTE_c;

    /**
     * The origin of this curve (the location of the leading edge; gives the curve
     * location).
     */
    private Point origin;

    /**
     * The chord length of this curve (gives the curve scale).
     */
    private Parameter<Length> chord;

    /**
     * The X-direction for the curve (gives the curve direction).
     */
    private Vector<Dimensionless> xhat;

    /**
     * The Y-direction for the curve (gives the curve direction).
     */
    private Vector<Dimensionless> yhat;

    /**
     * A parametric mapping between the input "s" and the "s" used in calculations.
     */
    private NurbsCurve pmap;
    
    /**
     * The minimum bounding point for this subrange curve.
     */
    private Point _boundsMin;

    /**
     * The maximum bounding point for this subrange curve.
     */
    private Point _boundsMax;

    /**
     * Return a new planar {@link BasicCSTCurve} instance using the supplied input
     * parameters.
     *
     * @param Cf       The class function for this BasicCSTCurve. May not be null.
     * @param Sf       The shape function for this BasicCSTCurve. May not be null.
     * @param offsetTE The trailing edge offset from the chord-line for this
     *                 BasicCSTCurve. May not be null.
     * @param origin   The origin of the leading edge of this curve (gives the curve
     *                 location). May not be null.
     * @param chord    The chord length for this curve (gives the curve scale). May not be
     *                 null.
     * @param xhat     The chord-wise direction of this curve (the curve X-axis
     *                 direction). May not be null.
     * @param yhat     The direction perpendicular to xhat for this curve (the curve
     *                 Y-axis direction). May not be null.
     * @return A new BasicCSTCurve instance
     */
    public static BasicCSTCurve newInstance(CSTClassFunction Cf, CSTShapeFunction Sf, Parameter<Length> offsetTE,
            Point origin, Parameter<Length> chord, GeomVector<Dimensionless> xhat, GeomVector<Dimensionless> yhat) {
        requireNonNull(Cf);
        requireNonNull(Sf);
        requireNonNull(offsetTE);
        
        int numDims = GeomUtil.maxPhyDimension(origin, xhat, yhat);
        if (numDims < 2)
            throw new DimensionException(
                    MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"),
                            "input origin or normal", numDims));

        //  Convert everything to common units and dimensions.
        Unit<Length> unit = origin.getUnit();
        origin = origin.toDimension(numDims);
        chord = chord.to(unit);
        xhat = xhat.toDimension(numDims);
        yhat = yhat.toDimension(numDims);

        //  Create an instance of this object and store all the inputs.
        BasicCSTCurve obj = FACTORY.object();
        obj.pmap = CurveFactory.createLine(Point.valueOf(0), Point.valueOf(1));
        obj.Cf = Cf;
        obj.Sf = Sf;
        obj.offsetTE_c = offsetTE.divide(chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT);
        obj.origin = origin;
        obj.chord = chord;
        obj.xhat = xhat.toUnitVector();
        obj.yhat = yhat.toUnitVector();

        //  Compute and store the bounds min/max values.
        obj.calcBoundsMinMax();

        return obj;
    }

    /**
     * Return a new 2D, planar {@link BasicCSTCurve} instance using the supplied input
     * parameters with the origin at (0,0), a chord length of 1.0 m and with the chord
     * aligned with the X-axis.
     *
     * @param Cf     The class function for this BasicCSTCurve. May not be null.
     * @param Sf     The shape function for this BasicCSTCurve. May not be null.
     * @param dZTE_c The trailing edge offset from the chord-line to chord length ratio
     *               for this BasicCSTCurve.
     * @return A new BasicCSTCurve instance
     */
    public static BasicCSTCurve newInstance(CSTClassFunction Cf, CSTShapeFunction Sf, double dZTE_c) {
        Parameter<Length> chord = Parameter.valueOf(1, SI.METER);
        Parameter<Length> dZTE = chord.times(dZTE_c);
        BasicCSTCurve crv = BasicCSTCurve.newInstance(Cf, Sf, dZTE, Point.newInstance(2), chord,
                Vector.valueOf(1, 0), Vector.valueOf(0, 1));
        return crv;
    }

    /**
     * Return the class function associated with this CST curve.
     *
     * @return The class function associated with this CST curve.
     */
    @Override
    public CSTClassFunction getCFunction() {
        return Cf;
    }

    /**
     * Return the shape function associated with this CST curve.
     *
     * @return The shape function associated with this CST curve.
     */
    @Override
    public CSTShapeFunction getSFunction() {
        return Sf;
    }

    /**
     * Return the offset of the trailing edge of the curve from the chord line.
     *
     * @return The offset of the trailing edge of the curve from the chord line.
     */
    @Override
    public Parameter<Length> getOffsetTE() {
        return chord.times(offsetTE_c);
    }

    /**
     * Return the chord length of the curve. This defines the scale or size of the curve.
     *
     * @return The chord length of the curve.
     */
    @Override
    public Parameter<Length> getChord() {
        return chord;
    }

    /**
     * Return the origin of the leading edge of this curve. This defines the s=0 location
     * of the curve in space.
     *
     * @return The origin of the leading edge of this curve.
     */
    @Override
    public Point getOrigin() {
        return origin;
    }

    /**
     * Return the chord-wise direction of this curve (the curve X-direction).
     *
     * @return The chord-wise direction of this curve.
     */
    @Override
    public Vector<Dimensionless> getXhat() {
        return xhat;
    }

    /**
     * Return the in-plane direction of this curve perpendicular to the chord-wise
     * direction (the curve Y-direction).
     *
     * @return The in-plane direction of this curve perpendicular to Xhat.
     */
    @Override
    public Vector<Dimensionless> getYhat() {
        return yhat;
    }

    /**
     * Return a copy of this CSTCurve with the class function changed to the specified
     * class function.
     *
     * @param Cf The new class function for this BasicCSTCurve. May not be null.
     * @return A copy of this CSTCurve with the class function changed to the specified
     *         class function.
     */
    public BasicCSTCurve changeCFunction(CSTClassFunction Cf) {
        requireNonNull(Cf);
        BasicCSTCurve obj = copy();
        obj.Cf = Cf;
        obj.calcBoundsMinMax();
        return obj;
    }

    /**
     * Return a copy of this CSTCurve with the shape function changed to the specified
     * shape function.
     *
     * @param Sf The new shape function for this BasicCSTCurve. May not be null.
     * @return A copy of this CSTCurve with the shape function changed to the specified
     *         shape function.
     */
    public BasicCSTCurve changeSFunction(CSTShapeFunction Sf) {
        requireNonNull(Sf);
        BasicCSTCurve obj = copy();
        obj.Sf = Sf;
        obj.calcBoundsMinMax();
        return obj;
    }

    /**
     * Return a copy of this CSTCurve with the shape function Bernstein Polynomial
     * coefficients changed to the specified coefficients.
     *
     * @param order  The order of the shape function to use.
     * @param Acoefs The new coefficients of the Bernstein Polynomial used to construct
     *               the shape function. If more than "order" coefficients are provided
     *               then the additional coefficients are ignored. May not be null.
     * @return A copy of this CSTCurve with the shape function Bernstein Polynomial
     *         coefficients changed to the specified coefficients.
     */
    public BasicCSTCurve changeBFCoefficients(int order, double... Acoefs) {
        CSTShapeFunction Sfunc = CSTShapeFunction.newInstance(order, requireNonNull(Acoefs));
        return changeSFunction(Sfunc);
    }

    /**
     * Return a copy of this CSTCurve with the trailing edge offset from the chord-line
     * changed to the specified value.
     *
     * @param offsetTE The new trailing edge offset from the chord-line for this
     *                 BasicCSTCurve. May not be null.
     * @return A copy of this CSTCurve with the trailing edge offset changed to the
     *         specified value.
     */
    public BasicCSTCurve changeOffsetTE(Parameter<Length> offsetTE) {
        BasicCSTCurve obj = copy();
        obj.offsetTE_c = offsetTE.divide(obj.chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT);
        obj.calcBoundsMinMax();
        return obj;
    }
    
    /**
     * Return a copy of this CSTCurve with the trailing edge offset from the chord-line to
     * chord ratio changed to the specified value.
     *
     * @param offsetTE_c The new trailing edge offset from the chord-line to chord ratio
     *                   for this BasicCSTCurve.
     * @return A copy of this CSTCurve with the trailing edge offset changed to the
     *         specified value.
     */
    public BasicCSTCurve changeOffsetTE(double offsetTE_c) {
        BasicCSTCurve obj = copy();
        obj.offsetTE_c = offsetTE_c;
        obj.calcBoundsMinMax();
        return obj;
    }
    
    /**
     * Return a copy of this CSTCurve with the curve origin or leading-edge point changed
     * to the specified value.
     *
     * @param origin The new origin of the leading edge of this curve (gives the curve
     *               location). May not be null.
     * @return A copy of this CSTCurve with the origin changed to the specified value.
     */
    public BasicCSTCurve changeOrigin(Point origin) {
        requireNonNull(origin);
        BasicCSTCurve obj = copy();
        obj.origin = origin;
        obj.calcBoundsMinMax();
        return obj;
    }
    
    /**
     * Return a copy of this CSTCurve with the chord length changed to the specified
     * value.
     *
     * @param chord The chord length for this curve (gives the curve scale). May not be null.
     * @return A copy of this CSTCurve with the chord length changed to the specified
     *         value.
     */
    public BasicCSTCurve changeChord(Parameter<Length> chord) {
        requireNonNull(chord);
        BasicCSTCurve obj = copy();
        obj.chord = chord;
        obj.calcBoundsMinMax();
        return obj;
    }
    
    /**
     * Returns the number of physical dimensions of the geometry element.
     *
     * @return The number of physical dimensions of the geometry element.
     */
    @Override
    public int getPhyDimension() {
        return origin.getPhyDimension();
    }

    /**
     * Return the order of this CST curve.
     *
     * @return order of this curve (degree + 1)
     *
     */
    public int getOrder() {
        return Sf.getOrder();
    }

    /**
     * Return the value of this CST curve function at the specified parametric location.
     */
    private double getValue(double s) {
        double S = Sf.getValue(s);
        double C = Cf.getValue(s);
        return C * S + offsetTE_c * s;
    }

    /**
     * Calculate a point on the curve for the given parametric distance along the curve,
     * <code>p(s)</code>.
     *
     * @param s parametric distance to calculate a point for (0.0 to 1.0 inclusive).
     * @return the calculated point
     */
    @Override
    public Point getRealPoint(double s) {
        validateParameter(s);

        StackContext.enter();
        try {
            //  Map "s" from parameter space to CST space.
            s = pmap.getRealPoint(s).getValue(Point.X);

            double ch = chord.getValue(SI.METER);
            double aValue = ch * s;
            double bValue = ch * getValue(s);
            Point p = Point.valueOf(xhat.times(aValue).plus(yhat.times(bValue))).to(getUnit());
            p = p.plus(origin);

            return StackContext.outerCopy(p);
        } finally {
            StackContext.exit();
        }
    }

    /**
     * Calculate all the derivatives from <code>0</code> to <code>grade</code> with
     * respect to parametric distance on the curve for the given parametric distance along
     * the curve, <code>d^{grade}p(s)/d^{grade}s</code>.
     * <p>
     * Example:<br>
     * 1st derivative (grade = 1), this returns <code>[p(s), dp(s)/ds]</code>;<br>
     * 2nd derivative (grade = 2), this returns <code>[p(s), dp(s)/ds, d^2p(s)/d^2s]</code>; etc.
     * </p>
     *
     * @param s     Parametric distance to calculate derivatives for (0.0 to 1.0
     *              inclusive).
     * @param grade The maximum grade to calculate the derivatives for (1=1st derivative,
     *              2=2nd derivative, etc)
     * @return A list of derivatives up to the specified grade of the curve at the
     *         specified parametric position.
     * @throws IllegalArgumentException if the grade is &lt; 0.
     */
    @Override
    public List<Vector<Length>> getSDerivatives(double s, int grade) {
        validateParameter(s);
        s = roundParNearEnds(s);

        if (grade < 0)
            throw new IllegalArgumentException(RESOURCES.getString("gradeLTZeroErr"));

        //  CST Terms:  1: C(s)*S(s), 2: s*offsetTE_c
        Vector[] ders = Vector.allocateArray(grade + 1);   //  Created outside of StackContext.
        StackContext.enter();
        try {
            //  Map "s" from user space to calculation space.
            s = pmap.getRealPoint(s).getValue(Point.X);

            //  Get derivatives of the class and shape functions.
            double[] Cders = Cf.getDerivatives(s, grade);
            double[] Sders = Sf.getDerivatives(s, grade);

            //  Calculate components of the 1st term derivatives using Leibniz's Rule.
            double[] term1 = leibnizRule(Cders, Sders, grade);

            GeomVector<Length> o = origin.toGeomVector();
            Unit<Length> unit = getUnit();
            for (int n = 0; n <= grade; ++n) {
                //  Get the derivatives of the 1st term.
                double v = term1[n];

                //  Add in the derivatives of the 2nd term (and of "s").
                double ds = 0;  //  Derivate of "s".
                if (n == 0) {
                    v += s * offsetTE_c;
                    ds = s;
                } else if (n == 1) {
                    v += offsetTE_c;
                    ds = 1;
                }
                //  2nd term derivatives are 0 for all higher grades.

                Vector<Dimensionless> vhat = xhat.times(ds).plus(yhat.times(v));
                Vector<Length> vec = vhat.times(chord).to(unit);
                if (n == 0)
                    vec = vec.plus(o);
                ders[n] = StackContext.outerCopy(vec);
            }
        } finally {
            StackContext.exit();
        }

        //  Convert to a list of Vector objets.
        FastTable<Vector<Length>> output = FastTable.newInstance();
        Point o = Point.valueOf(ders[0]);
        for (int i = 0; i <= grade; ++i) {
            Vector<Length> v = ders[i];
            v.setOrigin(o);
            output.add(v);
        }

        Vector.recycleArray(ders);

        return output;
    }

    /**
     * Compute the derivatives of the product of two functions using the Leibniz General
     * Rule for Differentiation: <code>d^n(f(s)*g(s))/ds^n</code>.
     *
     * @param f     The derivatives of the function "f" from 0 up to grade.
     * @param g     The derivatives of the function "g" from 0 up to grade.
     * @param grade The grade to which the derivative terms are to be calculated.
     * @return An array of the derivatives, from 0 to grade, of the product of two
     *         functions: d^n(f(s)*g(s))/ds^n. The returned array was allocated using
     *         javolution.context.ArrayFactory.DOUBLES_FACTORY, could be longer than the
     *         number of derivatives requested, and could be recycled by the user when no
     *         longer needed.
     */
    private static double[] leibnizRule(double[] f, double[] g, int grade) {
        //  Reference: https://en.wikipedia.org/wiki/General_Leibniz_rule

        int gradeP1 = grade + 1;
        BinomialCoef bin = BinomialCoef.newInstance(gradeP1);
        double[] fgn = ArrayFactory.DOUBLES_FACTORY.array(gradeP1);
        for (int n = 0; n < gradeP1; ++n) {
            for (int k = n; k >= 0; --k) {
                double bCoef = bin.get(n, k);
                double fk = f[k];
                double gnmk = g[n - k];
                double der = fk * gnmk * bCoef;
                fgn[n] += der;
            }
        }
        BinomialCoef.recycle(bin);

        return fgn;
    }

    /**
     * Return the coordinate point representing the minimum bounding box corner of this
     * geometry element (e.g.: min X, min Y, min Z).
     *
     * @return The minimum bounding box coordinate for this geometry element.
     * @throws IndexOutOfBoundsException if this list contains no geometry.
     */
    @Override
    public Point getBoundsMin() {
        return _boundsMin;
    }

    /**
     * Return the coordinate point representing the maximum bounding box corner (e.g.: max
     * X, max Y, max Z).
     *
     * @return The maximum bounding box coordinate for this geometry element.
     * @throws IndexOutOfBoundsException if this list contains no elements.
     */
    @Override
    public Point getBoundsMax() {
        return _boundsMax;
    }

    private static final int NPTS = 51;

    /**
     * Calculate the minimum & maximum bounding box corner points of this geometry
     * element. This method may ONLY be called at the time that an instance of this object
     * is created by a "newInstance()" method.
     */
    private void calcBoundsMinMax() {
        StackContext.enter();
        try {
            //  Space out some points on the curve.
            List<Double> spacing = GridSpacing.linear(NPTS);
            PointString<SubrangePoint> str = this.extractGrid(GridRule.PAR, spacing);

            //  Get the bounds of the grid of points.
            _boundsMin = StackContext.outerCopy(str.getBoundsMin());
            _boundsMax = StackContext.outerCopy(str.getBoundsMax());
        } finally {
            StackContext.exit();
        }
    }

    /**
     * Return a new curve that is identical to this one, but with the parameterization
     * reversed. This does not change the origin or axis directions, but only the curve
     * parameterization; what was s=0 will become s=1, etc.
     *
     * @return A new curve identical to this one, but with the parameterization reversed.
     */
    @Override
    public CSTCurve reverse() {
        BasicCSTCurve obj = copy();
        obj.pmap = this.pmap.reverse();
        return obj;
    }
    
    /**
     * Split this curve at the specified parametric position returning a list containing
     * two curves (a lower curve with smaller parametric positions than "s" and an upper
     * curve with larger parametric positions). The origin, chord length, and trailing
     * edge offset from chord-line will be identical for the the returned segments and
     * this curve.  The only change is in the parameterization.
     *
     * @param s The parametric position where this curve should be split (must not be 0 or
     *          1!).
     * @return A list containing two curves: 0 == the lower curve, 1 == the upper curve.
     */
    @Override
    public GeomList<CSTCurve> splitAt(double s) {
        validateParameter(s);
        if (parNearEnds(s, TOL_S))
            throw new IllegalArgumentException(
                    MessageFormat.format(RESOURCES.getString("canNotSplitAtEnds"), "curve"));
        
        GeomList<CSTCurve> output = GeomList.newInstance();
        
        //  Split the parametric mapping at "s".
        GeomList<NurbsCurve> p_split = pmap.splitAt(s);
        
        //  Create a "lower" curve segment.
        BasicCSTCurve cl = copy();
        cl.pmap = p_split.getFirst();
        cl.calcBoundsMinMax();
        BasicCSTCurve cu = copy();
        cu.pmap = p_split.getLast();
        cu.calcBoundsMinMax();
        
        output.add(cl);
        output.add(cu);
        
        return output;
    }
    
    /**
     * Return <code>true</code> if this curve contains valid and finite numerical
     * components. A value of <code>false</code> will be returned if any of the coordinate
     * values are NaN or Inf.
     *
     * @return true if this element contains valid and finite numerical values.
     */
    @Override
    public boolean isValid() {
        return origin.isValid() && xhat.isValid() && yhat.isValid()
                && !chord.isInfinite() && !chord.isNaN()
                && !Double.isInfinite(offsetTE_c) && !Double.isNaN(offsetTE_c);
    }

    /**
     * Returns the unit in which this curve is stated.
     *
     * @return The unit that this curves points are stated in.
     */
    @Override
    public Unit<Length> getUnit() {
        return origin.getUnit();
    }

    /**
     * Returns the equivalent to this curve but stated in the specified unit.
     *
     * @param unit the length unit of the curve to be returned. May not be null.
     * @return an equivalent to this curve but stated in the specified unit.
     * @throws ConversionException if the the input unit is not a length unit.
     */
    @Override
    public BasicCSTCurve to(Unit<Length> unit) throws ConversionException {
        if (unit.equals(getUnit()))
            return this;
        BasicCSTCurve crv = newInstance(Cf, Sf, getOffsetTE(), origin.to(unit), chord.to(unit), xhat, yhat);
        return copyState(crv);
    }

    /**
     * Return the equivalent of this curve converted to the specified number of physical
     * dimensions. If the number of dimensions is greater than this element, then zeros
     * are added to the additional dimensions. If the number of dimensions is less than
     * this element, then the extra dimensions are simply dropped (truncated). If the new
     * dimensions are the same as the dimension of this element, then this element is
     * simply returned.
     *
     * @param newDim The dimension of the curve to return.
     * @return The equivalent of this curve converted to the new dimensions.
     */
    @Override
    public BasicCSTCurve toDimension(int newDim) {
        if (getPhyDimension() == newDim)
            return this;
        BasicCSTCurve crv = newInstance(Cf, Sf, getOffsetTE(), origin.toDimension(newDim), chord,
                xhat.toDimension(newDim), yhat.toDimension(newDim));
        return copyState(crv);
    }

    /**
     * Return a NURBS curve representation of this curve to within the specified
     * tolerance.
     *
     * @param tol The greatest possible difference between this curve and the NURBS
     *            representation returned. May not be null.
     * @return A NURBS curve that represents this curve to within the specified tolerance.
     */
    @Override
    public NurbsCurve toNurbs(Parameter<Length> tol) {
        StackContext.enter();
        try {
            //  Can this curve be represented exactly by a NURBS curve?
            //  Only integer N2 (either 0 or 1) can be represented exactly as NURBS.
            if (MathTools.isApproxEqual(Cf.getN2(), 1.0)) {
                if (MathTools.isApproxEqual(Cf.getN1(), 0.5)) {
                    BasicNurbsCurve curve = bluntAirfoil2Bezier();
                    curve = matchParameterization(curve);
                    copyState(curve);
                    return StackContext.outerCopy(curve);
                }
                if (MathTools.isApproxEqual(Cf.getN1(), 1.0)) {
                    BasicNurbsCurve curve = sharpAirfoil2Bezier();
                    curve = matchParameterization(curve);
                    copyState(curve);
                    return StackContext.outerCopy(curve);
                }
                
            } else if (MathTools.isApproxEqual(Cf.getN2(), 0)) {
                if (MathTools.isApproxEqual(Cf.getN1(), 1.0)) {
                    BasicNurbsCurve curve = cone2Bezier();
                    curve = matchParameterization(curve);
                    copyState(curve);
                    return StackContext.outerCopy(curve);
                }
                if (MathTools.isApproxEqual(Cf.getN1(), 0.0)) {
                    BasicNurbsCurve curve = rectangular2Bezier();
                    curve = matchParameterization(curve);
                    copyState(curve);
                    return StackContext.outerCopy(curve);
                }
            }

            //  Fall-back to gridding points onto the curve and then fitting a curve through them.
            
            //  Grid some points onto the curve.
            PointString<SubrangePoint> str = this.gridToTolerance(requireNonNull(tol));

            //  Fit a cubic NURBS curve to the points.
            int deg = Sf.getOrder() - 1;
            if (deg < 3)
                deg = 3;
            if (str.size() <= deg)
                deg = str.size() - 1;
            BasicNurbsCurve curve = CurveFactory.fitPoints(deg, str);
            copyState(curve);   //  Copy over the super-class state for this curve to the new one.

            return StackContext.outerCopy(curve);

        } finally {
            StackContext.exit();
        }
    }

    /**
     * Modify the input NurbsCurve to have the same parameterization as this CST curve.
     * 
     * @param crv The NURBS curve to be modified.
     * @return The modified NURBS curve.
     */
    private BasicNurbsCurve matchParameterization(BasicNurbsCurve crv) {
        //  Get the limits of this curve's parameter mapping.
        double s0 = pmap.getRealPoint(0).getValue(0);
        double s1 = pmap.getRealPoint(1).getValue(0);

        //  Is the parameterization reversed?
        boolean reversed = s1 < s0;

        if (!reversed) {
            if (!parNearStart(s0, TOL_S)) {
                //  s0 is not at the start of the curve.
                
                //  Find s0 on the NurbsCurve and split the NurbsCurve there.
                SubrangePoint p = crv.getClosest(this.getRealPoint(0), s0, TOL_S);
                double s_split = p.getParPosition().getValue(0);
                GeomList<NurbsCurve> lst = crv.splitAt(s_split);
                crv = (BasicNurbsCurve)lst.getLast();
            }
            if (!parNearEnd(s1, TOL_S)) {
                //  s1 is not at the end of the curve.
                
                //  Find s1 on the NurbsCurve and split the NurbsCurve there.
                SubrangePoint p = crv.getClosest(this.getRealPoint(1), s1, TOL_S);
                double s_split = p.getParPosition().getValue(0);
                GeomList<NurbsCurve> lst = crv.splitAt(s_split);
                crv = (BasicNurbsCurve)lst.getFirst();
            }

        } else {
            //  The parameterization is reversed.
            if (!parNearStart(s1, TOL_S)) {
                //  s01 is not at the start of the curve.
                
                //  Find s1 on the NurbsCurve and split the NurbsCurve there.
                SubrangePoint p = crv.getClosest(this.getRealPoint(1), s1, TOL_S);
                double s_split = p.getParPosition().getValue(0);
                GeomList<NurbsCurve> lst = crv.splitAt(s_split);
                crv = (BasicNurbsCurve)lst.getLast();
            }
            if (!parNearEnd(s0, TOL_S)) {
                //  s0 is not at the end of the curve.
                
                //  Find s0 on the NurbsCurve and split the NurbsCurve there.
                SubrangePoint p = crv.getClosest(this.getRealPoint(0), s0, TOL_S);
                double s_split = p.getParPosition().getValue(0);
                GeomList<NurbsCurve> lst = crv.splitAt(s_split);
                crv = (BasicNurbsCurve)lst.getFirst();
            }
            crv = (BasicNurbsCurve)crv.reverse();
        }

        return crv;
    }
    
    /**
     * An exact conversion to a Bezier curve is possible for a blunt-nosed airfoil CST curve.
     * 
     * @return A NURBS curve that represents this airfoil curve exactly.
     */
    private BasicNurbsCurve bluntAirfoil2Bezier() {
        //  Ref.: Marshall, D.D., "Creating Exact Bezier Representations of CST Shapes",
        //  21st AIAA Computational Fluid Dynamics Conference Proceedings: San Diego, CA, June 24, 2013.
        
        int n = Sf.getOrder() - 1;      //  Degree of the shape function.
        int m = 2*n + 3;                //  Degree of output NURBS curve.
        
        //  Convert the Bezier curve shape function to a monomial of the form: S(t) = sum_{i=0}^{n} ai*t^i
        BinomialCoef bin = BinomialCoef.newInstance(m + 1);   //  Get a set of binomial coefficients.
        double[] a = shape2Monomial(Sf.getBasisFunction(), bin);
        
        //  Create the "c" coefficients for the new monomial using Equation 13 from the reference.
        //  xi = s^2 due to parameter substitution (see below), so c[2] = 1.
        double[] c = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        c[0] = 0;
        c[1] = 0;
        c[2] = 1;
        for (int i=3; i <= m; ++i)
            c[i] = 0;
        
        //  Create the "b" coefficients for the new monomial using Equation 11 from the reference.
        //  zeta = sqrt(t)*(1-t)*sum_{i=0}^{n}{ai*t^i} + t*DzetaTE
        //      With parameter substitution:  t ==> s^2:
        //  zeta = s*[a0 + sum_{i=1}^{n}{(ai - a_(i-1))*s^(2*i) - an*s^(2*(n+1))] + s^2*DzetaTE; or simplifying:
        //  zeta = a0*s + DzetaTE*s^2 + sum_{i=1}^{n}{ (ai - a_(i-1)) * s^(2*i + 1) } - an*s^m = sum_{i=0}^{n}{bi*t^i}
        double[] b = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        for (int i = 0; i <= m; ++i) {
            if (i == 1)
                b[i] = a[0];
            else if (i == 2)
                b[i] = offsetTE_c;
            else if (MathTools.even(i))
                b[i] = 0;
            else if (i == m)
                b[i] = -a[n];
            else {
                int i2 = (i - 1) / 2;
                b[i] = a[i2] - a[i2 - 1];
            }
        }

        //  Form the new Bezier curve from the newly computed monomial coefficients.
        BasicNurbsCurve newCrv = monomial2Bezier(c,b, m, bin);
        
        return newCrv;
    }
    
    /**
     * An exact conversion to a Bezier curve is possible for a sharp-nosed airfoil CST curve.
     * 
     * @return A NURBS curve that represents this airfoil curve exactly.
     */
    private BasicNurbsCurve sharpAirfoil2Bezier() {
        //  Ref.: Marshall, D.D., "Creating Exact Bezier Representations of CST Shapes",
        //  21st AIAA Computational Fluid Dynamics Conference Proceedings: San Diego, CA, June 24, 2013.
        
        int n = Sf.getOrder() - 1;      //  Degree of the shape function.
        int m = n + 2;                  //  Degree of output NURBS curve.
        
        //  Convert the Bezier curve shape function to a monomial of the form: S(t) = sum_{i=0}^{n} ai*t^i
        BinomialCoef bin = BinomialCoef.newInstance(m + 1);   //  Get a set of binomial coefficients.
        double[] a = shape2Monomial(Sf.getBasisFunction(), bin);
        
        //  Create the "c" coefficients for the new monomial.
        //  xi = t, to c[1] = 1.
        double[] c = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        c[0] = 0;
        c[1] = 1;
        for (int i=2; i <= m; ++i)
            c[i] = 0;
        
        //  Create the "b" coefficients for the new monomial.
        //  zeta = t*(1 - t)*sum_{i=0}^{n}{ai*t^i} + t*DzetaTE (Eqn. 14 with xi = t)
        //  zeta = t*[a0 + sum_{i=1}^{n}{(ai - a_(i-1))*t^i - an*s^(n+1)] + t*DzetaTE; simplifying:
        //  zeta = (a0 + DzetaTE)*t + sum_{i=1}^{n}{ (ai - a_(i-1))*t^(i+1) } - an*t^m = sum_{i=0}^{n}{bi*t^i}
        double[] b = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        for (int i = 0; i <= m; ++i) {
            if (i == 0)
                b[i] = 0;
            else if (i == 1)
                b[i] = a[0] + offsetTE_c;
            else if (i == m)
                b[i] = -a[n];
            else {
                int i2 = i - 1;
                b[i] = a[i2] - a[i2 - 1];
            }
        }

        //  Form the new Bezier curve from the newly computed monomial coefficients.
        BasicNurbsCurve newCrv = monomial2Bezier(c,b, m, bin);
        
        return newCrv;
    }
    
    /**
     * An exact conversion to a Bezier curve is possible for a cone or wedge shaped CST curve.
     * 
     * @return A NURBS curve that represents this CST curve exactly.
     */
    private BasicNurbsCurve cone2Bezier() {
        //  Ref.: Marshall, D.D., "Creating Exact Bezier Representations of CST Shapes",
        //  21st AIAA Computational Fluid Dynamics Conference Proceedings: San Diego, CA, June 24, 2013.
        
        int n = Sf.getOrder() - 1;      //  Degree of the shape function.
        int m = n + 1;                  //  Degree of output NURBS curve.
        
        //  Convert the Bezier curve shape function to a monomial of the form: S(t) = sum_{i=0}^{n} ai*t^i
        BinomialCoef bin = BinomialCoef.newInstance(m + 1);   //  Get a set of binomial coefficients.
        double[] a = shape2Monomial(Sf.getBasisFunction(), bin);
        
        //  Create the "c" coefficients for the new monomial.
        //  xi = t, to c[1] = 1.
        double[] c = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        c[0] = 0;
        c[1] = 1;
        for (int i=2; i <= m; ++i)
            c[i] = 0;
        
        //  Create the "b" coefficients for the new monomial.
        //  zeta = t*sum_{i=0}^{n}{ai*t^i} + t*DzetaTE (Eqn. 15 with xi = t)
        //  zeta = (a0 + DzetaTE)*t + sum_{i=1}^{n}{ai*t^(i+1)} + an*t^m = sum_{i=0}^{n}{bi*t^i}
        double[] b = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        for (int i = 0; i <= m; ++i) {
            if (i == 0)
                b[i] = 0;
            else if (i == 1)
                b[i] = a[0] + offsetTE_c;
            else if (i == m)
                b[i] = a[n];
            else
                b[i] = a[i];
        }

        //  Form the new Bezier curve from the newly computed monomial coefficients.
        BasicNurbsCurve newCrv = monomial2Bezier(c,b, m, bin);
        
        return newCrv;
    }
    
    /**
     * An exact conversion to a Bezier curve is possible for a rectangular shaped CST curve.
     * 
     * @return A NURBS curve that represents this CST curve exactly.
     */
    private BasicNurbsCurve rectangular2Bezier() {
        //  Ref.: Marshall, D.D., "Creating Exact Bezier Representations of CST Shapes",
        //  21st AIAA Computational Fluid Dynamics Conference Proceedings: San Diego, CA, June 24, 2013.
        
        int n = Sf.getOrder() - 1;      //  Degree of the shape function.
        int m = n;                      //  Degree of output NURBS curve.
        
        //  Convert the Bezier curve shape function to a monomial of the form: S(t) = sum_{i=0}^{n} ai*t^i
        BinomialCoef bin = BinomialCoef.newInstance(m + 1);   //  Get a set of binomial coefficients.
        double[] a = shape2Monomial(Sf.getBasisFunction(), bin);
        
        //  Create the "c" coefficients for the new monomial.
        //  xi = t, to c[1] = 1.
        double[] c = ArrayFactory.DOUBLES_FACTORY.array(m + 1);
        c[0] = 0;
        c[1] = 1;
        for (int i=2; i <= m; ++i)
            c[i] = 0;
        
        //  Create the "b" coefficients for the new monomial.
        //  zeta = sum_{i=0}^{n}{ai*t^i} + t*DzetaTE (Eqn. 16 with xi = t)
        //  zeta = a0 + (a1 + DzetaTE)*t + sum_{i=2}^{n}{ai*t^i} = sum_{i=0}^{n}{bi*t^i}
        double[] b = a;
        b[1] += offsetTE_c;

        //  Form the new Bezier curve from the newly computed monomial coefficients.
        BasicNurbsCurve newCrv = monomial2Bezier(c,b, m, bin);
        
        return newCrv;
    }
    
    /**
     * Method that converts the 1D shape function Bezier curve to a monomial
     * representation of the form: <code>S(t) = sum_{i=0}^{n}{ai*t^i}</code>.
     *
     * @param S   The Bezier curve that is the basis for the shape function.
     * @param bin A set of binomial coefficients of at least the S.getOrder() in size.
     * @return The coefficients of a monomial representation of the shape function.
     */
    private static double[] shape2Monomial(BasicNurbsCurve S, BinomialCoef bin) {
        
        //  Convert the Bezier curve shape function to a monomial of the form:
        //      S(t) = sum_{i=0}^{n} {ai*t^i}
        //  Equation 7 from reference: aj = sum_{i=0}^{j}{ (n,j) * (j,i) * (-1)^(j-i) * pi }
        int n = S.getDegree();
        List<ControlPoint> cps = S.getControlPoints();
        double[] a = ArrayFactory.DOUBLES_FACTORY.array(n + 1);
        for (int j = 0; j <= n; ++j) {
            a[j] = 0;
            for (int i = 0; i <= j; ++i) {
                double Bnj = bin.get(n, j);
                double Bji = bin.get(j, i);
                double sign = (MathTools.even(j - i) ? 1.0 : -1.0);   //  = (-1)^(j-i)
                double pi = cps.get(i).getPoint().getValue(Point.X);
                a[j] += Bnj * Bji * sign * pi;
            }
        }
        
        return a;
    }
    
    /**
     * Form a 2D Bezier curve (as a NURBS curve) from the input list of monomial
     * coefficients of degree "m".
     *
     * @param c   An array of monomial coefficients (c.length > m) for the "X" dimension.
     * @param b   An array of monomial coefficients (b.length > m) for the "Y" dimension.
     * @param m   The degree of the monomial equation.
     * @param bin A set of binomial coefficients of at least size m+1.
     * @return A new 2D Bezier curve representation of the input monomial equation
     *         coefficients.
     */
    private BasicNurbsCurve monomial2Bezier(double[] c, double[] b, int m, BinomialCoef bin) {
        
        //  Form the new control points using Equation 12 from the reference:
        //  qveci = 1/(m,n)*sum_{j=0}^{i}{ (m-j, i-j) * dvecj }
        //  dvecj = [cj, bj]
        FastTable<ControlPoint> cpLst = FastTable.newInstance();
        for (int i = 0; i <= m; ++i) {
            double invBmi = 1. / bin.get(m, i);
            Point qi = Point.valueOf(0, 0);
            for (int j = 0; j <= i; ++j) {
                double B = invBmi * bin.get(m - j, i - j);
                Point qij = Point.valueOf(c[j]*B, b[j]*B);
                qi = qi.plus(qij);
            }
            double ch = chord.getValue(SI.METER);
            double qixh = ch*qi.getValue(0);    //  qi in the xhat direction.
            double qiyh = ch*qi.getValue(1);    //  qi in the yhat direction.
            Point cpp = Point.valueOf(xhat.times(qixh).plus(yhat.times(qiyh))).to(getUnit());
            ControlPoint cpi = ControlPoint.valueOf(cpp, 1);
            cpLst.add(cpi);
        }
        
        //  Create a Bezier curve knot vector.
        KnotVector bezierKV = KnotVector.bezierKnotVector(m);
        
        //  Construct the new Bezier curve as a NURBS curve.
        BasicNurbsCurve newCrv = BasicNurbsCurve.newInstance(cpLst, bezierKV);
        
        return newCrv;
    }
    
    /**
     * Returns a copy of this {@link BasicCSTCurve} instance
     * {@link javolution.context.AllocatorContext allocated} by the calling thread
     * (possibly on the stack).
     *
     * @return an identical and independent copy of this curve.
     */
    @Override
    public BasicCSTCurve copy() {
        return copyOf(this);
    }

    /**
     * Return a copy of this object with any transformations or subranges removed
     * (applied).
     *
     * @return A copy of this object with any transformations or subranges removed
     *         (applied).
     */
    @Override
    public BasicCSTCurve copyToReal() {
        return copy();
    }

    /**
     * Return an immutable version of this CST curve. This implementation simply returns
     * this BasicCSTCurve instance.
     *
     * @return an immutable version of this curve.
     */
    @Override
    public BasicCSTCurve immutable() {
        return this;
    }

    /**
     * Compares this BasicCSTCurve against the specified object for strict equality (same
     * values and same units).
     *
     * @param obj the object to compare with.
     * @return <code>true</code> if this point is identical to that point;
     *         <code>false</code> otherwise.
     */
    @Override
    public boolean equals(Object obj) {
        if (this == obj)
            return true;
        if ((obj == null) || (obj.getClass() != this.getClass()))
            return false;

        BasicCSTCurve that = (BasicCSTCurve)obj;
        return this.offsetTE_c == that.offsetTE_c
                && this.Cf.equals(that.Cf)
                && this.Sf.equals(that.Sf)
                && this.chord.equals(that.chord)
                && this.origin.equals(that.origin)
                && this.xhat.equals(that.xhat)
                && this.yhat.equals(that.yhat)
                && super.equals(obj);
    }

    /**
     * Returns the hash code for this parameter.
     *
     * @return the hash code value.
     */
    @Override
    public int hashCode() {
        return 31*super.hashCode() 
                + 31*makeVarCode(offsetTE_c) 
                + Objects.hash(Cf, Sf, chord, origin, xhat, yhat);
    }

    private static int makeVarCode(double value) {
        long bits = Double.doubleToLongBits(value);
        int var_code = (int)(bits ^ (bits >>> 32));
        return var_code;
    }

    /**
     * Recycles a <code>BasicCSTCurve</code> instance immediately (on the stack when
     * executing in a StackContext).
     *
     * @param instance The instance to be recycled.
     */
    public static void recycle(BasicCSTCurve instance) {
        FACTORY.recycle(instance);
    }

    /**
     * Holds the default XML representation for this object.
     */
    @SuppressWarnings("FieldNameHidesFieldInSuperclass")
    protected static final XMLFormat<BasicCSTCurve> XML = new XMLFormat<BasicCSTCurve>(BasicCSTCurve.class) {

        @Override
        public BasicCSTCurve newInstance(Class<BasicCSTCurve> cls, XMLFormat.InputElement xml) throws XMLStreamException {
            return FACTORY.object();
        }

        @Override
        public void read(XMLFormat.InputElement xml, BasicCSTCurve obj) throws XMLStreamException {
            double offsetTE_c = xml.getAttribute("offsetTE_c", 0.0);

            CSTCurve.XML.read(xml, obj);     // Call parent read.
            
            NurbsCurve pmap = xml.get("pmap");
            CSTClassFunction Cf = xml.get("Cf");
            CSTShapeFunction Sf = xml.get("Sf");
            Parameter<Length> chord = xml.get("Chord");
            Point origin = xml.get("Origin");
            Vector<Dimensionless> xhat = xml.get("xhat");
            Vector<Dimensionless> yhat = xml.get("yhat");

            obj.pmap = pmap;
            obj.Cf = Cf;
            obj.Sf = Sf;
            obj.offsetTE_c = offsetTE_c;
            obj.chord = chord;
            obj.origin = origin;
            obj.xhat = xhat;
            obj.yhat = yhat;
            obj.calcBoundsMinMax();

        }

        @Override
        public void write(BasicCSTCurve obj, XMLFormat.OutputElement xml) throws XMLStreamException {
            xml.setAttribute("offsetTE_c", obj.offsetTE_c);

            CSTCurve.XML.write(obj, xml);    // Call parent write.
            
            xml.add(obj.pmap, "pmap");
            xml.add(obj.Cf, "Cf");
            xml.add(obj.Sf, "Sf");
            xml.add(obj.chord, "Chord");
            xml.add(obj.origin, "Origin");
            xml.add(obj.xhat, "xhat");
            xml.add(obj.yhat, "yhat");
        }
    };

    ///////////////////////
    // Factory creation. //
    ///////////////////////
    private BasicCSTCurve() { }

    @SuppressWarnings("unchecked")
    private static final ObjectFactory<BasicCSTCurve> FACTORY = new ObjectFactory<BasicCSTCurve>() {
        @Override
        protected BasicCSTCurve create() {
            return new BasicCSTCurve();
        }

        @Override
        protected void cleanup(BasicCSTCurve obj) {
            obj.reset();
            obj.pmap = null;
            obj.Cf = null;
            obj.Sf = null;
            obj.chord = null;
            obj.origin = null;
            obj.xhat = null;
            obj.yhat = null;
            obj._boundsMax = null;
            obj._boundsMin = null;
        }
    };

    @SuppressWarnings("unchecked")
    private static BasicCSTCurve copyOf(BasicCSTCurve original) {
        BasicCSTCurve obj = FACTORY.object();
        obj.pmap = original.pmap.copy();
        obj.Cf = original.Cf.copy();
        obj.Sf = original.Sf.copy();
        obj.offsetTE_c = original.offsetTE_c;
        obj.origin = original.origin.copy();
        obj.chord = original.chord.copy();
        obj.xhat = original.xhat.copy();
        obj.yhat = original.yhat.copy();
        obj._boundsMax = original._boundsMax.copy();
        obj._boundsMin = original._boundsMin.copy();
        original.copyState(obj);
        return obj;
    }

}