#  Stinger length and width calculation for shaker
#  from Harris and Bush, Journal of Sound and Vibration 334 (2015) 255–269.
#  https://doi.org/10.1016/j.jsv.2014.09.015
#

from pylab import *

E 		= 205.0 * 1e9   	# Youngs modulus in Pa 
rho 	= 7830   			# density in kg/m^3
sig_e 	= 515.0 * 1e6 		# endurance limit in Pa

M_plate = 3.31   			# mass of experimental plate.  All masses in kg.
M_connect_plate = 0.520		# mass of the connecting plate below the experimental plate
M_bearing = 1.877			# mass of the sliding part of the square air bearing
M_top_fixture = 0.229		# mass of fixture at top end of stinger 
M_bot_fixture = 0.149		# mass of fixture at bottom end of stinger
M_shaker = 0.454			# mass of the shaker armature

M 		= M_plate + M_connect_plate + M_bearing + M_top_fixture		# mass of payload supported by rod in kg
m 		= M + M_bot_fixture + M_shaker				# mass of entire payload above and below the rod

g      	= 9.8				# gravitational acceleration m/s^2

lam1sq = 22.4  				# frequency parameter for loaded clamped beam (see Harris)

###################################################################################
#    These are the proposed design limits and the operating point

gam 	= 8 * g 			# max ampltude of acceleration
fm 		= 50				# max test frequency in Hz

d_op_mm = 1.65			# operating value of d in mm   (diameter of stinger)
L_op_mm = 45.0       	# operating value of L in mm   (length of stinger)

###################################################################################

beta 	= sqrt(1 - (g/(gam+g)))

d_mm = arange(0.001, 3.0, 0.01)  # range of rod diameters in d_mm (avoiding zero)

d = d_mm / 1000.0  # rod diameters in m

d_e = sqrt((4 * M * gam)/(pi * sig_e))  # Eqn 3 of Harris minimum "endurance" diameter d > d_e
d_e_mm = 1000.0*d_e                     # d_e in mm

L_a = (E * d**2)/(16 * pi * M * fm**2)  # Eqn 6  "axial resonance"  L < L_a

L_b = (d**2/4)*sqrt((pi**3 * E)/(M*(gam+g)))  # Eqn 9 "buckling" L < L_b

L_l = (E/rho)**(1/4) * sqrt((beta*lam1sq*d)/(8*pi*fm))  # Eqn 13 "lateral resonance" L < L_l

#  put operating point into right units

d_op_m = d_op_mm/1000.0  # operating diameter in m
L_op_m  = L_op_mm/1000.0  # operating length in m

# plot a version of Harris' Fig 7a

plot(d_mm, L_a)  					# Eqn 6  "axial resonance"  L < L_a
plot(d_mm, L_b)  					# Eqn 9 "buckling" L < L_b
plot(d_mm, L_l)  					# Eqn 13 "lateral resonance" L < L_l
plot([d_e_mm, d_e_mm], [0.0,0.5])  	# endurance  diameter d > d_e

plot(d_op_mm,L_op_m, 'r+')  # plot the operating point

xlim(0.0,3.0)
ylim(0.0, 0.3)

text(0.1,0.275,'(a)', fontsize=18)

legend(('axial resonance','buckling','lateral resonance', 'endurance', 'operating point'), loc='upper right')

xlabel('diameter d (mm)')
ylabel('Length L (m)')

savefig('L_vs_d.pdf')

show()

# plot a version of Fig 7b

inv_k_L_a = (16 * L_a**3)/(3 * pi * E * d**4)  # Eqn 14
inv_k_L_b = (16 * L_b**3)/(3 * pi * E * d**4)
inv_k_L_l = (16 * L_l**3)/(3 * pi * E * d**4)

inv_k_L_op = (16 * L_op_m**3)/(3 * pi * E * d_op_m**4)  # operating point

semilogy(d_mm, inv_k_L_a)
semilogy(d_mm, inv_k_L_b)
semilogy(d_mm, inv_k_L_l)
semilogy([d_e_mm, d_e_mm], [0.00001,0.1])  # endurance

semilogy(d_op_mm,inv_k_L_op,'r+') # plot the operating point

xlim(0.0,3.0)
ylim(0.00001,0.1)

text(0.1,0.04,'(b)', fontsize=18)

xlabel('diameter d (mm)')
ylabel('1/k (m/N)')

legend(('axial resonance','buckling','lateral resonance', 'endurance', 'operating point'), loc='lower right')

savefig('invk_vs_d.pdf')

show()

# plot a version of Fig 7c

deg_rad = 180.0/pi  # convert radians into degrees

inv_kappa_L_a = deg_rad * (16 * L_a)/(pi * E * d**4)   # Eqn 15
inv_kappa_L_b = deg_rad * (16 * L_b)/(pi * E * d**4)
inv_kappa_L_l = deg_rad * (16 * L_l)/(pi * E * d**4)

inv_kappa_L_op = deg_rad * (16 * L_op_m)/(pi * E * d_op_m**4) # operating point 

semilogy(d_mm, inv_kappa_L_a)
semilogy(d_mm, inv_kappa_L_b)
semilogy(d_mm, inv_kappa_L_l)
semilogy([d_e_mm, d_e_mm], [1.0,1000.0])  # endurnace

semilogy(d_op_mm,inv_kappa_L_op,'r+') # plot the operating point

xlim(0.0,3.0)
ylim(1.0,1000.0)

text(0.06,300.0,'(c)', fontsize=18)

xlabel('diameter d (mm)')
ylabel(r'1/kappa (deg/Nm)')

legend(('axial resonance','buckling','lateral resonance', 'endurance', 'operating point'), loc='upper right')

savefig('invkappa_vs_d.pdf')

show()