コンピュータ、ソフトウェアに関すること。忘れないようにね。

アポロニウス【問題5】の解答

問題5の解答

【問題5】1つの円と2本の直線に接する円

幾何的に作図をしようと思ったのですが,うまい方法が見つからず,解析的に作図をすることにしました。その結果,ソースがかなりプログラムになってしまいました。わかりにくいかとは思いますが,ご容赦ください。

(1) 2本の直線が平行になる場合

(2) 2本の直線が平行にならない場合



(1) 2本の直線が平行になる場合


ソース

import geometry;
unitsize(1cm);

/**
 * 2直線が平行な場合
 *    この2直線の距離が求める接円の直径となる。
 *  そこで,接円の中心をx軸上とし,与えられる円の中心をy軸上とした。
 */

// 2直線と円を与える
real r = 3, R = 2, Y = 2.5;

circle cG = circle((0,Y), R);
line l1 = line((-R,r),(R,r));  // y=r
line l2 = line((-R,-r),(R,-r));// y=-r

// 接円の中心を求める
real x = sqrt((R+r)^2-Y^2);
circle c1 = circle((x,0),r);
circle c2 = circle((-x,0),r);

draw(l1);draw(l2);draw(cG);
draw(c1, red);
dot("$C_1$",c1.C);
draw(c2, red);
dot("$C_2$",c2.C);


(2) 2本の直線が平行にならない場合


ソース

import geometry;
size(0,10cm);

/**
 * 2直線が平行でない場合
 *    この2直線の交点を原点とする。
 *  接円の中心は、2等分線上にある。これをx軸とする。
 *  与えられる2直線は、y=mx , y=-mx で表される。
 *
 *    接円があまり大きくならないように注意しましょう。
 */

// 2直線と円を与える
/**
 *  円が2直線の上側にある
 *real m = 1, R = .5, X = 1, Y = 2;
 */
/**
 *  円が直線のひとつと接している
 *real m = 1, R = sqrt(2)/2, X = 2, Y = 1;
 */
/**
 *  円が直線のひとつと2点で交わる
 *real m = 1, R = 1, X = 2, Y = 1;
 */
/**
 *  円が2直線の右側にある
 *real m = 1, R = .5, X = 1, Y = 0;
 */
/**
 *  円が2直線でできる4つの領域にまたがる
 *real m = 1, R = 1, X = 0, Y =0;
 */

real m = 1, R = 1, X = 0, Y =0;

real sin = m/sqrt(1+m^2);
real cos = 1/sqrt(1+m^2);
// 2直線と円を描く
circle cG = circle((X,Y), R);
line l1 = line((0,0),(1,m));  // y=mx
line l2 = line((0,0),(1,-m)); // y=-mx
draw(l1);draw(l2);draw(cG);

// 判別式
real discriminant(real R, real X, real Y, real sin, real cos){
  real D = (X+R*sin)^2 - cos^2*(X^2+Y^2-R^2);

  return D;
}

/**
 *  接円の中心を求める
 */
real [] zcal(real X, real Y, real sin, real cos){
  real[] x={0,0,0,0};
  if(Y==0){
    x[0] = (X-R)/(1-sin);
    x[1] = (X+R)/(1-sin);
    x[2] = (X-R)/(1+sin);
    x[3] = (X+R)/(1+sin);
  }else{
    real D = discriminant(R, X, Y, sin, cos);
    if(D>=0){
      x[0] = (X+R*sin - sqrt(D))/(cos^2);
      x[1] = (X+R*sin + sqrt(D))/(cos^2);
    }
    real D = discriminant(-R, X, Y, sin, cos);
    if(D>=0){
      x[2] = (X-R*sin - sqrt(D))/(cos^2);
      x[3] = (X-R*sin + sqrt(D))/(cos^2);
    }
  }
  return x;
}
/**
 *  接円を描画する
 */
void drawC(real R, real X, real Y, real sin, real cos, int i=0){
  circle c1,c2;
  real[] z = zcal(X,Y,sin,cos);
  if(i==0){
    if(z[0]!=0){
      c1 = circle((z[0],0),z[0]*sin);
      draw(c1, red);
      dot(format("$C_{%i1}$",i),c1.C);
      if(z[1]!=z[0]){
	c2 = circle((z[1],0),z[1]*sin);
	draw(c2, red);
	dot(format("$C_{%i2}$",i),c2.C);
      }
    }
    if(z[2]!=0){
      c1 = circle((z[2],0),z[2]*sin);
      draw(c1, red);
      dot(format("$C_{%i3}$",i),c1.C);
      if(z[3]!=z[2]){
	c2 = circle((z[3],0),z[3]*sin);
	draw(c2, red);
	dot(format("$C_{%i4}$",i),c2.C);
      }
    }
  }
  z = zcal(Y,X,cos,sin);
  if(z[0]!=0){
    i=1;
    c1 = circle((0,z[0]),z[0]*cos);
    draw(c1, blue);
    dot(format("$C_{%i1}$",i),c1.C);
    if(z[1]!=z[0]){
      c2 = circle((0,z[1]),z[1]*cos);
      draw(c2, blue);
      dot(format("$C_{%i2}$",i),c2.C);
    }
  }
  if(z[2]!=0){
    i=1;
    c1 = circle((0,z[2]),z[2]*cos);
    draw(c1, blue);
    dot(format("$C_{%i3}$",i),c1.C);
    if(z[3]!=z[2]){
      c2 = circle((0,z[3]),z[3]*cos);
      draw(c2, blue);
      dot(format("$C_{%i4}$",i),c2.C);
    }
  }
  
}

real D = discriminant(R, X, Y, sin, cos);
if(D >= 0){
  drawC(R, X, Y, sin, cos);
}else if(D < 0){
  drawC(R, X, Y, sin, cos, 1);
}