function [xw,w]=q_ossinf(M,omega,q)
% Quadrature of oscillating function over infinite interval
%   xw=q_ossinf(M,omega,q)
%   [x,w]=q_ossinf(M,omega,q)
% x = xw(:,1) = nodes, w = xw(:,2) = weights 
% sum(w.*f(x))  approximates  integral_0^inf f(t) sin(omega*t+q*pi) dt.
% In particular, use q=0.5 for  integral_0^inf f(t) cos(omega*t) dt.
% Method: double exponential transformation with original step size pi/M. 
% The number of points in the formula is approximately 5*M.
% Source: Ooura and Mori, JCAM 112 (1999) 229-241
% WARNING: this routine may deliver some inaccurate weights due to numerical
%   instability if M is large or q is close but not equal to an integer.  
%   Read the comments in the source for details. 


magic=8.835; h=pi/M; n=floor(magic/h); n0=floor(magic/h-q); 
t=h*(-n0-q:n); dt=ones(1,length(t));
% Perturb a possible zero entry to avoid a NaN.  The correct value
%   will be inserted later.  Numerical instability arises when t is small.
nz=find(t==0); t(nz)=1e-8; 

beta=0.25; alpha=beta/sqrt(1+M*log(1+M)/(4*pi));

% Formulas mathematically equivalent to the original
% gamma=(alpha+beta)/2; delta=(alpha-beta)/2;
% phi(t) = t*exp(v)/(2*sinh(v)) where v=v(t),
%   v(t) = t+gamma*sinh(t)+2*delta*(sinh(t))^2, 
%  v'(t) = 1+gamma*cosh(t)+delta*sinh(t)

% Formulated in terms of sinh for numerical stability
het=exp(t/2); et=het.*het; dx=dt.*(-2-alpha./et-beta*et); 
x=-2*(t+(alpha./et+beta).*het.*sinh(t/2));
% Again
hex=exp(x/2); ex=hex.*hex; dy=-dx.*ex; y=2*hex.*sinh(-x/2);
phi=t./y; dphi=dt./y-phi.*(dy./y);  
% This is the crunch.  When t(k) is small, this formula for dphi(k)
%   suffers from cancellation.

% Correct the perturbed entry
phi(nz)=1/(2+alpha+beta); dphi(nz)=((alpha-beta)*phi(nz).^2+1)/2;

A=M*phi; x=A/omega; w=dphi.*sin(A+q*pi)*(M/omega); 
% Anti-roundoff refinement: (4.7) in paper
m=find(abs(ex)<0.1); w0=w(m);
w(m)=(-1).^m.*dphi(m).*sin(A(m).*ex(m))*(M/omega); 
% Possible sign change
if sum(w(m(1:2)).*w0(1:2))<0, w(m)=-w(m); end
w=h*w;

% remove NaN's
k=isnan(w); w(k)=[]; x(k)=[];
% Remove underflows
k=find(w==0 | x==0); w(k)=[]; x(k)=[];  
 
xw=[x;w]'; if nargout==2, xw=x; end