import numpy as np
np.set_printoptions(precision=2)

def load_getdp_integral(filename):
    return np.loadtxt(filename)[1]+1j*np.loadtxt(filename)[2]

p_max   = 3
Tmatrix = np.zeros((2*p_max,2*p_max),dtype=complex)

for k1 in range(2*p_max):
    for k2 in range(2*p_max):
        if k1<p_max and k2<p_max: 
            Tmatrix[k1,k2] = load_getdp_integral('resT/fhM_pe%gpo%g.dat'%(k1+1,k2+1))
        if k1<p_max and k2>=p_max:
            Tmatrix[k1,k2] = load_getdp_integral('resT/feM_pe%gpo%g.dat'%(k1+1,k2+1-p_max))
        if k1>=p_max and k2<p_max:
            Tmatrix[k1,k2] = load_getdp_integral('resT/fhN_pe%gpo%g.dat'%(k1+1-p_max,k2+1))
        if k1>=p_max and k2>=p_max:
            Tmatrix[k1,k2] = load_getdp_integral('resT/feN_pe%gpo%g.dat'%(k1+1-p_max,k2+1-p_max))

print('|2*Tmatrix+1|\n',np.abs(2*Tmatrix+1))
print('Re Tmatrix\n',Tmatrix.real)
print('Im Tmatrix\n',Tmatrix.imag)