Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: quadl_error
Date: Mon, 10 Dec 2012 12:14:06 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 102
Message-ID: <ka4jme$g91$1@newscl01ah.mathworks.com>
References: <ka4d4g$pj0$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-03-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1355141646 16673 172.30.248.48 (10 Dec 2012 12:14:06 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Mon, 10 Dec 2012 12:14:06 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3799640
Xref: news.mathworks.com comp.soft-sys.matlab:784530

"george veropoulos" <veropgr@yahoo.gr> wrote in message <ka4d4g$pj0$1@newscl01ah.mathworks.com>...
> 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;

Your function does not seem to depend on the independent variable z.
Thus an integration with respect to z should be easy ...

Best wishes
Torsten.