Рассчитать долготу дня

[RN3QGA]

  Все понял, Сказал Хорек - какие суслики!

[Tau]

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

[RN3QGA]

  Ох, Tau лучше в мою шкуру не надо, ну её.
Сложности просто в том, что я не в том возрасте видимо увлекся этим надо сказать интересным занятием., поэтому пробелы сплошные, но это ничего я ведь тока учусь.
Спасибо, что не оставляете и не отфутболиваете. Если будет необходимость я обязательно Вам напишу и попрошу о помощи, а просто отнимать драгоценное время которого всегда нехватает, или  по пустякам забивать голову я просто не имею права.
  С уважением RN3QGA.

[RN3QGA]

  Для луны, что-то в 8 пункте застрял. Правда и еще есть вопросы но по порядку.

1. H (высота луны)  в пункте 7. на  7 октября 2010 г. 0ч:00м у меня выходит =56,1961751298217
    Я думаю где то в расчете косяк поэтому с вашего позволения приведу кусок кода в прицепе и ЕХЕ для наглядности в
    нем есть и другие вычисленные значения.
2. Ну и пункт 8. честно говоря не совсем понял, а точнее в коде я отметил это комментарием.
    Поправьте пожалуйста, а если нужно подробней-  насколько смогу опишу.
С уважением RN3QGA.
 

[Tau]

1) Time:= Hour + Minute/100;  делить на 60 (патамушта в часе 60 минут), а не на 100. Для полуночи неважно, в ином случае вылезет ерунда.

2) Аргумент тригонометрических функций (sin,cos,tg) всегда должен быть в радианах.
Здесь у вас S в градусах:
Xg:=x1*cos(S)+Y1*sin(S);
Yg:=(-X1)*sin(S)+Y1*cos(S);

Аналогично здесь широта и звездное время:
Xn:=(ER*Czv+Hvus)*cos(Sl)*cos(fshir);
Yn:=(ER*Czv+Hvus)*sin(Sl)*cos(fshir);
Zn:= (ER*Wzv+Hvus)*sin(fshir);

Аналогично и здесь все аргументы должны быть в радианах:
LunHH:=Pil*ArcSin(sin(fshir)*sin(Btpc)+cos(fshir)*cos(Btpc)*cos(sl-Atpc)+0.85)

3) Czv1:= sqr(1+0.00335*(0.00335-2)*sin(2*fshir))
Здесь должен быть квадрат синуса, а не синус двойного аргумента.

4) Czv1:= sqr(1+0.00335*(0.00335-2)*sin(2*fshir))
sqr - это возведение в квадрат, а нужно извлечь корень, то есть использовать паскалевскую функцию sqrt

5) В длинных выражениях при переносе строки не ставьте минусы и в конце строки и в начале следующей,  языки программирования тупые - они два минуса считают плюсом по закону арифметики.
Правильно таК:
w = 5.13 * Sin(q * (93.3 + 483202.03 * T)) +
       0.28 * Sin(q * (228.2 + 960400.87 * T)) -
       0.28 * Sin(q * (318.3 + 6003.18 * T)) -
       0.17 * Sin(q * (217.6 - 407332.2 * T))
L = 218.32 + 481267.883 * T + 6.29 * Sin(q * (134.9 + 477198.85 * T)) -
      1.27 * Sin(q * (259.2 - 413335.38 * T)) + 0.66 * Sin(q * (235.7 + 890534.23 * T)) +
      0.21 * Sin(q * (269.9 + 954397.7 * T)) - 0.19 * Sin(q * (357.5 + 35999.05 * T)) -
      0.11 * Sin(q * (186.6 + 966404.05 * T))

6) Арктангенс от двух аргументов вы уже встречали не раз . Когда долготу и широту Солнца и Луны вычисляли, то вы раскладывали вектор на сферические углы. Это универсальная процедура. Любой вектор раскладывается на углы также.
Соответственно, топоцентрические координаты:
Atpc:=ArcTan2(Yt,Xt);
Btpc:=ArcTan2(Zt,sqrt(Xt*Xt+Yt*Yt));
и перевести их из радиан в градусы.

7) К вопросу о пункте 8.
Запустили цикл, вычисляете сначала Tk и MD1, а потом вызываете процедуру TheLunCalc, чтобы она вычислила высоту Луны на момент времени MD1. Проблема в том, что входными параметрами в TheLunCalc вы сделали Year, Month, Day, Hour, Minute. В цикле удобнее работать с равномерной шкалой времени, например с юлианской датой. В этой задаче вам юлианскую дату нужно вычислить один раз - ту дату, от которой начинаются вычисления,  вот и вычислите ее в главном теле программы до запуска цикла. Потом в цикле вы просто наращиваете  юлианскую дату путем прибавления к ней шага по времени (MD1:=MD+(Tk/24)). Эту юлианскую дату и отправляйте в TheLunCalc в качестве MD, соответственно в самой процедуре вычисления начнутся с  координат Луны, а не с вычисления юлианской даты.

[RN3QGA]

  Да, странно даже самому. Особенно не понятно почему я в часе определил 100 минут,  во лапух .  Вроде же все очевидно, ну кроме 8. пункта, там, я и запутался с этими входными данными, ладненько как переделаю отпишусь.   
С уважением RN3QGA.

[RN3QGA]

  Так, решил пока без 8 пункта разобраться с вычислениями.
Все вроде переделал в блоке вычислений, выкинул пока вычисление Юл. даты поставил фиксированное
MD=55476.0; и время Time:=0;  у меня высота равна 42,4900629969896. у Вас на 0.00, H=-25.800079
В прицепке исправленный Мной вариант. Видимо, где то еще нужно исправить, но я уже запутался честно говоря.

[Tau]

1) Atpc:=ArcTan2(Yt,Xt);
  Atpc:=DegtoRad(Atpc);
  Btpc:=ArcTan2(Zt,sqrt((Xt*Xt)+(Yt*Yt)));
  Btpc:=DegtoRad(Btpc);

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

2)  LunHH:=Pil*ArcSin(sin(fshir)*sin(Btpc)+cos(fshir)*cos(Btpc)*cos(sl-Atpc)+0.85);
Здесь аргумент арксинуса закрывается раньше, 0.85 вне арксинуса, скобочку переставьте. Даже странно, что у вас в этом месте программа не вылетает из-за того, что аргумент арксинуса больше 1.

Сравниваем:
Atpc=3.153429 радиан
Btpc=-0.111419 радиан
Если у вас не так,  гоните экзешник с выводом на печать цифири.

[RN3QGA]

  Вот он, наконец то, -  это Ваш "ПОСТ"       (Настоящие программисты вообще предпочитают работать только с радианами)  это же точно про меня, ну я думаю Вы все поняли , ну , а то, что (2) программа не Вылетает, это Вы не поверите: Вы, даже не 20-100 человек который удивляется, что у меня некоторые вещи работают  против всех правил.  .
Поправлю отпишусь.

[RN3QGA]

  Вот то все, что наваял см. прицеп.  

[Tau]

В переписанном на VB виде ваша процедура работает близко к правде, а в экзешнике на стадии вычисления звездного времени у вас происходит непонятно что. Покажите функции DegToRad и DegToGrad.

[RN3QGA]

  С личкой понял. не обращаем внимания.
Функции DegtoRad и DegtoGrad это встроенные функции Delphi  перевода из радиан в градусы и наоборот. При вычислении терминатора и положения на карте солнца и луны я их использовал чтоб лишка не писать, все работает ОК. Я наверное чуть попозже (сегодня арбайтн) пошагово каждое вычисление сделаю и прицеплю с комментами, и самому будет полезно. 
С уважением RN3QGA.

[Tau]

Ой  , в справочнике по функциям Дельфи DegToGrad - функция перевода градусов в грады. Вместо умножения угла на 180/pi эта функция делит угол на 0.9.
А радианы в градусы переводит RadToDeg

[RN3QGA]

Оп, это моя опечаточка виноват исправлюсь.

[RN3QGA]

Вот переделал, но точно как у Вас все равно не выходит. Вот в ЕХЕ вариант на 16.00 
хотя для 0 часов вроде все совпало, тока при этом Atpc=-3.153429 радиан (со знаком -),
Btpc=-0.111419 радиан. это совпало точно.

Код: [Выделить]

procedure TheLunCalc;
 const
  A1 = 24110.54841;
  a2 = 8640184.812 ;
  A3 = 0.093104 ;
  a4 = 0.0000062 ;
var
  MD,T,W,L,X,Y,Z,e,e2,X1,Y1,Z1,So,ww,AT,Ss,S,Xg,Yg,Zg,Ao,Fo,F,AD:Double;   
  A,Xl,Yl,X0,Y0:integer;
  q1,q2,q3,ao1,ao2,ao3,ao4,q,p,ER,r:Double;
  i,ii: Integer;

  Fshir,Ldolg,Sl:Double; 
  Hvus:integer;           
  Xn,Yn,Zn:double;     
  Fzv,Wzv,Czv,Czv1,Wxz:Double;

    Xt,Yt,Zt:Double;
 
    Atpc,Btpc,ftpc,htpc,LunHH,PiL:Double;

 begin

  q:= (pi / 180);
 
  MD:=55476.0;
  Time:=16;

  T := (MD+Time/24 - 51544.5) / 36525;

  L := 218.32 + 481267.883 * T +6.29 * Sin(q * (134.9 + 477198.85 * T))-
      1.27 * Sin(q * (259.2 - 413335.38 * T))+ 0.66 * Sin(q * (235.7 + 890534.23 * T))+
      0.21 * Sin(q * (269.9 + 954397.7 * T))-0.19 * Sin(q * (357.5 + 35999.05 * T))-
      0.11 * Sin(q * (186.6 + 966404.05 * T));


   W := 5.13 * Sin(q * (93.3 + 483202.03 * T)) +
       0.28 * Sin(q * (228.2 + 960400.87 * T)) -
       0.28 * Sin(q * (318.3 + 6003.18 * T)) -
       0.17 * Sin(q * (217.6 - 407332.2 * T));

    p := 0.9508 + 0.0518 * Cos(q * (134.9 + 477198.85 * T)) +
       0.0095 * Cos(q * (259.2 - 413335.38 * T)) +
       0.0078 * Cos(q * (235.7 + 890534.23 * T)) +
       0.0028 * Cos(q * (269.9 + 954397.7 * T));

    r := 1 / Sin(p * q);

    ER := 6378.136;

    X:= r * ER * Cos(q*L) * Cos(q*W);
    Y := r * ER * Sin(q*L) * Cos(q*W);
    Z := r * ER * Sin(q*W);

 
  e2 := 84381.488 - 46.815 * T - 0.00059 * Sqr(T) + 0.001813 * T * Sqr(T);
  e  := e2/3600;
  e  := (e*pi)/180 ;// в радианы

  X1 := X;
  Y1 := Y * cos(e) - Z * sin(e);
  Z1 := Y * sin(e) + Z * cos(e);


  So := A1+a2*T+a3*sqr(T)-a4*T*sqr(T);
  ww := 1.0027379093508;
  AT := (Time/24)*86400*ww;
  Ss := So+AT;
  S  := Ss/240;

  S  := (S*pi)/180;// звездное время в радианы
 
  Xg:=x1*cos(S)+Y1*sin(S);
  Yg:=(-X1)*sin(S)+Y1*cos(S);
  Zg:=Z1;

  Ao:=ArcTan2(yg,xg);

  Fo:=ArcTan2(Zg,Hypot(Xg, Yg));

  LunX := Ao+180;
  LunY := Fo;

  fshir:=60;     
  Ldolg:=30;     
  Hvus:=0;       
  fzv:=0.00335; 
  Wxz:=1-0.00335;

  Fshir:=(Fshir*pi)/180; //Широту в радианы
  Ldolg:=(Ldolg*pi)/180; //Долготу в радианы


  Sl:=S+Ldolg;         

  Form1.ListBox1.Items[2]:=FloatToStr(Sl);

  Czv1:= sqrt(1+0.00335*(0.00335-2)*(sqr(sin(fshir)))); // Вот здесь может я чегой то напортачил ?

  Czv:=(1/czv1);
  Wzv:=Czv*(Wxz*Wxz);

  Xn:=(ER*Czv+Hvus)*cos(Sl)*cos(fshir);
  Yn:=(ER*Czv+Hvus)*sin(Sl)*cos(fshir);
  Zn:=(ER*Wzv+Hvus)*sin(fshir);
 
  Xt:=X1-Xn;
  Yt:=Y1-Yn;
  Zt:=Z1-Zn;
 
  Atpc:=ArcTan2(Yt,Xt);
  Btpc:=ArcTan2(Zt,sqrt((Xt*Xt)+(Yt*Yt)));


  Pil:=180/pi;

  LunHH:=Pil*ArcSin(sin(fshir)*sin(Btpc)+cos(fshir)*cos(Btpc)*cos(sl-Atpc))+0.85;

  end;
Лишние комменты и вывод на печать здесь убрал. В прицепе ЕХЕ по этому коду с выводом на печать.

[Tau]

При вычислении координат Луны переменная времени T должна быть выражена в юлианских датах на момент расчета:
T := (MD+Time/24 - 51544.5) / 36525;
L := 218.32 + 481267.883  + 481267.883 +6.29 * Sin(q * (134.9 + 477198.85 * T))- ....
А при вычислении звездного гринвического времени на полночь даты переменная T должна быть выражена в юлианских датах на момент полуночи.
T := (MD - 51544.5) / 36525;
So := A1+a2*T+a3*sqr(T)-a4*T*sqr(T);
а уж потом прибавляется огрызок времени AT, прошедший от полуночи..

Лучше все-таки передавать в TheLunCalc юлианскую дату момента расчета MD = MD(полуночи)+Time/24 и вычислять аргументы T так, как в моем опусе в п.1 и п.2 с использованием операции округления INT. Но вам виднее.

Atpc=-3.153429 радиан

А не atpc=-3.129757?

[RN3QGA]

  Да уж кому кому, а мне точно не видней.
atpc=-3.129757 - точно , у меня по ходу действительно конфуз 
Щас попробую разобраться. отпишусь минут через 10-20, так во чего навоял
начало совпадает дальше хоть убейте : в прицепке посмотрите

Код: [Выделить]

procedure TheLunCalc; //Âû÷èñëåíèå âîñõîäà è çàõîäà ëóíû
 const
  A1 = 24110.54841;
  a2 = 8640184.812 ;
  A3 = 0.093104 ;
  a4 = 0.0000062 ;
var
  T,W,L,X,Y,Z,e,e2,X1,Y1,Z1,So,ww,AT,Ss,S,Xg,Yg,Zg,Ao,Fo,F,AD:Double;   //Îñíîâíûå ïåðåìåííûå
  A,Xl,Yl,X0,Y0:integer;
  ao1,ao2,ao3,ao4,q,p,ER,r,ah:Double;
  i,ii: Integer; //Âñïîìîãàòåëüíûé ïåðåìåííûå

   //ïåðåìåííûå äëÿ ãåîöåíòð. âåêòîðà íàáëþäàòåëÿ

  Fshir,Ldolg,Sl:Double;  //øèðîòà, äîëãîòà,ìåñòíîå çâåçäíîå âðåìÿ äëÿ ìåñòà
  Hvus:integer;           //Âûñîòà íà óðîâíåì ãåîèäà
  Xn,Yn,Zn:double;       //êîîðäèíàòû íà âûõîäå
  Fzv,Wzv,Czv,Czv1,Wxz:Double;//Âñïîìîãàòåëüíûå äëÿ ðàñ÷åòà

   // ïåðåìåííûå äëÿ òîïîöåíòðè÷åñêîãî âåêòîðà ëóíû

    Xt,Yt,Zt:Double;

  // ïåðåìåííûå äëÿ òîïîöåíòðè÷åñêîãî âåêòîðà ëóíû è ðàññ÷åòà âûñîòû
    Atpc,Btpc,ftpc,htpc,LunHH,PiL:Double;

 begin

  q:= (pi / 180);

   //Âû÷èñëÿåì åäèíè÷íûé ãåîöåíòðè÷åñêèé âåêòîð Ëóíû íà äàòó MD è âðåìÿ UT  â ýêëèïòè÷åñêîé ñèñòåìå êîîðäèíàò


  T := (MD - 51544.5) / 36525;

  L := 218.32 + 481267.883 * T +6.29 * Sin(q * (134.9 + 477198.85 * T))-
      1.27 * Sin(q * (259.2 - 413335.38 * T))+ 0.66 * Sin(q * (235.7 + 890534.23 * T))+
      0.21 * Sin(q * (269.9 + 954397.7 * T))-0.19 * Sin(q * (357.5 + 35999.05 * T))-
      0.11 * Sin(q * (186.6 + 966404.05 * T));


   W := 5.13 * Sin(q * (93.3 + 483202.03 * T)) +
       0.28 * Sin(q * (228.2 + 960400.87 * T)) -
       0.28 * Sin(q * (318.3 + 6003.18 * T)) -
       0.17 * Sin(q * (217.6 - 407332.2 * T));

    p := 0.9508 + 0.0518 * Cos(q * (134.9 + 477198.85 * T)) +
       0.0095 * Cos(q * (259.2 - 413335.38 * T)) +
       0.0078 * Cos(q * (235.7 + 890534.23 * T)) +
       0.0028 * Cos(q * (269.9 + 954397.7 * T));

    r := 1 / Sin(p * q);

    ER := 6378.136;

    X:= r * ER * Cos(q*L) * Cos(q*W);
    Y := r * ER * Sin(q*L) * Cos(q*W);
    Z := r * ER * Sin(q*W);


  //Âû÷åñëåíèå âåêòîðà Ëóíû â ýêâàòîðèàëüíîé ñèñòåìå êîîðäèíàò
  e2 := 84381.488 - 46.815 * T - 0.00059 * Sqr(T) + 0.001813 * T * Sqr(T);
  e  := e2/3600;
  e  := (e*pi)/180 ;//èç ãðàä. â ðàäèàíû

  X1 := X;
  Y1 := Y * cos(e) - Z * sin(e);
  Z1 := Y * sin(e) + Z * cos(e);

  //Âû÷èñëåíèå ñðåäíåãî çâåçäíîãî âðåìåíè íà äàòó MD è âðåìÿ UT.

  T := (Int(MD) - 51544.5) / 36525;
  So := A1+a2*T+a3*sqr(T)-a4*T*sqr(T);
  ww := 1.002737909350795;
  AT := (MD-int(MD))*86400*ww;
  Ss := So+AT;
  S  := Ss/240;
  S  := (S*pi)/180;//âðåìÿ èç ãðàäóñîâ â ðàäèàíû

  //Ïåðåâîä âåêòîðà Ëóíû â ãðèíâè÷ ñèñòåìó êîîðäèíàò
  Xg:=x1*cos(S)+Y1*sin(S);
  Yg:=(-X1)*sin(S)+Y1*cos(S);
  Zg:=Z1;

  Ao:=ArcTan2(yg,xg);

  Fo:=ArcTan2(Zg,Hypot(Xg, Yg));

  //êîîðäèíàòû ëóíû íà êàðòå

  LunX := Ao+180;
  LunY := Fo;

  // Áëîê âûâîäà íà ýêðàí ïðîìåæóòî÷íûõ äàííûõ
  Form1.ListBox1.Items[0]:=FloatToStr(Lunx)+'      //ïîëîæåíèå ëóíû íà êàðòå îñü Õ';
  Form1.ListBox1.Items[1]:=FloatToStr(LunY)+'      //ïîëîæåíèå ëóíû íà êàðòå îñü Y';

  // Âû÷èñëåíèå ãåîöåíòðè÷åñêîãî âåêòîðà íàáëþäàòåëÿ íà MD

  fshir:=60;      //ïåðåìåííàÿ øèðîòû
  Ldolg:=30;      //ïåðåìåííàÿ äîëãîòû
  Hvus:=0;        // ïåðåìåííàÿ âûñîòû ãåîèäà
  fzv:=0.00345709;   // ñæàòèå çåìëè óñðåäíåííàÿ îò (1/298,26)
  Wxz:=1-0.00345709;

  Fshir:=(Fshir*pi)/180; //øèðîòà â ðàäèàíû
  Ldolg:=(Ldolg*pi)/180; //äîëãîòà â ðàäèàíû


  Sl:=S+Ldolg;         // ìåñòíîå çâåçäíîå âðåìÿ



  Form1.ListBox1.Items[2]:=FloatToStr(Sl)+'     // ìåñòíîå çâåçäíîå âðåìÿ';

 

  Czv1:= sqrt(1+0.00345709*(0.00335-2)*(sqr(sin(fshir)))); //??? Óòî÷íèòü

  Czv:=(1/czv1);
  Wzv:=Czv*(Wxz*Wxz);


  Xn:=(ER*Czv+Hvus)*cos(Sl)*cos(fshir);
  Yn:=(ER*Czv+Hvus)*sin(Sl)*cos(fshir);
  Zn:=(ER*Wzv+Hvus)*sin(fshir);

  Form1.ListBox1.Items[3]:=FloatToStr(Xn)+'     //âåêòîð íàáëþäàòåëÿ Õ íà MD';
  Form1.ListBox1.Items[4]:=FloatToStr(Yn)+'     //âåêòîð íàáëþäàòåëÿ Y íà MD';
  Form1.ListBox1.Items[5]:=FloatToStr(Zn)+'     //âåêòîð íàáëþäàòåëÿ Z íà MD';


  // ÂÛ×ÈÑËÅÍÈÅ ÒÎÏÎÖÅÍÒÐÈ×ÅÑÊÎÃÎ ÂÅÊÒÎÐÀ ËÓÍÛ

  Xt:=X1-Xn;
  Yt:=Y1-Yn;
  Zt:=Z1-Zn;


  Form1.ListBox1.Items[6]:=FloatToStr(Xt)+'         //òîïîöåíòðè÷åñêèé âåêòîð ëóíû Õ';
  Form1.ListBox1.Items[7]:=FloatToStr(Yt)+'         //òîïîöåíòðè÷åñêèé âåêòîð ëóíû Y';
  Form1.ListBox1.Items[8]:=FloatToStr(Zt)+'         //òîïîöåíòðè÷åñêèé âåêòîð ëóíû Z';

  //Ðàçëîæåíèå òîïîöåíòðè÷åñêîãî Âåêòîðà íà òîïîöåíòðè÷åñêèå êîîðäèíàòû

  Atpc:=ArcTan2(Yt,Xt);
  Btpc:=ArcTan2(Zt,sqrt((Xt*Xt)+(Yt*Yt)));


  //Âûñîòà Ëóíû íàä ãîðèçîíòîì
  Pil:=180/pi;

  LunHH:=Pil*ArcSin(sin(fshir)*sin(Btpc)+cos(fshir)*cos(Btpc)*cos(sl-Atpc))+0.85;

  LunHHH:=LunHH;

  Form1.ListBox1.Items[9]:=FloatToStr(Atpc)+'          //Òîïîöåíòð. êîîðäèíàòû';
  Form1.ListBox1.Items[10]:=FloatToStr(Btpc)+'         //Òîïîöåíòð. êîîðäèíàòû';


  Form1.ListBox1.Items[11]:=FloatToStr(LunHHH)+'         //Âûñîòà ëóíû íàä ãîðèçîíòîì';
  Form1.ListBox1.Items[12]:='_____________________________________________________';
  Form1.ListBox1.Items[13]:='_____________________________________________________';

  end;

procedure TForm1.Button1Click(Sender: TObject);
var
 k:integer;
 Tk:double;
begin
  for k := 0 to N do
 begin
   Tk:=k*Atx ;
   MD:=MD+(Tk/24);
   TheLunCalc;
   Form1.ListBox1.Items[14+k]:=FloatTostr(Tk)+'                //'+FloatTostr(LunHHH);// ïå÷àòü äëÿ íàãëÿäíîñòè è ïðîâåðêè
 end;
   
end;

procedure TForm1.FormCreate(Sender: TObject);
 begin
 // Âû÷èñëÿåì ìîäèôèöèðîâàííóþ Þëèàíñêóþ äàòó íà ìîìåíò ðàñ÷åòà

  Year:=2010;
  Month:=10;
  Day:=7;
  Hour:=0;
  Minute:=0;

   Atx:=0.5;
   N:=round(24/Atx);

  q1:=int((Month+9)/12);
  q2:=int(7*((Year+q1)/4));
  q3:=int((275*Month)/9)+Day;
  Time:=Hour + Minute/60;
  MD:=367*Year-q2+q3-678987+(Time/24);


end;

[Tau]

Смотрите, что у вас присходит:
k= 0  Tk= 0  Md= 55476
k= 1  Tk= .5  Md= 55476.0208333333
k= 2  Tk= 1  Md= 55476.0625
k= 3  Tk= 1.5  Md= 55476.125
k= 4  Tk= 2  Md= 55476.2083333333
k= 5  Tk= 2.5  Md= 55476.3125
k= 6  Tk= 3  Md= 55476.4375
k= 7  Tk= 3.5  Md= 55476.5833333333
k= 8  Tk= 4  Md= 55476.75
k= 9  Tk= 4.5  Md= 55476.9375
k= 10  Tk= 5  Md= 55477.1458333333
За 10 шагов (5 часов) вы прошли целые сутки (MD стал не 55476, а 55477)
На каждом k-ом шаге цикла вы прибавляете к предыдущей юлианской дате  k шагов по времени.

Запомните начальный момент времени в процедуре FormCreate
MD:=367*Year-q2+q3-678987+(Time/24);
MD0:=MD;

А в процедуре нажатия на кнопку прибавляйте k шагов (Tk/24) не к предыдущей дате, а к начальной
MD:=MD0+(Tk/24);

atpc=-3.129757

Угол - это циклическая переменная, кратная 360 градусам. То есть 20 градусов это все равно что 20+360=380 градусов и все равно что 20-360 =-340 градусов. Плюс-минус 360 градусов к углу сам угол не меняют. Соответственно, если вы прибавите к -3.129757 величину в 2*pi (360 градусов в радианах), то получите 3.153429 радиан.

[RN3QGA]

  Ох, я не знаю  чтоб я без Вас делал Tau. Все заработало, теперь осталось разобраться с моментами пересечения Луной линии горизонта.
Да кстати H_k-1 просто H или опять приводить к чему нибудь надо Ну в смысле положительно или отрицательного? Уж простите лучше сразу спрошу, что опять не забивать голову Вам.

[Tau]

H - высота Луны в градусах. Положительное H - Луна над горизонтом, отрицательное H - Луна под горизонтом. Ничего больше переводить не нужно.  Запоминаете в массиве все H и потом в отдельном цикле ищете, где у соседних H меняется знак, и вычисляете время пересечения Луной горизонта.