Untitled

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 15 18:11:39 2024

@author: Danny
"""

import numpy as np
from numpy.polynomial import Polynomial
from numpy.polynomial.polynomial import polyval
import matplotlib.pyplot as plt


def zero_cross(ar):
    H1t=np.sign(ar)
    H1s=np.abs(H1t[:-1]-H1t[1:])
    return np.nonzero(H1s)[0]

def Power_Spectrum(t,x):
    # t is the time array x is the signal array
    T=t[-1] # T is the last element of the time array
    dt=t[1]-t[0]
    df=1/T # 1/T is the freq resolution
    X=np.fft.fft(x) # calculate the fast fourier transform of the signal
    P=np.abs(X) # find the power spectrum
    f=np.arange(0,1/dt,df )
    halfN=t.size//2 # only the first half of the power spectrum is of interest
    return f[:halfN],P[:halfN]
#input
f_in=float(input('Enter frequncy (Hz): '))


#your code
#(1)
# Function to simulate an ideal spring
def ideal(x0, v0, k):
    w0 = np.sqrt(k/m)
    A, f = [], []
    
    for x1 in x0:
        x = np.zeros_like(t)
        v = np.zeros_like(t)
        x[0], v[0] = x1, v0
        
        for i in range(1, len(t)):
            v[i] = v[i-1] - (w0**2) * x[i-1] * dt1
            x[i] = x[i-1] + v[i] * dt1
        
        plt.plot(t, x, label=f"x0 = {x1} m")
        A.append(np.max(x))
        
        zero_crossings = zero_cross(x)
        if len(zero_crossings) > 1:
            f.append(1 / (t[zero_crossings[1]] - t[zero_crossings[0]]))

    plt.title("Ideal Spring Motion")
    plt.xlabel("Time (s)")
    plt.ylabel("Position (m)")
    plt.legend()
    plt.show()

    plt.plot(A, f, '-o')
    plt.title("Amplitude vs Frequency (Ideal Spring)")
    plt.xlabel("Amplitude (m)")
    plt.ylabel("Frequency (Hz)")
    plt.show()

    f_spectrum, P_spectrum = Power_Spectrum(t, x)
    plt.plot(f_spectrum, P_spectrum)
    plt.title('Power Spectrum (Ideal Spring)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power')
    plt.grid(True)
    plt.show()
#(2)
# Function to simulate a non-ideal spring
def not_ideal(x0, v0, k):
    global poly
    w0 = np.sqrt(k/m)
    A, f = [], []

    for x1 in x0:
        x = np.zeros_like(t)
        v = np.zeros_like(t)
        x[0], v[0] = x1, v0
        
        for i in range(1, len(t)):
            v[i] = v[i-1] - (w0**2) * np.power(x[i-1], 3) * dt1
            x[i] = x[i-1] + v[i] * dt1
        
        plt.plot(t, x, label=f"x0 = {x1} m")
        A.append(np.max(x))

        zero_crossings = zero_cross(x)
        if len(zero_crossings) > 1:
            f.append(1 / (t[zero_crossings[2]] - t[zero_crossings[0]]))
    
    plt.title("Non-Ideal Spring Motion")
    plt.xlabel("Time (s)")
    plt.ylabel("Position (m)")
    plt.legend()
    plt.show()

    poly = Polynomial.fit(f, A, deg=3)
    plt.plot(A, f, '-o')
    plt.title("Amplitude vs Frequency (Non-Ideal Spring)")
    plt.xlabel("Amplitude (m)")
    plt.ylabel("Frequency (Hz)")
    plt.show()

    f_spectrum, P_spectrum = Power_Spectrum(t, x)
    plt.plot(f_spectrum, P_spectrum)
    plt.title('Power Spectrum (Non-Ideal Spring)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power')
    plt.grid(True)
    plt.show()

#parameters
dt1 = 1e-2
t = np.arange(0, 10, dt1)
K_ideal = 20
K_not_ideal = 5
m = 1
x0 = [1, 1.2, 1.4, 1.6, 1.8, 2]

# Run ideal spring and not_ideal
ideal(x0, 0, K_ideal)
not_ideal(x0, 0, K_not_ideal)

# Calculate output for given input frequency
f_in = float(input('Enter frequency (Hz): '))
poly_coef = poly.convert().coef 
A_of_f = polyval(f_in, poly_coef)
print(f'{A_of_f:.4g}')


#(3)
#parameters
x0 = float(input("Initial position: "))
v0 = float(input("Initial velocity: "))
t_input = float(input("Time: "))
dt = 1e-3
t_array = np.arange(0, 10, dt)
m = 2
k = 40
b = 5

gamma = b / (2 * m)
w0 = np.sqrt(k / m)
w = np.sqrt(w0**2 - gamma**2)

# Initialize arrays for position and velocity
v = np.zeros_like(t_array)
x = np.zeros_like(t_array)

# Initial conditions
v[0] = v0
x[0] = x0

# Numerical integration using Euler's method
for i in range(len(t_array) - 1):
    v[i + 1] = v[i] + (-2 * gamma * v[i] - w0**2 * x[i]) * dt
    x[i + 1] = x[i] + v[i + 1] * dt

# Find velocity at the input time
t_index = np.argmin(np.abs(t_array - t_input))
print(f"Velocity at t={t_input}s: {v[t_index]:.5f} m/s")

# Analytical solution
alpha = -gamma
A = x0 / np.cos(np.arctan((-x0 * gamma - v0) / (x0 * w0)))
phi = np.arctan((-x0 * gamma - v0) / (x0 * w0))
x_theory = A * np.exp(alpha * t_array) * np.cos(w0 * t_array + phi)

# Plot position as a function of time
plt.figure(figsize=(10, 6))
plt.plot(t_array, x, label='Numerical')
plt.plot(t_array, x_theory, '--', label='Analytical')
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Position vs Time')
plt.legend()
plt.grid(True)
plt.show()

# Plot Power Spectrum
f, P = Power_Spectrum(t_array, x)
plt.plot(f, P)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.title('Power Spectrum')
plt.grid(True)
plt.show()