import argparse
import astropy.io as aio
import galsim
import galsim.roman
import numpy as np
from astropy import constants as const
from astropy import units as u
from .utils.trapz import trapz
[docs]
fNuRef = 3.631e-23 * (u.W / u.m**2) / u.Hz # W/m^2/Hz
[docs]
def makeParser():
"""Input file parser."""
parser = argparse.ArgumentParser(description="Input Files for PSF Simulation")
parser.add_argument("wcsFileName", help="Pathname to Fits File for WCS")
parser.add_argument("starCat", help="Pathname to Star Catalog")
parser.add_argument("SCA", help="SCA Number")
parser.add_argument(
"-r",
"--random",
dest="randomPos",
action="store_true",
help="Use random positions for stars",
default=False,
)
parser.add_argument(
"-b",
"--blackbody",
dest="blackBody",
action="store_true",
help="Each star is assumed to be 5000 K blackbody for testing purposes",
default=False,
)
return parser
[docs]
def sedBB(w, T):
"""Blackbody distribution"""
return (
(8 * np.pi * const.h * const.c**2 / w**5) * 1 / (np.exp(const.h * const.c / (w * const.k_B * T)) - 1)
).decompose()
if __name__ == "__main__":
args = parser.parse_args()
# Read RA,Dec from star catalog
if ".fits" not in args.starCat:
raise Exception("Star Catalog should be a .fits file")
cat = galsim.Catalog(args.starCat)
degrees = galsim.AngleUnit(np.pi / 180)
readImage = galsim.fits.read(file_name=args.wcsFileName, hdu=1)
mybounds = readImage.bounds
mywcs = readImage.wcs
scaNum = int(args.SCA)
if scaNum < 10:
effAreaTable = aio.ascii.read(
f"Roman_effarea_tables_20240327/Roman_effarea_v8_SCA0{scaNum}_20240301.ecsv"
)
else:
effAreaTable = aio.ascii.read(
f"Roman_effarea_tables_20240327/Roman_effarea_v8_SCA{scaNum}_20240301.ecsv"
)
mirrorDiameter = 2.37 * u.m
geomArea = np.pi * mirrorDiameter**2 / 4
transmissionCurve = effAreaTable["F184"] * u.m**2 / geomArea
tExp = 120 * u.s
outImage = galsim.Image(wcs=mywcs, bounds=mybounds)
roman_bandpasses = galsim.roman.getBandpasses()
# The following PSF will be used if not using random positions for stars. Otherwise thePSF will be
# slightly different depending on the SCA position of the object's image. Will modify this later.
psf = galsim.roman.getPSF(scaNum, "F184", wcs=mywcs, wavelength=roman_bandpasses["F184"])
source = galsim.Convolve([psf, galsim.DeltaFunction(flux=1.0)])
for i in range(cat.nobjects):
if not args.randomPos:
ra = cat.get(i, "RAJ2000") * degrees
dec = cat.get(i, "DECJ2000") * degrees
mag = cat.get(i, "H")
worldCenter = galsim.CelestialCoord(ra=ra, dec=dec)
imageCenter = mywcs.posToImage(worldCenter)
else:
x = np.random.random_sample() * mybounds.getXMax()
y = np.random.random_sample() * mybounds.getYMax()
mag = cat.get(i, "H")
imageCenter = galsim.PositionD(x=x, y=y)
pos_SCA = galsim.PositionD(x=x - (mybounds.getXMax() / 2.0), y=y - (mybounds.getYMax() / 2.0))
psf = galsim.roman.getPSF(
scaNum, "F184", SCA_pos=pos_SCA, wcs=mywcs, wavelength=roman_bandpasses["F184"]
)
source = galsim.Convolve([psf, galsim.DeltaFunction(flux=1.0)])
if args.blackBody:
# bb = BlackBody(temperature = 5000*u.K)
wav = np.arange(0.400, 2.600, 0.001) * u.um
fluxUnnorm = sedBB(wav, 5000 * u.K)
fLambdaRef = fNuRef * const.c / wav**2
norm = (
10 ** (-0.4 * mag)
* trapz(fLambdaRef * transmissionCurve * wav, x=wav)
/ trapz(fluxUnnorm * transmissionCurve * wav, x=wav)
)
flux = norm * fluxUnnorm
nPhotQ = trapz(flux * effAreaTable["F184"] * u.m**2 * wav * tExp / (const.h * const.c), x=wav)
nPhotQ = nPhotQ.decompose()
# nPhotQ = nPhotQ.to(1/(u.s*u.cm**2))
nPhot = nPhotQ.value
print(nPhot)
else:
nPhot = 300000 # Need to update as a function of SED for a given star
geomAreaCM = geomArea.to(u.cm**2).value
source.drawImage(outImage, method="phot", n_photons=nPhot, center=imageCenter, add_to_image=True)
print("Image Drawn!")
outImage.write("outImage.fits")
print("Image written to outImage.fits")