#
# A demonstration of the capabilities of the MVEXPLICIT function.
#
#					Gershon, May 2002
#

CoerceToBezSrf = function( Mv, UMin, UMax, VMin, VMax ): MvBzr: SrfBzr:
    MvBzr = coerce( MV, bezier_type ):
    SrfBzr = coerce( coerce( MvBzr, surface_type ), E3 )
						* rotx( -90 ) * roty( -90 ):
    SrfBzr = sregion( sregion( SrfBzr, row, UMin, UMax ), col, VMin, VMax ):
    return = SrfBzr;

#
# Simple quadratic surfaces:
#
Parab = CoerceToBezSrf( mvexplicit( "A^2 + B^2 + 0.5" ), -1, 1, -1, 1 );
color( Parab, yellow );
Saddl1 = CoerceToBezSrf( mvexplicit( "A^2 - B^2" ), -1, 1, -1, 1 );
color( Saddl1, green );
Saddl2 = CoerceToBezSrf( mvexplicit( "A * B - 1.5" ), -1, 1, -1, 1 );
color( Saddl2, cyan );

interact( list( axes, Parab, Saddl1, Saddl2 ) * sc( 0.35 ) );
save( "mvexplc1", list( axes, Parab, Saddl1, Saddl2 ) );

free( Parab );
free( Saddl1 );
free( Saddl2 );

#
# Monkey saddle.
#

Monkey = CoerceToBezSrf( mvexplicit( "A^3 - 3 * A * B^2" ), -1, 1, -1, 1 );
color( Monkey, yellow );

Pln = ruledSrf( ctlpt( E3, 0, -5, -5 ) + ctlpt( E3, 0, -5,  5 ),
		ctlpt( E3, 0,  5, -5 ) + ctlpt( E3, 0,  5,  5 ) );

ICrv = iritstate( "intercrv", true );
IntCrvs = list( Monkey * ( Pln * rz( 90 ) ),
		Monkey * ( Pln * rz( 90 + 60 ) ),
		Monkey * ( Pln * rz( 90 + 120 ) ) );
adwidth( IntCrvs, 2 );
color( IntCrvs, red );
ICrv = iritstate( "intercrv", ICrv );
free( ICrv );
free( Pln );

interact( list( axes, IntCrvs, Monkey ) * sc( 0.35 ) );
save( "mvexplc2", list( axes, IntCrvs, Monkey ) );

free( Monkey );
free( IntCrvs );

#
# Enneper minimal surface.
#

Enneper = ffmerge( list( mvexplicit( "A - A^3 / 3 + A * B^2" ),
			 mvexplicit( "    B^3 / 3 - A^2 * B - B" ),
			 mvexplicit( "A^2 - B^2" ) ),
		   E3 );
Enneper = coerce( coerce( Enneper, bezier_type ), surface_type );
Enneper = sregion( sregion( Enneper, row,  -1, 1 ), col, -1, 1 );

b = bbox( smean( Enneper, false ) );
printf( "Mean curvature sqaure computed for the Enneper surface spans %.14f to %.14f\\n",
	list( nth( b, 1 ), nth( b, 2 ) ) );
free( b );

save( "mvexplc3", list( axes, Enneper ) );

free( Enneper );

#
# Solve non linear equations.
#

# A^2 + B^2 = 1
M1 = coerce( mvexplicit( "A^2 + B^2 - 1" ), bezier_type );
M1 = mregion( mregion( M1, 0, -3, 3 ), 1, -3, 3 );

# (A+1)^2 + B^2 = 1
M2 = coerce( mvexplicit( "A^2 + 2 * A + 1 + B^2 - 1" ), bezier_type );
M2 = mregion( mregion( M2, 0, -3, 3 ), 1, -3, 3 );

z = mzero( list( M1, M2 ), 0.01, 1e-6 ) * sc( 6 ) * tx( -3 ) * ty( -3 );
color( z, yellow );

interact( list( circle( vector( 0, 0, 0 ), 1 ),
		circle( vector( -1, 0, 0 ), 1 ),
		z ) * sc( 0.6 ) );


# A^2 + B^2 = 1
M1 = coerce( mvexplicit( "A^2 + B^2 - 1" ), bezier_type );
M1 = mregion( mregion( M1, 0, -3, 3 ), 1, -3, 3 );

# 4A^2 + B^2/4 = 1
M2 = coerce( mvexplicit( "4 * A^2 + B^2 / 4 - 1" ), bezier_type );
M2 = mregion( mregion( M2, 0, -3, 3 ), 1, -3, 3 );

z = mzero( list( M1, M2 ), 0.01, 1e-6 ) * sc( 6 ) * tx( -3 ) * ty( -3 );
color( z, yellow );
free( M1 );
free( M2 );

interact( list( circle( vector( 0, 0, 0 ), 1 ),
		circle( vector( 0, 0, 0 ), 1 ) * sx( 0.5 ) * sy( 2 ),
		z ) * sc( 0.6 ) );

#
# Intersect the three surfaces we started with - paraboloid and two saddles.
#
Parab = CoerceToBezSrf( mvexplicit( "A^2 + B^2 - 0.5" ), -1, 1, -1, 1 );
color( Parab, yellow );
Saddl1 = CoerceToBezSrf( mvexplicit( "A^2 - B^2" ), -1, 1, -1, 1 );
color( Saddl1, green );
Saddl2 = CoerceToBezSrf( mvexplicit( "A * B" ), -1, 1, -1, 1 );
color( Saddl2, cyan );

Pln = ruledSrf( ctlpt( E3, 0, -1, -1 ) + ctlpt( E3, 0, -1,  1 ),
		ctlpt( E3, 0,  1, -1 ) + ctlpt( E3, 0,  1,  1 ) ) * ry( 90 );
color( Pln, red );
attrib( Pln, "transp", 0.95 );

# Parab = Saddl1
M12 = coerce( mvexplicit( "A^2 + B^2 - 0.5 - A^2 + B^2" ), bezier_type );
M12 = mregion( mregion( M12, 0, -1, 1 ), 1, -1, 1 );
# Parab = Saddl1
M23 = coerce( mvexplicit( "A^2 - B^2 - A * B" ), bezier_type );
M23 = mregion( mregion( M23, 0, -1, 1 ), 1, -1, 1 );

# Intersection of the two saddles:
z = mzero( list( M23 ), 0.01, 1e-6 ) * sc( 2 ) * tx( -1 ) * ty( -1 );
color( z, magenta );
interact( list( Saddl1, Saddl2, Pln, z ) * sc( 0.35 ) );

# Intersection of the three surfaces:
z = mzero( list( M12, M23 ), 0.01, 1e-6 ) * sc( 2 ) * tx( -1 ) * ty( -1 );
color( z, magenta );
interact( list( Parab, Saddl1, Saddl2, Pln, z ) * sc( 0.35 ) );
save( "mvexplc4", list( Parab, Saddl1, Saddl2, Pln, z ) );

free( M12 );
free( M23 );
free( z );

free( Pln );
free( Parab );
free( Saddl1 );
free( Saddl2 );
