Notes/Prog/4_rb

Kepler

# генерирует html-таблицу координат x,y тела, брошенного под углом к горизонту в центральном поле.
# usage: ruby kepler.rb > yo-ho.html.
# таблицу нужно открывать calc-ом и строить диаграмму.

# at home
view_wth=1366	# a viewport width!
view_hth=683	# a viewport height!

# at work
#view_wth=1905	# a viewport width!
#view_hth=977	# a viewport height!

# -- initialization values block --------------------

  v0=848.6	# initializing velocity, m/s;
  angle=20.95	# a angle of the begin of the traectory, deg;
  r0=500000.0	# a planet radius, m;
  r1_0=300000.0	# a satellite radius, m;
  x1_0=20*r0	# satellite coordinates in 'r0' measure
  y1_0=-1*r0

  density=5500	# a planet density, kg/m^3;
  GM=6.67*4/3*(Math::PI)*r0**3*density/100000000000
  GM1=6.67*4/3*(Math::PI)*r1_0**3*density/100000000000
  
  dt=0.4	# time step, s;
  N=1227500	# a number of a steps;

  v0x=v0*(Math.cos(angle*Math::PI/180))
  v0y=v0*(Math.sin(angle*Math::PI/180))
  x=0
  y=r0
  r=r0
  r1=Math.sqrt((x1_0-x)**2+(y1_0-y)**2)
  ax=0 + GM1*x1_0/r1**3
  ay=-GM/r0**2 + GM1*(y1_0-r0)/r1**3

  xmax=x
  xmin=x
  ymax=y
  ymin=y

  xlist=[x]
  ylist=[y]



# -- END of initialization values block --------------------

N.times do


  vx=v0x+ax*dt
  vy=v0y+ay*dt

  x=x+(vx+v0x)*dt/2
  y=y+(vy+v0y)*dt/2

  v0x=vx
  v0y=vy

  r=Math.sqrt(x**2+y**2)
  r1=Math.sqrt((x1_0-x)**2+(y1_0-y)**2)
  ax= -GM*x/r**3 + GM1*(x1_0-x)/r1**3
  ay= -GM*y/r**3 + GM1*(y1_0-y)/r1**3

  xlist.push x		# list of x-coordinates
  ylist.push y

# min-max calculation for traectory drowing
  if x > xmax
	xmax = x end
  if x < xmin
	xmin = x end
  if y > ymax
	ymax = y end
  if y < ymin
	ymin = y end

end	# -- END of general solving block --------------

# picture scaling 

padding=10

  view_wth=view_wth-2*padding
  view_hth=view_hth-2*padding

  xscale=view_wth/(xmax-xmin)
  yscale=view_hth/(ymax-ymin)

  if xscale < yscale
	  scale = xscale
    else  scale = yscale
  end

  xcenter=-xmin*scale + padding
  ycenter=ymax*scale + padding
  x1center=(x1_0-xmin)*scale + padding
  y1center=(ymax-y1_0)*scale + padding

# -- svg-output block ---------------

svg_output=Proc.new do
  
  puts "<!-- it's generated by the 'kepler.rb' -->"
  puts '<style> svg{background:#bbb; width:100%; height: 100%; stroke:black; fill:white}
   *{margin:0}</style>'
  puts
  puts '<svg><polyline points="'
  i=0
  while i<N+1 do
    xx=(xlist[i]-xmin)*scale + padding
    yy=(ymax-ylist[i])*scale + padding
    puts xx.to_s + ',' + yy.to_s + ' '
    i+=1
  end
  puts '"/>
	<circle cx=' + xcenter.to_s + ' cy=' + ycenter.to_s + 
	' r=' + (r0*scale).to_s + ' fill=rgba(0,0,200,0.2) stroke=#353 stroke-width="3"/>
	<circle cx=' + x1center.to_s + ' cy=' + y1center.to_s + 
	' r=' + (r1_0*scale).to_s + ' fill=rgba(100,50,0,0.2) stroke=#533 stroke-width="3"/>
	<line x1="' + padding.to_s + '" y1="' + ycenter.to_s + 
	'" x2="' + (view_wth+padding).to_s + '" y2="' + ycenter.to_s + '" style="stroke:#060;"/>
	<line y1="' + padding.to_s + '" x1="' + xcenter.to_s + 
	'" y2="' + (view_hth+padding).to_s + '" x2="' + xcenter.to_s + '" style="stroke:#060;"/>'
end	# -- END of svg-output block --------------

# -- html-output block --------------

html_output=Proc.new do
  puts "<!-- it's generated by the 'kepler.rb' -->"
  puts '<!-- open this document by the LibreOffice calc & drow a chart! -->'
  puts
  puts '<table>'
  i=0
  while i<N+1 do
    # short output form
    puts i.to_s + '<tr><td>' + xlist[i].to_s + '</td><td>' + ylist[i].to_s + '</td></tr>'
    i+=1
  end
  puts '</table>'
end	# -- END of html-output block --------------

svg_output.call