#/***********************************************************
#    contour.rb -- 等高線
#***********************************************************/
require "plotter.rb"

def evaluate(x,  y)
    f = x * x + 4 * y * y;  # $f(x, y) = x^2 + 4y^2$ 
    fx = 2 * x;             # $x$ で微分したもの: $f_x(x, y) = 2x$ 
    fy = 8 * y;             # $y$ で微分したもの: $f_y(x, y) = 8y$ 
    return f, fx, fy
end

def sq(x); x * x; end    # $x^2$ 
EPS = 0.001

def newton(fc, x, y)
    begin 
        f, fx, fy = evaluate(x, y)
        grad2 = sq(fx) + sq(fy)
        if (grad2 < 1e-10)
            x = 1e30;  return x, y, fx, fy, grad2
        end
        t = (fc - f) / grad2.to_f
        x += t * fx;  y += t * fy
    end while (t * t * grad2 > EPS * EPS)
    return x, y, fx, fy, grad2
end

def contour(fc,  x,  y,  step) # 等高線を描く 
    i = 0
    while (true)
        x, y, fx, fy, grad2 = newton(fc, x, y)
        if (x.abs + y.abs > 1e10); return; end
        if (i == 0);  move(x, y);  x0 = x;  y0 = y
        else; draw(x, y)
        end
        if (i > 2 && sq(x - x0) + sq(y - y0) < sq(step))
            break
        end
        t = step / Math::sqrt(grad2).to_f
        x += fy * t;  y += -fx * t
    i += 1
    end
    draw(x0, y0)
end

$gplot = IO.popen("gnuplot -persist", "w")
$gplot.puts "set multiplot"
$gplot.puts "set noautoscale"
$gplot.puts "set xrange[-3:3]"
$gplot.puts "set yrange[-2:2]"
$gplot.puts 'plot "-" w l'
while (true) 
    printf("f(x, y) = ");                      # 等高線の関数値 
    fc = gets
    if (fc == nil || fc == "\n")
        $gplot.puts "end"
        $gplot.puts "exit"
        exit 0
    end
    fc = fc.to_f
    printf("initial x = ");  x = gets.to_f # $x$ の初期値 
    printf("initial y = ");  y = gets.to_f # $y$ の初期値 
    printf("step = ");  step = gets.to_f   # ステップサイズ 
    contour(fc, x, y, step)
    $gplot.puts "end"
    $gplot.puts "replot"
end
gr_off

exit 0