MATLAB: How to integrate an array properly in matlab

arraysintegration

Hello, I am trying to finish a m-file to find the inductance in 2 coils. I have finished the program to the point of integration. I believe the problem is the integration of an array. I have tried various methods such as int, trapz, and quad but all of these seem to be returning an error. I am not sure as to whether I am implementing the commands wrong or rather I have a bad equation. Here is my code
% This program finds mutual inductance for TET coil system %
% Using Neumann's definition %
% M = sqrt(ap*as)*(1/(2*pi))*(int(((cos(theta)-(d/as)*((cos(psi)*cos(phi))-(sin(psi)*sin(phi)*cos(theta))))/(R^(3/2)))*f,phi,0,2*pi))
%
% Increment angles
%psi = [0:1:90];
%theta = 0:1:90;
%phi = 0:1:90;
psi = input('Enter psi value in degrees \n')
theta = input('Enter theta value in degrees \n')
phi = input('Enter phi value in degrees \n')
% Arguments
%dim = 0:2*pi
h = input('Enter h value in mm \n')
ap = input('Enter ap value in mm \n')
as = sqrt((ap^2)-(h^2))
delta = h/ap
alpha = as/ap
d = h.*tan(phi)
Ra = (1-(cos(phi).*cos(phi).*sin(theta).*sin(theta)))
Rb = ((2)*(d/as)).*((sin(psi).*sin(phi))-(cos(psi).*cos(phi).*cos(theta)))
Rc = ((d.^2)/(as.^2))
R = sqrt(Ra+Rb+Rc)
z = delta-(alpha*sin(theta)*cos(phi))
%kprime_2 = (((1-(alpha.*R)).^2)+z.^2)/(((1+(alpha.*R)).^2+z.^2))
kprime_2a = ((1-(alpha.*R)).^2)+z.^2
kprime_2b = ((1+(alpha.*R)).^2)+z.^2
kprime_2 = kprime_2a./kprime_2b
f = -0.011*(log(kprime_2))-0.0021
integrand = ((cos(theta)-(d/as).*((cos(psi).*cos(phi))-
(sin(psi).*sin(phi).*cos(theta))))/(R.^(3/2))).*f
%integrand1 = double(int(((cos(theta)-(d/as).*((cos(psi).*cos(phi))-(sin(psi).*sin(phi).*cos(theta))))./(R.^(3/2))).*f,phi,0,2.*pi))
%integrand = double(integrand);
stuff = trapz(phi,integrand)
M = sqrt(ap.*as).*(1/(2*pi)).*stuff
I am setting psi and phi to 0 and setting theta to 0:10:90. H is usually 3 and ap is usuall 6. This gives me different error messages for each method of integration I use.
Any help would be appreciated. Thanks.

Best Answer

Assuming you are trying to integrate over theta, a good approach is to create a function for your integrand like this:
function y = inductanceIntegrand(theta,psi,phi,h,ap)
as = sqrt((ap^2)-(h^2));
delta = h/ap;
alpha = as/ap;
d = h.*tan(phi);
Ra = (1-(cos(phi).*cos(phi).*sin(theta).*sin(theta)));
Rb = ((2)*(d/as)).*((sin(psi).*sin(phi))-(cos(psi).*cos(phi).*cos(theta)));
Rc = ((d.^2)/(as.^2));
R = sqrt(Ra+Rb+Rc);
z = delta-(alpha*sin(theta).*cos(phi));
kprime_2a = ((1-(alpha.*R)).^2)+z.^2;
kprime_2b = ((1+(alpha.*R)).^2)+z.^2;
kprime_2 = kprime_2a./kprime_2b;
f = -0.011*(log(kprime_2))-0.0021;
y = sqrt(ap.*as).*(1/(2*pi)).*((cos(theta)-(d/as).*((cos(psi).*cos(phi))-(sin(psi).*sin(phi).*cos(theta))))/(R.^(3/2))).*f;
Then you can create an anonymous function that is a parameter of theta only, and integrate:
psi = 0;
phi = 0;
h = 3;
ap = 6;
f = @(theta) inductanceIntegrand(theta,psi,phi,h,ap);
M = quadl(f,0,2*pi);
disp(M)
-4.5074e-04
If you're integrating over phi, just change the argument of f to phi and provide a value for theta.
Related Question