Why the function I wrote doesnt plot ?

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

Answer 1

[Walter Roberson]

You are overwriting all of E_theta every time through the loops, so at the end it is a scalar.

[Burak]

so how can I correct it , I am new at matlab.I want to plot according to R.

[Walter Roberson]

Look on the line right above your assignment to E_theta to see how to save a value each iteration of the nested "for" loops.

[Burak]

still didn't get it , which value , by the way thanks for your concern.I am still trying to solve the problem.

[Walter Roberson]

E_theta(n,thetad) = ....

Answer 2

[Youssef Khmou]

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))

[Youssef Khmou]

NaN values ate still produced in a b and c, try to change the equation

[Burak]

thanks for help , I will try to write it another ways but the equation is exactly what it is.

[Burak]

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