计算太阳在地平线以下/上方X度的时间 [英] calculate the time when the sun is X degrees below/above the Horizon
问题描述
我想知道太阳在地平线以下/上方X度的时间.
例如，我想找到太阳在地平线以下19.75度的时间.我认为这与date_sunrise/date_sunset函数中的最高点有关，但我不确定.
提前谢谢！
-
收集所需的太阳星历表数据
采取1个小时的步骤，以所需的地理位置获取太阳在方位角坐标中的位置.可以使用找到的方程式，也可以使用某些WEB服务，例如:
- JPL地平线不喜欢这一点，因为它们使用了一些奇怪的输出参考帧与我的测量值不符，但很可能会在此过程中发生错误……
-
Presov天文台这是我的最爱(但它在斯洛伐克)输出可以很容易地复制到矿井发动机上，并且输出结果与矿井的观测，计算，估计和测量结果相对应.只需填写:
- 地理位置(Miesto pozorovania)
- 日期，时间(Dátumačaspozorovania)
- 在左下方:时间间隔[天]，间隔步长[天]
- 单击太阳(Slnko)，月亮(Mesiac)，行星(Planety)的按钮
那里有很多这样的页面，但是总是检查它们是否输出正确的数据.我使用开普勒定律/方程构成行星运动(精度较低，但对于地球太阳应该可以).如今，引擎改用重力模型了(但是随着时间的流逝，它会变得很不稳定)
-
将数据作为多段线(方位角，高度，时间)上的3D点进行处理
-
现在仅在数据中找到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.
[Edit1]开普勒定律
如果您要采用这种方式，请在此处.您将需要地球的轨道和旋转参数.这些是从雷星历表引擎* .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，纹理，环和大气参数.所以:
- 让Sun在笛卡尔坐标中定位
(0,0,0) -
根据开普勒定律计算地球位置
(x,y,z)然后太阳在地心坐标中为
(-x,-y,-z) -
通过每日轮换，进动，章动回滚
(-x,-y,-z) -> (x',y',z') - 计算地理位置
(North,East,High(Up))的 - 将
(x',y',z')转换为NEH本地坐标(xx,yy,zz) -
计算:
azimut=atanxy(-xx,-yy) height=atanxy(sqrt((xx*xx)+(yy*yy)),-zz)就是这样
NEH框架
这是我的日心星体星体位置计算:
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)=inverse(rep)*(gx,gy,gz)
local是地球的坐标系和全局Sun的坐标系
I want to know what is the time when the sun is X degrees below/above the Horizon.
For example, I want to find the time when the sun is 19.75 degrees below the horizon. I think it has something to do with the zenith in the function date_sunrise/date_sunset but I'm not sure.
Thanks in advance!
collect Sun ephemerids data for day you need
take 1 hour steps and get Suns position in azimuthal coordinates for the geo-location you need. Either use equations you found or use some WEB service like:
- JPL Horizons don't like this one as they use some weird output reference frames that does not correspond to my measurements, but it is more likely I transform something wrong along the way ...
Presov observatory this is mine favorite (but it is in Slovak) output is easily copyable to mine engines and the output is corresponding with mine observations, computations, estimations and measurements. Just fill in the:
- geo-location (Miesto pozorovania)
- date,time (Dátum a čas pozorovania)
- in the bottom from left: interval[days],interval step[days]
- click on the button for Sun(Slnko), Moon(Mesiac),Planets(Planety)
there are many such pages out there just look but always check if they outputting correct data. I use Kepler's laws/equations form planetary motions (lower precision but for Earth-Sun should be OK). Nowadays engines use gravity model instead (but that is unstable with longer times from epoch)
handle the data as set of 3D points along polyline (azimut,height,time)
now just find in the data 2 points
one below desired angle and the next above desired angle. Booth points must be neighboring. If any point is on the desired angle then you already have the solution so stop
interpolate the height angle crossing time
so if desired height angle is
band wanted timetthen:a0,a1are azimutal anglesb0,b1are height anglest0,t1are times
then just solve this linear system:
b=b0+(b1-b0)*u t=t0+(t1-t0)*uso if I did not make some silly mistake:
t=t0+((t1-t0)*(b-b0)/(b1-b0))
[Notes]
if you do not need too high precision (and usage above 100 years) and geo-location is fixed then you can table whole year and use this data periodically. This way you will not need step 1 on runtime.
[Edit1] Kepler's law
if you want to go this way look here. You will need orbital and rotation parameters of Earth. These are extracted from mine ephemeris engine *.ini for solar system:
[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
and here are the explanations:
[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
Ignore the is_star, textures, rings and atmosphere parameters. So:
- get Sun to position
(0,0,0)in cartesian coordinates compute Earth position
(x,y,z)from Kepler's lawSun is then
(-x,-y,-z)in geocentric coordinatesrotate back by dayly rotation,precession,nutation
(-x,-y,-z) -> (x',y',z')- compute
NEHframe for your geolocation(North,East,High(Up)) - convert
(x',y',z')to NEH local coordinates(xx,yy,zz) compute:
azimut=atanxy(-xx,-yy) height=atanxy(sqrt((xx*xx)+(yy*yy)),-zz)and that is it
here is mine Heliocentric astro body position computation:
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 class is quite complicated (something like GLM) the only thing you need from it are local rotations (all other stuff is basic). this is how 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);
repis the input and output matrixangis the rotation angle [rad]
now you can use rep to convert to/from Earth local coordinate system
- LCS to GCS
(l2g) ... (gx,gy,gz)=rep*(lx,ly,lz) - GCS to LCS
(g2l) ... (lx,ly,lz)=inverse(rep)*(gx,gy,gz)
local is Earth's coordinate system and global Sun's coordinate system
这篇关于计算太阳在地平线以下/上方X度的时间的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
