Previous ToC Next

6. 「とんでる力学」2004年 11月分

6.1. プログラムと解説

赤木: さて、問題の 11月分よね。

学生S: です。えーと、その、今回スキップ、というわけにいきません?

赤木: いかないわよね。てゆーか、なんで?

学生S: つまりですね、本文のほうみてもらえばわからないように、こ の回で使ったプログラム、作者が大昔に書いたものなんですよね、、、

赤木: えーと、まあ、そこは後でいいから君が書いたところは?

学生S: えーと、、、これですね


  # plot.rb : postprocessing of three body integration

  require 'vector.rb'
  include Math
  class Body
    attr_accessor :mass, :pos, :vel, :acc, :pot, :withring
    def initialize
      @pos=[]
      @vel=[]
      @withring=nil
    end
    def froms(s)
      a=s.split.collect{|x| x.to_f}
      p a
      n=(a.size-1)/2
      @mass = a[0]
      @pos  = a[1..n].to_v
      @vel  = a[(n+1)..(2*n)].to_v
      p self
    end

    def convertfromelements(a,e)
      @pos[0]=a*(1+e)
      @pos[1]=0
      @vel[0]=0
      @vel[1]= sqrt(2*(1/@pos[0]-0.5/a))
      @pos=@pos.to_v
      @vel=@vel.to_v
    end
    def fromvectortocontactelements(mass,pos,vel)
      r=sqrt(pos*pos)
      be = -mass/r+0.5*(vel*vel)
      # be = -mass/2a, a = -mass/2be
      a = -mass/(be*2)
      vcirc = sqrt(mass/a)
      jmax = a*vcirc
      j = pos[0]*vel[1]-pos[1]*vel[0]
      vbyj = vel*j
      vbyh=[vbyj[1],-vbyj[0]].to_v
      e = vbyh-pos*(1/r)
      eabs=sqrt(e*e)
      argp = atan2(-e[1],-e[0])
      p e
      p eabs
      p a
      p argp
      [e,eabs,a,argp]
    end

    def contactelements(ref)
      fromvectortocontactelements(@mass+ref.mass,
                                @pos-ref.pos,@vel-ref.vel)
    end
    def origincontactelements
      p @mass
      p @pos
      p @vel
      fromvectortocontactelements(1.0,@pos,@vel)
    end
      #
  # relation between J and e
  #
  # assume a=1,j=1
  # ra = (1+e)
  # E= -1/r+v^2/2
  # va = sqrt(1/ra - 1/2a)

    def accumulate_force(otherbody)
      dr= @pos-otherbody.pos
      r2 = dr*dr
      r = sqrt(r2)
      r3 = r2*r
      mfact=@mass/r3
      @acc += dr*(-otherbody.mass/r3)
      otherbody.acc += dr*(@mass/r3)
      @pot -= otherbody.mass/r
      otherbody.pot -= @mass/r
    end
    def forcefromring
      dummyp = Body.new
      dummyp.mass=@ringmass/@ringndivide
      dummyp.acc=@acc
      dummyp.pot=@pot
      dt = 2*Math::PI/@ringndivide
      for i in 0..(@ringndivide-1)
        theta = dt*i
        dummyp.pos[0]=@ringradius*cos(theta)
        dummyp.pos[1]=@ringradius*sin(theta)
        accumulate_force(dummyp)
      end
    end

    def ringinit(radius,mass,ndivide)
      @ringradius,@ringmass,@ringndivide = radius,mass,ndivide
      @withring=1
    end
    def makenull
      null = ([0.0]*(@pos.size)).to_v
    end
  end

  class Nbody
    attr_accessor :ba,  :time
    def initialize
      @ba=[]
    end

    def calculate_force
      @ba.each{|x|
        x.pot=0; x.acc = x.makenull
        x.forcefromring if x.withring
      }
      for i in 0..(@ba.size-2)
        for j in (i+1)..(@ba.size-1)
        @ba[i].accumulate_force(@ba[j])
        end
      end
    end
    def leapfrog(first_call,dt)
      if first_call
        calculate_force
        first_call=nil
      end
      @ba.each{|p| p.vel += p.acc*0.5*dt}
      @ba.each{|p| p.pos += p.vel*dt}
      calculate_force
      @ba.each{|p| p.vel += p.acc*0.5*dt}
    end

    def yoshida6b(first_call,dt)
      d =   [0.784513610477560e0,
        0.235573213359357e0,
        -1.17767998417887e0,
        1.31518632068391e0]
      for i in 0..2 do leapfrog(first_call,dt*d[i]) end
      leapfrog(first_call,dt*d[3])
      for i in 0..2 do leapfrog(first_call,dt*d[2-i]) end
    end
    def calculate_energy
      ke=pe=0
      @ba.each{|p|
        ke += 0.5*p.mass*(p.vel*p.vel)
        pe += p.pot*p.mass
      }
      pe *= 0.5
      [ke+pe,ke,pe]
    end

    def cmadjust
      cmp=@ba[0].makenull
      cmv=@ba[0].makenull
      cmm=0
      @ba.each{|p|
        mass = p.mass
        cmp += p.pos*mass
        cmv += p.vel*mass
        cmm += mass
      }
      cmp *= 1/cmm
      cmv *= 1/cmm
      @ba.each{|p|
        p.pos -= cmp
        p.vel -= cmv
      }
    end
    def plot
      @ba.each{|p|
        relocate(*p.pos)
        dot
      }
    end

    def run(dt,tend,interval,plotendtime=0)
      steps = ((tend+dt*0.5)/dt).to_i
      calculate_force
      e0= calculate_energy
      print "Initial energy = #{e0[0]}\n"
      i=0
      @time=0
      first_call = 1
      steps.times{
        yoshida6b(first_call,dt)
        @time += dt
        i+= 1
        if i==interval or @time < plotendtime
        i=0
        yield
        end
      }
      e= calculate_energy
      print "Final energy = #{e[0]}, de = #{e0[0]-e[0]}\n"
    end
    def xyplot
      @ba.each{|p|
        reloc p.pos[0]  p.pos[1]
        dot
      }
    end

    def taplot
      p "taplot called"
      @ba.each{|p|
        p p
        a=p.origincontactelements
        p a
        reloc @time a[2]
        dot
      }
    end
    def teplot
      p = @ba[1]
      a=p.origincontactelements
      reloc @time a[1]
      dot
    end
    def teplot2
      for i in 0..(@ba.size-1)
        p = @ba[i]
        a=p.origincontactelements
        reloc @time a[1]
        color i+1
        dot
      end
      color 1
    end
    def teplot3(i)
      if i > 0
        i=i-1
        p = @ba[i]
        a=p.origincontactelements
        reloc @time a[1]
        color i+1
        dot
        color 1
      end
    end

  end
  def vcirc(r)
    sqrt(1.0/r)
  end

  def plot(deltar,nbody,tend)
    integrator=open("|solar3","r+");
    mass = 1e-5
    r0=1-deltar
    r1=1
    r2=1+deltar
    mass0=0
    dtplot = 8
    dtplot = tend/1000 if dtplot < tend/1000
    integrator.puts <<-EOD
    0 0 #{tend} 0.1 #{dtplot} #{nbody}
    #{mass} 1e-4  #{r0} 0 0    0 #{vcirc(r0)}        0
    #{mass} 1e-4  #{r1} 0 0    0 #{vcirc(r1)}        0
    #{mass} 1e-4  #{r2} 0 0    0 #{vcirc(r2)}        0
    EOD
    b = Nbody.new
    while s = integrator.gets
      sa = s.split;
      if sa[0] == "singlemass"
        nbody = sa[2].to_i
        b.ba = Array.new(nbody){Body.new}
        p b
        p b.ba
        for i in 0..(nbody-1)
        b.ba[i].mass = sa[3+i].to_f
        end
      end
      if sa[0] == "single"
        b.time = sa[1].to_f
        for i in 0..(nbody-1)
        for j in 0..2
          b.ba[i].pos[j] = sa[2+i*6+j].to_f
          b.ba[i].vel[j] = sa[5+i*6+j].to_f
        end
        end
        b.xyplot
      end
      if sa[0] == "steps"
        print("t=", b.time, " ",s)
      end
    end
  end

  def plot2(deltar,nbody,tend,plotfuncstring)
    integrator=open("|solar3","r+");
    mass = 1e-5
    r0=1-deltar
    r1=1
    r2=1+deltar
    mass0=0
    dtplot = 32
    dtplot = tend/1000 if dtplot < tend/1000
    integrator.puts <<-EOD
    0 0 #{tend} 0.05 #{dtplot} #{nbody}
    #{mass} 1e-4  #{r0} 0 0    0 #{vcirc(r0)}        0
    #{mass} 1e-4  #{r1} 0 0    0 #{vcirc(r1)}        0
    #{mass} 1e-4  #{r2} 0 0    0 #{vcirc(r2)}        0
    EOD
    b = Nbody.new
    while s = integrator.gets
      sa = s.split;
      if sa[0] == "singlemass"
        nbody = sa[2].to_i
        b.ba = Array.new(nbody){Body.new}
        p b
        p b.ba
        for i in 0..(nbody-1)
        b.ba[i].mass = sa[3+i].to_f
        end
      end
      if sa[0] == "single"
        b.time = sa[1].to_f
        for i in 0..(nbody-1)
        for j in 0..2
          b.ba[i].pos[j] = sa[2+i*6+j].to_f
          b.ba[i].vel[j] = sa[5+i*6+j].to_f
        end
        b.ba[i].pos = b.ba[i].pos.to_v
        b.ba[i].vel = b.ba[i].vel.to_v
        end
        eval(plotfuncstring)
      end
      if sa[0] == "steps"
        print("t=", b.time, " ",s)
      end
    end
  end

  unless File.exist?("november-fig1.ps")
    printer november-fig1.ps/vcps
    square
    expand 1.0
    limit -1.3 1.3  -1.3 1.3
    viewport 0 0.45 0.55 1
    box BCTSV  BCNTSV
    ylabel y
    expand 0.5
    ptype 8 0
    plot(0.06,3,5000)

    limit -1.3 1.3  -1.3 1.3
    viewport 1.0 2.0 0 1
    expand 1
    box BCTSV  BCTSV
    expand 0.5
    ptype 8 0
    plot(0.08,3,50000)
    limit -1.3 1.3  -1.3 1.3
    viewport -1 0 -1 0
    expand 1
    box BCNTSV  BCNTSV
    xlabel x
    ylabel y
    expand 0.5
    ptype 8 0
    plot(0.08,3,60000)

    limit -1.3 1.3  -1.3 1.3
    viewport 1.0 2.0 0 1
    expand 1
    box BCNTSV  BCTSV
    xlabel x
    expand 0.5
    ptype 8 0
    plot(0.1,3,900000)
  end
  unless File.exist?("november-fig1b.ps")
    printer november-fig1b.ps/vcps
    square
    expand 1.0

    limit -1.3 1.3  -1.3 1.3
    viewport 0 0.45 0.55 1
    box BCTSV  BCNTSV
    ylabel y
    expand 0.5
    ptype 8 0
    plot(0.06,2,5000)
    limit -1.3 1.3  -1.3 1.3
    viewport 1.0 2.0 0 1
    expand 1
    box BCTSV  BCTSV
    expand 0.5
    ptype 8 0
    plot(0.07,2,50000)

    limit -1.3 1.3  -1.3 1.3
    viewport -1 0 -1 0
    expand 1
    box BCNTSV  BCNTSV
    xlabel x
    ylabel y
    expand 0.5
    ptype 8 0
    plot(0.08,2,500000)
    limit -1.3 1.3  -1.3 1.3
    viewport 1.0 2.0 0 1
    expand 1
    box BCNTSV  BCTSV
    xlabel x
    expand 0.5
    ptype 8 0
    plot(0.1,2,900000)
  end

  unless File.exist?("november-fig2.ps")
    printer november-fig2.ps/vcps
    square
    expand 1.0
    limit 0 2000 0.9 1.1
    viewport 0 1  0.75 1
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.06,3,2000,"b.taplot")

    expand 1.0
    limit 0 60000 0.88 1.12
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.08,3,60000,"b.taplot")
    expand 1.0
    limit 0 900000 0.86 1.14
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    xlabel t
    expand 0.5
    ptype 8 0
    plot2(0.1,3,900000,"b.taplot")

  end
  unless File.exist?("november-fig2b.ps")
    printer november-fig2b.ps/vcps
    square

    expand 1.0
    limit 0 2000 0.9 1.1
    viewport 0 1  0.75 1
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.06,2,2000,"b.taplot")
    expand 1.0
    limit 0 60000 0.88 1.12
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.08,2,60000,"b.taplot")

    expand 1.0
    limit 0 900000 0.86 1.14
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    xlabel t
    expand 0.5
    ptype 8 0
    plot2(0.1,2,900000,"b.taplot")
  end

  unless File.exist?("november-fig2c.ps")
    printer november-fig2c.ps/vcps
    square
    expand 1.0
    limit 0 4000 0.9 1.1
    viewport 0 1  0.8 1
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.07,3,4000,"b.taplot")

    expand 1.0
    limit 0 60000 0.88 1.12
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.08,3,60000,"b.taplot")
    expand 1.0
    limit 0 100000 0.86 1.14
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    xlabel t
    expand 0.5
    ptype 8 0
    plot2(0.09,3,100000,"b.taplot")

    expand 1.0
    limit 0 900000 0.86 1.14
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    xlabel t
    expand 0.5
    ptype 8 0
    plot2(0.1,3,900000,"b.taplot")
    expand 1.0
    limit 0 900000 0.86 1.14
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    xlabel t
    expand 0.5
    ptype 8 0
    plot2(0.11,3,1800000,"b.taplot")

  end
  unless File.exist?("november-fig3.ps")
    printer november-fig3.ps/vcps
    square

    expand 1.0
    viewport 0.1 1.1 0.55 1
    limit 0 60000 -0.001 0.02
    box BCTSV  BCNTSV
    expand 1.5
    reloc -8000 0.012
    label e
    ptype 8 0
    expand 0.5
    plot2(0.08,3,60000,"b.teplot")
    expand 1
    expand 1.0
    viewport 0 1 -1 0
    limit 0 60000 -0.001 0.02
    box BCNTSV  BCNTSV
    expand 1.5
    reloc -8000 0.012
    label e
    xlabel t
    ptype 8 0
    expand 0.5
    plot2(0.08,2,60000,"b.teplot")
    expand 1
  end

  unless File.exist?("november-fig3b.ps")
    printer november-fig3b.ps/vcps
    square
    expand 1.0
    viewport 0.1 1.1 0.55 1
    limit 0 60000 -0.001 0.02
    box BCTSV  BCNTSV
    expand 1.5
    reloc -8000 0.012
    label e
    ptype 8 0
    expand 0.5
    plot2(0.08,3,60000,"b.teplot2")
    expand 1

    expand 1.0
    viewport 0 1 -1 0
    limit 0 60000 -0.001 0.02
    box BCNTSV  BCNTSV
    expand 1.5
    reloc -8000 0.012
    label e
    xlabel t
    ptype 8 0
    expand 0.5
    plot2(0.08,2,60000,"b.teplot2")
    expand 1
  end
  unless File.exist?("november-fig3c.ps")
    printer november-fig3c.ps/vcps
    term
    square

    expand 1.0
    viewport 0.1 1.1 0.55 1
    limit 0 900000 -0.001 0.02
    box BCTSV  BCNTSV
    expand 1.5
    reloc -120000 0.012
    label e
    ptype 8 0
    expand 0.5
    plot2(0.1,3,900000,"b.teplot2")
    expand 1
    expand 1.0
    viewport 0 1 -1 0
    limit 0 900000 -0.001 0.02
    box BCNTSV  BCNTSV
    expand 1.5
    reloc -120000 0.012
    label e
    xlabel t
    ptype 8 0
    expand 0.5
    plot2(0.1,2,900000,"b.teplot2")
    expand 1
  end

赤木: えーと、ちょっと待って、本文で使った図は3枚だけよね?

学生S: ですね。 november-fig1.ps と、後は 2 と 3c です。

赤木: その辺、使ってないのを消してみて。

学生S: ちょっと待って下さい。こんな感じかな?


  # plot.rb : postprocessing of three body integration

  require 'vector.rb'
  include Math
  class Body
    attr_accessor :mass, :pos, :vel, :acc, :pot, :withring
    def initialize
      @pos=[]
      @vel=[]
      @withring=nil
    end
    def froms(s)
      a=s.split.collect{|x| x.to_f}
      p a
      n=(a.size-1)/2
      @mass = a[0]
      @pos  = a[1..n].to_v
      @vel  = a[(n+1)..(2*n)].to_v
      p self
    end

    def convertfromelements(a,e)
      @pos[0]=a*(1+e)
      @pos[1]=0
      @vel[0]=0
      @vel[1]= sqrt(2*(1/@pos[0]-0.5/a))
      @pos=@pos.to_v
      @vel=@vel.to_v
    end
    def fromvectortocontactelements(mass,pos,vel)
      r=sqrt(pos*pos)
      be = -mass/r+0.5*(vel*vel)
      # be = -mass/2a, a = -mass/2be
      a = -mass/(be*2)
      vcirc = sqrt(mass/a)
      jmax = a*vcirc
      j = pos[0]*vel[1]-pos[1]*vel[0]
      vbyj = vel*j
      vbyh=[vbyj[1],-vbyj[0]].to_v
      e = vbyh-pos*(1/r)
      eabs=sqrt(e*e)
      argp = atan2(-e[1],-e[0])
      p e
      p eabs
      p a
      p argp
      [e,eabs,a,argp]
    end

    def contactelements(ref)
      fromvectortocontactelements(@mass+ref.mass,
                                @pos-ref.pos,@vel-ref.vel)
    end
    def origincontactelements
      p @mass
      p @pos
      p @vel
      fromvectortocontactelements(1.0,@pos,@vel)
    end
  end

  class Nbody
    attr_accessor :ba,  :time
    def initialize
      @ba=[]
    end

    def plot
      @ba.each{|p|
        relocate(*p.pos)
        dot
      }
    end
    def xyplot
      @ba.each{|p|
        reloc p.pos[0]  p.pos[1]
        dot
      }
    end

    def taplot
      p "taplot called"
      @ba.each{|p|
        p p
        a=p.origincontactelements
        p a
        reloc @time a[2]
        dot
      }
    end
    def teplot2
      for i in 0..(@ba.size-1)
        p = @ba[i]
        a=p.origincontactelements
        reloc @time a[1]
        color i+1
        dot
      end
      color 1
    end

  end
  def vcirc(r)
    sqrt(1.0/r)
  end

  def plot(deltar,nbody,tend)
    integrator=open("|solar3","r+");
    mass = 1e-5
    r0=1-deltar
    r1=1
    r2=1+deltar
    mass0=0
    dtplot = 8
    dtplot = tend/1000 if dtplot < tend/1000
    integrator.puts <<-EOD
    0 0 #{tend} 0.1 #{dtplot} #{nbody}
    #{mass} 1e-4  #{r0} 0 0    0 #{vcirc(r0)}        0
    #{mass} 1e-4  #{r1} 0 0    0 #{vcirc(r1)}        0
    #{mass} 1e-4  #{r2} 0 0    0 #{vcirc(r2)}        0
    EOD
    b = Nbody.new
    while s = integrator.gets
      sa = s.split;
      if sa[0] == "singlemass"
        nbody = sa[2].to_i
        b.ba = Array.new(nbody){Body.new}
        p b
        p b.ba
        for i in 0..(nbody-1)
        b.ba[i].mass = sa[3+i].to_f
        end
      end
      if sa[0] == "single"
        b.time = sa[1].to_f
        for i in 0..(nbody-1)
        for j in 0..2
          b.ba[i].pos[j] = sa[2+i*6+j].to_f
          b.ba[i].vel[j] = sa[5+i*6+j].to_f
        end
        end
        b.xyplot
      end
      if sa[0] == "steps"
        print("t=", b.time, " ",s)
      end
    end
  end

  def plot2(deltar,nbody,tend,plotfuncstring)
    integrator=open("|solar3","r+");
    mass = 1e-5
    r0=1-deltar
    r1=1
    r2=1+deltar
    mass0=0
    dtplot = 32
    dtplot = tend/1000 if dtplot < tend/1000
    integrator.puts <<-EOD
    0 0 #{tend} 0.05 #{dtplot} #{nbody}
    #{mass} 1e-4  #{r0} 0 0    0 #{vcirc(r0)}        0
    #{mass} 1e-4  #{r1} 0 0    0 #{vcirc(r1)}        0
    #{mass} 1e-4  #{r2} 0 0    0 #{vcirc(r2)}        0
    EOD
    b = Nbody.new
    while s = integrator.gets
      sa = s.split;
      if sa[0] == "singlemass"
        nbody = sa[2].to_i
        b.ba = Array.new(nbody){Body.new}
        p b
        p b.ba
        for i in 0..(nbody-1)
        b.ba[i].mass = sa[3+i].to_f
        end
      end
      if sa[0] == "single"
        b.time = sa[1].to_f
        for i in 0..(nbody-1)
        for j in 0..2
          b.ba[i].pos[j] = sa[2+i*6+j].to_f
          b.ba[i].vel[j] = sa[5+i*6+j].to_f
        end
        b.ba[i].pos = b.ba[i].pos.to_v
        b.ba[i].vel = b.ba[i].vel.to_v
        end
        eval(plotfuncstring)
      end
      if sa[0] == "steps"
        print("t=", b.time, " ",s)
      end
    end
  end

  unless File.exist?("november-fig1.ps")
    printer november-fig1.ps/vcps
    square
    expand 1.0
    limit -1.3 1.3  -1.3 1.3
    viewport 0 0.45 0.55 1
    box BCTSV  BCNTSV
    ylabel y
    expand 0.5
    ptype 8 0
    plot(0.06,3,5000)

    limit -1.3 1.3  -1.3 1.3
    viewport 1.0 2.0 0 1
    expand 1
    box BCTSV  BCTSV
    expand 0.5
    ptype 8 0
    plot(0.08,3,50000)
    limit -1.3 1.3  -1.3 1.3
    viewport -1 0 -1 0
    expand 1
    box BCNTSV  BCNTSV
    xlabel x
    ylabel y
    expand 0.5
    ptype 8 0
    plot(0.08,3,60000)

    limit -1.3 1.3  -1.3 1.3
    viewport 1.0 2.0 0 1
    expand 1
    box BCNTSV  BCTSV
    xlabel x
    expand 0.5
    ptype 8 0
    plot(0.1,3,900000)
    pgend
    psfix
  end
  unless File.exist?("november-fig2.ps")
    printer november-fig2.ps/vcps
    square

    expand 1.0
    limit 0 2000 0.9 1.1
    viewport 0 1  0.75 1
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.06,3,2000,"b.taplot")
    expand 1.0
    limit 0 60000 0.88 1.12
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    expand 0.5
    ptype 8 0
    plot2(0.08,3,60000,"b.taplot")

    expand 1.0
    limit 0 900000 0.86 1.14
    viewport 0 1  -1.25 -0.25
    box BCNTSV  BCNTSV
    ylabel a
    xlabel t
    expand 0.5
    ptype 8 0
    plot2(0.1,3,900000,"b.taplot")
    pgend
    psfix
  end

  unless File.exist?("november-fig3c.ps")
    printer november-fig3c.ps/vcps
    term
    square
    expand 1.0
    viewport 0.1 1.1 0.55 1
    limit 0 900000 -0.001 0.02
    box BCTSV  BCNTSV
    expand 1.5
    reloc -120000 0.012
    label e
    ptype 8 0
    expand 0.5
    plot2(0.1,3,900000,"b.teplot2")
    expand 1

    expand 1.0
    viewport 0 1 -1 0
    limit 0 900000 -0.001 0.02
    box BCNTSV  BCNTSV
    expand 1.5
    reloc -120000 0.012
    label e
    xlabel t
    ptype 8 0
    expand 0.5
    plot2(0.1,2,900000,"b.teplot2")
    expand 1
    pgend
    psfix
  end

赤木: 長さ半分ね。ちょっとはいいかしら。最初の図を書いてるところは?

学生S: これですね。

 unless File.exist?("november-fig1.ps")
   printer november-fig1.ps/vcps
   square
   expand 1.0

   limit -1.3 1.3  -1.3 1.3
   viewport 0 0.45 0.55 1
   box BCTSV  BCNTSV
   ylabel y
   expand 0.5
   ptype 8 0
   plot(0.06,3,5000)
   limit -1.3 1.3  -1.3 1.3
   viewport 1.0 2.0 0 1
   expand 1
   box BCTSV  BCTSV
   expand 0.5
   ptype 8 0
   plot(0.08,3,50000)

   limit -1.3 1.3  -1.3 1.3
   viewport -1 0 -1 0
   expand 1
   box BCNTSV  BCNTSV
   xlabel x
   ylabel y
   expand 0.5
   ptype 8 0
   plot(0.08,3,60000)
   limit -1.3 1.3  -1.3 1.3
   viewport 1.0 2.0 0 1
   expand 1
   box BCNTSV  BCTSV
   xlabel x
   expand 0.5
   ptype 8 0
   plot(0.1,3,900000)
   pgend
   psfix
 end

枠書いたあと plot っていうのを呼んでるだけみたいですね。

赤木: そうね。 plot の中身は?

学生S: これです:

 def plot(deltar,nbody,tend)
   integrator=open("|solar3","r+");
   mass = 1e-5
   r0=1-deltar
   r1=1
   r2=1+deltar
   mass0=0
   dtplot = 8
   dtplot = tend/1000 if dtplot < tend/1000
   integrator.puts <<-EOD
   0 0 #{tend} 0.1 #{dtplot} #{nbody}
   #{mass} 1e-4  #{r0} 0 0    0 #{vcirc(r0)}        0
   #{mass} 1e-4  #{r1} 0 0    0 #{vcirc(r1)}        0
   #{mass} 1e-4  #{r2} 0 0    0 #{vcirc(r2)}        0
   EOD

   b = Nbody.new
   while s = integrator.gets
     sa = s.split;
     if sa[0] == "singlemass"
       nbody = sa[2].to_i
       b.ba = Array.new(nbody){Body.new}
       p b
       p b.ba
       for i in 0..(nbody-1)
        b.ba[i].mass = sa[3+i].to_f
       end
     end
     if sa[0] == "single"
       b.time = sa[1].to_f
       for i in 0..(nbody-1)
        for j in 0..2
          b.ba[i].pos[j] = sa[2+i*6+j].to_f
          b.ba[i].vel[j] = sa[5+i*6+j].to_f
        end
       end
       b.xyplot
     end
     if sa[0] == "steps"
       print("t=", b.time, " ",s)
     end
   end
 end

赤木: 最初の

   integrator=open("|solar3","r+");

がわからないわよね。

学生S: これはサブプロセスを作るというやつです。 solar3 という別 の、ここでは作者が書いた C のプログラムが実行されて、その入力と出力が integrator の書き込みと読み出しになるものです。 UNIX のコマンドでのパイプと同 様ですが、双方向なのがちょっと違います。

赤木: あ、だからその後の

   integrator.puts <<-EOD
   0 0 #{tend} 0.1 #{dtplot} #{nbody}
   #{mass} 1e-4  #{r0} 0 0    0 #{vcirc(r0)}        0
   #{mass} 1e-4  #{r1} 0 0    0 #{vcirc(r1)}        0
   #{mass} 1e-4  #{r2} 0 0    0 #{vcirc(r2)}        0
   EOD

で solar3 の入力を与えてるわけね。なるほど。

学生S: そうです。その後の while では solar3 に出力を受け取って グラフを書くわけです。 solar3 の出力は

  E, |DE/E0|,  PE, VR: -1.50403e-05 0.000000e+00 -3.00764e-05    -0.499931
 am:                     0                     0    2.999098985582e-05
 singlemass nb= 3   1e-05 1e-05 1e-05
 single 0  0.931199 0.128523 0 -0.140592 1.02172 0 0.992198 0.124675 0 -0.124674 0.992199 0 1.05303 0.121146 0 -0.111438 0.964934 0
 steps =       0 E, |DE/E0|,  PE, VR: -1.50403e-05 1.135362e-13 -3.00763e-05    -0.499929
 single 8  -0.802115 0.485464 0 -0.531287 -0.884667 0 -0.211358 0.977229 0 -0.977314 -0.211439 0 0.511369 0.929056 0 -0.853098 0.46786 0
 steps =     123 E, |DE/E0|,  PE, VR: -1.50403e-05 1.485657e-13 -3.00974e-05    -0.500279
 single 16  0.402918 -0.8415 0 0.936139 0.445985 0 -0.933565 -0.356606 0 0.357279 -0.934403 0 -0.511405 0.929818 0 -0.853205 -0.466523 0
 steps =     251 E, |DE/E0|,  PE, VR: -1.50403e-05 1.280662e-13 -3.01484e-05    -0.501125
 single 24  0.137574 0.929347 0 -1.01865 0.150961 0 0.495079 -0.867864 0 0.86919 0.495659 0 -1.06268 0.0716562 0 -0.0692701 -0.966385 0

で、僕のプログラムでは singlemass で始まる行の 4, 5, 6 列目が粒子の質 量、 single で始まる行の 3 列目以降が粒子の位置、速度がその順番でとなっ てます。だから、 solar3 はそういうふうに出力するんだと思います。

赤木: 読んでから何をするかっていうと、 xyplot ね。これは?

学生S: これですね。

  def xyplot
    @ba.each{|p|
      reloc p.pos[0]  p.pos[1]
      dot
    }
  end

それぞれの粒子の xy 座標で点を打つだけです。

赤木: 他の図も同じようなものかしら?

学生S: えーと、そうですね、 plot2 ってのになってますね。これは なんか怪しい文字列を引数にしてて、さっきの xyplot の代わりにその引数で

      eval(plotfuncstring)

してます。

赤木: eval って?

学生S: その文字列をプログラムとして実行するものです。 Ruby はイ ンタプリタだから、そういうこともできるわけです。

赤木: で、図2では taplot、図3では teplot2 というのを呼んでるわけね。

学生S: です。これらは軌道要素に変換して時間を横軸にしてプロット するだけです。

赤木: まあ、君が書いたところはそれでいいわ。 solar3 は?

学生S: すみません、これはちょっとまたの機会に、、、

赤木: んー、そうね、他の章先にやったほうがいいわね。

Previous ToC Next