Commit 8f63aba4 authored by Wei Chengliang's avatar Wei Chengliang
Browse files

add tools to getPSF

parent 2cd2bbb9
Loading
Loading
Loading
Loading

tools/getPSF.py

0 → 100644
+82 −0
Original line number Diff line number Diff line
import os
import numpy as np
import ObservationSim.PSF.PSFInterp as PSFInterp
from ObservationSim.Instrument import Chip, Filter, FilterParam
import yaml
import galsim
import astropy.io.fits as fitsio

# Setup PATH
SIMPATH = "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C8/testRun_FGS"
config_filename= SIMPATH+"/config_C6_fits.yaml"
cat_filename   = SIMPATH+"/MSC_00000000/MSC_10106100000000_chip_40_filt_FGS.cat"

# Read cat file
catFn = open(cat_filename,"r")
line = catFn.readline()
print(cat_filename,'\n',line)
imgPos = []
chipID = -1
for line in catFn:
    line = line.strip()
    columns = line.split()
    
    if chipID == -1:
        chipID = int(columns[1])
    else:
        assert chipID == int(columns[1])
    ximg = float(columns[3])
    yimg = float(columns[4])
    imgPos.append([ximg, yimg])
imgPos = np.array(imgPos)
nobj = imgPos.shape[0]
print('chipID, nobj::', chipID, nobj)

# Read config file
with open(config_filename, "r") as stream:
    try:
        config = yaml.safe_load(stream)
        for key, value in config.items():
            print (key + " : " + str(value))
    except yaml.YAMLError as exc:
        print(exc)

# Setup Chip
chip = Chip(chipID=chipID, config=config)

for iobj in range(nobj):
    print("\nget psf for iobj-", iobj, '\t', 'bandpass:', end=" ", flush=True)
    # Setup Position
    x, y = imgPos[iobj, :] # try get the PSF at some location (1234, 1234) on the chip
    pos_img = galsim.PositionD(x, y)
    
    # Setup sub-bandpass
    # (There are 4 sub-bandpasses for each PSF sample)
    filter_param = FilterParam()
    filter_id, filter_type = chip.getChipFilter()
    filt = Filter(
        filter_id=filter_id, 
        filter_type=filter_type, 
        filter_param=filter_param, 
        ccd_bandpass=chip.effCurve)
    bandpass_list = filt.bandpass_sub_list
    for i in range(len(bandpass_list)):
        print(i, end=" ", flush=True)
        bandpass = bandpass_list[i] # say you want to access the PSF for the sub-bandpass at the blue end for that chip
        
        # Get corresponding PSF model
        psf_model = PSFInterp(chip=chip, npsf=100, PSF_data_file=config["psf_setting"]["psf_dir"])
        psf  = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, galsimGSObject=False)
        
        if True:
            fn = "psf_{:}.{:}.{:}.fits".format(chipID, iobj, i)
            if fn != None:
                if os.path.exists(fn):
                    os.remove(fn)
            hdu = fitsio.PrimaryHDU()
            hdu.data = psf
            hdu.header.set('pixScale', 5)
            hdu.writeto(fn)