Logo Search packages:      
Sourcecode: sagcad version File versions  Download package

Ellipse.c

/* ====================================================================
 * ===  Copyright (C) 1998-2007 Yutaka Sagiya. All rights reserved. ===
 * ====================================================================
 * 
 *    Project              : SagCAD
 *    Description          : CAD/CAM
 *    Source               : Ellipse.c
 * 
 *    ----------------------------------
 * 
 *    License              : GNU General Public License (GPL)
 *    Copyright            : (C) 1998-2007 by Yutaka Sagiya
 *    email                : kappa@a6s.highway.ne.jp
 *                         : yutaka@sagiya.com
 *    Begin                : 2003/04/22
 *    Last                 : 2007/11/08
 * ====================================================================
 */

#ifdef HAVE_CONFIG_H
#  include "config.h"
#endif

#include <gtk/gtk.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "MemoryLeak.h"
#include "types.h"
#include "List_cad.h"
#include "List_Ellipse.h"
//#include "CopyFunc.h"
#include "culcfunc.h"

#define _ELLIPSE_
#include "Ellipse.h"





/* -------------------------------------------------------------------
 * 長軸 a を求める
 * 
 */
double get_ellipse_a(ELLIPSE ellipse)
{
      struct RtnDat rtn;

      rtn.sx[1] = ellipse.cx;
      rtn.sy[1] = ellipse.cy;
      rtn.ex[1] = ellipse.cx + ellipse.dx;
      rtn.ey[1] = ellipse.cy + ellipse.dy;
      pp(&rtn);
      return rtn.l;
}


/* -------------------------------------------------------------------
 * 短軸 b を求める
 * 
 */
double get_ellipse_b(ELLIPSE ellipse)
{
      return (get_ellipse_a(ellipse) * ellipse.k);
}


/* -------------------------------------------------------------------
 * 楕円の傾きを求める
 * 
 */
double get_ellipse_angle(ELLIPSE ellipse)
{
      struct RtnDat rtn;


      rtn.sx[1] = ellipse.cx;
      rtn.sy[1] = ellipse.cy;
      rtn.ex[1] = ellipse.cx + ellipse.dx;
      rtn.ey[1] = ellipse.cy + ellipse.dy;
      la(&rtn);
      angle_check(&rtn.angle);
      return rtn.angle;
}


/* -------------------------------------------------------------------
 * 楕円の始点と終点を求める
 * 
 */
int get_ellipse_se_point(ELLIPSE ellipse, SAG_POINT *s_point, SAG_POINT *e_point)
{
      double RAD, DX, DY;




      /* 始点 */
      RAD = degrad(get_ellipse_angle(ellipse));

      s_point->x = get_ellipse_a(ellipse) * cos(degrad(ellipse.sa)) + ellipse.cx;
      s_point->y = get_ellipse_b(ellipse) * sin(degrad(ellipse.sa)) + ellipse.cy;

      DX = s_point->x - ellipse.cx;
      DY = s_point->y - ellipse.cy;
      s_point->x = sg((DX * cos(RAD) - DY * sin(RAD)) + ellipse.cx, calcu_digits);
      s_point->y = sg((DY * cos(RAD) + DX * sin(RAD)) + ellipse.cy, calcu_digits);


      /* 終点 */
      RAD = degrad(get_ellipse_angle(ellipse));

      e_point->x = get_ellipse_a(ellipse) * cos(degrad(ellipse.ea)) + ellipse.cx;
      e_point->y = get_ellipse_b(ellipse) * sin(degrad(ellipse.ea)) + ellipse.cy;

      DX = e_point->x - ellipse.cx;
      DY = e_point->y - ellipse.cy;
      e_point->x = sg((DX * cos(RAD) - DY * sin(RAD)) + ellipse.cx, calcu_digits);
      e_point->y = sg((DY * cos(RAD) + DX * sin(RAD)) + ellipse.cy, calcu_digits);

      return 1;
}


/* -------------------------------------------------------------------
 * 楕円の焦点を求める
 * 
 */
int get_ellipse_focus(ELLIPSE ellipse, SAG_POINT *focus1, SAG_POINT *focus2)
{
      double a, b, angle;


      a = get_ellipse_a(ellipse);
      b = get_ellipse_b(ellipse);
      angle = get_ellipse_angle(ellipse);

      /* 焦点を求める */
      focus1->x = -1 * (sqrt(a*a - b*b)) + ellipse.cx;
      focus1->y = 0 + ellipse.cy;
      focus2->x = sqrt(a*a - b*b) + ellipse.cx;
      focus2->y = 0 + ellipse.cy;

      rotation(&focus1->x, &focus1->y, ellipse.cx, ellipse.cy, angle);
      rotation(&focus2->x, &focus2->y, ellipse.cx, ellipse.cy, angle);

      return 1;
}





/* -------------------------------------------------------------------
 * 楕円と角度で楕円上の点を求める
 * 
 */
int get_ellipse_point_from_angle (ELLIPSE ellipse, double angle, SAG_POINT *point)
{
//#define GET_ELLIPSE_POINT_FROM_ANGLE
      struct RtnDat rtn;
      SAG_POINT point1, point2;
      double angle1, angle2;
      ELLIPSE eps;
      int debug = 0;


#ifdef GET_ELLIPSE_POINT_FROM_ANGLE
      debug = 1;
#endif


      if (sg(angle, 3) == 360) angle = 0;
      if (debug > 0) g_print ("get_ellipse_point_from_angle() : in\n");
      eps = ellipse;
      eps.dx = get_ellipse_a(ellipse);
      eps.dy = 0;


      rtn.sx[1] = ellipse.cx;
      rtn.sy[1] = ellipse.cy;
      rtn.angle = angle;
      rtn.l = 50;
      pap(&rtn);
      rtn.type = 0;

      lep(&rtn, eps, 10);

      point1.x = rtn.sx[2];
      point1.y = rtn.sy[2];
      point2.x = rtn.sx[3];
      point2.y = rtn.sy[3];

      if (debug > 0) g_print ("    p1 (%f,%f)     p2 (%f,%f)\n", point1.x, point1.y, point2.x, point2.y);

      
      rtn.sx[1] = ellipse.cx;
      rtn.sy[1] = ellipse.cy;
      rtn.ex[1] = point1.x;
      rtn.ey[1] = point1.y;
      la(&rtn);
      angle1 = rtn.angle;

      rtn.sx[1] = ellipse.cx;
      rtn.sy[1] = ellipse.cy;
      rtn.ex[1] = point2.x;
      rtn.ey[1] = point2.y;
      la(&rtn);
      angle2 = rtn.angle;

      if (debug > 0) g_print ("    angle1 = %f   angle2 = %f   angle = %f\n", angle1, angle2, angle);

      if ((angle1 + 1) > angle && (angle1 - 1) < angle) {
            rotation(&point1.x, &point1.y, ellipse.cx, ellipse.cy, get_ellipse_angle(ellipse));
            point->x = point1.x;
            point->y = point1.y;
            if (debug > 0) g_print ("    (%f,%f)\n", point->x, point->y);
            return 1;
      }
      else if ((angle2 + 1) > angle && (angle2 - 1) < angle) {
            rotation(&point2.x, &point2.y, ellipse.cx, ellipse.cy, get_ellipse_angle(ellipse));
            point->x = point2.x;
            point->y = point2.y;
            if (debug > 0) g_print ("    (%f,%f)\n", point->x, point->y);
            return 1;
      }



      if (debug > 0) g_print ("get_ellipse_point_from_angle() : out\n");
      return 0;
}





/* -------------------------------------------------------------------
 * 楕円と点で楕円上の角度を求める
 * 
 */
double get_ellipse_angle_from_point (ELLIPSE ellipse, SAG_POINT point)
{
      struct RtnDat rtn;


      rotation(&point.x, &point.y, ellipse.cx, ellipse.cy, (-1 * get_ellipse_angle(ellipse)));

      /* -------------------------------------------------------------------
       * la : 直線の角度
       */
      rtn.sx[1] = ellipse.cx;
      rtn.sy[1] = ellipse.cy;
      rtn.ex[1] = point.x;
      rtn.ey[1] = point.y;
      la(&rtn);
      angle_check(&rtn.angle);

      return rtn.angle;
}





/* -------------------------------------------------------
 * 点 (px, py) は、楕円上の点であるかの判定
 * 
 */
int point_on_ellipse (ELLIPSE *ellipse, double px, double py)
{
//#define POINT_ON_ELLIPSE
      struct RtnDat rtn;
      double angle = 0;
      double pic_angle;
      int debug = 0;


#ifdef POINT_ON_ELLIPSE
      debug = 1;
      if (debug > 0) g_print("List_Ellipse.c : point_on_ellipse() : in\n");
#endif


      /* 点が楕円上にあるか */
      if (point_in_ellipse(ellipse, px, py, 0.000001) != 3) {
            return 0;
      }


      /* 開始角・終了角で判定 */
      /* ピック点の楕円に対する角度を求める */
      /* -------------------------------------------------------------------
       * la : [8] 直線の角度
       */
      rtn.sx[1] = ellipse->cx;
      rtn.sy[1] = ellipse->cy;
      rtn.ex[1] = px;
      rtn.ey[1] = py;
      la(&rtn);
      pic_angle = rtn.angle - angle;
      angle_check(&pic_angle);


      if (debug > 0) {
            g_print("Ellipse.c : point_on_ellipse() : pic_angle = %f  Ellipse sa = %f  ea = %f\n", 
                         pic_angle, ellipse->sa, ellipse->ea);
      }


      /* SA < EA */
      if (ellipse->sa < ellipse->ea) {
            if (debug > 0) g_print("Ellipse.c : point_on_ellipse() : sa (%f) < ea (%f)\n", ellipse->sa, ellipse->ea);

            if (pic_angle > ellipse->sa && pic_angle < ellipse->ea) {
                  if (debug > 0) g_print("Ellipse.c : point_on_ellipse() : sa < pic (%f) < ea ---> OK !!\n", pic_angle);
                  return 1;
            }
      }

      /* SA > EA */
      else if (ellipse->sa > ellipse->ea) {
            if (debug > 0) g_print("Ellipse.c : point_on_ellipse() : sa (%f) > ea (%f)\n", ellipse->sa, ellipse->ea);

            if (pic_angle > ellipse->sa || pic_angle < ellipse->ea) {
                  if (debug > 0) g_print("Ellipse.c : point_on_ellipse() : pic (%f) > sa  pic (%f) < ea ---> OK !!\n", 
                              pic_angle, pic_angle);
                  return 1;
            }
      }
      else {
            return 0;
      }

      return 0;
}





/* -------------------------------------------------------
 * 点が楕円の中にあるか外にあるかの判定
 * 
 * 
 * 
 * value : 判定範囲(誤差)    0.001 とか 0.5 とか
 * 
 * return : 1 : IN
 *          2 : OUT
 *          3 : ON
 */
int point_in_ellipse (ELLIPSE *ellipse, double px, double py, double value)
{
//#define POINT_IN_ELLIPSE
      struct RtnDat rtn;
      double a, b, angle;
      double l1, l2;
      SAG_POINT focus1, focus2;
      int debug = 0;


#ifdef POINT_IN_ELLIPSE
      debug = 1;
      if (debug > 0) g_print("Ellipse.c : point_in_ellipse() : in\n");
#endif



      /* -----------------------------------------------------
       * 楕円の要素を求める
       */
      /* 中心から長軸の位置の角度と距離を求める */
      angle = get_ellipse_angle(*ellipse);
      a = get_ellipse_a(*ellipse);
      b = get_ellipse_b(*ellipse);
      get_ellipse_focus(*ellipse, &focus1, &focus2);

      if (debug > 0) {
            g_print("Ellipse.c : point_in_ellipse() : a = %f  b = %f  angle = %f  焦点 (%f,%f) (%f,%f)\n", 
                         a, b, angle, focus1.x, focus1.y, focus2.x, focus2.y);
      }



      /* (px,py) と 焦点の距離の関係で、(1/pow(10, compa_digits)) で判定 */
      /* -------------------------------------------------------------------
       * pp : [9] 2点間の距離
       */
      rtn.sx[1] = focus1.x;
      rtn.sy[1] = focus1.y;
      rtn.ex[1] = px;
      rtn.ey[1] = py;
      pp(&rtn);
      l1 = rtn.l;

      rtn.sx[1] = focus2.x;
      rtn.sy[1] = focus2.y;
      rtn.ex[1] = px;
      rtn.ey[1] = py;
      pp(&rtn);
      l2 = rtn.l;

      if (debug > 0) g_print("Ellipse.c : point_in_ellipse() : X10000 : (l1 + l2) = %f   (2 X a) = %f   value = %f\n", 10000*(l1 + l2), 10000*(2*a), 10000*value);



      /* (l1 + l2) = 2a なら Ellipse の円周上の点 */
      if ( (l1 + l2) + (value) > (2*a) && (l1 + l2) - (value) < (2*a) ) {
            if (debug > 0) g_print("Ellipse.c : point_in_ellipse() : on Ellipse\n");
            return 3;
      }

      /* Ellipse の中の点 */
      else if (l1 + l2 < (2*a)) {
            if (debug > 0) g_print("Ellipse.c : point_in_ellipse() : in Ellipse\n");
            return 1;
      }

      /* Ellipse の外の点 */
      else if (l1 + l2 > (2*a)) {
            if (debug > 0) g_print("Ellipse.c : point_in_ellipse() : out Ellipse\n");
            return 2;
      }

      return 0;
}





/* -------------------------------------------------------------------
 * 点から楕円の接線の交点
 * 
 * in :
 *      Point (a->sx[1],a->sy[1])
 * 
 * out :
 *      type 0,1,2
 *      Point 1 (a->ex[1],a->ey[1])
 *      Point 2 (a->ex[2],a->ey[2])
 * 
 */
int pel(struct RtnDat *rtn, ELLIPSE ellipse)
{
      int Ret;
      ELLIPSE eps;
      CAD cad;
      double a, b, A, B, C, X1, Y1, X2, Y2;


      /* -----------------------------------------------------
       * 点の位置を確認
       */
      eps = ellipse;
      Ret = point_in_ellipse (&eps, rtn->sx[1], rtn->sy[1], 0.0000001);
      /* in */
      if (Ret == 1) {
            rtn->type = 0;
            return 0;
      }
      /* out */
      else if (Ret == 2) {
            rtn->type = 2;
      }
      /* out */
      else if (Ret == 3) {
            rtn->type = 1;
            rtn->ex[1] = rtn->sx[1];
            rtn->ey[1] = rtn->sy[1];
            return 1;
      }



      /* -----------------------------------------------------
       * まず、楕円の中心が (0,0) で、角度が 0 になるような座標に変換。
       */
      /* 点 */
      init_cad(&cad);
      cad.sx = rtn->sx[1];
      cad.sy = rtn->sy[1];

      parallel(&cad.sx, &cad.sy, -ellipse.cx, -ellipse.cy);
      rotation(&cad.sx, &cad.sy, 0, 0, (-1 * get_ellipse_angle(ellipse)));



      /* 交点を求める */
      /* 楕円の式と接線の式から */
      a = get_ellipse_a(eps);
      b = get_ellipse_b(eps);

      A = (1/(a*a)) + (b * b * cad.sx * cad.sx) / (a * a * a * a * cad.sy * cad.sy);
      B = -1 * ((2 * b * b * cad.sx) / (a * a * cad.sy * cad.sy));
      C = ((b * b) / (cad.sy * cad.sy)) - 1;

      X1 = (-B + sqrt((B * B) - (4 * A * C))) / (2 * A);
      X2 = (-B - sqrt((B * B) - (4 * A * C))) / (2 * A);

      Y1 = ( (b*b) * (1 - ((X1 * cad.sx)/(a*a))) ) / cad.sy;
      Y2 = ( (b*b) * (1 - ((X2 * cad.sx)/(a*a))) ) / cad.sy;


      /* -----------------------------------------------------
       * 最後に、交点を元の座標に戻す
       */
      rtn->type = 2;
      rtn->ex[1] = X1;
      rtn->ey[1] = Y1;
      rtn->ex[2] = X2;
      rtn->ey[2] = Y2;

      rotation(&rtn->ex[1], &rtn->ey[1], 0, 0, get_ellipse_angle(ellipse));
      parallel(&rtn->ex[1], &rtn->ey[1], ellipse.cx, ellipse.cy);
      rotation(&rtn->ex[2], &rtn->ey[2], 0, 0, get_ellipse_angle(ellipse));
      parallel(&rtn->ex[2], &rtn->ey[2], ellipse.cx, ellipse.cy);

      return 2;
}





/* -------------------------------------------------------------------
 * 楕円と線の交点
 * lep
 * RtnDat
 *   線   : (sx[1],sy[1]) - (ex[1],ey[1])
 *   楕円 : 
 * RtnDat
 *   type   0 , 1 , 2
 *   点 1  (sx[2],sy[2])
 *   点 2  (sx[3],sy[3])
 */
int lep(struct RtnDat *rtn, ELLIPSE ellipse, int dgt)
{
//#define LEP
      ELLIPSE eps;
      CAD cad;
      double a, b, m, n, angle;
      double A, B, C, D, X1, Y1, X2, Y2;
      int debug = 0;


#ifdef LEP
      debug = 1;
#endif


      if (debug > 0) g_print("ellipse : (%f,%f),(%f,%f),%f, %f-%f\n", ellipse.cx, ellipse.cy, 
                                                                                                      ellipse.dx, ellipse.dy, 
                                                                                                      ellipse.k, 
                                                                                                      ellipse.sa, ellipse.ea);
      if (debug > 0) g_print("rtn : (%f,%f)-(%f,%f)\n", rtn->sx[1], rtn->sy[1], rtn->ex[1], rtn->ey[1]);





      /* -----------------------------------------------------
       * まず、楕円の中心が (0,0) で、角度が 0 になるような座標に変換。
       */
      /* 点 */
      init_cad(&cad);
      cad.sx = rtn->sx[1];
      cad.sy = rtn->sy[1];
      cad.ex = rtn->ex[1];
      cad.ey = rtn->ey[1];

      parallel(&cad.sx, &cad.sy, -ellipse.cx, -ellipse.cy);
      rotation(&cad.sx, &cad.sy, 0, 0, (-1 * get_ellipse_angle(ellipse)));
      parallel(&cad.ex, &cad.ey, -ellipse.cx, -ellipse.cy);
      rotation(&cad.ex, &cad.ey, 0, 0, (-1 * get_ellipse_angle(ellipse)));



      eps = ellipse;
      a = get_ellipse_a(eps);
      b = get_ellipse_b(eps);


      /* y = mx + n  */

      /* 垂直 */
      if (cad.ex > cad.sx - (1/pow(10, dgt)) && cad.ex < cad.sx + (1/pow(10, dgt))) {
            if (debug > 0) g_print("垂直\n");

            /* -a で接する */
            if (cad.sx > -a - (1/pow(10, dgt)) && cad.sx < -a + (1/pow(10, dgt))) {
                  if (debug > 0) g_print("-a で接する\n");
                  rtn->type = 1;
                  X1 = -a;
                  Y1 = 0;
            }

            /* a で接する */
            else if (cad.sx > a - (1/pow(10, dgt)) && cad.sx < a + (1/pow(10, dgt))) {
                  if (debug > 0) g_print("a で接する\n");
                  rtn->type = 1;
                  X1 = a;
                  Y1 = 0;
            }

            /* 重なる */
            else if (cad.sx > -a && cad.sx < a) {
                  if (debug > 0) g_print("重なる\n");
                  rtn->type = 2;

                  angle = acos(cad.sx / a);
                  X1 = cad.sx;
                  Y1 = b * sin(angle);
                  X2 = cad.sx;
                  Y2 = -Y1;
            }

            else {
                  if (debug > 0) g_print("交点無し\n");
                  rtn->type = 0;
                  return 0;
            }
      }



      /* 水平 */
      else if (cad.ey > cad.sy - (1/pow(10, dgt)) && cad.ey < cad.sy + (1/pow(10, dgt))) {
            if (debug > 0) g_print("水平\n");

            /* -b で接する */
            if (cad.sy > -b - (1/pow(10, dgt)) && cad.sy < -b + (1/pow(10, dgt))) {
                  rtn->type = 1;
                  X1 = 0;
                  Y1 = -b;
            }

            /* b で接する */
            else if (cad.sy > b - (1/pow(10, dgt)) && cad.sy < b + (1/pow(10, dgt))) {
                  rtn->type = 1;
                  X1 = 0;
                  Y1 = b;
            }

            /* 重なる */
            else if (cad.sy > -b && cad.sy < b) {
                  rtn->type = 2;
                  
                  angle = asin(cad.sy / b);
                  X1 = a * cos(angle);
                  Y1 = cad.sy;
                  X2 = -X1;
                  Y2 = Y1;
            }
            
            else {
                  rtn->type = 0;
                  return 0;
            }
      }



      /* 普通の線 */
      else {
            m = (cad.ey - cad.sy) / (cad.ex - cad.sx);
            n = cad.sy - (m * cad.sx);

            A = (1/(a*a)) + ((m*m)/(b*b));
            B = (2*m*n) / (b*b);
            C = ((n*n) / (b*b)) - 1;

            D = (B*B) - (4*A*C);


            /* 線と楕円は接している(交点が1個) */
            if (D < (1/pow(10, dgt)) && D > -(1/pow(10, dgt)) ) {
                  if (debug > 0) g_print("接している D = %f\n", D);
                  rtn->type = 1;

                  X1 = (-B / (2*A));
                  Y1 = m * X1 + n;
            }

            /* 重なっている(交点が2個) */
            else if (D > 0) {
                  if (debug > 0) g_print("重なっている D = %f\n", D);
                  rtn->type = 2;

                  X1 = (-B-sqrt(B*B-4*A*C)) / (2*A);
                  Y1 = m * X1 + n;
                  X2 = (-B+sqrt(B*B-4*A*C)) / (2*A);
                  Y2 = m * X2 + n;
            }

            /* 離れている(交点が0個) */
            else if (D < 0) {
                  if (debug > 0) g_print("離れている D = %f\n", D);
                  rtn->type = 0;
                  return 0;
            }

            else {
                  if (debug > 0) g_print("アレ D = %f\n", D);
                  return -1;
            }
      }





      if (rtn->type == 1) {
            rotation(&X1, &Y1, 0, 0, get_ellipse_angle(ellipse));
            parallel(&X1, &Y1, ellipse.cx, ellipse.cy);
            rtn->sx[2] = X1;
            rtn->sy[2] = Y1;
      }

      if (rtn->type == 2) {
            rotation(&X1, &Y1, 0, 0, get_ellipse_angle(ellipse));
            parallel(&X1, &Y1, ellipse.cx, ellipse.cy);
            rotation(&X2, &Y2, 0, 0, get_ellipse_angle(ellipse));
            parallel(&X2, &Y2, ellipse.cx, ellipse.cy);
            rtn->sx[2] = X1;
            rtn->sy[2] = Y1;
            rtn->sx[3] = X2;
            rtn->sy[3] = Y2;
      }

      return rtn->type;
}





/* -------------------------------------------------------------------
 * 楕円上の点(angle)での接線を求める
 * 
 */
int ellipse_tangent(struct RtnDat *rtn, ELLIPSE ellipse, double angle)
{
      double a, b;


      /* 楕円上の点を angle から求める */
      /* 中心 (0,0) で傾き 0 度の楕円の座標を求める */
      a = get_ellipse_a(ellipse);
      b = get_ellipse_b(ellipse);

      if (angle > -0.0001 && angle < 0.0001) {
            rtn->sx[1] = a;
            rtn->sy[1] = 0.;
            rtn->ex[1] = a;
            rtn->ey[1] = 100.;
            //g_print("ellipse_tangent() : rtn->ey[1] = %f\n", rtn->ey[1]);
      }
      else if (angle > 180 - 0.0001 && angle < 180 + 0.0001) {
            rtn->sx[1] = -a;
            rtn->sy[1] = 0.;
            rtn->ex[1] = -a;
            rtn->ey[1] = -100.;
            //g_print("ellipse_tangent() : rtn->ey[1] = %f\n", rtn->ey[1]);
      }
      else {
            rtn->sx[1] = get_ellipse_a(ellipse) * cos(degrad(angle));
            rtn->sy[1] = get_ellipse_b(ellipse) * sin(degrad(angle));

            rtn->ex[1] = rtn->sx[1] + 100;
            rtn->ey[1] = (1 - (rtn->sx[1] * rtn->ex[1]) / (a*a)) * ((b*b) / rtn->sy[1]);
      }


      /* -----------------------------------------------------
       * 最後に、交点を元の座標に戻す
       */
      rtn->type = 1;
      if (get_ellipse_angle(ellipse) != 0) {
            rotation(&rtn->sx[1], &rtn->sy[1], 0, 0, get_ellipse_angle(ellipse));
      }
      parallel(&rtn->sx[1], &rtn->sy[1], ellipse.cx, ellipse.cy);
      if (get_ellipse_angle(ellipse) != 0) {
            rotation(&rtn->ex[1], &rtn->ey[1], 0, 0, get_ellipse_angle(ellipse));
      }
      parallel(&rtn->ex[1], &rtn->ey[1], ellipse.cx, ellipse.cy);

      return 1;
}





/* -------------------------------------------------------------------
 * cc_check : 2円の位置関係を調べる
 * 
 *     : rtndat
 *     :     円 1      (cx[1] , cy[1]) ,r[1]
 *     :     円 2      (cx[2] , cy[2]) ,r[2]
 * 
 * return = type : 0 エラー
 *                 1 2円が離れている
 *                 2 2円が接している
 *                 3 2円が重なっている
 *                 4 2円が重なって接している
 *                 5 2円が重なって離れている
 * 
 */
int cc_check(struct RtnDat *rtn)
{
//#define CC_CHECK
      struct RtnDat pph;
      double dist,r1pr2,r1mr2;
      int debug = 0;


#ifdef CC_CHECK
      debug = 1;
#endif

      pph.sx[1] = rtn->cx[1];
      pph.sy[1] = rtn->cy[1];
      pph.ex[1] = rtn->cx[2];
      pph.ey[1] = rtn->cy[2];
      pp(&pph);

      dist = sg(pph.l, compa_digits);
      r1pr2 = sg((rtn->r[1] + rtn->r[2]), compa_digits);
      r1mr2 = sg((rtn->r[1] - rtn->r[2]), compa_digits);


      if (debug > 1) {
            /* 2円が離れている。接線4本 */
            g_print("type 4 : dist > r1pr2 = %d\n", dist > r1pr2);
            /* 2円が接している。接線3本 */
            g_print("type 3 : dist == r1pr2 = %d\n", dist == r1pr2);
            /* 2円が重なっている。接線2本 */
            g_print("type 2 : dist < r1pr2 && dist > r1mr2 = %d\n", (dist < r1pr2 && dist > r1mr2));
            g_print("  type 2-1 : dist < r1pr2 = %d\n", (dist < r1pr2));
            g_print("  type 2-2 : dist > r1mr2 = %d\n", (dist > r1mr2));
            /* 2円が重なって接している。接線1本 */
            g_print("type 1 : dist == r1mr2 = %d\n", (dist == r1mr2));
            /* 2円が重なって離れている。接線 0 本 */
            g_print("type 0 : (dist < r1mr2 && sg(rtn->l, compa_digits) > 0) = %d\n", 
                        (dist < r1mr2 && sg(rtn->l, compa_digits) > 0));
            g_print("  type 0-1 : (dist < r1mr2) = %d\n", (dist < r1mr2));
            g_print("  type 0-2 : (sg(rtn->l, compa_digits) > 0) = %d\n", (sg(rtn->l, compa_digits) > 0));
      }


      /* 2円が離れている。接線4本 */
      if (dist > r1pr2) {
            if (debug > 0) g_print("type 1 : 2円が離れている\n");
            rtn->type = 1;
            return 1;
      }

      /* 2円が接している。接線3本 */
      else if (dist == r1pr2) {
            if (debug > 0) g_print("type 2 : 2円が接している\n");
            rtn->type = 2;
            return 2;
      }

      /* 2円が重なっている。接線2本 */
      else if (dist < r1pr2 && dist > r1mr2) {
            if (debug > 0) g_print("type 2 : 2円が重なっている\n");
            rtn->type = 3;
            return 3;
      }

      /* 2円が重なって接している。接線1本 */
      else if (dist == r1mr2) {
            if (debug > 0) g_print("type 4 : 2円が重なって接している\n");
            rtn->type = 4;
            return 4;
      }

      /* 2円が重なって離れている。接線 0 本 */
      else if (dist < r1mr2 && sg(rtn->l, compa_digits) > 0) {
            if (debug > 0) g_print("type 5 : 2円が重なって離れている\n");
            rtn->type = 5;
            return 5;
      }

      return 0;
}





/* -------------------------------------------------------
 * 楕円と楕円(円)の交点 (MAIN)
 * 
 * 楕円は4個もありえる
 * 
 * rtn->type : 0 交点無し
 *             1 交点1個  (rtn->sx[1],rtn->sy[1])
 *             2 交点2個  (rtn->sx[2],rtn->sy[2])
 *             3 交点3個  (rtn->sx[3],rtn->sy[3])
 *             4 交点4個  (rtn->sx[4],rtn->sy[4])
 * 
 *          楕円の位置関係も調べたい。
 *  return : 0 : エラー
 *           1 : 離れている
 *           2 : 並んで接している
 *           3 : 重なっている
 *           4 : 重なって接している
 *           5 : 重なって離れている
 */
int ellipse_on_ellipse(struct RtnDat *rtn, ELLIPSE ellipseA, ELLIPSE ellipseB)
{
//#define ELLIPSE_ON_ELLIPSE
      int i, j, Ret, Ret1, Ret2;
      double value;
      int A_data[360];
      int B_data[360];
      int AinB = 0, BinA = 0, start_frag = 0;
      double sa[5], ea[5];
      double sa1, ea1, sa2, ea2;
      SAG_POINT point1, point2;
      int BeforeRet = 0;
      int debug = 0;


#ifdef ELLIPSE_ON_ELLIPSE
      debug = 1;
#endif

      rtn->type = 0;
      if (debug > 0) g_print ("ellipse_on_ellipse() : in\n");


      /* -----------------------------------------------------
       * データの収集
       * 
       */
      /* 楕円Aを1度ごとにデータを収集 */
      value = (2 * PI * get_ellipse_a(ellipseA)) / (360*2);
      if (debug > 0) g_print ("楕円Aの解析\n");
      for (i = 0 ; i < 360 ; i++) {
            get_ellipse_point_from_angle (ellipseA, (double)i, &point1);
            Ret = point_in_ellipse (&ellipseB, point1.x, point1.y, value);
            if (Ret == 3) Ret = 1;
            if (debug > 0) {
                  g_print ("%d", Ret);
                  if (i == 89 || i == 179 || i == 269 || i == 359) g_print("\n");
            }
            A_data[i] = Ret;
      }

      /* 楕円Bを1度ごとにデータを収集 */
      value = (2 * PI * get_ellipse_a(ellipseB)) / (360*2);
      if (debug > 0) g_print ("楕円Bの解析\n");
      for (i = 0 ; i < 360 ; i++) {
            get_ellipse_point_from_angle (ellipseB, (double)i, &point1);
            Ret = point_in_ellipse (&ellipseA, point1.x, point1.y, value);
            if (Ret == 3) Ret = 1;
            if (debug > 0) {
                  g_print ("%d", Ret);
                  if (i == 89 || i == 179 || i == 269 || i == 359) g_print("\n");
            }
            B_data[i] = Ret;
      }


      /* 楕円Bの中心は楕円Aの中 */
      Ret = point_in_ellipse (&ellipseA, ellipseB.cx, ellipseB.cy, 0.0000001);
      if (Ret == 1) BinA = 1;
      /* 楕円Aの中心は楕円Bの中 */
      Ret = point_in_ellipse (&ellipseB, ellipseA.cx, ellipseA.cy, 0.0000001);
      if (Ret == 1) AinB = 1;





      /* -----------------------------------------------------
       * SA(out->in) or EA(in->out) 情報を解析 
       * 
       */
      for (i = 0 ; i < 5 ; i++) {
            sa[i] = 0;
            ea[i] = 0;
      }
      if (A_data[0] == 1) {
            start_frag = 1;
      }
      j = 0;
      BeforeRet = 0;
      for (i = 0 ; i < 360 ; i++) {
            /* out */
            if (A_data[i] == 2) {
                  if (BeforeRet == 1) {
                        ea[j] = (double) i;
                        j++;
                  }
                  BeforeRet = 2;
            }
            /* in or on */
            else if (A_data[i] == 1) {
                  if (BeforeRet == 2) {
                        sa[j] = (double) i-1;
                  }
                  BeforeRet = 1;
            }
      }
      if (start_frag == 1) {
            sa[0] = sa[j];
      }

      if (debug > 0) {
            for (i = 0 ; i < j ; i++) {
                  g_print ("sa[%d] = %f  ea[%d] = %f\n", i, sa[i], i, ea[i]);
            }
      }





      /* -----------------------------------------------------
       * SA(out->in) or EA(in->out) 情報から交点を探索 
       * 
       */
      for (i = 0 ; i < j ; i++) {
            /*  */
            if (sa[i] == 0 && ea[i] == 0) {
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa[i] == 0 && ea[i] == 0 で終了\n");
                  return 0;
            }

            /* SA < EA */
            else if (sa[i] < ea[i]) {
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa[i] < ea[i]\n");

                  /* 探索範囲を初期設定 */
                  sa1 = sa[i];
                  ea1 = ((ea[i] - sa[i]) / 2) + sa[i];
                  sa2 = ((ea[i] - sa[i]) / 2) + sa[i];
                  ea2 = ea[i];

                  /* 領域の最初 */
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa1 = %f  ea1 = %f\n", sa1, ea1);
                  Ret1 = ellipse_on_ellipse_point(ellipseA, ellipseB, &sa1, &ea1, ((ea[i] - sa[i]) / 10), 0.0000001);
                  if (Ret1 == 1) {
                        rtn->type++;
                        get_ellipse_point_from_angle (ellipseA, ((ea1 + sa1) / 2), &point1);
                        rtn->sx[rtn->type] = point1.x;
                        rtn->sy[rtn->type] = point1.y;
                        if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : (%f,%f)\n", point1.x, point1.y);
                  }

                  /* 領域の最後 */
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa2 = %f  ea2 = %f\n", sa2, ea2);
                  Ret2 = ellipse_on_ellipse_point(ellipseA, ellipseB, &sa2, &ea2, ((ea[i] - sa[i]) / 10), 0.0000001);
                  if (Ret2 == 1) {
                        rtn->type++;
                        get_ellipse_point_from_angle (ellipseA, ((ea2 + sa2) / 2), &point2);
                        rtn->sx[rtn->type] = point2.x;
                        rtn->sy[rtn->type] = point2.y;
                        if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : (%f,%f)\n", point2.x, point2.y);
                  }

                  if (Ret1 == 1 && Ret2 == 1) {
                        if (sg(point1.x, 4) == sg(point2.x, 4) && sg(point1.y, 4) == sg(point2.y, 4)) {
                              rtn->type--;
                              rtn->sx[rtn->type] = sg(point1.x, 4);
                              rtn->sy[rtn->type] = sg(point1.y, 4);
                        }
                  }
            }



            /* EA < SA */
            else if (sa[i] > ea[i]) {
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa[i] > ea[i]\n");
                  /* 探索範囲を初期設定 */
                  sa1 = sa[i];
                  ea1 = sa[i] + (360-sa[i]+ea[i])/2;
                  if (ea1 <= 0) ea1 = ea1 + 360;

                  sa2 = sa[i] + (360-sa[i]+ea[i])/2;
                  if (sa2 >= 360) sa2 = sa2 - 360;
                  if (sa2 == 360) sa2 = 0;
                  ea2 = ea[i];


                  /* 領域の最初 */
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa1 = %f  ea1 = %f\n", sa1, ea1);
                  Ret1 = ellipse_on_ellipse_point(ellipseA, ellipseB, &sa1, &ea1, ((360-sa[i]+ea[i])/20), 0.0000001);
                  if (Ret1 == 1) {
                        rtn->type++;
                        get_ellipse_point_from_angle (ellipseA, ((ea1 + sa1) / 2), &point1);
                        rtn->sx[rtn->type] = point1.x;
                        rtn->sy[rtn->type] = point1.y;
                        if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : (%f,%f)\n", point1.x, point1.y);
                  }

                  /* 領域の最後 */
                  if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : sa2 = %f  ea2 = %f\n", sa2, ea2);
                  Ret2 = ellipse_on_ellipse_point(ellipseA, ellipseB, &sa2, &ea2, ((360-sa[i]+ea[i])/20), 0.0000001);
                  if (Ret2 == 1) {
                        rtn->type++;
                        get_ellipse_point_from_angle (ellipseA, ((ea2 + sa2) / 2), &point2);
                        rtn->sx[rtn->type] = point2.x;
                        rtn->sy[rtn->type] = point2.y;
                        if (debug > 0) g_print("Ellipse.c : ellipse_on_ellipse_rect() : (%f,%f)\n", point2.x, point2.y);
                  }

                  if (Ret1 == 1 && Ret2 == 1) {
                        if (sg(point1.x, 4) == sg(point2.x, 4) && sg(point1.y, 4) == sg(point2.y, 4)) {
                              rtn->type--;
                              rtn->sx[rtn->type] = sg(point1.x, 4);
                              rtn->sy[rtn->type] = sg(point1.y, 4);
                        }
                  }
            }
      }

      /* 離れている */
      if (rtn->type == 0 && BinA == 0 && AinB == 0) {
            if (debug > 0) g_print ("ellipse_on_ellipse() : out 1\n");
            return 1;
      }
      /* 並んで接している */
      else if (rtn->type == 1 && BinA == 0 && AinB == 0) {
            if (debug > 0) g_print ("ellipse_on_ellipse() : out 2\n");
            return 2;
      }
      /* 重なっている */
      else if (rtn->type > 1) {
            if (debug > 0) g_print ("ellipse_on_ellipse() : out 3\n");
            return 3;
      }
      /* 重なって接している */
      else if (rtn->type == 1 && (BinA == 1 || AinB == 1)) {
            if (debug > 0) g_print ("ellipse_on_ellipse() : out 4\n");
            return 4;
      }
      /* 重なって離れている */
      else if (rtn->type == 0 && (BinA == 1 || AinB == 1)) {
            if (debug > 0) g_print ("ellipse_on_ellipse() : out 5\n");
            return 5;
      }

      if (debug > 0) g_print ("ellipse_on_ellipse() : out 0\n");
      return 0;
}





/* -------------------------------------------------------------------
 * 楕円Aと楕円Bの交点となる、楕円Aの角度の範囲 
 * 
 */
int ellipse_on_ellipse_point(ELLIPSE ellipseA, 
                                           ELLIPSE ellipseB, 
                                           double *sa, 
                                           double *ea, 
                                           double angle_degits, 
                                           double last_degits 
                                           )
{
      double i, j = 0;
      int Ret = 0, BeforeRet = 0, End = 0;
      double dgt, value;
      double min = 0, max = 0;
      SAG_POINT point;
      int debug = 0;


//    value = (2 * PI * get_ellipse_a(ellipseA) * angle_degits) / (360*2);
      value = last_degits/5000;
      if (debug > 0) g_print("ellipse_on_ellipse_point() : sa = %f   ea = %f   value = %e   angle_degits = %e   last_degits = %e\n", *sa, *ea, value, angle_degits, last_degits);




      /* SA < EA */
      if (*sa < *ea) {
            if (debug > 0) g_print("ellipse_on_ellipse_point() : sa < ea\n");

            min = *sa;
            max = *ea;

            for (i = *sa - angle_degits ; i < *ea + angle_degits ; i = i + angle_degits) {
                  if (End != 1) {
                        get_ellipse_point_from_angle (ellipseA, (double)i, &point);
                        Ret = point_in_ellipse (&ellipseB, point.x, point.y, value);
                        if (Ret == 3) Ret = 1;

                        if (debug > 0) g_print("ellipse_on_ellipse_point() : i = %f   Ret = %d\n", i, Ret);
                        j++;

                        /* out */
                        if (Ret == 2) {
                              /* start in */
                              if (BeforeRet == 0) {
                                    //min = (double) i - angle_degits;
                              }
                              /* out -> out */
                              else if(BeforeRet == 2) {
                                    min = (double) i;
                              }
                              /* in -> out */
                              else if (BeforeRet == 1) {
                                    max = (double) i;
                                    End = 1;
                              }
                              BeforeRet = 2;
                        }
                        /* in */
                        else if (Ret == 1) {
                              /* start in */
                              if (BeforeRet == 0) {
                                    //min = (double) i - angle_degits;
                              }
                              /* in -> in */
                              else if(BeforeRet == 1) {
                                    min = (double) i;
                              }
                              /* out -> in */
                              else if (BeforeRet == 2) {
                                    max = (double) i;
                                    End = 1;
                              }
                              BeforeRet = 1;
                        }
                  }
            }

            *sa = min;
            *ea = max;
            dgt = (max - min);

            if (debug > 0) g_print("ellipse_on_ellipse_point() : j = %f   X10000 : max[%f] - min[%f] = dgt[%f]\n\n\n", j, 10000*max, 10000*min, 10000*dgt);

            if (dgt < last_degits) {
//                if (debug > 0) g_print("dgt < last_degits\n");
                  return 1;  // Ans (max + min) / 2
            }
            else {
//                if (debug > 0) g_print("dgt > last_degits\n");
                  return ellipse_on_ellipse_point(ellipseA, ellipseB, sa, ea, (dgt / 10), last_degits);
            }
      }





      /* EA < SA */
      else if (*ea < *sa) {
            if (debug > 0) g_print("ellipse_on_ellipse_point() : sa > ea\n");

            min = *sa;
            max = *ea;

            for (i = *ea - angle_degits ; i < 360 + angle_degits ; i = i + angle_degits) {
                  if (End != 1) {
                        get_ellipse_point_from_angle (ellipseA, (double)i, &point);
                        Ret = point_in_ellipse (&ellipseB, point.x, point.y, value);
                        if (Ret == 3) Ret = 1;
                        
                        if (debug > 0) g_print("ellipse_on_ellipse_point() : i = %f   Ret = %d\n", i, Ret);
                        j++;
      //                if (debug > 0) g_print("%d", Ret);

                        /* out */
                        if (Ret == 2) {
                              /* start in */
                              if (BeforeRet == 0) {
                                    //min = (double) i - angle_degits;
                              }
                              /* out -> out */
                              else if(BeforeRet == 2) {
                                    min = (double) i;
                              }
                              /* in -> out */
                              else if (BeforeRet == 1) {
                                    max = (double) i;
                                    End = 1;
                              }
                              BeforeRet = 2;
                        }
                        /* in */
                        else if (Ret == 1) {
                              /* start in */
                              if (BeforeRet == 0) {
                                    //min = (double) i - angle_degits;
                              }
                              /* in -> in */
                              else if(BeforeRet == 1) {
                                    min = (double) i;
                              }
                              /* out -> in */
                              else if (BeforeRet == 2) {
                                    max = (double) i;
                                    End = 1;
                              }
                              BeforeRet = 1;
                        }
                  }
            }

            //if (max == *sa) {
                  for (i = 0 - angle_degits ; i < *sa + angle_degits + angle_degits ; i = i + angle_degits) {
                        if (End != 1) {
                              get_ellipse_point_from_angle (ellipseA, (double)i, &point);
                              Ret = point_in_ellipse (&ellipseB, point.x, point.y, value);
                              if (Ret == 3) Ret = 1;

                              if (debug > 0) g_print("ellipse_on_ellipse_point() : i = %f   Ret = %d\n", i, Ret);
                              j++;
//                            if (debug > 0) g_print("%d", Ret);

                              /* out */
                              if (Ret == 2) {
                                    /* start in */
                                    if (BeforeRet == 0) {
                                          //min = (double) i - angle_degits;
                                    }
                                    /* out -> out */
                                    else if(BeforeRet == 2) {
                                          min = (double) i;
                                    }
                                    /* in -> out */
                                    else if (BeforeRet == 1) {
                                          max = (double) i;
                                          End = 1;
                                    }
                                    BeforeRet = 2;
                              }
                              /* in */
                              else if (Ret == 1) {
                                    /* start in */
                                    if (BeforeRet == 0) {
                                          //min = (double) i - angle_degits;
                                    }
                                    /* in -> in */
                                    else if(BeforeRet == 1) {
                                          min = (double) i;
                                    }
                                    /* out -> in */
                                    else if (BeforeRet == 2) {
                                          max = (double) i;
                                          End = 1;
                                    }
                                    BeforeRet = 1;
                              }
                        }
                  }
            //}
            


            *sa = min;
            *ea = max;
            if (min < max) {
                  dgt = (max - min);
            }
            else if (min > max) {
                  dgt = ((360-min) + max);
            }
            else dgt = 0;
            
            //dgt = max - min;
            if (debug > 0) g_print("ellipse_on_ellipse_point() : j = %f   X10000 : max[%f] - min[%f] = dgt[%f]\n\n\n", j, 10000*max, 10000*min, 10000*dgt);

            if (dgt < last_degits) {
                  return 1;  // Ans (max + min) / 2
            }
            else {
                  return ellipse_on_ellipse_point(ellipseA, ellipseB, sa, ea, dgt / 10, last_degits);
            }
      }
      return 0;
}





/* -------------------------------------------------------
 * 楕円から楕円(円)の接線 (Main)
 * 
 * 
 * return : 0 : 接線なし
 *          1 : 
 *          2 : 
 *          3 : 
 *          4 : 
 * 
 */
int ellipse_to_ellipse(struct RtnDat *rtn, ELLIPSE ellipseA, ELLIPSE ellipseB)
{
      int Ret, i, j = 0, BeforeRet = 0;
      int A_data[360];
      double sp[5], ep[5], AnsAngle;
      struct RtnDat lo_rtn;
      SAG_POINT point;
      int debug = 0;


      /* -----------------------------------------------------
       * データの収集
       * 
       * 楕円Aの接線を 1 度ごとに求め、その接線と
       * 楕円Bの位置関係を A_data[] の配列に入れる。
       */
      /* 楕円Aを1度ごとにデータを収集 */
      if (debug > 0) g_print ("楕円Aの解析\n");
      for (i = 0 ; i < 360 ; i++) {
            ellipse_tangent(rtn, ellipseA, (double)i);
            lep(rtn, ellipseB, 4);
            Ret = rtn->type;
            if (Ret == 1) Ret = 2;
            if (debug > 0) {
                  g_print ("%d", Ret);
                  if (i == 89 || i == 179 || i == 269 || i == 359) g_print("\n");
            }
            A_data[i] = Ret;
      }



      /* -----------------------------------------------------
       * 情報を解析 
       * 
       */
      /* 接線があると思われる角度の範囲を格納する配列を初期化 */
      for (i = 0 ; i < 5 ; i++) {
            sp[i] = 0;
            ep[i] = 0;
      }

      /* 接線があると思われる角度の範囲を解析 */
      j = 0;
      BeforeRet = 1;
      for (i = 0 ; i < 360 ; i++) {
            /* out -> in */
            if (A_data[i] == 2) {
                  if (BeforeRet == 0) {
                        if (debug > 0) g_print ("i = %d  sp = %d  ep = %d\n", i, i-2, i);
                        sp[j] = (double)(i - 2.);
                        ep[j] = (double)i;
                        j++;
                  }
                  BeforeRet = 2;
            }
            /* in -> out */
            else if (A_data[i] == 0) {
                  if (BeforeRet == 2) {
                        if (debug > 0) g_print ("i = %d  sp = %d  ep = %d\n", i, i-2, i);
                        sp[j] = (double)(i - 2.);
                        ep[j] = (double)i;
                        j++;
                  }
                  BeforeRet = 0;
            }
      }

      /* 接線があると思われる角度の範囲のデータを確認 */
      if (debug > 0) {
            for (i = 0 ; i < j ; i++) {
                  g_print ("sp[%d] = %f  ep[%d] = %f\n", i, sp[i], i, ep[i]);
            }
      }



      /* -----------------------------------------------------
       * 接線を求める
       * 
       * 接線があると思われる角度の範囲を格納している配列の各データから
       * さらに精度を細かくして接線を求め、交点を求める。
       */
      rtn->type = 0;
      for (i = 0 ; i < j ; i++) {
            point.frag = 0;
            Ret = ellipse_to_ellipse_search (ellipseA, ellipseB, &sp[i], &ep[i], &point, 0.1, 0.0000001);
            if (Ret == 1) {
                  rtn->type++;
                  AnsAngle = (sp[i] + ep[i]) / 2;
                  ellipse_tangent(&lo_rtn, ellipseA, AnsAngle);
                  rtn->sx[rtn->type] = lo_rtn.sx[1];
                  rtn->sy[rtn->type] = lo_rtn.sy[1];

                  lep(&lo_rtn, ellipseB, 10);
                  if (lo_rtn.type == 1) {
                        rtn->ex[rtn->type] = lo_rtn.sx[2];
                        rtn->ey[rtn->type] = lo_rtn.sy[2];
                  }
                  else if (lo_rtn.type == 2) {
                        rtn->ex[rtn->type] = (lo_rtn.sx[2] + lo_rtn.sx[3]) / 2;
                        rtn->ey[rtn->type] = (lo_rtn.sy[2] + lo_rtn.sy[3]) / 2;
                  }
                  else if (point.frag == 1) {
                        rtn->ex[rtn->type] = point.x;
                        rtn->ey[rtn->type] = point.y;
                  }
                  else {
                        rtn->type--;
                  }
            }
      }

      return rtn->type;
}





/* -------------------------------------------------------
 * 楕円Aの sp ep から楕円B(円)の接線を探す (再帰関数)
 * 
 * 
 * 楕円Aの sp から ep の角度の間をステップ angle_degits で接線を求め、
 * 
 * 楕円Bの中から外へ抜ける区間の sp ep の範囲をもっと狭く
 * または
 * 楕円Bの外から中へ入る区間の sp ep の範囲をもっと狭く
 * して、再帰呼び出しにてまた ellipse_to_ellipse_search() を呼ぶ。
 * 
 * 終了条件は、sp ep の範囲が last_degits よりも小さくなったときか、
 * sp ep の範囲が変わらないとき。
 * 
 */
int ellipse_to_ellipse_search (ELLIPSE ellipseA, 
                                             ELLIPSE ellipseB, 
                                             double *sp, 
                                             double *ep, 
                                             SAG_POINT *point, 
                                             double angle_degits, 
                                             double last_degits 
                                             )
{
      double i, dgt;
      int Ret = 0, BeforeRet = 1;
      struct RtnDat rtn;
      int debug = 0;
      double osp, oep;


      osp = *sp;
      oep = *ep;

      if (debug > 0) g_print ("sp = %f   ep = %f   angle_degits = %f\n", *sp, *ep, angle_degits);



      /* -----------------------------------------------------
       * 楕円Aの接線が、楕円Bの中から外、または外から中に
       * 切り替わる範囲を探す。
       */
      BeforeRet = 1;
      for (i = *sp ; i < *ep + angle_degits; i = i + angle_degits) {
            ellipse_tangent(&rtn, ellipseA, i);
            lep(&rtn, ellipseB, 10);
            if (rtn.type == 1) {
                  point->frag = 1;
                  point->x = rtn.sx[2];
                  point->y = rtn.sy[2];
                  if (debug > 0) g_print("point : Ret = %d   i = %f\n", Ret, i);
            }
            Ret = rtn.type;
            /* 交点が 1 個でも 2 個として考える */
            if (Ret == 1) Ret = 2;
            if (debug > 0) g_print ("i = %f   Ret = %d\n", i, Ret);

            /* out -> in */
            if (Ret == 2) {
                  if (BeforeRet == 0) {
                        if (debug > 0) g_print ("    i = %f  sp = %f  ep = %f\n", i, i - angle_degits, i + angle_degits);
                        *sp = i - angle_degits;
                        *ep = i + angle_degits;
                  }
                  BeforeRet = 2;
            }
            /* in -> out */
            else if (Ret == 0) {
                  if (BeforeRet == 2) {
                        if (debug > 0) g_print ("    i = %f  sp = %f  ep = %f\n", i, i - angle_degits, i + angle_degits);
                        *sp = i - angle_degits;
                        *ep = i + angle_degits;
                  }
                  BeforeRet = 0;
            }
      }



      /* -----------------------------------------------------
       * 求めた範囲によって次の処理をわける。
       */
      /* 誤差を求める */
      dgt = (*ep - *sp);
      if (debug > 0) g_print("    X10000 : ep[%f] - sp[%f] = dgt[%f]\n\n\n", 10000*(*ep), 10000*(*sp), 10000*dgt);

      /* sp ep が変化無しなら接線無しとして終了 */
      if (osp == *sp && oep == *ep) {
            if (debug > 0) g_print("osp == *sp && oep == *ep によってキャンセル\n");
            return 0;
      }

      /* 誤差が last_degits より小さければ終了 */
      else if (dgt < last_degits) {
            //if (debug > 0) g_print("dgt < last_degits\n");
            if (point->frag != 1) {
                  for (i = *sp ; i < *ep + (angle_degits / 10); i = i + (angle_degits / 10)) {
                        ellipse_tangent(&rtn, ellipseA, i);
                        lep(&rtn, ellipseB, 10);
                        if (rtn.type == 1) {
                              point->frag = 1;
                              point->x = rtn.sx[2];
                              point->y = rtn.sy[2];
                        }
                        else {
                              if (debug > 0) g_print("ellipse_to_ellipse_search() : 1 : の中では交点が出せなかった。\n");
                        }
                  }
                  if (point->frag != 1) {
                        ellipse_tangent(&rtn, ellipseA, (*ep + *sp) / 2);
                        lep(&rtn, ellipseB, 4);
                        if (rtn.type == 1) {
                              point->frag = 1;
                              point->x = rtn.sx[2];
                              point->y = rtn.sy[2];
                        }
                        else {
                              if (debug > 0) g_print("ellipse_to_ellipse_search() : 2 : の中では交点が出せなかった。\n");
                        }
                  }
            }
            return 1;  // Ans (max + min) / 2
      }

      /* 誤差が last_degits より大きければもっと細かく探すために ellipse_to_ellipse_search() を呼ぶ */
      else {
            //if (debug > 0) g_print("dgt > last_degits\n");
            return ellipse_to_ellipse_search(ellipseA, ellipseB, sp, ep, point, (dgt / 10), last_degits);
      }

      return 0;
}





/* ====================================================================
 * ===  Copyright (C) 1998-2007 Yutaka Sagiya. All rights reserved. ===
 * ====================================================================
 *    Project              : SagCAD
 *    Source               : Ellipse.c
 * ====================================================================
 */

Generated by  Doxygen 1.6.0   Back to index