我想知道太阳在地平线以下/之上 X 度的时间是什么时候。
例如，我想找到太阳位于地平线以下 19.75 度的时间。我认为这与函数中的最高点有关date_sunrise/date_sunset但我不确定。
提前致谢！
采取 1 小时步长，获取太阳在方位角坐标中的位置，以获得您所需的地理位置。使用您找到的方程式或使用一些 WEB 服务，例如：
喷气推进实验室地平线不喜欢这个，因为他们使用了一些奇怪的输出参考框架，这些参考框架与我的测量不符，但更有可能的是我一路上转换了一些错误的东西......
普雷绍夫天文台这是我最喜欢的（但它是斯洛伐克语）输出可以轻松复制到矿井引擎，并且输出与矿井观察、计算、估计和测量相对应。只需填写：
有很多这样的页面，只是看看，但总是检查它们是否输出正确的数据。我使用开普勒定律/方程形成行星运动（精度较低，但对于地球-太阳应该没问题）。如今，引擎使用重力模型（但随着时间的推移，它不稳定）
将数据处理为沿折线的 3D 点集（方位角、高度、时间）
现在只要在数据中找到2点
一个低于所需角度，下一个高于所需角度。展位点必须相邻。如果任何点位于所需的角度，那么您已经有了解决方案，所以停止
所以如果所需的高度角是b想要时间t then:
a0,a1是方位角b0,b1是高度角t0,t1是次然后求解这个线性系统：
b=b0+(b1-b0)*u
t=t0+(t1-t0)*u
所以如果我没有犯一些愚蠢的错误：
t=t0+((t1-t0)*(b-b0)/(b1-b0))
[Notes]
如果你不需要太高的精度（以及上面的用法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')
compute NEH您的地理位置框架(North,East,High(Up))
convert (x',y',z')到 NEH 局部坐标(xx,yy,zz)
compute:
azimut=atanxy(-xx,-yy)
height=atanxy(sqrt((xx*xx)+(yy*yy)),-zz)
就是这样
这是我的日心天体身体位置计算：
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 works:
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与地球局部坐标系进行转换
(l2g) ... (gx,gy,gz)=rep*(lx,ly,lz)
(g2l) ... (lx,ly,lz)=inverse(rep)*(gx,gy,gz)
local 是地球坐标系，global 是太阳坐标系