Here is the latest version of the entire code, as suggested by @David Goodmanson. The value still isn't the same sadly, but seems much closer to the corect result.
Any advice would be much appreciated!
clc; clear all;
%% self inductance of the coils
%https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1669896
%all dimensions in inches
a = 1.476;
b = 0.315;
c = 0.945;
N12 = 27; %total number of turns
mu0 = 4*pi*10^(-7); %permeability of free space in [H/m]
%Wheeler's formula for multi-layer coil
%in microhenries
L_Wheeler = 0.8*a.^2*N12^2/(6*a+9*b+10*c);
disp('L_Wheeler = '); disp(L_Wheeler);
%% Mutual inductance calculations based on https://www.researchgate.net/publication/265332426_Mutual_inductance_calculation_between_misalignment_coils_for_wireless_power_transfer_of_energy
%everything in centimeters
R1 = 7.8232;
R2 = 11.7729;
a_axial = 14.2748;
b_axial = 2.413;
hp = 1.397;
hs = 4.1529;
zc = 7.366;
yc = 30.988;
xc = 0;
N1 = 1142;
N2 = 516;
N = 2;
n = N;
S = 2;
m = S;
h = 0;
q = h;
p = 0;
g = p;
Rp = @(h) R1 + sum(hp*h/(2*N+1));
Rs = @(q) R2 + sum(hs*q/(2*n+1));
alfa = @(q,h) Rs(q)/ Rp(h);
%initial coefficients of the plane of the secondary coil
a = 0;
b = 1;
c = 0.001;
l = @(a, c) sqrt(a^2+c^2);
L = @(a, b, c) sqrt(b^2 + l(a, c)^2);
x = @(p, xc) xc + b_axial*a*p/(2*m+1);
y = @(p, yc) yc + b_axial*b.*p./(2.*m+1);
z = @(zc, g, p) zc + a_axial.*g./(2.*S+1) + b_axial*c.*p./(2.*m+1);
beta = @(p, xc, h) x(p,xc)./Rp(h);
gamma = @(p, yc,h) y(p, yc)./Rp(h);
delta = @(zc, g, p, h) z(zc, g, p)./Rp(h);
phi = 0;
p1 = @(p, yc, h, a, c) gamma(p,yc, h) .*c/l(a, c);
p2 = @(p, xc, yc, h, a, b, c) (beta(p, xc, h).*l(a, c)^2+gamma(p, yc, h)*a*b)./(l(a,c)*L(a, b, c));
p3 = @(a, b, c) a*c/L(a, b, c);
p4 = @(a,b,c,xc,yc,zc,g,p,h) (gamma(p,yc,h)*l(a,c)^2-beta(p, xc, h)*a-delta(zc, g, p, h)*b*c)/(l(a,c)*L(a,b,c));
p5 = @(a,b,c,xc,zc,g,p,h) (delta(zc,g,p,h)*a-beta(p, xc, h)*c)/l(a,c);
l1 = @(a, b, c) 1 - (b^2*c^2)/(l(a,c)^2*L(a,b,c)^2);
l2 = @(a,b,c) c^2/l(a, c)^2;
l3 = @(a,b,c) a*b*c/(l(a,c)^2*L(a,b,c));
q1 = @(a,b,c,xc,yc,p,h) (gamma(p,yc,h)*l(a,c)^2-beta(p,xc,h)*a*b)/(l(a,c)*L(a,b,c));
q2 = @(a,c,xc,p,h) (-beta(p,xc,h)*c)/l(a,c);
A0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi) 1 + alfa(q,h)^2 + ...
beta(p,xc,h).^2 + gamma(p,yc,h).^2 + delta(zc,g,p,h).^2 + ...
2*alfa(q,h)*(p4(a,b,c,xc,yc,zc,g,p,h).*cos(phi) + p5(a,b,c,xc,zc,g,p,h).*sin(phi));
V0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi) sqrt(beta(p,xc,h).^2 + gamma(p,yc,h).^2 + ...
alfa(q,h)^2*(l1(a,b,c)*cos(phi).^2+l2(a,b,c)*sin(phi).^2 + l3(a,b,c)*sin(2*phi)) + ...
2*alfa(q,h)*(q1(a,b,c,xc,yc,p,h)*cos(phi) + q2(a,c,xc,p,h)*sin(phi)));
k_squared = @(a,b,c,xc,yc,zc,q,g,p,h,phi) 4*V0(a,b,c,xc,yc,zc,q,g,p,h,phi)./(A0(a,b,c,xc,yc,zc,q,g,p,h,phi)+2*V0(a,b,c,xc,yc,zc,q,g,p,h,phi));
k = sqrt(k_squared(a,b,c,xc,yc,zc,q,g,p,h,phi));
[K,E] = ellipke(k);
Psi = @(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E) (1 - k_squared(a,b,c,xc,yc,zc,q,g,p,h,phi)./2).*K-E;
int_M0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E) (p1(p, yc, h, a, c).*cos(phi) + ...
p2(p, xc, yc, h, a, b, c).*sin(phi) + p3(a, b, c)).*Psi(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E)./(k.*sqrt(V0(a,b,c,xc,yc,zc,q,g,p,h,phi).^3));
M0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E) mu0*Rs(q)./pi .* integral(@(phi) int_M0(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E), 0, 2*pi, 'ArrayValued', true);
disp('M0 = '); disp(M0(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E));
%% 1 test on example 1 from the paper
M = @(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E) N1*N2 .* M0(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E)./(2*n+1)./(2*N+1)./(2*S+1)./(2*m+1);
M_temp=0;
for i = -n:1:n
q = i;
for ii = -m:1:m
p = ii;
for iii = -N:1:N
h = iii;
for iv = -S:1:S
g = iv;
k = sqrt(abs(k_squared(a,b,c,xc,yc,zc,q,g,p,h,phi)));
[K,E] = ellipke(k);
M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
end
end
end
end
disp('M_temp = '); disp(M_temp);
disp('M = '); disp(M(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E));