print(`This is Project 1.txt`);
print(`by Rachel Adelman, Sophie Droppa, Anna Janik, and Sydney Yao`);

Help := proc()
    if nargs = 0 then
        print("The available procedures are: NBH, JAC, ORB, SSg, Beddington, EvalNBH, BeddingtonEQ, PhaseDiagE");

    elif nargs = 1 and args[1] = NBH then
        print(`NBH(r,a,K,H,P) : Nicholson-Bailey equation`);
        print(`Using parameters r, a, and K.  K is the prey carrying capacity, a is the strength of interaction between the predator and prey, and r is the reproductive rate of the prey.`);
        print(`Try:`);
        print(`NBH(1.2, 0.1, 50, 3, 3)`);
  elif nargs = 1 and args[1] = JAC then
        print(`JAC(r,a,K,H,P) : Jacobian of the Nicholson-Bailey equation`);
        print(`Using parameters r, a, and K.  K is the prey carrying capacity, a is the strength of interaction between the predator and prey, and r is the reproductive rate of the prey.`);
        print(`Try:`);
        print(`JAC(1.2, 0.1, 50, H, P)`);
elif nargs = 1 and args[1] = ORB then
        print(`ORB(f,x,x0,K1,K2)`);
        print(`Inputs a transformation F in the list of variables x with initial point pt, outputs the trajectory of the discrete dynamical system (i.e. solutions of the difference equation): x(n)=F(x(n-1)) with x(0)=x0 from n=K1 to n=K2.`);
        print(`Try:`);
        print(`ORB([5/2*x*(1-x)],[x], [0.5], 1000,1010);`);
    elif nargs = 1 and args[1] = SSg then
        print(`SSg(T,x)`);
        print(`Given a transormation T in the list of variables, finds all the steady-states (in floating point).`);
        print(`Try:`);
        print(`SSg([(1+x)/(3+y),(3+x)/(1+y)],[x,y]);`);

    elif nargs = 1 and args[1] = Beddington then
        print(`Beddington(r,a,K,N,P): Beddington equation`);
        print(`Evaluates the Beddington equating using parameters r, a, and K.  K is the prey carrying capacity, a is the strength of interaction between the predator and prey, and r is the reproductive rate of the prey.`);
        print(`Try:`);
        print(`Beddington(1.2, 0.1, 50, 3, 3);`);

    elif nargs = 1 and args[1] = EvalNBH then
        print(`EvalNBH(T,r,a,K,H,P)`);
        print(`After finding the density functions H and P using NBH, find the eigenvalues of the Jacobian and determine if the system is stable using r, a and K.`);
        print(`Try:`);
        print(`EvalNBH(NBH(1.2,0.1,50,3,3),1.2,0.1,50,H,P);`);
    elif nargs = 1 and args[1] = BeddingtonEQ then
        print(`BeddingtonEQ(r,a,K,N,P)`);
        print(`With the prey carrying capacity, K, the interaction strength, a, and the prey's reproductive rater, this function defines the Beddington equations for plotting.`);
        print(`Try:`);
        print(`BeddingtonEQ(1.2, 0.1, 50, 3, 3);`);
   elif nargs = 1 and args[1] = PhaseDiagE then
        print(`PhaseDiagE((F,x,x0,A)`);
        print(`Takes an input of a transformation F and a pair of variables x and their intial conditions, x0 and creates a phase diagram`);
        print(`Try:`);
        print(`PhaseDiagE([y,-x],[x,y],[0,1],  10);`);
    else
        print("Invalid call to Help.");
    fi;
end proc:

with(linalg):

SSg := proc(T, x)
    local eq, var, i;

    if nops(T) <> nops(x) then
        RETURN(FAIL);
    fi;

    eq := { seq(T[i] = x[i], i=1..nops(T)) };

    var := solve(eq, {op(x)});

    return evalf( [ seq( subs(var, x[i]), i=1..nops(x) ) ], 10 );
end proc:


JAC := proc(r, a, K, H, P)
    local Hnext, Pnext, dHdH, dHdP, dPdH, dPdP, J;

    Hnext := H*exp(r*(1 - H/K) - a*P);
    Pnext := a*H*(1 - exp(-a*P));

    dHdH := diff(Hnext, H);
    dHdP := diff(Hnext, P);
    dPdH := diff(Pnext, H);
    dPdP := diff(Pnext, P);

	 Matrix([[dHdH, dHdP],
                 [dPdH, dPdP]]);
end proc:

NBH := proc(r, a, K, H, P);
	[H*exp(r*(1 - H/K) - a*P), a*H*(1 - exp(-a*P))];
end proc:

EvalNBH:= proc(T,r,a,K,H,P) local N,M,L,B;
	N:= SSg(T,[H,P]);
	M:=JAC(r,a,K,H,P);
	L:= subs({H=N[1],P=N[2]},M);
	B:= evalf(eigenvalues(L));
	if abs(Re(B[1])) < 1 and abs(Re(B[2])) < 1 then
            print("the equation is stable");
        else
            print("the equation is not stable");
        fi;
	return(B);

end proc:
	

ORB:=proc(F,x,x0,K1,K2) local x1,i,L,i1,i2:

	if not (type(F,list) and type(x,list) and type(x0,list) and nops(F)=nops(x) and  nops(x)=nops(x0) and type(K1, integer) and type(K2,integer) and K1>=0 and K1<=K2) then
	 print(`bad input`):
	  RETURN(FAIL):
	fi:

	x1:=x0:

	for i from 0 to K1-1 do
	x1:=evalf([seq(subs( {seq(x[i2]=x1[i2],i2=1..nops(x))},F[i1]),i1=1..nops(F))]):
	od:

	L:=[x1]:

	for i from K1 to K2-1 do

	x1:=evalf([seq(subs( {seq(x[i2]=x1[i2],i2=1..nops(x))},F[i1]),i1=1..nops(F))]):
	 L:=[op(L),expand(x1)]:
	 od:

	L:
	evalf(L,10):

end:

Beddington := proc(r, a, K, H, P) local q, phi, sol, sol_num, l;

    q := H/K;
    phi := r*(1 - q)/(1 - exp(-r*(1 - q)));

    sol := solve(l^2 - l*(1 - r + phi) +
                 (1 - r*q)*phi + r^2*q*(1 - q) = 0, l);

    sol := [sol];
    sol_num := evalf(sol);

    if abs(Re(sol_num[1])) < 1 and abs(Re(sol_num[2])) < 1 then
        print("the equation is stable");
    else
        print("the equation is not stable");
    fi;

    return sol_num;
end proc:

BeddingtonEQ := proc(r, a, K, H, P) 
	[H*exp(r*(1 - H/K) - a*P), a*H*(1 - exp(-a*P))]; 
end proc:

PhaseDiagE:=proc(F,x,x0,A) local sol,t,i1,X,F1:

	if not (type(F,list) and type(x,list) and nops(x)=2 and type(x0,list) and nops(F)=nops(x) and nops(F)=nops(x0)  and type(A,numeric) and A>0  ) then
 		print(`bad input`):
 		RETURN(FAIL):
	fi:

	F1:=subs({seq(x[i1]=X[i1](t),i1=1..nops(x))},F):
	sol:=dsolve({seq(diff(X[i1](t),t)=F1[i1],i1=1..nops(x)), seq(X[i1](0)=x0[i1],i1=1..nops(x0))}):

	plot([subs(sol,X[1](t)),subs(sol,X[2](t)),t=0..A]):

end: