Как расчитываются...-II (Луна)

[an1440]

Теперь хотелось бы поставить несколько вопросов по расчетам, связанным с Луной.

1. Как пишут в книгах, "фаза измеряется отношением площади освещенной части видимого диска к общей площади". То есть, для Луны она определяется как разностью эклиптических долгот луны и Солнца, так и эклиптической широтой Луны.
То, что в календаре (который висит в соседней комнате ) называется "первой четвертью" - это действительно фаза, или это то, что скорее можно назвать "возрастом" (т.е. без привлечения к расчету широты)?

2. Вычислить положение (Луны, Солнца и т.д.) на заданную дату можно однократным решением уравнений. А как решается обратная задача - вычислить дату/время для заданного положения? Компьютер, конечно, позволяет сделать это "варварским" способом - приближением от обратного, т.е. "наматыванием/разматыванием" времени до получения заданного положения. Но ведь это способ действительно варварский? Тем более, я подозреваю что такие вычисления астрономы делали еще и до появления компьютеров.
Как это делается на самом деле?

[Анатолий Волчков]

Фазы Луны рассчитываются по разности геоцентрических эклиптических долгот Луны и Солнца.

Учет широты не имеет смысла, так как в этом случае для разных мест на Земле были бы разные фазы. Широты Луны достигает 5 градусов, а разность топоцентрических широт Луны на полюсах Земли может достигать почти половины этой величины. Да и новолуние или полнолуние в таком случае были бы лишь в моменты затмений, а не ежемесячно.


Что касается вычисления моментов, когда Луна достигает определенного положения, то эта задача вообще может не иметь решения.

Во-первых, видимое положение Луны зависит от места наблюдателя на Земле, все же Луна довольно близко и параллакс ее большой.

Во-вторых, в заданном месте на небе Луна может никогда не появиться в пределах разумного периода времени. Кроме того, вблизи некоторых точек Луна может появляться неоднократно, и требуются критерии отбора момента из серии моментов.

Да и не ясно, какова должна быть структура файла для решения этой задачи. Ведь в последнее время теория движения Луны строится как численная модель с корректировкой реальными наблюдениями. Например, это файлы DE200/LE200. Файлы определяют барицентрические положения Луны по аргументу времени. А как строить файлы для решения обратной задачи? Ответ совершенно неясен. Неужели для ряда моментов времени строить трехмерную матрицу положений, в которой первая координата соответствует датам из некоторого интервала времени, а две других – положению на Земле. Каков будет объем такого массива? Думаю, что простой поиск в таком массиве займет времени не меньше, чем итерационный способ нахождения одного момента времени. Да и не нужен такой массив при использовании компьютеров.

Похоже, что на построение такой теории потребуется машинного времени больше, чем на все мыслимые вычисления подобных данных за время использования указанной теории. Мало кому требуются такие вычисления.

[an1440]

...в разных местах Земли разные фазы...
Так вот почему огурцы в Китае сажают раньше, чем у нас...
------------------------------------------------------------------
...движения Луны строится как численная модель с корректировкой реальными ...

  Ну - это уже "высшая" астрономия. Я, конечно, имел в виду совсем другое: как определяют например, момент полнолуния в ноябре 2017 года.
  Или как проверить, действительно ли  в пятницу 03 апреля  0033 года "в шестом же часу настала тьма по всей земле и продолжалась до часа девятого" (Марка 15,33).  (Хотя, прости Господи - такие вещи проверять нельзя. В них положено просто верить)

тут же поправляюсь - второй вопрос в ТАКОЙ формулировке - это как раз ПРЯМАЯ задача. Я хотел сказать: как "найти" например ближайшее затмение (лунное, в данном случае) после 1 янв 0033 (не прибегая к помощи очевидцев).

[Анатолий Волчков]

Найти момент ближайшего затмения можно либо с применением методики вычисления обстоятельств затмений, описанных в книгах, например у Мейеса в «Астрономические формулы для калькуляторов» (книга вроде есть в продаже в интернет-магазине http://edurss.ru/cgi-bin/db.pl?cp=&page=Catalog&id=〈=Ru&blang=ru&list=25&sort=name), либо применив программу-планетарий, показывающую затмения. Я, например, отыскиваю заранее неизвестные моменты затмений на десяток лет минут за двадцать.

[an1440]

Попробовал найти все полнолуния путем прокручивания времени с уточненным приближением вблизи каждых 180 град.
Да. Способ надо сказать варварский, хотя результат все же и получен.
Расчет на 100 лет занял более 120 секунд (пентиум-166), процедуры  MINI_SUN/MINI_MOON обрабатывались около 140 тысяч раз каждая.
Все-таки полнолуния считали как-то еще и до компъютеров...

Кто, кстати, автор упомянутых процедур? Если потребуется использовать - надо ведь будет сослаться?

Как далеко можно заглядывать  от базовых констант J2000,0 ? ( не с точностью до астронома, но и не до огурцесажальщика - в последнем случае, я полагаю, достаточно функции FRAC и средней продолжительности синодического месяца?)

[Анатолий Волчков]

Так зачем нужно ссылаться на автора программы, можно ведь просто сослаться на сайт.

Насчет того, как глубоко во времени можно использовать программу. Думаю, что на период, на котором имеются сведения о древних затмениях, можно использовать. Ведь среднее движение Луны и скорость замедления вращения Земли выводились с учетом сведений о древних затмениях. Сам я этим не интересовался, но статьи о вычислениях моментов старых затмений встречаются на сайтах, посвященных критике «Новой хронологии» Фоменко.

Этими программами я не пользовался, но случайно заметил, что в процедуре вычисления координат Луны нет данных о расстоянии до Луны. А это делает практически невозможным использование процедур для расчета моментов восходов-заходов, обстоятельств затмений, особенно солнечных. Можно, конечно, применить среднее расстояние до Луны, но точность такого приближения низкая. Почему автор поленился вставить вычисление расстояния – не ясно.

[an1440]

...невозможным использование процедур для расчета моментов восходов-заходов...  - то есть, параллакс нельзя учесть?

[Анатолий Волчков]

Могу предложить процедуру, вычисляющую координаты Луны с погрешностью в несколько секунд дуги. Процедура сделана на основе книги Мееса «Астрономические формулы для калькуляторов» (фактически это «урезанная» теория Брауна).

При использовании процедуры для вычисления координат на дату необходимо задать JdRez = Jd. Можно вычислять и на фиксированную эпоху, например при вычислении положений Луны относительно каталожных положений звезд, для чего задать JdRez = равноденствию каталога.

Точности процедуры вполне достаточно для расчета обстоятельств затмений, предвычисления моментов покрытий Луной объектов на небе и видимости Луны в разных пунктах Земли.

Ввиду того, что форум не позволяет публиковать длинные сообщения, разбил сообщений на части.


  Procedure MOON ( JD, JDREZ : DOUBLE; Var L, B, R, RA, DECL : DOUBLE );

{ MOON  HEOCENTRIC COORDINATES

   JD              - TIME MOMENT
   JDREZ      - RESULTING COORDINTES SYSTEM
   L, B           - ECLIPTICAL COORDINATES
        R                - DISTANCE IN AU ( 149.59787E+6 km )
        RA, DECL - EQUATORIAL COORDINATES
}

Var

   PAR, M, N,
   PI2, T2, GR,
   EPS, SEPS, CEPS, G, G1, G2, G3, G4, G5, DL, DZ, DR, F,
   SB, CB, TB, V, U,  X, Y, Z,
   T, TE, TCP, P, PP, DT, AA, BB, CP, C1, SL, CL,
   L1, M1, D, E, OM, OM1, OM2, BV, SA : DOUBLE ;

 Function Atan2 ( S, C : Double ) : Double  ;

 Const
   Pi    : Double = 3.14159265359 ;
   Pi2   : Double = Pi * 2.0 ;
   Pid2  : Double = Pi / 2.0 ;
   Pi270 : Double = Pid2 * 3.0 ;

 Var
   A  : Double ;

 Begin

          If ( Abs ( C ) > 0.0 )
           Then
             Begin
               If ( Abs ( C ) > Abs ( S ) )
                 Then
                   Begin
                     A := ArcTan ( S / C ) ;
                   End
                 Else
                   Begin
                     If ( ( C / S ) < 0.0 )
                       Then  A := Pi270 - ArcTan ( C / S )
                       Else  A := Pid2  - ArcTan ( C / S ) ;
                     End ;

               If ( C < 0.0 )  Then  A := A + Pi  ;
               If ( A < 0.0 )  Then  A := A + Pi2 ;
               If ( A > Pi2 )  Then  A := A - Pi2 ;

             End
           Else
             Begin
               If ( S > 0.0 )
                 Then  A := Pi / 2.0
                 Else  A := Pi270    ;
             End ;

         ATan2 := A ;

 End ;

[Анатолий Волчков]

 Begin

   CP  := 206264.806247;

   PI2 := 3.14159265359 * 2.0;
   GR  := PI2 / 360.0 ;

   T  := ( JD - 2415020.0  ) / 36525.0 ;
   T2 := T * T ;
   TB := T ;

       { фундаментальные аргументы }

   L1 := 270.4342  + 481267.8831  * T -0.0011  * T2 ;
   M  := 358.4758  +  35999.0498  * T -0.0002  * T2 ;
   M1 := 296.1046  + 477198.8491  * T +0.0092  * T2 ;
   D  := 350.7375  + 445267.1142  * T -0.0014  * T2 ;
   F  :=  11.2509  + 483202.0251  * T -0.0032  * T2 ;

   L1 := L1 - Trunc ( L1 / 360.0 ) * 360.0;
   M  := M  - Trunc ( M  / 360.0 ) * 360.0;
   M1 := M1 - Trunc ( M1 / 360.0 ) * 360.0;
   D  := D  - Trunc ( D  / 360.0 ) * 360.0;
   F  := F  - Trunc ( F  / 360.0 ) * 360.0;

   OM := 259.1833  -   1934.1420  * T +0.0021  * T2 ;
   SA := ABS ( OM );
   If ( SA > 360.0 ) Then
     Begin
       BV := OM / SA ;
       OM := ( SA - Trunc ( SA / 360.0 ) * 360.0 ) * BV;
     End;

      {дополнительные члены}

   SA := SIN ( (51.2 + 20.2 * T ) * GR );

   M  := M  -0.0018 * SA;
   M1 := M1 +0.0008 * SA;
   D  := D  +0.0020 * SA;

   BV := 0.0040 * SIN ( ( 346.560 + 132.870 * T ) * GR );

   L1 := L1 + BV;
   M1 := M1 + BV;
   D  := D  + BV;
   F  := F  + BV;   

   SA := SIN ( OM * GR );

   L1 := L1 + 0.0020 * SA;
   M1 := M1 + 0.0025 * SA;
   D  := D  + 0.0020 * SA;
   F  := F  - 0.0247 * SA;

   F  := F  - 0.0043 * SIN ( ( OM + 275.05 - 2.3 * T ) * GR );

   E  := 1.0 - 0.002495 * T - 0.0000075 * T * T;


   L := L1;

   L1 := L1 * GR;
   M  := M  * GR;
   M1 := M1 * GR;
   D  := D  * GR;
   F  := F  * GR;
   OM := OM * GR;


      {долгота, широта и параллакс Луны (видимые на дату}

   L := L +  6.2888 * SIN ( M1 );
   L := L +  1.2740 * SIN ( 2.0 * D - M1 );
   L := L +  0.6583 * SIN ( 2.0 * D );
   L := L +  0.2136 * SIN ( 2.0 * M1 );
   L := L -  0.1856 * SIN ( M ) * E;
   L := L -  0.1143 * SIN ( 2.0 * F );
   L := L +  0.0588 * SIN ( 2.0 * ( D - M1 ) );
   L := L +  0.0572 * SIN ( 2.0 * D - M - M1 ) * E;
   L := L +  0.0533 * SIN ( 2.0 * D + M1 );
   L := L +  0.0459 * SIN ( 2.0 * D - M ) * E;
   L := L +  0.0410 * SIN ( M1 - M ) * E;
   L := L -  0.0347 * SIN ( D );
   L := L -  0.0305 * SIN ( M + M1 ) * E;
   L := L +  0.0153 * SIN ( 2.0 * ( D - F ) );
   L := L -  0.0125 * SIN ( 2.0 * F + M1 );
   L := L -  0.0110 * SIN ( 2.0 * F - M1 );
   L := L +  0.0107 * SIN ( 4.0 * D - M1 );
   L := L +  0.0100 * SIN ( 3.0 * M1 );
   L := L +  0.0085 * SIN ( 4.0 * D - 2.0 * M1 );
   L := L -  0.0079 * SIN ( M - M1 + 2.0 * D );
   L := L -  0.0068 * SIN ( 2.0 * D + M );
   L := L +  0.0052 * SIN ( M1 - D );
   L := L +  0.0050 * SIN ( M + D );
   L := L +  0.0040 * SIN ( M1 - M + 2.0 * D );
   L := L +  0.0040 * SIN ( 2.0 * ( M1 + D ) );
   L := L +  0.0039 * SIN ( 4.0 * D );
   L := L +  0.0037 * SIN ( 2.0 * D - 3.0 * M1 );
   L := L +  0.0027 * SIN ( 2.0 * M1 - M );
   L := L +  0.0026 * SIN ( M1 - 2.0 * F - 2.0 * D );
   L := L +  0.0024 * SIN ( 2.0 * D - M - 2.0 * M1 );
   L := L -  0.0023 * SIN ( M1 + D );
   L := L +  0.0022 * SIN ( 2.0 * D - 2.0 * M );
   L := L -  0.0021 * SIN ( 2.0 * M1 + M );
   L := L -  0.0021 * SIN ( 2.0 * M );
   L := L +  0.0021 * SIN ( 2.0 * D - M1 - 2.0 * M );
   L := L -  0.0018 * SIN ( M1 + 2.0 * D - 2.0 * F );
   L := L -  0.0016 * SIN ( 2.0 * F + 2.0 * D );
   L := L +  0.0012 * SIN ( 4.0 * D - M - M1 );
   L := L -  0.0011 * SIN ( 2.0 * M1 + 2.0 * F );
   L := L +  0.0009 * SIN ( M1 - 3.0 * D );
   L := L -  0.0008 * SIN ( M + M1 + 2.0 * D );
   L := L +  0.0008 * SIN ( 4.0 * D - M - 2.0 * M1 );
   L := L +  0.0007 * SIN ( M1 - 2.0 * M );
   L := L +  0.0007 * SIN ( M1 - 2.0 * M - 2.0 * D );
   L := L +  0.0007 * SIN ( M  - 2.0 * M1 + 2.0 * D );
   L := L +  0.0006 * SIN ( 2.0 * D - M - 2.0 * F );
   L := L +  0.0006 * SIN ( M1 + 4.0 * D );
   L := L +  0.0005 * SIN ( 4.0 * M1 );
   L := L +  0.0005 * SIN ( 4.0 * D - M );
   L := L +  0.0005 * SIN ( 2.0 * M1 - D );


   B :=  0.0;

   B := B +  5.128 * SIN ( F );
   B := B +  0.281 * SIN ( M1 + F );
   B := B +  0.278 * SIN ( M1 - F );
   B := B +  0.173 * SIN ( 2.0 * D - F );
   B := B +  0.055 * SIN ( 2.0 * D + F - M1 );
   B := B +  0.046 * SIN ( 2.0 * D - F - M1 );
   B := B +  0.033 * SIN ( 2.0 * D + F );
   B := B +  0.017 * SIN ( 2.0 * M1 + F );
   B := B +  0.009 * SIN ( 2.0 * D + M1 - F );
   B := B +  0.009 * SIN ( 2.0 * M1 - F );
   B := B +  0.008 * SIN ( 2.0 * D - M - F );
   B := B +  0.004 * SIN ( 2.0 * D - F - 2.0 * M1 );
   B := B +  0.004 * SIN ( 2.0 * D + F + M1 );
   B := B +  0.003 * SIN ( F - M - 2.0 * D );
   B := B +  0.002 * SIN ( 2.0 * D + F - M - M1 );
   B := B +  0.002 * SIN ( 2.0 * D + F - M );
   B := B +  0.002 * SIN ( 2.0 * D - F - M - M1 );
   B := B +  0.002 * SIN ( F - M + M1 );
   B := B +  0.002 * SIN ( 4.0 * D - F - M1 );
   B := B -  0.002 * SIN ( F + M );
   B := B -  0.002 * SIN ( 3.0 * F );
   B := B +  0.002 * SIN ( M1 - M - F );
   B := B -  0.001 * SIN ( F + D );
   B := B -  0.001 * SIN ( F + M + M1 );
   B := B +  0.001 * SIN ( F - M - M1 );
   B := B +  0.001 * SIN ( F - M );
   B := B +  0.001 * SIN ( F - D );
   B := B +  0.001 * SIN ( F + 3.0 * M1 );
   B := B +  0.001 * SIN ( 4.0 * D - F );
   B := B +  0.001 * SIN ( F + 4.0 * D - M1 );
   B := B +  0.001 * SIN ( M1 - 3.0 * F );
   B := B +  0.001 * SIN ( F + 4.0 * D - 2.0 * M1 );
   B := B +  0.001 * SIN ( 2.0 * D - 3.0 * F );
   B := B +  0.001 * SIN ( 2.0 * D + 2.0 * M1 - F );

   OM1 := 0.0004664 * COS ( OM );
   OM2 := 0.0000754 * COS ( OM + 275.05 - 2.30 * T );

   B := B * ( 1.0 - OM1 - OM2 );

   PAR := 0.950724;

   PAR := PAR + 0.051818 * COS ( M1 );
   PAR := PAR + 0.009531 * COS ( 2.0 * D - M1 );
   PAR := PAR + 0.007843 * COS ( 2.0 * D );
   PAR := PAR + 0.002824 * COS ( 2.0 * M1 );
   PAR := PAR + 0.000857 * COS ( 2.0 * D + M1 );
   PAR := PAR + 0.000533 * COS ( 2.0 * D - M );
   PAR := PAR + 0.000401 * COS ( 2.0 * D - M - M1 );
   PAR := PAR + 0.000320 * COS ( M1 - M );
   PAR := PAR - 0.000271 * COS ( D );
   PAR := PAR - 0.000264 * COS ( M + M1 );
   PAR := PAR - 0.000198 * COS ( 2.0 * F - M1 );
   PAR := PAR + 0.000173 * COS ( 3.0 * M1 );
   PAR := PAR + 0.000167 * COS ( 4.0 * D - M1 );
   PAR := PAR - 0.000111 * COS ( M );
   PAR := PAR + 0.000103 * COS ( 4.0 * D - 2.0 * M1 );
   PAR := PAR - 0.000084 * COS ( 2.0 * M1 - 2.0 * D );
   PAR := PAR - 0.000083 * COS ( 2.0 * D + M );
   PAR := PAR + 0.000079 * COS ( 2.0 * D + 2.0 * M1 );
   PAR := PAR + 0.000072 * COS ( 4.0 * D );
   PAR := PAR + 0.000064 * COS ( 2.0 * D - M + M1 );
   PAR := PAR - 0.000063 * COS ( 2.0 * D + M - M1 );
   PAR := PAR + 0.000041 * COS ( M + D );
   PAR := PAR + 0.000035 * COS ( 2.0 * M1 - M );
   PAR := PAR - 0.000033 * COS ( 3.0 * M1 - 2.0 * D );
   PAR := PAR - 0.000030 * COS ( M1 + D );
   PAR := PAR - 0.000029 * COS ( 2.0 * F - 2.0 * D );
   PAR := PAR - 0.000029 * COS ( 2.0 * M1 + M );
   PAR := PAR + 0.000026 * COS ( 2.0 * D - 2.0 * M );
   PAR := PAR - 0.000023 * COS ( 2.0 * F - 2.0 * D + M1 );
   PAR := PAR + 0.000019 * COS ( 4.0 * D - M - M1 );

   L   := L   * GR;
   B   := B   * GR;
   PAR := PAR * GR;

   R := 6378.14 / SIN ( PAR ) / 149.59787E+6 ;

   SB := SIN ( B );
   CB := COS ( B );

      {перевод TB -> TE}

   TE  := ( JDREZ - 2415020.0 ) / 36525.0;

   TCP := ( TB + TE ) / 2.0;
   P   := ( 5025.6 + 2.2 * TCP ) / CP;
   PP  :=   3.036015 +  0.0159358 * TCP;

   DT := TB - TE;
   AA := P * DT;
   BB := 47.1 / CP * DT;
   C1 := L + PI2 / 2.0 -  0.5 * AA - PP;

   L := L - AA + BB * COS ( C1 ) * SB/CB;
   B := B -      BB * SIN ( C1 );

   SB := SIN ( B );
   CB := COS ( B );
   TB := SB / CB  ;
   L  := L - Trunc ( L / PI2 ) * PI2 ;

   SL := SIN ( L );
   CL := COS ( L );

       T    := ( JDREZ - 2415020.0 ) / 36525.0;
       EPS  := ( 23.0 * 3600.0 + 27.0 * 60.0 + 8.260 - 46.845 * T  ) / CP ;
       SEPS := SIN ( EPS ) ;
       CEPS := COS ( EPS ) ;

   X := R * CB *   CL;
   Y := R * CB * ( SL * CEPS - TB * SEPS );
   Z := R * CB * ( SL * SEPS + TB * CEPS );

       RA   := ATan2 ( Y, X ) ;
       Decl := Atan2 ( Z, Sqrt ( Sqr ( X ) + Sqr ( Y ) )  ) ;
       If ( Decl > Pi )  Then
          Decl := Decl - Pi2 ;

End;

[an1440]

Все!!!!!
Ошибку нашел. Оказалась она не в программе, а на стадии задания исходных данных - вместо JDREZ=JD затесалось IDREZ=ID...

Радианы-то само собой; а "бред" - в смысле что почти до самого конца все идет нормально, а в самом конце - все сбивается. Как теперь оказалось, спотыкаясь об злополучный idrez...

[an1440]

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

[Sleepwalker]

Вопрос к знатокам:
Вот в старкальке для Луны, да и для всех других планет, экваториальные координаты выдаются как геоцентрические так и топоцентрические. Насколько я понял для перевода в горизонтальные координаты нужны именно экваториальные топоцентрические, так ?
Геоцентрические я считать научился, а дальше туман...может кто подскажет?

[ZUDA32.EXE]

Гео- и топоцентрические координаты отличаются на величину параллакса. Считается, что для технических расчетов эту разницу достаточно учитывать только для Луны (ну и, разумеется, для того, что запущено на небо человеком), т.к. для остальных светил она слишком мала.