subroutine position( 1 n,tjd,icoord,glon,glat,ht, 2 aname, 3 xequ,yequ,zequ,vxequ,vyequ,vzequ, 4 ra,dec,rah,ram,ras,decd,decm,decs, 5 xecl,yecl,zecl,vxecl,vyecl,vzecl, 6 lam,bet,lamd,lamm,lams,betd,betm,bets, 7 xalt,yalt,zalt,vxalt,vyalt,vzalt, 8 alt,azi,altd,altm,alts,azid,azim,azis, 9 dist,drdt,dsun,dsundt, & semi_maj,qper,xecen,lan,xinc,aop,xmanom,period,tjdp) implicit none save c..this routine drives the novas routines, jpl planet ephemeris, c..and jpl minor body data derive position and velocity data c..on some 15,000 objects. c..input: c..n = body number 1=sun 2=moon 3=mercury ... 11=pluto c.. 12-309 are comets, 309-15097 are asteroids c..tjd = ut julian date where data is desired c..icoord = origin of coordinate system 1=earth 2=sun 3=location on earth c..glon = observers geographic longitude in decimal degrees c.. positive west, negative east, of greenwhich c..glat = observers geographic latitude in decimal degrees c.. positive north, negative south, of equator c..ht = observers height above mean sea level in meters c..output: c..aname = name of the body number c..xequ = x coordinate of the rectangular equatorial position vector in au c..yequ = y coordinate of the rectangular equatorial position vector in au c..zequ = z coordinate of the rectangular equatorial position vector in au c..vxequ = x coordinate of the rectangular equatorial velocity vector in au c..vyequ = y coordinate of the rectangular equatorial velocity vector in au c..vzequ = z coordinate of the rectangular equatorial velocity vector in au c..ra = right ascension in decimal hours c..rah = right ascension hours c..ram = right ascension minutes c..ras = right ascension seconds c..dec = declination in decimal degrees c..decd = declination degrees c..decm = declination minutes c..decs = declination seconds c..xecl = x coordinate of the rectangular ecliptic position vector in au c..yecl = y coordinate of the rectangular ecliptic position vector in au c..zecl = z coordinate of the rectangular ecliptic position vector in au c..vxecl = x coordinate of the rectangular ecliptic velocity vector in au c..vyecl = y coordinate of the rectangular ecliptic velocity vector in au c..vzecl = z coordinate of the rectangular ecliptic velocity vector in au c..lam = ecliptic longitude in decimal degrees c..lamd = ecliptic longitude degrees c..lamm = ecliptic longitude minutes c..lams = ecliptic longitude seconds c..bet = ecliptic latitude in decimal degrees c..betd = ecliptic latitude in degrees c..betm = ecliptic latitude in minutes c..bets = ecliptic latitude in seconds c..xalt = x coordinate of the rectangular horizon position vector in au c..yalt = y coordinate of the rectangular horizon position vector in au c..zalt = z coordinate of the rectangular horizon position vector in au c..vxalt = x coordinate of the rectangular horizon velocity vector in au c..vyalt = y coordinate of the rectangular horizon velocity vector in au c..vzalt = z coordinate of the rectangular horizon velocity vector in au c..alt = altitude in decimal degrees c..altd = altitude degrees c..altm = altitude minutes c..alts = altitude seconds c..alt = azimuth in decimal degrees c..altd = azimuth in degrees c..altm = azimuth in minutes c..alts = azimuth in seconds c..dist = distance from the earth to the object in au c..drdt = rate of change in the earth-object distance in au/day c..dsun = distance from sun to object in au c..dsundt = rate of change in the sun-object distance in au/day c..semi_maj = semi major axis in au c..qper = perihelion distance in au c..xecen = eccentricity, dimensionless c..lan = longitude of the ascending node in decimal degrees c..xinc = inclination in decimal degrees c..aop = angle of perhelion in decimal degrees c..xmanom = mean anomaly in decimal degrees c..period = obital period in days c..tjdp = julian day of perihelion c..declare the pass character*40 aname integer n,icoord double precision rhour,glon,glat,ht,tjd, 1 xequ,yequ,zequ,vxequ,vyequ,vzequ, 2 ra,dec,rah,ram,ras,decd,decm,decs, 3 xecl,yecl,zecl,vxecl,vyecl,vzecl, 4 lam,bet,lamd,lamm,lams,betd,betm,bets, 5 xalt,yalt,zalt,vxalt,vyalt,vzalt, 5 alt,azi,altd,altm,alts,azid,azim,azis, 6 dist,drdt,dsun,dsundt, 7 semi_maj,qper,xecen,lan,xinc,aop,xmanom, 8 period,tjdp c..for converting between coordinate systems double precision pos(3),vel(3),pos_geo(3),vel_geo(3), 1 pos_sun(3),vel_sun(3),pos2(3),vel2(3), 2 gdlon,oblm,oblt,eqeq,dpsi,deps,ca,sa c..for getting the geodetic latitude external geodetic integer niter double precision lo,hi,tol,zbrent,gdlat,geodetic parameter (tol = 1.0d-12) double precision glat_com,ht_com common /gedc1/ glat_com,ht_com c.. integer irefr double precision x,y,zd,az,rar,decr,posvel(6) c..for the novas body positions integer nearth,nsun,m parameter (nearth = 3, nsun = 10) c..for converting between ut and et double precision tjde,secdif c..gravitational constants double precision kgauss,kearth,mu,gm parameter (kgauss = 0.01720209895d0, 1 kearth = 0.07436680d0) c.for the masses of the planets in solar masses, iau 1976 values integer ifirst,i double precision muplan(9) c..for the names and data structures include 'ephem.dek' data ifirst/0/ c..initialize the planet masses in solar masses c..jpl 1996 reference values c..each is the sum of the parent body + any moons + any rings if (ifirst .eq. 0) then ifirst = 1 muplan(1) = 1.0d0/6023600.0d0 muplan(2) = 1.0d0/408523.71d0 c muplan(3) = 1.0d0/328900.56d0 c muplan(3) = 1.0d0/332946.038d0 muplan(3) = 1.0d0/331500.0d0 muplan(4) = 1.0d0/3098708.0d0 muplan(5) = 1.0d0/1047.3486d0 muplan(6) = 1.0d0/3497.898d0 muplan(7) = 1.0d0/22902.98d0 muplan(8) = 1.0d0/19412.24d0 muplan(9) = 1.0d0/1.35e8 end if c..initialize the output aname = ' ' xequ = 0.0d0 yequ = 0.0d0 zequ = 0.0d0 vxequ = 0.0d0 vyequ = 0.0d0 vzequ = 0.0d0 ra = 0.0d0 rah = 0.0d0 ram = 0.0d0 ras = 0.0d0 dec = 0.0d0 decd = 0.0d0 decm = 0.0d0 decs = 0.0d0 xecl = 0.0d0 yecl = 0.0d0 zecl = 0.0d0 vxecl = 0.0d0 vyecl = 0.0d0 vzecl = 0.0d0 lam = 0.0d0 lamd = 0.0d0 lamm = 0.0d0 lams = 0.0d0 bet = 0.0d0 betd = 0.0d0 betm = 0.0d0 bets = 0.0d0 xalt = 0.0d0 yalt = 0.0d0 zalt = 0.0d0 vxalt = 0.0d0 vyalt = 0.0d0 vzalt = 0.0d0 alt = 0.0d0 altd = 0.0d0 altm = 0.0d0 alts = 0.0d0 azi = 0.0d0 azid = 0.0d0 azim = 0.0d0 azis = 0.0d0 dist = 0.0d0 drdt = 0.0d0 dsun = 0.0d0 dsundt = 0.0d0 semi_maj = 0.0d0 qper = 0.0d0 xecen = 0.0d0 lan = 0.0d0 xinc = 0.0d0 aop = 0.0d0 xmanom = 0.0d0 tjdp = 0.0d0 period = 0.0d0 c..convert between the input body code and the c..body code used by the novas routines if (n .eq. 1) then m = 10 else if (n .eq. 2) then m = 11 else if (n .ge. 3 .and. n .le. 11) then m = n - 2 else m = n endif c..get the ephemeris julian date call deltat(tjd,tjde,secdif) c..for geocentric, topocentric or heliocentric origins c..get the geocentric, equatorial, apparent state vector and angular position call applan (tjde,m,nearth,pos,vel,ra,dec,dist,drdt) c if (dist .eq. 0.0) return do i=1,3 pos_geo(i) = pos(i) vel_geo(i) = vel(i) enddo c..for topographic coordinates c..the conventional geocentric longitude and latitude is input, but c..the geodetic longitude and latitude is needed, so perform the conversion. c..then get the topgraphic, equatorial apparent state vector if (icoord .eq. 3) then if (glon .gt. 0.0) gdlon = 360.0d0 - glon if (glon .lt. 0.0) gdlon = abs(glon) glat_com = glat ht_com = ht lo = 0.8d0 * glat_com hi = 1.2d0 * glat_com gdlat = zbrent(geodetic,lo,hi,tol,niter) call tpplan (tjde,gdlon,gdlat,ht,pos,vel,ra,dec,dist,drdt) end if c..get the heliocentric equatorial, apparent state vector of the sun call applan2(tjd,m,nearth,pos_sun,vel_sun,x,y,dsun,dsundt) c..for heliocentric coordinates if (icoord .eq. 2) then do i=1,3 pos(i) = pos_sun(i) vel(i) = vel_sun(i) end do ra = x dec = y end if c..form the name of the planet aname = zname(n) c..set the output equatorial rectangular coordinates xequ = pos(1) yequ = pos(2) zequ = pos(3) vxequ = vel(1) vyequ = vel(2) vzequ = vel(3) c..this next section converts equatorial to ecliptic coordinates c..get the obliquity of the ecliptic, and rotate about the x-axis call etilt (tjde,oblm,oblt,eqeq,dpsi,deps) call rotax(pos,oblt,pos2) call rotax(vel,oblt,vel2) xecl = pos2(1) yecl = pos2(2) zecl = pos2(3) vxecl = vel2(1) vyecl = vel2(2) vzecl = vel2(3) call angles(pos2,lam,bet) lam = lam*15.0d0 c..get the heliocentric, ecliptic, osculatory orbital elements c..for the moon, return the geocentric ecliptic orbital elements if (n .ne. 1) then if (n .eq. 2) then c call pleph(tjde,10,13,posvel) c do i=1,3 c pos_geo(i) = posvel(i) c vel_geo(i) = posvel(i+3) c enddo call rotax(pos_geo,oblt,pos2) call rotax(vel_geo,oblt,vel2) else if (n .ne. 2) then call rotax(pos_sun,oblt,pos2) call rotax(vel_sun,oblt,vel2) end if mu = 1.0d0 if (n .ge. 3 .and. n.le. 11) mu = 1.0d0 + muplan(n-2) gm = kgauss*kgauss * mu if (n .eq. 2) gm = gm*muplan(3) call state_to_kep(gm,tjde,pos2,vel2, 1 semi_maj,qper,xecen,lan,xinc,aop,xmanom,period,tjdp) end if c..this next section converts equatorial to horizon coordinates c..note the pass is the ut julian date tjd and not the ephemeris c..julian date tjde, since the ut julian date is to be used c..when computing greenwhich sidereal time, which the call c..to requ_to_hor does. c..remember to convert azimuths from hours to degrees at the end if (icoord .eq. 3) then call requ_to_rhor(xequ,yequ,zequ,tjd,glon,glat, 1 xalt,yalt,zalt) call requ_to_rhor(vxequ,vyequ,vzequ,tjd,glon,glat, 1 vxalt,vyalt,vzalt) pos2(1) = xalt pos2(2) = yalt pos2(3) = zalt call angles(pos2,azi,alt) azi = azi*15.0d0 endif c..account for refraction c x = 0.0d0 c y = 0.0d0 c irefr = 0 c call zdaz (tjd,x,y,gdlon,glat,ht,ra,dec,irefr, c 1 zd,az,rar,decr) c c write(6,222) 90.0d0-zd,az,rar,decr c 222 format(1x,1p4e16.8) c ra = rar c dec = decr c..this next section converts all the decimal angles into hour:min:sec format c..convert the right ascension and declination call dh2h(ra,rah,ram,ras) call dh2h(dec,decd,decm,decs) c..convert the ecliptic longitude and latitude call dh2h(lam,lamd,lamm,lams) call dh2h(bet,betd,betm,bets) c..convert the alt-horizon coordinates call dh2h(alt,altd,altm,alts) call dh2h(azi,azid,azim,azis) return end c..do all the includes here. makes it easier to maintain include 'initialize_routine.f' include 'jpl_routines.f' include 'novas_routines.f' include 'fxt_routines.f' include 'marsat_routines.f' include 'lieske_routines.f' include 'satsat_routines.f' include 'gust86_routines.f' c..for silly compilers that don't have intrinsic trig functions in degrees include 'trig_degrees.f'