1568-1644 Solar eclipses, make 1617(8) and 1620(12) Lunar eclipses the same as Chinese records , C++ (gcc)

1568-1644 Solar eclipses, make 1617(8) and 1620(12) Lunar eclipses the same as Chinese records

// C++ Solar eclipse program with some Lunar eclipse(1568-1644)=>adjusted MinSun(printed) with Zhuanzhong_C = 27.554368
// Author: Mr Lam/ Yee Din(林貽典) email : yeedlam@yahoo.com
// Date : Nov. 4, 2014. Copyright reserved, please contact the author if you want to use it.
 
#include <iostream>
using namespace std;
#include <cmath>
 
const int start_year = 1568; // enter start_year
const int end_year = 1644; // enter end_year
 
const double Tropical_year = 365.2425;
const int Base_year = 1281;
//===================================================================
// 4 Paramters for adjustments
const double Qiying_P = 55.06, Runying_P = 20.2050;
const double Jiaoying_P = 26.0388-0.10243;   
const double Zhuanying_P = 13.0205-0.785192;
//-------------------------------
// Constants for Moon to determine solar eclipse purpose
const double Shuoshi_C = 29.530593, Zhuanzhong_C = 27.554368;      // 27.5546;
const double Jiao_eclipse_C = 2.318369, Jiaozhong_C = 27.212224;
const double half_Zhuanzhong = Zhuanzhong_C/2;
//-------------------------------
// Constant for Moon positions : Lead or Lag
const double MoonRef_C = 1.976225;                                 //  1.975993;
const double MoonMeanSpeed = 13.36878311;                          // 13.36875;
const double MinSun_drift=0.023452;        //  MinSun=MinSun_drift*(calculated_year-Base_year);
//  RefSun_ = half_year-LeapDay;
//  RefSun=RefSun_- MinSun;               //  MinSun added back to LunStep 18 as: LEclNetSun_=LEclNetSun_+ MinSun
//===================================================================
 
typedef struct{
int year; int month; int day;
double EclipseFrac; double EclipseBegin; double MaxPhase; double EclipseEnd;
double Winter_solitice;
int SunStatus; int MoonStatus;
double observed_time; double HalfDayFan;
double ReqSun; double ReqMoon;
double half_dayoff;
 
double MoonSp;
double noon_offset;
 
double ConjuncT24;
double SunX; double MoonX; double SunMoonXS;
} EclipseTime;
 
EclipseTime Shoushi[100];
 
 EclipseTime Theo[]
= {
   {1568,  9, 21,     0.109955,  11.1367,   11.7366,   12.3428},  // from eclipse.org.uk
   {1572,  6, 10,     0.916,      6.66,      7.93,      9.93},
   {1575,  5, 10,     0.701,     13.05,     14.29,     15.45},
   {1582,  6, 20,     0.723024,  12.2055,   13.6179,   14.9154},  // from eclipse.org.uk
   {1583, 12, 14,     0.825,     11.23,     13.08,     14.77},
   {1587, 10,  2,     0.735831,  11.3223,   12.6436,   13.9343},  // from eclipse.org.uk
   {1589,  2, 15,     0.301344,   6.4499,    7.1961,    7.99914}, // from eclipse.org.uk
   {1590,  7, 31,     0.430839,  16.0803,   17.1059,   18.0372},  // from eclipse.org.uk
   {1594,  5, 20,     0.319,      9.83,     10.84,     11.9},
   {1596,  9, 22,     0.949,      9.87,     11.17,     12.5},
 
   {1603,  5, 11,     0.903,      7.4,       8.59,      9.93},
   {1604,  4, 29,     0.078077,  15.8827,   16.4079,   16.907},   // from eclipse.org.uk
   {1607,  2, 26,     999,       17.32,     18.26,     999},      // no data from eclipse.org.uk
   {1609, 12, 26,     0.044430,  16.197,    16.594,    16.9753},  // from eclipse.org.uk but no apparent noon
   {1610, 12, 15,     0.607,     14.61,     16.07,     17.32},
   {1612, 12, 99,     999,         999,      999,      999},      // no data from eclipse.org.uk
   {1615,  3, 29,     0.889903,  15.1624,   16.4118,   17.5463},  // from eclipse.org.uk
   {1617,  8,  1,     0.470019,  19.07,     19.9041,   20.6909},  // eclipse_fraction, eclipse_maximum, eclipse_end from eclipse.org.uk
   {1621,  5, 21,     0.711,     17.27,     18.23,     19.19},
   {1629,  6, 21,     0.166,     10.98,     11.71,     12.45},
 
   {1631, 10, 25,     0.108,     12.34,     12.91,     13.48},
   {1634,  3, 29,     0.81964,    8.77022,   9.94358,  11.1899 },  // from eclipse.org.uk
   {1636,  8,  1,     0.632763,   8.37958,   9.62653,  10.9643},   // from eclipse.org.uk
   {1637,  1, 26,     0.199,     12.71,     13.68,     14.58},
   {1638,  1, 15,     0.656,     11.65,     13.32,     14.84},
   {1641, 11,  3,     0.827,     12.65,     14.04,     15.33},
   {1643,  3, 20,     0.53,       7.83,      8.85,      9.94},
   {1644,  9,  1,     0.303,     10.93,     12.14,     13.34},
 
   {99999, 0,  0,     0,        0,        0,        0}
  };
 
 EclipseTime Observe[]
= {
   {1568,  9, 21,     0.109955,  11.1367,   11.7366,   12.3428},  // from eclipse.org.uk
   {1572,  6, 10,     0.916,      6.66,      7.93,      9.93},
   {1575,  5, 10,     0.701,     13.05,     14.29,     15.45},
   {1582,  6, 20,     0.723024,  12.2055,   13.6179,   14.9154},  // from eclipse.org.uk
   {1583, 12, 14,     0.825,     11.23,     13.08,     14.77},
   {1587, 10,  2,     0.735831,  11.3223,   12.6436,   13.9343},  // from eclipse.org.uk
   {1589,  2, 15,     0.301344,   6.4499,    7.1961,    7.99914}, // from eclipse.org.uk
   {1590,  7, 31,     0.430839,  16.0803,   17.1059,   18.0372},  // from eclipse.org.uk
   {1594,  5, 20,     0.319,      9.83,     10.84,     11.9},
   {1596,  9, 22,     0.949,      9.87,     11.17,     12.5},
 
   {1603,  5, 11,     0.903,      7.4,       8.59,      9.93},
   {1604,  4, 29,     0.078077,  15.8827,   16.4079,   16.907},   // from eclipse.org.uk
   {1607,  2, 26,     999,       17.32,     18.26,     999},      // no data from eclipse.org.uk
   {1609, 12, 26,     0.044430,  16.197,    16.594,    16.9753},  // from eclipse.org.uk but no apparent noon
   {1610, 12, 15,     0.607,     14.61,     16.07,     17.32},
   {1612, 12, 99,     999,         999,      999,      999},      // no data from eclipse.org.uk
   {1615,  3, 29,     0.889903,  15.1624,   16.4118,   17.5463},  // from eclipse.org.uk
   {1617,  8,  1,     0.470019,  19.07,     19.9041,   20.6909},  // eclipse_fraction, eclipse_maximum, eclipse_end from eclipse.org.uk
   {1621,  5, 21,     0.8,       17.24,     9999,      9999},
   {1629,  6, 21,     0.166,     10.98,     11.71,     12.45},
 
   {1631, 10, 25,     0.108,     12.34,     12.91,     13.48},
   {1634,  3, 29,     0.81964,    8.77022,   9.94358,  11.1899 },  // from eclipse.org.uk
   {1636,  8,  1,     0.632763,   8.37958,   9.62653,  10.9643},   // from eclipse.org.uk
   {1637,  1, 26,     0.199,     12.71,     13.68,     14.58},
   {1638,  1, 15,     0.656,     11.65,     13.32,     14.84},
   {1641, 11,  3,     0.827,     12.65,     14.04,     15.33},
   {1643,  3, 20,     0.53,       7.83,      8.85,      9.94},
   {1644,  9,  1,     0.303,     10.93,     12.14,     13.34},
 
   {99999, 0,  0,     0,        0,        0,        0}
  };
 
enum Status{LagStart, LagEnd, LeadStart, LeadEnd};
struct Status_type{char *sun, *moon;};
 
class Ref_data
{
  private:
    Status SunStatus, MoonStatus;
 
  protected:
    int calculated_year, MonthNum;
    double SumDay, Winter_solitice, LeapDay;
    double RefEclipse, RefSun, RefMoonDay, RefMoon;
    double ReqEclipse, ReqSun, ReqMoonDay, ReqMoon;
    double UseSun, SunAdd, SunSum, SunX;
    double MoonSum, SoD;
    double moon_x12_2;
    int FrontEclipseSign;
 
  public:
    Ref_data(int m,int n){calculated_year = m;  MonthNum = n;}
    ~Ref_data(){ }
 
void check_SunStatus(double half_year,double ReqSun, double& UseSun, int& SunStatus)
    {
      SunStatus=LagStart;
      UseSun = 0;
 
      if (ReqSun < half_year)
      {
        SunStatus = LagStart;
        UseSun = ReqSun;
        if (ReqSun > 93.712025){UseSun = half_year - UseSun; SunStatus = LagEnd;}
      }
//  推初末限:　視盈曆在盈初縮末限已下，縮曆在縮初盈末限已下，各為初。已上用減半歲周為末。
     if (ReqSun > half_year)
     {
        SunStatus = LeadStart;
        UseSun = ReqSun - half_year ;
        if (UseSun > 88.909225){UseSun = half_year - UseSun; SunStatus = LeadEnd;}
     }
   }//end check_SunStatus
 
void show_SunStatus(int SunStatus)
    {
    cout << "SunStatus = ";
    switch (SunStatus)
    {
      case LagStart:cout << "SunLagStart\n"; break;
      case LagEnd:cout << "SunLagEnd\n"; break;
      case LeadStart:cout << "SunLeadStart"; break;
      case LeadEnd:cout << "SunLeadEnd"; break;
      default: cout << "no SunStatus"; break;
    }
    cout << endl;
    }//end show_SunStatus
 
void get_SunX(int SunStatus,double UseSun,int IUseSun,double& SunAdd,double& SunSum,double& SunX)
    {
    SunAdd=SunSum=SunX=0;
    if (fmod(SunStatus, 3.0)==0) // SunStaus=LagStart(0) & LeadEnd(3)
    {
    SunAdd = 484.8473-(4.4362 + 0.0162*(IUseSun - 1)/2)*IUseSun ;
    SunSum = (487.06-(2.21 + 0.0027*IUseSun )*IUseSun )*IUseSun ;
    }
    if (fmod(SunStatus, 3.0) > 0) // SunStatus=LagEnd(1) & LeadStart(2)
    {
    SunAdd = 510.8569-(4.9386 + 0.0186*(IUseSun - 1)/2)*IUseSun ;
    SunSum =(513.32-(2.46 + 0.0031*IUseSun)*IUseSun)*IUseSun;
    }
    SunX =((UseSun - IUseSun)*SunAdd + SunSum)/10000.0;
    }// end get_SunX
//--------------------------------------------------------------
 void get_LSunX(int LSunStatus,double LUseSun, double& LSunX)
 {
    double LSunSum=0;
 
    if (fmod(LSunStatus, 3.0) == 0){// for LSunStaus = LagStart(0) & LeadEnd(3)
                                    LSunSum = (487.06-(2.21 + 0.0027*LUseSun)*LUseSun)*LUseSun ;
                                   }
    if (fmod(LSunStatus, 3.0) > 0){// for LSunStatus = LagEnd(1) & LeadStart(2)
                                    LSunSum =(513.32-(2.46 + 0.0031*LUseSun)*LUseSun)*LUseSun;
                                  }
    LSunX=LSunSum/10000.0;
 }// end get_LSunX
//--------------------------------------------------------------
void get_UseMoon(double ReqMoon,double& UseMoon,double& SoD,double& moon_x12_2,double& MoonSolar,
                 int& MoonStatus, double& MoonCh, double& MoonAdd)
    {
    double moon_84_check;
    if (ReqMoon > half_Zhuanzhong){UseMoon = ReqMoon - half_Zhuanzhong; MoonStatus = LagStart;}
    if (ReqMoon < half_Zhuanzhong){UseMoon = ReqMoon; MoonStatus = LeadStart;}
 
//  推遲疾限:　各置遲疾曆，以日轉限乘之，即得限數。(日轉限一十二限二十。)
    moon_x12_2=UseMoon*12.2;
    SoD=int(moon_x12_2)*0.082008;
 
    moon_84_check=moon_x12_2;
    if (moon_x12_2 > 84){
                         if (MoonStatus==LeadStart){MoonStatus=LeadEnd;}
                         if (MoonStatus==LagStart){MoonStatus=LagEnd;}
                         moon_84_check = 167 - moon_84_check;
                        }
    MoonCh = moon_84_check;
    MoonSolar = UseMoon - SoD;
    if (MoonSolar < 0){
//  求遲疾差:　置遲疾曆，以立成日率減之，(如不及減，則退一位。)餘以其下損益分乘之為實，八百二十分為法除之，得數以加其下遲疾積，即遲疾差。
                        SoD=SoD-0.082008;
                        MoonSolar=UseMoon-SoD;
                        if (fmod(MoonStatus, 2.0)==0){MoonCh=int(MoonCh);}
                        if (fmod(MoonStatus, 2.0) > 0){MoonCh=167-int(moon_x12_2);}
                       }
    MoonAdd = 11.081575-(0.05815 + 0.00195*(MoonCh - 1)/2)*MoonCh;
    }// end get_UseMoon
 
void get_MoonX(int MoonStatus, double MoonCh, double MoonAdd,double MoonSolar,
               double& MoonSum, double& MoonX, double& MoonSp)
    {
    MoonSum=0, MoonX=0, MoonSp=0;
    switch(MoonStatus)
     {
      case LagStart:
           MoonSp=0.9855+(0.05815+0.00195*(MoonCh-1)/2)*MoonCh/100;
           MoonSum=(0.1111-(0.000281+0.00000325*MoonCh)*MoonCh)*MoonCh;
           MoonX=MoonSum+MoonSolar*MoonAdd/820;
           break;
 
      case LagEnd:
           MoonSp=1.2071-(0.05815+0.00195*(MoonCh-1)/2)*MoonCh/100;
           MoonCh=MoonCh+1;
           MoonSum=(0.1111-(0.000281+0.00000325*MoonCh)*MoonCh)*MoonCh;
           MoonX=MoonSum-MoonSolar*MoonAdd/820;
           break;
 
      case LeadStart:
           MoonSp=1.2071-(0.05815+0.00195*(MoonCh-1)/2)*MoonCh/100;
           MoonSum=(0.1111-(0.000281+0.00000325*MoonCh)*MoonCh)*MoonCh;
           MoonX=MoonSum+MoonSolar*MoonAdd/820;
           break;
 
      case LeadEnd:
           MoonSp=0.9855+(0.05815+0.00195*(MoonCh-1)/2)*MoonCh/100;
           MoonCh=MoonCh+1;
           MoonSum=(0.1111-(0.000281+0.00000325*MoonCh)*MoonCh)*MoonCh;
           MoonX=MoonSum-MoonSolar*MoonAdd/820;
           break;
 
      default:
           break;
     }// end switch(MoonStatus)
    }//  end get_MoonX
 //----------------------------------------------------
 
void get_LMoonX(int LMoonStatus,double LMoonCh_more,double& LMoonX,double& LMoonSp,double& LMoonCh_)
  {
    LMoonX=0,LMoonSp=0;
    if (fmod(LMoonStatus,2) > 0){LMoonCh_more = LMoonCh_more + 1;}
 
    LMoonSp=0.9855 +(0.05815 + 0.00195*(LMoonCh_more-1)/2)*LMoonCh_more/100;
    if(fmod(LMoonStatus, 3) > 0){LMoonSp=1.2071-(0.05815 + 0.00195*(LMoonCh_more-1)/2)*LMoonCh_more/100;}
 
    LMoonX =(0.1111-(0.000281 + 0.00000325*LMoonCh_more)*LMoonCh_more)*LMoonCh_more;
    LMoonCh_=LMoonCh_more;
  } // end get_LMoonX
//-----------------------------------------------------------
 
void get_SunMoonT(int SunStatus,int MoonStatus,double MoonX,double SunX,double MoonSp,
                  double& SunXS, double& MoonXS, double& SunMoonT)
//  推加減差:　視經朔弦望下所得盈縮差、遲疾差，以盈遇遲、縮遇疾為同相併，盈遇疾、縮遇遲為異相較，各以八百二十分乘之為實，
    {
    SunXS=SunX;   if (SunStatus < LeadStart){SunXS = -SunX;}
    MoonXS=MoonX; if( MoonStatus > LagEnd){MoonXS = -MoonX;}
 
    SunMoonT =(MoonXS+SunXS)*0.082/MoonSp;
    } // end get_SunMoonT
//--------------------------------------------------------
 
void print_SolarEclipseTable(int p)
    {
        int index=p;
    printf("Reference:\n");
    printf("(1)Accuracy of Solar Eclipse Calculations, 1644-1785, with the Chinese Shoushi Calandar(Table 2)\n");
    printf("http://articles.adsabs.harvard.edu//full/1999JHA....30..161L/0000161.000.html\n\n");
 
    printf("(2)Accuracy of Solar Eclipse Observations made by Jesuit Astronomers in China\n");
    printf("http://adsabs.harvard.edu/full/1995JHA....26..227S\n");
    printf("Observation data from Table 2\n");
    printf("Observation error data from Table 3\n");
    printf("Theorectical=Observation-Observation error\n\n");
 
    printf("\t Shoushi(Observation)\t\tRemark : 999 or 9999 = no data\n");
    printf("\t-)Theorectical(Obsevation-Shoushi)\n");
    printf("\t--------------\n");
    printf("Shoushi-Theorectical(Observation-Theorectical)\n");
    printf("==================================================================================================================\n");
    printf("N  year(month)\t  EclipseFraction\tEclipseBegin\t\tMaxPhase\t\tEclipseEnd\n");
    printf("      month/day\n");
    printf("==================================================================================================================\n");
    int record;
    for (record=0; record < index; record ++)
    {
      if(Shoushi[record].year==0){Shoushi[record].year=Observe[record].year;}
      printf("%2d %4d(%1d)",record,Shoushi[record].year, Shoushi[record].month);
      printf(" \t   %5.4f", Shoushi[record].EclipseFrac);
      printf("(%5.3f)\t",Observe[record].EclipseFrac);
      printf("%5.4f", Shoushi[record].EclipseBegin);
      printf("(%5.3f)   \t", Observe[record].EclipseBegin);
      printf("%5.4f", Shoushi[record].MaxPhase);
      printf("(%5.3f)   \t", Observe[record].MaxPhase);
      printf("%5.4f", Shoushi[record].EclipseEnd);
      printf("(%5.3f)", Observe[record].EclipseEnd);
      printf("\n");
      printf("    %4d/%1d",Theo[record].month, Theo[record].day);
 
      printf("\t-) ");
      printf("%5.3f", Theo[record].EclipseFrac);
      printf("(%5.4f)   \t", Observe[record].EclipseFrac-Shoushi[record].EclipseFrac);
 
      printf("%5.3f", Theo[record].EclipseBegin);
      printf("(%5.4f)   \t",Observe[record].EclipseBegin-Shoushi[record].EclipseBegin);
 
      printf("%5.3f", Theo[record].MaxPhase);
      printf("(%5.4f)   \t",Observe[record].MaxPhase-Shoushi[record].MaxPhase);
      printf("%5.3f", Theo[record].EclipseEnd);
      printf("(%5.4f)",Observe[record].EclipseEnd-Shoushi[record].EclipseEnd);
      printf("\n");
      printf("\t\t----------------------------------------------------------------------------------------\n");
      printf("\t\t   ");
      printf("%5.4f", Shoushi[record].EclipseFrac-Theo[record].EclipseFrac);
      printf("(%5.3f)   \t", Observe[record].EclipseFrac-Theo[record].EclipseFrac);
 
      printf("%5.4f",Shoushi[record].EclipseBegin-Theo[record].EclipseBegin);
      printf("(%5.3f)   \t",Observe[record].EclipseBegin-Theo[record].EclipseBegin);
 
      printf("%5.4f",Shoushi[record].MaxPhase-Theo[record].MaxPhase);
      printf("(%5.3f)   \t",Observe[record].MaxPhase-Theo[record].MaxPhase);
 
      printf("%5.4f",Shoushi[record].EclipseEnd-Theo[record].EclipseEnd);
      printf("(%5.3f)",Observe[record].EclipseEnd-Theo[record].EclipseEnd);
      printf("\n\n");
    }; // end for loop
   }  // end print_SolarEclipseTable
  };// end class Ref_data
 
void print_CalculatedRecords_1617()
{
 printf("\t\t《大明神宗顯皇帝實錄卷之五百四十七》\n");
 printf("\t\t\t\t邢雲路獻七政真數..(望交泛分:13.845663)\n");
 printf("\t\t以此法布算今歲七月十六日戊寅夜望月食,推得是月:\n");
 printf("\t\t\t望交泛分一十三日八十四刻五十六分六十三秒(13.845663)\n");
 printf("\t\t\t陰歷交前度一度二十二分零三秒九十五微(1.220395)\n");
 printf("\t\t\t月食分一十三分五十九秒六十九微(1.35969)\n");
 printf("\t\t\t定用分六刻六十八分九十五秒二十一微(0.06689521)\n");
 printf("\t\t\t既內分二刻八十六分六十五秒\n");
 printf("\t\t\t既外分三刻八十一分七十九秒\n");
 printf("\t\t\t初虧分六刻四十三分二十秒八十九微(0.064[3]2089,應修正為:復圓分-2*定用分=0.19801131-0.06689521*2=0.064[2]2089)\n");
 printf("\t\t\t食既分一十刻二十五分零四十微\n");
 printf("\t\t\t食甚分一十三刻一十一分六十六秒一十微(0.1311[6]610,應修正為:復圓分-定用分=0.19801131-0.06689521=0.1311[1]610)\n");
 printf("\t\t\t生光分一十五刻九十八分三十一秒八十微\n");
 printf("\t\t\t復圓分一十九刻八十分一十一秒三十一微(0.19801131)\n\n");
 
 printf("\t\t\t\t   初虧丑初二刻(1.48h<=>0.064[2]2089*24=1.54130136,仍然算是丑初二刻,所以該記錄是可靠的)\n");
 printf("\t\t\t食既丑正一刻,食甚寅初二刻(3.48h? <=>0.1311[1]610*24=3.1467864h=>應修正為寅初一刻)\n");
 printf("\t\t\t生光寅初三刻,復明寅正三刻(4.72h<=>0.19801131*24=4.75227144h=寅正三刻,所以該記錄亦是可靠的)\n");
}// end  print_CalculatedRecords_1617
 
void print_CalculatedRecords_1620()
{
 printf("\t\t\t《大明熹宗哲皇帝實錄卷之七》\n");
 printf("\t\t\t\t按新法推泰昌元年庚申歲十一月十六日己丑夜望月食:\n");
 printf("\t\t\t\t\t初虧漏下二百七十三籌五十三分,計九十七刻二十四分(0.9724)\n");
 printf("\t\t\t\t\t食既漏下二百九十四籌六十分,計八十八分\n");
 printf("\t\t\t\t\t生光漏下三百三十二籌四十分,計七刻四十分\n");
 printf("\t\t\t\t\t復圓漏下三百五十三籌四十八分,計一十一刻(0.11)\n");
 printf("\t\t\t\t各以發斂求之,得:\n");
 printf("\t\t\t\t\t初虧夜子初一刻(23.24h)\n");
 printf("\t\t\t\t\t食既子正初刻\n");
 printf("\t\t\t\t\t食甚子正四刻(0.96h)\n");
 printf("\t\t\t\t\t生光丑初三刻\n");
 printf("\t\t\t\t\t復圓丑正二刻(2.48h,但見有改自三刻的痕跡)\n");
 printf("\t\t\t\t\t月食一十四分九十九秒(1.499)\n");
 printf("\t\t\t\t\t食甚月離黃道畢宿一十四度三十分一十秒(14.3010)\n");
 printf("\t\t\t\t而授時則推:\n");
 printf("\t\t\t\t\t初虧子正一刻(0.24h)\n");
 printf("\t\t\t\t\t食既丑初一刻\n");
 printf("\t\t\t\t\t食甚丑正初刻(2.00h?)\n");
 printf("\t\t\t\t\t生光丑正三刻\n");
 printf("\t\t\t\t\t復圓寅初三刻(3.72h?)\n");
 printf("\t\t\t\t\t月食一十三分三十一秒(1.331)\n");
 printf("\t\t\t\t\t月離黃道畢宿一十四度三十四分四十四秒(14.3444),與天不合\n");
}   // end print_CalculatedRecords_1620
 
void print_CalculatedRecords_1623()
{
 printf("\t\t\t《大明熹宗哲皇帝實錄卷之七》\n");
 printf("\t\t\t\t以新法推天啟三年癸亥歲九月十五壬寅夜望月食:\n");
 printf("\t\t\t\t\t初虧申正三刻\n");
 printf("\t\t\t\t\t食甚酉初三刻(17.72?,該記錄不可靠)\n");
 printf("\t\t\t\t\t復圓戍初初刻\n");
 printf("\t\t\t\t\t月食五分三十八秒\n");
 printf("\t\t\t\t\t月未出巳食五分二十五秒,月已出見食一十三秒,不及一分不救\n");
 printf("\t\t\t\t而授時則推月食:\n");
 printf("\t\t\t\t\t初虧酉初一刻\n");
 printf("\t\t\t\t\t食甚酉正二刻(18.48h,該記錄可靠)\n");
 printf("\t\t\t\t\t復圓戌正一刻\n");
 printf("\t\t\t\t\t月食七分八秒\n");
 printf("\t\t\t\t\t月未出已食一分九十八秒,月已出見食五分一十秒,與天不合\n");
}   // end print_CalculatedRecords_1623
 
void print_FindingRecords_1621()
{
 printf("\t\t《大明熹宗哲皇帝實錄卷之七》\n");
 printf("\t\t原任陝西按察使邢雲路獻七政真數.....復以此法推天啟元年辛酉歲四月壬申朔日食:\n");
 printf("\t\t\t\t初虧申正一刻(16.24h)\n");
 printf("\t\t\t\t食甚酉初刻(17.00h)\n");
 printf("\t\t\t\t復圓酉正初刻(18.00h)\n");
 printf("\t\t\t\t日食一分八十六秒(0.186)\n");
 printf("\t\t\t\t復圓日在天未入地,食不及三分不救(以上記錄缺乏其它數據可供證驗,是否可靠,甚有疑問?)\n");
 printf("\t\t而授時則推是日食:\n");
 printf("\t\t\t\t初虧申正三刻(16.72h)\n");
 printf("\t\t\t\t食甚酉正初刻(18.00h\n");
 printf("\t\t\t\t復圓戍初刻(19.00h)\n");
 printf("\t\t\t\t日食三分九十一秒(0.391)\n");
 printf("\t\t\t\t日未入已復光三分一十一秒,日已入未復光八十秒,與天不合(以上授時數據可按計算方法證驗)\n\n");
 printf("\t\t《明熹宗悊皇帝實錄卷之十六》\n");
 printf("\t\t庚申原任陝西按察使閒住邢雲路奏:\n");
 printf("\t\t\t\t臣本年春言四月朔日食歷數,臨期策驗分刻乃有不合\n");
 printf("\t\t\t\t臣始悟歷元玄法原自無差,而其差者,乃今大統歷沿郭守敬授時歷之誤也\n");
 printf("\t\t\t\t蓋日食者,月食之也,故正歷一準於月食,而月之分數與限度定法,皆其綱領最大者\n");
 printf("\t\t\t\t授時定月食分三十分及十五分,月食限十三度五分,定限八十七分\n");
 printf("\t\t\t\t此守敬率南北日官以儀器測驗于天,自丙子至庚辰又經甲午,再驗幾二十年,自謂推步之精,可永久無弊\n");
 printf("\t\t\t\t以故大統歷沿用之,臣今所步歷亦用之,及四月朔食測驗,方見其差,乃悟以異乘同\n");
 printf("\t\t\t\t除另變其法,變月食分三十分為三十五分(30->35)\n");
 printf("\t\t\t\t         變十五分為十七分五十秒\n");
 printf("\t\t\t\t         變定限十三度五分為十三度十二分五十秒(13.05->13.1250)\n");
 printf("\t\t\t\t         變定法八十七分為七十五分(8.7->7.5)\n");
 printf("\t\t\t\t諸數既變,乃依舊歷元玄法布歷,一一方合,然後知歷元法原自無差,而其差者,大統相沿授時法誤之也,故必如臣所推:\n");
 printf("\t\t\t\t如五千七百四十為日食成限,以四千九百二十為月食成限,以三十五分為月食分,而黜郭守敬之三十分,方為正法,下禮部\n\n");
 
 printf("\t\t\t\t  LunFraction=(13.05 - abs(SexLunEclA_))/8.7\n");
 printf("\t\t\t月食食分=定限(13.05->13.1250) - 視入陰陽曆/定法(8.7->7.5)\n\n");
 
 printf("\t\t\t\tLHalfEclT=(sqrt((30-LunFraction*10)*LunFraction*10)*49.20)/(LEclM_SunSp*10000)\n");
 printf("\t\t\t\t定用分=(sqrt(三十分(30->35)-月食分)*月食分)*49.20)/(定限行度*10000)\n");                   
}   // end print_FindingRecords_1621
 
void print_LunStep_13_Extra(int k)
{
    int calculated_year=k;
    if (calculated_year==1617){
                               printf("\t\t\t\t\t\t\t\t\t\t\t邢雲路食甚分 0.13111610*24=3.1467864h\n");
                               printf("\t\t\t\t\t\t\t\t\t\t\t\t\t   理論值:3.219\n");
                              }
    if (calculated_year==1620){
          printf("   邢雲路..按新法推..初虧..計九十七刻二十四分..復圓..計一十一刻=>食甚分=0.11-(0.11+(1-0.9724))/2=0.0412");
          printf("=>食甚=0.0412*24=0.9888h=子正四刻,與記錄一致\n");
          printf("\t\t\t\t\t\t\t\t\t\t\t\t\t    理論值:1.17h \n");
                             }
    if (calculated_year==1623){
                               printf("\t\t\t\t\t\t\t\t\t\t\t  而授時則推...食甚酉正二刻(18.48h),該記錄可靠\n");
                               printf("\t\t\t\t\t\t\t\t\t\t\t\t\t   理論值:19.23h \n");
                              }
}   // print_LunStep_13_Extra
 
void print_LunStep_23_Extra(int k)
{
    int calculated_year=k;
    if (calculated_year == 1617){
                                 printf("\t\t\t\t\t邢雲路獻七政真數..以此法布算..月食分一十三分五十九秒六十九微(1.35969)\n");
                                 printf("\t\t\t\t\t\t\t\t\t\t\t 理論值:1.403\n");        
                            }
    if (calculated_year == 1620){
                                 printf("\t\t\t\t\t\t\t    授時則推..月食一十三分三十一秒(1.331)\n");
                                 printf("\t\t\t\t\t\t邢雲路.....按新法推....月食一十四分九十九秒(1.499)\n");
                                 printf("\t\t\t\t\t\t\t\t\t\t   理論值:1.6\n");
                                }
    if (calculated_year == 1623){
                                 printf("\t\t\t\t以新法推天啟三年癸亥歲九月十五壬寅夜望月食...月食五分三十八秒(0.538)\n");
                                 printf("\t\t\t\t\t\t\t\t\t\t\t理論值:0.6\n");
                                }
}   // end  print_LunStep_23_Extra
 
void print_LunStep_25_Extra(int k)
{
    int calculated_year=k;
    if (calculated_year == 1617)
    {
        printf("\t\t\t\t\t(邢雲路的初虧*24=1.54130136h=>算是丑初二刻\t   4.75227144h=寅正三刻,兩項記錄都可靠)\n");
        printf("\t\t\t\t\t\t理論值: 1.49h      \t\t   4.95h\n\n");
                  
double  SexLunEclA_1617=13.05-1.35969*8.7; // 計算出原有的 SexLunEclA_1617(不分陰陽曆)               
double  NewFraction=(13.125-1.220395)/7.5;              
                                           
double  LEclM_SunSp_1617=(sqrt((30-13.5969)*13.5969)*49.20)/(0.06689521*10000);  // 計算出原有的 LEclM_SunSp_1617
double  NewHalfT=(sqrt((35-NewFraction*10)*NewFraction*10)*49.20)/(LEclM_SunSp_1617*10000);
                  
double  NewBeginT=0.13111610-NewHalfT;
double  NewBeginT24=NewBeginT*24;                  
                  
double  NewEndT=0.13111610+NewHalfT;
double  NewEndT24=NewEndT*24;
        printf("\t\t\t\t邢雲路原有的陰陽曆 SexLunEclA_1617(不分陰陽)=13.05-1.35969*8.7=%1.5f,與記錄所指的 1.220395 基本一致\n\n",SexLunEclA_1617);
 
        printf("\t\t\t\t邢雲路變法後的食分 NewFraction=(13.125-1.220395)/7.5=%1.5f,大於理論值\n",NewFraction);
        printf("\t\t\t\t\t\t\t\t\t     理論值=1.403\n");
        printf("\t\t\t\t\t\t\t\t   邢雲路原有的食分: 1.35969\n\n");
 
        printf("\t\t  邢雲路原有的定限行度 LEclM_SunSp_1617=(sqrt((30-13.5969)*13.5969)*49.20)/(0.06689521*10000)=%1.5f\n",LEclM_SunSp_1617);                    
        printf("\t\t  邢雲路變法後改大的定用分 NewHalfT=(sqrt((35-NewFraction*10)*NewFraction*10)*49.20)/(LEclM_SunSp_1617*10000)=%1.5f",NewHalfT);
        printf("\n");
        printf("\t\t  邢雲路變法後的 NewBeginT=0.13111610-NewHalfT=%1.5f\t NewEndT=0.13111610+NewHalfT=%1.5f\n\n",NewBeginT,NewEndT);    
        printf("\t\t\t變法後的初虧=NewBeginT*24=%1.5f\t\t食甚不變=0.13111610*24=%1.5f\t\t復圓=NewEndT*24=%1.5f\n",NewBeginT24,0.13111610*24, NewEndT24);
        printf("\t\t\t\t     理論值: 初虧=1.49h\t\t\t\t  食甚=3.219h\t\t\t   復圓=4.95h\n");                   
        printf("\t\t\t\t 邢雲路原有的初虧=1.54130136h\t\t\t\t\t\t\t   復圓=4.75227144h\n");                                             
    }
    if (calculated_year == 1620)
    {
        printf("\t\t (邢雲路...按新法推..初虧=0.9724*24 = 23.3376h\t\t\t\t復圓=0.11*24=2.64h)\n");
        printf("\t\t\t\t\t理論值:初虧= 23.28h\t\t\t\t\t復圓=3.05h\n\n");
//-----------
double  SexLunEclA_1620=13.05-1.499*8.7; // 計算出原有的 SexLunEclA_1620(不分陰陽曆)                
double  NewFraction=(13.125-SexLunEclA_1620)/7.5;              
 
double  HalfT_1620=(0.11+(1-0.9724))/2;
double  LEclM_SunSp_1620=(sqrt((30-14.99)*14.99)*49.20)/(HalfT_1620*10000);  // 計算出原有的 LEclM_SunSp_1620
double  NewHalfT=(sqrt((35-NewFraction*10)*NewFraction*10)*49.20)/(LEclM_SunSp_1620*10000);
                  
double  NewBeginT=1+0.0412-NewHalfT;
double  NewBeginT24=NewBeginT*24;                  
                  
double  NewEndT=0.0412+NewHalfT;
double  NewEndT24=NewEndT*24;
        printf("\t\t\t\t邢雲路原有的陰陽曆 SexLunEclA_1620(不分陰陽)=13.05-1.499*8.7=%1.5f\n",SexLunEclA_1620);
        printf("\t\t\t\t邢雲路變法後的食分 NewFraction=(13.125-SexLunEclA_1620)/7.5=%1.5f\n\n",NewFraction);
 
        printf("\t\t\t\t邢雲路原有的定用分 HalfT_1620=(0.11+(1-0.9724))/2=%1.5f\n",HalfT_1620);
        printf("\t\t\t\t邢雲路原有的定限行度 LEclM_SunSp_1620=(sqrt((30-14.99)*14.99)*49.20)/(HalfT_1620*10000)=%1.5f\n",LEclM_SunSp_1620);                    
        printf("\t\t  邢雲路變法後的定用分 NewHalfT=(sqrt((35-NewFraction*10)*NewFraction*10)*49.20)/(LEclM_SunSp_1620*10000)=%1.5f\n",NewHalfT);
        printf("\t\t\t\t邢雲路變法後的 NewBeginT=1+0.0412-NewHalfT=%1.5f\tNewEndT=0.0412+NewHalfT=%1.5f\n\n",NewBeginT,NewEndT);
 
        printf("\t\t邢雲路變法後的初虧=NewBeginT*24=%1.5f\t\t食甚不變=0.0412*24=%1.5f\t     復圓=NewEndT*24=%1.5f\n",NewBeginT24,0.0412*24,NewEndT24);
        printf("\t\t\t\t  理論值:初虧= 23.28h\t\t\t     食甚=1.17h\t\t\t\t復圓=3.05h\n");
        printf("\t\t    邢雲路原有的初虧=0.9724*24=%1.5f\t\t\t\t\t\t\t復圓=0.11*24=%1.5fh\n",0.9724*24,0.11*24); 
    }
    if (calculated_year == 1623)
    {
        printf("    以新法推天啟三年癸亥歲九月十五壬寅夜望月食[初虧申正三刻](16.72)不可靠, 因為這應該是授時的計算值\n");
        printf("\t\t\t     而授時則推月食[初虧酉初一刻](17.24)不可靠, 因為與計算值不符合,很可能是新法推的初虧\n");
        printf("\t\t\t 以新法推天啟三年..[初虧申正三刻](16.72h ?)\t ....  復圓[戍初初刻](19h ?)\n");
        printf("\t\t\t\t\t\t理論值:17.7h\t\t\t\t     20.79h\n\n");
    }
}   //  end print_LunStep_25_Extra
 //---------------------------------------------------------------------
    int main()
    {
    cout << "\tSolar eclipse time is in 24 hours\n"
         << "\tBeijing local noon time is 12 as reference\n"
         << "\tmonth means total months from winter solitice\n"
         << "\tReqMDay=ReqMoonDay" << endl;
    cout << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt"
         << "\tt" << "\tt" << "\tt" << endl
         << " N year(month)\n";
 
    int index = 0; int calculated_year, MonthNum;
    for (calculated_year = start_year; calculated_year - 1< end_year; calculated_year ++)
    {
    for (MonthNum = 0; MonthNum < 14; MonthNum ++)
    {
//--------------
//     Step 1:
 
       Ref_data year(calculated_year,MonthNum);
double YearDays=Tropical_year-(calculated_year-Base_year)*0.0001/100;  //0.0001/100
double half_year=YearDays/2;
double SumDay=(calculated_year-Base_year)*YearDays;    
double Winter_solitice=fmod((SumDay + Qiying_P + int(abs(SumDay))*60), 60);
 
//     推天正閏餘:　置中積，加閏應，滿朔策去之，至不滿之數，為天正閏餘。
double LeapDay = fmod((SumDay + Runying_P + int(abs(SumDay))*Shuoshi_C), Shuoshi_C);
 
//--------------
//     Step 2:
 
//     推天正盈縮:　置半歲周，內減其年閏餘全分，餘為所求天正縮曆。
double RefSun_ = half_year-LeapDay;
double MinSun=MinSun_drift*(calculated_year-Base_year);
double RefSun=RefSun_- MinSun;
        
//     推天正經朔:　置冬至，減閏餘，遇不及減，加紀法減之，為天正經朔。
double RefMoonDay = Winter_solitice-LeapDay;
       if (RefMoonDay < 0){RefMoonDay = RefMoonDay + 60;}
 
//     推天正遲疾:　置中積，加轉應，減去其年閏餘全分，餘滿轉終去之，即天正入轉。        
double RefMoon = fmod((SumDay + Zhuanying_P - LeapDay + int(abs(SumDay))*Zhuanzhong_C), Zhuanzhong_C);
 
//     推天正入交:　置中積，減閏餘，加交應，滿交終去之，即天正入交汎日。        
double RefEclipse = fmod((SumDay + Jiaoying_P - LeapDay + int(abs(SumDay))*Jiaozhong_C), Jiaozhong_C);
 
//--------------
//     Step 3:
 
//     推各月盈縮曆:　置天正縮曆，加二朔策，去半歲周，即得正月經朔下盈曆。
double ReqSun = RefSun + Shuoshi_C*MonthNum;
       ReqSun = fmod(ReqSun, YearDays);            
             
//     推各月經朔:　置天正經朔，加二朔策，滿紀法去之，即得正月經朔。
double ReqMoonDay = fmod((RefMoonDay + Shuoshi_C*MonthNum), 60);
 
//     推各月遲疾曆:　置天正經朔遲疾曆，加二轉差，得正月經朔下遲疾曆。
double ReqMoon = fmod(( RefMoon + MoonRef_C*MonthNum), Zhuanzhong_C);
 
//     推各月入交:　置天正經朔入交汎日加二交差，得正月經朔下入交汎日。
double ReqEclipse = fmod((RefEclipse + Jiao_eclipse_C*MonthNum), Jiaozhong_C);
 
//--------------
//     Step 4:
//     推交常度　置有食之朔入交汎分，以月平行度乘之，即得。
double EclAngle = ReqEclipse * MoonMeanSpeed;
 
//--------------
//     Step 5:
//     如滿半歲周去之交縮，滿半歲周又去之即復交盈。(quoted from Step 3 for 推各月盈縮曆:)
double ChSun = ReqSun;
       if (ReqSun > half_year){ChSun = ReqSun - half_year;}
 
//--------------
//     Step 6:
 
//     推初末限:　視盈曆在盈初縮末限已下，縮曆在縮初盈末限已下，各為初。已上用減半歲周為末。
       double UseSun; int SunStatus;
       year.check_SunStatus(half_year, ReqSun, UseSun, SunStatus);
 
//--------------
//     Step 7:
 
int    IUseSun = int(UseSun);
double SunAdd, SunSum, SunX;
       year.get_SunX(SunStatus,UseSun,IUseSun,SunAdd,SunSum, SunX);
             
//--------------
//     Step 8:
 
double EclA_S; int FrontEclipseSign;
 
    if (SunStatus < LeadStart){EclA_S = EclAngle - SunX;}
    if (SunStatus > LagEnd){EclA_S = EclAngle + SunX;}
 
    if (EclA_S < 7){FrontEclipseSign = 1;}
    if (EclA_S > 342){FrontEclipseSign = 1;}
    if (EclA_S > 175){if (EclA_S < 202){FrontEclipseSign = -1;}}
//     Eclipse conditions:
//     342 < EclA_S ----> FrontEclipseSign = 1
//     202 < EclA_S < 342 ----> No eclipse
//     175 < EclA_S < 202 ----> FrontEclipseSign = -1
//       7 < EclA_S < 175 ----> No eclipse
//     EclA_S < 7 ----> FrontEclipseSign = 1
 
//--------------
//     Step 9:
 
int    MoonStatus = 0;double UseMoon = 0; double MoonCh,MoonAdd, moon_x12_2,SoD, MoonSolar;
       year.get_UseMoon(ReqMoon,UseMoon,SoD,moon_x12_2,MoonSolar,MoonStatus,MoonCh,MoonAdd);
 
//--------------
//     Step 10:
 
double MoonSum, MoonX, MoonSp;
       year.get_MoonX(MoonStatus,MoonCh,MoonAdd,MoonSolar,MoonSum,MoonX,MoonSp);
 
//---------------
//     Step 11:
 
double SunXS, MoonXS, SunMoonT;
       year.get_SunMoonT(SunStatus, MoonStatus, MoonX,SunX, MoonSp,SunXS,MoonXS, SunMoonT);

double SunMoonT24 = SunMoonT * 24;
 
//---------------
//     Step 12:

double Conjunc0 = ReqMoonDay - int(ReqMoonDay);
double Conjunc0_24 = Conjunc0 * 24;
 
double ConjuncT = ReqMoonDay + SunMoonT;
       if (ConjuncT < 0){ConjuncT = 1 + ConjuncT;}
 
//double ConjuncTDay = int(ConjuncT);
double ConjuncTHour = ConjuncT - int(ConjuncT);
double ConjuncT24 = ConjuncTHour * 24;
 
//------------------
//     Step 13:
 
double time_from_noon =(ConjuncT-int(ConjuncT))-0.5;
double observed_time =(0.5-abs(time_from_noon))* time_from_noon/0.96;
double MaxPhase = ConjuncTHour + observed_time;
 
double TNoon = time_from_noon * 24;
double SeeT = observed_time * 24 * 60;
double MaxPhase24 = MaxPhase * 24;
 
//------------------
//     Step 14:
 
//     Calculate half day hours
double Sunrise_time
     = 2068.3 +(0.122928-(IUseSun -53)*0.000608)*(IUseSun - 1)/2 *IUseSun + 0.06*IUseSun ;
 
double Sunset_time = 10000 - Sunrise_time;
 
double HalfDayFan = Sunset_time - 5000;
 
       if (SunStatus == LagEnd){HalfDayFan = 5000 - HalfDayFan;}
       if (SunStatus == LeadStart){HalfDayFan = 5000 - HalfDayFan;}
 
//----------------
//     Step 15:
 
double Net_time_from_noon = time_from_noon + observed_time;
       Net_time_from_noon = abs(Net_time_from_noon*10000);
 
//----------------
//     Step 16:
 
//     推食甚入盈縮曆:　置原得盈縮曆，加入定朔大餘及食甚定分，即得。
double EclRef = ReqSun + int(ConjuncT)+ MaxPhase;
double EclSun = EclRef - ReqMoonDay;
double EclUseSun;
int    EclSunStatus;
       year.check_SunStatus(half_year,EclSun, EclUseSun, EclSunStatus);
 
//----------------
//     Step 17:
 
//     推食甚盈縮差:　依步氣朔求之。
double Ecl_SunAdd, Ecl_SunSum,Ecl_SunX;
int    IEclUseSun = int(EclUseSun);
       year.get_SunX(EclSunStatus, EclUseSun,IEclUseSun, Ecl_SunAdd, Ecl_SunSum,Ecl_SunX);
 
//-----------------
//     Step 18:
 
//     推食甚入盈縮曆行定度:　置食甚入盈縮曆，以盈縮差，盈加縮減之，即得。
double NetSun = 0;
       if (SunStatus < LeadStart){NetSun = EclSun - Ecl_SunX;}
       if (SunStatus > LagEnd){NetSun = EclSun + Ecl_SunX;}
//------------------
//     Step 19:
 
double NS_Ref;
int NetSunStatus;
//     推南北汎差:　視食甚入盈縮曆行定度，在周天象限已下為初限(NS_Ref)，已上與半歲周相減為末限(NS_Ref)。
//               以初末限自之，如一千八百七十度而一，得數，置四度四十六分減之，餘為南北汎差(NS_X)。
       year.check_SunStatus(half_year,NetSun, NS_Ref, NetSunStatus);
double NS_X = 4.46 -((NS_Ref*NS_Ref)/1870);
 
//     推南北定差:　置南北汎差，以距午定分乘之，如半晝分而一，以減汎差，餘為南北定差。若汎差數少，即反減之。
double NS_Sh = NS_X-(NS_X* Net_time_from_noon)/HalfDayFan;
double NS_ = NS_Sh;
 
//------------------
//     Step 20:
 
//     推東西汎差:　置半歲周，減去食甚入盈縮曆行定度，餘以食甚入盈縮曆行定度乘之，以一千八百七十除之為度，即東西汎差(EW_X)。
double EW_Ref = NetSun;
       if (EW_Ref > half_year){EW_Ref = EW_Ref - half_year;}
 
double EW_X =((half_year - EW_Ref)*EW_Ref)/1870;
 
//     推東西定差:　置東西汎差，以距午定分乘之，如二千五百度而一，視得數在東西汎差以下，即為東西定差。
//     若在汎差已上，倍汎差減之，餘為定差。
double EW_H = (EW_X * Net_time_from_noon)/2500;
 
double EW_Sh = EW_H;
       if( EW_H > EW_X ){EW_Sh = 2 * EW_X - EW_H;}
double EW_ = EW_Sh;
 
//------------------
//     Step 21:
 
//     盈曆中前，縮曆中後者，正交減，中交加。盈曆中後，縮曆中前者，正交加，中交減。
int    AfterNoonSign = 1;
       if (time_from_noon < 0){AfterNoonSign = -1;}
 
       switch(SunStatus)
       {
        case LagStart:
                     NS_Sh = FrontEclipseSign * NS_Sh;
                     EW_Sh = -AfterNoonSign * FrontEclipseSign * EW_Sh; break;
 
        case LagEnd:
                     NS_Sh = -FrontEclipseSign * NS_Sh;
                     EW_Sh = -AfterNoonSign * FrontEclipseSign * EW_Sh; break;
 
        case LeadStart:
                     NS_Sh = -FrontEclipseSign * NS_Sh;
                     EW_Sh = AfterNoonSign * FrontEclipseSign * EW_Sh; break;
 
        case LeadEnd:
                     NS_Sh = FrontEclipseSign * NS_Sh;
                     EW_Sh = AfterNoonSign * FrontEclipseSign * EW_Sh; break;
        default:
                     break;
       }
//----------------
//     Step 22:
 
double EclRefA = 0;
double EclipseFrac = 0;
double NetEclA_S = EclA_S;
//     推正交中交定限度:　視日食在正交者置正交度，在中交者置中交度，以南北東西二定差加減之，即得。
       if(FrontEclipseSign == -1){EclRefA = 188.05 + NS_Sh + EW_Sh;}
       if (FrontEclipseSign == 1){EclRefA = 357.64 + NS_Sh + EW_Sh;}
       if (EclA_S < 7){NetEclA_S = EclA_S + 363.7934196;}
 
//-----------------
//     Step 23:
 
double PosA = NetEclA_S - EclRefA;
double Ecl_Pos = FrontEclipseSign * PosA;
 
int    EclipseRef = 6;
       if (Ecl_Pos < 0){EclipseRef = 8;}
 
       EclipseFrac = (EclipseRef - abs(PosA))/EclipseRef;
 
//-----------------
//     Step 24:
 
double EclM = ReqMoon + SunMoonT;
int    EclMoonStatus;
double EclM_Used,EclSoD,EclM_x12_2,EclMoonSolar,EclM_Ch,EclMoonAdd;
       year.get_UseMoon(EclM,EclM_Used,EclSoD,EclM_x12_2,EclMoonSolar,EclMoonStatus,EclM_Ch,EclMoonAdd);
 
double EclMoonSum,EclMoonX,EclM_Sp;
       year.get_MoonX(EclMoonStatus,EclM_Ch,EclMoonAdd,EclMoonSolar, EclMoonSum, EclMoonX, EclM_Sp);
 
//     定限行度:以定限,取立成內行度，遲用遲，疾用疾，內減日行分八分二十秒，得之。
double EclM_SunSp = EclM_Sp - 0.082;
 
//------------------
//     Step 25:
double HalfEclT =(sqrt((20 - EclipseFrac*10)*EclipseFrac*10)*57.40)/EclM_SunSp;
 
double EclBegin_hour = MaxPhase - HalfEclT/10000;
double EclBegin24 = EclBegin_hour * 24;
 
double EclEnd_hour = MaxPhase + HalfEclT/10000;
double EclEnd24 = EclEnd_hour * 24;
 
//MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
if ((calculated_year == 1605)||(calculated_year == 1617)||(calculated_year == 1620)
    ||(calculated_year == 1621)||(calculated_year == 1623))
{
//     LunStep 3:
double LunSun = fmod(ReqSun + Shuoshi_C/2, YearDays);
double LunEcl = fmod(ReqEclipse + Shuoshi_C/2, Jiaozhong_C);
double LunDay = fmod(ReqMoonDay + Shuoshi_C/2, 60);
double LunMoon = fmod(ReqMoon + Shuoshi_C/2, Zhuanzhong_C);
//-------------
//     LunStep 4:
//     推交常度:　置望下入交汎分，乘月平行，如日食法。
double LunEclA = LunEcl * MoonMeanSpeed;
 
//-------------
//     LunStep 5:
double LunChSun = LunSun;
       if (LunSun > half_year){LunChSun = LunSun - half_year;}
 
//-------------
//     LunStep 6:
double LUseSun; int LSunStatus;
year.check_SunStatus(half_year,LunSun, LUseSun, LSunStatus);
 
//-------------
//     LunStep 7:
double LSunX = 0;
year.get_LSunX(LSunStatus,LUseSun,LSunX);
 
//-----------------------------
//     LunStep 8:
//     推交定度:　置交常度，以望下盈縮差盈加縮減之即得。不及減者，加交終度減之。
double LSunXS = -LSunX;
       if (LSunStatus > LagEnd){LSunXS = LSunX;}
double LunEclA_S = LunEclA + LSunXS;
 
//-----------------------------
//     LunStep 9:
double LUseMoon,LSoD,Lmoon_x12_2,LMoonSolar,LMoonCh,LMoonAdd; int LMoonStatus;
year.get_UseMoon(LunMoon,LUseMoon,LSoD,Lmoon_x12_2,LMoonSolar,LMoonStatus, LMoonCh,LMoonAdd);
 
//-----------------------------
//     LunStep 10:
double LMoonX, LMoonSp,LMoonCh_;
double LMoonCh_more = LMoonCh+0.08306;
year.get_LMoonX(LMoonStatus,LMoonCh_more,LMoonX, LMoonSp,LMoonCh_);
 
//-----------------------------
//     LunStep 11:
 
double LMoonXS = -LMoonX;
if (LMoonStatus > LagEnd){LMoonXS = LMoonX;}
 
double LSunMoonX = LSunXS - LMoonXS;
double LSunMoonT=LSunMoonX*0.082/LMoonSp;
 
//------------------------------
//     LunStep 12:
double LConjuncT = LunDay + LSunMoonT;
if (LConjuncT < 0){LConjuncT = 1 + LConjuncT;}
 
int LConjuncTDay=int(LConjuncT);
double LConjuncTime=LConjuncT - int(LConjuncT);
 
//------------------------------
//     LunStep 13:
//     月食：视定望分在日周四分之一已下，为卯前； 已上，覆减半周，为卯后；
//         在四分之三已下，减去半周，为酉前；已上，覆减日周，为酉后。
//         以卯酉前后分自乘，退二位，如四百七十八而一，为时差；
//         子前以减，子后以加，皆加减定望分，为食甚定分；各依发敛求之，即食甚辰刻。
double LUseTime;
       LUseTime = LConjuncTime;
       if(LConjuncTime > 0.25){LUseTime = 0.5 - LConjuncTime;}
       if(LConjuncTime > 0.5){LUseTime = LConjuncTime - 0.5;}
       if(LConjuncTime > 0.75){LUseTime = 1 - LConjuncTime;}
double LunSeeT = LUseTime*LUseTime/4.78;
       if(LConjuncTime > 0.5){LunSeeT = -LunSeeT;}
double LunMaxT = LConjuncTime + LunSeeT;
double LunMaxT24 = LunMaxT * 24;
 
//--------------------------------------------------
//     LunStep 16:
//     推食甚入盈縮曆:　置原得盈縮曆，加入定朔大餘及食甚定分，即得。
double LunEclRef = LunSun + int(LConjuncT)+ LunMaxT;
double LunEclSun = LunEclRef - LunDay;
double LEcl_ChSun = LunEclSun;
       if (LunEclSun > half_year){LEcl_ChSun = LunEclSun - half_year;}
double LEclUseSun;
int    LEclSunStatus;
       year.check_SunStatus(half_year,LunEclSun, LEclUseSun, LEclSunStatus);
 
//-----------------
//     LunStep 17:
 
//     推食甚盈縮差:　依步氣朔求之。
double LEcl_SunX = 0;
year.get_LSunX(LEclSunStatus,LEclUseSun,LEcl_SunX);
 
//-----------------
//     LunStep 18:
 
//     推食甚入盈縮曆行定度:　置食甚入盈縮曆，以盈縮差，盈加縮減之，即得。
double LEcl_SunX_=-LEcl_SunX;
if (LEclSunStatus > LagEnd){LEcl_SunX_=LEcl_SunX;}
double LEclNetSun = LEcl_ChSun + LEcl_SunX_;       
 
double LEclNetSun_= LEclNetSun;
       if(LEclSunStatus > LagEnd){LEclNetSun_=LEclNetSun+182.62875;}

LEclNetSun_=LEclNetSun_+ MinSun;  
    
int total_year=calculated_year-Base_year;
double Added_drift=0.015*total_year;
double Ref_drift=10-Added_drift;
double Ref_drift_check=fmod(SumDay+315.1075,365.2575)-305.1075;
double Next_Ref_drift=Ref_drift_check-4.3445;
 
double Sun_drift=Next_Ref_drift/1.0849;
double Net_drift=4 + Sun_drift;
double LEclNetSun_All=LEclNetSun_+ Net_drift;
 
double LEclNetSun_From=0;
       if(calculated_year==1605){LEclNetSun_From=LEclNetSun_All-277.9775;}
       if(calculated_year==1617){LEclNetSun_From=LEclNetSun_All-51.08;}
       if(calculated_year==1620){LEclNetSun_From=LEclNetSun_All-160.8175;}
       if(calculated_year==1623){LEclNetSun_From=LEclNetSun_All-103.6975;}
//------------------
//     LunStep 22:
 
//求月食入陰陽曆去交前後度
// 視交定度，在交中度(181.8967)已下，為陽曆；已上，減去交中，為陰曆。
// 視入陰陽曆，在後准十五度半已下，為交後度；
// 前准一百六十六度三十九分六十八秒已上，覆減交中，余為交前度及分。
 
double SexLunEclA=LunEclA_S;
if (LunEclA_S > 181.8967){SexLunEclA = 181.8967 - LunEclA_S;}
 
int LunEclFront=-1;
    if(abs(SexLunEclA)<15.5){LunEclFront=1;}
 
double SexLunEclA_= SexLunEclA;
if (LunEclA_S > 166.3968){SexLunEclA_=181.8967-LunEclA_S;LunEclFront=-1;}
if (SexLunEclA < -166.3968){SexLunEclA_=-181.8967-SexLunEclA;LunEclFront=-1;}
 
//-------------------------------------------------
//     LunStep 23:
//求月食分秒
// 視去交前後度，
// 不用南北東西差者。用減食限，不及減者不食。
// 余如定法而一，為月食之分秒。
 
double LunFraction=(13.05-abs(SexLunEclA_))/8.7;
//------------------
//     LunStep 24:
 
double LEclM = LunMoon+LSunMoonT;
int    LEclMoonStatus;
double LEclM_Used,LEclSoD,LEclM_x12_2,LEclMoonSolar,LEclM_Ch,LEclMoonAdd;
       year.get_UseMoon(LEclM,LEclM_Used,LEclSoD,LEclM_x12_2,LEclMoonSolar,LEclMoonStatus,LEclM_Ch,LEclMoonAdd);
 
double LEclMoonX,LEclM_Sp,LEclM_Ch_;
double LEclM_Ch_more=LEclM_Ch + 0.08306;
       year.get_LMoonX(LEclMoonStatus,LEclM_Ch_more,LEclMoonX, LEclM_Sp,LEclM_Ch_);
 
//     定限行度:以定限,取立成內行度，遲用遲，疾用疾，內減日行分八分二十秒，得之。
double LEclM_SunSp = LEclM_Sp - 0.082;
 
//     LunStep 25:
double LHalfEclT=(sqrt((30-LunFraction*10)*LunFraction*10)*49.20)/(LEclM_SunSp*10000);  // Shoushi uses 57.40
double LEclBeginT=LunMaxT - LHalfEclT;
    if(LEclBeginT<0){LEclBeginT=1+LEclBeginT;}
double LEclBeginT24=LEclBeginT*24;
double LEclEndT=LunMaxT+LHalfEclT;
double LEclEndT24=LEclEndT*24;
 
if (((calculated_year==1620)&(MonthNum==13))||(calculated_year==1621))
{
printf("%2d %4d(%2d) LunEcl=%6.4f\tLunMoon=%6.4f\tLunFraction=%6.4f\tLunMaxT24=%6.4fh\tLEclNetSun_=%6.4f\n",
            index,calculated_year,MonthNum,LunEcl,LunMoon,LunFraction,LunMaxT24,LEclNetSun_);
}
//------------------------------------------------
if (((calculated_year==1605)&(MonthNum==3))
||((calculated_year==1617)&(MonthNum==8))
||((calculated_year==1620)&(MonthNum==12))
||((calculated_year==1623)&(MonthNum==10))
||((calculated_year==1621)&(MonthNum==6)))// lunar(about 15 days earlier) and solar elipse occurred in 1621(6)
{
printf("%2d %4d(%2d) ReqEcl(交泛)=%9.6f\tLunEcl(食月望交泛)=%10.7f\n", index,calculated_year,MonthNum,ReqEclipse,LunEcl);
if (calculated_year==1617){print_CalculatedRecords_1617();}
if (calculated_year==1620){print_CalculatedRecords_1620();}
if (calculated_year==1623){print_CalculatedRecords_1623();}
 
cout << endl;
printf("\tSumDay=%8.3f\t[LunStep 3]: LunDay(食月經望)=%10.7f\t\tLunSun=%5.5f\tLunMoon=%5.5f\n",SumDay,LunDay,LunSun,LunMoon);
printf("\tWinter_solitice=%4.3f\t[LunStep 4]: LunEclA(交常度)=%10.8f\n", Winter_solitice,LunEclA);
printf("\tLeapDay=%10.7f\t[LunStep 5]: LunChSun(食月盈縮曆)=%10.8f\n",LeapDay,LunChSun);
printf("\tSumDay + Q=%8.3f\t[LunStep 6]: LUseSun=%10.8f(%1d)\n", SumDay + Qiying_P,LUseSun, LSunStatus);
printf("\tRefMoonDay=%9.6f\t[LunStep 7]: LSunX(盈縮差)=%11.9f\n",RefMoonDay, LSunX);
printf("\tRefSun=%9.6f\t[LunStep 8]: LSunXS=%11.9f \tLunEclA_S(交定度)=%11.9f\n",RefSun,LSunXS,LunEclA_S);
printf("\tRefMoon=%9.6f\t[LunStep 9]: LUseMoon(食月望遲疾曆)=%10.7f(%1d) \tLmoon_x12_2(遲疾限)=%10.7f\tLMoonCh=%5.5f\n",
                                                         RefMoon,LUseMoon,LMoonStatus,Lmoon_x12_2, LMoonCh);
printf("\tRefEclipse=%9.6f\t[LunStep 10]:(+0.08306=>)LMoonCh_=%10.7f\tLMoonX(遲疾差)=%10.8f\tLMoonSp=%9.7f\n",
                                                    RefEclipse,LMoonCh_, LMoonX, LMoonSp);
if (calculated_year==1605)
{
  printf("\t\t\t\t\t\t古今律厯考卷五十(787-563)原計算推遲末限=總限(168.08306限)減遲限(138.3628779限)=29.7201821限\n");
  printf("\t\t\t\t\t\t有別大統以168限或167限減遲限計算,因此 LMoonCh_ 加回0.08306\n");
  printf("\t\t\t\t\t\t該遲行度以29限整數計算,得出1.1823。現改回29.7201821限計算,於是 LMoonSp=1.1814954\n\n");
}
printf("\t\t\t\t[LunStep 11]:LMoonXS=%10.8f\tLSunMoonX=%10.8f \tLSunMoonT(加減差)=%10.8f\n",LMoonXS,LSunMoonX,LSunMoonT);
printf("\tReqMoonDay=%9.6f\t[LunStep 12]:LConjuncT(食月定望)=%10.8f\tLConjuncTDay=%2d\t\tLConjuncTime=%10.8f\n",
          ReqSun,LConjuncT,LConjuncTDay,LConjuncTime);
printf("\tReqSun=%9.6f\t[LunStep 13]:LunSeeT(時差)=%2.7f\tLunMaxT(食甚定分)=%2.8f  =>\tLunMaxT24=%6.7fh",
          UseSun, LunSeeT,LunMaxT,LunMaxT24);
    if(calculated_year==1620){printf("[但丑正初刻略遲]");}
printf("\n");
printf("\tReqMoon=%6.6f\n",ReqMoon);  print_LunStep_13_Extra(calculated_year); 
printf("\n");
printf("\t\t\t\t[LunStep 16]:LunEclSun=%8.5f \tLEcl_ChSun(食甚入盈縮曆)=%10.8f \tLEclUseSun=%10.8f(%1d)\n",
          LunEclSun,LEcl_ChSun,LEclUseSun,LEclSunStatus);
printf("\t\t\t\t[LunStep 17]:LEcl_SunX(食甚盈縮差)=%2.8f\n",LEcl_SunX);
printf("\t\t\t\t[LunStep 18]:LEcl_SunX_=%2.8f \tLEclNetSun(食甚入盈縮曆定度)=%2.8f(%1d)\n",LEcl_SunX_,LEclNetSun,LEclSunStatus);
if(calculated_year == 1605){printf("\t\t\t\t\t     原計算=103.57495347(食甚入盈縮曆)+2.33579615(食甚盈縮差)=105.91074962\n");
                            printf("\t\t\t\t\t\t     但古今律厯考卷五十(787-565)給出錯誤食甚入盈曆定度=105.93074962\n\n");}
printf("\t\t\t\t\t\tLEclNetSun_= LEclNetSun\n");
printf("\t\t\t\t\t\tif(LEclSunStatus > LagEnd){LEclNetSun_=LEclNetSun+182.62875}\n");
printf("\t\t\t\t\t\t\tAs LEclSunStatus=%1d\tLEclNetSun_=%10.8f\n\n",LEclSunStatus,LEclNetSun_);
 
printf("\t\t\t\t\t\ttotal_year=calculated_year-Base_year=%3d => Added_drift=0.015*total_year=%7.5f\n",total_year,Added_drift);
printf("\t\t\t\t\t\tRef_drift(天正加時赤道度)=10-Added_drift=%7.5f(邢雲路計算方法)\n",10-Added_drift);
printf("\t\t\t\t\t\t  (授時曆計算方法:Ref_drift_check=fmod(SumDay+315.1075,365.2575)-305.1075=%7.5f)\n",Ref_drift_check);
printf("\t\t\t\t\t\tNext_Ref_drift=Ref_drift_check-4.3445=%7.5f\n\n",Next_Ref_drift);
 
printf("\t\t\t\t\t\tSun_drift=Next_Ref_drift/1.0849=%7.5f\n",Sun_drift);
printf("\t\t\t\t\t\tNet_drift=4+Sun_drift=%7.5f\n",Net_drift);
printf("\t\t\t\t\t\tLEclNetSun_All=LEclNetSun_ + Net_drift=%7.5f\n",LEclNetSun_All);
 
if(calculated_year==1605){printf("\t\t\t\t\t\tLEclNetSun_From(食甚月躔黃道軫宿)=LEclNetSun_All-277.9775(翼宿)=%7.5f\n",LEclNetSun_From);}
if(calculated_year==1617) {printf("\t\t\t\t\t\tLEclNetSun_From(食甚月躔黃道虚宿)=LEclNetSun_All-51.08(女宿)=%7.5f\n",LEclNetSun_From);}
if(calculated_year==1620)
                {printf("\t\t\t\t\tLEclNetSun_From(食甚月躔黃道畢宿)=LEclNetSun_All-160.8175(昴宿)=%7.5f",LEclNetSun_From);
                 printf("(幾乎=14.3444)\n");
                 printf("\t\t\t\t\t\t\t而授時則推....月離黃道畢宿一十四度三十四分四十四秒(14.3444 與計算值符合=>記錄可靠)\n");
                 printf("\t\t\t\t\t\t  邢雲路...按新法推...食甚月離黃道畢宿一十四度三十分一十秒(14.3010)\n");
                 printf("\t\t\t\t\t\t按日行每天一度,估計授時比新法遲(14.3444-14.3010)*24=1.0416h\n");
                 printf("\t\t\t\t\t\t與食甚遲的時刻比較,差距亦大致相同 LunMaxT24-0.9888h=%1.7fh,表示該些記錄亦可靠\n",LunMaxT24-0.9888);}
if(calculated_year==1623)
                {printf("\t\t\t\t\t\tLEclNetSun_From(食甚月躔黃道奎宿)=LEclNetSun_All-103.6975(壁宿)=%7.5f\n",LEclNetSun_From);}
 
printf("\t\t\t\t[LunStep 22]:\n");
      printf("\t\t\t\t\t求交常交定度:(refer to LunStep 4 and LunStep 8)\n");
      printf("\t\t\t\t\t\t置經朔望入交泛日及分秒(LunEcl)，以月平行度(%10.8f)乘之，為交常度(LunEclA)；\n",MoonMeanSpeed);
      printf("\t\t\t\t\t\t以盈縮差(LSunXS)盈加縮減之，為交定度(LunEclA_S)。\n\n");
 
      printf("\t\t\t\t\t求月食入陰陽曆去交前後度:\n");
      printf("\t\t\t\t\t\t視交定度，在交中度(181.8967)已下，為陽曆；已上，減去交中，為陰曆。\n");
      printf("\t\t\t\t\t\t視入陰陽曆，在後準十五度半已下，為交後度；\n");
      printf("\t\t\t\t\t\t前準一百六十六度三十九分六十八秒已上，覆減交中，餘為交前度及分。\n\n");
      printf("\t\t\t\t\t\t交中度=181.8967\n");
      printf("\t\t\t\t\t\tSexLunEclA=LunEclA_S=%11.9f\n",LunEclA_S);
      printf("\t\t\t\t\t\tif (LunEclA_S > 181.8967){SexLunEclA=181.8967-LunEclA_S}\n");
      printf("\t\t\t\t\t\tResult:SexLunEclA(陰陽曆,負號[-]表示陰曆;正數表示陽曆)=%11.9f\n\n",SexLunEclA);
 
      printf("\t\t\t\t\t\tLunEclFront=-1(負數[-1]表示交前度);\n");
      printf("\t\t\t\t\t\tif(abs(SexLunEclA)<15.5){LunEclFront=1(正數[1]表示交後度);}\n\n");
 
      printf("\t\t\t\t\t\tSexLunEclA_=SexLunEclA=%11.9f\n",SexLunEclA);
      printf("\t\t\t\t\t\tif (LunEclA_S > 166.3968){SexLunEclA_=181.8967-LunEclA_S; LunEclFront=-1}\n");
      printf("\t\t\t\t\t\tif (SexLunEclA < -166.3968){SexLunEclA_=-181.8967-SextLunEclA; LunEclFront=-1}\n");
      printf("\t\t\t\t\t\tResult: SexLunEclA_(LunEclFront) = 月食入陰陽曆(交前後度)= %11.9f(%1d)\n", SexLunEclA_,LunEclFront);
      if (calculated_year == 1617) { printf("\t\t\t\t\t\t邢雲路獻七政真數..以此法布算..陰歷交前度一度二十二分零三秒九十五微=>-1.220395(-1)\n");}
 
printf("\t\t\t\t[LunStep 23]:\n");
                             printf("\t\t\t\t\t求月食分秒:\n");
                             printf("\t\t\t\t\t\t視去交前後度，\n");
                             printf("\t\t\t\t\t\t不用南北東西差者。用減食限，不及減者不食。\n");
                             printf("\t\t\t\t\t\t餘如定法而一，為月食之分秒。\n");
      printf("\t\t\t\t\t\tLunFraction=(13.05-abs(SexLunEclA_))/8.7=%11.9f\n", LunFraction);
      print_LunStep_23_Extra(calculated_year);
 
printf("\t\t\t\t[LunStep 24]:\n");
      printf("\t\t\t\t\tLEclM=%2.8f(%1d)\tLEclM_Used(定入遲疾曆)=%2.8f\n",LEclM, LEclMoonStatus,LEclM_Used);
      printf("\t\t\t\t\tLEclM_x12_2(定入遲疾限)=%2.8f\tLEclM_Ch=%2.8f\tLEclM_Sp=%2.8f\n",LEclM_x12_2,LEclM_Ch,LEclM_Sp);
      printf("\t\t\t\t\tLEclM_SunSp(定限行度)=%2.5f\t    LHalfEclT(定用分)=%2.8f",LEclM_SunSp,LHalfEclT);
      printf("(邢雲路用 4920,不是 5740)\n");
      if (calculated_year == 1617)
      {
          printf("\t\t\t\t\t邢雲路獻七政真數..以此法布算..定用分六刻六十八分九十五秒二十一微(0.06689521)\n");
      }
 
printf("\t\t\t\t[LunStep 25]:\n");
      printf("\t\t\t\t\tLEclBeginT(初虧)=%2.6f\t   LEclEndT(復圓)=%2.6f\n",LEclBeginT,LEclEndT);
      if (calculated_year == 1617){printf("\t\t\t\t\t   (邢雲路的初虧: 0.064[2]2089\t\t    復圓: 0.19801131)\n\n");}
      if (calculated_year == 1620){printf("\t\t邢雲路...按新法推..初虧..計九十七刻二十四分(0.9724)\t\t複圓....計一十一刻(0.11 )\n\n");}
      printf("\t\t\t\t\tLEclBeginT24=%6.8fh",LEclBeginT24);if(calculated_year==1620){printf("[初虧子正一刻]是可靠的");}
      printf("\tLEclEndT24=%6.8fh",LEclEndT24);if(calculated_year==1620){printf("[復圓寅初三刻]略為稍遲");}
      printf("\n");                                                             
      print_LunStep_25_Extra(calculated_year);
cout << "----------------------------------------------------------\n";
 }
} //-----------------
       if (EclipseFrac > 0)
       {
         if (calculated_year >= Theo[index].year)
         {
           if (calculated_year > Theo[index].year) { Shoushi[index].EclipseFrac = -99; index = index + 1;}
 
         if (calculated_year == Theo[index].year)
         {
           if (MonthNum >= Theo[index].month)
           {
             Shoushi[index].year = calculated_year;
             Shoushi[index].month = MonthNum;
 
             Shoushi[index].EclipseFrac = EclipseFrac;
             Shoushi[index].EclipseBegin = EclBegin24;
             Shoushi[index].MaxPhase = MaxPhase24;
             Shoushi[index].EclipseEnd = EclEnd24;
 
             Shoushi[index].Winter_solitice = Winter_solitice;
             Shoushi[index].SunStatus = SunStatus;
 
             Shoushi[index].observed_time = observed_time;
             Shoushi[index].HalfDayFan = HalfDayFan;
             Shoushi[index].ReqSun = ReqSun;
             Shoushi[index].ReqMoon = ReqMoon;
 
             Shoushi[index].MoonStatus = MoonStatus;
             Shoushi[index].ConjuncT24 = ConjuncT24;
             Shoushi[index].observed_time = observed_time;
 
             index = index + 1;
           }
         }
             if (EclipseFrac > 0)
             {
                 if (calculated_year == 1621){print_FindingRecords_1621();}
 
    printf("%2d %4d(%1d)", index,calculated_year,MonthNum);
//-------------------------- Next Row for Step [1] to [6]-------------------------
    printf("\t\t\t\t[2]MinSun=%2.6f\n",MinSun);             
    printf("\t[1](365.)SumDay(WinSol)LeapDay\t[2]RefS(RefMDay)RefM(RefEcl)");
    printf("\t[3]ReqS(ReqMDay)ReqM(ReqEcl)\t[4](EclA)[5]ChSun[6]UseSun(S)\n");
    printf("\t(%5.6f)%5.2f(%5.2f)%5.2f\t  %5.2f(%5.2f)%5.2f(%5.2f)",
           YearDays-365,SumDay,Winter_solitice,LeapDay,RefSun,RefMoonDay,RefMoon,RefEclipse);
    printf("\t   %4.2f(%5.2f)%5.2f(%5.2f)\t   (%5.2f)  %5.2f  %5.2f(%1d)\n\n",
           ReqSun,ReqMoonDay,ReqMoon,ReqEclipse,EclAngle,ChSun,UseSun,SunStatus);
//-------------------------- Next Row for Step [7] to [11]-------------------------
    printf("\t[7]Sun_Add/Sum(X)[8]Ft(EclA_S)\t[9]UseM/SoD(M12_2)MSo(M)");
    printf("\t(MoonCh)MAdd[10]Sum(X)Speed\t[11]SunXS/MoonXS/SunMoonT\n");
    printf("\t%5.2f/%5.2f(%4.2f)%1d(%4.2f)",
             SunAdd,SunSum,SunX,FrontEclipseSign,EclA_S);
    printf("\t %5.2f/%5.2f(%5.2f)%1.2f(%1d)\t(%5.2f)%5.2f",
           UseMoon,SoD,moon_x12_2,MoonSolar,MoonStatus,MoonCh,MoonAdd);
    printf("  %1.2f(%1.2f)%1.4f\t   %1.3f/%1.3f/%5.3f(%5.3fh)\n\n",
           MoonSum,MoonX,MoonSp,SunXS,MoonXS,SunMoonT, SunMoonT24);
//-------------------------- Next Row for Step [12] to [16] ------------------------
    printf("\t[12]Conjunc0(T)  [13]TNoon\tSeeT(MaxPhase) [14]Srise/Sset");
    printf("\tHalfDayFan[15]NetTNoon\t\t[16]Ecl)Ref/Sun(S)UseSun\n");
    printf("\t  %1.4fh(%1.4fh)%5.2fh \t  %5.2fm(%5.4fh)  %4.0f/%4.0f",
           Conjunc0_24, ConjuncT24, TNoon,SeeT,MaxPhase24,Sunrise_time,Sunset_time);
    printf("\t  %5.2f\t%5.2f\t\t   %5.2f/%5.2f(%1d)%5.2f\n\n",
           HalfDayFan,Net_time_from_noon,EclRef,EclSun,EclSunStatus,EclUseSun);
//---------------------------- Next Row for Step [17] to [21]-----------------------
    printf("\t[17]Ecl_)SAdd/SSum(SunX)\t[18]NSun(S)[19]NS_R/NS_X(NS_)");
    printf("\t[20]EW_R/EW_X/EW_H(EW_)\t\t[21]A12_Sgn/NS_Sh/EW_Sh\n");
    printf("\t  %5.2f/%5.2f(%4.2f)\t\t  %5.2f(%1d)  %1.2f/%1.2f/(%1.2f)",
           Ecl_SunAdd,Ecl_SunSum,Ecl_SunX,NetSun,NetSunStatus,NS_Ref, NS_X, NS_);
    printf("\t %1.2f/%1.2f/%1.2f(%1.2f)",EW_Ref,EW_X,EW_H,EW_);
    printf("\t\t\t %1d/%5.2f/%5.2f\n\n",AfterNoonSign,NS_Sh,EW_Sh);
//---------------------------- Next Row for Step [22] to [25]-----------------------
    printf("\t[22]EclRefA/NetEclA_S/PosA\t[23]EclPos(Ref)Frac");
    printf("\t\t[24]EclM(M)_Ch(_Sp)(_SunSp)\t[25]HalfEclT(EclBegin)(EclEnd)\n");
    printf("\t  %5.2f/%5.2f/%5.2f\t\t   %5.2f(%1d)%5.4f",
           EclRefA, NetEclA_S,PosA,Ecl_Pos, EclipseRef,EclipseFrac);
    printf("\t\t %5.2f(%1d)%5.2f(%1.2f)(%1.2f)\t   %5.2f(%5.4fh)(%5.4fh)\n\n",
           EclM,EclMoonStatus,EclM_Ch,EclM_Sp,EclM_SunSp,HalfEclT,EclBegin24,EclEnd24);
    printf("----------------------------------------------\n");
    }
   }
  }
   }//Start_End for loop end(1)
  }//Start_End for loop end(2)
 
     Ref_data Eclipse_Comparison(calculated_year,MonthNum);
     int index_end=index;
     Eclipse_Comparison.print_SolarEclipseTable(index_end);
    return(0);
 }// main end


run \| edit \| history \| help	0