home *** CD-ROM | disk | FTP | other *** search
/ Dream 52 / Amiga_Dream_52.iso / Linux / Divers / freedraft.tar.gz / freedraft.tar / FREEdraft-050298 / GEOMLIB2D / geomlib2d.cpp < prev    next >
C/C++ Source or Header  |  1998-04-27  |  41KB  |  1,455 lines
  1.  
  2. // Copyright (C) 1997  Cliff Johnson                                       //
  3. // Copyright (C) 1998  Cliff Johnson                                       //
  4. //                                                                         //
  5. // This program is free software; you can redistribute it and/or           //
  6. // modify it under the terms of the GNU  General Public                    //
  7. // License as published by the Free Software Foundation; either            //
  8. // version 2 of the License, or (at your option) any later version.        //
  9. //                                                                         //
  10. // This software is distributed in the hope that it will be useful,        //
  11. // but WITHOUT ANY WARRANTY; without even the implied warranty of          //
  12. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU       //
  13. // General Public License for more details.                                //
  14. //                                                                         //
  15. // You should have received a copy of the GNU General Public License       //
  16. // along with this software (see COPYING.LIB); if not, write to the        //
  17. // Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. //
  18.  
  19.  
  20. //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21.  
  22. //
  23. //    FreeEngineer 2D Drafting Library Functions - libFE2D
  24. //
  25. //
  26. //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  27.  
  28. #include <math.h>
  29. #include <vector>
  30.  
  31. #include "geomlib2d.h"
  32. #include "point.h"
  33. #include "line.h"
  34. #include "segment.h"
  35. #include "circle.h"
  36. #include "arc.h"
  37. #include "geomexception.h"
  38. #include "pick.h"
  39. #include "pickgeom.h"
  40. #include "geom_enum.h"
  41. //#include <menuhandler_enum.h>
  42.  
  43. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  44. //    GL2D_PointByConstraints(Pick,Pick)
  45. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  46.  
  47. Point GL2D_PointByConstraints(const Pick& pk1, const Pick& pk2) 
  48.     throw (GeomException)
  49. {
  50.  
  51. // point point
  52.     if( pk1.Type()==POINT && pk2.Type()==POINT)
  53.     {
  54.         
  55.         Point p1 = GetPointFromPick(pk1);
  56.         Point p2 = GetPointFromPick(pk2);
  57.         return Point(p1 + 0.5 * (p2 - p1));
  58.     }
  59. // point vs. line/segment 
  60.     if( ( pk1.Type()==POINT && ( pk2.Type()==LINE || pk2.Type()==SEGMENT     ) ) 
  61.     ||  ( pk2.Type()==POINT && ( pk1.Type()==LINE || pk1.Type()==SEGMENT     ) ) ) 
  62.     {
  63. // swap around if necessary to get point first
  64.         Pick lpk1 = pk1;
  65.         Pick lpk2 = pk2;
  66.         if( pk2.Type()==POINT)
  67.         {
  68.             Pick lpktmp = lpk2;
  69.             lpk2 = lpk1;
  70.             lpk1 = lpktmp;
  71.         }
  72. // get the point and line
  73.         Point p = GetPointFromPick(lpk1);
  74.         Line l = GetLineFromPick(lpk2);
  75.  
  76. // project point onto line
  77.         return l.Project(p);
  78.     
  79.     }
  80. // point vs. circle/arc 
  81.     if( ( pk1.Type()==POINT && ( pk2.Type()==CIRCLE || pk2.Type()==ARC     ) ) 
  82.     ||  ( pk2.Type()==POINT && ( pk1.Type()==CIRCLE || pk1.Type()==ARC     ) ) ) 
  83.     {
  84. // swap around if necessary to get point first
  85.         Pick lpk1 = pk1;
  86.         Pick lpk2 = pk2;
  87.         if( pk2.Type()==POINT)
  88.         {
  89.             Pick lpktmp = lpk2;
  90.             lpk2 = lpk1;
  91.             lpk1 = lpktmp;
  92.         }
  93. // get the point and circle 
  94.         Point p = GetPointFromPick(lpk1);
  95.         Circle c = GetCircleFromPick(lpk2);
  96.  
  97. // project point onto circle
  98.         Point px1 = c.Project(p);
  99. // find opposite solution for circle
  100.         Point px2 = ( 2 * c.Center()) - px1 ;
  101.  
  102. // return solution closest to pick 2 point
  103.         if (Distance(px1,lpk2.GetPoint()) < Distance(px2,lpk2.GetPoint()))
  104.             return px1;
  105.         return px2;
  106.     }
  107.  
  108. // line/segment vs. line/segment 
  109.     if(    ( pk1.Type()==LINE || pk1.Type()==SEGMENT     ) &&
  110.              ( pk2.Type()==LINE || pk2.Type()==SEGMENT     ))
  111.     {
  112.         Line l1 = GetLineFromPick(pk1);
  113.         Line l2 = GetLineFromPick(pk2);
  114.  
  115. // solve the system
  116.         return  GL2D_SolveLineLineIntersection(l1,l2);
  117.     }
  118.  
  119.  
  120.  
  121. // line/segment vs. circle/arc
  122.  
  123.     if((     ( pk1.Type()==LINE || pk1.Type()==SEGMENT     ) &&
  124.                  ( pk2.Type()==CIRCLE || pk2.Type()==ARC    )    ) ||
  125.                (     ( pk1.Type()==CIRCLE || pk1.Type()==ARC    ) &&
  126.                  ( pk2.Type()==LINE || pk2.Type()==SEGMENT     )     )    ) 
  127.     {   
  128.  
  129. // swap around if necessary to get line/segment first
  130.         Pick lpk1 = pk1;
  131.         Pick lpk2 = pk2;
  132.         if( pk2.Type()==LINE || pk2.Type()==SEGMENT)
  133.         {
  134.             Pick lpktmp = lpk2;
  135.             lpk2 = lpk1;
  136.             lpk1 = lpktmp;
  137.         }
  138. // get the line
  139.         Line l = GetLineFromPick(lpk1);
  140.  
  141. // get the circle;
  142.         Circle c = GetCircleFromPick(lpk2);
  143.  
  144. // get the solutions
  145.         return GL2D_SolveLineCircleIntersection(l,c,pk2.GetPoint());
  146.     }
  147. // circle/arc vs. circle/arc
  148.     if(( pk1.Type()==CIRCLE || pk1.Type()==ARC     ) &&
  149.              ( pk2.Type()==CIRCLE || pk2.Type()==ARC     ))
  150.     { 
  151.  
  152. // in case of arcs, get circles
  153.  
  154.         Circle c1 = GetCircleFromPick(pk1);
  155.         Circle c2 = GetCircleFromPick(pk2);
  156. // solve
  157.         return GL2D_SolveCircleCircleIntersection(c1,c2,pk2.GetPoint());
  158.     }  
  159.     
  160.     throw GeomException("GL2D_PointByConstraints() : a Pick is not any expected type");
  161. }
  162.     
  163. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  164.  
  165. Line GL2D_LineParallel(const Line& l, const Pick& pk)
  166.     throw (GeomException)
  167. {
  168.     
  169.     switch( pk.Type())
  170.     {
  171.     // Line/Segment, Point        - Parallel to Line/Segment, passing throgh point
  172.     case POINT:
  173.         {
  174.         Point p = GetPointFromPick(pk);
  175.         return Line(p,l.Direction());
  176.         break;
  177.         }
  178.      // Circle/arc    - tangent to Circle/arc
  179.     case CIRCLE:
  180.     case ARC:
  181.         {
  182.         Circle c = GetCircleFromPick(pk);
  183.         
  184.         Point R = Cross(l.Direction(),Point(0,0,1)) * c.Radius();
  185.         Point p1 = c.Center() + R;
  186.         Point p2 = c.Center() - R;
  187.  
  188.     //choose the solution closest to the selected point.
  189.  
  190.         if(Distance( pk.GetPoint(), p1 ) < Distance( pk.GetPoint(), p2))
  191.         {
  192.             return Line(p1,l.Direction());
  193.         }
  194.         else
  195.         {
  196.             return Line(p2,l.Direction());
  197.         }
  198.  
  199.         break;
  200.         }
  201.     }
  202.     throw GeomException("GL2D_LineParallel() : a Pick is not any expected type");
  203. }
  204. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  205. //
  206. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  207.  
  208. Line GL2D_LineByConstraints(const Pick& pk1, const Pick& pk2)
  209.     throw (GeomException)
  210. {
  211.  
  212.     // determine which case to solve (lets get this out of the way up front)
  213.  
  214.     // POINT-POINT
  215.     if(pk1.Type() == POINT && pk2.Type() == POINT)
  216.     {
  217.         Point p1  = GetPointFromPick(pk1);
  218.         Point p2  = GetPointFromPick(pk2);
  219.         if (p1 == p2 )
  220.             throw GeomException("GeomLib2d::GL2D_LineByConstraints() : Points are equal");
  221.         return Line(p1, p2 - p1);
  222.     }
  223.     
  224.     // POINT-LINE
  225.     else if(( pk1.Type() == POINT && (pk2.Type() == LINE || pk2.Type() == SEGMENT )) ||
  226.     ((pk1.Type() == LINE || pk1.Type() == SEGMENT) && pk2.Type() == POINT ))
  227.     {
  228.  
  229. // local copies
  230.         Pick lpk1 = pk1;
  231.         Pick lpk2 = pk2;
  232.  
  233. // reorder picks - point first
  234.         if(lpk2.Type()==POINT)
  235.         {
  236.             Pick tmp = lpk1;
  237.             lpk1 = lpk2;
  238.             lpk2 = tmp;
  239.         }
  240.  
  241.     // get point 
  242.         Point p = GetPointFromPick(lpk1);
  243.     // get line
  244.         Line l = GetLineFromPick(lpk2);
  245.     // parallel to line, passing through point
  246.         return Line(p,l.Direction());
  247.     }
  248.         
  249.     // POINT-CIRCLE
  250.     else if(( pk1.Type() == POINT && (pk2.Type() == CIRCLE || pk2.Type() == ARC )) ||
  251.     ((pk1.Type() == CIRCLE || pk1.Type() == ARC) && pk2.Type() == POINT ))
  252.     {
  253. // local copies
  254.         Pick lpk1 = pk1;
  255.         Pick lpk2 = pk2;
  256.  
  257. // reorder picks - point first
  258.         if(lpk2.Type()==POINT)
  259.         {
  260.             Pick tmp = lpk1;
  261.             lpk1 = lpk2;
  262.             lpk2 = tmp;
  263.         }
  264.  
  265.     // get geometry
  266.         Point p = GetPointFromPick(lpk1);
  267.         Circle c = GetCircleFromPick(lpk2);
  268.  
  269.         return GL2D_SolveLineTangentPointCircle(p,c,lpk2.GetPoint());
  270.     }
  271.     
  272.     // LINE-LINE
  273.     else if((pk1.Type() == LINE || pk1.Type() == SEGMENT ) && 
  274.                 (pk2.Type() == LINE || pk2.Type() == SEGMENT )  )
  275.     {
  276.         Line l1 = GetLineFromPick(pk1);
  277.         Line l2 = GetLineFromPick(pk2);
  278.  
  279.         return GL2D_SolveLineConstrainedLineLine(l1,l2,pk1.GetPoint(), pk2.GetPoint());
  280.     }
  281.  
  282.     // LINE-CIRCLE
  283.     else if(( ( pk1.Type() == LINE   || pk1.Type() == SEGMENT) && 
  284.           ( pk2.Type() == CIRCLE || pk2.Type() == ARC    )     ) ||
  285.             ( (pk1.Type() == CIRCLE  || pk1.Type() == ARC    ) && 
  286.               (pk2.Type() == LINE  || pk2.Type() == SEGMENT  )     ))
  287.     {
  288. // local copies
  289.         Pick lpk1 = pk1;
  290.         Pick lpk2 = pk2;
  291.  
  292. // reorder picks - line /segment first
  293.         if(lpk2.Type()==LINE || lpk2.Type() == SEGMENT)
  294.         {
  295.             Pick tmp = lpk1;
  296.             lpk1 = lpk2;
  297.             lpk2 = tmp;
  298.         }
  299.  
  300. // get line 
  301.         Line l = GetLineFromPick(lpk1);
  302.         Circle c = GetCircleFromPick(lpk2);
  303.  
  304.         return GL2D_LineParallel(l,lpk2);
  305.     }
  306.  
  307.     // CIRCLE-CIRCLE
  308.     else if((pk1.Type() == CIRCLE || pk1.Type() == ARC ) && 
  309.             (pk2.Type() == CIRCLE || pk2.Type() == ARC )  )
  310.     {
  311.         Circle c1 = GetCircleFromPick(pk1);
  312.         Circle c2 = GetCircleFromPick(pk2);
  313.  
  314.         return GL2D_SolveLineTangentCircleCircle(c1,c2,pk1.GetPoint(),pk2.GetPoint());
  315.     }
  316.  
  317.     throw GeomException("GeomLib2d::GL2D_LineByConstraints() : syntax error");
  318. }    
  319.  
  320.  
  321. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  322. //
  323. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  324. // 
  325. // Line Orthoginal(const Pick& pk1 , const Pick& pk2 ) // Constraint, Subject
  326. // {
  327. //     // Line/Segment, Point        - Ortho to Line/Segment, passing throgh point
  328. //     // Circle/arc, Point        - Ortho to circle/arc, passing through point
  329. //     // Line/Segment, Circle/arc    - Ortho to Line/Segment, tangent to Circle/arc
  330. //     // Circle/arc, Circle/arc     - Ortho to circle/arc, tangent to circle/arc
  331. // }
  332. // 
  333. // Circles:
  334. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  335. //
  336. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  337. // 
  338. // Circle GL2D_CircleCenter(const Point& p, double r){ // solve }
  339. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  340. //
  341. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  342.  
  343. Circle GL2D_CircleCenter(const Pick& pk1, const Pick& pk2) 
  344.     throw (GeomException)
  345. {
  346. // pk1 must be a  Point
  347. // pk2 can be  Point,Line or Circle, segment or arc
  348.  
  349.  
  350.     if(pk1.Type() != POINT)
  351.         throw GeomException("GL2D_CircleCenter() : Pick 1 is not a point");
  352.  
  353.     Point ctr = GetPointFromPick(pk1);
  354.     Point origin1,origin2;
  355.  
  356.     switch(pk2.Type())
  357.     {
  358.     // Passing another point
  359.     case POINT:
  360.         {
  361.         origin1 = GetPointFromPick(pk2);
  362.         if(ctr == origin1)
  363.             throw GeomException("GL2D_CircleCenter() : center and origin points are equal");
  364.         return Circle(ctr,origin1);
  365.         }
  366.     // Tangent to a line
  367.     case LINE:
  368.     case SEGMENT:
  369.         {
  370.         Line l = GetLineFromPick(pk2);
  371.         origin1 = l.Project(ctr);
  372.         if(ctr == origin1)
  373.             throw GeomException("GL2D_CircleCenter() : Line passes through Circle center");
  374.         return Circle(ctr,origin1);
  375.         }
  376.     // Tangent to a Circle, arc
  377.     case CIRCLE:
  378.     case ARC:
  379.         {
  380.         Circle c = GetCircleFromPick(pk2);
  381.     // selected center is center of selected circle
  382.         if(c.Center() == ctr)
  383.             throw GeomException("GL2D_CircleCenter() : selected Point is center of selected Circle");
  384.     // point is tangent to circle
  385.         if(fabs(c.Radius() - Distance(ctr,c.Center())) < Point::NullDist)
  386.             throw GeomException("GL2D_CircleCenter() : selected Point is tangent to selected Circle");
  387. // find two solutions
  388.         origin1 = c.Project(ctr); // projection
  389.         origin2 =  ( 2 * c.Center())  - origin1  ;
  390.         if( origin1.Distance(pk2.GetPoint()) < origin2.Distance(pk2.GetPoint()) )
  391. // return solution closest to pick point
  392.             return Circle(ctr,origin1);
  393.         else
  394.             return Circle(ctr,origin2);
  395.         }
  396.     default:
  397.         throw GeomException("GL2D_CircleCenter() : second Pick is not one of the expected types");
  398.     }
  399. }
  400. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  401. //
  402. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  403. // 
  404. // Circle CircleParallel(const Pick& pk, double d)
  405. // {
  406. //     //Circle/Arc
  407. // }
  408. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  409. //
  410. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  411. // 
  412. // Circle CircleParallel(const Pick& c, const Pick& t) // parallel to c,  tangent to t
  413. // {
  414. //     // t is Point
  415. //     // t is Line/Segment
  416. //     // t is Circle/Arc
  417. // }
  418. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  419. //
  420. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  421. // 
  422. // Circle CircleConstraint(const Pick&, const Pick&, const Pick&)
  423. // {
  424. //     // Point-Point-Point
  425. //     // Point-Point-Line
  426. //     // Point-Point-Circle
  427. // 
  428. //     // Point-Line-Line        PLL
  429. //     // Point-Line-Circle        PLC
  430. //     // Point-Circle-Circle        PCC
  431. //     // Line-Line-Line        LLL
  432. //     // Line-Line-Circle        LLC
  433. //     // Line-Circle-Circle        LCC
  434. //     // Circle-Circle-Circle        CCC
  435. // }
  436. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  437. //
  438. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  439.  
  440. // private "subroutines"
  441.  
  442. Point GL2D_SolveLineLineIntersection(const Line& l1, const Line& l2)
  443.     throw (GeomException)
  444. {
  445.     double a1,b1,c1;
  446.     double a2,b2,c2;
  447.  
  448. //cerr << endl << " GL2D_SolveLineLineIntersection( " << endl;
  449. //cerr << l1 << " , " << l2 << " ) " << endl;
  450. // get directions
  451.  
  452.     Vector v1 = l1.Direction();
  453.     Vector v2 = l2.Direction();
  454.  
  455.  
  456. //cerr << "v1 = " << v1 << " v2 = " << v2 << endl;
  457. //+ + + + + + + + + + + + + + 
  458. // check for degenerate cases
  459. //+ + + + + + + + + + + + + + 
  460.  
  461. // parallel lines
  462.  
  463.     a1 = Angle(v1,v2);
  464. //cerr << " a1 = " << a1 << endl;
  465.     a2 = Angle(v1,-v2);
  466. //cerr << " a2 = " << a2 << endl;
  467.     if(a1 < Point::NullAngle || a2 < Point::NullAngle)
  468.     {
  469. //cerr << "throwing exception" << endl;
  470.         throw GeomException("GL2D_SolveLineLineIntersection() : parallel lines");
  471.     }
  472.  
  473. //cerr << " lines are NOT parallel " << endl;
  474. //+ + ++ + + + + + + + + + + 
  475. // Formulate explicit standard form  ax + by + c = 0 
  476. // reference: faux and pratt, chapter 1
  477. //+ + + + + + + + + + + + + + 
  478.  
  479. // account for vertical lines as well
  480.  
  481.     int lineflag1=0,lineflag2=0;
  482.  
  483. // line 1
  484.     lineflag1 = GL2D_NormalizedExplicitForm(a1,b1,c1,l1);
  485.  
  486. // line 2
  487.     lineflag2 = GL2D_NormalizedExplicitForm(a2,b2,c2,l2);
  488.  
  489. // vertical line 1
  490.     
  491.     Point p;
  492.  
  493.     // check the results of the last operations 
  494.  
  495.     if(lineflag1 == LINE_IS_VERTICAL)
  496.  
  497.     // ax + by + c = 0:   y = -(ax + c ) /b
  498.     // line 1 vertical:  y value on line 2 =>   c = -x  --> y = -(ax -x)/b =  (-ax + x)/b
  499.     //                              = x(1-a)/b 
  500.  
  501.     {
  502.     // solve line 2 at  x = -c1
  503.         // by = -a*c1 + c2
  504.  
  505.         double x = -c1/a1;
  506.  
  507.         p = Point(  x  , -(a2*x+c2)/b2 );
  508.     }
  509.  
  510. // vertical line 2
  511.  
  512.     else if(lineflag2 == LINE_IS_VERTICAL)
  513.     {
  514.         // solve line 1 at Y = -c2; 
  515.         double x = -c2/a2;
  516.         p = Point(  x  , -(a1*x+c1)/b1 );
  517.     }
  518.  
  519. // general case 
  520.  
  521.  
  522.     else
  523.     {
  524.     //f & p eq. 1.8
  525.         p = Point((b1*c2 - b2*c1)/(a1*b2 - a2*b1),(c1*a2 - c2*a1)/(a1*b2 - a2*b1));
  526.     }
  527.  
  528.     return Point(p);
  529. }
  530.  
  531. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  532. //
  533. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  534.  
  535.  
  536. int GL2D_NormalizedExplicitForm(double& a, double& b, double& c,const Line& l)
  537.     throw (GeomException)
  538. {
  539.  
  540.     double x1,y1,x2,y2;
  541.  
  542.     Point p1 = l.Origin();
  543.     Point p2 = p1 + l.Direction();
  544.     x1 = p1[0];
  545.     y1 = p1[1];
  546.     x2 = p2[0];
  547.     y2 = p2[1];
  548.     
  549.  
  550.     if ( fabs(x2 - x1)  < Point::NullDist )
  551.     {
  552.         a=1;
  553.         b=0;
  554.         c=-x1;
  555.         if(c > 0){ a = -a; c = -c; }
  556.         return LINE_IS_VERTICAL;
  557.     }
  558.     
  559.     a = (y2 - y1) / (x2 - x1 ); 
  560.     b = -1;
  561.     c = ( x2 * y1 - x1 * y2 ) / ( x2 - x1) ;   
  562.  
  563.  
  564.     // normalize into standard form
  565.  
  566.  
  567.     if( fabs(a) < Point::NullDist&& fabs(b) < Point::NullDist)
  568.         throw GeomException("GL2D_NormalizedExplicitForm() : zero denomenator");
  569.  
  570.     if(c > 0){ a = -a; b = -b; c = -c; }     // 1st rule: c < 0, multiply by -1 
  571.  
  572.     // now scale so that a**2 + b**2 = 1
  573.  
  574.     double s = sqrt(1/ (a * a + b * b));
  575.     a = a * s;
  576.     b = b * s;
  577.     c = c * s;
  578.  
  579.     return 0;
  580. }
  581.  
  582. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  583. //
  584. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  585.  
  586. // rotate a vector in the plane
  587.         
  588. Vector GL2D_Rotate(const Vector& v ,double theta) throw(GeomException)
  589. {
  590.     // rotate v by angle theta
  591.  
  592.     return Vector(v[0]*cos(theta)-v[1]*sin(theta), v[0]*sin(theta)+v[1]*cos(theta));
  593. }
  594. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  595.  
  596. Line GL2D_SolveLineTangentPointCircle(const Point& p, const Circle& c, const Point& pkpt)
  597.     throw (GeomException)
  598. {
  599.  
  600. // check for degenerate cases
  601.         
  602.     Point ctr = c.Center();
  603.     double r = c.Radius();
  604.     double d = Distance(ctr, p);
  605.     // point inside circle
  606.     if(d < r)throw GeomException("GL2D_SolveLineTangentPointCircle : point is inside circle");
  607.     
  608. // general solution
  609.     double theta = acos( r / Distance(c.Center(),p) );
  610.     Vector v = p - c.Center();
  611.     v = v.Normal();            // should be safe from exception :)
  612.     v = v * r;
  613.  
  614.     Point s1 = c.Center() + GL2D_Rotate(v,theta);
  615.     Point s2 = c.Center() + GL2D_Rotate(v,-theta);
  616.  
  617. // return solution closest to pick point.
  618.  
  619.     if(Distance(s1,pkpt) < Distance(s2,pkpt) ) return Line(s1, p - s1);
  620.     return Line(s2, p - s2);
  621. }
  622.  
  623. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  624. Line GL2D_SolveLineConstrainedLineLine(const Line&  l1,const Line& l2, const Point& p1, const Point& p2)
  625.     throw (GeomException)
  626. {
  627.  
  628. //cerr << " GL2D_SolveLineConstrainedLineLine( " << endl;
  629. //cerr << l1 << " , " << endl;
  630. //cerr << l2 << " , " << endl;
  631. //cerr << p1 << " , " << endl;
  632. //cerr << p2 << " ) " << endl;
  633.  
  634. // general solution
  635. // intersect the lines, get the point
  636.  
  637.     Point p;
  638.     try 
  639.     {
  640. //cerr << "trying" << endl;
  641.         p = GL2D_SolveLineLineIntersection(l1,l2);
  642.     }
  643.     catch(GeomException & ge) 
  644.     {
  645. //cerr << "catching" << endl;
  646.         return Line(l1.Origin() + 0.5 * (l2.Origin() - l1.Origin()), l1.Direction());
  647.     }
  648.  
  649. // otherwise
  650.     // find both solutions
  651.  
  652.     Line lx1(p,l1.Direction() + l2.Direction());    // remember, line directions are 
  653.       Line lx2(p,l1.Direction() - l2.Direction());     // unit vectors so I can do this!  
  654.  
  655. // solution is closest to midpoint between selection points
  656.     p = p1 + 0.5 * ( p2 - p1 );
  657.     if(Distance(p,lx1.Project(p)) < Distance(p,lx2.Project(p)))return lx1;
  658.     else return lx2;
  659. }
  660. //+++++++++++++++++++++++++++++++++++++++++++++++++++++
  661. Line GL2D_SolveLineTangentCircleCircle(const Circle&  c1, const Circle& c2, const Point&   p1, const Point&  p2)
  662.     throw (GeomException)
  663. {
  664.  
  665.     // Degenerate cases
  666.     // circles have same center - no solution
  667.     double dc = Distance(c1.Center(), c2.Center());    // do these only once cause they
  668.     double r1 = c1.Radius();            // use square roots
  669.     double r2 = c2.Radius();
  670.  
  671. // make local copies    ( must make locals so I can swap them if necessary)
  672.     Circle cx1(c1);
  673.     Circle cx2(c2);
  674.     Point px1(p1);
  675.     Point px2(p2);
  676.  
  677. // assure larger circle is cx1    
  678.     if(r1 < r2)
  679.     {
  680.         cx1 = c2;
  681.         cx2 = c1;
  682.         double tmp = r2;
  683.         r2 = r1;
  684.         r1 = tmp;
  685.         Point ptmp = px2;
  686.         px2 = px1;
  687.         px1 = ptmp;
  688.     }
  689.  
  690.     if( dc < Point::NullDist )
  691.         throw GeomException("GL2D_SolveLineTangentCircleCircle() : coincident Circle centers");
  692.  
  693.     // one circle completly inside the other - no solution
  694.  
  695.     else if( (dc+r2-r1) < -Point::NullDist)
  696.         throw GeomException("GL2D_SolveLineTangentCircleCircle() : Circle inside Circle  ");
  697.  
  698.     // one circle inside another and tangent to it - one solution
  699.  
  700.     else if( fabs(dc+r2-r1) < Point::NullDist)
  701.     {
  702.         Vector v = cx2.Center() - cx1.Center();
  703.         v = v.Normal() * cx2.Radius();          //should be exception safe 
  704.         return Line( cx2.Center() + v, Cross(v,Point(0,0,1)));
  705.     }
  706.  
  707.  
  708. // first find classical - outside solutions
  709.  
  710.     Vector v = cx2.Center() - cx1.Center();
  711.  
  712.     Vector v1= v.Normal() * cx1.Radius();    // v normalized to larger radius
  713.     Vector v2= v.Normal() * cx2.Radius();    // v normalized to smaller radius
  714.  
  715.     double theta = acos((r1-r2)/dc);    // calculate critical angle
  716.  
  717.     Point p11 = cx1.Center() + GL2D_Rotate(v1,theta);    // p11 - circle 1 solution 2
  718.     Point p12 = cx1.Center() + GL2D_Rotate(v1,-theta);   // p12 - circle 1 solution 2
  719.  
  720.     Point p21 = cx2.Center() + GL2D_Rotate(v2,theta);    // p21 - circle 2 solution 1
  721.     Point p22 = cx2.Center() + GL2D_Rotate(v2,-theta);   // p22 - circle 2 solution 2
  722.  
  723. // if only two solutions, solve now
  724.     if((dc - r2) < (r1 - Point::NullDist))
  725.     {
  726.         if(Distance(p1,p11) < Distance(p1,p12))
  727.         {
  728.             return Line(p11,p21-p11);
  729.         }
  730.         return Line(p12,p22-p12);
  731.     }
  732.  
  733. // now find possible third solution ---- circles are tangent ----
  734.     if( fabs(dc - r2 - r1)  < Point::NullDist) 
  735.     {
  736.         Point p13 = cx1.Center() + v1;    
  737.         Point p23 = Cross(v1,Point(0,0,1));
  738.  
  739.         // choose closest solution based on c1 picks
  740.         double d1 = Distance(p11,p1);
  741.         double d2 = Distance(p12,p1);
  742.         double d3 = Distance(p13,p1);
  743.         if(d1 < d2 && d1 < d3)return Line(p11,p21-p11);
  744.         else if(d2 < d3)return Line(p12,p22-p12);
  745.         return Line(p13,p23-p13);
  746.     }
  747.  
  748. // distinct non-intersecting, non-tangent, circles - four solutions
  749.     theta = acos((r1+r2)/dc);    // calculate secondary critical angle
  750.  
  751.     Point p13 = cx1.Center() + GL2D_Rotate(v1,theta);    // circle 1 solution 3
  752.     Point p14 = cx1.Center() + GL2D_Rotate(v1,-theta);   // circle 1 solution 4
  753.  
  754.     Point p23 = cx2.Center() + GL2D_Rotate(v2,-(M_PI-theta)); // circle 2 solution 3
  755.     Point p24 = cx2.Center() + GL2D_Rotate(v2,M_PI-theta);    // circle 2 solution 4
  756.  
  757. // now return correct solution ... hmmmm
  758.  
  759.     double d[4];
  760.     d[0] = Distance(p11,px1) + Distance(p21,px2);
  761.     d[1] = Distance(p12,px1) + Distance(p22,px2);
  762.     d[2] = Distance(p13,px1) + Distance(p23,px2);
  763.     d[3] = Distance(p14,px1) + Distance(p24,px2);
  764.  
  765.     double dmin = d[0];
  766.     int idx = 0;
  767.     for(int i=1;i<4;i++)
  768.     {
  769.         if(d[i] < dmin)
  770.         {
  771.             dmin = d[i];
  772.             idx = i;
  773.         }
  774.     }
  775.  
  776.     switch (idx)
  777.     {
  778.         case 0:
  779.         // outside solution 1
  780.         {
  781.             return Line(p11,p21-p11);
  782.         }
  783.         case 1:
  784.         // outside solution 2
  785.         {
  786.             return Line(p12,p22-p12);
  787.         }
  788.         case 2:
  789.         // inside solution 3
  790.         {
  791.             return Line(p13,p23-p13);
  792.         }
  793.     }
  794.     // inside solution 4
  795.     return Line(p14,p24-p14);        
  796. }
  797. //***********************************************************************************
  798.  
  799. Arc GL2D_SolveArc3Points(const Point& p1, const Point& p2, const Point& p3) throw (GeomException)
  800. {
  801.     Circle c = GL2D_SolveCircle3Points(p1,p2,p3);
  802.     return Arc(c.Center(),p1,p3);
  803. }
  804.  
  805. //***********************************************************************************
  806. Circle GL2D_SolveCircle3Points(const Point& p1, const Point& p2, const Point& p3)
  807.     throw (GeomException)
  808. {
  809.  
  810.     Point px1(p1);
  811.     Point px2(p2);
  812.     Point px3(p3);
  813.  
  814. // error conditions
  815. // points are equal
  816.     if(px1 == px2 || px1 == px3 || px2 == px3)
  817.         throw GeomException("GL2D_SolveArc3Points() : some Points are equal");
  818. // points are co-linear
  819.     
  820.     Vector v1 = px2 - px1;
  821.         Vector v2 = px3 - px1;
  822.     if(Angle(v1,v2) < Point::NullAngle)
  823.         throw GeomException("GL2D_SolveArc3Points() : Points are colinear");
  824.  
  825. // check for the sense of rotation  - reverse the points if necessary
  826.         Vector n = Cross(v1,v2);
  827.         if(n[2] < 0)
  828.         {
  829.         Point tmpx = px3;
  830.         px3 = px1;
  831.         px1 = tmpx;
  832.         v1 = px2 - px1;
  833.         v2 = px3 - px1;
  834.         }
  835. // find the cener from the 3 points 
  836.     Point m1 = px1 + 0.5 * v1; 
  837.     Point m2 = px1 + 0.5 * v2;
  838.  
  839.     Vector d1 = Cross(v1,Vector(0,0,1));
  840.     Vector d2 = Cross(v2,Vector(0,0,1));
  841.  
  842.     Line l1(m1,d1);
  843.     Line l2(m2,d2);
  844.  
  845.     Point c = GL2D_SolveLineLineIntersection(l1,l2);
  846.  
  847.     return Circle(c,px1);
  848. }
  849. //***********************************************************************************
  850.  
  851. Point GL2D_SolveLineCircleIntersection(const Line& l, const Circle& c, const Point& pickloc) 
  852.     throw (GeomException)
  853. {
  854.  
  855. // copy the line and circle
  856.     Line lx(l);
  857.     Circle cx(c);
  858.  
  859. // get distance from center to line
  860.     Point center = cx.Center();
  861.     double d = lx.Distance(cx.Center());
  862.  
  863. // line tangent to circle
  864.     if(fabs(d - cx.Radius()) < Point::NullDist)
  865.     {
  866.         return lx.Project(cx.Center());
  867.     }
  868.  
  869. // line is outside circle
  870.     if(d > cx.Radius())
  871.         throw GeomException("GL2D_SolveLineCircleIntersection() : no intersection between Line and Circle");
  872.  
  873. // two solution cases
  874.     Point p1,p2;
  875. // passing through center
  876.     if(fabs(d) < Point::NullDist)
  877.     {
  878.         Vector vx = lx.Direction();
  879.         vx = vx.Normal();
  880.         p1 = cx.Center() + (cx.Radius() * vx);
  881.         p2 = cx.Center() - (cx.Radius() * vx);
  882.     }
  883.     else
  884.     {
  885.         Point px = lx.Project(cx.Center());
  886.         Vector v1 = px - cx.Center();
  887.         v1 = cx.Radius() * v1.Normal();
  888.         double theta = acos(d/cx.Radius());
  889.         p1 = cx.Center() + GL2D_Rotate(v1,theta);
  890.         p2 = cx.Center() + GL2D_Rotate(v1,-theta);
  891.     }
  892.  
  893. // take closet of 2 solutions to 2nd pick point 
  894.     if (pickloc.Distance(p1) < pickloc.Distance(p2)) return p1;
  895.     return p2;
  896. }
  897.  
  898.     
  899. //***********************************************************************************
  900. Point GL2D_SolveCircleCircleIntersection(const Circle& c1,const Circle& c2,const Point& pickloc)
  901.     throw (GeomException)
  902. {
  903.  
  904. // make local copies    ( must make locals so I can swap them if necessary)
  905.     Circle cx1(c1);
  906.     Circle cx2(c2);
  907.  
  908. // assure larger circle is cx1    
  909.     if(cx1.Radius() < cx2.Radius())
  910.     {
  911.         Circle tmp = cx2;
  912.         cx2 = cx1;
  913.         cx1 = tmp;
  914.     }
  915.  
  916. // circle centers distance
  917.     double dc = Distance(cx1.Center(), cx2.Center());
  918.  
  919.     // common centers or one circle completly inside the other - no solution
  920.     // dc > r1 + r2 --- no solution
  921.  
  922.     if( dc < Point::NullDist || 
  923.           (dc+ cx2.Radius()-cx1.Radius()) < -Point::NullDist ||
  924.       (dc- cx1.Radius()-cx2.Radius()) > Point::NullDist )
  925.           throw GeomException("GL2D_SolveCircleCircleIntersection() : no intersections between Circles");
  926.  
  927.     // one circle tangent to  another - one solution
  928.     if( fabs( dc - ( cx1.Radius() + cx2.Radius() ) ) < Point::NullDist ||
  929.         fabs( dc - ( cx1.Radius() - cx2.Radius() ) ) < Point::NullDist )
  930.     {
  931.         Vector v = cx2.Center() - cx1.Center();
  932.         v = v.Normal();
  933.         return cx1.Center() + v * cx1.Radius();
  934.     }
  935.  
  936.     // classical two solutions
  937.  
  938. // solve for trianlge height using Heron's formula
  939.     double s = 0.5 * ( dc + cx1.Radius() + cx2.Radius());
  940.     double h = 2 * sqrt( s * ( s - cx1.Radius() ) * ( s - cx2.Radius()) * ( s - dc ) ) /dc;
  941.  
  942. // now solve for solution points
  943.     Vector v = cx2.Center() - cx1.Center();
  944.     v = v.Normal();
  945.     double theta = asin( h / cx1.Radius());
  946.  
  947.     Point p1 = cx1.Center() +  cx1.Radius() * GL2D_Rotate(v,theta) ;
  948.     Point p2 = cx1.Center() +  cx1.Radius() * GL2D_Rotate(v,-theta) ;
  949.  
  950. // take closet of 2 solutions to 2nd pick point 
  951.     if (pickloc.Distance(p1) < pickloc.Distance(p2)) return p1;
  952.     return p2;
  953. }
  954.  
  955. //***********************************************************************************
  956. Line GL2D_SolveLineTangentCircleParallelLine(const Circle& c, const Line& l, const Point& p)
  957.     throw (GeomException)
  958. {
  959. // degenerate cases:
  960. // selection point at circle center - no way to correctly interpret.
  961.     if(p.Distance(c.Center()) < Point::NullDist )
  962.         throw GeomException("GL2D_SolveLineTangentCircleParallelLine() : selection point is at Circle center");
  963.  
  964. // line passes through circle center
  965.     if(l.Distance(c.Center()) < Point::NullDist )
  966.         throw GeomException("GL2D_SolveLineTangentCircleParallelLine() : Line passes through Circle center");
  967.  
  968.     Point px;
  969. // unit vector from center parallel to line
  970.     Vector v = l.Project(c.Center()) - c.Center();
  971.     v = v.Normal();
  972. // vector from center to selection point
  973.     Vector vx = p - c.Center();
  974.  
  975.     if(( v * vx ) > 0 )
  976.         return Line(c.Center() + ( c.Radius() * v ), l.Direction());
  977.     else
  978.         return Line(c.Center() - ( c.Radius() * v ), l.Direction());
  979. }
  980.  
  981. //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  982. Point GL2D_PointByModifier(const Pick& pk,int modifier) throw(GeomException)
  983. {
  984.     int type = pk.Type();
  985.  
  986.     switch(modifier)
  987.     {
  988.         case GL2D_MIDPOINT:
  989.         {    
  990.             switch(type)
  991.             {
  992.                 case POINT:
  993.                 {
  994.                     return GetPointFromPick(pk);
  995.                 }
  996.                 case SEGMENT:
  997.                 {
  998.                     Segment s = GetSegmentFromPick(pk);
  999.                     return s.U(0.5);
  1000.                 }
  1001.                 case ARC:
  1002.                 {
  1003.                     Arc a = GetArcFromPick(pk);
  1004.                     return a.U(0.5);
  1005.                 }
  1006.             }
  1007.             throw GeomException("GL2D_PointByModifier() : midpoint is not a valid modifier for the selected type");
  1008.             break;
  1009.         }
  1010.         case GL2D_ENDPOINT:
  1011.         {
  1012.             switch(type)
  1013.             {
  1014.                 case POINT:
  1015.                 {    
  1016.                     return GetPointFromPick(pk);
  1017.                 }
  1018.                 case SEGMENT:
  1019.                 {
  1020.                     Segment s = GetSegmentFromPick(pk);
  1021.                     Point p0 = s.U(0.0);
  1022.                     Point p1 = s.U(1.0);
  1023.                     if(Distance(p0,pk.GetPoint()) < Distance(p1,pk.GetPoint()) )
  1024.                         return p0;
  1025.                     return p1;
  1026.                 }
  1027.                 case ARC:
  1028.                 {
  1029.                     Arc a = GetArcFromPick(pk);
  1030.                     Point p0 = a.U(0.0);
  1031.                     Point p1 = a.U(1.0);
  1032.                     if(Distance(p0,pk.GetPoint()) < Distance(p1,pk.GetPoint()) )
  1033.                         return p0;
  1034.                     return p1;
  1035.                 }    
  1036.             }
  1037.             throw GeomException("GL2D_PointByModifier() : endpoint is not a valid modifier for the selected type");
  1038.             break;
  1039.         }
  1040.         case GL2D_CENTER:
  1041.         {
  1042.             switch(type)
  1043.             {
  1044.                 case POINT:
  1045.                 {
  1046.                     return GetPointFromPick(pk);
  1047.                 }
  1048.                 case CIRCLE:
  1049.                 {
  1050.                     Circle c = GetCircleFromPick(pk);
  1051.                     return c.Center();
  1052.                 }
  1053.                 case ARC:
  1054.                 {
  1055.                     Arc a = GetArcFromPick(pk);
  1056.                     return a.Center();
  1057.                 }
  1058.             }
  1059.             throw GeomException("GL2D_PointByModifier() : center is not a valid modifier for the selected type");
  1060.             break;
  1061.         }
  1062.         case GL2D_NEAREST:
  1063.         {
  1064.             switch(type)
  1065.             {
  1066.                 case POINT:
  1067.                 {
  1068.                     return Point();
  1069.                 }
  1070.                 case LINE:
  1071.                 {
  1072.                     return Point();
  1073.                 }
  1074.                 case CIRCLE:
  1075.                 {
  1076.                     return Point();
  1077.                 }    
  1078.                 case SEGMENT:
  1079.                 {
  1080.                     return Point();
  1081.                 }
  1082.                 case ARC:
  1083.                 {
  1084.                     return Point();
  1085.                 }
  1086.             }
  1087.             throw GeomException("GL2D_PointByModifier() : endpoint is not a valid modifier for the selected type");
  1088.             break;
  1089.         }
  1090.         break;
  1091.     }
  1092.     throw GeomException("GL2D_PointByModifier() : invalid modifier value");
  1093. }
  1094.  
  1095.  
  1096. Circle  GL2D_CircleRadiusAndConstraints(double r, const Pick& pk1, const Pick& pk2) 
  1097.     throw (GeomException)
  1098. {
  1099. // cases
  1100. // 1 line 1 circle - minimizing total distance...
  1101. // 2 circle minimize total distance
  1102.  
  1103. // 2 points - right handed circle from p1 to p2
  1104.     if(pk1.Type() == POINT && pk2.Type() == POINT)
  1105.     {
  1106.         // failures : points are the same, distance between points is more than 2r 
  1107.         Point p1 = GetPointFromPick(pk1);
  1108.         Point p2 = GetPointFromPick(pk2);
  1109.         return GL2D_SolveCircleTanPointPointRadius(p1,p2,r);
  1110.     }
  1111. // 1 point 1 line - return solution closest to where line was picked
  1112. // 1 point 1 circle - solution closest to where circle was picked
  1113.     else if(pk1.Type()==POINT || pk2.Type() == POINT)
  1114.     {
  1115.  
  1116. // point vs.line/segment
  1117.         if(pk1.Type()==LINE || pk1.Type()==SEGMENT || pk2.Type()==LINE || pk2.Type()==SEGMENT)
  1118.         {
  1119. // get the data
  1120.             Point p;
  1121.             Line l;
  1122.             Point pkloc;
  1123.             if(pk1.Type() == POINT)
  1124.             {
  1125.                 p = GetPointFromPick(pk1);                
  1126.                 l = GetLineFromPick(pk2);
  1127.                 pkloc = pk2.GetPoint();
  1128.             }
  1129.             else
  1130.             {    
  1131.                 p = GetPointFromPick(pk2);                
  1132.                 l = GetLineFromPick(pk1);
  1133.                 pkloc = pk1.GetPoint();
  1134.             }
  1135.             return GL2D_SolveCircleTanPointLineRadius(p,l,r,pkloc);
  1136.         }
  1137.         else if(pk1.Type()==CIRCLE || pk1.Type()==ARC || pk2.Type()==CIRCLE || pk2.Type()==ARC)
  1138.         {
  1139.         // get circle and point
  1140.             Point p,px;
  1141.             Circle c,coffset,cx1,cx2;
  1142.             Point pkloc;
  1143.             if(pk1.Type() == POINT)
  1144.             {
  1145.                 p = GetPointFromPick(pk1);                
  1146.                 c = GetCircleFromPick(pk2);
  1147.                 pkloc = pk2.GetPoint();
  1148.             }
  1149.             else
  1150.             {    
  1151.                 p = GetPointFromPick(pk2);                
  1152.                 c = GetCircleFromPick(pk1);
  1153.                 pkloc = pk1.GetPoint();
  1154.             }
  1155.             return GL2D_SolveCircleTanPointCircleRadius(p,c,r,pkloc);
  1156.         }
  1157.     }
  1158. // 2 lines - take midpoint of line selection location, offset lines on that side by r, find intersection 
  1159.     else if((pk1.Type() == LINE || pk1.Type() == SEGMENT) && ( pk2.Type() == LINE || pk2.Type() ==SEGMENT))
  1160.     {
  1161.         Line l1,l2;
  1162.         Point mpt;
  1163.  
  1164.         l1 = GetLineFromPick(pk1);
  1165.         l2 = GetLineFromPick(pk2);
  1166.         mpt = pk1.GetPoint() + 0.5 * (pk2.GetPoint() - pk1.GetPoint());
  1167.  
  1168.         return GL2D_SolveCircleTanLineLineRadius(l1,l2,r,mpt);
  1169.     }
  1170. // circle/circle
  1171.       else if((pk1.Type() == CIRCLE || pk1.Type() == ARC ) && (pk2.Type() == CIRCLE || pk2.Type() == ARC))
  1172.       {
  1173.         Circle c1,c2;
  1174.         Point pkloc1, pkloc2;
  1175.         c1 = GetCircleFromPick(pk1);
  1176.         c2 = GetCircleFromPick(pk2);
  1177.         pkloc1 = pk1.GetPoint();
  1178.         pkloc2 = pk2.GetPoint();
  1179.         return GL2D_SolveCircleTanCircleCircleRadius(c1,c2,r,pkloc1,pkloc2);
  1180.       }
  1181.       else if(((pk1.Type() == CIRCLE || pk1.Type() == ARC ) && (pk2.Type() == LINE || pk2.Type() == SEGMENT)) ||
  1182.               ((pk2.Type() == CIRCLE || pk2.Type() == ARC ) && (pk1.Type() == LINE || pk1.Type() == SEGMENT)))
  1183.     {
  1184.         Line l;
  1185.         Circle c;
  1186.         Point pkloc_c, pkloc_l;
  1187.         if(pk1.Type() == LINE || pk1.Type() == SEGMENT)
  1188.         {
  1189.             l = GetLineFromPick(pk1);                
  1190.             pkloc_l = pk1.GetPoint();
  1191.             c = GetCircleFromPick(pk2);
  1192.             pkloc_c = pk2.GetPoint();
  1193.         }
  1194.         else
  1195.         {    
  1196.             l = GetLineFromPick(pk2);                
  1197.             pkloc_l = pk2.GetPoint();
  1198.             c = GetCircleFromPick(pk1);
  1199.             pkloc_c = pk1.GetPoint();
  1200.         }
  1201.         return GL2D_SolveCircleTanLineCircleRadius(l,c,r,pkloc_l,pkloc_c);
  1202.     }
  1203.     throw GeomException("GL2D_CircleByRadiusAndConstraints() : combination not implemented");
  1204. }
  1205. //=======================================================================================
  1206.  
  1207. Circle GL2D_SolveCircleTanPointCircleRadius(const Point& p, const Circle& c, const double& r, const Point& pkloc)
  1208.     throw (GeomException)
  1209. {
  1210.     Circle coffset;
  1211.     Point px;
  1212.     double d;
  1213.     Circle cx1,cx2;
  1214.  
  1215. // get distance from point to circle
  1216.     d = c.Distance(p);
  1217.  
  1218. // if point is on circle -> no solution (questionable...) 
  1219.     if(d < Point::NullDist)
  1220.         throw GeomException("GL2D_CircleTanPointCircleRadius() : Point on circle -> no solution");
  1221. // if distance from circle to point is > 2r -> no solution
  1222.     if(d > (2 * r + Point::NullDist)) 
  1223.         throw GeomException("GL2D_CircleTanPointCircleRadius() :: Point too far from circle -> no solution");
  1224. // if point is at circle center, infinite solutions
  1225.     if(p == c.Center()) 
  1226.         throw GeomException("GL2D_CircleTanPointCircleRadius() : Point is at circle center, infinite solutions");
  1227.     d = p.Distance(c.Center());
  1228. // is point inside circle?
  1229.     if( d < c.Radius())
  1230.     {
  1231. // if point is inside circle and r > circle radius -> no solution possible
  1232.         if(r > c.Radius())
  1233.             throw GeomException("GL2D_CircleTanPointCircleRadius() : Point inside circle and r > circle radius -> no solution");
  1234. // otherwise 1 of 2 solutions
  1235.         coffset = Circle(c.Center(), c.Radius() - r);
  1236.         px = GL2D_SolveCircleCircleIntersection(coffset,Circle(p,r),pkloc);
  1237.         return Circle(px,r);
  1238.     }
  1239. // point is outside circle - there are from 2 to 4 solutions remaining...
  1240. // is there a possible enclosing solution?
  1241. // if the diameter of the requested circle is less than or equal to the 
  1242. // distance to the point + the radius of the existing circle...
  1243.     bool multi = false;
  1244.     if((2 * r - d - c.Radius()  ) > -(Point::NullDist) )
  1245.     {
  1246. // yup
  1247.         coffset = Circle(c.Center(), fabs(c.Radius() - r)); // interior offset circle
  1248.     // pkloc for solution selection should be opposite side of circle
  1249.         Line lx(p,c.Center() - p );
  1250.         Point pkmirror = 2 * lx.Project(pkloc) - pkloc;
  1251.         px = GL2D_SolveCircleCircleIntersection(coffset,Circle(p,r),pkmirror);
  1252.         cx1 = Circle(px,r);
  1253.         multi = true;
  1254.     }
  1255.     // the 'near' solution
  1256.     coffset = Circle(c.Center(),c.Radius() + r);
  1257.     px = GL2D_SolveCircleCircleIntersection(coffset,Circle(p,r),pkloc);
  1258.     // return solution nearest to circle pick point
  1259.     cx2 = Circle(px,r);
  1260.     if(multi && cx1.Distance(pkloc) < cx2.Distance(pkloc)) return cx1;
  1261.     return cx2;
  1262. }
  1263.  
  1264. //=======================================================================================
  1265. Circle GL2D_SolveCircleTanPointLineRadius(const Point& p, const Line& l, const double& r, const Point& pkloc) 
  1266.     throw (GeomException)
  1267. {
  1268. // point on line - ambiguous - no solution ( could choose from side of line, ehhh...)
  1269.     if(l.Distance(p) < Point::NullDist)
  1270.         throw GeomException("GL2D_CircleByTanPointLineRadius() : Ambiguous solution > point is on line");
  1271.     if(l.Distance(p) > 2 * r + Point::NullDist)
  1272.         throw GeomException("GL2D_CircleByTanPointLineRadius() : no solution > point is more than 2r from line");
  1273. // solve the case
  1274. // offset the line 
  1275.     Vector v = p - l.Project(p);
  1276.     v = v.Normal();
  1277.     Line offset(l.Origin() + r * v,l.Direction());
  1278. // circle around r
  1279.     Circle cx(p,r);
  1280.     Point pc = GL2D_SolveLineCircleIntersection(offset,cx,pkloc);
  1281.     return Circle(pc,r);
  1282. }
  1283. //============================================================================================
  1284. Circle GL2D_SolveCircleTanLineLineRadius(const Line& l1, const Line& l2, const double& r, const Point& pkloc) 
  1285.     throw (GeomException)
  1286. {
  1287.     Line lx1,lx2;
  1288.     Vector v;
  1289.     Point px;
  1290.  
  1291.     if(Angle(l1.Direction(), l2.Direction()) < Point::NullAngle)
  1292.         throw GeomException("GL2D_CircleTanLineLineRadius() : lines are parallel -> no solution");
  1293. // midpoint between picks
  1294. // offset first line
  1295.     v = pkloc - l1.Project(pkloc);
  1296.     v = v.Normal();
  1297.     lx1 = Line(l1.Origin() + (r * v), l1.Direction());
  1298. // offset second line
  1299.     v = pkloc - l2.Project(pkloc);
  1300.     v = v.Normal();
  1301.     lx2 = Line(l2.Origin() + (r * v), l2.Direction());
  1302.  
  1303. // get intersection
  1304.     px = GL2D_SolveLineLineIntersection(lx1,lx2);
  1305.     return Circle(px,r);
  1306. }
  1307. //============================================================================================
  1308. Circle GL2D_SolveCircleTanPointPointRadius(const Point& p1, const Point& p2, const double& r)
  1309.     throw (GeomException)
  1310. {
  1311.     if(p1 == p2 || p1.Distance(p2) > (2 * r )) 
  1312.         throw GeomException("GL2D_CircleTanPointPointRadius() : No solution for the selected points");
  1313.     Circle c1(p1,r);
  1314.     Circle c2(p2,r);
  1315.     Vector v = p2 - p1;
  1316.     v = GL2D_Rotate(v,90);
  1317.     Point px = GL2D_SolveCircleCircleIntersection(c1,c2,p2+v);
  1318.     return Circle(px,r);
  1319. }
  1320. //============================================================================================
  1321. Circle GL2D_SolveCircleTanCircleCircleRadius(const Circle& c1, const Circle& c2, const double& r, const Point& pkloc1, const Point& pkloc2)
  1322.     throw (GeomException)
  1323. {
  1324. // try this :   consider the Rsmall, radius of the smaller of the two circles Csmall, and Psmall it's
  1325. // center point. find the centers of (up to) 4 circles :
  1326. // passing Psmall, tangent Cbig + Rsmall, radius r + Rsmall
  1327. // passing Psmall, tangent Cbig + Rsmall, radius r - Rsmall
  1328. // passing Psmall, tangent Cbig - Rsmall, radius r + Rsmall
  1329. // passing Psmall, tangent Cbig - Rsmall, radius r - Rsmall
  1330. // constitute 4 solution circles from these centers and choose the best to return
  1331.     vector<Circle> sols;
  1332.     Circle cx,czzz;
  1333.     double Rsmall,Rzzz;
  1334.     Point Psmall,pkbig,pkzzz;
  1335.     Circle Csmall,Cbig;
  1336.     Vector v;
  1337.  
  1338.     if(c1.Radius() < c2.Radius())
  1339.     {
  1340.         Csmall = c1;
  1341.         Cbig = c2;
  1342.         pkbig = pkloc2;
  1343.     }
  1344.     else
  1345.     {
  1346.         Csmall = c2;
  1347.         Cbig = c1;
  1348.         pkbig = pkloc1;
  1349.     }
  1350.     Psmall = Csmall.Center();
  1351.     Rsmall = Csmall.Radius();
  1352.     v = pkbig - Cbig.Center();    // vector from center of Cbig to pkpoint on Cbig, normalized
  1353.     v = v.Normal();
  1354. // passing Psmall, tangent Cbig + Rsmall, radius r + Rsmall
  1355.     czzz = Circle(Cbig.Center(), Cbig.Radius() + Rsmall);
  1356.     Rzzz = r + Rsmall;
  1357.     pkzzz = pkbig + Rsmall * v;
  1358.     try
  1359.     {
  1360.         cx = GL2D_SolveCircleTanPointCircleRadius(Psmall,czzz,Rzzz,pkzzz);
  1361.         sols.push_back(Circle(cx.Center(),r));
  1362.     }
  1363.     catch (GeomException& ge)
  1364.     {
  1365.     }
  1366. // passing Psmall, tangent Cbig + Rsmall, radius r - Rsmall
  1367.     czzz = Circle(Cbig.Center(), Cbig.Radius() + Rsmall);
  1368.     Rzzz = r - Rsmall;
  1369.     try
  1370.     {
  1371.         cx = GL2D_SolveCircleTanPointCircleRadius(Psmall,czzz,Rzzz,pkzzz);
  1372.         sols.push_back(Circle(cx.Center(),r));
  1373.     }
  1374.     catch (GeomException& ge)
  1375.     {
  1376.     }
  1377. // passing Psmall, tangent Cbig - Rsmall, radius r + Rsmall
  1378.     czzz = Circle(Cbig.Center(), Cbig.Radius() - Rsmall);
  1379.     Rzzz = r + Rsmall;
  1380.     pkzzz = pkbig - Rsmall * v;
  1381.     try
  1382.     {
  1383.         cx = GL2D_SolveCircleTanPointCircleRadius(Psmall,czzz,Rzzz,pkzzz);
  1384.         sols.push_back(Circle(cx.Center(),r));
  1385.     }
  1386.     catch (GeomException& ge)
  1387.     {
  1388.     }
  1389. // passing Psmall, tangent Cbig - Rsmall, radius r - Rsmall
  1390.     czzz = Circle(Cbig.Center(), Cbig.Radius() - Rsmall);
  1391.     Rzzz = r - Rsmall;
  1392.     try
  1393.     {
  1394.         cx = GL2D_SolveCircleTanPointCircleRadius(Psmall,czzz,Rzzz,pkzzz);
  1395.         sols.push_back(Circle(cx.Center(),r));
  1396.     }
  1397.     catch (GeomException& ge)
  1398.     {
  1399.     }
  1400. // figure out the best of the solution points
  1401.     if(sols.size() == 0) throw GeomException("Gl2D_SolveCircleTanCircleCircleRadius() : no solutions found");
  1402.     double minDist = 100000.;
  1403.     vector<Circle>::iterator minIx;
  1404.     for(vector<Circle>::iterator ix = sols.begin(); ix != sols.end() ; ix++)
  1405.     {
  1406.         double totDist = ix->Distance(pkloc1) + ix->Distance(pkloc2);
  1407.         if(totDist < minDist)
  1408.         {
  1409.             minDist = totDist;
  1410.             minIx = ix;
  1411.         }
  1412.     }
  1413.     return Circle(*minIx);
  1414.             
  1415.  
  1416. //    throw GeomException("GL2D_SolveCircleTanCircleCircleRadius() : function not implemented");
  1417. }
  1418. //============================================================================================
  1419. Circle GL2D_SolveCircleTanLineCircleRadius(const Line& l, const Circle& c,const double& r, const Point& pkloc_l, const Point& pkloc_c)
  1420.     throw (GeomException)
  1421. {
  1422.     Circle cx,cs1,cs2;
  1423.     Point ps1,ps2;
  1424.     bool multi = false;
  1425.     double d = l.Distance(c.Center());
  1426. // no solution if... circle is more than r from line
  1427.     if(( d - c.Radius() -r ) > Point::NullDist)
  1428.         throw GeomException("GL2D_SolveCircleTanLineCircleRadius() : No solution -> circle too far from line");
  1429. // 2 solutions if line intersects circle, 4 otherwise
  1430. // both solutions require offsets of line and circle
  1431. // offset the line on the circle selection side
  1432.     Vector v = pkloc_c - l.Project(pkloc_c);
  1433.     v = v.Normal();
  1434.     Line lx(l.Origin() + r * v , l.Direction());
  1435. // offset the circle on the outside
  1436.     cx = Circle(c.Center(), c.Radius() + r);
  1437.  
  1438. // get the first solution ( this automatically resolves from 2 solutions based on pkloc_l )
  1439.     ps1 = GL2D_SolveLineCircleIntersection(lx,cx,pkloc_l);
  1440.     cs1 = Circle(ps1,r);
  1441.  
  1442. // check for 4 solutions - if the circle does not intersect the line AND the circle is not 
  1443. // too far from the line, there are an extra 2 solutions
  1444.     if(d > c.Radius() &&  (d + c.Radius()) < r )    
  1445.     {
  1446.         cx = Circle(c.Center(), fabs(c.Radius() - r));
  1447.         ps2 = GL2D_SolveLineCircleIntersection(lx,cx,pkloc_l);
  1448.         cs2 = Circle(ps2,r);
  1449.         multi = true;
  1450.     }
  1451. // return the correct solution - closest to pkloc_c 
  1452.     if(multi  && cs2.Distance(pkloc_c) < cs1.Distance(pkloc_c)) return cs2;
  1453.     return cs1;
  1454. }
  1455.