Advertisement
frithbuck

Planetarium Project

Apr 4th, 2023
672
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 7.99 KB | Science | 0 0
  1. MODULE ORBITALTRIG
  2.  
  3. REAL X, RADEG
  4. PARAMETER(RADEG=57.2957795130823)
  5.  
  6. CONTAINS
  7.  
  8. FUNCTION SIND(X)
  9. SIND = SIN( X * (1.0/RADEG) )
  10. END FUNCTION
  11.  
  12. FUNCTION ATAND(X)
  13. ATAND = RADEG * ATAN(X)
  14. END FUNCTION
  15.  
  16. FUNCTION ATAN2D(Y,X)
  17. ATAN2D = RADEG * ATAN2(Y,X)
  18. END FUNCTION
  19.  
  20. FUNCTION REV(X)
  21. REV = X - AINT(X/360.0)*360.0
  22. IF (REV.LT.0.0)  REV = REV + 360.0
  23. END FUNCTION
  24.  
  25. FUNCTION CBRT(X)
  26. IF (X.GE.0.0) THEN
  27. CBRT = X ** (1.0/3.0)
  28. ELSE
  29. CBRT = -((-X)**(1.0/3.0))
  30. ENDIF
  31. END FUNCTION
  32.  
  33. FUNCTION COSD(X)
  34. COSD = COS( X * (1.0/RADEG))
  35. END FUNCTION
  36.  
  37. FUNCTION TAND(X)
  38. TAND = TAN(X * (1.0/RADEG))
  39. END FUNCTION
  40.  
  41. FUNCTION ASIND(X)
  42. ASIND = RADEG * ASIN (X)
  43. END FUNCTION
  44.  
  45. FUNCTION ACOSD(X)
  46. ACOSD = RADEG * ACOS (X)
  47. END FUNCTION
  48.  
  49. END MODULE ORBITALTRIG
  50.    
  51. PROGRAM PLANETARIUM
  52. USE ORBITALTRIG
  53. IMPLICIT NONE
  54. REAL, PARAMETER :: PI = 3.14159263585
  55. !TIME PARAMETERS
  56. INTEGER YEAR, MONTH, DATE !USER INPUT YEAR, MONTH, AND DAY
  57. INTEGER HOUR, MINUT, SECON, ZONE !USER LOCAL TIME (UNIT=HOUR), USER LOCAL TIME (UNIT=MINUTE), USER LOCAL TIME (UNIT=SECOND), USER TIMEZONE (UTC OFFSET) (UNIT=HOUR)
  58. REAL EPOCH !TIME SINCE MIDNIGHT JAN 1ST 2000 (UNIT=DAYS)
  59. REAL UT !GREENWICH MEAN TIME
  60. REAL GMST0 !SIDERIAL TIME AT GREENWICH MERIDIAN AT 00:00
  61. REAL SIDTIME !LOCAL SIDERIAL TIME
  62. !OBSERVER TERRESTRIAL COORDINATES
  63. REAL OLON !OBSERVER LONGITUDE
  64. REAL OLAT !OBSERVER LATITUDE'
  65. !THESE ARE FOR CONVERTING DATE_AND_TIME OUTPUT TO FEED EPOCH
  66. INTEGER V(8), YYYY, MM, DD, HH, MINS, SS, MSS, GMT, DAYS
  67. EQUIVALENCE (V(1),YYYY), (V(2),MM), (V(3),DD), &
  68. (V(4),GMT),  (V(5),HH), (V(6),MINS), &
  69. (V(7),SS), (V(8),MSS)
  70. !RECTANGULAR COORDINATES FOR THE ELIPTIC PLANE
  71. REAL XECLIP
  72. REAL YECLIP
  73. REAL ZECLIP
  74. !RECTANGULAR COORDINATES FOR THE EQUATORIAL PLANE
  75. REAL XEQUAT
  76. REAL YEQUAT
  77. REAL ZEQUAT
  78. !RECTANGULAR COORDINATES FOR OBSERVER
  79. REAL XHOR
  80. REAL YHOR
  81. REAL ZHOR
  82. REAL AZI
  83. REAL ALT
  84. !PLACEHOLDERS BETWEEN AXIS ROTATIONS
  85. REAL XHOR1
  86. REAL YHOR1
  87. REAL ZHOR1
  88. !POLAR COORDINATES
  89. REAL RIASC !RIGHT ASCENSION (UNIT=DEGREES)
  90. REAL DECLIN !DECLINATION (UNIT=DEGREES)
  91. !DERIVED PARAMETERS:
  92. REAL PERDIS !PERIAPSIS DISTANCE (UNIT=AU)
  93. REAL APODIS !APOAPSIS DISTANCE (UNIT=AU)
  94. REAL ORBDAY !ORBITAL PERIOD (UNIT=DAYS)
  95. REAL DAYMOT !DAILY MOTION (UNIT=DEGREES PER DAY)
  96. REAL LASPER !TIME SINCE LAST PERIAPSIS (UNIT=DAYS)
  97. REAL TRUANO !TRUE ANOMALY (UNIT=DEGREES)
  98. REAL MENANO !MEAN ANOMALY (UNIT=DEGREES)
  99. REAL MENLON !MEAN LONGITUDE (UNIT=DEGREES)
  100. REAL OBLECL !OBLEQUITY OF THE ECLIPTIC
  101. REAL ECCANO !ECCENTRIC ANOMALY
  102. CHARACTER :: OBJTYP !OBJECT TYPE (PLANET, STAR, COMET, ASTEROID, MOON)
  103. CHARACTER STAR, PLANET, ASTEROID, COMET, MOON
  104. REAL :: HELDIS !HELIOCENTRIC DISTANCE (UNIT=AU)
  105. REAL :: SOLMAS !SOLAR MASSES (UNITLESS)
  106. REAL HORAN ! HOUR ANGLE
  107. !KEPLERIAN ORBITAL PARAMETERS:
  108. REAL :: ARGPER !ARGUMENT OF PERIAPSIS (BIG OMEGA) (UNIT=DEGREES)
  109. REAL :: SMAJAX !SEMI-MAJOR AXIS (SMALL A) (UNIT=AU)
  110. REAL :: ECCEN !ECCENTRICITY (SMALL E) (UNITLESS)
  111. REAL :: INCLI !INCLINATION (UNIT=DEGREES)
  112. REAL :: RAOAN !RIGHT ASCENSION OF ASCENDING NODE (UNIT=DEGREES)
  113. REAL :: PERIDAY !TIME AT PERIAPSIS (BIG T) (UNIT=DAYS)
  114.  
  115. CHARACTER (LEN=15) OBJNAME
  116. CHARACTER (LEN=15) LOCNAME
  117. CHARACTER YESNO
  118.  
  119. WRITE (*,*) 'USE CURRENT DATE AND TIME? <Y/N>'
  120. READ (*,*) YESNO
  121. IF (YESNO.EQ.'N') THEN
  122. WRITE (*,*) 'ENTER DATE <YYYY MM DD>'
  123. READ (*,*) YEAR,MONTH,DATE
  124. EPOCH = 367*YEAR-7*(YEAR+(MONTH+9)/12)/4-3*((YEAR+(MONTH-9)/7)/100+1)/4+275*MONTH/9+DATE-730515
  125. WRITE (*,*) 'ENTER LOCAL TIME <HH MM SS>'
  126. READ (*,*) HOUR, MINUT, SECON
  127. WRITE (*,*) 'ENTER UTC OFFSET, PST=-8, EST=-5, GMT=0, EET=2'
  128. READ (*,*) ZONE
  129. EPOCH = EPOCH + ((HOUR+ZONE) / 24.0) + (MINUT / 1440) + (SECON / 86400.0)
  130. WRITE (*,'(A, TR3, I2, A, I2, A, I5)') 'THE JULIAN DAY NUMBER FOR: ', DATE, '/', MONTH, '/', YEAR
  131. WRITE (*,"(T28, A, TR1, I2.1, A, I2.1, A, I2.2, A, I4.2, A)") 'AT', HOUR, ':', MINUT, ':', SECON,' UTC',ZONE,' IS:'
  132. WRITE (*,"(T28, F15.5)") EPOCH
  133. ELSE
  134. CALL DATE_AND_TIME (VALUES=V)
  135. DAYS = 367*YYYY - 7 * (YYYY + (MM+9)/12 ) / 4 - 3 * ( (YYYY + (MM-9)/7 ) / 100 + 1 ) / 4 + 275*MM/9 + DD - 730515
  136. EPOCH = DAYS + ((HH + GMT)/ 24.0) + (MINS / 1440) + (SS / 86400.0) + (MSS / 86400000.0) + GMT/24.0
  137. WRITE (*,'(A, TR3, I2.2, A, I2.2, A, I4.4)') 'THE JULIAN DAY NUMBER FOR: ', DD, '/', MM, '/', YYYY
  138. WRITE (*,"(T28, A, TR1, I2.2, A, I2.2, A, I2.2, A)") 'AT', HH, ':', MINS, ':', SS,' UTC+0 IS:'
  139. WRITE (*,"(T28, F15.10)") EPOCH
  140. END IF
  141. UT = ZONE + HOUR + (MINUT/60) + (SECON/3600)
  142. WRITE (*,'(A,T25,A)') 'CHOOSE A LOCATION:', 'OPTION'
  143. WRITE (*,'(A,T25,A)') 'SAN FRANCISCO', '{1}'
  144. WRITE (*,'(A,T25,A)') 'NEW YORK', '{2}'
  145. WRITE (*,'(A,T25,A)') 'LONDON', '{3}'
  146. WRITE (*,'(A,T25,A)') 'MOSCOW', '{4}'
  147. WRITE (*,'(A,T25,A)') 'CAIRO', '{5}'
  148. WRITE (*,'(A,T25,A)') 'NEW DELHI', '{6}'
  149. WRITE (*,'(A,T25,A)') 'HONG KONG', '{7}'
  150. WRITE (*,'(A,T25,A)') 'TOKYO', '{8}'
  151. WRITE (*,'(A,T25,A)') 'HONOLULU', '{9}'
  152. WRITE (*,'(A,T25,A)') 'SYDNEY', '{10}'
  153. WRITE (*,'(A,T25,A)') 'CAPE TOWN', '{11}'
  154. WRITE (*,'(A,T25,A)') 'BUENOSAIRES', '{12}'
  155. WRITE (*,'(A,T25,A)') '<ENTER COORDINATES>', '{13}'
  156. READ (*,*) LOCNAME
  157. IF (LOCNAME.EQ.'1') THEN
  158. OLAT = 37.7775
  159. OLON = -122.416389
  160. ELSE IF (LOCNAME.EQ.'2') THEN
  161. OLAT = 40.712778
  162. OLON = -74.006111
  163. ELSE IF (LOCNAME.EQ.'3') THEN
  164. OLAT = 51.507222
  165. OLON = -0.1275
  166. ELSE IF (LOCNAME.EQ.'4') THEN
  167. OLAT = 55.755833
  168. OLON = 37.617222
  169. ELSE IF (LOCNAME.EQ.'5') THEN
  170. OLAT = 30.044444
  171. OLON = 31.235833
  172. ELSE IF (LOCNAME.EQ.'6') THEN
  173. OLAT = 28.613895
  174. OLON = 77.209006
  175. ELSE IF (LOCNAME.EQ.'7') THEN
  176. OLAT = 22.280792
  177. OLON = 114.165578
  178. ELSE IF (LOCNAME.EQ.'8') THEN
  179. OLAT = 35.689722
  180. OLON = 139.692222
  181. ELSE IF (LOCNAME.EQ.'9') THEN
  182. OLAT = 21.306944
  183. OLON = -157.858333
  184. ELSE IF (LOCNAME.EQ.'10') THEN
  185. OLAT = -33.867778
  186. OLON = 151.21
  187. ELSE IF (LOCNAME.EQ.'11') THEN
  188. OLAT = -33.925278
  189. OLON = 18.423889
  190. ELSE IF (LOCNAME.EQ.'12') THEN
  191. OLAT = -34.603333
  192. OLON = -58.381667
  193. ELSE IF  (LOCNAME.EQ.'13') THEN
  194. WRITE (*,*) 'ENTER LATITUDE <DDD MM SS>'
  195. READ (*,*) HOUR, MINUT, SECON
  196. OLAT = HOUR + (MINUT/60) + (SECON/3600)
  197. WRITE (*,*) 'ENTER LONGITUDE <DDD MM SS>'
  198. READ (*,*) HOUR, MINUT, SECON
  199. OLON = HOUR + (MINUT/60) + (SECON/3600)
  200. ELSE
  201. WRITE (*,*) 'INVALID ENTRY'
  202. WRITE (*,*) 'ENTER LATITUDE <DDD MM SS>'
  203. READ (*,*) HOUR, MINUT, SECON
  204. OLAT = HOUR + (MINUT/60) + (SECON/3600)
  205. WRITE (*,*) 'ENTER LONGITUDE <DDD MM SS>'
  206. READ (*,*) HOUR, MINUT, SECON
  207. OLON = HOUR + (MINUT/60) + (SECON/3600)
  208. END IF
  209. WRITE (*,*) LOCNAME
  210. WRITE (*,*) 'LON: ', OLON
  211. WRITE (*,*) 'LAT: ', OLAT
  212. ARGPER = 282.9404 + (4.70935E-5   * EPOCH) ! ARGUMENT OF PERIHELION
  213. SMAJAX = 1.000000 ! MEAN DISTANCE, I.E. SEMIMAJOR AXIS
  214. ECCEN = 0.016709 - (1.151E-9 * EPOCH) ! ECCENTRICITY  
  215. MENANO = REV(356.0470 + (0.9856002585 * EPOCH)) ! MEAN ANOMALY
  216. OBLECL = 23.4393 - (0.0000003563 * EPOCH) ! OBLEQUITY OF ECLIPTIC
  217. MENLON = ARGPER + MENANO
  218.  
  219. ECCANO = MENANO + (180/PI) * ECCEN * SIND(MENANO) * (1 + ECCEN * COSD(MENANO))
  220. XECLIP = COSD(ECCANO) - ECCEN
  221. YECLIP = SIND(ECCANO) * SQRT(1 - ECCEN*ECCEN)
  222. HELDIS = SQRT(XECLIP*XECLIP + YECLIP*YECLIP)
  223. TRUANO = ATAN2D(YECLIP, XECLIP)
  224. RIASC = TRUANO + ARGPER
  225.  
  226. XECLIP = HELDIS * COSD(RIASC)
  227. YECLIP = HELDIS * SIND(RIASC)
  228. ZECLIP = 0.0
  229.  
  230. XEQUAT = XECLIP
  231. YEQUAT = YECLIP * COSD(OBLECL) - ZECLIP * SIND(OBLECL)
  232. ZEQUAT = YECLIP * SIND(OBLECL) + ZECLIP * COSD(OBLECL)
  233.  
  234. HELDIS = SQRT(XEQUAT*XEQUAT + YEQUAT*YEQUAT + ZEQUAT*ZEQUAT )
  235. RIASC = ATAN2D( YEQUAT, XEQUAT )
  236. DECLIN = ATAN2D(ZEQUAT, SQRT( XEQUAT*XEQUAT + YEQUAT*YEQUAT))
  237.  
  238. GMST0 = (MENLON/15.0 + 12.0) - AINT((MENLON/15.0 + 12.0)/24)*24
  239. IF (GMST0.LT.0.0)  GMST0 = GMST0 + 24
  240. SIDTIME = GMST0 + UT + OLON/15.0
  241. HORAN = (SIDTIME*15) - RIASC
  242.  
  243. XHOR1 = COSD(HORAN) * COSD(DECLIN)
  244. YHOR1 = SIND(HORAN) * COSD(DECLIN)
  245. ZHOR1 = SIND(DECLIN)
  246.  
  247. XHOR = XHOR1 * COSD(90 - OLAT) - ZHOR1 * SIND(90 - OLAT)
  248. YHOR = YHOR1
  249. ZHOR = XHOR1 * SIND(90 - OLAT) + ZHOR1 * COSD(90 - OLAT)
  250. AZI  = ATAN2D( YHOR, XHOR ) + 180
  251. ALT = ASIND( ZHOR )
  252.  
  253. WRITE (*,*) 'POSITION OF SUN: '
  254. WRITE (*,'(A, T18, F12.7, A)') 'AZIMUTH ', AZI, '°'
  255. WRITE (*,'(A, T18, F12.7, A)') 'ALTITUDE: ' , ALT, '°'
  256. WRITE (*,'(A, T18, F12.7, A)') 'DISTANCE: ' , HELDIS, ' AU'
  257.  
  258. END PROGRAM PLANETARIUM
Tags: Astronomy
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement