%matplotlib inline
%reload_ext autoreload
%autoreload 2
import os
import json
import fpfs
import galsim
import matplotlib.pylab as plt
import astropy.io.fits as pyfits
from configparser import ConfigParser
from astropy.visualization import simple_norm

def prepare_psf(pscale, seeing, psf_type, outdir):
    psffname = os.path.join(outdir, "psf-%d.fits" % (seeing * 100))
    if psf_type.lower() == "moffat":
        psfInt = galsim.Moffat(
            beta=3.5, fwhm=seeing, trunc=seeing * 4.0
        ).shear(e1=0.02, e2=-0.02)
    else:
        raise ValueError("Only support moffat PSF.")
    psfImg = psfInt.drawImage(nx=45, ny=45, scale=pscale)
    psfImg.write(psffname)
    return psfInt

Id = 1
config_name = 'config_sim_gal.ini'
cparser = ConfigParser()
cparser.read(config_name)
simname = cparser.get("simulation", "sim_name")
imgdir = cparser.get("simulation", "img_dir")
infname = cparser.get("simulation", "input_name")
scale = cparser.getfloat("survey", "pixel_scale")
image_nx = cparser.getint("survey", "image_nx")
image_ny = cparser.getint("survey", "image_ny")
assert image_ny == image_nx, "'image_nx' must equals 'image_ny'!"
if "basic" in simname or "small" in simname:
    assert image_nx % 256 == 0, "'image_nx' must be divisible by 256 ."
outdir = os.path.join(imgdir, simname)
if not os.path.exists(outdir):
    os.makedirs(outdir, exist_ok=True)

seeing = cparser.getfloat("survey", "psf_fwhm")
psfInt = prepare_psf(scale, seeing, psf_type="moffat", outdir=outdir)

glist = []
if cparser.getboolean("distortion", "test_g1"):
    glist.append("g1")
if cparser.getboolean("distortion", "test_g2"):
    glist.append("g2")
if len(glist) > 0:
    zlist = json.loads(cparser.get("distortion", "shear_z_list"))
    pendList = ["%s-%s" % (i1, i2) for i1 in glist for i2 in zlist]
else:
    # this is for non-distorted image simulation
    pendList = ["g1-1111"]
shear_value = cparser.getfloat("distortion", "shear_value")



p2List = ["0000", "2222"]
p1List = ["g1"]
pendList = ["%s-%s" % (i1, i2) for i1 in p1List for i2 in p2List]
pp = pendList[0]

sim_img = fpfs.simutil.make_isolate_sim(
    gal_type="mixed",
    psf_obj=psfInt,
    gname=pp,
    seed=Id,
    ny=image_ny//2,
    nx=image_nx,
    scale=scale,
    do_shift=False,
    shear_value=shear_value,
    nrot_per_gal=2,
)[0]
plt.imshow(sim_img,aspect='equal',cmap='RdYlBu_r',origin='lower',interpolation='None',\
             norm=simple_norm(sim_img,'asinh',asinh_a=0.1,min_cut=-0.01,max_cut=0.8))

2024/10/10 20:03:58 ---  Creating Mixed galaxy profiles
2024/10/10 20:03:58 ---  Processing for g1-0000, and shear is -0.03.
WARNING: AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead. [warnings]
2024/10/10 20:04:00 ---  AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead.
WARNING: AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead. [warnings]
2024/10/10 20:04:00 ---  AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead.

<matplotlib.image.AxesImage at 0x7fbcd1caf790>

sim_img2 = fpfs.simutil.make_isolate_sim(
    gal_type="mixed",
    psf_obj=psfInt,
    gname=pp,
    seed=Id,
    ny=image_ny//2,
    nx=image_nx,
    scale=scale,
    do_shift=False,
    shear_value=shear_value,
    nrot_per_gal=2,
    sim_method="mc",
)[0]
plt.imshow(sim_img2,aspect='equal',cmap='RdYlBu_r',origin='lower',interpolation='None',\
             norm=simple_norm(sim_img,'asinh',asinh_a=0.1,min_cut=-0.01,max_cut=0.8))

2024/10/10 20:04:00 ---  Creating Mixed galaxy profiles
2024/10/10 20:04:00 ---  Processing for g1-0000, and shear is -0.03.
WARNING: AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead. [warnings]
2024/10/10 20:04:01 ---  AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead.
WARNING: AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead. [warnings]
2024/10/10 20:04:01 ---  AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead.

<matplotlib.image.AxesImage at 0x7fbcd1b6b280>

plt.imshow(sim_img*0.70+sim_img2*0.3,aspect='equal',cmap='RdYlBu_r',origin='lower',interpolation='None',\
             norm=simple_norm(sim_img,'asinh',asinh_a=0.1,min_cut=-0.01,max_cut=0.8))

WARNING: AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead. [warnings]
2024/10/10 20:04:01 ---  AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead.
WARNING: AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead. [warnings]
2024/10/10 20:04:01 ---  AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead.

<matplotlib.image.AxesImage at 0x7fbccfabc4c0>