#include "SSQ.h"

#include <math.h>

#include "const.h"

#include "eph.h"

SSQ::SSQ()

{

std::string suoS = "", qiS = "";

// 619-01-21开始16598个朔日修正表 d0=1947168

suoS = "EqoFscDcrFpmEsF2DfFideFelFpFfFfFiaipqti1ksttikptikqckstekqttgkqttgkqteksttikptikq2fjstgjqttjkqttgkqt";

suoS += "ekstfkptikq2tijstgjiFkirFsAeACoFsiDaDiADc1AFbBfgdfikijFifegF1FhaikgFag1E2btaieeibggiffdeigFfqDfaiBkF";

suoS += "1kEaikhkigeidhhdiegcFfakF1ggkidbiaedksaFffckekidhhdhdikcikiakicjF1deedFhFccgicdekgiFbiaikcfi1kbFibef";

suoS += "gEgFdcFkFeFkdcfkF1kfkcickEiFkDacFiEfbiaejcFfffkhkdgkaiei1ehigikhdFikfckF1dhhdikcfgjikhfjicjicgiehdik";

suoS += "cikggcifgiejF1jkieFhegikggcikFegiegkfjebhigikggcikdgkaFkijcfkcikfkcifikiggkaeeigefkcdfcfkhkdgkegieid";

suoS += "hijcFfakhfgeidieidiegikhfkfckfcjbdehdikggikgkfkicjicjF1dbidikFiggcifgiejkiegkigcdiegfggcikdbgfgefjF1";

suoS += "kfegikggcikdgFkeeijcfkcikfkekcikdgkabhkFikaffcfkhkdgkegbiaekfkiakicjhfgqdq2fkiakgkfkhfkfcjiekgFebicg";

suoS += "gbedF1jikejbbbiakgbgkacgiejkijjgigfiakggfggcibFifjefjF1kfekdgjcibFeFkijcfkfhkfkeaieigekgbhkfikidfcje";

suoS += "aibgekgdkiffiffkiakF1jhbakgdki1dj1ikfkicjicjieeFkgdkicggkighdF1jfgkgfgbdkicggfggkidFkiekgijkeigfiski";

suoS += "ggfaidheigF1jekijcikickiggkidhhdbgcfkFikikhkigeidieFikggikhkffaffijhidhhakgdkhkijF1kiakF1kfheakgdkif";

suoS += "iggkigicjiejkieedikgdfcggkigieeiejfgkgkigbgikicggkiaideeijkefjeijikhkiggkiaidheigcikaikffikijgkiahi1";

suoS += "hhdikgjfifaakekighie1hiaikggikhkffakicjhiahaikggikhkijF1kfejfeFhidikggiffiggkigicjiekgieeigikggiffig";

suoS += "gkidheigkgfjkeigiegikifiggkidhedeijcfkFikikhkiggkidhh1ehigcikaffkhkiggkidhh1hhigikekfiFkFikcidhh1hit";

suoS += "cikggikhkfkicjicghiediaikggikhkijbjfejfeFhaikggifikiggkigiejkikgkgieeigikggiffiggkigieeigekijcijikgg";

suoS += "ifikiggkideedeijkefkfckikhkiggkidhh1ehijcikaffkhkiggkidhh1hhigikhkikFikfckcidhh1hiaikgjikhfjicjicgie";

suoS += "hdikcikggifikigiejfejkieFhegikggifikiggfghigkfjeijkhigikggifikiggkigieeijcijcikfksikifikiggkidehdeij";

suoS += "cfdckikhkiggkhghh1ehijikifffffkhsFngErD1pAfBoDd1BlEtFqA2AqoEpDqElAEsEeB2BmADlDkqBtC1FnEpDqnEmFsFsAFn";

suoS += "llBbFmDsDiCtDmAB2BmtCgpEplCpAEiBiEoFqFtEqsDcCnFtADnFlEgdkEgmEtEsCtDmADqFtAFrAtEcCqAE1BoFqC1F1DrFtBmF";

suoS += "tAC2ACnFaoCgADcADcCcFfoFtDlAFgmFqBq2bpEoAEmkqnEeCtAE1bAEqgDfFfCrgEcBrACfAAABqAAB1AAClEnFeCtCgAADqDoB";

suoS += "mtAAACbFiAAADsEtBqAB2FsDqpFqEmFsCeDtFlCeDtoEpClEqAAFrAFoCgFmFsFqEnAEcCqFeCtFtEnAEeFtAAEkFnErAABbFkAD";

suoS += "nAAeCtFeAfBoAEpFtAABtFqAApDcCGJ";

//1645-09-23开始7567个节气修正表

qiS = "FrcFs22AFsckF2tsDtFqEtF1posFdFgiFseFtmelpsEfhkF2anmelpFlF1ikrotcnEqEq2FfqmcDsrFor22FgFrcgDscFs22FgEe";

qiS += "FtE2sfFs22sCoEsaF2tsD1FpeE2eFsssEciFsFnmelpFcFhkF2tcnEqEpFgkrotcnEqrEtFermcDsrE222FgBmcmr22DaEfnaF22";

qiS += "2sD1FpeForeF2tssEfiFpEoeFssD1iFstEqFppDgFstcnEqEpFg11FscnEqrAoAF2ClAEsDmDtCtBaDlAFbAEpAAAAAD2FgBiBqo";

qiS += "BbnBaBoAAAAAAAEgDqAdBqAFrBaBoACdAAf1AACgAAAeBbCamDgEifAE2AABa1C1BgFdiAAACoCeE1ADiEifDaAEqAAFe1AcFbcA";

qiS += "AAAAF1iFaAAACpACmFmAAAAAAAACrDaAAADG0";

SB = jieya(suoS); //定朔修正表解压

QB = jieya(qiS); //定气修正表解压

//朔直线拟合参数

long double suoKBTmp[] = {

1457698.231017,29.53067166, // -721-12-17 h=0.00032 古历·春秋

1546082.512234,29.53085106, // -479-12-11 h=0.00053 古历·战国

1640640.735300,29.53060000, // -221-10-31 h=0.01010 古历·秦汉

1642472.151543,29.53085439, // -216-11-04 h=0.00040 古历·秦汉

1683430.509300,29.53086148, // -104-12-25 h=0.00313 汉书·律历志(太初历)平气平朔

1752148.041079,29.53085097, // 85-02-13 h=0.00049 后汉书·律历志(四分历)

1807665.420323,29.53059851, // 237-02-12 h=0.00033 晋书·律历志(景初历)

1883618.114100,29.53060000, // 445-01-24 h=0.00030 宋书·律历志(何承天元嘉历)

1907360.704700,29.53060000, // 510-01-26 h=0.00030 宋书·律历志(祖冲之大明历)

1936596.224900,29.53060000, // 590-02-10 h=0.01010 随书·律历志(开皇历)

1939135.675300,29.53060000, // 597-01-24 h=0.00890 随书·律历志(大业历)

1947168.00// 619-01-21

};

suoKB = new std::vector<long double>(suoKBTmp, suoKBTmp + sizeof(suoKBTmp) / sizeof(long double));

long double qiKBTmp[] = {

1640650.479938,15.21842500, // -221-11-09 h=0.01709 古历·秦汉

1642476.703182,15.21874996, // -216-11-09 h=0.01557 古历·秦汉

1683430.515601,15.218750011, // -104-12-25 h=0.01560 汉书·律历志(太初历)平气平朔 回归年=365.25000

1752157.640664,15.218749978, // 85-02-23 h=0.01559 后汉书·律历志(四分历) 回归年=365.25000

1807675.003759,15.218620279, // 237-02-22 h=0.00010 晋书·律历志(景初历) 回归年=365.24689

1883627.765182,15.218612292, // 445-02-03 h=0.00026 宋书·律历志(何承天元嘉历) 回归年=365.24670

1907369.128100,15.218449176, // 510-02-03 h=0.00027 宋书·律历志(祖冲之大明历) 回归年=365.24278

1936603.140413,15.218425000, // 590-02-17 h=0.00149 随书·律历志(开皇历) 回归年=365.24220

1939145.524180,15.218466998, // 597-02-03 h=0.00121 随书·律历志(大业历) 回归年=365.24321

1947180.798300,15.218524844, // 619-02-03 h=0.00052 新唐书·历志(戊寅元历)平气定朔 回归年=365.24460

1964362.041824,15.218533526, // 666-02-17 h=0.00059 新唐书·历志(麟德历) 回归年=365.24480

1987372.340971,15.218513908, // 729-02-16 h=0.00096 新唐书·历志(大衍历,至德历) 回归年=365.24433

1999653.819126,15.218530782, // 762-10-03 h=0.00093 新唐书·历志(五纪历) 回归年=365.24474

2007445.469786,15.218535181, // 784-02-01 h=0.00059 新唐书·历志(正元历,观象历) 回归年=365.24484

2021324.917146,15.218526248, // 822-02-01 h=0.00022 新唐书·历志(宣明历) 回归年=365.24463

2047257.232342,15.218519654, // 893-01-31 h=0.00015 新唐书·历志(崇玄历) 回归年=365.24447

2070282.898213,15.218425000, // 956-02-16 h=0.00149 旧五代·历志(钦天历) 回归年=365.24220

2073204.872850,15.218515221, // 964-02-16 h=0.00166 宋史·律历志(应天历) 回归年=365.24437

2080144.500926,15.218530782, // 983-02-16 h=0.00093 宋史·律历志(乾元历) 回归年=365.24474

2086703.688963,15.218523776, // 1001-01-31 h=0.00067 宋史·律历志(仪天历,崇天历) 回归年=365.24457

2110033.182763,15.218425000, // 1064-12-15 h=0.00669 宋史·律历志(明天历) 回归年=365.24220

2111190.300888,15.218425000, // 1068-02-15 h=0.00149 宋史·律历志(崇天历) 回归年=365.24220

2113731.271005,15.218515671, // 1075-01-30 h=0.00038 李锐补修(奉元历) 回归年=365.24438

2120670.840263,15.218425000, // 1094-01-30 h=0.00149 宋史·律历志 回归年=365.24220

2123973.309063,15.218425000, // 1103-02-14 h=0.00669 李锐补修(占天历) 回归年=365.24220

2125068.997336,15.218477932, // 1106-02-14 h=0.00056 宋史·律历志(纪元历) 回归年=365.24347

2136026.312633,15.218472436, // 1136-02-14 h=0.00088 宋史·律历志(统元历,乾道历,淳熙历) 回归年=365.24334

2156099.495538,15.218425000, // 1191-01-29 h=0.00149 宋史·律历志(会元历) 回归年=365.24220

2159021.324663,15.218425000, // 1199-01-29 h=0.00149 宋史·律历志(统天历) 回归年=365.24220

2162308.575254,15.218461742, // 1208-01-30 h=0.00146 宋史·律历志(开禧历) 回归年=365.24308

2178485.706538,15.218425000, // 1252-05-15 h=0.04606 淳祐历 回归年=365.24220

2178759.662849,15.218445786, // 1253-02-13 h=0.00231 会天历 回归年=365.24270

2185334.020800,15.218425000, // 1271-02-13 h=0.00520 宋史·律历志(成天历) 回归年=365.24220

2187525.481425,15.218425000, // 1277-02-12 h=0.00520 本天历 回归年=365.24220

2188621.191481,15.218437494, // 1280-02-13 h=0.00015 元史·历志(郭守敬授时历) 回归年=365.24250

2322147.76// 1645-09-21

};

qiKB = new std::vector<long double>(qiKBTmp, qiKBTmp + sizeof(qiKBTmp) / sizeof(long double));

}

SSQ::~SSQ()

{

delete suoKB;

delete qiKB;

}

void str_replace(std::string & str, const std::string strsrc, const std::string strdst)

{

std::string::size_type pos = 0;//位置

std::string::size_type srclen = strsrc.size();//要替换的字符串大小

std::string::size_type dstlen = strdst.size();//目标字符串大小

while ((pos = str.find(strsrc, pos)) != std::string::npos)

{

str.replace(pos, srclen, strdst);

pos += dstlen;

}

}

std::string SSQ::jieya(std::string s) { //气朔解压缩

std::string o = "0000000000", o2 = o + o;

str_replace(s, "J", "00");

str_replace(s, "I", "000");

str_replace(s, "H", "0000");

str_replace(s, "G", "00000");

str_replace(s, "t", "02");

str_replace(s, "s", "002");

str_replace(s, "r", "0002");

str_replace(s, "q", "00002");

str_replace(s, "p", "000002");

str_replace(s, "o", "0000002");

str_replace(s, "n", "00000002");

str_replace(s, "m", "000000002");

str_replace(s, "l", "0000000002");

str_replace(s, "k", "01");

str_replace(s, "j", "0101");

str_replace(s, "i", "001");

str_replace(s, "h", "001001");

str_replace(s, "g", "0001");

str_replace(s, "f", "00001");

str_replace(s, "e", "000001");

str_replace(s, "d", "0000001");

str_replace(s, "c", "00000001");

str_replace(s, "b", "000000001");

str_replace(s, "a", "0000000001");

str_replace(s, "A", o2 + o2 + o2);

str_replace(s, "B", o2 + o2 + o);

str_replace(s, "C", o2 + o2);

str_replace(s, "D", o2 + o);

str_replace(s, "E", o2);

str_replace(s, "F", o);

return s;

}

int SSQ::calc(long double jd, QSType qs)

{

jd += 2451545;

int i;

long double D;

std::string n;

std::vector<long double> B = *suoKB;

long double pc = 14;

//如果查的是气朔

if (qs == QType)

{

B = *qiKB, pc = 7;

}

long double f1 = B[0] - pc, f2 = B[B.size() - 1] - pc, f3 = 2436935;

if (jd < f1 || jd >= f3)

{

//平气朔表中首个之前，使用现代天文算法。1960.1.1以后，使用现代天文算法 (这一部分调用了qi_high和so_high,所以需星历表支持)

if (qs == QType)

{

return floor(qi_high(floor((jd + pc - 2451259) / 365.2422 * 24) * PI / 12) + 0.5); //2451259是1999.3.21,太阳视黄经为0,春分.定气计算

}

else

{

return floor(so_high(floor((jd + pc - 2451551) / 29.5306) * PI * 2) + 0.5); //2451551是2000.1.7的那个朔日,黄经差为0.定朔计算

}

}

if (jd >= f1 && jd < f2) { //平气或平朔

for (i = 0; i < B.size(); i += 2) if (jd + pc < B[i + 2]) break;

D = B[i] + B[i + 1] * floor((jd + pc - B[i]) / B[i + 1]);

D = floor(D + 0.5);

if (D == 1683460) D++; //如果使用太初历计算-103年1月24日的朔日,结果得到的是23日,这里修正为24日(实历)。修正后仍不影响-103的无中置闰。如果使用秦汉历，得到的是24日，本行D不会被执行。

return D - 2451545;

}

if (jd >= f2 && jd < f3) { //定气或定朔

if (qs == QType) {

D = floor(qi_low(floor((jd + pc - 2451259) / 365.2422 * 24) * PI / 12) + 0.5); //2451259是1999.3.21,太阳视黄经为0,春分.定气计算

n = QB.substr(floor((jd - f2) / 365.2422 * 24), 1); //找定气修正值

}

else {

D = floor(so_low(floor((jd + pc - 2451551) / 29.5306) * PI * 2) + 0.5); //2451551是2000.1.7的那个朔日,黄经差为0.定朔计算

n = SB.substr(floor((jd - f2) / 29.5306), 1); //找定朔修正值

}

if (n == "1") return D + 1;

if (n == "2") return D - 1;

return D;

}

return 0;

}

long double SSQ::qi_high(long double W)

{

long double t = XL::S_aLon_t2(W) * 36525;

t = t - dt_T(t) + 8.0 / 24.0;

long double v = fmod(t + 0.5, 1) * 86400;

if (v < 1200 || v >86400 - 1200) t = XL::S_aLon_t(W) * 36525 - dt_T(t) + 8.0 / 24.0;

return t;

}

long double SSQ::so_high(long double W)

{ //较高精度朔

long double t = XL::MS_aLon_t2(W) * 36525;

t = t - dt_T(t) + 8.0 / 24.0;

long double v = fmod(t + 0.5, 1) * 86400;

if (v < 1800 || v >86400 - 1800) t = XL::MS_aLon_t(W) * 36525 - dt_T(t) + 8.0 / 24.0;

return t;

}

long double SSQ::so_low(long double W) { //低精度定朔计算,在2000年至600，误差在2小时以内(仍比古代日历精准很多)

long double v = 7771.37714500204;

long double t = (W + 1.08472) / v, L;

t -= (-0.0000331*t*t

+ 0.10976 *cos(0.785 + 8328.6914*t)

+ 0.02224 *cos(0.187 + 7214.0629*t)

- 0.03342 *cos(4.669 + 628.3076*t)) / v

+ (32 * (t + 1.8)*(t + 1.8) - 20) / 86400 / 36525;

return t * 36525 + 8.0 / 24.0;

}

long double SSQ::qi_low(long double W) { //最大误差小于30分钟，平均5分

long double t, L, v = 628.3319653318;

t = (W - 4.895062166) / v; //第一次估算,误差2天以内

t -= (53 * t*t + 334116 * cos(4.67 + 628.307585*t) + 2061 * cos(2.678 + 628.3076*t)*t) / v / 10000000; //第二次估算,误差2小时以内

L = 48950621.66 + 6283319653.318*t + 53 * t*t //平黄经

+ 334166 * cos(4.669257 + 628.307585*t) //地球椭圆轨道级数展开

+ 3489 * cos(4.6261 + 1256.61517*t) //地球椭圆轨道级数展开

+ 2060.6 * cos(2.67823 + 628.307585*t) * t //一次泊松项

- 994 - 834 * sin(2.1824 - 33.75705*t); //光行差与章动修正

t -= (L / 10000000 - W) / 628.332 + (32 * (t + 1.8)*(t + 1.8) - 20) / 86400 / 36525;

return t * 36525 + 8.0 / 24.0;

}

//农历排月序计算,可定出农历,有效范围：两个冬至之间(冬至一 <= d < 冬至二)

void SSQ::calcY(int jd) {

std::vector<long double>& A = ZQ;

std::vector<int>& B = HS; //中气表,日月合朔表(整日)

int i, k;

long double W, w;

//该年的气

W = int2((jd - 355 + 183) / 365.2422)*365.2422 + 355; //355是2000.12冬至,得到较靠近jd的冬至估计值

if (calc(W, QType) > jd)

{

W -= 365.2422;

}

//25个节气时刻(北京时间),从冬至开始到下一个冬至以后;

A.clear();

for (i = 0; i < 25; i++)

{

int t = calc(W + 15.2184*i, QType);

A.push_back(t);

}

ZQ_pe1 = calc(W - 15.2, QType); ZQ_pe2 = calc(W - 30.4, QType); //补算二气,确保一年中所有月份的“气”全部被计算在内

//今年"首朔"的日月黄经差w

w = calc(A[0], SType); //求较靠近冬至的朔日

if (w > A[0])

{

w -= 29.53;

}

//该年所有朔,包含14个月的始末

B.clear();

for (i = 0; i < 15; i++)

{

B.push_back( calc(w + 29.5306*i, SType) );

}

//月大小

leap = 0;

dx.clear();

ym.clear();

for (i = 0; i < 14; i++) {

dx.push_back( HS[i + 1] - HS[i] ); //月大小

ym.push_back(i); //月序初始化

}

//-721年至-104年的后九月及月建问题,与朔有关，与气无关

int YY = int2((ZQ[0] + 10 + 180) / 365.2422) + 2000; //确定年份

if (YY >= -721 && YY <= -104) {

int ns[12];

int yy;

for (i = 0; i < 3; i++) {

yy = YY + i - 1;

//颁行历年首, 闰月名称, 月建

if (yy >= -721)

{

ns[i] = calc(1457698 - J2000 + int2(0.342 + (yy + 721)*12.368422)*29.5306, SType);

//ns[i + 3] = '十三';

ns[i + 6] = 2; //春秋历,ly为-722.12.17

}

if (yy >= -479)

{

ns[i] = calc(1546083 - J2000 + int2(0.500 + (yy + 479)*12.368422)*29.5306, SType);

//ns[i + 3] = '十三';

ns[i + 6] = 2; //战国历,ly为-480.12.11

}

if (yy >= -220) {

ns[i] = calc(1640641 - J2000 + int2(0.866 + (yy + 220)*12.369000)*29.5306, SType);

// ns[i + 3] = '后九';

ns[i + 6] = 11; //秦汉历,ly为-221.10.31

}

}

int nn, f1;

for (i = 0; i < 14; i++) {

for (nn = 2; nn >= 0; nn--) if (HS[i] >= ns[nn]) break;

f1 = int2((HS[i] - ns[nn] + 15) / 29.5306); //该月积数

if (f1 < 12) ym[i] = (f1 + ns[nn + 6]) % 12; else ym[i] = ns[nn + 3];

}

return;

}

//无中气置闰法确定闰月,(气朔结合法,数据源需有冬至开始的的气和朔)

if (B[13] <= A[24]) { //第13月的月末没有超过冬至(不含冬至),说明今年含有13个月

for (i = 1; B[i + 1] > A[2 * i] && i < 13; i++); //在13个月中找第1个没有中气的月份

leap = i;

for (; i < 14; i++) ym[i]--;

}

//名称转换(月建别名)

for (i = 0; i < 14; i++) {

int Dm = HS[i] + J2000, v2 = ym[i]; //Dm初一的儒略日,v2为月建序号

int mc = v2 % 12; //月建对应的默认月名称：建子十一,建丑十二,建寅为正……

if (Dm >= 1724360 && Dm <= 1729794) mc = (v2 + 1) % 12; // 8.01.15至 23.12.02 建子为十二,其它顺推

else if (Dm >= 1807724 && Dm <= 1808699) mc = (v2 + 1) % 12; //237.04.12至239.12.13 建子为十二,其它顺推

else if (Dm >= 1999349 && Dm <= 1999467) mc = (v2 + 2) % 12; //761.12.02至762.03.30 建子为正月,其它顺推

else if (Dm >= 1973067 && Dm <= 1977052) { if (v2 % 12 == 0) mc = 2; if (v2 == 2) mc = 2; } //689.12.18至700.11.15 建子为正月,建寅为一月,其它不变

if (Dm == 1729794 || Dm == 1808699)

{

mc = 12;

// mc = '拾贰'; //239.12.13及23.12.02均为十二月,为避免两个连续十二月，此处改名

}

ym[i] = mc;

}

}