| 1 |
SUBROUTINE ORBIT (OBLIQ,ECCN,OMEGT,DAY,SDIST,SIND,COSD,LAMBDA) 8201. |
| 2 |
C**** 8202. |
| 3 |
C**** ORBIT receives the orbital parameters and time of year, and 8203. |
| 4 |
C**** returns the distance from the sun and its declination angle. 8204. |
| 5 |
C**** The reference for the following caculations is: V.M.Blanco 8205. |
| 6 |
C**** and S.W.McCuskey, 1961, "Basic Physics of the Solar System", 8206. |
| 7 |
C**** pages 135 - 151. 8207. |
| 8 |
C**** 8208. |
| 9 |
C**** Program authors: Gary L. Russell and Robert J. Suozzo, 12/13/85 8209. |
| 10 |
C**** 8210. |
| 11 |
C**** All computations are in double-precision; 8211. |
| 12 |
C**** but the arguments are single-precision. 8212. |
| 13 |
C**** Input: OBLIQ = latitude of tropics in degrees 8213. |
| 14 |
C**** ECCEN = eccentricity of the orbital ellipse 8214. |
| 15 |
C**** OMEGT = angle from vernal equinox to perihelion in degrees 8215. |
| 16 |
C**** DAY = day of the year in days; 0 = Jan 1, hour 0 8216. |
| 17 |
C**** 8217. |
| 18 |
C**** Constants: EDAYPY = Earth days per year = 365 8218. |
| 19 |
C**** VERQNX = occurence of vernal equinox = day 79 = Mar 21 8219. |
| 20 |
C**** 8220. |
| 21 |
C**** Intermediate quantities: 8221. |
| 22 |
C**** PERIHE = perihelion during the year in temporal radians 8222. |
| 23 |
C**** MA = mean anomaly in temporal radians = 2J DAY/365 - PERIHE8223. |
| 24 |
C**** EA = eccentric anomaly in radians 8224. |
| 25 |
C**** TA = true anomaly in radians 8225. |
| 26 |
C**** BSEMI = semi minor axis in units of the semi major axis 8226. |
| 27 |
C**** GREENW = longitude of Greenwich in the Earth's reference frame 8227. |
| 28 |
C**** 8228. |
| 29 |
C**** Output: DIST = distance to the sun in units of the semi major axis8229. |
| 30 |
C**** SDIST = square of DIST 8229.5 |
| 31 |
C**** SIND = sine of the declination angle 8230. |
| 32 |
C**** COSD = cosine of the declination angle 8231. |
| 33 |
C**** LAMBDA = sun longitude in Earth's rotating reference frame 8232. |
| 34 |
C**** 8233. |
| 35 |
IMPLICIT REAL*8 (A-H,O-Z) 8234. |
| 36 |
REAL*8 MA 8235. |
| 37 |
C REAL*4 SIND,COSD,SDIST,LAMBDA,OBLIQ,ECCN,OMEGT,DAY 8236. |
| 38 |
C**** 8237. |
| 39 |
PI = 3.14159265358979D0 8238. |
| 40 |
EDAYPY = 365. 8239. |
| 41 |
VERQNX = 79. 8240. |
| 42 |
OMEGA=OMEGT*(PI/180.D0) 8241. |
| 43 |
DOBLIQ=OBLIQ*(PI/180.D0) 8242. |
| 44 |
ECCEN=ECCN 8243. |
| 45 |
C**** 8244. |
| 46 |
C**** Determine time of perihelion using Kepler's equation: 8245. |
| 47 |
C**** PERIHE-VERQNX = OMEGA - ECCEN sin(OMEGA) 8246. |
| 48 |
C**** 8247. |
| 49 |
PERIHE = OMEGA-ECCEN*SIN(OMEGA)+VERQNX*2.*PI/365. 8248. |
| 50 |
C PERIHE = DMOD(PERIHE,2.*PI) 8249. |
| 51 |
MA = 2.*PI*DAY/365.-PERIHE 8250. |
| 52 |
MA = DMOD(MA,2.*PI) 8251. |
| 53 |
C**** 8252. |
| 54 |
C**** Numerically solve Kepler's equation: MA = EA - ECCEN sin(EA) 8253. |
| 55 |
C**** 8254. |
| 56 |
EA = MA+ECCEN*(SIN(MA)+ECCEN*SIN(2.*MA)/2.) 8255. |
| 57 |
110 DEA = (MA-EA+ECCEN*SIN(MA))/(1.-ECCEN*COS(EA)) 8256. |
| 58 |
EA = EA+DEA 8257. |
| 59 |
IF (DABS(DEA).GT.1.D-8) GO TO 110 8258. |
| 60 |
C**** 8259. |
| 61 |
C**** Calculate the distance to the sun and the true anomaly 8260. |
| 62 |
C**** 8261. |
| 63 |
BSEMI = DSQRT(1.-ECCEN*ECCEN) 8262. |
| 64 |
COSEA = COS(EA) 8263. |
| 65 |
SINEA = SIN(EA) 8264. |
| 66 |
SDIST = (1.-ECCEN*COSEA)*(1.-ECCEN*COSEA) 8265. |
| 67 |
TA = DATAN2(SINEA*BSEMI,COSEA-ECCEN) 8266. |
| 68 |
C**** 8267. |
| 69 |
C**** Change the reference frame to be the Earth's equatorial plane 8268. |
| 70 |
C**** with the Earth at the center and the positive x axis parallel to 8269. |
| 71 |
C**** the ray from the sun to the Earth were it at vernal equinox. 8270. |
| 72 |
C**** The distance from the current Earth to that ray (or x axis) is: 8271. |
| 73 |
C**** DIST sin(TA+OMEGA). The sun is located at: 8272. |
| 74 |
C**** 8273. |
| 75 |
C**** SUN = (-DIST cos(TA+OMEGA), 8274. |
| 76 |
C**** -DIST sin(TA+OMEGA) cos(OBLIQ), 8275. |
| 77 |
C**** DIST sin(TA+OMEGA) sin(OBLIQ)) 8276. |
| 78 |
C**** SIND = sin(TA+OMEGA) sin(OBLIQ) 8277. |
| 79 |
C**** COSD = sqrt(1-SIND**2) 8278. |
| 80 |
C**** LAMBDA = atan[tan(TA+OMEGA) cos(OBLIQ)] - GREENW 8279. |
| 81 |
C**** GREENW = 2*3.14159 DAY (EDAYPY-1)/EDAYPY 8280. |
| 82 |
C**** 8281. |
| 83 |
SINDD = SIN(TA+OMEGA)*SIN(DOBLIQ) 8282. |
| 84 |
COSD = DSQRT(1.-SINDD*SINDD) 8283. |
| 85 |
SIND = SINDD 8284. |
| 86 |
C GREENW = 2.*PI*(DAY-VERQNX)*(EDAYPY+1.)/EDAYPY 8285. |
| 87 |
C SUNX = -COS(TA+OMEGA) 8286. |
| 88 |
C SUNY = -SIN(TA+OMEGA)*COS(DOBLIQ) 8287. |
| 89 |
C LAMBDA = DATAN2(SUNY,SUNX)-GREENW 8288. |
| 90 |
C LAMBDA = DMOD(LAMBDA,2.*PI) 8289. |
| 91 |
C**** 8290. |
| 92 |
RETURN 8291. |
| 93 |
END 8292. |
Parent Directory
| 
Revision Log
| 
Revision Graph