diff --git a/config/config_surveysim.yaml b/config/config_surveysim.yaml index bfceb9d5e2ed1701b002c7654a9f6c158317d4a0..b803be5bae094970d9aca941590c3d58784ac61b 100644 --- a/config/config_surveysim.yaml +++ b/config/config_surveysim.yaml @@ -1,5 +1,5 @@ out_path_param: - wokr_dir: /Users/zhangxin/Work/SurveyPlan/CSST_Survey # define the folder of output file + work_dir: /Users/zhangxin/Work/SurveyPlan/CSST_Survey # define the folder of output file out_filename_prefix: test_2461876_1 # output file prefix, 3 files: prefix_plan.json, prefix_info.dat, prefix_log.log # input_survey_data ---- whether to use input survey data; if yes, set the file name (including path), otherwise set null diff --git a/config/config_surveysim_fdw.yaml b/config/config_surveysim_fdw.yaml index ca4fc348e40496a1588ad1c7d81a5ae09ce43ca1..bdb16fe1a8a6d75e2fa89404998f3f46a0bccf1f 100644 --- a/config/config_surveysim_fdw.yaml +++ b/config/config_surveysim_fdw.yaml @@ -1,5 +1,5 @@ out_path_param: - wokr_dir: D:/Projects/CSST/dev/csst-planning/upload/sim/output # define the folder of output file + work_dir: D:/Projects/CSST/dev/csst-planning/upload/sim/output # define the folder of output file out_filename_prefix: test_2461876_1 # output file prefix, 3 files: prefix_plan.json, prefix_info.dat, prefix_log.log # input_survey_data ---- whether to use input survey data; if yes, set the file name (including path), otherwise set null diff --git a/run_survey_sim.py b/run_survey_sim.py index cb11b07f3c8160c596bb3842984b73d6da283be5..6727fe1122f7d2fafcbbca814908805b6677bfbe 100644 --- a/run_survey_sim.py +++ b/run_survey_sim.py @@ -158,6 +158,6 @@ if __name__ == '__main__': if config_data.get("input_survey_data") is not None: input_sim_data = config_data.get("input_survey_data") sim_from_input(jsonfn=input_sim_data, - outDir=config_data["out_path_param"].get("wokr_dir"), info_out_filename_pre=config_data["out_path_param"].get("out_filename_prefix")) + outDir=config_data["out_path_param"].get("work_dir"), info_out_filename_pre=config_data["out_path_param"].get("out_filename_prefix")) else: run_survey_sim(config=config_data) diff --git a/survey_sim/config/infoOutput.py b/survey_sim/config/infoOutput.py index 07b526846b766f696173139d41ea41a2b2182e06..fbcd4181f182c6ffc8a6f6257ba7e20d23edbd97 100644 --- a/survey_sim/config/infoOutput.py +++ b/survey_sim/config/infoOutput.py @@ -183,7 +183,7 @@ class InfoOutput(object): data['simresult']['is_in_sun'] = isInSunSide data['simresult']['cmg'] = cmg data['simresult']['charge'] = energy - data['simresult']['skyid'] = sky_ID + data['simresult']['skyid'] = int(sky_ID) # data['simresult']['board_angle'] = d1['boardAngle'] data['simresult']['quad_0'] = sat_q[0] data['simresult']['quad_i'] = sat_q[1] diff --git a/survey_sim/constraints/surveyConstraint.py b/survey_sim/constraints/surveyConstraint.py index 29c545239ae35be983678132d2cc052c93e4dc4f..9867b2fb5a71950608b6e921fa18c276cd343cf8 100644 --- a/survey_sim/constraints/surveyConstraint.py +++ b/survey_sim/constraints/surveyConstraint.py @@ -5,13 +5,13 @@ class surveyConstraint(object): def __init__( self, sunAngle=50, - moonAgle=40, + moonAngle=40, earthLLimbAngle_limit=80, earthLLimbAngle_plus=70, earthDLimgAngle=30, ): self.sun_los_angle_cos = math.cos(sunAngle * math.pi / 180.0) - self.moon_los_angle_cos = math.cos(moonAgle * math.pi / 180.0) + self.moon_los_angle_cos = math.cos(moonAngle * math.pi / 180.0) # 约为70°是根据400km轨道，地球平均半径计算的，切线和卫星矢量夹角 self.satZenith_angle_light_min_cos = math.cos( (180 - 70 - earthLLimbAngle_limit) * math.pi / 180.0 diff --git a/survey_sim/satOrbit/locate_sat.py b/survey_sim/satOrbit/locate_sat.py index aa765c7b18a3c2b5fd9b8ab638609f6808ddfa98..088eb5281bf48c6396175f943241789a33a134cf 100644 --- a/survey_sim/satOrbit/locate_sat.py +++ b/survey_sim/satOrbit/locate_sat.py @@ -68,6 +68,59 @@ def loadSatOrbitDat(datDir=None): return oData +def findTimeIndex(time, time_col, startId=0): + """ + 查找时间在轨道数据中的索引位置（使用二分查找优化，并利用 startId） + + 参数: + time: 目标时间 + time_col: 时间列数组 + startId: 起始搜索索引（用于优化连续调用） + + 返回: + int: 找到的索引，如果未找到返回 -1 + """ + orbDataLen = len(time_col) + + if orbDataLen < 2: + return -1 + + # 如果 startId 有效，先检查 startId 附近是否包含目标时间 + if startId > 0 and startId < orbDataLen - 1: + # 检查 startId 附近的时间范围（前后各检查一定范围） + search_range = 100 # 检查范围 + start_check = max(0, startId - search_range) + end_check = min(orbDataLen - 1, startId + search_range) + + # 检查 startId 附近的时间是否包含目标时间 + t_start = time_col[start_check].item() + t_end = time_col[end_check].item() + + if t_start <= time <= t_end: + # 在 startId 附近使用线性搜索（因为范围小，线性搜索更快） + for i in range(start_check, end_check): + if i >= orbDataLen - 1: + break + t1 = time_col[i].item() + t2 = time_col[i + 1].item() + if time >= t1 and time < t2: + return i + + # 如果 startId 附近没找到，使用全局二分查找 + # 使用 np.searchsorted 进行二分查找 + idx = np.searchsorted(time_col, time, side='right') + + # 检查是否找到有效的时间区间 + if idx > 0 and idx < orbDataLen: + i = idx - 1 + t1 = time_col[i].item() + t2 = time_col[i + 1].item() + if time >= t1 and time < t2: + return i + + return -1 + + def locateSat(time=2459766.0, OrbitData=None, startId=0, orbDataLen=-1): if orbDataLen <= 0 and OrbitData is not None: @@ -75,78 +128,83 @@ def locateSat(time=2459766.0, OrbitData=None, startId=0, orbDataLen=-1): satPos = np.zeros(3, dtype=np.float64) satVel = np.zeros(3, dtype=np.float64) nextSid = startId - for i in np.arange(startId, orbDataLen - 1, 1): - t1 = OrbitData[i, 0] - t2 = OrbitData[i + 1, 0] - if time >= t1 and time < t2: - if (t2 - t1) > 130.0 / 86400.0: # should be 120s, for error, set 120+10s - print("Warning: time interval is not 120s!!!!!!!!!") - # return satPos-1, 0 - # break - x1 = OrbitData[i, 1] - y1 = OrbitData[i, 2] - z1 = OrbitData[i, 3] - x2 = OrbitData[i + 1, 1] - y2 = OrbitData[i + 1, 2] - z2 = OrbitData[i + 1, 3] - - line1 = np.sqrt(x1 * x1 + y1 * y1 + z1 * z1) - line2 = np.sqrt(x2 * x2 + y2 * y2 + z2 * z2) - theta = np.arccos((x1 * x2 + y1 * y2 + z1 * z2) / (line1 * line2)) - theta1 = (time - t1) / (t2 - t1) * theta - theta2 = theta - theta1 - line = (t2 - time) / (t2 - t1) * line1 + \ - (time - t1) / (t2 - t1) * line2 - # // double ef = sin(theta2)/sin(theta1); - - x0 = np.sin(theta2) * x1 / line1 + np.sin(theta1) * x2 / line2 - y0 = np.sin(theta2) * y1 / line1 + np.sin(theta1) * y2 / line2 - z0 = np.sin(theta2) * z1 / line1 + np.sin(theta1) * z2 / line2 - line_ = np.sqrt(x0 * x0 + y0 * y0 + z0 * z0) - - satPos[0] = x0 * line / line_ - satPos[1] = y0 * line / line_ - satPos[2] = z0 * line / line_ - - xdot1 = OrbitData[i, 4] - ydot1 = OrbitData[i, 5] - zdot1 = OrbitData[i, 6] - - xdot2 = OrbitData[i + 1, 4] - ydot2 = OrbitData[i + 1, 5] - zdot2 = OrbitData[i + 1, 6] - - ldot1 = np.sqrt(xdot1 * xdot1 + ydot1 * ydot1 + zdot1 * zdot1) - ldot2 = np.sqrt(xdot2 * xdot2 + ydot2 * ydot2 + zdot2 * zdot2) - thetadot = math.acos( - (xdot1 * xdot2 + ydot1 * ydot2 + zdot1 * zdot2) / (ldot1 * ldot2) - ) - - ldot = (t2 - time) / (t2 - t1) * ldot1 + \ - (time - t1) / (t2 - t1) * ldot2 - - if time != t1: - thetadot1 = (time - t1) / (t2 - t1) * thetadot - thetadot2 = thetadot - thetadot1 - efdot = math.sin(thetadot2) / math.sin(thetadot1) - - xdot0 = efdot * xdot1 / ldot1 + xdot2 / ldot2 - ydot0 = efdot * ydot1 / ldot1 + ydot2 / ldot2 - zdot0 = efdot * zdot1 / ldot1 + zdot2 / ldot2 - else: - xdot0 = xdot1 / ldot1 - ydot0 = ydot1 / ldot1 - zdot0 = zdot1 / ldot1 - - ldot_ = np.sqrt(xdot0 * xdot0 + ydot0 * ydot0 + zdot0 * zdot0) - - satVel[0] = xdot0 * ldot / ldot_ - satVel[1] = ydot0 * ldot / ldot_ - satVel[2] = zdot0 * ldot / ldot_ + # 使用独立的查找函数 + time_col = OrbitData[:, 0] + i = findTimeIndex(time, time_col, startId) + + # 如果找到匹配的时间区间，进行插值计算 + if i >= 0: + t1 = OrbitData[i, 0].item() + t2 = OrbitData[i + 1, 0].item() + if (t2 - t1) > 130.0 / 86400.0: # should be 120s, for error, set 120+10s + print("Warning: time interval is not 120s!!!!!!!!!") + # return satPos-1, 0 + # break + + x1 = OrbitData[i, 1] + y1 = OrbitData[i, 2] + z1 = OrbitData[i, 3] + x2 = OrbitData[i + 1, 1] + y2 = OrbitData[i + 1, 2] + z2 = OrbitData[i + 1, 3] + + line1 = np.sqrt(x1 * x1 + y1 * y1 + z1 * z1) + line2 = np.sqrt(x2 * x2 + y2 * y2 + z2 * z2) + theta = np.arccos((x1 * x2 + y1 * y2 + z1 * z2) / (line1 * line2)) + theta1 = (time - t1) / (t2 - t1) * theta + theta2 = theta - theta1 + line = (t2 - time) / (t2 - t1) * line1 + \ + (time - t1) / (t2 - t1) * line2 + # // double ef = sin(theta2)/sin(theta1); + + x0 = np.sin(theta2) * x1 / line1 + np.sin(theta1) * x2 / line2 + y0 = np.sin(theta2) * y1 / line1 + np.sin(theta1) * y2 / line2 + z0 = np.sin(theta2) * z1 / line1 + np.sin(theta1) * z2 / line2 + line_ = np.sqrt(x0 * x0 + y0 * y0 + z0 * z0) + + satPos[0] = x0 * line / line_ + satPos[1] = y0 * line / line_ + satPos[2] = z0 * line / line_ + + xdot1 = OrbitData[i, 4] + ydot1 = OrbitData[i, 5] + zdot1 = OrbitData[i, 6] + + xdot2 = OrbitData[i + 1, 4] + ydot2 = OrbitData[i + 1, 5] + zdot2 = OrbitData[i + 1, 6] + + ldot1 = np.sqrt(xdot1 * xdot1 + ydot1 * ydot1 + zdot1 * zdot1) + ldot2 = np.sqrt(xdot2 * xdot2 + ydot2 * ydot2 + zdot2 * zdot2) + thetadot = math.acos( + (xdot1 * xdot2 + ydot1 * ydot2 + zdot1 * zdot2) / (ldot1 * ldot2) + ) + + ldot = (t2 - time) / (t2 - t1) * ldot1 + \ + (time - t1) / (t2 - t1) * ldot2 + + if time != t1: + thetadot1 = (time - t1) / (t2 - t1) * thetadot + thetadot2 = thetadot - thetadot1 + efdot = math.sin(thetadot2) / math.sin(thetadot1) + + xdot0 = efdot * xdot1 / ldot1 + xdot2 / ldot2 + ydot0 = efdot * ydot1 / ldot1 + ydot2 / ldot2 + zdot0 = efdot * zdot1 / ldot1 + zdot2 / ldot2 + else: + xdot0 = xdot1 / ldot1 + ydot0 = ydot1 / ldot1 + zdot0 = zdot1 / ldot1 + + ldot_ = np.sqrt(xdot0 * xdot0 + ydot0 * ydot0 + zdot0 * zdot0) + + satVel[0] = xdot0 * ldot / ldot_ + satVel[1] = ydot0 * ldot / ldot_ + satVel[2] = zdot0 * ldot / ldot_ + + nextSid = i - nextSid = i - break return satPos, satVel, nextSid diff --git a/survey_sim/strategies/SurveySim_MSC_mpi.py b/survey_sim/strategies/SurveySim_MSC_mpi.py index 11f889ff0660442c0a52ef8fb8b705aee8af1754..c5e885ce3b3daac98aeb376a2b108cc96f5acc7f 100644 --- a/survey_sim/strategies/SurveySim_MSC_mpi.py +++ b/survey_sim/strategies/SurveySim_MSC_mpi.py @@ -267,13 +267,13 @@ def calulateSearchPointWeight(skymap, sky_attr_tab, idx, sun_grcs_car, sun_ecl_ class SurveySim(object): - def __init__(self, time_start=2459767.0, time_end=2463416.0, bettery_level=97200., cmgCons=None, dir=None, outFn=None, obsStartNum=1, config=None, beta_flag=True): + def __init__(self, time_start=2459767.0, time_end=2463416.0, bettery_level=97200., cmgCons=None, dir=None, outFn=None, obsStartNum=1, config=None, beta_flag=True, orbit_data=None): if config is not None: self.plan_startTime = config["sim_param"].get("start_time") self.plan_endTime = config["sim_param"].get("end_time") self.bettery_level = config["sim_param"].get("bettery_level") - self.outDir = config["out_path_param"].get("wokr_dir") + self.outDir = config["out_path_param"].get("work_dir") self.info_out_filename_pre = config["out_path_param"].get( "out_filename_prefix") self.config = config @@ -298,10 +298,11 @@ class SurveySim(object): self.initSkyMap(surveyCons=self.surveyCons) # print("memory analysis 2 : %f" % ( # (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024))) - tmpData = loadSatOrbitDat() - ids1 = tmpData[:, 0] > self.plan_startTime - 0.2 - # ids2 = tmpData[:, 0] < self.plan_endTime + 2 - self.orbitDat = tmpData[ids1] + if orbit_data is None: + orbit_data = loadSatOrbitDat() + ids1 = orbit_data[:, 0] > self.plan_startTime - 0.2 + # ids2 = orbit_data[:, 0] < self.plan_endTime + 2 + self.orbitDat = orbit_data[ids1] # self.orbitDat = loadSatOrbitDat() # print("memory analysis 3 : %f" % ( # (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)))