# To unbundle, sh this file
echo Attic.cc 1>&2
cat >Attic.cc <<'End of Attic.cc'
#include "Attic.hxx"
#include <stdlib.h>
#include <stdio.h>

Attic::Attic( int m ) {
 halfMax = m/2;
 max = 2*halfMax;
 if( max <= 0 ) {
  printf( "*** ERROR ***                                     \n" );
  printf( "*** ERROR *** Attic::Attic()                      \n" );
  printf( "*** ERROR *** Negative mass anomaly!!             \n" );
  printf( "*** ERROR ***                                     \n" );
  exit(0);
 }
 upper = max - 1;
 stuff = new Junk* [ max ];
 count = 0;
}

Attic::~Attic() {
 delete [ max ] stuff;   // ??? Will not work well!!
}

void Attic::store( Junk* x ) {
 int i;

 if( count > upper ) {
  if( ++upper == max ) {
   upper = halfMax - 1;
   count = 0;
  }
  else {
   upper = max - 1;
  }
  /* DGN printf( "\n*** GARBAGE COLLECTION ***\n" ); */
  /* DGN printf( "    Deleting from %d to %d\n", count, upper ); */
  for( i = count; i <= upper; i++ ) {
   /* DGN printf( "%d  %x\n", i, stuff[i] ); */
   delete stuff[i];
  }
 }

 stuff[ count++ ] = x;
}
End of Attic.cc
echo dlist.cc 1>&2
cat >dlist.cc <<'End of dlist.cc'
/*************************************************************************
**************************************************************************
**************************************************************************
******                                                                ****
******  DLIST: Doubly linked lists. Based on slist, found in          ****
******         Stroustrup's C++ Programming Language                  ****
******                                                                ****
******  Copyright (c) 1990  Universities Research Association, Inc.   ****
******                All Rights Reserved                             ****
******                                                                ****
******  Author:    Leo Michelotti                                     ****
******                                                                ****
******             Fermilab                                           ****
******             P.O.Box 500                                        ****
******             Mail Stop 345                                      ****
******             Batavia, IL   60510                                ****
******                                                                ****
******             Phone: (708) 840 4956                              ****
******             Email: MICHELOTTI@FNALAD                           ****
******                    MICHELOTTI@ADCALC.FNAL.GOV                  ****
******                    ALMOND::MICHELOTTI                          ****
******                                                                ****
******  Release    Version 1.0 : May 20, 1990                         ****
******  Dates:                                                        ****
**************************************************************************
**************************************************************************
******                                                                ****
******  This material resulted from work developed under a            ****
******  Government Contract and is subject to the following           ****
******  license: The Government retains a paid-up, nonexclusive,      ****
******  irrevocable worldwide license to reproduce, prepare           ****
******  derivative works, perform publicly and display publicly       ****
******  by or for the Government, including the right to distribute   ****
******  to other Government contractors.  Neither the United          ****
******  States nor the United States Department of Energy, nor        ****
******  any of their employees, makes any warranty, express or        ****
******  implied, or assumes any legal liability or responsibility     ****
******  for the accuracy, completeness or usefulness of any           ****
******  information, apparatus, product, or process disclosed, or     ****
******  represents that its use would not infringe privately owned    ****
******  rights.                                                       ****
******                                                                ****
******  These files are made avaliable for use by specific            ****
******  individuals or a specific group.  They are not to be          ****
******  reproduced or redistributed.                                  ****
******                                                                ****
******  These files are provided "as is" with no warranties of any    ****
******  kind including the warranties of design,                      ****
******  merchantibility and fitness for a particular purpose,         ****
******  or arising from a course of dealing, usage or trade           ****
******  practice.                                                     ****
******                                                                ****
******  These files are provided with no support or obligation on     ****
******  the part of U.R.A. or Fermilab to assist in their use,        ****
******  correction, modification or enhancement.                      ****
******                                                                ****
******  Neither U.R.A. nor Fermilab shall be liable with              ****
******  respect to the infringement of copyrights, trade              ****
******  secrets or any patents by these files or any part             ****
******  thereof.                                                      ****
******                                                                ****
******  In no event will U.R.A., Fermilab, or the author(s) be        ****
******  liable for any lost revenue or profits or other               ****
******  special, indirect and consequential damages, even if          ****
******  they have been advised of the possibility of such             ****
******  damages.                                                      ****
******                                                                ****
**************************************************************************
******                                                                ****
******  Sep 18, 1992:                                                 ****
******  Copy constructor dlist::dlist( dlist& x ) installed.          ****
******                                                                ****
******  Oct 8, 1992:                                                  ****
******  dlist::putAbove and dlist::putBelow installed.                ****
******                                                                ****
******  Nov 16,1992                                                   ****
******  dlist::putBelow corrected to handle the "last" element.       ****
******                                                                ****
**************************************************************************
*************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include "dlist.hxx"

dlist::dlist( dlist& x )
{
ent w;
last = 0;
dlist_iterator getNext( x );
while( w = getNext() ) append( w );
}

void dlink::putAbove( dlink* a ) {
a -> prev = prev;
a -> next = this;
prev -> next = a;
prev = a;
}

dlist& dlist::operator=( dlist& x )
{
dlist_iterator getNext( x );
ent p;
clear();
while( p = getNext() ) append( p );
return *this;
}

void dlink::putBelow( dlink* a ) {
a -> prev = this;
a -> next = next;
next -> prev = a;
next = a;
}

ent dlist::remove( dlink* p ) {
ent a;

// NOTE:  This method assumes that *p is actually in the dlist;
//        it does not double-check to see that this is true.
//        This method purposely does not use dlist_iterator

// If p is the first element in the list ...
if ( p == last -> next ) {
  a = get();
  return a;
  }

// If p is the last element in the list ...
if ( p == last ) {
  // .. and there is only one element ..
  if( last->prev == last ) {
    a = get();
    return a;
    }
  // .. or there is more then one element ..
  last = p->prev;
  }

// If p is neither the first nor the last element
// in the list, or p is the last element in a list
// with more than one element ..
( p -> prev ) -> next  =  p -> next;
( p -> next ) -> prev  =  p -> prev;
a = p -> e;
delete p;
return a;

}


char dlist::remove( ent a )
{
dlist_traversor getNext( *this );
dlink* p;
ent w;
while( p = getNext() )  if( p->e == a ) {
                          w = remove( p );
                          return 1;
                          }
return 0;
} // End function char dlist::remove( ent a )


ent dlist::remove( int n )
{
dlist_traversor getNext( *this );
dlink* p;
ent w;
int i = 0;
while( p = getNext() )  if( i++ == n ) {
                          w = remove( p );
                          return w;
                          }
return 0;
} // End function char dlist::remove( ent a )


dlist dlist::remove( ent a, ent b )
{
dlist  ret;
dlink* p;
dlink* pr;
ent    w;
ent    wd;
ent    y;
dlist_traversor getNext ( *this );

do { p = getNext();
     w = p->e;
   } while ( p != 0 && w != a && w != b );

if( p == 0 ) {
  printf( "\n*** ERROR ***                              \n" );
  printf(   "*** ERROR *** dlist::remove( ent, ent )    \n" );
  printf(   "*** ERROR *** Reached the end of the list  \n" );
  printf(   "*** ERROR *** without finding markers.     \n" );
  printf(   "*** ERROR ***                              \n" );
  exit(1);
}

if( w == a ) y = b;
else  {      y = a;
             if( w != b ) {
               printf( "\n*** IMPOSSIBLE ERROR!\n" );
               exit(1);
             }
}

p = getNext();
w = ( pr = p )->e;
if( w == y || p == 0 ) {
  printf( "\n*** WARNING ***                                     \n" );
  printf(   "*** WARNING *** dlist::remove( ent, ent )           \n" );
  printf(   "*** WARNING *** No elements available for removal.  \n" );
  printf(   "*** WARNING ***                                     \n" );
  return *this;
}

do {
  p = getNext();
  wd = remove( pr );
  ret.append( w );
  w = ( pr = p )->e;
} while ( ( p != 0 ) && ( w != y ) );

if( p == 0 ) {
  printf( "\n*** WARNING ***                                \n" );
  printf(   "*** WARNING *** dlist::remove( ent, ent )      \n" );
  printf(   "*** WARNING *** Reached the end of the list.   \n" );
  printf(   "*** WARNING ***                                \n" );
}

return ret;
} // End function dlist dlist::remove( ent a, ent b )


dlist dlist::remove( int x, int y )     // ??? This is a kludged version. Rewrite!!
{
ent   a;
ent   b;
dlist c;
a = operator[]( x );
b = operator[]( y );
return remove( a, b );
}

void dlist::insert( ent a) {
if( last ) {
  last->next = new dlink( a, last, last->next );
  ((last->next)->next)->prev = last->next;
  }
else {
  last = new dlink( a, 0, 0 );
  last->next = last;
  last->prev = last;
  }
}

void dlist::append( ent a ) {
if( last ) {
  last = last->next = new dlink( a, last, last->next );
  (last->next)->prev = last;
  }
else {
  last = new dlink( a, 0, 0 );
  last->next = last;
  last->prev = last;
  }
}


char dlist::putAbove( ent a, ent b )    // Untested???
{
char   found;
dlink* w;
dlink* z;
dlist_traversor getNext( *this );
found = 0;
while( w = getNext() )
  if( w->e == a ) {
    found = 1;
    break;
    }
if( !found ) return 0;
z = new dlink( b, 0, 0 );
w->putAbove( z );
return 1;
}


char dlist::putBelow( ent a, ent b )    // Untested???
{
char   found;
dlink* w;
dlink* z;
dlist_traversor getNext( *this );
found = 0;
while( w = getNext() )
  if( w->e == a ) {
    found = 1;
    break;
    }
if( !found ) return 0;
z = new dlink( b, 0, 0 );
w->putBelow( z );
if( w == last ) last = z;
return 1;
}


ent dlist::get() {
if ( last == 0 ) return 0;   // In case the list is empty.
dlink* f = last->next;
ent r = f->e;
if ( f == last )
  last = 0;
else {
  last->next = f->next;
  (last->next)->prev = last;
  }
delete f;
return r;
}


ent dlist::operator[]( int n )
{
dlist_iterator getNext( *this );
int i = 0;
ent p = 0;
while( ( i++ <= n ) && ( p = getNext() ) ) ;
return p;
}


int dlist::size()
{
dlist_iterator getNext( *this );
int i = 0;
ent p;
while( p = getNext() ) i++;
return i;
}


void dlist::clear() {
dlink* l = last;
if( l == 0 ) return;
do {
  dlink* ll = l;
  l = l->next;
  delete ll;
  }
while( l != last );
last = 0;
}

void dlist::zap(){
dlink* l;
dlink* ll;
l = last;
if( l == 0 ) return;
do {
  ll = l;
  l = l->next;
  delete ll->e;
  delete ll;
  }
while( l != last );
delete l->e;
delete l;
last = 0;
}

ent dlist_iterator::operator()() {
ent ret = ce ? ( ce = ce->next )->e : 0;
if( ce == cs->last ) ce = 0;
return ret;
}

ent dlist_looper::operator()() {
ent ret = ce;
if( ret != 0 ) {
  ret = ( ce = ce->next )->e;
  if( ce == cs->last ) ce = 0;
  }
else {
  ce = cs->last;
  }
return ret;
}

dlink* dlist_traversor::operator()() {
dlink* ret;
if( ce != 0 ) ret = ( ce = ce->next );
else          ret = 0;
if( ce == cs->last ) ce = 0;
return ret;
}

ent dlist_reverseIterator::operator()() {
if( ce == 0 ) return 0;
ent ret = ce->e;
if( ( ce = ce->prev ) == cs->last ) ce = 0;
return ret;
}

/*
 *   Pool functions written by Michael Allen
 *   February, 1993
 */

#ifdef POOLED

_dlink_Shell* _dlink_Pool::top = 0;
_dlink_Pool::_dlink_Pool()
{
  register int i = _dlink_BLOCK_COUNT;
  _dlink_Shell* p = slot;
  _dlink_Shell* q;
  while (--i)
    {
      q = p++;
      q->next = p;
    }
  p->next = top;
  top = slot;
}

#endif   POOLED
End of dlist.cc
echo menu.cc 1>&2
cat >menu.cc <<'End of menu.cc'
/* -------------------------------------------------- 
		Required for use of menu tools....... */

#include <stdio.h>

#include "menu.hxx"

/* --------------------------------------------------- */

ttMenu::ttMenu() {
  MENUFUNC menuNullFunc = 0;
  for( int i = 0; i < MAXMENUCHOICES; i++ ) menuFunction[i] = menuNullFunc;
  promptMessage = "\0";
  numChoices = 0;
}


void ttMenu::go() {
  int menuChoice = 0;
  this->print();
  printf( "\n%s", promptMessage );
  scanf( "%d", &menuChoice );
  getchar();
  while( ( menuChoice >= 0 )&& ( menuChoice <= numChoices ) ) {
    if( menuChoice == 0 ) this->print();
    else (*menuFunction[menuChoice])();
    printf( "\n%s", promptMessage );
    scanf( "%d", &menuChoice );
    getchar();
    }
}


void ttMenu::print() {
int i = 0;
printf( "\n" );
while( ++i <= numChoices ) printf( "%2d) %s\n", i, menuMessage[i] );
printf( "\n" );
}


void ttMenu::setItem( char* message, MENUFUNC functionName ) {
  if( numChoices < MAXMENUCHOICES - 2 ) {
    menuMessage[ ++numChoices ] = message;
    menuFunction[ numChoices ] = functionName;
    }
}


void ttMenu::setOperation( int choice, MENUFUNC functionName ) {
  if( ( choice > 0 ) && ( choice <= numChoices ) )
    menuFunction[ choice ] = functionName;
}


void ttMenu::setMessage( int choice, char* message ) {
  if( ( choice > 0 ) && ( choice <= numChoices ) )
    menuMessage[ choice ] = message;
}


void ttMenu::setPrompt( char* message ) {
  promptMessage = message;
}


void menuNullFunc() {
return;
}
End of menu.cc
echo nexcom.cc 1>&2
cat >nexcom.cc <<'End of nexcom.cc'
char nexcom( int n, int k, int* r ) {
int t, h, i;
static char firstEntry = 1;
static char allDone = 0;

if( allDone ) {
  firstEntry = 1;
  allDone = 0;
  return 0;
  }

if( firstEntry ) {
  r[0] = n;
  if( k > 1 ) {
    for( i = 1; i < k; i++ ) r[i] = 0;
    firstEntry = 0;
    }
  else allDone = 1;
  return 1;
  }

h = 0;
while( r[h++] == 0 ) ;
h--;

t = r[h];
r[h] = 0;
r[0] = --t;
r[++h]++;

allDone = ( r[k-1] == n );

return 1;

}
End of nexcom.cc
echo slist.cc 1>&2
cat >slist.cc <<'End of slist.cc'
/* 
Implementation of the class slist, as detailed in Section 7.3.2
of The C++ Programming Language, by Stroustrup.
*/

#include <stdio.h>
#include <stdlib.h>
#include "slist.hxx"


void slist::insert( ent a) {
if( last )
  last->next = new slink( a, last->next );
else {
  last = new slink( a, 0 );
  last->next = last;
  }
}

void slist::append( ent a ) {
if( last )
  last = last->next = new slink( a, last->next );
else {
  last = new slink( a, 0 );
  last->next = last;
  }
}

ent slist::get() {
if ( last == 0 ) return 0;   // In case the list is empty.
// if ( last == 0 ) { 
//   printf( "ERROR:SLIST> get from empty slist." );
//   exit(0);
//   }
slink* f = last->next;
ent r = f->e;
if ( f == last )
  last = 0;
else
  last->next = f->next;
delete f;
return r;
}

slist& slist::operator=( slist& x )
{
slist_iterator getNext( x );
ent p;
clear();
while( p = getNext() ) append( p );
return *this;
}

ent slist::operator[]( int n )
{
slist_iterator getNext( *this );
int i = 0;
ent p = 0;
while( ( i++ <= n ) && ( p = getNext() ) ) ;
return p;
}

int slist::size()
{
slist_iterator getNext( *this );
int i = 0;
ent p;
while( p = getNext() ) i++;
return i;
}


void slist::clear() {
slink* l = last;
if( l == 0 ) return;
do {
  slink* ll = l;
  l = l->next;
  delete ll;
  } 
while( l != last );
last = 0;
}


ent slist_iterator::operator()() {
ent ret = ce ? ( ce = ce->next )->e : 0;
if( ce == cs->last ) ce = 0;
return ret;
}
End of slist.cc
echo binomcoeff_recursive.c 1>&2
cat >binomcoeff_recursive.c <<'End of binomcoeff_recursive.c'
int bcfRec( n, m ) 
int n;
int m;
{
if( ( n < 0 ) || ( m < 0 )  || ( n < m )  ) return 0;
if( ( n < 2 ) || ( m == 0 ) || ( m == n ) ) return 1;
return bcfRec( n-1, m ) + bcfRec( n-1, m-1);
}
End of binomcoeff_recursive.c
echo gcd.c 1>&2
cat >gcd.c <<'End of gcd.c'
/*
    Function returns the greatest common divisor of x and y.
*/
int gcd( x, y )
int x;
int y;
{
int larger, smaller;
int remainder;

if( x < 0 ) x = -x;
if( y < 0 ) y = -y;

if( x == 0 ) {
  if( y == 0 ) return 1;
  else return y;
  }

if( y == 0 ) return x;

if ( x < y ) {
  larger = y;
  smaller = x;
  }
else {
  larger = x;
  smaller = y;
  }

remainder = larger%smaller;

while( remainder != 0 ) {
  larger = smaller;
  smaller = remainder;
  remainder = larger%smaller;
  }

return smaller;

}

End of gcd.c
echo mod2pi.c 1>&2
cat >mod2pi.c <<'End of mod2pi.c'
#define twopi 6.2831853072

double mod2pi( x )
double x;
{
if ( x == 0.0 ) return x;
while ( x < 0.0 ) x+= twopi;
while ( x >= twopi ) x -= twopi;
return x;
}
End of mod2pi.c
echo nearint.c 1>&2
cat >nearint.c <<'End of nearint.c'
#include <math.h>
#include <stdio.h>

int nearestInteger( x )
double x;
{
double lo, hi;

lo = floor(x);
hi = ceil(x);
if( x - lo >= hi - x )
  return (int) hi;
else
  return (int) lo;
}
End of nearint.c
echo parsef.c 1>&2
cat >parsef.c <<'End of parsef.c'
#include <stdio.h>

#define term_char_size 80   /* Number of characters in a terminal line.    */
#define maxchar 1300        /* Max number of characters in a file line.    */
#define maxtokens 50        /* Max number of tokens in a parse list.        */
#define max_token_size 25   /* Max number of characters in a token. */

typedef char token[max_token_size];
struct array_of_tokens
 {
    token piece[maxtokens];
 };

struct array_of_tokens parsef( line, numitems )
char *line;
int  *numitems;
{
int counter, place, length;
struct array_of_tokens item;

  place = 0;
  counter = 0;
  length = strspn( &line[place], " \t" );
  place = place + length;
  while( sscanf( &line[place], "%s", &(item.piece[counter]) ) != EOF 
         && counter < maxtokens )
  {
    length = strlen( &(item.piece[counter]) );
    place = place + length;
    length = strspn( &line[place], " \t" );
    place = place + length;
    counter = counter + 1;
  };  

*numitems = counter;
return item;

}
End of parsef.c
echo strequals.c 1>&2
cat >strequals.c <<'End of strequals.c'
#include <stdio.h>

strEquals( str_1, str_2 )
char *str_1, *str_2;
{
int answer;

answer = ( strlen(str_1) == strlen(str_2) ) && 
         ( strcmp( str_1, str_2 ) == 0 );

return answer;

}
End of strequals.c
echo threeVector.cc 1>&2
cat >threeVector.cc <<'End of threeVector.cc'
#include "threeVector.hxx"

extern threeVector THRV_buff;

threeVector::threeVector()
{
for ( int i = 0; i < 3; i++ ) comp[i] = 0.0;
}

threeVector::threeVector( double x, double y, double z )
{
comp[0] = x;
comp[1] = y;
comp[2] = z;
}

threeVector::threeVector( double* x )
{
for ( int i = 0; i < 3; i++ ) comp[i] = x[i];
}

threeVector::threeVector( threeVector& x )
// ??? Put in garbage collection!!!
{
for ( int i = 0; i < 3; i++ ) comp[i] = x.comp[i];
}

threeVector::~threeVector()
{
}

void threeVector::set( double x, double y, double z )
{
comp[0] = x;
comp[1] = y;
comp[2] = z;
}

void threeVector::set( double* x )
{
for ( int i = 0; i < 3; i++ ) comp[i] = x[i];
}

double threeVector::operator* ( threeVector& x )
{
double u;
int i;
u = 0.0;
for ( i = 0; i < 3; i++ ) u += comp[i] * x.comp[i];
return u;
}

threeVector threeVector::operator^ ( threeVector& x )
{
THRV_buff.comp[ 0 ] = comp[ 1 ] * x.comp[ 2 ] - comp[ 2 ] * x.comp[ 1 ];
THRV_buff.comp[ 1 ] = comp[ 2 ] * x.comp[ 0 ] - comp[ 0 ] * x.comp[ 2 ];
THRV_buff.comp[ 2 ] = comp[ 0 ] * x.comp[ 1 ] - comp[ 1 ] * x.comp[ 0 ];
return THRV_buff;
}

threeVector threeVector::operator= ( threeVector& x )
{
for ( int i = 0; i < 3; i++ ) comp[i] = x.comp[i];
return *this;
}

char threeVector::operator== ( threeVector& x )
{
return ( comp[0] == x.comp[0] ) && 
       ( comp[1] == x.comp[1] ) && 
       ( comp[2] == x.comp[2] );
}

char threeVector::operator!= ( threeVector& x )
{
return !( operator==( x ) );
}

threeVector threeVector::operator+ ( threeVector& x )
{
for ( int i = 0; i < 3; i++ ) 
  THRV_buff.comp[i] = comp[i] + x.comp[i];
return THRV_buff;
}

void threeVector::operator+= ( threeVector& x )
{
for ( int i = 0; i < 3; i++ ) comp[i] += x.comp[i];
}

threeVector threeVector::operator- ( threeVector& x )
{
for ( int i = 0; i < 3; i++ ) 
  THRV_buff.comp[i] = comp[i] - x.comp[i];
return THRV_buff;
}

threeVector operator- ( threeVector& x )
{
for ( int i = 0; i < 3; i++ ) 
  THRV_buff.comp[i] = - x.comp[i];
return THRV_buff;
}

double threeVector::operator[] ( int i )
{
if ( ( 0 <= i ) && ( i <= 2 ) ) return comp[i];
else                            return 0.0;
}

double& threeVector::operator() ( int i )
{
static double dummy;
dummy = 0.0;
if ( ( 0 <= i ) && ( i <= 2 ) ) return comp[i];
else                            return dummy;
}

double threeVector::norm ()
{
double x;
x = 0.0;
for ( int i = 0; i < 3; i++ ) x += comp[i]*comp[i];
return sqrt(x);
}

threeVector threeVector::unit ()
{
double x;
int    i;
x = 0.0;
for ( i = 0; i < 3; i++ ) x += comp[i]*comp[i];
x = sqrt(x);
for ( i = 0; i < 3; i++ ) THRV_buff.comp[i] = comp[i] / x;
return THRV_buff;
}

void threeVector::rotate ( threeVector& v, double theta )
{
double c, s;
threeVector e, u;

e = unit();
c = cos( theta );
s = sin( theta );
u = ( c*v ) +
    ( s*( e^v) ) +
    ( ( ( 1.0 - c )*(e*v) )*e );
for ( int i = 0; i < 3; i++ ) v.comp[i] = u.comp[i];
}

threeVector operator* ( double c, threeVector& v )
{
threeVector u;
for ( int i = 0; i < 3; i++ ) u.comp[i] = c * v.comp[i];
return u;
}

threeVector operator* ( threeVector& v, double c )
{
threeVector u;
for ( int i = 0; i < 3; i++ ) u.comp[i] = c * v.comp[i];
return u;
}

threeVector threeVector::operator/ ( double x ) {
 threeVector u;
 for( int i = 0; i < 3; i++ ) u.comp[i] = comp[i] / x;
 return u;
}

threeVector threeVector::operator*=( double x ) {
 for( int i = 0; i < 3; i++ ) comp[i] *= x;
 return *this;
}

threeVector threeVector::operator/=( double x ) {
 for( int i = 0; i < 3; i++ ) comp[i] /= x;
 return *this;
}

End of threeVector.cc
echo twoVector.cc 1>&2
cat >twoVector.cc <<'End of twoVector.cc'
#include "twoVector.hxx"

extern twoVector TWOV_buff;

twoVector::twoVector()
{
for ( int i = 0; i < 2; i++ ) comp[i] = 0.0;
}

twoVector::twoVector( double x, double y )
{
comp[0] = x;
comp[1] = y;
}

twoVector::twoVector( double* x )
{
for ( int i = 0; i < 2; i++ ) comp[i] = x[i];
}

twoVector::twoVector( twoVector& x )
{
for ( int i = 0; i < 2; i++ ) comp[i] = x.comp[i];
}

twoVector::~twoVector()
{
}

void twoVector::set( double x, double y )
{
comp[0] = x;
comp[1] = y;
}

void twoVector::set( double* x )
{
for ( int i = 0; i < 2; i++ ) comp[i] = x[i];
}

double twoVector::operator* ( twoVector& x )
{
double u;
int i;
u = 0.0;
for ( i = 0; i < 2; i++ ) u += comp[i] * x.comp[i];
return u;
}

twoVector twoVector::operator= ( twoVector& x )
{
for ( int i = 0; i < 2; i++ ) comp[i] = x.comp[i];
return *this;
}

char twoVector::operator== ( twoVector& x )
{
return ( comp[0] == x.comp[0] ) && 
       ( comp[1] == x.comp[1] );
}

char twoVector::operator!= ( twoVector& x )
{
return !( operator==( x ) );
}

twoVector twoVector::operator+ ( twoVector& x )
{
for ( int i = 0; i < 2; i++ ) 
  TWOV_buff.comp[i] = comp[i] + x.comp[i];
return TWOV_buff;
}

void twoVector::operator+= ( twoVector& x )
{
for ( int i = 0; i < 2; i++ ) comp[i] += x.comp[i];
}

twoVector twoVector::operator- ( twoVector& x )
{
for ( int i = 0; i < 2; i++ ) 
  TWOV_buff.comp[i] = comp[i] - x.comp[i];
return TWOV_buff;
}

twoVector operator- ( twoVector& x )
{
for ( int i = 0; i < 2; i++ ) 
  TWOV_buff.comp[i] = - x.comp[i];
return TWOV_buff;
}

double twoVector::operator[] ( int i )
{
if ( ( 0 <= i ) && ( i <= 1 ) ) return comp[i];
else                            return 0.0;
}

double& twoVector::operator() ( int i )
{
static double dummy;
dummy = 0.0;
if ( ( 0 <= i ) && ( i <= 1 ) ) return comp[i];
else                            return dummy;
}

double twoVector::norm ()
{
double x;
x = 0.0;
for ( int i = 0; i < 2; i++ ) x += comp[i]*comp[i];
return sqrt(x);
}

twoVector twoVector::unit ()
{
double x;
int    i;
x = 0.0;
for ( i = 0; i < 2; i++ ) x += comp[i]*comp[i];
x = sqrt(x);
for ( i = 0; i < 2; i++ ) TWOV_buff.comp[i] = comp[i] / x;
return TWOV_buff;
}

twoVector twoVector::rotate ( double theta )
{
double c, s;
c = cos( theta );
s = sin( theta );
TWOV_buff.set( c*comp[0] - s*comp[1], 
               s*comp[0] + c*comp[1] );
return TWOV_buff;
}

twoVector operator* ( double c, twoVector& v )
{
twoVector u;
for ( int i = 0; i < 2; i++ ) u.comp[i] = c * v.comp[i];
return u;
}

twoVector operator* ( twoVector& v, double c )
{
twoVector u;
for ( int i = 0; i < 2; i++ ) u.comp[i] = c * v.comp[i];
return u;
}

twoVector twoVector::operator/ ( double x ) {
 twoVector u;
 for( int i = 0; i < 3; i++ ) u.comp[i] = comp[i] / x;
 return u;
}

twoVector twoVector::operator*=( double x ) {
 for( int i = 0; i < 3; i++ ) comp[i] *= x;
 return *this;
}

twoVector twoVector::operator/=( double x ) {
 for( int i = 0; i < 3; i++ ) comp[i] /= x;
 return *this;
}

End of twoVector.cc
