home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / CPM / TURBOM2 / M2ASTRO.LBR / ASTROLIB.MZD / ASTROLIB.MOD
Text File  |  2000-06-30  |  6KB  |  282 lines
  1. IMPLEMENTATION MODULE AstroLib;
  2.  
  3. (* Jim Lill  August, 1987 *)
  4.  
  5. FROM LongMath IMPORT Cos, Sin;
  6. FROM LongTrig IMPORT Arccos, Arcsin, Arctan2, Tan, Rad, Deg;
  7.  
  8. CONST Epoch0 = 1980;
  9.       mwdec = 27.4D0;     (* Milky Way north Galactic Pole *)
  10.       mwra = 192.25D0;    (* 1950 Coordinates *)
  11.          D24 = 24.0D0;
  12.          D60 = 60.0D0;
  13.           D0 = 0.0D0;
  14.          D15 = 15.0D0;
  15.  
  16. VAR rmwdec,rmwra,oe, Result : LONGREAL;
  17.  
  18. (***************************************)
  19.  
  20. PROCEDURE DLONG(dnum:LONGREAL):LONGREAL;
  21.  
  22. (* 8/16/87 *)
  23.  
  24. BEGIN
  25.   Result:= DOUBLE(LONG(dnum));
  26.   RETURN Result
  27. END DLONG;
  28.  
  29. PROCEDURE JD(yr,mo,day:CARDINAL; hrs:LONGREAL):LONGREAL;
  30.  
  31. (* 8/15/87 *)
  32.  
  33. VAR  D, C, d : LONGREAL;
  34.      A, B : INTEGER;
  35.  
  36. BEGIN
  37.   IF mo < 3 THEN
  38.     mo:= mo + 12;
  39.     yr:= yr - 1
  40.   END;
  41.   d:= DOUBLE(day) + hrs/D24;
  42.   A:= INT(yr DIV 100);
  43.   B:= (A DIV 4) + 2 - A;
  44.   C:= DLONG(3.6525D2 * DOUBLE(yr));
  45.   D:= DLONG(30.6001D0 * (DOUBLE(mo) + 1.0D0));
  46.   Result:= D + C + d + 1.7209945D6 + DOUBLE(B);
  47.   RETURN Result
  48. END JD;
  49.  
  50. PROCEDURE HRS(hr,min,sec:CARDINAL):LONGREAL;
  51.  
  52. (* 8/14/87 *)
  53.  
  54. VAR h,m,s : LONGREAL;
  55.  
  56. BEGIN
  57.   h:= DOUBLE(hr);
  58.   m:= DOUBLE(min);
  59.   s:= DOUBLE(sec);
  60.   Result:= h + (m/D60) + (s/(D60 * D60));
  61.   RETURN Result
  62. END HRS;
  63.  
  64. PROCEDURE HMS(VAR hr,min,sec:CARDINAL; hrs:LONGREAL);
  65.  
  66. (* 8/14/87 *)
  67.  
  68. VAR r : LONGREAL;
  69.  
  70. BEGIN
  71.   hr:= CARD(hrs);
  72.   r:= D60 * (hrs - DOUBLE(hr));
  73.   min:= CARD(r);
  74.   r:= D60 * (r - DOUBLE(min));
  75.   sec:= CARD(r)
  76. END HMS;
  77.  
  78. PROCEDURE ND(yr:CARDINAL; jd:LONGREAL):LONGREAL;
  79.  
  80. (* 8/16/87 *)
  81.  
  82. BEGIN
  83.   Result:= jd - JD(yr,1,0,D0);
  84.   RETURN Result
  85. END ND;
  86.  
  87. PROCEDURE TC(yr:CARDINAL):LONGREAL;
  88.  
  89. (* 8/16/87 *)
  90.  
  91. BEGIN
  92.   Result:= (JD(yr,1,0,D0) - JD(1900,1,0,12.0D0))/36.525D3;
  93.   RETURN Result
  94. END TC;
  95.  
  96. PROCEDURE GST(nd,tc,hrs:LONGREAL; yr:CARDINAL):LONGREAL;
  97.  
  98. (* 8/15/87 *)
  99.  
  100. CONST Ka = 6.57098D-2;
  101.       Kc = 1.002738D0;
  102.       Kd = 9.97270D-1;
  103.  
  104. VAR R,S,U : LONGREAL;
  105.  
  106. BEGIN
  107.   R:= 6.64606560D0 + (2400.0512620D0 * tc) + (0.00002581D0 * tc * tc);
  108.   U:= R - DOUBLE(24 * (yr - 1900));
  109.   Result:= (nd * Ka) - (D24 - U) +  (DOUBLE(hrs) * Kc);
  110.   IF Result > D24
  111.     THEN Result:= Result - D24
  112.   END;
  113.   IF Result < D0
  114.     THEN Result:= Result + D24
  115.   END;
  116.   RETURN Result
  117. END GST;
  118.  
  119. PROCEDURE LST(gst, long:LONGREAL):LONGREAL;
  120.  
  121. (* 8/15/87 *)
  122.  
  123. BEGIN
  124.   long:= long / D15;          (* offset in hours *)
  125.   Result:= gst - DOUBLE(long);   (* W long only *)
  126.   IF Result > D24
  127.     THEN Result:= Result - D24
  128.   END;
  129.   IF Result < D0
  130.     THEN Result:= Result + D24
  131.   END;
  132.   RETURN Result
  133. END LST;
  134.  
  135. PROCEDURE LHA(ra,lst:LONGREAL):LONGREAL;
  136.  
  137. (* 8/15/87 *)
  138.  
  139. BEGIN
  140.   Result:= lst - ra;
  141.   IF Result < D0
  142.     THEN Result:= Result + D24
  143.   END;
  144.   RETURN Result
  145. END LHA;
  146.  
  147. PROCEDURE RA(lha,lst:LONGREAL):LONGREAL;
  148.  
  149. (* 8/16/87 *)
  150.  
  151. BEGIN
  152.   Result:= lst - lha;
  153.   IF Result < D0
  154.     THEN Result:= Result + D24
  155.   END;
  156.   RETURN Result
  157. END RA;
  158.  
  159. PROCEDURE OE(tc:LONGREAL):LONGREAL;
  160.  
  161. (* 8/16/87 *)
  162.  
  163. BEGIN
  164.   Result:= 46.845D0 + (0.0059D0 * tc * tc) - (0.00181D0 * tc * tc * tc);
  165.   Result:= Result / 3.6D3;
  166.   Result:= Rad(23.452294D0 - Result);
  167.   RETURN Result
  168. END OE;
  169.  
  170. PROCEDURE Equ2Hor(lha,dec,lat:LONGREAL; VAR az,alt:LONGREAL);
  171.  
  172. (* 8/16/87 *)
  173.  
  174. VAR Ap, H, ca, sa : LONGREAL;
  175.  
  176. BEGIN
  177.   H:= Rad(lha * D15);
  178.   dec:= Rad(dec);
  179.   lat:= Rad(lat);
  180.   sa:= (Sin(dec) * Sin(lat)) + (Cos(dec) * Cos(lat) * Cos(H));
  181.   alt:= Arcsin(sa);
  182.   ca:= (Sin(dec) - (Sin(lat) * Sin(alt))) / (Cos(lat) * Cos(alt));
  183.   Ap:= Deg(Arccos(ca));
  184.   IF Sin(H) < D0 THEN
  185.     az:= Ap
  186.   ELSE
  187.     az:= 360.0D0  - Ap
  188.   END;
  189.   alt:= Deg(alt)
  190. END Equ2Hor;
  191.  
  192. PROCEDURE Hor2Equ(az,alt,lat:LONGREAL; VAR lha,dec:LONGREAL);
  193.  
  194. (* 8/16/87
  195.  
  196. This is the opposite of Equ2Hor literally! *)
  197.  
  198. VAR l1, d1 : LONGREAL;
  199.  
  200. BEGIN
  201.   az:= az / D15;
  202.   Equ2Hor(az,alt,lat,l1,d1);
  203.   lha:= l1 / D15;              (* convert to hours *)
  204.   dec:= d1
  205. END Hor2Equ;
  206.  
  207. PROCEDURE Ecl2Equ(elong,elat:LONGREAL; VAR ra,dec:LONGREAL);
  208.  
  209. (* 8/16/87 *)
  210.  
  211. VAR  sd, x, y : LONGREAL;
  212.  
  213. BEGIN
  214.   elong:= Rad(elong);
  215.   elat:= Rad(elat);
  216.   sd:= (Sin(elat) * Cos(oe)) + (Cos(elat) * Sin(oe) * Sin(elong));
  217.   dec:= Deg(Arcsin(sd));
  218.   y:= (Sin(elong) * Cos(oe)) - (Tan(elat) * Sin(oe));
  219.   x:= Cos(elong);
  220.   ra:= Arctan2(x,y);
  221.   ra:= Deg(ra) / D15
  222. END Ecl2Equ;
  223.  
  224. PROCEDURE Equ2Ecl(ra,dec:LONGREAL; VAR elong,elat:LONGREAL);
  225.  
  226. (* 8/16/87 *)
  227.  
  228. VAR sb, x, y : LONGREAL;
  229.  
  230. BEGIN
  231.   dec:= Rad(dec);
  232.   ra:= Rad(ra * D15);
  233.   sb:= (Sin(dec) * Cos(oe)) - (Cos(dec) * Sin(oe) * Sin(ra));
  234.   WRITELN("sb= ",sb);
  235.   elat:= Deg(Arcsin(sb));
  236.   y:= (Sin(ra) * Cos(oe)) + (Tan(dec) * Sin(oe));
  237.   x:= Cos(ra);
  238.   elong:=Deg(Arctan2(x,y))
  239. END Equ2Ecl;
  240.  
  241. PROCEDURE Equ2Gal(ra,dec:LONGREAL; VAR l,b:LONGREAL);
  242.  
  243. (* 8/18/87 *)
  244.  
  245.  
  246. VAR sb,x,y : LONGREAL;
  247.  
  248. BEGIN
  249.   dec:= Rad(dec);
  250.   ra:= Rad(ra * D15);
  251.   sb:= (Cos(dec) * Cos(rmwdec) * Cos(ra - rmwra)) +
  252.        (Sin(dec) * Sin(rmwdec));
  253.   b:= Arcsin(sb);
  254.   y:= Sin(dec) - (Sin(b) * Sin(rmwdec));
  255.   x:= Cos(dec) * Sin(ra - rmwra) * Cos(rmwdec);
  256.   l:= Arctan2(x,y);
  257.   b:= Deg(b);
  258.   l:= Deg(l) + 33.0D0
  259. END Equ2Gal;
  260.  
  261. PROCEDURE Gal2Equ(l,b:LONGREAL; VAR ra,dec:LONGREAL);
  262.  
  263. (* 8/17/87 *)
  264.  
  265. VAR sd,x,y : LONGREAL;
  266.  
  267. BEGIN
  268.   l:= Rad(l - 33.0D0);
  269.   b:= Rad(b);
  270.   sd:= (Cos(b) * Cos(rmwdec) * Sin(l)) + (Sin(b) * Sin(rmwdec));
  271.   dec:= Deg(Arcsin(sd));
  272.   y:= Cos(b) * Cos(l);
  273.   x:= (Sin(b) * Cos(rmwdec)) - (Cos(b) * Sin(rmwdec) * Sin(l));
  274.   ra:= (Deg(Arctan2(x,y) + rmwra)) / D15
  275. END Gal2Equ;
  276.  
  277. BEGIN  (* Initialization *)
  278.   JD0:= JD(Epoch0,0,0,D0);
  279.   oe:= OE(TC(Epoch0));            (*obliquity of the ecliptic, Epoch0 *)
  280.   rmwdec:= Rad(mwdec);
  281.   rmwra:= Rad(mwra);
  282. END AstroLib.