计算太阳在地平线以下/以上X度的时间

我想知道什么时候太阳在地平线以下X度以上。

例如，我想find太阳在地平线以下19.75度的时间。 我认为这与函数date_sunrise / date_sunset但我不确定。

提前致谢！

  1. 收集您需要的一天的Sun ephemerids数据

    采取1小时的步骤，并获得所需的地理位置的方位angular坐标太阳位置。 您可以使用您find的公式或使用一些WEB服务：

    • JPL Horizo​​ns不喜欢这个，因为他们使用了一些奇怪的输出参考框架，这些参考框架与我的测量结果并不相符，但是更有可能在这个过程中改变了一些错误。

    • Presov天文台这是我最喜欢的（但它是在斯洛伐克）输出很容易复制到我的引擎和输出是相应的我的观测，计算，估计和测量。 只需填写：

      • 地理位置（Miesto pozorovania）
      • date，时间（Dátumačaspozorovania）
      • 在左边的底部：区间[days]，区间步[days]
      • 点击button为太阳（Slnko），月亮（Mesiac），行星（Planety）

    这里有很多这样的页面，但总是检查它们是否输出正确的数据。 我使用开普勒定律/方程forms的行星运动（精度较低，但地球 – 太阳应该没问题）。 现在引擎使用重力模型（但是这是不稳定的，从时代开始的时间更长）

  2. 将数据作为一组沿着折线（方位angular，高度，时间）的3D点处理

  3. 现在在数据中find2点

    一个低于所需的angular度，另一个低于所需的angular度。 展位点必须是相邻的。 如果任何一点是在所需的angular度，那么你已经有解决scheme，所以停止

  4. 内插高度angular穿越时间

    在这里输入图像描述

    所以如果期望的高度angular是b并且想要时间t那么：

    • a0 ， a1是方位angular
    • b0 ， b1是高度angular
    • t0 ， t1是次数

    那么就解决这个线性系统吧： 

     b=b0+(b1-b0)*u t=t0+(t1-t0)*u

    所以如果我没有犯一些愚蠢的错误： 

     t=t0+((t1-t0)*(b-b0)/(b1-b0))

[笔记]

如果您不需要太高的精确度（并且使用时间超过100年），并且地理位置是固定的，那么您可以整年进行分析并定期使用这些数据。 这样你就不需要运行时的第一步。

开普勒定律

如果你想这样看这里 。 你将需要地球的轨道和旋转参数。 这些是从太阳系的地雷星历引擎* .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，纹理，环和大气参数。 所以：

  1. 让Sun在笛卡尔坐标中定位(0,0,0)

  2. 根据开普勒定律计算地球位置(x,y,z)

    太阳然后(-x,-y,-z)以地心坐标表示

  3. 向后旋转，每日旋转，岁差，章动(-x,-y,-z) -> (x',y',z')

  4. 为您的地理位置计算NEH帧(North,East,High(Up))
  5. 将(x',y',z')成NEH本地坐标(xx,yy,zz)

  6. 计算： 

     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 (ea>=1.0) ea=0; c0=sqrt((1.0-ea)/(1.0+ea)); // some helper constants computation c1=sqrt((1.0+ea)/(1.0-ea)); c2=aa*(1-ea*ea); sO=sin(Oa); cO=cos(Oa); si=sin(-ia); ci=cos(-ia); b=aa*sqrt(1.0-ea); Ma-=oa; // correction Ma=Ma-pi2*floor(Ma/pi2); E=Ma; for (q=0;q<20;q++) E=M.a+ea*sin(E); // Kepler's equation V=2.0*atan(c1*tan(E/2.0)); r=c2/(1.0+ea*cos(V)); pos[0]=r*cos(V+oa-Oa); // now just compute heliocentric position along ecliptic ellipse pos[1]=r*sin(V+oa-Oa); // 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是input和输出matrix
  • ang是旋转angular度[rad]

现在你可以用rep来转换地球本地坐标系

  • LCS到GCS (l2g) ... (gx,gy,gz)=rep*(lx,ly,lz)
  • GCS到LCS (g2l) ... (lx,ly,lz)=inverse(rep)*(gx,gy,gz)

局部是地球坐标系和全球太阳坐标系

    Interesting Posts

    如何绘制天空图？

    太阳的位置，一天的时间，经度和纬度