Integration HELP Mie Scattering
4 次查看(过去 30 天)
显示 更早的评论
Dear friends, i am really in trable, i am trying to find the specific attenuation due to rain but i have problem in my code and i am asking you your help. My problem is in the last command which is an integration operation, and i am sure you will do your best for helping me. The link below is for the equations and it's followed by the code.
radius = 1;
nMax = 40; % maximum mode number
No=8*10^3; %m^ -4
R1=140; %rainfall rate=1mm/hr
AS1=8.2*R1^-0.21;
N1D=No*exp(-AS1*radius*1e-3);
% mode numbers
mode = 1:nMax;
frequency = 6e9;
% speed of light
c = 299792458.0;
lambda = c / ( frequency ) ;
for n=1:10;
n2 = (2*n+1);
% radian frequency
w = 2.0*pi*frequency;
% wavenumber
k = w/c;
% conversion factor between cartesian and spherical Bessel/Hankel function
s = sqrt(0.5*pi/(k*radius));
% compute spherical bessel, hankel functions
[J(mode)] = besselj(mode + 1/2, k*radius); J = J*s;
[H(mode)] = besselh(mode + 1/2, 2, k*radius); H = H*s;
[J2(mode)] = besselj(mode + 1/2 - 1, k*radius); J2 = J2*s;
[H2(mode)] = besselh(mode + 1/2 - 1, 2, k*radius); H2 = H2*s;
% derivatives of spherical bessel and hankel functions % Recurrence relationship, Abramowitz and Stegun Page 361
kaJ1P(mode) = (k*radius*J2 - mode .* J );
kaH1P(mode) = (k*radius*H2 - mode .* H );
% Ruck, et. al. (3.2-1)
An = -((1i).^mode) .* ( J ./ H ) .* (2*mode + 1) ./ (mode.*(mode + 1));
% Ruck, et. al. (3.2-2), using derivatives of bessel functions
Bn = ((1i).^(mode+1)) .* (kaJ1P ./ kaH1P) .* (2*mode + 1) ./ (mode.*(mode + 1));
Qt = (lambda^2/2*pi)*sum(n2).*sum(real(An + Bn));
end
Co1=4.343*No;
==============================================================
NOW SHOULD BE THE INTEGRATION OPERATION
0 个评论
回答(1 个)
Mike Hosea
2013-8-5
Since I don't know much about this application area, I'm not in position to work on the code here. I would approach a problem like this in stages. The first stage is to write a function that you're integrating, with all parameters being represented as inputs to that
Step 1: Write a function that takes all your parameters as inputs and produces the value of the integrand for those inputs. That is to say, if I wanted to integrate f(x,a,b)*g(x,c) over x, I would write a function with the signature fun(x,a,b,c). In this case it will be fun(r,m,lambda).
Step 2: If possible make this function "vectorized" in the integration variable. Usually this means using .* instead of *, ./ instead of /, and .^ instead of ^. You will not need to vectorize it in all the variables, as these other variables will be fixed scalars during the integration. Sometimes vectorization is not practical, in which case you can skip this step, but you will need to add "'ArrayValued',true" to the argument list of INTEGRAL in step 3.
Step 3: Define the integral as an anonymous function. It looks like the integral you want is an antiderivative, meaning that you probably want to fix the left-hand limit of integration and to consider the right-hand limit to be r, so I'll write it this way. Just realize that you could easily make the limits variable. If the integrand is fun(r,m,lambda), then you would write
kscalar = @(r,m,lambda)integral(@(r1)fun(r1,m,lambda),a,r);
If you skipped step 2, then
kscalar = @(r,m,lambda)integral(@(r1)fun(r1,m,lambda),a,r,'ArrayValued',true);
Step 4: Vectorize the integral as needed for plotting, tabulating, or whatever.
k = @(r,m,lambda)arrayfun(kscalar,r,m,lambda);
With this definition, all inputs need to be the same size, so if you wanted scalar expansion on one input, e.g. r, you'd need to write something like k(r*ones(size(m)),m,lambda) where you use k.
Not sure if that will help or not, but the bottom line is at first ignore the integral. First focus on writing a MATLAB function that accepts all the variables and parameters (and at first, only scalar values of these), and then once you have that working, move towards integrating. There's no sense in trying to integrate anything until you have coded up a function that can correctly evaluate the integrand. -- Mike
5 个评论
Mike Hosea
2013-8-5
编辑:Mike Hosea
2013-8-5
Your reference defines k as an integral, but not a definite integral. This is the so-called Fundamental Theorem of Calculus. Here k is the anti-derivative of a function of r with m and lambda as parameters. So k is a function of r with m and lambda as parameters. I am doing numerical integration, not symbolic, so the only way I can represent an anti-derivative is as a definite integral with a variable upper bound, e.g. ln(t) = integral(@(x)1./x,1,t). For example, in MATLAB
>> lns = @(t)integral(@(x)1./x,1,t);
>> lns(3)
ans =
1.0986
>> log(3)
ans =
1.0986
But this function only works for scalar inputs
>> lns(1:10)
Error using integral (line 85)
A and B must be floating point scalars.
Error in @(t)integral(@(x)1./x,1,t)
To make my function more useful in MATLAB, I need to vectorize it:
>> ln = @(t)arrayfun(lns,t);
>> ln(1:4)
ans =
0 0.69315 1.0986 1.3863
To use the kscalar or k functions, you provide the arguments to evaluate them, e.g., k(1,2,3) should return a number corresponding to the value of r=1, m=2, and lambda=3. If r should have been integrated out, then the antiderivative in your reference would need to be a definite integral with fixed limits.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Special Functions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!