#!/usr/bin/octave

% This program, called EULER.m, has been published by Abate and Whitt
% in "The Fourier series method for inverting transforms of
% probability distributions", Queueing Systems, 10, pp. 5-88, and in
% "Numerical Inversion of Laplace Transforms of Probability
% Distributions", ORSA Journal on Computing 7 (1995), pp 36-43. A&W
% wrote it originally in Ubasic. This is a port to Octave, a GNU version
% of Matlab, written by Nicky van Foreest, foreest@math.utwente.nl
% Date: 08/21/01

% I tried to stay as close as possible to the published code to
% facilitate the reading of this program. Doubtless there will be
% numerous improvements to make by exploiting specific functionality
% of Octave. I did not opt to include these in the present
% program (save the definition of the vector C---I was too lazy to type
% in the numbers; and the definition of the variable M). Easy of 
% comparison with the published version had my  highest priority.

M=11;
SU=zeros(1,M+2);
C=zeros(1,M+1);
C(1)=1;
for k=[2:M+1]
  C(k)=(M-(k-2))*C(k-1)/(k-1);
endfor

T=input("TIME=");
A=19.1;
Ntr=15;
U=exp(A/2.)/T;
X=A/(2.*T);
H=pi/T;

% The function fnRf has to be defined before we can invoke it
function Rfs=fnRf(X,Y)
  S=X+i*Y;
  Rho=0.75;
  Mean=1.;
  Gs=1./sqrt(1.+2.*S);
  Gse=(1.-Gs)/(Mean*S);
  Fs=(1.-Gse)/(S*(1.-Rho*Gse));
  Rfs=real(Fs);
endfunction 

Sum=fnRf(X,0)/2.;
for N=[1:Ntr]
  Y=N*H;
  Sum=Sum+(-1)^N*fnRf(X,Y);
endfor

SU(1)=Sum;
for K=[1:M+1]
  N=Ntr+K;
  Y=N*H;
  SU(K+1)=SU(K)+(-1)^N*fnRf(X,Y);
endfor  

Avgsu=0.;
Avgsu1=0.;
for J=[1:M+1]
  Avgsu=Avgsu+C(J)*SU(J);
  Avgsu1=Avgsu1+C(J)*SU(J+1);
endfor
Fun=U*Avgsu/2048;
Fun1=U*Avgsu1/2048;

Errt=abs(Fun-Fun1)/2.;

printf("\nTime=%f\tFUNCTION=%10.9f\n",T,Fun1);
printf("Truncation Error Estimate=%e\n",Errt);