| [f,a,omega,integr,epsabs,limlst,limit,maxp1,resultf,abserr,neval,ier,rslst,erlst,ierlst,lst,alist,blist,rlist,elist,iord,nnlog,chebmo]=dqawfe(f,a,omega,integr,epsabs,limlst,limit,maxp1,resultf,abserr,neval,ier,rslst,erlst,ierlst,lst,alist,blist,rlist,elis |
function [f,a,omega,integr,epsabs,limlst,limit,maxp1,resultf,abserr,neval,ier,rslst,erlst,ierlst,lst,alist,blist,rlist,elist,iord,nnlog,chebmo]=dqawfe(f,a,omega,integr,epsabs,limlst,limit,maxp1,resultf,abserr,neval,ier,rslst,erlst,ierlst,lst,alist,blist,rlist,elist,iord,nnlog,chebmo);
%***BEGIN PROLOGUE DQAWFE
%***PURPOSE The routine calculates an approximation result to a
% given Fourier integral
% I = Integral of F(X)*W(X) over (A,INFINITY)
% where W(X)=COS(OMEGA*X) or W(X)=SIN(OMEGA*X),
% hopefully satisfying following claim for accuracy
% ABS(I-RESULT).LE.EPSABS.
%***LIBRARY SLATEC (QUADPACK)
%***CATEGORY H2A3A1
%***TYPE doubleprecision (QAWFE-S, DQAWFE-D)
%***KEYWORDS AUTOMATIC INTEGRATOR, CONVERGENCE ACCELERATION,
% FOURIER INTEGRALS, INTEGRATION BETWEEN ZEROS, QUADPACK,
% QUADRATURE, SPECIAL-PURPOSE INTEGRAL
%***AUTHOR Piessens, Robert
% Applied Mathematics and Programming Division
% K. U. Leuven
% de Doncker, Elise
% Applied Mathematics and Programming Division
% K. U. Leuven
%***DESCRIPTION
%
% Computation of Fourier integrals
% Standard fortran subroutine
% doubleprecision version
%
% PARAMETERS
% ON ENTRY
% F - doubleprecision
% function subprogram defining the integrand
% function F(X). The actual name for F needs to
% be declared E X T E R N A L in the driver program.
%
% A - doubleprecision
% Lower limit of integration
%
% OMEGA - doubleprecision
% Parameter in the WEIGHT function
%
% INTEGR - Integer
% Indicates which WEIGHT function is used
% INTEGR = 1 W(X) = COS(OMEGA*X)
% INTEGR = 2 W(X) = SIN(OMEGA*X)
% If INTEGR.NE.1.AND.INTEGR.NE.2, the routine will
% end with IER = 6.
%
% EPSABS - doubleprecision
% absolute accuracy requested, EPSABS.GT.0
% If EPSABS.LE.0, the routine will end with IER = 6.
%
% LIMLST - Integer
% LIMLST gives an upper bound on the number of
% cycles, LIMLST.GE.1.
% If LIMLST.LT.3, the routine will end with IER = 6.
%
% LIMIT - Integer
% Gives an upper bound on the number of subintervals
% allowed in the partition of each cycle, LIMIT.GE.1
% each cycle, LIMIT.GE.1.
%
% MAXP1 - Integer
% Gives an upper bound on the number of
% Chebyshev moments which can be stored, I.E.
% for the intervals of lengths ABS(B-A)*2**(-L),
% L=0,1, ..., MAXP1-2, MAXP1.GE.1
%
% ON RETURN
% RESULTF - doubleprecision
% Approximation to the integral X
%
% ABSERR - doubleprecision
% Estimate of the modulus of the absolute error,
% which should equal or exceed ABS(I-RESULT)
%
% NEVAL - Integer
% Number of integrand evaluations
%
% IER - IER = 0 Normal and reliable termination of
% the routine. It is assumed that the
% requested accuracy has been achieved.
% IER.GT.0 Abnormal termination of the routine. The
% estimates for integral and error are less
% reliable. It is assumed that the requested
% accuracy has not been achieved.
% ERROR MESSAGES
% If OMEGA.NE.0
% IER = 1 Maximum number of cycles allowed
% Has been achieved., i.e. of subintervals
% (A+(K-1)C,A+KC) where
% C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
% for K = 1, 2, ..., LST.
% One can allow more cycles by increasing
% the value of LIMLST (and taking the
% according dimension adjustments into
% account).
% Examine the array IWORK which contains
% the error flags on the cycles, in order to
% look for eventual local integration
% difficulties. If the position of a local
% difficulty can be determined (e.g.
% SINGULARITY, DISCONTINUITY within the
% interval) one will probably gain from
% splitting up the interval at this point
% and calling appropriate integrators on
% the subranges.
% = 4 The extrapolation table constructed for
% convergence acceleration of the series
% formed by the integral contributions over
% the cycles, does not converge to within
% the requested accuracy. As in the case of
% IER = 1, it is advised to examine the
% array IWORK which contains the error
% flags on the cycles.
% = 6 The input is invalid because
% (INTEGR.NE.1 AND INTEGR.NE.2) or
% EPSABS.LE.0 or LIMLST.LT.3.
% RESULTF, ABSERR, NEVAL, LST are set
% to zero.
% = 7 Bad integrand behaviour occurs within one
% or more of the cycles. Location and type
% of the difficulty involved can be
% determined from the vector IERLST. Here
% LST is the number of cycles actually
% needed (see below).
% IERLST(K) = 1 The maximum number of
% subdivisions (= LIMIT) has
% been achieved on the K th
% cycle.
% = 2 Occurrence of roundoff error
% is detected and prevents the
% tolerance imposed on the
% K th cycle, from being
% achieved.
% = 3 Extremely bad integrand
% behaviour occurs at some
% points of the K th cycle.
% = 4 The integration procedure
% over the K th cycle does
% not converge (to within the
% required accuracy) due to
% roundoff in the
% extrapolation procedure
% invoked on this cycle. It
% is assumed that the result
% on this interval is the
% best which can be obtained.
% = 5 The integral over the K th
% cycle is probably divergent
% or slowly convergent. It
% must be noted that
% divergence can occur with
% any other value of
% IERLST(K).
% If OMEGA = 0 and INTEGR = 1,
% The integral is calculated by means of DQAGIE
% and IER = IERLST(1) (with meaning as described
% for IERLST(K), K = 1).
%
% RSLST - doubleprecision
% Vector of dimension at least LIMLST
% RSLST(K) contains the integral contribution
% over the interval (A+(K-1)C,A+KC) where
% C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
% K = 1, 2, ..., LST.
% Note that, if OMEGA = 0, RSLST(1) contains
% the value of the integral over (A,INFINITY).
%
% ERLST - doubleprecision
% Vector of dimension at least LIMLST
% ERLST(K) contains the error estimate corresponding
% with RSLST(K).
%
% IERLST - Integer
% Vector of dimension at least LIMLST
% IERLST(K) contains the error flag corresponding
% with RSLST(K). For the meaning of the local error
% flags see description of output parameter IER.
%
% LST - Integer
% Number of subintervals needed for the integration
% If OMEGA = 0 then LST is set to 1.
%
% ALIST, BLIST, RLIST, ELIST - doubleprecision
% vector of dimension at least LIMIT,
%
% IORD, NNLOG - Integer
% Vector of dimension at least LIMIT, providing
% space for the quantities needed in the subdivision
% process of each cycle
%
% CHEBMO - doubleprecision
% Array of dimension at least (MAXP1,25), providing
% space for the Chebyshev moments needed within the
% cycles
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, DQAGIE, DQAWOE, DQELG
%***REVISION HISTORY (YYMMDD)
% 800101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 891009 Removed unreferenced variable. (WRB)
% 891009 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
%***end PROLOGUE DQAWFE
%
persistent abseps c1 c2 correc cyclef dl drl ep eps epsa errsum fact firstCall gt ktmin l last ll momcom nev nres numrl2 p p1 pi psum res3la reseps uflow ; if isempty(firstCall),firstCall=1;end;
if isempty(abseps), abseps=0; end;
if isempty(correc), correc=0; end;
if isempty(cyclef), cyclef=0; end;
if isempty(c1), c1=0; end;
if isempty(c2), c2=0; end;
if isempty(dl), dl=0; end;
if isempty(drl), drl=0; end;
if isempty(ep), ep=0; end;
if isempty(eps), eps=0; end;
if isempty(epsa), epsa=0; end;
if isempty(errsum), errsum=0; end;
if isempty(fact), fact=0; end;
if isempty(p), p=0; end;
if isempty(pi), pi=0; end;
if isempty(p1), p1=0; end;
if isempty(psum), psum=zeros(1,52); end;
if isempty(reseps), reseps=0; end;
if isempty(res3la), res3la=zeros(1,3); end;
if isempty(uflow), uflow=0; end;
if isempty(ktmin), ktmin=0; end;
if isempty(l), l=0; end;
if isempty(last), last=0; end;
if isempty(ll), ll=0; end;
if isempty(momcom), momcom=0; end;
if isempty(nev), nev=0; end;
if isempty(nres), nres=0; end;
if isempty(numrl2), numrl2=0; end;
if isempty(gt), gt=0; end;
%
alist_shape=size(alist);alist=reshape(alist,1,[]);
blist_shape=size(blist);blist=reshape(blist,1,[]);
chebmo_orig=chebmo;chebmo_shape=[maxp1,25];chebmo=reshape([chebmo_orig(1:min(prod(chebmo_shape),numel(chebmo_orig))),zeros(1,max(0,prod(chebmo_shape)-numel(chebmo_orig)))],chebmo_shape);
elist_shape=size(elist);elist=reshape(elist,1,[]);
erlst_shape=size(erlst);erlst=reshape(erlst,1,[]);
ierlst_shape=size(ierlst);ierlst=reshape(ierlst,1,[]);
iord_shape=size(iord);iord=reshape(iord,1,[]);
nnlog_shape=size(nnlog);nnlog=reshape(nnlog,1,[]);
rlist_shape=size(rlist);rlist=reshape(rlist,1,[]);
rslst_shape=size(rslst);rslst=reshape(rslst,1,[]);
%
%
%
% THE DIMENSION OF PSUM IS DETERMINED BY THE VALUE OF
% LIMEXP IN SUBROUTINE DQELG (PSUM MUST BE OF DIMENSION
% (LIMEXP+2) AT LEAST).
%
% LIST OF MAJOR VARIABLES
% -----------------------
%
% C1, C2 - END POINTS OF SUBINTERVAL (OF LENGTH CYCLEF)
% CYCLEF - (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA)
% PSUM - VECTOR OF DIMENSION AT LEAST (LIMEXP+2)
% (SEE ROUTINE DQELG)
% PSUM CONTAINS THE PART OF THE EPSILON TABLE
% WHICH IS STILL NEEDED FOR FURTHER COMPUTATIONS.
% EACH ELEMENT OF PSUM IS A PARTIAL SUM OF THE
% SERIES WHICH SHOULD SUM TO THE VALUE OF THE
% INTEGRAL.
% ERRSUM - SUM OF ERROR ESTIMATES OVER THE SUBINTERVALS,
% CALCULATED CUMULATIVELY
% EPSA - ABSOLUTE TOLERANCE REQUESTED OVER CURRENT
% SUBINTERVAL
% CHEBMO - ARRAY CONTAINING THE MODIFIED CHEBYSHEV
% MOMENTS (SEE ALSO ROUTINE DQC25F)
%
if firstCall, p=[0.9d+00]; end;
if firstCall, pi=[3.14159265358979323846264338327950d0]; end;
firstCall=0;
%
% TEST ON VALIDITY OF PARAMETERS
% ------------------------------
%
%***FIRST EXECUTABLE STATEMENT DQAWFE
resultf = 0.0d+00;
abserr = 0.0d+00;
neval = 0;
lst = 0;
ier = 0;
if((integr~=1 && integr~=2) || epsabs<=0.0d+00 ||limlst<3 )
ier = 6;
end;
gt=0;
if( ier~=6 )
while (1);
if( omega~=0.0d+00 )
%
% INITIALIZATIONS
% ---------------
%
l = fix(abs(omega));
dl = 2.*l + 1;
cyclef = dl.*pi./abs(omega);
ier = 0;
ktmin = 0;
neval = 0;
numrl2 = 0;
nres = 0;
c1 = a;
c2 = cyclef + a;
p1 = 0.1d+01 - p;
[uflow ]=d1mach(1);
eps = epsabs;
if( epsabs>uflow./p1 )
eps = epsabs.*p1;
end;
ep = eps;
fact = 0.1d+01;
correc = 0.0d+00;
abserr = 0.0d+00;
errsum = 0.0d+00;
%
% MAIN DO-LOOP
% ------------
%
for lst = 1 : limlst;
%
% INTEGRATE OVER CURRENT SUBINTERVAL.
%
epsa = eps.*fact;
[f,c1,c2,omega,integr,epsa,dumvar7,limit,lst,maxp1,rslst(lst),erlst(lst),nev,ierlst(lst),last,alist,blist,rlist,elist,iord,nnlog,momcom,chebmo]=dqawoe(f,c1,c2,omega,integr,epsa,0.0d+00,limit,lst,maxp1,rslst(lst),erlst(lst),nev,ierlst(lst),last,alist,blist,rlist,elist,iord,nnlog,momcom,chebmo);
neval = fix(neval + nev);
fact = fact.*p;
errsum = errsum + erlst(lst);
drl = 0.5d+02.*abs(rslst(lst));
%
% TEST ON ACCURACY WITH PARTIAL SUM
%
if((errsum+drl)<=epsabs && lst>=6 )
gt=1;
break;
end;
correc = max(correc,erlst(lst));
if( ierlst(lst)~=0 )
eps = max(ep,correc.*p1);
end;
if( ierlst(lst)~=0 )
ier = 7;
end;
if( ier==7 &&(errsum+drl)<=correc.*0.1d+02 && lst>5 )
gt=1;
break;
end;
numrl2 = fix(numrl2 + 1);
if( lst>1 )
psum(numrl2) = psum(ll) + rslst(lst);
if( lst~=2 )
%
% TEST ON MAXIMUM NUMBER OF SUBINTERVALS
%
if( lst==limlst )
ier = 1;
end;
%
% PERFORM NEW EXTRAPOLATION
%
[numrl2,psum,reseps,abseps,res3la,nres]=dqelg(numrl2,psum,reseps,abseps,res3la,nres);
%
% TEST WHETHER EXTRAPOLATED RESULT IS INFLUENCED BY ROUNDOFF
%
ktmin = fix(ktmin + 1);
if( ktmin>=15 && abserr<=0.1d-02.*(errsum+drl) )
ier = 4;
end;
if( abseps<=abserr || lst==3 )
abserr = abseps;
resultf = reseps;
ktmin = 0;
%
% IF IER IS NOT 0, CHECK WHETHER DIRECT RESULT (PARTIAL SUM)
% OR EXTRAPOLATED RESULT YIELDS THE BEST INTEGRAL
% APPROXIMATION
%
if((abserr+0.1d+02.*correc)<=epsabs ||(abserr<=epsabs && 0.1d+02.*correc>=epsabs) )
break;
end;
end;
if( ier~=0 && ier~=7 )
break;
end;
end;
else;
psum(1) = rslst(1);
end;
ll = fix(numrl2);
c1 = c2;
c2 = c2 + cyclef;
end;
if(gt==1)
break;
end;
%
% SET FINAL RESULT AND ERROR ESTIMATE
% -----------------------------------
%
abserr = abserr + 0.1d+02.*correc;
if( ier==0 )
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
chebmo_orig(1:prod(chebmo_shape))=chebmo;chebmo=chebmo_orig;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
erlst_shape=zeros(erlst_shape);erlst_shape(:)=erlst(1:numel(erlst_shape));erlst=erlst_shape;
ierlst_shape=zeros(ierlst_shape);ierlst_shape(:)=ierlst(1:numel(ierlst_shape));ierlst=ierlst_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
nnlog_shape=zeros(nnlog_shape);nnlog_shape(:)=nnlog(1:numel(nnlog_shape));nnlog=nnlog_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
rslst_shape=zeros(rslst_shape);rslst_shape(:)=rslst(1:numel(rslst_shape));rslst=rslst_shape;
return;
end;
if( resultf==0.0d+00 || psum(numrl2)==0.0d+00 )
if( abserr>errsum )
break;
end;
if( psum(numrl2)==0.0d+00 )
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
chebmo_orig(1:prod(chebmo_shape))=chebmo;chebmo=chebmo_orig;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
erlst_shape=zeros(erlst_shape);erlst_shape(:)=erlst(1:numel(erlst_shape));erlst=erlst_shape;
ierlst_shape=zeros(ierlst_shape);ierlst_shape(:)=ierlst(1:numel(ierlst_shape));ierlst=ierlst_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
nnlog_shape=zeros(nnlog_shape);nnlog_shape(:)=nnlog(1:numel(nnlog_shape));nnlog=nnlog_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
rslst_shape=zeros(rslst_shape);rslst_shape(:)=rslst(1:numel(rslst_shape));rslst=rslst_shape;
return;
end;
end;
if( abserr./abs(resultf)<=(errsum+drl)./abs(psum(numrl2)) )
if( ier>=1 && ier~=7 )
abserr = abserr + drl;
end;
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
chebmo_orig(1:prod(chebmo_shape))=chebmo;chebmo=chebmo_orig;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
erlst_shape=zeros(erlst_shape);erlst_shape(:)=erlst(1:numel(erlst_shape));erlst=erlst_shape;
ierlst_shape=zeros(ierlst_shape);ierlst_shape(:)=ierlst(1:numel(ierlst_shape));ierlst=ierlst_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
nnlog_shape=zeros(nnlog_shape);nnlog_shape(:)=nnlog(1:numel(nnlog_shape));nnlog=nnlog_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
rslst_shape=zeros(rslst_shape);rslst_shape(:)=rslst(1:numel(rslst_shape));rslst=rslst_shape;
return;
end;
else;
%
% INTEGRATION BY DQAGIE IF OMEGA IS ZERO
% --------------------------------------
%
if( integr==1 )
[f,a,dumvar3,epsabs,dumvar5,limit,resultf,abserr,neval,ier,alist,blist,rlist,elist,iord,last]=dqagie(f,a,1,epsabs,0.0d+00,limit,resultf,abserr,neval,ier,alist,blist,rlist,elist,iord,last);
end;
rslst(1) = resultf;
erlst(1) = abserr;
ierlst(1) = fix(ier);
lst = 1;
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
chebmo_orig(1:prod(chebmo_shape))=chebmo;chebmo=chebmo_orig;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
erlst_shape=zeros(erlst_shape);erlst_shape(:)=erlst(1:numel(erlst_shape));erlst=erlst_shape;
ierlst_shape=zeros(ierlst_shape);ierlst_shape(:)=ierlst(1:numel(ierlst_shape));ierlst=ierlst_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
nnlog_shape=zeros(nnlog_shape);nnlog_shape(:)=nnlog(1:numel(nnlog_shape));nnlog=nnlog_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
rslst_shape=zeros(rslst_shape);rslst_shape(:)=rslst(1:numel(rslst_shape));rslst=rslst_shape;
return;
end;
break;
end;
resultf = psum(numrl2);
abserr = errsum + drl;
end;
alist_shape=zeros(alist_shape);alist_shape(:)=alist(1:numel(alist_shape));alist=alist_shape;
blist_shape=zeros(blist_shape);blist_shape(:)=blist(1:numel(blist_shape));blist=blist_shape;
chebmo_orig(1:prod(chebmo_shape))=chebmo;chebmo=chebmo_orig;
elist_shape=zeros(elist_shape);elist_shape(:)=elist(1:numel(elist_shape));elist=elist_shape;
erlst_shape=zeros(erlst_shape);erlst_shape(:)=erlst(1:numel(erlst_shape));erlst=erlst_shape;
ierlst_shape=zeros(ierlst_shape);ierlst_shape(:)=ierlst(1:numel(ierlst_shape));ierlst=ierlst_shape;
iord_shape=zeros(iord_shape);iord_shape(:)=iord(1:numel(iord_shape));iord=iord_shape;
nnlog_shape=zeros(nnlog_shape);nnlog_shape(:)=nnlog(1:numel(nnlog_shape));nnlog=nnlog_shape;
rlist_shape=zeros(rlist_shape);rlist_shape(:)=rlist(1:numel(rlist_shape));rlist=rlist_shape;
rslst_shape=zeros(rslst_shape);rslst_shape(:)=rslst(1:numel(rslst_shape));rslst=rslst_shape;
end %subroutine dqawfe
%DECK DQAWF
|
|
