C===========================================================================
C
C				ROTA
C
C===========================================================================
C
	character yn*1,month(12)*3,tel(6)*6,dweek(7)*3,line*80,tmp*180,
     +html*8,obs(6,8)*180,inst(6,8)*80,fast(6,8)*80,save*7,
     +months(12)*9
	dimension idmon(12),night(6,8),num(6),nt(6),nnt(6,8)
	real*8 djulian
C
	data tel/'1.9-m ','1.0-m ','0.75-m','0.5-m ','IRSF  ','SALT  '/
	data dweek /'Tue','Wed','Thu','Fri','Sat','Sun','Mon'/
	data month /'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug',
     +'Sep','Oct','Nov','Dec'/
	data months /'January','February','March','April','May','June',
     +'July','August','September','October','November','December'/
cdata idmon /31,29,31,30,31,30,31,31,30,31,30,31/
C*****
C   Feb=29 for 2000
       data idmon /31,28,31,30,31,30,31,31,30,31,30,31/
C*****
C
	write(*,'(''*************************************************************************************''/
     +            ''This program generates a Web page for one month of the Sutherland observing rota.''/
     +            ''The output is a file of the form Mmm.html where Mmm is the month. This is in the ''/
     +            ''markup language used on the Web. The program also produces a file Mmm.inp which ''/
     +            ''is a transcript of what you have entered on the keyboard. Modify this file and type''/
     +            ''    Newrota < Mmm.inp''/
     +            ''to generate a new version of the page.  If a note on an instrument configuration ''/
     +            ''is required, place an asterisk after the name (eg: RPCS*). The program will ask you''/
     +            ''for the name of the file containing the note.  The note itself is written as a normal''/
     +            ''text file. Convert it into HTML by typing''/
     +            ''    remarks filename''/
     +            ''The web page is now ready.''//)') 
C
	write(*,'(''When there is more than one observer, use / to separate the names. The support ''/
     +            ''astronomer can be included in the list using square brackets, e.g.,''/
     +            ''    Bloggs/Smith/[Laney]/NA-Marang''//)')
C IRSF,SALT
        write(*,'(''Add IRSF observers, SALT astronomers, SALT operators in the same way as for other''/
     +            ''telescopes, marking SA by an asterisk, e.g.,''/
     +            ''    *Encarni/Sandisa/Lucian''/
     +            ''************************************************************************************''//)')
C
C
	open(unit=2,file='rota.inp',status='unknown')
	write(*,'(''Year (eg: 1997)? '',$)')
	read(*,*) jyear
	write(2,1014) jyear
	iyear = jyear - 1900
C
10	write(*,'(''Month (1-12)? '',$)')
	read(*,*) mon
	write(2,1014) mon
        mon0 = mon
C
	html(1:8) = month(mon)//'.html'
	save(1:7) = month(mon)//'.inp'
C
C
	open(unit=1,file=html,status='unknown')
C
	m1 = mon - 1
	if(m1.lt.1) m1 = 12
	m2 = mon + 1
	if(m2.gt.12) m2 = 1
	write(1,2000) months(mon),jyear,month(m1),months(m1),month(m2),
     +months(m2)
2000	format('<HEAD><TITLE>SAAO Rota for ',a,i5,'</TITLE></HEAD>'/
     +'<BODY BGCOLOR="#FFFFFF">'/
     +'Back to <A HREF= "http://www.saao.ac.za/observing/rota/"> Rota</A>'/
     +': <A HREF= "',a,'.html"> ',a,'</A>'/
     +': <A HREF= "',a,'.html"> ',a,'</A><P>')
C
	write(1,1009) months(mon),jyear
1009    format('<CENTER><H2>SAAO Observing rota for ',a,i5,
     *   '</H2></CENTER>')
        write(1,1010)tel
1010	format('<TABLE BORDER=2 WIDTH=100% CELLSPACING=1 
     +  CELLPADDING=1>'/
     +'<TR><TH colspan=2>Date</TH><TH>M %</TH><TH>',
     +A6,'</TH><TH>',A6,'</TH><TH>',A6,'</TH><TH>',A6,'</TH>
     +<TH>',A6,'</TH><TH>',A6,'</TH><TH>Duty Astr</TH>')
C
	write(*,'(''Starting Tuesday date (eg: 28/3)? '',$)')
	read(*,1000) line
1000	format(a)
	k = lnblnk(line)
	write(2,1000) line(1:k)
	k = index(line,'/')
	line(k:k) = ' '
	read(line,*) iday,mon
	iday = iday - 1
C
C	
20	do i =1,6
	 num(i) = 0
	 do j=1,8
	  fast(i,j)(1:1) = ' '
	 enddo
	enddo
C

	do it = 1,6
	 ntot = 0
21	 write(*,1012) tel(it)
1012	 format('Enter data for ',a,':')
C
	 write(*,'(''Name (Start with * if an observer, or SALT astronomer)? '',$)')
	 read(*,1000) line
	 k = lnblnk(line)
         write(*,*)k
	 write(2,1000) line(1:k)
C
	 num(it) = num(it) + 1
	 n1 = lnblnk(line)
	 k = 0
         kobs = 0
C
	 do i = 1,n1
          if(line(i:i).eq.'*')then
            tmp(k+1:k+8) = '<strong>'
            k = k+8
            kobs = 1
          else
	    if(line(i:i).ne.'/') then
	       k = k + 1
	       tmp(k:k) = line(i:i)
	    else
               if(kobs.eq.1)then
                tmp(k+1:k+9) = '</strong>'
                kobs = 0
                k = k + 9
               end if
	       tmp(k+1:k+4) = '<BR>'
	       k = k + 4
	    endif
          end if
c        write(*,'(i4,i4,2x,a)')k,kobs,tmp(1:k)
	 enddo
           if(kobs.eq.1)then
            tmp(k+1:k+9) = '</strong>'
            kobs = 0
            k = k + 9
           end if
	 tmp(k+1:180) = ' '
C
	 nob = num(it)
	 obs(it,nob) = tmp
C
	 write(*,'(''Instrument? (blank for SALT, Sirius for IRSF)'',$)')
	 read(*,1000) line
	 k = lnblnk(line)
         n1 = 0
         do j=1,k
         if(line(j:j).ne.'/')then
           n1 = n1 + 1
           tmp(n1:n1) = line(j:j)
         else
           tmp(n1+1:n1+4) = '<BR>'
           n1 = n1 + 4
         end if
         end do
         tmp(n1+1:80) = ' '
         k = n1
	 write(2,1000) tmp(1:k)
	 if(index(tmp,'*').gt.0) then
	  write(*,'(''Filename for special remarks? '',$)')
	  read(*,1000) fast(it,nob)
	  k = lnblnk(fast(it,nob))
	  write(2,1000) fast(it,nob)(1:k)
	 endif
	 inst(it,nob) = tmp
C
	 write(*,'(''No of nights allocated [7]? '',$)')
	 read(*,1000) line
	 nts = 7
	 if(line(1:5).ne.'   ') read(line,*) nts
	 write(2,1014) nts
1014	 format(i5)
	 nnt(it,nob) = nts
	 if(nob.eq.1) then
	  night(it,nob) = 1
	 else
	  night(it,nob) = ntot + 1
	 endif
	 ntot = ntot + nts
	 if(mod(ntot,7).ne.0) goto 21
C
	enddo
C
	n = 0
	do i=1,6
	 nt(i) = 1
	enddo
C
	do i = 1,7
	 iday = iday + 1
	 if(iday.gt.idmon(mon)) then
	  mon = mon + 1
	  if(mon.gt.12) mon = 1
	  iday = 1
	 endif
C
C
	 id = iday
	 call moonphase(id, mon, iyear, 0, 0, 0,iphase,djulian)
	 write(1,1020) dweek(i),iday,iphase
1020	 format('<TR><TD>',a3,'</td><td align=center>',i02,
     *                       '</TD><TD ALIGN=CENTER>',i3,'</TD>')

	        if(i.eq.1)then
        write(*,'('' Type of week (B, G, D)? '')')
	  read(*,1000) line
c  n1 = lnblnk(line)
	   write(2,1000) line(1:n1)
C
	  k = 0
	  do j = 1,n1
	   if(line(j:j).ne.'/') then
	    k = k + 1
	    tmp(k:k) = line(j:j)
	   else
	    tmp(k+1:k+4) = '<BR>'
	    k = k + 4
	   endif
	  enddo
	  tmp(k+1:80) = ' '
c  write(1,1024) tmp(1:k)
1024	  format('<TD ALIGN=CENTER ROWSPAN=7>',a,'</TD>')
          end if
C
	 do it = 1,6
	  if(nt(it).le.num(it)) then
 	   k = nt(it)
	   if(night(it,k).eq.i) then
	    k1 = lnblnk(obs(it,k))
	    k2 = lnblnk(inst(it,k))
	    if(fast(it,k)(1:1).eq.' ') then
	     write(1,1021) nnt(it,k),obs(it,k)(1:k1),inst(it,k)(1:k2)
1021	     format('<TD ALIGN=CENTER ROWSPAN=',i2,'>',a,'<BR><BR>',
     +a,'</TD>')
	    else
	     nf = lnblnk(fast(it,k))
	     write(1,1022) nnt(it,k),obs(it,k)(1:k1),fast(it,k)(1:nf),
     +inst(it,k)(1:k2-1)
1022	     format('<TD ALIGN=CENTER ROWSPAN=',i2,'>',a,'<BR><BR>
     +<A HREF= "http://rota.saao.ac.za/',a,'.html">',a,'</A></TD>')
	    endif
C
	    nt(it) = nt(it) + 1
	   endif
	  endif
	 enddo
C
	 if(i.eq.1) then
	  write(*,'(''Duty Astronomer? '',$)')
	  read(*,1000) line
C
	  n1 = lnblnk(line)
	  write(2,1000) line(1:n1)
C
	  k = 0
	  do j = 1,n1
	   if(line(j:j).ne.'/') then
	    k = k + 1
	    tmp(k:k) = line(j:j)
	   else
	    tmp(k+1:k+4) = '<BR>'
	    k = k + 4
	   endif
	  enddo
	  tmp(k+1:80) = ' '
	  write(1,1025) tmp(1:k)
1025	  format('<TD ALIGN=CENTER ROWSPAN=7>',a,'</TD>')
	 endif
C
	 write(1,1031)
1031	 format('</TR>')
C
	enddo
C
	write(*,'(''Finished? '',$)')
	read(*,1000) yn
	write(2,1000) yn
	if(.not.(yn.eq.'y'.or.yn.eq.'Y')) goto 20
C
	call FDATE(line)
	k = lnblnk(line)
	write(1,1040) month(m1),months(m1),month(m2),months(m2),
     + month(mon0),line(1:k)
1040	format('</TABLE>'/
     +'The percentage of moon illuminated at local midnight is given '/
     +'in the third column.<br>Abbreviations: B, G and D: Bright,'/
     +'grey and dark weeks. NA: Night assistant.<br>'/ 
     +'Names in <strong>square brackets</strong> are visitor support'/
     +'astronomers, in <strong>bold</strong>'/
     +' are actual observers.<br>'/
     +'<font color="#00dd00"><strong>Dates '/
     +'of Public Holidays</strong>'/
     +'</font> marked in <strong>bold</strong>.<HR><P>'/
     +'Back to <A HREF= "http://www.saao.ac.za/observing/rota/">Rota</A>'/
     +': <A HREF= "',a,'.html"> ',a,'</A>'/
     +': <A HREF= "',a,'.html"> ',a,'</A>'/
     +': Print this <A HREF= "',a,'.ps"> version</A>'/
     +'<P>'/'Last updated: ',a/'</BODY>')
C
	close(1)
	call system('mv rota.inp '//save)
	stop
	end
C===========================================================================
C
C				ASTR
C
C===========================================================================
C
C
C**********************************************************************
C                            MOONPHASE
C
C A program to calculate the altitude of the Sun and Moon at any
C date and time.
C
C Written by:  Dr. L.A. Balona
C              South African Astronomical Observatory
C              P.O. Box 9
C              Observatory 7935
C              Cape Town
C              South Africa
C              Tel: (021)47-0025
C**********************************************************************/
C
	subroutine moonphase(iday, month, iyear, ihr, mn, isc,
     + iphase, djulian)
	implicit real*8 (a-h, o-z)
	parameter (frad = 0.017453292, twopi = 6.283185307)
	parameter (alat = -32.3795,along = -4994.6)
	parameter (time_zone = -2.0)
C
	std_time = 3600.0*ihr + 60.0*mn + isc
	ut = std_time + 3600.0*time_zone
C
	if ( ut.lt.0.0) then
	 ut = ut + 86400.0 
	 iday = iday - 1
	endif
C
	if ( ut .gt. 86400.0 ) then
	 ut = ut - 86400.0 
	 iday = iday + 1
	endif
C
	td = tdays(iday, month, iyear)
	fjd = ut/86400.00 + 0.5		 !   fractional julian iday 
	Tcen = (td + 0.5 )/36525.00	 ! No of centuries since 1900 Jan 0 
C
	sidt0 = 23925.836 + Tcen*(8640184.542 + 0.0929*Tcen)
	k = sidt0/86400.0
	sidt0 = sidt0 - 86400.0*k
	djulian = td + fjd + 2415020.000
C
C	Tcen = No of centuries since 1900
	Tcen = (djulian  - 2415020.000)/36525.00
C
	call sun (Tcen,ut,sidt,along,alat,slong,sun_ra,sun_dec,
     +sidt0 )
C
	call moon(alat, djulian, slong, sun_ra, sun_dec, ra_moon,
     +dec_moon, sidt, iphase)
C
	return
	end
C
C**********************************************************************
C                            TDAYS
C                Calculates no of days since 1900
C**********************************************************************
C
	function tdays(iday, month, iyear)
	implicit real*8 (a-h, o-z)
C
	ig = iyear !  Formula adopted from HP calculater handbook
	if ( month .lt. 3) ig = iyear - 1
	jf = month + 1
	if ( month .lt. 3 ) jf = month + 13 
C
	nt1 = 365.25*ig 
	nt2 = 30.6001*jf 
	tdays = nt1 + nt2 + iday - 64
	return
	end
C
C
C**********************************************************************
C                           SUN
C  Calculates sun's ra, dec and along and returns the sidereal time
C  sdt in seconds.
C
C  This algorithm is based to a large extent on  "Methods and formulae for
C  computations for Navigation and Astronomy", by Yallop (1978)
C  NAO technical note.
C**********************************************************************/
C
	subroutine sun (Tcen,ut,sdt,along,alat,slong,sun_ra,sun_dec,
     +sidt0 )
	implicit real*8 (a-h, o-z)
	parameter (frad=0.017453292)
C
	sunL = 279.697 + 36000.769 * Tcen	! Sun's mean along in degrees
	itemp = sunL/360.0 
	sunL = frad*(sunL - 360.0*itemp)	! Suns mean long in radians
	sunM = 358.476 + 35999.050 * Tcen 	! sunMean anomaly in degrees
	itemp = sunM/360.0
	sunM = frad*(sunM - 360.0*itemp) 	! Sun's mean anomaly in radians
	eps = frad*(23.452 - 0.013*Tcen) 	! Obliquity of ecliptic in radians 
C
C	Sun's true long:
	slong = sunL + frad*((1.919 - 0.005*Tcen)*sin(sunM) 
     ++ 0.020*sin(2.0*sunM))
	sun_ra = atan( tan(slong) * cos(eps) ) 	! Sun's RA 
	sun_dec = asin( sin(slong) * sin(eps) ) ! Sun's Dec
	if ( (slong.gt.1.570796).and.(slong.lt.4.71238898) )
     + sun_ra = sun_ra + 3.141592654		! Puts RA in right quadrant
	sun_nu = slong - 4.90822926 - 0.030*Tcen 	! True anomaly
	r =  0.9997210/(1.0 + 0.016751*cos(sun_nu)) ! Sun's distance in AU 
C
C	Now work out the Siderial time in seconds 
	r = sidt0  - along + 1.00273791*ut 	! Calculate LSIDT
	if ( r .lt. 0.0 ) r = r + 86400.0 
	k = r/86400.0
	sdt = r - k*86400.0
C
	return
	end
C
C*********************************************************************
C                               MOON
C                      Gets moon's altitude, phase
C**********************************************************************
C
	subroutine moon (alat, djulian, slong, sun_ra,
     +sun_dec, ra_moon, dec_moon, sidt, ip)
	implicit real*8 (a-h, o-z)
	parameter (frad = 0.017453292,twopi = 6.283185307)
C
	call moon_pos ( ra_moon,dec_moon,amlong,alat,djulian,
     +sidt )
C
	alt = 90.5 - zenith ( ra_moon, dec_moon, ha, alat, sidt ) 
	phs = amlong - slong 
C
	do while ( phs .gt. twopi )
	 phs = phs - twopi
	enddo
	do while ( phs .lt.0.0 )
	 phs = phs + twopi
	enddo
C
	if ( phs .gt. 0.5*twopi )  phs = twopi - phs
	ip = 200.0*phs/twopi 
C
	return
	end
C
C**********************************************************************
C                            MOON_POS
C                Calculates position of the moon given JD.
C**********************************************************************
C
	subroutine moon_pos ( ra, dec, amlong, alat, djulian,
     +sidt )
	implicit real*8 (a-h, o-z)
	parameter (frad = 0.017453292, twopi = 6.283185307)
C
	t = (djulian - 2451545.00)/36525.000
	am  =  218.32 + 481267.883*t
	am = am + 6.29 * sin ( residue (134.9 + 477198.85*t) )
	am = am - 1.27 * sin ( residue (259.2 - 413335.38*t) ) 
	am = am + 0.66 * sin ( residue (235.7 + 890534.23*t) ) 
	am = am + 0.21 * sin ( residue (269.9 + 954397.70*t) ) 
	am = am - 0.19 * sin ( residue (357.5 +  35999.05*t) ) 
	am = am - 0.11 * sin ( residue (186.6 + 966404.05*t) ) 
C
	bet  =  5.13 * sin ( residue ( 93.3 + 483202.03*t) ) 
	bet = bet +  0.28 * sin ( residue (228.2 + 960400.87*t) ) 
	bet = bet - 0.28 * sin ( residue (318.3 +   6003.18*t) ) 
	bet = bet - 0.17 * sin ( residue (217.6 - 407332.20*t) ) 
C
	am = residue ( am ) 
	bet = residue ( bet ) 
 	amlong = am 
C
	sl = cos(bet)*cos(am) 
	sm = 0.9175*cos(bet)*sin(am) - 0.3978*sin(bet) 
	sn = 0.3978*cos(bet)*sin(am) + 0.9175*sin(bet) 
C
C	Geocentric tan(RA) = sm/sl, sin(DEC) = sn.
C	Now obtain topocentric posn 
C
	am = frad*alat 
	bet = 7.272205e-5*sidt 
	x = 60.0*sl - cos (am)*cos (bet) 
	y = 60.0*sm - cos (am)*sin (bet) 
	z = 60.0*sn - sin (am) 
	r = sqrt ( x*x + y*y + z*z ) 
 	dec = asin ( z/r ) 
C
	if ( abs ( x) .lt. 1.0e-5 ) then
	 z = 1.0E+5 
	else
	 z = y/x 
	endif
C
	r = atan ( z ) 
	if ( x .lt. 0.0 ) r = r + 0.5*twopi 
	if ( r .lt. 0.0 ) r = r + twopi 
 	ra = r 
C
	return
	end

C
C**********************************************************************
C                            RESIDUE
C    Returns an angle in the range 0 - 360 and converts to radians.
C**********************************************************************
C
	function residue ( tht )
	implicit real*8 (a-h, o-z)
	parameter (frad = 0.017453292)
C
  	i = tht/360.0 
	residue = frad*(tht - 360.0*i) 
C
	return
	end
C
C**********************************************************************
C                             ZENITH
C             Returns the zenith distance in degrees
C**********************************************************************
C
	function zenith( ra, dec, hr_ang, alat, sidt )
	implicit real*8 (a-h, o-z)
	parameter (frad = 0.017453292)
C
	ha = 7.272205e-5*sidt - ra 
	cosz =  sin(frad*alat)*sin(dec) 
     +        + cos(frad*alat)*cos(dec)*cos(ha) 
	hr_ang = ha/frad 
	zenith = acos(cosz)/frad 
C
	return
	end
C

