C++ Solar eclipse program(Shoushi integrated) 1644 - 1785 from Ideone( Date: August 14, 2014 ), C++ (gcc)

C++ Solar eclipse program(Shoushi integrated) 1644 - 1785 from Ideone( Date: August 14, 2014 )

// C++ Solar eclipse program(Shoushi integrated) input as start_year and end_year
    // author : Mr Lam/ Yee Din email : yeedlam@yahoo.com
    // Copyright reserved if you want to use the program, please contact the author
    // August 7, 2012
     
    #include <iostream>
    using namespace std;
    #include <cmath> // or math.h
     
    const int start_year = 1644; // enter start_year
    const int end_year = 1785; // enter end_year
     
    const double Tropical_year = 365.2425;
    const double half_year = Tropical_year/2;
    const int Base_year = 1281;
    const double pi = 3.14159265;
    //===========================================================
    // 4 Paramters for adjustments
    const double Qiying_P = 55.06, Runying_P = 20.2050;
    const double Jiaoying_P = 26.0388, Zhuanying_P = 13.0205;
    //===========================================================
    // Constants for Moon to determine solar eclipse purpose
    const double Shuoshi_C = 29.530593, Zhuanzhong_C = 27.5546;
    const double Jiaozhong_C = 27.212224;
    const double Jiao_eclipse_C = 2.318369;
    const double half_Zhuanzhong = Zhuanzhong_C/2;
    //===================================================
    // Constant for Moon positions : Lead or Lag
    const double Moon_Ref_C = 1.975993;
    const double MoonMeanSpeed = 13.36875;
    //===================================================
     
    typedef struct{
    int year; int month; int day;
    double eclipseFrac; double eclipseBegin; double eclipseMax; double eclipseEnd;
    double eclipseBeginS; double eclipseEndS;
    double Winter_solitice;
    int SunStatus; int MoonStatus;
    int SunStatusS; int MoonStatusS;
    double observed_time; double half_day_hours;
    double ReqSun; double ReqMoon;
    double half_dayoff;
    double Gon;
     
    double EclipseFracS;
    double MoonSpeedFactor;
    double noon_offset;
     
    double MDayT_hour;
    double observed_timeS;
    double SunXS; double MoonXS; double SunMoonXS;
    } EclipseTime;
     
    EclipseTime eclipse_shoushi[100];
     
    EclipseTime eclipse_theo[]
    = {
    {1644, 9, 1, 0.303, 10.941, 12.152, 13.354},
    {1648, 6, 21, 0.982, 6.169, 7.213, 8.374},
    {1650, 10, 25, 0.904, 10.267, 11.603, 12.991},
    {1657, 6, 12, 0.611, 4.301, 5.143, 6.049},
    {1658, 6, 11, 0.53, 9.084, 10.501, 12.049},
    {1665, 1, 16, 0.908, 14.822, 16.267, 17.518},
    {1666, 7, 2, 0.988, 15.563, 16.764, 17.835},
    {1669, 4, 30, 0.538, 13.117, 14.386, 15.578},
    {1671, 9, 3, 0.122, 16.426, 16.968, 17.491},
    {1676, 6, 11, 0.092382, 18.3401, 18.779, 19.207}, //eclipseFrac & eclipseBegin from www.eclipse.org.uk
     
    {1681, 9, 12, 0.383, 8.201, 9.225, 10.333},
    {1685, 11, 26, 0.264, 14.663, 15.478, 16.247},
    {1688, 4, 30, 0.96, 8.182, 9.291, 10.495},
    {1690, 9, 3, 0.216, 6.787, 7.476, 8.206},
    {1691, 2, 28, 0.375, 12.289, 13.526, 14.681},
    {1692, 2, 17, 0.532, 11.276, 12.787, 14.223},
    {1695, 12, 6, 0.973, 15.347, 16.518, 17.570},
    {1697, 4, 21, 0.996, 7.795, 8.899, 10.109},
    {1704, 11, 27, 0.577, 12.241, 13.624, 14.938},
    {1706, 5, 12, 0.574, 18.261, 19.074, 19.839},
     
    {1708, 9, 14, 0.541, 16.541, 17.465, 18.328},
    {1709, 9, 4, 0.523, 5.859, 6.681, 7.56},
    {1712, 7, 4, 0.626, 3.784, 4.684, 5.659},
    {1715, 5, 3, 0.604, 18.146, 18.955, 19.715},
    {1719, 2, 19, 0.741, 14.586, 15.969, 17.186},
    {1720, 8, 4, 0.736, 10.691, 12.225, 13.717},
    {1721, 7, 24, 0.39, 17.33, 18.202, 19.012},
    {1730, 7, 15, 0.829, 11.097, 12.85, 14.604},
    {1731, 12, 29, 0.867, 6.3417, 7.409, 8.602},// eclipseBegin from eclipse.org.uk
    {1735, 10, 16, 0.825, 7.696, 8.869, 10.168},
     
    {1742, 6, 3, 0.702, 6.593, 7.512, 8.511},
    {1745, 4, 2, 0.082, 10.852, 11.496, 12.143},
    {1746, 3, 22, 0.765, 9.315, 10.832, 12.43},
    {1747, 8, 16, 0.249, 17.001, 17.7, 18.361},
    {1751, 5, 25, 0.477, 6.714, 7.543, 8.442},
    {1758, 12, 30, 0.9, 15.005, 16.265, 17.3861}, // eclipseEnd from eclipse.org.uk
    {1760, 6, 13, 0.985, 16.312, 17.343, 18.284},
    {1762, 10, 17, 0.756, 16.68, 9999, 9999},// eclipse.org.uk has no data on Beijing
    {1763, 10, 7, 0.734, 6.0403, 6.933, 7.891},// eclipseBegin from eclipse.org.uk
    {1769, 6, 4, 0.333, 16.984, 17.686, 18.348},
     
    {1770, 5, 25, 0.414, 7.52, 8.394, 9.339},
    {1773, 3, 23, 0.415, 13.015, 14.39, 15.644},
    {1774, 9, 16, 0.387, 7.259, 8.251, 9.337},
    {1775, 8, 26, 0.468, 11.405, 12.923, 14.334},
    {1776, 1, 21, 0.145, 9.526, 10.197, 10.889},
    {1784, 8, 16, 0.187, 5.574, 6.27, 7.019},
    {1785, 8, 5, 0.459, 6.634, 7.691, 8.861},
     
    {99999, 0, 0, 0, 0, 0, 0}
    };
     
    enum Status{LagStart, LagEnd, LeadStart, LeadEnd};
    struct Status_type{char *sun; char *moon;};
     
    class Ref_data
    {
    private:
    Status SunStatus, MoonStatus;
     
    protected:
    int calculated_year; int total_months;
     
    double total_days; double Winter_solitice; double LeapDay;
    double RefEcl; double RefSun; double RefMoonDay; double RefMoon;
    double ReqEcl; double ReqSun; double ReqMoonDay; double ReqMoon;
    double sun_days; int intSun_days;
    double sun_add_minus; double sun_sum; double SunX;
    double moon_add_minus; double moon_sum; double moon_difference; double solar_days;
    double SolarEclipse_angle; int Front;
     
    public:
    Ref_data(int m,int n)
    {
    calculated_year = m;total_months = n;
     
    RefEcl = 0; RefSun = 0; RefMoonDay = 0; RefMoon = 0;
    ReqEcl = 0; ReqSun = 0; ReqMoonDay = 0; ReqMoon = 0;
     
    total_days = 0; LeapDay = 0; Winter_solitice = 0; sun_days = 0;
     
    sun_add_minus = 0; sun_sum = 0; SunX = 0;
    }
    ~Ref_data(){}
     
    void get_basic(double& total_days, double& Winter_solitice, double& LeapDay)
    {
    total_days = (calculated_year - Base_year)* Tropical_year;
    Winter_solitice = fmod((total_days + Qiying_P + int(abs(total_days))* 60), 60);
    LeapDay = fmod((total_days + Runying_P + int(abs(total_days)) * Shuoshi_C), Shuoshi_C);
    }
     
    void get_ref(double total_days, double Winter_solitice, double LeapDay,
    double& BaseSun,double& RefEcl,double& RefSun,double& RefMoonDay,double& RefMoon,double& MinSun)
    {
    RefMoonDay = Winter_solitice - LeapDay;
    if (RefMoonDay < 0){RefMoonDay = RefMoonDay + 60;}
     
    BaseSun = Tropical_year/2 - LeapDay;
     
    RefMoon = fmod((total_days + Zhuanying_P - LeapDay + int(abs(total_days))* Zhuanzhong_C),Zhuanzhong_C);
     
    RefEcl = fmod((total_days + Jiaoying_P - LeapDay + int(abs(total_days))*Jiaozhong_C), Jiaozhong_C);
     
    //--------------------------------------------------------------------------------------
    // Shoushi drift problems referring to base year 1628 from Adam Schall 200 eclipse table
    // calculated data from http://ideone.com/BlkMI
    // for MoonDay = -(0.005499/347) => approx -0.0000158473
    // for sun = -(5.99089/347) => -0.0172648
    // for moon = (0.162147/347) => 0.000467281
    // for eclipse = (0.0122486/347) => 0.0000352987
    //--------------------------------------------------------------------------------------
    // slightly modified according to data record in Chongzhen Lishu(page 788-34) as :
     
    RefMoonDay = -0.0081 + RefMoonDay + 0.01041 - 0.01041/347*(calculated_year - 1281);
    if (RefMoonDay < 0){RefMoonDay = RefMoonDay + 60;}
     
    MinSun = 6.572 - 6.1419 + 6.1419/347 *(calculated_year - 1281);
    RefSun = BaseSun - MinSun;
     
    RefMoon = 0.0402 + RefMoon + 0.16309/347 *(calculated_year - 1281);
    RefEcl = RefEcl + 0.01735/347*(calculated_year - 1281);
    }
     
    void get_required(double RefEcl, double RefSun, double RefMoonDay, double RefMoon,
    double& ReqEcl, double& ReqSun, double& ReqMoonDay, double& ReqMoon)
    {
    ReqEcl = fmod(((Jiao_eclipse_C * total_months)+ RefEcl),Jiaozhong_C);
    ReqSun = RefSun + (Shuoshi_C * total_months);
    ReqSun = fmod(ReqSun, Tropical_year);
    ReqMoonDay = fmod((RefMoonDay + Shuoshi_C * total_months),60);
    ReqMoon = fmod(( RefMoon + Moon_Ref_C * total_months),Zhuanzhong_C);
    }
     
    void checkSunStatus(double ReqSun, double& sun_days, int& SunStatus)
    {
    SunStatus = LagStart; sun_days = 0;
     
    if (ReqSun < half_year)
    {
    SunStatus = LagStart; sun_days = ReqSun;
    if (ReqSun > 93.712025){sun_days = half_year - sun_days; SunStatus = LagEnd;}
    }
     
    if (ReqSun > half_year )
    {
    SunStatus = LeadStart; sun_days = ReqSun - half_year ;
    if (sun_days > 88.909225){sun_days = half_year - sun_days; SunStatus = LeadEnd;}
    }
    }
     
    void calculateSunOffset(int SunStatus, double sun_days, int intSun_days,
    double& sun_add_minus, double& sun_sum, double& SunX)
    {
    sun_add_minus = sun_sum = SunX = 0;
     
    if (fmod(SunStatus, 3.0) == 0) // SunStaus = LagStart & LeadEnd
    {
    sun_add_minus = 484.8473 -(4.4362 + 0.0162 *(intSun_days - 1)/2)* intSun_days;
    sun_sum = (487.06 -(2.21 + 0.0027 * intSun_days)* intSun_days )* intSun_days;
    }
     
    if (fmod(SunStatus, 3.0) > 0) // SunStatus = LagEnd & LeadStart
    {
    sun_add_minus = 510.8569 -(4.9386 + 0.0186 *(intSun_days - 1)/2)*intSun_days;
    sun_sum = (513.32 -(2.46 + 0.0031 * intSun_days)*intSun_days)*intSun_days;
    }
    SunX =((sun_days - intSun_days)*sun_add_minus + sun_sum)/10000;
    }
     
    // Calculate Solar eclipse angle
     
    void get_SolarEclipse_angle(int SunStatus,double ReqEcl_angle,double SunX,double& SolarEclipse_angle,int& Front)
    {
    if (SunStatus < LeadStart){SolarEclipse_angle = ReqEcl_angle - abs(SunX);}
    if (SunStatus > LagEnd){SolarEclipse_angle = ReqEcl_angle + abs(SunX);}
     
    if (SolarEclipse_angle < 7){Front = 1;}
    if (SolarEclipse_angle > 342){Front = 1;}
     
    if(SolarEclipse_angle > 202)
    {
    if (SolarEclipse_angle < 342){
    // cout << endl;
    // cout << "==============================" << endl;
    // cout << " no solar eclipse observed" << endl;
    // cout << "==============================" << endl;
    // cout << endl;
    }
    }
     
    if(SolarEclipse_angle < 175)
    {
    if (SolarEclipse_angle > 7){
    // cout << endl;
    // cout << "=================================" << endl;
    // cout << " no solar eclipse observed" << endl;
    // cout << "=================================" << endl;
    // cout << endl;
    }
    }
     
    if (SolarEclipse_angle > 175)
    {
    if ( SolarEclipse_angle < 202){Front = -1;}
    }
    // cout << endl;
    // cout << "=================" << endl;
    // cout << "Front = " << Front << endl;
    // cout << "=================" << endl;
    // cout << endl;
    }
     
    void get_calculatedMoon_data(double ReqMoon, double& CalculatedMoon, double& solar_days,
    double& moon_x12_2, double& moon_minus_solar_days,
    int& MoonStatus, double& MoonLeadLag, double& moon_add_minus)
    {
    double moon_84_check;
    if (ReqMoon > half_Zhuanzhong){CalculatedMoon = ReqMoon - half_Zhuanzhong;MoonStatus = LagStart;}
    if (ReqMoon < half_Zhuanzhong){CalculatedMoon = ReqMoon;MoonStatus = LeadStart;}
     
    moon_x12_2 = CalculatedMoon * 12.2;
    solar_days = 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;
    }
    MoonLeadLag = moon_84_check;
     
    moon_minus_solar_days = CalculatedMoon - solar_days;
     
    if( moon_minus_solar_days < 0 )
    {
    solar_days = solar_days - 0.082008;
    moon_minus_solar_days = CalculatedMoon - solar_days;
     
    if (fmod(MoonStatus, 2.0) == 0){MoonLeadLag = int(MoonLeadLag);}
    if (fmod(MoonStatus, 2.0) > 0){MoonLeadLag = 167 - int(moon_x12_2);}
    }
    moon_add_minus = 11.081575 -(0.05815 + 0.00195 *(MoonLeadLag-1)/2)* MoonLeadLag;
    }
     
     
    void getMoon_adjusted_data(int MoonStatus, double MoonLeadLag, double moon_add_minus,
    double moon_minus_solar_days,
    double& moon_sum, double& MoonX, double& MoonSpeedFactor)
    {
    moon_sum = 0, MoonX = 0, MoonSpeedFactor = 0;
     
    switch(MoonStatus)
    {
    case LagStart:
    MoonSpeedFactor = 0.9855 +(0.05815 + 0.00195 *(MoonLeadLag-1)/2)* MoonLeadLag/100;
    moon_sum =(0.1111 -(0.000281 + 0.00000325 *MoonLeadLag)*MoonLeadLag)*MoonLeadLag;
    MoonX = moon_sum + moon_minus_solar_days * moon_add_minus/820;
     
    break;
     
    case LagEnd:
    MoonSpeedFactor = 1.2071 -(0.05815 + 0.00195 *(MoonLeadLag-1)/2)*MoonLeadLag/100;
    MoonLeadLag = MoonLeadLag + 1;
    moon_sum =(0.1111 -(0.000281 + 0.00000325 * MoonLeadLag)*MoonLeadLag)*MoonLeadLag;
    MoonX = moon_sum - moon_minus_solar_days * moon_add_minus/820;
     
    break;
     
    case LeadStart:
    MoonSpeedFactor = 1.2071 -(0.05815 + 0.00195 *(MoonLeadLag-1)/2)*MoonLeadLag/100;
    moon_sum =(0.1111 -(0.000281 + 0.00000325 *MoonLeadLag)*MoonLeadLag)*MoonLeadLag;
    MoonX = moon_sum + moon_minus_solar_days * moon_add_minus/820;
     
    break;
     
    case LeadEnd:
    MoonSpeedFactor = 0.9855 +(0.05815 + 0.00195*(MoonLeadLag-1)/2)*MoonLeadLag/100;
    MoonLeadLag = MoonLeadLag + 1;
    moon_sum =(0.1111 -(0.000281 + 0.00000325 * MoonLeadLag)*MoonLeadLag)*MoonLeadLag;
    MoonX = moon_sum - moon_minus_solar_days * moon_add_minus/820;
     
    break;
    default:
    cout << " no MoonStatus";
     
    break;
    }
    }
     
    };
     
    class base2 {
    protected:
    int k;
    public:
    base2(int x) {k=x;}
    ~base2() {}
    };
     
    class derived: public Ref_data, public base2 {int j;
    public:
    derived(int m, int n, int x, int z): Ref_data(m,n), base2(x){j=z;}
    ~derived() {}
     
    void show() {cout << calculated_year << " " << k << " " << j << "\n";}
    };
     
    int main()
    {
    cout << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt" << "\tt"
    << "\tt" << "\tt" << "\tt" << endl;
     
    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\n"
    << "year(month)\n";
     
    int index = 0;int Pindex = 0;int calculated_year, total_months;
     
    for (calculated_year = start_year; calculated_year - 1< end_year; calculated_year ++)
    {
    for (total_months = 0; total_months < 14; total_months ++)
    {
    Ref_data year(calculated_year,total_months);
     
    double total_days, Winter_solitice, LeapDay;
    year.get_basic(total_days, Winter_solitice, LeapDay);
     
    double BaseSun, RefEcl, RefSun, RefMoonDay, RefMoon, MinSun;
    year.get_ref(total_days,Winter_solitice,LeapDay,BaseSun,RefEcl,RefSun,RefMoonDay,RefMoon,MinSun);
     
    double ReqEcl, ReqEcl_angle, ReqSun, ReqMoonDay, ReqMoon;
    year.get_required(RefEcl,RefSun,RefMoonDay,RefMoon,ReqEcl,ReqSun,ReqMoonDay,ReqMoon);
     
    ReqEcl_angle = ReqEcl * MoonMeanSpeed;
     
    double ReqSun_days = ReqSun;
    if (ReqSun > half_year){ReqSun_days = ReqSun - half_year;};
     
    double sun_days;int SunStatus, intSun_days;
    year.checkSunStatus(ReqSun,sun_days, SunStatus);
     
    intSun_days = int(sun_days);
     
    double sun_add_minus, sun_sum, SunX;
    year.calculateSunOffset(SunStatus, sun_days, intSun_days, sun_add_minus, sun_sum, SunX);
     
    //-----------------
    double ReqSun360 = fmod(ReqSun * 360/Tropical_year + 180, 360);
    double e = 0.0167;
    double SunX_Th = 180/pi *(2 *e* sin(ReqSun360 *pi/180)+ 1.25 * e*e * sin(2*ReqSun360 *pi/180));
    double SFactor = MoonMeanSpeed * 0.082 + 0.082; // = 1.1782375
     
    double SunX360 = SunX_Th * SFactor;
    SunX = Tropical_year/360* SunX360;
    //-----------------
    double MinSun360 = MinSun * 360/Tropical_year;
    double RSunYel = ReqSun360 + SunX_Th + MinSun360 + 0.098;
     
    // 求正交平均
    double PosExtra = -570 *SunX_Th/6973;
    double PosX360 = PosExtra * SFactor;
    double PosX = PosX360 * Tropical_year/360;
     
    //^^^^^^^^^^^^^^(Convert Shoushi ReqEcl as SunPos)^^^^^^^^^^
    double SunPos = (ReqEcl - Jiaozhong_C/2)*360/Jiaozhong_C - PosExtra + SunX_Th;
    if (SunPos < 0){SunPos = SunPos + 360;}
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     
    double UsePos = RSunYel - SunPos;
     
    // 求正交實均
    double SunPosReal = 180/pi *atan(56.0/59 *tan(SunPos * pi/180));
    if ((SunPos > 90)&(SunPos < 270)){SunPosReal = SunPosReal + 180;}
    if (SunPos > 270){SunPosReal = SunPosReal + 360;}
     
    // 正交實均
    double PosRealAdd = SunPos - SunPosReal;
     
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    double PosRX360 = PosRealAdd * SFactor;
    double PosRX = PosRX360 * Tropical_year/360;
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     
    // 求正交實行
    double PosReal = UsePos + PosRealAdd;
    double MoonPos = RSunYel - PosReal; // same as MoonPos = SunPos - PosRealAdd
    if (MoonPos < 0){MoonPos = MoonPos + 360;}
    //----------
    double YelWh = 5 + 17.0/60 + 20.0/3600;
    double YMoonPos = 180/pi * atan(cos(YelWh * pi/180)* tan(MoonPos * pi/180));
    if ((MoonPos > 90)&(MoonPos < 270)){YMoonPos = YMoonPos + 180;}
     
    double YelD = YMoonPos - MoonPos;
    if(YelD < -350){YelD = YelD + 360;}
    //------
    double SolarEclipse_angle;int Front = 0;
    year.get_SolarEclipse_angle(SunStatus,ReqEcl_angle,SunX, SolarEclipse_angle, Front);
    if (SolarEclipse_angle < 7){SolarEclipse_angle = SolarEclipse_angle + 363.7934196;}
     
    double SolEA = SolarEclipse_angle - PosX - PosRX;
    // if(NS_f < 0){SolEA = SolarEclipse_angle - NS_Hu/60;}
     
    double SolEA360 = fmod(SolEA + 363.7934196/2, 363.7934196)* 360/363.7934196;
     
    double LatYel = 180/pi *asin(sin(YelWh *pi/180)*sin(MoonPos *pi/180));
    LatYel = 180/pi *asin(sin(YelWh *pi/180)*sin(SolEA360 *pi/180));
     
    //RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
    // Convert Shoushi ReqMoon to consider Ee_o
    //---------------------------------------------------------------
    double OReqMoon = ReqMoon;
    double OReqMoon360 = fmod(OReqMoon * 360/Zhuanzhong_C + 180, 360);
    if (OReqMoon360 < 0){OReqMoon360 = fmod(360 + OReqMoon360, 360);}
    //------------------------------
    // apply Earth Center Offset for Moon adjustment
    double Ee_o = 0;
    if (cos(OReqMoon360 * pi/180)< 0){Ee_o = 0;}
    double y_scale = 1;
    double x_cir = cos(OReqMoon360 * pi/180);
    double y_cir = sin(OReqMoon360 * pi/180);
    double r_elip = 1/sqrt(x_cir*x_cir + y_cir*y_cir/(y_scale*y_scale));
    double x_o = r_elip * x_cir;
    double y_o = r_elip * y_cir;
     
    double EMoon360 = 180/pi *atan(y_o/(x_o - Ee_o));
    if ((OReqMoon360 > 90) & (OReqMoon360 < 270)){EMoon360 = 180 + EMoon360;}
    if (OReqMoon360 > 270){EMoon360 = 360 + EMoon360;}
    if ((OReqMoon360 < 90) & (EMoon360 < 0)){EMoon360 = 180 + EMoon360;}
    //------------------------------
    double MOff = EMoon360 - OReqMoon360;
    double EMoon = fmod((EMoon360 + 180)*Zhuanzhong_C/360, Zhuanzhong_C);
    double MF = 0;
    ReqMoon = EMoon;
    double ReqMoon360 = fmod(ReqMoon * 360/Zhuanzhong_C + 180, 360);
    //RRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
     
    int MoonStatus = 0;double CalculatedMoon = 0, MoonLeadLag, moon_add_minus, solar_days;
    double moon_x12_2, moon_minus_solar_days;
    year.get_calculatedMoon_data(ReqMoon,CalculatedMoon,solar_days,moon_x12_2, moon_minus_solar_days,
    MoonStatus, MoonLeadLag, moon_add_minus);
     
    double moon_sum, MoonX, MoonSpeedFactor;
    year.getMoon_adjusted_data(MoonStatus, MoonLeadLag, moon_add_minus, moon_minus_solar_days,
    moon_sum, MoonX, MoonSpeedFactor);
    //--------------
    double SunXS = -abs(SunX);
    if (SunStatus > LagEnd){SunXS = abs(SunX);}
     
    double MoonXS = -abs(MoonX);
    if (MoonStatus > LagEnd){MoonXS = abs(MoonX);}
     
    double SunMoonX = SunXS - MoonXS;
    double SunMoonT = SunMoonX * 0.082/MoonSpeedFactor;
    //--------------
    double RSunEnd = RSunYel + SunMoonT * 360/Tropical_year;
     
    double MoonT24 = SunMoonT * 24;
     
    //=========================================
    // Use HouBian steps to simulate MidMoon
    //=========================================
     
    const double Qiying_HB = 32.12254;
    double SumDay_HB = (calculated_year - 1723)* 365.24233442;
    double Winter_HB = fmod((SumDay_HB + Qiying_HB +(int)(abs(SumDay_HB))*60),60);
    double UseDays = SumDay_HB +(Qiying_HB -(int)Qiying_HB)-(Winter_HB -(int)Winter_HB);
     
    double DayMoon = fmod(UseDays* 47435.0234086, 1296000)/3600;
    double MoonYearBase = DayMoon + 176 + 27.0/60 + 48.0/3600 + 53.0/21600;
     
    const double Shuoying_HB = 15.12633;
    const double Shuoshi_HB = 29.53059053;
    double days_for_moon = UseDays - Shuoying_HB;
    double Remainder_moon = fmod(days_for_moon, Shuoshi_HB);
    if (Remainder_moon < 0){Remainder_moon = Shuoshi_HB + Remainder_moon;}
     
    double SumMonth_HB = total_months - 1;
    if (SumMonth_HB < 0){SumMonth_HB = 0;}
    if (calculated_year == 1719){SumMonth_HB = SumMonth_HB - 1;}
    if (calculated_year == 1643){SumMonth_HB = SumMonth_HB - 1;}
     
    double Next_moon_days = Shuoshi_HB - Remainder_moon;
    double MoonDay = (Next_moon_days + SumMonth_HB * Shuoshi_HB)*47435.0234086/3600;
    double MidMoon = fmod(MoonYearBase + MoonDay, 360);
     
    //=========================================
    // Use HouBian steps to simulate MaxMoon
    //=========================================
    // 求積日最高平行
    double BaseMaxMoon = fmod(UseDays * 401.070226, 1296000)/3600;
    // 求最高年根
    double MaxYearBase = BaseMaxMoon +(241 + 15.0/60 + 45.0/3600 + 38.0/216000);
    // 求最高日數
    double MaxDay = (Next_moon_days + SumMonth_HB * Shuoshi_HB)* 401.070226/3600;
    // 求最高平行
    double MaxMoon =fmod(MaxYearBase + MaxDay, 360);
    if(MaxMoon < 0){MaxMoon = 360 + MaxMoon;}
    //-----------------------------------------
    double ReqMoon_HB = MidMoon - MaxMoon;
    if (ReqMoon_HB < 0){ReqMoon_HB = 360 + ReqMoon_HB;}
    //-----------------------------------------
     
    double MoonX360 = MoonXS * 360/Tropical_year;
    double MoonX_Th = MoonX360 + SunX_Th - SunX360;
    double UseRMoon = MidMoon + MoonX_Th;
    double RMoonEnd = fmod(UseRMoon + SunMoonT *(MoonSpeedFactor/0.082 + 1)*360/Tropical_year, 360);
    //------------
    double SunMoonXS = SunMoonX;
    int SunStatusS = SunStatus;
    int MoonStatusS = MoonStatus;
    //------------
     
    double MDay0 = ReqMoonDay - (int)ReqMoonDay;
    double MDay0_24 = MDay0 * 24;
     
    double MDayT = MDay0 + SunMoonT;
    if(MDayT < 0){MDayT = 1 + MDayT;}
    if(MDayT > 1){MDayT = MDayT - 1;}
    double MDayT_day = int(MDayT);
    double MDayT_hour = MDayT - int(MDayT);
    double MDayT_24hour = MDayT_hour * 24;
     
    //****************************
    // Calculate half day hours
     
    double sun_rise_time
    = 2068.3 +(((0.122928 -((intSun_days-53)* 0.000608))*(intSun_days-1))/2)*intSun_days + 0.06 * intSun_days;
     
    double sun_set_time = 10000 - sun_rise_time;
     
    double half_day_hours = sun_set_time - 5000;
    if (SunStatus == LagEnd){half_day_hours = 5000 - half_day_hours;}
    if (SunStatus == LeadStart){half_day_hours = 5000 - half_day_hours;}
    //*****************
    //----------------------
    double time_from_noon = (MDayT - int(MDayT)) - 0.5;
    double observed_timeS = (0.5 - abs(time_from_noon))* (time_from_noon)/0.96;
    //----------------------
     
    int SeeTimeTable[][19]
    = {
    //h_deg-120 -105 -90 -75 -60 -45 -30 -15 0 15 30 45 60 75 90 105 120 <--
    //Gon h20 h19 h18 h17 h16 h15 h14 h13 h12 h11 h10 h9 h8 h7 h6 h5 h4 Gon
    //--------------------------------------------------------------------------------------------
    {9, -61, -72, -81, -85, -85, -78, -61, -34, -0, 34, 61, 78, 85, 85, 81, 72, 61, 9},
    {10, -51, -61, -68, -70, -66, -54, -32, -2, 33, 63, 84, 95, 99, 96, 90, 80, 67, 8},
    {11, -40, -49, -55, -55, -50, -37, -14, 15, 48, 75, 93, 103, 106, 103, 96, 86, 74, 7},
    {0, -33, -42, -48, -50, -46, -35, -15, 13, 44, 72, 92, 103, 108, 106, 100, 91, 80, 6},
    {1, -34, -44, -52, -55, -54, -45, -28, -3, 29, 60, 84, 98, 105, 105, 100, 92, 81, 5},
    {2, -45, -56, -64, -69, -69, -62, -46, -20, 13, 46, 73, 90, 97, 98, 94, 86, 75, 4},
    {3, -61, -72, -81, -85, -85, -78, -61, -34, 0, 34, 61, 78, 85, 85, 81, 72, 61, 3},
    //---------------------------------------------------------------------------------------------
    //Gon h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16 h17 h18 h19 h20 Gon
    };
    //-->
     
    double degree = fmod(RSunEnd + 270, 360);
    double Gon = degree/30.0;
    double Hour = MDayT_hour * 24;
     
    int LRow = 9;
    int RRow = 9;
    if (fmod(Gon + 3,12.0) < 6){LRow = (int)(fmod((int)Gon + 3, 12.0));}
    if (fmod(Gon + 3,12.0) > 5){RRow = (int)(fmod(9 - (int)Gon, 12.0));}
     
    int TableTime = 0;
    int LHour = 0;
    double t0 = 0;
    double t1 = 0;
    double SeeT = 0;
     
    if (LRow < 6)
    {
    LHour = (int)Hour - 3;
    TableTime = SeeTimeTable[LRow][LHour];
    t0 = TableTime +(SeeTimeTable[LRow][LHour+1] - TableTime)*(Hour -(int)Hour);
     
    TableTime = SeeTimeTable[LRow+1][LHour];
    t1 = TableTime +(SeeTimeTable[LRow+1][LHour+1] - TableTime)*(Hour -(int)Hour);
    }
     
    int RHour = 0;
    if (LRow > 5)
    {
    if (RRow < 7)
    {
    RHour = 21 - (int)Hour;
    TableTime = SeeTimeTable[RRow][RHour];
    t0 = TableTime +(SeeTimeTable[RRow][RHour-1] - TableTime)*(Hour -(int)Hour);
    t0 = -t0;
     
    TableTime = SeeTimeTable[RRow-1][RHour];
    t1 = TableTime +(SeeTimeTable[RRow-1][RHour-1] - TableTime)*(Hour -(int)Hour);
    t1 = -t1;
    }
    }
    SeeT = t0 + (t1 - t0)*(Gon -(int)Gon);
    double SeeHour = SeeT/60;
    double observed_time = SeeHour/24;
    double eclipseMax_hour = MDayT_hour + observed_time;
     
    //-------------------
    double eclipseMax_24hour = eclipseMax_hour * 24;
     
    // cout << "===============================================" << endl;
    // cout << " eclipseMax_24hour = " << eclipseMax_24hour << endl;
    // cout << "===============================================" << endl;
    // cout << endl;
    //********************
     
    double time_difference_from_noon = time_from_noon + observed_time;
    time_difference_from_noon = abs(time_difference_from_noon)* 10000;
     
    double MaxPhase_lead_lag = ReqSun + int(MDayT) + eclipseMax_hour;
     
    double MaxPhaseSun_ref = MaxPhase_lead_lag - ReqMoonDay;
     
    double MaxPhaseSun_days;
    year.checkSunStatus(MaxPhaseSun_ref, MaxPhaseSun_days, SunStatus);
     
    intSun_days = int(MaxPhaseSun_days);
     
    double MaxPhaseSun_add_minus, MaxPhaseSun_sum, MaxPhaseSunX;
    year.calculateSunOffset(SunStatus, MaxPhaseSun_days, intSun_days,MaxPhaseSun_add_minus, MaxPhaseSun_sum, MaxPhaseSunX);
     
    //*************
    double MaxPhaseSun_adjust = 0;
    if (SunStatus < LeadStart){MaxPhaseSun_adjust = MaxPhaseSun_ref - MaxPhaseSunX;}
    if (SunStatus > LagEnd){MaxPhaseSun_adjust = MaxPhaseSun_ref + MaxPhaseSunX;}
     
    //*************
    double MaxPhaseSun_NS;
    year.checkSunStatus(MaxPhaseSun_adjust, MaxPhaseSun_NS, SunStatus);
     
    double NS_difference = 4.46 -((MaxPhaseSun_NS * MaxPhaseSun_NS)/1870);
    double NS_Sh = NS_difference -(NS_difference * time_difference_from_noon)/half_day_hours;
     
    double MaxPhaseSun_EW = MaxPhaseSun_adjust;
    if(MaxPhaseSun_EW > half_year){MaxPhaseSun_EW = MaxPhaseSun_EW - half_year;}
     
    double EW_difference = ((half_year - MaxPhaseSun_EW)* MaxPhaseSun_EW)/1870;
     
    double EW_level = (EW_difference * time_difference_from_noon)/2500;
    double EW_Sh = EW_level;
    if(EW_level > EW_difference ){EW_Sh = 2 * EW_difference - EW_level;}
     
    int AfterNoonSign;
    AfterNoonSign = 1;
    if (time_from_noon < 0){AfterNoonSign = -1;}
     
    switch(SunStatus)
    {
    case LagStart:
    NS_Sh = Front * NS_Sh;
    EW_Sh = -AfterNoonSign * Front * EW_Sh;
    break;
     
    case LagEnd:
    NS_Sh = -Front * NS_Sh;
    EW_Sh = -AfterNoonSign * Front * EW_Sh;
    break;
     
    case LeadStart:
    NS_Sh = -Front * NS_Sh;
    EW_Sh = AfterNoonSign * Front * EW_Sh;
    break;
     
    case LeadEnd:
    NS_Sh = Front * NS_Sh;
    EW_Sh = AfterNoonSign * Front * EW_Sh;
    break;
    default:
    cout << " no SunStatus" << endl;
    }
    //---------------------------
    double eclipse_ref_angle = 0;
    if (Front == -1){eclipse_ref_angle = 188.05 + NS_Sh + EW_Sh;}
    if (Front == 1){eclipse_ref_angle = 357.64 + NS_Sh + EW_Sh;}
     
    double postiveEclipse_angle = SolarEclipse_angle - eclipse_ref_angle;
    double SexEcl = Front * postiveEclipse_angle;
     
    double EclipseFracS = 0;
    double EclipseRef = 6;
    if (SexEcl < 0){EclipseRef = 8;}
     
    EclipseFracS = (EclipseRef - abs(postiveEclipse_angle))/EclipseRef;
    //**************
     
    int Lat_LongTable[][19]
    = {
    //h_deg-120 -105 -90 -75 -60 -45 -30 -15 0 15 30 45 60 75 90 105 120 <--
    //Gon h20 h19 h18 h17 h16 h15 h14 h13 h12 h11 h10 h9 h8 h7 h6 h5 h4 Gon
    //--------------------------------------------------------------------------------------------
    {9, 10, 15, 20, 25, 27, 30, 35, 40, 42, 40, 35, 30, 27, 25, 20, 15, 10, 9},//Lat(min)
    {9, 0, 2, 0, 47, 47, 41, 28, 26, 25, 27, 28, 41, 47, 47, 0, 2, 1, 9},//Lat(sec)
    {9, 30, 35, 40, 41, 42, 39, 31, 17, 0, 17, 31, 39, 42, 41, 40, 35, 30, 9},//Long(min)
    {9, 0, 0, 0, 20, 17, 30, 58, 6, 0, 6, 58, 30, 17, 21, 0, 0, 0, 9},//Long(sec)
     
    {10, 32, 33, 34, 35, 36, 40, 41, 41, 38, 32, 26, 20, 16, 14, 14, 13, 12, 8},
    {10, 0, 0, 40, 33, 38, 26, 29, 24, 25, 36, 40, 51, 48, 39, 0, 0, 0, 8},
    {10, 0, 0, 0, 35, 33, 30, 19, 3, 13, 30, 40, 46, 46, 46, 47, 48, 49, 8},
    {10, 0, 0, 0, 7, 27, 8, 49, 35, 10, 43, 37, 30, 42, 37, 0, 0, 0, 8},
     
    {11, 36, 37, 38, 39, 40, 41, 39, 36, 29, 21, 16, 11, 9, 8, 7, 6, 5, 7},
    {11, 0, 0, 26, 37, 39, 29, 34, 30, 29, 44, 44, 53, 57, 48, 51, 0, 0, 7},
    {11, 32, 31, 30, 30, 29, 23, 12, 1, 21, 35, 43, 47, 49, 49, 48, 47, 46, 7},
    {11, 0, 0, 13, 14, 16, 22, 29, 27, 15, 2, 1, 50, 44, 45, 47, 0, 0, 7},
     
    {0, 39, 40, 41, 41, 40, 38, 35, 28, 21, 13, 9, 7, 6, 5, 5, 5, 4, 6},
    {0, 0, 0, 27, 28, 30, 31, 25, 36, 41, 56, 53, 55, 54, 57, 53, 0, 0, 6},
    {0, 28, 29, 30, 30, 28, 21, 13, 0, 15, 32, 41, 47, 49, 49, 48, 48, 47, 6},
    {0, 0, 0, 3, 5, 4, 11, 13, 36, 32, 7, 57, 51, 44, 47, 47, 0, 0, 6},
     
    {1, 38, 39, 41, 40, 38, 33, 27, 21, 13, 8, 7, 6, 8, 8, 9, 10, 11, 5},
    {1, 0, 47, 26, 28, 21, 30, 36, 40, 55, 46, 52, 56, 48, 52, 27, 40, 0, 5},
    {1, 26, 27, 28, 30, 31, 27, 19, 7, 10, 26, 37, 45, 49, 49, 48, 48, 47, 5},
    {1, 0, 33, 26, 15, 9, 7, 11, 13, 23, 11, 5, 49, 46, 57, 53, 13, 0, 5},
     
    {2, 35, 35, 36, 34, 30, 25, 18, 10, 9, 6, 7, 10, 13, 16, 17, 16, 17, 4},
    {2, 17, 7, 33, 31, 30, 39, 54, 18, 50, 57, 52, 48, 47, 10, 11, 46, 0, 4},
    {2, 27, 28, 29, 31, 34, 34, 27, 13, 3, 18, 34, 43, 46, 47, 46, 45, 45, 4},
    {2, 41, 53, 52, 54, 36, 12, 22, 29, 20, 20, 38, 22, 27, 27, 34, 34, 0, 4},
     
    {3, 26, 27, 27, 26, 21, 17, 13, 8, 6, 8, 13, 17, 21, 26, 27, 27, 26, 3},
    {3, 40, 15, 55, 42, 48, 45, 45, 52, 56, 52, 45, 45, 48, 42, 55, 15, 40, 3},
    {3, 39, 40, 41, 42, 43, 39, 31, 18, 0, 18, 31, 39, 43, 42, 41, 40, 39, 3},
    {3, 13, 20, 27, 29, 35, 29, 37, 21, 0, 38, 37, 29, 25, 39, 27, 20, 13, 3},
     
    //---------------------------------------------------------------------------------------------
    //Gon h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16 h17 h18 h19 h20 Gon
    //-->
    };
     
    enum TRow{Lat_min, Lat_sec, Long_min, Long_sec};
     
    LRow = 9;
    RRow = 9;
    if (fmod(Gon + 3,12.0) < 6){LRow = (int)(fmod((int)Gon + 3, 12.0));}
    if (fmod(Gon + 3,12.0) > 5){RRow = (int)(fmod(9 - (int)Gon, 12.0));}
     
    double NS_0 = 0;
    double NS_1 = 0;
    double NS_01 = 0;
    double NS_N0 = 0;
    double NS_N1 = 0;
    double NS_N01 = 0;
     
    double EW_0 = 0;
    double EW_1 = 0;
    double EW_01 = 0;
    double EW_N0 = 0;
    double EW_N1 = 0;
    double EW_N01 = 0;
     
    LHour = 0;
    double NS_Hu = 0;
    double EW_Hu = 0;
     
    if (LRow < 6)
    {
    LHour = (int)Hour - 3;
    NS_0 = Lat_LongTable[LRow *4+Lat_min][LHour] + Lat_LongTable[LRow*4+Lat_sec][LHour]/60.0;
    NS_1 = Lat_LongTable[LRow *4+Lat_min][LHour+1] + Lat_LongTable[LRow*4+Lat_sec][LHour+1]/60.0;
    NS_01 = NS_0 + (NS_1 - NS_0)*(Hour -(int)Hour);
     
    NS_N0 = Lat_LongTable[(LRow+1)*4+Lat_min][LHour] + Lat_LongTable[(LRow+1)*4+Lat_sec][LHour]/60.0;
    NS_N1 = Lat_LongTable[(LRow+1)*4+Lat_min][LHour+1] + Lat_LongTable[(LRow+1)*4+Lat_sec][LHour+1]/60.0;
    NS_N01 = NS_N0 + (NS_N1 - NS_N0)*(Hour -(int)Hour);
    //--------------
    EW_0 = Lat_LongTable[LRow *4+Long_min][LHour] + Lat_LongTable[LRow*4+Long_sec][LHour]/60.0;
    EW_1 = Lat_LongTable[LRow *4+Long_min][LHour+1] + Lat_LongTable[LRow*4+Long_sec][LHour+1]/60.0;
    EW_01 = EW_0 + (EW_1 - EW_0)*(Hour -(int)Hour);
     
    EW_N0 = Lat_LongTable[(LRow+1)*4+Long_min][LHour] + Lat_LongTable[(LRow+1)*4+Long_sec][LHour]/60.0;
    EW_N1 = Lat_LongTable[(LRow+1)*4+Long_min][LHour+1] + Lat_LongTable[(LRow+1)*4+Long_sec][LHour+1]/60.0;
    EW_N01 = EW_N0 + (EW_N1 - EW_N0)*(Hour -(int)Hour);
    }
     
    RHour = 0;
    if (LRow > 5)
    {
    if (RRow < 7)
    {
    RHour = 21 - (int)Hour;
    NS_0 = Lat_LongTable[RRow *4+Lat_min][RHour] + Lat_LongTable[RRow*4+Lat_sec][RHour]/60.0;
    NS_1 = Lat_LongTable[RRow *4+Lat_min][RHour-1] + Lat_LongTable[RRow*4+Lat_sec][RHour-1]/60.0;
    NS_01 = NS_0 + (NS_1 - NS_0)*(Hour -(int)Hour);
     
    NS_N0 = Lat_LongTable[(RRow-1)*4+Lat_min][RHour] + Lat_LongTable[(RRow-1)*4+Lat_sec][RHour]/60.0;
    NS_N1 = Lat_LongTable[(RRow-1)*4+Lat_min][RHour-1] + Lat_LongTable[(RRow-1)*4+Lat_sec][RHour-1]/60.0;
    NS_N01 = NS_N0 + (NS_N1 - NS_N0)*(Hour -(int)Hour);
    //---------------------------------
    EW_0 = Lat_LongTable[RRow *4+Long_min][RHour] + Lat_LongTable[RRow*4+Long_sec][RHour]/60.0;
    EW_1 = Lat_LongTable[RRow *4+Long_min][RHour-1] + Lat_LongTable[RRow*4+Long_sec][RHour-1]/60.0;
    EW_01 = EW_0 + (EW_1 - EW_0)*(Hour -(int)Hour);
     
    EW_N0 = Lat_LongTable[(RRow-1)*4+Long_min][RHour] + Lat_LongTable[(RRow-1)*4+Long_sec][RHour]/60.0;
    EW_N1 = Lat_LongTable[(RRow-1)*4+Long_min][RHour-1] + Lat_LongTable[(RRow-1)*4+Long_sec][RHour+1]/60.0;
    EW_N01 = EW_N0 + (EW_N1 - EW_N0)*(Hour -(int)Hour);
    }
    }
    NS_Hu = NS_01 + (NS_N01 - NS_01)*(Gon -(int)Gon);
    EW_Hu = EW_01 + (EW_N01 - EW_01)*(Gon -(int)Gon);
    //=============
     
    // 求最高平均
    double MaxExtra = 1196 * SunX_Th/6973;
    // 求用最高
    double UseMax = MaxMoon + MaxExtra;
     
    // 求日距月最高
    double MaxSMoon = RSunYel - UseMax;
    if (MaxSMoon < 0){MaxSMoon = MaxSMoon + 360;}
     
    //mmmmmmmmmm
    // 所夾之角
    double UseAngle = 180 - fmod(MaxSMoon *2, 360);
    // 最高實均
    double MaxRealExtra
    = 180/pi *atan((117315 *sin(UseAngle *pi/180))/(550505 - 117315 *cos(UseAngle *pi/180)));
    //mmmmmmmmmm
    //==========
    // From HouBian
    // ------------
    // 求本天心距地數(本時两心差)
    double C_Earth = 117315.0/1000000 *sin(MaxSMoon *2*pi/180)/sin(MaxRealExtra *pi/180);
    //=========
    // From HouBian Step 3 :
    // 推地平高下差及日月視徑第三
     
    // 求太陰距地
    double DistCos = cos(ReqMoon360 * pi/180)* 2 * C_Earth;
    double DistSin = sin(ReqMoon360 * pi/180)* 2 * C_Earth;
    double AddDistCos = 20 - DistCos;
    double DistRef = DistSin * DistSin/AddDistCos;
    double CalRef = (DistRef + AddDistCos)/2;
     
    // 月距地心數(太陰距地)
    double MoonToEarth = 20 - CalRef;
     
    //求日距地心數
    // 太陽實引
    double ReqRealSun = fmod(ReqSun360 + SunX360, 360);
    // 分股
    double SideCos = cos(ReqRealSun * pi/180)* 0.338;
    // 勾
    double SideSin = sin(ReqRealSun * pi/180)* 0.338;
    // 勾弦和
    double AddSideCos = 20 + SideCos;
    // 勾弦較
    double SideRef = SideSin * SideSin/AddSideCos;
    // 弦
    double BaseRef =(SideRef + AddSideCos)/2;
    // 日距地心數
    double SunToEarth = 20 - BaseRef;
     
    //=======
    // From HouBian Step 2 :
     
    // 一小時太陰白道實行
    double HourWhiteMoon = MoonSpeedFactor/0.082 + 1; // simulate HouBian
    double HourRSunYel = 1; // simulate HouBian
     
    // 求斜距交角差
    double XA = 180/pi*atan(HourRSunYel*sin(YelWh*pi/180)/(HourWhiteMoon-HourRSunYel*cos(YelWh*pi/180)));
     
    // 求斜距黃道交角
    double XYelA = YelWh + XA;
    // 求食甚實緯
    double MaxLat = LatYel * cos(XYelA *pi/180);
    //^^^^^^^^^^
    double MaxMoonX =(3450.0 * 10/MoonToEarth)/60;
    double HighLow = MaxMoonX - 10.0/60;
     
    double SunSeeSemiDiameter =(966 * 10/SunToEarth)/60;
    double SunRSemiDiameter = SunSeeSemiDiameter - 15.0/60;
    double MoonSeeSemiDiameter =(940.5 * 10/MoonToEarth)/60;
    double NetDiameter = SunRSemiDiameter + MoonSeeSemiDiameter;
    //^^^^^^^^^
    double MaxSunYel = RSunEnd + observed_time/24;
    double MaxPhaseUseT = eclipseMax_hour;
    //=========
    // Step 4 : 推食甚太陽黃赤經緯宿度及黃赤二經交角第四
    //-----------------
    // 求食甚太陽赤道經度
    double SprYel = MaxSunYel - 90;
    double WinYel = MaxSunYel;
    double MaxSunRed = 180/pi*atan(cos((23+29.0/60)*pi/180)*tan(SprYel*pi/180));
     
    double SprRed = MaxSunRed;
    if ((SprYel > 90)&(SprYel < 270)){SprRed = 180 + MaxSunRed;}
    if (SprYel > 270){SprRed = 360 - MaxSunRed;}
    //-----------------
    // 求食甚太陽赤道緯度
    double ESunLat = 180/pi*asin(sin((23+29.0/60)*pi/180)*sin(SprYel*pi/180));
    // 求太陽距北極
    double SunNPole = 90 - ESunLat;
     
    // 求黃赤二經交角
    double YelRedXA = 90-180/pi*atan(tan((90-(23 + 29.0/60))*pi/180)/cos(SprYel*pi/180));
    if (cos(SprYel*pi/180)< 0){YelRedXA=YelRedXA-180;}
     
    // 求黄白二經交角
    double YelWhXA = XYelA;
    if (cos(MoonPos * pi/180) < 0){YelWhXA = -XYelA;}
     
    // 求赤白二經交角
    double RedWhXA = YelRedXA + YelWhXA;
     
    //=========
    // Step 5 : 又法 page(792-248)推食甚近時第五
     
    // 求用時太陽距午赤道度
    double T_12Red =(MaxPhaseUseT -(int)MaxPhaseUseT - 0.5)* 24 * 15;
    //-----------------
    // 求用時赤經高弧交角
     
    // 北極距天頂 = 90 -(Beijing Lat)
    double NPoleV = 50.08333333;
    // 距極分邊
    double T_12X=180/pi*atan(cos(T_12Red*pi/180)*tan(NPoleV*pi/180));
    // 距日分邊
    double T_12NP = SunNPole - T_12X;
    // 垂弧之正切線
    double T_12Arc = tan(T_12Red*pi/180)*sin(T_12X*pi/180);
    // 用時赤經高弧交角
    double URedXA = 180/pi*atan(T_12Arc/sin(T_12NP*pi/180));
    //---------------
    // 求用時太陽距天頂
    double UseSunHL = 180/pi*asin(sin(NPoleV*pi/180)*sin(T_12Red*pi/180)/sin(URedXA*pi/180));
     
    // 求用時白經高弧交角
    double UWhXA = URedXA + RedWhXA;
    // 求用時高下差
    double UseHL = HighLow*sin(UseSunHL*pi/180);
    // 求用時東西差
    double UseEW = sin(UWhXA*pi/180)*UseHL;
    double EW_f = UseEW/60;
    // 求用時南北差
    double UseNS = -cos(UWhXA*pi/180)*UseHL;
    double NS_f = UseNS/60;
     
    // 求用時視緯
    double UseSeeLat = MaxLat*60 + UseNS;
    //^^^^^^^^^^^^^^^
    double ESeeLat = UseSeeLat/60;
    double SunHalfD = SunRSemiDiameter/60;
    double addD = NetDiameter/60;
     
    double eclipseFrac =(addD - abs(ESeeLat))/(2* SunHalfD);
    //***************
    double RefEA_Sh = 0;
    if (Front == -1){RefEA_Sh = 188.05;}
    if (Front == 1){RefEA_Sh = 357.64;}
    double TestEA = SolarEclipse_angle - RefEA_Sh - PosX - PosRX;
    double NetEA = TestEA -(NS_Sh + EW_Sh);
    double SexEA = NetEA;
    if (Front == -1){SexEA = -1.0 * NetEA;}
     
    double EARef = 6;
    if (SexEA < 0){EARef = 8;}
    EclipseFracS = (EARef - abs(SexEA))/EARef;
    //***************
    // Calculate Eclipse begin time and end time
     
    double AdjMoonSpeedFactor = 0;double AdjMoonLeadLag;int AdjMoonStatus = 0;
     
    double AdjMoon = ReqMoon + SunMoonT;
     
    year.get_calculatedMoon_data(AdjMoon, CalculatedMoon, solar_days, moon_x12_2,
    moon_minus_solar_days, AdjMoonStatus, AdjMoonLeadLag, moon_add_minus);
     
    year.getMoon_adjusted_data(AdjMoonStatus, AdjMoonLeadLag, moon_add_minus,
    moon_minus_solar_days, moon_sum, MoonX, AdjMoonSpeedFactor);
     
    double relativeMoonSun_speed_factor = AdjMoonSpeedFactor - 0.082;
    //---------------------------------
    double relative_time_difference
    = (sqrt((20 - eclipseFrac * 10)* eclipseFrac * 10)* 57.40)/relativeMoonSun_speed_factor;
     
    // Hui Hui method here
    // relative_time_difference
    //= 24*60*sqrt(addD*addD - ESeeLat*ESeeLat)/AdjMoonSpeedFactor;
     
    double eclipseBegin_hour = eclipseMax_hour - relative_time_difference/10000;
    double eclipseBegin_24hour = eclipseBegin_hour * 24;
     
    double eclipseEnd_hour = eclipseMax_hour + relative_time_difference/10000;
    double eclipseEnd_24hour = eclipseEnd_hour * 24;
    //----------------------------------
    double MidEclipse_time
    = (sqrt((20 - EclipseFracS * 10)* EclipseFracS * 10)* 57.40)/relativeMoonSun_speed_factor;
     
    double eclipseBegin_hourS = eclipseMax_hour - MidEclipse_time/10000;
    double eclipseBegin_24hourS = eclipseBegin_hourS * 24;
     
    double eclipseEnd_hourS = eclipseMax_hour + MidEclipse_time/10000;
    double eclipseEnd_24hourS = eclipseEnd_hourS * 24;
    //------------------------------------
    int Print = 0;
    if (eclipseFrac < 0){eclipseFrac = EclipseFracS;}
    if (eclipseFrac > 0)
    {
    if (calculated_year >= eclipse_theo[index].year)
    {
    if (calculated_year > eclipse_theo[index].year)
    {eclipse_shoushi[index].eclipseFrac = -99;index = index + 1;}
     
    if (calculated_year == eclipse_theo[index].year)
    {
    if (total_months >= eclipse_theo[index].month)
    {
    Pindex = index;
    eclipse_shoushi[index].year = calculated_year;
    eclipse_shoushi[index].month = total_months;
     
    eclipse_shoushi[index].eclipseFrac = eclipseFrac;
    eclipse_shoushi[index].MoonSpeedFactor = MoonSpeedFactor;
    eclipse_shoushi[index].EclipseFracS = EclipseFracS;
    eclipse_shoushi[index].SunXS = SunXS;
    eclipse_shoushi[index].MoonXS = MoonXS;
     
    eclipse_shoushi[index].eclipseBegin = eclipseBegin_24hour;
    eclipse_shoushi[index].eclipseBeginS = eclipseBegin_24hourS;
    eclipse_shoushi[index].eclipseMax = eclipseMax_24hour;
    eclipse_shoushi[index].eclipseEnd = eclipseEnd_24hour;
    eclipse_shoushi[index].eclipseEndS = eclipseEnd_24hourS;
     
    eclipse_shoushi[index].Winter_solitice = Winter_solitice;
    eclipse_shoushi[index].SunStatus = SunStatus;
     
    eclipse_shoushi[index].observed_time = observed_time;
    eclipse_shoushi[index].half_day_hours = half_day_hours;
    eclipse_shoushi[index].ReqSun = ReqSun;
    eclipse_shoushi[index].ReqMoon = ReqMoon;
     
    eclipse_shoushi[index].Gon = Gon;
    eclipse_shoushi[index].MoonStatus = MoonStatus;
    eclipse_shoushi[index].SunMoonXS = SunMoonXS;
    eclipse_shoushi[index].MDayT_hour = MDayT_hour;
    eclipse_shoushi[index].observed_timeS = observed_timeS;
     
    eclipse_shoushi[index].SunStatusS = SunStatusS;
    eclipse_shoushi[index].MoonStatusS = MoonStatusS;
     
    index = index + 1;
    Print = 1;
    }
    }
    }
     
    if (Print == 1){//if_Print
    if (calculated_year == eclipse_theo[Pindex].year)
    {
    if (total_months >= eclipse_theo[Pindex].month)
    {
    if (calculated_year <= end_year)
    {
    printf("---------\n");
    printf("%4d(%1d) \t",calculated_year,total_months);
    printf("SprRed\tWinYel\tT_12Red\tYelRedXA+YelWhXA=RedWhXA[+]URedXA=UWhXA\tHighLow*sin(UseSunHL)=UseHL\n");
    printf("\t\t%3.2f\t%3.2f\t%3.2f\t%3.2f + %3.2f = %3.2f[+]%3.2f=%3.2f\t%3.2f*sin(%3.2f)=%3.2f\n",
    SprRed,WinYel,T_12Red,YelRedXA,YelWhXA,RedWhXA,URedXA,UWhXA,HighLow,UseSunHL,UseHL);
    printf("MidMoon-MaxMoon\tLeapDay ReqMDay\tReqSun()360\tReqMoon()360\tReqEcl");
    printf(" SunPos+UsePos=RSunYel UseRMoon RSunEnd/RMoonEnd\n");
    printf("%3.2f-%3.2f\t%3.2f\t%3.2f\t%3.2f(%1d)%3.2f",
    MidMoon,MaxMoon,LeapDay,ReqMoonDay, ReqSun,SunStatus,ReqSun360);
    printf("\t%3.2f(%1d)%3.2f\t%3.2f", ReqMoon,MoonStatus,ReqMoon360, ReqEcl);
    printf("\t%3.2f+%3.2f=%3.2f\t%3.2f\t%3.2f/%3.2f\n", SunPos,UsePos,RSunYel,UseRMoon, RSunEnd,RMoonEnd);
    printf("=%3.2f\n",ReqMoon_HB);
    printf("MinSun(360)\tT_12X\tT_12NP\tT_12Arc\tOReqMoon(360)\tSun_Th/SunXS/MoonXS\tEe_o\tEMoon360/MOff/MF\n");
    printf("%3.2f(%3.2f)\t%3.2f\t%3.2f\t%3.2f\t%3.2f(%3.2f)",MinSun,MinSun360,T_12X,T_12NP,T_12Arc,OReqMoon,OReqMoon360);
    printf("\t%3.2f/%3.2f/%3.2f \t%4.3f\t%3.2f/%3.2f/%3.3f\n",SunX_Th,SunXS,MoonXS,Ee_o,EMoon360, MOff,MF);
     
    printf("Eclipse)Frac\tMaxPh(MDayT)\tTestEA\tSexEA\n");
    printf("\t%4.3f\t%4.3f(%4.3f)",eclipseFrac,eclipseMax_24hour, MDayT_24hour);
    if (MoonLeadLag >= 82){cout << "(>)";}
    printf("\t%3.2f\t%3.2f\n",TestEA,SexEA);
    printf("MDay0_24 MoonT24 SeeT(min)\tXA+YelWh=XYelA\tLatYel\tMoonPos/YelD\tSolEA(360)\n");
    printf("%3.3f\t%3.3f\t%3.2f(%3.2f)\t%3.2f+%2.2f=%3.2f\t%3.2f\t%3.2f/%3.3f\t%3.2f(%3.2f)\n",
    MDay0_24,MoonT24,SeeHour,SeeT,XA,YelWh,XYelA,LatYel,MoonPos,YelD,SolEA,SolEA360);
    }
     
    printf("\tNS_Sh/EW_Sh\tNS_Hu/EW_Hu\tESunLat\tEW_f\tMaxLat+NS_f=ESeeLat\taddD\n");
    printf("\t%3.2f/%3.2f\t%3.2f/%3.2f\t%3.2f\t%3.2f\t%3.2f + %3.2f=%3.3f\t%3.3f\n\n",
    NS_Sh,EW_Sh, NS_Hu,EW_Hu,ESunLat, EW_f,MaxLat,NS_f,ESeeLat,addD);
    }
    }
    }//if_Print end
    }// if_1 end
    }
    }
    cout << endl;
    cout << "========================================\n"
    << "\tShoushi calendar method referring to Chonzhen Lishu 200 year eclipse table for drift problem\n"
    << "\t\tdrift rate = offset/347 and 347 = 1628 - 1281(base year for Shoushi)\n"
    << "\t\tobserved_time is integrated with theorectical HuiHui calendar method table\n"
    << "\t\teclipse fraction is integrated from HouBian for comparison with Shoushi\n"
    << endl
     
    << "\tData record in Chongzhen Lishu(page 788-34):\n"
    << "\t\t65 fan late => 0.0065\n"
    << "\t\t6 days late => more than 6\n"
    << "\t\t2 degree late => (2/365.2425)* 27.5546 => more than 0.150883864\n\n"
     
    << "\tShoushi drift from Adam Schall 200 eclipse table for year 1628\n"
    << "\t\tcalculated data from http://ideone.com/BlkMI\n"
    << "\t\tfor moon_days = -(0.005499/347) => approx -0.0000158473\n"
    << "\t\tfor sun = -(5.99089/347) => -0.0172648\n"
    << "\t\tfor moon = (0.162147/347) => 0.000467281\n"
    << "\t\tfor eclipse = (0.0122486/347) => 0.0000352987\n"
    << "-----------------------------\n"
    << endl
    << " slightly modified according to data record in Chongzhen Lishu(page 788-34) as :\n"
    << "RefMoonDay = -0.0081 + RefMoonDay + 0.01041 - 0.01041/347 *(calculated_year - 1281)\n"
    << endl
    << "MinSun = 6.572 - 6.1419 +(6.1419/347)*(calculated_year - 1281)\n"
    << "RefSun = BaseSun - MinSun\n"
    << endl
    << "RefMoon = 0.0402 + RefMoon + 0.16309/347 *(calculated_year - 1281)\n"
    << "RefEcl = RefEcl + 0.01735/347*(calculated_year - 1281)\n"
    << endl
    << "----------------------------\n"
    << "\t\tInverted Kepler elliptic model is implemented as :\n"
    << "\tReqSun360 = fmod(ReqSun * 360/Tropical_year + 180, 360)\n"
    << "\te = 0.0167\n"
    << "\tSunX_Th = 180/pi*(2 *e* sin(ReqSun360 *pi/180)+ 1.25* e*e * sin(2*ReqSun360 *pi/180))\n"
    << "\tSFactor = 13.36875 * 0.082 + 0.082// = 1.1782375\n"
    << "\tSunX360 = SFactor * SunX_Th\n"
    << "\tSunX = Tropical_year/360* SunX360\n"
    << "\t=================================================\n"
    << "\tRSunYel = ReqSun360 + SunX_Th + MinSun360 + 0.098\n"
    << "\t=================================================\n"
    << endl
    << "----------------------\n"
    << "\t\tConvert Shoushi ReqMoon to consider Ee_o\n"
    << "\tOReqMoon = ReqMoon\n"
    << "\tOReqMoon360 = fmod(OReqMoon * 360/Zhuanzhong_C + 180, 360)\n"
    << "\tif (OReqMoon360 < 0){OReqMoon360 = fmod(360 + OReqMoon360, 360);}\n\n"
     
    << "\t\tapply Earth Center Offset for Moon adjustment\n"
    << "\tEe_o = 0\n"
    << "\tif (cos(OReqMoon360 * pi/180)< 0){Ee_o = 0}\n"
    << "\ty_scale = 1\n"
    << "\tx_cir = cos(OReqMoon360 * pi/180)\n"
    << "\ty_cir = sin(OReqMoon360 * pi/180)\n"
    << "\tr_elip = 1/sqrt(x_cir*x_cir + y_cir*y_cir/(y_scale*y_scale))\n"
    << "\tx_o = r_elip * x_cir\n"
    << "\ty_o = r_elip * y_cir\n"
     
    << "\t\tEMoon360 = 180/pi *atan(y_o/(x_o - Ee_o))\n"
    << "\t\tif ((OReqMoon360 > 90)&(OReqMoon360 < 270)){EMoon360 = 180 + EMoon360;}\n"
    << "\t\tif (OReq_moon360 > 270){EMoon360 = 360 + EMoon360;}\n"
    << "----------------------\n"
    << "\tMOff = EMoon360 - OReqMoon360\n"
    << "\tEMoon = fmod((EMoon360 + 180)*Zhuanzhong_C/360, Zhuanzhong_C)\n\n"
    << "\tReqMoon = EMoon\n\n"
     
    << "\tMDay0 = ReqMoonDay -(int)ReqMoonDay\n\n"
     
    << "\t\tHuiHui observed time table is implemented as :\n"
    << "\tdegree = fmod(RSunEnd + 270, 360)\n"
    << "\tGon = degree/30.0\n"
    << "\tHour = MDayT_hour * 24\n"
    << endl
    << "int SeeTimeTable[][19]\n"
    << " = {\n"
    << "h_deg -120 -105 -90 -75 -60 -45 -30 -15 0 15 30 45 60 75 90 105 120 <--\n"
    << "Gon h20 h19 h18 h17 h16 h15 h14 h13 h12 h11 h10 h9 h8 h7 h6 h5 h4 Gon\n"
    << "-----------------------------------------------------------------------------------------------\n"
    << "{9, -61, -72, -81, -85, -85, -78, -61, -34, -0, 34, 61, 78, 85, 85, 81, 72, 61, 9},\n"
    << "{10, -51, -61, -68, -70, -66, -54, -32, -2, 33, 63, 84, 95, 99, 96, 90, 80, 67, 8},\n"
    << "{11, -40, -49, -55, -55, -50, -37, -14, 15, 48, 75, 93, 103, 106, 103, 96, 86, 74, 7},\n"
    << "{0, -33, -42, -48, -50, -46, -35, -15, 13, 44, 72, 92, 103, 108, 106, 100, 91, 80, 6},\n"
    << "{1, -34, -44, -52, -55, -54, -45, -28, -3, 29, 60, 84, 98, 105, 105, 100, 92, 81, 5},\n"
    << "{2, -45, -56, -64, -69, -69, -62, -46, -20, 13, 46, 73, 90, 97, 98, 94, 86, 75, 4},\n"
    << "{3, -61, -72, -81, -85, -85, -78, -61, -34, 0, 34, 61, 78, 85, 85, 81, 72, 61, 3},\n"
    <<"-----------------------------------------------------------------------------------------------\n"
    << "Gon h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16 h17 h18 h19 h20 Gon\n"
    << "-->\n" << endl
     
    << "\tobserved_time = SeeT/(24*60)\n"
    << endl
    << "\teclipseMax_hour = MDayT_hour + observed_time\n\n"
    << "\tt\tt\tt\tt\tt\tt\tt\tt\tt\tt\tt\tt\tt\tt" << endl
    << "======================================================================================================================\n"
    << "Shoushi(year)(months)"
    << "\t\teclipseFraction \teclipseBegin\t\tMaxPhase\t\teclipseEnd" << endl;
    cout << "Theorectical(month/day)"
    << "\t\t\t(EclipseFracS)\t(eclipseBeginS)\t(observed_time x 24 x 60)\t(eclipse_endS)" << endl;
    cout << " half_day_hours/ReqSun/ReqMoon"
    << "\t\t\t\t\t\t(observed_timeS x 24 x 60)" << endl
    << " SunXS(SunStatusS)/MoonSpeedFactor" << endl
    << " MoonXS(MoonStatusS)/SunMoonXS" << endl
    << " Gon/MDayT" << endl;
    cout << "___________________________________________________________"
    << "___________________________________________________________" << endl
    << endl;
     
    int record;
    for (record = 0; record < index; record ++)
    {
    printf("%4d(%1d)",eclipse_shoushi[record].year,eclipse_shoushi[record].month);
    printf(" %4.2f",eclipse_shoushi[record].half_day_hours);
    printf("/%2.1f/%2.1f", eclipse_shoushi[record].ReqSun, eclipse_shoushi[record].ReqMoon);
    printf("\t %5.4f(%5.4f)",eclipse_shoushi[record].eclipseFrac, eclipse_shoushi[record].EclipseFracS);
    cout << "\t " << eclipse_shoushi[record].eclipseBegin << "(" << eclipse_shoushi[record].eclipseBeginS << ")"
    << "\t " << eclipse_shoushi[record].eclipseMax << "(" << eclipse_shoushi[record].observed_time *24*60 << ")\t "
    << eclipse_shoushi[record].eclipseEnd << "(" << eclipse_shoushi[record].eclipseEndS << ")\n";
     
    cout << eclipse_theo[record].month << "/" << eclipse_theo[record].day;
    printf("\t%4.3f(%1d)/", eclipse_shoushi[record].SunXS, eclipse_shoushi[record].SunStatusS );
    printf("%4.3f ", eclipse_shoushi[record].MoonSpeedFactor);
     
    printf("\t-)%5.4f", eclipse_theo[record].eclipseFrac);
    printf("\t\t-)%5.4f", eclipse_theo[record].eclipseBegin);
    printf("\t\t-)%5.4f", eclipse_theo[record].eclipseMax);
    printf("(%5.4f)", eclipse_shoushi[record].observed_timeS * 24 * 60);
    printf(" \t-)%6.4f\n", eclipse_theo[record].eclipseEnd);
     
    printf("\t%4.3f(%1d)/", eclipse_shoushi[record].MoonXS, eclipse_shoushi[record].MoonStatusS);
    printf("%4.3f ", eclipse_shoushi[record].SunMoonXS);
    cout << "\t--------------------------------------------------------------------------------------" << endl;
    printf("\t%3.2f/%3.2f", eclipse_shoushi[record].Gon, eclipse_shoushi[record].MDayT_hour * 24);
     
    printf("\t\t %5.4f",eclipse_shoushi[record].eclipseFrac - eclipse_theo[record].eclipseFrac);
    printf("(%5.4f)", eclipse_shoushi[record].EclipseFracS - eclipse_theo[record].eclipseFrac);
    printf("\t %5.4f", eclipse_shoushi[record].eclipseBegin - eclipse_theo[record].eclipseBegin);
    printf("(%5.4f)\t ",eclipse_shoushi[record].eclipseBeginS - eclipse_theo[record].eclipseBegin);
    printf("%5.4f\t\t", eclipse_shoushi[record].eclipseMax - eclipse_theo[record].eclipseMax);
    printf("%5.4f", eclipse_shoushi[record].eclipseEnd - eclipse_theo[record].eclipseEnd);
    printf("(%5.4f)\n\n", eclipse_shoushi[record].eclipseEndS - eclipse_theo[record].eclipseEnd);
    };
     
    cout << endl;
    return 0;
    }
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    // 求正交平均
    // PosExtra = -570 *SunX_Th/6973;
    //======================================
    // Convert Shoushi ReqEcl as SunPos
    //======================================
    // 求日距正交
    // SunPos = (ReqEcl - Jiaozhong_C/2)* 360/Jiaozhong_C - PosExtra + SunX_Th;
    // if (SunPos < 0){SunPos = SunPos + 360;}
     
    // 求用正交
    // UsePos = RSunYel - SunPos;
     
    // MoonPos = RSunYel - PosReal
     
    // LatYel = 180/pi *asin(sin(YelWh *pi/180)*sin(MoonPos *pi/180))
    // Here use LatYel = 180/pi *asin(sin(YelWh *pi/180)*sin(SolEA360 *pi/180));
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    // ESeeLat = UseSeeLat/60;
    // SunHalfD = SunRSemiDiameter/60;
    // addD = NetDiameter/60;
     
    // eclipseFrac =(addD - abs(ESeeLat))/(2* SunHalfD);
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    // PosRX360 = PosRealAdd * SFactor;
    // PosRX = PosRX360 * Tropical_year/360;
     
    // SolEA = SolarEclipse_angle - PosX - PosRX
    // SolEA360 = fmod(SolEA + 363.7934196/2, 363.7934196)*360/363.7934196
    //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    // EclipseFracS modified to include PosX and PosRX as :
    //--------------------------------------------
    // if (Front == -1){RefEA_Sh = 188.05}
    // if (Front == 1){RefEA_Sh = 357.64}
    // TestEA = SolarEclipse_angle - RefEA_Sh - PosX - PosRX
    // NetEA = TestEA -(NS_Sh + EW_Sh)
    // SexEA = NetEA
    // if (Front == -1){SexEA = -1*NetEA}
     
    // EARef = 6;
    // if (SexEA < 0){EARef = 8}
    // EclipseFracS = (EARef - abs(SexEA))/EARef


run \| edit \| history \| help	0