# General solution to Part 1
t,K1,K2,w = var('t,K1,K2,w')
y = function ('y', t)
desolve ( diff (y,t ,2) + 25/2* y - 10, y)
# K2*cos(5/2*sqrt(2)*t) + K1*sin(5/2*sqrt(2)*t) + 4/5

# When y(0)=0 and y'(0)=1
w = 3.4
K1 = 4/5
K2 = (1-((0.1*w)/(w**2+25/2)))/(5*(2**(1/2))/2)


y1 = K2*cos(5/2*sqrt(2)*t) + K1*sin(5/2*sqrt(2)*t) + 4/5 + (0.1/(-w**2 + 25/2))*sin(w*t)
# plot(y1,(t,0,100))

q = 25/2
b = 0
w = 3.535
# A = 1/sqrt((q - w^2)^2 + b^2*w^2)
A = 1/sqrt((q - w^2)^2 + b^2*w^2);A
# As w approaches 5/2*sqrt(2), amplitude increases.
# For w = 1,    A = 2/23
# For w = 2,    A = 2/17
# For w = 3,    A = 2/7
# For w = 3.5,  A = 4
# For w = 3.52, A = 9.12408759124085
# For w = 4,    A = 2/7

# General solution to Part 2, varying value of b
b = 0.01
b1 = desolve ( diff (y,t ,2) + b* diff (y,t)+ 25/2* y - 10, y, ivar =t);b1

b = 2
b2 = desolve ( diff (y,t ,2) + b* diff (y,t)+ 25/2* y - 10, y, ivar =t);b2

b = 3
b3 = desolve ( diff (y,t ,2) + b* diff (y,t)+ 25/2* y - 10, y, ivar =t);b3

b = 4
b4 = desolve ( diff (y,t ,2) + b* diff (y,t)+ 25/2* y - 10, y, ivar =t);b4

b = 5
b5 = desolve ( diff (y,t ,2) + b* diff (y,t)+ 25/2* y - 10, y, ivar =t);b5

# General solution to Part 2.....
b,w,k1,k2 = var('b,w,k1,k2')
k1 = 1
k2 = 1
y_2 = k1*e^((-b - sqrt(b**2 - 50))*t/2) + k2*e^((-b - sqrt(b**2 - 50))*t/2) + 4/5 +((.1*b*w)/(-w^4-(b*w)^2 -(625/4)))*cos(w*t) + ((.1*(-w^2 + 25/2))/(w^4 + (b*w)^2+(625/4)))*sin(w*t)
# As b ->0, max amplitude increases and the value pf w where A reaches its maximum moves to the right.

# b = 2
A_1 = 1/sqrt((25/2 - w^2)^2 + (2^2)*(w^2))
# b = 1
A_2 = 1/sqrt((25/2 - w^2)^2 + (1^2)*(w^2))
# b = 0.1
A_3 = 1/sqrt((25/2 - w^2)^2 + (0.1^2)*(w^2))

plot(A_1,(w,0,5)) + plot(A_2,(w,0,5)) + plot(A_3,(w,0,5))

# General solution to Part 3 when y>=0

k1,k2 = var('k1,k2')
k1 = 0.427326508
k2 = -0.0684388889
y_3_a = k2*cos(sqrt(87)*sqrt(5)*t/5) + k1*sin(sqrt(87)*sqrt(5)*t/5) + 50/87 + (175/2452)*sin(4*t) - (5/2452)*cos(4*t)
dy1 = y_3_a.diff(t);dy1
ddy1 = y_3_a.diff(t).diff(t);ddy1
#plot(y_3_a,(t,0,100))

#General solution to part 3 when y < 0

k1,k2 = var('k1,k2')
k1 = -.028567
k2 = -.0003264879
y_3_b = (k2*cos(127/200*sqrt(31)*t) + k1*sin(127/200*sqrt(31)*t))*e^(-1/200*t) + 4/5 - 0.0285676973*sin(4*t) - 0.0003264879689*cos(4*t)
dy2 = y_3_b.diff(t);dy2
ddy2 =y_3_b.diff(t).diff(t);ddy2
#plot(y_3_b,(t,0,100))

# Search for large amplitude solution (given new initial conditions) Part 4



# If c1,c2 are equal to 175/2452 and -5/2452 "beating" occurs in the general equation for part 3 when y>=0
c1, c2 = var('c1,c2')
c1 = 175/2452
c2 = -5/2452
y_4 = c2*cos(sqrt(87)*sqrt(5)*t/5) + c1*sin(sqrt(87)*sqrt(5)*t/5) + 50/87 + (175/2452)*sin(4*t) - (5/2452)*cos(4*t)

plot(y_4,(t,0,100))

# Define a function F that takes t, y, and v values and gives you dy/dt and dv/dt from part 3.
# Make sure the code is indented as below so that Sage/Python recognizes it as being inside the function


def F(t,y):
    if y<0:
        dy = dy1
        dv = ddy1
    else:
        dy = dy2
        dv = ddy2
    return dy,dv


# Define a function eulersmethod that takes as input, a function, starting and ending times, y0, v0, and step size.
def eulersmethod(t0, tf, y0, v0, dt, F):
    # initialize time and derivative variables

    # initialize data storage variables
    tdata = [t0]
    ydata = [y0]
    vdata = [v0]
    dydata = [F(t0,y0)[0]]
    dvdata = [F(t0,y0)[1]]
    while t0 < tf:
        # update t, y, and dy/dt
        t1 = t0 + dt
        y1 = y0 + dt*F(t0,y0)[0]
        v1 = v0 + dt*F(t0,y0)[1]
        dy1 = F(t1,y1)[0]
        dv1 = F(t1,y1)[1]
        # store updated values of t, V, and dV/dt
        tdata = tdata+[t1]
        ydata = ydata+[y1]
        vdata = vdata+[v1]
        dydata = dydata+[dy1]
        dvdata = dvdata+[dv1]
        # relabel for next loop iteration
        t0 = t1
        y0 = y1
        v0 = v1
        dy0 = dy1
        dv0 = dv1

    return tdata, ydata, vdata, dydata, dvdata

# Here's an example of calling Euler's Method for t = 0 to 100, y0=1, v0=0, and step size .01. The results are placed in EM as 5 vectors - t_k, y_k, v_k, f(y_k, v_k), and g(y_k, v_k). To call any of the vectors, put [n] after EM where n = 0 to 4.  For example, EM[1] is the vector of y values.

EM = eulersmethod(0, 100, 1, 0, .1, F);EM[1]

# You should be able to graph with these vectors.  Something like the following (but you may need to play with it a little:
#P1=plot(EM[1],EM[0])
#show(P1)