How to integrate an array properly in matlab

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.

[Andrew Newell]

Are you trying to integrate over theta or phi?

[Alexander Lampe]

Hey, I tested your script with your values psi, phi... But then "phi" and "integrand" are both scalar. Then MATLAB expects, that the second Input off trapz() is the dimension, so it needs to be an integer. Write trapz(phi,integrand,1) or rather trapz(phi,integrand,dim) and it will work. But why you use trapz() for scalar values?

Alex

[Andrew Matthews]

Hello, thanks for the help. When I tried that earlier it did not work. Thank you for your time.

[Andrew Matthews]

Whenever I use 1 my integral becomes 0. It was not supposed to be 0. Would you know the reasoning? Thank you.

[Andrew Matthews]

I just realized, it is because my program is filling in the angles before integration. Because of this I am getting 0.

Answer 1

[Andrew Newell]

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.

[Andrew Matthews]

Hello again. This makes sense, and I realized that is why my code was not working, it was integrating an answered problem. However, I tried the way you have it specified but now it is saying that I have stuff undefined, such as the variable as. I see you actually had an answer for M. How was that obtained? I am just getting errors. Thank you for your time. Here is what I used:

function y = inductanceIntegrand(phi,psi,theta,h,ap)

% Arguments

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;

and then

psi = 0;

theta = 0:10:90;

h = 3;

ap = 6;

f = @(phi) inductanceIntegrand(phi,psi,theta,h,ap);

M = quadl(f,0,2*pi);

disp(M)

Thanks, any help would be appreciated.

[Andrew Matthews]

When I run the code the way you have it written it says 'requires more input arguments to run'. And when I have variables for input when I run the second m file it says the output and input vectors must be the same length. What do I do from here? Thank you.

[Andrew Matthews]

I just realize my error. I was trying to run to man variables at a time. Thank you all, you have been a great help.