g program riset implicit none save c..rise, transit, set times of objects c..this version uses a bombproof, but slower, search method c..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 rise,set,above character*3 amon(12) integer jj,i,nz,iy,im,id double precision tjd0,zoneday,tjd1,tstep,xhour,yhour, 1 sin_h0,yminus,y0,yplus, 2 xe,ye,zero1,zero2, 3 utrise,utriseh,utrisem,utrises, 4 utset,utseth,utsetm,utsets,xh, 5 altmax,utcul,utculh,utculm,utculs, 6 azirise,aziset parameter (tstep = 0.1d0) 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) 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..for ndays days do jj=1,ndays tjd1 = tjd0 + dble(jj-1) c..initialize, tstep sets the number of hours between time points utrise = 0.0d0 utset = 0.0d0 xhour = tstep tjd = tjd1 c..get the position of the body 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) 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..still initializing yminus = sind(alt) - sin_h0 above = .false. if (yminus .gt. 0.0) above = .true. rise = .false. set = .false. altmax = -100.0d0 utcul = -1.0d0 azirise = -100.0d0 aziset = -100.0d0 c..look in tstep increments do i=1,100000 yhour = xhour tjd = tjd1 + yhour/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) y0 = sind(alt) - sin_h0 c..if this is the largest altitude, capture the transit if (alt .gt. altmax) then altmax = alt utcul = yhour end if yhour = xhour + tstep tjd = tjd1 + yhour/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) yplus = sind(alt) - sin_h0 c..if this is the largest altitude, capture the transit if (alt .gt. altmax) then altmax = alt utcul = yhour end if c..find the parabola through the yminus, y0, and yplus call quad(yminus,y0,yplus,xe,ye,zero1,zero2,nz) c..if there is only one zero crossing if (nz .eq. 1) then if (yminus .lt. 0.0) then utrise = xhour + zero1 * tstep rise = .true. azirise = azi else utset = xhour + zero1 * tstep set = .true. aziset = azi end if c..if there are two zero crossings else if (nz .eq. 2) then if (ye .lt. 0.0) then utrise = xhour + zero2 * tstep utset = xhour + zero1 * tstep azirise = azi aziset = azi else utrise = xhour + zero1 * tstep utset = xhour + zero2 * tstep azirise = azi aziset = azi end if rise = .true. set = .true. end if c..set up for the next loop through yminus = yplus xhour = xhour + 2.0d0 * tstep if (xhour.gt. 24.0) goto 20 end do c..start the output 20 continue c..form the calander date call caldat(tjd1+zoneday,iy,im,id,xh) write(adate(jj),123) id,amon(im),iy 123 format(i2.2,a,i4.4) c..set the transit time if (utcul .ne. -1.0d0) then call dh2h(utcul,utculh,utculm,utculs) write(transit_string(jj),03) int(utculh), 1 nint(utculm+utculs/60.0d0),altmax else write(transit_string(jj),04) end if c..if there was a rise or set time found if (rise .or. set) then if (rise) then call dh2h(utrise,utriseh,utrisem,utrises) write(rise_string(jj),03) int(utriseh), 1 nint(utrisem+utrises/60.0d0),azirise else write(rise_string(jj),04) end if if (set) then call dh2h(utset,utseth,utsetm,utsets) write(set_string(jj),03) int(utseth), 1 nint(utsetm+utsets/60.0d0),aziset else write(set_string(jj),04) end if c..no rise or set was found, c..so its either always above or below the horizon else if (above) then write(rise_string(jj),05) 'always' write(transit_string(jj),05) 'above' write(set_string(jj),05) 'horizon' else write(rise_string(jj),05) 'always' write(transit_string(jj),05) 'below' write(set_string(jj),05) 'horizon' end if endif c..end of day loop enddo return end subroutine quad(yminus,y0,yplus,xe,ye,zero1,zero2,nz) implicit none save c..fits a parabola through 3 points c..input: c..yminus = left value c..y0 = middle value c..yplus = right value c..output: c..xe = extreme x value of the parabola c..ye = extreme y value of the parabola c..zero1 = first root, if there is one c..zero2 = second root, if there is one c..nz = number of roots c..declare the pass integer nz double precision yminus,y0,yplus,xe,ye,zero1,zero2 c..local variables double precision a,b,c,dis,dx nz = 0 a = 0.5d0*(yminus + yplus) - y0 b = 0.5d0 * (yplus - yminus) c = y0 xe = -0.5d0*b/a ye = (a*xe + b) * xe + c dis = b*b - 4.0d0*a*c if (dis .ge. 0.0) then dx = 0.5d0 * sqrt(dis)/abs(a) zero1 = xe - dx zero2 = xe + dx if (abs(zero1) .le. 1.0) nz = nz + 1 if (abs(zero2) .le. 1.0) nz = nz + 1 if (zero1 .lt. -1.0) zero1 = zero2 end if return end include 'position_routine.f' include 'glue.f'