вычислить время, когда солнце находится на X градусов ниже / выше горизонта

1. собирать данные эфемерид Солнца за день, в котором вы нуждаетесь

   сделайте шаг в 1 час и получите положение Солнца в азимутальных координатах для нужного вам географического местоположения. Либо используйте найденные вами уравнения, либо используйте некоторые службы WEB, например:

   • JPL Horizons не нравится этот, поскольку они используют некоторые странные исходные кадры отсчета, которые не соответствуют моим измерениям, но, скорее всего, я превращаю что-то не так на этом пути …

   • Пресовская обсерватория – это моя любимая (но она на словацком языке), ее продукция легко копируется в шахтные двигатели, а выход соответствует шахтным наблюдениям, вычислениям, оценкам и измерениям. Просто заполните:

     • гео-местоположение (Miesto pozorovania)
     • дата, время (Dátum a čas pozorovania)
     • в нижней части слева: интервал [дней], интервал шаг [дней]
     • нажмите на кнопку для Sun (Slnko), Moon (Mesiac), Planets (Planety)

   есть много таких страниц, которые просто смотрят, но всегда проверяют, выдают ли они правильные данные. Я использую законы Келлера / уравнения, формирующие планетарные движения (более низкая точность, но для Земли-Солнца должно быть ОК). В настоящее время в двигателях используется гравитационная модель (но это неустойчиво с более длительным периодом времени)

2. обрабатывать данные как набор трехмерных точек вдоль полилинии (азимут, высота, время)

3. теперь просто найдите в данных 2 балла

   один ниже желаемого угла и следующий выше желаемый угол. Точки будки должны быть соседними. Если какая-либо точка находится на нужном угле, то у вас уже есть решение, поэтому остановитесь

4. интерполировать время пересечения угла высоты

   

   поэтому, если желаемый угол роста равен 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, текстуры, кольца и атмосферу. Так:

1. получить Солнце в положение (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 – матрица ввода и вывода
• ang – угол поворота [рад]

теперь вы можете использовать rep для преобразования в / из локальной системы координат Земли

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

локальная система координат Земли и глобальная система координат Солнца