| [f,a,omega,integr,epsabs,result,abserr,neval,ier,limlst,lst,leniw,maxp1,lenw,iwork,work]=dqawf(f,a,omega,integr,epsabs,result,abserr,neval,ier,limlst,lst,leniw,maxp1,lenw,iwork,work); |
function [f,a,omega,integr,epsabs,result,abserr,neval,ier,limlst,lst,leniw,maxp1,lenw,iwork,work]=dqawf(f,a,omega,integr,epsabs,result,abserr,neval,ier,limlst,lst,leniw,maxp1,lenw,iwork,work);
%***BEGIN PROLOGUE DQAWF
%***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 (QAWF-S, DQAWF-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 integrand WEIGHT function
%
% INTEGR - Integer
% Indicates which of the WEIGHT functions 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.
%
% ON RETURN
% RESULT - doubleprecision
% Approximation to the integral
%
% 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 - Integer
% 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.1 or
% LENIW.LT.(LIMLST+2) or MAXP1.LT.1 or
% LENW.LT.(LENIW*2+MAXP1*25).
% RESULT, 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 first LST elements of
% vector IWORK. Here LST is the number of
% cycles actually needed (see below).
% IWORK(K) = 1 The maximum number of
% subdivisions (=(LENIW-LIMLST)
% /2) 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
% on this cycle.
% = 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
% IWORK(K).
% If OMEGA = 0 and INTEGR = 1,
% The integral is calculated by means of DQAGIE,
% and IER = IWORK(1) (with meaning as described
% for IWORK(K),K = 1).
%
% DIMENSIONING PARAMETERS
% LIMLST - Integer
% LIMLST gives an upper bound on the number of
% cycles, LIMLST.GE.3.
% If LIMLST.LT.3, the routine will end with IER = 6.
%
% LST - Integer
% On return, LST indicates the number of cycles
% actually needed for the integration.
% If OMEGA = 0, then LST is set to 1.
%
% LENIW - Integer
% Dimensioning parameter for IWORK. On entry,
% (LENIW-LIMLST)/2 equals the maximum number of
% subintervals allowed in the partition of each
% cycle, LENIW.GE.(LIMLST+2).
% If LENIW.LT.(LIMLST+2), the routine will end with
% IER = 6.
%
% MAXP1 - Integer
% MAXP1 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.
% If MAXP1.LT.1, the routine will end with IER = 6.
% LENW - Integer
% Dimensioning parameter for WORK
% LENW must be at least LENIW*2+MAXP1*25.
% If LENW.LT.(LENIW*2+MAXP1*25), the routine will
% end with IER = 6.
%
% WORK ARRAYS
% IWORK - Integer
% Vector of dimension at least LENIW
% On return, IWORK(K) FOR K = 1, 2, ..., LST
% contain the error flags on the cycles.
%
% WORK - doubleprecision
% Vector of dimension at least
% On return,
% WORK(1), ..., WORK(LST) contain the integral
% approximations over the cycles,
% WORK(LIMLST+1), ..., WORK(LIMLST+LST) contain
% the error estimates over the cycles.
% further elements of WORK have no specific
% meaning for the user.
%
%***REFERENCES (NONE)
%***ROUTINES CALLED DQAWFE, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800101 DATE WRITTEN
% 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)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
%***end PROLOGUE DQAWF
%
persistent l1 l2 l3 l4 l5 l6 limit ll2 lvl ;
if isempty(limit), limit=0; end;
if isempty(ll2), ll2=0; end;
if isempty(lvl), lvl=0; end;
if isempty(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(l3), l3=0; end;
if isempty(l4), l4=0; end;
if isempty(l5), l5=0; end;
if isempty(l6), l6=0; end;
%
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%
%
% CHECK VALIDITY OF LIMLST, LENIW, MAXP1 AND LENW.
%
%***FIRST EXECUTABLE STATEMENT DQAWF
ier = 6;
neval = 0;
result = 0.0d+00;
abserr = 0.0d+00;
if( limlst>=3 && leniw>=(limlst+2) && maxp1>=1 &&lenw>=(leniw.*2+maxp1.*25) )
%
% PREPARE CALL FOR DQAWFE
%
limit =fix(fix((leniw-limlst)./2));
l1 = fix(limlst + 1);
l2 = fix(limlst + l1);
l3 = fix(limit + l2);
l4 = fix(limit + l3);
l5 = fix(limit + l4);
l6 = fix(limit + l5);
ll2 = fix(limit + l1);
[f,a,omega,integr,epsabs,limlst,limit,maxp1,result,abserr,neval,ier,dumvar13,dumvar14,idumvar13,lst,dumvar17,dumvar18,dumvar19,dumvar20,idumvar14,dumvar22,dumvar23]=dqawfe(f,a,omega,integr,epsabs,limlst,limit,maxp1,result,abserr,neval,ier,work(sub2ind(size(work),max(1,1)):end),work(sub2ind(size(work),max(l1,1)):end),iwork(sub2ind(size(work),max(1,1)):end),lst,work(sub2ind(size(work),max(l2,1)):end),work(sub2ind(size(work),max(l3,1)):end),work(sub2ind(size(work),max(l4,1)):end),work(sub2ind(size(work),max(l5,1)):end),iwork(sub2ind(size(work),max(l1,1)):end),iwork(sub2ind(size(iwork),max(ll2,1)):end),work(l6:end)); dumvar13i=find((work(sub2ind(size(work),max(1,1)):end))~=(dumvar13));dumvar14i=find((work(sub2ind(size(work),max(l1,1)):end))~=(dumvar14));dumvar17i=find((work(sub2ind(size(work),max(l2,1)):end))~=(dumvar17));dumvar18i=find((work(sub2ind(size(work),max(l3,1)):end))~=(dumvar18));dumvar19i=find((work(sub2ind(size(work),max(l4,1)):end))~=(dumvar19));dumvar20i=find((work(sub2ind(size(work),max(l5,1)):end))~=(dumvar20));dumvar22i=find((iwork(sub2ind(size(iwork),max(ll2,1)):end))~=(dumvar22));dumvar23i=find((work(l6:end))~=(dumvar23)); work(1-1+dumvar13i)=dumvar13(dumvar13i); work(l1-1+dumvar14i)=dumvar14(dumvar14i); work(l2-1+dumvar17i)=dumvar17(dumvar17i); work(l3-1+dumvar18i)=dumvar18(dumvar18i); work(l4-1+dumvar19i)=dumvar19(dumvar19i); work(l5-1+dumvar20i)=dumvar20(dumvar20i); iwork(ll2-1+dumvar22i)=dumvar22(dumvar22i); work(l6-1+dumvar23i)=dumvar23(dumvar23i);
%
% CALL ERROR HANDLER IF NECESSARY
%
lvl = 0;
end;
if( ier==6 )
lvl = 1;
end;
if( ier~=0 )
[dumvar1,dumvar2,dumvar3,ier,lvl]=xermsg('SLATEC','DQAWF','ABNORMAL RETURN',ier,lvl);
end;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK DQAWOE
|
|
