4. 「とんでる力学」9 月分

4.1. プログラムと解説

赤木: というわけで、今月もプログラムとかの解説ね。

学生S: えーと、あれ、今月数値積分とかしてないですよ。図2は解析的な計算でグラフ書いただけだし、後はお絵かきツールで書いてますから。

赤木: あ、そうだったわね。折角始めたから、図2のプログラム 出して。

学生S: これですけど、、、

  include Math
  printer september-fig2.ps/vcps
  term
  def orbit_plot(x0,kappa)
    reloc 1+x0 0
    dtheta= PI/200;
    theta=dtheta;
    while theta < 2*PI+dtheta/2
      p theta
      x=cos(theta)
      y=sin(theta)
      epx = x0*cos(kappa*theta)
      epy = -2*x0/kappa*sin(kappa*theta)

      xp = x + x*epx-y*epy
      yp = y + epx*y*epy*x

      draw xp yp
      print(theta," ", xp," ",  yp," ",  x," ",  y," ",  epx," ", epy)
      theta += dtheta
    end
  end

  square
  expand 1.5
  lw 2
  limit -1.2 1.2 -1.2 1.2
  box

  lw 2
  lt 1
  x0=0.1
  orbit_plot(x0,1)
  lt 2
  orbit_plot(x0,2)
  lt 3
  orbit_plot(x0,4)
  pgend
  psfix

赤木: なんか力任せね。でも、確かにあんまり解説という話ではないわね。角度 を変えて座標を計算して線でつなぐだけだものね。 x, y がベースの円運動、 epx, epy が周転円上の座標ね。

学生S: そうです。足す時には周転円上の座標のほうを回転するので、行列を掛ける格好になってます。

赤木: まあ、これはそんなものね、、、

学生S: はい。というわけで今月はこれくらいで勘弁して下さい。来月のもあるし。

赤木: 来月のは、でも、来月書くのよね？

学生S: そうですけど、なんか作者が忙しいんですよ、ここんとこ。我々は所詮作者の分身というか創作物というかそういうものですから、作者に時間がないと我々に与えられる時間もなくなるわけで。