diff --git a/ME_2046/HW4/Figure_1.png b/ME_2046/HW4/Figure_1.png new file mode 100644 index 0000000..9a05c53 Binary files /dev/null and b/ME_2046/HW4/Figure_1.png differ diff --git a/ME_2046/HW4/S25HW4.pdf b/ME_2046/HW4/S25HW4.pdf new file mode 100644 index 0000000..2792b80 Binary files /dev/null and b/ME_2046/HW4/S25HW4.pdf differ diff --git a/ME_2046/HW4/hw4.py b/ME_2046/HW4/hw4.py new file mode 100644 index 0000000..0707180 --- /dev/null +++ b/ME_2046/HW4/hw4.py @@ -0,0 +1,130 @@ +import sympy as sm +import numpy as np +import matplotlib.pyplot as plt + +sm.init_printing() + +print("PROBLEM 1:") +print("part a:") +s, z, t = sm.symbols("s z t") + +J, B, RHO, RADIUS, L, RESIST, T = sm.symbols( + "J, B, rho, r, L, R, T", Real=True, Positive=True +) +K_T, K_B = sm.symbols("K_T, K_b") + +# loop gain +Loop_Gain = K_T * 1 / (J * s + B) + +# transfer function from motor current command to linear displacement (ignoring feedback) +X_over_I = Loop_Gain * RHO / s * RADIUS +sm.pprint(X_over_I) +X_over_I = X_over_I.expand().simplify() +sm.pprint(X_over_I) + +print("part b:") +# recall ZOH_eq +# G(z) = (1-z**-1) Z{L**-1{G(s)/s}} +print("Finding inverse laplace of G/s") +G_s = X_over_I +sm.pprint(G_s / s) +g_t = sm.inverse_laplace_transform(G_s / s, s, t) +sm.pprint(g_t.expand()) + +# make z substitution subs + +G_z = (1 - z**-1) * ( + K_T * T * RADIUS * RHO / B * (z / (z - 1) ** 2) + - J * K_T * RADIUS * RHO / B**2 + + J * K_T * RADIUS * RHO / B**2 * (z / (z - sm.exp(-B / J * T))) +) + +sm.pprint(G_z) + +# part c was on the board +print("part c") +F = sm.Matrix([[0, RHO], [0, -B / J]]) +G = sm.Matrix([0, K_T / J]) + +# part d +print("part d") +A = sm.exp(F * T) +B = sm.integrate(A, (T, 0, T)) +sm.pprint(F) +sm.pprint(G) +sm.pprint(A) +sm.pprint(B) +print(sm.latex(A)) +print(sm.latex(B)) + +########################################################a +print("PROBLEM 4:") +print("bilinear transform") +T = 0.2 # s + +G_s = 10 / (s + 10) +G_z = G_s.subs(s, 2 * (z - 1) / T / (z + 1)) +sm.pprint(G_z) +sm.pprint(G_z.simplify()) + +print("bilinear with warping") +omega_star = 10 # rad/s +G_z_warped = G_s.subs(s, (omega_star / np.tan(omega_star * T / 2)) * (z - 1) / (z + 1)) +sm.pprint(G_z_warped) +sm.pprint(G_z_warped.simplify()) + +# Create a frequency vector from 0 to 10 rad/s (avoid exactly 0 to prevent log(0) issues) +omega = np.linspace(0, 10, 500) + +# Continuous-time frequency response (evaluate at s = j*omega) +G_s_response = 10 / (10 + 1j * omega) + +# For the discrete transfer functions, substitute z = exp(j*omega*T) +z_vals = np.exp(1j * omega * T) + +# Convert symbolic expressions into functions that accept numpy arrays +G_z_func = sm.lambdify(z, G_z, "numpy") +G_z_warped_func = sm.lambdify(z, G_z_warped, "numpy") + +G_z_response = G_z_func(z_vals) +G_z_warped_response = G_z_warped_func(z_vals) + + +# Helper function to compute magnitude (in dB) and phase (in degrees) +def compute_mag_phase(response): + mag = 20 * np.log10(np.abs(response)) + phase = np.angle(response, deg=True) + return mag, phase + + +# Compute responses +mag_s, phase_s = compute_mag_phase(G_s_response) +mag_z, phase_z = compute_mag_phase(G_z_response) +mag_z_warped, phase_z_warped = compute_mag_phase(G_z_warped_response) + +# Plot the Bode plots +plt.figure(figsize=(10, 8)) + +# Magnitude plot +plt.subplot(2, 1, 1) +plt.semilogx(omega, mag_s, label="Continuous: G(s)") +plt.semilogx(omega, mag_z, label="Bilinear: G(z)") +plt.semilogx(omega, mag_z_warped, label="Warped: G(z)_warped") +plt.title("Bode Magnitude Plot") +plt.ylabel("Magnitude (dB)") +plt.legend() +plt.grid(True) + +# Phase plot +plt.subplot(2, 1, 2) +plt.semilogx(omega, phase_s, label="Continuous: G(s)") +plt.semilogx(omega, phase_z, label="Bilinear: G(z)") +plt.semilogx(omega, phase_z_warped, label="Warped: G(z)_warped") +plt.title("Bode Phase Plot") +plt.xlabel("Frequency (rad/s)") +plt.ylabel("Phase (degrees)") +plt.legend() +plt.grid(True) + +plt.tight_layout() +plt.show() diff --git a/ME_2046/Project/ME2046 - Project Description.pdf b/ME_2046/Project/ME2046 - Project Description.pdf new file mode 100644 index 0000000..2400c42 Binary files /dev/null and b/ME_2046/Project/ME2046 - Project Description.pdf differ diff --git a/NUCE_2113/lab6/exercise6_4.png b/NUCE_2113/lab6/exercise6_4.png new file mode 100644 index 0000000..e73b5c4 Binary files /dev/null and b/NUCE_2113/lab6/exercise6_4.png differ diff --git a/NUCE_2113/lab6/quick_maths.py b/NUCE_2113/lab6/quick_maths.py new file mode 100644 index 0000000..ae499a2 --- /dev/null +++ b/NUCE_2113/lab6/quick_maths.py @@ -0,0 +1,61 @@ +print("Experiment 6.1\n") +measured_activity = 5.2e-6 # curies +elapsed_time = 193.92 # s +sigma_u1 = 35331 # counts +background_time = 193.94 # s +sigma_b = 558 # counts + +R = sigma_u1 / elapsed_time - sigma_b / background_time +print(f"R: {R:.3e} counts / s") +print(f"R: {R/3.7e10*1e6:.3e} microcurie") +eps_ip = 0.101 # from figure 6.2 +f = 0.6617 + +d = 100 # mm +r = 38 / 2 # mm + +G = r**2 / 4 / d**2 + +A = R / eps_ip / G / f + +print(f"G: {G:.3e}") +print(f"A: {A:.3e} decays / s") +print(f"A: {A/3.7e10*1e6:.3e} microcurie") + +A = 4.5 * 3.7e10 / 1e6 +eps_ip_est = R / A / G / f +print(f"eps_ips: {eps_ip_est:.3e}") + +sigma_u1 = [10122, 11763, 14464, 17476, 22073, 28268, 38084] +sigma_b = 130 +d = [100, 90, 80, 70, 60, 50, 40] +t = 60 # s + +import numpy as np + +n = len(sigma_u1) +R = np.empty(n) +G = np.empty(n) +eps_ip = np.empty(n) + +for i, value in enumerate(sigma_u1): + print(i, value) + R[i] = (value - sigma_b) / t + G[i] = r**2 / 4 / d[i] ** 2 + eps_ip[i] = R[i] / A / G[i] / f + +print("R") +print(R.transpose()) +print("G") +print(G.transpose()) +print("eps_ip") +print(eps_ip.transpose()) + +import matplotlib.pyplot as plt + +plt.plot(d, eps_ip, ".b") +plt.xlabel("Distance [mm]") +plt.ylabel("$\epsilon_{ip}$") +plt.grid("both") +plt.savefig("exercise6_4.png") +plt.show()