# 构造Josephson junction
@everywhere using SharedArrays, LinearAlgebra,Distributed,DelimitedFiles,Printf,Arpack
# --------------------------------------
@everywhere function pauli()
    s0 = zeros(ComplexF64,2,2)
    s1 = zeros(ComplexF64,2,2)
    s2 = zeros(ComplexF64,2,2)
    s3 = zeros(ComplexF64,2,2)
    #----
    s0[1,1] = 1
    s0[2,2] = 1
    #----
    s1[1,2] = 1
    s1[2,1] = 1
    #----
    s2[1,2] = -im
    s2[2,1] = im
    #-----
    s3[1,1] = 1
    s3[2,2] = -1
    #-----
    return s0,s1,s2,s3
end
#---------------------------------------------------------------
@everywhere function boundary(xn::Int64)
    bry = zeros(Int64,2,xn)
    for i0 in 1:xn
        bry[1,i0] = i0 + 1
        if i0 == xn
            bry[1,i0] = bry[1,i0] - xn
        end
        bry[2,i0] = i0 - 1
        if i0 == 1
            bry[2,i0] = bry[2,i0] + xn
        end
    end
    return bry
end
#-------------------------------------------------------------------
@everywhere function junction(ix::Int64,phil::Float64,phir::Float64,lpos::Int64,rpos::Int64)
    # 两侧超导位相存在差
    if ix <= lpos
        return exp(im*phil)
    elseif ix >= rpos
        return exp(im*phir)
    else
        return 0 # 正常态区域不存在电子配对，设置为零
    end
end
#-----------------------------------------------------------------------------------------
@everywhere function orpar1(ix::Int64,phil::Float64,phir::Float64,lpos::Int64,rpos::Int64)::Matrix{ComplexF64}
    # dx^2-y^2 pairing
    phi = junction(ix,phil,phir,lpos,rpos)
    g1 = zeros(ComplexF64,4,4)
    g1[1,2] = -im*phi
    g1[2,1] = conj(g1[1,2])
    return g1
end
#-----------------------------------------------------------------------------------------
@everywhere function matset(phil::Float64,phir::Float64,lpos::Int64,rpos::Int64)
    t0::Float64 = 1.0
    mu::Float64 = 1.0
    d0::Float64 = 1.0
    hn::Int64 = 2
    xn::Int64 = 100  # 开边界格点数量
    N::Int64 = xn*hn
    ham = zeros(ComplexF64,N,N)
    s0,sx,sy,sz = pauli()
    bry = boundary(xn)
    #------------------------
    for i0 in 1:xn
        g1 = orpar1(i0,phil,phir,lpos,rpos) # 需要对序参量进行修正
        for i1 in 0:hn - 1,i2 in 0:hn - 1
            ham[i0 + i1*xn,i0 + i2*xn] = -mu*sz[i1 + 1,i2 + 1]
            if i0 != xn
                ham[i0 + i1*xn,bry[1,i0] + i2*xn] = -t0*sz[i1 + 1,i2 + 1] + d0/(2*im)*g1[i1 + 1,i2 + 1]
            end
            if i0 != 1
                ham[i0 + i1*xn,bry[2,i0] + i2*xn] = -t0*sz[i1 + 1,i2 + 1] - d0/(2*im)*g1[i1 + 1,i2 + 1]
            end
        end
    end
    #-----------------------------------
    # 修正边界hopping
    for ix in 1:xn
        i0 = ix
        for i1 in 0:hn -1,i2 in 0:hn - 1 
            if ix == lpos
                ham[i0 + xn * i1,bry[1,i0] + xn * i2] = -t0*sz[i1 + 1,i2 + 1]
            elseif ix == rpos
                ham[i0 + xn * i1,bry[2,i0] + xn * i2] = -t0*sz[i1 + 1,i2 + 1]
            end 
        end
    end
    #-----------------
    if ~ishermitian(ham)
        f1 = open("hermi.dat","w")
        for m1 in 1:N,m2 in 1:N
            if ham[m1,m2] != conj(ham[m2,m1])
                # println("(",m1,",",m2,")",ham[m1,m2],ham[m2,m1])
                writedlm(f1,[real(m1) real(m2) ham[m1,m2] ham[m2, m1]],"\t")
            end
        end
        close(f1)
    end
    #-----------------------------------------
    if ishermitian(ham)
        temp2 = (a->(@sprintf "%3.1f" a)).(phir/pi)
        fx1 = "eigval-phase-" * temp2 * ".dat"
        #f1 = open(fx1,"w")
        val,vec = eigen(ham)
        # val,vec = eigs(ham,nev = 50,maxiter = 30,which = :SM)
        ind = (a->(@sprintf "%5.2f" a)).(range(1,length(val),length = length(val)))
        val2 = (a->(@sprintf "%15.8f" a)).(sort(map(real,val)))
        #writedlm(f1,[ind val2],"\t")
        #close(f1)
    else
        println("Hamiltonian is not hermitian")
        # break
    end
    return sort(map(real,val))
end
#--------------------------------------------------------------------------------
@everywhere function phase()
    philist = 0:0.1:4
    vallist = SharedArray(zeros(Float64,length(philist),200))
    @sync @distributed for i0 in 1:length(philist)
        vallist[i0,:] = matset(0.0,philist[i0]*pi,50,51)
    end
    fx1 = "KitaevChain-short.dat"
    f1 = open(fx1,"w")
    x0 = (a->(@sprintf "%15.8f" a)).(philist)
    y0 = (a->(@sprintf "%15.8f" a)).(vallist)
    writedlm(f1,[x0 y0],"\t")
end
#-----------------------------------------------------------------------------------
@everywhere function current()
    kbT::Float64 = 0.001
    dphi::Float64 = 0.1
    philist = 0:dphi:4
    vallist = SharedArray(zeros(Float64,length(philist),200))
    @sync @distributed for i0 in 1:length(philist)
        vallist[i0,:] = matset(0.0,philist[i0]*pi,50,51)
    end
    len1 = length(vallist[1,:]) # 只需要能量为正的本征值
    I0 = []
    phi = []
    for i0 in 1:length(philist) - 1 # loop for phase
        re1 = 0
        for i1 in Int(len1/2) + 1:len1 # loop for eigvals
            re1 += tanh(1/kbT*vallist[i0,i1]/2)*(vallist[i0 + 1,i1] - vallist[i0,i1])/(dphi*pi)
        end
        append!(I0,re1)
        append!(phi,philist[i0])
    end
    fx1 = "KitaevChain-current.dat"
    f1 = open(fx1,"w")
    x0 = (a->(@sprintf "%15.8f" a)).(phi)
    y0 = (a->(@sprintf "%15.8f" a)).(I0)
    writedlm(f1,[x0 y0],"\t")
    # writedlm(f1,[y0],"\t")
    close(f1)
end
#------------------------------------------------------------------------------------
@time phase()
@time current()