csst_mci_sim.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 18 09:38:35 2021


@author: yan, zhaojun"""


"""
The CSST MCI Image Simulator
=============================================

This file contains an image simulator for the CSST MCI.

The approximate sequence of events in the simulator is as follows:

      #. Read in a configuration file, which defines for example,
         detector characteristics (bias, dark and readout noise, gain,
         plate scale and pixel scale, exposure time etc.).
      #. Read in another file containing charge trap definitions (for CTI modelling).
      #. Read in a file defining the cosmic rays (trail lengths and cumulative distributions).
      

      #. Loop over the number of exposures to co-add and for each object in the object catalog:



      #. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [optional].

      #. Add cosmic ray tracks onto the CCD with random positions but known distribution [optional].
      #. Apply detector charge bleeding in column direction [optional].
      #. Add constant dark current  [optional].

      #. Add photon (Poisson) noise [optional]
      #. Add cosmetic defects from an input file [optional].
      #. Add overscan regions in the serial direction [optional].
      #. Apply the CDM03 radiation damage model [optional].
      #. Apply non-linearity model to the pixel data [optional].
      #. Add readout noise selected [optional].
      #. Convert from electrons to ADUs using a given gain factor.
      #. Add a given bias level and discretise the counts (the output is going to be in 16bit unsigned integers).
      #. Finally the simulated image is converted to a FITS file, a WCS is assigned
         and the output is saved to the current working directory.

.. Warning:: The code is still work in progress and new features are being added.
             The code has been tested, but nevertheless bugs may be lurking in corners, so
             please report any weird or inconsistent simulations to the author.


Contact Information: zhaojunyan@shao.ac.cn
    
-------

"""
import galsim
import pandas as pd
from scipy.integrate import simps

from datetime import datetime, timedelta
from astropy.time import Time
from astropy.coordinates import get_sun

import scipy.io as sio
import cmath

import numpy as np

from tqdm import tqdm

from astropy.table import Table

from astropy.wcs import WCS as WCS

from  astropy.io import fits

from astropy import units as u

import os, sys, math 

import configparser as ConfigParser

from optparse import OptionParser

from matplotlib import pyplot as plt 

from scipy import ndimage

sys.path.append('./csst_mci_sim')
##sys.path.append('..')

from CTI import CTI
from support import logger as lg
from support import cosmicrays

from support import shao

from support import sed 
#from shao import onOrbitObsPosition

from support import MCIinstrumentModel

from joblib import Parallel, delayed
#import multiprocessing

from astropy.coordinates import SkyCoord

from scipy import interpolate

from scipy.signal import fftconvolve

import pandas

# sys.path.append('../MCI_inputData/TianCe')
# sys.path.append('../MCI_inputData/SED_Code')




 

# import jax
#import jax.numpy as jnp

# from jax import config
# config.update("jax_enable_x64", True)

# os.environ['CUDA_VISIBLE_DEVICES'] = '0'
# devices = jax.local_devices()


# set the folder
FOLDER ='../'


#####################################################################################
###############################################################################
import ctypes

import astropy.coordinates as coord

from scipy.interpolate import interp1d

#from sky_bkg import sky_bkg


filterPivotWave = {'nuv':2875.5,'u':3629.6,'g':4808.4,'r':6178.2, 'i':7609.0, 'z':9012.9,'y':9627.9}
filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6}
# filterCCD = {'nuv':'UV0','u':'UV0','g':'Astro_MB','r':'Astro_MB', 'i':'Basic_NIR', 'z':'Basic_NIR','y':'Basic_NIR'}
# bandRange = {'nuv':[2504.0,3230.0],'u':[3190.0,4039.0],'g':[3989.0,5498.0],'r':[5438.0,6956.0], 'i':[6886.0,8469.0],
#              'z':[8379.0,10855.0],'y':[9217.0, 10900.0], 'GU':[2550, 4000],'GV':[4000, 6200],'GI':[6200,10000]}
# Instrument_dir = '/Users/linlin/Data/csst/straylightsim-master/Instrument/'
# SpecOrder = ['-2','-1','0','1','2']
#
# filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8}


def transRaDec2D(ra, dec):
    # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz.
    x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795)
    y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795)
    z1 = np.sin(dec / 57.2957795)
    return np.array([x1, y1, z1])
###############################################################################

def flux2ill(wave, flux):

    # erg/s/cm^2/A/arcsec^2 to W/m^2

    # 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A
    # 1 sr = 1 rad^2 = 4.25452e10 arcsec^2
    # 1 J/s = 1 W
    # 1 J = 10^7 erg

    # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
    flux1 = flux / (1/4.25452e10)
    # convert erg/s/cm^2/A/sr to W/m^2/sr/um
    flux2 = flux1 * 10

    # 对接收面积积分，输出单位 W/m^2/nm
    D = 2   # meter
    f = 28  # meter
    flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3

    # # # u_lambda: W/m^2/nm
    # # obj0:
    # V = 73 * 1e-8  # W/m^2/nm
    # Es = V * np.pi * D**2 / 4 / f**2 / 10**3
    # c1 = 3.7418e-16
    # c2 = 1.44e-2
    # t = 5700
    # # wave需要代入meter.
    # wave0 = np.arange(380, 780)  # nm
    # delta_lamba0 = 1             # nm
    # u_lambda = c1 / (wave0*1e-9)**5 / (np.exp(c2 / (wave0*1e-9) / t) - 1)
    # f_lambda = (u_lambda / u_lambda[wave0 == 500]) * Es
    # E0 = np.sum(f_lambda * 1)
    # # plt.plot(wave, u_lambda)
    # # plt.show()

    # 对波长积分
    f = interp1d(wave, flux3)
    wave_interp = np.arange(3800, 7800)
    flux3_interp = f(wave_interp)
    # 输出单位 W/m^2
    delta_lamba = 0.1   # nm
    E = np.sum(flux3_interp * delta_lamba)

    # pdb.set_trace()

    return E

################################################################
def ill2flux(E):

    # use template from sky_bkg (background_spec_hst.dat)
    filename = self.information['dir_path']+'MCI_inputData/refs/background_spec_hst.dat'
    cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#')
    wave0 = cat_spec[0].values       # A
    spec0 = cat_spec[2].values      # erg/s/cm^2/A/arcsec^2

    # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
    flux1 = spec0 / (1/4.25452e10)
    # convert erg/s/cm^2/A/sr to W/m^2/sr/um
    flux2 = flux1 * 10

    # 对接收面积积分，输出单位 W/m^2/nm
    D = 2   # meter
    f = 28  # meter, 焦距，转换关系来源于王维notes.
    flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3

    f = interp1d(wave0, flux3)
    wave_range = np.arange(3800, 7800)
    flux3_mean = f(wave_range)
    delta_lamba = 0.1   # nm
    E0 = np.sum(flux3_mean * delta_lamba)

    factor = E / E0
    spec_scaled = factor * spec0

    return wave0, spec_scaled

##############################################################

##########################################################
def zodiacal(ra, dec, time):
    """
        For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998.

    :param ra: RA in unit of degree, ICRS frame
    :param dec: DEC in unit of degree, ICRS frame
    :param time: the specified string that in ISO format i.e., yyyy-mm-dd.
    :return:
        wave_A: wavelength of the zodical spectrum
        spec_mjy: flux of the zodical spectrum, in unit of MJy/sr
        spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr

    """

    # get solar position
    dt = datetime.fromisoformat(time)
    #jd = julian.to_jd(dt, fmt='jd')
    jd = time2jd(dt)
    t = Time(jd, format='jd', scale='utc')

    astro_sun = get_sun(t)
    ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg

    radec_sun = SkyCoord(ra=ra_sun*u.degree, dec=dec_sun*u.degree, frame='gcrs')
    lb_sun = radec_sun.transform_to('geocentrictrueecliptic')

    # get offsets between the target and sun.
    radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
    lb_obj = radec_obj.transform_to('geocentrictrueecliptic')

    beta = abs(lb_obj.lat.degree)
    lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree)

    # interpolated zodical surface brightness at 0.5 um
    zodi = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/zodi_map.dat', sep='\s+', header=None, comment='#')
    beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
    lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45,
                          60, 75, 90, 105, 120, 135, 150, 165, 180])
    xx, yy = np.meshgrid(beta_angle, lamda_angle)
    #xx, yy = np.meshgrid(beta_angle, lamda_angle,indexing='ij', sparse=True)
    
    f = interpolate.interp2d(xx, yy, zodi, kind='linear')
    #f = interpolate.RegularGridInterpolator((xx, yy), zodi, method='linear')
    
    zodi_obj = f(beta, lamda)       # 

    # read the zodical spectrum in the ecliptic
    cat_spec = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
    wave = cat_spec[0].values       # A
    spec0 = cat_spec[1].values      # 
    zodi_norm = 252                 # 

    spec = spec0 * (zodi_obj / zodi_norm) * 1e-8  # 

    # convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr
    wave_A = wave                                               # A
    #spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6      # MJy/sr
    spec_erg = spec * 0.1                                       # erg/s/cm^2/A/sr
    spec_erg2 = spec_erg / 4.25452e10                           # erg/s/cm^2/A/arcsec^2
    
    # self.zodiacal_wave=wave_A    # in A
    
    # self.zodiacal_flux=spec_erg2

    return  wave_A, spec_erg2
    ###################################################################################
    
    

#from astropy import units as u
#from astropy.coordinates import SkyCoord


def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):

    ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180
    dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180
    radec_sat = SkyCoord(ra=ra_sat*u.degree, dec=dec_sat*u.degree, frame='gcrs')
    lb_sat = radec_sat.transform_to('geocentrictrueecliptic')

    # get the obj location
    radec_obj = SkyCoord(ra=ra_obj*u.degree, dec=dec_obj*u.degree, frame='gcrs')
    lb_obj = radec_obj.transform_to('geocentrictrueecliptic')

    # calculate the angle between sub-satellite point and the earth side
    earth_radius = 6371     # km
    sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2)
    angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180

    # calculate the angle between satellite position and the target position
    angle_b = lb_sat.separation(lb_obj)

    # calculat the earth angle
    angle = 180 - angle_a - angle_b.degree

    return angle


    
    #################################################3
    
def MCIinformation():
    """
    Returns a dictionary describing MCI CCD. The following information is provided (id: value - reference)::

        
         dob: 0 - CDM03 (Short et al. 2010)
         fwc: 90000 - CCD spec EUCL-EST-RS-6-002 (for CDM03)      
         rdose: 30000000000.0 - derived (above the PLM requirement)
         sfwc: 730000.0 - CDM03 (Short et al. 2010), see also the CCD spec EUCL-EST-RS-6-002       
         st: 5e-06 - CDM03 (Short et al. 2010)        
         svg: 1e-10 - CDM03 (Short et al. 2010)
         t: 0.01024 - CDM03 (Short et al. 2010)
         trapfile: cdm_euclid.dat - CDM03 (derived, refitted to CCD204 data)
         vg: 6e-11 - CDM03 (Short et al. 2010)
         vth: 11680000.0 - CDM03 (Short et al. 2010)



    :return: instrument model parameters
    :rtype: dict
    """

    #########################################################################################################
    out=dict()
    out.update({'dob' : 0, 'rdose' : 8.0e9,
                'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat',
                'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11,
                'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10})

    return out
#######################################

def CCDnonLinearityModel(data, beta=6e-7):
    """   

    The non-linearity is modelled based on the results presented. 
    :param data: data to which the non-linearity model is being applied to
    :type data: ndarray

    :return: input data after conversion with the non-linearity model
    :rtype: float or ndarray
    """
    out = data-beta*data**2
    
    return out
#############################################################33333

class StrayLight(object):

    def __init__(self, jtime = 2460843., sat = np.array([0,0,0]), radec = np.array([0,0])):
        self.jtime = jtime
        self.sat = sat
        self.equator = coord.SkyCoord(radec[0]*u.degree, radec[1]*u.degree,frame='icrs')
        self.ecliptic = self.equator.transform_to('barycentrictrueecliptic')
        self.pointing = transRaDec2D(radec[0], radec[1])
        self.slcdll = ctypes.CDLL(self.information['dir_path']+'MCI_inputData/refs/libstraylight.so') #dylib

        self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double),
                                         ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
                                         ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]

        self.slcdll.PointSource.argtypes = [ctypes.c_double ,ctypes.POINTER(ctypes.c_double),
                                           ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double),
                                           ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]

        self.slcdll.EarthShine.argtypes = [ctypes.c_double ,ctypes.POINTER(ctypes.c_double),
                                          ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
                                          ctypes.POINTER(ctypes.c_double)]

        self.slcdll.Zodiacal.argtypes = [ctypes.c_double ,ctypes.POINTER(ctypes.c_double),
                                        ctypes.POINTER(ctypes.c_double)]
        self.slcdll.ComposeY.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
                                       ctypes.POINTER(ctypes.c_double)]
        self.slcdll.Init.argtypes=[ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p]
        self.deFn = "../MCI_inputData/refs/DE405"
        self.PSTFn = "../MCI_inputData/refs/PST"
        self.RFn = "../MCI_inputData/refs/R"
        self.ZolFn = "../MCI_inputData/refs/Zodiacal"
        self.brightStarTabFn = "../MCI_inputData/refs/BrightGaia_with_csst_mag"

        self.slcdll.Init(str.encode(self.deFn),str.encode(self.PSTFn),str.encode(self.RFn),str.encode(self.ZolFn))


    def caculateStarLightFilter(self, filter='i'):

        sat = (ctypes.c_double*3)()
        sat[:] = self.sat
        ob = (ctypes.c_double*3)()
        ob[:] = self.pointing

        py1 = (ctypes.c_double*3)()
        py2 = (ctypes.c_double*3)()
        self.slcdll.ComposeY(ob, py1, py2)

        star_e1 = (ctypes.c_double*7)()
        self.slcdll.PointSource(self.jtime, sat, ob, py1, star_e1, str.encode(self.brightStarTabFn))
        star_e2 = (ctypes.c_double*7)()
        self.slcdll.PointSource(self.jtime, sat, ob, py2, star_e2, str.encode(self.brightStarTabFn))

        band_star_e1 = star_e1[:][filterIndex[filter]]
        band_star_e2 = star_e2[:][filterIndex[filter]]

        return max(band_star_e1, band_star_e2)

    def caculateEarthShineFilter(self, filter='i'):
        sat = (ctypes.c_double*3)()
        sat[:] = self.sat
        ob = (ctypes.c_double*3)()
        ob[:] = self.pointing

        py1 = (ctypes.c_double*3)()
        py2 = (ctypes.c_double*3)()
        self.slcdll.ComposeY(ob, py1, py2)

        earth_e1 = (ctypes.c_double*7)()
        self.slcdll.EarthShine(self.jtime, sat, ob, py1, earth_e1)          # e[7]代表7个波段的照度
        
        earth_e2 = (ctypes.c_double*7)()
        self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2)

        band_earth_e1 = earth_e1[:][filterIndex[filter]]
        band_earth_e2 = earth_e2[:][filterIndex[filter]]

        return max(band_earth_e1, band_earth_e2)
    
###############################################################################  
####################################################################################


def processArgs(printHelp=False):
    """
    Processes command line arguments.
    """
    parser = OptionParser()

    parser.add_option('-c', '--configfile', dest='configfile',
                      help="Name of the configuration file", metavar="string")
    
    parser.add_option('-s', '--section', dest='section',
                      help="Name of the section of the config file [SCIENCE]", metavar="string")
    
    parser.add_option('-q', '--quadrant', dest='quadrant', help='CCD quadrant to simulate [0, 1, 2, 3]',
                      metavar='int')
    
    parser.add_option('-x', '--xCCD', dest='xCCD', help='CCD number in X-direction within the FPA matrix',
                      metavar='int')
    
    parser.add_option('-y', '--yCCD', dest='yCCD', help='CCD number in Y-direction within the FPA matrix',
                      metavar='int')
    
    parser.add_option('-d', '--debug', dest='debug', action='store_true',
                      help='Debugging mode on')
    
    parser.add_option('-t', '--test', dest='test', action='store_true',
                      help='Run unittest')
    
    parser.add_option('-f', '--fixed', dest='fixed', help='Use a fixed seed for the random number generators',
                      metavar='int')
    if printHelp:
        parser.print_help()
    else:
        return parser.parse_args()


###############################################################################

def make_c_coor(fov, step):
    """
    Draw the mesh grids for a fov * fov box with given step .
    """
    
    nc=int(fov/step)
    
    bs=fov
    ds = bs / nc
    xx01 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds
    xx02 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds
    xg2, xg1 = np.meshgrid(xx01, xx02)
    return xg1, xg2

##############################################################################    

###############################################################################

def str2time(strTime):
    if len(strTime)>20:#
        msec=int(float('0.'+strTime[20:])*1000000) # microsecond
    str2=strTime[0:19]+' '+str(msec)
    return datetime.strptime(str2,'%Y %m %d %H %M %S %f')
    
#datetime to mjd
def time2mjd(dateT):
    t0=datetime(1858,11,17,0,0,0,0)   #
    mjd=(dateT-t0).days
    mjd_s=dateT.hour*3600.0+dateT.minute*60.0+dateT.second+dateT.microsecond/1000000.0
    return mjd+mjd_s/86400.0

def time2jd(dateT):
    t0=datetime(1858,11,17,0,0,0,0)  # 
    mjd=(dateT-t0).days
    mjd_s=dateT.hour*3600.0+dateT.minute*60.0+dateT.second+dateT.microsecond/1000000.0
    return mjd+mjd_s/86400.0++2400000.5
 
#mjd to datetime
def mjd2time(mjd):
    t0=datetime(1858,11,17,0,0,0,0) # 
    return t0+timedelta(days=mjd)

def jd2time(jd):
    mjd=jd2mjd(jd)
    return mjd2time(mjd)
    
#mjd to jd
def mjd2jd(mjd):
    return mjd+2400000.5

def jd2mjd(jd):
    return jd-2400000.5

def dt2hmd(dt):
    ## dt is datetime 
    hour=dt.hour
    minute=dt.minute
    second=dt.second
    if hour<10:
        str_h='0'+str(hour)
    else:
        str_h=str(hour)
        
    if minute<10:
        str_m='0'+str(minute)
    else:
        str_m=str(minute)
    
    if second<10:
        str_d='0'+str(second+dt.microsecond*1e-6)
    else:
        str_d=str(second+dt.microsecond*1e-6)        
    
    return str_h+':'+str_m+':'+str_d

###############################################################################


def deg2HMS(ra0, dec0):
    '''convert deg to ra's HMS and  dec's DMS'''
    c = SkyCoord(ra=ra0*u.degree, dec=dec0*u.degree)
    ss=c.to_string('hmsdms')
    print(ss)
    for k in range(5,len(ss)-1):

        if ss[k]=='s':
            s_idx=k
            
        if ss[k]=='d':
            d_idx=k
        
        if ss[k]=='m':
            m_idx=k
            
    if c.ra.hms.s<10 and  c.ra.hms.s>0.1:
        temp=str(c.ra.hms.s)
        ra_ss='0'+temp[:3]
    if c.ra.hms.s>10:
        temp=str(c.ra.hms.s)
        ra_ss=temp[:4]
        
    if c.ra.hms.s<0.1:
        temp=str(c.ra.hms.s)
        ra_ss='00.0'
        
    dms_d=ss[s_idx+2:d_idx]
    dms_m=ss[d_idx+1:m_idx]
    dms_s=ss[m_idx+1:m_idx+3]
    
    hhmmss=ss[0:2]+ss[3:5]+ra_ss
    return hhmmss+dms_d+dms_m+dms_s

################################################################################
def cut_radius(rcutstar, mstar, mag):
    return rcutstar * 10**(0.4*2*(mstar-mag)/3.125)

def core_radius(rcorestar, mstar, mag):
    return rcorestar * 10**(0.4*(mstar-mag)/2)

def v_disp(sigstar, mstar, mag):
    return sigstar * 10**(0.4*(mstar-mag)/3.57)

###############################################################################
#####################################################################################
#### distortion function 

def distortField(ra, dec, ch):
  ###% ra ,dec are the idea position in arcsec , and the center position is ra=0, dec=0

    '''MCI geometry distortion effect in channel G R and I'''
    distortA=dict()
    
    distortB=dict()      
    
    distortA['g']=np.array([0.000586224851462036,0.999778507355332,-0.000377391397200834,-4.40403573019393e-06,0.000147444528630966,-1.93281825050268e-06,5.73520238714447e-07,-1.05194725549731e-08,7.55124671251098e-07,-3.85623216175251e-09,-1.36254798567168e-11,1.36983654024573e-09,-4.53104347730230e-11,1.50641840169747e-09,-1.05226890727366e-11,1.06471556810228e-12,-1.42299485385165e-16,5.90008722878028e-12,2.36749977009538e-17,4.11061448730916e-12,3.39287277469805e-17])
    distortB['g']=np.array([1.00113878750680,-0.000495107770623867,0.999162436209357,9.67138672907941e-05,-3.88808001621361e-06,0.000267940195644359,-3.03504629416955e-09,7.59508802931395e-07,-1.13952684052500e-08,9.58672893217047e-07,2.96397490595616e-10,-3.01506961685930e-11,2.22570315950441e-09,-4.26077028191289e-11,2.07336026666512e-09,1.71450079597168e-16,2.94455108664049e-12,8.18246797239659e-17,8.32235684990184e-12,1.05203957047210e-17,5.58541337682320e-12])
    
    distortA['r']=np.array([0.000530712822898294,0.999814397089242,-0.000366783811357609,-1.08650906916023e-05,0.000146500108063480,-4.48994291741769e-06,5.66992328511032e-07,-2.10826363791224e-08,7.51050898843171e-07,-7.74459106063614e-09,-5.78712436084704e-11,1.36389366137393e-09,-9.50057131576629e-11,1.49666815592076e-09,-2.16074939825761e-11,2.07124539578673e-12,-1.32787329448528e-13,5.78542850850562e-12,-5.82328830750271e-14,4.04881815088026e-12,2.46095156355282e-15])
    distortB['r']=np.array([1.00110972320988,-0.000483253233767592,0.999181377618578,9.60316928622784e-05,-9.06574382424422e-06,0.000266147076486786,-6.64513480754565e-09,7.53794491639207e-07,-2.25340150955077e-08,9.53032549996219e-07,2.93844965469408e-10,-6.80521155195455e-11,2.21272371816576e-09,-9.13097728847483e-11,2.06225316601308e-09,9.41096324034730e-14,3.00253083795091e-12,-1.40935245750903e-13,8.27122551593984e-12,-2.41500808188601e-13,5.52074803508046e-12])
    
    distortA['i']=np.array([0.000532645905256854,0.999813271238129,-0.000366606839338804,-1.06782246147986e-05,0.000146529811926289,-4.54488847969004e-06,5.68028930846942e-07,-2.11055358580630e-08,7.51187628288423e-07,-7.93885390157909e-09,-5.63104333653064e-11,1.36559362569725e-09,-9.64353839395137e-11,1.50231201562648e-09,-2.04159956020294e-11,1.91408535007684e-12,-7.46369323634455e-14,5.86071470639700e-12,-1.65328163262914e-13,4.07648238374705e-12,3.40880674871652e-14])
    distortB['i']=np.array([1.00111276400037,-0.000484310769918440,0.999180286636823,9.61240720951134e-05,-8.76526511019577e-06,0.000266192433565878,-6.59462433108999e-09,7.54391611102465e-07,-2.30957980169170e-08,9.53612393886641e-07,2.95558631849358e-10,-7.08307092409029e-11,2.21952307845185e-09,-9.03816603147003e-11,2.06450141393973e-09,-3.77836686311826e-14,3.13358046393416e-12,-5.20417845240954e-14,8.24487152802918e-12,-1.17846469609206e-13,5.35830020655191e-12])
    
    x=ra/0.05*10/1000  # convert ra in arcsec to mm
    
    y=dec/0.05*10/1000 # convert ra in arcsec to mm
    
    dd=np.array([1, x, y, x*x, x*y, y*y, x**3, x**2*y, x*y**2, y**3, x**4, x**3*y, x**2*y**2, x*y**3, y**4,  x**5, x**4*y, x**3*y**2, x**2*y**3 , x*y**4, y**5])
    
    xc= np.dot(distortA[ch],dd)  
    yc= np.dot(distortB[ch],dd)  
    
    distortRa =xc*1000/10*0.05-0.00265
    
    distortDec=yc*1000/10*0.05-5.0055
    
    return distortRa, distortDec

#####################################################################################
def world_to_pixel(sra,
                   sdec,
                   rotsky,
                   tra,
                   tdec,
                   x_center_pixel,
                   y_center_pixel,
                   pixelsize=0.05):
    """_summary_

    Parameters
    ----------
    sra : np.array
         star ra such as np.array([22,23,24]),unit:deg
    sdec : np.array
         stat dec such as np.array([10,20,30]),unit:deg
    rotsky : float
        rotation angel of the telescope,unit:deg
    tra : float
        telescope pointing ra,such as 222.2 deg
    rdec : float
        telescope pointing dec,such as 33.3 deg
    x_center_pixel : float
        x refer point of MCI,usually the center of the CCD,which start from (0,0)
    y_center_pixel : float
        y refer point of MCI,usually the center of the CCD,which start from(0,0)
    pixelsize : float
        pixelsize for MCI ccd, default :0.05 arcsec / pixel

    Returns
    -------
    pixel_position:list with two np.array
        such as [array([124.16937605,  99.        ,  99.        ]),
                array([ 97.52378661,  99.        , 100.50014483])]  
    """
    theta_r = rotsky / 180 * np.pi
    w = WCS(naxis=2)
    w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1]  # pixel from (1, 1)
    w.wcs.cd = np.array([[
        -np.cos(-theta_r) * pixelsize / 3600,
        np.sin(-theta_r) * pixelsize / 3600
    ],
                         [
                             np.sin(-theta_r) * pixelsize / 3600,
                             np.cos(-theta_r) * pixelsize / 3600
                         ]])
    w.wcs.crval = [tra, tdec]
    w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
    pixel_position = w.all_world2pix(sra, sdec, 0)
    return pixel_position

###############################################################################

def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec):
# ''' center_ra: telescope pointing in ra direction, unit: degree
#     center_dec: telescope pointing in dec direction, unit: degree    
#     rotTel:  telescope totation in degree    
#     rotSky:  telescope rotation in degree    
#     sRa:  source ra in arcsec
#     sDec: source dec in arcsec

    boresight = galsim.CelestialCoord(ra=-center_ra*galsim.degrees, dec=center_dec*galsim.degrees)

    q = rotTel - rotSky    
    cq, sq = np.cos(q), np.sin(q)
    jac    = [cq, -sq, sq, cq]
    affine = galsim.AffineTransform(*jac) 
    wcs    = galsim.TanWCS(affine, boresight, units=galsim.arcsec)    
    objCoord  = galsim.CelestialCoord(ra =-sRa*galsim.arcsec, dec=sDec*galsim.arcsec)
    field_pos = wcs.toImage(objCoord)    
    return -field_pos.x, field_pos.y  ## return the field position

#####################################################################################

def krebin(a, sample):
    """ Fast Rebinning with flux conservation
    New shape must be an integer divisor of the current shape.
    This algorithm is much faster than rebin_array
    Parameters
    ----------
    a : array_like
        input array
    sample :   rebinned sample 2 or 3 or 4 or other int value
    """
    # Klaus P's fastrebin from web
    size=a.shape
    shape=np.zeros(2,dtype=int)
    shape[0]=int(size[0]/sample)
    shape[1]=int(size[1]/sample)
    
    sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1]
    
    return a.reshape(sh).sum(-1).sum(1)

#####################################################################


def float2char(a, z):
    ### a is input float value
    ### transfer float a to chars and save z bis
    b=str(a)
    n=len(b)
  
    if z>=n:
        return b
    else:
        return b[:z]
#########################################################

def centroid(data):
    '''
    calculate the centroid of the input two-dimentional image data

    Parameters
    ----------
    data : input image.

    Returns
    -------
    cx: the centroid column number,  in horizontal direction definet in python image show 
    cy: the centroid row number  ,   in vertical direction 

    '''
    ###
    from scipy import ndimage
    cy,cx=ndimage.center_of_mass(data)
    return cx, cy

############################################################################### 

def fftrange(n, dtype=float):
    """FFT-aligned coordinate grid for n samples."""
    return np.arange(-n//2, -n//2+n, dtype=dtype)      
  
###############################################################################    
       
       
def dft2(ary, Q, samples, shift=None):
    """Compute the two dimensional Discrete Fourier Transform of a matrix.
    
    Parameters
    ----------
    ary : `numpy.ndarray`
        an array, 2D, real or complex.  Not fftshifted.
    Q : `float`
        oversampling / padding factor to mimic an FFT.  If Q=2, Nyquist sampled
    samples : `int` or `Iterable`
        number of samples in the output plane.
        If an int, used for both dimensions.  If an iterable, used for each dim
    shift : `float`, optional
        shift of the output domain, as a frequency.  Same broadcast
        rules apply as with samples.
    
    Returns
    -------
    `numpy.ndarray`
        2D array containing the shifted transform.
        Equivalent to ifftshift(fft2(fftshift(ary))) modulo output
        sampling/grid differences
    
    """
    
    # this is for dtype stabilization
    Q = float(Q)  ##  Q=lambda*f/D/pixelsize 

    n, m = ary.shape ###ary maybe is the pupil function
    N, M = samples,samples
    
    X, Y, U, V = (fftrange(n) for n in (m, n, M, N))

    a = 1 / Q
    
    ###################################################################           
    Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T)
    Ein_fwd  = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U))

    #############################################################################    
    
    out = Eout_fwd @ ary @ Ein_fwd   #### ary is the input pupil function 
    
    return out
##############################################################################
    
def idft2(ary, Q, samples, shift=None):
    """Compute the two dimensional inverse Discrete Fourier Transform of a matrix.
    
    Parameters
    ----------
    ary : `numpy.ndarray`
        an array, 2D, real or complex.  Not fftshifted.
    Q : `float`
        oversampling / padding factor to mimic an FFT.  If Q=2, Nyquist sampled
    samples : `int` or `Iterable`
        number of samples in the output plane.
        If an int, used for both dimensions.  If an iterable, used for each dim
    shift : `float`, optional
        shift of the output domain, as a frequency.  Same broadcast
        rules apply as with samples.
    
    Returns
    -------
    `numpy.ndarray`
        2D array containing the shifted transform.
        Equivalent to ifftshift(ifft2(fftshift(ary))) modulo output
        sampling/grid differences
    
    """
    # this is for dtype stabilization
    Q = float(Q)  ##  Q=lambda*f/D/pixelsize 

    n, m = ary.shape ###ary maybe is the pupil function
    N, M = samples
    
    X, Y, U, V = (fftrange(n) for n in (m, n, M, N))

    ###############################################################################
    nm = n*m
    NM = N*M
    r = NM/nm
    a = 1 / Q
    
    ###################################################################           
    # Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T)
    # Ein_fwd  = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U))
    Eout_rev = np.exp(1j * 2 * np.pi * a / n * np.outer(Y, V).T) * (1/r)
    Ein_rev  = np.exp(1j * 2 * np.pi * a / m * np.outer(X, U)) * (1/nm)

    ########################################################################### 
    
    out = Eout_rev @ ary @ Ein_rev
    
    return out

###############################################################################

###############################################################################
def opd2psf(opd,cwave,oversampling):
    '''
    calculate  psf from opd data with defined wavelength

    Parameters
    ----------
    opd : TYPE
        the input opd data.
    cwave : TYPE
        the defined wavelength.

    '''

    D=2;
    focLength=41.253              # % MCI focal length in meters
    pixel=10/oversampling         # % pixel size in microns 
    
    opd=opd*1e9; ## convert opd unit of meter to nm

    zoomopd=ndimage.zoom(opd, 2, order=1)

    m,n=np.shape(zoomopd);
    clambda=cwave;
 
    wfe=zoomopd/clambda;   #%% the main telescope aberration ,in wavelength;   
        
    Q=clambda*1e-3*focLength/D/pixel        
    
    pupil=np.zeros((m,m));

    pupil[abs(zoomopd)>0]=1;  
    #idx=pupil>0;
     
    Nt=512  
    
    phase=2*np.pi*wfe; #% wavefront phase in radians     
           
    #%generalized pupil function of channel1;        
    pk=pupil*np.exp(cmath.sqrt(-1)*(phase))   
    
    pf=dft2(pk, Q, Nt)
    
    pf2 = abs(pf)**2
    psfout=pf2/pf2.sum()
    
    cx,cy=centroid(psfout)    
    
    psf=ndimage.shift(psfout,[Nt/2-1-cy, Nt/2-1-cx],order=1, mode='nearest' )  
    
    return psf
########################################################################

def  cal_PSF_new(channel,wfetimes,oversampling):
    # filed postion: Ra=ys1,  Dec=ys2
    
    ############ use opd cal PSF on the defined field;
    #% psf interpolation test 20210207
    # xfmin=-0.07 # % the filed minmum  value in the x direction in degrees    
    # xfmax= 0.07  #% the filed maxmum  value in the x direction in degrees    
    # yfmin=-0.35  #% the filed minmum  value in the y direction in degrees    
    # yfmax=-0.49 # % the filed maxmum  value in the y direction in degrees
    
    cfx=  0.2*3600;  #   % field center in x derection in arcsec;    
    cfy= -0.45*3600; #   % field center in y derection in arcsec;
    
    #wavelength=[255,337,419,501,583,665,747,829,911,1000]
    
    cwave=dict()
    
    if channel == 'g':
        cwave[channel]=475
        waven=4
    if channel == 'r':
        cwave[channel]=625
        waven=6
    if channel == 'i':
        cwave[channel]=776
        waven=7
        
    n=np.arange(1,101,1)
    
    ffx= 0.1+ 0.02222222*((n-1)%10)
    ffy=-0.35-0.02222222*np.floor((n-0.1)/10)

    psf=dict()
    fx=dict()
    fy=dict()
    fx[channel]=ffx*3600-cfx   
     
    fy[channel]=ffy*3600-cfy
    
    Nt=512