c	**************************************************************
	program twilight
c	version 1.0   [11 August 1992]
c	version 2.0   [30 September 1992]
c	Neil D. Tyson & Roy R. Gal
c	'An Exposure Guide for Taking
c		Twilight Flat Fields with Large Format CCDs
c       Astronomical Journal 105, 1206 (May 1993)

c---------------------------------------------------------------------
c	SEQUENCE...
c	Test exposure gives primary exposure
c	Primary exposure gives SS (= Stotal/So)
c	SS gives next exposure.
c	That exposure gives the next exposure...
c
c	SS = (a**t1/lna) - 1/lna   ...t1 = exposure1 since t0 = 0
c
c       t2 = ln(a**(t1 + trot) + lna * (SS))/lna
c	     exposure2 = t2 - (t1 + trot)
c
c	tn+1 = ln(a**(tn + trot) * lna * (SS))/lna
c	     exposuren = tn+1 - (tn + trot)
c----------------------------------------------------------------------

	implicit none

	integer n, know
	real expo(12), rot(13)
	real tau, k, twi,  ccdrot
	real lna, a, norm
	data rot/.25, .5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6/
	data expo/.25, .5, .75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6/

	open (10, file = 'twilight.tab')

	call twilite(twi, tau, norm)
	call getrot(know, ccdrot, rot)
	call setcon(k, tau, a, lna)
	call title

c	Loop through CCD readouts.  Bail out if read out time is known.
	do n=1,13
	    call tabhead(rot, tau, n)
	    call compute(expo, a, lna, rot, n)
	    write(6, *) ' '
	    write(10,*) ' '
	    if(know.eq.1) stop '****End Exposure Table*****'
	enddo

	stop '****End Exposure Tables*****'
	end

c	********************************************************

	subroutine twilite(twi, tau, norm)
c	get coefficient of twilight duration *or* twilight duration
c	from the user

	implicit none
	real twi, tau, norm

25	continue
	write(6,*) ' '
2	write(6,*) 'Twilight duration [in minutes], or Twilight coefficient:'
	read(5,*,err=2) twi
	
c	input > 60 is assumed to be the twilight duration in minutes
c	input < 20 assumed to be the twilight duration coefficient 'tau'
	if (twi.gt.20.and.twi.le.60) print *,'***invalid input***'
	if (twi.gt.20.and.twi.le.60) goto 25
	if (twi.le.20) tau = twi
	if (twi.gt.60) tau = twi/80.

c	duration of twilight for normalization, in minutes
	norm = 80.

	return
	end

c	********************************************************

	subroutine getrot(know, ccdrot, rot)
c	get the CCD read out time from the user, if known

	implicit none
	real ccdrot, rot(13)
	integer know

	write(6,*) ' '
1	write(6,*) 'Do you know your CCD read out time? [no:0, yes:1]'
	read(5,*,err=1) know
	write(6,*) ' '
	   if (know.eq.1) then
3		write(6,*) 'Read out time in seconds:'
		read(5,*,err=3) ccdrot
		rot(1) = ccdrot / 60.
	   endif

	return
	end

c	********************************************************

	subroutine setcon(k, tau, a, lna)
c	set the known constants

	implicit none
	real k, tau, a, lna

c	Slope of the  Log (S)-vs-time  curve  [units: 1/min]
c	for 80 minutes of Astronomical Twilight
	k = 0.09391

c	collect and rename terms
	a = 10.**(k/tau)
	lna = alog(a)

	return
	end

c	********************************************************

	subroutine title

	write(6,100)
	write(10,100)
100	format(42x,'****** EXPOSURE SEQUENCES FOR TWILIGHT FLATS ******'/)

	return
	end

c	********************************************************

	subroutine tabhead(rot, tau, n)
c	output table and column heads

	implicit none
	integer n
	real rot(13), tau

c	header lines
	write(6, 101) rot(n), tau
	write(10, 101) rot(n), tau
101	format(' CCD Readout time: ',f5.2,' minutes',
     &   24x,'[All entries in seconds]',41x,'Tau: ' f5.2)
	write(6,102)
	write(10,102)
102	format('   (1)   (2)   (3)   (4)   (5)   (6)   (7)   (8)   (9)',
     &	'  (10)  (11)  (12)  (13)  (14)  (15)  (16)  (17)  (18)',
     &  '  (19)  (20)  (21)  (22)')
	write(6,103)
	write(10,103)
103	format(' -----------------------------------------------------',
     &  '------------------------------------------------------'
     &  '------------------------')

	return
	end

c	********************************************************

	subroutine compute(expo, a, lna, rot, n)
c	compute exposure sequences and write to output file

	implicit none
	integer i, j, m, n
	real expo(12), rot(13), seq(23)
	real a, lna, ss, t, tnext

c	Loop through exposure sequence
	  do i = 1, 12
	    seq(1) = expo(i)

c	determine SS from first exposure (assume t0 = 0)
	    ss = (a**seq(1))/lna - 1./lna
	    t = seq(1)

c	Loop to compute for 22 exposures per R.O.T
	      do j = 1, 22
		tnext = (alog(a**(t + rot(n)) + lna * ss))/lna
		seq(j+1) = tnext - ( t + rot(n) )
		t = tnext

c		convert to seconds for output
		seq(j) = seq(j) * 60.
	      enddo

	    write (6,104) (seq(m), m=1,22)
	    write (10,104) (seq(m), m=1,22)
104	    format(22(f6.1))
	  enddo

	return
	end

c	********************************************************