require 'gsl'
require 'nyaplot'
require 'nyaplot3d'
require 'nyaplot_utils'

alpha = 10.82
beta = 14.286
a = 1.3
b = 0.1
d = 0.2
func = lambda {|x| return -b*Math.sin((Math::PI*x)/(2*a)+d)}

ode_func = Proc.new { |t, y, dydt, mu|
  dydt[0] = alpha*(y[1] - func.call(y[0]))
  dydt[1] = y[0] - y[1] + y[2]
  dydt[2] = -beta*y[1]
}

solver = GSL::Odeiv::Solver.alloc(GSL::Odeiv::Step::RK8PD, [1e-11, 0.0], ode_func, 3)

t = 0.0       # initial time
t1 = 250.0    # end time
h = 1e-11      # initial step
y = GSL::Vector.alloc([0, 0.49899, 0.2]) # initial value
mat = []

while t < t1
  t, h, status = solver.apply(t, t1, h, y)
  break if status != GSL::SUCCESS
  mat.push([t, y[0],y[1],y[2]])
end

mat = mat.transpose
df = Nyaplot::DataFrame.new({t: mat[0], x: mat[1], y: mat[2], z: mat[3]})
df

color = Nyaplot::Colors.GnBu

fill = color.to_a[4]

plot = Nyaplot::Plot.new
line = plot.add_with_df(df, :line, :x, :y)
line.color(fill)
line.stroke_width(1)
plot

plot = Nyaplot::Plot3D.new
plot.add_with_df(df, :line, :x, :y, :z)
plot

a=1.013; b=-0.021;c=0.019;d=0.96;e=0;f=0.01;g=1;u=0.05;i=0.05
ode_func = Proc.new { |t, y, dydt, mu|
  dydt[0] = (y[0]-a*y[1])*Math.cos(y[2])-b*y[1]*Math.sin(y[2])
  dydt[1] = (y[0]+c*y[1])*Math.sin(y[2])+d*y[1]*Math.cos(y[2])
  dydt[2] = e + f*y[2] + g*Math.atan(((1-u)*y[1]) / (1-i)*y[0])
}

solver = GSL::Odeiv::Solver.alloc(GSL::Odeiv::Step::RK8PD, [1e-10, 0.0], ode_func, 3)

t = 0.0
t1 = 800.0
h = 1e-10
y = GSL::Vector.alloc([0.9,1,1])
mat = []

while t < t1
  t, h, status = solver.apply(t, t1, h, y)
  break if status != GSL::SUCCESS
  mat.push([y[0],y[1],y[2]])
end

mat = mat.transpose
df = Nyaplot::DataFrame.new({x: mat[0], y: mat[1].map{|v| -v}, z:mat[2]})

color = Nyaplot::Colors.GnBu

plot = Nyaplot::Plot.new
line = plot.add_with_df(df, :line, :y, :x)
line.stroke_width(0.3)
line.color(color.to_a[5])
plot

plot.export_contents 'context.svg'

plot = Nyaplot::Plot3D.new
plot.add_with_df(df, :line, :x, :y, :z)
plot