class SmoothChain {
  // number of points to be inserted between two adjacent control points
  // NOTE: Points are interpolated evenly in terms of t.
  //       The actual distance is usually not even due to higher order
  //       terms. Only if the coefficients of higher terms (t^2, t^3) are
  //       zero and the coefficient of linear term is non-zero,
  //       interpolated points are evenly separated.
  public var n_interp( __getInterpolateNum, __setInterpolateNum ):Int;
  public var curves( __getCurves, null ):Array< CubicCurve >;
  public var controls( __getControls, null ):Array< WMCPoint >;

  public function new( ?ops:Array< Dynamic > = null,
                       ?ni:Int = 3 ) {
    n_interp = ni;
    controls = new Array< WMCPoint >();
    curves = new Array< CubicCurve >();

    initialize( ops );
  }

  public function clear():Void {
    curves = [];
    controls = [];
  }

  public function initialize( ?ops:Array< Dynamic > = null ):Void {
    if ( ops == null ) return; // nothing to do
    if ( ops.length < 4 ) {    // unable to build chain
      trace( "WMChain::initialize: at least 4 control points are necessary to build a chain" );
      return;
    }
    genFaceDirections( ops );
    genCurves( ops );
    genPoints( ops );
  }

  private function genCurves( pos:Array< Dynamic > ):Void {
    curves = [];
    var num:Int = pos.length;
    // generate first curve
    curves.push( new CubicCurve( null, pos[0].p, pos[1].p, pos[2].p ) );
    for ( i in 0 ... num - 3 ) {
      curves.push( new CubicCurve( pos[i].p, pos[i+1].p, pos[i+2].p, pos[i+3].p ) );
    }
    // last one
    curves.push( new CubicCurve( pos[num-3].p, pos[num-2].p, pos[num-1].p ) );
  }

  private function genPoints( pos:Array< Dynamic > ):Void {
    controls = [];
    var num:Int = curves.length;
    controls.push( new WMCPoint( curves[0].getVal( 0.0 ), pos[0].dir ) );
    for ( i in 0 ... num ) {
      var ni:Int = n_interp;
      var d0:Point3D = pos[i].dir;
      var d1:Point3D = pos[i+1].dir;
      if ( pos[i].n != null && pos[i].n > 0 ) ni = pos[i].n;
      var fn:Float = 1.0 / ni;
      for ( j in 1 ... ni + 1 ) {
        var d0j:Point3D = Point3D.getMultiply( d0, ( ni - j ) / ni );
        var d1j:Point3D = Point3D.getMultiply( d1, j / ni );
        controls.push( new WMCPoint( curves[i].getVal( j * fn ),
                       Point3D.getAdd( d0j, d1j ) ) );
      }
    }
  }

  private function genFaceDirections( pos:Array< Dynamic > ):Void {
    var num:Int = pos.length;
    var prev_face:Point3D = null;
    for ( i in 1 ... num - 1 ) {
      var v0:Point3D = Point3D.getSub( pos[i].p, pos[i-1].p );
      var v1:Point3D = Point3D.getSub( pos[i].p, pos[i+1].p );
      if ( pos[i].dir != null ) {
        pos[i].dir = Point3D.fromString( pos[i].dir );
      } else {
        pos[i].dir = Point3D.getCross( v0, v1 );
      }
      pos[i].dir.normalize();
      if ( prev_face != null ) {
        if ( v0.getAngle( v1 ) > 0.945 * Math.PI ) {
          pos[i].dir = prev_face.clone();
          continue;
        } else {
          if ( pos[i].dir.getAngle( prev_face ) > 0.5 * Math.PI ) pos[i].dir.multiply( -1.0 );
        }
      }
      prev_face = pos[i].dir.clone();
    }
    if ( pos[0].dir == null ) pos[0].dir = pos[1].dir.clone();
    if ( pos[num-1].dir == null ) pos[num-1].dir = pos[num-2].dir.clone();
  }

  public function genCoil( width:Float,
                           quality:Int,
                           from:Int,
                           to:Int ):Tube3D {
    return( new Tube3D( getPartOfControls( from, to ), width, quality ) );
  }

  private function __LagrangeInterpolationDenoms( id0:Int,
                                                  id1:Int,
                                                  id2:Int,
                                                  id3:Int ):Array< Float > {
    var ret:Array< Float > = new Array< Float >();
    ret.push( ( id0 - id1 ) * ( id0 - id2 ) * ( id0 - id3 ) );
    ret.push( ( id1 - id0 ) * ( id1 - id2 ) * ( id1 - id3 ) );
    ret.push( ( id2 - id0 ) * ( id2 - id1 ) * ( id2 - id3 ) );
    ret.push( ( id3 - id0 ) * ( id3 - id1 ) * ( id3 - id2 ) );
    return( ret );
  }

  public function genRibbon( width:Float,
                             thickness:Float,
                             from:Int,
                             to:Int,
                             ?smoothing:Bool = false ):Ribbon3D {
    if ( smoothing ) {
      // count number
      var mynum:Int = 0;
      for ( i in from ... to ) {
        if ( i >= controls.length ) break;
        mynum++;
      }
      // recalculate positions, facedirs
      if ( mynum > 8 ) {
        // use Lagrange interpolation using four control points
        // first, select four control points
        var mynum5:Int = Std.int( mynum / 5 );
        var id0:Int = 0;
        var id1:Int = mynum5;
        var id2:Int = mynum - mynum5 - 1;
        var id3:Int = mynum - 1;
        // calculate denominators
        var denoms:Array< Float > = __LagrangeInterpolationDenoms( id0, id1, id2, id3 );
        var mycs:Array< WMCPoint > = new Array< WMCPoint >();
        // reassign positions, where facedirs are not changed
        var fn:Float = 1.0 / ( mynum - 1 );
        for ( i in 0 ... mynum ) {
          var c0:Point3D = Point3D.getMultiply( controls[from+id0].pos, ( i - id1 ) * ( i - id2 ) * ( i - id3 ) / denoms[0] );
          var c1:Point3D = Point3D.getMultiply( controls[from+id1].pos, ( i - id0 ) * ( i - id2 ) * ( i - id3 ) / denoms[1] );
          var c2:Point3D = Point3D.getMultiply( controls[from+id2].pos, ( i - id0 ) * ( i - id1 ) * ( i - id3 ) / denoms[2] );
          var c3:Point3D = Point3D.getMultiply( controls[from+id3].pos, ( i - id0 ) * ( i - id1 ) * ( i - id2 ) / denoms[3] );
          c0.add( c1 );
          c0.add( c2 );
          c0.add( c3 );
          mycs.push( new WMCPoint( c0, controls[from+i].facedir ) );
        }
        return( genRibbonMain( mycs, width, thickness, 0, mynum ) );
      } else {
        // using qudratic function
        // first, select (almost) evenly separated three points
        var mynum2:Int = Std.int( mynum / 2 );
        var id0:Int = from;
        var id1:Int = from + mynum2;
        var id2:Int = from + mynum - 1;
        // regenerate smooth curve with quadratic function
        var r10 = Point3D.getSub( controls[id1].pos, controls[id0].pos );
        var r20 = Point3D.getSub( controls[id2].pos, controls[id0].pos );
        var b:Point3D = Point3D.getSub( Point3D.getMultiply( r10, 4.0 ), r20 );
        var a:Point3D = Point3D.getSub( r20, b );
        var c:Point3D = controls[id0].pos.clone();
        var mycs:Array< WMCPoint > = new Array< WMCPoint >();
        // reassign positions, where facedirs are not changed
        var fn:Float = 1.0 / ( mynum - 1 );
        for ( i in 0 ... mynum ) {
          var ratio:Float = fn * i;
          var x2:Point3D = Point3D.getMultiply( a, ratio * ratio );
          var x10:Point3D = Point3D.getAdd( Point3D.getMultiply( b, ratio ), c );
          mycs.push( new WMCPoint( Point3D.getAdd( x10, x2 ), controls[from+i].facedir ) );
        }
        return( genRibbonMain( mycs, width, thickness, 0, mynum ) );
      }
    } else {
      return( genRibbonMain( controls, width, thickness, from, to ) );
    }
  }

  private function genRibbonMain( cps:Array< WMCPoint >,
                                  width:Float,
                                  thickness:Float,
                                  from:Int,
                                  to:Int ):Ribbon3D {
    if ( thickness > 0.0 ) {
      var _p0_0:Array< Point3D > = new Array< Point3D >();
      var _p0_1:Array< Point3D > = new Array< Point3D >();
      var _p1_0:Array< Point3D > = new Array< Point3D >();
      var _p1_1:Array< Point3D > = new Array< Point3D >();
      var vec0:Point3D = new Point3D();
      var vec1:Point3D = new Point3D();
      for ( i in from ... to ) {
        if ( i >= cps.length ) break;
        var work:Point3D = Point3D.getMultiply( cps[i].facedir, width );
        if ( i != cps.length - 1 ) {
          vec0 = Point3D.getSub( cps[i+1].pos, cps[i].pos );
        }
        var newvec:Point3D = Point3D.getCross( vec0, cps[i].facedir );
        if ( i != 0 ) {
          if ( newvec.dot( vec1 ) < 0 ) newvec.multiply( -1.0 );
        }
        newvec.multiply( 0.5 * thickness );
        vec1 = newvec;
        var base0:Point3D = Point3D.getAdd( cps[i].pos, work );
        var base1:Point3D = Point3D.getSub( cps[i].pos, work );
        _p0_0.push( Point3D.getAdd( base0, newvec ) );
        _p0_1.push( Point3D.getSub( base0, newvec ) );
        _p1_0.push( Point3D.getAdd( base1, newvec ) );
        _p1_1.push( Point3D.getSub( base1, newvec ) );
      }
      return( new Kishimen( _p0_0, _p0_1, _p1_0, _p1_1 ) );
    } else {
      var _p0:Array< Point3D > = new Array< Point3D >();
      var _p1:Array< Point3D > = new Array< Point3D >();
      for ( i in from ... to ) {
        if ( i >= cps.length ) break;
        var work:Point3D = Point3D.getMultiply( cps[i].facedir, width );
        _p0.push( Point3D.getAdd( cps[i].pos, work ) );
        _p1.push( Point3D.getSub( cps[i].pos, work ) );
      }
      return( new Ribbon3D( _p0, _p1 ) );
    }
  }

  private function getPartOfControls( from:Int,
                                      to:Int ):Array< Point3D > {
    var ret:Array< Point3D > = new Array< Point3D >();
    var rangeF:Int = Std.int( Math.max( from, 0 ) );
    var rangeT:Int = Std.int( Math.min( to, controls.length ) );
    for ( i in rangeF ... rangeT ) ret.push( controls[i].pos );
    return( ret );
  }

  public function getPositions():Array< Point3D > {
    var ret:Array< Point3D > = new Array< Point3D >();
    for ( cont in controls ) ret.push( cont.pos );
    return( ret );
  }

  public function getCurve( i:Int ):CubicCurve {
    return( curves[i] );
  }

  public function __getInterpolateNum():Int { return( n_interp ); }
  public function __getCurves():Array< CubicCurve > { return( curves ); }
  public function __getControls():Array< WMCPoint > { return( controls ); }
  public function __setInterpolateNum( n:Int ):Int {
    n_interp = n;
    return( n );
  }
}
