二十四節気を計算してみた

ふと　二十四節気クロックを作ってみたいと　思いました。暦は季節の移り変わりに即したものがよかろうという考えからですが。

節気計算の式はネットにある様なので容易かと思ったのですが、結構厄介でした。参考にしたのは、下記

にあるC++コードからお借りしましたが、これはこれで略式の様で今年の立春 を計算したところ　本当の2/3　23:59 立春時刻　にはギリギリで誤差が出てしまって2/4が立春になってしまいました。。。　これと同様の式を用いたと思われる節気サイトがありますが、同じく立春が2/4 00:00 に計算されているサイトがある様でした。

いまいち精確な計算法が分からず、、暫定的ですがこの式の結果に補正を加えることでどうにか節気としては正確に求められる様にしました。具体的にはこの式から求めたψに角度で 0.0072度程度を加えただけというアバウトな方法です。

次に、節気クロックでは　「立春5日、雨水前11日」　みたいな表示をしたいのですが、節気と節気間の日数は15日だったり16日だったりまちまち（当たり前ですね）　「今日の節気は何」　は分かっても　何日目なのか、黄経の角度だけから計算しても誤差が発生してしまい、立春の３日目が二日続いたり変な事も発生。。。

結局今日の黄経から一日ずつ遡って節気から何日目かを知り、一日づつ進めて次の節気日まで勘定するという方法で対処しました。

以下サンプルプログラムです。のちほどM5Stackに実装したいと思います。

sekkiday.c

/*

   Sekki Day Calc

  cc sekkiday.c -lm -o sekkiday

  2021/2/10   calc with julian today

*/

#include <stdio.h>
#include <math.h>
#include <stdbool.h>

#define deg2rad(deg)  ((deg) / 180.0 * 3.14159265358979323846)

#define  daydeg (360/365.24219)

#define false 0
#define true  1

double julian(int nYear, int nMonth, double dDay);

double getlamda(int yr , int month , double aday);

// new funcs
typedef const int cint;
typedef const double cdouble;

double JC(double);
double CLongit( cdouble aJC);

int yr=2000;
int mn=4;
int dy=1;
double ddy=1.0;
double ddy1=1.0;

int sekki;
int sday;

int sekki2;
int sday2;

double dsday;

//double  hosei = 0.386874;     //  from eco.mtk.nao.ac.jp
double  hosei = 0.0072;
//double  hosei = 0.102753;


int main()
{
   int i,n;
   int hoseideg;
   double psi,psi0,psi1,psit,psin;
   double jday,jday1;
   double jc;

   yr = 2021;    // Year
   mn = 3;       // Month
   dy = 10;      //  Day



       ddy  = (double)dy-(9.0/24.0) + 0.99930;  // 23:59
       ddy  = (double)dy-(9.0/24.0) + 0.9999930;  // 23:59:58


       jday= julian(yr,mn,ddy);   // today julian day


       jc =  JC(jday);
       psit  = CLongit(jc)+hosei;

       sekki = (int)(psit / 15);
       sekki2 = sekki+1;
       psi0 = sekki*15;          //  today's sekki degree
       psin = sekki2*15;         // next sekki degree

      sday=1;
      while(sday<=16){

         jday1=jday-1.0*sday;
         jc =  JC(jday1);
         psi1  = CLongit(jc)+hosei;
        if(sekki==0 && psi1 >= 345){
          psi1 = psi1-360.0;
         }

         printf("%d ,  %f ,  %f , %f  \n",sday,psi1,psi0,psit);

         if(psi1 < psi0){
            break;
          }
         sday++;
      }
      printf("%d/%d : sekki=%d , sday=%d \n",mn,dy,sekki,sday);

      sday2=1;

    while(sday2<=16){
       jday1=jday+1.0*sday2;
       jc=JC(jday1);
       psi1 = CLongit(jc)+hosei;
       if(sekki2 == 24 && psi1 <= 15){
       psi1 = psi1+360;
       }
       printf("%d ,  %f ,  %f , %f  \n",sday2,psi1,psin,psit);
       if(psi1 > psin){
          break;
       }
     sday2++;
     }

     printf("%d/%d : sekki2=%d , sday2=%d \n",mn,dy,sekki2,sday2);
  return sekki;
}


double julian(int nYear, int nMonth, double dDay)
{
    double da, db;
    double dJD;

    if (nMonth < 3){
        nMonth += 12;
        nYear--;
    }

    // Calc. body
    da = floor(nYear/100.0);
    db = 2.0 - da + floor(da / 4.0);
    dJD = floor(365.25 * nYear) + floor(30.6001 * (nMonth + 1)) + dDay + db +

1720994.5;
//  dJD = floor(365.25 * nYear) + floor(30.6001 * (nMonth + 1)) + dDay + db +

1721088.5;

    if (dJD < 2299160.5){
        dJD -= db;
    }
    return (dJD);
}

/*
          psi :  psi - koukei

*/
#define _rpd   (3.14159 / 180.0)
typedef struct {
bool T;
double P, Q, R;
} SCoef;

#define PCOEF_MAX (18)

const SCoef PCoef[PCOEF_MAX] = {
{ false, 1.9147, 35999.05, 267.52, }, // 1
{ false, 0.0200, 71998.10, 265.10, }, // 2
{ false, 0.0020, 32964.00, 158.00, }, // 3
{ false, 0.0018, 19.00, 159.00, }, // 4
{ false, 0.0018, 445267.00, 208.00, }, // 5
{ false, 0.0015, 45038.00, 254.00, }, // 6
{ false, 0.0013, 22519.00, 352.00, }, // 7
{ false, 0.0007, 65929.00, 45.00, }, // 8
{ false, 0.0007, 3035.00, 110.00, }, // 9
{ false, 0.0007, 9038.00, 64.00, }, // 10
{ false, 0.0006, 33718.00, 316.00, }, // 11
{ false, 0.0005, 155.00, 118.00, }, // 12
{ false, 0.0005, 2281.00, 221.00, }, // 13
{ false, 0.0004, 29930.00, 48.00, }, // 14
{ false, 0.0004, 31557.00, 161.00, }, // 15
{ true, -0.0048, 35999.00, 268.00, }, // 16*
{ false, 0.0048, 1934.00, 145.00, }, // 17
{ false, -0.0004, 72002.00, 111.00, }, // 18
};
//

double CLongit( cdouble aJC)
{

//
//    aJC : julian year
//
        bool IsAP = true;

        double T = aJC;
        double psi = 0.0, qtr, p;
        int icmax = (IsAP ? PCOEF_MAX : PCOEF_MAX - 2);
        for ( int ic = icmax; ic >= 1; ic-- ) {
                qtr = fmod( PCoef[ic - 1].Q * T + PCoef[ic - 1].R, 360.0 ) * _rpd;
                p = PCoef[ic - 1].P;

                if ( PCoef[ic - 1].T ) p *= T;
                psi += p * cos( qtr );
        }

        if ( IsAP ) psi -= 0.00569;
        psi += 36000.7695 * T + 280.4602;
        psi = fmod(psi,360.0);
        return ( psi );
}

/*
   julian day to julian seiki;

*/
double JC(double jday)
{
    double jc;
    jc = (jday - 2451545.0)/36525;
   return jc;
}