SUBROUTINE SOLSYS VERSION 3. THIS SUBROUTINE PROVIDES THE POSITION AND VELOCITY OF THE EARTH AT EPOCH TJD BY EVALUATING A CLOSED-FORM THEORY WITHOUT REFERENCE TO AN EXTERNAL FILE. THIS ROUTINE CAN ALSO PROVIDE THE POSITION AND VELOCITY OF THE SUN.
TJD = TDB JULIAN DATE OF DESIRED EPOCH (IN)
M = BODY IDENTIFICATION NUMBER (IN)
SET M=0 OR M=1 FOR THE SUN
SET M=2 OR M=3 FOR THE EARTH
K = ORIGIN SELECTION CODE (IN)
SET K=0 FOR ORIGIN AT SOLAR SYSTEM BARYCENTER
SET K=1 FOR ORIGIN AT CENTER OF SUN
POS = POSITION VECTOR, EQUATORIAL RECTANGULAR
COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX
OF J2000.0, COMPONENTS IN AU (OUT)
VEL = VELOCITY VECTOR, EQUATORIAL RECTANGULAR
COORDINATES, REFERRED TO MEAN EQUATOR AND EQUINOX
OF J2000.0, COMPONENTS IN AU/DAY (OUT)
IERR = ERROR INDICATOR (OUT)
IERR=0 MEANS EVERYTHING OK
IERR=1 MEANS TJD BEFORE FIRST ALLOWED DATE
IERR=2 MEANS TJD AFTER LAST ALLOWED DATE
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| double precision | :: | tjd | ||||
| integer | :: | m | ||||
| integer | :: | k | ||||
| double precision, | dimension(3) | :: | pos | |||
| double precision, | dimension(3) | :: | vel | |||
| integer | :: | ierr |
subroutine solsys (tjd,m,k,pos,vel,ierr) use novas_module, only: preces double precision tjd,pos,vel,pi,twopi,t0,obl,el,c,p,tlast, & pm,pa,pe,pj,po,pw,pl,pn, & tmass,se,ce,si,ci,sn,cn,sw,cw,p1,p2,p3,q1,q2,q3,roote,a,b, & qjd,e,mlon,ma,u,sinu,cosu,anr,pplan,vplan,f,pbary,vbary, & dfloat,dabs,dmod,dsin,dcos,dsqrt dimension pos(3), vel(3), el(21), c(13), p(3,3), & pm(4), pa(4), pe(4), pj(4), po(4), pw(4), pl(4), pn(4), & a(3,4), b(3,4), pplan(3), vplan(3), pbary(3), vbary(3) save parameter ( pi = 3.14159265358979324d0 ) parameter ( twopi = 2.d0 * pi ) parameter ( t0 = 2451545.0d0 ) parameter ( obl = 23.43927944d0 ) ! T0 = TDB JULIAN DATE OF EPOCH J2000.0 ! OBL = OBLIQUITY OF ECLIPTIC AT EPOCH J2000.0 data el, c, p / 43*0.d0 /, tlast / 0.d0 / ! ARRAYS BELOW CONTAIN MASSES AND ORBITAL ELEMENTS OF THE FOUR ! LARGEST PLANETS (SEE EXPLANATORY SUPPLEMENT (1992), P. 316) ! WITH ANGLES IN RADIANS ! THIS DATA USED FOR BARYCENTER COMPUTATIONS ONLY ! JUPITER SATURN URANUS NEPTUNE data pm / 1047.349d+0, 3497.898d+0, 22903.0d+0, 19412.2d+0 / ! data pa / 5.203363d+0, 9.537070d+0, 19.191264d+0, 30.068963d+0 / ! data pe / 0.048393d+0, 0.054151d+0, 0.047168d+0, 0.008586d+0 / ! data pj / 0.022782d+0, 0.043362d+0, 0.013437d+0, 0.030878d+0 / ! data po / 1.755036d+0, 1.984702d+0, 1.295556d+0, 2.298977d+0 / ! data pw / 0.257503d+0, 1.613242d+0, 2.983889d+0, 0.784898d+0 / ! data pl / 0.600470d+0, 0.871693d+0, 5.466933d+0, 5.321160d+0 / ! data pn / 1.450138d-3, 5.841727d-4, 2.047497d-4, 1.043891d-4 / ! if ( tlast < 1.d0 ) then ! first time computations ! mass of sun plus four inner planets tmass = 1.d0 + 5.977d-6 se = dsin ( obl * pi / 180.d0 ) ce = dcos ( obl * pi / 180.d0 ) do i = 1, 4 tmass = tmass + 1.d0 / pm(i) ! compute sine and cosine of orbital angles si = dsin ( pj(i) ) ci = dcos ( pj(i) ) sn = dsin ( po(i) ) cn = dcos ( po(i) ) sw = dsin ( pw(i) - po(i) ) cw = dcos ( pw(i) - po(i) ) ! compute p and q vectors (see brouwer & clemence (1961), ! methods of celestial mechanics, pp. 35-36.) p1 = cw * cn - sw * sn * ci p2 = ( cw * sn + sw * cn * ci ) * ce - sw * si * se p3 = ( cw * sn + sw * cn * ci ) * se + sw * si * ce q1 = -sw * cn - cw * sn * ci q2 = ( -sw * sn + cw * cn * ci ) * ce - cw * si * se q3 = ( -sw * sn + cw * cn * ci ) * se + cw * si * ce roote = dsqrt ( 1.d0 - pe(i)**2 ) a(1,i) = pa(i) * p1 a(2,i) = pa(i) * p2 a(3,i) = pa(i) * p3 b(1,i) = pa(i) * roote * q1 b(2,i) = pa(i) * roote * q2 b(3,i) = pa(i) * roote * q3 end do tlast = 1.d0 end if ierr = 0 ! VALID DATES ARE WITHIN 3 CENTURIES OF J2000, ALTHOUGH RESULTS ! DETERIORATE GRADUALLY if ( tjd < 2340000.5d0 ) ierr = 1 if ( tjd > 2560000.5d0 ) ierr = 2 if ( ierr /= 0 ) return if ( m >= 2 ) then ! FORM HELIOCENTRIC COORDINATES OF EARTH ! VELOCITIES ARE OBTAINED FROM CRUDE NUMERICAL DIFFERENTIATION do i = 1, 3 qjd = tjd + dfloat(i-2) * 0.1d0 ! SUBROUTINE SUN COMPUTES EARTH-SUN VECTOR call sun ( qjd, el, c ) call preces ( qjd, c(11), t0, pos ) p(i,1) = -pos(1) p(i,2) = -pos(2) p(i,3) = -pos(3) end do do j=1,3 pos(j) = p(2,j) vel(j) = ( p(3,j) - p(1,j) ) / 0.2d0 end do if ( k >= 1 ) return else ! FORM HELIOCENTRIC COORDINATES OF SUN do j=1,3 pos(j) = 0.d0 vel(j) = 0.d0 end do if ( k >= 1 ) return end if ! IF K=0, MOVE ORIGIN TO SOLAR SYSTEM BARYCENTER ! SOLAR SYSTEM BARYCENTER COORDINATES COMPUTED FROM KEPLERIAN ! APPROXIMATIONS OF THE COORDINATES OF THE FOUR LARGEST PLANETS if ( dabs ( tjd - tlast ) >= 1.d-6 ) then do j = 1, 3 pbary(j) = 0.d0 vbary(j) = 0.d0 end do ! THE FOLLOWING LOOP CYCLES ONCE FOR EACH OF THE FOUR LARGE PLANETS do i = 1, 4 ! COMPUTE MEAN LONGITUDE, MEAN ANOMALY, AND ECCENTRIC ANOMOLY e = pe(i) mlon = pl(i) + pn(i) * ( tjd - t0 ) ma = dmod ( mlon - pw(i), twopi ) u = ma + e * dsin ( ma ) + 0.5d0 * e * e * dsin ( 2.d0 * ma ) sinu = dsin ( u ) cosu = dcos ( u ) ! COMPUTE VELOCITY FACTOR anr = pn(i) / ( 1.d0 - e * cosu ) ! COMPUTE PLANET'S POSITION AND VELOCITY WRT EQ & EQ J2000 pplan(1) = a(1,i) * ( cosu - e ) + b(1,i) * sinu pplan(2) = a(2,i) * ( cosu - e ) + b(2,i) * sinu pplan(3) = a(3,i) * ( cosu - e ) + b(3,i) * sinu vplan(1) = anr * ( -a(1,i) * sinu + b(1,i) * cosu ) vplan(2) = anr * ( -a(2,i) * sinu + b(2,i) * cosu ) vplan(3) = anr * ( -a(3,i) * sinu + b(3,i) * cosu ) ! COMPUTE MASS FACTOR AND ADD IN TO TOTAL DISPLACEMENT f = 1.d0 / ( pm(i) * tmass ) pbary(1) = pbary(1) + pplan(1) * f pbary(2) = pbary(2) + pplan(2) * f pbary(3) = pbary(3) + pplan(3) * f vbary(1) = vbary(1) + vplan(1) * f vbary(2) = vbary(2) + vplan(2) * f vbary(3) = vbary(3) + vplan(3) * f end do tlast = tjd end if do j=1,3 pos(j) = pos(j) - pbary(j) vel(j) = vel(j) - vbary(j) end do end subroutine solsys