c  Copyright 1999 by Robert J. Nemiroff and J. Bruce Rafert (MTU).
c  Archived as http://ascl.net/nemiroff990513.html
c  Disclaimer: no known flaws but not guaranteed accurate.
c  consky.for: computes telescope times for sky observation projects.
c  Updated Noise calculation on 2000 April 6, should not affect paper.
c  See http://xxx.lanl.gov/abs/astro-ph/9807096 .
       program consky
       implicit real*8 (a-z)
       parameter (nB=20)
       real*8 mapp(nB)
       real*8 sig10(nB), sig10m(nB), sig30(nB), sig30m(nB)
       real*8 sig50(nB), sig50m(nB), sig70(nB), sig70m(nB)
       real*8 sig90(nB), sig90m(nB)
       real*8 Nstar(nB)
       real*8 Npixel(nB), Npoint(nB), Ntel(nB)
       real*8 Data_epoch(nB), Data_conti(nB)
       real*8 Omega_pixel(nB), f(nB), r(nB)
       real*8 tepoch(nB), tcontinuous(nB)
c                                  Hardwire telescope parameters.
       nCCD=4096.0d0
       alpha=25.0d0
       beta=16.0d0
       Ntelmin=4.0d0
       b=4.0d0
       p=1.0d-5
       msky=21.0d0
       fratio=2.0d0
       eCCD=0.5d0
       Ddark=1.0*(p/1.0d-5)**2
       SNRmin=100.0d0
       Rreadout=1.0d5
       enight=0.25d0
       pi=dacos(-1.0d0)
c                                        Open star density files.
       open (1, file='conden.out', status='old')
       open (2, file='mn.out', status='old')
c                                        Read in star density by JB.
         do 100 i=5, nB 
         read (1,20) mapp(i), sig90m(i)
         read (1,30)          sig70m(i)
         read (1,30)          sig50m(i)
         read (1,30)          sig30m(i)
         read (1,30)          sig10m(i)
         read (1,30)           sig10(i)
         read (1,30)           sig30(i)
         read (1,30)           sig50(i)
         read (1,30)           sig70(i)
         read (1,20) mtest,    sig90(i)
  20     format (1x, f5.1, 14x, e11.4)
  30     format (20x, e11.4)
c
         if (mtest.ne.mapp(i)) stop ' Read magnitude error.'
c                                       Surface density of stars at B.
         sig90(i)=(sig90(i) + sig90m(i))/2.0d0
         sig70(i)=(sig70(i) + sig70m(i))/2.0d0
         sig50(i)=(sig50(i) + sig50m(i))/2.0d0
         sig30(i)=(sig30(i) + sig30m(i))/2.0d0
         sig10(i)=(sig10(i) + sig10m(i))/2.0d0
c                                       Number of pixels --> 25 sig10
         Npixel(i)=alpha*sig10(i)*4.0d0*pi
  100    continue 
c                                       Npixel --> 25 Nstar
         do 150, i=1, nB
         read (2,*) mapp(i), Nstar(i)
c        write (*,*) mapp(i), Nstar(i)
         Npixel(i)=alpha*Nstar(i)
  150    continue
c                                       Compute tel, data limits.
       do 200 i=1, nB
       Npoint(i)=Npixel(i)/nCCD**2/Ntelmin
       if (Npoint(i).lt.1.0d0) Npoint(i)=1.0d0
       Ntel(i)=Npixel(i)/nCCD**2
       if (Ntel(i).lt.4.0d0) Ntel(i)=4.0d0
c                                       Epochal Data.
       Data_epoch(i)=b*Ntelmin*Npoint(i)*nCCD**2
       Data_conti(i)=b*Ntel(i)*nCCD**2
  200  continue
c                                       Compute Omega_pixel, f.
         do 300 i=1, nB
         Omega_tel=4.0d0*pi/Ntel(i)
         Omega_pixel(i)=Omega_tel/nCCD**2
         f(i)=p/sqrt(Omega_pixel(i))
         diameter=f(i)/fratio
         r(i)=diameter/2.0d0
  300    continue 
c
       write (*,*) ' mapp, tepoch, tcontinuous '
c                                       Loop over star brightness.
       do 800 i=1, nB
       Istar=5.9d5*10.0d0**(-mapp(i)/2.5d0)
       Isky= 8.5d5*10.0d0**(-msky/2.5d0)
       Nexposure=0.99d0
       Cbackground=0.0d0
       Csignal=0.0d0
       tm=0.0d0
c                                       Artificial integration.
c                                       One exposure, one second at a time.
        do 500 t=1.0d0, 1.0d8 
        Cstar=Istar*t*4.0d0*(r(i)*100.0d0)**2*eCCD
        Csky =Isky *t*4.0d0*(r(i)*100.0d0)**2*eCCD*
     2        Omega_pixel(i)*206265.0d0**2
        Cdark=Ddark*t
        Cread=0.0d0
        Cbackground=Csky + Cdark + Cread
        Cnoise=sqrt(Cstar+Cbackground)
        Csignal=Cstar
        SNR=Csignal/Cnoise
        if (SNR.ge.SNRmin) goto 700 
        Ctotal=Csignal + Cbackground
        if (Ctotal.ge.2.0d0**beta) goto 550
  500   continue
c
       write (*,*) ' Error! Pixel never saturated. '
  550  texposure=t
c      write (*,*) mapp(i), t, SNR
c                                      Now generalize to many exposures.
         do 600 Nexposure=1.0d0, 1.0d8, 0.001d0 
         Cback=Nexposure*Cbackground
         Csignal=Nexposure*Cstar
         Cnoise=sqrt(Csignal+Cback)
         SNR=Csignal/Cnoise
c        write (*,*) mapp(i), Nexposure, SNR
         if (SNR.ge.SNRmin) goto 700
  600    continue
c
       write (*,*) ' Error! Signal/Noise never saturated. '
c
  700     tm=t*Nexposure
          tdark=tm
          treadout=int(Nexposure+1)*b*nCCD**2/Rreadout  
          t1=tm + tdark + 2.0d0*treadout
          tepoch(i)=int(Npoint(i)+1)*t1/enight
          tcontinuous(i)=t1/enight
c
       Repoch=Data_epoch(i)/tepoch(i)*(60.0d0*60.0d0*24.0)
       write (*,*) mapp(i), t1, tepoch(i)
c      write (*,*) mapp(i), tepoch(i), tcontinuous(i)
  800  continue
c
  999  stop
       end