module pub
    implicit none
    integer N,ne,kn
    parameter(N = 4,ne = 100,kn = 200)
    real,parameter::pi = 3.14159265358
    complex,parameter::im = (0.0,1.0)  
    complex ham(N,N)
    complex gf0(N,N)     ! bare green's function
    complex tmat(N,N)   ! T-matrix
    complex::V(N,N) = 0.0
    real t1,t2,t3,t4,del0,mu,eta,v0
    !-----Package needed----------
    integer info
    end module pub
!=================================================================
    program sol
    use pub 
    t1 = -1
    t2 = 1.3
    t3 = -0.85
    t4 = -0.85
    del0 = 0.1
    mu = 1.54
    v0 = 0.4
    eta = 0.05  ! quasiparticle dumping
    !-------------
    !   impurity matrix
    V(1,1) = V0
    V(2,2) = -V0
    V(3,3) = V0
    V(4,4) = -V0
    !-------------
    ! call main(0.07)
    call spectrum(0.07)
    stop
    end program sol
!==================================================================
    subroutine spectrum(omg)
    !   纯净系统的谱函数计算 
    use pub
    integer m1,m2
    real kx,ky,omg
    open(11,file="spec.dat")
    do m1 = -kn,kn
        do m2 = -kn,kn
            kx = m1*pi/kn
            ky = m2*pi/kn
            call gf_clean(kx,ky,omg)
            write(11,*)kx/pi,ky/pi,-aimag((gf0(1,1) + gf0(3,3))/pi)
        end do
    end do
    close(11)
    return
    end subroutine spectrum
!==================================================================
    subroutine main(omg)
    use pub
    integer m1,m2,m3,m4
    real qx,qy,kx,ky,omg
    complex gf1(N,N),gf2(N,N),re1,rho0
    call scmat(omg) !    得到 T 矩阵
    open(10,file='scat.dat')
    do m1 = -ne,ne
        qx = m1*pi/ne
        do m2 = -ne,ne
            qy = m2*pi/ne
            rho0 = 0
            do m3 = -ne,ne
                kx = pi*m3/ne
                do m4 = -ne,ne
                    ky = pi*m4/ne
                    call gf_clean(kx,ky,omg)
                    gf1 = gf0
                    call gf_clean(kx + qx,ky + qy,omg)
                    gf2 =gf0
                    gf2 = matmul(matmul(gf1,tmat),gf2)  ! 计算杂质系统的格林函数  
                    re1 = gf2(1,1) - conjg(gf2(1,1)) + gf2(3,3) - conjg(gf2(3,3))    ! 态密度计算(这里需要具体参考文章，或者是你想计算什么态密度)
                    rho0 = rho0 + re1*(1.0/ne)**2
                end do
            end do
            rho0 = rho0/(2.0*pi)**3*im
            write(10,*)qx/pi,qy/pi,real(rho0)
        end do
    end do
    close(10)
    return
    end subroutine main
!==================================================================
    subroutine scmat(omg)
    !   加入杂质后利用杂质矩阵和纯净系统的GF计算散射T-matrix
    use pub
    real kx,ky,omg
    integer m1,m2,m3,m4
    complex ga(N,N),temp1(N,N),temp2(N,N) ! green's function FT 后的结果
    !   green's function 傅里叶变化
    do m1 = -ne,ne
        do m2 = -ne,ne
            kx = m1*pi/ne
            ky = m2*pi/ne
            call matset(kx,ky)
            call gf_clean(kx,ky,omg)
            do m3 = 1,N
                do m4 = 1,N
                    ga(m3,m4) = ga(m3,m4) + gf0(m3,m4)*1.0/ne
                end do
            end do
        end do
    end do
    ga = ga/(2.0*pi)**2
    !-----------------------------------------------
    !   T-matrix compute
!    tmat = matmul(1 - matmul(V,ga),V)

    temp1 = 1 - matmul(V,ga)
    call inv(temp1,temp2)
    tmat = matmul(temp2,V)
    return
    end subroutine scmat
!==========================================================================
    subroutine gf_clean(kx,ky,omg)
    !   未加杂质时的格林函数
    use pub
    real kx,ky,omg
    integer m1,m2
    complex::temp1(N,N) = 0
    call matset(kx,ky)
    do m1 = 1,N
        do m2 = 1,N
            temp1(m1,m2) = omg - im*eta - ham(m1,m2)
        end do
    end do
    call inv(temp1,gf0)     ! 求逆的结果保存在gf0中,由哈密顿量得到其对应的格林函数
    return
    end subroutine gf_clean
!==========================================================================
    subroutine matset(kx,ky)
    use pub
    real kx,ky
    ham = 0

    ham(1,1) = (-2*t1*cos(kx) - 2*t2*cos(ky) - 4*t3*cos(kx)*cos(ky)) - mu
    ham(2,2) = -(-2*t1*cos(kx) - 2*t2*cos(ky) - 4*t3*cos(kx)*cos(ky)) + mu
    ham(3,3) = (-2*t1*cos(ky)-2*t2*cos(kx)-4*t3*cos(kx)*cos(ky))-mu
    ham(4,4) = -(-2*t1*cos(ky)-2*t2*cos(kx)-4*t3*cos(kx)*cos(ky)) + mu

    ham(1,2) = del0*cos(kx)*cos(ky)
    ham(2,1) = del0*cos(kx)*cos(ky)
    ham(1,3) = -4*t4*sin(kx)*sin(ky)
    ham(3,1) = -4*t4*sin(kx)*sin(ky)
    ham(1,4) = 0.0
    ham(4,1) = 0.0

    ham(2,3) = 0.0
    ham(3,2) = 0.0
    ham(2,4) = -(-4*t4*sin(kx)*sin(ky))
    ham(4,2) = -(-4*t4*sin(kx)*sin(ky))

    ham(3,4) = del0*cos(kx)*cos(ky)
    ham(4,3) = del0*cos(kx)*cos(ky)

    return
    end subroutine matset
!=======================矩阵求逆====================================
    subroutine inv(matin,matout)
    use pub
    complex,intent(in) :: matin(N,N)
    complex:: matout(size(matin,1),size(matin,2))
    real:: work2(size(matin,1))            ! work2 array for LAPACK
    integer::ipiv(size(matin,1))     ! pivot indices
    ! Store matin in matout to prevent it from being overwritten by LAPACK
    matout = matin
    ! SGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call CGETRF(N,N,matout,N,ipiv,info)
    if (info.ne.0) stop 'Matrix is numerically singular!'
    ! SGETRI computes the inverse of a matrix using the LU factorization
    ! computed by SGETRF.
    call CGETRI(N,matout,N,ipiv,work2,N,info)
    if (info.ne.0) stop 'Matrix inversion failed!'
    return
    end subroutine inv