Search code examples
matlabanonymous-functionintegral

Handle function implicitly accounting for independent variables


I have that

clc, clear all, close all
tic
k1 =  1E-02:0.1:1E+02;
k2 =  1E-02:0.1:1E+02;
k3 =  1E-02:0.1:1E+02;
k = sqrt(k1.^2+k2.^2+k3.^2);
c = 1.476;
gamma = 3.9;
colors = {'b'};
Ek = (1.453*k.^4)./((1 + k.^2).^(17/6));
E = @(k) (1.453*k.^4)./((1 + k.^2).^(17/6));
E_int = zeros(1,numel(k));
E_int(1) = 1.5;

for i = 2:numel(k)
    if k(i) < 400
        E_int(i) = E_int(i-1) - integral(E,k(i-1),k(i));
    elseif k(i) > 400
        E_int(i) = 2.180/(k(i)^(2/3));
    end %end if
end %end i


beta = (c*gamma)./(k.*sqrt(E_int));
figure
plot(k,beta,colors{1})

count = 0;
%F_11 = zeros(1,numel(k1));
F_33 = zeros(1,numel(k1));

Afterwards, I should calculate F_33 as

for i = 1:numel(k1)   
    count = count + 1;    
    phi_33 = @(k2,k3) (1.453./(4.*pi)).*(((k1(i)^2+k2.^2+(k3 + beta(i).*k1(i)).^2).^2)./((k1(i)^2+k2.^2+k3.^2).^2)).*((k1(i)^2+k2.^2)./((1+k1(i)^2+k2.^2+(k3+beta(i).*k1(i)).^2).^(17/6)));
    F_33(count) = 4*integral2(phi_33,0,1000,0,1000);    
end

Now let's come to my question. I know from a paper that:

k = sqrt(k1.^2+k2.^2+k3.^2);
k30 = k3 + beta.*k1;
k0 = sqrt(k1.^2+k2.^2+k30.^2);
E_k0 = 1.453.*(k0.^4./((1+k0.^2).^(17/6)));

Therefore the expression for phi_33 would result in

phi_33 = (E_k0./(4*pi.*(k.^4))).*(k1.^2+k2.^2);

The question is: how can I make use of this final expression insted of the long one I'm using at the moment (within the for loop)?

The last expression for phi_33 is easier to handle (especially because of reckless mistakes in writing the former) and it would "pass by reference" (k2,k3), which are the independent variables.

Any hint is more than welcome.

Best regards, fpe


Solution

  • If I understand you correctly you want to use the new expression in exactly the same way as the old one-liner. You just want to divide your function phi33 into parts because of readability.

    You could do this by placing the expression in a separate function taking all values needed for the calculation. Using your old expression exactly this would look something like this:

    function phi_33 = phi_33_old(k1,k2,k3,beta,i)
    
    phi_33 = (1.453./(4.*pi)).*(((k1(i)^2+k2.^2+(k3 + beta(i).*k1(i)).^2).^2)./((k1(i)^2+k2.^2+k3.^2).^2)).*((k1(i)^2+k2.^2)./((1+k1(i)^2+k2.^2+(k3+beta(i).*k1(i)).^2).^(17/6)));
    
    end
    

    This function could then be called inside your for-loop like this.

    phi_33_test = @(k2,k3) phi_33_old(k1,k2,k3,beta,i);
    

    Using the same style, a new function could be defined as follows.

    function phi_33 = phi_33_new(k1,k2,k3,beta,i)
    
    k = sqrt(k1.^2+k2.^2+k3.^2);
    k30 = k3 + beta.*k1;
    k0 = sqrt(k1.^2+k2.^2+k30.^2);
    E_k0 = 1.453.*(k0.^4./((1+k0.^2).^(17/6)));
    phi_33_allValues = (E_k0./(4*pi.*(k.^4))).*(k1.^2+k2.^2);
    phi_33 = phi_33_allValues(i);
    
    end
    

    Note that here all values of phi_33 are calculated and then the ith value is selected. It is written in this way only to show the similarity to the old case. This new function can now be called inside the for-loop in the same way as the old one.

    phi_33 = @(k2,k3) phi_33_new(k1,k2,k3,beta,i);