diff --git a/src/bsp.cpp b/src/bsp.cpp
index 527c47c..7ff7191 100644
--- a/src/bsp.cpp
+++ b/src/bsp.cpp
@@ -36,7 +36,7 @@ SBsp3 *SBsp3::FromMesh(SMesh *m) {
     return bsp3;
 }
 
-Vector SBsp3::IntersectionWith(Vector a, Vector b) {
+Vector SBsp3::IntersectionWith(Vector a, Vector b) const {
     double da = a.Dot(n) - d;
     double db = b.Dot(n) - d;
     ssassert(da*db < 0, "Expected segment to intersect BSP node");
@@ -401,7 +401,7 @@ void SBsp3::Insert(STriangle *tr, SMesh *instead) {
     return;
 }
 
-void SBsp3::GenerateInPaintOrder(SMesh *m) {
+void SBsp3::GenerateInPaintOrder(SMesh *m) const {
 
     // Doesn't matter which branch we take if the normal has zero z
     // component, so don't need a separate case for that.
@@ -411,7 +411,7 @@ void SBsp3::GenerateInPaintOrder(SMesh *m) {
         if(neg) neg->GenerateInPaintOrder(m);
     }
 
-    SBsp3 *flip = this;
+    const SBsp3 *flip = this;
     while(flip) {
         m->AddTriangle(&(flip->tri));
         flip = flip->more;
@@ -424,7 +424,7 @@ void SBsp3::GenerateInPaintOrder(SMesh *m) {
     }
 }
 
-void SBsp3::DebugDraw() {
+void SBsp3::DebugDraw() const {
 
     if(pos) pos->DebugDraw();
     Vector norm = tri.Normal();
@@ -468,7 +468,7 @@ void SBsp3::DebugDraw() {
 
 /////////////////////////////////
 
-Vector SBsp2::IntersectionWith(Vector a, Vector b) {
+Vector SBsp2::IntersectionWith(Vector a, Vector b) const {
     double da = a.Dot(no) - d;
     double db = b.Dot(no) - d;
     ssassert(da*db < 0, "Expected segment to intersect BSP node");
@@ -665,7 +665,7 @@ void SBsp2::InsertTriangle(STriangle *tr, SMesh *m, SBsp3 *bsp3) {
     return;
 }
 
-void SBsp2::DebugDraw(Vector n, double d) {
+void SBsp2::DebugDraw(Vector n, double d) const {
     ssassert(fabs((edge.a).Dot(n) - d) < LENGTH_EPS, "Endpoint too close to BSP node plane");
     ssassert(fabs((edge.b).Dot(n) - d) < LENGTH_EPS, "Endpoint too close to BSP node plane");
 
diff --git a/src/constraint.cpp b/src/constraint.cpp
index 23bb126..59bbad4 100644
--- a/src/constraint.cpp
+++ b/src/constraint.cpp
@@ -6,7 +6,7 @@
 //-----------------------------------------------------------------------------
 #include "solvespace.h"
 
-std::string Constraint::DescriptionString() {
+std::string Constraint::DescriptionString() const {
     const char *s;
     switch(type) {
         case POINTS_COINCIDENT:     s = "pts-coincident"; break;
diff --git a/src/constrainteq.cpp b/src/constrainteq.cpp
index db5facb..c139cea 100644
--- a/src/constrainteq.cpp
+++ b/src/constrainteq.cpp
@@ -9,7 +9,7 @@
 
 const hConstraint ConstraintBase::NO_CONSTRAINT = { 0 };
 
-bool ConstraintBase::HasLabel() {
+bool ConstraintBase::HasLabel() const {
     switch(type) {
         case PT_LINE_DISTANCE:
         case PT_PLANE_DISTANCE:
@@ -193,7 +193,7 @@ void ConstraintBase::ModifyToSatisfy() {
     }
 }
 
-void ConstraintBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index)
+void ConstraintBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) const
 {
     Equation eq;
     eq.e = expr;
@@ -201,12 +201,12 @@ void ConstraintBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index)
     l->Add(&eq);
 }
 
-void ConstraintBase::Generate(IdList<Equation,hEquation> *l) {
+void ConstraintBase::Generate(IdList<Equation,hEquation> *l) const {
     if(!reference) {
         GenerateReal(l);
     }
 }
-void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) {
+void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) const {
     Expr *exA = Expr::From(valA);
 
     switch(type) {
diff --git a/src/drawconstraint.cpp b/src/drawconstraint.cpp
index df795bc..3cdd1f0 100644
--- a/src/drawconstraint.cpp
+++ b/src/drawconstraint.cpp
@@ -40,7 +40,7 @@ static void LineCallback(void *fndata, Vector a, Vector b)
     c->LineDrawOrGetDistance(a, b);
 }
 
-std::string Constraint::Label() {
+std::string Constraint::Label() const {
     std::string result;
     if(type == ANGLE) {
         if(valA == floor(valA)) {
@@ -104,7 +104,7 @@ void Constraint::DoLabel(Vector ref, Vector *labelPos, Vector gr, Vector gu) {
 
 
     if(dogd.drawing) {
-        ssglWriteTextRefCenter(s, th, ref, gr, gu, LineCallback, this);
+        ssglWriteTextRefCenter(s, th, ref, gr, gu, LineCallback, (void *)this);
     } else {
         double l = swidth/2 - sheight/2;
         l = max(l, 5/SS.GW.scale);
@@ -461,7 +461,7 @@ void Constraint::DoArcForAngle(Vector a0, Vector da, Vector b0, Vector db,
         Vector trans =
             (*ref).Plus(gu.ScaledBy(-1.5*ssglStrCapHeight(Style::DefaultTextHeight())));
         ssglWriteTextRefCenter("angle between skew lines", Style::DefaultTextHeight(),
-            trans, gr.WithMagnitude(px), gu.WithMagnitude(px), LineCallback, this);
+            trans, gr.WithMagnitude(px), gu.WithMagnitude(px), LineCallback, (void *)this);
     }
 }
 
@@ -921,7 +921,7 @@ void Constraint::DrawOrGetDistance(Vector *labelPos) {
 
             if(dogd.drawing) {
                 ssglWriteTextRefCenter("T", Style::DefaultTextHeight(),
-                    textAt, u, v, LineCallback, this);
+                    textAt, u, v, LineCallback, (void *)this);
             } else {
                 dogd.refp = textAt;
                 Point2d ref = SS.GW.ProjectPoint(dogd.refp);
@@ -1105,7 +1105,7 @@ s:
                                     (type == AT_MIDPOINT) ? "M" : NULL));
 
                     ssglWriteTextRefCenter(s, Style::DefaultTextHeight(),
-                        m.Plus(offset), r, u, LineCallback, this);
+                        m.Plus(offset), r, u, LineCallback, (void *)this);
                 } else {
                     dogd.refp = m.Plus(offset);
                     Point2d ref = SS.GW.ProjectPoint(dogd.refp);
@@ -1220,7 +1220,7 @@ void Constraint::GetEdges(SEdgeList *sel) {
     dogd.sel = NULL;
 }
 
-bool Constraint::IsStylable() {
+bool Constraint::IsStylable() const {
     if(type == COMMENT) return true;
     return false;
 }
@@ -1230,7 +1230,7 @@ hStyle Constraint::GetStyle() const {
     return { Style::CONSTRAINT };
 }
 
-bool Constraint::HasLabel() {
+bool Constraint::HasLabel() const {
     switch(type) {
         case COMMENT:
         case PT_PT_DISTANCE:
diff --git a/src/drawentity.cpp b/src/drawentity.cpp
index 875d45b..6de97ce 100644
--- a/src/drawentity.cpp
+++ b/src/drawentity.cpp
@@ -7,7 +7,7 @@
 //-----------------------------------------------------------------------------
 #include "solvespace.h"
 
-std::string Entity::DescriptionString() {
+std::string Entity::DescriptionString() const {
     if(h.isFromRequest()) {
         Request *r = SK.GetRequest(h.request());
         return r->DescriptionString();
@@ -137,7 +137,6 @@ void Entity::GenerateEdges(SEdgeList *el, bool includingConstruction) {
         }
         lv.Clear();
     }
-
 }
 
 SBezierList *Entity::GetOrGenerateBezierCurves() {
@@ -211,14 +210,14 @@ Vector Entity::GetReferencePos() {
     return dogd.refp;
 }
 
-bool Entity::IsStylable() {
+bool Entity::IsStylable() const {
     if(IsPoint()) return false;
     if(IsWorkplane()) return false;
     if(IsNormal()) return false;
     return true;
 }
 
-bool Entity::IsVisible() {
+bool Entity::IsVisible() const {
     Group *g = SK.GetGroup(group);
 
     if(g->h.v == Group::HGROUP_REFERENCES.v && IsNormal()) {
@@ -271,7 +270,7 @@ void Entity::CalculateNumerical(bool forExport) {
     }
 }
 
-bool Entity::PointIsFromReferences() {
+bool Entity::PointIsFromReferences() const {
     return h.request().IsFromReferences();
 }
 
@@ -280,7 +279,7 @@ bool Entity::PointIsFromReferences() {
 // routine for periodic splines (in a loop) or open splines (with specified
 // end tangents).
 //-----------------------------------------------------------------------------
-void Entity::ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) {
+void Entity::ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) const {
     static const int MAX_N = BandedMatrix::MAX_UNKNOWNS;
     int ep = extraPoints;
 
@@ -413,7 +412,7 @@ void Entity::ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) {
     }
 }
 
-void Entity::GenerateBezierCurves(SBezierList *sbl) {
+void Entity::GenerateBezierCurves(SBezierList *sbl) const {
     SBezier sb;
 
     int i = sbl->l.n;
diff --git a/src/dsc.h b/src/dsc.h
index 91ecf4b..0ab01fd 100644
--- a/src/dsc.h
+++ b/src/dsc.h
@@ -27,23 +27,23 @@ public:
     static Quaternion From(Vector u, Vector v);
     static Quaternion From(Vector axis, double dtheta);
 
-    Quaternion Plus(Quaternion b);
-    Quaternion Minus(Quaternion b);
-    Quaternion ScaledBy(double s);
-    double Magnitude();
-    Quaternion WithMagnitude(double s);
+    Quaternion Plus(Quaternion b) const;
+    Quaternion Minus(Quaternion b) const;
+    Quaternion ScaledBy(double s) const;
+    double Magnitude() const;
+    Quaternion WithMagnitude(double s) const;
 
     // Call a rotation matrix [ u' v' n' ]'; this returns the first and
     // second rows, where that matrix is generated by this quaternion
-    Vector RotationU();
-    Vector RotationV();
-    Vector RotationN();
-    Vector Rotate(Vector p);
+    Vector RotationU() const;
+    Vector RotationV() const;
+    Vector RotationN() const;
+    Vector Rotate(Vector p) const;
 
-    Quaternion ToThe(double p);
-    Quaternion Inverse();
-    Quaternion Times(Quaternion b);
-    Quaternion Mirror();
+    Quaternion ToThe(double p) const;
+    Quaternion Inverse() const;
+    Quaternion Times(Quaternion b) const;
+    Quaternion Mirror() const;
 };
 
 class Vector {
@@ -68,43 +68,43 @@ public:
                                          Vector pb, Vector db,
                                          double *ta, double *tb);
 
-    double Element(int i);
-    bool Equals(Vector v, double tol=LENGTH_EPS);
-    bool EqualsExactly(Vector v);
-    Vector Plus(Vector b);
-    Vector Minus(Vector b);
-    Vector Negated();
-    Vector Cross(Vector b);
-    double DirectionCosineWith(Vector b);
-    double Dot(Vector b);
-    Vector Normal(int which);
-    Vector RotatedAbout(Vector orig, Vector axis, double theta);
-    Vector RotatedAbout(Vector axis, double theta);
-    Vector DotInToCsys(Vector u, Vector v, Vector n);
-    Vector ScaleOutOfCsys(Vector u, Vector v, Vector n);
-    double DistanceToLine(Vector p0, Vector dp);
-    bool OnLineSegment(Vector a, Vector b, double tol=LENGTH_EPS);
-    Vector ClosestPointOnLine(Vector p0, Vector dp);
-    double Magnitude();
-    double MagSquared();
-    Vector WithMagnitude(double s);
-    Vector ScaledBy(double s);
-    Vector ProjectInto(hEntity wrkpl);
-    Vector ProjectVectorInto(hEntity wrkpl);
-    double DivPivoting(Vector delta);
-    Vector ClosestOrtho();
-    void MakeMaxMin(Vector *maxv, Vector *minv);
-    Vector ClampWithin(double minv, double maxv);
+    double Element(int i) const;
+    bool Equals(Vector v, double tol=LENGTH_EPS) const;
+    bool EqualsExactly(Vector v) const;
+    Vector Plus(Vector b) const;
+    Vector Minus(Vector b) const;
+    Vector Negated() const;
+    Vector Cross(Vector b) const;
+    double DirectionCosineWith(Vector b) const;
+    double Dot(Vector b) const;
+    Vector Normal(int which) const;
+    Vector RotatedAbout(Vector orig, Vector axis, double theta) const;
+    Vector RotatedAbout(Vector axis, double theta) const;
+    Vector DotInToCsys(Vector u, Vector v, Vector n) const;
+    Vector ScaleOutOfCsys(Vector u, Vector v, Vector n) const;
+    double DistanceToLine(Vector p0, Vector dp) const;
+    bool OnLineSegment(Vector a, Vector b, double tol=LENGTH_EPS) const;
+    Vector ClosestPointOnLine(Vector p0, Vector deltal) const;
+    double Magnitude() const;
+    double MagSquared() const;
+    Vector WithMagnitude(double s) const;
+    Vector ScaledBy(double s) const;
+    Vector ProjectInto(hEntity wrkpl) const;
+    Vector ProjectVectorInto(hEntity wrkpl) const;
+    double DivPivoting(Vector delta) const;
+    Vector ClosestOrtho() const;
+    void MakeMaxMin(Vector *maxv, Vector *minv) const;
+    Vector ClampWithin(double minv, double maxv) const;
     static bool BoundingBoxesDisjoint(Vector amax, Vector amin,
                                       Vector bmax, Vector bmin);
     static bool BoundingBoxIntersectsLine(Vector amax, Vector amin,
                                           Vector p0, Vector p1, bool segment);
-    bool OutsideAndNotOn(Vector maxv, Vector minv);
+    bool OutsideAndNotOn(Vector maxv, Vector minv) const;
     Vector InPerspective(Vector u, Vector v, Vector n,
-                         Vector origin, double cameraTan);
-    Point2d Project2d(Vector u, Vector v);
-    Point2d ProjectXy();
-    Vector4 Project4d();
+                         Vector origin, double cameraTan) const;
+    Point2d Project2d(Vector u, Vector v) const;
+    Point2d ProjectXy() const;
+    Vector4 Project4d() const;
 };
 
 class Vector4 {
@@ -115,10 +115,10 @@ public:
     static Vector4 From(double w, Vector v3);
     static Vector4 Blend(Vector4 a, Vector4 b, double t);
 
-    Vector4 Plus(Vector4 b);
-    Vector4 Minus(Vector4 b);
-    Vector4 ScaledBy(double s);
-    Vector PerspectiveProject();
+    Vector4 Plus(Vector4 b) const;
+    Vector4 Minus(Vector4 b) const;
+    Vector4 ScaledBy(double s) const;
+    Vector PerspectiveProject() const;
 };
 
 class Point2d {
@@ -165,12 +165,12 @@ public:
         }
     }
 
-    void Add(T *t) {
+    void Add(const T *t) {
         AllocForOneMore();
         new(&elem[n++]) T(*t);
     }
 
-    void AddToBeginning(T *t) {
+    void AddToBeginning(const T *t) {
         AllocForOneMore();
         new(&elem[n]) T();
         std::move_backward(elem, elem + 1, elem + n + 1);
@@ -181,11 +181,19 @@ public:
     T *First() {
         return (n == 0) ? NULL : &(elem[0]);
     }
+    const T *First() const {
+        return (n == 0) ? NULL : &(elem[0]);
+    }
     T *NextAfter(T *prev) {
         if(!prev) return NULL;
         if(prev - elem == (n - 1)) return NULL;
         return prev + 1;
     }
+    const T *NextAfter(const T *prev) const {
+        if(!prev) return NULL;
+        if(prev - elem == (n - 1)) return NULL;
+        return prev + 1;
+    }
 
     void ClearTags() {
         int i;
@@ -515,12 +523,12 @@ public:
 
     static BBox From(const Vector &p0, const Vector &p1);
 
-    Vector GetOrigin();
-    Vector GetExtents();
+    Vector GetOrigin() const;
+    Vector GetExtents() const;
 
     void Include(const Vector &v, double r = 0.0);
-    bool Overlaps(BBox &b1);
-    bool Contains(const Point2d &p);
+    bool Overlaps(const BBox &b1) const;
+    bool Contains(const Point2d &p) const;
 };
 
 #endif
diff --git a/src/entity.cpp b/src/entity.cpp
index 660628f..ebd8e85 100644
--- a/src/entity.cpp
+++ b/src/entity.cpp
@@ -10,7 +10,7 @@
 const hEntity  EntityBase::FREE_IN_3D = { 0 };
 const hEntity  EntityBase::NO_ENTITY = { 0 };
 
-bool EntityBase::HasVector() {
+bool EntityBase::HasVector() const {
     switch(type) {
         case LINE_SEGMENT:
         case NORMAL_IN_3D:
@@ -25,7 +25,7 @@ bool EntityBase::HasVector() {
     }
 }
 
-ExprVector EntityBase::VectorGetExprs() {
+ExprVector EntityBase::VectorGetExprs() const {
     switch(type) {
         case LINE_SEGMENT:
             return (SK.GetEntity(point[0])->PointGetExprs()).Minus(
@@ -42,7 +42,7 @@ ExprVector EntityBase::VectorGetExprs() {
     }
 }
 
-Vector EntityBase::VectorGetNum() {
+Vector EntityBase::VectorGetNum() const {
     switch(type) {
         case LINE_SEGMENT:
             return (SK.GetEntity(point[0])->PointGetNum()).Minus(
@@ -59,7 +59,7 @@ Vector EntityBase::VectorGetNum() {
     }
 }
 
-Vector EntityBase::VectorGetRefPoint() {
+Vector EntityBase::VectorGetRefPoint() const {
     switch(type) {
         case LINE_SEGMENT:
             return ((SK.GetEntity(point[0])->PointGetNum()).Plus(
@@ -76,7 +76,7 @@ Vector EntityBase::VectorGetRefPoint() {
     }
 }
 
-Vector EntityBase::VectorGetStartPoint() {
+Vector EntityBase::VectorGetStartPoint() const {
     switch(type) {
         case LINE_SEGMENT:
             return SK.GetEntity(point[1])->PointGetNum();
@@ -92,11 +92,11 @@ Vector EntityBase::VectorGetStartPoint() {
     }
 }
 
-bool EntityBase::IsCircle() {
+bool EntityBase::IsCircle() const {
     return (type == CIRCLE) || (type == ARC_OF_CIRCLE);
 }
 
-Expr *EntityBase::CircleGetRadiusExpr() {
+Expr *EntityBase::CircleGetRadiusExpr() const {
     if(type == CIRCLE) {
         return SK.GetEntity(distance)->DistanceGetExpr();
     } else if(type == ARC_OF_CIRCLE) {
@@ -104,7 +104,7 @@ Expr *EntityBase::CircleGetRadiusExpr() {
     } else ssassert(false, "Unexpected entity type");
 }
 
-double EntityBase::CircleGetRadiusNum() {
+double EntityBase::CircleGetRadiusNum() const {
     if(type == CIRCLE) {
         return SK.GetEntity(distance)->DistanceGetNum();
     } else if(type == ARC_OF_CIRCLE) {
@@ -114,7 +114,7 @@ double EntityBase::CircleGetRadiusNum() {
     } else ssassert(false, "Unexpected entity type");
 }
 
-void EntityBase::ArcGetAngles(double *thetaa, double *thetab, double *dtheta) {
+void EntityBase::ArcGetAngles(double *thetaa, double *thetab, double *dtheta) const {
     ssassert(type == ARC_OF_CIRCLE, "Unexpected entity type");
 
     Quaternion q = Normal()->NormalGetNum();
@@ -137,46 +137,46 @@ void EntityBase::ArcGetAngles(double *thetaa, double *thetab, double *dtheta) {
     while(*dtheta > (2*PI)) *dtheta -= 2*PI;
 }
 
-Vector EntityBase::CubicGetStartNum() {
+Vector EntityBase::CubicGetStartNum() const {
     return SK.GetEntity(point[0])->PointGetNum();
 }
-Vector EntityBase::CubicGetFinishNum() {
+Vector EntityBase::CubicGetFinishNum() const {
     return SK.GetEntity(point[3+extraPoints])->PointGetNum();
 }
-ExprVector EntityBase::CubicGetStartTangentExprs() {
+ExprVector EntityBase::CubicGetStartTangentExprs() const {
     ExprVector pon  = SK.GetEntity(point[0])->PointGetExprs(),
                poff = SK.GetEntity(point[1])->PointGetExprs();
     return (pon.Minus(poff));
 }
-ExprVector EntityBase::CubicGetFinishTangentExprs() {
+ExprVector EntityBase::CubicGetFinishTangentExprs() const {
     ExprVector pon  = SK.GetEntity(point[3+extraPoints])->PointGetExprs(),
                poff = SK.GetEntity(point[2+extraPoints])->PointGetExprs();
     return (pon.Minus(poff));
 }
-Vector EntityBase::CubicGetStartTangentNum() {
+Vector EntityBase::CubicGetStartTangentNum() const {
     Vector pon  = SK.GetEntity(point[0])->PointGetNum(),
            poff = SK.GetEntity(point[1])->PointGetNum();
     return (pon.Minus(poff));
 }
-Vector EntityBase::CubicGetFinishTangentNum() {
+Vector EntityBase::CubicGetFinishTangentNum() const {
     Vector pon  = SK.GetEntity(point[3+extraPoints])->PointGetNum(),
            poff = SK.GetEntity(point[2+extraPoints])->PointGetNum();
     return (pon.Minus(poff));
 }
 
-bool EntityBase::IsWorkplane() {
+bool EntityBase::IsWorkplane() const {
     return (type == WORKPLANE);
 }
 
-ExprVector EntityBase::WorkplaneGetOffsetExprs() {
+ExprVector EntityBase::WorkplaneGetOffsetExprs() const {
     return SK.GetEntity(point[0])->PointGetExprs();
 }
 
-Vector EntityBase::WorkplaneGetOffset() {
+Vector EntityBase::WorkplaneGetOffset() const {
     return SK.GetEntity(point[0])->PointGetNum();
 }
 
-void EntityBase::WorkplaneGetPlaneExprs(ExprVector *n, Expr **dn) {
+void EntityBase::WorkplaneGetPlaneExprs(ExprVector *n, Expr **dn) const {
     if(type == WORKPLANE) {
         *n = Normal()->NormalExprsN();
 
@@ -188,18 +188,18 @@ void EntityBase::WorkplaneGetPlaneExprs(ExprVector *n, Expr **dn) {
     } else ssassert(false, "Unexpected entity type");
 }
 
-bool EntityBase::IsDistance() {
+bool EntityBase::IsDistance() const {
     return (type == DISTANCE) ||
            (type == DISTANCE_N_COPY);
 }
-double EntityBase::DistanceGetNum() {
+double EntityBase::DistanceGetNum() const {
     if(type == DISTANCE) {
         return SK.GetParam(param[0])->val;
     } else if(type == DISTANCE_N_COPY) {
         return numDistance;
     } else ssassert(false, "Unexpected entity type");
 }
-Expr *EntityBase::DistanceGetExpr() {
+Expr *EntityBase::DistanceGetExpr() const {
     if(type == DISTANCE) {
         return Expr::From(param[0]);
     } else if(type == DISTANCE_N_COPY) {
@@ -214,11 +214,11 @@ void EntityBase::DistanceForceTo(double v) {
     } else ssassert(false, "Unexpected entity type");
 }
 
-EntityBase *EntityBase::Normal() {
+EntityBase *EntityBase::Normal() const {
     return SK.GetEntity(normal);
 }
 
-bool EntityBase::IsPoint() {
+bool EntityBase::IsPoint() const {
     switch(type) {
         case POINT_IN_3D:
         case POINT_IN_2D:
@@ -233,7 +233,7 @@ bool EntityBase::IsPoint() {
     }
 }
 
-bool EntityBase::IsNormal() {
+bool EntityBase::IsNormal() const {
     switch(type) {
         case NORMAL_IN_3D:
         case NORMAL_IN_2D:
@@ -246,7 +246,7 @@ bool EntityBase::IsNormal() {
     }
 }
 
-Quaternion EntityBase::NormalGetNum() {
+Quaternion EntityBase::NormalGetNum() const {
     Quaternion q;
     switch(type) {
         case NORMAL_IN_3D:
@@ -310,27 +310,27 @@ void EntityBase::NormalForceTo(Quaternion q) {
     }
 }
 
-Vector EntityBase::NormalU() {
+Vector EntityBase::NormalU() const {
     return NormalGetNum().RotationU();
 }
-Vector EntityBase::NormalV() {
+Vector EntityBase::NormalV() const {
     return NormalGetNum().RotationV();
 }
-Vector EntityBase::NormalN() {
+Vector EntityBase::NormalN() const {
     return NormalGetNum().RotationN();
 }
 
-ExprVector EntityBase::NormalExprsU() {
+ExprVector EntityBase::NormalExprsU() const {
     return NormalGetExprs().RotationU();
 }
-ExprVector EntityBase::NormalExprsV() {
+ExprVector EntityBase::NormalExprsV() const {
     return NormalGetExprs().RotationV();
 }
-ExprVector EntityBase::NormalExprsN() {
+ExprVector EntityBase::NormalExprsN() const {
     return NormalGetExprs().RotationN();
 }
 
-ExprQuaternion EntityBase::NormalGetExprs() {
+ExprQuaternion EntityBase::NormalGetExprs() const {
     ExprQuaternion q;
     switch(type) {
         case NORMAL_IN_3D:
@@ -430,7 +430,7 @@ void EntityBase::PointForceTo(Vector p) {
     }
 }
 
-Vector EntityBase::PointGetNum() {
+Vector EntityBase::PointGetNum() const {
     Vector p;
     switch(type) {
         case POINT_IN_3D:
@@ -479,7 +479,7 @@ Vector EntityBase::PointGetNum() {
     return p;
 }
 
-ExprVector EntityBase::PointGetExprs() {
+ExprVector EntityBase::PointGetExprs() const {
     ExprVector r;
     switch(type) {
         case POINT_IN_3D:
@@ -528,7 +528,7 @@ ExprVector EntityBase::PointGetExprs() {
     return r;
 }
 
-void EntityBase::PointGetExprsInWorkplane(hEntity wrkpl, Expr **u, Expr **v) {
+void EntityBase::PointGetExprsInWorkplane(hEntity wrkpl, Expr **u, Expr **v) const {
     if(type == POINT_IN_2D && workplane.v == wrkpl.v) {
         // They want our coordinates in the form that we've written them,
         // very nice.
@@ -559,7 +559,7 @@ void EntityBase::PointForceQuaternionTo(Quaternion q) {
     SK.GetParam(param[6])->val = q.vz;
 }
 
-Quaternion EntityBase::GetAxisAngleQuaternion(int param0) {
+Quaternion EntityBase::GetAxisAngleQuaternion(int param0) const {
     Quaternion q;
     double theta = timesApplied*SK.GetParam(param[param0+0])->val;
     double s = sin(theta), c = cos(theta);
@@ -570,7 +570,7 @@ Quaternion EntityBase::GetAxisAngleQuaternion(int param0) {
     return q;
 }
 
-ExprQuaternion EntityBase::GetAxisAngleQuaternionExprs(int param0) {
+ExprQuaternion EntityBase::GetAxisAngleQuaternionExprs(int param0) const {
     ExprQuaternion q;
 
     Expr *theta = Expr::From(timesApplied)->Times(
@@ -583,7 +583,7 @@ ExprQuaternion EntityBase::GetAxisAngleQuaternionExprs(int param0) {
     return q;
 }
 
-Quaternion EntityBase::PointGetQuaternion() {
+Quaternion EntityBase::PointGetQuaternion() const {
     Quaternion q;
 
     if(type == POINT_N_ROT_AA) {
@@ -595,7 +595,7 @@ Quaternion EntityBase::PointGetQuaternion() {
     return q;
 }
 
-bool EntityBase::IsFace() {
+bool EntityBase::IsFace() const {
     switch(type) {
         case FACE_NORMAL_PT:
         case FACE_XPROD:
@@ -608,7 +608,7 @@ bool EntityBase::IsFace() {
     }
 }
 
-ExprVector EntityBase::FaceGetNormalExprs() {
+ExprVector EntityBase::FaceGetNormalExprs() const {
     ExprVector r;
     if(type == FACE_NORMAL_PT) {
         Vector v = Vector::From(numNormal.vx, numNormal.vy, numNormal.vz);
@@ -637,7 +637,7 @@ ExprVector EntityBase::FaceGetNormalExprs() {
     return r;
 }
 
-Vector EntityBase::FaceGetNormalNum() {
+Vector EntityBase::FaceGetNormalNum() const {
     Vector r;
     if(type == FACE_NORMAL_PT) {
         r = Vector::From(numNormal.vx, numNormal.vy, numNormal.vz);
@@ -660,7 +660,7 @@ Vector EntityBase::FaceGetNormalNum() {
     return r.WithMagnitude(1);
 }
 
-ExprVector EntityBase::FaceGetPointExprs() {
+ExprVector EntityBase::FaceGetPointExprs() const {
     ExprVector r;
     if(type == FACE_NORMAL_PT) {
         r = SK.GetEntity(point[0])->PointGetExprs();
@@ -689,7 +689,7 @@ ExprVector EntityBase::FaceGetPointExprs() {
     return r;
 }
 
-Vector EntityBase::FaceGetPointNum() {
+Vector EntityBase::FaceGetPointNum() const {
     Vector r;
     if(type == FACE_NORMAL_PT) {
         r = SK.GetEntity(point[0])->PointGetNum();
@@ -714,12 +714,12 @@ Vector EntityBase::FaceGetPointNum() {
     return r;
 }
 
-bool EntityBase::HasEndpoints() {
+bool EntityBase::HasEndpoints() const {
     return (type == LINE_SEGMENT) ||
            (type == CUBIC) ||
            (type == ARC_OF_CIRCLE);
 }
-Vector EntityBase::EndpointStart() {
+Vector EntityBase::EndpointStart() const {
     if(type == LINE_SEGMENT) {
         return SK.GetEntity(point[0])->PointGetNum();
     } else if(type == CUBIC) {
@@ -728,7 +728,7 @@ Vector EntityBase::EndpointStart() {
         return SK.GetEntity(point[1])->PointGetNum();
     } else ssassert(false, "Unexpected entity type");
 }
-Vector EntityBase::EndpointFinish() {
+Vector EntityBase::EndpointFinish() const {
     if(type == LINE_SEGMENT) {
         return SK.GetEntity(point[1])->PointGetNum();
     } else if(type == CUBIC) {
@@ -738,14 +738,14 @@ Vector EntityBase::EndpointFinish() {
     } else ssassert(false, "Unexpected entity type");
 }
 
-void EntityBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) {
+void EntityBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) const {
     Equation eq;
     eq.e = expr;
     eq.h = h.equation(index);
     l->Add(&eq);
 }
 
-void EntityBase::GenerateEquations(IdList<Equation,hEquation> *l) {
+void EntityBase::GenerateEquations(IdList<Equation,hEquation> *l) const {
     switch(type) {
         case NORMAL_IN_3D: {
             ExprQuaternion q = NormalGetExprs();
diff --git a/src/export.cpp b/src/export.cpp
index f38f0fc..427de60 100644
--- a/src/export.cpp
+++ b/src/export.cpp
@@ -235,9 +235,9 @@ void SolveSpaceUI::ExportWireframeCurves(SEdgeList *sel, SBezierList *sbl,
 }
 
 void SolveSpaceUI::ExportLinesAndMesh(SEdgeList *sel, SBezierList *sbl, SMesh *sm,
-                                    Vector u, Vector v, Vector n,
-                                        Vector origin, double cameraTan,
-                                    VectorFileWriter *out)
+                                      Vector u, Vector v, Vector n,
+                                      Vector origin, double cameraTan,
+                                      VectorFileWriter *out)
 {
     double s = 1.0 / SS.exportScale;
 
diff --git a/src/expr.cpp b/src/expr.cpp
index 2f51000..08691f0 100644
--- a/src/expr.cpp
+++ b/src/expr.cpp
@@ -37,7 +37,7 @@ ExprVector ExprVector::From(double x, double y, double z) {
     return ve;
 }
 
-ExprVector ExprVector::Minus(ExprVector b) {
+ExprVector ExprVector::Minus(ExprVector b) const {
     ExprVector r;
     r.x = x->Minus(b.x);
     r.y = y->Minus(b.y);
@@ -45,7 +45,7 @@ ExprVector ExprVector::Minus(ExprVector b) {
     return r;
 }
 
-ExprVector ExprVector::Plus(ExprVector b) {
+ExprVector ExprVector::Plus(ExprVector b) const {
     ExprVector r;
     r.x = x->Plus(b.x);
     r.y = y->Plus(b.y);
@@ -53,7 +53,7 @@ ExprVector ExprVector::Plus(ExprVector b) {
     return r;
 }
 
-Expr *ExprVector::Dot(ExprVector b) {
+Expr *ExprVector::Dot(ExprVector b) const {
     Expr *r;
     r =         x->Times(b.x);
     r = r->Plus(y->Times(b.y));
@@ -61,7 +61,7 @@ Expr *ExprVector::Dot(ExprVector b) {
     return r;
 }
 
-ExprVector ExprVector::Cross(ExprVector b) {
+ExprVector ExprVector::Cross(ExprVector b) const {
     ExprVector r;
     r.x = (y->Times(b.z))->Minus(z->Times(b.y));
     r.y = (z->Times(b.x))->Minus(x->Times(b.z));
@@ -69,7 +69,7 @@ ExprVector ExprVector::Cross(ExprVector b) {
     return r;
 }
 
-ExprVector ExprVector::ScaledBy(Expr *s) {
+ExprVector ExprVector::ScaledBy(Expr *s) const {
     ExprVector r;
     r.x = x->Times(s);
     r.y = y->Times(s);
@@ -77,12 +77,12 @@ ExprVector ExprVector::ScaledBy(Expr *s) {
     return r;
 }
 
-ExprVector ExprVector::WithMagnitude(Expr *s) {
+ExprVector ExprVector::WithMagnitude(Expr *s) const {
     Expr *m = Magnitude();
     return ScaledBy(s->Div(m));
 }
 
-Expr *ExprVector::Magnitude() {
+Expr *ExprVector::Magnitude() const {
     Expr *r;
     r =         x->Square();
     r = r->Plus(y->Square());
@@ -90,7 +90,7 @@ Expr *ExprVector::Magnitude() {
     return r->Sqrt();
 }
 
-Vector ExprVector::Eval() {
+Vector ExprVector::Eval() const {
     Vector r;
     r.x = x->Eval();
     r.y = y->Eval();
@@ -126,7 +126,7 @@ ExprQuaternion ExprQuaternion::From(Quaternion qn) {
     return qe;
 }
 
-ExprVector ExprQuaternion::RotationU() {
+ExprVector ExprQuaternion::RotationU() const {
     ExprVector u;
     Expr *two = Expr::From(2);
 
@@ -144,7 +144,7 @@ ExprVector ExprQuaternion::RotationU() {
     return u;
 }
 
-ExprVector ExprQuaternion::RotationV() {
+ExprVector ExprQuaternion::RotationV() const {
     ExprVector v;
     Expr *two = Expr::From(2);
 
@@ -162,7 +162,7 @@ ExprVector ExprQuaternion::RotationV() {
     return v;
 }
 
-ExprVector ExprQuaternion::RotationN() {
+ExprVector ExprQuaternion::RotationN() const {
     ExprVector n;
     Expr *two = Expr::From(2);
 
@@ -180,14 +180,14 @@ ExprVector ExprQuaternion::RotationN() {
     return n;
 }
 
-ExprVector ExprQuaternion::Rotate(ExprVector p) {
+ExprVector ExprQuaternion::Rotate(ExprVector p) const {
     // Express the point in the new basis
     return (RotationU().ScaledBy(p.x)).Plus(
             RotationV().ScaledBy(p.y)).Plus(
             RotationN().ScaledBy(p.z));
 }
 
-ExprQuaternion ExprQuaternion::Times(ExprQuaternion b) {
+ExprQuaternion ExprQuaternion::Times(ExprQuaternion b) const {
     Expr *sa = w, *sb = b.w;
     ExprVector va = { vx, vy, vz };
     ExprVector vb = { b.vx, b.vy, b.vz };
@@ -203,7 +203,7 @@ ExprQuaternion ExprQuaternion::Times(ExprQuaternion b) {
     return r;
 }
 
-Expr *ExprQuaternion::Magnitude() {
+Expr *ExprQuaternion::Magnitude() const {
     return ((w ->Square())->Plus(
             (vx->Square())->Plus(
             (vy->Square())->Plus(
@@ -262,7 +262,7 @@ Expr *Expr::AnyOp(int newOp, Expr *b) {
     return r;
 }
 
-int Expr::Children() {
+int Expr::Children() const {
     switch(op) {
         case PARAM:
         case PARAM_PTR:
@@ -288,7 +288,7 @@ int Expr::Children() {
     }
 }
 
-int Expr::Nodes() {
+int Expr::Nodes() const {
     switch(Children()) {
         case 0: return 1;
         case 1: return 1 + a->Nodes();
@@ -297,7 +297,7 @@ int Expr::Nodes() {
     }
 }
 
-Expr *Expr::DeepCopy() {
+Expr *Expr::DeepCopy() const {
     Expr *n = AllocExpr();
     *n = *this;
     int c = n->Children();
@@ -307,7 +307,7 @@ Expr *Expr::DeepCopy() {
 }
 
 Expr *Expr::DeepCopyWithParamsAsPointers(IdList<Param,hParam> *firstTry,
-    IdList<Param,hParam> *thenTry)
+    IdList<Param,hParam> *thenTry) const
 {
     Expr *n = AllocExpr();
     if(op == PARAM) {
@@ -333,7 +333,7 @@ Expr *Expr::DeepCopyWithParamsAsPointers(IdList<Param,hParam> *firstTry,
     return n;
 }
 
-double Expr::Eval() {
+double Expr::Eval() const {
     switch(op) {
         case PARAM:         return SK.GetParam(parh)->val;
         case PARAM_PTR:     return parp->val;
@@ -357,7 +357,7 @@ double Expr::Eval() {
     }
 }
 
-Expr *Expr::PartialWrt(hParam p) {
+Expr *Expr::PartialWrt(hParam p) const {
     Expr *da, *db;
 
     switch(op) {
@@ -400,7 +400,7 @@ Expr *Expr::PartialWrt(hParam p) {
     }
 }
 
-uint64_t Expr::ParamsUsed() {
+uint64_t Expr::ParamsUsed() const {
     uint64_t r = 0;
     if(op == PARAM)     r |= ((uint64_t)1 << (parh.v % 61));
     if(op == PARAM_PTR) r |= ((uint64_t)1 << (parp->h.v % 61));
@@ -411,7 +411,7 @@ uint64_t Expr::ParamsUsed() {
     return r;
 }
 
-bool Expr::DependsOn(hParam p) {
+bool Expr::DependsOn(hParam p) const {
     if(op == PARAM)     return (parh.v    == p.v);
     if(op == PARAM_PTR) return (parp->h.v == p.v);
 
@@ -510,7 +510,7 @@ void Expr::Substitute(hParam oldh, hParam newh) {
 //-----------------------------------------------------------------------------
 const hParam Expr::NO_PARAMS       = { 0 };
 const hParam Expr::MULTIPLE_PARAMS = { 1 };
-hParam Expr::ReferencedParams(ParamList *pl) {
+hParam Expr::ReferencedParams(ParamList *pl) const {
     if(op == PARAM) {
         if(pl->FindByIdNoOops(parh)) {
             return parh;
@@ -546,7 +546,7 @@ hParam Expr::ReferencedParams(ParamList *pl) {
 // Routines to pretty-print an expression. Mostly for debugging.
 //-----------------------------------------------------------------------------
 
-std::string Expr::Print() {
+std::string Expr::Print() const {
 
     char c;
     switch(op) {
diff --git a/src/expr.h b/src/expr.h
index 42222ac..6b4bfc3 100644
--- a/src/expr.h
+++ b/src/expr.h
@@ -80,33 +80,33 @@ public:
     inline Expr *ASin  () { return AnyOp(ASIN,   NULL); }
     inline Expr *ACos  () { return AnyOp(ACOS,   NULL); }
 
-    Expr *PartialWrt(hParam p);
-    double Eval();
-    uint64_t ParamsUsed();
-    bool DependsOn(hParam p);
+    Expr *PartialWrt(hParam p) const;
+    double Eval() const;
+    uint64_t ParamsUsed() const;
+    bool DependsOn(hParam p) const;
     static bool Tol(double a, double b);
     Expr *FoldConstants();
     void Substitute(hParam oldh, hParam newh);
 
     static const hParam NO_PARAMS, MULTIPLE_PARAMS;
-    hParam ReferencedParams(ParamList *pl);
+    hParam ReferencedParams(ParamList *pl) const;
 
     void ParamsToPointers();
 
-    std::string Print();
+    std::string Print() const;
 
     // number of child nodes: 0 (e.g. constant), 1 (sqrt), or 2 (+)
-    int Children();
+    int Children() const;
     // total number of nodes in the tree
-    int Nodes();
+    int Nodes() const;
 
     // Make a simple copy
-    Expr *DeepCopy();
+    Expr *DeepCopy() const;
     // Make a copy, with the parameters (usually referenced by hParam)
     // resolved to pointers to the actual value. This speeds things up
     // considerably.
     Expr *DeepCopyWithParamsAsPointers(IdList<Param,hParam> *firstTry,
-        IdList<Param,hParam> *thenTry);
+        IdList<Param,hParam> *thenTry) const;
 
     static Expr *From(const char *in, bool popUpError);
     static void  Lex(const char *in);
@@ -136,15 +136,15 @@ public:
     static ExprVector From(hParam x, hParam y, hParam z);
     static ExprVector From(double x, double y, double z);
 
-    ExprVector Plus(ExprVector b);
-    ExprVector Minus(ExprVector b);
-    Expr *Dot(ExprVector b);
-    ExprVector Cross(ExprVector b);
-    ExprVector ScaledBy(Expr *s);
-    ExprVector WithMagnitude(Expr *s);
-    Expr *Magnitude();
+    ExprVector Plus(ExprVector b) const;
+    ExprVector Minus(ExprVector b) const;
+    Expr *Dot(ExprVector b) const;
+    ExprVector Cross(ExprVector b) const;
+    ExprVector ScaledBy(Expr *s) const;
+    ExprVector WithMagnitude(Expr *s) const;
+    Expr *Magnitude() const;
 
-    Vector Eval();
+    Vector Eval() const;
 };
 
 class ExprQuaternion {
@@ -155,14 +155,14 @@ public:
     static ExprQuaternion From(Quaternion qn);
     static ExprQuaternion From(hParam w, hParam vx, hParam vy, hParam vz);
 
-    ExprVector RotationU();
-    ExprVector RotationV();
-    ExprVector RotationN();
+    ExprVector RotationU() const;
+    ExprVector RotationV() const;
+    ExprVector RotationN() const;
 
-    ExprVector Rotate(ExprVector p);
-    ExprQuaternion Times(ExprQuaternion b);
+    ExprVector Rotate(ExprVector p) const;
+    ExprQuaternion Times(ExprQuaternion b) const;
 
-    Expr *Magnitude();
+    Expr *Magnitude() const;
 };
 
 #endif
diff --git a/src/mesh.cpp b/src/mesh.cpp
index fe30f94..d04eb1e 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -33,12 +33,12 @@ void SMesh::AddTriangle(STriMeta meta, Vector a, Vector b, Vector c) {
     t.c = c;
     AddTriangle(&t);
 }
-void SMesh::AddTriangle(STriangle *st) {
+void SMesh::AddTriangle(const STriangle *st) {
     if(st->meta.color.alpha != 255) isTransparent = true;
     l.Add(st);
 }
 
-void SMesh::DoBounding(Vector v, Vector *vmax, Vector *vmin) {
+void SMesh::DoBounding(Vector v, Vector *vmax, Vector *vmin) const {
     vmax->x = max(vmax->x, v.x);
     vmax->y = max(vmax->y, v.y);
     vmax->z = max(vmax->z, v.z);
@@ -47,7 +47,7 @@ void SMesh::DoBounding(Vector v, Vector *vmax, Vector *vmin) {
     vmin->y = min(vmin->y, v.y);
     vmin->z = min(vmin->z, v.z);
 }
-void SMesh::GetBounding(Vector *vmax, Vector *vmin) {
+void SMesh::GetBounding(Vector *vmax, Vector *vmin) const {
     int i;
     *vmin = Vector::From( 1e12,  1e12,  1e12);
     *vmax = Vector::From(-1e12, -1e12, -1e12);
@@ -302,8 +302,8 @@ void SMesh::MakeFromAssemblyOf(SMesh *a, SMesh *b) {
     MakeFromCopyOf(b);
 }
 
-void SMesh::MakeFromTransformationOf(SMesh *a,
-                                      Vector trans, Quaternion q, double scale)
+void SMesh::MakeFromTransformationOf(SMesh *a, Vector trans,
+                                     Quaternion q, double scale)
 {
     STriangle *tr;
     for(tr = a->l.First(); tr; tr = a->l.NextAfter(tr)) {
@@ -322,11 +322,11 @@ void SMesh::MakeFromTransformationOf(SMesh *a,
     }
 }
 
-bool SMesh::IsEmpty() {
+bool SMesh::IsEmpty() const {
     return (l.n == 0);
 }
 
-uint32_t SMesh::FirstIntersectionWith(Point2d mp) {
+uint32_t SMesh::FirstIntersectionWith(Point2d mp) const {
     Vector p0 = Vector::From(mp.x, mp.y, 0);
     Vector gn = Vector::From(0, 0, 1);
 
@@ -487,7 +487,7 @@ leaf:
     return ret;
 }
 
-void SKdNode::ClearTags() {
+void SKdNode::ClearTags() const {
     if(gt && lt) {
         gt->ClearTags();
         lt->ClearTags();
@@ -524,7 +524,7 @@ void SKdNode::AddTriangle(STriangle *tr) {
     }
 }
 
-void SKdNode::MakeMeshInto(SMesh *m) {
+void SKdNode::MakeMeshInto(SMesh *m) const {
     if(gt) gt->MakeMeshInto(m);
     if(lt) lt->MakeMeshInto(m);
 
@@ -537,7 +537,7 @@ void SKdNode::MakeMeshInto(SMesh *m) {
     }
 }
 
-void SKdNode::ListTrianglesInto(std::vector<STriangle *> *tl) {
+void SKdNode::ListTrianglesInto(std::vector<STriangle *> *tl) const {
     if(gt) gt->ListTrianglesInto(tl);
     if(lt) lt->ListTrianglesInto(tl);
 
@@ -654,7 +654,7 @@ void SKdNode::SnapToMesh(SMesh *m) {
 // them for occlusion. Keep only the visible segments. sel is both our input
 // and our output.
 //-----------------------------------------------------------------------------
-void SKdNode::SplitLinesAgainstTriangle(SEdgeList *sel, STriangle *tr, bool removeHidden) {
+void SKdNode::SplitLinesAgainstTriangle(SEdgeList *sel, STriangle *tr, bool removeHidden) const {
     SEdgeList seln = {};
 
     Vector tn = tr->Normal().WithMagnitude(1);
@@ -767,7 +767,7 @@ void SKdNode::SplitLinesAgainstTriangle(SEdgeList *sel, STriangle *tr, bool remo
 // Given an edge orig, occlusion test it against our mesh. We output an edge
 // list in sel, containing the visible portions of that edge.
 //-----------------------------------------------------------------------------
-void SKdNode::OcclusionTestLine(SEdge orig, SEdgeList *sel, int cnt, bool removeHidden) {
+void SKdNode::OcclusionTestLine(SEdge orig, SEdgeList *sel, int cnt, bool removeHidden) const {
     if(gt && lt) {
         double ac = (orig.a).Element(which),
                bc = (orig.b).Element(which);
@@ -808,7 +808,7 @@ void SKdNode::OcclusionTestLine(SEdge orig, SEdgeList *sel, int cnt, bool remove
 // with a triangle in the mesh, otherwise not.
 //-----------------------------------------------------------------------------
 void SKdNode::FindEdgeOn(Vector a, Vector b, int cnt, bool coplanarIsInter,
-                         EdgeOnInfo *info)
+                         EdgeOnInfo *info) const
 {
     if(gt && lt) {
         double ac = a.Element(which),
@@ -925,7 +925,7 @@ static bool CheckAndAddTrianglePair(std::set<std::pair<STriangle *, STriangle *>
 //      a triangle from a different face)
 //-----------------------------------------------------------------------------
 void SKdNode::MakeCertainEdgesInto(SEdgeList *sel, int how, bool coplanarIsInter,
-                                   bool *inter, bool *leaky, int auxA)
+                                   bool *inter, bool *leaky, int auxA) const
 {
     if(inter) *inter = false;
     if(leaky) *leaky = false;
@@ -1022,7 +1022,7 @@ void SKdNode::MakeCertainEdgesInto(SEdgeList *sel, int how, bool coplanarIsInter
     }
 }
 
-void SKdNode::MakeOutlinesInto(SOutlineList *sol)
+void SKdNode::MakeOutlinesInto(SOutlineList *sol) const
 {
     std::vector<STriangle *> tris;
     ClearTags();
diff --git a/src/polygon.cpp b/src/polygon.cpp
index f035f10..59bca0b 100644
--- a/src/polygon.cpp
+++ b/src/polygon.cpp
@@ -5,12 +5,12 @@
 //-----------------------------------------------------------------------------
 #include "solvespace.h"
 
-Vector STriangle::Normal() {
+Vector STriangle::Normal() const {
     Vector ab = b.Minus(a), bc = c.Minus(b);
     return ab.Cross(bc);
 }
 
-double STriangle::MinAltitude() {
+double STriangle::MinAltitude() const {
     double altA = a.DistanceToLine(b, c.Minus(b)),
            altB = b.DistanceToLine(c, a.Minus(c)),
            altC = c.DistanceToLine(a, b.Minus(a));
@@ -18,7 +18,7 @@ double STriangle::MinAltitude() {
     return min(altA, min(altB, altC));
 }
 
-bool STriangle::ContainsPoint(Vector p) {
+bool STriangle::ContainsPoint(Vector p) const {
     Vector n = Normal();
     if(MinAltitude() < LENGTH_EPS) {
         // shouldn't happen; zero-area triangle
@@ -27,7 +27,7 @@ bool STriangle::ContainsPoint(Vector p) {
     return ContainsPointProjd(n.WithMagnitude(1), p);
 }
 
-bool STriangle::ContainsPointProjd(Vector n, Vector p) {
+bool STriangle::ContainsPointProjd(Vector n, Vector p) const {
     Vector ab = b.Minus(a), bc = c.Minus(b), ca = a.Minus(c);
 
     Vector no_ab = n.Cross(ab);
@@ -63,7 +63,7 @@ SEdge SEdge::From(Vector a, Vector b) {
     return se;
 }
 
-bool SEdge::EdgeCrosses(Vector ea, Vector eb, Vector *ppi, SPointList *spl) {
+bool SEdge::EdgeCrosses(Vector ea, Vector eb, Vector *ppi, SPointList *spl) const {
     Vector d = eb.Minus(ea);
     double t_eps = LENGTH_EPS/d.Magnitude();
 
@@ -155,7 +155,7 @@ void SEdgeList::AddEdge(Vector a, Vector b, int auxA, int auxB) {
 }
 
 bool SEdgeList::AssembleContour(Vector first, Vector last, SContour *dest,
-                                SEdge *errorAt, bool keepDir)
+                                SEdge *errorAt, bool keepDir) const
 {
     int i;
 
@@ -195,7 +195,7 @@ bool SEdgeList::AssembleContour(Vector first, Vector last, SContour *dest,
     return true;
 }
 
-bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir) {
+bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir) const {
     dest->Clear();
 
     bool allClosed = true;
@@ -234,11 +234,10 @@ bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir) {
 // If pi is not NULL, then a crossing is returned in that.
 //-----------------------------------------------------------------------------
 int SEdgeList::AnyEdgeCrossings(Vector a, Vector b,
-                                Vector *ppi, SPointList *spl)
+                                Vector *ppi, SPointList *spl) const
 {
     int cnt = 0;
-    SEdge *se;
-    for(se = l.First(); se; se = l.NextAfter(se)) {
+    for(const SEdge *se = l.First(); se; se = l.NextAfter(se)) {
         if(se->EdgeCrosses(a, b, ppi, spl)) {
             cnt++;
         }
@@ -250,16 +249,14 @@ int SEdgeList::AnyEdgeCrossings(Vector a, Vector b,
 // Returns true if the intersecting edge list contains an edge that shares
 // an endpoint with one of our edges.
 //-----------------------------------------------------------------------------
-bool SEdgeList::ContainsEdgeFrom(SEdgeList *sel) {
-    SEdge *se;
-    for(se = l.First(); se; se = l.NextAfter(se)) {
+bool SEdgeList::ContainsEdgeFrom(const SEdgeList *sel) const {
+    for(const SEdge *se = l.First(); se; se = l.NextAfter(se)) {
         if(sel->ContainsEdge(se)) return true;
     }
     return false;
 }
-bool SEdgeList::ContainsEdge(SEdge *set) {
-    SEdge *se;
-    for(se = l.First(); se; se = l.NextAfter(se)) {
+bool SEdgeList::ContainsEdge(const SEdge *set) const {
+    for(const SEdge *se = l.First(); se; se = l.NextAfter(se)) {
         if((se->a).Equals(set->a) && (se->b).Equals(set->b)) return true;
         if((se->b).Equals(set->a) && (se->a).Equals(set->b)) return true;
     }
@@ -395,7 +392,7 @@ SKdNodeEdges *SKdNodeEdges::From(SEdgeLl *sell) {
 }
 
 int SKdNodeEdges::AnyEdgeCrossings(Vector a, Vector b, int cnt,
-        Vector *pi, SPointList *spl)
+        Vector *pi, SPointList *spl) const
 {
     int inters = 0;
     if(gt && lt) {
@@ -462,11 +459,11 @@ void SPointList::Clear() {
     l.Clear();
 }
 
-bool SPointList::ContainsPoint(Vector pt) {
+bool SPointList::ContainsPoint(Vector pt) const {
     return (IndexForPoint(pt) >= 0);
 }
 
-int SPointList::IndexForPoint(Vector pt) {
+int SPointList::IndexForPoint(Vector pt) const {
     int i;
     for(i = 0; i < l.n; i++) {
         SPoint *p = &(l.elem[i]);
@@ -506,31 +503,29 @@ void SContour::AddPoint(Vector p) {
     l.Add(&sp);
 }
 
-void SContour::MakeEdgesInto(SEdgeList *el) {
+void SContour::MakeEdgesInto(SEdgeList *el) const {
     int i;
     for(i = 0; i < (l.n - 1); i++) {
         el->AddEdge(l.elem[i].p, l.elem[i+1].p);
     }
 }
 
-void SContour::CopyInto(SContour *dest) {
-    SPoint *sp;
-    for(sp = l.First(); sp; sp = l.NextAfter(sp)) {
+void SContour::CopyInto(SContour *dest) const {
+    for(const SPoint *sp = l.First(); sp; sp = l.NextAfter(sp)) {
         dest->AddPoint(sp->p);
     }
 }
 
 void SContour::FindPointWithMinX() {
-    SPoint *sp;
     xminPt = Vector::From(1e10, 1e10, 1e10);
-    for(sp = l.First(); sp; sp = l.NextAfter(sp)) {
+    for(const SPoint *sp = l.First(); sp; sp = l.NextAfter(sp)) {
         if(sp->p.x < xminPt.x) {
             xminPt = sp->p;
         }
     }
 }
 
-Vector SContour::ComputeNormal() {
+Vector SContour::ComputeNormal() const {
     Vector n = Vector::From(0, 0, 0);
 
     for(int i = 0; i < l.n - 2; i++) {
@@ -544,12 +539,12 @@ Vector SContour::ComputeNormal() {
     return n.WithMagnitude(1);
 }
 
-Vector SContour::AnyEdgeMidpoint() {
+Vector SContour::AnyEdgeMidpoint() const {
     ssassert(l.n >= 2, "Need two points to find a midpoint");
     return ((l.elem[0].p).Plus(l.elem[1].p)).ScaledBy(0.5);
 }
 
-bool SContour::IsClockwiseProjdToNormal(Vector n) {
+bool SContour::IsClockwiseProjdToNormal(Vector n) const {
     // Degenerate things might happen as we draw; doesn't really matter
     // what we do then.
     if(n.Magnitude() < 0.01) return true;
@@ -557,7 +552,7 @@ bool SContour::IsClockwiseProjdToNormal(Vector n) {
     return (SignedAreaProjdToNormal(n) < 0);
 }
 
-double SContour::SignedAreaProjdToNormal(Vector n) {
+double SContour::SignedAreaProjdToNormal(Vector n) const {
     // An arbitrary 2d coordinate system that has n as its normal
     Vector u = n.Normal(0);
     Vector v = n.Normal(1);
@@ -574,7 +569,7 @@ double SContour::SignedAreaProjdToNormal(Vector n) {
     return area;
 }
 
-bool SContour::ContainsPointProjdToNormal(Vector n, Vector p) {
+bool SContour::ContainsPointProjdToNormal(Vector n, Vector p) const {
     Vector u = n.Normal(0);
     Vector v = n.Normal(1);
 
@@ -618,7 +613,7 @@ void SPolygon::AddEmptyContour() {
     l.Add(&c);
 }
 
-void SPolygon::MakeEdgesInto(SEdgeList *el) {
+void SPolygon::MakeEdgesInto(SEdgeList *el) const {
     int i;
     for(i = 0; i < l.n; i++) {
         (l.elem[i]).MakeEdgesInto(el);
@@ -630,22 +625,21 @@ Vector SPolygon::ComputeNormal() {
     return (l.elem[0]).ComputeNormal();
 }
 
-double SPolygon::SignedArea() {
-    SContour *sc;
+double SPolygon::SignedArea() const {
     double area = 0;
     // This returns the true area only if the contours are all oriented
     // correctly, with the holes backwards from the outer contours.
-    for(sc = l.First(); sc; sc = l.NextAfter(sc)) {
+    for(const SContour *sc = l.First(); sc; sc = l.NextAfter(sc)) {
         area += sc->SignedAreaProjdToNormal(normal);
     }
     return area;
 }
 
-bool SPolygon::ContainsPoint(Vector p) {
+bool SPolygon::ContainsPoint(Vector p) const {
     return (WindingNumberForPoint(p) % 2) == 1;
 }
 
-int SPolygon::WindingNumberForPoint(Vector p) {
+int SPolygon::WindingNumberForPoint(Vector p) const {
     int winding = 0;
     int i;
     for(i = 0; i < l.n; i++) {
@@ -690,17 +684,17 @@ void SPolygon::FixContourDirections() {
     }
 }
 
-bool SPolygon::IsEmpty() {
+bool SPolygon::IsEmpty() const {
     if(l.n == 0 || l.elem[0].l.n == 0) return true;
     return false;
 }
 
-Vector SPolygon::AnyPoint() {
+Vector SPolygon::AnyPoint() const {
     ssassert(!IsEmpty(), "Need at least one point");
     return l.elem[0].l.elem[0].p;
 }
 
-bool SPolygon::SelfIntersecting(Vector *intersectsAt) {
+bool SPolygon::SelfIntersecting(Vector *intersectsAt) const {
     SEdgeList el = {};
     MakeEdgesInto(&el);
     SKdNodeEdges *kdtree = SKdNodeEdges::From(&el);
@@ -728,7 +722,7 @@ bool SPolygon::SelfIntersecting(Vector *intersectsAt) {
 // polygon is in the xy plane, and the contours all go in the right direction
 // with respect to normal (0, 0, -1).
 //-----------------------------------------------------------------------------
-void SPolygon::OffsetInto(SPolygon *dest, double r) {
+void SPolygon::OffsetInto(SPolygon *dest, double r) const {
     int i;
     dest->Clear();
     for(i = 0; i < l.n; i++) {
@@ -781,7 +775,7 @@ static bool IntersectionOfLines(double x0A, double y0A, double dxA, double dyA,
 
     return true;
 }
-void SContour::OffsetInto(SContour *dest, double r) {
+void SContour::OffsetInto(SContour *dest, double r) const {
     int i;
 
     for(i = 0; i < l.n; i++) {
diff --git a/src/polygon.h b/src/polygon.h
index 5b4084f..737594a 100644
--- a/src/polygon.h
+++ b/src/polygon.h
@@ -22,7 +22,7 @@ public:
     Vector a, b;
 
     static SEdge From(Vector a, Vector b);
-    bool EdgeCrosses(Vector a, Vector b, Vector *pi=NULL, SPointList *spl=NULL);
+    bool EdgeCrosses(Vector a, Vector b, Vector *pi=NULL, SPointList *spl=NULL) const;
 };
 
 class SEdgeList {
@@ -31,13 +31,13 @@ public:
 
     void Clear();
     void AddEdge(Vector a, Vector b, int auxA=0, int auxB=0);
-    bool AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir=false);
+    bool AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir=false) const;
     bool AssembleContour(Vector first, Vector last, SContour *dest,
-                            SEdge *errorAt, bool keepDir);
+                            SEdge *errorAt, bool keepDir) const;
     int AnyEdgeCrossings(Vector a, Vector b,
-        Vector *pi=NULL, SPointList *spl=NULL);
-    bool ContainsEdgeFrom(SEdgeList *sel);
-    bool ContainsEdge(SEdge *se);
+        Vector *pi=NULL, SPointList *spl=NULL) const;
+    bool ContainsEdgeFrom(const SEdgeList *sel) const;
+    bool ContainsEdge(const SEdge *se) const;
     void CullExtraneousEdges();
     void MergeCollinearSegments(Vector a, Vector b);
 };
@@ -68,7 +68,7 @@ public:
     static SKdNodeEdges *From(SEdgeLl *sell);
     static SKdNodeEdges *Alloc();
     int AnyEdgeCrossings(Vector a, Vector b, int cnt,
-        Vector *pi=NULL, SPointList *spl=NULL);
+        Vector *pi=NULL, SPointList *spl=NULL) const;
 };
 
 class SPoint {
@@ -91,8 +91,8 @@ public:
     List<SPoint>    l;
 
     void Clear();
-    bool ContainsPoint(Vector pt);
-    int IndexForPoint(Vector pt);
+    bool ContainsPoint(Vector pt) const;
+    int IndexForPoint(Vector pt) const;
     void IncrementTagFor(Vector pt);
     void Add(Vector pt);
 };
@@ -105,18 +105,18 @@ public:
     List<SPoint>    l;
 
     void AddPoint(Vector p);
-    void MakeEdgesInto(SEdgeList *el);
+    void MakeEdgesInto(SEdgeList *el) const;
     void Reverse();
-    Vector ComputeNormal();
-    double SignedAreaProjdToNormal(Vector n);
-    bool IsClockwiseProjdToNormal(Vector n);
-    bool ContainsPointProjdToNormal(Vector n, Vector p);
-    void OffsetInto(SContour *dest, double r);
-    void CopyInto(SContour *dest);
+    Vector ComputeNormal() const;
+    double SignedAreaProjdToNormal(Vector n) const;
+    bool IsClockwiseProjdToNormal(Vector n) const;
+    bool ContainsPointProjdToNormal(Vector n, Vector p) const;
+    void OffsetInto(SContour *dest, double r) const;
+    void CopyInto(SContour *dest) const;
     void FindPointWithMinX();
-    Vector AnyEdgeMidpoint();
+    Vector AnyEdgeMidpoint() const;
 
-    bool IsEar(int bp, double scaledEps);
+    bool IsEar(int bp, double scaledEps) const;
     bool BridgeToContour(SContour *sc, SEdgeList *el, List<Vector> *vl);
     void ClipEarInto(SMesh *m, int bp, double scaledEps);
     void UvTriangulateInto(SMesh *m, SSurface *srf);
@@ -134,16 +134,16 @@ public:
 
     Vector ComputeNormal();
     void AddEmptyContour();
-    int WindingNumberForPoint(Vector p);
-    double SignedArea();
-    bool ContainsPoint(Vector p);
-    void MakeEdgesInto(SEdgeList *el);
+    int WindingNumberForPoint(Vector p) const;
+    double SignedArea() const;
+    bool ContainsPoint(Vector p) const;
+    void MakeEdgesInto(SEdgeList *el) const;
     void FixContourDirections();
     void Clear();
-    bool SelfIntersecting(Vector *intersectsAt);
-    bool IsEmpty();
-    Vector AnyPoint();
-    void OffsetInto(SPolygon *dest, double r);
+    bool SelfIntersecting(Vector *intersectsAt) const;
+    bool IsEmpty() const;
+    Vector AnyPoint() const;
+    void OffsetInto(SPolygon *dest, double r) const;
     void UvTriangulateInto(SMesh *m, SSurface *srf);
     void UvGridTriangulateInto(SMesh *m, SSurface *srf);
 };
@@ -156,12 +156,12 @@ public:
     Vector      an, bn, cn;
 
     static STriangle From(STriMeta meta, Vector a, Vector b, Vector c);
-    Vector Normal();
+    Vector Normal() const;
     void FlipNormal();
-    double MinAltitude();
-    int WindingNumberForPoint(Vector p);
-    bool ContainsPoint(Vector p);
-    bool ContainsPointProjd(Vector n, Vector p);
+    double MinAltitude() const;
+    int WindingNumberForPoint(Vector p) const;
+    bool ContainsPoint(Vector p) const;
+    bool ContainsPointProjd(Vector n, Vector p) const;
 };
 
 class SBsp2 {
@@ -180,12 +180,13 @@ public:
     enum { POS = 100, NEG = 101, COPLANAR = 200 };
     void InsertTriangleHow(int how, STriangle *tr, SMesh *m, SBsp3 *bsp3);
     void InsertTriangle(STriangle *tr, SMesh *m, SBsp3 *bsp3);
-    Vector IntersectionWith(Vector a, Vector b);
+    Vector IntersectionWith(Vector a, Vector b) const;
     void InsertEdge(SEdge *nedge, Vector nnp, Vector out);
-    static SBsp2 *InsertOrCreateEdge(SBsp2 *where, SEdge *nedge, Vector nnp, Vector out);
+    static SBsp2 *InsertOrCreateEdge(SBsp2 *where, SEdge *nedge,
+                                     Vector nnp, Vector out);
     static SBsp2 *Alloc();
 
-    void DebugDraw(Vector n, double d);
+    void DebugDraw(Vector n, double d) const;
 };
 
 class SBsp3 {
@@ -204,7 +205,7 @@ public:
     static SBsp3 *Alloc();
     static SBsp3 *FromMesh(SMesh *m);
 
-    Vector IntersectionWith(Vector a, Vector b);
+    Vector IntersectionWith(Vector a, Vector b) const;
 
     enum { POS = 100, NEG = 101, COPLANAR = 200 };
     void InsertHow(int how, STriangle *str, SMesh *instead);
@@ -217,9 +218,9 @@ public:
 
     void InsertInPlane(bool pos2, STriangle *tr, SMesh *m);
 
-    void GenerateInPaintOrder(SMesh *m);
+    void GenerateInPaintOrder(SMesh *m) const;
 
-    void DebugDraw();
+    void DebugDraw() const;
 };
 
 class SMesh {
@@ -232,11 +233,12 @@ public:
     bool    isTransparent;
 
     void Clear();
-    void AddTriangle(STriangle *st);
+    void AddTriangle(const STriangle *st);
     void AddTriangle(STriMeta meta, Vector a, Vector b, Vector c);
-    void AddTriangle(STriMeta meta, Vector n, Vector a, Vector b, Vector c);
-    void DoBounding(Vector v, Vector *vmax, Vector *vmin);
-    void GetBounding(Vector *vmax, Vector *vmin);
+    void AddTriangle(STriMeta meta, Vector n,
+                     Vector a, Vector b, Vector c);
+    void DoBounding(Vector v, Vector *vmax, Vector *vmin) const;
+    void GetBounding(Vector *vmax, Vector *vmin) const;
 
     void Simplify(int start);
 
@@ -245,17 +247,17 @@ public:
     void MakeFromDifferenceOf(SMesh *a, SMesh *b);
 
     void MakeFromCopyOf(SMesh *a);
-    void MakeFromTransformationOf(SMesh *a,
-                                    Vector trans, Quaternion q, double scale);
+    void MakeFromTransformationOf(SMesh *a, Vector trans,
+                                  Quaternion q, double scale);
     void MakeFromAssemblyOf(SMesh *a, SMesh *b);
 
     void MakeEdgesInPlaneInto(SEdgeList *sel, Vector n, double d);
     void MakeCertainEdgesAndOutlinesInto(SEdgeList *sel, SOutlineList *sol, int type);
 
-    bool IsEmpty();
+    bool IsEmpty() const;
     void RemapFaces(Group *g, int remap);
 
-    uint32_t FirstIntersectionWith(Point2d mp);
+    uint32_t FirstIntersectionWith(Point2d mp) const;
 };
 
 // A linked list of triangles
@@ -310,11 +312,11 @@ public:
     static SKdNode *From(STriangleLl *tll);
 
     void AddTriangle(STriangle *tr);
-    void MakeMeshInto(SMesh *m);
-    void ListTrianglesInto(std::vector<STriangle *> *tl);
-    void ClearTags();
+    void MakeMeshInto(SMesh *m) const;
+    void ListTrianglesInto(std::vector<STriangle *> *tl) const;
+    void ClearTags() const;
 
-    void FindEdgeOn(Vector a, Vector b, int cnt, bool coplanarIsInter, EdgeOnInfo *info);
+    void FindEdgeOn(Vector a, Vector b, int cnt, bool coplanarIsInter, EdgeOnInfo *info) const;
     enum {
         NAKED_OR_SELF_INTER_EDGES  = 100,
         SELF_INTER_EDGES           = 200,
@@ -323,11 +325,11 @@ public:
         SHARP_EDGES                = 500,
     };
     void MakeCertainEdgesInto(SEdgeList *sel, int how, bool coplanarIsInter,
-                              bool *inter, bool *leaky, int auxA=0);
-    void MakeOutlinesInto(SOutlineList *sel);
+                              bool *inter, bool *leaky, int auxA = 0) const;
+    void MakeOutlinesInto(SOutlineList *sel) const;
 
-    void OcclusionTestLine(SEdge orig, SEdgeList *sel, int cnt, bool removeHidden);
-    void SplitLinesAgainstTriangle(SEdgeList *sel, STriangle *tr, bool removeHidden);
+    void OcclusionTestLine(SEdge orig, SEdgeList *sel, int cnt, bool removeHidden) const;
+    void SplitLinesAgainstTriangle(SEdgeList *sel, STriangle *tr, bool removeHidden) const;
 
     void SnapToMesh(SMesh *m);
     void SnapToVertex(Vector v, SMesh *extras);
diff --git a/src/request.cpp b/src/request.cpp
index aa3a4c8..7c43930 100644
--- a/src/request.cpp
+++ b/src/request.cpp
@@ -83,7 +83,7 @@ int EntReqTable::GetRequestForEntity(int ent) {
 
 
 void Request::Generate(IdList<Entity,hEntity> *entity,
-                       IdList<Param,hParam> *param)
+                       IdList<Param,hParam> *param) const
 {
     int points = 0;
     int et = 0;
@@ -169,7 +169,7 @@ void Request::Generate(IdList<Entity,hEntity> *entity,
     if(et) entity->Add(&e);
 }
 
-std::string Request::DescriptionString() {
+std::string Request::DescriptionString() const {
     const char *s;
     if(h.v == Request::HREQUEST_REFERENCE_XY.v) {
         s = "#XY";
@@ -184,7 +184,7 @@ std::string Request::DescriptionString() {
     return ssprintf("r%03x-%s", h.v, s);
 }
 
-int Request::IndexOfPoint(hEntity he) {
+int Request::IndexOfPoint(hEntity he) const {
     if(type == DATUM_POINT) {
         return (he.v == h.entity(0).v) ? 0 : -1;
     }
diff --git a/src/sketch.h b/src/sketch.h
index 8fc76fc..d545b20 100644
--- a/src/sketch.h
+++ b/src/sketch.h
@@ -21,7 +21,6 @@ class Param;
 class Equation;
 class Style;
 
-
 // All of the hWhatever handles are a 32-bit ID, that is used to represent
 // some data structure in the sketch.
 class hGroup {
@@ -29,19 +28,19 @@ public:
     // bits 15: 0   -- group index
     uint32_t v;
 
-    inline hEntity entity(int i);
-    inline hParam param(int i);
-    inline hEquation equation(int i);
+    inline hEntity entity(int i) const;
+    inline hParam param(int i) const;
+    inline hEquation equation(int i) const;
 };
 class hRequest {
 public:
     // bits 15: 0   -- request index
     uint32_t v;
 
-    inline hEntity entity(int i);
-    inline hParam param(int i);
+    inline hEntity entity(int i) const;
+    inline hParam param(int i) const;
 
-    inline bool IsFromReferences();
+    inline bool IsFromReferences() const;
 };
 class hEntity {
 public:
@@ -49,10 +48,10 @@ public:
     //      31:16   -- request index
     uint32_t v;
 
-    inline bool isFromRequest();
-    inline hRequest request();
-    inline hGroup group();
-    inline hEquation equation(int i);
+    inline bool isFromRequest() const;
+    inline hRequest request() const;
+    inline hGroup group() const;
+    inline hEquation equation(int i) const;
 };
 class hParam {
 public:
@@ -60,7 +59,7 @@ public:
     //      31:16   -- request index
     uint32_t v;
 
-    inline hRequest request();
+    inline hRequest request() const;
 };
 
 class hStyle {
@@ -68,7 +67,6 @@ public:
     uint32_t v;
 };
 
-
 class EntityId {
 public:
     uint32_t v;     // entity ID, starting from 0
@@ -301,10 +299,10 @@ public:
     std::string font;
 
     static hParam AddParam(ParamList *param, hParam hp);
-    void Generate(EntityList *entity, ParamList *param);
+    void Generate(EntityList *entity, ParamList *param) const;
 
-    std::string DescriptionString();
-    int IndexOfPoint(hEntity he);
+    std::string DescriptionString() const;
+    int IndexOfPoint(hEntity he) const;
 
     void Clear() {}
 };
@@ -378,75 +376,75 @@ public:
     // times to apply the transformation.
     int timesApplied;
 
-    Quaternion GetAxisAngleQuaternion(int param0);
-    ExprQuaternion GetAxisAngleQuaternionExprs(int param0);
+    Quaternion GetAxisAngleQuaternion(int param0) const;
+    ExprQuaternion GetAxisAngleQuaternionExprs(int param0) const;
 
-    bool IsCircle();
-    Expr *CircleGetRadiusExpr();
-    double CircleGetRadiusNum();
-    void ArcGetAngles(double *thetaa, double *thetab, double *dtheta);
+    bool IsCircle() const;
+    Expr *CircleGetRadiusExpr() const;
+    double CircleGetRadiusNum() const;
+    void ArcGetAngles(double *thetaa, double *thetab, double *dtheta) const;
 
-    bool HasVector();
-    ExprVector VectorGetExprs();
-    Vector VectorGetNum();
-    Vector VectorGetRefPoint();
-    Vector VectorGetStartPoint();
+    bool HasVector() const;
+    ExprVector VectorGetExprs() const;
+    Vector VectorGetNum() const;
+    Vector VectorGetRefPoint() const;
+    Vector VectorGetStartPoint() const;
 
     // For distances
-    bool IsDistance();
-    double DistanceGetNum();
-    Expr *DistanceGetExpr();
+    bool IsDistance() const;
+    double DistanceGetNum() const;
+    Expr *DistanceGetExpr() const;
     void DistanceForceTo(double v);
 
-    bool IsWorkplane();
+    bool IsWorkplane() const;
     // The plane is points P such that P dot (xn, yn, zn) - d = 0
-    void WorkplaneGetPlaneExprs(ExprVector *n, Expr **d);
-    ExprVector WorkplaneGetOffsetExprs();
-    Vector WorkplaneGetOffset();
-    EntityBase *Normal();
+    void WorkplaneGetPlaneExprs(ExprVector *n, Expr **d) const;
+    ExprVector WorkplaneGetOffsetExprs() const;
+    Vector WorkplaneGetOffset() const;
+    EntityBase *Normal() const;
 
-    bool IsFace();
-    ExprVector FaceGetNormalExprs();
-    Vector FaceGetNormalNum();
-    ExprVector FaceGetPointExprs();
-    Vector FaceGetPointNum();
+    bool IsFace() const;
+    ExprVector FaceGetNormalExprs() const;
+    Vector FaceGetNormalNum() const;
+    ExprVector FaceGetPointExprs() const;
+    Vector FaceGetPointNum() const;
 
-    bool IsPoint();
+    bool IsPoint() const;
     // Applies for any of the point types
-    Vector PointGetNum();
-    ExprVector PointGetExprs();
-    void PointGetExprsInWorkplane(hEntity wrkpl, Expr **u, Expr **v);
+    Vector PointGetNum() const;
+    ExprVector PointGetExprs() const;
+    void PointGetExprsInWorkplane(hEntity wrkpl, Expr **u, Expr **v) const;
     void PointForceTo(Vector v);
     // These apply only the POINT_N_ROT_TRANS, which has an assoc rotation
-    Quaternion PointGetQuaternion();
+    Quaternion PointGetQuaternion() const;
     void PointForceQuaternionTo(Quaternion q);
 
-    bool IsNormal();
+    bool IsNormal() const;
     // Applies for any of the normal types
-    Quaternion NormalGetNum();
-    ExprQuaternion NormalGetExprs();
+    Quaternion NormalGetNum() const;
+    ExprQuaternion NormalGetExprs() const;
     void NormalForceTo(Quaternion q);
 
-    Vector NormalU();
-    Vector NormalV();
-    Vector NormalN();
-    ExprVector NormalExprsU();
-    ExprVector NormalExprsV();
-    ExprVector NormalExprsN();
+    Vector NormalU() const;
+    Vector NormalV() const;
+    Vector NormalN() const;
+    ExprVector NormalExprsU() const;
+    ExprVector NormalExprsV() const;
+    ExprVector NormalExprsN() const;
 
-    Vector CubicGetStartNum();
-    Vector CubicGetFinishNum();
-    ExprVector CubicGetStartTangentExprs();
-    ExprVector CubicGetFinishTangentExprs();
-    Vector CubicGetStartTangentNum();
-    Vector CubicGetFinishTangentNum();
+    Vector CubicGetStartNum() const;
+    Vector CubicGetFinishNum() const;
+    ExprVector CubicGetStartTangentExprs() const;
+    ExprVector CubicGetFinishTangentExprs() const;
+    Vector CubicGetStartTangentNum() const;
+    Vector CubicGetFinishTangentNum() const;
 
-    bool HasEndpoints();
-    Vector EndpointStart();
-    Vector EndpointFinish();
+    bool HasEndpoints() const;
+    Vector EndpointStart() const;
+    Vector EndpointFinish() const;
 
-    void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index);
-    void GenerateEquations(IdList<Equation,hEquation> *l);
+    void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) const;
+    void GenerateEquations(IdList<Equation,hEquation> *l) const;
 
     void Clear() {}
 };
@@ -499,15 +497,16 @@ public:
         int         stippleType;
         int         data;
     } dogd; // state for drawing or getting distance (for hit testing)
+
     void LineDrawOrGetDistance(Vector a, Vector b, bool maybeFat=false, int userData = -1);
     void DrawOrGetDistance();
 
-    bool IsStylable();
-    bool IsVisible();
-    bool PointIsFromReferences();
+    bool IsStylable() const;
+    bool IsVisible() const;
+    bool PointIsFromReferences() const;
 
-    void ComputeInterpolatingSpline(SBezierList *sbl, bool periodic);
-    void GenerateBezierCurves(SBezierList *sbl);
+    void ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) const;
+    void GenerateBezierCurves(SBezierList *sbl) const;
     void GenerateEdges(SEdgeList *el, bool includingConstruction=false);
 
     static void DrawAll(bool drawAsHidden);
@@ -517,7 +516,7 @@ public:
 
     void CalculateNumerical(bool forExport);
 
-    std::string DescriptionString();
+    std::string DescriptionString() const;
 
     SBezierList *GetOrGenerateBezierCurves();
     SEdgeList *GetOrGenerateEdges();
@@ -575,7 +574,7 @@ class hConstraint {
 public:
     uint32_t v;
 
-    inline hEquation equation(int i);
+    inline hEquation equation(int i) const;
 };
 
 class ConstraintBase {
@@ -643,13 +642,13 @@ public:
     bool        reference;  // a ref dimension, that generates no eqs
     std::string comment;    // since comments are represented as constraints
 
-    bool HasLabel();
+    bool HasLabel() const;
 
-    void Generate(IdList<Equation,hEquation> *l);
-    void GenerateReal(IdList<Equation,hEquation> *l);
+    void Generate(IdList<Equation,hEquation> *l) const;
+    void GenerateReal(IdList<Equation,hEquation> *l) const;
     // Some helpers when generating symbolic constraint equations
     void ModifyToSatisfy();
-    void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index);
+    void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) const;
     static Expr *DirectionCosine(hEntity wrkpl, ExprVector ae, ExprVector be);
     static Expr *Distance(hEntity workplane, hEntity pa, hEntity pb);
     static Expr *PointLineDistance(hEntity workplane, hEntity pt, hEntity ln);
@@ -685,14 +684,14 @@ public:
     Vector GetReferencePos();
     void Draw();
     void GetEdges(SEdgeList *sel);
-    bool IsStylable();
+    bool IsStylable() const;
     hStyle GetStyle() const;
-    bool HasLabel();
+    bool HasLabel() const;
 
     void LineDrawOrGetDistance(Vector a, Vector b);
     bool IsVisible() const;
     void DrawOrGetDistance(Vector *labelPos);
-    std::string Label();
+    std::string Label() const;
     bool DoLineExtend(Vector p0, Vector p1, Vector pt, double salient);
     void DoArcForAngle(Vector a0, Vector da, Vector b0, Vector db,
                         Vector offset, Vector *ref, bool trim);
@@ -706,7 +705,7 @@ public:
     void DoEqualLenTicks(Vector a, Vector b, Vector gn);
     void DoEqualRadiusTicks(hEntity he);
 
-    std::string DescriptionString();
+    std::string DescriptionString() const;
 
     static hConstraint AddConstraint(Constraint *c, bool rememberForUndo);
     static hConstraint AddConstraint(Constraint *c);
@@ -724,8 +723,8 @@ class hEquation {
 public:
     uint32_t v;
 
-    inline bool isFromConstraint();
-    inline hConstraint constraint();
+    inline bool isFromConstraint() const;
+    inline hConstraint constraint() const;
 };
 
 class Equation {
@@ -850,49 +849,49 @@ public:
     static int PatternType(hStyle hs);
     static double StippleScaleMm(hStyle hs);
 
-    std::string DescriptionString();
+    std::string DescriptionString() const;
 
     void Clear() {}
 };
 
 
-inline hEntity hGroup::entity(int i)
+inline hEntity hGroup::entity(int i) const
     { hEntity r; r.v = 0x80000000 | (v << 16) | (uint32_t)i; return r; }
-inline hParam hGroup::param(int i)
+inline hParam hGroup::param(int i) const
     { hParam r; r.v = 0x80000000 | (v << 16) | (uint32_t)i; return r; }
-inline hEquation hGroup::equation(int i)
+inline hEquation hGroup::equation(int i) const
     { hEquation r; r.v = (v << 16) | 0x80000000 | (uint32_t)i; return r; }
 
-inline bool hRequest::IsFromReferences() {
+inline bool hRequest::IsFromReferences() const {
     if(v == Request::HREQUEST_REFERENCE_XY.v) return true;
     if(v == Request::HREQUEST_REFERENCE_YZ.v) return true;
     if(v == Request::HREQUEST_REFERENCE_ZX.v) return true;
     return false;
 }
-inline hEntity hRequest::entity(int i)
+inline hEntity hRequest::entity(int i) const
     { hEntity r; r.v = (v << 16) | (uint32_t)i; return r; }
-inline hParam hRequest::param(int i)
+inline hParam hRequest::param(int i) const
     { hParam r; r.v = (v << 16) | (uint32_t)i; return r; }
 
-inline bool hEntity::isFromRequest()
+inline bool hEntity::isFromRequest() const
     { if(v & 0x80000000) return false; else return true; }
-inline hRequest hEntity::request()
+inline hRequest hEntity::request() const
     { hRequest r; r.v = (v >> 16); return r; }
-inline hGroup hEntity::group()
+inline hGroup hEntity::group() const
     { hGroup r; r.v = (v >> 16) & 0x3fff; return r; }
-inline hEquation hEntity::equation(int i)
+inline hEquation hEntity::equation(int i) const
     { hEquation r; r.v = v | 0x40000000; return r; }
 
-inline hRequest hParam::request()
+inline hRequest hParam::request() const
     { hRequest r; r.v = (v >> 16); return r; }
 
 
-inline hEquation hConstraint::equation(int i)
+inline hEquation hConstraint::equation(int i) const
     { hEquation r; r.v = (v << 16) | (uint32_t)i; return r; }
 
-inline bool hEquation::isFromConstraint()
+inline bool hEquation::isFromConstraint() const
     { if(v & 0xc0000000) return false; else return true; }
-inline hConstraint hEquation::constraint()
+inline hConstraint hEquation::constraint() const
     { hConstraint r; r.v = (v >> 16); return r; }
 
 // The format for entities stored on the clipboard.
diff --git a/src/solvespace.h b/src/solvespace.h
index de687ea..3c29344 100644
--- a/src/solvespace.h
+++ b/src/solvespace.h
@@ -521,7 +521,7 @@ public:
     void BezierAsPwl(SBezier *sb);
     void BezierAsNonrationalCubic(SBezier *sb, int depth=0);
 
-    virtual void StartPath( RgbaColor strokeRgb, double lineWidth,
+    virtual void StartPath(RgbaColor strokeRgb, double lineWidth,
                             bool filled, RgbaColor fillRgb, hStyle hs) = 0;
     virtual void FinishPath(RgbaColor strokeRgb, double lineWidth,
                             bool filled, RgbaColor fillRgb, hStyle hs) = 0;
@@ -856,8 +856,9 @@ public:
     void ExportWireframeCurves(SEdgeList *sel, SBezierList *sbl,
                                VectorFileWriter *out);
     void ExportLinesAndMesh(SEdgeList *sel, SBezierList *sbl, SMesh *sm,
-                            Vector u, Vector v, Vector n, Vector origin,
-                                double cameraTan,
+                            Vector u, Vector v,
+                            Vector n, Vector origin,
+                            double cameraTan,
                             VectorFileWriter *out);
 
     static void MenuAnalyze(int id);
@@ -974,7 +975,7 @@ void ImportDwg(const std::string &file);
 extern SolveSpaceUI SS;
 extern Sketch SK;
 
-};
+}
 
 #ifndef __OBJC__
 using namespace SolveSpace;
diff --git a/src/srf/boolean.cpp b/src/srf/boolean.cpp
index 5cc8152..ff84553 100644
--- a/src/srf/boolean.cpp
+++ b/src/srf/boolean.cpp
@@ -36,13 +36,13 @@ static int ByTAlongLine(const void *av, const void *bv)
     return (ta > tb) ? 1 : -1;
 }
 SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
-                                    SSurface *srfA, SSurface *srfB)
+                                    SSurface *srfA, SSurface *srfB) const
 {
     SCurve ret;
     ret = *this;
     ret.pts = {};
 
-    SCurvePt *p = pts.First();
+    const SCurvePt *p = pts.First();
     ssassert(p != NULL, "Cannot split an empty curve");
     SCurvePt prev = *p;
     ret.pts.Add(p);
@@ -364,7 +364,7 @@ void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv,
         ClosestPointTo(*pt, &muv);
     } else if(!sc->isExact) {
         SSurface *trimmedA = sc->GetSurfaceA(sha, shb),
-                 *trimmedB = sc->GetSurfaceB(sha, shb);
+                       *trimmedB = sc->GetSurfaceB(sha, shb);
         *pt = trimmedA->ClosestPointOnThisAndSurface(trimmedB, *pt);
         ClosestPointTo(*pt, &muv);
     }
@@ -602,9 +602,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
     return ret;
 }
 
-void SShell::CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into,
-                                        int type)
-{
+void SShell::CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into, int type) {
     SSurface *ss;
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
         SSurface ssn;
@@ -716,7 +714,7 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
     SCurve *sc;
     for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) {
         SSurface *srfA = sc->GetSurfaceA(a, b),
-                 *srfB = sc->GetSurfaceB(a, b);
+                       *srfB = sc->GetSurfaceB(a, b);
 
         sc->RemoveShortSegments(srfA, srfB);
     }
@@ -785,6 +783,7 @@ static int ByLength(const void *av, const void *bv)
     // stability for the normals.
     return (la < lb) ? 1 : -1;
 }
+
 SBspUv *SBspUv::From(SEdgeList *el, SSurface *srf) {
     SEdgeList work = {};
 
@@ -812,7 +811,8 @@ SBspUv *SBspUv::From(SEdgeList *el, SSurface *srf) {
 // time we care about exact correctness is when we're very close to the line,
 // which is when the linearization is accurate.
 //-----------------------------------------------------------------------------
-void SBspUv::ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf) {
+
+void SBspUv::ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf) const {
     Vector tu, tv;
     srf->TangentsAt(pt->x, pt->y, &tu, &tv);
     double mu = tu.Magnitude(), mv = tv.Magnitude();
@@ -821,8 +821,9 @@ void SBspUv::ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf) {
     a ->x *= mu; a ->y *= mv;
     b ->x *= mu; b ->y *= mv;
 }
+
 double SBspUv::ScaledSignedDistanceToLine(Point2d pt, Point2d a, Point2d b,
-                                            SSurface *srf)
+                                          SSurface *srf) const
 {
     ScalePoints(&pt, &a, &b, srf);
 
@@ -831,15 +832,16 @@ double SBspUv::ScaledSignedDistanceToLine(Point2d pt, Point2d a, Point2d b,
 
     return pt.Dot(n) - d;
 }
+
 double SBspUv::ScaledDistanceToLine(Point2d pt, Point2d a, Point2d b, bool seg,
-                                        SSurface *srf)
+                                    SSurface *srf) const
 {
     ScalePoints(&pt, &a, &b, srf);
 
     return pt.DistanceToLine(a, b, seg);
 }
 
-SBspUv *SBspUv::InsertOrCreateEdge(SBspUv *where, const Point2d &ea, const Point2d &eb, SSurface *srf) {
+SBspUv *SBspUv::InsertOrCreateEdge(SBspUv *where, Point2d ea, Point2d eb, SSurface *srf) {
     if(where == NULL) {
         SBspUv *ret = Alloc();
         ret->a = ea;
@@ -896,12 +898,11 @@ void SBspUv::InsertEdge(Point2d ea, Point2d eb, SSurface *srf) {
     return;
 }
 
-int SBspUv::ClassifyPoint(Point2d p, Point2d eb, SSurface *srf) {
-
+int SBspUv::ClassifyPoint(Point2d p, Point2d eb, SSurface *srf) const {
     double dp = ScaledSignedDistanceToLine(p, a, b, srf);
 
     if(fabs(dp) < LENGTH_EPS) {
-        SBspUv *f = this;
+        const SBspUv *f = this;
         while(f) {
             Point2d ba = (f->b).Minus(f->a);
             if(ScaledDistanceToLine(p, f->a, ba, true, srf) < LENGTH_EPS) {
@@ -931,7 +932,7 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb, SSurface *srf) {
     }
 }
 
-int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) {
+int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) const {
     int ret = ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb, srf);
     if(ret == EDGE_OTHER) {
         // Perhaps the edge is tangent at its midpoint (and we screwed up
@@ -942,7 +943,7 @@ int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) {
     return ret;
 }
 
-double SBspUv::MinimumDistanceToEdge(Point2d p, SSurface *srf) {
+double SBspUv::MinimumDistanceToEdge(Point2d p, SSurface *srf) const {
 
     double dn = (neg) ? neg->MinimumDistanceToEdge(p, srf) : VERY_POSITIVE;
     double dp = (pos) ? pos->MinimumDistanceToEdge(p, srf) : VERY_POSITIVE;
diff --git a/src/srf/curve.cpp b/src/srf/curve.cpp
index ae38c19..50f2da9 100644
--- a/src/srf/curve.cpp
+++ b/src/srf/curve.cpp
@@ -60,11 +60,11 @@ SBezier SBezier::From(Vector p0, Vector p1, Vector p2, Vector p3) {
                          p3.Project4d());
 }
 
-Vector SBezier::Start() {
+Vector SBezier::Start() const {
     return ctrl[0];
 }
 
-Vector SBezier::Finish() {
+Vector SBezier::Finish() const {
     return ctrl[deg];
 }
 
@@ -84,7 +84,7 @@ void SBezier::ScaleSelfBy(double s) {
 }
 
 void SBezier::GetBoundingProjd(Vector u, Vector orig,
-                               double *umin, double *umax)
+                               double *umin, double *umax) const
 {
     int i;
     for(i = 0; i <= deg; i++) {
@@ -94,7 +94,7 @@ void SBezier::GetBoundingProjd(Vector u, Vector orig,
     }
 }
 
-SBezier SBezier::TransformedBy(Vector t, Quaternion q, double scale) {
+SBezier SBezier::TransformedBy(Vector t, Quaternion q, double scale) const {
     SBezier ret = *this;
     int i;
     for(i = 0; i <= deg; i++) {
@@ -108,7 +108,7 @@ SBezier SBezier::TransformedBy(Vector t, Quaternion q, double scale) {
 // Does this curve lie entirely within the specified plane? It does if all
 // the control points lie in that plane.
 //-----------------------------------------------------------------------------
-bool SBezier::IsInPlane(Vector n, double d) {
+bool SBezier::IsInPlane(Vector n, double d) const {
     int i;
     for(i = 0; i <= deg; i++) {
         if(fabs((ctrl[i]).Dot(n) - d) > LENGTH_EPS) {
@@ -122,7 +122,7 @@ bool SBezier::IsInPlane(Vector n, double d) {
 // Is this Bezier exactly the arc of a circle, projected along the specified
 // axis? If yes, return that circle's center and radius.
 //-----------------------------------------------------------------------------
-bool SBezier::IsCircle(Vector axis, Vector *center, double *r) {
+bool SBezier::IsCircle(Vector axis, Vector *center, double *r) const {
     if(deg != 2) return false;
 
     if(ctrl[1].DistanceToLine(ctrl[0], ctrl[2].Minus(ctrl[0])) < LENGTH_EPS) {
@@ -170,7 +170,7 @@ bool SBezier::IsCircle(Vector axis, Vector *center, double *r) {
     return true;
 }
 
-bool SBezier::IsRational() {
+bool SBezier::IsRational() const {
     int i;
     for(i = 0; i <= deg; i++) {
         if(fabs(weight[i] - 1) > LENGTH_EPS) return true;
@@ -183,7 +183,7 @@ bool SBezier::IsRational() {
 // the new weights as required.
 //-----------------------------------------------------------------------------
 SBezier SBezier::InPerspective(Vector u, Vector v, Vector n,
-                               Vector origin, double cameraTan)
+                               Vector origin, double cameraTan) const
 {
     Quaternion q = Quaternion::From(u, v);
     q = q.Inverse();
@@ -206,7 +206,7 @@ SBezier SBezier::InPerspective(Vector u, Vector v, Vector n,
     return ret;
 }
 
-bool SBezier::Equals(SBezier *b) {
+bool SBezier::Equals(SBezier *b) const {
     // We just test of identical degree and control points, even though two
     // curves could still be coincident (even sharing endpoints).
     if(deg != b->deg) return false;
@@ -262,15 +262,14 @@ void SBezierList::CullIdenticalBeziers() {
 // curves. So this will screw up on tangencies and stuff, but otherwise should
 // be fine.
 //-----------------------------------------------------------------------------
-void SBezierList::AllIntersectionsWith(SBezierList *sblb, SPointList *spl) {
-    SBezier *sba, *sbb;
-    for(sba = l.First(); sba; sba = l.NextAfter(sba)) {
-        for(sbb = sblb->l.First(); sbb; sbb = sblb->l.NextAfter(sbb)) {
+void SBezierList::AllIntersectionsWith(SBezierList *sblb, SPointList *spl) const {
+    for(const SBezier *sba = l.First(); sba; sba = l.NextAfter(sba)) {
+        for(const SBezier *sbb = sblb->l.First(); sbb; sbb = sblb->l.NextAfter(sbb)) {
             sbb->AllIntersectionsWith(sba, spl);
         }
     }
 }
-void SBezier::AllIntersectionsWith(SBezier *sbb, SPointList *spl) {
+void SBezier::AllIntersectionsWith(const SBezier *sbb, SPointList *spl) const {
     SPointList splRaw = {};
     SEdgeList sea, seb;
     sea = {};
@@ -304,12 +303,11 @@ void SBezier::AllIntersectionsWith(SBezier *sbb, SPointList *spl) {
 // Returns true if all the curves are coplanar, otherwise false.
 //-----------------------------------------------------------------------------
 bool SBezierList::GetPlaneContainingBeziers(Vector *p, Vector *u, Vector *v,
-                        Vector *notCoplanarAt)
+                        Vector *notCoplanarAt) const
 {
     Vector pt, ptFar, ptOffLine, dp, n;
     double farMax, offLineMax;
     int i;
-    SBezier *sb;
 
     // Get any point on any Bezier; or an arbitrary point if list is empty.
     if(l.n > 0) {
@@ -321,7 +319,7 @@ bool SBezierList::GetPlaneContainingBeziers(Vector *p, Vector *u, Vector *v,
 
     // Get the point farthest from our arbitrary point.
     farMax = VERY_NEGATIVE;
-    for(sb = l.First(); sb; sb = l.NextAfter(sb)) {
+    for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {
         for(i = 0; i <= sb->deg; i++) {
             double m = (pt.Minus(sb->ctrl[i])).Magnitude();
             if(m > farMax) {
@@ -341,7 +339,7 @@ bool SBezierList::GetPlaneContainingBeziers(Vector *p, Vector *u, Vector *v,
     // Get the point farthest from the line between pt and ptFar
     dp = ptFar.Minus(pt);
     offLineMax = VERY_NEGATIVE;
-    for(sb = l.First(); sb; sb = l.NextAfter(sb)) {
+    for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {
         for(i = 0; i <= sb->deg; i++) {
             double m = (sb->ctrl[i]).DistanceToLine(pt, dp);
             if(m > offLineMax) {
@@ -369,7 +367,7 @@ bool SBezierList::GetPlaneContainingBeziers(Vector *p, Vector *u, Vector *v,
     n = u->Cross(*v);
     n = n.WithMagnitude(1);
     double d = p->Dot(n);
-    for(sb = l.First(); sb; sb = l.NextAfter(sb)) {
+    for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {
         for(i = 0; i <= sb->deg; i++) {
             if(fabs(n.Dot(sb->ctrl[i]) - d) > LENGTH_EPS) {
                 if(notCoplanarAt) *notCoplanarAt = sb->ctrl[i];
@@ -454,17 +452,15 @@ void SBezierLoop::Reverse() {
 }
 
 void SBezierLoop::GetBoundingProjd(Vector u, Vector orig,
-                                   double *umin, double *umax)
+                                   double *umin, double *umax) const
 {
-    SBezier *sb;
-    for(sb = l.First(); sb; sb = l.NextAfter(sb)) {
+    for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {
         sb->GetBoundingProjd(u, orig, umin, umax);
     }
 }
 
-void SBezierLoop::MakePwlInto(SContour *sc, double chordTol) {
-    SBezier *sb;
-    for(sb = l.First(); sb; sb = l.NextAfter(sb)) {
+void SBezierLoop::MakePwlInto(SContour *sc, double chordTol) const {
+    for(const SBezier *sb = l.First(); sb; sb = l.NextAfter(sb)) {
         sb->MakePwlInto(sc, chordTol);
         // Avoid double points at join between Beziers; except that
         // first and last points should be identical.
@@ -478,7 +474,7 @@ void SBezierLoop::MakePwlInto(SContour *sc, double chordTol) {
     }
 }
 
-bool SBezierLoop::IsClosed() {
+bool SBezierLoop::IsClosed() const {
     if(l.n < 1) return false;
     Vector s = l.elem[0].Start(),
            f = l.elem[l.n-1].Finish();
@@ -533,10 +529,9 @@ SBezierLoopSet SBezierLoopSet::From(SBezierList *sbl, SPolygon *poly,
 }
 
 void SBezierLoopSet::GetBoundingProjd(Vector u, Vector orig,
-                                      double *umin, double *umax)
+                                      double *umin, double *umax) const
 {
-    SBezierLoop *sbl;
-    for(sbl = l.First(); sbl; sbl = l.NextAfter(sbl)) {
+    for(const SBezierLoop *sbl = l.First(); sbl; sbl = l.NextAfter(sbl)) {
         sbl->GetBoundingProjd(u, orig, umin, umax);
     }
 }
@@ -545,9 +540,8 @@ void SBezierLoopSet::GetBoundingProjd(Vector u, Vector orig,
 // Convert all the Beziers into piecewise linear form, and assemble that into
 // a polygon, one contour per loop.
 //-----------------------------------------------------------------------------
-void SBezierLoopSet::MakePwlInto(SPolygon *sp) {
-    SBezierLoop *sbl;
-    for(sbl = l.First(); sbl; sbl = l.NextAfter(sbl)) {
+void SBezierLoopSet::MakePwlInto(SPolygon *sp) const {
+    for(const SBezierLoop *sbl = l.First(); sbl; sbl = l.NextAfter(sbl)) {
         sp->AddEmptyContour();
         sbl->MakePwlInto(&(sp->l.elem[sp->l.n - 1]));
     }
@@ -732,8 +726,8 @@ void SBezierLoopSetSet::Clear() {
     l.Clear();
 }
 
-SCurve SCurve::FromTransformationOf(SCurve *a,
-                                          Vector t, Quaternion q, double scale)
+SCurve SCurve::FromTransformationOf(SCurve *a, Vector t,
+                                    Quaternion q, double scale)
 {
     SCurve ret = {};
 
@@ -757,7 +751,7 @@ void SCurve::Clear() {
     pts.Clear();
 }
 
-SSurface *SCurve::GetSurfaceA(SShell *a, SShell *b) {
+SSurface *SCurve::GetSurfaceA(SShell *a, SShell *b) const {
     if(source == FROM_A) {
         return a->surface.FindById(surfA);
     } else if(source == FROM_B) {
@@ -767,7 +761,7 @@ SSurface *SCurve::GetSurfaceA(SShell *a, SShell *b) {
     } else ssassert(false, "Unexpected curve source");
 }
 
-SSurface *SCurve::GetSurfaceB(SShell *a, SShell *b) {
+SSurface *SCurve::GetSurfaceB(SShell *a, SShell *b) const {
     if(source == FROM_A) {
         return a->surface.FindById(surfB);
     } else if(source == FROM_B) {
diff --git a/src/srf/ratpoly.cpp b/src/srf/ratpoly.cpp
index bfbb951..7060d20 100644
--- a/src/srf/ratpoly.cpp
+++ b/src/srf/ratpoly.cpp
@@ -93,7 +93,7 @@ double SolveSpace::BernsteinDerivative(int k, int deg, double t)
     ssassert(false, "Unexpected degree of spline");
 }
 
-Vector SBezier::PointAt(double t) {
+Vector SBezier::PointAt(double t) const {
     Vector pt = Vector::From(0, 0, 0);
     double d = 0;
 
@@ -107,7 +107,7 @@ Vector SBezier::PointAt(double t) {
     return pt;
 }
 
-Vector SBezier::TangentAt(double t) {
+Vector SBezier::TangentAt(double t) const {
     Vector pt = Vector::From(0, 0, 0), pt_p = Vector::From(0, 0, 0);
     double d = 0, d_p = 0;
 
@@ -130,7 +130,7 @@ Vector SBezier::TangentAt(double t) {
     return ret;
 }
 
-void SBezier::ClosestPointTo(Vector p, double *t, bool converge) {
+void SBezier::ClosestPointTo(Vector p, double *t, bool converge) const {
     int i;
     double minDist = VERY_POSITIVE;
     *t = 0;
@@ -162,7 +162,7 @@ void SBezier::ClosestPointTo(Vector p, double *t, bool converge) {
     }
 }
 
-bool SBezier::PointOnThisAndCurve(SBezier *sbb, Vector *p) {
+bool SBezier::PointOnThisAndCurve(const SBezier *sbb, Vector *p) const {
     double ta, tb;
     this->ClosestPointTo(*p, &ta, false);
     sbb ->ClosestPointTo(*p, &tb, false);
@@ -187,7 +187,7 @@ bool SBezier::PointOnThisAndCurve(SBezier *sbb, Vector *p) {
     return false;
 }
 
-void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) {
+void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) const {
     Vector4 ct[4];
     int i;
     for(i = 0; i <= deg; i++) {
@@ -226,7 +226,7 @@ void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) {
     }
 }
 
-void SBezier::MakePwlInto(SEdgeList *sel, double chordTol) {
+void SBezier::MakePwlInto(SEdgeList *sel, double chordTol) const {
     List<Vector> lv = {};
     MakePwlInto(&lv, chordTol);
     int i;
@@ -235,7 +235,7 @@ void SBezier::MakePwlInto(SEdgeList *sel, double chordTol) {
     }
     lv.Clear();
 }
-void SBezier::MakePwlInto(List<SCurvePt> *l, double chordTol) {
+void SBezier::MakePwlInto(List<SCurvePt> *l, double chordTol) const {
     List<Vector> lv = {};
     MakePwlInto(&lv, chordTol);
     int i;
@@ -248,7 +248,7 @@ void SBezier::MakePwlInto(List<SCurvePt> *l, double chordTol) {
     }
     lv.Clear();
 }
-void SBezier::MakePwlInto(SContour *sc, double chordTol) {
+void SBezier::MakePwlInto(SContour *sc, double chordTol) const {
     List<Vector> lv = {};
     MakePwlInto(&lv, chordTol);
     int i;
@@ -257,7 +257,7 @@ void SBezier::MakePwlInto(SContour *sc, double chordTol) {
     }
     lv.Clear();
 }
-void SBezier::MakePwlInto(List<Vector> *l, double chordTol) {
+void SBezier::MakePwlInto(List<Vector> *l, double chordTol) const {
     if(EXACT(chordTol == 0)) {
         // Use the default chord tolerance.
         chordTol = SS.ChordTolMm();
@@ -273,7 +273,7 @@ void SBezier::MakePwlInto(List<Vector> *l, double chordTol) {
         MakePwlInitialWorker(l, 0.5, 1.0, chordTol);
     }
 }
-void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb, double chordTol)
+void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb, double chordTol) const
 {
     Vector pa = PointAt(ta);
     Vector pb = PointAt(tb);
@@ -291,7 +291,7 @@ void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb, double chordT
         MakePwlWorker(l, tm, tb, chordTol);
     }
 }
-void SBezier::MakePwlInitialWorker(List<Vector> *l, double ta, double tb, double chordTol)
+void SBezier::MakePwlInitialWorker(List<Vector> *l, double ta, double tb, double chordTol) const
 {
     Vector pa = PointAt(ta);
     Vector pb = PointAt(tb);
@@ -322,10 +322,10 @@ void SBezier::MakePwlInitialWorker(List<Vector> *l, double ta, double tb, double
     }
 }
 
-Vector SSurface::PointAt(Point2d puv) {
+Vector SSurface::PointAt(Point2d puv) const {
     return PointAt(puv.x, puv.y);
 }
-Vector SSurface::PointAt(double u, double v) {
+Vector SSurface::PointAt(double u, double v) const {
     Vector num = Vector::From(0, 0, 0);
     double den = 0;
 
@@ -343,7 +343,7 @@ Vector SSurface::PointAt(double u, double v) {
     return num;
 }
 
-void SSurface::TangentsAt(double u, double v, Vector *tu, Vector *tv) {
+void SSurface::TangentsAt(double u, double v, Vector *tu, Vector *tv) const {
     Vector num   = Vector::From(0, 0, 0),
            num_u = Vector::From(0, 0, 0),
            num_v = Vector::From(0, 0, 0);
@@ -377,10 +377,11 @@ void SSurface::TangentsAt(double u, double v, Vector *tu, Vector *tv) {
     *tv = tv->ScaledBy(1.0/(den*den));
 }
 
-Vector SSurface::NormalAt(Point2d puv) {
+Vector SSurface::NormalAt(Point2d puv) const {
     return NormalAt(puv.x, puv.y);
 }
-Vector SSurface::NormalAt(double u, double v) {
+
+Vector SSurface::NormalAt(double u, double v) const {
     Vector tu, tv;
     TangentsAt(u, v, &tu, &tv);
     return tu.Cross(tv);
@@ -389,6 +390,7 @@ Vector SSurface::NormalAt(double u, double v) {
 void SSurface::ClosestPointTo(Vector p, Point2d *puv, bool converge) {
     ClosestPointTo(p, &(puv->x), &(puv->y), converge);
 }
+
 void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
     // A few special cases first; when control points are coincident the
     // derivative goes to zero at the conrol points, and would result in
@@ -455,7 +457,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
     }
 }
 
-bool SSurface::ClosestPointNewton(Vector p, double *u, double *v, bool converge)
+bool SSurface::ClosestPointNewton(Vector p, double *u, double *v, bool converge) const
 {
     // Initial guess is in u, v; refine by Newton iteration.
     Vector p0 = Vector::From(0, 0, 0);
@@ -488,7 +490,7 @@ bool SSurface::ClosestPointNewton(Vector p, double *u, double *v, bool converge)
     return false;
 }
 
-bool SSurface::PointIntersectingLine(Vector p0, Vector p1, double *u, double *v)
+bool SSurface::PointIntersectingLine(Vector p0, Vector p1, double *u, double *v) const
 {
     int i;
     for(i = 0; i < 15; i++) {
@@ -564,8 +566,7 @@ Vector SSurface::ClosestPointOnThisAndSurface(SSurface *srf2, Vector p) {
            ((srf[1])->PointAt(puv[1]))).ScaledBy(0.5);
 }
 
-void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2,
-                                                double *up, double *vp)
+void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2, double *up, double *vp)
 {
     double u[3] = { *up, 0, 0 }, v[3] = { *vp, 0, 0 };
     SSurface *srf[3] = { this, s1, s2 };
diff --git a/src/srf/raycast.cpp b/src/srf/raycast.cpp
index 17d7742..75c1809 100644
--- a/src/srf/raycast.cpp
+++ b/src/srf/raycast.cpp
@@ -16,7 +16,7 @@ const double SShell::DOTP_TOL = 1e-5;
 extern int FLAG;
 
 
-double SSurface::DepartureFromCoplanar() {
+double SSurface::DepartureFromCoplanar() const {
     int i, j;
     int ia, ja, ib = 0, jb = 0, ic = 0, jc = 0;
     double best;
@@ -392,7 +392,7 @@ void SShell::AllPointsIntersecting(Vector a, Vector b,
 
 
 int SShell::ClassifyRegion(Vector edge_n, Vector inter_surf_n,
-                           Vector edge_surf_n)
+                           Vector edge_surf_n) const
 {
     double dot = inter_surf_n.DirectionCosineWith(edge_n);
     if(fabs(dot) < DOTP_TOL) {
diff --git a/src/srf/surface.cpp b/src/srf/surface.cpp
index e1d1c4d..c702165 100644
--- a/src/srf/surface.cpp
+++ b/src/srf/surface.cpp
@@ -24,7 +24,7 @@ SSurface SSurface::FromExtrusionOf(SBezier *sb, Vector t0, Vector t1) {
     return ret;
 }
 
-bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) {
+bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) const {
     int i;
 
     if(degn != 1) return false;
@@ -52,7 +52,7 @@ bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) {
 }
 
 bool SSurface::IsCylinder(Vector *axis, Vector *center, double *r,
-                            Vector *start, Vector *finish)
+                            Vector *start, Vector *finish) const
 {
     SBezier sb;
     if(!IsExtrusion(&sb, axis)) return false;
@@ -175,7 +175,7 @@ SSurface SSurface::FromTransformationOf(SSurface *a,
     return ret;
 }
 
-void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) {
+void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) const {
     *ptMax = Vector::From(VERY_NEGATIVE, VERY_NEGATIVE, VERY_NEGATIVE);
     *ptMin = Vector::From(VERY_POSITIVE, VERY_POSITIVE, VERY_POSITIVE);
 
@@ -187,7 +187,7 @@ void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) {
     }
 }
 
-bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) {
+bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) const {
     Vector amax, amin;
     GetAxisAlignedBounding(&amax, &amin);
     if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, segment)) {
@@ -275,8 +275,7 @@ void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags,
 // by taking the cross product of the surface normals. We choose the direction
 // of this tangent so that its dot product with dir is positive.
 //-----------------------------------------------------------------------------
-Vector SSurface::ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB,
-                                       Vector dir)
+Vector SSurface::ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB, Vector dir)
 {
     Point2d puva, puvb;
     srfA->ClosestPointTo(p, &puva);
@@ -295,8 +294,7 @@ Vector SSurface::ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB,
 // add its exact form to sbl. Otherwise, add its piecewise linearization to
 // sel.
 //-----------------------------------------------------------------------------
-void SSurface::MakeSectionEdgesInto(SShell *shell,
-                                    SEdgeList *sel, SBezierList *sbl)
+void SSurface::MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl)
 {
     STrimBy *stb;
     for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
@@ -324,8 +322,8 @@ void SSurface::MakeSectionEdgesInto(SShell *shell,
             sbl->l.Add(&keep_bef);
         } else if(sbl && !sel && !sc->isExact) {
             // We must approximate this trim curve, as piecewise cubic sections.
-            SSurface *srfA = shell->surface.FindById(sc->surfA),
-                     *srfB = shell->surface.FindById(sc->surfB);
+            SSurface *srfA = shell->surface.FindById(sc->surfA);
+            SSurface *srfB = shell->surface.FindById(sc->surfB);
 
             Vector s = stb->backwards ? stb->finish : stb->start,
                    f = stb->backwards ? stb->start : stb->finish;
@@ -842,8 +840,7 @@ void SShell::MakeEdgesInto(SEdgeList *sel) {
     }
 }
 
-void SShell::MakeSectionEdgesInto(Vector n, double d,
-                                 SEdgeList *sel, SBezierList *sbl)
+void SShell::MakeSectionEdgesInto(Vector n, double d, SEdgeList *sel, SBezierList *sbl)
 {
     SSurface *s;
     for(s = surface.First(); s; s = surface.NextAfter(s)) {
@@ -860,7 +857,7 @@ void SShell::TriangulateInto(SMesh *sm) {
     }
 }
 
-bool SShell::IsEmpty() {
+bool SShell::IsEmpty() const {
     return (surface.n == 0);
 }
 
diff --git a/src/srf/surface.h b/src/srf/surface.h
index 22a226b..27a71f0 100644
--- a/src/srf/surface.h
+++ b/src/srf/surface.h
@@ -39,17 +39,17 @@ public:
     static SBspUv *Alloc();
     static SBspUv *From(SEdgeList *el, SSurface *srf);
 
-    void ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf);
+    void ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf) const;
     double ScaledSignedDistanceToLine(Point2d pt, Point2d a, Point2d b,
-        SSurface *srf);
+        SSurface *srf) const;
     double ScaledDistanceToLine(Point2d pt, Point2d a, Point2d b, bool seg,
-        SSurface *srf);
+        SSurface *srf) const;
 
     void InsertEdge(Point2d a, Point2d b, SSurface *srf);
-    static SBspUv *InsertOrCreateEdge(SBspUv *where, const Point2d &ea, const Point2d &eb, SSurface *srf);
-    int ClassifyPoint(Point2d p, Point2d eb, SSurface *srf);
-    int ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf);
-    double MinimumDistanceToEdge(Point2d p, SSurface *srf);
+    static SBspUv *InsertOrCreateEdge(SBspUv *where, Point2d ea, Point2d eb, SSurface *srf);
+    int ClassifyPoint(Point2d p, Point2d eb, SSurface *srf) const;
+    int ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) const;
+    double MinimumDistanceToEdge(Point2d p, SSurface *srf) const;
 };
 
 // Now the data structures to represent a shell of trimmed rational polynomial
@@ -80,33 +80,33 @@ public:
     double          weight[4];
     uint32_t        entity;
 
-    Vector PointAt(double t);
-    Vector TangentAt(double t);
-    void ClosestPointTo(Vector p, double *t, bool converge=true);
-    void SplitAt(double t, SBezier *bef, SBezier *aft);
-    bool PointOnThisAndCurve(SBezier *sbb, Vector *p);
+    Vector PointAt(double t) const;
+    Vector TangentAt(double t) const;
+    void ClosestPointTo(Vector p, double *t, bool converge=true) const;
+    void SplitAt(double t, SBezier *bef, SBezier *aft) const;
+    bool PointOnThisAndCurve(const SBezier *sbb, Vector *p) const;
 
-    Vector Start();
-    Vector Finish();
-    bool Equals(SBezier *b);
-    void MakePwlInto(SEdgeList *sel, double chordTol=0);
-    void MakePwlInto(List<SCurvePt> *l, double chordTol=0);
-    void MakePwlInto(SContour *sc, double chordTol=0);
-    void MakePwlInto(List<Vector> *l, double chordTol=0);
-    void MakePwlWorker(List<Vector> *l, double ta, double tb, double chordTol);
-    void MakePwlInitialWorker(List<Vector> *l, double ta, double tb, double chordTol);
+    Vector Start() const;
+    Vector Finish() const;
+    bool Equals(SBezier *b) const;
+    void MakePwlInto(SEdgeList *sel, double chordTol=0) const;
+    void MakePwlInto(List<SCurvePt> *l, double chordTol=0) const;
+    void MakePwlInto(SContour *sc, double chordTol=0) const;
+    void MakePwlInto(List<Vector> *l, double chordTol=0) const;
+    void MakePwlWorker(List<Vector> *l, double ta, double tb, double chordTol) const;
+    void MakePwlInitialWorker(List<Vector> *l, double ta, double tb, double chordTol) const;
 
-    void AllIntersectionsWith(SBezier *sbb, SPointList *spl);
-    void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax);
+    void AllIntersectionsWith(const SBezier *sbb, SPointList *spl) const;
+    void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax) const;
     void Reverse();
 
-    bool IsInPlane(Vector n, double d);
-    bool IsCircle(Vector axis, Vector *center, double *r);
-    bool IsRational();
+    bool IsInPlane(Vector n, double d) const;
+    bool IsCircle(Vector axis, Vector *center, double *r) const;
+    bool IsRational() const;
 
-    SBezier TransformedBy(Vector t, Quaternion q, double scale);
+    SBezier TransformedBy(Vector t, Quaternion q, double scale) const;
     SBezier InPerspective(Vector u, Vector v, Vector n,
-                          Vector origin, double cameraTan);
+                          Vector origin, double cameraTan) const;
     void ScaleSelfBy(double s);
 
     static SBezier From(Vector p0, Vector p1, Vector p2, Vector p3);
@@ -124,9 +124,9 @@ public:
     void Clear();
     void ScaleSelfBy(double s);
     void CullIdenticalBeziers();
-    void AllIntersectionsWith(SBezierList *sblb, SPointList *spl);
+    void AllIntersectionsWith(SBezierList *sblb, SPointList *spl) const;
     bool GetPlaneContainingBeziers(Vector *p, Vector *u, Vector *v,
-                                        Vector *notCoplanarAt);
+                                        Vector *notCoplanarAt) const;
 };
 
 class SBezierLoop {
@@ -135,10 +135,10 @@ public:
     List<SBezier>   l;
 
     inline void Clear() { l.Clear(); }
-    bool IsClosed();
+    bool IsClosed() const;
     void Reverse();
-    void MakePwlInto(SContour *sc, double chordTol=0);
-    void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax);
+    void MakePwlInto(SContour *sc, double chordTol=0) const;
+    void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax) const;
 
     static SBezierLoop FromCurves(SBezierList *spcl,
                                   bool *allClosed, SEdge *errorAt);
@@ -155,8 +155,8 @@ public:
                                bool *allClosed, SEdge *errorAt,
                                SBezierList *openContours);
 
-    void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax);
-    void MakePwlInto(SPolygon *sp);
+    void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax) const;
+    void MakePwlInto(SPolygon *sp) const;
     void Clear();
 };
 
@@ -204,13 +204,13 @@ public:
     hSSurface       surfA;
     hSSurface       surfB;
 
-    static SCurve FromTransformationOf(SCurve *a, Vector t, Quaternion q,
-                                        double scale);
+    static SCurve FromTransformationOf(SCurve *a, Vector t,
+                                       Quaternion q, double scale);
     SCurve MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
-                                SSurface *srfA, SSurface *srfB);
+                                SSurface *srfA, SSurface *srfB) const;
     void RemoveShortSegments(SSurface *srfA, SSurface *srfB);
-    SSurface *GetSurfaceA(SShell *a, SShell *b);
-    SSurface *GetSurfaceB(SShell *a, SShell *b);
+    SSurface *GetSurfaceA(SShell *a, SShell *b) const;
+    SSurface *GetSurfaceB(SShell *a, SShell *b) const;
 
     void Clear();
 };
@@ -302,35 +302,35 @@ public:
     void CopyRowOrCol(bool row, int this_ij, SSurface *src, int src_ij);
     void BlendRowOrCol(bool row, int this_ij, SSurface *a, int a_ij,
                                               SSurface *b, int b_ij);
-    double DepartureFromCoplanar();
+    double DepartureFromCoplanar() const;
     void SplitInHalf(bool byU, SSurface *sa, SSurface *sb);
     void AllPointsIntersecting(Vector a, Vector b,
-                                    List<SInter> *l,
-                                    bool seg, bool trimmed, bool inclTangent);
+                               List<SInter> *l,
+                               bool seg, bool trimmed, bool inclTangent);
     void AllPointsIntersectingUntrimmed(Vector a, Vector b,
-                                            int *cnt, int *level,
-                                            List<Inter> *l, bool segment,
-                                            SSurface *sorig);
+                                        int *cnt, int *level,
+                                        List<Inter> *l, bool segment,
+                                        SSurface *sorig);
 
     void ClosestPointTo(Vector p, Point2d *puv, bool converge=true);
     void ClosestPointTo(Vector p, double *u, double *v, bool converge=true);
-    bool ClosestPointNewton(Vector p, double *u, double *v, bool converge=true);
+    bool ClosestPointNewton(Vector p, double *u, double *v, bool converge=true) const;
 
-    bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v);
+    bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v) const;
     Vector ClosestPointOnThisAndSurface(SSurface *srf2, Vector p);
     void PointOnSurfaces(SSurface *s1, SSurface *s2, double *u, double *v);
-    Vector PointAt(double u, double v);
-    Vector PointAt(Point2d puv);
-    void TangentsAt(double u, double v, Vector *tu, Vector *tv);
-    Vector NormalAt(Point2d puv);
-    Vector NormalAt(double u, double v);
-    bool LineEntirelyOutsideBbox(Vector a, Vector b, bool segment);
-    void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin);
-    bool CoincidentWithPlane(Vector n, double d);
-    bool CoincidentWith(SSurface *ss, bool sameNormal);
-    bool IsExtrusion(SBezier *of, Vector *along);
+    Vector PointAt(double u, double v) const;
+    Vector PointAt(Point2d puv) const;
+    void TangentsAt(double u, double v, Vector *tu, Vector *tv) const;
+    Vector NormalAt(Point2d puv) const;
+    Vector NormalAt(double u, double v) const;
+    bool LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) const;
+    void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) const;
+    bool CoincidentWithPlane(Vector n, double d) const;
+    bool CoincidentWith(SSurface *ss, bool sameNormal) const;
+    bool IsExtrusion(SBezier *of, Vector *along) const;
     bool IsCylinder(Vector *axis, Vector *center, double *r,
-                        Vector *start, Vector *finish);
+                        Vector *start, Vector *finish) const;
 
     void TriangulateInto(SShell *shell, SMesh *sm);
 
@@ -341,16 +341,16 @@ public:
     };
     void MakeTrimEdgesInto(SEdgeList *sel, int flags, SCurve *sc, STrimBy *stb);
     void MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags,
-            SShell *useCurvesFrom=NULL);
+                       SShell *useCurvesFrom=NULL);
 
     Vector ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB,
                                 Vector dir);
     void MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl);
     void MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom);
-    double ChordToleranceForEdge(Vector a, Vector b);
+    double ChordToleranceForEdge(Vector a, Vector b) const;
     void MakeTriangulationGridInto(List<double> *l, double vs, double vf,
-                                    bool swapped);
-    Vector PointAtMaybeSwapped(double u, double v, bool swapped);
+                                    bool swapped) const;
+    Vector PointAtMaybeSwapped(double u, double v, bool swapped) const;
 
     void Reverse();
     void Clear();
@@ -377,8 +377,7 @@ public:
     };
     void MakeFromBoolean(SShell *a, SShell *b, int type);
     void CopyCurvesSplitAgainst(bool opA, SShell *agnst, SShell *into);
-    void CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into,
-                                    int type);
+    void CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into, int type);
     void MakeIntersectionCurvesAgainst(SShell *against, SShell *into);
     void MakeClassifyingBsps(SShell *useCurvesFrom);
     void AllPointsIntersecting(Vector a, Vector b, List<SInter> *il,
@@ -398,23 +397,24 @@ public:
         COINC_OPP  = 400
     };
     static const double DOTP_TOL;
-    int ClassifyRegion(Vector edge_n, Vector inter_surf_n, Vector edge_surf_n);
+    int ClassifyRegion(Vector edge_n, Vector inter_surf_n,
+                       Vector edge_surf_n) const;
+
     bool ClassifyEdge(int *indir, int *outdir,
                       Vector ea, Vector eb,
-                      Vector p,
-                      Vector edge_n_in, Vector edge_n_out, Vector surf_n);
+                      Vector p, Vector edge_n_in,
+                      Vector edge_n_out, Vector surf_n);
 
     void MakeFromCopyOf(SShell *a);
     void MakeFromTransformationOf(SShell *a,
-                                    Vector trans, Quaternion q, double scale);
+                                  Vector trans, Quaternion q, double scale);
     void MakeFromAssemblyOf(SShell *a, SShell *b);
     void MergeCoincidentSurfaces();
 
     void TriangulateInto(SMesh *sm);
     void MakeEdgesInto(SEdgeList *sel);
-    void MakeSectionEdgesInto(Vector n, double d,
-                                SEdgeList *sel, SBezierList *sbl);
-    bool IsEmpty();
+    void MakeSectionEdgesInto(Vector n, double d, SEdgeList *sel, SBezierList *sbl);
+    bool IsEmpty() const;
     void RemapFaces(Group *g, int remap);
     void Clear();
 };
diff --git a/src/srf/surfinter.cpp b/src/srf/surfinter.cpp
index 998ce4c..641bcbc 100644
--- a/src/srf/surfinter.cpp
+++ b/src/srf/surfinter.cpp
@@ -10,7 +10,7 @@
 extern int FLAG;
 
 void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
-                            SShell *agnstA, SShell *agnstB, SShell *into)
+                                         SShell *agnstA, SShell *agnstB, SShell *into)
 {
     SCurve sc = {};
     // Important to keep the order of (surfA, surfB) consistent; when we later
@@ -460,7 +460,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
 // Are two surfaces coincident, with the same (or with opposite) normals?
 // Currently handles planes only.
 //-----------------------------------------------------------------------------
-bool SSurface::CoincidentWith(SSurface *ss, bool sameNormal) {
+bool SSurface::CoincidentWith(SSurface *ss, bool sameNormal) const {
     if(degm != 1 || degn != 1) return false;
     if(ss->degm != 1 || ss->degn != 1) return false;
 
@@ -480,7 +480,7 @@ bool SSurface::CoincidentWith(SSurface *ss, bool sameNormal) {
     return true;
 }
 
-bool SSurface::CoincidentWithPlane(Vector n, double d) {
+bool SSurface::CoincidentWithPlane(Vector n, double d) const {
     if(degm != 1 || degn != 1) return false;
     if(fabs(n.Dot(ctrl[0][0]) - d) > LENGTH_EPS) return false;
     if(fabs(n.Dot(ctrl[0][1]) - d) > LENGTH_EPS) return false;
diff --git a/src/srf/triangulate.cpp b/src/srf/triangulate.cpp
index 816eb3e..7473978 100644
--- a/src/srf/triangulate.cpp
+++ b/src/srf/triangulate.cpp
@@ -213,7 +213,7 @@ haveEdge:
     return true;
 }
 
-bool SContour::IsEar(int bp, double scaledEps) {
+bool SContour::IsEar(int bp, double scaledEps) const {
     int ap = WRAP(bp-1, l.n),
         cp = WRAP(bp+1, l.n);
 
@@ -359,7 +359,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
     ClipEarInto(m, 0, scaledEps); // add the last triangle
 }
 
-double SSurface::ChordToleranceForEdge(Vector a, Vector b) {
+double SSurface::ChordToleranceForEdge(Vector a, Vector b) const {
     Vector as = PointAt(a.x, a.y), bs = PointAt(b.x, b.y);
 
     double worst = VERY_NEGATIVE;
@@ -374,7 +374,7 @@ double SSurface::ChordToleranceForEdge(Vector a, Vector b) {
     return sqrt(worst);
 }
 
-Vector SSurface::PointAtMaybeSwapped(double u, double v, bool swapped) {
+Vector SSurface::PointAtMaybeSwapped(double u, double v, bool swapped) const {
     if(swapped) {
         return PointAt(v, u);
     } else {
@@ -383,7 +383,7 @@ Vector SSurface::PointAtMaybeSwapped(double u, double v, bool swapped) {
 }
 
 void SSurface::MakeTriangulationGridInto(List<double> *l, double vs, double vf,
-                                         bool swapped)
+                                         bool swapped) const
 {
     double worst = 0;
 
diff --git a/src/style.cpp b/src/style.cpp
index c9f3ea0..2f33721 100644
--- a/src/style.cpp
+++ b/src/style.cpp
@@ -337,7 +337,7 @@ double Style::StippleScaleMm(hStyle hs) {
     return 1.0;
 }
 
-std::string Style::DescriptionString() {
+std::string Style::DescriptionString() const {
     if(name.empty()) {
         return ssprintf("s%03x-(unnamed)", h.v);
     } else {
diff --git a/src/util.cpp b/src/util.cpp
index 9b13085..dec49d5 100644
--- a/src/util.cpp
+++ b/src/util.cpp
@@ -279,7 +279,7 @@ Quaternion Quaternion::From(Vector u, Vector v)
     return q.WithMagnitude(1);
 }
 
-Quaternion Quaternion::Plus(Quaternion b) {
+Quaternion Quaternion::Plus(Quaternion b) const {
     Quaternion q;
     q.w  = w  + b.w;
     q.vx = vx + b.vx;
@@ -288,7 +288,7 @@ Quaternion Quaternion::Plus(Quaternion b) {
     return q;
 }
 
-Quaternion Quaternion::Minus(Quaternion b) {
+Quaternion Quaternion::Minus(Quaternion b) const {
     Quaternion q;
     q.w  = w  - b.w;
     q.vx = vx - b.vx;
@@ -297,7 +297,7 @@ Quaternion Quaternion::Minus(Quaternion b) {
     return q;
 }
 
-Quaternion Quaternion::ScaledBy(double s) {
+Quaternion Quaternion::ScaledBy(double s) const {
     Quaternion q;
     q.w  = w*s;
     q.vx = vx*s;
@@ -306,15 +306,15 @@ Quaternion Quaternion::ScaledBy(double s) {
     return q;
 }
 
-double Quaternion::Magnitude() {
+double Quaternion::Magnitude() const {
     return sqrt(w*w + vx*vx + vy*vy + vz*vz);
 }
 
-Quaternion Quaternion::WithMagnitude(double s) {
+Quaternion Quaternion::WithMagnitude(double s) const {
     return ScaledBy(s/Magnitude());
 }
 
-Vector Quaternion::RotationU() {
+Vector Quaternion::RotationU() const {
     Vector v;
     v.x = w*w + vx*vx - vy*vy - vz*vz;
     v.y = 2*w *vz + 2*vx*vy;
@@ -322,7 +322,7 @@ Vector Quaternion::RotationU() {
     return v;
 }
 
-Vector Quaternion::RotationV() {
+Vector Quaternion::RotationV() const {
     Vector v;
     v.x = 2*vx*vy - 2*w*vz;
     v.y = w*w - vx*vx + vy*vy - vz*vz;
@@ -330,7 +330,7 @@ Vector Quaternion::RotationV() {
     return v;
 }
 
-Vector Quaternion::RotationN() {
+Vector Quaternion::RotationN() const {
     Vector v;
     v.x = 2*w*vy + 2*vx*vz;
     v.y = 2*vy*vz - 2*w*vx;
@@ -338,14 +338,14 @@ Vector Quaternion::RotationN() {
     return v;
 }
 
-Vector Quaternion::Rotate(Vector p) {
+Vector Quaternion::Rotate(Vector p) const {
     // Express the point in the new basis
     return (RotationU().ScaledBy(p.x)).Plus(
             RotationV().ScaledBy(p.y)).Plus(
             RotationN().ScaledBy(p.z));
 }
 
-Quaternion Quaternion::Inverse() {
+Quaternion Quaternion::Inverse() const {
     Quaternion r;
     r.w = w;
     r.vx = -vx;
@@ -354,7 +354,7 @@ Quaternion Quaternion::Inverse() {
     return r.WithMagnitude(1); // not that the normalize should be reqd
 }
 
-Quaternion Quaternion::ToThe(double p) {
+Quaternion Quaternion::ToThe(double p) const {
     // Avoid division by zero, or arccos of something not in its domain
     if(w >= (1 - 1e-6)) {
         return From(1, 0, 0, 0);
@@ -374,7 +374,7 @@ Quaternion Quaternion::ToThe(double p) {
     return r;
 }
 
-Quaternion Quaternion::Times(Quaternion b) {
+Quaternion Quaternion::Times(Quaternion b) const {
     double sa = w, sb = b.w;
     Vector va = { vx, vy, vz };
     Vector vb = { b.vx, b.vy, b.vz };
@@ -390,7 +390,7 @@ Quaternion Quaternion::Times(Quaternion b) {
     return r;
 }
 
-Quaternion Quaternion::Mirror() {
+Quaternion Quaternion::Mirror() const {
     Vector u = RotationU(),
            v = RotationV();
     u = u.ScaledBy(-1);
@@ -413,7 +413,7 @@ Vector Vector::From(hParam x, hParam y, hParam z) {
     return v;
 }
 
-double Vector::Element(int i) {
+double Vector::Element(int i) const {
     switch(i) {
         case 0: return x;
         case 1: return y;
@@ -422,7 +422,7 @@ double Vector::Element(int i) {
     }
 }
 
-bool Vector::Equals(Vector v, double tol) {
+bool Vector::Equals(Vector v, double tol) const {
     // Quick axis-aligned tests before going further
     double dx = v.x - x; if(dx < -tol || dx > tol) return false;
     double dy = v.y - y; if(dy < -tol || dy > tol) return false;
@@ -431,13 +431,13 @@ bool Vector::Equals(Vector v, double tol) {
     return (this->Minus(v)).MagSquared() < tol*tol;
 }
 
-bool Vector::EqualsExactly(Vector v) {
+bool Vector::EqualsExactly(Vector v) const {
     return EXACT(x == v.x &&
                  y == v.y &&
                  z == v.z);
 }
 
-Vector Vector::Plus(Vector b) {
+Vector Vector::Plus(Vector b) const {
     Vector r;
 
     r.x = x + b.x;
@@ -447,7 +447,7 @@ Vector Vector::Plus(Vector b) {
     return r;
 }
 
-Vector Vector::Minus(Vector b) {
+Vector Vector::Minus(Vector b) const {
     Vector r;
 
     r.x = x - b.x;
@@ -457,7 +457,7 @@ Vector Vector::Minus(Vector b) {
     return r;
 }
 
-Vector Vector::Negated() {
+Vector Vector::Negated() const {
     Vector r;
 
     r.x = -x;
@@ -467,7 +467,7 @@ Vector Vector::Negated() {
     return r;
 }
 
-Vector Vector::Cross(Vector b) {
+Vector Vector::Cross(Vector b) const {
     Vector r;
 
     r.x = -(z*b.y) + (y*b.z);
@@ -477,17 +477,17 @@ Vector Vector::Cross(Vector b) {
     return r;
 }
 
-double Vector::Dot(Vector b) {
+double Vector::Dot(Vector b) const {
     return (x*b.x + y*b.y + z*b.z);
 }
 
-double Vector::DirectionCosineWith(Vector b) {
+double Vector::DirectionCosineWith(Vector b) const {
     Vector a = this->WithMagnitude(1);
     b = b.WithMagnitude(1);
     return a.Dot(b);
 }
 
-Vector Vector::Normal(int which) {
+Vector Vector::Normal(int which) const {
     Vector n;
 
     // Arbitrarily choose one vector that's normal to us, pivoting
@@ -521,13 +521,13 @@ Vector Vector::Normal(int which) {
     return n;
 }
 
-Vector Vector::RotatedAbout(Vector orig, Vector axis, double theta) {
+Vector Vector::RotatedAbout(Vector orig, Vector axis, double theta) const {
     Vector r = this->Minus(orig);
     r = r.RotatedAbout(axis, theta);
     return r.Plus(orig);
 }
 
-Vector Vector::RotatedAbout(Vector axis, double theta) {
+Vector Vector::RotatedAbout(Vector axis, double theta) const {
     double c = cos(theta);
     double s = sin(theta);
 
@@ -550,7 +550,7 @@ Vector Vector::RotatedAbout(Vector axis, double theta) {
     return r;
 }
 
-Vector Vector::DotInToCsys(Vector u, Vector v, Vector n) {
+Vector Vector::DotInToCsys(Vector u, Vector v, Vector n) const {
     Vector r = {
         this->Dot(u),
         this->Dot(v),
@@ -559,7 +559,7 @@ Vector Vector::DotInToCsys(Vector u, Vector v, Vector n) {
     return r;
 }
 
-Vector Vector::ScaleOutOfCsys(Vector u, Vector v, Vector n) {
+Vector Vector::ScaleOutOfCsys(Vector u, Vector v, Vector n) const {
     Vector r = u.ScaledBy(x).Plus(
                v.ScaledBy(y).Plus(
                n.ScaledBy(z)));
@@ -567,7 +567,7 @@ Vector Vector::ScaleOutOfCsys(Vector u, Vector v, Vector n) {
 }
 
 Vector Vector::InPerspective(Vector u, Vector v, Vector n,
-                             Vector origin, double cameraTan)
+                             Vector origin, double cameraTan) const
 {
     Vector r = this->Minus(origin);
     r = r.DotInToCsys(u, v, n);
@@ -579,12 +579,12 @@ Vector Vector::InPerspective(Vector u, Vector v, Vector n,
     return r;
 }
 
-double Vector::DistanceToLine(Vector p0, Vector dp) {
+double Vector::DistanceToLine(Vector p0, Vector dp) const {
     double m = dp.Magnitude();
     return ((this->Minus(p0)).Cross(dp)).Magnitude() / m;
 }
 
-bool Vector::OnLineSegment(Vector a, Vector b, double tol) {
+bool Vector::OnLineSegment(Vector a, Vector b, double tol) const {
     if(this->Equals(a, tol) || this->Equals(b, tol)) return true;
 
     Vector d = b.Minus(a);
@@ -600,7 +600,7 @@ bool Vector::OnLineSegment(Vector a, Vector b, double tol) {
     return true;
 }
 
-Vector Vector::ClosestPointOnLine(Vector p0, Vector dp) {
+Vector Vector::ClosestPointOnLine(Vector p0, Vector dp) const {
     dp = dp.WithMagnitude(1);
     // this, p0, and (p0+dp) define a plane; the min distance is in
     // that plane, so calculate its normal
@@ -614,15 +614,15 @@ Vector Vector::ClosestPointOnLine(Vector p0, Vector dp) {
     return this->Plus(n.WithMagnitude(d));
 }
 
-double Vector::MagSquared() {
+double Vector::MagSquared() const {
     return x*x + y*y + z*z;
 }
 
-double Vector::Magnitude() {
+double Vector::Magnitude() const {
     return sqrt(x*x + y*y + z*z);
 }
 
-Vector Vector::ScaledBy(double v) {
+Vector Vector::ScaledBy(double v) const {
     Vector r;
 
     r.x = x * v;
@@ -632,7 +632,7 @@ Vector Vector::ScaledBy(double v) {
     return r;
 }
 
-Vector Vector::WithMagnitude(double v) {
+Vector Vector::WithMagnitude(double v) const {
     double m = Magnitude();
     if(EXACT(m == 0)) {
         // We can do a zero vector with zero magnitude, but not any other cases.
@@ -645,7 +645,7 @@ Vector Vector::WithMagnitude(double v) {
     }
 }
 
-Vector Vector::ProjectVectorInto(hEntity wrkpl) {
+Vector Vector::ProjectVectorInto(hEntity wrkpl) const {
     EntityBase *w = SK.GetEntity(wrkpl);
     Vector u = w->Normal()->NormalU();
     Vector v = w->Normal()->NormalV();
@@ -656,7 +656,7 @@ Vector Vector::ProjectVectorInto(hEntity wrkpl) {
     return (u.ScaledBy(up)).Plus(v.ScaledBy(vp));
 }
 
-Vector Vector::ProjectInto(hEntity wrkpl) {
+Vector Vector::ProjectInto(hEntity wrkpl) const {
     EntityBase *w = SK.GetEntity(wrkpl);
     Vector p0 = w->WorkplaneGetOffset();
 
@@ -665,25 +665,25 @@ Vector Vector::ProjectInto(hEntity wrkpl) {
     return p0.Plus(f.ProjectVectorInto(wrkpl));
 }
 
-Point2d Vector::Project2d(Vector u, Vector v) {
+Point2d Vector::Project2d(Vector u, Vector v) const {
     Point2d p;
     p.x = this->Dot(u);
     p.y = this->Dot(v);
     return p;
 }
 
-Point2d Vector::ProjectXy() {
+Point2d Vector::ProjectXy() const {
     Point2d p;
     p.x = x;
     p.y = y;
     return p;
 }
 
-Vector4 Vector::Project4d() {
+Vector4 Vector::Project4d() const {
     return Vector4::From(1, x, y, z);
 }
 
-double Vector::DivPivoting(Vector delta) {
+double Vector::DivPivoting(Vector delta) const {
     double mx = fabs(delta.x), my = fabs(delta.y), mz = fabs(delta.z);
 
     if(mx > my && mx > mz) {
@@ -695,7 +695,7 @@ double Vector::DivPivoting(Vector delta) {
     }
 }
 
-Vector Vector::ClosestOrtho() {
+Vector Vector::ClosestOrtho() const {
     double mx = fabs(x), my = fabs(y), mz = fabs(z);
 
     if(mx > my && mx > mz) {
@@ -707,7 +707,7 @@ Vector Vector::ClosestOrtho() {
     }
 }
 
-Vector Vector::ClampWithin(double minv, double maxv) {
+Vector Vector::ClampWithin(double minv, double maxv) const {
     Vector ret = *this;
 
     if(ret.x < minv) ret.x = minv;
@@ -721,7 +721,7 @@ Vector Vector::ClampWithin(double minv, double maxv) {
     return ret;
 }
 
-void Vector::MakeMaxMin(Vector *maxv, Vector *minv) {
+void Vector::MakeMaxMin(Vector *maxv, Vector *minv) const {
     maxv->x = max(maxv->x, x);
     maxv->y = max(maxv->y, y);
     maxv->z = max(maxv->z, z);
@@ -731,7 +731,7 @@ void Vector::MakeMaxMin(Vector *maxv, Vector *minv) {
     minv->z = min(minv->z, z);
 }
 
-bool Vector::OutsideAndNotOn(Vector maxv, Vector minv) {
+bool Vector::OutsideAndNotOn(Vector maxv, Vector minv) const {
     return (x > maxv.x + LENGTH_EPS) || (x < minv.x - LENGTH_EPS) ||
            (y > maxv.y + LENGTH_EPS) || (y < minv.y - LENGTH_EPS) ||
            (z > maxv.z + LENGTH_EPS) || (z < minv.z - LENGTH_EPS);
@@ -918,19 +918,19 @@ Vector4 Vector4::Blend(Vector4 a, Vector4 b, double t) {
     return (a.ScaledBy(1 - t)).Plus(b.ScaledBy(t));
 }
 
-Vector4 Vector4::Plus(Vector4 b) {
+Vector4 Vector4::Plus(Vector4 b) const {
     return Vector4::From(w + b.w, x + b.x, y + b.y, z + b.z);
 }
 
-Vector4 Vector4::Minus(Vector4 b) {
+Vector4 Vector4::Minus(Vector4 b) const {
     return Vector4::From(w - b.w, x - b.x, y - b.y, z - b.z);
 }
 
-Vector4 Vector4::ScaledBy(double s) {
+Vector4 Vector4::ScaledBy(double s) const {
     return Vector4::From(w*s, x*s, y*s, z*s);
 }
 
-Vector Vector4::PerspectiveProject() {
+Vector Vector4::PerspectiveProject() const {
     return Vector::From(x / w, y / w, z / w);
 }
 
@@ -1036,8 +1036,8 @@ BBox BBox::From(const Vector &p0, const Vector &p1) {
     return bbox;
 }
 
-Vector BBox::GetOrigin() { return minp.Plus(maxp.Minus(minp)).ScaledBy(0.5); }
-Vector BBox::GetExtents() { return maxp.Minus(minp).ScaledBy(0.5); }
+Vector BBox::GetOrigin() const { return minp.Plus(maxp.Minus(minp)).ScaledBy(0.5); }
+Vector BBox::GetExtents() const { return maxp.Minus(minp).ScaledBy(0.5); }
 
 void BBox::Include(const Vector &v, double r) {
     minp.x = min(minp.x, v.x - r);
@@ -1049,7 +1049,7 @@ void BBox::Include(const Vector &v, double r) {
     maxp.z = max(maxp.z, v.z + r);
 }
 
-bool BBox::Overlaps(BBox &b1) {
+bool BBox::Overlaps(const BBox &b1) const {
 
     Vector t = b1.GetOrigin().Minus(GetOrigin());
     Vector e = b1.GetExtents().Plus(GetExtents());
@@ -1057,6 +1057,6 @@ bool BBox::Overlaps(BBox &b1) {
     return fabs(t.x) < e.x && fabs(t.y) < e.y && fabs(t.z) < e.z;
 }
 
-bool BBox::Contains(const Point2d &p) {
+bool BBox::Contains(const Point2d &p) const {
     return p.x >= minp.x && p.y >= minp.y && p.x <= maxp.x && p.y <= maxp.y;
 }