php – 计算太阳低于/高于地平线X度的时间
作者：互联网
我想知道太阳低于/高于地平线X度的时间是什么时候.
例如,我想找到太阳在地平线以下19.75度的时间.我认为它与函数date_sunrise / date_sunset中的天顶有关,但我不确定.
提前致谢！
解决方法:
>收集所需日期的太阳星历表数据
步行1小时,在方位坐标处获得太阳位置,以获得所需的地理位置.使用您找到的方程式或使用某些WEB服务,例如：
> JPL Horizons不喜欢这个,因为他们使用了一些与我的测量不对应的奇怪的输出参考帧,但更可能的是我在路上改变了一些错误…
> Presov observatory这是我最喜欢的(但它在斯洛伐克)输出很容易复制到矿井引擎,输出与矿山观测,计算,估算和测量相对应.只需填写：
>地理位置(Miesto pozorovania)
>日期,时间(Dátumačaspozorovania)
>从左下方：间隔[天],间隔步[天]
>点击太阳(Slnko),月亮(Mesiac),行星(Planety)按钮
有许多这样的页面只是看,但总是检查他们是否输出正确的数据.我使用开普勒的定律/方程形成行星运动(精度较低,但对于地球 – 太阳应该没问题).现在发动机使用重力模型代替(但是从纪元开始的时间越长,它就越不稳定)
>将数据作为沿折线的三维点集合处理(方位角,高度,时间)
>现在只需在数据中找到2分
一个低于所需角度,另一个高于所需角度.展位必须相邻.如果任何一点在所需的角度,那么你已经有了解决方案,所以停止
>插入高度角交叉时间
所以如果需要的高度角是b,想要的时间t那么：
> a0,a1是方位角
> b0,b1是高度角
> t0,t1是次
然后就解决这个线性系统：
b=b0+(b1-b0)*u
t=t0+(t1-t0)*u
所以,如果我没有犯一些愚蠢的错误：
t=t0+((t1-t0)*(b-b0)/(b1-b0))
[笔记]
如果你不需要太高的精度(并且使用超过100年)并且地理位置是固定的,那么你可以整年制表并定期使用这些数据.这样,您就不需要在运行时执行第1步.
[编辑1]开普勒定律
如果你想这样看here.你将需要地球的轨道和旋转参数.这些是从太阳系的矿山星历引擎* .ini中提取的：
[Earth]
txr_map=Earth_Map.jpg
txr_nor=Earth_Normal.jpg
txr_clouds=Earth_Cloud.jpg
txr_lights=Earth_Light.jpg
txr_ring_map=
txr_ring_alpha=
is_star=0
mother=Sun
re=6378141.2
rp=6356754.79506139
r0=-1
r1=-1
ha=60000
vd=250000
B0r=0.1981
B0g=0.4656
B0b=0.8625
B0a=0.75
t0=-0.0833333333333333 ; this means 1.1.2000 00:00:00 UT
a=149597896927.617
da=-0.122872993839836
e=0.01673163
de=-1.00232717316906E-9
i=-9.48516635288838E-6
di=-6.38963964003634E-9
O=-0.004695
dO=-1.15274665428334E-7
o=1.79646842620403
do=1.51932094052745E-7
M =1.7464
dM =0.0172021242603194
ddM=0
rota0 =3.0707963267949
rotda =6.30038738085328
prea0 =1.5707963267949
preda =-6.68704522111755E-7
prei =0.409124584728753
predi =0
nuta =0
nutda =0
nutia =0
nutdia=0
nutii =0
nutdii=0
以下是解释：
[Name] [string id] object ID name
txr_map [filename] surface texture
txr_nor [filename] surface normal/bump texture
txr_clouds [filename] cloud blend texture (white cloud, black clear sky)
txr_lights [filename] night surface texture
txr_ring_map [filename] rings color texture
txr_ring_alpha [filename] rings alpha texture (alpha0 transparent, alpha1 solid)
is_star [0/1] is star ?
mother [string] "" or owner object name
re [m] equator radius
rp [m] polar radius
r0 [m] -1 or rings inner radius
r1 [m] -1 or rings outer radius
ha [m] 0 or atmosphere thickness
vd [m] -1 or atmosphere view depth
B0r <0,1> star R light or atmosphere color
B0g <0,1> star G light or atmosphere color
B0b <0,1> star B light or atmosphere color
B0a <0,1> overglow of star below horizont
t0 [day] t0 time the parameters are taken after 1.1.2000 00:00:00
a [m] a main semiaxis
da [m/day] a change in time
e [-] e eccentricity
de [-/day] e change in time
i [rad] i inclination
di [rad/day] i change in time
O [rad] O (node n) position of inclination axis
dO [rad/day] O node shift (pi2/T)
o [rad] o perihelium (no change in inclination position)
do [rad/day] o perihelium shift (pi2/T)
M [rad] M rotation around owner position in t0
dM [rad/day] dM orbital rotation (pi2/draconic month)
ddM0 [rad/day^2] dM change in time
rota0 [rad] rota0 rotation around self axis position in t0
rotda [rad/day] rotda mean rotation around self axis
prea0 [rad] prea rotation axis position in t0
preda [rad/day] preda precession rotation (pi2/Platonic year)
prei [rad] prei equator inclination to ecliptic
predi [rad/day] prei change in time
nuta [rad] nuta angle position on nutation ellipse
nutda [rad/day] nutation rotation (pi2/T)
nutia [rad] nutia nutation (of rotation axis) ellipse semiaxis axis in ecliptic plane
nutdia [rad/day] nutia change in time
nutii [rad] nutii nutation (of rotation axis) ellipse semiaxis axis in rotation axis direction
nutdii [rad/day] nutii change in time
忽略is_star,纹理,环和大气参数.所以：
>让太阳以笛卡尔坐标定位(0,0,0)
>根据开普勒定律计算地球位置(x,y,z)
然后,太阳是地心坐标中的(-x,-y,-z)
>通过每日旋转,进动,章动(-x,-y,-z)旋转 – > (X ‘Y’,Z’)
>为您的地理位置计算NEH帧(北,东,高(上))
>将(x’,y’,z’)转换为NEH局部坐标(xx,yy,zz)
>计算：
azimut=atanxy(-xx,-yy)
height=atanxy(sqrt((xx*xx)+(yy*yy)),-zz)
就是这样
这是我的Heliocentric astro体位计算：
void astro_body::compute(double t)
{
// t is time in days after 1.1.2000 00:00:00
// double pos[3] is output heliocentric position [m]
// reper rep is output heliocentric position [m] and orientation transform matrix (mine class) where Z is rotation axis (North pole) and X is long=0,lat=0
rot_a.compute(t); // compute actual value for orbital parameters changing in time
pre_a.compute(t); // the actual parameter is in XXX.a you can ignore this part
pre_i.compute(t);
nut_a.compute(t);
nut_ia.compute(t);
nut_ii.compute(t);
// pre_a=pre_a0+(pre_da.a*dt)+(nut_ia*cos(nut_a)); // some old legacy dead code
// pre_i=pre_i0+(pre_di.a*sin(pre_e))+(nut_ii*sin(nut_a));
rep.reset(); // rep is the transform matrix representing body coordinate system (orientation and position)
rep.lrotz(pre_a.a); // local rotation around reps Z axis by pre_a.a [rad] angle
rep.lroty(pre_i.a);
rep.lrotx(nut_ia.a*cos(nut_a.a));
rep.lroty(nut_ii.a*sin(nut_a.a));
rep.lrotz(rot_a.a);
a.compute(t); // the same as above can ignore this part
e.compute(t);
i.compute(t);
O.compute(t);
o.compute(t);
M.compute(t);
M.compute(t);
double c0,c1,c2,sO,si,cO,ci,b; // trajectory constants
double x,y;
int q;
if (e.a>=1.0) e.a=0;
c0=sqrt((1.0-e.a)/(1.0+e.a)); // some helper constants computation
c1=sqrt((1.0+e.a)/(1.0-e.a));
c2=a.a*(1-e.a*e.a);
sO=sin(O.a);
cO=cos(O.a);
si=sin(-i.a);
ci=cos(-i.a);
b=a.a*sqrt(1.0-e.a);
M.a-=o.a; // correction
M.a=M.a-pi2*floor(M.a/pi2);
E=M.a;
for (q=0;q<20;q++) E=M.a+e.a*sin(E); // Kepler's equation
V=2.0*atan(c1*tan(E/2.0));
r=c2/(1.0+e.a*cos(V));
pos[0]=r*cos(V+o.a-O.a); // now just compute heliocentric position along ecliptic ellipse
pos[1]=r*sin(V+o.a-O.a); // and then rotate by inclination
pos[2]=-pos[1]*si;
pos[1]=+pos[1]*ci;
x=pos[0]; y=pos[1];
pos[0]=x*cO-y*sO;
pos[1]=x*sO+y*cO;
if ((mother>=0)&&(tab!=NULL)) vector_add(pos,pos,tab[mother].pos); // if satelite like Moon add owners position
rep.gpos_set(pos); // set the global position to transform matrix also
}
//---------------------------------------------------------------------------
reper类非常复杂(类似GLM)你需要的唯一东西是局部旋转(所有其他东西都是基本的).这就是lrotx的工作原理：
double c=cos(ang),s=sin(ang);
double rot[16],inv[16]; // rot is the rotation around x transform matrix
rot=(1, 0, 0, 0,
0, c,-s, 0,
0, s, c, 0,
0, 0, 0, 1);
inv=inverse(rep); // inverse is inverse matrix 4x4
inv=inv*rot
rep=inverse(inv);
> rep是输入和输出矩阵
> ang是旋转角度[rad]
现在您可以使用rep转换为/从地球局部坐标系转换
> LCS到GCS(l2g)…(gx,gy,gz)= rep *(lx,ly,lz)
> GCS到LCS(g2l)……(lx,ly,lz)=逆(rep)*(gx,gy,gz)
本地是地球的坐标系和全球太阳的坐标系
标签：php,astronomy 来源： https://codeday.me/bug/20190917/1808856.html