Math Forum Discussions

george veropoulos

Posts: 20
Registered: 6/26/10

quadl_error
Posted: Dec 10, 2012 5:22 AM

Plain Text

Dear friends

i m try to find of integral (between 0, s =0.03 of a simpe function (name is current ). I use the quadl

when i run my code i have the message

??? Attempted to access y(13); index out of bounds because numel(y)=1.

Error in ==> quadl at 78
if ~isfinite(y(13))

Error in ==> test_propability at 4
g=quadl(@current,0,s)

ANY HELP
THANK YOU IN ADVANCE

function f=current(z)
global E0;
global L;
global s;
lambda=1;
s=0.03.*lambda;
p=0.3.*lambda;
%s is the separation distance
%p is pinth of helix
a=((s./2).^2+(p./(2.*pi)).^2).^(-1./2);
thetae=pi./2;
phip=pi./2;
thetap=pi./4;
%Lz
Lz=15.*lambda;
L=2.*pi.*Lz./(p.*a);
k=2.*pi./lambda;
theta0=3.*pi./4;
%E0 in mV
E0=1000;
%normalize factor
Einc=E0.*exp(j.*k.*cos(theta0).*Lz./2);
%lossless line
gamma=j.*k;
Zc=50;
Zs=Zc;
Zl=Zc;
%polarization angle
%

ex=sin(thetae).*sin(thetap);
ey=-sin(thetae).*cos(thetap).*cos(phip)-cos(thetae).*sin(phip);
ez=-sin(thetae).*cos(thetap).*sin(phip)+cos(thetae).*cos(phip);
%k parameters
kwx=-k.*cos(thetap);
kwy=-k.*sin(thetap).*cos(phip);
kwz=-k.*sin(thetap).*sin(phip);
% d factor
D=cosh(gamma.*L).*(Zs+Zl)+sinh(gamma.*L).*(Zc+Zs.*Zl./Zc);
kx=(kwx.*a.*p)./(2.*pi);
ky=(kwy.*a.*p)./(2.*pi);
kz=(kwz.*a.*p)./(2.*pi);
num1=(gamma+j.*kz).^2+a.^2;
num2=(gamma-j.*kz).^2+a.^2;
F1=( a.*exp(gamma.*L)-( a.*cos(a.*L)+ (gamma+j.*kz).*sin(a.*L) ).*exp(-j.*kz.*L) )./num1;
F2=( a.*exp(-gamma.*L)-( a.*cos(a.*L)+ (-gamma+j.*kz).*sin(a.*L) ).*exp(-j.*kz.*L) )./num2;
K1=( (gamma+j.*kz).*exp(gamma.*L)+( a.*sin(a.*L)-(gamma+j.*kz).*cos(a.*L) ).*exp(-j.*kz.*L) )./num1;
K2=( (-gamma+j.*kz).*exp(-gamma.*L)+( a.*sin(a.*L)+(gamma-j.*kz).*cos(a.*L) )*exp(-j.*kz.*L) )./num2;
%F, K
Fplus=F1+F2;
Fminus=F1-F2;
Kplus=K1+K2;
Kminus=K1-K2;
Fplusstar=F1.*exp(-gamma.*L)+F2.*exp(gamma.*L);
Fminusstar=-F1.*exp(-gamma.*L)+F2.*exp(gamma.*L);
Kplusstar=K1.*exp(-gamma.*L)+K2.*exp(gamma.*L);
Kminusstar=-K1.*exp(-gamma.*L)+K2.*exp(gamma.*L);
I10=(a.*ex+j.*ky.*ez).*(Fplus+(Zl./Zc).*Fminus)-(a.*ey-j.*kx.*ez).*(Kplus+(Zl./Zc).*Kminus);
I20=2.*(ex.*cos(a.*L)+ey.*sin(a.*L)).*exp(-j.*kz.*L)-2.*ex.*(cosh(gamma.*L)+(Zl./Zc).*sinh(gamma.*L));
I0=-((E0.*s)./(2.*D)).*(I10+I20);
%current and voltage
I1L=(a.*ex+j.*ky.*ez).*(Fplusstar+(Zs./Zc).*Fminusstar)-(a.*ey-j.*kx.*ez).*(Kplusstar+(Zs./Zc).*Kminusstar);
I2L=2.*(cosh(gamma.*L)+(Zs./Zc).*sinh(gamma.*L) ).*(ex.*cos(a.*L)+ey.*sin(a.*L) ).*exp(-j.*kz.*L)-2.*ex;
%
P=(E0.*s)./(2.*D);
IL=-P.*(I1L+I2L);
zeta=IL./(E0.*P);
xi=E0.*zeta;
xi_prime=zeta./(P.*zeta);
%distributin amplitude
g=amplitude(xi_prime)
f=(1./(P.*zeta)).*g;

function amp=amplitude(z)
%uniform distribution of em field
global E0;
amp=1./E0;