计算太阳在地平线以下/以上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那么:

      • a0a1是方位angular
      • b0b1是高度angular
      • t0t1是次数

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

       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)

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