#
# Computation of flecnodal curves on freeform surfaces.
#
#			Gershon Elber, August 2003
#

EvalOneFlecnodalPts = function( Srf, Pts ): i: Pt:
    return = nil():
    for ( i = 1, 1, sizeof( Pts ),
 	Pt = nth( Pts, i ):
	snoc( seval( Srf, coord( Pt, 1 ), coord( Pt, 2 ) ), return ) );

EvalFlecnodalPts = function( Srfs, Pts ): i:
    if ( sizeof( Pts ) == 1,
	return = EvalOneFlecnodalPts( Srfs, Pts ),
	return = nil():
	for ( i = 1, 1, sizeof( Pts ),
	    snoc( EvalOneFlecnodalPts( Srfs, nth( Pts, i ) ),
	          return ) ) );

EvalOneFlecnodalVecs = function( Srf, Pts, Size ): i: Pt: Vec: DuSrf: DvSrf:
    return = nil():
    DuSrf = sderive( Srf, col ):
    DvSrf = sderive( Srf, row ):
    for ( i = 1, 1, sizeof( Pts ),
 	Pt = nth( Pts, i ):
	Vec = coerce( seval( DuSrf, coord( Pt, 1 ),
			            coord( Pt, 2 ) ), vector_type )
		    * coord( Pt, 3 ) +
	      coerce( seval( DvSrf, coord( Pt, 1 ),
				    coord( Pt, 2 ) ), vector_type )
		    * coord( Pt, 4 ):
        #printf( "mag = %f\\n", list( sqrt( sqr( coord( Pt, 3 ) ) +
	#				   sqr( coord( Pt, 4 ) ) ) ) ):
	Vec = normalize( Vec ) * Size:
	Pt = seval( Srf, coord( Pt, 1 ), coord( Pt, 2 ) ):
	snoc( Pt + coerce( coerce( Pt, point_type ) + Vec, e3 ),
	      return ) );

EvalFlecnodalVecs = function( Srfs, Pts, Size ): i:
    if ( sizeof( Pts ) == 1,
	return = EvalOneFlecnodalVecs( Srfs, Pts, Size ),
	return = nil():
	for ( i = 1, 1, sizeof( Pts ),
	    snoc( EvalOneFlecnodalVecs( Srfs, nth( Pts, i ), Size ),
	          return ) ) );

save_res = resolution;
resolution = 40;

#############################################################################

s1 = sbezier(
    list( list( ctlpt( E3, -1, -1, 0 ),
                ctlpt( E3, -0.635, -0.743, 2.04 ),
                ctlpt( E3, 0.333, -1, -0.8 ),
                ctlpt( E3, 1.83, -1.1, 0.675 ) ),
          list( ctlpt( E3, -1, -0.333, 2.8 ),
                ctlpt( E3, -0.333, -0.333, -4 ),
                ctlpt( E3, 0.333, -0.333, 4 ),
                ctlpt( E3, 0.846, -0.305, -4.72 ) ),
          list( ctlpt( E3, -1, 0.333, 1.92 ),
                ctlpt( E3, -0.333, 0.333, 2.4 ),
                ctlpt( E3, 0.333, 0.333, 3.2 ),
                ctlpt( E3, 1, 0.333, 1.6 ) ),
          list( ctlpt( E3, -1.14, 1.21, -1.49 ),
                ctlpt( E3, -0.333, 1, 2.8 ),
                ctlpt( E3, 0.333, 1, 2 ),
                ctlpt( E3, 0.741, 1.47, -0.622 ) ) ) );
attrib( s1, "transp", 0.5 );

Parab = sparabolc( s1, true );
color( Parab, green );

f = SFlecnodal( s1, 0.05, -1e-6, 1e-1, 3 );

PtE3 = EvalFlecnodalPts( s1, f );
color( PtE3, yellow );
VecE3 = EvalFlecnodalVecs( s1, f, 0.1 );
color( VecE3, red );

interact( list( s1, Parab, PtE3, VecE3 ) );

save( "flecndl1", list( s1, Parab, PtE3, VecE3 ) );

#############################################################################

s1 = sbezier(
    list( list( ctlpt( E3, -1, -1, 0 ),
                ctlpt( E3, -0.667, -1, 0.02 ),
                ctlpt( E2, 0, -1 ),
                ctlpt( E3, 0.667, -1, -0.02 ),
                ctlpt( E2, 1, -1 ) ),
          list( ctlpt( E3, -1, -0.667, 0.07 ),
                ctlpt( E3, -0.494, -0.262, 6 ),
                ctlpt( E2, 0, -0.667 ),
                ctlpt( E3, 0.667, -0.667, 0.005 ),
                ctlpt( E3, 1, -0.667, -0.07 ) ),
          list( ctlpt( E3, -1, 0, 0.118 ),
                ctlpt( E3, -0.667, 0, 0.039 ),
                ctlpt( E3, 0.0671, 0.0801, -0.833 ),
                ctlpt( E3, 0.667, 0, 0.075 ),
                ctlpt( E3, 1, 0, -0.03 ) ),
          list( ctlpt( E3, -1, 0.667, 0.048 ),
                ctlpt( E3, -0.667, 0.667, 0.089 ),
                ctlpt( E3, 0, 0.667, 0.13 ),
                ctlpt( E3, 0.765, 0.682, 5 ),
                ctlpt( E3, 1, 0.667, 0.04 ) ),
          list( ctlpt( E2, -1, 1 ),
                ctlpt( E3, -0.667, 1, 0.07 ),
                ctlpt( E3, 0, 1, 0.12 ),
                ctlpt( E3, 0.667, 1, 0.05 ),
                ctlpt( E2, 1, 1 ) ) ) );
attrib( s1, "transp", 0.5 );

Parab = sparabolc( s1, true );
color( Parab, green );

f = SFlecnodal( s1, 0.05, -1e-6, 1e-1, 3 );

PtE3 = EvalFlecnodalPts( s1, f );
color( PtE3, yellow );
VecE3 = EvalFlecnodalVecs( s1, f, 0.1 );
color( VecE3, red );

interact( list( s1, Parab, PtE3, VecE3 ) );

save( "flecndl2", list( s1, Parab, PtE3, VecE3 ) );

#############################################################################

s1 = sbspline( 4, 4,
    list( list( ctlpt( E3, -1, -1, 0 ),
                ctlpt( E2, -0.6667, -1 ),
                ctlpt( E2, 0, -1 ),
                ctlpt( E2, 0.6667, -1 ),
                ctlpt( E2, 1, -1 ) ),
          list( ctlpt( E3, -1.026, -0.1631, 0.8864 ),
                ctlpt( E3, -0.4817, -0.5642, 0.176 ),
                ctlpt( E3, -0.03526, -0.595, 0.1037 ),
                ctlpt( E3, 0.6046, -0.05832, 1.267 ),
                ctlpt( E3, 1.01, -0.1999, 0.929 ) ),
          list( ctlpt( E3, -1.083, 0.377, 0.8322 ),
                ctlpt( E3, -0.6685, 0.5344, 1.272 ),
                ctlpt( E1, 0 ),
                ctlpt( E3, 0.6425, 0.5309, 1.208 ),
                ctlpt( E3, 1.045, 0.41, 0.8736 ) ),
          list( ctlpt( E3, -1.108, 0.8457, 0.5472 ),
                ctlpt( E3, -0.6899, 0.8424, 0.4805 ),
                ctlpt( E2, 0, 0.6667 ),
                ctlpt( E2, 0.6667, 0.6667 ),
                ctlpt( E2, 1, 0.6667 ) ),
          list( ctlpt( E3, -1.125, 0.8487, -0.1662 ),
                ctlpt( E3, -0.4435, 1.007, -0.05021 ),
                ctlpt( E2, 0, 1 ),
                ctlpt( E2, 0.6667, 1 ),
                ctlpt( E2, 1, 1 ) ) ),
    list( list( kv_open ),
          list( kv_open ) ) );
attrib( s1, "transp", 0.5 );

Parab = sparabolc( s1, true );
color( Parab, green );

f = SFlecnodal( s1, 0.05, -1e-6, 1e-1, 3 );

PtE3 = EvalFlecnodalPts( s1, f );
color( PtE3, yellow );
VecE3 = EvalFlecnodalVecs( s1, f, 0.1 );
color( VecE3, red );

view( list( s1, Parab, PtE3, VecE3 ), 1 );

save( "flecndl3", list( s1, Parab, PtE3, VecE3 ) );

#############################################################################

resolution = save_res;
free( s1 );
free( Parab );
free( f );
free( PtE3 );
free( VecE3 );