home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / math / daylite / calculat.pas next >
Encoding:
Pascal/Delphi Source File  |  1989-07-17  |  4.5 KB  |  153 lines
  1. unit Calculate;
  2.  
  3. interface
  4.  
  5.   uses Globals;
  6.  
  7.   type
  8.  
  9.     SunTime = record Dawn,Dusk:Float; end;
  10.  
  11.     Daytime = record
  12.       Astronomical,
  13.       Nautical,
  14.       Civil,
  15.       Actual:SunTime;
  16.       end;
  17.  
  18.   procedure Phenomina(lat,long,zone:float; y:longint;m,d:integer;
  19.                       var Day:DayTime);
  20.  
  21.  
  22. implementation
  23.  
  24.   const
  25.         Sunup    = pi/2;
  26.         Sundown  = pi*3/2;
  27.  
  28.  
  29.   function Normalize(z:float):float;
  30.     begin
  31.       z := 2*pi * Frac(z/(2*pi));
  32.       if z<0 then
  33.         z := z + 2*pi;
  34.       Normalize := z;
  35.     end;
  36.  
  37.   procedure CalcSolar(SunAt,lat,long:float; y:longint; m,d:integer;
  38.                   var ApproxTime,
  39.                       SolarRightAscension, SolarDeclination:float);
  40.     var
  41.       SolarMeanAnomaly,
  42.       SolarTrueLong
  43.         :Float;
  44.       n :integer;
  45.     begin
  46.       ApproxTime := DayOfYear(y,m,d) + (SunAt+long)/(2*pi);
  47.       SolarMeanAnomaly := ApproxTime * 0.017202 - 0.0574039;
  48.       SolarTrueLong := Normalize(
  49.                         SolarMeanAnomaly
  50.                         + 0.0334405 * sin(SolarMeanAnomaly)
  51.                         + 3.49066E-4 * sin(2*SolarMeanAnomaly)
  52.                         + 4.93289
  53.                        );
  54.       {Quadrant Determination}
  55.         if Frac(SolarTrueLong*2/pi) = 0.0 then
  56.           SolarTrueLong := SolarTrueLong + 4.84814e-6;
  57.         n := 2;
  58.         if SolarTrueLong <= pi/2 then
  59.           n := 0
  60.         else if SolarTrueLong <= pi*3/2 then
  61.           n := 1;
  62.       SolarRightAscension := ArcTan(0.91746 * Tan(SolarTrueLong));
  63.       {Quadrant Adjustment}
  64.         if n = 1 then
  65.           SolarRightAscension := SolarRightAscension + pi
  66.         else if n = 2 then
  67.           SolarRightAscension := SolarRightAscension + 2*pi;
  68.       {Solar Declination}
  69.         SolarDeclination := 0.39782 * Sin(SolarTrueLong);
  70.         SolarDeclination := ArcTan(SolarDeclination
  71.                                   / Sqrt(1-Sqr(SolarDeclination))
  72.                                   );
  73.     end;
  74.  
  75.   procedure ConvertCoords(r,q,lat:float; var s:float);
  76.     begin
  77.       s := (r-(sin(q)*sin(lat))) / (cos(q)*cos(lat));
  78.     end;
  79.  
  80.   function SettingAdjustment(var SolarDeclination:float):float;
  81.     begin
  82.       SettingAdjustment := pi/2 - ArcTan( SolarDeclination
  83.                                         / Sqrt(1-sqr(SolarDeclination))
  84.                                         );
  85.     end;
  86.  
  87.   function RisingAdjustment(var SolarDeclination:float):float;
  88.     begin
  89.       RisingAdjustment := 2*pi - SettingAdjustment(SolarDeclination);
  90.     end;
  91.  
  92.   procedure GetTime(long,zone:float; y,m,d: integer;
  93.                     ApproxTime,SolarRightAscension, SolarDeclination:float;
  94.                 var t:float);
  95.     begin
  96.       t := {Loval Apparent Time}
  97.              SolarDeclination + SolarRightAscension
  98.                  - ApproxTime * 2*pi/365.2425 - 1.73364
  99.            {Universal Time}
  100.                + long
  101.            {Wall Clock Time}
  102.                - zone;
  103.       t := Normalize(t) * 24{hours}/(2*pi{radians});
  104.     end;
  105.  
  106.  
  107.   procedure Phenomina(lat,long,zone:float; y:longint;m,d:integer;
  108.                       var Day:DayTime);
  109.     var time,ascension,declination:float;
  110.     procedure RisingTime(r:float; var t:SunTime);
  111.       var s:float;
  112.           i:integer;
  113.       begin;
  114.         ConvertCoords(r,declination,lat,s);
  115.         if abs(s)>1.0 then
  116.           t.Dawn := 0
  117.         else
  118.           GetTime(long,zone,y,m,d,time,ascension,RisingAdjustment(s),t.Dawn);
  119.       end;
  120.     procedure SettingTime(r:float; var t:SunTime);
  121.       var s:float;
  122.           i:integer;
  123.       begin;
  124.         ConvertCoords(r,declination,lat,s);
  125.         if abs(s)>1.0 then
  126.           t.Dusk := 0
  127.         else
  128.           GetTime(long,zone,y,m,d,time,ascension,SettingAdjustment(s),t.Dusk);
  129.       end;
  130.     const
  131.        Astro   = -0.309017;
  132.        Nautic  = -0.207912;
  133.        Civil   = -0.104528;
  134.        Actual  = -0.0145439;
  135.     begin
  136.       lat  := Radians(lat);
  137.       long := Radians(long);
  138.       zone := Radians(zone)*15.0; {15 degrees/Time Zone}
  139.  
  140.       CalcSolar(Sunup,lat,long, y,m,d, time,ascension,declination);
  141.       RisingTime(Astro  ,Day.Astronomical);
  142.       RisingTime(Nautic ,Day.Nautical);
  143.       RisingTime(Civil  ,Day.Civil);
  144.       RisingTime(Actual ,Day.Actual);
  145.  
  146.       CalcSolar(Sundown,lat,long, y,m,d, time,ascension,declination);
  147.       SettingTime(Actual ,Day.Actual);
  148.       SettingTime(Civil  ,Day.Civil);
  149.       SettingTime(Nautic ,Day.Nautical);
  150.       SettingTime(Astro  ,Day.Astronomical);
  151.     end;
  152.  
  153. end.