Commit 627689db authored by Zhang Xin's avatar Zhang Xin
Browse files

adjust observe num; add some sky obs

parent 95624292

survey_sim/skyField/skyMap.py

+6 −1
Original line number Diff line number Diff line
@@ -210,6 +210,8 @@ class SkyMap(object): 
                    # self.skymap[key].setCoordinateValue_icrs()

                    # self.skymap[key].setCoordinateValue_cartesian()

                    self.skymap[key].isObserve = 1

                    self.skymap[key].targetCoverNum = 1

                    self.skymap[key].maxCoverNum = 2

                    # self.skymap_id.append(key)

                else:

                    del self.skymap[key]
@@ -231,7 +233,7 @@ class SkyMap(object): 
                coor.down_neighbour = -1

        gc.collect()



    def addSkyUnit(self, lonList=None, latList=None, idBase=100000000, obsTyp=0):

    def addSkyUnit(self, lonList=None, latList=None, idBase=100000000, obsTyp=0, targetCoverNum=2, maxCoverNum=2):

        if lonList is None or latList is None:

            print("Add sky unit failed!!!!!!")

            return
@@ -244,6 +246,9 @@ class SkyMap(object): 
                ecl_lat=lonlat[1], ecl_lon=lonlat[0], id=idBase+i)

            coor.obsType = obsTyp

            self.skymap[id] = coor

            self.isObserve = 1  # 0表示不可观测； 1表示可观测

            self.targetCoverNum = targetCoverNum  # 该天区需要观测次数

            self.maxCoverNum = maxCoverNum





# figure()

survey_sim/strategies/SurveySim_MSC_mpi.py

+76 −41
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@ 
Author: Zhang Xin zhangx@bao.ac.cn

Date: 2024-09-19 15:53:27

LastEditors: Zhang Xin zhangx@bao.ac.cn

LastEditTime: 2025-07-09 10:11:22

LastEditTime: 2025-07-22 15:45:23

FilePath: /CSST_Survey/survey_sim/strategies/SurveySim_MSC.py

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

'''
@@ -98,7 +98,8 @@ def fast_interpolate(x_points, y_points, query_points, kind='linear', fill_value 
def calulateSearchPointWeight_skyList_sigle(skymap=None, sky_attr_tab=None, sat_grcs_car=np.array([0, 0, 0]), sun_grcs_car=np.array([0, 0, 0]), sat_ecl_car=np.array([0, 0, 0]), sun_ecl_car=np.array([0, 0, 0]), moon_ecl_car=np.array([0, 0, 0]), dist_sun=1.5e8, lastPointing=np.array([0, 0, 0]), lastSatAtt=None, satOrbit_interp=None, curTime=2459766, endTime=2459766., battery_q=12000, subPoint_interp=None, curTime_yr=2023., lastPointingCar=np.array([0, 0, 0]), isInSunside=True, next_shadow_s=2459766., next_shadow_e=2459766., surveyCons=None, cmgCons=None, shadowTime=[123., 124.]):

    result_list = np.zeros([len(sky_attr_tab[:, 0]), 8], dtype=float)

    for i in np.arange(0, len(sky_attr_tab[:, 0]), 1):

        result_list[i] = calulateSearchPointWeight(skymap, sky_attr_tab, i, sun_grcs_car, sun_ecl_car, sat_ecl_car, curTime, endTime, battery_q, curTime_yr, lastPointingCar, isInSunside, next_shadow_s, next_shadow_e, surveyCons, shadowTime)

        result_list[i] = calulateSearchPointWeight(skymap, sky_attr_tab, i, sun_grcs_car, sun_ecl_car, sat_ecl_car, curTime,

                                                   endTime, battery_q, curTime_yr, lastPointingCar, isInSunside, next_shadow_s, next_shadow_e, surveyCons, shadowTime)



    valid_ids = result_list[:, 0] != MAX_VALUE

    result_list_ = result_list[valid_ids]
@@ -120,6 +121,7 @@ def calulateSearchPointWeight_skyList_sigle(skymap=None, sky_attr_tab=None, sat_ 
    # gc.collect()

    return min_result





def calulateSearchPointWeight(skymap, sky_attr_tab, idx,  sun_grcs_car, sun_ecl_car, sat_ecl_car, curTime, endTime, battery_q, curTime_yr, lastPointingCar, isInSunside, next_shadow_s, next_shadow_e, surveyCons, shadowTime):



    # idx = cuda.grid(1)
@@ -141,7 +143,7 @@ def calulateSearchPointWeight(skymap, sky_attr_tab, idx, sun_grcs_car, sun_ecl_ 
    obseveSkyUnit = skymap[s_id]

    # obseveSkyUnit.weight = 10000  # 初始权重

    # weight = obseveSkyUnit.weight

    weight = 1000  # 初始权重

    weight = 20000  # 初始权重



    tAngle = transAngle
@@ -153,7 +155,6 @@ def calulateSearchPointWeight(skymap, sky_attr_tab, idx, sun_grcs_car, sun_ecl_ 
    #     # self.local_fail_out[0] += 1.0

    #     return MAX_VALUE, s_id, MAX_VALUE, MAX_VALUE, MAX_VALUE, MAX_VALUE, MAX_VALUE, 0



 

    # sat_ecl_car_Atrans = sat_ecl_car_tTime



    # isObscureEarthFlag, angleValue = obscure_constraint.IsObscureByEarth(
@@ -229,7 +230,6 @@ def calulateSearchPointWeight(skymap, sky_attr_tab, idx, sun_grcs_car, sun_ecl_ 
        #     if panel_is_ok != 1:

        #         return MAX_VALUE, s_id, MAX_VALUE, MAX_VALUE, MAX_VALUE, MAX_VALUE, MAX_VALUE, 4





        if obscure_constraint.IsInSunSide(sat=sat_ecl_car_extTime, sun=sun_ecl_car) == 1:

            # panel_is_ok = (int)(

            #     value_sun_angle_to_nomal >= thres_for_sun_angle)
@@ -291,6 +291,7 @@ def calulateSearchPointWeight(skymap, sky_attr_tab, idx, sun_grcs_car, sun_ecl_ 


    return weight, s_id, tTime + exTime, exTime, tAngle, cmg_use, energyT1, 99





class SurveySim(object):

    def __init__(self, time_start=2459767.0, time_end=2463416.0, bettery_level=97200., cmgCons=None, dir=None, logName=None, outFn=None):
@@ -337,7 +338,8 @@ class SurveySim(object): 


    def saveSkyMap(self, ctime):

        skydat_path = resource_filename("survey_sim", "strategies/data")

        skyMap_save_data = os.path.join(skydat_path, "skyMap_Init" + str(np.round(ctime,6)) + ".pkl")

        skyMap_save_data = os.path.join(

            skydat_path, "skyMap_Init" + str(np.round(ctime, 6)) + ".pkl")

        smF = open(skyMap_save_data, "wb")

        pickle.dump(self.skymap, smF, protocol=pickle.HIGHEST_PROTOCOL)
@@ -387,7 +389,8 @@ class SurveySim(object): 
            outInfo = InfoOutput(

                dir=self.outDir, logger_filename=self.logger_filename, info_out_filename=self.info_out_filename)

        for start_end_t_beta in self.betaCons.MSC_time:

            if start_end_t_beta[1] - start_end_t_beta[0] < 0.0002055123903420137: ## 时间间隔小于300s的时间段跳过

            # 时间间隔小于300s的时间段跳过

            if start_end_t_beta[1] - start_end_t_beta[0] < 0.0002055123903420137:

                continue

            saaSeg = saa_constraint.SAA_constraint(

                startTime=start_end_t_beta[0], endTime=start_end_t_beta[1], orbitData=self.orbitDat)
@@ -520,12 +523,10 @@ class SurveySim(object): 
            #     curTime, self.orbitDat, self.eph

            # )





        sat, satVel, orbitDataId = locateSat1(

            time=curTime, OrbitData=self.orbitDat, startId=orbitDataId

        )





        if lastSatAtt is None or curPSky is None:

            lastSatAtt = satAttitude(curTime=curTime, sat=sat, q=None)

            sat_ecl_car = lastSatAtt.sat_ecl
@@ -640,7 +641,8 @@ class SurveySim(object): 
        else:

            # skymap_ids, value_sun_angle_to_nomals, cut_ids = obscure_constraint.IsObscureSolarPlane(

            #     skymap_ids=skymap_ids, sun_ecl_car=sun_ecl_car, dist_sun=dist_sun, skyMap=self.skymap)

            skymap_ids, value_sun_angle_to_nomals, cut_ids = obscure_constraint.IsObscureSolarPlane_batch(skymap_ids= skymap_ids, sun_ecl_car=sun_ecl_car, dist_sun=dist_sun, skyMap=self.skymap)

            skymap_ids, value_sun_angle_to_nomals, cut_ids = obscure_constraint.IsObscureSolarPlane_batch(

                skymap_ids=skymap_ids, sun_ecl_car=sun_ecl_car, dist_sun=dist_sun, skyMap=self.skymap)

            # print(np.sum(skymap_ids-skymap_ids1))



        pointings = pointings[cut_ids]
@@ -686,7 +688,8 @@ class SurveySim(object): 
        after_expTime_JD = after_tTime_JD + expTimes/86400.



        cmg_uses = self.cmgCons.get_cmg_use_array(tAngle_array)

        cmg_total_uses = self.cmgCons.computeCMGTotal_newComein_array(curTime, tTime_array, expTimes, cmg_uses)

        cmg_total_uses = self.cmgCons.computeCMGTotal_newComein_array(

            curTime, tTime_array, expTimes, cmg_uses)



        ids_select1 = cmg_total_uses <= self.surveyCons.CMG_THRES
@@ -698,7 +701,8 @@ class SurveySim(object): 
            time_list, satInfo_list[:, 2], after_tTime_JD[ids_select1], kind='linear', fill_value=np.nan)



        isObscureE_array = obscure_constraint.IsObscureByEarth_array(

            sat=np.stack((sat_ecl_x_tTime, sat_ecl_y_tTime, sat_ecl_z_tTime),axis=1),

            sat=np.stack((sat_ecl_x_tTime, sat_ecl_y_tTime,

                         sat_ecl_z_tTime), axis=1),

            sun=sun_ecl_car,

            p=pointings[:, 0:3][ids_select1],

            constr=self.surveyCons,
@@ -714,7 +718,8 @@ class SurveySim(object): 
            time_list, satInfo_list[:, 2], after_expTime_JD[ids_select1][ids_select2], kind='linear', fill_value=np.nan)



        isObscureE_array_ = obscure_constraint.IsObscureByEarth_array(

            sat=np.stack((sat_ecl_x_expTime, sat_ecl_y_expTime, sat_ecl_z_expTime),axis=1),

            sat=np.stack((sat_ecl_x_expTime, sat_ecl_y_expTime,

                         sat_ecl_z_expTime), axis=1),

            sun=sun_ecl_car,

            p=pointings[:, 0:3][ids_select1][ids_select2],

            constr=self.surveyCons,
@@ -724,7 +729,8 @@ class SurveySim(object): 
        sky_attr_tab = np.stack(

            (skymap_ids[ids_select1][ids_select2], value_sun_angle_to_nomals[ids_select1][ids_select2], tAngle_array[ids_select1][ids_select2], tTime_array[ids_select1][ids_select2], expTimes[ids_select1][ids_select2], after_tTime_JD[ids_select1][ids_select2], after_expTime_JD[ids_select1][ids_select2], sat_ecl_x_tTime[ids_select2], sat_ecl_y_tTime[ids_select2], sat_ecl_z_tTime[ids_select2], sat_ecl_x_expTime, sat_ecl_y_expTime, sat_ecl_z_expTime, pointings[:, 0][ids_select1][ids_select2], pointings[:, 1][ids_select1][ids_select2], pointings[:, 2][ids_select1][ids_select2], pointings[:, 3][ids_select1][ids_select2], pointings[:, 4][ids_select1][ids_select2], cmg_uses[ids_select1][ids_select2], cmg_total_uses[ids_select1][ids_select2]), axis=1)



        ids_select3 = (isObscureE_array_ == 0) & (after_tTime_JD[ids_select1][ids_select2]<=endTime) & (after_expTime_JD[ids_select1][ids_select2]<=endTime)

        ids_select3 = (isObscureE_array_ == 0) & (after_tTime_JD[ids_select1][ids_select2] <= endTime) & (

            after_expTime_JD[ids_select1][ids_select2] <= endTime)



        sky_attr_tab = sky_attr_tab[ids_select3]

        # print("Time:", curTime, isInSunSide,
@@ -832,7 +838,7 @@ class SurveySim(object): 
                battery_q = self.surveyCons.BATTERY_LOW

                print("ERROR:-----Energy not balance")

            self.bettery_level = battery_q

            if returnTime - lastObsEndTime > 0.003472222222222222: ##空转300s后指向指到天顶，默认没有消耗cmg，在期间慢慢转过去

            if returnTime - lastObsEndTime > 0.003472222222222222:  # 空转300s后指向指到天顶，默认没有消耗cmg，在期间慢慢转过去

                return returnTime, returnTime, None, None, None, None, None, orbitDataId, isInSunSide, None, None, None, self.cmgCons.total_cmg, self.bettery_level, lastObsEndTime

            else:

                return returnTime, returnTime, curPSky, None, None, None, None, orbitDataId, isInSunSide, None, None, lastSatAtt, self.cmgCons.total_cmg, self.bettery_level, lastObsEndTime
@@ -903,16 +909,45 @@ class SurveySim(object): 
        return weight





def addSkyList_inputRADEC(radecs):

    from astropy.coordinates import SkyCoord

    from astropy import units as u

    ra = radecs[:, 0] * u.deg      # 赤经 (Right Ascension)

    dec = radecs[:, 1] * u.deg     # 赤纬 (Declination)

    distance = 1 * u.AU  # 距离（假设单位是 AU，如果是 km 需转换）



    # 创建赤道坐标

    equatorial_coord = SkyCoord(

        ra=ra, dec=dec, distance=distance, frame='icrs')



    # 转换为黄道坐标（默认使用 'mean' 黄道坐标系）

    ecliptic_coord = equatorial_coord.transform_to('geocentricmeanecliptic')



    return ecliptic_coord.lon.value, ecliptic_coord.lat.value





if __name__ == '__main__':

    if len(sys.argv) < 6:

        print("Error: input parameter error, must be : start time, end time, out dir, logName, outfile name")

    print(sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5])



    ra = np.array([227.242, 132.8726, 15.0183, 100.4065, 210.02, 260.0684, 114.1066,

                  158.7706, 37.389, 202.0091, 186.7454, 247.7722, 152.1146, 189.5908, 168.3627])

    dec = np.array([67.2221, 63.1335, -33.7186, -50.9593, 14.5135, 57.9185, -57.9991,

                   51.9479, -79.3089, 33.5521, 23.9069, 12.7852, 12.3059, -40.9043, 22.1529])

    radecs = np.column_stack((ra, dec))



    ecl_lons, ecl_lats = addSkyList_inputRADEC(radecs)



    # dir="/Users/zhangxin/Work/SurveyPlan/CSST_Survey", logName="MSC_survey_2459767_log.dat", outFn="MSC_survey_2459767_log_outInfo.dat"

    # ss = SurveySim(time_start=2459767, time_end=2459767+10,

    #                bettery_level=97200., cmgCons=None, dir=sys.argv[1], logName=sys.argv[2], outFn=sys.argv[3])

    ss = SurveySim(time_start=float(sys.argv[1]), time_end=float(sys.argv[2]),

                   bettery_level=97200., cmgCons=None, dir=sys.argv[3], logName=sys.argv[4], outFn=sys.argv[5])



    ss.skymap.addSkyUnit(lonList=ecl_lons, latList=ecl_lats,

                         idBase=100000000, obsTyp=9, targetCoverNum=2, maxCoverNum=2)



    t1 = time.time()

    tt = ss.do_observe_plan(jump_time_unObs=150)

    # ss.saveSkyMap(time_end)

survey_sim/strategies/_utils_cuda.py

+15 −3
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@ 
Author: Zhang Xin zhangx@bao.ac.cn

Date: 2024-10-14 13:54:24

LastEditors: Zhang Xin zhangx@bao.ac.cn

LastEditTime: 2025-06-17 09:36:22

LastEditTime: 2025-07-22 14:53:03

FilePath: /CSST_Survey/survey_sim/strategies/_utils.py

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

'''
@@ -125,7 +125,7 @@ def calculateEnergyInSunback(battery_q=9000, curTime=2459777., exEndTime=2459777 
# 根据CMG调整权重





def setWeight_CMG(tAngle=1.0, curTime_yr=0, weight=0.0, surveyCons=None):

def setWeight_CMG(tAngle=1.0, curTime_yr=0, weight=0.0, surveyCons=None, obs_type=9):

    # double angles[8] = {5,10,20,35,45,75,90,180};

    angles = np.array([5, 10, 20, 35, 45, 75, 90, 180])
@@ -140,6 +140,10 @@ def setWeight_CMG(tAngle=1.0, curTime_yr=0, weight=0.0, surveyCons=None): 
                    weight = weight - 50.0 * np.exp(

                        -(tAngle - 5) * (tAngle - 5) / 5

                    )

                if obs_type == 9:

                    weight = weight - 50.0 * np.exp(

                        -(tAngle - 5) * (tAngle - 5) / 5

                    )

            elif i == 1:

                weight = weight / (-2 * tAngle + 45)

                # if obseveSkyUnit.obsType == 1:
@@ -149,8 +153,16 @@ def setWeight_CMG(tAngle=1.0, curTime_yr=0, weight=0.0, surveyCons=None): 
                    weight = weight - 50.0 * np.exp(

                        -(tAngle - 5) * (tAngle - 5) / 5

                    )

                if obs_type == 9:

                    weight = weight - 50.0 * np.exp(

                        -(tAngle - 5) * (tAngle - 5) / 5

                    )

            elif i == 2:

                weight = weight / (-tAngle + 35.0)

                if obs_type == 9:

                    weight = weight - 50.0 * np.exp(

                        -(tAngle - 5) * (tAngle - 5) / 5

                    )

            elif i == 3:

                weight = weight / (-0.333333 * tAngle + 21.666667)

            elif i == 4: