hystersis.py (3079B)
1 import numpy as np 2 import matplotlib.pyplot as plt 3 from matplotlib import animation 4 5 fs = 48000 * 16 6 T = 1/fs # sample interval 7 M_s = 350000 # Jiles book 8 a = 2.2e4 #adjustable parameter 9 alpha = 1.6e-3 10 k = 27.0e3 #Coercivity 11 c = 1.7e-1 #1.867e-1 12 #M_r = 0.5e6 (A/m) 13 #H_c (coercivity) = 25000 - 30000 A/m (Jiles book) 14 15 # Langevin function 16 def L (x): 17 if (abs (x) > 10 ** -4): 18 return (1 / np.tanh (x)) - (1/x) 19 else: 20 return (x / 3) 21 22 # Langevin derivative 23 def L_d (x): 24 if (abs(x) > 10 ** -4): 25 return (1 / x ** 2) - (1 / np.tanh (x)) ** 2 + 1 26 else: 27 return (1 / 3) 28 29 # trapezoidal rule derivative 30 def deriv (x_n, x_n1, xDeriv_n1): 31 return ((2 / T) * (x_n - x_n1)) - xDeriv_n1 32 33 # dM/dt or "non-linear function" 34 def f (M, H, H_d): 35 Q = (H + alpha * M) / a 36 M_diff = M_s * L (Q) - M 37 delta = 1 if H_d > 0 else -1 38 delta_M = 1 if np.sign (delta) == np.sign (M_diff) else 0 39 L_prime = L_d (Q) 40 41 denominator = 1 - c * alpha * (M_s / a) * L_prime 42 43 t1_num = (1 - c) * delta_M * M_diff 44 t1_den = (1 - c) * delta * k - alpha * M_diff 45 t1 = (t1_num / t1_den) * H_d 46 47 t2 = c * (M_s / a) * H_d * L_prime 48 49 return (t1 + t2) / denominator 50 51 def M_n (M_n1, k1, k2, k3, k4): 52 return M_n1 + (k1 / 6) + (k2 / 3) + (k3 / 3) + (k4 / 6) 53 54 #input signal 55 t = np.linspace (0, 1, fs) 56 #t = 1 57 H_in = (1e5) * np.sin (2 * np.pi * 30000 * t) 58 freq = 2000 59 #H_in = np.concatenate ((5e2 * np.sin (2 * np.pi * freq * t[0:fs*5]), 1e3 * np.sin (2 * np.pi * freq * t[fs*5:fs*10]), \ 60 # 3e3 * np.sin (2 * np.pi * freq * t[fs*10:fs*15]), 5e3 * np.sin (2 * np.pi * freq * t[fs*15:fs*20]), \ 61 # 1e4 * np.sin (2 * np.pi * freq * t[fs*20:fs*25]), 3e4 * np.sin (2 * np.pi * freq * t[fs*25:fs*30]), \ 62 # 5e4 * np.sin (2 * np.pi * freq * t[fs*30:fs*35]), 1e5 * np.sin (2 * np.pi * freq * t[fs*35:fs*40]), \ 63 # 3e5 * np.sin (2 * np.pi * freq * t[fs*40:fs*45]), 5e5 * np.sin (2 * np.pi * freq * t[fs*45:fs*50]))) 64 #H_in= 5e5 * np.array ([1]) 65 # plt.plot (t, H_in) 66 # plt.show() 67 68 M_out = np.zeros (fs) 69 M_n1 = 0 70 H_n1 = 0 71 H_d_n1 = 0 72 73 n = 0 74 percent = 0 75 for H in H_in: 76 H_d = deriv (H, H_n1, H_d_n1) 77 #print (H_d) 78 79 k1 = T * f (M_n1, H_n1, H_d_n1) 80 #print (k1) 81 k2 = T * f (M_n1 + k1/2, (H + H_n1) / 2, (H_d + H_d_n1) / 2) 82 #print (k2) 83 k3 = T * f (M_n1 + k2/2, (H + H_n1) / 2, (H_d + H_d_n1) / 2) 84 #print (k3) 85 k4 = T * f (M_n1 + k3, H, H_d) 86 #print (k4) 87 88 M = M_n (M_n1, k1, k2, k3, k4) 89 #print (M) 90 91 M_n1 = M 92 H_n1 = H 93 H_d_n1 = H_d 94 95 M_out[n] = M 96 n += 1 97 98 curPercent = (int) (n / (fs * 50) * 100) 99 if (curPercent > percent + 4): 100 percent = curPercent 101 print (str (percent) + "% completed") 102 103 MH = plt.figure (1) 104 plt.plot (H_in[0:20000] / 1000, M_out[0:20000] / M_s) 105 plt.xlabel ("Magnetic Field (A/m)") 106 plt.ylabel ("Tape Magnetisation (A/m)") 107 plt.title ("Simulated Digital Tape Magnetization Hysteresis Loop") 108 #MH.show() 109 110 #Mt = plt.figure (2) 111 #plt.plot (t[0:20000], M_out[0:20000] / M_s) 112 plt.show()