#!/usr/bin/env julia

#
# This example is taken from Phys. Rev. D 106, L051502 (2022)
# Resonance peak + Continuum spectral function
#

using Random
using Printf
using ACFlow

heaviside(x::AbstractFloat) = ifelse(x < 0, zero(x), one(x))

function Î¾(Ï‰, ð�‘€, Î”)
    return 1.0 / ( 1.0 + exp( (ð�‘€ ^ 2.0 - Ï‰ ^ 2.0) / (Ï‰ * Î”) ) )
end

function Ï�áµ£(Ï‰, ð��¶, ð�‘€, Î“)
    return ð��¶ * Ï‰ ^ 2.0 / ( (Ï‰ ^ 2.0 / ð�‘€ / Î“ - ð�‘€ / Î“) ^ 2.0 + 1.0 )
end

function Ï�â‚œ(Ï‰, ð��¶, ð�‘€, Î²)
    ð�‘ƒâ‚� = ð��¶ * 3.0 * Ï‰ ^ 2.0 / 8.0 / Ï€
    ð�‘ƒâ‚‚ = heaviside(Ï‰ ^ 2.0 - 4.0 * ð�‘€ ^ 2.0)
    ð�‘ƒâ‚ƒ = tanh(Ï‰ * Î² / 4.0)
    ð�‘ƒâ‚„ = sqrt(abs(1.0 - (2.0 * ð�‘€ / Ï‰) ^ 2.0))
    ð�‘ƒâ‚… = 2.0 + (2.0 * ð�‘€ / Ï‰) ^ 2.0
    return ð�‘ƒâ‚� * ð�‘ƒâ‚‚ * ð�‘ƒâ‚ƒ * ð�‘ƒâ‚„ * ð�‘ƒâ‚…
end

# Setup parameters
wmin = +0.0  # Left boundary
wmax = +4.0  # Right boundary
nmesh = 2001 # Number of real-frequency points
niw  = 50    # Number of Matsubara frequencies
ntau = 501   # Number of imaginary time points
beta = 50.0  # Inverse temperature
ð�‘€áµ£   = 0.10  # Parameters for model
ð�‘€â‚œ   = 0.05
ð��¶áµ£   = 2.00
ð��¶â‚œ   = 0.00
Î“    = 0.06
Î”    = 1.00

#
# For true spectrum
#

# Real frequency mesh
rmesh = collect(LinRange(wmin, wmax, nmesh))

# Spectral function
image = similar(rmesh)
#
rmesh[1] = 1e-8 # To avoid NaN
for i in eachindex(rmesh)
    Ï�â‚� = Î¾(rmesh[i], ð�‘€áµ£, Î“) * Ï�áµ£(rmesh[i], ð��¶áµ£, ð�‘€áµ£, Î“) * (1.0 - Î¾(rmesh[i], ð�‘€áµ£ + Î“, Î“))
    Ï�â‚‚ = Î¾(rmesh[i], ð�‘€â‚œ + Î“, Î“) * Ï�â‚œ(rmesh[i], ð��¶â‚œ, ð�‘€â‚œ, beta)
    image[i] = Ï�â‚� + Ï�â‚‚
end
#
image = image ./ rmesh ./ rmesh ./ rmesh
rmesh[1] = 0.0

# Write spectral function
open("image.data", "w") do fout
    for i in eachindex(image)
        @printf(fout, "%20.16f %20.16f\n", rmesh[i], image[i])
    end
end

#
# For Matsubara frequency data
#

# Matsubara frequency mesh
iw = Ï€ / beta * (2.0 * collect(0:niw-1) .+ 0.0)

# Noise
seed = rand(1:100000000)
rng = MersenneTwister(seed)
noise_ampl = 1.0e-4
noise = randn(rng, F64, niw) * noise_ampl

# Kernel function
kernel = -2.0 * reshape(rmesh .^ 2.0, (1,nmesh)) ./
         (reshape(iw .^ 2.0, (niw,1)) .+ reshape(rmesh .^ 2.0, (1,nmesh)))
kernel[1,1] = -2.0

# Build green's function
KA = kernel .* reshape(image, (1,nmesh))
chiw = zeros(F64, niw)
for i in eachindex(chiw)
    chiw[i] = trapz(rmesh, KA[i,:]) + noise[i]
end

# Build error
err = ones(F64, niw) * noise_ampl

# Write green's function
open("chiw.data", "w") do fout
    for i in eachindex(chiw)
        z = chiw[i]
        @printf(fout, "%20.16f %20.16f %20.16f\n", iw[i], z, err[i])
    end
end

#
# For imaginary time data
#

# Imaginary time mesh
tmesh = collect(LinRange(0, beta, ntau))

# Noise
seed = rand(1:100000000)
rng = MersenneTwister(seed)
noise_ampl = 1.0e-4
noise = randn(rng, F64, ntau) * noise_ampl

# Build green's function
chit = zeros(F64, ntau)
for i = 1:ntau
    tw = exp.(-tmesh[i] * rmesh)
    bw = exp.(-beta * rmesh)
    btw = exp.(-(beta - tmesh[i]) * rmesh)
    K = rmesh .* (tw .+ btw) ./ (1.0 .- bw)
    K[1] = 2.0 / beta
    global KA = K .* image
    chit[i] = trapz(rmesh, KA) + noise[i]
end

# Build error
err = ones(F64, ntau) * noise_ampl

# Write green's function
open("chit.data", "w") do fout
    for i in eachindex(chit)
        @printf(fout, "%16.12f %16.12f %16.12f\n", tmesh[i], chit[i], err[i])
    end
end
