response_convolution() [2/4]

template<typename T >

void cntr::response_convolution	(	int	tstp,
		CPLX &	cc,
		GREEN_TSTP &	W,
		int	a1,
		int	a2,
		function< T > &	f,
		int	b1,
		int	b2,
		int	SolveOrder,
		T	beta,
		T	h
	)

Definition at line 68 of file cntr_response_convolution_impl.hpp.

                                                                       {
    int ntau = W.ntau();
    T dtau = beta / ntau;
    int sizew = W.size1(), idxw = a1 * sizew + a2, sw = sizew * sizew;
    int sizef = f.size1(), idxf = b1 * sizef + b2, sf = sizef * sizef;
    CPLX res;

    assert(W.tstp() >= tstp);
    assert(tstp == -1 || tstp >= SolveOrder);
    assert(SolveOrder <= ntau);
    assert(b1 >= 0 && b1 <= sizef-1);
    assert(b2 >= 0 && b2 <= sizef-1);
    assert(a1 >= 0 && a1 <= sizew-1);
    assert(a2 >= 0 && a2 <= sizew-1);

    if (tstp == -1) {
        response_integrate<T>(ntau, dtau, res, W.matptr(0), sw, idxw, integration::I<T>(SolveOrder));
        cc = res * f.ptr(-1)[idxf];
    } else {
        response_integrate<T>(ntau, dtau, res, W.tvptr(0), sw, idxw, integration::I<T>(SolveOrder));
        cc = res * (CPLX(0, -1.0) * f.ptr(-1)[idxf]);
        response_integrate<T>(tstp, h, res, W.retptr(0), sw, idxw, f.ptr(0), sf, idxf,
                              integration::I<T>(SolveOrder));
        cc += res;
    }
}