/*  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 <stdio.h>
#include <math.h>
#include <complex.h>

void InvLapTrans(complex<double> (*vs)(complex<double> s),double tstart,
    double tend,double tdelta)
{
    double t,vt;
    int count;
static  const complex<double> z[5] = {
        complex<double>(11.83009373916819,1.593753005885813),
        complex<double>(11.22085377939519,4.792964167565670),
        complex<double>(9.933383722175002,8.033106334266296),
        complex<double>(7.781146264464616,11.36889164904993),
        complex<double>(4.234522494797000,14.95704378128156)};    
static  const complex<double> Kp[5] = {
        complex<double>(16286.62368050479,-139074.7115516051),
        complex<double>(-28178.11171305163,74357.58237274176),
        complex<double>(14629.74025233142,-19181.80818501836),
        complex<double>(-2870.418161032078,1674.109484084304),
        complex<double>(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;
    }
}