call_dayperi.F

路径

LMDZ.MARS\libf\phymars\call_dayperi.F

所属目录/模块

libf\phymars

文件定位

该文件定义 call_dayperi_mod 模块和一个轨道日历换算例程 call_dayperi。例程输入近日点太阳经度 Lsperi、轨道偏心率 e_elips 和火星年长度 year_day，输出从北半球春分 Ls=0 起算的近日点日期 dayperi（sol）。

确认调用方是 1D 初始化模块 init_testphys1d_mod.F90：它从配置读取 Lsperi（单位为度）和 eccentric，把 Lsperi 转成弧度后调用 call_dayperi(Lsperi,eccentric,peri_day,year_day)，再把结果写入 planete_h 的 peri_day。3D 主路径中也有 planete_h:lsp2solp 可根据 lsperi 计算近日点 sol，因此本文件主要服务于 1D/testphys 初始化的轨道参数换算。

定义的符号

依赖的模块

调用的关键例程

输入

输出

共享状态与副作用

• 不读写 module 变量或 common block。
• 不写文件、不输出日志、不分配内存。
• 通过 dayperi 形参把结果返回给调用方。
• 形参未声明 intent，但源码注释明确 Lsperi、e_elips、year_day 为输入，dayperi 为输出。

核心逻辑

1. 用 pi = 2.*asin(1.) 计算圆周率。
2. 计算椭圆轨道变换中的两个中间量：
   • x1 = sqrt((1-e_elips)/(1+e_elips))
   • x2 = e_elips*sqrt(1-e_elips**2)
3. 用源码公式计算从 Ls=0 到近日点的日期：

dayperi = 0.5*(year_day/pi) *
          ( 2*atan(x1*tan(0.5*Lsperi))
            - x2*sin(Lsperi)/(1+e_elips*cos(Lsperi)) )

4. 若结果为负，则加一个火星年：dayperi = dayperi + year_day，把日期放回 [0, year_day) 的年度范围。

伪代码

call_dayperi(Lsperi, e_elips, dayperi, year_day):
    pi = 2 * asin(1)
    x1 = sqrt((1 - e_elips) / (1 + e_elips))
    x2 = e_elips * sqrt(1 - e_elips**2)

    anomaly_part = 2 * atan(x1 * tan(0.5 * Lsperi))
    correction = x2 * sin(Lsperi) / (1 + e_elips * cos(Lsperi))
    dayperi = 0.5 * (year_day / pi) * (anomaly_part - correction)

    if dayperi < 0:
        dayperi = dayperi + year_day

参与的主题流程

写法特点

• .F 文件内使用 module/contains，代码为固定格式风格。
• 无 implicit none 之外的模块依赖，公式完全由输入标量和 intrinsic 组成。
• 未使用 intent(in/out) 标注，复现接口时应按源码注释理解输入输出方向。
• 输入 Lsperi 必须是弧度；init_testphys1d_mod.F90 在调用前显式做 Lsperi = Lsperi*pi/180.。
• 只在 dayperi < 0 时加 year_day，没有对 dayperi >= year_day 的情形做取模。

复现要点

• 不要把 Lsperi 的配置单位（度）直接传入本例程；调用前必须转弧度。
• e_elips 应满足 0 <= e_elips <= 1。1D 初始化在调用前检查 eccentric 范围。
• 公式使用普通单精度 real，不是 double precision。
• atan(x1*tan(0.5*Lsperi)) 含有 tan 的象限/分支行为；复现时应保持源码表达式，不要随意替换成不同象限约定的 atan2，除非重新核对全年范围。
• 负值结果只加一次 year_day；对超出常规范围的 Lsperi 输入，源码没有完整归一化。

待确认

• init_testphys1d_mod.F90 中检查 Lsperi 范围时源码写成 if (eccentric < 0. .or. eccentric > 360.)，看起来应检查 Lsperi 而不是 eccentric；本页只记录 call_dayperi 行为，不把该上游疑点写成事实。
• 本公式与 planete_h:lsp2solp 在所有 Lsperi 象限上是否完全一致，尚未逐值比对。
• 对 Lsperi 超出 [0,2*pi] 的输入，当前例程没有归一化，调用方应保证范围。

相关页面

• testphys1d-config：1D 配置和初始化流程，包含轨道参数输入。
• planete_h.md：保存 year_day、peri_day、lsperi、e_elips 等轨道状态。
• solarlong.md：用 peri_day 把 sol 转换为太阳经度。
• orbite：用 lsperi、偏心率和太阳经度计算日火距离与太阳赤纬。