#
# Some routines to test bspline curves/surfaces.
#

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

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

s45 = sin( pi / 4 );
cbsp = list( ctlpt( P2, 1.0, 1.0, 0.0 ),
	     ctlpt( P2, s45, s45, s45 ),
	     ctlpt( P2, 1.0, 0.0, 1.0 ) );
free( s45 );
sbsp = list ( list( ctlpt( E3, 0.1, 0.0, 1.0 ),
	            ctlpt( E3, 0.3, 1.0, 0.0 ),
	            ctlpt( E3, 0.0, 2.0, 1.0 ) ),
              list( ctlpt( E3, 1.1, 0.0, 0.0 ),
	            ctlpt( E3, 1.3, 1.0, 2.0 ),
	            ctlpt( E3, 1.0, 2.0, 0.0 ) ),
              list( ctlpt( E3, 2.1, 0.0, 2.0 ),
	            ctlpt( E3, 2.3, 1.0, 0.0 ),
	            ctlpt( E3, 2.0, 2.0, 2.0 ) ),
              list( ctlpt( E3, 3.1, 0.0, 0.0 ),
	            ctlpt( E3, 3.3, 1.0, 2.0 ),
	            ctlpt( E3, 3.0, 2.0, 0.0 ) ),
              list( ctlpt( E3, 4.1, 0.0, 1.0 ),
	            ctlpt( E3, 4.3, 1.0, 0.0 ),
	            ctlpt( E3, 4.0, 2.0, 1.0 ) ) );

cb = cbspline( 3, cbsp, list( 0, 0, 0, 1, 1, 1 ) )
	* scale( vector( 0.7, 1.4, 1.0 ) );
free( cbsp );
color( cb, red );
sb = sbspline( 3, 3, sbsp, list( list( 1, 1, 1, 2, 2, 2 ),
				 list( 3, 3, 3, 4, 5, 6, 6, 6 ) ) );
color( sb, red );

fforder( cb );
ffmsize( cb );
ffkntvec( cb );
ffctlpts( cb );
fforder( sb );
ffmsize( sb );
ffkntvec( sb );
ffctlpts( sb );

if ( display == on, interact( list( axes, cb, sb ) ):
		    viewstate( "DSrfMesh" ):
		    pause():
		    viewstate( "DSrfMesh" ):
		    pause() );

#
# Float end condition conversion to open end condition using subdivision.
#
sbsp = list ( list( ctlpt( E3, 0.1, 0.0, 1.0 ),
	            ctlpt( E3, 0.3, 1.0, 0.0 ),
	            ctlpt( E3, 0.0, 2.0, 1.0 ) ),
              list( ctlpt( E3, 1.1, 0.0, 0.0 ),
	            ctlpt( E3, 1.3, 1.0, 2.0 ),
	            ctlpt( E3, 1.0, 2.0, 0.0 ) ),
              list( ctlpt( E3, 2.1, 0.0, 2.0 ),
	            ctlpt( E3, 2.3, 1.0, 0.0 ),
	            ctlpt( E3, 2.0, 2.0, 2.0 ) ),
              list( ctlpt( E3, 3.1, 0.0, 0.0 ),
	            ctlpt( E3, 3.3, 1.0, 2.0 ),
	            ctlpt( E3, 3.0, 2.0, 0.0 ) ),
              list( ctlpt( E3, 4.1, 0.0, 1.0 ),
	            ctlpt( E3, 4.3, 1.0, 0.0 ),
	            ctlpt( E3, 4.0, 2.0, 1.0 ) ) );
sb = sbspline( 3, 4, sbsp, list( list( 1, 2, 3, 4, 5, 6 ),
				 list( 1, 2, 3, 4, 5, 6, 7, 8, 9 ) ) );
color(sb, red);
iritState("DumpLevel", 9);
sbOpenRow = nth(sdivide(nth(sdivide( sb, row, 4), 2), row, 6), 1);
sbOpenCol = nth(sdivide(nth(sdivide( sbOpenRow, col, 3), 2), col, 4), 1);
color(sbOpenCol,cyan);
if ( display == on, view( list( sb, sbOpenCol ), on ):
		    viewstate( "DSrfMesh" ):
		    pause():
		    viewstate( "DSrfMesh" ):
		    pause() );
free(sb);
free(sbOpenRow);
free(sbOpenCol);

#
# Curve refinement.
#
cb_ref = crefine( cb, false, list( 0.25, 0.5, 0.75 ) );
color( cb_ref, yellow );
if ( display == on, interact( list( axes, cb, cb_ref ) ) );
free( cb_ref );

#
# Knot substitution (One internal knot is moved from 0.1 to 0.9).
#
cb_ref = crefine( cb, false, list( 0.5 ) );
cb_all = list( axes );
for ( t = 0.1, 0.1, 0.9,
	cb1 = crefine( cb_ref, true, list( 0, 0, 0, t, 1, 1, 1 ) ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
if ( display == on, interact( cb_all ) );
free( cb_ref );
free( cb_all );

#
# Curve subdivision.
#
cb_lst = cdivide( cb, 0.5 );
cb1 = nth( cb_lst, 1 );
color( cb1, green );
cb2 = nth( cb_lst, 2 );
color( cb2, yellow );
free( cb_lst );
if ( display == on, interact( list( axes, cb, cb1, cb2 ) ) );
free( cb1 );
free( cb2 );

#
# Region from curve.
#
cbr1 = cregion( cb, 0.3, 0.6 );
color( cbr1, yellow );
cbr2 = cregion( cb, 0.5, 1.0 );
color( cbr2, green );
cbr3 = cregion( cb, 0.3, 0.0 );
color( cbr3, blue );
if ( display == on, interact( list( cb, cbr1, cbr2, cbr3 ) ) );
free( cbr1 ); free( cbr2 ); free( cbr3 );

#
# Surface subdivision and merging.
#
sb = sbspline( 3, 3, sbsp, list( list( 1, 1, 1, 2, 2, 2 ),
				 list( 3, 3, 3, 4, 5, 6, 6, 6 ) ) );
color( sb, red );
sb_lst = sdivide( sb, COL, 1.4 );
sb1 = nth( sb_lst, 1 );
color( sb1, green );
sb2 = nth( sb_lst, 2 );
color( sb2, yellow );
free( sb_lst );
if ( display == on, interact( list( axes, sb, sb1, sb2 ) ) );
sbm = smerge( sb1, sb2, COL, 1 );
if ( display == on, interact( list( axes, sbm, sb1, sb2 ) ) );
sbm = smerge( sb1 * trans( vector( 0.0, -0.5, 0.0 ) ),
	      sb2 * trans( vector( 0.0,  0.5, 0.0 ) ), COL, 0 );
if ( display == on, interact( list( axes, sbm, sb1, sb2 ) ) );
free( sb1 );
free( sb2 );
free( sbsp );

sb_lst = sdivide( sb, ROW, 4.8 );
sb1 = nth( sb_lst, 1 );
color( sb1, green );
sb2 = nth( sb_lst, 2 );
color( sb2, yellow );
free( sb_lst );
if ( display == on, interact( list( axes, sb, sb1, sb2 ) ) );
sbm = smerge( sb1, sb2, ROW, 1 );
if ( display == on, interact( list( axes, sbm, sb1, sb2 ) ) );
sbm = smerge( sb1 * trans(vector( -0.5, 0.0, 0.0 ) ),
	      sb2 * trans(vector(  0.5, 0.0, 0.0 ) ), ROW, 0 );
if ( display == on, interact( list( axes, sbm, sb1, sb2 ) ) );
free( sbm );
free( sb1 );
free( sb2 );

#
# Region from surface.
#
sbr1 = sregion( sb, COL, 1.3, 1.6 );
color( sbr1, yellow );
sbr2 = sregion( sb, COL, 1.8, 2.0 );
color( sbr2, green );
sbr3 = sregion( sb, ROW, 4.0, 5.0 );
color( sbr3, blue );
interact( list( sb, sbr1, sbr2, sbr3 ) );
free( sbr1 ); free( sbr2 ); free( sbr3 );

#
# Derivative and normal curves/surfaces
#
dcb = cderive( cb );
if ( display == on, interact( list( axes, dcb, cb ) ) );
dsb1 = sderive( sb, ROW );
color( dsb1, magenta );
dsb2 = sderive( sb, COL );
color( dsb2, yellow );
if ( display == on, interact( list( axes, dsb1, dsb2, sb ) ) );
ncb = cnrmlcrv( cb );
color( ncb, yellow );
if ( display == on, interact( list( axes, ncb, cb ) ) );
nsb = snrmlsrf( sb );
color( nsb, yellow );
if ( display == on, interact( list( axes, nsb, sb ) ) );
free( dcb );
free( dsb1 );
free( dsb2 );
free( ncb );
free( nsb );


#
# Iso curves extraction from surface.
#
cb_all = list( axes );
snoc( sb, cb_all );
for ( t = 1.1, 0.1, 1.9,
	cb1 = csurface( sb, COL, t ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
for ( t = 3.1, 0.2, 5.9,
	cb1 = csurface( sb, ROW, t ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
if ( display == on, interact( cb_all ) );
free( cb_all );


#
# curves extraction from surface mesh. Note curves may be not on the surface.
#
cb_all = list( axes );
snoc( sb, cb_all );
for ( t = 0, 1, 2,
	cb1 = cmesh( sb, COL, t ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
for ( t = 0, 1, 4,
	cb1 = cmesh( sb, ROW, t ):
	color( cb1, green ):
	snoc( cb1, cb_all )
    );
free( t );
if ( display == on, interact( cb_all ) );
free( cb_all );

#
# convert into polygons/polylines (using default resolution).
#
iritState("DumpLevel", 1);
resolution;

p = gpolyline( list( sb, cb ), off );
if ( display == on, interact( list( p, axes ) ) );

p = gpolygon( sb, on );
if ( display == on, viewstate("DrawVNrml"):
		    interact( list( p, axes ) ):
		    viewstate("DrawVNrml") );

#
# reverse surface ( flip normals ).
#
q = gpolygon( -sb, on );
if ( display == on, viewstate("DrawVNrml"):
		    interact( list( q, axes ) ):
		    viewstate("DrawVNrml") );

free(p);
free(q);

#
# Offset approximation by translation of srf/crv in normal direction.
#
cbo = offset(cb, 0.1, 0.1, off);
if ( display == on, interact( list( axes, cb, cbo ) ) );
free(cbo);

sbo = offset(sb, 0.2, 0.1, off);
if ( display == on, interact( list( axes, sb, sbo ) ) );
free(sbo);

#
# Surface and Curve evaluation.
#
ceval( cb, 0.0 );
ceval( cb, 0.1 );
ceval( cb, 0.3 );
ceval( cb, 0.5 );
ceval( cb, 0.9 );
ceval( cb, 1.0 );

pause();

seval( sb, 1.0, 3.0 );
seval( sb, 1.1, 3.0 );
seval( sb, 1.3, 3.0 );
seval( sb, 1.5, 3.5 );
seval( sb, 1.9, 3.1 );
seval( sb, 1.0, 4.0 );
seval( sb, 1.5, 4.0 );

pause();

#
#
# Surface and Curve tangents.
#
ctangent( cb, 0.0 );
ctangent( cb, 0.1 );
ctangent( cb, 0.3 );
ctangent( cb, 0.5 );
ctangent( cb, 0.9 );
ctangent( cb, 1.0 );

pause();

stangent( sb, ROW, 1.0, 3.0 );
stangent( sb, COL, 1.1, 3.0 );
stangent( sb, ROW, 1.3, 3.0 );
stangent( sb, COL, 1.5, 3.5 );
stangent( sb, ROW, 1.9, 3.1 );
stangent( sb, COL, 1.0, 4.0 );
stangent( sb, COL, 1.5, 4.0 );

#
#
# Surface and Curve normals.
#
cnormal( cb, 0.0 );
cnormal( cb, 0.1 );
cnormal( cb, 0.3 );
cnormal( cb, 0.5 );
cnormal( cb, 0.9 );
cnormal( cb, 1.0 );

snormal( sb, 1.0, 3.0 );
snormal( sb, 1.1, 3.0 );
snormal( sb, 1.3, 3.0 );
snormal( sb, 1.5, 3.5 );
snormal( sb, 1.9, 3.1 );
snormal( sb, 1.0, 4.0 );
snormal( sb, 1.5, 4.0 );

#
# Compute moments:
#
a = circle( vector( 0, 0, 0 ), 1 );
a = cregion( a, 0, 1 );
moment( a, 0 );
moment( a, 1 );

a = cregion( a, 0, 1 ) * rz( 45 );
moment( a, 0 );
moment( a, 1 );
free( a );

pause();

#
# save("cb", cb);
# save("sb", sb);
#
# cb1 = load("cb.crv");
# sb1 = load("sb.srf");
#
# save("cb1", cb1);
# save("sb1", sb1);
#

free( cb1 );
free( cb );
free( sb );
free( display );