Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- MODULE ORBITALTRIG
- REAL X, RADEG
- PARAMETER(RADEG=57.2957795130823)
- CONTAINS
- FUNCTION SIND(X)
- SIND = SIN( X * (1.0/RADEG) )
- END FUNCTION
- FUNCTION ATAND(X)
- ATAND = RADEG * ATAN(X)
- END FUNCTION
- FUNCTION ATAN2D(Y,X)
- ATAN2D = RADEG * ATAN2(Y,X)
- END FUNCTION
- FUNCTION REV(X)
- REV = X - AINT(X/360.0)*360.0
- IF (REV.LT.0.0) REV = REV + 360.0
- END FUNCTION
- FUNCTION CBRT(X)
- IF (X.GE.0.0) THEN
- CBRT = X ** (1.0/3.0)
- ELSE
- CBRT = -((-X)**(1.0/3.0))
- ENDIF
- END FUNCTION
- FUNCTION COSD(X)
- COSD = COS( X * (1.0/RADEG))
- END FUNCTION
- FUNCTION TAND(X)
- TAND = TAN(X * (1.0/RADEG))
- END FUNCTION
- FUNCTION ASIND(X)
- ASIND = RADEG * ASIN (X)
- END FUNCTION
- FUNCTION ACOSD(X)
- ACOSD = RADEG * ACOS (X)
- END FUNCTION
- END MODULE ORBITALTRIG
- PROGRAM PLANETARIUM
- USE ORBITALTRIG
- IMPLICIT NONE
- REAL, PARAMETER :: PI = 3.14159263585
- !TIME PARAMETERS
- INTEGER YEAR, MONTH, DATE !USER INPUT YEAR, MONTH, AND DAY
- 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)
- REAL EPOCH !TIME SINCE MIDNIGHT JAN 1ST 2000 (UNIT=DAYS)
- REAL UT !GREENWICH MEAN TIME
- REAL GMST0 !SIDERIAL TIME AT GREENWICH MERIDIAN AT 00:00
- REAL SIDTIME !LOCAL SIDERIAL TIME
- !OBSERVER TERRESTRIAL COORDINATES
- REAL OLON !OBSERVER LONGITUDE
- REAL OLAT !OBSERVER LATITUDE'
- !THESE ARE FOR CONVERTING DATE_AND_TIME OUTPUT TO FEED EPOCH
- INTEGER V(8), YYYY, MM, DD, HH, MINS, SS, MSS, GMT, DAYS
- EQUIVALENCE (V(1),YYYY), (V(2),MM), (V(3),DD), &
- (V(4),GMT), (V(5),HH), (V(6),MINS), &
- (V(7),SS), (V(8),MSS)
- !RECTANGULAR COORDINATES FOR THE ELIPTIC PLANE
- REAL XECLIP
- REAL YECLIP
- REAL ZECLIP
- !RECTANGULAR COORDINATES FOR THE EQUATORIAL PLANE
- REAL XEQUAT
- REAL YEQUAT
- REAL ZEQUAT
- !RECTANGULAR COORDINATES FOR OBSERVER
- REAL XHOR
- REAL YHOR
- REAL ZHOR
- REAL AZI
- REAL ALT
- !PLACEHOLDERS BETWEEN AXIS ROTATIONS
- REAL XHOR1
- REAL YHOR1
- REAL ZHOR1
- !POLAR COORDINATES
- REAL RIASC !RIGHT ASCENSION (UNIT=DEGREES)
- REAL DECLIN !DECLINATION (UNIT=DEGREES)
- !DERIVED PARAMETERS:
- REAL PERDIS !PERIAPSIS DISTANCE (UNIT=AU)
- REAL APODIS !APOAPSIS DISTANCE (UNIT=AU)
- REAL ORBDAY !ORBITAL PERIOD (UNIT=DAYS)
- REAL DAYMOT !DAILY MOTION (UNIT=DEGREES PER DAY)
- REAL LASPER !TIME SINCE LAST PERIAPSIS (UNIT=DAYS)
- REAL TRUANO !TRUE ANOMALY (UNIT=DEGREES)
- REAL MENANO !MEAN ANOMALY (UNIT=DEGREES)
- REAL MENLON !MEAN LONGITUDE (UNIT=DEGREES)
- REAL OBLECL !OBLEQUITY OF THE ECLIPTIC
- REAL ECCANO !ECCENTRIC ANOMALY
- CHARACTER :: OBJTYP !OBJECT TYPE (PLANET, STAR, COMET, ASTEROID, MOON)
- CHARACTER STAR, PLANET, ASTEROID, COMET, MOON
- REAL :: HELDIS !HELIOCENTRIC DISTANCE (UNIT=AU)
- REAL :: SOLMAS !SOLAR MASSES (UNITLESS)
- REAL HORAN ! HOUR ANGLE
- !KEPLERIAN ORBITAL PARAMETERS:
- REAL :: ARGPER !ARGUMENT OF PERIAPSIS (BIG OMEGA) (UNIT=DEGREES)
- REAL :: SMAJAX !SEMI-MAJOR AXIS (SMALL A) (UNIT=AU)
- REAL :: ECCEN !ECCENTRICITY (SMALL E) (UNITLESS)
- REAL :: INCLI !INCLINATION (UNIT=DEGREES)
- REAL :: RAOAN !RIGHT ASCENSION OF ASCENDING NODE (UNIT=DEGREES)
- REAL :: PERIDAY !TIME AT PERIAPSIS (BIG T) (UNIT=DAYS)
- CHARACTER (LEN=15) OBJNAME
- CHARACTER (LEN=15) LOCNAME
- CHARACTER YESNO
- WRITE (*,*) 'USE CURRENT DATE AND TIME? <Y/N>'
- READ (*,*) YESNO
- IF (YESNO.EQ.'N') THEN
- WRITE (*,*) 'ENTER DATE <YYYY MM DD>'
- READ (*,*) YEAR,MONTH,DATE
- EPOCH = 367*YEAR-7*(YEAR+(MONTH+9)/12)/4-3*((YEAR+(MONTH-9)/7)/100+1)/4+275*MONTH/9+DATE-730515
- WRITE (*,*) 'ENTER LOCAL TIME <HH MM SS>'
- READ (*,*) HOUR, MINUT, SECON
- WRITE (*,*) 'ENTER UTC OFFSET, PST=-8, EST=-5, GMT=0, EET=2'
- READ (*,*) ZONE
- EPOCH = EPOCH + ((HOUR+ZONE) / 24.0) + (MINUT / 1440) + (SECON / 86400.0)
- WRITE (*,'(A, TR3, I2, A, I2, A, I5)') 'THE JULIAN DAY NUMBER FOR: ', DATE, '/', MONTH, '/', YEAR
- WRITE (*,"(T28, A, TR1, I2.1, A, I2.1, A, I2.2, A, I4.2, A)") 'AT', HOUR, ':', MINUT, ':', SECON,' UTC',ZONE,' IS:'
- WRITE (*,"(T28, F15.5)") EPOCH
- ELSE
- CALL DATE_AND_TIME (VALUES=V)
- DAYS = 367*YYYY - 7 * (YYYY + (MM+9)/12 ) / 4 - 3 * ( (YYYY + (MM-9)/7 ) / 100 + 1 ) / 4 + 275*MM/9 + DD - 730515
- EPOCH = DAYS + ((HH + GMT)/ 24.0) + (MINS / 1440) + (SS / 86400.0) + (MSS / 86400000.0) + GMT/24.0
- WRITE (*,'(A, TR3, I2.2, A, I2.2, A, I4.4)') 'THE JULIAN DAY NUMBER FOR: ', DD, '/', MM, '/', YYYY
- WRITE (*,"(T28, A, TR1, I2.2, A, I2.2, A, I2.2, A)") 'AT', HH, ':', MINS, ':', SS,' UTC+0 IS:'
- WRITE (*,"(T28, F15.10)") EPOCH
- END IF
- UT = ZONE + HOUR + (MINUT/60) + (SECON/3600)
- WRITE (*,'(A,T25,A)') 'CHOOSE A LOCATION:', 'OPTION'
- WRITE (*,'(A,T25,A)') 'SAN FRANCISCO', '{1}'
- WRITE (*,'(A,T25,A)') 'NEW YORK', '{2}'
- WRITE (*,'(A,T25,A)') 'LONDON', '{3}'
- WRITE (*,'(A,T25,A)') 'MOSCOW', '{4}'
- WRITE (*,'(A,T25,A)') 'CAIRO', '{5}'
- WRITE (*,'(A,T25,A)') 'NEW DELHI', '{6}'
- WRITE (*,'(A,T25,A)') 'HONG KONG', '{7}'
- WRITE (*,'(A,T25,A)') 'TOKYO', '{8}'
- WRITE (*,'(A,T25,A)') 'HONOLULU', '{9}'
- WRITE (*,'(A,T25,A)') 'SYDNEY', '{10}'
- WRITE (*,'(A,T25,A)') 'CAPE TOWN', '{11}'
- WRITE (*,'(A,T25,A)') 'BUENOSAIRES', '{12}'
- WRITE (*,'(A,T25,A)') '<ENTER COORDINATES>', '{13}'
- READ (*,*) LOCNAME
- IF (LOCNAME.EQ.'1') THEN
- OLAT = 37.7775
- OLON = -122.416389
- ELSE IF (LOCNAME.EQ.'2') THEN
- OLAT = 40.712778
- OLON = -74.006111
- ELSE IF (LOCNAME.EQ.'3') THEN
- OLAT = 51.507222
- OLON = -0.1275
- ELSE IF (LOCNAME.EQ.'4') THEN
- OLAT = 55.755833
- OLON = 37.617222
- ELSE IF (LOCNAME.EQ.'5') THEN
- OLAT = 30.044444
- OLON = 31.235833
- ELSE IF (LOCNAME.EQ.'6') THEN
- OLAT = 28.613895
- OLON = 77.209006
- ELSE IF (LOCNAME.EQ.'7') THEN
- OLAT = 22.280792
- OLON = 114.165578
- ELSE IF (LOCNAME.EQ.'8') THEN
- OLAT = 35.689722
- OLON = 139.692222
- ELSE IF (LOCNAME.EQ.'9') THEN
- OLAT = 21.306944
- OLON = -157.858333
- ELSE IF (LOCNAME.EQ.'10') THEN
- OLAT = -33.867778
- OLON = 151.21
- ELSE IF (LOCNAME.EQ.'11') THEN
- OLAT = -33.925278
- OLON = 18.423889
- ELSE IF (LOCNAME.EQ.'12') THEN
- OLAT = -34.603333
- OLON = -58.381667
- ELSE IF (LOCNAME.EQ.'13') THEN
- WRITE (*,*) 'ENTER LATITUDE <DDD MM SS>'
- READ (*,*) HOUR, MINUT, SECON
- OLAT = HOUR + (MINUT/60) + (SECON/3600)
- WRITE (*,*) 'ENTER LONGITUDE <DDD MM SS>'
- READ (*,*) HOUR, MINUT, SECON
- OLON = HOUR + (MINUT/60) + (SECON/3600)
- ELSE
- WRITE (*,*) 'INVALID ENTRY'
- WRITE (*,*) 'ENTER LATITUDE <DDD MM SS>'
- READ (*,*) HOUR, MINUT, SECON
- OLAT = HOUR + (MINUT/60) + (SECON/3600)
- WRITE (*,*) 'ENTER LONGITUDE <DDD MM SS>'
- READ (*,*) HOUR, MINUT, SECON
- OLON = HOUR + (MINUT/60) + (SECON/3600)
- END IF
- WRITE (*,*) LOCNAME
- WRITE (*,*) 'LON: ', OLON
- WRITE (*,*) 'LAT: ', OLAT
- ARGPER = 282.9404 + (4.70935E-5 * EPOCH) ! ARGUMENT OF PERIHELION
- SMAJAX = 1.000000 ! MEAN DISTANCE, I.E. SEMIMAJOR AXIS
- ECCEN = 0.016709 - (1.151E-9 * EPOCH) ! ECCENTRICITY
- MENANO = REV(356.0470 + (0.9856002585 * EPOCH)) ! MEAN ANOMALY
- OBLECL = 23.4393 - (0.0000003563 * EPOCH) ! OBLEQUITY OF ECLIPTIC
- MENLON = ARGPER + MENANO
- ECCANO = MENANO + (180/PI) * ECCEN * SIND(MENANO) * (1 + ECCEN * COSD(MENANO))
- XECLIP = COSD(ECCANO) - ECCEN
- YECLIP = SIND(ECCANO) * SQRT(1 - ECCEN*ECCEN)
- HELDIS = SQRT(XECLIP*XECLIP + YECLIP*YECLIP)
- TRUANO = ATAN2D(YECLIP, XECLIP)
- RIASC = TRUANO + ARGPER
- XECLIP = HELDIS * COSD(RIASC)
- YECLIP = HELDIS * SIND(RIASC)
- ZECLIP = 0.0
- XEQUAT = XECLIP
- YEQUAT = YECLIP * COSD(OBLECL) - ZECLIP * SIND(OBLECL)
- ZEQUAT = YECLIP * SIND(OBLECL) + ZECLIP * COSD(OBLECL)
- HELDIS = SQRT(XEQUAT*XEQUAT + YEQUAT*YEQUAT + ZEQUAT*ZEQUAT )
- RIASC = ATAN2D( YEQUAT, XEQUAT )
- DECLIN = ATAN2D(ZEQUAT, SQRT( XEQUAT*XEQUAT + YEQUAT*YEQUAT))
- GMST0 = (MENLON/15.0 + 12.0) - AINT((MENLON/15.0 + 12.0)/24)*24
- IF (GMST0.LT.0.0) GMST0 = GMST0 + 24
- SIDTIME = GMST0 + UT + OLON/15.0
- HORAN = (SIDTIME*15) - RIASC
- XHOR1 = COSD(HORAN) * COSD(DECLIN)
- YHOR1 = SIND(HORAN) * COSD(DECLIN)
- ZHOR1 = SIND(DECLIN)
- XHOR = XHOR1 * COSD(90 - OLAT) - ZHOR1 * SIND(90 - OLAT)
- YHOR = YHOR1
- ZHOR = XHOR1 * SIND(90 - OLAT) + ZHOR1 * COSD(90 - OLAT)
- AZI = ATAN2D( YHOR, XHOR ) + 180
- ALT = ASIND( ZHOR )
- WRITE (*,*) 'POSITION OF SUN: '
- WRITE (*,'(A, T18, F12.7, A)') 'AZIMUTH ', AZI, '°'
- WRITE (*,'(A, T18, F12.7, A)') 'ALTITUDE: ' , ALT, '°'
- WRITE (*,'(A, T18, F12.7, A)') 'DISTANCE: ' , HELDIS, ' AU'
- END PROGRAM PLANETARIUM
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement