Add a constraint to make line length equal arc length. That's quite

tricky; can't just use the dot product, since that blows up when
you cross pi radians. A gear shift approach, use either sin or cos,
same kind of thing as the 3d-parallel constraint.

And report a NaN constraint as unconverged, of course.

[git-p4: depot-paths = "//depot/solvespace/": change = 1890]

constraint.cpp

@@ -36,6 +36,7 @@ char *Constraint::DescriptionString(void) {
         case PERPENDICULAR:     s = "perpendicular"; break;
         case EQUAL_RADIUS:      s = "eq-radius"; break;
         case EQUAL_ANGLE:       s = "eq-angle"; break;
+        case EQUAL_LINE_ARC_LEN:s = "eq-line-len-arc-len"; break;
         case COMMENT:           s = "comment"; break;
         default:                s = "???"; break;
     }
@@ -239,6 +240,15 @@ void Constraint::MenuConstrain(int id) {
                 c.type = EQUAL_RADIUS;
                 c.entityA = gs.entity[0];
                 c.entityB = gs.entity[1];
+            } else if(gs.arcs == 1 && gs.lineSegments == 1 && gs.n == 2) {
+                c.type = EQUAL_LINE_ARC_LEN;
+                if(SS.GetEntity(gs.entity[0])->type == Entity::ARC_OF_CIRCLE) {
+                    c.entityA = gs.entity[1];
+                    c.entityB = gs.entity[0];
+                } else {
+                    c.entityA = gs.entity[0];
+                    c.entityB = gs.entity[1];
+                }
             } else {
                 Error("Bad selection for equal length / radius constraint. "
                       "This constraint can apply to:\r\n\r\n"
@@ -253,7 +263,9 @@ void Constraint::MenuConstrain(int id) {
                               "(equal angle between A,B and C,D)\r\n"
                       "    * three line segments or normals "
                               "(equal angle between A,B and B,C)\r\n"
-                      "    * two circles or arcs (equal radius)\r\n");
+                      "    * two circles or arcs (equal radius)\r\n"
+                      "    * a line segment and an arc "
+                              "(line segment length equals arc length)\r\n");
                 return;
             }
             if(c.type == EQUAL_ANGLE) {

constrainteq.cpp

@@ -249,6 +249,49 @@ void Constraint::GenerateReal(IdList<Equation,hEquation> *l) {
             break;
         }
 
+        case EQUAL_LINE_ARC_LEN: {
+            Entity *line = SS.GetEntity(entityA),
+                   *arc  = SS.GetEntity(entityB);
+
+            // Get the line length
+            ExprVector l0 = SS.GetEntity(line->point[0])->PointGetExprs(),
+                       l1 = SS.GetEntity(line->point[1])->PointGetExprs();
+            Expr *ll = (l1.Minus(l0)).Magnitude();
+
+            // And get the arc radius, and the cosine of its angle
+            Entity *ao = SS.GetEntity(arc->point[0]),
+                   *as = SS.GetEntity(arc->point[1]),
+                   *af = SS.GetEntity(arc->point[2]);
+
+            ExprVector aos = (as->PointGetExprs()).Minus(ao->PointGetExprs()),
+                       aof = (af->PointGetExprs()).Minus(ao->PointGetExprs());
+            Expr *r = aof.Magnitude();
+
+            ExprVector n = arc->Normal()->NormalExprsN();
+            ExprVector u = aos.WithMagnitude(Expr::From(1.0));
+            ExprVector v = n.Cross(u);
+            // so in our new csys, we start at (1, 0, 0)
+            Expr *costheta = aof.Dot(u)->Div(r);
+            Expr *sintheta = aof.Dot(v)->Div(r);
+
+            double thetas, thetaf, dtheta;
+            arc->ArcGetAngles(&thetas, &thetaf, &dtheta);
+            Expr *theta;
+            if(dtheta < 3*PI/4) {
+                theta = costheta->ACos();
+            } else if(dtheta < 5*PI/4) {
+                // As the angle crosses pi, cos theta is not invertible;
+                // so use the sine to stop blowing up
+                theta = Expr::From(PI)->Minus(sintheta->ASin());
+            } else {
+                theta = (Expr::From(2*PI))->Minus(costheta->ACos());
+            }
+
+            // And write the equation; r*theta = L
+            AddEq(l, (r->Times(theta))->Minus(ll), 0);
+            break;
+        }
+
         case POINTS_COINCIDENT: {
             Entity *a = SS.GetEntity(ptA);
             Entity *b = SS.GetEntity(ptB);

drawconstraint.cpp

@@ -104,6 +104,30 @@ void Constraint::DoEqualLenTicks(Vector a, Vector b, Vector gn) {
     LineDrawOrGetDistance(m.Minus(n), m.Plus(n));
 }
 
+void Constraint::DoEqualRadiusTicks(hEntity he) {
+    Entity *circ = SS.GetEntity(he);
+
+    Vector center = SS.GetEntity(circ->point[0])->PointGetNum();
+    double r = circ->CircleGetRadiusNum();
+    Quaternion q = circ->Normal()->NormalGetNum();
+    Vector u = q.RotationU(), v = q.RotationV();
+
+    double theta;
+    if(circ->type == Entity::CIRCLE) {
+        theta = PI/2;
+    } else if(circ->type == Entity::ARC_OF_CIRCLE) {
+        double thetaa, thetab, dtheta;
+        circ->ArcGetAngles(&thetaa, &thetab, &dtheta);
+        theta = thetaa + dtheta/2;
+    } else oops();
+
+    Vector d = u.ScaledBy(cos(theta)).Plus(v.ScaledBy(sin(theta)));
+    d = d.ScaledBy(r);
+    Vector p = center.Plus(d);
+    Vector tick = d.WithMagnitude(10/SS.GW.scale);
+    LineDrawOrGetDistance(p.Plus(tick), p.Minus(tick));
+}
+
 void Constraint::DoArcForAngle(Vector a0, Vector da, Vector b0, Vector db,
                                    Vector offset, Vector *ref)
 {
@@ -505,30 +529,22 @@ void Constraint::DrawOrGetDistance(Vector *labelPos) {
 
         case EQUAL_RADIUS: {
             for(int i = 0; i < 2; i++) {
-                Entity *circ = SS.GetEntity(i == 0 ? entityA : entityB);
+                DoEqualRadiusTicks(i == 0 ? entityA : entityB);
-                Vector center = SS.GetEntity(circ->point[0])->PointGetNum();
-                double r = circ->CircleGetRadiusNum();
-                Quaternion q = circ->Normal()->NormalGetNum();
-                Vector u = q.RotationU(), v = q.RotationV();
-
-                double theta;
-                if(circ->type == Entity::CIRCLE) {
-                    theta = PI/2;
-                } else if(circ->type == Entity::ARC_OF_CIRCLE) {
-                    double thetaa, thetab, dtheta;
-                    circ->ArcGetAngles(&thetaa, &thetab, &dtheta);
-                    theta = thetaa + dtheta/2;
-                } else oops();
-
-                Vector d = u.ScaledBy(cos(theta)).Plus(v.ScaledBy(sin(theta)));
-                d = d.ScaledBy(r);
-                Vector p = center.Plus(d);
-                Vector tick = d.WithMagnitude(10/SS.GW.scale);
-                LineDrawOrGetDistance(p.Plus(tick), p.Minus(tick));
             }
             break;
         }
 
+        case EQUAL_LINE_ARC_LEN: {
+            Entity *line = SS.GetEntity(entityA);
+            DoEqualLenTicks(
+                SS.GetEntity(line->point[0])->PointGetNum(),
+                SS.GetEntity(line->point[1])->PointGetNum(),
+                gn);
+
+            DoEqualRadiusTicks(entityB);
+            break;
+        }
+
         case LENGTH_RATIO:
         case EQUAL_LENGTH_LINES: {
             Vector a, b;

expr.cpp

@@ -243,6 +243,8 @@ int Expr::Children(void) {
         case SQUARE:
         case SIN:
         case COS:
+        case ASIN:
+        case ACOS:
             return 1;
 
         default: oops();
@@ -312,6 +314,8 @@ double Expr::Eval(void) {
         case SQUARE:        { double r = a->Eval(); return r*r; }
         case SIN:           return sin(a->Eval());
         case COS:           return cos(a->Eval());
+        case ACOS:          return acos(a->Eval());
+        case ASIN:          return asin(a->Eval());
 
         default: oops();
     }
@@ -349,6 +353,13 @@ Expr *Expr::PartialWrt(hParam p) {
         case SIN:       return (a->Cos())->Times(a->PartialWrt(p));
         case COS:       return ((a->Sin())->Times(a->PartialWrt(p)))->Negate();
 
+        case ASIN:
+            return (From(1)->Div((From(1)->Minus(a->Square()))->Sqrt()))
+                        ->Times(a->PartialWrt(p));
+        case ACOS:
+            return (From(-1)->Div((From(1)->Minus(a->Square()))->Sqrt()))
+                        ->Times(a->PartialWrt(p));
+
         default: oops();
     }
 }
@@ -421,6 +432,8 @@ Expr *Expr::FoldConstants(void) {
         case NEGATE:
         case SIN:
         case COS:
+        case ASIN:
+        case ACOS:
             if(n->a->op == CONSTANT) {
                 double nv = n->Eval();
                 n->op = CONSTANT;
@@ -526,6 +539,8 @@ p:
         case SQUARE:    App("(square "); a->PrintW(); App(")"); break;
         case SIN:       App("(sin "); a->PrintW(); App(")"); break;
         case COS:       App("(cos "); a->PrintW(); App(")"); break;
+        case ASIN:      App("(asin "); a->PrintW(); App(")"); break;
+        case ACOS:      App("(acos "); a->PrintW(); App(")"); break;
 
         default: oops();
     }

expr.h

@@ -29,6 +29,8 @@ public:
     static const int SQUARE         = 106;
     static const int SIN            = 107;
     static const int COS            = 108;
+    static const int ASIN           = 109;
+    static const int ACOS           = 110;
 
     // Special helpers for when we're parsing an expression from text.
     // Initially, literals (like a constant number) appear in the same
@@ -70,6 +72,8 @@ public:
     inline Expr *Square(void) { return AnyOp(SQUARE, NULL); }
     inline Expr *Sin   (void) { return AnyOp(SIN,    NULL); }
     inline Expr *Cos   (void) { return AnyOp(COS,    NULL); }
+    inline Expr *ASin  (void) { return AnyOp(ASIN,   NULL); }
+    inline Expr *ACos  (void) { return AnyOp(ACOS,   NULL); }
 
     Expr *PartialWrt(hParam p);
     double Eval(void);

sketch.h

@@ -458,6 +458,7 @@ public:
     static const int EQ_LEN_PT_LINE_D   =  52;
     static const int EQ_PT_LN_DISTANCES =  53;
     static const int EQUAL_ANGLE        =  54;
+    static const int EQUAL_LINE_ARC_LEN =  55;
     static const int SYMMETRIC          =  60;
     static const int SYMMETRIC_HORIZ    =  61;
     static const int SYMMETRIC_VERT     =  62;
@@ -527,6 +528,7 @@ public:
     void DoLabel(Vector ref, Vector *labelPos, Vector gr, Vector gu);
     void DoProjectedPoint(Vector *p);
     void DoEqualLenTicks(Vector a, Vector b, Vector gn);
+    void DoEqualRadiusTicks(hEntity he);
 
     double GetDistance(Point2d mp);
     Vector GetLabelPos(void);

system.cpp

@@ -516,7 +516,7 @@ didnt_converge:
     (g->solved.remove).Clear();
     SS.constraint.ClearTags();
     for(i = 0; i < eq.n; i++) {
-        if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
+        if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE || isnan(mat.B.num[i])) {
             // This constraint is unsatisfied.
             if(!mat.eq[i].isFromConstraint()) continue;
 

textwin.cpp

@@ -343,6 +343,9 @@ void TextWindow::DescribeSelection(void) {
                 double r = e->CircleGetRadiusNum();
                 Printf(true, "   diameter =  %Fi%s", SS.MmToString(r*2));
                 Printf(false, "     radius =  %Fi%s", SS.MmToString(r));
+                double thetas, thetaf, dtheta;
+                e->ArcGetAngles(&thetas, &thetaf, &dtheta);
+                Printf(false, "    arc len =  %Fi%s", SS.MmToString(dtheta*r));
                 break;
             }
             case Entity::CIRCLE: {