1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
| import numpy as np from numba import njit
@njit def compute_integral(L_y, L_z, l_0, n_y, n_z, n_mu, n_phi): y_vals = np.linspace(0, L_y, n_y)[1:-1] z_vals = np.linspace(0, L_z, n_z)[1:-1] mu_vals = np.linspace(0, 1, n_mu)[1:-1] dy = y_vals[1] - y_vals[0] dz = z_vals[1] - z_vals[0] dmu = mu_vals[1] - mu_vals[0]
integral = 0.0 for i in range(len(y_vals)): for j in range(len(z_vals)): for k in range(len(mu_vals)): phi_1 = np.arctan((y_vals[i]) / (z_vals[j])) phi_2 = np.pi / 2 + np.arctan((L_z - z_vals[j]) / (y_vals[i])) phi_3 = np.pi + np.arctan((L_y - y_vals[i]) / (L_z - z_vals[j])) phi_vals1 = np.linspace(phi_1, phi_2, n_phi)[1:-1] phi_vals2 = np.linspace(phi_2, phi_3, n_phi)[1:-1]
dphi1 = phi_vals1[1] - phi_vals1[0] dphi2 = phi_vals2[1] - phi_vals2[0] for l in range(len(phi_vals1)): term1 = np.exp(- y_vals[i] / (np.sin(phi_vals1[l]) * np.sqrt(1 - mu_vals[k]**2) * l_0)) * mu_vals[k]**2 integral += term1 * dphi1 for l in range(len(phi_vals2)): term2 = np.exp((L_z - z_vals[j]) / (np.cos(phi_vals2[l]) * np.sqrt(1 - mu_vals[k]**2) * l_0)) * mu_vals[k]**2 integral += term2 * dphi2
return integral
k = []
for L in [0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]: L_y = L L_z = L L_x = 10.0 l_0 = 1.0
n_y = 100 n_z = 100 n_mu = 100 n_phi = 100
integral = compute_integral(L_y, L_z, l_0, n_y, n_z, n_mu, n_phi) integral *= (3 / (np.pi * L_y * L_z)) * (L_y / n_y) * (L_z / n_z) * (1 / n_mu)
G_rec_inv = 1 - integral G = 1 / G_rec_inv F = 1 + 4 * l_0 / (3 * L_x)
k.append(1 / (F + G - 1))
|