Merge branch 'master' into 'fandongwei'

"""

Author: Zhang Xin zhangx@bao.ac.cn

Date: 2024-07-02 13:55:19

LastEditors: Zhang Xin zhangx@bao.ac.cn

LastEditTime: 2024-09-26 15:43:32

FilePath: /CSST_Survey/setup.py

Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE

"""

from setuptools import setup, find_packages

    license="NAOC",

    packages=find_packages(),

    # python_requires=">=3.9",

    install_requires=[

        "numpy",

        "jax"

    ],

    install_requires=["numpy", "jax"],

    # requires=['numpy', 'scipy', 'astropy', 'drizzlepac', 'stwcs'],

    # long_description=read('README.rst'),

    # classifiers=[

    #     'Topic :: Scientific/Engineering :: Astronomy',

    # ],

    package_data={

        # "survey_sim.ephemeris.data": ["*.dylib", "*.so", "jpl.405"],

        # "survey_sim.satOrbit.data": ["*"],

        "survey_sim.strategies.data": ["*"],

        "survey_sim.ephemeris": ["data/**/*"],

        "survey_sim.satOrbit": ["data/**/*"],

    },

from astropy import units as u

from astropy.time import Time

import math

import numpy as np

"""

description: 

param {*} t: time, jd data

return {*} [lonitude, latitude]

"""

PI_180 = math.pi / 180

M_1_PI = 1 / math.pi

def getSatSubpoint(satPos=None, t=2459766.0):

    longitude = earth_location.lon.deg

    return [longitude, latitude]

def calculateAngle(ra1, dec1, ra2, dec2):

    # double x1, y1, z1, x2, y2, z2, angle, cosValue;

    x1 = math.cos(dec1 * PI_180) * math.cos(ra1 * PI_180)

    y1 = math.cos(dec1 * PI_180) * math.sin(ra1 * PI_180)

    z1 = math.sin(dec1 * PI_180)

    x2 = math.cos(dec2 * PI_180) * math.cos(ra2 * PI_180)

    y2 = math.cos(dec2 * PI_180) * math.sin(ra2 * PI_180)

    z2 = math.sin(dec2 * PI_180)

    cosValue = x1 * x2 + y1 * y2 + z1 * z2

    if cosValue < -1.0:

        return 180.0

    elif cosValue > 1.0:

        return 0.0

    angle = math.acos(x1 * x2 + y1 * y2 + z1 * z2)

    return angle * 180 * M_1_PI

    # 返回值为角度（以度数为单位，不是以弧度为单位）

def getAngle132(x1, y1, z1, x2, y2, z2, x3, y3, z3):

    x11 = x1 - x3

    y11 = y1 - y3

    z11 = z1 - z3

    x22 = x2 - x3

    y22 = y2 - y3

    z22 = z2 - z3

    tt = np.sqrt(

        (x11 * x11 + y11 * y11 + z11 * z11) * (x22 * x22 + y22 * y22 + z22 * z22)

    )

    cosValue = (x11 * x22 + y11 * y22 + z11 * z22) / tt

    cosValue = np.clip(cosValue, -1, 1)

    angle = math.acos(cosValue)

    return angle * 180 * M_1_PI

# 根据指定的转动轴和转动角度生成旋转矩阵

def GenRotationMatrix(u=np.array([0, 0, 0]), angle_deg=10.0):

    theta = angle_deg * PI_180

    cos_theta = math.cos(theta)

    sin_theta = math.sin(theta)

    One_cos_theta = 1 - cos_theta

    u_norm = np.linalg.norm(u)

    u = u / u_norm

    ux = u[0]

    uy = u[1]

    uz = u[2]

    if np.fabs(u_norm - 1.0) > 1e-5:

        print("in GenRotationMatrix: u_norm differs too much from 1.0!")

    R = np.zeros([3, 3])

    R[0, 0] = cos_theta + ux * ux * One_cos_theta

    R[0, 1] = ux * uy * One_cos_theta - uz * sin_theta

    R[0, 2] = ux * uz * One_cos_theta + uy * sin_theta

    R[1, 0] = uy * ux * One_cos_theta + uz * sin_theta

    R[1, 1] = cos_theta + uy * uy * One_cos_theta

    R[1, 2] = uy * uz * One_cos_theta - ux * sin_theta

    R[2, 0] = uz * ux * One_cos_theta - uy * sin_theta

    R[2, 1] = uz * uy * One_cos_theta + ux * sin_theta

    R[2, 2] = cos_theta + uz * uz * One_cos_theta

    return R

def rodrigues_rotation_formula(axis, theta):

    """根据旋转轴向量和旋转角度生成旋转矩阵, 与GenRotationMatrix()结果一样"""

    theta = theta * PI_180

    axis = axis / np.linalg.norm(axis)  # 确保轴向量归一化

    cos_theta, sin_theta = np.cos(theta), np.sin(theta)

    # 构造斜对称矩阵

    S = np.array(

        [[0, -axis[2], axis[1]], [axis[2], 0, -axis[0]], [-axis[1], axis[0], 0]]

    )

    # 构造旋转矩阵

    R = cos_theta * np.eye(3) + (1 - cos_theta) * np.outer(axis, axis) + sin_theta * S

    return R

def get_rotation_axis(mat):

    """根据旋转轴旋转矩阵,生成旋转轴"""

    assert np.allclose(np.dot(mat.T, mat), np.eye(3)), "R must be orthogonal"

    assert np.isclose(np.linalg.det(mat), 1.0), "R must have a determinant of 1"

    axis = np.zeros(3)

    # axis[0] = mat[7]-mat[5];

    # axis[1] = mat[2]-mat[6];

    # axis[2] = mat[3]-mat[1];

    axis[0] = mat[2, 1] - mat[1, 2]

    axis[1] = mat[0, 2] - mat[2, 0]

    axis[2] = mat[1, 0] - mat[0, 1]

    # axis_norm = np.sqrt(axis[0]*axis[0]+axis[1]*axis[1]+axis[2]*axis[2]);

    axis_norm = np.linalg.norm(axis)

    axis[0] /= axis_norm

    axis[1] /= axis_norm

    axis[2] /= axis_norm

    return axis

def get_RotationAngleFromMatrix(mat):

    assert np.allclose(np.dot(mat.T, mat), np.eye(3)), "R must be orthogonal"

    assert np.isclose(np.linalg.det(mat), 1.0), "R must have a determinant of 1"

    trace = np.trace(mat)

    cos_theta = np.clip((trace - 1) / 2, -1, 1)

    angle = 180.0 / math.pi * math.acos(cos_theta)

    return angle

# def getRotationAngleForAttitude()

# //  ==================================================================

# //  原先所采取的旋转方式在某些情况下给出的旋转角度不是最合理的。例如在靠近黄极的时候，

# //  如果两个指向的经度相差接近于180度，利用原先的旋转方式得到旋转角度就接近于180度，

# //  但实际上可以通过更小角度的转动来完成指向的移动。

# //

# //  注意：新版本返回的角度单位是“度数”,不再使用“弧度”！

# //

# //  另外，目前还没有优化这个函数。

def Get_RotationAngle(

    ra_old=60.0,

    dec_old=-40.0,

    ra_new=65.0,

    dec_new=-43.0,

):

    angle_rot = 0

    rot_axis = np.zeros(3)

    dec_old = 90 - dec_old

    dec_new = 90 - dec_new

    sin_ra_old = math.sin(ra_old * PI_180)

    cos_ra_old = math.cos(ra_old * PI_180)

    sin_dec_old = math.sin(dec_old * PI_180)

    cos_dec_old = math.cos(dec_old * PI_180)

    sin_ra_new = math.sin(ra_new * PI_180)

    cos_ra_new = math.cos(ra_new * PI_180)

    sin_dec_new = math.sin(dec_new * PI_180)

    cos_dec_new = math.cos(dec_new * PI_180)

    #   2021-07-03:这里采用的是数学上标准的球坐标与直角坐标之间的转换，与天文上的定义不一样，所以在使用该

    #              函数时，需要对dec做一个变换: "dec" --> "90-dec"

    # TODO 为什么不用天文的坐标系？？？？？？？？？？？？

    p_old = np.array([sin_dec_old * cos_ra_old, sin_dec_old * sin_ra_old, cos_dec_old])

    p_new = np.array([sin_dec_new * cos_ra_new, sin_dec_new * sin_ra_new, cos_dec_new])

    p_tmp = np.ones(3) * -999

    # ====================================================================

    # step 1:绕z轴旋转，将旧指向旋转到新指向所在的子午圈，并计算出旋转后的“临时”指向矢量

    delta_alpha = 0.0

    delta_alpha1 = 0.0

    delta_alpha2 = 0.0

    if (

        np.fabs(ra_new - ra_old) < 1e-5

    ):  # 这个情况下可以近似认为两个指向处在同一个子午圈内

        if np.fabs(dec_new - dec_old) < 1e-5:

            # // printf("==> getting rotation axis : 1\n");

            rot_axis[0] = 0

            rot_axis[1] = 0

            rot_axis[2] = 1

            return (

                angle_rot,

                rot_axis,

                0,

            )  # 两个指向相同的情况下不需要进行任何转动,可直接返回

        delta_alpha = 0

        delta_alpha1 = 0  # 等价于不做任何转动。

        delta_alpha2 = 180  # 这种情况应该予以排除，因为会导致帆板面的指向反转，造成无法接受太阳光照进行发电。

        p_tmp[0] = p_old[0]

        p_tmp[1] = p_old[1]

        p_tmp[2] = p_old[2]

        alpha_n = np.array([0, 0, 1])

        Rz = GenRotationMatrix(alpha_n, 0)

    else:

        # alpha_old,alpha_new是在XY平面内的单位向量

        alpha_old = np.array([cos_ra_old, sin_ra_old, 0])

        alpha_new = [cos_ra_new, sin_ra_new, 0]

        # 通过ＸＹ平面内的两个（单位长度）指向的叉乘获取旋转轴以及相应的旋转方向(用于调整黄经)

        alpha_n = np.cross(alpha_old, alpha_new)

        alpha_n_norm = np.linalg.norm(alpha_n)

        if np.fabs(alpha_n_norm) < 1e-6:  # 原先是1e-8,导致在台式机上可能出现NaN问题

            # 	说明 alpha_old[] 与 alpha_new[] 指向相反的方向，就直接绕z轴旋转180度即可????? 应该不需要转，要不帆板就转了

            alpha_n[0] = 0

            alpha_n[1] = 0

            alpha_n[2] = 1

            delta_alpha = 0.0

            # delta_alpha = 180.

            Rz = GenRotationMatrix(alpha_n, delta_alpha)

        else:

            alpha_n[0] /= alpha_n_norm

            alpha_n[1] /= alpha_n_norm

            alpha_n[2] /= alpha_n_norm

            # 根据alpha_old,alpha_new来确定如何绕z-轴旋转

            # TODO cosval 难道不是经度的夹角？？？？？？？？

            cosval = np.dot(alpha_old, alpha_new)

            # csst_test( fabs(cosval) < -1.0 || fabs(cosval) > 1.0, errmsg, "cosval is out of range [-1,1]!" );

            # 确保不会出现数值计算错误（之前在台式机上运行仿真时，总在这个模块内出现问题）

            delta_alpha1 = math.acos(cosval) * 180 / math.pi

            # 这个值始终是大于0的

            delta_alpha2 = 180 - delta_alpha1  # 这个值始终也是大于0的

            #  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            #  因为旋转角超过90度时帆板面的指向会“反转”，导致不能接受太阳光照，因此只考虑

            #  旋转角小于90度的情况。

            #  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            if np.fabs(delta_alpha1) <= 90:

                delta_alpha = delta_alpha1

                # 转动的方向由转动轴的指向决定，因此传入的角度值都是正数

                Rz = GenRotationMatrix(alpha_n, delta_alpha1)

            elif np.fabs(delta_alpha2) <= 90:

                delta_alpha = delta_alpha2

                # 转动的方向由转动轴的指向决定，因此传入的角度值都是正数

                # 这里需要注意的是需要将alpha_n反号

                alpha_n[0] *= -1

                alpha_n[1] *= -1

                alpha_n[2] *= -1

                Rz = GenRotationMatrix(alpha_n, delta_alpha2)

            else:

                return -999, -999, -999

        # MatrixVecProduct(Rz,p_old,p_tmp);

        p_tmp = np.dot(Rz, p_old)

    #  ====================================================================

    #  step 2:根据 p_tmp 和 p_new 的叉乘来计算第二次旋转的旋转轴；虽然确定了旋转轴，

    #  但是依旧存在两种旋转方式。

    axis = np.zeros(3)

    # 首先还是先判断“临时指向”与目标指向是否共线(包括同向与反向)

    p_tmp_dot_p_new = np.dot(p_tmp, p_new)

    if np.fabs(p_tmp_dot_p_new - 1) < 1e-8:  # p_tmp与p_new重合

        angle_rot = delta_alpha

        rot_axis = get_rotation_axis(Rz)

        return angle_rot, rot_axis, Rz

    if np.fabs(p_tmp_dot_p_new + 1) < 1e-8:  # p_tmp与p_new反向

        #  旋转轴位于XY平面内，因此可以直接根据ra_new+pi/2直接计算出来旋转轴；此时

        #  旋转轴的指向不影响结果

        axis[0] = math.cos(0.5 * math.pi + ra_new * PI_180)

        axis[1] = math.sin(0.5 * math.pi + ra_new * PI_180)

        axis[2] = 0

        Rn = GenRotationMatrix(axis, 180)

        R = np.dot(Rn, Rz)  # 此处需要注意旋转矩阵乘积的顺序

        # MatrixMultiplication(Rn,Rz);  //  此处需要注意旋转矩阵乘积的顺序

        angle_rot = get_RotationAngleFromMatrix(R)

        # printf("==> getting rotation axis : 3\n");

        get_rotation_axis(R, rot_axis)

        return angle_rot, rot_axis, R

    #  给存贮绕n-轴的两个旋转矩阵分配内存（n轴由两个矢量的叉积得到）

    cosval = np.clip(np.dot(p_tmp, p_new), -1, 1)

    delta_alpha1 = math.acos(cosval) * 180 / math.pi  # 这个值始终是大于0的

    # delta_alpha2 = 180-delta_alpha1; //这个值始终也是大于0的，但表示反方向旋转

    axis = np.cross(p_tmp, p_new)

    axis_norm = np.linalg.norm(axis)

    axis[0] /= axis_norm

    axis[1] /= axis_norm

    axis[2] /= axis_norm

    R1 = GenRotationMatrix(axis, delta_alpha1)

    R = np.dot(R1, Rz)

    angle_rot = get_RotationAngleFromMatrix(R)

    rot_axis = get_rotation_axis(R)

    return angle_rot, rot_axis, R

# 计算下来貌似没有考虑到姿态,实际上计算的就是p1和p2的夹角？？？？？？？？？？？

def Get_RotationAngle_old(ra_old, dec_old, ra_new, dec_new):

    dec_old = 90 - dec_old

    dec_new = 90 - dec_new

    p1 = np.zeros(3)

    p2 = np.zeros(3)

    p1[0] = math.sin(dec_old * PI_180) * math.cos(ra_old * PI_180)

    p1[1] = math.sin(dec_old * PI_180) * math.sin(ra_old * PI_180)

    p1[2] = math.cos(dec_old * PI_180)

    p2[0] = math.sin(dec_new * PI_180) * math.cos(ra_new * PI_180)

    p2[1] = math.sin(dec_new * PI_180) * math.sin(ra_new * PI_180)

    p2[2] = math.cos(dec_new * PI_180)

    ppNorm = np.cross(p1, p2)

    pn1 = np.cross(p1, ppNorm)

    pn2 = np.cross(p2, ppNorm)

    angle = getAngle132(pn1[0], pn1[1], pn1[2], pn2[0], pn2[1], pn2[2], 0, 0, 0)

    R = GenRotationMatrix(u=ppNorm, angle_deg=angle)

    return angle, ppNorm, R

from survey_sim.satOrbit import locateSat

import numpy as np

import math

from survey_sim.ephemeris import locate_sun

import matplotlib.pyplot as plt

# def get_betaAngle(time = 2459799, orbitData = None):

class beta_time_constraint(object):

    def __init__(

        self,

        startTime=2459766.0,

        endTime=2463416.0,

        orbitData=None,

        ephLib=None,

        surveyCosntraint=None,

    ):

        self.startTime = startTime

        self.endTime = endTime

        self.orbitData = orbitData

        self.ephLib = ephLib

        self.surveyCosntraint = surveyCosntraint

        self.beta_time_seg = self.get_beta_time()

    def get_beta_time(self):

        if self.orbitData is None:

            return

        orbDataLen = len(self.orbitData)

        beta_time_i = 0

        beta_time_tot = 0

        i_start = 0

        i_end = 0

        if self.startTime < self.orbitData[0, 0]:

            print("ERROR: start Time is not in the range of orbit data!!!!!!!!!!!!")

            return

        if self.endTime > self.orbitData[orbDataLen - 1, 0]:

            print("ERROR: end Time is not in the range of orbit data!!!!!!!!!!!!")

            return

        for i in np.arange(0, orbDataLen, 1):

            t = self.orbitData[i, 0]

            if self.startTime < t:

                i_start = i

                break

        for i in np.arange(i_start, orbDataLen, 1):

            t = self.orbitData[i, 0]

            if self.endTime < t:

                i_end = i - 1

                break

        if i_end < i_start:

            print(

                "ERROR: the time between start and end is tooooo short, it must be large than the data interval (about 120s)"

            )

            return

        tmp_orbit_len = i_end + 1 - i_start + 2

        tmp_orbitData = np.zeros([tmp_orbit_len, 4])

        sat_start, _, nid = locateSat(

            time=self.startTime, OrbitData=self.orbitData, orbDataLen=orbDataLen

        )

        sat_end, _, _ = locateSat(

            time=self.endTime,

            OrbitData=self.orbitData,

            startId=nid,

            orbDataLen=orbDataLen,

        )

        tmp_orbitData[0, :] = np.array(

            [self.startTime, sat_start[0], sat_start[1], sat_start[2]]

        )

        tmp_orbitData[-1, :] = np.array(

            [self.endTime, sat_end[0], sat_end[1], sat_end[2]]

        )

        tmp_orbitData[1:-1] = self.orbitData[i_start : i_end + 1, 0:4]

        orbitSegStartIds = [0]

        for i in np.arange(1, tmp_orbit_len, 1):

            t1 = tmp_orbitData[i, 0]

            t2 = tmp_orbitData[i - 1, 0]

            if t1 - t2 > 0.5:

                orbitSegStartIds.append(i)

        beta_time_seg = []

        segNum = len(orbitSegStartIds)

        for k in np.arange(len(orbitSegStartIds)):

            seg_start_i = orbitSegStartIds[k]

            if k + 1 < segNum:

                seg_end_i = orbitSegStartIds[k + 1]

            else:

                seg_end_i = tmp_orbit_len

            in_beta = 0

            time_seg_start = 0

            time_seg_end = 0

            for i in np.arange(seg_start_i, seg_end_i - 1, 1):

                t1 = tmp_orbitData[i, 0]

                t2 = tmp_orbitData[i + 1, 0]

                if t2 - t1 < 1.1574074074074074e-08:  # 1ms

                    continue

                curTime = 0.5 * (t1 + t2)

                # // double sunAngle = getSun2OrbitAngle1(infp, curTime,orbitData, orbitDataNum, Orbit_File_Num);

                # double a[3],b[3],normalVect[3],sun[3];

                a = np.zeros(3)

                b = np.zeros(3)

                a[0] = tmp_orbitData[i, 1]

                a[1] = tmp_orbitData[i, 2]

                a[2] = tmp_orbitData[i, 3]

                b[0] = tmp_orbitData[i + 1, 1]

                b[1] = tmp_orbitData[i + 1, 2]

                b[2] = tmp_orbitData[i + 1, 3]

                normalVect = np.cross(a, b)

                sun = locate_sun(curTime, ephLib=self.ephLib)

                pointMul = (

                    sun[0] * normalVect[0]

                    + sun[1] * normalVect[1]

                    + sun[2] * normalVect[2]

                )

                modSun = np.sqrt(sun[0] * sun[0] + sun[1] * sun[1] + sun[2] * sun[2])

                modNormal = np.sqrt(

                    normalVect[0] * normalVect[0]

                    + normalVect[1] * normalVect[1]

                    + normalVect[2] * normalVect[2]

                )

                sunAngle = math.acos(pointMul / (modSun * modNormal)) * 57.29577951

                sunAngle = 90 - sunAngle

                # print(t1, t2, (t2 - t1) * 86400, in_beta, sunAngle)

                if np.fabs(sunAngle) < self.surveyCosntraint.BETA_ANGLE:

                    if in_beta == 0:

                        time_seg_start = t1

                        in_beta = 1

                    time_seg_end = t1

                else:

                    # elif np.fabs(sunAngle) >= surveyCosntraint.BETA_ANGLE:

                    if in_beta == 1:

                        time_seg_end = t1

                        in_beta = 0

                        beta_time_seg.append([time_seg_start, time_seg_end])

            if in_beta == 1:

                beta_time_seg.append([time_seg_start, time_seg_end])

        beta_time_seg = np.array(beta_time_seg)

        beta_time = beta_time_seg[:, 1] - beta_time_seg[:, 0]

        d_ids = np.where(beta_time < 0.5)

        d_ids_flat = np.hstack((d_ids[0] * 2, d_ids[0] * 2 + 1))

        beta_time_seg_d_flat = np.delete(beta_time_seg, d_ids_flat)

        self.beta_time_seg = beta_time_seg_d_flat.reshape(

            beta_time_seg.shape[0] - d_ids[0].shape[0], 2

        )

        # return beta_time_seg_d

    def get_survey_time_segment_MSC(self):

        # if self.beta_time_seg is None:

        #     self.get_beta_time()

        MSC_time = np.hstack((np.array(self.startTime), self.beta_time_seg.flatten()))

        # beta_time_seg_flat = self.beta_time_seg.flatten()

        # MSC_time = np.stack((np.array([self.startTime]), beta_time_seg_flat))

        MSC_time = np.hstack((MSC_time, np.array(self.endTime)))

        self.MSC_time = MSC_time.reshape(self.beta_time_seg.shape[0] + 1, 2)

        # return MSC_time

    def get_beta_useTime(self):

        time_seg = self.beta_time_seg[:, 1] - self.beta_time_seg[:, 0]

        return np.sum(time_seg)

    def plotFig(self):

        plt.figure()

        h = 1

        for xx in self.MSC_time:

            plt.plot((xx - self.startTime) / 365.25, [h, h], "b")

            h = h + 1

        h = 1.5

        for xx in self.beta_time_seg:

            plt.plot((xx - self.startTime) / 365.25, [h, h], "r")

            h = h + 1

        plt.yticks([])

        plt.show()

'''

Author: Zhang Xin zhangx@bao.ac.cn

Date: 2020-06-17 17:03:15

LastEditors: Zhang Xin zhangx@bao.ac.cn

LastEditTime: 2024-06-12 10:54:54

FilePath: /survey-dev/Users/zhangxin/Work/SurveyPlan/CSST_Survey/csst_survey_sim/constraints/sun_constraint.py

Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE

'''

'''

Author: Zhang Xin zhangx@bao.ac.cn

Date: 2020-06-17 17:03:15

LastEditors: Zhang Xin zhangx@bao.ac.cn

LastEditTime: 2024-06-12 10:54:54

FilePath: /survey-dev/Users/zhangxin/Work/SurveyPlan/CSST_Survey/csst_survey_sim/constraints/sun_constraint.py

Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE

'''