#
# Some routines to test surface/curve operators
#

#
# 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;

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

s45 = sin( pi / 4 );

#  
# Circular constructors
#

circ = circle( vector( 0.25, 0.5, 0.5 ), 1.5 );
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 ) );

dumplvl = 9; pdomain( circ ); dumplvl = 1;
if ( display == on, interact( list( axes, circ, arc1, arc2 ) ) );

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

#
# 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 ) );
dumplvl = 9; pdomain( ruled ); dumplvl = 1;
if ( display == on, interact( list( 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 ) ) );
if ( display == on, interact( list( axes, skewcyl ) ) );
free( skewcyl );

skew2cyl = ruledsrf( circ*rotx( 20 ),
		     circ*rotx( -20 )* trans( vector( 0.0, 0.0, 1.0 ) ) );
if ( display == on, interact( list( axes, skew2cyl ) ) );
free( skew2cyl );

drawctlpt = 0;
free( arc1 ); free( arc2 ); free( arc3 ); free( circ );

#
# Curve constructors.
#
crv1 = ctlpt( E3, 0.5, 0.0, 1.0 ) +
       ctlpt( E3, 0.5, 0.0, 1.2 ) +
       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 ) ):
		    viewstate("FinrAprx"):
		    pause(1):
		    viewstate("CrsrAprx") );

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( E2,  0, 2 ),
		      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( P3,  1.0,  1.0, 0.1,  0.0 ),
		       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 ) ) );

srf1 = sbezier( list( list( ctlpt( E2, 0.0, 0.0 ),
		            ctlpt( E2, 0.0, 1.0 ),
		            ctlpt( E2, 0.0, 2.0 ) ),
	              list( ctlpt( E2, 1.0, 0.0 ),
		            ctlpt( E2, 1.0, 1.0 ),
		            ctlpt( E2, 1.0, 2.0 ) ) ) );
srf2 = 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_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, interact( list( 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 ) ) );

#
# Create floating end condition examples (although barely used!).
#
drawctlpt = red;
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( crv1f, on ):
		    viewstate("DSrfMesh"):
		    pause( 1 ):
		    viewstate("DSrfMesh"):
		    pause( 1 ) );

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( srf1f, on ):
		    viewstate("DSrfMesh"):
		    pause( 1 ):
		    viewstate("DSrfMesh"):
		    pause( 1 ) );

free(crv1f);
free(srf1f);

#
# 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, interact( 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 )
    );
if ( display == on, interact( sb_all ) );
free( sb_all ); free( sb ); free( sb1 );

#
# Curve order raise:
#
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 ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
if ( display == on, view( cb_all, on ):
		    viewstate("DSrfMesh"):
		    pause( 1 ):
		    viewstate("DSrfMesh"):
		    pause( 1 ) );
free( cb_all ); free( cb ); free( cb1 );

drawctlpt = 0;

#
# 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 ) );

if ( display == on, interact( list( axes, s, cbzr ) ) );
free( cbzr ); free( s );

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

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

resolution = 10;
ps = gpolygon( s );
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, interact( list( sb, 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, interact( list( sp, halfcirc ) ) );

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 = surfrev( gcross );
if ( display == on, interact( list( axes, glass, gcross ) ) );

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

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.5, 1.5, 2.5, 3.5 ) );

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, 0.2, off );
if ( display == on, interact( list( 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 = sweepsrf( ccross, arc1, scalecrv, off );
if ( display == on, interact( list( axes, s ) ) );

circ = circle( vector( 0.0, 0.0, 0.0 ), 1.0 );

s = sweepsrf( circ, circ, 0.25, off );
if ( display == on, interact( list( axes, s ) ) );

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

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

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

s = sweepsrf( ccross, sweep_axis, 0.2, off );
if ( display == on, interact( list( axes, s ) ) );

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

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

free( ccross );
free( cross );
free( circ );
free( arc1 );

#
# 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, interact( list( axes, s, cbzr1, cbzr2, cbzr3, cbzr4 ) ) );

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, s, cbzr1, cbzr2, cbzr3, cbzr4 ) ) );

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, s, cbzr1, cbzr2, cbsp3, cbsp4 ) ) );

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, interact( list( 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, interact( list( 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 );
color( s, green );
if ( display == on, interact( list( s, c1, c2, c3, c4 ) ) );
s = sfromcrvs( list( c1, c2, c3, c4 ), 4 );
color( s, green );
if ( display == on, interact( list( s, c1, c2, c3, c4 ) ) );
free( c1 ); free( c2 ); free( c3 ); free( c4 );

#
# Offset and adaptive, with global tolerance, aoffset.
#
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 ) );
c1 = offset(cpawn, 0.5);
color(c1, magenta);
c2 = aoffset(cpawn, 0.5, 0.05, false);
color(c2, green);
c3 = aoffset(cpawn, 0.5, 0.05, true);
color(c3, yellow);
if ( display == on, interact( list( c1, c2, c3, cpawn ) ) );
free(c1);
free(c2);
free(c3);

#
# 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 ) );

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

yzeros = czeros( cb, 0.001, 2 );
pt_yzeros = nil();
pt = nil();
for ( i = 1, 1, listsize( 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, listsize( 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, listsize( 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 );
pt_inflect = nil();
pt = nil();
for ( i = 1, 1, listsize( inflect ),
	pt = ceval( cb, nth( inflect, i ) ):
	snoc( pt, pt_inflect )
    );
if ( display == on, interact( list( axes, cb, pt_inflect ) ) );
free(inflect);
free(pt_inflect);

#
# 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 );
pt_crvtr = nil();
pt = nil();
for ( ( i = 1 ), 1, listsize( crvtr ),
	(
	  ( pt = ceval( crv, nth( crvtr, i ) ) ):
	  snoc( pt, pt_crvtr )
	)
    );
if ( display == on, interact( list( crv, 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 );
pt_crvtr = nil();
pt = nil();
for ( ( i = 1 ), 1, listsize( crvtr ),
	(
	  ( pt = ceval( crv, nth( crvtr, i ) ) ):
	  snoc( pt, pt_crvtr )
	)
    );
if ( display == on, interact( list( crv, pt_crvtr ) ) );

display = off;

resolution = save_res;
drawctlpt = 0;