When I run the following script below, I get the error Too many input arguments (I have figured out that this specifically happens when calculating intNH(5,2) in the second for loop). Basically what is happening is that fun2 for i=5 equals 0 (because the relevant G element is 0) and MATLAB cannot numerically integrate a function equal to zero (or that's how I've interpreted this).
Does anyone have any pointers about how I can overcome this? I have tried some if statements saying that if the function is equal to 0 then intNH = 0, else numerically integrate as normal, but more errors keep happening when I try this!
%%Calculation of dH/dt for mode m=1 for the entire sphere, NH and SH
clear all
%%Radius of photosphere
R = 6.957*(10^5); %In km
%%Call in spherical harmonic coefficients, change the 535 figure as more
%%data is added to the spreadsheets
G(:,1) = xlsread('G Coefficients.xls', 'C3:C535');
G(:,2) = xlsread('G Coefficients.xls', 'E3:E535');
G(:,3) = xlsread('G Coefficients.xls', 'H3:H535');
G(:,4) = xlsread('G Coefficients.xls', 'L3:L535');
G(:,5) = xlsread('G Coefficients.xls', 'Q3:Q535');
G(:,6) = xlsread('G Coefficients.xls', 'W3:W535');
G(:,7) = xlsread('G Coefficients.xls', 'AD3:AD535');
G(:,8) = xlsread('G Coefficients.xls', 'AL3:AL535');
G(:,9) = xlsread('G Coefficients.xls', 'AU3:AU535');
%%Set function v which always remains the same
nhztoradperday = 2*pi*86400*(10^(-9));
a = 460.7*nhztoradperday;
b = -62.69*nhztoradperday;
c = -67.13*nhztoradperday;
syms t %t short for theta here
V = a + (b*cos(t)^2) + (c*cos(t)^4);
zero = @(t) 0*t;
%%Set B and A functions within separate matrices
for i = 1:length(G)
B(i,1) = G(i,1)*cos(t);
B(i,2) = 0.5*G(i,2)*(3*(cos(t)^2)-1);
B(i,3) = 0.5*G(i,3)*(5*(cos(t)^3)-3*cos(t));
B(i,4) = 1/8*G(i,4)*(35*(cos(t)^4)-30*(cos(t)^2)+3);
B(i,5) = 1/8*G(i,5)*(63*(cos(t)^5)-70*(cos(t)^3)+15*cos(t));
B(i,6) = 1/16*G(i,6)*(231*(cos(t)^6)-315*(cos(t)^4)+105*(cos(t)^2)-5);
B(i,7) = 1/16*G(i,7)*(429*(cos(t)^7)-693*(cos(t)^5)+315*(cos(t)^3)-35*cos(t));
B(i,8) = 1/128*G(i,8)*(6435*(cos(t)^8)-12012*(cos(t)^6)+6930*(cos(t)^4)-1260*(cos(t)^2)+35);
B(i,9) = 1/128*G(i,9)*(12155*(cos(t)^9)-25740*(cos(t)^7)+18018*(cos(t)^5)-4620*(cos(t)^3)+315*cos(t));
A(i,1) = 0.5*R*G(i,1)*sin(t);
A(i,2) = 0.5*R*G(i,2)*csc(t)*(cos(t)-(cos(t)^3));
A(i,3) = 0.5*R*G(i,3)*csc(t)*(1.5*(cos(t)^2)-1.25*(cos(t)^4)-0.25);
A(i,4) = 1/8*R*G(i,4)*csc(t)*(-7*(cos(t)^5)+10*(cos(t)^3)-3*cos(t));
A(i,5) = 1/8*R*G(i,5)*csc(t)*(-21/2*(cos(t)^6)+35/2*(cos(t)^4)-15/2*(cos(t)^2)+0.5);
A(i,6) = 1/16*R*G(i,6)*csc(t)*(-33*(cos(t)^7)+63*(cos(t)^5)-35*(cos(t)^3)+5*cos(t));
A(i,7) = 1/16*R*G(i,7)*csc(t)*(-429/8*(cos(t)^8)+231/2*(cos(t)^6)-315/4*(cos(t)^4)+35/2*(cos(t)^2)-5/8);
A(i,8) = 1/128*R*G(i,8)*csc(t)*(-715*(cos(t)^9)+1716*(cos(t)^7)-1386*(cos(t)^5)+420*(cos(t)^3)-35*cos(t));
A(i,9) = 1/128*R*G(i,9)*csc(t)*(-2431/2*(cos(t)^10)+6435/2*(cos(t)^8)-3003*(cos(t)^6)+1155*(cos(t)^4)-315/2*(cos(t)^2)+7/2);
end
%Turn vector V from a sym to a function and compute V at equator
Vfunc = matlabFunction(V);
Veq = Vfunc(0);
for i=1:length(G)
fun(1) = 0; fun(2) = 0; fun(3) = 0; fun(4) = 0; fun(5) = 0; fun(6) = 0;
fun(7) = 0; fun(8) = 0; fun(9) = 0;
%Create functions for integration
fun(1) = matlabFunction(B(i,1).*A(i,1).*(V-Veq).*sin(t));
fun(2) = matlabFunction(B(i,2).*A(i,2).*(V-Veq).*sin(t));
fun(3) = matlabFunction(B(i,3).*A(i,3).*(V-Veq).*sin(t));
fun(4) = matlabFunction(B(i,4).*A(i,4).*(V-Veq).*sin(t));
fun(5) = matlabFunction(B(i,5).*A(i,5).*(V-Veq).*sin(t));
fun(6) = matlabFunction(B(i,6).*A(i,6).*(V-Veq).*sin(t));
fun(7) = matlabFunction(B(i,7).*A(i,7).*(V-Veq).*sin(t));
fun(8) = matlabFunction(B(i,8).*A(i,8).*(V-Veq).*sin(t));
fun(9) = matlabFunction(B(i,9).*A(i,9).*(V-Veq).*sin(t));
%Compute numerical integral for each order and store values
%Northern Hemisphere
intNH(i,1) = integral(fun1,0,pi/2);
intNH(i,2) = integral(fun2,0,pi/2);
intNH(i,3) = integral(fun3,0,pi/2);
intNH(i,4) = integral(fun4,0,pi/2);
intNH(i,5) = integral(fun5,0,pi/2);
intNH(i,6) = integral(fun6,0,pi/2);
intNH(i,7) = integral(fun7,0,pi/2);
intNH(i,8) = integral(fun8,0,pi/2);
intNH(i,9) = integral(fun9,0,pi/2);
%Southern Hemisphere
intSH(i,1) = int(fun1,t,pi/2,pi);
intSH(i,2) = int(fun2,t,pi/2,pi);
intSH(i,3) = int(fun3,t,pi/2,pi);
intSH(i,4) = int(fun4,t,pi/2,pi);
intSH(i,5) = int(fun5,t,pi/2,pi);
intSH(i,6) = int(fun6,t,pi/2,pi);
intSH(i,7) = int(fun7,t,pi/2,pi);
intSH(i,8) = int(fun8,t,pi/2,pi);
intSH(i,9) = int(fun9,t,pi/2,pi);
%Whole Sun
intSun(i,1) = int(fun1,t,0,pi);
intSun(i,2) = int(fun2,t,0,pi);
intSun(i,3) = int(fun3,t,0,pi);
intSun(i,4) = int(fun4,t,0,pi);
intSun(i,5) = int(fun5,t,0,pi);
intSun(i,6) = int(fun6,t,0,pi);
intSun(i,7) = int(fun7,t,0,pi);
intSun(i,8) = int(fun8,t,0,pi);
intSun(i,9) = int(fun9,t,0,pi);
%Sum over all orders to get a value for dH/dt
NH(i) = sum(dintNH(i,:));
SH(i) = sum(intSH(i,:));
Sun(i) = sum(intSun(i,:));
end
x = linspace(1643,2175,533);
figure(1)
plot(x,NH)
xlabel('Carrington Rotation')
ylabel('Magnetic helicity transport')
title('Northern Hemisphere magnetic helicity transport via twisting motions on the boundary vs Carrington Rotations')
figure(2)
plot(x,SH)
xlabel('Carrington Rotation')
ylabel('Magnetic helicity transport')
title('Southern Hemisphere magnetic helicity transport via twisting motions on the boundary vs Carrington Rotations')
figure(3)
plot(x,Sun)
xlabel('Carrington Rotation')
ylabel('Magnetic helicity transport')
title('Whole sun magnetic helicity transport via twisting motions on the boundary vs Carrington Rotations')
There are many problems with your code. All this magic you're doing with syms and setting function-valued arrays manually is entirely unnecessary, and it's hard to tell where the exact problem is coming from (if you put this code into a function, at least we'd know the line where the error comes from...).
You should just use anonymous functions, just like the one called zero in your code. All you have are analytically defined functions and numbers. You can only store function (handles) in cell objects, though, not in arrays. So fun{k} and friends should be assigned for this. But if you rewrite your code to work with numerical values (which is actually optimal for MATLAB too), you'll end up a happier person.
For instance, your first loop is equivalent to this:
B{1} = @(t) G(:,1)*cos(t);
B{2} = @(t) 0.5*G(:,2)*(3*(cos(t)^2)-1);
B{3} = @(t) 0.5*G(:,3)*(5*(cos(t)^3)-3*cos(t));
B{4} = @(t) 1/8*G(:,4)*(35*(cos(t)^4)-30*(cos(t)^2)+3);
B{5} = @(t) 1/8*G(:,5)*(63*(cos(t)^5)-70*(cos(t)^3)+15*cos(t));
B{6} = @(t) 1/16*G(:,6)*(231*(cos(t)^6)-315*(cos(t)^4)+105*(cos(t)^2)-5);
B{7} = @(t) 1/16*G(:,7)*(429*(cos(t)^7)-693*(cos(t)^5)+315*(cos(t)^3)-35*cos(t));
B{8} = @(t) 1/128*G(:,8)*(6435*(cos(t)^8)-12012*(cos(t)^6)+6930*(cos(t)^4)-1260*(cos(t)^2)+35);
B{9} = @(t) 1/128*G(:,9)*(12155*(cos(t)^9)-25740*(cos(t)^7)+18018*(cos(t)^5)-4620*(cos(t)^3)+315*cos(t));
A{1} = @(t) 0.5*R*G(:,1)*sin(t);
A{2} = @(t) 0.5*R*G(:,2)*csc(t)*(cos(t)-(cos(t)^3));
A{3} = @(t) 0.5*R*G(:,3)*csc(t)*(1.5*(cos(t)^2)-1.25*(cos(t)^4)-0.25);
A{4} = @(t) 1/8*R*G(:,4)*csc(t)*(-7*(cos(t)^5)+10*(cos(t)^3)-3*cos(t));
A{5} = @(t) 1/8*R*G(:,5)*csc(t)*(-21/2*(cos(t)^6)+35/2*(cos(t)^4)-15/2*(cos(t)^2)+0.5);
A{6} = @(t) 1/16*R*G(:,6)*csc(t)*(-33*(cos(t)^7)+63*(cos(t)^5)-35*(cos(t)^3)+5*cos(t));
A{7} = @(t) 1/16*R*G(:,7)*csc(t)*(-429/8*(cos(t)^8)+231/2*(cos(t)^6)-315/4*(cos(t)^4)+35/2*(cos(t)^2)-5/8);
A{8} = @(t) 1/128*R*G(:,8)*csc(t)*(-715*(cos(t)^9)+1716*(cos(t)^7)-1386*(cos(t)^5)+420*(cos(t)^3)-35*cos(t));
A{9} = @(t) 1/128*R*G(:,9)*csc(t)*(-2431/2*(cos(t)^10)+6435/2*(cos(t)^8)-3003*(cos(t)^6)+1155*(cos(t)^4)-315/2*(cos(t)^2)+7/2);
These are all cells containing array-valued function handles. Each of these functions gets a single value t and outputs an array of length length(G).
Later, you can do
V = @(t) a + (b*cos(t)^2) + (c*cos(t)^4);
Veq = V(0);
% todo: pre-allocate fun
intNH = zeros(length(G),9);
for k=1:9
fun{k} = @(t) B{k}(t).*A{k}(t).*(V(t)-Veq).*sin(t);
intNH(:,k) = integral(fun{k},0,pi/2,'arrayvalued',true);
end
and so on.