string vec3dString( const vec3d v )
{
    return "(" + v.x() + ", " + v.y() + ", " + v.z() + ")";
}

int main()
{
    int failure = 0;

    // The checks below lean on two invariants rather than stored numbers: a
    // surface normal is a unit vector, and a point pushed one unit along that
    // normal projects back to the point it came from at a distance of one.
    double tol = 1e-6;

    string gid = AddGeom( "POD", "" );

    Update();

    int surf_indx = 0;
    double u = 0.12345;
    double w = 0.67890;

    vec3d pnt = CompPnt01( gid, surf_indx, u, w );
    Print( "Point: " + vec3dString( pnt ) );

    vec3d norm = CompNorm01( gid, surf_indx, u, w );
    Print( "Norm: " + vec3dString( norm ) );

    if ( !closeTo( norm.mag(), 1.0, tol ) )
    {
        Print( "FAIL CompNorm01 is not a unit vector, mag " + norm.mag() );
        failure++;
    }

    vec3d tanu = CompTanU01( gid, surf_indx, u, w );
    Print( "Tan U: " + vec3dString( tanu ) );

    vec3d tanw = CompTanW01( gid, surf_indx, u, w );
    Print( "Tan W: " + vec3dString( tanw ) );

    double k1, k2, ka, kg;
    CompCurvature01( gid, surf_indx, u, w, k1, k2, ka, kg );
    Print( "Curvature : k1 " + k1 + " k2 " + k2 + " ka " + ka + " kg " + kg );


    double uout, wout;

    int surf_indx_out;

    // Offset point from surface
    pnt = pnt + norm;

    double d = ProjPnt01( gid, surf_indx, pnt, uout, wout );

    Print( "Dist " + d + " u " + uout + " w " + wout );

    if ( !closeTo( d, 1.0, tol ) || !closeTo( uout, u, tol ) || !closeTo( wout, w, tol ) )
    {
        Print( "FAIL ProjPnt01 gave d " + d + " u " + uout + " w " + wout );
        failure++;
    }

    d =ProjPnt01I( gid, pnt, surf_indx_out, uout, wout );

    Print( "Dist " + d + " u " + uout + " w " + wout + " surf_index " + surf_indx_out );

    if ( !closeTo( d, 1.0, tol ) || !closeTo( uout, u, tol ) || !closeTo( wout, w, tol ) || surf_indx_out != surf_indx )
    {
        Print( "FAIL ProjPnt01I gave d " + d + " u " + uout + " w " + wout + " surf " + surf_indx_out );
        failure++;
    }

    d = ProjPnt01Guess( gid, surf_indx, pnt, u + 0.1, w + 0.1, uout, wout );

    Print( "Dist " + d + " u " + uout + " w " + wout );

    if ( !closeTo( d, 1.0, tol ) || !closeTo( uout, u, tol ) || !closeTo( wout, w, tol ) )
    {
        Print( "FAIL ProjPnt01Guess gave d " + d + " u " + uout + " w " + wout );
        failure++;
    }

    array<double> uvec;
    array<double> wvec;

    int n = 5;

    uvec.resize( n );
    wvec.resize( n );

    for( int i = 0 ; i < n ; i++ )
    {
        uvec[i] = (i+1)*1.0/(n+1);
        wvec[i] = (n-i)*1.0/(n+1);
    }

    array< vec3d > ptvec = CompVecPnt01( gid, 0, uvec, wvec );
    array< vec3d > normvec = CompVecNorm01( gid, 0, uvec, wvec );

    array<double> k1vec;
    array<double> k2vec;
    array<double> kavec;
    array<double> kgvec;

    CompVecCurvature01( gid, 0, uvec, wvec, k1vec, k2vec, kavec, kgvec );

    Print( "Surface query:" );

    for( int i = 0 ; i < n ; i++ )
    {
        Print( "u: " + uvec[i] + " w: " + wvec[i] + " pnt: " + vec3dString( ptvec[i] ) + " norm: " + vec3dString( normvec[i] ) );
        Print( "    k1: " + k1vec[i] + " k2: " + k2vec[i]  + " ka: " + kavec[i]  + " kg: " + kgvec[i] );

        if ( !closeTo( normvec[i].mag(), 1.0, tol ) )
        {
            Print( "FAIL CompVecNorm01 point " + i + " is not a unit vector, mag " + normvec[i].mag() );
            failure++;
        }

        ptvec[i] = ptvec[i] + normvec[i];
    }

    array<double> uoutv;
    array<double> woutv;
    array<double> doutv;

    ProjVecPnt01( gid, 0, ptvec, uoutv, woutv, doutv );

    Print( "Surface projection:" );

    for( int i = 0 ; i < n ; i++ )
    {
        Print( "u: " + uoutv[i] + " w: " + woutv[i] + " d: " + doutv[i] );

        if ( !closeTo( doutv[i], 1.0, tol ) || !closeTo( uoutv[i], uvec[i], tol ) || !closeTo( woutv[i], wvec[i], tol ) )
        {
            Print( "FAIL ProjVecPnt01 point " + i + " did not project back to its origin" );
            failure++;
        }
    }

    array<double> u0v;
    array<double> w0v;

    u0v.resize( n );
    w0v.resize( n );

    for( int i = 0 ; i < n ; i++ )
    {
        u0v[i] = uvec[i] + 0.01234;
        w0v[i] = wvec[i] - 0.05678;
    }

    ProjVecPnt01Guess( gid, 0, ptvec, u0v,  w0v,  uoutv, woutv, doutv );

    Print( "Surface projection with guesses:" );

    for( int i = 0 ; i < n ; i++ )
    {
        Print( "u: " + uoutv[i] + " w: " + woutv[i] + " d: " + doutv[i] );

        if ( !closeTo( doutv[i], 1.0, tol ) || !closeTo( uoutv[i], uvec[i], tol ) || !closeTo( woutv[i], wvec[i], tol ) )
        {
            Print( "FAIL ProjVecPnt01Guess point " + i + " did not project back to its origin" );
            failure++;
        }
    }


    array<double> utess;
    array<double> wtess;
    GetUWTess01( gid, surf_indx, utess, wtess );

    if ( utess.size() == 0 || wtess.size() == 0 )
    {
        Print( "FAIL GetUWTess01 returned empty tessellation arrays" );
        failure++;
    }

    for( uint i = 0 ; i < utess.size() ; i++ )
    {
        Print( "utess: " + utess[i] );
    }

    for( uint j = 0 ; j < wtess.size() ; j++ )
    {
        Print( "wtess: " + wtess[j] );
    }

    //==== Check For API Errors ====//
    while ( GetNumTotalErrors() > 0 )
    {
        ErrorObj err = PopLastError();
        Print( err.GetErrorString() );
        failure++;
    }

    return failure;
}
