program riset implicit none save c..rise, transit, set times of objects c..this version uses a faster iteration method, which can be fooled c..for unusual cases, to get the rise, transit, and set times c..declare locals character*1 asign integer i,iwant,ilo,ihi,ioption c..for the rise_and_set routine integer nmax,imax,ndays parameter (nmax=1000, imax = 11) character*13 rise(nmax,imax), 1 transit(nmax,imax), 2 set(nmax,imax) character*9 adate(nmax) character*40 aname(imax) integer j,ibody,iyear,imonth,iday,ifirst double precision glon,glat,height,zone c..for picking which ephemeris to use integer which_eph common /dop/ which_eph c..popular formats 01 format( 1 1x,a,1pe11.3,t27,a,1pe11.3,/, 2 1x,a,a,i3.2,a,/) 04 format(1x,t22,10(a,' ')) 05 format(1x,t2,a,t17,a,t25,a,t32,a,t40,a,t47,a,t55,a) 06 format(1x,a13,a15,a15,a15) c..initialize the name database which_eph = 405 call init_names(0) c..get the input vector 10 write(6,*) 'a body over lots of days = 1' write(6,*) 'solar system over a few days = 2 =>' read(5,*) ioption if (ioption .ne. 1 .and. ioption .ne. 2) goto 10 if (ioption .eq. 1) then write(6,*) 'give body number =>' read(5,*) iwant end if write(6,*) 'give iyear imonth iday', 1 ' longitude latitude height timezone ndays =>' read(5,*) iyear,imonth,iday,glon,glat,height,zone,ndays c..write a header write(6,*) asign='+' if (zone .le. 0.0) asign ='-' write(6,01) 1 'longitude =',glon, 2 'latitude =',glat, 3 'all times in local time = UT ',asign,int(abs(zone)),' hours' c..for the sun, moon and eight planets ilo = 1 ihi = 11 if (ioption .eq. 1) then ilo = iwant ihi = iwant end if do ibody = ilo,ihi if (ibody .ne. 5) then call rise_and_set( 1 ibody,iyear,imonth,iday,glon,glat,height,zone,ndays, 2 aname(ibody),adate, 3 rise(1,ibody),transit(1,ibody),set(1,ibody)) end if enddo c..various writes depending on what one wants c..for one body over a period of time if (ioption .eq. 1) then write(6,06) aname(iwant) write(6,05) 'body','rise','azi','trans','alt','set','azi' do j=1,ndays do i=ilo,ihi if (i .ne. 5) then write(6,06) adate(j)//' ', 1 rise(j,i),transit(j,i),set(j,i) end if enddo enddo write(6,*) c..for all bodies over a limited time span else if (ioption .eq. 2) then do j=1,ndays ifirst = 0 do i=1,11 if (i .ne. 5) then if (ifirst .eq. 0) then ifirst = 1 write(6,04) adate(j) write(6,05) 'body','rise','azi','trans','alt','set','azi' endif write(6,06) aname(i),rise(j,i),transit(j,i),set(j,i) end if enddo write(6,*) enddo end if c..back for another goto 10 end subroutine rise_and_set( 1 ibody,iyear,imonth,iday,glon,glat,height,zone,ndays, 2 aname,adate,rise_string,transit_string,set_string) implicit none save c..this routine computes the local rise and set times for c..any body supported by routine position. c..input: c..ibody = body number 1=sun 2=moon 3=mercury ... 11=pluto c.. 12-309 are comets, 309-15097 are asteroids c..iyear = year c..imonth = month c..iday = day c..glon = observers geographic longitude in decimal degrees c..glat = observers geographic latitude in decimal degrees c..height = observers height above mean sea level in meters c..zone = local time zone (+ east, - west on greenwhich) c.. eastern standard time is -5.0 c.. central standard time is -6.0 c.. mountain standard time is -7.0 c.. pacific standard time is -8.0 c.. add one hour to daylight savings time is on. c..ndays = number of days to output, dimensions of output arrays c..output: c..anem = name of the object, character*40 string c..adate = date of event c..rise_string = local rise time of object c..transit_string = local transit (highest altitude) time of object c..set_string = local set time of object c..declare the pass integer ndays character*13 rise_string(ndays), 1 transit_string(ndays), 2 set_string(ndays) character*9 adate(ndays) character*40 aname integer ibody,iyear,imonth,iday double precision glon,glat,height,zone c..for getting the positional data integer icoord parameter (icoord = 3) double precision 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, 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..local variables logical ok,circum,never character*3 amon(12) integer jj,i,nz,iy,im,id,icase double precision tjd0,zoneday,tjd1, 1 sin_h0,sid,ra_0h,dec_0h,ra_24h,dec_24h, 2 tjdh,tjdl,gst,lst_0h,zt0,sda,lst,dtau,dzt,zt, 3 utrise,utriseh,utrisem,utrises, 4 utset,utseth,utsetm,utsets,xh, 5 utcul,utculh,utculm,utculs,mod24,tol parameter (tol = 1.0d-3) c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..data statements data amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..format statements for the output 03 format(i2.2,':',i2.2,' (',0p1f5.1,')') 04 format('------------') 05 format(a) 123 format(i2.2,a,i4.4) 222 format(1x,1p6e14.6) c..get the julian date at 0 hour ut c..then correct the julian date for the time zone difference call juldat(iyear,imonth,iday,0.0d0,tjd0) zoneday = zone/24.0d0 tjd0 = tjd0 - zoneday c..solar to sidereal time conversion sid = 1.002737909350795D0 c..altitude for the sun c..parallax is already taken into account in the novas routines. c..all we need is the corrections for the radius of the sun and refraction c..here i use average values if (ibody .eq. 1) then sin_h0 = sind( -0.5d0 * 0.533d0 - 34.0d0/60.0d0) c..altitude for the moon c..parallax is already taken into account in the novas routines. c..all we need is the corrections for the radius of the sun and refraction c..here i use average values else if (ibody .eq. 2) then sin_h0 = sind(-0.5d0 * 0.5181d0 - 34.0d0/60.0d0) c..altitude for planets and stars c..parallax is already taken into account in the novas routines. c..all we need here is refraction else sin_h0 = sind(-34.0d0/60.0d0) end if c..for ndays days do jj=1,ndays tjd1 = tjd0 + dble(jj-1) c..form the output calander date string call caldat(tjd1+zoneday,iy,im,id,xh) write(adate(jj),123) id,amon(im),iy c..get the position of the body at 0h local time tjd = tjd1 call position( 1 ibody,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) ra_0h = 15.0d0 * ra dec_0h = dec c..get the local sidereal time in hours at 0h local time tjdh = int(tjd) tjdl = tjd - tjdh call sidtim(tjdh,tjdl,1,gst) lst_0h = gst - glon/15.0d0 if (lst_0h .lt. 24.0) lst_0h = lst_0h + 24.0d0 c..get the position of the body at 24h local time tjd = tjd1 + 1.0d0 call position( 1 ibody,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) ra_24h = 15.0d0 * ra dec_24h = dec c..generate continuous right ascension in case of jumps between 0h and 24h if (ra_0h - ra_24h .gt. 180.0) ra_24h = ra_24h + 360.0d0 if (ra_0h - ra_24h .lt. -180.0) ra_0h = ra_0h + 360.0d0 c..for each case, rising transit, setting do icase=1,3 c..initialize for the iteration zt0 = 12.0d0 ok = .false. circum = .false. never = .false. c..start the iteration loop do i=1,10 c..linear interpolation of planetary position c ra = ra_0h + (zt0/24.0) * (ra_24h - ra_0h) c dec = dec_0h + (zt0/24.0) * (dec_24h - dec_0h) c write(6,222) ra,dec c..no interpolation, one position call per iteration tjd = tjd1 + zt0/24.0d0 call position( 1 ibody,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) ra = ra * 15.0d0 c write(6,222) ra,dec c..possible altitude angle sda = (sin_h0 - sind(dec)*sind(glat))/(cosd(dec)*cosd(glat)) c..check if its circumpolar or always above the horizon if (abs(sda) .lt. 1.0) then sda = acosd(sda) ok = .true. circum = .false. never = .false. else ok = .false. if (glat .ge. 0.0) then if (dec .ge. 90.0-glat) then circum = .true. never = .false. else circum = .false. never = .true. end if else if (dec .lt. -90.0-glat) then circum = .true. never = .false. else circum = .false. never = .true. end if end if endif c..if its ok, get sidereal time at universal time zt0 if (ok) then lst = lst_0h + zt0*sid if (icase .eq. 1) then dtau = (lst - ra/15.0d0) + sda/15.0d0 else if (icase .eq. 2) then dtau = (lst - ra/15.0d0) else if (icase .eq. 3) then dtau = (lst - ra/15.0d0) - sda/15.0d0 end if dzt = (mod24(dtau+12.0d0) - 12.0d0)/sid zt = zt0 - dzt zt0 = zt end if c..check convergence of iteration loop if (abs(dzt) .le. tol .or. .not.ok) goto 100 enddo write(6,*) write(6,*) 'did not converge',icase write(6,*) stop 'error in routine rise_and_set' c..start the output 100 continue c..form the rise time string if (icase .eq. 1) then if (ok) then if (zt .gt. 24.0d0) then write(rise_string(jj),04) else call dh2h(zt,utriseh,utrisem,utrises) write(rise_string(jj),03) int(utriseh), 1 nint(utrisem+utrises/60.0d0),azi end if else if (circum) then write(rise_string(jj),05) 'circum' else if (never) then write(rise_string(jj),05) 'below' endif c..form the transit time string else if (icase .eq. 2) then if (ok) then call dh2h(zt,utculh,utculm,utculs) write(transit_string(jj),03) int(utculh), 1 nint(utculm+utculs/60.0d0),alt else if (circum) then write(transit_string(jj),05) 'circum' else if (never) then write(transit_string(jj),05) 'below' endif c..form the setting time string else if (icase .eq. 3) then if (ok) then if (zt .lt. 0.0) then write(set_string(jj),04) else call dh2h(zt,utseth,utsetm,utsets) write(set_string(jj),03) int(utseth), 1 nint(utsetm+utsets/60.0d0),azi end if else if (circum) then write(set_string(jj),05) 'circum' else if (never) then write(set_string(jj),05) 'below' endif endif c..go back for another case of rising, transit or setting enddo c..end of day loop enddo return end include 'position_routine.f' include 'glue.f'