# Code needed for drawing phase plane plots:
#
# See the web: https://cran.r-project.org/web/packages/phaseR/phaseR.pdf
#
if( !require('phaseR') ){
    install.packages('phaseR')
}
library(phaseR)

# Here is the example from the document above (showing what one can put on a plot);
#
if( FALSE ){
    logistic.flowField <- flowField(logistic, x.lim = c(0, 5), y.lim = c(-1, 3),
                                    parameters = c(1, 2), points = 21, system = "one.dim", add = FALSE)
    logistic.nullclines <- nullclines(logistic, x.lim = c(0, 5), y.lim = c(-1, 3),
                                      parameters = c(1, 2), system = "one.dim")
    logistic.trajectory <- trajectory(logistic, y0 = c(-0.5, 0.5, 1.5, 2.5), t.end = 5,
                                      parameters = c(1, 2), system = "one.dim")
}

#
# All of the functions in this set are of the same form:
#
# y' = a + b y
#
# for different values of a and b
#
my_yprime = function(t, y, parameters){
    y_prime = parameters[1] + parameters[2] * y
    list( y_prime )
}


#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_1_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( 3, -2 )

logistic.flowField <- flowField(my_yprime, x.lim = c(0, 5), y.lim = c(-1, 3),
                                parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

logistic.nullclines <- nullclines(my_yprime, x.lim = c(0, 5), y.lim = c(-1, 3),
                                  parameters = diff_eq_params, system ="one.dim")

logistic.trajectory <- trajectory(my_yprime, y0 = c(-0.5, 0.5, 2.5), t.end = 5,
                                  parameters = diff_eq_params, system = "one.dim")

#dev.off()


#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_2_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( -3, 2 )

logistic.flowField <- flowField(my_yprime, x.lim = c(0, 5), y.lim = c(-1, 3),
                                parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

logistic.nullclines <- nullclines(my_yprime, x.lim = c(0, 5), y.lim = c(-1, 3),
                                  parameters = diff_eq_params, system ="one.dim")

logistic.trajectory <- trajectory(my_yprime, y0 = c(-0.5, 0.5, 2.5), t.end = 5,
                                  parameters = diff_eq_params, system = "one.dim")

#dev.off()


#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_3_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( 3, 2 )

logistic.flowField <- flowField(my_yprime, x.lim = c(0, 5), y.lim = c(-3, 1),
                                parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

logistic.nullclines <- nullclines(my_yprime, x.lim = c(0, 5), y.lim = c(-3, 1),
                                  parameters = diff_eq_params, system ="one.dim")

logistic.trajectory <- trajectory(my_yprime, y0 = c(-2.0, -1.0, 0.5), t.end = 5,
                                  parameters = diff_eq_params, system = "one.dim")

#dev.off()


#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_4_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( -1, -2 )

logistic.flowField <- flowField(my_yprime, x.lim = c(0, 5), y.lim = c(-2, 1),
                                parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

logistic.nullclines <- nullclines(my_yprime, x.lim = c(0, 5), y.lim = c(-2, 1),
                                  parameters = diff_eq_params, system ="one.dim")

logistic.trajectory <- trajectory(my_yprime, y0 = c(-1.0, 0.0, 0.5), t.end = 5,
                                  parameters = diff_eq_params, system = "one.dim")

#dev.off()


#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_5_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( 1, 2 )

logistic.flowField <- flowField(my_yprime, x.lim = c(0, 5), y.lim = c(-2, 1),
                                parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

logistic.nullclines <- nullclines(my_yprime, x.lim = c(0, 5), y.lim = c(-2, 1),
                                  parameters = diff_eq_params, system ="one.dim")

logistic.trajectory <- trajectory(my_yprime, y0 = c(-1.0, 0.0, 0.5), t.end = 5,
                                  parameters = diff_eq_params, system = "one.dim")

#dev.off()


#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_6_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( 2, 1 )

logistic.flowField <- flowField(my_yprime, x.lim = c(0, 5), y.lim = c(-4, 0),
                                parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

logistic.nullclines <- nullclines(my_yprime, x.lim = c(0, 5), y.lim = c(-4, 0),
                                  parameters = diff_eq_params, system ="one.dim")

logistic.trajectory <- trajectory(my_yprime, y0 = c(-3.0, -1.0), t.end = 5,
                                  parameters = diff_eq_params, system = "one.dim")

#dev.off()