mv vector out from tksao

7 files changed, 1929 insertions, 0 deletions

diff --git a/vector/fuzzy.h b/vector/fuzzy.h
new file mode 100644
index 0000000..d7f5f1b
--- /dev/null
+++ b/vector/fuzzy.h
@@ -0,0 +1,38 @@
+// Copyright (C) 1999-2018
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#ifndef __fuzzy_h__
+#define __fuzzy_h__
+
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+using namespace std;
+
+#include <float.h>
+
+inline void tzero(double* ff, const double epsilon= DBL_EPSILON)
+{if (*ff>=-epsilon && *ff<=epsilon) *ff = 0;}
+
+inline bool teq(const double f1, const double f2,
+		const double epsilon= DBL_EPSILON)
+{return f1-f2 >= -epsilon && f1-f2 <= epsilon;}
+
+inline bool tlt(const double f1, const double f2,
+		const double epsilon= DBL_EPSILON)
+{return f1-f2 < -epsilon;}
+
+inline bool tle(const double f1, const double f2,
+		const double epsilon= DBL_EPSILON)
+{return f1-f2 <= -epsilon;  }
+
+inline bool tgt(const double f1, const double f2, 
+		const double epsilon= DBL_EPSILON)
+{return f1-f2 > epsilon;}
+
+inline bool tge(const double f1, const double f2, 
+		const double epsilon= DBL_EPSILON)
+{return f1-f2 >= epsilon;}
+
+#endif
diff --git a/vector/vector.C b/vector/vector.C
new file mode 100644
index 0000000..88949f1
--- /dev/null
+++ b/vector/vector.C
@@ -0,0 +1,408 @@
+// Copyright (C) 1999-2018
+// 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
+
+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;
+}
+
+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/vector/vector.h b/vector/vector.h
new file mode 100644
index 0000000..783336e
--- /dev/null
+++ b/vector/vector.h
@@ -0,0 +1,361 @@
+// Copyright (C) 1999-2018
+// 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&);
+
+  const double& operator[](int i) const {return v[i];} // return element
+  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 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/vector/vector3d.C b/vector/vector3d.C
new file mode 100644
index 0000000..78575ae
--- /dev/null
+++ b/vector/vector3d.C
@@ -0,0 +1,468 @@
+// Copyright (C) 1999-2018
+// 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;
+}
+
+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/vector/vector3d.h b/vector/vector3d.h
new file mode 100644
index 0000000..3d77d86
--- /dev/null
+++ b/vector/vector3d.h
@@ -0,0 +1,395 @@
+// Copyright (C) 1999-2018
+// 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;}
+
+  const double& operator[](int i) const {return v[i];} // return element
+  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;}
+};
+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/vector/vectorstr.C b/vector/vectorstr.C
new file mode 100644
index 0000000..7505fc6
--- /dev/null
+++ b/vector/vectorstr.C
@@ -0,0 +1,205 @@
+// Copyright (C) 1999-2018
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+using namespace std;
+
+#include "vectorstr.h"
+#include "vector.h"
+#include "vector3d.h"
+
+static char* dupstr(const char* str)
+{
+  char* copy;
+  if (str) {
+    copy=new char[strlen(str)+1];
+    strcpy(copy,str);
+  }
+  else
+    copy=NULL;
+
+  return copy;
+}
+
+// VectorStr
+
+VectorStr::~VectorStr()
+{
+  if (c[0])
+    delete [] c[0];
+  if (c[1])
+    delete [] c[1];
+}
+
+VectorStr::VectorStr(const char* aa, const char* bb)
+{
+  c[0] = dupstr(aa);
+  c[1] = dupstr(bb);
+}
+
+VectorStr::VectorStr(const VectorStr& vv)
+{
+  c[0] = dupstr(vv.c[0]);
+  c[1] = dupstr(vv.c[1]);
+}
+
+VectorStr::VectorStr(const Vector& vv)
+{
+  ostringstream str0;
+  str0 << vv[0];
+  c[0] = dupstr(str0.str().c_str());
+
+  ostringstream str1;
+  str1 << vv[1];
+  c[1] = dupstr(str1.str().c_str());
+}
+
+VectorStr& VectorStr::operator=(const VectorStr& vv)
+{
+  if (c[0])
+    delete [] c[0];
+  c[0]=dupstr(vv.c[0]);
+
+  if (c[1])
+    delete [] c[1];
+  c[1]=dupstr(vv.c[1]);
+
+  return *this;
+}
+
+VectorStr& VectorStr::operator=(const Vector& vv)
+{
+  if (c[0])
+    delete [] c[0];
+  ostringstream str0;
+  str0 << vv[0];
+  c[0] = dupstr(str0.str().c_str());
+
+  if (c[1])
+    delete [] c[1];
+  ostringstream str1;
+  str1 << vv[1];
+  c[1] = dupstr(str1.str().c_str());
+
+  return *this;
+}
+
+ostream& operator<<(ostream& os, const VectorStr& vv)
+{
+  unsigned char sep = (unsigned char)os.iword(Vector::separator);
+  if (!sep)
+    sep = ' ';
+
+  unsigned char unit = (unsigned char)os.iword(Vector::unit);
+  if (!unit)
+    os << vv.c[0] << sep << vv.c[1];
+  else
+    os << vv.c[0] << unit << sep << vv.c[1] << unit;
+
+  // reset unit
+  os.iword(Vector::unit) = '\0';
+
+  return os;
+}
+
+// VectorStr3d
+
+VectorStr3d::~VectorStr3d()
+{
+  if (c[0])
+    delete [] c[0];
+  if (c[1])
+    delete [] c[1];
+  if (c[2])
+    delete [] c[2];
+}
+
+VectorStr3d::VectorStr3d(const char* aa, const char* bb, const char* cc)
+{
+  c[0] = dupstr(aa);
+  c[1] = dupstr(bb);
+  c[2] = dupstr(cc);
+}
+
+VectorStr3d::VectorStr3d(const VectorStr3d& vv)
+{
+  c[0] = dupstr(vv.c[0]);
+  c[1] = dupstr(vv.c[1]);
+  c[2] = dupstr(vv.c[2]);
+}
+
+VectorStr3d::VectorStr3d(const Vector3d& vv)
+{
+  ostringstream str0;
+  str0 << vv[0];
+  c[0] = dupstr(str0.str().c_str());
+
+  ostringstream str1;
+  str1 << vv[1];
+  c[1] = dupstr(str1.str().c_str());
+
+  ostringstream str2;
+  str2 << vv[2];
+  c[2] = dupstr(str2.str().c_str());
+}
+
+VectorStr3d& VectorStr3d::operator=(const VectorStr3d& vv)
+{
+  if (c[0])
+    delete [] c[0];
+  c[0]=dupstr(vv.c[0]);
+
+  if (c[1])
+    delete [] c[1];
+  c[1]=dupstr(vv.c[1]);
+
+  if (c[2])
+    delete [] c[2];
+  c[2]=dupstr(vv.c[2]);
+
+  return *this;
+}
+
+VectorStr3d& VectorStr3d::operator=(const Vector3d& vv)
+{
+  if (c[0])
+    delete [] c[0];
+  ostringstream str0;
+  str0 << vv[0];
+  c[0] = dupstr(str0.str().c_str());
+
+  if (c[1])
+    delete [] c[1];
+  ostringstream str1;
+  str1 << vv[1];
+  c[1] = dupstr(str1.str().c_str());
+
+  if (c[2])
+    delete [] c[2];
+  ostringstream str2;
+  str2 << vv[2];
+  c[2] = dupstr(str2.str().c_str());
+
+  return *this;
+}
+
+ostream& operator<<(ostream& os, const VectorStr3d& vv)
+{
+  unsigned char sep = (unsigned char)os.iword(Vector::separator);
+  if (!sep)
+    sep = ' ';
+
+  unsigned char unit = (unsigned char)os.iword(Vector::unit);
+  if (!unit)
+    os << vv.c[0] << sep << vv.c[1] << sep << vv.c[2];
+  else
+    os << vv.c[0] << unit << sep << vv.c[1] << unit << sep << vv.c[2] << unit;
+
+  // reset unit
+  os.iword(Vector::unit) = '\0';
+
+  return os;
+}
diff --git a/vector/vectorstr.h b/vector/vectorstr.h
new file mode 100644
index 0000000..71181fb
--- /dev/null
+++ b/vector/vectorstr.h
@@ -0,0 +1,54 @@
+// Copyright (C) 1999-2018
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#ifndef __vectorstr_h__
+#define __vectorstr_h__
+
+#include "vector.h"
+
+#include <iostream>
+using namespace std;
+
+class VectorStr {
+ public:
+  char* c[2];
+
+ public:
+  VectorStr() {c[0]=NULL; c[1]=NULL;}
+  ~VectorStr();
+  VectorStr(const char*, const char*);
+  VectorStr(const VectorStr&);
+  VectorStr(const Vector&);
+  VectorStr& operator=(const VectorStr&);
+  VectorStr& operator=(const Vector&);
+
+  const char* operator[](int i) const {return c[i];} // return element
+  char** cc() {return c;} // return vector
+};
+
+ostream& operator<<(ostream&, const VectorStr&);
+
+class VectorStr3d {
+ public:
+  char* c[3];
+
+ public:
+  VectorStr3d() {c[0]=NULL; c[1]=NULL; c[2]=NULL;}
+  ~VectorStr3d();
+  VectorStr3d(const char*, const char*, const char*);
+  VectorStr3d(const VectorStr3d&);
+  VectorStr3d(const Vector3d&);
+  VectorStr3d& operator=(const VectorStr3d&);
+  VectorStr3d& operator=(const Vector3d&);
+
+  const char* operator[](int i) const {return c[i];} // return element
+  char** cc() {return c;} // return vector
+};
+
+ostream& operator<<(ostream&, const VectorStr3d&);
+#endif
+
+
+
+