unit sunriseset;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls, Math;             ///remember to add Math

type
  TForm1 = class(TForm)
    Edit1: TEdit;
    Edit2: TEdit;
    Edit3: TEdit;
    Edit4: TEdit;
    Edit5: TEdit;
    Button1: TButton;
    Button2: TButton;
    Edit9: TEdit;
    Edit10: TEdit;
    Edit11: TEdit;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label5: TLabel;
    Label6: TLabel;
    Label7: TLabel;
    Label8: TLabel;
    Label9: TLabel;
    Label10: TLabel;
    Label11: TLabel;
    Label12: TLabel;
    Label13: TLabel;
    ListBox1: TListBox;
    procedure Button2Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.DFM}

procedure TForm1.Button2Click(Sender: TObject);
var
d : real ;
T : real ;
K1 : real ;
K2 : real ;
K3 : real ;
K4 : real ;
DegToRad, K, N, V, T1, T2 : real;
Tb, Kb1, Kb2, Kb3, Nb, Vb, T1b, T2b : real;
year1 : Integer ;
day1 : integer;
month1:integer;
c:real;
leap : integer;
D9:integer;
daynum:integer;
d1:integer;
cnt:integer;
d2:real;
tb1,tb2:real;
Pday,Pmonth,PD1:integer;
f9:integer;
SRA,SSA,SRB,SSB,DAY,DN:string;

begin
V := strtofloat(Edit1.Text) ; //'longitude A station
Vb := strtofloat(Edit10.Text);  //longitude B station
N := strtofloat(Edit2.Text); //'latitude  A station
Nb := strtofloat(Edit11.Text);//'latitude B station
c := strtofloat(Edit3.Text) ;
day1 := strtoint(Edit5.Text) ;
month1 := strtoint(Edit4.Text);
year1 := strtoint(Edit9.Text) ;

//clear the current sunrise/sunset times
ListBox1.Clear;

//'routine to calculate day of year number

leap := 0;
If (year1 / 4) = Int(year1 / 4) Then
leap := 1 ;

   Case month1 of
   1 : D9 := 0;
   2 : D9 := 31;
   3 : D9 := 59;
   4 : D9 := 90;
   5 : D9 := 120;
   6 : D9 := 151;
   7 : D9 := 181;
   8:  D9 := 212;
   9 : D9 := 243;
   10 : D9 := 273;
   11 : D9 := 304;
   12 : D9 := 334;
   End;
daynum := day1 + D9;
If month1 > 2 Then
daynum := leap + daynum;

K := 3.1415927 / 180;
DegToRad := 3.1415927 / 180;

d1 := daynum ;

For cnt := 0 To 31 do//'number of calculations
 begin
     d2 := 23.45 * Sin((360 / 365 * (284) + d1) * DegToRad);
      d := -23.5 * Cos((10 + 0.988 * (day1 + 30.33 * (month1 - 1))) * DegToRad);
        //'declination of the sun
           d := (d + d2) / 2 ;// 'average two methods of determining sun's declination
     //'this is a simplification!  Results in sunset/sunrise times plus/minus several
      //minutes

//' correction factors (each is a segment of the equation of time curve)
{
If d1 < 107 Then
K4 := 15 * Sin((107 - d1) * 180 * K / 117)
Else If d1 < 160 Then
K4 := -3.5 * Sin((160 - d1) * 180 * K / 53)
Else If d1 < 247 Then
K4 := 8 * Sin((247 - d1) * 180 * K / 87)
Else If d1 >= 247 Then
K4 := -14.3 * Sin(355 - d1) * 180 * K / 108
}  //above ifs or below case statements work equally well
 case  d1 of
        1..107:K4 := 15 * Sin((107 - d1) * 180 * K / 117);
        108..160 :K4 := -3.5 * Sin((160 - d1) * 180 * K / 53);
        161..246 :K4 := 8 * Sin((247 - d1) * 180 * K / 87);
        247..366 :K4 := -14.3 * Sin(355 - d1) * 180 * K / 108 ;
 end;

K4 := K4 / 60;
 
//' time  for station A
T := Tan(N * K) * Tan(d * K);
If (T * T) > 1 Then  exit;
T := (1.5708 - Arctan(T / Sqrt(1 - (T * T)) )) / (15 * K);
//'correction
K1 := Sin(N * K) / Cos(d * K);
K1 := 1.5708 - Arctan(K1 / Sqrt(1 - (K1 * K1)));
K2 := Sin(0.835608 * K) / Sin(K1);
K2 := Arctan(K2 / Sqrt(1 - (K2 * K2)));
K2 := K2 / K ;
K3 := 240 * K2 / (3600 * Cos(d * K));
//'calculate sunrise/sunset times
T1 := V / 15 + T - K3 + K4;
If T1 < 0 Then T1 := T1 + 24;
If T1 > 24 Then T1 := T1 - 24;
T2 := V / 15 - T + K3 + K4;
If T2 < 0 Then T2 := T2 + 24 ;
If T2 > 24 Then  T2 := T2 - 24 ;
T1 := Int(T1) * 100 + Int((T1 - Int(T1)) * 60);
T2 := Int(T2) * 100 + Int((T2 - Int(T2)) * 60);

//'calculate second station B
Tb := Tan(Nb * K) * Tan(d * K);
If Tb * Tb > 1 Then exit;
Tb := (1.5708 - ArcTan(Tb / Sqrt(1 - (Tb *Tb)  ))) / (15 * K);
//'correction
Kb1 := Sin(Nb * K) / Cos(d * K);
Kb1 := 1.5708 - ArcTan(Kb1 / Sqrt(1 - (Kb1 * Kb1)));
Kb2 := Sin(0.835608 * K) / Sin(Kb1) ;
Kb2 := ArcTan(Kb2 / Sqrt(1 - (Kb2 *Kb2)));
Kb2 := Kb2 / K ;
Kb3 := 240 * Kb2 / (3600 * Cos(d * K));
//'calculate sunrise/sunset times
Tb1 := Vb / 15 + Tb - Kb3 + K4;
If Tb1 < 0 Then Tb1 := Tb1 + 24;
If Tb1 > 24 Then Tb1 := Tb1 - 24;
Tb2 := Vb / 15 - Tb + Kb3 + K4;
If Tb2 < 0 Then Tb2 := Tb2 + 24 ;
If Tb2 > 24 Then Tb2 := Tb2 - 24;
Tb1 := Int(Tb1) * 100 + Int((Tb1 - Int(Tb1)) * 60);
Tb2 := Int(Tb2) * 100 + Int((Tb2 - Int(Tb2)) * 60);

 Pday := day1 ;
 Pmonth := month1;
 pd1 := d1 ;

Case Pmonth of
   1: If day1 > 31 Then
    begin
    month1 := month1 + 1;
    Pmonth := month1;
    day1 := 1 ;
    Pday := day1;
    End;
   3:If day1 > 31 Then
    begin
    month1 := month1 + 1;
    Pmonth := month1;
    day1 := 1 ;
    Pday := day1 ;
    End;
   5:If day1 > 31 Then
   begin
    month1 := month1 + 1 ;
    Pmonth := month1 ;
    day1 := 1 ;
    Pday := day1;
    end;

   7: If day1 > 31 Then
   begin
    month1 := month1 + 1 ;
    Pmonth := month1 ;
    day1 := 1 ;
    Pday := day1 ;
    End;
   8: If day1 > 31 Then
   begin
    month1 := month1 + 1 ;
    Pmonth := month1 ;
    day1 := 1 ;
    Pday := day1 ;
    End;
   10: If day1 > 31 Then
   begin
    month1 := month1 + 1 ;
    Pmonth := month1 ;
    day1 := 1 ;
    Pday := day1 ;
    End;
    12: If day1 > 31 Then
    begin
    month1 := month1 - 11;
    Pmonth := month1 ;
    day1 := 1 ;
    Pday := day1;
    End;
   4: If day1 > 30 Then
    begin
    month1 := month1 + 1;
    Pmonth := month1;
    day1 := 1;
    Pday := day1 ;
    End;
   6 :If day1 > 30 Then
   begin
    month1 := month1 + 1;
    Pmonth := month1;
    day1 := 1 ;
    Pday := day1;
    End;

   9:If day1 > 30 Then
   begin
    month1 := month1 + 1;
    Pmonth := month1;
    day1 := 1;
    Pday := day1;
    End;

   11: If day1 > 30 Then
   begin
    month1 := month1 + 1;
    Pmonth := month1;
    day1 := 1;
    Pday := day1;
    End;

   2 :If (day1 > 29) And (f9 = 1) Then
   begin
    month1 := 3;
    Pmonth := month1;
    day1 := 1;
    Pday := day1 ;
    end
  Else If (day1 > 28) And (f9 = 0) Then
  begin
    month1 := 3;
    Pmonth := month1;
    day1 := 1;
    Pday := day1;
  End;

Else

End;

SRA := floattostr(T1);        //routines to format output
SSA := floattostr(T2);
SRB := floattostr(Tb1);
SSB := floattostr(Tb2);
day := floattostr(Pmonth)+ '/'+ floattostr( Pday) ;
DN := floattostr(pd1);

if length(DAY)= 3 then
DAY := DAY+'    '
else if length(Day)= 4 then
DAY := DAY +'  '
else if length(DAY)=5 then
DAY := DAY;

if length(SRA)=1 then
SRA := '000'+ SRA
else if length(SRA)=2 then
SRA := '00'+ SRA
else if length(SRA)=3 then
SRA :='0'+ SRA ;

if length(SSA)=1 then
SSA := '000'+ SSA
else if length (SSA)=2 then
SSA := '00'+ SSA
else if length(SSA)=3 then
SSA :='0'+ SSA;

if length(SRB)=1 then
SRB := '000'+ SRB
else if length( SRB)=2 then
SRB := '00'+ SRB
else if length(SRB)=3 then
SRB :='0'+ SRB;
if length(SSB)=1 then
SSB := '000'+ SSB
else if length(SSB)=2 then
SSB := '00'+ SSB
else if length(SSB)=3 then
SSB :='0'+ SSB;
ListBox1.items.add   (SRA +'          '+SSA  +'          '+DAY +'          '+SRB +'          '+SSB +'          '+'   '+ floattostr( pd1) );

d1 := d1 + 1 ;// 'increment the day number
day1 := day1 + 1;//   'increment the day
if d1>366 then exit;   //only a years worth
Edit3.Text := inttostr(d1 - 1 );

end;
end;
end.
