using Winston

function doPlot(V_in, k_l, k_p, K_m)
    
    dt    = 0.2
    tlast = 1000
    
    iterations = integer(round(tlast/dt))
    ATPall = zeros(Float64, (iterations, 1))
    Gall = zeros(Float64, (iterations, 1))

    ATP = 4
    G = 3
    for i = 1:iterations
        ATPall[i] = ATP
        Gall[i] = G

        dATPdt = 2 * k_l * G * ATP - (k_p * ATP) / (ATP + K_m)
        dGdt = V_in - k_l * G * ATP

        ATP = ATP + dATPdt * dt
        G = G + dGdt * dt
    end

    tall = dt*(0:iterations-1)'
    figure()
    hold(true)
    plot(tall, ATPall, color="red")
    plot(tall, Gall, color="blue")
    hold(false)
    xlabel("Time")

    println("max(G) = ", maximum(Gall))
    println("max(ATP) = ", maximum(ATPall))
    # max(G) = 20.97673324771467
    # max(ATP) = 18.65978074314273
    title("")
end

doPlot(0.36, 0.02, 6, 12.0)

function doXY(V_in, k_l, k_p, K_m)
    
    dt    = 0.2
    tlast = 1000
    
    iterations = integer(round(tlast/dt))
    ATPall = zeros(Float64, (iterations, 1))
    Gall = zeros(Float64, (iterations, 1))

    ATP = 4
    G = 3
    for i = 1:iterations
        ATPall[i] = ATP
        Gall[i] = G

        dATPdt = 2 * k_l * G * ATP - (k_p * ATP) / (ATP + K_m)
        dGdt = V_in - k_l * G * ATP

        ATP = ATP + dATPdt * dt
        G = G + dGdt * dt
    end

    tall = dt*(0:iterations-1)'
    figure()
    hold(true)
    plot(Gall, ATPall, color="red")
    hold(false)
    xlabel("G")
    ylabel("ATP")
end

function testStability(V_in, k_l, k_p, K_m)
    
    dt    = 0.05
    tlast = 2000
    
    iterations = integer(round(tlast/dt))
    ATPall = zeros(Float64, (iterations, 1))
    Gall = zeros(Float64, (iterations, 1))

    ATP = 4
    G = 3
    for i = 1:iterations
        ATPall[i] = ATP
        Gall[i] = G

        dATPdt = 2 * k_l * G * ATP - (k_p * ATP) / (ATP + K_m)
        dGdt = V_in - k_l * G * ATP

        ATP = ATP + dATPdt * dt
        G = G + dGdt * dt
    end

    #deltaATP = maximum(ATPall[length(ATPall)/2:end]) - minimum(ATPall[length(ATPall)/2:end])
    #deltaG = maximum(Gall[length(Gall)/2:end]) - minimum(Gall[length(Gall)/2:end])
    #return (deltaATP, deltaG)
    return (maximum(ATPall[length(ATPall)/2:end]),
            minimum(ATPall[length(ATPall)/2:end]),
            maximum(Gall[length(Gall)/2:end]),
            minimum(Gall[length(Gall)/2:end]), )
end

function doBistable()
    x = 0.1:0.05:1.6
    varATP = zeros(length(x))
    varG = zeros(length(x))
    ATP_h = zeros(length(x))
    ATP_l = zeros(length(x))
    G_h = zeros(length(x))
    G_l = zeros(length(x))
    for i=1:length(x)
        (a_h, a_l, g_h, g_l) = testStability(x[i], 0.02, 6, 12.0)
        varATP[i] = a_h - a_l
        varG[i] = g_h - g_l
        ATP_h[i] = a_h
        ATP_l[i] = a_l
        G_h[i] = g_h
        G_l[i] = g_l
        println((a_h, g_h))
    end
    figure()
    hold(true)
    plot(x, varG, color="blue")
    plot(x, varATP, color="red")
    xlabel("V_in")
    ylabel("Stability")
    hold(false)
    title("Stability plot")

    figure()
    hold(true)
    plot(x, G_h, color="blue")
    plot(x, G_l, color="blue")
    plot(x, ATP_h, color="red")
    plot(x, ATP_l, color="red")
    xlabel("V_in")
    ylabel("Value")
    hold(false)
    title("Values plot")
end