#
# Some routines to test surface/curve operators
#

#
# Set states.
#

IProd = iritstate("InterpProd", off);

#
# Set display to on to view some results, off to view nothing.
#
display = on;

#
# Control the surface to polygons subdivison resolution, and isolines gen.
#
save_res = resolution;
save_mat = view_mat;

if ( machine == msdos, resolution = 5, resolution = 20 );

s45 = sin( pi / 4 );

#  
# Circular constructors
#

circ = circle( vector( 0.25, 0.5, 0.5 ), 1.5 );
circ = creparam( circ, 0.0, 1.0 );
arc1 = arc( vector( 0.0, 0.0, 0.0 ),
	    vector( 0.5, 2.0, 0.0 ),
	    vector( 1.0, 0.0, 0.0 ) );
arc2 = arc( vector( 0.0, 0.0, 0.0 ),
	    vector( 0.0, 2.0, 0.5 ),
	    vector( 0.0, 0.0, 1.0 ) );

DLevel = iritState( "DumpLevel", 255 );

dummy = iritState( "DumpLevel", DLevel );

if ( display == on, view_mat = save_mat * sc( 0.5 ):
		    interact( list( view_mat, axes, circ, arc1, arc2 ) ) );

#
# Piecewise linear approximation to a curve using ceval:
#
cb_all = nil();
for ( t = 0.0, 0.05, 1.0,
	cb = ceval( circ, t ):
	snoc( cb, cb_all )
    );
free( t );
cb_crv = creparam( cbspline( 2, cb_all, list( KV_OPEN ) ), 0.0, 2.0 );
color( cb_crv, yellow );
if ( display == on, interact( list( cb_crv, circ, axes ) ) );
free( cb_crv );
free( cb_all );
free( cb );

#
# Power basis curve construction
#

c = cpower( list( ctlpt( E3, 0, 1, 0 ),
		  ctlpt( E3, 1, 0, 0 ),
		  ctlpt( E3, 0, 0, 1 ) ) );
printf( "Power/Bezier curve coercion test = %d\\n",
	list( c == coerce( coerce( c, bezier_type ), power_type ) ) );

if ( display == on, interact( list( axes, coerce( c, bezier_type ) ) ) );

s = spower( list( list( ctlpt( E3, 1, 0, 1 ),
		        ctlpt( E3, 0, 1, 1 ) ),
		  list( ctlpt( E3, 0, 0, 1 ),
		        ctlpt( E3, 0, 0, 1 ) ) ) );
printf( "Power/Bezier surface coercion test = %d\\n",
	list( s == coerce( coerce( s, bezier_type ), power_type ) ) );
if ( display == on, interact( list( axes, coerce( s, bezier_type ) ) ) );

s = spower( list( list( ctlpt( E3, 0, 1, 0 ),
		        ctlpt( E3, 1, 0, 0 ),
		        ctlpt( E3, 0, 0, 1 ) ),
		  list( ctlpt( E3, 1, 1, 0 ),
		        ctlpt( E3, 1, 0, 1 ),
		        ctlpt( E3, 0, 1, 1 ) ),
		  list( ctlpt( E3, 0, -2, 0 ),
		        ctlpt( E3, 2, 0, 0 ),
		        ctlpt( E3, 0, 0, 2 ) ),
		  list( ctlpt( E3, -3, 0, 0 ),
		        ctlpt( E3, 0, -1, 0 ),
		        ctlpt( E3, 0, 0, -3 ) ) ) );
printf( "Power/Bezier surface coercion test = %d\\n",
	list( s == coerce( coerce( s, bezier_type ), power_type ) ) );
if ( display == on, view_mat = save_mat * sc( 0.5 ):
		    interact( list( view_mat, axes,
				    coerce( s, bezier_type ) ) ) );
#
# Reparametrization examples.
# 

c1 = cbspline( 3,
               list( ctlpt( E2, 1.7, 0.0 ),
                     ctlpt( E2, 0.7, 0.7 ),
                     ctlpt( E2, 1.7, 0.3 ),
                     ctlpt( E2, 1.5, 0.8 ),
                     ctlpt( E2, 1.6, 1.0 ) ),
               list( KV_OPEN ) );
color( c1, red );

c2 = creparam( c1, param_chord, param_chord );
color( c2, yellow );

c3 = creparam( c1, param_centrip, param_centrip );
color( c3, green );

if ( display == on, view_mat = save_mat * sc( 0.5 ):
		    interact( list( view_mat, axes, c1, c2, c3 ) ) );

s1 = coerce( s, bspline_type );
s1 = srefine( srefine( s1, row, false, list( 0.3, 0.6 ) ),	
	      col, false, list( 0.25, 0.5, 0.75 ) );
color( s1, red );

s2 = sreparam( sreparam( s1, col, param_chord, param_chord ),
	       row, param_chord, param_chord );
color( s2, yellow );

s3 = sreparam( sreparam( s1, col, param_centrip, param_centrip ),
	       row, param_centrip, param_centrip );
color( s3, green );

if ( display == on, view_mat = save_mat * sc( 0.5 ):
		    interact( list( view_mat, axes, s1, s2, s3 ) ) );

save( "reparam", list( c1, c2, c3, s1, s2, s3 ) );

free( c1 );
free( c2 );
free( c3 );
free( s1 );
free( s2 );
free( s3 );

#
# Monkey saddle (u, v, u^3 - 3v^2*u):
#
s = spower( list( list( ctlpt( E3, 0, 0, 0 ),
		        ctlpt( E3, 1, 0, 0 ),
		        ctlpt( E3, 0, 0, 0 ),
		        ctlpt( E3, 0, 0, 1 ) ),
		  list( ctlpt( E3, 0, 1, 0 ),
		        ctlpt( E3, 0, 0, 0 ),
		        ctlpt( E3, 0, 0, 0 ),
		        ctlpt( E3, 0, 0, 0 ) ),
		  list( ctlpt( E3, 0, 0, 0 ),
		        ctlpt( E3, 0, 0, -3 ),
		        ctlpt( E3, 0, 0, 0 ),
		        ctlpt( E3, 0, 0, 0 ) ) ) );
Monkey = sregion( sregion( coerce( s, bezier_type ),
			   row, -1, 1 ),
		  col, -1, 1 );
if ( display == on, view_mat = save_mat * sc( 0.3 ):
		    interact( list( view_mat, axes, Monkey ) ) );

free( c );
free( s );
free( Monkey );

#
# Ruled surface constructor examples.
#
arc3 = arc( vector( 0.0,  0.0, 1.0 ),
	    vector( 0.5, -0.2, 1.0 ),
	    vector( 1.0,  0.0, 1.0 ) );
ruled = ruledsrf( arc3,
		  ctlpt( E2, 0.0, 0.0 ) + ctlpt( E2, 1.0, 0.0 ) );

dummy = iritState( "DumpLevel", 255 );

dummy = iritState( "DumpLevel", DLevel );
free( DLevel );

if ( display == on, view_mat = save_mat * ty( -0.3 ):
		    interact( list( view_mat, axes, ruled ) ) );
free( ruled );

circ = circle( vector( 0.0, 0.0, 0.0 ), 0.25 );
cyl = ruledsrf( circ, circ*trans( vector( 0.0, 0.0, 1.0 ) ) );
if ( display == on, interact( list( axes, cyl ) ) );
free( cyl );

skewcyl = ruledsrf( circ, circ * trans( vector( 0.2, 0.0, 1.0 ) ) );
skewcylMesh = ffmesh( skewcyl );
color( skewcylMesh, red );
if ( display == on, interact( list( axes, skewcyl, skewcylMesh ) ) );
free( skewcylMesh );
free( skewcyl );

skew2cyl = ruledsrf( circ * rotx( 20 ),
		     circ * rotx( -20 ) * trans( vector( 0.0, 0.0, 1.0 ) ) );
skewcyl2Mesh = ffmesh( skew2cyl );
color( skewcyl2Mesh, red );
if ( display == on, interact( list( axes, skew2cyl, skewcyl2Mesh ) ) );
free( skew2cyl );
free( skewcyl2Mesh );
free( arc1 );
free( arc2 );
free( arc3 );
free( circ );

#
# Degeneracy prevention in ruled surfaces.
#

c1 = cbspline( 3,
               list( ctlpt( E3, 1.7, 0.0, 0 ),
                     ctlpt( E3, 0.7, 0.7, 0 ),
                     ctlpt( E3, 1.7, 0.3, 0 ),
                     ctlpt( E3, 1.5, 0.8, 0 ),
                     ctlpt( E3, 1.6, 1.0, 0 ) ),
               list( KV_OPEN ) );
color( c1, red );
c2 = cbspline( 3,
               list( ctlpt( E3,  0.7, 0.0, 0 ),
                     ctlpt( E3, -0.7, 0.2, 0 ),
                     ctlpt( E3,  0.7, 0.5, 0 ),
                     ctlpt( E3, -0.7, 0.7, 0 ),
                     ctlpt( E3,  0.7, 1.0, 0 ) ) ,
               list( KV_OPEN ) );
color( c2, red );
s1 = ruledsrf( c1, c2 );
color( s1, yellow );
if ( display == on, view_mat = trans( vector( -0.7, -0.5, 0.0 ) ):
		    interact( list( view_mat, c1, c2, s1 ) ) );

c2a = ffmatch( c2, c1, 50, 100, 2, false, 1 );
s2 = ruledsrf( c2, c2a );
color( s2, yellow );
if ( display == on, interact( list( c1, c2, s2 ) ):
		    view_mat = save_mat:
		    viewobj( view_mat ) );

free( s1 );
free( s2 );
free( c1 );
free( c2 );
free( c2a );

#
# Curve constructors.
#
crv1 = ctlpt( E3, 0.5, 0.0, 1.0 ) +
       ctlpt( P3, 1.0, 0.5, 0.0, 1.2 ) +  # Note ctlpts can be different type.
       ctlpt( E3, 1.0, 0.0, 1.2 );
crv2 = crv1 +
       arc( vector( 1.0,  0.0, 0.75 ),
	    vector( 0.75, 0.0, 0.7 ),
	    vector( 0.5,  0.0, 0.85 ) ) +
       arc( vector( 0.5,  0.0, 0.75 ),
	    vector( 0.75, 0.0, 0.8 ),
	    vector( 1.0,  0.0, 0.65 ) );
crv3 = crv2 +
       crv2 * trans( vector( -0.5, 0.15, -0.5 ) ) +
       crv2 * trans( vector( -1.0, 0.3,  -1.0 ) );

if ( display == on, interact( list( axes, crv3 ) ) );

free( crv1 );
free( crv2 );
free( crv3 );

cross = arc( vector( 0.2, 0.0, 0.0 ),
	     vector( 0.2, 0.2, 0.0 ),
	     vector( 0.0, 0.2, 0.0 ) ) +
	arc( vector( 0.0, 0.4, 0.0 ),
	     vector( 0.1, 0.4, 0.0 ),
	     vector( 0.1, 0.5, 0.0 ) ) +
	arc( vector( 0.8, 0.5, 0.0 ),
	     vector( 0.8, 0.3, 0.0 ),
	     vector( 1.0, 0.3, 0.0 ) ) +
	arc( vector( 1.0, 0.1, 0.0 ),
	     vector( 0.9, 0.1, 0.0 ),
	     vector( 0.9, 0.0, 0.0 ) ) +
	ctlpt( E2, 0.2, 0.0 );
if ( display == on, interact( list( axes, cross ) ) );

#
# Curves and surfaces convertions and compatibility.
#
crv1 = cbezier( list( ctlpt( E2, -1, 0 ),
		      ctlpt( E3,  0, 2, 0 ),
		      ctlpt( E2,  1, 0 ) ) );
crv2 = cbspline( 4,
		 list( ctlpt( P3,  1.0,  0.0, 0.1,  1.0 ),
		       ctlpt( P3,  s45, -s45, 1.0,  s45 ),
		       ctlpt( P2,  1.0,  1.0, 0.1 ),
		       ctlpt( P3,  s45, -s45, 1.0, -s45 ),
		       ctlpt( P3,  1.0,  0.0, 0.1, -1.0 ) ),
		 list( KV_OPEN ) );
color( crv1, yellow );
color( crv2, yellow );

crv1bsp = bzr2bsp( crv1 );
color( crv1bsp, green );
crv2bzr = bsp2bzr( crv2 );
color( crv2bzr, green );
if ( display == on, interact( list( crv1bsp, crv1, crv2bzr, crv2 ) ) );
free( crv1bsp );
free( crv2bzr );

if ( display == on, interact( list( crv1, crv2 ) ) );
ffcompat( crv1, crv2 );
if ( display == on, interact( list( crv1, crv2 ) ) );
free( crv1 );
free( crv2 );

srf1 = sbezier( list( list( ctlpt( E2, 0.0, 0.0 ),
		            ctlpt( E2, 0.0, 1.0 ),
		            ctlpt( E3, 0.0, 2.0, 0.5 ) ),
	              list( ctlpt( E2, 1.0, 0.0 ),
		            ctlpt( E2, 1.0, 1.0 ),
		            ctlpt( E3, 1.0, 2.0, -0.5 ) ) ) );
srf2 = sbspline( 2, 4,
		 list( list( ctlpt( E3, 0.0, 0.0, 1.0 ),
		             ctlpt( E2, 0.0, 1.0 ),
		             ctlpt( E3, 0.0, 2.0, 1.0 ) ),
        	       list( ctlpt( E2, 1.0, 0.0 ),
		             ctlpt( E3, 1.0, 1.0, 2.0 ),
		             ctlpt( E2, 1.0, 2.0 ) ),
        	       list( ctlpt( E3, 2.0, 0.0, 2.0 ),
	        	     ctlpt( E2, 2.0, 1.0 ),
		             ctlpt( E3, 2.0, 2.0, 2.0 ) ),
        	       list( ctlpt( E2, 3.0, 0.0 ),
	        	     ctlpt( E3, 3.0, 1.0, 2.0 ),
		             ctlpt( E2, 3.0, 2.0 ) ),
        	       list( ctlpt( E3, 4.0, 0.0, 1.0 ),
	        	     ctlpt( E2, 4.0, 1.0 ),
		             ctlpt( E3, 4.0, 2.0, 1.0 ) ) ),
		 list( list( KV_OPEN ),
		       list( KV_OPEN ) ) );

color( srf1, yellow );
color( srf2, yellow );

srf1bsp = bzr2bsp( srf1 );
color( srf1bsp, green );
srf2bzr = bsp2bzr( srf2 );
color( srf2bzr, green );

if ( display == on, view_mat = save_mat * sc( 0.45 ) 
					* trans( vector( 0.35, 0.5, -0.3 ) ):
		    interact( list( view_mat,
				    srf1bsp, srf1, srf2bzr, srf2 ) ) );
free( srf1bsp );
free( srf2bzr );

if ( display == on, interact( list( srf1, srf2 ) ) );
ffcompat( srf1, srf2 );
if ( display == on, interact( list( srf1, srf2 ) ) );
free( srf1 );
free( srf2 );

#
# Create floating end condition examples (although barely used!).
#
viewstate( "DSrfMesh", on );
crv1f = cbspline( 4, 
		  list( ctlpt( E2, 0.0, 0.0 ),
		 	ctlpt( E2, 1.0, 1.0 ),
		        ctlpt( E2, 0.0, 2.0 ),
	                ctlpt( E2, 1.0, 0.0 ) ),
		  list( KV_FLOAT ) );
color( crv1f, yellow );

if ( display == on, view_mat = trans( vector( -0.5, -0.8, 0.0 ) ) * sc( 0.8 ):
		    view( list( view_mat, crv1f ), on ):
		    viewstate( "DSrfMesh", 0 ):
		    pause():
		    viewstate( "DSrfMesh", 1 ):
		    pause() );

srf1f = sbspline( 2, 4,
		  list( list( ctlpt( E3, 0.0, 0.0, 1.0 ),
		              ctlpt( E3, 0.0, 1.0, 0.0 ),
		              ctlpt( E3, 0.0, 2.0, 1.0 ) ),
        	        list( ctlpt( E3, 1.0, 0.0, 0.0 ),
		              ctlpt( E3, 1.0, 1.0, 2.0 ),
		              ctlpt( E3, 1.0, 2.0, 0.0 ) ),
        	        list( ctlpt( E3, 2.0, 0.0, 2.0 ),
	        	      ctlpt( E3, 2.0, 1.0, 0.0 ),
		              ctlpt( E3, 2.0, 2.0, 2.0 ) ),
        	        list( ctlpt( E3, 3.0, 0.0, 0.0 ),
	        	      ctlpt( E3, 3.0, 1.0, 2.0 ),
		              ctlpt( E3, 3.0, 2.0, 0.0 ) ),
        	        list( ctlpt( E3, 4.0, 0.0, 1.0 ),
	        	      ctlpt( E3, 4.0, 1.0, 0.0 ),
		              ctlpt( E3, 4.0, 2.0, 1.0 ) ) ),
		  list( list( KV_FLOAT ),
		        list( KV_FLOAT ) ) );
color( srf1f, yellow );
if ( display == on, view_mat = save_mat * trans( vector( 0.6, 0.6, 0.0 ) )
					* sc( 0.4 ):
		    view( list( view_mat, srf1f ), on ):
		    viewstate("DSrfMesh", 1 ):
		    pause():
		    viewstate("DSrfMesh", 0 ):
		    pause() );

free( crv1f );
free( srf1f );
viewstate( "DSrfMesh", 0 );

#
# Direct control points manipulation.
#
cb = cbezier( list( ctlpt( E3, -1, 0, 0 ),
		    ctlpt( E3,  0, 0, 0 ),
		    ctlpt( E3,  1, 0, 0 ) ) );

cb_all = list( axes );
for ( z = -0.9, 0.3, 0.9,
	cb1 = ceditpt( cb, ctlpt( E3, 0, 0, z ), 1 ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
if ( display == on, view_mat = save_mat:
		    interact( list( view_mat, cb_all ) ) );
free( cb_all );
free( cb );
free( cb1 );

sb = ruledSrf( cbezier( list( ctlpt( E3, -0.5, -0.5, 0.0 ),
			      ctlpt( E3,  0.5, -0.5, 0.0 ) ) ),
	       cbezier( list( ctlpt( E3, -0.5,  0.5, 0.0 ),
			      ctlpt( E3,  0.5,  0.5, 0.0 ) ) ) );
sb = sraise( sraise( sb, row, 3 ),
	     col, 3 );
sb = srefine( srefine( sb, row, false, list( 0.333, 0.667 ) ),
	      col, false, list( 0.333, 0.667 ) );
sb_all = list( axes );
for ( z = -0.9, 0.6, 0.9,
	sb1 = seditpt( sb, ctlpt( E3, 0, 0, z ), 2, 2 ):
	color( sb1, green ):
	snoc( sb1, sb_all )
    );
free( z );
if ( display == on, interact( sb_all ) );
free( sb_all );
free( sb );
free( sb1 );

#
# Curve order raise/reduce:
#
cb = cbezier( list( ctlpt( E2, -1, 0 ),
		    ctlpt( E2,  0, 2 ),
		    ctlpt( E2,  1, 0 ) ) );
cb_all = list( axes, cb );
for ( o = 4, 1, 8,
	cb1 = craise( cb, o ) * tz( 0.1 * (o - 3) ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
for ( o = 7, -1, 3,
	cb1 = creduce( cb1, o ) * tz( 0.1 ):
	color( cb1, yellow ):
	snoc( cb1, cb_all )
    );
free( o );
if ( display == on, view_mat = save_mat * sc( 0.7 ) * tx( -0.3 ):
		    view( list( view_mat, cb_all ), on ):
		    viewstate("DSrfMesh", 1 ):
		    pause():
		    viewstate("DSrfMesh", 0 ):
		    pause() );
free( cb_all );
free( cb );
free( cb1 );
viewstate("DSrfMesh", 0 ):

#
# comaprison tests
#

x = cregion( circle( vector( 1, 2, 0 ), 1 ) * sy( 0.5 ), 0, 2.7 );

y1 = x * rz( 45 ) * sc( 0.5 ) * tx( 0.5 ) * ty( -1.3333 );
y2 = x * sc( 2.5 ) * rz( -45 ) * tx( -1/3 ) * ty( pi );
y3 = x * sc( -0.5 );

CrvRigidSim = list( FfRigidSim( x, y1, 1e-6 ),
                    FfRigidSim( x, y2, 1e-6 ),
		    FfRigidSim( x, y3, 1e-6 ) );

x = surfRev( x * rx( 90 ) );

y1 = x * rz( -11.5 ) * sc( 11.5 ) * tx( -3.5 ) * ty( 1.3333 );
y2 = x * sc( 0.35 ) * rz( -45 ) * tx( 1/3 ) * ty( sqrt( 2 ) );
y3 = x * rz( 1 );

SrfRigidSim = list( FfRigidSim( x, y1, 1e-6 ),
                    FfRigidSim( x, y2, 1e-6 ),
		    FfRigidSim( x, y3, 1e-6 ) );

save("ffrgdcmp", list( CrvRigidSim, SrfRigidSim ) );

free( CrvRigidSim );
free( SrfRigidSim );
free( x );
free( y1 );
free( y2 );
free( y3 );

#
# Extrusion examples.
#
cbzr = cbezier( list( ctlpt( E2, 0.0, 0.0 ),
		      ctlpt( E2, 1.0, 0.0 ),
		      ctlpt( E2, 1.0, 1.0 ) ) );
color( cbzr, white );

s = extrude( cbzr, vector( 0.0, 0.0, 1.0 ), 0 );
if ( display == on, view_mat = save_mat:
		    interact( list( view_mat, axes, s, cbzr ) ) );
free( cbzr );
free( s );

s = extrude( cross, vector( 0.0, 0.0, 1.0 ), 0 );
if ( display == on, interact( list( axes, s ) ) );

s = extrude( cross, vector( 0.1, 0.2, 1.0 ), 0 );
if ( display == on, interact( list( axes, s ) ) );

resolution = 10;
ps = gpolygon( s, off );
if ( display == on, interact( list( axes, ps ) ) );
free( ps );
free( s );

if ( machine == msdos, resolution = 5, resolution = 20 );

#
# Srf of revolution examples
#

cbzr = cbezier( list( ctlpt( E3, 0.0, 0.0, 0.0 ),
		      ctlpt( E3, 1.0, 0.0, 0.0 ),
		      ctlpt( E3, 1.0, 0.0, 1.0 ) ) );
color( cbzr, white );

sb = surfrev( cbzr );
if ( display == on, view_mat = save_mat * ty( -0.5 ):
		    interact( list( view_mat, sb, cbzr ) ) );
free( sb );
free( cbzr );

halfcirc = cbspline( 3,
		     list( ctlpt( P3, 1.0,  0.0, 0.0,  1.0 ),
		           ctlpt( P3, s45, -s45, 0.0,  s45 ),
		           ctlpt( P3, 1.0, -1.0, 0.0,  0.0 ),
		           ctlpt( P3, s45, -s45, 0.0, -s45 ),
			   ctlpt( P3, 1.0,  0.0, 0.0, -1.0 ) ),
		     list( 0, 0, 0, 1, 1, 2, 2, 2 ) );
color( halfcirc, white );

sp = surfrev( halfcirc );

if ( display == on, view_mat = save_mat:
		    interact( list( view_mat, sp, halfcirc ) ) );
free( halfcirc );
free( sp );

gcross = cbspline( 3,
		   list( ctlpt( E3, 0.3, 0.0, 0.0 ),
		         ctlpt( E3, 0.3, 0.0, 0.05 ),
		         ctlpt( E3, 0.1, 0.0, 0.05 ),
		         ctlpt( E3, 0.1, 0.0, 0.4 ),
		         ctlpt( E3, 0.5, 0.0, 0.4 ),
		         ctlpt( E3, 0.6, 0.0, 0.8 ) ),
		   list( 0, 0, 0, 1, 2, 3, 4, 4, 4 ) );
color( gcross, white );
glass = surfprev( gcross );
color( glass, red );
if ( display == on, view_mat = save_mat * ty( -0.3 ):
		    interact( list( view_mat, axes, gcross, glass ) ) );
free( gcross );

#
# Ray surface intersection.
#

RayOrigin = point( 2, 0.1, 0.3 );
RayDir = vector( -4, 0, 0 );

RayLine = coerce( RayOrigin, E3 ) + coerce( RayOrigin + RayDir, E3 );
color( RayLine, magenta );
adwidth( RayLine, 2 );

resolution = 5;
InterPt = srinter( glass, RayOrigin, RayDir );
InterPtE3 = seval( glass, coord( InterPt, 0 ), coord( InterPt, 1 ) );
color( InterPtE3, cyan );
adwidth( InterPtE3, 3 );

if ( display == on, interact( list( InterPtE3, RayLine, glass, axes ) ) );

resolution = 80;
InterPt = srinter( glass, RayOrigin, RayDir );
InterPtE3 = seval( glass, coord( InterPt, 0 ), coord( InterPt, 1 ) );
color( InterPtE3, cyan );
adwidth( InterPtE3, 3 );

if ( display == on, interact( list( InterPtE3, RayLine, glass, axes ) ) );

free( RayOrigin );
free( RayDir );
free( RayLine );
free( InterPt );
free( InterPtE3 );

resolution = save_res;

#
# Contouring example.
#

resolution = 50;
cntrs = nil();
for ( i = 0, 1, 5,
   cntrs = cntrs + list( contour( glass,
				  plane( 1, 1, 1, ( i - 4 ) / 3  ) ) ) );
color( cntrs, cyan );
adwidth( cntrs, 4 );

if ( display == on, interact( list( axes, glass, cntrs ) ) );
free( cntrs );

#
# Silhouette extraction using contouring of normal field surface.
#

resolution = 30;
nglass = snrmlsrf( glass ) * vector( 1, 1, 1 );

sils = contour( nglass, plane( 1, 0, 0, 0 ), glass );
color( sils, cyan );
adwidth( sils, 3 );

if ( display == on, interact( list( axes, glass, sils ) ) );

# Silhouettes using the silhouette function:
sils = silhouette( glass, vector( 1, -2, 1 ), true );
color( sils, cyan );
adwidth( sils, 3 );

if ( display == on,
     interact( list( axes, glass, sils, vector( 1, -2, 1 ) ) ) );

for ( i = 10, 10, 180,
    sils = contour( nglass, plane( 1, 0, 0, 0 ), glass,
		   list( vector( -1, 1, 0 ), i ) ):
    if ( display == on, view( list( sils, glass ), 1 ) ):
    miliSleep( 50 ) );

free( nglass );
free( sils );

# Polar silhouettes using the polar silhouette function:
sils = polarsil( glass * tx( 0.15 ) * ty( 0.15 ), vector( 0, 0, 1 ), true );
color( sils, cyan );
adwidth( sils, 3 );

if ( display == on,
     interact( list( axes, glass * tx( 0.15 ) * ty( 0.15 ), sils ) ) );

# Isoclines:
resolution = 15;
sils = silhouette( glass, vector( 1, -2, 1 ), true );
color( sils, yellow );
isocs = list( sils,
	      isocline( glass, vector( 1, -2, 1 ), 75, true, 0 ),
	      isocline( glass, vector( 1, -2, 1 ), 60, true, 0 ) );
color( isocs, white );
adwidth( isocs, 3 );

if ( display == on, interact( list( axes, glass, isocs,
				    vector( 1, -2, 1 ) ) ) );

free( sils );
free( isocs );

# Silhouette, given a view point:
silh = function( srf, viewPt ): Zset:
    Zset = SymbDProd( srf * scale( vector( -1, -1, -1 ) ) * trans( viewPt ),
		      sNrmlSrf( srf ) ):
    return = contour( Zset, plane( 1, 0, 0, 0 ), srf );

ViewPoint = point( 0.5, 0.5, 0.5 );
color( ViewPoint, green );
adwidth( ViewPoint, 3 );

sils = silh( glass, coerce( ViewPoint, vector_type ) );
color( sils, cyan );
adwidth( sils, 3 );

if ( display == on, interact( list( axes, glass, sils, ViewPoint ) ) );
free( sils );
free( ViewPoint );

free( glass );

#
# Sweep examples.
#
ccross = -cross * trans( vector( -0.5, -0.25, 0.0 ) ) * sc( 0.2 );

sweep_axis = crefine( cbspline( 3,
			        list( ctlpt( E3,  0.0,  1.0, 1.0 ),
				      ctlpt( E3,  1.0,  1.0, 1.0 ),
				      ctlpt( E3,  1.0, -1.0, 0.0 ),
				      ctlpt( E3, -1.0, -1.0, 0.0 ),
				      ctlpt( E3, -1.0,  1.0, 1.0 ),
				      ctlpt( E3,  0.0,  1.0, 1.0 ) ),
			        list( KV_OPEN ) ),
		      false,
		      list( 0.125, 0.375, 0.625, 0.875 ) );

arc1 = arc( vector( 1.0, 0.0, 0.0 ),
	    vector( 0.0, 0.0, 0.0 ),
	    vector( 0.0, 1.0, 0.0 ) );
s = sweepsrf( ccross, arc1, off );
if ( display == on, view_mat = save_mat * ty( 0.2 ):
		    interact( list( view_mat, axes, s ) ) );

arc1 = arc( vector( -1.0, 0.0, 0.0 ),
	    vector( 0.0, 0.0, 0.1 ),
	    vector( 1.0, 0.0, 0.0 ) );
arc1 = crefine( arc1, false, list( 0.25, 0.5, 0.75 ) );
scalecrv = cbezier( list( ctlpt( E2, 0.0, 0.1 ),
	       	          ctlpt( E2, 1.0, 0.5 ),
	       	          ctlpt( E2, 2.0, 0.1 ) ) );
s = swpsclsrf( ccross * sc( 6 ), arc1, scalecrv, off, 2 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( list( ccross * sc( 6 ),
		     ccross * sx( 6 ),
		     ccross * sy( 6 ),
		     ccross * sx( 6 ),
		     ccross * sc( 6 ) ),
	       arc1, scalecrv, off, 3 );
if ( display == on, interact( list( axes, s ) ) );

circ = -pcircle( vector( 0.0, 0.0, 0.0 ), 1.0 );

s = sweepsrf( circ * scale( vector( 0.25, 0.25, 1.0 ) ), circ, off );

if ( display == on, view_mat = save_mat * sc( 0.75 ):
		    interact( list( view_mat, axes, s ) ) );

resolution = 10;
ps = gpolygon( s, off );
if ( display == on, interact( list( axes, ps ) ) );
ps = gpolygon( s, on );
if ( display == on, interact( list( axes, ps ) ) );
free( ps );

if ( machine == msdos, resolution = 5, resolution = 20 );

s = swpsclsrf( circ, circ, scalecrv, off, 3 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( ccross, sweep_axis, 1.0, off, 1 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( ccross, sweep_axis, 1.0, vector( 0.0, 0.0, 1.0 ), 1 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( ccross, circle( vector( 0.0, 0.0, 0.0 ), 1.0 ), 1.0,
	       circle( vector( 0.0, 0.0, 1.0 ), -1.0 ), 3 );
if ( display == on, interact( list( axes, s ) ) );

free( sweep_axis );
free( ccross );
free( circ );
free( arc1 );
view(axes,1);

Cross = arc( vector( -0.11, -0.1,  0.0 ),
             vector( -0.1,  -0.1,  0.0 ),
             vector( -0.1,  -0.11, 0.0 ) ) +
        arc( vector(  0.1,  -0.11, 0.0 ),
             vector(  0.1,  -0.1,  0.0 ),
             vector(  0.11, -0.1,  0.0 ) ) +
        arc( vector(  0.11,  0.1,  0.0 ),
             vector(  0.1,   0.1,  0.0 ),
             vector(  0.1,   0.11, 0.0 ) ) +
        arc( vector( -0.1,   0.11, 0.0 ),
             vector( -0.1,   0.1,  0.0 ),
             vector( -0.11,  0.1,  0.0 ) ) +
        ctlpt( E2, -0.11, -0.1 );
scaleCrv = cbspline( 3,
                     list( ctlpt( E2, 0.05, 1.0 ),
                           ctlpt( E2, 0.1,  0.0 ),
                           ctlpt( E2, 0.2,  2.0 ),
                           ctlpt( E2, 0.3,  0.0 ),
                           ctlpt( E2, 0.4,  2.0 ),
                           ctlpt( E2, 0.5,  0.0 ),
                           ctlpt( E2, 0.6,  2.0 ),
                           ctlpt( E2, 0.7,  0.0 ),
                           ctlpt( E2, 0.8,  2.0 ),
                           ctlpt( E2, 0.85, 1.0 ) ),
                     list( KV_OPEN ) );
Axis = -pcircle( vector( 0, 0, 0 ), 1 );
Frame = -pcircle( vector( 0, 0, 0 ), 1 ) *
        rotx( 90 ) * trans( vector( 1.5, 0.0, 0.0 ) );

s = swpsclsrf( Cross, Axis, scaleCrv, off, 2 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( Cross, Axis, 1.0, off, 2 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( Cross, Axis, 1.0, Frame, 7 );
if ( display == on, interact( list( axes, s ) ) );

#
# Periodic curves and surfaces.
#
cPer = cbspline( 4,
		 list( ctlpt( E2, -1, -1 ),
		       ctlpt( E2, -1,  1 ),
		       ctlpt( E2,  1,  1 ),
		       ctlpt( E2,  1, -1 ) ),
		 list( KV_PERIODIC ) );
color( cPer, red );
cPer2 = crefine( cPer, false, list( 0.0625, 0.125, 0.375 ) );
color( cPer2, green );
cPer3 = crefine( cPer2, false, list( 0.875, 0.875, 0.875 ) );
color( cPer3, yellow );
if ( display == on, interact( list( cPer, cPer2, cPer3 ) ) );

s = extrude( cPer, vector( 0.0, 0.0, 1.0 ), 0 );
color( s, red );
s2 = srefine( srefine( srefine( s, col, false, list( 0.125, 0.125, 0.125 ) ),
	               row, false, list( 0.25, 0.5 ) ),
	      col, false, list( 0.625, 0.75, 0.875 ) );
color( s2, yellow );
if ( display == on, view_mat = save_mat * sc( 0.75 ) * ty( -0.2 ):
		    interact( list( view_mat, axes, s, s2, cper ) ) );

s = swpsclsrf( cPer, Axis, 0.1, off, 0 );
color( s, red );
s2 = srefine( srefine( srefine( s, col, false, list( 0.125, 0.125, 0.125 ) ),
	               row, false, list( 0.25, 0.5 ) ),
	      col, false, list( 0.625, 0.75, 0.875 ) );
color( s2, yellow );
if ( display == on, interact( list( axes, s, s2, cper ) ) );

s = swpsclsrf( cPer * scale( vector( 0.2, 0.2, 0.2 ) ),
	       Axis, ScaleCrv, off, 2 );
color( s, red );
s2 = srefine( srefine( srefine( s, col, false, list( 0.125, 0.125, 0.125 ) ),
	               row, false, list( 0.25, 0.5 ) ),
	      col, false, list( 0.625, 0.75, 0.875 ) );
color( s2, yellow );
if ( display == on, interact( list( axes, s, s2, cper ) ) );

free( Axis );
free( ScaleCrv );
free( cPer );
free( cPer2 );
free( cPer3 );
free( s );
free( s2 );

#
# Sweep with varying cross section
#

Cross = -pcircle( vector( 0.0, 0.0, 0.0 ), 1.0 ) * rz( -90 );

Cross2 = cbspline( 4,
	list( ctlpt( E3, 1.0,   0.0,  0 ),
	      ctlpt( E3, 1.0,   0.2,  0 ),
	      ctlpt( E3, 0.6,   0.12, 0 ),
	      ctlpt( E3, 0.5,   0.2,  0 ),
	      ctlpt( E3, 0.8,   0.6,  0 ),
	      ctlpt( E3, 0.6,   0.8,  0 ),
	      ctlpt( E3, 0.2,   0.5,  0 ),
	      ctlpt( E3, 0.12,  0.6,  0 ),
	      ctlpt( E3, 0.2,   1.0,  0 ),
	      ctlpt( E3, 0.0,   1.0,  0 ) ),
	list( kv_open ) );
Cross2 = Cross2 + Cross2 * rz( 90 ) + Cross2 * rz( 180 ) + Cross2 * rz( 270 );

Axis = cbspline( 4,
	         list( ctlpt( E3,  0.0,  1.0, 0.0 ),
		       ctlpt( E3,  1.0,  1.0, 1.0 ),
		       ctlpt( E3,  1.0, -1.0, 0.0 ),
		       ctlpt( E3, -1.0, -1.0, 0.0 ) ),
	         list( kv_periodic ) );

s = swpsclsrf( Cross2, Axis, 0.2, off, 3 );
if ( display == on, interact( list( axes, s ) ) );

s = swpsclsrf( list( Cross,
		     Cross2,
		     Cross2,
		     Cross,
		     Cross,
		     Cross2,
		     Cross2,
		     Cross,
		     Cross,
		     Cross2,
		     Cross2,
		     Cross,
		     Cross ),
	       Axis, 0.2, off, 3 );
if ( display == on, interact( list( axes, s ) ) );

free( Cross );
free( Cross2 );
free( Frame );
free( Axis );
free( s );

#
# Boolean sum examples.
#
cbzr1 = cbezier( list( ctlpt( E3, -0.2,  0.1,  0.5 ),
		       ctlpt( E3,  0.0,  0.5,  1.0 ),
		       ctlpt( E3,  0.1,  1.0, -0.2 ) ) );
cbzr2 = cbezier( list( ctlpt( E3,  1.0,  0.0, -0.3 ),
		       ctlpt( E3,  0.8,  0.5, -1.0 ),
		       ctlpt( E3,  1.0,  1.0,  0.2 ) ) );
cbzr3 = cbezier( list( ctlpt( E3, -0.2,  0.1,  0.5 ),
		       ctlpt( E3,  0.5,  0.0, -1.0 ),
		       ctlpt( E3,  1.0,  0.0, -0.3 ) ) );
cbzr4 = cbezier( list( ctlpt( E3,  0.1,  1.0, -0.2 ),
		       ctlpt( E3,  0.5,  1.0,  1.0 ),
		       ctlpt( E3,  1.0,  1.0,  0.2 ) ) );

s = boolsum( cbzr1, cbzr2, cbzr3, cbzr4 );
color( s, green );
if ( display == on, view_mat = save_mat * ty( 0.25 ):
		    interact( list( view_mat, axes,
				    cbzr1, cbzr2, cbzr3, cbzr4, s ) ) );

cbzr1 = cbezier( list( ctlpt( E3, -0.1,  0.1,  0.2 ),
		       ctlpt( E3,  0.0,  0.5,  1.0 ),
		       ctlpt( E3,  0.1,  1.0,  0.2 ) ) );
cbzr2 = cbezier( list( ctlpt( E3,  1.0,  0.2, -0.1 ),
		       ctlpt( E3,  1.0,  0.5, -1.0 ),
		       ctlpt( E3,  1.1,  1.1,  0.1 ) ) );
cbzr3 = cbezier( list( ctlpt( E3, -0.1,  0.1,  0.2 ),
		       ctlpt( E3,  0.2,  0.1, -1.0 ),
		       ctlpt( E3,  0.4,  0.0,  2.0 ),
		       ctlpt( E3,  0.5, -0.1, -1.0 ),
		       ctlpt( E3,  1.0,  0.2, -0.1 ) ) );
cbzr4 = cbezier( list( ctlpt( E3,  0.1,  1.0,  0.2 ),
		       ctlpt( E3,  0.5,  0.8,  1.0 ),
		       ctlpt( E3,  0.7,  0.9, -2.0 ),
		       ctlpt( E3,  0.8,  1.0,  1.0 ),
		       ctlpt( E3,  1.1,  1.1,  0.1 ) ) );
s = boolsum( cbzr1, cbzr2, cbzr3, cbzr4 );
color( s, green );
if ( display == on, interact( list( axes, cbzr1, cbzr2, cbzr3, cbzr4, s ) ) );

cbzr1 = cbezier( list( ctlpt( E3, 0.1, 0.1, 0.1 ),
		       ctlpt( E3, 0.0, 0.5, 1.0 ),
		       ctlpt( E3, 0.4, 1.0, 0.4 ) ) );
cbzr2 = cbezier( list( ctlpt( E3, 1.0, 0.2, 0.2 ),
		       ctlpt( E3, 1.0, 0.5, -1.0 ),
		       ctlpt( E3, 1.0, 1.0, 0.3 ) ) );
cbsp3 = cbspline( 4,
		  list( ctlpt( E3, 0.1,  0.1, 0.1 ),
		        ctlpt( E3, 0.25, 0.0, -1.0 ),
		        ctlpt( E3, 0.5,  0.0, 2.0 ),
		        ctlpt( E3, 0.75, 0.0, -1.0 ),
		        ctlpt( E3, 1.0,  0.2, 0.2 ) ),
		  list( KV_OPEN ) );
cbsp4 = cbspline( 4,
		  list( ctlpt( E3, 0.4,  1.0, 0.4 ),
		        ctlpt( E3, 0.25, 1.0, 1.0 ),
		        ctlpt( E3, 0.5,  1.0, -2.0 ),
		        ctlpt( E3, 0.75, 1.0, 1.0 ),
		        ctlpt( E3, 1.0,  1.0, 0.3 ) ),
		  list( KV_OPEN ) );
s = boolsum( cbzr1, cbzr2, cbsp3, cbsp4 );
color( s, green );
if ( display == on, interact( list( axes, cbzr1, cbzr2, cbsp3, cbsp4, s ) ) );

free( cbzr1 );
free( cbzr2 );
free( cbzr3 );
free( cbzr4 );
free( cbsp3 );
free( cbsp4 );

#
# Boolean one examples.
#
s = boolone( circle( vector( 0.0, 0.0, 0.0 ), 1.0 ) );
color( s, green );
if ( display == on, view_mat = save_mat:
		    interact( list( view_mat, axes, s ) ) );

c1 = ctlpt( E3, 0, 0, 0.5 ) +
     ctlpt( E3, 1, 0, 0 ) +
     ctlpt( E3, 1, 1, 0.5 ) +
     ctlpt( E3, 0, 1, 0 ) +
     ctlpt( E3, 0, 0, 0.5 );
s = boolone( c1 );
color( s, green );
if ( display == on, view_mat = save_mat * rx( -10 ) * ry( 5 ):
		    interact( list( view_mat, axes, s, c1 ) ) );

#
# Surface from curves constructor.
#
c1 = cbspline( 3,
	       list( ctlpt( E3, 0.0, 0.0, 0.0 ),
		     ctlpt( E3, 1.0, 0.0, 0.0 ),
		     ctlpt( E3, 1.0, 1.0, 0.0 ) ),
	       list( KV_OPEN ) );
c2 = cbspline( 3,
	       list( ctlpt( E3, 0.0, 0.0, 1.0 ),
		     ctlpt( E3, 1.0, 0.0, 1.0 ),
		     ctlpt( E3, 1.0, 2.0, 1.0 ) ),
	       list( KV_OPEN ) );
c3 = cbspline( 3,
	       list( ctlpt( E3, 0.0, 0.0, 1.5 ),
		     ctlpt( E3, 2.0, 0.0, 1.5 ),
		     ctlpt( E3, 1.0, 0.5, 1.5 ),
		     ctlpt( E3, 1.0, 1.0, 1.5 ) ),
	       list( KV_OPEN ) );
c4 = cbspline( 3,
	       list( ctlpt( E3, 0.0, 0.0, 2.5 ),
		     ctlpt( E3, 1.0, 0.0, 2.5 ),
		     ctlpt( E3, 1.0, 1.0, 2.5 ) ),
	       list( KV_OPEN ) );

s = sfromcrvs( list( c1, c2, c3, c4 ), 2, KV_OPEN );
color( s, green );

if ( display == on, view_mat = save_mat * sc( 0.6 ) * ty( -0.3 ):
		    interact( list( view_mat, c1, c2, c3, c4, s ) ) );
s = sfromcrvs( list( c1, c2, c3, c4 ), 4, KV_OPEN );
color( s, green );
if ( display == on, interact( list( c1, c2, c3, c4, s ) ) );
free( c1 );
free( c2 );
free( c3 );
free( c4 );
free( s );

#
# Offset and adaptive/variable offset, with global tolerance.
#
cpawn = cbspline( 4,
	list( ctlpt( E2, 0.95, 0.05 ),
	      ctlpt( E2, 0.95, 0.76 ),
	      ctlpt( E2, 0.30, 1.52 ),
	      ctlpt( E2, 0.30, 1.90 ),
	      ctlpt( E2, 0.50, 2.09 ),
	      ctlpt( E2, 0.72, 2.24 ),
	      ctlpt( E2, 0.72, 2.32 ),
	      ctlpt( E2, 0.38, 2.50 ),
	      ctlpt( E2, 0.42, 2.70 ),
	      ctlpt( E2, 0.57, 2.81 ),
	      ctlpt( E2, 0.57, 3.42 ),
	      ctlpt( E2, 0.19, 3.57 ),
	      ctlpt( E2, 0.00, 3.57 ) ),
	 list( KV_OPEN ) );
color( cpawn, white );
c1 = offset( cpawn, 0.5, 0.05, off );
color( c1, magenta );
c2 = aoffset( cpawn, 0.5, 0.05, false, off );
color( c2, green );
c3 = aoffset( cpawn, 0.5, 0.005, true, off );
color( c3, yellow );
c4 = loffset( cpawn, 0.5, 100, 20, 4 );
color( c4, cyan );
#c5 = moffset( cpawn, 0.5, 10 );
#color( c5, red );
if ( display == on, view_mat = sc( 0.45 ) * tx( -0.3 ) * ty( -0.9 ):
		    interact( list( view_mat, c1, c2, c3, c4, cpawn ) ) );
save( "offset1", list( c1, c2, c3, c4, cpawn ) );

# Try variable offset...

voff = cbspline( 3, list( ctlpt( E1,  0.25 ),
			  ctlpt( E1, -0.25 ),
			  ctlpt( E1,  0.25 ),
			  ctlpt( E1, -0.25 ),
			  ctlpt( E1,  0.25 ),
			  ctlpt( E1, -0.25 ),
			  ctlpt( E1,  0.25 ) ),
		 list( kv_float ) ) * tx( 0.125 );
voff = coerce( voff, kv_open );

c1 = offset( cpawn, voff, 0.05, off );
color( c1, magenta );
c2 = offset( cpawn, voff * tx( -0.125 ), 0.05, on );
color( c2, green );
c3 = offset( cpawn, voff * tx( -0.25 ), 0.05, on );
color( c3, yellow );
c4 = aoffset( cpawn, voff, 0.02, false, off );
color( c4, cyan );
c5 = aoffset( cpawn, voff * sx( -1 ), 0.02, false, on );
color( c5, red );

if ( display == on, view_mat = sc( 0.45 ) * tx( -0.3 ) * ty( -0.9 ):
		    interact( list( view_mat, c1, c2, c3, c4, c5, cpawn ) ) );
save( "offset2", list( c1, c2, c3, c4, c5, cpawn ) );

free( voff );
free( cpawn );
free( c1 );
free( c2 );
free( c3 );
free( c4 );
free( c5 );

#
# Zero set of curves and inflection points
#
cbsp = list( ctlpt( E2, 1.0, 1.0 ),
	     ctlpt( E2,-1.0, 1.0 ),
	     ctlpt( E2,-1.0,-1.0 ),
	     ctlpt( E2, 1.5, 1.0 ),
	     ctlpt( E2,-1.0,-1.4 ),
	     ctlpt( E2, 1.0, 0.1 ) );
cb = cbspline( 4, cbsp, list( KV_OPEN ) );
free( cbsp );

xzeros = czeros( cb, 0.001, 1 );
pt_xzeros = nil();
pt = nil();
for ( i = 1, 1, sizeof( xzeros ),
	pt = ceval( cb, nth( xzeros, i ) ):
	snoc( pt, pt_xzeros )
    );
if ( display == on, view_mat = save_mat:
		    interact( list( view_mat, axes, cb, pt_xzeros ) ) );
free( xzeros );
free( pt_xzeros );

yzeros = czeros( cb, 0.001, 2 );
pt_yzeros = nil();
pt = nil();
for ( i = 1, 1, sizeof( yzeros ),
	pt = ceval( cb, nth( yzeros, i ) ):
	snoc( pt, pt_yzeros )
    );
if ( display == on, interact( list( axes, cb, pt_yzeros ) ) );
free( yzeros );
free( pt_yzeros );

xextremes = cextremes( cb, 0.001, 1 );
pt_xextremes = nil();
pt = nil();
for ( i = 1, 1, sizeof( xextremes ),
	pt = ceval( cb, nth( xextremes, i ) ):
	snoc( pt, pt_xextremes )
    );
if ( display == on, interact( list( axes, cb, pt_xextremes ) ) );
free( xextremes );
free( pt_xextremes );

yextremes = cextremes( cb, 0.001, 2 );
pt_yextremes = nil();
pt = nil();
for ( i = 1, 1, sizeof( yextremes ),
	pt = ceval( cb, nth( yextremes, i ) ):
	snoc( pt, pt_yextremes )
    );
if ( display == on, interact( list( axes, cb, pt_yextremes ) ) );
free( yextremes );
free( pt_yextremes );

inflect = cinflect( cb, 0.001, 2 );
pt_inflect = nil();
pt = nil();
for ( i = 1, 1, sizeof( inflect ),
	pt = ceval( cb, nth( inflect, i ) ):
	snoc( pt, pt_inflect )
    );
if ( display == on, interact( list( axes, cb, pt_inflect ) ) );
free( inflect );

Crvs = cinflect( cb, 0.001, 3 );
for ( i = 1, 1, sizeof( Crvs ),
    attrib( nref( Crvs, i ), "rgb", GenRandomColor() ) );

if ( display == on, interact( list( view_mat, Crvs, pt_inflect ) ) );

free( pt_inflect );
free( Crvs );
free( cb );

#
# Computation of extremum of curvature for curves
#
crv = cbezier( list( ctlpt( E3, -0.5,  0.5,  0.5 ),
		     ctlpt( E3, -0.5, -0.6,  0.5 ),
		     ctlpt( E3,  0.0,  1.0, -1.0 ),
		     ctlpt( E3,  0.5, -0.5,  0.5 ) ) );

crvtr = ccrvtr( crv, 0.001, 2 );
pt_crvtr = nil();
pt = nil();
for ( ( i = 1 ), 1, sizeof( crvtr ),
	(
	  ( pt = ceval( crv, nth( crvtr, i ) ) ):
	  snoc( pt, pt_crvtr )
	)
    );
if ( display == on, interact( list( crv, pt_crvtr ) ) );

Crvs = ccrvtr( crv, 0.001, 3 );
for ( i = 1, 1, sizeof( Crvs ),
    attrib( nref( Crvs, i ), "rgb", GenRandomColor() ) );

if ( display == on, interact( list( view_mat, Crvs, pt_crvtr ) ) );

crv = cbspline( 4,
		list( ctlpt( E2, -1.0,  0.5 ),
		      ctlpt( E2,  0.0,  0.5 ),
		      ctlpt( E2,  0.0,  0.0 ),
		      ctlpt( E2, -1.0, -1.0 ),
		      ctlpt( E2,  1.0, -1.0 ),
		      ctlpt( E2,  1.0, -0.7 ),
		      ctlpt( E2,  0.0,  1.0 ),
		      ctlpt( E2,  1.0,  1.0 ) ),
		list( KV_OPEN ) );

crvtr = ccrvtr( crv, 0.001, 2 );
pt_crvtr = nil();
pt = nil();
for ( ( i = 1 ), 1, sizeof( crvtr ),
	(
	  ( pt = ceval( crv, nth( crvtr, i ) ) ):
	  snoc( pt, pt_crvtr )
	)
    );
free( pt );
free( crvtr );
if ( display == on, view_mat = save_mat * sc( 0.8 ) * tx( 0.05 ):
		    interact( list( view_mat, crv, pt_crvtr ) ) );

Crvs = ccrvtr( crv, 0.001, 3 );
for ( i = 1, 1, sizeof( Crvs ),
    attrib( nref( Crvs, i ), "rgb", GenRandomColor() ) );

if ( display == on, interact( list( view_mat, Crvs, pt_crvtr ) ) );

free( Crvs );
free( i );
free( pt_crvtr );

#
# Computation of evolute curves
#
pl = nil();
for ( x = 0, 1, 20,
	snoc( point( cos(x * Pi / 10), sin(x * Pi / 10), 0 ), pl )
    );
free( x );

#
# Evolute of a very crude approximation of a circle
#
crv = cinterp( pl, 3, 4, PARAM_UNIFORM, false );
cev = evolute( crv );
color( cev, red );
color( crv, green );
if ( display == on, view_mat = save_mat * sc( 0.25 ) * tx( -0.3 ) * ty( -0.3 ):
		    interact( list( view_mat, crv, cev ) ) );

#
# Better approximation of a circle (evolute should degen. to a point if exact)
#
crv = cinterp( pl, 3, 8, PARAM_UNIFORM, false );
cev = evolute( crv );
color( cev, red );
color( crv, green );
if ( display == on, view_mat = save_mat:
		    interact( list( view_mat, crv, cev ) ) );
free( pl );

#
# And an evolute of a space curve.
#
crv = cbspline( 3,
		list( ctlpt( E3, -1.0,  0.1,  0.2 ),
		      ctlpt( E3, -0.1,  1.0,  0.1 ),
		      ctlpt( E3,  0.1,  0.1,  1.0 ),
		      ctlpt( E3,  1.0,  0.1,  0.1 ),
		      ctlpt( E3,  0.1,  1.0,  0.2 ) ),
		list( KV_OPEN ) );
cev = evolute( crv );
color( cev, red );
color( crv, green );
if ( display == on, view_mat = save_mat * sc( 0.5 ) * ty( -0.3 ):
		    interact( list( view_mat, crv, cev ) ) );

free( cev );
free( crv );

free( s45 );
free( display );

resolution = save_res;
view_mat = save_mat;

IProd = iritstate("InterpProd", IProd );
free( IProd );