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 ********************************************************