Solving Dependant Integrals in MATlab

  • Context: MATLAB 
  • Thread starter Thread starter Zandman
  • Start date Start date
  • Tags Tags
    Integrals Matlab
Join the discussion
Registration is free. Ask a follow-up in this thread, or start your own.
1 reply · 3K views
Zandman
Messages
2
Reaction score
0
Hi i need to create a MATlab m file solving the following function for 0 to 90 degrees of \theta_0 and for any function F(\theta^').

[PLAIN]http://rogercortesi.com/eqn/tempimagedir/eqn9903.png

[edit] dx should be d\theta^', sorry about that.

I managed to do it in MATlab using symbolic maths toolkit and the int() function, but the script breaks for more complex F(\theta^'). I thus need ways of going about it using trapz() or any other recommendable numerical solutions. But I'm really finding it hard to wrap my head around where to start. Any ideas?

Oh just for interest the formula is used for calculating a reflector antenna's aperture efficiency. Where \theta_0 is the subtended angle or f/d ratio (Focus point), and F(\theta^') is the feed used normally in my case, (sin(x/2))^n, (sin(x))^n or 10^((nx^2)/10)
 
Last edited by a moderator:
Physics news on Phys.org
Okay I played around a bit and this is what I've ended up. But there is still problems with the amplitude. Any help will still be appreciated. The aperture efficiency for the given problem below should be between 0.8 and 0.83 but the results are out with a factor bigger than 3000. The shap is correct but the amplitude is just out by a factor.

Code:
clc
clear all

theta = 0:89;   %Angle range

for n= 2:2:8    %Specify range of n
   
    % calculate values for all angles of the integral function
    for i=1:length(theta)   
        angle = theta(i)*pi/180;
        g =(2*(n+1))*((cos(angle))^n);  %The feed 
        int(i) = sqrt(g)*tan(angle/2);  % The integral
    end
    
    %Set initial value of iterative integral sum
    res(1)=int(1);
    for j=2:length(theta)
        res(j) = res(j-1)+int(j);   %calculate and store integral sum for each case
    end

    res = abs(res);     %apply absolute value
    res = res.^2;       %square integral

    %calculate resulting values of complete function
    for k=1:length(theta) 
        res(k) = res(k)*(cot(k*pi/360))^2;
    end

plot(res)
hold on;
grid on;

end