Initial commit

6 files changed, 2785 insertions, 0 deletions

diff --git a/tksao/vector/vector.C b/tksao/vector/vector.C
new file mode 100644
index 0000000..b7c3340
--- /dev/null
+++ b/tksao/vector/vector.C
@@ -0,0 +1,413 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include <tk.h>
+
+#include "vector.h"
+#include "vector3d.h"
+#include "fuzzy.h"
+
+// Vector::manip
+
+int Vector::separator = ios_base::xalloc();
+int Vector::unit = ios_base::xalloc();
+
+Vector::Vector(const Vector3d& a)
+{
+  v[0]=a.v[0];
+  v[1]=a.v[1];
+  v[2]=1;
+}
+
+Vector& Vector::operator=(const Vector3d& a)
+{
+  v[0]=a.v[0];
+  v[1]=a.v[1];
+  v[2]=1;
+  return *this;
+}
+
+Vector& Vector::clip(const BBox& bb)
+{
+  // restrict vector by bbox
+  Vector ll(bb.ll);
+  Vector ur(bb.ur);
+  if (v[0]<ll[0])
+    v[0]=ll[0];
+  if (v[0]>ur[0])
+    v[0]=ur[0];
+  if (v[1]<ll[1])
+    v[1]=ll[1];
+  if (v[1]>ur[1])
+    v[1]=ur[1];
+  return *this;
+}
+
+Vector Vector::TkCanvasPs(void* canvas)
+{
+  return Vector(v[0], Tk_CanvasPsY((Tk_Canvas)canvas, v[1]));
+}
+
+ostream& operator<<(ostream& os, const Vector& v)
+{
+  unsigned char sep = (unsigned char)os.iword(Vector::separator);
+  if (!sep)
+    sep = ' ';
+
+  unsigned char unit = (unsigned char)os.iword(Vector::unit);
+  if (!unit)
+    os << v.v[0] << sep << v.v[1];
+  else
+    os << v.v[0] << unit << sep << v.v[1] << unit;
+
+  // reset unit
+  os.iword(Vector::unit) = '\0';
+
+  return os;
+}
+
+istream& operator>>(istream& s, Vector& v)
+{
+  s >> v.v[0] >> v.v[1];
+  return s;
+}
+
+// Vertex
+
+ostream& operator<<(ostream& os, const Vertex& v)
+{
+  os << v.vector;
+  return os;
+}
+
+// Matrix
+
+Matrix& Matrix::operator*=(const Matrix& a)
+{
+  Matrix r;
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      r.m[i][j] = 
+	m[i][0]*a.m[0][j] + 
+	m[i][1]*a.m[1][j] + 
+	m[i][2]*a.m[2][j];
+
+  return *this=r;
+}
+
+Matrix Matrix::invert()
+{
+  Matrix cc = this->cofactor();
+  Matrix aa = cc.adjoint();
+  double dd = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + m[0][2]*aa.m[2][0];
+
+  Matrix rr;
+  for (int ii=0; ii<3; ii++ )
+    for (int jj=0; jj<3; jj++)
+      rr.m[ii][jj] = aa.m[ii][jj]/dd;
+
+  return rr;
+}
+
+Matrix Matrix::cofactor()
+{
+  Matrix rr;
+
+  rr.m[0][0] = +(m[1][1]*m[2][2]-m[1][2]*m[2][1]);
+  rr.m[0][1] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]);
+  rr.m[0][2] = +(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
+
+  rr.m[1][0] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]);
+  rr.m[1][1] = +(m[0][0]*m[2][2]-m[0][2]*m[2][0]);
+  rr.m[1][2] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]);
+
+  rr.m[2][0] = +(m[0][1]*m[1][2]-m[0][2]*m[1][1]);
+  rr.m[2][1] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]);
+  rr.m[2][2] = +(m[0][0]*m[1][1]-m[0][1]*m[1][0]);
+
+  return rr;
+}
+
+double Matrix::det()
+{
+  return 
+    + m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])
+    - m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])
+    + m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
+}
+
+Matrix Matrix::adjoint()
+{
+  Matrix rr;
+  for (int ii=0; ii<3; ii++)
+    for (int jj=0; jj<3; jj++)
+      rr.m[jj][ii] = m[ii][jj];
+
+  return rr;
+}
+
+ostream& operator<<(ostream& os, const Matrix& m)
+{
+  os << ' ';
+  for (int i=0; i<3; i++)
+    for (int j=0; j<2; j++)
+      os << m.m[i][j] << ' ';
+
+  return os;
+}
+
+istream& operator>>(istream& s, Matrix& m)
+{
+  for (int i=0; i<3; i++ )
+    for (int j=0; j<2; j++)
+      s >> m.m[i][j];
+
+  return s;
+}
+
+// Translate
+
+ostream& operator<<(ostream& os, const Translate& m)
+{
+  os << ' ' << m.m[2][0] << ' ' << m.m[2][1] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, Translate& m)
+{
+  s >> m.m[2][0] >> m.m[2][1];
+  return s;
+}
+
+// Scale
+
+ostream& operator<<(ostream& os, const Scale& m)
+{
+  os << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, Scale& m)
+{
+  s >> m.m[0][0] >> m.m[1][1];
+  return s;
+}
+
+// Rotate
+
+Rotate::Rotate(double a) : Matrix()
+{
+  // note: signs reverse for X-Windows (origin is upper left)
+  m[0][0] = cos(a);
+  m[0][1] = -sin(a);
+  m[1][0] = sin(a);
+  m[1][1] = cos(a);
+
+  // this fixes a problem with numbers too small and tring to invert the matrix
+  tzero(&m[0][0]);
+  tzero(&m[0][1]);
+  tzero(&m[1][0]);
+  tzero(&m[1][1]);
+}
+
+ostream& operator<<(ostream& os, const Rotate& m)
+{
+  os << ' ' << m.m[0][0] << ' ' << m.m[0][1]
+    << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, Rotate& m)
+{
+  s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1];
+  return s;
+}
+
+// BBox
+
+BBox::BBox(double a, double b, double c, double d)
+{
+  // we want a 'positive' box
+  ll.v[0] = a < c ? a : c;
+  ll.v[1] = b < d ? b : d;
+  ur.v[0] = a < c ? c : a;
+  ur.v[1] = b < d ? d : b;
+}
+
+BBox::BBox(const Vector& l, const Vector& h)
+{
+  // we want a 'positive' box
+  ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0];
+  ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1];
+  ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0];
+  ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1];
+}
+
+int BBox::isIn(const Vector& v) const
+{
+  return !(v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || 
+	   v.v[0] > ur.v[0] || v.v[1] > ur.v[1]);
+}
+
+int BBox::isIn(const BBox& bb) const
+{
+  // return 0 if outside, > 0 if intersection
+  // = 4 if inside
+  BBox b = bb;
+  return isIn(b.ll) + isIn(b.ur) + isIn(b.ul()) + isIn(b.lr());
+}
+
+BBox& BBox::bound(const Vector& v)
+{
+  if (v.v[0] < ll[0])
+    ll[0] = v.v[0];
+  if (v.v[1] < ll[1])
+    ll[1] = v.v[1];
+    
+  if (v.v[0] > ur[0])
+    ur[0] = v.v[0];
+  if (v.v[1] > ur[1])
+    ur[1] = v.v[1];
+
+  return *this;
+}
+
+BBox& BBox::bound(BBox b)
+{
+  this->bound(b.ll);
+  this->bound(b.lr());
+  this->bound(b.ur);
+  this->bound(b.ul());
+
+  return *this;
+}
+
+BBox intersect(const BBox& a, const BBox& b)
+{
+  // test for obvious
+  int ab = a.isIn(b);
+  int ba = b.isIn(a);
+
+  // no intersection?
+  if (ab==0 && ba == 0) {
+    // maybe they are just crossed, check the centers
+    int abc = a.isIn(((BBox&)b).center());
+    int bac = b.isIn(((BBox&)a).center());
+    if (abc==0 && bac==0)
+      return BBox();
+  }
+
+  if (ab == 4) // b is inside a
+    return b; 
+  if (ba == 4) // a is inside b
+    return a;
+
+  // else, there seems to be some overlap
+  BBox r;
+  r.ll.v[0] = (a.ll.v[0] > b.ll.v[0]) ? a.ll.v[0] : b.ll.v[0];
+  r.ll.v[1] = (a.ll.v[1] > b.ll.v[1]) ? a.ll.v[1] : b.ll.v[1];
+  r.ur.v[0] = (a.ur.v[0] < b.ur.v[0]) ? a.ur.v[0] : b.ur.v[0];
+  r.ur.v[1] = (a.ur.v[1] < b.ur.v[1]) ? a.ur.v[1] : b.ur.v[1];
+
+  return r;
+}
+
+ostream& operator<<(ostream& os, const BBox& b)
+{
+  os << b.ll << ' ' << b.ur;
+  return os;
+}
+
+// Cohen–Sutherland clipping algorithm
+// bounded diagonally by (0,0), and (width, height)
+
+const int INSIDE = 0; // 0000
+const int LEFT = 1;   // 0001
+const int RIGHT = 2;  // 0010
+const int BOTTOM = 4; // 0100
+const int TOP = 8;    // 1000
+
+static int calcOutCode(Vector vv, int width, int height)
+{
+  int code = INSIDE;
+
+  if (vv[0] < 0)           // to the left of clip window
+    code |= LEFT;
+  else if (vv[0] > width)  // to the right of clip window
+    code |= RIGHT;
+  if (vv[1] < 0)           // below the clip window
+    code |= BOTTOM;
+  else if (vv[1] > height) // above the clip window
+    code |= TOP;
+
+  return code;
+}
+
+int clip(Vector* v0, Vector* v1, int width, int height) 
+{
+  int outcode0 = calcOutCode(*v0, width, height);
+  int outcode1 = calcOutCode(*v1, width, height);
+  int accept = false;
+
+  while (true) {
+    if (!(outcode0 | outcode1)) { 
+      // Bitwise OR is 0. Trivially accept and get out of loop
+      accept = true;
+      break;
+    } 
+    else if (outcode0 & outcode1) {
+      // Bitwise AND is not 0. Trivially reject and get out of loop
+      break;
+    } 
+    else {
+      // At least one endpoint is outside the clip rectangle; pick it.
+      int outcodeOut = outcode0 ? outcode0 : outcode1;
+
+      // Now find the intersection point
+      // y = y0 + slope * (x - x0)
+      // x = x0 + (1 / slope) * (y - y0)
+      double x, y;
+      if (outcodeOut & TOP) {
+	// point is above the clip rectangle
+	x = (*v0)[0] + ((*v1)[0]-(*v0)[0]) * 
+	  (height-(*v0)[1]) / ((*v1)[1]-(*v0)[1]);
+	y = height;
+      } 
+      else if (outcodeOut & BOTTOM) {
+	// point is below the clip rectangle
+	x = (*v0)[0] + ((*v1)[0]-(*v0)[0]) * 
+	  (0-(*v0)[1]) / ((*v1)[1]-(*v0)[1]);
+	y = 0;
+      } 
+      else if (outcodeOut & RIGHT) {
+	// point is to the right of clip rectangle
+	y = (*v0)[1] + ((*v1)[1]-(*v0)[1]) * 
+	  (width-(*v0)[0]) / ((*v1)[0]-(*v0)[0]);
+	x = width;
+      } 
+      else if (outcodeOut & LEFT) {
+	// point is to the left of clip rectangle
+	y = (*v0)[1] + ((*v1)[1]-(*v0)[1]) * 
+	  (0-(*v0)[0]) / ((*v1)[0]-(*v0)[0]);
+	x = 0;
+      }
+
+      // Now we move outside point to intersection point to clip
+      // and get ready for next pass.
+      if (outcodeOut == outcode0) {
+	(*v0)[0] = x;
+	(*v0)[1] = y;
+	outcode0 = calcOutCode(*v0, width, height);
+      } 
+      else {
+	(*v1)[0] = x;
+	(*v1)[1] = y;
+	outcode1 = calcOutCode(*v1, width, height);
+      }
+    }
+  }
+
+  return accept;
+}
+
diff --git a/tksao/vector/vector.h b/tksao/vector/vector.h
new file mode 100644
index 0000000..da0932b
--- /dev/null
+++ b/tksao/vector/vector.h
@@ -0,0 +1,372 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#ifndef __vector_h__
+#define __vector_h__
+
+#include <math.h>
+#include <float.h>
+
+#include <iostream>
+using namespace std;
+
+class Vector3d;
+class Matrix;
+class BBox;
+
+class Vector {
+ public:
+  static int separator;
+  static int unit;
+
+  double v[3];
+
+ public:
+  Vector() 
+    {v[0]=0; v[1]=0; v[2]=1;}
+  Vector(double* f) 
+    {v[0]=f[0]; v[1]=f[1]; v[2]=1;}
+  Vector(double x, double y) 
+    {v[0]=x; v[1]=y; v[2]=1;}
+
+  Vector(const Vector& a) 
+    {v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];}
+  Vector& operator=(const Vector& a)
+    {v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2]; return *this;}
+
+  Vector(const Vector3d&);
+  Vector& operator=(const Vector3d&);
+
+  double& operator[](int i) 
+  {return v[i];} // return element
+  double* vv() 
+  {return v;} // return vector
+
+  Vector& origin()
+    {v[0]=0; v[1]=0; v[2]=1; return *this;}
+
+  Vector& operator+=(const Vector& a) // addition
+    {v[0]+=a.v[0]; v[1]+=a.v[1]; return *this;}
+  Vector& operator-=(const Vector& a) // subtraction
+    {v[0]-=a.v[0]; v[1]-=a.v[1]; return *this;}
+  Vector& operator*=(double f) // scalar multipy
+    {v[0]*=f; v[1]*=f; return *this;}
+  Vector& operator/=(double f) // scalar division
+    {v[0]/=f; v[1]/=f; return *this;}
+  Vector& operator*=(const Matrix& m); // vector multipy
+
+  Vector abs() 
+  {return Vector(fabs(v[0]),fabs(v[1]));}
+  double angle() 
+  {return atan2(v[1],v[0]);}
+  Vector ceil() 
+  {return Vector(::ceil(v[0]),::ceil(v[1]));}
+  Vector floor() 
+  {return Vector(::floor(v[0]),::floor(v[1]));}
+  Vector invert() 
+  {return Vector(1/v[0],1/v[1]);}
+  double length() 
+  {return sqrt(v[0]*v[0]+v[1]*v[1]);}
+  double area()
+  {return v[0]*v[1];}
+  Vector round() 
+  {return Vector((int)(v[0]+.5),(int)(v[1]+.5));}
+  Vector normalize() 
+  {double d = sqrt(v[0]*v[0]+v[1]*v[1]);
+    return d ? Vector(v[0]/d,v[1]/d) : Vector();}
+  // restrict vector by bbox
+  Vector& clip(const BBox&);
+
+  Vector TkCanvasPs(void* canvas);
+};
+
+// This I/O manipulator flips the number stored in iword between 0 and 1
+/*
+inline ios_base& setfoo(ios_base& os)
+{
+  os.iword(Vector::foo) = !os.iword(Vector::foo);
+  return os;
+}
+*/
+
+// Vector separator(int)
+
+struct _Setseparator {int _M_n;};
+inline _Setseparator setseparator(int __n)
+{
+  _Setseparator __x;
+  __x._M_n = __n;
+  return __x;
+}
+
+template<class _CharT, class _Traits>
+  inline basic_ostream<_CharT,_Traits>& 
+  operator<<(basic_ostream<_CharT,_Traits>& __os, _Setseparator __f)
+{ 
+  __os.iword(Vector::separator) = __f._M_n;
+  return __os; 
+}
+
+// Vector unit(int)
+
+struct _Setunit {int _M_n;};
+inline _Setunit setunit(int __n)
+{
+  _Setunit __x;
+  __x._M_n = __n;
+  return __x;
+}
+
+template<class _CharT, class _Traits>
+  inline basic_ostream<_CharT,_Traits>& 
+  operator<<(basic_ostream<_CharT,_Traits>& __os, _Setunit __f)
+{ 
+  __os.iword(Vector::unit) = __f._M_n;
+  return __os; 
+}
+
+ostream& operator<<(ostream&, const Vector&);
+istream& operator>>(istream&, Vector&);
+
+inline Vector operator-(const Vector& a)
+{return Vector(-a.v[0],-a.v[1]);}
+inline Vector operator+(const Vector& a, const Vector& b) 
+{return Vector(a) +=b;}
+inline Vector operator-(const Vector& a, const Vector& b) 
+{return Vector(a) -=b;}
+inline Vector operator*(const Vector& a, double b)
+{return Vector(a) *=b;}
+inline Vector operator/(const Vector& a, double b)
+{return Vector(a) /=b;}
+inline Vector operator*(const Vector& v, const Matrix& m) 
+{return Vector(v) *=m;}
+inline double operator*(const Vector& a, const Vector& b) // dot product
+{double r =0; r+=a.v[0]*b.v[0]; r+=a.v[1]*b.v[1]; return r;}
+
+class Vertex {
+public:
+  Vector vector;
+
+private:
+  Vertex* next_;
+  Vertex* previous_;
+
+public:
+  Vertex() 
+    {next_=NULL; previous_=NULL;}
+  Vertex(double x, double y) 
+    {vector=Vector(x,y); next_=NULL; previous_=NULL;}
+  Vertex(const Vector& a) 
+    {vector=a; next_=NULL; previous_=NULL;}
+
+  Vertex(const Vertex& a) 
+    {vector=a.vector; next_=a.next_; previous_=a.previous_;}
+  Vertex& operator=(const Vertex& a)
+    {vector=a.vector; next_=a.next_; previous_=a.previous_; return *this;}
+  
+  Vertex* next() 
+  {return next_;}
+  Vertex* previous()
+  {return previous_;}
+  void setNext(Vertex* v)
+  {next_=v;}
+  void setPrevious(Vertex* v)
+  {previous_=v;}
+};
+ostream& operator<<(ostream&, const Vertex&);
+
+class Matrix {
+ public:
+  double m[3][3];
+
+ public:
+  Matrix()
+    { m[0][0]=1; m[0][1]=0; m[0][2]=0; 
+      m[1][0]=0; m[1][1]=1; m[1][2]=0; 
+      m[2][0]=0; m[2][1]=0; m[2][2]=1;}
+  Matrix(double a, double b, 
+	 double c, double d, 
+	 double e, double f)
+    { m[0][0]=a; m[0][1]=b; m[0][2]=0; 
+      m[1][0]=c; m[1][1]=d; m[1][2]=0; 
+      m[2][0]=e; m[2][1]=f; m[2][2]=1;}
+  Matrix(double a, double b, double c, 
+	 double d, double e, double f, 
+	 double g, double h, double i)
+    { m[0][0]=a; m[0][1]=b; m[0][2]=c; 
+      m[1][0]=d; m[1][1]=e; m[1][2]=f; 
+      m[2][0]=g; m[2][1]=h; m[2][2]=i;}
+
+  Matrix(const Matrix& a)
+    { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; m[0][2]=a.m[0][2]; 
+      m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1]; m[1][2]=a.m[1][2]; 
+      m[2][0]=a.m[2][0]; m[2][1]=a.m[2][1]; m[2][2]=a.m[2][2];}
+  Matrix& operator=(const Matrix& a)
+    { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; m[0][2]=a.m[0][2]; 
+      m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1]; m[1][2]=a.m[1][2]; 
+      m[2][0]=a.m[2][0]; m[2][1]=a.m[2][1]; m[2][2]=a.m[2][2]; return *this;}
+  
+  double matrix(int i, int j) // return element
+  {return m[i][j];}
+  Vector operator[](int i) // return row
+  {return Vector(m[i]);}
+  double* mm() const // return matrix
+  {return (double*)m;}
+
+  Matrix& identity()
+    { m[0][0]=1; m[0][1]=0; m[0][2]=0; 
+      m[1][0]=0; m[1][1]=1; m[1][2]=0; 
+      m[2][0]=0; m[2][1]=0; m[2][2]=1; return *this;}
+  Matrix& operator*=(const Matrix&); // matrix multiply
+
+  Matrix invert();
+  Matrix cofactor();
+  Matrix adjoint();
+  double det();
+};
+ostream& operator<<(ostream&, const Matrix&);
+istream& operator>>(istream&, Matrix&);
+
+inline Matrix operator*(const Matrix& a, const Matrix& b) 
+{return Matrix(a) *= b;}
+inline Vector& Vector::operator*=(const Matrix& m)
+{
+  double vv[3];
+  double* mm = (double*)(m.m);
+  vv[0] = v[0]*mm[0] + v[1]*mm[3] + v[2]*mm[6]; 
+  vv[1] = v[0]*mm[1] + v[1]*mm[4] + v[2]*mm[7]; 
+  vv[2] = v[0]*mm[2] + v[1]*mm[5] + v[2]*mm[8];
+  v[0] = vv[0];
+  v[1] = vv[1];
+  v[2] = vv[2];
+  return *this;
+}
+
+class Translate : public Matrix {
+public:
+  Translate() 
+    {};
+  Translate(double x, double y) 
+    {m[2][0]=x; m[2][1]=y;}
+  Translate(const Vector& v) 
+    {m[2][0]=v.v[0]; m[2][1]=v.v[1];}
+  Translate(const Matrix& a) 
+    {m[2][0] = a.m[2][0]; m[2][1] = a.m[2][1];}
+};
+ostream& operator<<(ostream&, const Translate&);
+istream& operator>>(istream&, Translate&);
+
+class Scale : public Matrix {
+public:
+  Scale() 
+    {};
+  Scale(double a) 
+    {m[0][0]=a; m[1][1]=a;}
+  Scale(double a, double b) 
+    {m[0][0]=a; m[1][1]=b;}
+  Scale(const Vector& v) 
+    {m[0][0]=v.v[0]; m[1][1]=v.v[1];}
+  Scale(const Matrix& a) 
+    {m[0][0] = a.m[0][0]; m[1][1] = a.m[1][1];}
+};
+ostream& operator<<(ostream&, const Scale&);
+istream& operator>>(istream&, Scale&);
+
+class FlipX : public Matrix {
+public:
+  FlipX() 
+    {m[0][0] = -1;}
+};
+
+class FlipY : public Matrix {
+public:
+  FlipY() 
+    {m[1][1] = -1;}
+};
+
+class FlipXY : public Matrix {
+public:
+  FlipXY() 
+    {m[0][0] = -1; m[1][1] = -1;}
+};
+
+class Rotate : public Matrix {
+public:
+  Rotate() 
+    {};
+  Rotate(double);
+  Rotate(double a, double b, double c, double d)
+    {m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d;}
+  Rotate(const Matrix& a)
+    { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; 
+      m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1];}
+};
+ostream& operator<<(ostream&, const Rotate&);
+istream& operator>>(istream&, Rotate&);
+
+class BBox {
+public:
+  Vector ll;
+  Vector ur;
+
+public:
+  BBox() 
+    {}
+  BBox(double w, double h) 
+    {ll.v[0] = 0; ll.v[1] = 0; ur.v[0] = w; ur.v[1] = h;}
+  BBox(const Vector& v) 
+    {ll=v; ur=v;}
+  BBox(double, double, double, double);
+  BBox(const Vector&, const Vector&);
+
+  BBox(const BBox& a)
+    {ll=a.ll; ur=a.ur;}
+  BBox& operator=(const BBox& a)
+    {ll=a.ll; ur=a.ur; return *this;}
+
+  Vector lr() {return Vector(ur[0],ll[1]);}
+  Vector ul() {return Vector(ll[0],ur[1]);}
+
+  BBox& operator+=(const Vector& v) // addition
+    {ll+=v; ur+=v; return *this;}
+  BBox& operator-=(const Vector& a) // subtraction
+    {ll-=a; ur-=a; return *this;}
+  BBox& operator*=(const Matrix& m) // multiply
+    {ll*=m; ur*=m; return *this;}
+
+  Vector center() 
+  {return (ur-ll)/2 + ll;}
+  Vector size() 
+  {return ur - ll;}
+  int isEmpty() const 
+  {Vector v = ur-ll; return (v[0]==0 && v[1]==0);}
+  int isIn(const Vector&) const;
+  int isIn(const BBox&) const;
+
+  BBox& expand(double a) 
+    {ll-=Vector(a,a); ur+=Vector(a,a); return *this;}
+  BBox& expand(const Vector& v) 
+    {ll-=v; ur+=v; return *this;}
+  BBox& shrink(double a) 
+    {ll+=Vector(a,a); ur-=Vector(a,a); return *this;}
+  BBox& shrink(const Vector& v) 
+    {ll+=v; ur-=v; return *this;}
+  BBox& bound(BBox);
+  BBox& bound(const Vector&);
+};
+ostream& operator<<(ostream&, const BBox&);
+
+inline BBox operator+(const BBox& b, const Vector& v) {return BBox(b) += v;}
+inline BBox operator-(const BBox& b, const Vector& v) {return BBox(b) -= v;}
+inline BBox operator*(const BBox& b, const Matrix& m) {return BBox(b) *= m;}
+
+BBox intersect(const BBox&, const BBox&);
+
+// clip line to box from origin to width,height
+int clip(Vector*, Vector*, int, int);
+#endif
+
+
+
+
diff --git a/tksao/vector/vector3d.C b/tksao/vector/vector3d.C
new file mode 100644
index 0000000..f0ed8d0
--- /dev/null
+++ b/tksao/vector/vector3d.C
@@ -0,0 +1,473 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include <tk.h>
+
+#include "vector3d.h"
+#include "vector.h"
+#include "fuzzy.h"
+
+// Vector3d
+
+Vector3d::Vector3d(const Vector& a)
+{
+  v[0]=a.v[0];
+  v[1]=a.v[1];
+  v[2]=0;
+  v[3]=1;
+}
+
+Vector3d::Vector3d(const Vector& a, double z)
+{
+  v[0]=a.v[0];
+  v[1]=a.v[1];
+  v[2]=z;
+  v[3]=1;
+}
+
+Vector3d& Vector3d::operator=(const Vector& a)
+{
+  v[0]=a.v[0];
+  v[1]=a.v[1];
+  v[2]=0;
+  v[3]=1;
+  return *this;
+}
+
+Vector Vector3d::TkCanvasPs(void* canvas)
+{
+  return Vector(v[0], Tk_CanvasPsY((Tk_Canvas)canvas, v[1]));
+}
+
+ostream& operator<<(ostream& os, const Vector3d& v)
+{
+  unsigned char sep = (unsigned char)os.iword(Vector::separator);
+  if (!sep)
+    sep = ' ';
+
+  unsigned char unit = (unsigned char)os.iword(Vector::unit);
+  if (!unit)
+    os << v.v[0] << sep << v.v[1] << sep << v.v[2];
+  else
+    os << v.v[0] << unit << v.v[1] << unit << v.v[2] << unit;
+
+  // reset unit
+  os.iword(Vector::unit) = '\0';
+
+  return os;
+}
+
+istream& operator>>(istream& s, Vector3d& v)
+{
+  s >> v.v[0] >> v.v[1] >> v.v[2];
+  return s;
+}
+
+// Vertex3d
+
+ostream& operator<<(ostream& os, const Vertex3d& v)
+{
+  os << v.vector;
+  return os;
+}
+
+// Matrix3d
+
+Matrix3d& Matrix3d::operator*=(const Matrix3d& a)
+{
+  Matrix3d r;
+  for (int ii=0; ii<4; ii++)
+    for (int jj=0; jj<4; jj++)
+      r.m[ii][jj] = 
+	m[ii][0]*a.m[0][jj] + 
+	m[ii][1]*a.m[1][jj] + 
+	m[ii][2]*a.m[2][jj] + 
+	m[ii][3]*a.m[3][jj];
+
+  return *this=r;
+}
+
+Matrix3d::Matrix3d(const Matrix& a)
+{ 
+  m[0][0]=a.m[0][0];
+  m[0][1]=a.m[0][1];
+  m[0][2]=0;
+  m[0][3]=0;
+
+  m[1][0]=a.m[1][0];
+  m[1][1]=a.m[1][1];
+  m[1][2]=0;
+  m[1][3]=0;
+
+  m[2][0]=0;
+  m[2][1]=0;
+  m[2][2]=1;
+  m[2][3]=0;
+  
+  m[3][0]=a.m[2][0];
+  m[3][1]=a.m[2][1];
+  m[3][2]=0;
+  m[3][3]=1;
+}
+
+Matrix3d Matrix3d::invert()
+{
+  Matrix3d cc = this->cofactor();
+  Matrix3d aa = cc.adjoint();
+  double dd = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + 
+    m[0][2]*aa.m[2][0] + m[0][3]*aa.m[3][0];
+
+  Matrix3d rr;
+  for (int ii=0; ii<4; ii++ )
+    for (int jj=0; jj<4; jj++)
+      rr.m[ii][jj] = aa.m[ii][jj]/dd;
+
+  return rr;
+}
+
+Matrix3d Matrix3d::cofactor()
+{
+  Matrix3d rr;
+
+  rr.m[0][0] = +det2d(m[1][1],m[1][2],m[1][3],
+		      m[2][1],m[2][2],m[2][3],
+		      m[3][1],m[3][2],m[3][3]);
+  rr.m[0][1] = -det2d(m[1][0],m[1][2],m[1][3],
+		      m[2][0],m[2][2],m[2][3],
+		      m[3][0],m[3][2],m[3][3]);
+  rr.m[0][2] = +det2d(m[1][0],m[1][1],m[1][3],
+		      m[2][0],m[2][1],m[2][3],
+		      m[3][0],m[3][1],m[3][3]);
+  rr.m[0][3] = -det2d(m[1][0],m[1][1],m[1][2],
+		      m[2][0],m[2][1],m[2][2],
+		      m[3][0],m[3][1],m[3][2]);
+
+  rr.m[1][0] = -det2d(m[0][1],m[0][2],m[0][3],
+		      m[2][1],m[2][2],m[2][3],
+		      m[3][1],m[3][2],m[3][3]);
+  rr.m[1][1] = +det2d(m[0][0],m[0][2],m[0][3],
+		      m[2][0],m[2][2],m[2][3],
+		      m[3][0],m[3][2],m[3][3]);
+  rr.m[1][2] = -det2d(m[0][0],m[0][1],m[0][3],
+		      m[2][0],m[2][1],m[2][3],
+		      m[3][0],m[3][1],m[3][3]);
+  rr.m[1][3] = +det2d(m[0][0],m[0][1],m[0][2],
+		      m[2][0],m[2][1],m[2][2],
+		      m[3][0],m[3][1],m[3][2]);
+
+  rr.m[2][0] = +det2d(m[0][1],m[0][2],m[0][3],
+		      m[1][1],m[1][2],m[1][3],
+		      m[3][1],m[3][2],m[3][3]);
+  rr.m[2][1] = -det2d(m[0][0],m[0][2],m[0][3],
+		      m[1][0],m[1][2],m[1][3],
+		      m[3][0],m[3][2],m[3][3]);
+  rr.m[2][2] = +det2d(m[0][0],m[0][1],m[0][3],
+		      m[1][0],m[1][1],m[1][3],
+		      m[3][0],m[3][1],m[3][3]);
+  rr.m[2][3] = -det2d(m[0][0],m[0][1],m[0][2],
+		      m[1][0],m[1][1],m[1][2],
+		      m[3][0],m[3][1],m[3][2]);
+
+  rr.m[3][0] = -det2d(m[0][1],m[0][2],m[0][3],
+		      m[1][1],m[1][2],m[1][3],
+		      m[2][1],m[2][2],m[2][3]);
+  rr.m[3][1] = +det2d(m[0][0],m[0][2],m[0][3],
+		      m[1][0],m[1][2],m[1][3],
+		      m[2][0],m[2][2],m[2][3]);
+  rr.m[3][2] = -det2d(m[0][0],m[0][1],m[0][3],
+		      m[1][0],m[1][1],m[1][3],
+		      m[2][0],m[2][1],m[2][3]);
+  rr.m[3][3] = +det2d(m[0][0],m[0][1],m[0][2],
+		      m[1][0],m[1][1],m[1][2],
+		      m[2][0],m[2][1],m[2][2]);
+
+  return rr;
+}
+
+Matrix3d Matrix3d::adjoint()
+{
+  Matrix3d rr;
+  for (int ii=0; ii<4; ii++)
+    for (int jj=0; jj<4; jj++)
+      rr.m[jj][ii] = m[ii][jj];
+
+  return rr;
+}
+
+double Matrix3d::det()
+{
+  Matrix3d cc = this->cofactor();
+  Matrix3d aa = cc.adjoint();
+  return m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] 
+    + m[0][2]*aa.m[2][0] + m[0][3]*aa.m[3][0];
+}
+
+void Matrix3d::dump()
+{
+  for (int ii=0; ii<4; ii++) {
+    for (int jj=0; jj<4; jj++)
+      cerr << m[ii][jj] << ' ';
+    cerr << endl;
+  }
+  cerr << endl;
+}
+
+ostream& operator<<(ostream& os, const Matrix3d& m)
+{
+  os << ' ';
+  for (int ii=0; ii<4; ii++)
+    for (int jj=0; jj<3; jj++)
+      os << m.m[ii][jj] << ' ';
+
+  return os;
+}
+
+istream& operator>>(istream& s, Matrix3d& m)
+{
+  for (int ii=0; ii<4; ii++ )
+    for (int jj=0; jj<3; jj++)
+      s >> m.m[ii][jj];
+
+  return s;
+}
+
+// Translate3d
+
+Translate3d::Translate3d(const Vector& v) 
+{
+  m[3][0]=v.v[0];
+  m[3][1]=v.v[1];
+  m[3][2]=0;
+}
+
+Translate3d::Translate3d(const Vector& v, double z) 
+{
+  m[3][0]=v.v[0]; 
+  m[3][1]=v.v[1]; 
+  m[3][2]=z;
+}
+
+ostream& operator<<(ostream& os, const Translate3d& m)
+{
+  os << ' ' << m.m[3][0] << ' ' << m.m[3][1] << ' ' << m.m[3][2] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, Translate3d& m)
+{
+  s >> m.m[3][0] >> m.m[3][1] >> m.m[3][2];
+  return s;
+}
+
+// Scale3d
+
+Scale3d::Scale3d(const Vector& v)
+{
+  m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=1;
+}
+
+Scale3d::Scale3d(const Vector& v, double c)
+{
+  m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=c;
+}
+
+ostream& operator<<(ostream& os, const Scale3d& m)
+{
+  os << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' ' << m.m[2][2] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, Scale3d& m)
+{
+  s >> m.m[0][0] >> m.m[1][1] >> m.m[2][2];
+  return s;
+}
+
+// RotateX3d
+
+RotateX3d::RotateX3d(double a) : Matrix3d()
+{
+  m[1][1] = cos(a);
+  m[1][2] = sin(a);
+  m[2][1] = -sin(a);
+  m[2][2] = cos(a);
+
+  // this fixes a problem with numbers too small and tring to invert the matrix
+  tzero(&m[1][1]);
+  tzero(&m[1][2]);
+  tzero(&m[2][1]);
+  tzero(&m[2][2]);
+}
+
+ostream& operator<<(ostream& os, const RotateX3d& m)
+{
+  os << ' ' << m.m[1][1] << ' ' << m.m[1][2]
+    << ' ' << m.m[2][1] << ' ' << m.m[2][2] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, RotateX3d& m)
+{
+  s >> m.m[1][1] >> m.m[1][2] >> m.m[2][1] >> m.m[2][2];
+  return s;
+}
+
+// RotateY3d
+
+RotateY3d::RotateY3d(double a) : Matrix3d()
+{
+  m[0][0] = cos(a);
+  m[0][2] = -sin(a);
+  m[2][0] = sin(a);
+  m[2][2] = cos(a);
+
+  // this fixes a problem with numbers too small and tring to invert the matrix
+  tzero(&m[0][0]);
+  tzero(&m[0][2]);
+  tzero(&m[2][0]);
+  tzero(&m[2][2]);
+}
+
+ostream& operator<<(ostream& os, const RotateY3d& m)
+{
+  os << ' ' << m.m[0][0] << ' ' << m.m[0][2]
+    << ' ' << m.m[2][0] << ' ' << m.m[2][2] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, RotateY3d& m)
+{
+  s >> m.m[0][0] >> m.m[0][2] >> m.m[2][0] >> m.m[2][2];
+  return s;
+}
+
+// RotateZ3d
+
+RotateZ3d::RotateZ3d(double a) : Matrix3d()
+{
+  m[0][0] = cos(a);
+  m[0][1] = sin(a);
+  m[1][0] = -sin(a);
+  m[1][1] = cos(a);
+
+  // this fixes a problem with numbers too small and tring to invert the matrix
+  tzero(&m[0][0]);
+  tzero(&m[0][1]);
+  tzero(&m[1][0]);
+  tzero(&m[1][1]);
+}
+
+ostream& operator<<(ostream& os, const RotateZ3d& m)
+{
+  os << ' ' << m.m[0][0] << ' ' << m.m[0][1]
+    << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' ';
+  return os;
+}
+
+istream& operator>>(istream& s, RotateZ3d& m)
+{
+  s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1];
+  return s;
+}
+
+// BBox3d
+
+BBox3d::BBox3d(double a, double b, double c, double d, double e, double f)
+{
+  // we want a 'positive' cube
+  ll.v[0] = a < d ? a : d;
+  ll.v[1] = b < e ? b : e;
+  ll.v[2] = c < f ? c : f;
+  ur.v[0] = a < d ? d : a;
+  ur.v[1] = b < e ? e : b;
+  ur.v[2] = c < f ? f : c;
+}
+
+BBox3d::BBox3d(const Vector3d& l, const Vector3d& h)
+{
+  // we want a 'positive' cube
+  ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0];
+  ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1];
+  ll.v[2] = l.v[2] < h.v[2] ? l.v[2] : h.v[2];
+  ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0];
+  ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1];
+  ur.v[2] = l.v[2] < h.v[2] ? h.v[2] : l.v[2];
+}
+
+int BBox3d::isIn(const Vector3d& v) const
+{
+  return !(v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || v.v[2] < ll.v[2] ||
+	   v.v[0] > ur.v[0] || v.v[1] > ur.v[1] || v.v[2] > ur.v[2]);
+}
+
+double BBox3d::volume()
+{
+  Vector3d aa =ur-ll;
+  return aa[0]*aa[1]*aa[2];
+}
+
+BBox3d& BBox3d::bound(const Vector3d& vv)
+{
+  // expand bbox3d to include vector3d
+  if (vv.v[0] < ll[0])
+    ll[0] = vv.v[0];
+  if (vv.v[1] < ll[1])
+    ll[1] = vv.v[1];
+  if (vv.v[2] < ll[2])
+    ll[2] = vv.v[2];
+    
+  if (vv.v[0] > ur[0])
+    ur[0] = vv.v[0];
+  if (vv.v[1] > ur[1])
+    ur[1] = vv.v[1];
+  if (vv.v[2] > ur[2])
+    ur[2] = vv.v[2];
+
+  return *this;
+}
+
+ostream& operator<<(ostream& os, const BBox3d& b)
+{
+  os << b.ll << ' ' << b.ur;
+  return os;
+}
+
+// WorldToView
+
+Matrix3d WorldToView3d(const Vector3d& cop, 
+		       const Vector3d& vpn, 
+		       const Vector3d& vup)
+{
+  Vector3d zv = ((Vector3d)vpn).normalize();
+  Vector3d xv = cross(zv,(Vector3d&)vup).normalize();
+  Vector3d yv = cross(xv,zv).normalize();
+
+  return Translate3d(-cop) * 
+    Matrix3d(xv[0],yv[0],zv[0],
+	     xv[1],yv[1],zv[1],
+	     xv[2],yv[2],zv[2],
+	         0,    0,    0);
+}
+
+Matrix3d WorldToView3d(const Vector3d& cop, 
+		       double head, double pitch, double bank)
+{
+  return Translate3d(-cop) *
+    RotateY3d(head) * 
+    RotateX3d(pitch) * 
+    RotateZ3d(bank) *
+    Scale3d(1,1,-1);
+}
+
+Matrix3d WorldToView3d(const Vector3d& cop, const Vector3d& vpn, double bank)
+{
+  Vector3d zv = -((Vector3d)vpn).normalize();
+  double l=sqrt(zv[0]*zv[0]+zv[2]*zv[2]);
+
+  return Translate3d(-cop) *
+    RotateY3d(zv[2]/l,zv[0]/l,-zv[0]/l,zv[2]/l) * 
+    RotateX3d(l,zv[1],-zv[1],l) * 
+    RotateZ3d(bank) *
+    Scale3d(1,1,-1);
+}
diff --git a/tksao/vector/vector3d.h b/tksao/vector/vector3d.h
new file mode 100644
index 0000000..03aeb58
--- /dev/null
+++ b/tksao/vector/vector3d.h
@@ -0,0 +1,395 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#ifndef __vector3d_h__
+#define __vector3d_h__
+
+#include <math.h>
+#include <float.h>
+
+#include <iostream>
+using namespace std;
+
+class Vector;
+class Matrix;
+class Matrix3d;
+
+class Vector3d {
+ public:
+  double v[4];
+
+ public:
+  Vector3d() 
+    {v[0]=0;v[1]=0;v[2]=0;v[3]=1;}
+  Vector3d(double* f) 
+    {v[0]=f[0]; v[1]=f[1]; v[2]=f[2]; v[3]=1;}
+  Vector3d(double x, double y) 
+    {v[0]=x;v[1]=y;v[2]=0;v[3]=1;}
+  Vector3d(double x, double y, double z) 
+    {v[0]=x;v[1]=y;v[2]=z;v[3]=1;}
+
+  Vector3d(const Vector&);
+  Vector3d(const Vector&, double);
+  Vector3d& operator=(const Vector&);
+  
+  Vector3d(const Vector3d& a) 
+    {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];v[3]=a.v[3];}
+  Vector3d& operator=(const Vector3d& a) 
+    {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];v[3]=a.v[3]; return *this;}
+
+  double& operator[](int i) {return v[i];} // return element
+  double* vv() {return v;} // return vector
+
+  Vector3d& operator+=(const Vector3d& a) // addition
+    {v[0]+=a.v[0]; v[1]+=a.v[1]; v[2]+=a.v[2]; return *this;}
+  Vector3d& operator-=(const Vector3d& a) // subtraction
+    {v[0]-=a.v[0]; v[1]-=a.v[1]; v[2]-=a.v[2]; return *this;}
+  Vector3d& operator*=(double f) // scalar multiply
+    {v[0]*=f; v[1]*=f; v[2]*=f; return *this;}
+  Vector3d& operator/=(double f) // scalar division
+    {v[0]/=f; v[1]/=f; v[2]/=f; return *this;}
+  Vector3d& operator*=(const Matrix3d&); // vector multiply
+
+  Vector3d abs()
+  {return Vector3d(fabs(v[0]),fabs(v[1]),fabs(v[2]));}
+  double angleX() 
+  {return atan2(v[2],v[1]);}
+  double angleY() 
+  {return atan2(v[0],v[2]);}
+  double angleZ() 
+  {return atan2(v[1],v[0]);}
+  Vector3d ceil()
+  {return Vector3d(::ceil(v[0]),::ceil(v[1]),::ceil(v[2]));}
+  Vector3d floor()
+  {return Vector3d(::floor(v[0]),::floor(v[1]),::floor(v[2]));}
+  Vector3d invert()
+  {return Vector3d(1/v[0],1/v[1],1/v[2]);}
+  double length()
+  {return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);}
+  Vector3d round()
+  {return Vector3d((int)(v[0]+.5),(int)(v[1]+.5),(int)(v[2]+.5));}
+  Vector3d normalize()
+  {double d = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+    return d ? Vector3d(v[0]/d,v[1]/d,v[2]/d) : Vector3d();}
+  Vector3d project()
+  {return (v[3]!=1) ? Vector3d(v[0]/v[3],v[1]/v[3],v[2]/v[3]) : *this;}
+  Vector TkCanvasPs(void* canvas);
+};
+ostream& operator<<(ostream&, const Vector3d&);
+istream& operator>>(istream&, Vector3d&);
+
+inline Vector3d operator-(const Vector3d& a) 
+{return Vector3d(-a.v[0],-a.v[1],-a.v[2]);}
+inline Vector3d operator+(const Vector3d& a, const Vector3d& b)
+{return Vector3d(a) +=b;}
+inline Vector3d operator-(const Vector3d& a, const Vector3d& b)
+{return Vector3d(a) -=b;}
+inline Vector3d operator*(const Vector3d& a, double b) 
+{return Vector3d(a) *=b;}
+inline Vector3d operator/(const Vector3d& a, double b)
+{return Vector3d(a) /=b;}
+inline Vector3d operator*(const Vector3d& v, const Matrix3d& m)
+{return Vector3d(v) *=m;}
+inline double operator*(const Vector3d& a, const Vector3d& b) // dot product
+{double r=0; r+=a.v[0]*b.v[0]; r+=a.v[1]*b.v[1]; r+=a.v[2]*b.v[2]; return r;}
+inline Vector3d cross(Vector3d& a, Vector3d& b) // cross product
+{return Vector3d(a[1]*b[2]-b[1]*a[2],a[2]*b[0]-b[2]*a[0],a[0]*b[1]-b[0]*a[1]);}
+
+class Vertex3d {
+public:
+  Vector3d vector;
+
+private:
+  Vertex3d* next_;
+  Vertex3d* previous_;
+
+public:
+  Vertex3d() 
+    {next_=NULL; previous_=NULL;}
+  Vertex3d(double x, double y, double z) 
+    {vector=Vector3d(x,y,z); next_=NULL; previous_=NULL;}
+  Vertex3d(const Vector3d& a) 
+    {vector=a; next_=NULL; previous_=NULL;}
+
+  Vertex3d(const Vertex3d& a) 
+    {vector=a.vector; next_=a.next_; previous_=a.previous_;}
+  Vertex3d& operator=(const Vertex3d& a)
+    {vector=a.vector; next_=a.next_; previous_=a.previous_; return *this;}
+
+  Vertex3d* next()
+  {return next_;}
+  Vertex3d* previous()
+  {return previous_;}
+  void setNext(Vertex3d* v)
+  {next_ = v;}
+  void setPrevious(Vertex3d* v)
+  {previous_ = v;}
+};
+ostream& operator<<(ostream&, const Vertex3d&);
+
+class Matrix3d {
+ public:
+  double m[4][4];
+
+ public:
+  Matrix3d()
+    { m[0][0]=1; m[0][1]=0; m[0][2]=0; m[0][3]=0;
+      m[1][0]=0; m[1][1]=1; m[1][2]=0; m[1][3]=0;
+      m[2][0]=0; m[2][1]=0; m[2][2]=1; m[2][3]=0;
+      m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; }
+
+  Matrix3d(double a, double b, double c,
+	   double d, double e, double f,
+	   double g, double h, double i,
+	   double j, double k, double l)
+    { m[0][0]=a; m[0][1]=b; m[0][2]=c; m[0][3]=0;
+      m[1][0]=d; m[1][1]=e; m[1][2]=f; m[1][3]=0;
+      m[2][0]=g; m[2][1]=h; m[2][2]=i; m[2][3]=0;
+      m[3][0]=j; m[3][1]=k; m[3][2]=l; m[3][3]=1; }
+
+  Matrix3d(Vector3d& x, Vector3d& y, Vector3d& z)
+    { m[0][0]=x[0]; m[0][1]=y[0]; m[0][2]=z[0]; m[0][3]=0;
+      m[1][0]=x[1]; m[1][1]=y[1]; m[1][2]=z[1]; m[1][3]=0;
+      m[2][0]=x[2]; m[2][1]=y[2]; m[2][2]=z[2]; m[2][3]=0;
+      m[3][0]=0;    m[3][1]=0;    m[3][2]=0;    m[3][3]=1; }
+
+  Matrix3d(const Matrix3d& a)
+    { m[0][0]=a.m[0][0];m[0][1]=a.m[0][1];m[0][2]=a.m[0][2];m[0][3]=a.m[0][3];
+      m[1][0]=a.m[1][0];m[1][1]=a.m[1][1];m[1][2]=a.m[1][2];m[1][3]=a.m[1][3];
+      m[2][0]=a.m[2][0];m[2][1]=a.m[2][1];m[2][2]=a.m[2][2];m[2][3]=a.m[2][3];
+      m[3][0]=a.m[3][0];m[3][1]=a.m[3][1];m[3][2]=a.m[3][2];m[3][3]=a.m[3][3]; }
+
+  Matrix3d& operator=(const Matrix3d& a)
+    { m[0][0]=a.m[0][0];m[0][1]=a.m[0][1];m[0][2]=a.m[0][2];m[0][3]=a.m[0][3];
+      m[1][0]=a.m[1][0];m[1][1]=a.m[1][1];m[1][2]=a.m[1][2];m[1][3]=a.m[1][3];
+      m[2][0]=a.m[2][0];m[2][1]=a.m[2][1];m[2][2]=a.m[2][2];m[2][3]=a.m[2][3];
+      m[3][0]=a.m[3][0];m[3][1]=a.m[3][1];m[3][2]=a.m[3][2];m[3][3]=a.m[3][3];
+      return *this;}
+  
+  Matrix3d(const Matrix& a);
+  double matrix(int i, int j)
+  {return m[i][j];} // return element
+  Vector3d operator[](int i)
+  {return Vector3d(m[i]);} // return row
+  double* mm()
+  {return (double*)m;} // return matrix
+
+  Matrix3d& identity()
+    { m[0][0]=1; m[0][1]=0; m[0][2]=0; m[0][3]=0;
+      m[1][0]=0; m[1][1]=1; m[1][2]=0; m[1][3]=0;
+      m[2][0]=0; m[2][1]=0; m[2][2]=1; m[2][3]=0;
+      m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; return *this;}
+  Matrix3d& operator*=(const Matrix3d&); // matrix multiply
+
+  Matrix3d invert();
+  Matrix3d cofactor();
+  Matrix3d adjoint();
+  double det();
+  double det2d(double& a, double& b, double& c, 
+	       double& d, double& e, double& f,
+	       double& g, double& h, double& i)
+  {return a*(e*i-f*h) - b*(d*i-f*g) + c*(d*h-e*g);}
+
+  void dump();
+};
+ostream& operator<<(ostream&, const Matrix3d&);
+istream& operator>>(istream&, Matrix3d&);
+
+inline Matrix3d operator*(const Matrix3d& a, const Matrix3d& b)
+{return Matrix3d(a) *= b;}
+inline Vector3d& Vector3d::operator*=(const Matrix3d& m)
+{
+  double vv[4];
+  double* mm = (double*)(m.m);
+  vv[0] = v[0]*mm[0] + v[1]*mm[4] + v[2]*mm[8]  + v[3]*mm[12]; 
+  vv[1] = v[0]*mm[1] + v[1]*mm[5] + v[2]*mm[9]  + v[3]*mm[13]; 
+  vv[2] = v[0]*mm[2] + v[1]*mm[6] + v[2]*mm[10] + v[3]*mm[14];
+  vv[3] = v[0]*mm[3] + v[1]*mm[7] + v[2]*mm[11] + v[3]*mm[15];
+  v[0] = vv[0];
+  v[1] = vv[1];
+  v[2] = vv[2];
+  v[3] = vv[3];
+  return *this;
+}
+
+class Translate3d : public Matrix3d {
+public:
+  Translate3d() {};
+  Translate3d(double x, double y, double z) 
+    {m[3][0]=x; m[3][1]=y; m[3][2]=z;}
+  Translate3d(const Vector3d& v) 
+    {m[3][0]=v.v[0]; m[3][1]=v.v[1]; m[3][2]=v.v[2];}
+  Translate3d(const Vector& v);
+  Translate3d(const Vector& v, double z);
+  Translate3d(const Matrix3d& a)
+    {m[3][0] = a.m[3][0]; m[3][1] = a.m[3][1]; m[3][2] = a.m[3][2];}
+};
+ostream& operator<<(ostream&, const Translate3d&);
+istream& operator>>(istream&, Translate3d&);
+
+class Scale3d : public Matrix3d {
+public:
+  Scale3d() {};
+  Scale3d(double a)
+    {m[0][0]=a; m[1][1]=a; m[2][2]=a;}
+  Scale3d(double a, double b)
+    {m[0][0]=a; m[1][1]=a; m[2][2]=b;}
+  Scale3d(double a, double b, double c)
+    {m[0][0]=a; m[1][1]=b; m[2][2]=c;}
+  Scale3d(const Vector& v);
+  Scale3d(const Vector& v, double c);
+  Scale3d(const Vector3d& v)
+    {m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=v.v[2];}
+  Scale3d(const Matrix3d& a) 
+    {m[0][0] = a.m[0][0]; m[1][1] = a.m[1][1]; m[2][2] = a.m[2][2];}
+};
+ostream& operator<<(ostream&, const Scale3d&);
+istream& operator>>(istream&, Scale3d&);
+
+class FlipX3d : public Matrix3d {
+public:
+  FlipX3d() 
+    {m[0][0] = -1;}
+};
+
+class FlipY3d : public Matrix3d {
+public:
+  FlipY3d() 
+    {m[1][1] = -1;}
+};
+
+class FlipZ3d : public Matrix3d {
+public:
+  FlipZ3d() 
+    {m[2][2] = -1;}
+};
+
+class FlipXY3d : public Matrix3d {
+public:
+  FlipXY3d() 
+    {m[0][0] = -1; m[1][1] = -1;}
+};
+
+class FlipXYZ3d : public Matrix3d {
+public:
+  FlipXYZ3d() 
+    {m[0][0] = -1; m[1][1] = -1; m[2][2] = -1;}
+};
+
+class RotateX3d : public Matrix3d {
+public:
+  RotateX3d()
+    {};
+  RotateX3d(double);
+  RotateX3d(double a, double b, double c, double d)
+    {m[1][1] = a; m[1][2] = b; m[2][1] = c; m[2][2] = d;}
+};
+ostream& operator<<(ostream&, const RotateX3d&);
+istream& operator>>(istream&, RotateX3d&);
+
+class RotateY3d : public Matrix3d {
+public:
+  RotateY3d()
+    {};
+  RotateY3d(double);
+  RotateY3d(double a, double b, double c, double d)
+    {m[0][0] = a; m[0][2] = b; m[2][0] = c; m[2][2] = d;}
+};
+ostream& operator<<(ostream&, const RotateY3d&);
+istream& operator>>(istream&, RotateY3d&);
+
+class RotateZ3d : public Matrix3d {
+public:
+  RotateZ3d()
+    {};
+  RotateZ3d(double);
+  RotateZ3d(double a, double b, double c, double d)
+    {m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d;}
+};
+ostream& operator<<(ostream&, const RotateZ3d&);
+istream& operator>>(istream&, RotateZ3d&);
+
+class Shear3d : public Matrix3d {
+ public:
+  Shear3d(double x, double y, double dist)
+    {m[2][0] = -x/dist; m[2][1] = -y/dist;}
+};
+
+class Perspective3d : public Matrix3d {
+ public:
+  Perspective3d(double front, double back, double dist)
+  { m[2][2] = back/(back-front); m[2][3] = 1; 
+    m[3][2] = -front*back/(dist*(back-front)); m[3][3] = 0;}
+};
+
+class BBox3d {
+ public:
+  Vector3d ll;
+  Vector3d ur;
+
+ public:
+  BBox3d()
+    {}
+  BBox3d(double w, double h, double d) 
+    {ll.v[0]=0; ll.v[1]=0; ll.v[2]=0; ur.v[0]=w; ur.v[1]=h; ur.v[2]=d;}
+  BBox3d(const Vector3d& v) 
+    {ll=v; ur=v;}
+  BBox3d(double, double, double, double, double, double);
+  BBox3d(const Vector3d&, const Vector3d&);
+
+  BBox3d(const BBox3d& a)
+    {ll=a.ll; ur=a.ur;}
+  BBox3d& operator=(const BBox3d& a)
+    {ll=a.ll; ur=a.ur; return *this;}
+
+  BBox3d& operator+=(const Vector3d& v) // addition
+    {ll+=v; ur+=v; return *this;}
+  BBox3d& operator-=(const Vector3d& a) // subtraction
+    {ll-=a; ur-=a; return *this;}
+  BBox3d& operator*=(const Matrix3d& m) // multiplication
+    {ll*=m; ur*=m; return *this;} 
+
+  double volume();
+  Vector3d center() 
+  {return (ur-ll)/2 + ll;}
+  Vector3d size() 
+  {return ur - ll;}
+  int isEmpty() const 
+  {Vector3d v = ur-ll; return (v[0]==0 && v[1]==0 && v[2]==0);}
+  int isIn(const Vector3d&) const;
+
+  BBox3d& expand(double a) 
+    {ll-=Vector3d(a,a,a); ur+=Vector3d(a,a,a); return *this;}
+  BBox3d& expand(const Vector3d& v) 
+    {ll-=v; ur+=v; return *this;}
+  BBox3d& shrink(double a) 
+    {ll+=Vector3d(a,a,a); ur-=Vector3d(a,a,a); return *this;}
+  BBox3d& shrink(const Vector3d& v) 
+    {ll+=v; ur-=v; return *this;}
+  // expand bbox3d to include vector3d
+  BBox3d& bound(const Vector3d&);
+};
+ostream& operator<<(ostream&, const BBox3d&);
+
+inline BBox3d operator+(const BBox3d& b, const Vector3d& v)
+{return BBox3d(b) += v;}
+inline BBox3d operator-(const BBox3d& b, const Vector3d& v)
+{return BBox3d(b) -= v;}
+inline BBox3d operator*(const BBox3d& b, const Matrix3d& m)
+{return BBox3d(b) *= m;}
+
+// WorldToView
+
+Matrix3d WorldToView3d(const Vector3d& cop, 
+		       const Vector3d& vpn, 
+		       const Vector3d& vup);
+Matrix3d WorldToView3d(const Vector3d& cop, 
+		       double head, double pitch, double bank);
+Matrix3d WorldToView3d(const Vector3d& cop, const Vector3d& vpn, double bank);
+
+#endif
+
+
+
+
diff --git a/tksao/vector/vectorold.C b/tksao/vector/vectorold.C
new file mode 100644
index 0000000..a0d9e4a
--- /dev/null
+++ b/tksao/vector/vectorold.C
@@ -0,0 +1,884 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include "vector.h"
+
+// Vector
+
+Vector operator-(const Vector& a)
+{
+  return Vector(-a.v[0], -a.v[1]);
+}
+
+Vector& Vector::operator+=(const Vector& a)
+{
+  v[0]+=a.v[0];
+  v[1]+=a.v[1];
+  return *this;
+}
+
+Vector operator+(const Vector& a, const Vector& b)
+{
+  Vector r=a;
+  r+=b;
+  return r;
+}
+
+Vector& Vector::operator-=(const Vector& a)
+{
+  v[0]-=a.v[0];
+  v[1]-=a.v[1];
+  return *this;
+}
+
+Vector operator-(const Vector& a, const Vector& b)
+{
+  Vector r=a;
+  r-=b;
+  return r;
+}
+
+Vector& Vector::operator*=(double f)
+{
+  v[0]*=f;
+  v[1]*=f;
+  return *this;
+}
+
+Vector operator*(const Vector& a, double b)
+{
+  Vector r=a;
+  r*=b;
+  return r;
+}
+
+Vector& Vector::operator/=(double f)
+{
+  v[0]/=f;
+  v[1]/=f;
+  return *this;
+}
+
+Vector operator/(const Vector& a, double b)
+{
+  Vector r=a;
+  r/=b;
+  return r;
+}
+
+double operator*(const Vector& a, const Vector& b) //  2D dot product
+{
+  double r = 0;
+  for (int i=0; i<2; i++)
+    r += a.v[i]*b.v[i];
+  return r;
+}
+
+Vector Vector::normalize()    // 2D normalize
+{
+  Vector r = *this;
+  double d = sqrt(v[0]*v[0]+v[1]*v[1]);
+  if (d) {
+    r[0] /= d;
+    r[1] /= d;
+  }
+  else {
+    r[0] = 0;
+    r[1] = 0;
+  }
+
+  return r;
+}
+
+Vector Vector::abs() // absolute value
+{
+  Vector r;
+  r[0] = fabs(v[0]);
+  r[1] = fabs(v[1]);
+  r[2] = 1;
+  return r;
+}
+
+double Vector::length()
+{
+  return sqrt(v[0]*v[0] + v[1]*v[1]);
+}
+
+double Vector::angle()
+{
+  return atan2(v[1],v[0]);
+}
+
+double Vector::avg()
+{
+  return (v[0]+v[1])/2;
+}
+
+double Vector::max()
+{
+  return v[0]>v[1]?v[0]:v[1];
+}
+
+double Vector::min()
+{
+  return v[0]<v[1]?v[0]:v[1];
+}
+
+Vector Vector::round()
+{
+  return Vector((int)(v[0]+.5), (int)(v[1]+.5));
+}
+
+Vector Vector::floor()
+{
+  return Vector(::floor(v[0]), ::floor(v[1]));
+}
+
+Vector Vector::ceil()
+{
+  return Vector(::ceil(v[0]), ::ceil(v[1]));
+}
+
+Vector Vector::invert()
+{
+  Vector r = *this;
+  r[0]=1/r[0];
+  r[1]=1/r[1];
+  r[2]=1;
+  return r;
+}
+
+Vector Vector::TkCanvasPs(Tk_Canvas canvas)
+{
+  return Vector(v[0], Tk_CanvasPsY(canvas, v[1]));
+}
+
+Vector& Vector::operator*=(const Matrix& m)
+{
+  Vector vv = *this;
+  for (int i=0; i<3; i++)
+    v[i] = vv.v[0]*m.m[0][i] + vv.v[1]*m.m[1][i] + vv.v[2]*m.m[2][i];
+
+  return *this;
+}
+
+Vector operator*(const Vector& v, const Matrix& m)
+{
+  Vector r;
+  for (int i=0; i<3; i++)
+    r.v[i] = v.v[0]*m.m[0][i] + v.v[1]*m.m[1][i] + v.v[2]*m.m[2][i];
+
+  return r;
+}
+
+ostream& operator<<(ostream& s, const Vector& v)
+{
+  s << ' ' << v.v[0] << ' ' << v.v[1] << ' ';
+  return s;
+}
+
+istream& operator>>(istream& s, Vector& v)
+{
+  s >> v.v[0] >> v.v[1];
+  return s;
+}
+
+Vector cross(Vector& a, Vector& b)
+{
+  return Vector(a[1]*b[2]-b[1]*a[2], a[2]*b[0]-b[2]*a[0], a[0]*b[1]-b[0]*a[1]);
+}
+
+// the following are not valid for 2D graphics:
+
+Vector& Vector::div(double f)
+{
+  v[0]/=f;
+  v[1]/=f;
+  v[2]/=f;
+  return *this;
+}
+
+Vector& Vector::minus(const Vector& a)
+{
+  v[0]-=a.v[0];
+  v[1]-=a.v[1];
+  v[2]-=a.v[2];
+  return *this;
+}
+
+Vector mult(const Vector& a, double f)
+{
+  Vector r = a;
+  r.v[0]*= f;
+  r.v[1]*= f;
+  r.v[2]*= f;
+  return r;
+}
+
+Vector& Vector::operator*=(const Vector& b)
+{
+  v[0]*=b.v[0];
+  v[1]*=b.v[1];
+  v[2]=1;
+  return *this;
+}
+
+Vector& Vector::operator/=(const Vector& b)
+{
+  v[0]/=b.v[0];
+  v[1]/=b.v[1];
+  v[2]=1;
+  return *this;
+}
+
+// Vertex
+
+Vertex::Vertex()
+{
+  next_ = NULL;
+  previous_ = NULL;
+}
+
+Vertex::Vertex(double x, double y)
+{
+  vector = Vector(x,y);
+  next_ = NULL;
+  previous_ = NULL;
+}
+
+Vertex::Vertex(const Vector& a)
+{
+  vector = a;
+  next_ = NULL;
+  previous_ = NULL;
+}
+
+ostream& operator<<(ostream& s, const Vertex& v)
+{
+  s << v.vector;
+  return s;
+}
+
+// Matrix
+
+Matrix::Matrix()
+{
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      m[i][j] = (i==j ? 1 : 0);
+}
+
+Matrix::Matrix(const Matrix& a)
+{
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      m[i][j] = a.m[i][j];
+}
+
+Matrix& Matrix::operator=(const Matrix& a)
+{
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      m[i][j] = a.m[i][j];
+
+  return *this;
+}
+
+Matrix::Matrix(const Vector& v0, const Vector& v1, const Vector& v2)
+{
+  int i;
+  for (i=0; i<3; i++)
+    m[0][i] = v0.v[i];
+  for (i=0; i<3; i++)
+    m[1][i] = v1.v[i];
+  for (i=0; i<3; i++)
+    m[2][i] = v2.v[i];
+}
+
+Matrix::Matrix(double a, double b, double c, double d, double e, double f)
+{
+  m[0][0]=a;
+  m[0][1]=b;
+  m[0][2]=0;
+  m[1][0]=c;
+  m[1][1]=d;
+  m[1][2]=0;
+  m[2][0]=e;
+  m[2][1]=f;
+  m[2][2]=1;
+}
+
+Matrix& Matrix::identity()
+{
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      m[i][j] = (i==j ? 1 : 0);
+
+  return *this;
+}
+
+int Matrix::isIdentity()
+{
+  return ((m[0][0]==1) && (m[0][1]==0) && (m[0][2]==0) &&
+	  (m[1][0]==0) && (m[1][1]==1) && (m[1][2]==0) &&
+	  (m[2][0]==0) && (m[2][1]==0) && (m[2][2]==1));
+}
+
+Matrix Matrix::abs()
+{
+  Matrix r;
+  r.m[0][0] = fabs(m[0][0]);
+  r.m[1][1] = fabs(m[1][1]);
+  return r;
+}
+
+Matrix& Matrix::operator*=(double a)
+{
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      m[i][j] *= a;
+
+  return *this;
+}
+
+Matrix& Matrix::operator/=(double a)
+{
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      m[i][j] /= a;
+
+  return *this;
+}
+
+Matrix& Matrix::operator*=(const Matrix& a)
+{
+  Matrix r;
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      r.m[i][j] = m[i][0]*a.m[0][j] + m[i][1]*a.m[1][j] + m[i][2]*a.m[2][j];
+
+  return *this=r;
+}
+
+Matrix Matrix::invert()
+{
+  // Method of Row Reduction-- Calculus and Analylic Geometry, A-13
+
+  Matrix r; // identy matrix
+
+  Vector m0 = m[0];
+  Vector m1 = m[1];
+  Vector m2 = m[2];
+
+  Vector r0 = r[0];
+  Vector r1 = r[1];
+  Vector r2 = r[2];
+
+  // check for 0 in m[0], swap rows 0 and 1 if so.
+  // this seems to be a problem only in the case of rotation of 90/270 degress
+
+  if ((m0[0]==0) || 
+      (m0[0]>0 && m0[0]<DBL_EPSILON) || 
+      (m0[0]<0 && m0[0]>-DBL_EPSILON)) {
+
+    Vector tm0 = m0;
+    Vector tr0 = r0;
+    m0 = m1;
+    r0 = r1;
+    m1 = tm0;
+    r1 = tr0;
+  }
+
+  // scale first row by 1/m0[0]
+
+  r0.div(m0[0]);
+  m0.div(m0[0]);
+
+  // zero out m1[0] and m2[0]
+
+  r1.minus(mult(r0,m1[0]));
+  m1.minus(mult(m0,m1[0]));
+  r2.minus(mult(r0,m2[0]));
+  m2.minus(mult(m0,m2[0]));
+  
+  // scale second row by 1/m1[1]
+
+  r1.div(m1[1]);
+  m1.div(m1[1]);
+
+  // zero out m0[1] and m2[1]
+
+  r0.minus(mult(r1,m0[1]));
+  m0.minus(mult(m1,m0[1]));
+  r2.minus(mult(r1,m2[1]));
+  m2.minus(mult(m1,m2[1]));
+
+  // scale third row by 1/m[2]
+
+  r2.div(m2[2]);
+  m2.div(m2[2]);
+
+  // and zero out m0[2] and m1[2]
+
+  r0.minus(mult(r2,m0[2]));
+  m0.minus(mult(m2,m0[2]));
+  r1.minus(mult(r2,m1[2]));
+  m2.minus(mult(m2,m1[2]));
+
+  return Matrix(r0,r1,r2);
+}
+
+Matrix Matrix::invert2()
+{
+  Matrix cc = this->cofactor();
+  Matrix aa = cc.adjoint();
+  double det = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + m[0][2]*aa.m[2][0];
+
+  Matrix rr;
+  for (int ii=0; ii<3; ii++ )
+    for (int jj=0; jj<3; jj++)
+      rr.m[ii][jj] = aa.m[ii][jj]/det;
+
+  return rr;
+}
+
+Matrix Matrix::cofactor()
+{
+  Matrix rr;
+
+  rr.m[0][0] = +(m[1][1]*m[2][2]-m[1][2]*m[2][1]);
+  rr.m[0][1] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]);
+  rr.m[0][2] = +(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
+
+  rr.m[1][0] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]);
+  rr.m[1][1] = +(m[0][0]*m[2][2]-m[0][2]*m[2][0]);
+  rr.m[1][2] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]);
+
+  rr.m[2][0] = +(m[0][1]*m[1][2]-m[0][2]*m[1][1]);
+  rr.m[2][1] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]);
+  rr.m[2][2] = +(m[0][0]*m[1][1]-m[0][1]*m[1][0]);
+
+  return rr;
+}
+
+Matrix Matrix::adjoint()
+{
+  Matrix rr;
+  for (int ii=0; ii<3; ii++)
+    for (int jj=0; jj<3; jj++)
+      rr.m[jj][ii] = m[ii][jj];
+
+  return rr;
+}
+
+Matrix operator*(const Matrix& a, double b)
+{
+  Matrix r=a;
+  return r*=b;
+}
+
+Matrix operator/(const Matrix& a, double b)
+{
+  Matrix r=a;
+  return r/=b;
+}
+
+Matrix operator*(const Matrix& a, const Matrix& b)
+{
+  Matrix r=a;
+  return r*=b;
+}
+
+Vector Matrix::operator[](int i)
+{
+  return Vector(m[i]);
+}
+
+ostream& operator<<(ostream& s, const Matrix& m)
+{
+  s << ' ';
+  for (int i=0; i<3; i++)
+    for (int j=0; j<2; j++)
+      s << m.m[i][j] << ' ';
+
+  return s;
+}
+
+istream& operator>>(istream& s, Matrix& m)
+{
+  for (int i=0; i<3; i++ )
+    for (int j=0; j<2; j++)
+      s >> m.m[i][j];
+
+  return s;
+}
+
+// Translate
+
+Translate::Translate(double x, double y) : Matrix()
+{
+  m[2][0]=x;
+  m[2][1]=y;
+}
+
+Translate::Translate(const Vector& v) : Matrix()
+{
+  m[2][0]=v.v[0];
+  m[2][1]=v.v[1];
+}
+
+Translate::Translate(const Matrix& a) : Matrix()
+{
+  m[2][0] = a.m[2][0];
+  m[2][1] = a.m[2][1];
+}
+
+Translate& Translate::operator=(const Matrix& a)
+{
+  m[2][0] = a.m[2][0];
+  m[2][1] = a.m[2][1];
+  return *this;
+}
+
+ostream& operator<<(ostream& s, const Translate& m)
+{
+  s << ' ' << m.m[2][0] << ' ' << m.m[2][1] << ' ';
+  return s;
+}
+
+istream& operator>>(istream& s, Translate& m)
+{
+  s >> m.m[2][0] >> m.m[2][1];
+  return s;
+}
+
+// Rotate
+
+Rotate::Rotate(double a) : Matrix()
+{
+  // note: signs reverse for X-Windows (origin is upper left)
+
+  m[0][0] = cos(a);
+  m[0][1] = -sin(a);
+  m[1][0] = sin(a);
+  m[1][1] = cos(a);
+
+  // this fixes a problem with numbers too small and tring to invert the matrix
+
+  if ((m[0][0]>0 && m[0][0]<DBL_EPSILON) || 
+      (m[0][0]<0 && m[0][0]>-DBL_EPSILON))
+    m[0][0] = 0;
+
+  if ((m[0][1]>0 && m[0][1]<DBL_EPSILON) ||
+      (m[0][1]<0 && m[0][1]>-DBL_EPSILON))
+    m[0][1] = 0;
+
+  if ((m[1][0]>0 && m[1][0]<DBL_EPSILON) ||
+      (m[1][0]<0 && m[1][0]>-DBL_EPSILON))
+    m[1][0] = 0;
+
+  if ((m[1][1]>0 && m[1][1]<DBL_EPSILON) ||
+      (m[1][1]<0 && m[1][1]>-DBL_EPSILON))
+    m[1][1] = 0;
+
+}
+
+Rotate::Rotate(double a, double b, double c, double d) : Matrix()
+{
+  m[0][0] = a;
+  m[0][1] = b;
+  m[1][0] = c;
+  m[1][1] = d;
+}
+
+Rotate::Rotate(const Matrix& a) : Matrix()
+{
+  for (int i=0; i<2; i++)
+    for (int j=0; j<2; j++)
+      m[i][j] = a.m[i][j];
+}
+
+Rotate& Rotate::operator=(const Matrix& a)
+{
+  for (int i=0; i<2; i++)
+    for (int j=0; j<2; j++)
+      m[i][j] = a.m[i][j];
+
+  return *this;
+}
+
+ostream& operator<<(ostream& s, const Rotate& m)
+{
+  s << ' ' << m.m[0][0] << ' ' << m.m[0][1]
+    << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' ';
+  return s;
+}
+
+istream& operator>>(istream& s, Rotate& m)
+{
+  s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1];
+  return s;
+}
+
+// Scale
+
+Scale::Scale(double a) : Matrix()
+{
+  m[0][0]=a;
+  m[1][1]=a;
+}
+
+Scale::Scale(double a, double b) : Matrix()
+{
+  m[0][0]=a;
+  m[1][1]=b;
+}
+
+Scale::Scale(const Vector& v) : Matrix()
+{
+  m[0][0]=v.v[0];
+  m[1][1]=v.v[1];
+}
+
+Scale::Scale(const Matrix& a) : Matrix()
+{
+  m[0][0] = a.m[0][0];
+  m[1][1] = a.m[1][1];
+}
+
+Scale& Scale::operator=(const Matrix& a)
+{
+  m[0][0] = a.m[0][0];
+  m[1][1] = a.m[1][1];
+  return *this;
+}
+
+ostream& operator<<(ostream& s, const Scale& m)
+{
+  s << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' ';
+  return s;
+}
+
+istream& operator>>(istream& s, Scale& m)
+{
+  s >> m.m[0][0] >> m.m[1][1];
+  return s;
+}
+
+// Flip
+
+FlipX::FlipX() : Matrix()
+{
+  m[0][0] = -1;
+}
+
+FlipY::FlipY() : Matrix()
+{
+  m[1][1] = -1;
+}
+
+FlipXY::FlipXY() : Matrix()
+{
+  m[0][0] = -1;
+  m[1][1] = -1;
+}
+
+// BBox
+
+BBox::BBox(double a, double b, double c, double d)
+{
+  // we want a 'positive' box
+
+  ll.v[0] = a < c ? a : c;
+  ll.v[1] = b < d ? b : d;
+  ur.v[0] = a < c ? c : a;
+  ur.v[1] = b < d ? d : b;
+}
+
+BBox::BBox(double w, double h)
+{
+  ll.v[0] = 0;
+  ll.v[1] = 0;
+  ur.v[0] = w;
+  ur.v[1] = h;
+}
+
+BBox::BBox(const Vector& l, const Vector& h)
+{
+  // we want a 'positive' box
+
+  ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0];
+  ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1];
+  ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0];
+  ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1];
+}
+
+BBox::BBox(const Vector& v)
+{
+  ll = v;
+  ur = v;
+}
+
+BBox& BBox::operator+=(const Vector& v)
+{
+  ll+=v;
+  ur+=v;
+  return *this;
+}
+
+BBox operator+(const BBox& b, const Vector& v)
+{
+  BBox r=b;
+  r+=v;
+  return r;
+}
+
+BBox& BBox::operator-=(const Vector& a)
+{
+  ll-=a;
+  ur-=a;
+  return *this;
+}
+
+BBox operator-(const BBox& b, const Vector& v)
+{
+  BBox r=b;
+  r-=v;
+  return r;
+}
+
+int BBox::isIn(const Vector& v) const
+{
+  if (v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || v.v[0] > ur.v[0] || v.v[1] > ur.v[1])
+    return 0;
+  else
+    return 1;
+}
+
+int BBox::isIn(const BBox& bb) const
+{
+  // return 0 if outside, > 0 if intersection
+  // = 4 if inside
+  BBox b = bb;
+  return isIn(b.ll) + isIn(b.ur) + isIn(b.ul()) + isIn(b.lr());
+}
+
+int BBox::isEmpty() const
+{
+  Vector v = ur-ll;
+  return (v[0]==0 && v[1]==0);
+}
+
+Vector BBox::center()
+{
+  return (ur-ll)/2 + ll;
+}
+
+Vector BBox::size()
+{
+  return ur - ll;
+}
+
+BBox& BBox::operator*=(const Matrix& m)
+{
+  ll *= m;
+  ur *= m;
+  return *this;
+}
+
+BBox operator*(const BBox& bb, const Matrix& m)
+{
+  BBox r = bb;
+  r*=m;
+  return r;
+}
+
+BBox& BBox::expand(double a)
+{
+  ll -= Vector(a,a);
+  ur += Vector(a,a);
+  return *this;
+}
+
+BBox& BBox::expand(const Vector& v)
+{
+  ll -= v;
+  ur += v;
+  return *this;
+}
+
+BBox& BBox::shrink(double a)
+{
+  ll += Vector(a,a);
+  ur -= Vector(a,a);
+  return *this;
+}
+
+BBox& BBox::shrink(const Vector& v)
+{
+  ll += v;
+  ur -= v;
+  return *this;
+}
+
+BBox& BBox::bound(const Vector& v)
+{
+  if (v.v[0] < ll[0])
+    ll[0] = v.v[0];
+  if (v.v[1] < ll[1])
+    ll[1] = v.v[1];
+    
+  if (v.v[0] > ur[0])
+    ur[0] = v.v[0];
+  if (v.v[1] > ur[1])
+    ur[1] = v.v[1];
+
+  return *this;
+}
+
+BBox& BBox::bound(BBox b)
+{
+  this->bound(b.ll);
+  this->bound(b.lr());
+  this->bound(b.ur);
+  this->bound(b.ul());
+
+  return *this;
+}
+
+BBox intersect(const BBox& a, const BBox& b)
+{
+  // test for obvious
+
+  int ab = a.isIn(b);
+  int ba = b.isIn(a);
+
+  // we are missing a case here! two bbox's that cross
+
+  if (ab==0 && ba == 0) // outside each other
+    return BBox();
+
+  if (ab == 4) // b is inside a
+    return b; 
+  if (ba == 4) // a is inside b
+    return a;
+
+  // else, there seems to be some overlap
+
+  BBox r;
+  r.ll.v[0] = (a.ll.v[0] > b.ll.v[0]) ? a.ll.v[0] : b.ll.v[0];
+  r.ll.v[1] = (a.ll.v[1] > b.ll.v[1]) ? a.ll.v[1] : b.ll.v[1];
+  r.ur.v[0] = (a.ur.v[0] < b.ur.v[0]) ? a.ur.v[0] : b.ur.v[0];
+  r.ur.v[1] = (a.ur.v[1] < b.ur.v[1]) ? a.ur.v[1] : b.ur.v[1];
+
+  return r;
+}
+
+ostream& operator<<(ostream& s, const BBox& b)
+{
+  s << b.ll << b.ur;
+  return s;
+}
+
+
diff --git a/tksao/vector/vectorold.h b/tksao/vector/vectorold.h
new file mode 100644
index 0000000..f41efa5
--- /dev/null
+++ b/tksao/vector/vectorold.h
@@ -0,0 +1,248 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#ifndef __vector_h__
+#define __vector_h__
+
+#include <math.h>
+#include <float.h>
+
+#include <iostream>
+using namespace std;
+
+#include <tk.h>
+
+class Vector;
+class Vertex;
+class Matrix;
+class Translate;
+class Rotate;
+class Scale;
+class FlipX;
+class FlipY;
+class FlipXY;
+class BBox;
+
+class Vector {
+  friend class Matrix;
+  friend class Translate;
+  friend class Scale;
+  friend class BBox;
+
+ public:
+  double v[3];
+
+ public:
+  Vector() {v[0]=0;v[1]=0;v[2]=1;} // constructor
+  Vector(double* f) {v[0]=f[0];v[1]=f[1];v[2]=f[2];} // constructor
+  Vector(double x, double y) {v[0]=x;v[1]=y;v[2]=1;} // constructor
+  Vector(double x, double y, double z) {v[0]=x;v[1]=y;v[2]=z;} // constructor
+  Vector(const Vector& a) {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];} // copy
+  Vector& operator=(const Vector& a) 
+    {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];return *this;}
+
+  double& operator[](int i) {return v[i];}
+  Vector abs();
+  double length();
+  double angle();
+  double avg();
+  double max();
+  double min();
+  Vector round();
+  Vector floor();
+  Vector ceil();
+  Vector invert();
+  Vector normalize();
+  Vector TkCanvasPs(Tk_Canvas);
+
+  friend Vector operator-(const Vector&);                // unary minus
+  friend Vector operator*(const Vector&, const Matrix&); // vector/matrix mult
+  friend double operator*(const Vector&, const Vector&); // dot product
+
+  friend ostream& operator<<(ostream&, const Vector&);
+  friend istream& operator>>(istream&, Vector&);
+
+  Vector& operator+=(const Vector&);   // addition
+  Vector& operator-=(const Vector&);   // subtraction
+  Vector& operator*=(double);          // scalar multipy
+  Vector& operator/=(double);          // scalar division
+  Vector& operator*=(const Matrix&);   // vector multiply
+
+  friend BBox intersect(const BBox&, const BBox&);
+
+  // the following are not valid for 2D graphics:
+  Vector& div(double);
+  Vector& minus(const Vector&);
+  friend Vector mult(const Vector&, double);
+  Vector& operator*=(const Vector&);
+  Vector& operator/=(const Vector&);
+};
+
+Vector operator+(const Vector&, const Vector&);    // addition
+Vector operator-(const Vector&, const Vector&);    // subtration
+Vector operator*(const Vector&, double);           // scalar multiply
+Vector operator/(const Vector&, double);           // scalar division
+
+Vector cross(Vector&, Vector&);
+
+class Vertex {
+public:
+  Vector vector;
+
+private:
+  Vertex* next_;
+  Vertex* previous_;
+
+public:
+  Vertex();
+  Vertex(double, double);
+  Vertex(const Vector&);
+  
+  Vertex* next() {return next_;}
+  Vertex* previous() {return previous_;}
+
+  void setNext(Vertex* v) {next_ = v;}
+  void setPrevious(Vertex* v) {previous_ = v;}
+
+  friend ostream& operator<<(ostream&, const Vertex&);
+};
+
+class Matrix {
+  friend class Vector;
+  friend class Translate;
+  friend class Rotate;
+  friend class Scale;
+
+protected:
+  double m[3][3];
+
+public:
+  Matrix();                               // constructor
+  Matrix(const Vector&, const Vector&, const Vector&);      // constructor
+  Matrix(double, double, double, double, double, double);
+  Matrix(const Matrix&);                  // copy constructor
+  Matrix& operator=(const Matrix&);       // assignment
+  
+  Vector operator[](int);                 // return row
+  double matrix(int i, int j) {return m[i][j];} // return element
+  Matrix& operator*=(double);             // scalar multiply
+  Matrix& operator/=(double);             // scalar division
+  Matrix& operator*=(const Matrix&);      // matrix multiply
+
+  Matrix invert();                        // matrix inverse
+
+  Matrix invert2();
+  Matrix cofactor();
+  Matrix adjoint();
+
+  Matrix abs();                           // returns abs with no translation
+  Matrix& identity();                        // sets to identity matrix;
+  int isIdentity();
+
+  double* mm() {return (double*)m;}
+
+  friend Vector operator*(const Vector&, const Matrix&); // vector/matrix mult
+  friend ostream& operator<<(ostream&, const Matrix&);
+  friend istream& operator>>(istream&, Matrix&);
+};
+
+Matrix operator*(const Matrix&, double);           // scalar multiply
+Matrix operator/(const Matrix&, double);           // scalar division
+Matrix operator*(const Matrix&, const Matrix&);    // matrix multiply
+
+class Translate : public Matrix {
+public:
+  Translate() {};                      // constructor
+  Translate(double, double);           // constructor
+  Translate(const Vector&);            // constructor
+  Translate(const Matrix&);            // copy constructor
+  Translate& operator=(const Matrix&); // assignment
+
+  friend ostream& operator<<(ostream&, const Translate&);
+  friend istream& operator>>(istream&, Translate&);
+};
+
+class Rotate : public Matrix {
+public:
+  Rotate() {};                         // constructor
+  Rotate(double);                      // constructor
+  Rotate(double, double, double, double); // constructor
+  Rotate(const Matrix&);               // copy constructor
+  Rotate& operator=(const Matrix&);    // assignment
+
+  friend ostream& operator<<(ostream&, const Rotate&);
+  friend istream& operator>>(istream&, Rotate&);
+};
+
+class Scale : public Matrix {
+public:
+  Scale() {};                          // constructor
+  Scale(double);                       // constructor
+  Scale(double, double);               // constructor
+  Scale(const Vector&);                // constructor
+  Scale(const Matrix&);                // copy constructor
+  Scale& operator=(const Matrix&);     // assignment
+
+  friend ostream& operator<<(ostream&, const Scale&);
+  friend istream& operator>>(istream&, Scale&);
+};
+
+class FlipX : public Matrix {
+public:
+  FlipX();                             // constructor
+};
+
+class FlipY : public Matrix {
+public:
+  FlipY();                             // constructor
+};
+
+class FlipXY : public Matrix {
+public:
+  FlipXY();                            // constructor
+};
+
+class BBox {
+public:
+  Vector ll;
+  Vector ur;
+
+public:
+  BBox() {}
+  BBox(const Vector&, const Vector&);
+  BBox(const Vector&);
+  BBox(double, double);
+  BBox(double, double, double, double);
+
+  BBox& operator+=(const Vector&);   // addition
+  BBox& operator-=(const Vector&);   // subtraction
+  BBox& operator*=(const Matrix&);   // multiplication
+
+  int isIn(const Vector&) const;
+  int isIn(const BBox&) const;
+  int isEmpty() const;
+  Vector center();
+  Vector size();
+  Vector lr() {return Vector(ur[0],ll[1]);}
+  Vector ul() {return Vector(ll[0],ur[1]);}
+  BBox& expand(double);
+  BBox& expand(const Vector&);
+  BBox& shrink(double);
+  BBox& shrink(const Vector&);
+  BBox& bound(const Vector&);
+  BBox& bound(BBox);
+
+  friend ostream& operator<<(ostream&, const BBox&);
+};
+
+BBox operator+(const BBox&, const Vector&);    // addition
+BBox operator-(const BBox&, const Vector&);    // subtraction
+BBox operator*(const BBox&, const Matrix&);
+BBox intersect(const BBox&, const BBox&);
+
+#endif
+
+
+
+