ot: analemma

[Ales Prochaska, Divesoft]

{ Vypocet polohy Slunce }

uses
  datim,graf0,m2;

type
  real         = extended;

const
  lat          = 50.0875;
  lon          = +14.4714;

{----------------------------------------------------------------------------}
function daynumber (Y,M,D : longint; UT : real) : real;
{ Plati pro roky 1901 - 2099                                                 }
{----------------------------------------------------------------------------}
begin
  daynumber:=367*Y-int(7*(Y+int((M+9)/12))/4)+int(275*M/9)+D-730530+UT/24;
  end;



{----------------------------------------------------------------------------}
function arcsin (x : real) : real;
{ ve stupnich                                                                }
{----------------------------------------------------------------------------}
begin
  arcsin:=arctan(x/sqrt(1-sqr(x)))/pi*180;
  end;



{----------------------------------------------------------------------------}
function sin (x : real) : real;
{ ve stupnich                                                                }
{----------------------------------------------------------------------------}
begin
  sin:=system.sin(x/180*pi);
  end;



{----------------------------------------------------------------------------}
function cos (x : real) : real;
{ ve stupnich                                                                }
{----------------------------------------------------------------------------}
begin
  cos:=system.cos(x/180*pi);
  end;



{----------------------------------------------------------------------------}
function tan (x : real) : real;
{ ve stupnich                                                                }
{----------------------------------------------------------------------------}
begin
  tan:=sin(x)/cos(x)
  end;



{----------------------------------------------------------------------------}
function arctan (x : real) : real;
{ ve stupnich                                                                }
{----------------------------------------------------------------------------}
begin
  arctan:=system.arctan(x)/pi*180;
  end;



{----------------------------------------------------------------------------}
function atan2 (y,x : real) : real;
{ ve stupnich                                                                }
{----------------------------------------------------------------------------}
begin
  if x>0 then atan2:=arctan(y/x)
  else if x<0 then atan2:=arctan(y/x)+180
  else if y>0 then atan2:=90
  else if y<0 then atan2:=-90
  else atan2:=0;
  end;



{----------------------------------------------------------------------------}
function norm (x : real) : real;
{ normalizuje uhel do intervalu (-180, 180>                                  }
{----------------------------------------------------------------------------}
begin
  { jak od blbyho (co byste taky cekali) }
  while x>180 do x:=x-360;
  while x<=-180 do x:=x+360;
  norm:=x;
  end;



{----------------------------------------------------------------------------}
procedure poloha     (    lat   : real;     { sirka                          }
                          lon   : real;     { delka                          }
                          year  : integer;  { rok                            }
                          month : integer;  { mesic                          }
                          day   : integer;  { den                            }
                          UT    : real;     { cas                            }
                      var az    : real;     { azimut                         }
                      var h     : real);    { vyska                          }
{ Vypocita azimut a vysku Slunce pro zadane misto v zadanem okamziku.        }
{ Uhly jsou ve stupnich a desetinach, cas v hodinach a desetinach.           }
{----------------------------------------------------------------------------}
var
  d            : real;         { day number                                  }
  lst          : real;         { local sidereal time (ve stupnich)           }
  oblecl       : real;         { obliquity of the ecliptic                   }
  w            : real;         { argument of perihelion                      }
  M            : real;         { mean anomaly (ve stupnich)                  }
  e            : real;         { eccentricity                                }
  EA           : real;         { eccentric anomaly (ve stupnich)             }
  A,B          : real;         { pracovni promenne                           }
  r            : real;         { distance                                    }
  v            : real;         { true anomaly                                }
  slon         : real;         { Sun's true longitude                        }
  sRA          : real;         { Sun's right ascension                       }
  sDec         : real;         { Sun's declination                           }
  t_sRA        : real;         { tan(sRA)                                    }
  s_sDec       : real;         { sin(sDec)                                   }
  HA           : real;         { hour angle                                  }

begin {poloha}
  d:=daynumber(year,month,day,UT);
  lst:=norm(98.9818+0.985647352*d+UT*15+lon);
  oblecl:=23.4393-3.563e-7*d;
  w:=282.9494+4.70935e-5*d;
  M:=356.0470+0.98566002585*d;
  e:=0.016709-1.151e-9*d;
  EA:=M+e*(180/pi)*sin(M)*(1+e*cos(M));
  A:=cos(EA)-e;
  B:=sqrt(1-e*e)*sin(EA);
  r:=sqrt(A*A+B*B);
  v:=atan2(B,A);
  slon:=v+w;
  t_sRA:=(sin(slon)*cos(oblecl))/cos(slon);
  sRA:=atan2(sin(slon)*cos(oblecl),cos(slon));
  s_sDec:=sin(oblecl)*sin(slon);
  sDec:=arcsin(s_sDec);
  HA:=lst-sRA;
  az:=atan2(sin(HA),cos(HA)*sin(lat)-tan(sDec)*cos(lat))+180;
  h:=arcsin(sin(lat)*sin(sDec)+cos(lat)*cos(sDec)*cos(HA));
  end {poloha};