genSurveyDFSByPointing.py

'''
Author: Zhang Xin zhangx@bao.ac.cn
Date: 2025-04-21 09:53:28
LastEditors: Zhang Xin zhangx@bao.ac.cn
LastEditTime: 2025-05-20 13:06:28
FilePath: /CSST_Survey/Users/zhangxin/Work/SurveyPlan/point/pointing_c9/genSurveyDFSByPointing.py
Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
'''
import json
from datetime import datetime
from astropy.table import Table, join  # type: ignore
import numpy as np  # type: ignore
from astropy.time import Time
from datetime import datetime, timedelta
import pytz
import sys


def calculateTransTime(transAngle=1.0, shuterTime=1.5):
    angleVStime = np.array([[1, 20, 45, 180], [80, 127, 196, 581]])

    if transAngle < angleVStime[0, 0]:
        tTime = 70
    elif transAngle == angleVStime[0, 0]:
        tTime = angleVStime[1, 0]
    else:
        for i in np.arange(1, 4, 1):
            if (transAngle > angleVStime[0, i-1] and transAngle <= angleVStime[0, i]):
                tTime = angleVStime[1, i-1] * ((transAngle - angleVStime[0, i])) / (((angleVStime[0, i-1] - angleVStime[0, i]))) + \
                    angleVStime[1, i] * ((transAngle - angleVStime[0, i-1])) / \
                    (((angleVStime[0, i] - angleVStime[0, i-1])))
                break

    return tTime + shuterTime * 2.0


def JDTimeToBeijingTime(jd_time=2459766.):
    # 使用 astropy 转换到UTC时间
    t = Time(jd_time, format='jd', scale='utc')
    utc_time = t.datetime  # 转换为 datetime 对象
    # 转换为北京时间（UTC+8）
    beijing_time = utc_time + timedelta(hours=8)
    # 使用 pytz 添加时区信息（可选）
    beijing_tz = pytz.timezone('Asia/Shanghai')
    beijing_time = beijing_tz.localize(beijing_time)
    # 格式化输出
    formatted_time = beijing_time.strftime("%Y-%m-%d %H:%M:%S")
    return formatted_time, beijing_time


def dataInit():
    # 定义数据字典
    data = {
        "id": 1593,
        "opid": "101000001593",
        "obs_id": "10100001593",
        "instrument": "MSC",
        "project_id": "None",
        "file_type": "SCI",
        "obs_type": "WIDE",
        "object": "m31",
        "mode": "OM1",
        "starttime": "2031-04-23 16:40:19",
        "endtime": "2031-04-23 16:43:37",
        "duration": 150.0,
        "params": {
            "target_ra": 146.4567,
            "target_dec": 54.0354,
            "T1": 150.0,
            "t0": 190.0,
            "T2": 22.0,
            "T4": 20,
            "exposure_start_time": "2031-04-23 16:41:04",
            "starttime": "2031-04-23 16:40:19",
            "exposure_duration": 150.0,
            "endtime": "2031-04-23 16:43:37"
        },
        "simresult": {
            "jd": 2462979.8618552,
            "decH": 54.0354,
            "raH": 146.4567,
            "sat_x": -4006.7755,
            "sat_y": 5432.4733,
            "sat_z": 321.34,
            "sun_x": 126733141.1541,
            "sun_y": 74307866.5619,
            "sun_z": 32209476.4337,
            "moon_x": 245820.5699,
            "moon_y": 284285.6117,
            "moon_z": 127374.2175,
            "deepArea": 0.0,
            "veryDeepArea": 0.0,
            "ca": -1.0,
            "cb": "0",
            "exposure": 150.0,
            "turnAngle": 0.95482,
            "isInSun": "1",
            "cmg": 0.6479,
            "charge": 97200.0,
            "betweenAngle": 19.1178,
            "ssaSec": 0.0,
            "skyid": "523664",
            "boardAngle": -0.6703,
            "quad_0": 1.0,
            "quad_i": 0.004,
            "quad_j": 0.0067,
            "quad_k": 0.003,
            "sky_type": "1",
            "id": 1220,
            "decCenter": 60.6576587868062,
            "raCenter": 182.834345163608,
            "decHigh": 60.7576587868062,
            "decLow": 60.5576587868062,
            "topLeftRa": 182.619403486662,
            "topRightRa": 183.049286840554,
            "bottomLeftRa": 182.620734042365,
            "bottomRightRa": 183.047956284851
        }
    }
    return data


if __name__ == "__main__":
    if len(sys.argv) < 3:
        print('input_file_name  output_file_name')

    fn1 = 'E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.obsid.csv'
    # fn2 = 'pointing_50_5_n.dat'
    fn2 = sys.argv[1]
    outFn = sys.argv[2]
    t1 = Table.read(fn1, format='ascii')
    t2 = Table.read(fn2, format='ascii')
    t1['obs_id'] = (t1['obs_id']//1e11*1e8+t1['obs_id'] % 1e8).astype(np.int64)

    t2.rename_column('id', 'obs_id')
    inputData = join(t1, t2, keys='obs_id')

    # dfn = 'pointing_50_5_n_obsid.dat'
    # inputData = Table.read(dfn, format='ascii')

    all_data = []
    for i in np.arange(len(inputData)):
        d1 = inputData[i]
        data = dataInit()
        data['id'] = int(d1['id'])
        data['obs_id'] = d1['obs_id']
        data['opid'] = (str(d1['obs_id'])[0:3]+'0'+str(d1['obs_id'])[3:])
        data['instrument'] = 'MSC'
        data['file_type'] = 'SCI'
        if d1['sky_type'] == 1:
            data['obs_type'] = 'WIDE'
        else:
            data['deep'] = 'WIDE'
        data['object'] = d1['skyid']

        data['starttime'] = JDTimeToBeijingTime(
            d1['jd'])[1].isoformat()
        data['endtime'] = JDTimeToBeijingTime(
            d1['jd']+d1['exposure']/86400.)[1].isoformat()
        data['duration'] = d1['exposure']

        # param
        data['params']['target_ra'] = d1['ra_2']
        data['params']['target_dec'] = d1['dec_2']
        data['params']['T1'] = d1['exposure']
        data['params']['t0'] = data['params']['T1']+40
        data['params']['T2'] = calculateTransTime(
            d1['turnAngle'])-20-1.5*2  # 快门打开和关闭的时间，分别是1.5s
        data['params']['T4'] = 20
        # exposure_time1, exposure_time = JDTimeToBeijingTime(d1['jd'])
        data['params']['exposure_start_time'] = JDTimeToBeijingTime(
            d1['jd'])[1].isoformat()
        data['params']['starttime'] = JDTimeToBeijingTime(
            d1['jd']-(1+1+1 + data['params']['T2'] + data['params']['T4'])/86400.)[1].isoformat()
        # exposure_time - \
        #     timedelta(seconds=1+1+1 + data['params']['T2'] + data['params']['T4'])
        data['params']['exposure_duration'] = d1['exposure']
        data['params']['endtime'] = JDTimeToBeijingTime(
            d1['jd']+(d1['exposure']+3)/86400.)[1].isoformat()  # 曝光完了，3s关快门，总控配电1s不计

        # simresult
        data['simresult']['jd'] = d1['jd']
        data['simresult']['decH'] = d1['decH']
        data['simresult']['raH'] = d1['raH']
        data['simresult']['sat_x'] = d1['sat_x_1']
        data['simresult']['sat_y'] = d1['sat_y_1']
        data['simresult']['sat_z'] = d1['sat_z_1']
        data['simresult']['sun_x'] = d1['sun_x_1']
        data['simresult']['sun_y'] = d1['sun_y_1']
        data['simresult']['sun_z'] = d1['sun_z_1']
        data['simresult']['moon_x'] = d1['moon_x_1']
        data['simresult']['moon_y'] = d1['moon_y_1']
        data['simresult']['moon_z'] = d1['moon_z_1']
        data['simresult']['deepArea'] = d1['deepArea']
        data['simresult']['veryDeepArea'] = d1['veryDeepArea']
        data['simresult']['ca'] = d1['ca']
        data['simresult']['cb'] = d1['cb']
        data['simresult']['exposure'] = d1['exposure']
        data['simresult']['turnAngle'] = d1['turnAngle']
        data['simresult']['isInSun'] = d1['isInSun']
        data['simresult']['cmg'] = d1['cmg']
        data['simresult']['charge'] = d1['charge']
        data['simresult']['skyid'] = d1['skyid']
        data['simresult']['boardAngle'] = d1['boardAngle']
        data['simresult']['quad_0'] = d1['quad_0']
        data['simresult']['quad_i'] = d1['quad_i']
        data['simresult']['quad_j'] = d1['quad_j']
        data['simresult']['quad_k'] = d1['quad_k']
        data['simresult']['sky_type'] = d1['sky_type']
        data['simresult']['id'] = d1['id']
        data['simresult']['decCenter'] = d1['dec_1']
        data['simresult']['raCenter'] = d1['ra_1']
        data['simresult']['topLeftRa'] = 0.
        data['simresult']['topRightRa'] = 0.
        data['simresult']['bottomLeftRa'] = 0.
        data['simresult']['bottomRightRa'] = 0.

        all_data.append(data)
        # break

    json_data = json.dumps(all_data, indent=4, default=str)

    # 打印JSON格式的字符串
    # print(json_data)

    with open(outFn, "w") as file:
        file.write(json_data)