c     EXPOSE - calculate the exposure time to reach a given S/N ratio 
c     for a given fixed photometry aperture and magnitude. Time is
c     calculated for all bands: UBVRI JHK (griz not yet supported; 
c     zero-points and nepers in place, but no good sky values for these
c     bands). Default sky values and CCD quantum efficiencies are used 
c     (see below) unless otherwise indicated and supplied. Zero-point 
c     is in Jy. To get to photons per band, the relation 
c              1 muJy = 15.1 photons/sec-m**2-neper 
c     is used. Nepers are defined here as dlambda/lambda.
c
c     Author: Matthew A. Bershady
c     
c     Modification History
c     Feb 16 1992   MAB  A quick hash.
c     Jan 18 1993   MAB  Add tr calculation: time when photon noise
c                        equals shot-noise
c     
      integer i
      real t(12),b(12),zp(12),neper(12),qe(12),nphot(12)
      real ge,mlim,diameter,area,snr,apert,rapert,rn,pscale
      real aa,bb,cc,tr(12)
      character band(12)*1
c
c     Default sky background is from Kron and Spinrad (KS), except for
c     griz bands which are faked (i.e. not given by KS). Default QE
c     is for CFHT SAIC3 CCD for UBVRI, and KPNO IRIM for JHK, and griz
c     is undefined (set to 0.5). Note that I find a total efficiency of 
c     8% for IRIM+4-meter, compared to 11% quoted for IRIM on 1.3 meter.
c
      data b / 22.4,22.9,22.0,21.2,19.9, 15.2,13.5,13.5,
     &         22.0,21.7,22.0,22.0 /
      data qe / 0.15,0.25,0.40,0.45,0.20, 0.80,0.60,0.75, 
     &          0.50,0.50,0.50,0.50 /
      data zp / 1810.,4260.,3640.,3080.,2550., 1600.,1080.,670.,
     &          3730.,4490.,4760.,4810. /
      data neper / 0.15,0.22,0.16,0.23,0.19, 0.16,0.23,0.23,
     &             0.14,0.14, 0.16,0.13 /
      data band / 'U','B','V','R','I','J','H','K','g','r','i','z' /
c
      call readparam ( mlim,snr,rapert,diameter,rn,pscale,ge,qe,b)
      area = 3.14159*(diameter/2.)**2
      apert = 3.14159* rapert**2
      do 30 i=1,12
         nphot(i) = 15.1*10**6 * zp(i) * neper(i)
 30   continue
c
      do 40 i=1,12
         sky = apert*ge*qe(i)*area*nphot(i)*10**(-0.4*b(i))
         obj = ge*qe(i)*area*nphot(i)*10**(-0.4*mlim)
         aa = obj**2
         bb = -1.*snr**2*(obj+sky)
         cc = -1.*snr**2*(apert/pscale**2)*rn**2
         t(i) = (-1.*bb+sqrt(bb**2-4.*aa*cc))/(2.*aa)
         tr(i) = (apert/pscale**2) * rn**2 / (sky+obj) 
 40   continue
c
      write (*,120) '#                 S/N = ',snr
      write (*,120) '# photometry aperture = ',apert,' (arcsec^2)'
      write (*,120) '#      telescope area = ',area,' (m^2)'
      write (*,120) '#          read noise = ',rn,' (e-/pix)'
      write (*,120) '#         plate scale = ',pscale,' (acrsec/pix)'
      write (*,115) '#',' ','photons/  '
      write (*,110) '#','Band','m^2-sec','sky','e','m(lim)',
     &             't(rn)[sec]','t(snr)[sec]'
      write (*,115) '#',' ','(m=0)    '
      do 50 i=1,5
         write (*,100) band(i),nphot(i),b(i),ge*qe(i),mlim,tr(i),t(i)
 50   continue
c
 100  format (a5,e12.3,f7.2,f8.3,f7.2,2e12.3)
 110  format (a1,a5,a10,2a7,a10,a12,a12)
 115  format (a1,a7,a11 )
 120  format (a24,f7.2,a14 )
c
      stop
      end
c ---------------------------------------------------------------------
      subroutine readparam (mlim,snr,rapert,diameter,rn,pscale,ge,qe,b)
c ---------------------------------------------------------------------
c
      real  mlim,snr,rapert,diameter,rn,pscale,ge,qe(12),b(12)
      character sfile*40,qefile*40
c
      write (*,*) 'Magnitude limit:'
      read (*,*) mlim
      write (*,*) 'S/N limit:'
      read (*,*) snr
      write (*,*) 'Photometry aperture (radius in arcsec):'
      read (*,*) rapert
      write (*,*) 'Effective Telescope Diameter (meters):'
      read (*,*) diameter
      write (*,*) 'Readnoise (e- rms):' 
      read (*,*) rn
      write (*,*) 'plate scale (arcsec/pix):' 
      read (*,*) pscale
      write (*,*) 'Wavelength independent throughput:' 
      read (*,*) ge
      write (*,*) 'QE file (list: UBVRI JHK griz):'
      read (*,*) qefile
      if ( qefile.ne.'none' ) then
         open (1, file=qefile)
         do 10 i=1,5
            read (1,*) qe(i)
 10      continue
         close (1)
      endif
      write (*,*) 'Sky file (list: UBVRI JHK griz):'
      read (*,*) sfile
      if ( sfile.ne.'none' ) then
         open (1, file=sfile)
         do 20 i=1,5
            read (1,*) b(i)
 20      continue
         close (1)
      endif
c
      return
      end