#!/usr/bin/ruby -w

=begin
/***************************************************************************
 *   Copyright (C) 2008, Paul Lutus                                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
=end

require 'gravity_ui'

PROGRAM_VERSION = "1.5"

# main class

class Gravity < GravityGlade

   Gravity::PlanetColors = [ Gdk::Color.parse("white"),Gdk::Color.parse("yellow"),
      Gdk::Color.parse("cyan"), Gdk::Color.new(128 * 256,128 * 256,255 * 256),
      Gdk::Color.parse("red"),Gdk::Color.parse("green"),
   Gdk::Color.parse("magenta"),Gdk::Color.parse("blue") ]

   Gravity::AnimStrings = [ "1 hour","2 hours","4 hours","8 hours","16 hours",
   "1 day","2 days","4 days","8 days","16 days","32 days","64 days","128 days","256 days" ]

   Gravity::CometStrings = [ "1","2","4","8","16","32","64","128","256","512" ]

   def initialize(path,xxx,name)
      super(path,xxx,name)
      get_widgets()
      @program_name = self.class.name
      @program_title = @program_name + " " + PROGRAM_VERSION
      GravityUI().set_title(@program_title)

      @total_time_hours = 0
      @initialized = false
      # this sets the sleep interval (units are seconds)
      @anim_time = 0.005
      @anim_thread = nil
      @anim_flag = false
      # initial drawing scale, change with mouse wheel
      @drawing_scale = 6e-12
      @rotx = -20
      @roty = 0
      @planet_list = []
      @time_step = nil
      @pixmap = nil
      @erase = true
      @anaglyph_mode = false
      @symmetric = false
      @nice = true
      @mouse_down = false

      @color_cyan = Gdk::Color.parse("cyan")
      @color_red = Gdk::Color.parse("red")
      @color_black = Gdk::Color.parse("black")
      @rotator = RotationMatrix.new
      @time_step_combobox_manager = ComboBoxManager.new(time_step_combobox(),AnimStrings,"16 hours")
      set_time_step(true)
      @comet_combobox_manager = ComboBoxManager.new(comet_combobox(),CometStrings,"16")
      solar_system_checkbutton.active = true
      comets_checkbutton.active = true
      # graphic_pane.signal_connect("expose_event") do
      #   draw_image
      # end
      @min_draw_radius = 5
      pixels_spinbutton.value = @min_draw_radius
      load_objects(true)
      @initialized = true
   end

   def get_widgets()
      @glade.widget_names.each do |name|
         # create accessor methods for each defined widget
         eval("def #{name}() return @glade.get_widget(\"#{name}\") end")
      end
   end

   def Gravity::message_dialog(window,message,inquiry = false)
      if inquiry
         dlg = Gtk::MessageDialog.new(nil,
         Gtk::MessageDialog::MODAL,
         Gtk::MessageDialog::QUESTION,
         Gtk::MessageDialog::BUTTONS_YES_NO,
         message)
      else # just an alert
         dlg = Gtk::MessageDialog.new(nil,
         Gtk::MessageDialog::MODAL,
         Gtk::MessageDialog::INFO,
         Gtk::MessageDialog::BUTTONS_OK,
         message)
      end
      dlg.set_title(window.class.name)
      response = dlg.run
      dlg.destroy
      return response == Gtk::MessageDialog::RESPONSE_YES || response == Gtk::MessageDialog::RESPONSE_OK
   end

   def close()
      Gtk.main_quit
   end


   def set_time_step(force = false)
      if(@initialized || force)
         s = @time_step_combobox_manager.active_string()
         v,units = s.split(" ")
         v = v.to_i
         v *= 3600 if units =~ /hour/i
         v *= 86400 if units =~ /day/i
         @time_step = v
      end
   end

   def show_status
      y = @total_time_hours
      h = y % 24
      y /= 24
      y *= 100
      d = (y % 36525) / 100
      y /= 36525
      str = sprintf("Elapsed time: %04dy %03dd %02dh",y,d,h)
      status_bar().text = str
   end

   def draw_planets(xsize,ysize,gc,anaglyph_flag = nil,td_color = nil)
      if(anaglyph_flag)
         gc.rgb_fg_color = td_color
      end
      i = 0
      @planet_list.each do |planet|
         v = planet.pos * @drawing_scale
         @rotator.rotate(v)
         @rotator.convert_3d_to_2d(v,anaglyph_flag)
         sxa = @x_screen_center + (v.x * @screen_scale)
         sya = @y_screen_center - (v.y * @screen_scale)
         if(sxa >= 0 && sxa < xsize && sya >= 0 && sya < ysize)
            # fake the sun's radius for aesthetics
            sr = (i == 0)?4e7:(planet.radius)
            s = Cart3.new(sr * @drawing_scale,0,-planet.pos.z * @drawing_scale);
            @rotator.convert_3d_to_2d(s)
            s.x *= @screen_scale * 100
            s.x = (s.x < @min_draw_radius)?@min_draw_radius:s.x
            sc = s.x/2
            unless(anaglyph_flag)
               col = PlanetColors[i % PlanetColors.size]
               gc.rgb_fg_color = col
            end
            @pixmap.draw_arc(gc,true,sxa-sc,sya-sc,s.x,s.x,0,23040)
         end
         i += 1
      end
   end

   def fill_block(gc,color,sx,sy,wx,wy)
      gc.fill = Gdk::GC::Fill::SOLID
      gc.rgb_fg_color = color
      @pixmap.draw_rectangle(gc,true,sx,sy,wx+2,wy+2)
   end

   def draw_image(erase = false)
      if(@initialized && @planet_list)
         show_status
         @rotator.populate_matrix(@rotx,@roty)
         alloc = graphic_pane().allocation
         xsize,ysize = alloc.width,alloc.height
         # create an off-screen pixmap for image drawing
         # whenever a change requires it
         if(@pixmap == nil || xsize != @old_xsize || ysize != @old_ysize)
            @pixmap = Gdk::Pixmap.new(graphic_pane().window,xsize,ysize,-1)
            @old_xsize = xsize
            @old_ysize = ysize
            @x_screen_center = xsize / 2
            @y_screen_center = ysize / 2
            @screen_scale = (@x_screen_center > @y_screen_center)?@y_screen_center:@x_screen_center

         end
         gc = Gdk::GC.new(@pixmap)
         # set image background to black
         fill_block(gc,@color_black,0,0,@old_xsize,@old_ysize) if @erase || erase
         if(@anaglyph_mode)
            # In 3D mode, let overlapping red & cyan lines appear white
            gc.function = Gdk::GC::OR
            # draw complete, rotated right-hand and left-hand
            # images in cyan and red for anaglyphic glasses
            draw_planets(xsize,ysize,gc,-1,@color_cyan) # right eye image
            draw_planets(xsize,ysize,gc,1,@color_red) # left eye image
            gc.function = Gdk::GC::COPY
         else
            # draw image once
            draw_planets(xsize,ysize,gc)
         end
         # pixpainter.end
         # move the completed image to the screen
         graphic_pane().window.draw_drawable(gc,@pixmap,0,0,0,0,@old_xsize,@old_ysize)
         @total_time_hours += @time_step / 3600
      end
   end

   def perform_orbit_calc
      OrbitalPhysics::process_planets(@planet_list,@time_step,@symmetric)
      draw_image
   end

   def toggle_animation
      if(@anim_flag)
         @anim_flag = false
         while(@anim_thread.alive?)
         end
      else
         @total_time_hours = 0
         @anim_flag = true
         @anim_thread = Thread.new {
            while @anim_flag
               # must have some thread "sleep" time
               # or GUI updates will stop
               sleep (@nice)?@anim_time:0.002
               perform_orbit_calc()
            end
         }
      end
      run_stop_button.label = ((@anim_thread == nil)?"Run":"Stop")
      draw_image(true)
   end

   def load_comets(force = false)
      if(@initialized || force)
         n = @comet_combobox_manager.active_string().to_i
         1.upto(n) do |i|
            name = "comet#{i}"
            ca = rand(360) # angle in x-z plane
            cr = rand(4e11) + 4e11 # distance from sun
            pos = Cart3.new(cr * Math.sin(ca * CommonCode::ToRad),0,cr * Math.cos(ca * CommonCode::ToRad))
            # comet initial velocity
            v = (rand(200) + 100) * 50.0
            v = (i % 2 == 1)?-v:v
            vel = Cart3.new(0,v,0)
            comet = Planet.new(name,1e3,pos,vel,1e9)
            @planet_list << comet
         end
      end
   end

   def load_orbital_data(data,sun_only = false)
      data = data.split("\n")
      data.shift # drop header line
      data.each do |line|
         fields = line.split(",")
         pos = Cart3.new(-fields[1].to_f,0,0)
         vel = Cart3.new(0,0,fields[4].to_f)
         planet = Planet.new(fields[0],fields[2].to_f,pos,vel,fields[3].to_f)
         @planet_list << planet
         break if sun_only
      end
   end

   def load_objects(force = false)
      @planet_list = []
      sun_only = !solar_system_checkbutton.active?
      load_orbital_data(SolarSystem::Data,sun_only)
      load_comets(force) if comets_checkbutton.active?
      draw_image(true)
   end

   def beep
      Gdk.beep
   end

   # close application cleanly

   def close(*x)
      Gtk.main_quit
   end

   # mouse events

   def mouseMoveEvent (e)
      if(@mouse_down)
         # rotate image by dragging mouse
         dx = (e.y - @mouse_press_x) / 2
         dy = (e.x - @mouse_press_y) / 2
         @rotx = @mouse_press_rx - dx
         @roty = @mouse_press_ry - dy
         draw_image(true)
      end
   end

   def mousePressEvent (e)
      # set up to control rotation
      # by dragging mouse
      @mouse_down = true
      @mouse_press_rx = @rotx
      @mouse_press_ry = @roty
      @mouse_press_x = e.y
      @mouse_press_y = e.x
      draw_image(true)
   end

   def mouseReleaseEvent (e)
      @mouse_down = false
      # set up to control rotation
      # by dragging mouse
      draw_image(true)
   end

   def wheelEvent (e)
      # change drawing scale using mouse wheel
      # get mouse wheel delta
      v = (e.direction == Gdk::EventScroll::DOWN)?-1:1
      v = v.to_f
      v = 1.0 + (v/10.0)
      @drawing_scale *= v
      draw_image(true)
   end

   # action handlers

   def on_GravityUI_delete_event(widget, arg0)
      close
   end

   def on_quit_button_clicked(widget)
      close
   end

   def on_run_stop_button_clicked(widget)
      toggle_animation
   end

   def on_eventbox1_button_press_event(widget, arg0)
      mousePressEvent(arg0)
   end

   def on_eventbox1_button_release_event(widget, arg0)
      mouseReleaseEvent(arg0)
   end

   def on_eventbox1_motion_notify_event(widget, arg0)
      mouseMoveEvent(arg0)
   end

   def on_eventbox1_scroll_event(widget, arg0)
      wheelEvent(arg0)
   end

   def on_time_step_combobox_changed(widget)
      set_time_step
   end
   def on_anaglyphic_checkbutton_toggled(widget)
      @anaglyph_mode = anaglyphic_checkbutton.active?
      draw_image(true)
   end
   def on_comet_combobox_changed(widget)
      comets_checkbutton.active = true
      load_objects
      draw_image(true)
   end
   def on_step_button_clicked(widget)
      toggle_animation if @anim_flag
      perform_orbit_calc
   end
   def on_pixels_spinbutton_value_changed(widget)
      @min_draw_radius = pixels_spinbutton.value
      draw_image(true)
   end
   def on_solar_system_checkbutton_toggled(widget)
      load_objects
   end
   def on_trails_checkbutton_toggled(widget)
      @erase = !trails_checkbutton.active?
   end

   def on_nice_checkbutton_toggled(widget)
      @nice = nice_checkbutton.active?
   end
   def on_comets_checkbutton_toggled(widget)
      load_objects
   end

   def on_graphic_pane_expose_event(widget, arg0)
      draw_image(true)
   end


end # class Gravity

class ComboBoxManager
   def initialize(box,list,default = nil)
      @box = box
      @hash = {}
      index = 0
      # a placeholder item is required to get around
      # a bug in the Glade designer that won't create
      # a sane combobox without it. So first,
      # remove the placeholder item
      @box.remove_text(0)
      list.each do |item|
         @box.append_text(item)
         @hash[item] = index
         index += 1
      end
      if (@hash[default])
         @box.set_active(@hash[default])
      else
         @box.set_active(0)
      end
   end
   def active_string()
      return @box.active_text
   end
end

# a class for routines and constants of common utility

class CommonCode
   CommonCode::ToRad = Math::PI / 180.0
   CommonCode::ToDeg = 180.0 / Math::PI
   def CommonCode::fmt_num(n)
      sprintf("%.2e",n)
   end
end

# solar system data class, all units mks

class SolarSystem
   SolarSystem::Data=<<-EOF
   "Name","OrbitRad","BodyRad","Mass","OrbitVel"
   "Sun",0,695000000,1.989E+030,0
   "Mercury",57900000000,2440000,3.33E+023,47900
   "Venus",108000000000,6050000,4.869E+024,35000
   "Earth",150000000000,6378140,5.976E+024,29800
   "Mars",227940000000,3397200,6.421E+023,24100
   "Jupiter",778330000000,71492000,1.9E+027,13100
   "Saturn",1429400000000,60268000,5.688E+026,9640
   "Uranus",2870990000000,25559000,8.686E+025,6810
   "Neptune",4504300000000,24746000,1.024E+026,5430
   "Pluto",5913520000000,1137000,1.27E+022,4740
   EOF
end

# a Cartesian 3D vector class
# with a number of important operator overrides

class Cart3
   attr_accessor :x,:y,:z
   def initialize(x = 0,y = 0,z = 0)
      if(x.class == self.class)
         @x = x.x
         @y = x.y
         @z = x.z
      else
         @x = x
         @y = y
         @z = z
      end
   end

   def -(e)
      Cart3.new(@x - e.x,@y - e.y,@z - e.z)
   end

   def +(e)
      Cart3.new(@x + e.x,@y + e.y,@z + e.z)
   end

   def *(e)
      if(e.class != self.class)
         # multiply by scalar
         Cart3.new(@x * e,@y * e,@z * e)
      else
         # multiply by vector
         Cart3.new(@x * e.x,@y * e.y,@z * e.z)
      end
   end

   def /(e)
      if(e.class != self.class)
         # divide by scalar
         Cart3.new(@x / e,@y / e,@z / e)
      else
         # divide by vector
         Cart3.new(@x / e.x,@y / e.y,@z / e.z)
      end
   end

   # sum of squares
   def sumsq
      @x*@x+@y*@y+@z*@z
   end

   def to_s
      "[#{CommonCode::fmt_num(@x)},#{CommonCode::fmt_num(@y)},#{CommonCode::fmt_num(@z)}]"
   end
end # class Cart3

=begin

Planet, a simple data class

name = string
radius = float
pos = Cart3
vel = Cart3
mass = float

=end

class Planet
   attr_accessor :name,:radius,:pos,:vel,:mass
   def initialize(name,radius,pos,vel,mass = 0)
      @name = name.gsub(/"/,"")
      @radius = radius
      @pos = pos
      @vel = vel
      @mass = mass
   end
   def to_s
      "#{@name},#{CommonCode::fmt_num(@radius)},#{@pos},#{@vel},#{CommonCode::fmt_num(@mass)}"
   end
end

# RotationMatrix performs 3D rotations and perspective

class RotationMatrix

   # perspective depth cue for 3D -> 2D transformation
   PerspectiveDepth = 16
   # empirical constant for anaglyphic rotation
   AnaglyphScale = 0.03

   RotationMatrix::ToRad = Math::PI / 180.0
   RotationMatrix::ToDeg = 180.0 / Math::PI

   # populate 3D matrix with values for x,y,z rotations

   def populate_matrix(xa,ya)

      # create trig values
      @sy = Math.sin(xa * RotationMatrix::ToRad);
      @cy = Math.cos(xa * RotationMatrix::ToRad);
      @sx = Math.sin(ya * RotationMatrix::ToRad);
      @cx = Math.cos(ya * RotationMatrix::ToRad);
   end

   # 3D -> 2D, add perspective cue,
   # perform anaglyph position change if specified

   def convert_3d_to_2d(v,anaglyph_flag = 0)
      v.x = (v.x * (PerspectiveDepth + v.z))/PerspectiveDepth
      v.x += v.z * anaglyph_flag * AnaglyphScale if anaglyph_flag
      v.y = (v.y * (PerspectiveDepth + v.z))/PerspectiveDepth
   end

   # rotate a 3D point using matrix values

   def rotate(v)
      # borrowed from my "Apple World" 1979
      hf = (v.x * @sx - v.z * @cx)
      py = v.y * @cy + @sy * hf
      px = v.x * @cx + v.z * @sx
      pz = -v.y * @sy + @cy * hf
      v.x = px; v.y = py; v.z = pz
   end
end # class RotationMatrix

=begin

Gravitational force f, Newtons, between two masses M and m:

f = G M m
   ------
    r^2

G = universal gravitational constant, 6.6742e-11 N m^2 / kg^2

M,m = masses of the two bodies, kilograms

r = radius (distance) between M and m, meters

acceleration, m/s, a for a force f and a mass m:

a = f/m

The shorthand version drops one of the masses
for a slight speed improvement, and goes
directly to acceleration a:

a = G M
   -----
    r^2

BUT the shorthand form assumes one
of the masses is infinite

In a numerical simulation using time slice dt:

velocity v += a * dt

position p += v * dt

All mks units

=end

class OrbitalPhysics

   # The all-important force of gravity

   OrbitalPhysics::G = 6.6742e-11 # N m^2 / kg^-2

   def OrbitalPhysics::compute_acceleration(pa,pb,dt)
      # don't compute self-gravitation
      if(pa != pb)
         radius = pa.pos - pb.pos
         sumsq = radius.sumsq
         hypot = Math.sqrt(sumsq)
         acc = -G * pb.mass / sumsq
         # this line assigns the acceleration to
         # the x,y,z velocity components along the
         # radius pa - pb without using trig functions
         pa.vel += radius / hypot * acc * dt
      end
   end

   def OrbitalPhysics::process_planets(planet_list,dt,symmetric = false)
      if(symmetric)
         # compute gravitation interactively for all bodies
         # very slow ... requires p^2 time
         planet_list.each do |p1|
            planet_list.each do |p2|
               compute_acceleration(p1,p2,dt)
            end
            p1.pos += p1.vel * dt
         end
      else
         # compute gravitation only wrt the sun
         # which is assumed to be first member of array
         sun = planet_list.first
         planet_list.each do |planet|
            compute_acceleration(planet,sun,dt)
            planet.pos += planet.vel * dt
         end
      end
   end

end # class OrbitalPhysics

# Main program
if __FILE__ == $0
   # Set values as your own application.
   PROG_PATH = "gravity.glade"
   PROG_NAME = "Gravity"
   Gravity.new(PROG_PATH, nil, PROG_NAME)
   Gtk.main
end