tyme4ts
tyme4ts copied to clipboard
6tail
有计划支持真太阳时间转换吗
export class SolarTimeUtil{
/**
* 标准时间发出地经度(角度表示,东经为正西经为负),北京时间的经度为+120度0分
*/
private J = 120.0;
/**
* 默认纬度(角度表示,北纬为正南纬为负),这里是中国标准时间发出地(陕西省渭南市蒲城县)
*/
private W = 35.0;
/** 统一东经为正 */
private lng = 120.0;
/** 北纬为正,南纬为负 */
private lat = 35.0;
private constructor(lng:number=120.0,lat:number=35.0) {
if (lng < -180.0 || lng > 180.0) {
throw new Error(`illegal longitude: ${lng}`);
}
if (lat < -90.0 || lat > 90.0) {
throw new Error(`illegal latitude: ${lat}`);
}
this.lng = lng;
this.lat = lat;
}
/**
* @param lng 经度 -180~180
* @param lat 纬度 -90~90
*/
static initLocation(lng:number=120.0,lat:number=35.0){
return new SolarTimeUtil(lng,lat);
}
/**
* 计算平太阳时
* @param solartime 公历时间
* @returns 平太阳时的公历时间 SolarTime
*/
getMeanSolarTime(solartime:SolarTime){
let spcjd = solartime.getJulianDay().getDay() ; // 將公历時间转换为儒略日 special jd,这里依然是标准时间,即this.J处的平太阳时
let deltPty = spcjd - (this.J - this.lng) * 4 / 60 / 24; // 计算地方平太阳时,每经度时差4分钟
let pty = JulianDay.fromJulianDay(deltPty).getSolarTime();
return pty;
}
/**
* 计算真太阳时
* @param solartime 公历时间
* @returns 真太阳时的公历时间 SolarTime
*/
getRealSolarTime(solartime:SolarTime):SolarTime{
let spcjd = solartime.getJulianDay().getDay(); // 將公历時间转换为儒略日 special jd,这里依然是标准时间,即this.J处的平太阳时
let realZty = this.zty(spcjd);
let zty = JulianDay.fromJulianDay(realZty).getSolarTime();
return zty;
}
/**
* 真太阳时模块,sn代表sin
*/
private sn(x:number) {
return Math.sin(x * 1.74532925199433E-02);
}
/**
* 真太阳时模块,cn代表cosine
*/
private cn(x:number) {
return Math.cos(x * 1.74532925199433E-02);
}
/**
* 真太阳时模块,返回小数部分(负数特殊) returns fractional part of a number
*/
private fpart(x:number) {
x = x - Math.floor(x);
if (x < 0) {
x = x + 1;
}
return x; //只取小数部份
}
/**
* 真太阳时模块,只取整数部份
*/
private ipart(x:number) {
if (x == 0) {
return 0;
}
return (x / Math.abs(x)) * Math.floor(Math.abs(x));
}
/**
* 真太阳时模块
* finds a parabola through three points and
* returns values of coordinates of extreme value (xe, ye) and
* zeros if any (z1, z2)
* assumes that the x values are -1, 0, +1
*/
private quad(ym:number, y0:number, yp:number) {
let nz = 0;
let A = 0.5 * (ym + yp) - y0;
let b = 0.5 * (yp - ym);
let c = y0;
let xe = -b / (2 * A); //x coord of symmetry line
let ye = (A * xe + b) * xe + c; //extreme value for y in interval
let dis = b * b - 4 * A * c; //discriminant
let z1:number = 0.0;
let z2:number = 0.0;
if (dis > 0) { //there are zeros
let dx = 0.5 * Math.sqrt(dis) / Math.abs(A);
z1 = xe - dx;
z2 = xe + dx;
if (Math.abs(z1) <= 1) {
nz = nz + 1;
} //This zero is in interval
if (Math.abs(z2) <= 1) {
nz = nz + 1;
} //This zero is in interval
if (z1 < -1) {
z1 = z2;
}
}
return [xe, ye, z1, z2, nz];
}
/**
* 真太阳时模块
* returns sine of the altitude of either the sun or the moon given the modified julian day of the UT
* @param float jd 公历時间的儒略日
* @param int LX 1月亮 2太阳日升日落 3太阳海上微光
*/
private sinalt(jd:number, LX:number):number {
let instant = jd - 2400001.0;
let t = (instant - 51544.5) / 36525; //减51544.5为相对2000年01月01日零点
let ra:number,dec:number;
if (LX == 1) {
let moon = this.moon(t);
ra = moon[0];
dec = moon[1];
} else {
let sun = this.sun(t);
ra = sun[0];
dec = sun[1];
}
let mjd0 = this.ipart(instant); //UT时间0点;returns the local sidereal time(计算观测地区的恒星时)开始
let ut = (instant - mjd0) * 24;
let t2 = (mjd0 - 51544.5) / 36525;
let gmst = 6.697374558 + 1.0027379093 * ut;
gmst = gmst + (8640184.812866 + (0.093104 - 0.0000062 * t2) * t2) * t2 / 3600;
let lmst = 24 * this.fpart((gmst + this.lng / 15) / 24); //结束
let tau = 15 * (lmst - ra); //hour angle of object
return this.sn(this.lat) * this.sn(dec) + this.cn(this.lat) * this.cn(dec) * this.cn(tau);
}
/**
* 真太阳时模块,关于太阳的
* Returns RA and DEC of Sun to roughly 1 arcmin for few hundred years either side of J2000.0
*/
private sun(t:number) {
let p2 = 2 * Math.PI;
let COSEPS = 0.91748;
let SINEPS = 0.39778;
let m = p2 * this.fpart(0.993133 + 99.997361 * t); //Mean anomaly
let dL = 6893 * Math.sin(m) + 72 * Math.sin(2 * m); //Eq centre
let L = p2 * this.fpart(0.7859453 + m / p2 + (6191.2 * t + dL) / 1296000);
//convert to RA and DEC - ecliptic latitude of Sun taken as zero
let sl = Math.sin(L);
let x = Math.cos(L);
let y = COSEPS * sl;
let Z = SINEPS * sl;
let rho = Math.sqrt(1 - Z * Z);
let dec = (360 / p2) * Math.atan(Z / rho);
let ra = (48 / p2) * Math.atan(y / (x + rho));
if (ra < 0) {
ra = ra + 24;
}
return [ra, dec];
}
/**
* 真太阳时模块,关于月球的,Returns RA and DEC of Moon to 5 arc min (ra) and 1 arc min (dec) for a few centuries either side of J2000.0
* Predicts rise and set times to within minutes for about 500 years in past - TDT and UT time diference may become significant for long times
*/
private moon(t:number) {
let p2 = 2 * Math.PI;
let ARC = 206264.8062;
let COSEPS = 0.91748;
let SINEPS = 0.39778;
let L0 = this.fpart(0.606433 + 1336.855225 * t); //mean long Moon in revs
let L = p2 * this.fpart(0.374897 + 1325.55241 * t); //mean anomaly of Moon
let LS = p2 * this.fpart(0.993133 + 99.997361 * t); //mean anomaly of Sun
let d = p2 * this.fpart(0.827361 + 1236.853086 * t); //diff longitude sun and moon
let F = p2 * this.fpart(0.259086 + 1342.227825 * t); //mean arg latitude
//longitude correction terms
let dL = 22640 * Math.sin(L) - 4586 * Math.sin(L - 2 * d);
dL = dL + 2370 * Math.sin(2 * d) + 769 * Math.sin(2 * L);
dL = dL - 668 * Math.sin(LS) - 412 * Math.sin(2 * F);
dL = dL - 212 * Math.sin(2 * L - 2 * d) - 206 * Math.sin(L + LS - 2 * d);
dL = dL + 192 * Math.sin(L + 2 * d) - 165 * Math.sin(LS - 2 * d);
dL = dL - 125 * Math.sin(d) - 110 * Math.sin(L + LS);
dL = dL + 148 * Math.sin(L - LS) - 55 * Math.sin(2 * F - 2 * d);
//latitude arguments
let S = F + (dL + 412 * Math.sin(2 * F) + 541 * Math.sin(LS)) / ARC;
let h = F - 2 * d;
//latitude correction terms
let N:number = -526 * Math.sin(h) + 44 * Math.sin(L + h) - 31 * Math.sin(h - L) - 23 * Math.sin(LS + h);
N = N + 11 * Math.sin(h - LS) - 25 * Math.sin(F - 2 * L) + 21 * Math.sin(F - L);
let lmoon = p2 * this.fpart(L0 + dL / 1296000); //Lat in rads
let bmoon = (18520 * Math.sin(S) + N) / ARC; //long in rads
//convert to equatorial coords using a fixed ecliptic
let CB = Math.cos(bmoon);
let x = CB * Math.cos(lmoon);
let V = CB * Math.sin(lmoon);
let C = Math.sin(bmoon);
let y = COSEPS * V - SINEPS * C;
let Z = SINEPS * V + COSEPS * C;
let rho = Math.sqrt(1 - Z * Z);
let dec = (360 / p2) * Math.atan(Z / rho); //算出月球的视赤纬(apparent declination)
let ra = (48 / p2) * Math.atan(y / (x + rho)); //算出月球的视赤经(apparent right ascension)
if (ra < 0) {
ra = ra + 24;
}
return [ra, dec];
}
/**
* 真太阳时模块,rise and set(升降计算) [升起时刻(真太阳时),落下时刻(真太阳时),真平太阳时差(仅类型2),升起时刻(标准时间,仅类型2),落下时刻(标准时间,仅类型2)]
* @param float jd 公历時间的儒略日值
* @param int LX 类型:1月亮;2太阳日升日落;3太阳海上微光
* @return array
*/
private risenset(jd:number, LX:number) {
let noon = Math.round(jd) - this.J / 360.0; //儒略日,中午12点,減去8小時時差
let sinho:number[] = [] //太阳盘面几何中心与理想地平面之间的夹角
sinho[1] = this.sn(8 / 60); //moonrise - average diameter used(月亮升降)
sinho[2] = this.sn(-50 / 60); //sunrise - classic value for refraction(太阳升降)
sinho[3] = this.sn(-12); //nautical twilight(海上微光)
let rise = 0; //是否有升起动作
let utrise:number|boolean = false; //升起的时间
let sett = 0; //是否有落下动作
let utset:number|boolean = false; //落下的时间
let hour = 1;
let ym = this.sinalt(noon + (hour - 1)/24, LX) - sinho[LX]; //See STEP 1 and 2 of Web page description.
let zero2 = 0; //两小时内是否进行了升起和落下两个动作(极地附近有这种情况,如1999年12月25日,经度0,纬度67.43,当天的太阳只有8分钟-_-)
let above = ym > 0 ? 1 : 0; //used later to classify non-risings 是否在地平线上方,用于判断极昼极夜
do {
//STEP 1 and STEP 3 of Web page description
let y0 = this.sinalt(noon + (hour + 0)/24, LX) - sinho[LX];
let yp = this.sinalt(noon + (hour + 1)/24, LX) - sinho[LX];
//STEP 4 of web page description
let quad = this.quad(ym, y0, yp);
let ye = quad[1];
let z1 = quad[2];
let z2 = quad[3];
let nz = quad[4];
switch (nz) { //cases depend on values of discriminant - inner part of STEP 4
case 0: //nothing - go to next time slot
break;
case 1: //simple rise / set event
if (ym < 0) { //must be a rising event
utrise = hour + z1;
rise = 1;
} else { //must be setting
utset = hour + z1;
sett = 1;
}
break;
case 2: //rises and sets within interval
if (ye < 0) { //minimum - so set then rise
utrise = hour + z2;
utset = hour + z1;
} else { //maximum - so rise then set
utrise = hour + z1;
utset = hour + z2;
}
rise = 1;
sett = 1;
zero2 = 1;
break;
}
ym = yp; //reuse the ordinate in the next interval
hour = hour + 2;
} while (!((hour == 25) || (rise * sett == 1))); //STEP 5 of Web page description - have we finished for this object?
if(utset !== false){ //注意这里转成了真太阳时
utset = Math.round(jd) - 0.5 + utset/24 - (this.J - this.lng) * 4 / 60 / 24;
}
if(utrise !== false){
utrise = Math.round(jd) - 0.5 + utrise/24 - (this.J - this.lng) * 4 / 60 / 24;
}
let dt = 0; //地方平太阳时 减 真太阳时 的差值,即"真平太阳时差换算表",单位为天
let tset = (LX == 2) ? utset as number : 0; //用于返回标准时间,关于月亮的必须先通过太阳升降获取到dt再转标准时间
let trise = (LX == 2) ? utrise as number : 0;
if((LX == 2) && (rise * sett == 1)){ //太阳相关,非极昼极夜且有升有落
while(tset < trise){ //太阳先落下再升起,时区与经度不匹配的情况下会出现此种情况,加一天修正
tset += 1;
}
dt = Math.round(jd) - (trise + (tset - trise) / 2); //单位为天.比较两者的中午12点(上午和下午是对称的)
tset = tset - dt + (this.J - this.lng) * 4 / 60 / 24; //真太阳时转标准时间
trise = trise - dt + (this.J - this.lng) * 4 / 60 / 24;
}
return [utrise, utset, dt, trise, tset];
}
/**
* 真太阳时模块,改编自 https://bieyu.com/ (月亮与太阳出没时间)
* 原理:用天文方法计算出太阳升起和落下时刻,中间则为当地正午(自创),与12点比较得到时差;与寿星万年历比较,两者相差在20秒内
* @param float jd 公历時间的儒略日值
* @param float J 经度,东经为正西经为负,注意西经60度38分转换方式是: -60 + -1 * 38/60
* @param float W 纬度,北纬为正南纬为负,太阳并不是严格从正东方升起,所以纬度也有影响,只是相对影响较小
*/
private zty(jd:number):number {
let dt = this.risenset(jd, 2)[2] as number;
return jd - (this.J - this.lng) * 4 / 60 / 24 + dt; //转地方平太阳时+修正
}
}
`export class SolarTimeUtil{
/** * 标准时间发出地经度(角度表示,东经为正西经为负),北京时间的经度为+120度0分 */ private J = 120.0; /** * 默认纬度(角度表示,北纬为正南纬为负),这里是中国标准时间发出地(陕西省渭南市蒲城县) */ private W = 35.0; /** 统一东经为正 */ private lng = 120.0; /** 北纬为正,南纬为负 */ private lat = 35.0; private constructor(lng:number=120.0,lat:number=35.0) { if (lng < -180.0 || lng > 180.0) { throw new Error(`illegal longitude: ${lng}`); } if (lat < -90.0 || lat > 90.0) { throw new Error(`illegal latitude: ${lat}`); } this.lng = lng; this.lat = lat; } /** * @param lng 经度 -180~180 * @param lat 纬度 -90~90 */ static initLocation(lng:number=120.0,lat:number=35.0){ return new SolarTimeUtil(lng,lat); } /** * 计算平太阳时 * @param solartime 公历时间 * @returns 平太阳时的公历时间 SolarTime */ getMeanSolarTime(solartime:SolarTime){ let spcjd = solartime.getJulianDay().getDay() ; // 將公历時间转换为儒略日 special jd,这里依然是标准时间,即this.J处的平太阳时 let deltPty = spcjd - (this.J - this.lng) * 4 / 60 / 24; // 计算地方平太阳时,每经度时差4分钟 let pty = JulianDay.fromJulianDay(deltPty).getSolarTime(); return pty; } /** * 计算真太阳时 * @param solartime 公历时间 * @returns 真太阳时的公历时间 SolarTime */ getRealSolarTime(solartime:SolarTime):SolarTime{ let spcjd = solartime.getJulianDay().getDay(); // 將公历時间转换为儒略日 special jd,这里依然是标准时间,即this.J处的平太阳时 let realZty = this.zty(spcjd); let zty = JulianDay.fromJulianDay(realZty).getSolarTime(); return zty; } /** * 真太阳时模块,sn代表sin */ private sn(x:number) { return Math.sin(x * 1.74532925199433E-02); } /** * 真太阳时模块,cn代表cosine */ private cn(x:number) { return Math.cos(x * 1.74532925199433E-02); } /** * 真太阳时模块,返回小数部分(负数特殊) returns fractional part of a number */ private fpart(x:number) { x = x - Math.floor(x); if (x < 0) { x = x + 1; } return x; //只取小数部份 } /** * 真太阳时模块,只取整数部份 */ private ipart(x:number) { if (x == 0) { return 0; } return (x / Math.abs(x)) * Math.floor(Math.abs(x)); } /** * 真太阳时模块 * finds a parabola through three points and * returns values of coordinates of extreme value (xe, ye) and * zeros if any (z1, z2) * assumes that the x values are -1, 0, +1 */ private quad(ym:number, y0:number, yp:number) { let nz = 0; let A = 0.5 * (ym + yp) - y0; let b = 0.5 * (yp - ym); let c = y0; let xe = -b / (2 * A); //x coord of symmetry line let ye = (A * xe + b) * xe + c; //extreme value for y in interval let dis = b * b - 4 * A * c; //discriminant let z1:number = 0.0; let z2:number = 0.0; if (dis > 0) { //there are zeros let dx = 0.5 * Math.sqrt(dis) / Math.abs(A); z1 = xe - dx; z2 = xe + dx; if (Math.abs(z1) <= 1) { nz = nz + 1; } //This zero is in interval if (Math.abs(z2) <= 1) { nz = nz + 1; } //This zero is in interval if (z1 < -1) { z1 = z2; } } return [xe, ye, z1, z2, nz]; } /** * 真太阳时模块 * returns sine of the altitude of either the sun or the moon given the modified julian day of the UT * @param float jd 公历時间的儒略日 * @param int LX 1月亮 2太阳日升日落 3太阳海上微光 */ private sinalt(jd:number, LX:number):number { let instant = jd - 2400001.0; let t = (instant - 51544.5) / 36525; //减51544.5为相对2000年01月01日零点 let ra:number,dec:number; if (LX == 1) { let moon = this.moon(t); ra = moon[0]; dec = moon[1]; } else { let sun = this.sun(t); ra = sun[0]; dec = sun[1]; } let mjd0 = this.ipart(instant); //UT时间0点;returns the local sidereal time(计算观测地区的恒星时)开始 let ut = (instant - mjd0) * 24; let t2 = (mjd0 - 51544.5) / 36525; let gmst = 6.697374558 + 1.0027379093 * ut; gmst = gmst + (8640184.812866 + (0.093104 - 0.0000062 * t2) * t2) * t2 / 3600; let lmst = 24 * this.fpart((gmst + this.lng / 15) / 24); //结束 let tau = 15 * (lmst - ra); //hour angle of object return this.sn(this.lat) * this.sn(dec) + this.cn(this.lat) * this.cn(dec) * this.cn(tau); } /** * 真太阳时模块,关于太阳的 * Returns RA and DEC of Sun to roughly 1 arcmin for few hundred years either side of J2000.0 */ private sun(t:number) { let p2 = 2 * Math.PI; let COSEPS = 0.91748; let SINEPS = 0.39778; let m = p2 * this.fpart(0.993133 + 99.997361 * t); //Mean anomaly let dL = 6893 * Math.sin(m) + 72 * Math.sin(2 * m); //Eq centre let L = p2 * this.fpart(0.7859453 + m / p2 + (6191.2 * t + dL) / 1296000); //convert to RA and DEC - ecliptic latitude of Sun taken as zero let sl = Math.sin(L); let x = Math.cos(L); let y = COSEPS * sl; let Z = SINEPS * sl; let rho = Math.sqrt(1 - Z * Z); let dec = (360 / p2) * Math.atan(Z / rho); let ra = (48 / p2) * Math.atan(y / (x + rho)); if (ra < 0) { ra = ra + 24; } return [ra, dec]; } /** * 真太阳时模块,关于月球的,Returns RA and DEC of Moon to 5 arc min (ra) and 1 arc min (dec) for a few centuries either side of J2000.0 * Predicts rise and set times to within minutes for about 500 years in past - TDT and UT time diference may become significant for long times */ private moon(t:number) { let p2 = 2 * Math.PI; let ARC = 206264.8062; let COSEPS = 0.91748; let SINEPS = 0.39778; let L0 = this.fpart(0.606433 + 1336.855225 * t); //mean long Moon in revs let L = p2 * this.fpart(0.374897 + 1325.55241 * t); //mean anomaly of Moon let LS = p2 * this.fpart(0.993133 + 99.997361 * t); //mean anomaly of Sun let d = p2 * this.fpart(0.827361 + 1236.853086 * t); //diff longitude sun and moon let F = p2 * this.fpart(0.259086 + 1342.227825 * t); //mean arg latitude //longitude correction terms let dL = 22640 * Math.sin(L) - 4586 * Math.sin(L - 2 * d); dL = dL + 2370 * Math.sin(2 * d) + 769 * Math.sin(2 * L); dL = dL - 668 * Math.sin(LS) - 412 * Math.sin(2 * F); dL = dL - 212 * Math.sin(2 * L - 2 * d) - 206 * Math.sin(L + LS - 2 * d); dL = dL + 192 * Math.sin(L + 2 * d) - 165 * Math.sin(LS - 2 * d); dL = dL - 125 * Math.sin(d) - 110 * Math.sin(L + LS); dL = dL + 148 * Math.sin(L - LS) - 55 * Math.sin(2 * F - 2 * d); //latitude arguments let S = F + (dL + 412 * Math.sin(2 * F) + 541 * Math.sin(LS)) / ARC; let h = F - 2 * d; //latitude correction terms let N:number = -526 * Math.sin(h) + 44 * Math.sin(L + h) - 31 * Math.sin(h - L) - 23 * Math.sin(LS + h); N = N + 11 * Math.sin(h - LS) - 25 * Math.sin(F - 2 * L) + 21 * Math.sin(F - L); let lmoon = p2 * this.fpart(L0 + dL / 1296000); //Lat in rads let bmoon = (18520 * Math.sin(S) + N) / ARC; //long in rads //convert to equatorial coords using a fixed ecliptic let CB = Math.cos(bmoon); let x = CB * Math.cos(lmoon); let V = CB * Math.sin(lmoon); let C = Math.sin(bmoon); let y = COSEPS * V - SINEPS * C; let Z = SINEPS * V + COSEPS * C; let rho = Math.sqrt(1 - Z * Z); let dec = (360 / p2) * Math.atan(Z / rho); //算出月球的视赤纬(apparent declination) let ra = (48 / p2) * Math.atan(y / (x + rho)); //算出月球的视赤经(apparent right ascension) if (ra < 0) { ra = ra + 24; } return [ra, dec]; } /** * 真太阳时模块,rise and set(升降计算) [升起时刻(真太阳时),落下时刻(真太阳时),真平太阳时差(仅类型2),升起时刻(标准时间,仅类型2),落下时刻(标准时间,仅类型2)] * @param float jd 公历時间的儒略日值 * @param int LX 类型:1月亮;2太阳日升日落;3太阳海上微光 * @return array */ private risenset(jd:number, LX:number) { let noon = Math.round(jd) - this.J / 360.0; //儒略日,中午12点,減去8小時時差 let sinho:number[] = [] //太阳盘面几何中心与理想地平面之间的夹角 sinho[1] = this.sn(8 / 60); //moonrise - average diameter used(月亮升降) sinho[2] = this.sn(-50 / 60); //sunrise - classic value for refraction(太阳升降) sinho[3] = this.sn(-12); //nautical twilight(海上微光) let rise = 0; //是否有升起动作 let utrise:number|boolean = false; //升起的时间 let sett = 0; //是否有落下动作 let utset:number|boolean = false; //落下的时间 let hour = 1; let ym = this.sinalt(noon + (hour - 1)/24, LX) - sinho[LX]; //See STEP 1 and 2 of Web page description. let zero2 = 0; //两小时内是否进行了升起和落下两个动作(极地附近有这种情况,如1999年12月25日,经度0,纬度67.43,当天的太阳只有8分钟-_-) let above = ym > 0 ? 1 : 0; //used later to classify non-risings 是否在地平线上方,用于判断极昼极夜 do { //STEP 1 and STEP 3 of Web page description let y0 = this.sinalt(noon + (hour + 0)/24, LX) - sinho[LX]; let yp = this.sinalt(noon + (hour + 1)/24, LX) - sinho[LX]; //STEP 4 of web page description let quad = this.quad(ym, y0, yp); let ye = quad[1]; let z1 = quad[2]; let z2 = quad[3]; let nz = quad[4]; switch (nz) { //cases depend on values of discriminant - inner part of STEP 4 case 0: //nothing - go to next time slot break; case 1: //simple rise / set event if (ym < 0) { //must be a rising event utrise = hour + z1; rise = 1; } else { //must be setting utset = hour + z1; sett = 1; } break; case 2: //rises and sets within interval if (ye < 0) { //minimum - so set then rise utrise = hour + z2; utset = hour + z1; } else { //maximum - so rise then set utrise = hour + z1; utset = hour + z2; } rise = 1; sett = 1; zero2 = 1; break; } ym = yp; //reuse the ordinate in the next interval hour = hour + 2; } while (!((hour == 25) || (rise * sett == 1))); //STEP 5 of Web page description - have we finished for this object? if(utset !== false){ //注意这里转成了真太阳时 utset = Math.round(jd) - 0.5 + utset/24 - (this.J - this.lng) * 4 / 60 / 24; } if(utrise !== false){ utrise = Math.round(jd) - 0.5 + utrise/24 - (this.J - this.lng) * 4 / 60 / 24; } let dt = 0; //地方平太阳时 减 真太阳时 的差值,即"真平太阳时差换算表",单位为天 let tset = (LX == 2) ? utset as number : 0; //用于返回标准时间,关于月亮的必须先通过太阳升降获取到dt再转标准时间 let trise = (LX == 2) ? utrise as number : 0; if((LX == 2) && (rise * sett == 1)){ //太阳相关,非极昼极夜且有升有落 while(tset < trise){ //太阳先落下再升起,时区与经度不匹配的情况下会出现此种情况,加一天修正 tset += 1; } dt = Math.round(jd) - (trise + (tset - trise) / 2); //单位为天.比较两者的中午12点(上午和下午是对称的) tset = tset - dt + (this.J - this.lng) * 4 / 60 / 24; //真太阳时转标准时间 trise = trise - dt + (this.J - this.lng) * 4 / 60 / 24; } return [utrise, utset, dt, trise, tset]; } /** * 真太阳时模块,改编自 https://bieyu.com/ (月亮与太阳出没时间) * 原理:用天文方法计算出太阳升起和落下时刻,中间则为当地正午(自创),与12点比较得到时差;与寿星万年历比较,两者相差在20秒内 * @param float jd 公历時间的儒略日值 * @param float J 经度,东经为正西经为负,注意西经60度38分转换方式是: -60 + -1 * 38/60 * @param float W 纬度,北纬为正南纬为负,太阳并不是严格从正东方升起,所以纬度也有影响,只是相对影响较小 */ private zty(jd:number):number { let dt = this.risenset(jd, 2)[2] as number; return jd - (this.J - this.lng) * 4 / 60 / 24 + dt; //转地方平太阳时+修正 }}`
似乎不太完善
export class SolarTimeUtil{ /** * 标准时间发出地经度(角度表示,东经为正西经为负),北京时间的经度为+120度0分 */ private J = 120.0; /** * 默认纬度(角度表示,北纬为正南纬为负),这里是中国标准时间发出地(陕西省渭南市蒲城县) */ private W = 35.0; /** 统一东经为正 */ private lng = 120.0; /** 北纬为正,南纬为负 */ private lat = 35.0; private constructor(lng:number=120.0,lat:number=35.0) { if (lng < -180.0 || lng > 180.0) { throw new Error(`illegal longitude: ${lng}`); } if (lat < -90.0 || lat > 90.0) { throw new Error(`illegal latitude: ${lat}`); } this.lng = lng; this.lat = lat; } /** * @param lng 经度 -180~180 * @param lat 纬度 -90~90 */ static initLocation(lng:number=120.0,lat:number=35.0){ return new SolarTimeUtil(lng,lat); } /** * 计算平太阳时 * @param solartime 公历时间 * @returns 平太阳时的公历时间 SolarTime */ getMeanSolarTime(solartime:SolarTime){ let spcjd = solartime.getJulianDay().getDay() ; // 將公历時间转换为儒略日 special jd,这里依然是标准时间,即this.J处的平太阳时 let deltPty = spcjd - (this.J - this.lng) * 4 / 60 / 24; // 计算地方平太阳时,每经度时差4分钟 let pty = JulianDay.fromJulianDay(deltPty).getSolarTime(); return pty; } /** * 计算真太阳时 * @param solartime 公历时间 * @returns 真太阳时的公历时间 SolarTime */ getRealSolarTime(solartime:SolarTime):SolarTime{ let spcjd = solartime.getJulianDay().getDay(); // 將公历時间转换为儒略日 special jd,这里依然是标准时间,即this.J处的平太阳时 let realZty = this.zty(spcjd); let zty = JulianDay.fromJulianDay(realZty).getSolarTime(); return zty; } /** * 真太阳时模块,sn代表sin */ private sn(x:number) { return Math.sin(x * 1.74532925199433E-02); } /** * 真太阳时模块,cn代表cosine */ private cn(x:number) { return Math.cos(x * 1.74532925199433E-02); } /** * 真太阳时模块,返回小数部分(负数特殊) returns fractional part of a number */ private fpart(x:number) { x = x - Math.floor(x); if (x < 0) { x = x + 1; } return x; //只取小数部份 } /** * 真太阳时模块,只取整数部份 */ private ipart(x:number) { if (x == 0) { return 0; } return (x / Math.abs(x)) * Math.floor(Math.abs(x)); } /** * 真太阳时模块 * finds a parabola through three points and * returns values of coordinates of extreme value (xe, ye) and * zeros if any (z1, z2) * assumes that the x values are -1, 0, +1 */ private quad(ym:number, y0:number, yp:number) { let nz = 0; let A = 0.5 * (ym + yp) - y0; let b = 0.5 * (yp - ym); let c = y0; let xe = -b / (2 * A); //x coord of symmetry line let ye = (A * xe + b) * xe + c; //extreme value for y in interval let dis = b * b - 4 * A * c; //discriminant let z1:number = 0.0; let z2:number = 0.0; if (dis > 0) { //there are zeros let dx = 0.5 * Math.sqrt(dis) / Math.abs(A); z1 = xe - dx; z2 = xe + dx; if (Math.abs(z1) <= 1) { nz = nz + 1; } //This zero is in interval if (Math.abs(z2) <= 1) { nz = nz + 1; } //This zero is in interval if (z1 < -1) { z1 = z2; } } return [xe, ye, z1, z2, nz]; } /** * 真太阳时模块 * returns sine of the altitude of either the sun or the moon given the modified julian day of the UT * @param float jd 公历時间的儒略日 * @param int LX 1月亮 2太阳日升日落 3太阳海上微光 */ private sinalt(jd:number, LX:number):number { let instant = jd - 2400001.0; let t = (instant - 51544.5) / 36525; //减51544.5为相对2000年01月01日零点 let ra:number,dec:number; if (LX == 1) { let moon = this.moon(t); ra = moon[0]; dec = moon[1]; } else { let sun = this.sun(t); ra = sun[0]; dec = sun[1]; } let mjd0 = this.ipart(instant); //UT时间0点;returns the local sidereal time(计算观测地区的恒星时)开始 let ut = (instant - mjd0) * 24; let t2 = (mjd0 - 51544.5) / 36525; let gmst = 6.697374558 + 1.0027379093 * ut; gmst = gmst + (8640184.812866 + (0.093104 - 0.0000062 * t2) * t2) * t2 / 3600; let lmst = 24 * this.fpart((gmst + this.lng / 15) / 24); //结束 let tau = 15 * (lmst - ra); //hour angle of object return this.sn(this.lat) * this.sn(dec) + this.cn(this.lat) * this.cn(dec) * this.cn(tau); } /** * 真太阳时模块,关于太阳的 * Returns RA and DEC of Sun to roughly 1 arcmin for few hundred years either side of J2000.0 */ private sun(t:number) { let p2 = 2 * Math.PI; let COSEPS = 0.91748; let SINEPS = 0.39778; let m = p2 * this.fpart(0.993133 + 99.997361 * t); //Mean anomaly let dL = 6893 * Math.sin(m) + 72 * Math.sin(2 * m); //Eq centre let L = p2 * this.fpart(0.7859453 + m / p2 + (6191.2 * t + dL) / 1296000); //convert to RA and DEC - ecliptic latitude of Sun taken as zero let sl = Math.sin(L); let x = Math.cos(L); let y = COSEPS * sl; let Z = SINEPS * sl; let rho = Math.sqrt(1 - Z * Z); let dec = (360 / p2) * Math.atan(Z / rho); let ra = (48 / p2) * Math.atan(y / (x + rho)); if (ra < 0) { ra = ra + 24; } return [ra, dec]; } /** * 真太阳时模块,关于月球的,Returns RA and DEC of Moon to 5 arc min (ra) and 1 arc min (dec) for a few centuries either side of J2000.0 * Predicts rise and set times to within minutes for about 500 years in past - TDT and UT time diference may become significant for long times */ private moon(t:number) { let p2 = 2 * Math.PI; let ARC = 206264.8062; let COSEPS = 0.91748; let SINEPS = 0.39778; let L0 = this.fpart(0.606433 + 1336.855225 * t); //mean long Moon in revs let L = p2 * this.fpart(0.374897 + 1325.55241 * t); //mean anomaly of Moon let LS = p2 * this.fpart(0.993133 + 99.997361 * t); //mean anomaly of Sun let d = p2 * this.fpart(0.827361 + 1236.853086 * t); //diff longitude sun and moon let F = p2 * this.fpart(0.259086 + 1342.227825 * t); //mean arg latitude //longitude correction terms let dL = 22640 * Math.sin(L) - 4586 * Math.sin(L - 2 * d); dL = dL + 2370 * Math.sin(2 * d) + 769 * Math.sin(2 * L); dL = dL - 668 * Math.sin(LS) - 412 * Math.sin(2 * F); dL = dL - 212 * Math.sin(2 * L - 2 * d) - 206 * Math.sin(L + LS - 2 * d); dL = dL + 192 * Math.sin(L + 2 * d) - 165 * Math.sin(LS - 2 * d); dL = dL - 125 * Math.sin(d) - 110 * Math.sin(L + LS); dL = dL + 148 * Math.sin(L - LS) - 55 * Math.sin(2 * F - 2 * d); //latitude arguments let S = F + (dL + 412 * Math.sin(2 * F) + 541 * Math.sin(LS)) / ARC; let h = F - 2 * d; //latitude correction terms let N:number = -526 * Math.sin(h) + 44 * Math.sin(L + h) - 31 * Math.sin(h - L) - 23 * Math.sin(LS + h); N = N + 11 * Math.sin(h - LS) - 25 * Math.sin(F - 2 * L) + 21 * Math.sin(F - L); let lmoon = p2 * this.fpart(L0 + dL / 1296000); //Lat in rads let bmoon = (18520 * Math.sin(S) + N) / ARC; //long in rads //convert to equatorial coords using a fixed ecliptic let CB = Math.cos(bmoon); let x = CB * Math.cos(lmoon); let V = CB * Math.sin(lmoon); let C = Math.sin(bmoon); let y = COSEPS * V - SINEPS * C; let Z = SINEPS * V + COSEPS * C; let rho = Math.sqrt(1 - Z * Z); let dec = (360 / p2) * Math.atan(Z / rho); //算出月球的视赤纬(apparent declination) let ra = (48 / p2) * Math.atan(y / (x + rho)); //算出月球的视赤经(apparent right ascension) if (ra < 0) { ra = ra + 24; } return [ra, dec]; } /** * 真太阳时模块,rise and set(升降计算) [升起时刻(真太阳时),落下时刻(真太阳时),真平太阳时差(仅类型2),升起时刻(标准时间,仅类型2),落下时刻(标准时间,仅类型2)] * @param float jd 公历時间的儒略日值 * @param int LX 类型:1月亮;2太阳日升日落;3太阳海上微光 * @return array */ private risenset(jd:number, LX:number) { let noon = Math.round(jd) - this.J / 360.0; //儒略日,中午12点,減去8小時時差 let sinho:number[] = [] //太阳盘面几何中心与理想地平面之间的夹角 sinho[1] = this.sn(8 / 60); //moonrise - average diameter used(月亮升降) sinho[2] = this.sn(-50 / 60); //sunrise - classic value for refraction(太阳升降) sinho[3] = this.sn(-12); //nautical twilight(海上微光) let rise = 0; //是否有升起动作 let utrise:number|boolean = false; //升起的时间 let sett = 0; //是否有落下动作 let utset:number|boolean = false; //落下的时间 let hour = 1; let ym = this.sinalt(noon + (hour - 1)/24, LX) - sinho[LX]; //See STEP 1 and 2 of Web page description. let zero2 = 0; //两小时内是否进行了升起和落下两个动作(极地附近有这种情况,如1999年12月25日,经度0,纬度67.43,当天的太阳只有8分钟-_-) let above = ym > 0 ? 1 : 0; //used later to classify non-risings 是否在地平线上方,用于判断极昼极夜 do { //STEP 1 and STEP 3 of Web page description let y0 = this.sinalt(noon + (hour + 0)/24, LX) - sinho[LX]; let yp = this.sinalt(noon + (hour + 1)/24, LX) - sinho[LX]; //STEP 4 of web page description let quad = this.quad(ym, y0, yp); let ye = quad[1]; let z1 = quad[2]; let z2 = quad[3]; let nz = quad[4]; switch (nz) { //cases depend on values of discriminant - inner part of STEP 4 case 0: //nothing - go to next time slot break; case 1: //simple rise / set event if (ym < 0) { //must be a rising event utrise = hour + z1; rise = 1; } else { //must be setting utset = hour + z1; sett = 1; } break; case 2: //rises and sets within interval if (ye < 0) { //minimum - so set then rise utrise = hour + z2; utset = hour + z1; } else { //maximum - so rise then set utrise = hour + z1; utset = hour + z2; } rise = 1; sett = 1; zero2 = 1; break; } ym = yp; //reuse the ordinate in the next interval hour = hour + 2; } while (!((hour == 25) || (rise * sett == 1))); //STEP 5 of Web page description - have we finished for this object? if(utset !== false){ //注意这里转成了真太阳时 utset = Math.round(jd) - 0.5 + utset/24 - (this.J - this.lng) * 4 / 60 / 24; } if(utrise !== false){ utrise = Math.round(jd) - 0.5 + utrise/24 - (this.J - this.lng) * 4 / 60 / 24; } let dt = 0; //地方平太阳时 减 真太阳时 的差值,即"真平太阳时差换算表",单位为天 let tset = (LX == 2) ? utset as number : 0; //用于返回标准时间,关于月亮的必须先通过太阳升降获取到dt再转标准时间 let trise = (LX == 2) ? utrise as number : 0; if((LX == 2) && (rise * sett == 1)){ //太阳相关,非极昼极夜且有升有落 while(tset < trise){ //太阳先落下再升起,时区与经度不匹配的情况下会出现此种情况,加一天修正 tset += 1; } dt = Math.round(jd) - (trise + (tset - trise) / 2); //单位为天.比较两者的中午12点(上午和下午是对称的) tset = tset - dt + (this.J - this.lng) * 4 / 60 / 24; //真太阳时转标准时间 trise = trise - dt + (this.J - this.lng) * 4 / 60 / 24; } return [utrise, utset, dt, trise, tset]; } /** * 真太阳时模块,改编自 https://bieyu.com/ (月亮与太阳出没时间) * 原理:用天文方法计算出太阳升起和落下时刻,中间则为当地正午(自创),与12点比较得到时差;与寿星万年历比较,两者相差在20秒内 * @param float jd 公历時间的儒略日值 * @param float J 经度,东经为正西经为负,注意西经60度38分转换方式是: -60 + -1 * 38/60 * @param float W 纬度,北纬为正南纬为负,太阳并不是严格从正东方升起,所以纬度也有影响,只是相对影响较小 */ private zty(jd:number):number { let dt = this.risenset(jd, 2)[2] as number; return jd - (this.J - this.lng) * 4 / 60 / 24 + dt; //转地方平太阳时+修正 } }似乎不太完善
有什么问题呢?
算法从这个项目提取修改的 https://github.com/hkargc/paipan