/* invlapt.cpp -- inverse laplace transform
* (C) 1998, C. Bond. All rights reserved.
* Based on FORTRAN algorithm given in "Computer methods for
* circuit analysis and design", by J. Vlach and K.Singhai.
*
* 'vs' is a function of 's' -- generally a rational function
* of polynomials in 's' representing the transfer function of
* a lumped parameter filter. It returns the
* complex value of the function, given
* the complex value of 's' as a parameter. 'vs' is not restricted
* to ratios of polynomials. Distributed networks, such as strip
* lines and transmission lines can also be solved. It is only
* necessary that the terminal behavior be expressible as a
* rational function in 's'.
*
* 'tstart', 'tend' and 'tdelta' are descriptors for the time
* interval during which the response function is to be
* determined. Do not start at t = 0, as a divide by 't' is
* present. Use a small value such as 10^-8 instead.
*
* 'vt' is the (real) value of the response at time 't'.
*/
#include
#include
#include
void InvLapTrans(complex (*vs)(complex s),double tstart,
double tend,double tdelta)
{
double t,vt;
int count;
static const complex z[5] = {
complex(11.83009373916819,1.593753005885813),
complex(11.22085377939519,4.792964167565670),
complex(9.933383722175002,8.033106334266296),
complex(7.781146264464616,11.36889164904993),
complex(4.234522494797000,14.95704378128156)};
static const complex Kp[5] = {
complex(16286.62368050479,-139074.7115516051),
complex(-28178.11171305163,74357.58237274176),
complex(14629.74025233142,-19181.80818501836),
complex(-2870.418161032078,1674.109484084304),
complex(132.1659412474876,17.47674798877164)};
t = tstart;
count = 0;
while (t <= tend) {
vt = 0.0;
for (int i=0;i<5;i++) {
vt = vt - real(vs(z[i] / t) * Kp[i]);
}
vt /= t;
printf("%d %lf\n",count,vt);
count++;
t += tdelta;
}
}