require "acs"

class Body

  attr_accessor :mass, :pos, :vel

  def initialize(mass = 0, pos = Vector[0,0,0], vel = Vector[0,0,0])
    @mass, @pos, @vel = mass, pos, vel
  end

end

class NBody

  attr_accessor :time, :body

  def initialize(n = 0)
    @body = []
    for i in 0...n
      @body[i] = Body.new
    end
  end

end

include Math

def frand(low, high)
  low + rand * (high - low)
end

def spherical(r)
  vector = Vector.new
  theta = acos(frand(-1, 1))
  phi = frand(0, 2*PI)
  vector[0] = r * sin( theta ) * cos( phi )
  vector[1] = r * sin( theta ) * sin( phi )
  vector[2] = r * cos( theta )
  vector
end  

def mkplummer(n, seed)
  if seed == 0
    srand
  else
    srand seed
  end
  nb = NBody.new(n)
  nb.body.each do |b|                                                        
    b.mass = 1.0/n                                                           
    radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0)                           
    b.pos = spherical(radius)                                                
    x = 0.0                                                                  
    y = 0.1                                                                  
    while y > x*x*(1.0-x*x)**3.5                                             
      x = frand(0,1)                                                         
      y = frand(0,0.1)                                                       
    end                                                                      
    velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)               
    b.vel = spherical(velocity)                                              
  end                                                                        
  STDERR.print "             actual seed used\t: ", srand, "\n"
  nb.acs_write
end

options_text= <<-END

  Description: Plummer's Model Builder
  Long description:
    This program creates an N-body realization of Plummer's Model.
    (c) 2004, Piet Hut and Jun Makino; see ACS at www.artcompsi.org

    The algorithm used is described in Aarseth, S., Henon, M., & Wielen, R.,
    Astron. Astroph. 37, 183 (1974).


  Short name:           -n
  Long name:            --n_particles
  Value type:           int
  Default value:        1
  Variable name:        n
  Print name:           N
  Description:          Number of particles
  Long description:
    Number of particles in a realization of Plummer's Model.

    Each particles is drawn at random from the Plummer distribution,
    and therefore there are no correlations between the particles.

    Physical Units are used in which G = M = a = 1, where
      G is the gravitational constant
      M is the total mass of the N-body system
      a is the structural length, with potential U(r) = GM/(r^2 + a^2)^{1/2}


  Short name:           -s
  Long name:            --seed
  Value type:           int
  Default value:        0
  Description:          pseudorandom number seed given
  Print name:           
  Variable name:      seed
  Long description:
    Seed for the pseudorandom number generator.  If a seed is given with
    value zero, a preudorandom number is chosen as the value of the seed.
    The seed value used is echoed separately from the seed value given,
    to allow the possibility to repeat the creation of an N-body realization.

      Example:

        |gravity> kali mkplummer1.rb -n 42 -s 0
        . . .
        pseudorandom number seed given  : 0
                     actual seed used   : 1087616341
        . . .
        |gravity> kali mkplummer1.rb -n 42 -s 1087616341
        . . .
        pseudorandom number seed given  : 1087616341
                     actual seed used   : 1087616341
        . . .


  END

c = parse_command_line(options_text)

mkplummer(c.n, c.seed)