iritState("EchoSource", false);
#
# Init file for IRIT solid modeller.
#
iritstate( "FloatFrmt", "%8.6lg" );

demo = procedure():
    include("demo.irt");

pause = procedure():t:
    printf("Press return to continue:", nil()):
    t = getline(string_type);

#
# Some simple functions.
#
min = function(x, y):
    if (x > y, return = y, return = x);
max = function(x, y):
    if (x > y, return = x, return = y);
sqr = function(x):
    return = x * x;
normalize = function(x):len:
    return = 1.0:
    if (thisobj("x") == vector_type || thisobj("x") == point_type,
	len = sqrt(x * x):
	return = point(coord(x, 0) / len,
		       coord(x, 1) / len,
		       coord(x, 2) / len)):
    if (thisobj("return") == numeric_type,
	printf("NORMALIZE: Can handle only vectors or points, found %8.6lDg\n",
	       list(x)));
midPoint = function(pt1, pt2):
    return = (pt1 + pt2) * 0.5;
InterpPoint = function(pt1, pt2, t):
    return = (pt1 * (1 - t) + pt2 * t);
edge2d = function( x1, y1, x2, y2 ):
    return = ctlpt( E2, x1, y1 ) + ctlpt( E2, x2, y2 );
edge3d = function( x1, y1, z1, x2, y2, z2 ):
    return = ctlpt( E3, x1, y1, z1 ) + ctlpt( E3, x2, y2, z2 );

#
# Sweep of circular cross section.
#
SwpCircSrf = function( AxisCrv, ScaleCrv, ScaleRefine ):
    return = SwpSclSrf( circle( vector( 0, 0, 0 ), 1 ),
			AxisCrv,
			ScaleCrv,
			off,
			ScaleRefine );

#
# Extractions of Control Polygon/Mesh/Points from a curve or a surface.
#
GetCtlPoints = function( Crv, Vecs ):i:p:
    return = nil():
    for ( i = 1, 1, sizeof( Crv ),
	if ( Vecs,
	     p = coerce( coord( Crv, i - 1 ), vector_type ),
	     p = coerce( coord( Crv, i - 1 ), point_type ) ):
	snoc( p, return )
    );
GetCtlPolygon = function( Crv ):
    return = poly( GetCtlPoints( Crv, false ), true );

GetCtlMeshPts = function( Srf, Vecs ):l:i:j:p:rsize:csize:
    return = nil():
    rsize = meshsize( Srf, ROW ):
    csize = meshsize( Srf, ROW ):
    for ( i = 1, 1, rsize,
        for ( j = 1, 1, csize,
	     if ( Vecs,
	          p = coerce( coord( Srf, (i - 1) * csize + j - 1 ),
			      vector_type ),
	          p = coerce( coord( Srf, (i - 1) * csize + j - 1 ),
			      point_type ) ):
	     snoc( p, return )
	)
    );
GetCtlMesh = function( Srf ):l:i:j:p:pl:rsize:csize:first:
    first = true:
    rsize = meshsize( Srf, ROW ):
    csize = meshsize( Srf, ROW ):
    for ( i = 1, 1, rsize,
	pl = nil():
        for ( j = 1, 1, csize,
	     p = coerce( coord( Srf, (i - 1) * csize + j - 1 ), vector_type ):
	     snoc( p, pl )
	):
	if ( first == true,
	     return = poly( pl, true ): first = false,
	     return = return + poly( pl, true ) )
    ):
    for ( j = 1, 1, csize,
	pl = nil():
        for ( i = 1, 1, rsize,
	     p = coerce( coord( Srf, (i - 1) * csize + j - 1 ), vector_type ):
	     snoc( p, pl )
	):
	return = return + poly( pl, true )
    );

#
# Approximate a (assumed to be) closed planar curve as a polygon
# with n vertices.
#
CnvrtCrvToPolygon = function( Crv, n ):ptl:pt:t:t0:t1:dt:
    ptl = nil():
    t0 = nth( pdomain( Crv ), 1 ):
    t1 = nth( pdomain( Crv ), 2 ):
    if ( n < 2, n = 2 ):
    dt = (t1 - t0) / (n - 1):
    for ( t = t0, dt, t1 - dt / 2,
	pt = ceval( Crv, t ):
	snoc( pt, ptl )
    ):
    return = poly( ptl, false );

#
# Primitives in freeform surface form.
#
planeSrf = function( x1, y1, x2, y2 ):
    return = ruledSrf( edge2d( x1, y1, x2, y1 ),
		       edge2d( x1, y2, x2, y2 ) );
sphereSrf = function( Radius ):s45:
    s45 = sin( 45 * pi / 180 ):
    return = surfRev( 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 ) ) ) *
	scale( vector( Radius, Radius, Radius ) );
torusSrf = function( MRadius, mRad ):
    return = SwpCircSrf( circle( vector( 0, 0, 0 ), MRadius ), mRad, 0 );
coneSrf = function( Height, Radius ):
    return = surfRev( ctlpt( E3, 0.0, 0.0, 0.0 ) +
		      ctlpt( E3, Radius, 0.0, 0.0 ) +
		      ctlpt( E3, 0.0, 0.0, Height ) );
cone2Srf = function( Height, Radius1, Radius2 ):
    return = surfRev( ctlpt( E3, 0.0, 0.0, 0.0 ) +
		      ctlpt( E3, Radius1, 0.0, 0.0 ) +
		      ctlpt( E3, Radius2, 0.0, Height ) +
		      ctlpt( E3, 0.0, 0.0, Height ) );
cylinSrf = function( Height, Radius ):
    return = surfRev( ctlpt( E3, 0.0, 0.0, 0.0 ) +
		      ctlpt( E3, Radius, 0.0, 0.0 ) +
		      ctlpt( E3, Radius, 0.0, Height ) +
		      ctlpt( E3, 0.0, 0.0, Height ) );
boxSrf = function( Width, Depth, Height ):
    return = list( coerce( planeSrf( 0, 0, Width, Depth ), e3 ),
		   coerce( planeSrf( 0, 0, Width, Depth ), e3 ) *
		       trans( vector( 0, 0, Height ) ),
		   coerce( planeSrf( 0, 0, Width, Height ), e3 ) *
		       rotx( 90 ),
		   coerce( planeSrf( 0, 0, Width, Height ), e3 ) *
		       rotx( 90 ) * trans( vector( 0, Depth, 0 ) ),
		   coerce( planeSrf( 0, 0, Height, Depth ), e3 ) *
		       roty( -90 ),
		   coerce( planeSrf( 0, 0, Height, Depth ), e3 ) *
		       roty( -90 ) * trans( vector( Width, 0, 0 ) ) );

#
# Transformations.
#
tx = function( r ):
    return = trans( vector( r, 0, 0 ) );
ty = function( r ):
    return = trans( vector( 0, r, 0 ) );
tz = function( r ):
    return = trans( vector( 0, 0, r ) );
sx = function( r ):
    return = scale( vector( r, 1, 1 ) );
sy = function( r ):
    return = scale( vector( 1, r, 1 ) );
sz = function( r ):
    return = scale( vector( 1, 1, r ) );
sc = function( s ):
    return = scale( vector( s, s, s ) );
rx = function( r ):
    return = rotx( r );
ry = function( r ):
    return = roty( r );
rz = function( r ):
    return = rotz( r );

RotVec2Z = function( w ):u:v: # Rotation Trans. of w dir to Z axis.
    if ( abs( coord( w, 0 ) ) > abs( coord( w, 1 ) ),
	 u = vector( 0, 1, 0 ),
	 u = vector( 1, 0, 0 ) ):
    w = normalize( w ):
    v = normalize( u ^ w ):
    u = normalize( w ^ v ):
    return =
	homomat( list( list( coord( u, 0 ), coord( v, 0 ), coord( w, 0 ), 0 ),
		       list( coord( u, 1 ), coord( v, 1 ), coord( w, 1 ), 0 ),
		       list( coord( u, 2 ), coord( v, 2 ), coord( w, 2 ), 0 ),
		       list( 0, 0, 0, 1 ) ) );
RotZ2Vec = function( w ): # Rotation Trans. of Z axis to W dir.
    return = RotVec2Z( w )^-1;

#
# Arrows.
#
arrow3d = function( Pt, Dir, Length, Width, HeadLength, HeadWidth ):
    return = list( cylinSrf( Length - HeadLength, Width / 2 )
			* RotZ2Vec( Dir )
			* trans( Pt ),
		   coneSrf( HeadLength, HeadWidth )
			* trans( vector( 0, 0, Length - HeadLength ) )
			* RotZ2Vec( Dir )
			* trans( Pt ) );

#
# Emulation of view, interact and other useful viewing functions using VIEWOBJ
# for the default display device.
#
beep = procedure():command_: # Make some noise.
    command_ = "BEEP":
    viewobj(command_);
viewclear = procedure():command_: # Clear the screen.
    command_ = "CLEAR":
    viewobj(command_);
viewdclear = procedure():command_: # Delayed clear screen.
    command_ = "DCLEAR":
    viewobj(command_);
viewdisc = procedure():command_: # Disconnect from this display device
    command_ = "DISCONNECT":
    viewobj(command_);
viewexit = procedure():command_: # Force Display device to exit.
    command_ = "EXIT":
    viewobj(command_);
viewmsave = procedure(name):command_: # Save viewing matrix under name.
    command_ = "MSAVE " + name:
    viewobj(command_);
viewremove = procedure(name):command_: # Remove an object from display.
    command_ = "REMOVE " + name:
    viewobj(command_);
viewanim = procedure(TMin, TMax, Dt):command_: # Animate a sequence.
    command_ = "ANIMATE " ^ TMin ^ " " ^ TMax ^ " " ^ Dt:
    viewobj(command_);
viewstate = procedure(state):command_: # Change state of display device.
    command_ = "STATE " + state:
    viewobj(command_);
view = procedure(none,clear): # Emulation of old VIEW command.
    if (clear != 0.0, viewdclear()):
    viewobj(none);
interact = procedure(none): # Emulation of old INTERACT command.
    viewdclear():
    viewobj(none):
    pause();

#
# Client communication helper functions.
#
clntclear = procedure(h):command_: # Clear the screen.
    command_ = "CLEAR":
    clntwrite(h, command_);
clntdclear = procedure(h):command_: # Delayed clear screen.
    command_ = "DCLEAR":
    clntwrite(h, command_);
clntdisc = procedure(h):command_: # Disconnect from this display device
    command_ = "DISCONNECT":
    clntwrite(h, command_);
clntexit = procedure(h):command_: # Force Display device to exit.
    command_ = "EXIT":
    clntwrite(h, command_);
clntgetobj = function(h, Name):command_: # Get an object from a client.
    command_ = "GETOBJ " + Name:
    clntwrite(h, command_):
    return = clntread(h, 100);
clntmsave = procedure(h, name):command_: # Save viewing matrix under name.
    command_ = "MSAVE " + name:
    clntwrite(h, command_);
clntremove = procedure(h, name):command_: # Remove an object from display.
    command_ = "REMOVE " + name:
    clntwrite(h, command_);
clntanim = procedure(h, TMin, TMax, Dt):command_: # Animate a sequence.
    command_ = "ANIMATE " ^ TMin ^ " " ^ TMax ^ " " ^ Dt:
    clntwrite(h, command_);
clntstate = procedure(h, state):command_: # Change state of display device.
    command_ = "STATE " + state:
    clntwrite(h, command_);
clntview = procedure(h, none, clear): # Emulation of old VIEW command.
    if (clear != 0.0, clntclear(h)):
    clntwrite(h, none);
cntrintr = procedure(none): # Emulation of old INTERACT command.
    clntclear(h):
    clntwrite(h, none):
    pause();

iritState("EchoSource", true);