使用DE历表预报二十四节气（Fortran77）

铺垫

希望你（不）是在天体测量这门课上接到相关使用历表的任务的。
太阳系历表，具体介绍找了一圈没找到。总之它记录了天体的位置、速度等信息。
二十四节气太阳所对应的黄经：

然后，经过一系列操作拿到了历表文件、IAU的SOFA程序，就可以开始了。
编译环境：Simply Fortran

代码

直接上主程序：

      program main
      implicit none
      integer year,month,day,hour,minute,J
      real second
      DOUBLE PRECISION RV(6),R(3),V(3),RVe(6),RVs(6),velocitys
!RV:record the data(position and velocity)
      DOUBLE PRECISION tJD(2),tJD0(2),JD,JD0,ET(2),tJDLTDB(2)
      DOUBLE PRECISION step(2),predict,Elong,Length
!tJD,tJD0 is julian date, devided in 2 parts
!JD=tJD(1)+tJD(2)
      DOUBLE PRECISION tTAI(2),tTT(2),tDTR,tTDB(2),tUTC(2)
      DOUBLE PRECISION RNPB(3,3),RM(3,3),Ro(3,3)
      double precision S,BM1,PPR(3),AU,c,alpha
      c=2.99792458D8
      AU=1.49597870699626200D11
      
！定义黄赤交角数值，并计算赤道坐标系转黄道坐标系的旋转矩阵      
      alpha=23.4+1.0/30.0
      Ro(1,1)=1D0
      Ro(1,2)=0D0
      Ro(1,3)=0D0
      Ro(2,1)=0D0
      Ro(2,2)=cos(dacos(-1d0)*alpha/180.0)
      Ro(2,3)=sin(dacos(-1d0)*alpha/180.0)
      Ro(3,1)=0d0
      Ro(3,2)=-sin(dacos(-1d0)*alpha/180.0)
      Ro(3,3)=cos(dacos(-1d0)*alpha/180.0)
!**********************************************************
！设置计算起点（其实不用定义时分秒）
      year=2020;month=2;day=1
      hour=0;minute=0;second=0D0
！预测接下来366天的      
      predict=366.0
      step(1)=1.0/24.0/60.0
      step(2)=1.0/24.0/60.0/60.0/2.22
      Length=360.0/365.25/24.0/60.0
！将格里高利转化为儒略日      
      call iau_cal2jd(year,month,day,tJD0(1),tJD0(2),J)
      JD0=tJD0(1)+tJD0(2)
      JD=JD0
      
!**********************************************************      
!Begin the cycle
      DO while(JD<(JD0+predict))
      tJD(1)=floor(JD-0.5D0)+0.5D0
      tJD(2)=JD-tJD(1)
！由于历表时间系统为TDB，所以需要转化时间系统，TT用做读取岁差章动矩阵      
!UTC-TAI
      call iau_UTCTAI(tJD(1),tJD(2),tTAI(1),tTAI(2),J)
!TAI-TT
      call iau_TAITT(tTAI(1),tTAI(2),tTT(1),tTT(2),J)
!TT-TDB
      call iau_TTTDB(tTT(1),tTT(2),tDTR,tTDB(1),tTDB(2),J)

!光行时修正，作用为找到8分钟前太阳的位置，也可以直接减去8分钟
      call LTT(tTDB(1),tTDB(2),tJDLTDB(1),tJDLTDB(2))
      
      ET(1)=tJDLTDB(1)
      ET(2)=tJDLTDB(2)
!**********************************************************      
!这点需要注意：下面使用的testeph.f在历表文件中
!DPLEPH(t,target,center,RV)读历表，取出太阳相对地球的位置和速度
      call PLEPH(ET(1)+ET(2),3,11,RV)
      R(1)=RV(1);R(2)=RV(2);R(3)=RV(3)
      V(1)=RV(4);V(2)=RV(5);V(3)=RV(6)
      
!Change R to the opposite direction（因为我们要求的是“视黄经”）
      R(1)=-R(1);R(2)=-R(2);R(3)=-R(3)
!**********************************************************  
!光行差修正，为了调动AB函数，改变速度的单位，计算单位向量以及洛伦兹量 AU/DAY-nc
      V(1)=V(1)*AU/86400.0/c
      V(2)=V(2)*AU/86400.0/c
      V(3)=V(3)*AU/86400.0/c
!Normalized
      S=dsqrt(R(1)**2+R(2)**2+R(3)**2)
      R(1)=R(1)/S
      R(2)=R(2)/S
      R(3)=R(3)/S
      
      velocitys=V(1)**2+V(2)**2+V(3)**2
      BM1=dsqrt(1D0-velocitys)
      call iau_AB(R,V,S,BM1,R)
      
！利用SOFA的程序读出进动章动修正矩阵      
!precession-nutation matrix;jindongzhangdong
      call iau_PNM06a(tTT(1),tTT(2),RNPB)
      
！对位矢计算新坐标系下的各分量     
      R=matmul(RNPB,R)
      R=matmul(Ro,R)
      
!求黄经，即用xy分量求极坐标的角度
      call cosp(R(1),R(2),Elong)
      
！把时间变成北京时间；Judge函数为判断黄经是否满足一定条件，满足条件则输出时间以及节气名      
      JD=JD+8.0/24.0
      call Judge(Elong,JD)
      JD=JD-8.0/24.0
      
！如果黄经离15度的倍数还远，则时间加step（1）+step（2）约等于1分钟；如果接近了，每次增量改为step（2）约等于半秒      
      if((15.0-mod(Elong+15.0,15.0)).GE.Length) then
      JD=JD+step(1)
      end if
    
      JD=JD+step(2)
      
      end DO
      
      end program

Judge子函数节选：

      function Judge(Elong,JD)
      DOUBLE PRECISION Elong,JD
      double precision jump,tol
      integer year4,month4,day4,hour4,minute4,second4
      Real*8  year3,month3,day3,hour3,minute3,second3
      REAL*8  A,B,C,D,E,G
      
！根据公式，把儒略日转化为格里高利日      
      G=FLOOR((FLOOR(JD+0.5D0)-2305447.5)/36524.25)
      A=FLOOR(JD+0.5D0)+10+G-FLOOR(G/4)
      B=A+1524
      C=FLOOR((B-122.1)/365.25)
      D=FLOOR(365.25*C)
      E=FLOOR((B-D)/30.6001D0)
      day3=B-D-FLOOR(30.6001D0*E)+(JD+0.5D0)-FLOOR(JD+0.5D0)
      day4=FLOOR(day3)
      month3=REAL(E)-1D0-12D0*FLOOR(E/14D0)
      month4=FLOOR(month3)
      year3=REAL(C)-4715D0-FLOOR((7D0+month3)/10D0)
      year4=FLOOR(year3)
      hour3=(day3-FLOOR(day3))*24
      hour4=FLOOR(hour3)
      minute3=((day3-FLOOR(day3))*24-FLOOR(hour3))*60
      minute4=FLOOR(minute3)
      second3=FLOOR((minute3-FLOOR(minute3))*60)
      second4=FLOOR(second3)
      
！tol为判断精度，精确度为1毫角秒；如果已经达到一个节气，则直接将日期+14.1，跳过那些一定不是某节气的日子。下面只是2个分支，还有23个，略   
      tol=1.0/3600.0/2.0
      jump=14.1
   
      IF(ABS(ELONG-0D0).LT.tol)THEN 
          WRITE(*,*) "春分",year4,"年",month4,"月",day4,"日",
     :hour4,"时",minute4,"分",second4,"秒"
          JD=JD+jump
      ELSE IF(ABS(ELONG-15.0).LT.tol)THEN
          WRITE(*,*) "清明" ,year4,"年",month4,"月",day4,"日",
     :hour4,"时",minute4,"分",second4,"秒"
          JD=JD+jump

后记

感谢队友们的支持与鼓励，我们小组最终完成了这项任务。
这样做完的结果与网上公布的2020年二十四节气相差七分钟左右，我试了一下，如果去掉光行差的修正，差距缩小为一分钟之内。道理说不清（摊手）我不是这专业的也不想管了。
这是我第一次用Fortran编程，也应该是我最后一次用Fortran，格式不好看还请见谅。
里面涉及到的子程序（光行时修正、直角坐标求黄经、判断节气）函数功能极简单，不赘述。（光行时应迭代求出修正后的时间，直接减去八分钟误差也没大到哪儿去）
最后放两张结果图：
（1）考虑光行差：

（2）不考虑光行差：