Why the function I wrote doesnt plot ?

1 次查看（过去 30 天）
Burak
Burak 2013-4-6
I wrote a function that calculates E(theta) for PEC sphere about electromagnetic plane wave scattering but it does not plot.Where is my fault ?thanks for your help.
format long
tic
N_cut=20;
eps0=(10^-9)/(36*pi);
mu0=4*pi*10^-7;
epsr1=1.;
epsr2=2.;
mur1=1.;
mur2=1.;
R=linspace(0,2,N_cut);
global phid;
global thetad;
eps1=epsr1*eps0;
eps2=epsr2*eps0;
mu1=mur1*mu0;
mu2=mur2*mu0;
freq=6*10^8;
omeg=2*pi*freq;
k=omeg*sqrt(eps1*mu1);
lambda=2*pi/k;
for phid=1:360
for thetad=1:180
phi(phid)=phid*180/pi;
theta(thetad)=thetad*180/pi;
end
end
for thetad=1:180
for n=1:N_cut
bes_kur(n,:)=sqrt(pi.*k.*R./2).*besselj(n+0.5,k.*R);
han_kur(n,:)=sqrt(pi.*k.*R./2).*besselh(n+0.5,2,k.*R);
bes_kur_der(n,:)=-n.*sqrt(pi.*k./(2.*R)).*besselj(n+0.5,k.*R)+k.*sqrt(pi.*k.*R./2).*besselj(n-0.5,k.*R);
han_kur_der(n,:)=-n.*sqrt(pi.*k./(2.*R)).*besselh(n+0.5,2,k.*R)+k.*sqrt(pi.*k.*R/2).*besselh(n-0.5,2,k.*R);
a(n,:)=(1i.^n).*(2.*n+1);
b(n,:)=-((1i.^n).*(2.*n+1).*bes_kur_der(n))./han_kur_der(n);
c(n,:)=-((1i.^n).*(2.*n+1).*bes_kur(n))./han_kur(n);
L1=legendre(n,cos(theta(thetad)));
L11=legendre(n-1,cos(theta(thetad)));
L2(n,thetad)=L1(2,:);
if n==1
L3(n,thetad)=0.;
else
L3(n,thetad)=L11(2,:);
end
L2_der(n,thetad)=(-(n+1).*L3(n,thetad)+n.*cos(theta(thetad)).*L2(n,thetad))/sin(theta(thetad));
E_theta=(1i^n).*((b(n).*sin(theta(thetad))).*L2_der(n,thetad)).*(c(n).*L3(n,thetad)./sin(theta(thetad)));
end
end
F_Etheta=sum(E_theta,2);
figure
plot(R,abs(F_Etheta))
hold

请先登录，再回答此问题。

回答（2 个）

Walter Roberson
Walter Roberson 2013-4-6
You are overwriting all of E_theta every time through the loops, so at the end it is a scalar.
  4 个评论
显示 2更早的评论
Burak
Burak 2013-4-6
编辑：Burak 2013-4-7
still didn't get it , which value , by the way thanks for your concern.I am still trying to solve the problem.

请先登录，再进行评论。

Youssef Khmou
Youssef Khmou 2013-4-7
hi, The function does not plot because the data contains NaN values in both Real and Imaginary parts , so try to change the code, the NaN values occurs in the loop, also the arguments for SIN /COS should be in Radian not Degree or use 'sind' and 'cosd' instead , i tried to change the code, try to verify the vectors 'a' 'b'and 'c' :
clear all,
format long
tic
N_cut=20;
eps0=(10^-9)/(36*pi);
mu0=4*pi*10^-7;
epsr1=1.;
epsr2=2.;
mur1=1.;
mur2=1.;
R=linspace(0,2,N_cut);
eps1=epsr1*eps0;
eps2=epsr2*eps0;
mu1=mur1*mu0;
mu2=mur2*mu0;
freq=6*10^8;
omeg=2*pi*freq;
k=omeg*sqrt(eps1*mu1);
lambda=2*pi/k;
thetad=1:0.5:180;
for t=1:length(thetad)
for n=1:N_cut
% THE PROBLEM IS IN THIS LOOP BES,BES',HAN, and HAN' contain
% NaN so you cant visualize the Result
% Temporal solution : Replace NaN with ZERO.
A=sqrt(pi.*k.*R./2).*besselh(n+0.5,2,k.*R);
A(isnan(A))=0;
bes_kur(n,:)=A;
B=sqrt(pi.*k.*R./2).*besselh(n+0.5,2,k.*R);
B(isnan(B))=0;
han_kur(n,:)=B;
C=-n.*sqrt(pi.*k./(2.*R)).*besselj(n+0.5,k.*R)+k.*sqrt(pi.*k.*R./2).*besselj(n-0.5,k.*R);
C(isnan(C))=0;
bes_kur_der(n,:)=C;
D=-n.*sqrt(pi.*k./(2.*R)).*besselh(n+0.5,2,k.*R)+k.*sqrt(pi.*k.*R/2).*besselh(n-0.5,2,k.*R);
D(isnan(D))=0;
han_kur_der(n,:)=D;
a(n)=(1i.^n).*(2.*n+1);
b(n)=(-((1i.^n).*(2.*n+1).*bes_kur_der(n)))./han_kur_der(n);
c(n)=-((1i.^n).*(2.*n+1).*bes_kur(n))./han_kur(n);
L1=legendre(n,cos(thetad(t)*pi/180));
L11=legendre(n-1,cos(thetad(t)*pi/180));
L2(n,t)=L1(2,:);
if n==1
L3(n,t)=0.;
else
L3(n,t)=L11(2,:);
end
L2_der(n,t)=(-(n+1).*L3(n,t)+n.*cos(thetad(t)*pi/180).*L2(n,t))/sin(thetad(t)*pi/180);
E_theta(n)=(1i^n).*((b(n).*sin(thetad(t)*pi/180)).*L2_der(n,t)).*(c(n).*L3(n,t)./sin(thetad(t)*pi/180));
end
end
%F_Etheta=sum(E_theta,2);
figure
plot(R,abs(E_theta))
  3 个评论
显示 1更早的评论
Burak
Burak 2013-4-7
thanks for help , I will try to write it another ways but the equation is exactly what it is.
Burak
Burak 2013-4-7
编辑：Burak 2013-4-7
hi again here is the functions that I have to calculate and plot http://i48.tinypic.com/axoarn.jpg , I changed my code like this and now I got this error Subscripted assignment dimension mismatch.
Error in ==> Untitled at 55 bes_kur_der(n,mg)=C;
How can I fix that? if true
clear all,
format long
tic
N_cut=20;
eps0=(10^-9)/(36*pi);
mu0=4*pi*10^-7;
epsr1=1.;
epsr2=2.;
mur1=1.;
mur2=1.;
eps1=epsr1*eps0;
eps2=epsr2*eps0;
mu1=mur1*mu0;
mu2=mur2*mu0;
freq=6*10^8;
omeg=2*pi*freq;
k=omeg*sqrt(eps1*mu1);
lambda=2*pi/k;
Mr=21;
Mtheta=21;
Mphi=21;
thetabegin=0.;
thetaend=2*pi;
deltatheta=(thetaend-thetabegin)/Mtheta;
phibegin=0.;
phiend=pi;
deltaphi=(phiend-phibegin)/Mphi;
Rbegin=0.;
Rend=0.4;
deltaR=(Rend-Rbegin)/Mr;
R=Rbegin:deltaR:Rend;
theta=thetabegin:deltatheta:thetaend;
phi=phibegin:deltaphi:phiend;
for mg=1:22
R(mg)=R(mg);
theta(mg)=theta(mg);
phi(mg)=phi(mg);
end
for mg=1:22
for n=1:N_cut
% THE PROBLEM IS IN THIS LOOP BES,BES',HAN, and HAN' contain
% NaN so you cant visualize the Result
% Temporal solution : Replace NaN with ZERO.
A=sqrt(pi.*k.*R(mg)./2).*besselh(n+0.5,2,k.*R(mg));
A(isnan(A))=0;
bes_kur(n,mg)=A;
B=sqrt(pi.*k.*R(mg)./2).*besselh(n+0.5,2,k.*R(mg));
B(isnan(B))=0;
han_kur(n,mg)=B;
C=-n.*sqrt(pi.*k./(2.*R(mg))).*besselj(n+0.5,k.*R)+k.*sqrt(pi.*k.*R(mg)./2).*besselj(n-0.5,k.*R(mg));
C(isnan(C))=0;
bes_kur_der(n,mg)=C;
D=-n.*sqrt(pi.*k./(2.*R(mg))).*besselh(n+0.5,2,k.*R(mg))+k.*sqrt(pi.*k.*R(mg)/2).*besselh(n-0.5,2,k.*R(mg));
D(isnan(D))=0;
han_kur_der(n,mg)=D;
a(n)=(1i.^n).*(2.*n+1);
b(n)=(-((1i.^n).*(2.*n+1).*bes_kur_der(n,mg)))./han_kur_der(n,mg);
c(n)=-((1i.^n).*(2.*n+1).*bes_kur(n,mg))./han_kur(n,mg);
end
end
for mg=1:22
for n=1:N_cut
L1=legendre(n,cos(theta(mg)));
L11=legendre(n-1,cos(theta(mg)));
L0(n,mg)=L1(1,:);
L2(n,mg)=L1(2,:);
if n==1
L3(n,mg)=0.;
else
L3(n,mg)=L11(2,:);
end
%L2_der(n,mg)=n*(cos(thetag(mg))*L2(n,mg)-L0(n,mg))/sin(thetag(mg));
% L2_der(n,mg)=-(n+1)*(cos(thetag(mg))*L2(n,mg)-L0(n,mg));
L2_der(n,mg)=(-(n+1)*L3(n,mg)+n*cos(theta(mg))*L2(n,mg))/sin(theta(mg));
E_theta(n,mg)=(1i/k*R(mg))*exp(-1i*k*R(mg)*cos(phi(mg))*(1i^n).*((b(n,R(mg)).*sin(theta(mg))).*L2_der(n,mg)).*(c(n,R(mg)).*L3(n,mg)./sin(theta(mg))));
end
end
F_Etheta=sum(E_theta,2);
figure
plot(R,abs(E_theta))
% code
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Types 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by