program within implicit none save c..two objects with a specified angular diatance c..declare character*40 aname1,aname2 character*3 amon(12) integer nmax parameter (nmax=1000) integer ibody1,ibody2, 1 iyear_beg,imonth_beg,iday_beg, 1 iyear_end,imonth_end,iday_end, 3 nfound,i, 2 iy,im,id double precision tstep,angsep,tjd(nmax),angfind(nmax),xh c..for picking which ephemeris to use integer which_eph common /dop/ which_eph c..data statements data amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..initialize the name database which_eph = 406 call init_names(0) 100 write(6,*) write(6,*) 'give the two body numbers =>' read(5,*) ibody1,ibody2 write(6,*) write(6,*) 'give the starting year month and day =>' read(5,*) iyear_beg,imonth_beg,iday_beg write(6,*) write(6,*) 'give the ending year month and day =>' read(5,*) iyear_end,imonth_end,iday_end write(6,*) write(6,*) 'give the tstep to search over' write(6,*) '1.0=1day 2.0d=2 days 0.5 = half days and so on =>' read(5,*) tstep write(6,*) write(6,*) 'give the angular seperation in decimal degrees =>' read(5,*) angsep c..get the dates on which these two objects are within the c..specified angular seperation call find_angsep( 1 ibody1,ibody2, 2 iyear_beg,imonth_beg,iday_beg, 3 iyear_end,imonth_end,iday_end, 4 tstep,angsep, 5 aname1,aname2, 6 tjd,angfind,nmax,nfound) c..report what was found if (nfound .eq. 0) then write(6,*) write(6,*) 'no dates were found for these two objects' write(6,*) 'within the specified limits and angular seperation' write(6,*) else write(6,*) do i=1,nfound call caldat(tjd(i),iy,im,id,xh) write(6,01) i,id,amon(im),iy,xh,aname1,aname2,angfind(i) 01 format(1x,i4,' on ',i2.2,a,i4.4,' +',0pf5.2,' hours UT ', 1 a10,a10,'are within',1pe11.3,' degrees') enddo end if c..go back for another search goto 100 end subroutine find_angsep( 1 ibody1,ibody2, 2 iyear_beg,imonth_beg,iday_beg, 3 iyear_end,imonth_end,iday_end, 4 tstep,angsep, 5 aname1,aname2, 6 tjd_found,angsep_found,nmax,nfound) implicit none save c..this routines returns the dates that two objects are within a c..specified angular seperation. c..input: c..ibody1 = body number of first object c..ibody2 = body number of second object c..iyear_beg = starting year of the window in which to search c..imonth_beg = starting month of the window in which to search c..iday_beg = starting day of the window in which to search c..tstep = time step in days to step through in search c..angsep = record the date if the angular seperation is at least angsep c..output: c..aname1 = name of the first object c..aname2 = name of the second object c..tjd_found = julian dates when seperation is at least angsep c..angsep_found = angular seperation on tjd_found c..nmax = maximum number of dates and seperations to record c..nfound = the number of dates and seperations found within the window c..declare the pass character*40 aname1,aname2 integer ibody1,ibody2, 1 iyear_beg,imonth_beg,iday_beg, 2 iyear_end,imonth_end,iday_end, 3 nmax,nfound double precision tstep,angsep,tjd_found(nmax), 1 angsep_found(nmax) c..for getting the positional data character*40 aname integer icoord double precision tjd,glon,glat,height, 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, 6 alt,azi,altd,altm,alts,azid,azim,azis, 7 dist,drdt,dsun,dsundt, 8 xmanom,semi_maj,qper,xecen,lan,xinc,aop, 9 period,tjdp c..assume an geocentric origin parameter (icoord = 1, 1 glon = 0.0d0, 2 glat = 0.0d0, 3 height = 0.0d0) c..local variables integer nstep,i double precision tjd_beg,tjd_end,ra1,dec1,ra2,dec2,xx,yy, 1 delra,havra,deldec,havdec,xh,s1,c1 c..initialize nfound = 0 c..get the julian dates of the start and end times call juldat(iyear_beg,imonth_beg,iday_beg,0.0d0,tjd_beg) call juldat(iyear_end,imonth_end,iday_end,0.0d0,tjd_end) c..how many steps to take depends on how fine to do the search nstep = (tjd_end - tjd_beg)/tstep + 1.0d0 c..for every time step in the search window do i=1,nstep tjd = tjd_beg + float(i-1)*tstep c..get the position of the first body call position( 1 ibody1,tjd,icoord,glon,glat,height, 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) c..convert right ascension to degrees aname1 = aname ra1 = ra * 15.0d0 dec1 = dec c..get the position of the second body call position( 1 ibody2,tjd,icoord,glon,glat,height, 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) c..convert right ascension to degrees aname2 = aname ra2 = ra * 15.0d0 dec2 = dec c..use this formulation for the angular seperation c..to avoid roundoff and truncation errors near 0 or 180 degrees delra = ra1 - ra2 havra = ( sind(0.5d0*delra) )**2 deldec = dec1 - dec2 havdec = ( sind(0.5d0*deldec) )**2 xx = havdec + cosd(dec1)*cosd(dec2)*havra s1 = sqrt(xx) c1 = sqrt(1.0d0 - xx) yy = 2.0d0*atan2d(s1,c1) c yy = 2.0d0*asind(sqrt(xx)) c..and see if we are within the users specified angular seperation if (yy .le. angsep) then c..if the number found is less than the storage allocated, the store it if (nfound+1 .le. nmax) then nfound = nfound + 1 tjd_found(nfound) = tjd angsep_found(nfound) = yy c..otherwise abort the search and return else return end if end if c..back for another julian date end do return end include 'position_routine.f'