/*  tdrsp.cpp -- time domain (step) response by inverse laplace transform.
 *  (C) 2000, C. Bond
 *
 *  Based on "Computer methods for circuit analysis and design,"
 *  by J. Vlach and K.Singhai. The authors do not make it clear
 *  that the 'stepping method' (implemented here) produces the
 *  step response of the circuit, not the impulse response.
 *
 *  Find: h(t) given H(s)
 *                                              
 *            b(0) + b(1)*s + b(2)*s^2 + ... + b(m-1)*s^(m-1)
 *  H(s) = -----------------------------------------------------
 *         a(0) + a(1)*s + a(2)*s^2 + ... + a(n-1)*s^(n-1) + s^n
 *
 *  n >= m,     n <= 15,    n >= 2
 *
 *  To find the impulse response of a filter, multiply the above
 *  Laplace transform by 's'. This is the same as shifting the
 *  numerator polynomial up one order and inserting a zero in the
 *  constant position.
 *
 *
 *  Third order Butterworth filter impulse response:
 *                   s
 *  H(s) = ------------------------
 *          1 + 2*s + 2*s^2 + s^3
 */
#include <complex.h>

static complex<double> qi[15][5];

int TDResp(double *a,double *b,int n,int m,double tend,
    double tdel,int kk,double *rsp)
{
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)};

    double x[16],xt[16],t,respon;
    complex<double> res[5], pk[5], sk, p, pp, resk, xn;
    int i,j,k,l,idx;

// Check input parameters
    if ((n < 2) || (m > n)) return 1;

    a[n] = 1.0;
// Initialize state vector.
    for (i=0;i<n;i++)
        x[i] = xt[i] = 0.0;

// Compute LU factorization at the 5 preselected frequencies.
    for (k=0;k<5;k++) {
        sk = z[k] / tdel;
        res[k] = Kp[k] / tdel;
        pk[k] = 1.0 / sk;
        qi[0][k] = a[0];
        for (j=1;j<n;j++)
            qi[j][k] = a[j] + qi[j-1][k] * pk[k];
        qi[n-1][k] = 1.0 / (qi[n-1][k] + sk);
    }
    t = 0.0;
    idx = 0;
// time domain computation starts here.
    while (t <= tend) {
        for (l=0;l<kk;l++) {
            for (k=0;k<5;k++) {
                p = pk[k];
                resk = res[k];
                pp = complex<double>(-1.0,0.0);
                for (j=0;j<n-1;j++)
                    pp += (qi[j][k] * x[j]);
                xn = (x[n-1] - p*pp) * qi[n-1][k];
                xt[n-1] -= real(resk*xn);
                for (i=0;i<n-1;i++) {
                    xn = (x[n-i-2] + xn) * p;
                    xt[n-i-2] -= real(xn * resk);
                }
            }
            for (i=0;i<n;i++) {
                x[i] = xt[i];
                xt[i] = 0.0;
            }
        }
// x now contains new state vector.
        respon = 0.0;
        for (i=0;i<m;i++)
            respon += (b[i]*x[i]);
// response computed from state vector at ouput points
        rsp[idx++] = respon;
        t += (kk*tdel);
    }
    return 0;
}