#!/usr/bin/env julia

#
# Well, this script will generate a pseudo-spectrum, which is then used
# to create a nonlinear mesh for the StochPX solver.
#
# Note that the parameters used to construct the pseudo-spectrum are just
# inspired by the analytic continuation results obtained by the MaxEnt
# and the NAC (using the Nevanlinna.jl package) methods.
#
# (1) Try to generate a pseudo-spectrum by `julia genspec.jl`.
#
# (2) Try to create a nonlinear mesh by `../../../util/gmesh.jl ac.toml`
#

using Random
using Printf
using ACFlow

# Setup parameters
wmin = +9.0  # Left boundary
wmax = +16.  # Right boundary
nmesh = 1001 # Number of real-frequency points
ğ�‘€â‚�   = 9.60  # Parameters for Gaussian mixture model
ğ�‘€â‚‚   = 11.5
Î“â‚�   = 0.01
Î“â‚‚   = 5.00
ğ��´â‚�   = 5.00
ğ��´â‚‚   = 1.80

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

# Spectral function
image = similar(rmesh)
#
for i in eachindex(rmesh)
    image[i] =            ğ��´â‚� * exp(-(rmesh[i] - ğ�‘€â‚�) ^ 2.0 / Î“â‚�)
    image[i] = image[i] + ğ��´â‚‚ * exp(-(rmesh[i] - ğ�‘€â‚‚) ^ 2.0 / Î“â‚‚)
end

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