! %W% (ESO-IPG) %G% %U%
! "@(#) $Id: efficiency.prg,v 1.14 2005/12/15 15:39:31 jpritcha Exp $"
!==========================================================================
!
!    efficiency.prg
!    O. Stahl, ESO, 8-Oct-1998
!
!    Compute FEROS efficiency + telescope + fibres
!
!    Syntax:  @@ efficiency [filenum] [standard] [fibre]
!
!    with 
!     filenum : the 4-digit number of the spectrum of the object, 
!    standard : the name of the MIDAS tabe containin the standard star.
!       fibre : Which fibre to measure the efficiency for.
!###########################################################################

def/par P1 ? I "Enter object frame number: "
def/par P2 ? I "       Standard star name: "
def/par P3 1 N "Fibre (1) : "
def/par P4 MIDASHOME:/calib/data/spec/ctio    CHAR "Spectral flux directory"

def/local N/i/1/1 0
def/local i/i/1/1 0
def/local j/i/1/1 0
def/local start/d/1/1 0.0
def/local effi/d/1/1 0.0
def/local step/d/1/1 0.0
def/local end/d/1/1 0.0
def/local phot/d/1/1 0.0
def/local npoints/d/1/1 0.0
def/local nstand/d/1/1 0.0
def/local diff/d/1/1 0.0
def/local num/i/1/1 0
def/local texp/r/1/1 0.0
def/local mergemtd/c/1/20 " " all

def/local numcoun/d/1/1 0.0
def/local wave/d/1/3 0.0 all
def/local eff1/d/1/1 0.0
def/local eff2/d/1/1 0.0

def/local EMPREF/c/1/2      " "
def/local EFFLINE/r/1/1     0.
def/local ordb/r/1/2        0.,0.
def/local mjd/d/1/1         0.
define/local EffArea/r/1/1  0.
define/local order/i/1/1    0

define/local OFLABEL/C/1/100 " "
define/local mtfits/c/1/4    "mt"
define/local xbin/i/1/1       1
define/local ybin/i/1/1       1

i = m$existk("ARCFILE")
if i .eq. 0 DEFINE/LOCAL ARCFILE/C/1/100 " " ? +lower_levels
ARCFILE = "fero{p1}.mt"
i = m$existd("fero{p1}.mt","ARCFILE")
if i .eq. 1 ARCFILE = "{fero{p1}.mt,ARCFILE}"
!! Output FITS & PS files in VLT compliant name format, i.e.:
!! r.<ARCFILE>.nnnn.fits
!! To set either:
!! In MIDAS
!!   write/key  UseAFDRSFITS/C/1/1 "Y"
!! or in .bashrc
!!   export FEROSUseAFDRSFITS="Y"
!! or in .cshrc or .tcshrc
!!   setenv FEROSUseAFDRSFITS="Y"
i = m$existk("UseAFDRSFITS")
if i .ne. 1 then
  DEFINE/LOCAL UseAFDRSFITS/C/1/1 " " ? +lower_levels
  UseAFDRSFITS = M$SYMBOL("FEROSUseAFDRSFITS")
endif
OFLABEL = "{P1}{P3}"
if "{ARCFILE}" .ne. "fero{p1}.mt" .and. "{UseAFDRSFITS}" .EQ. "Y" then
  mtfits = "fits"
  OFLABEL = ".{ARCFILE(1:29)}.{P3}"
endif

i = m$existd("fero{p1}.mt","ESO.DET.WIN1.BINX") * m$existd("fero{p1}.mt","ESO.DET.WIN1.BINY")
if i .eq. 1 then
  xbin = {fero{p1}.mt,ESO.DET.WIN1.BINX}
  ybin = {fero{p1}.mt,ESO.DET.WIN1.BINY}
endif

crea/tab efficiency 2 10000
crea/col efficiency :wave
crea/col efficiency :eff
crea/col efficiency :order

EMPREF = "bf"
if EXT_MODE(1:1) .EQ. "B" then
  EMPREF = "bo"
elseif EXT_MODE(1:1) .EQ. "O" then
  EMPREF = "bo"
endif

REBIN/FERO {EMPREF}{p1}ext {LINE_POS_TBL} rebinned I {REBIN_STEP},0

mergemtd = "{MERGE_MTD}"
MERGE/FERO rebinned eff 1,39 NOAPP
SET/FEROS MERGE_MTD={mergemtd}

write/out "Extinction correction for airmass = {{EMPREF}{p1}ext{P3},O_AIRM}"
i = 0
do order = {ORDER_FIRST} {ORDER_LAST} 1
  i = i+1
  write/out "Order {i} [{order}]: "
  copy/dd {EMPREF}{p1}ext{P3} *,3 eff{P3}{i}
  EXTINC/SPEC eff{P3}{i} ext{P3}{i} ? MID_EXTINCTION:atmoexan 
enddo


texp = {{EMPREF}{p1}ext{P3},O_TIME(7)}
write/out "---------------------------------------------------------------"
write/out "The following numbers have been extracted from the descriptors:"
write/out "Exposure time = {texp} seconds"
write/out "CCD gain      = {CCD_GAIN} e-/ADU"
write/out "---------------------------------------------------------------"

!!
!! Note: For the computation of phan below...
!!       For CTIO flux table, then #2 is fine (and it is the old
!!       setting). If you use MIDAS tables, then put #4 below. Basically,
!!       these must be the column nr. of the monochromatic magnitudes.
!!
mjd = {{EMPREF}{p1}ext{P3},O_TIME(4)}
!!! UTC 2002-Oct-01T12:00:00 == JD 2452549.00000 == MJD 52548.5
!!! UTC 2003-Oct-01T12:00:00 == JD 2452914.00000 == MJD 52913.5
if mjd .lt. 52548.5 then
  !! ESO-1.52m
  ! Compute the photons /A at the 1.5m telescope
  ! Area M1:      1.824 m^2 = pi*(1.524m/2)^2
  ! Area M2:      0.147 m^2 = pi*(0.432m/2)^2
  ! -> effective Area       = [Area M1] - [Area M2] = 1.678m^2
  ! Factor actually used historically... 1.5
  EffArea = 1.5

elseif mjd .lt. 52913.5 then
  !! MPG/ESO-2.20m with 3 Aluminium Mirrors (M1+M2+M3)
  ! Compute the photons /A at the 2.2m telescope
  ! Area M1:      3.801 m^2 = pi*(2.2m/2)^2
  ! Area M2:      0.618 m^2 = pi*(0.887m/2)^2
  ! -> effective Area       = [Area M1] - [Area M2] = 3.183m^2
  EffArea = 3.183

else
  !! MPG/ESO-2.20m with 2 Aluminium Mirrors (M1+M2) and High reflectivity M3
  ! Compute the photons /A at the 2.2m telescope
  ! Area M1:      3.801 m^2 = pi*(2.2m/2)^2
  ! Area M2:      0.618 m^2 = pi*(0.887m/2)^2
  ! -> effective Area       = [Area M1] - [Area M2] = 3.183m^2
  EffArea = 3.183

endif

!!

!! The flux to be compared with is obtained from tables of
!! Spectrophotometric standards, either from
!! provided at ftp://...
!! These are used instead of the MIDAS calib ones because then
!! everything is in a standard format...
!! The flux (f[\nu]) MUST be in units of milli-Jansky,
!! i.e. 1.0E-26 x erg/s/cm^2/Hz, and tabulated in the column :f_ny, or
!! in magnitudes (m[\nu]) and tabulated in the column :mag.
!
!     m[\nu]  = -2.5 log10[ f[\nu] ] + 16.410
!             = -2.5 log10[ f[\lambda] * \lambda^2 / c ] + 16.410
!
!! The factor 16.410 comes from the Hamuy et al calibration (48.590) which
!! is appropriate for f[\nu] in erg/s/cm^2/Hz.
!
!     m[\nu]  = -2.5 log10[ f[\nu] ] - 48.590
!             = -2.5 log10[ f[\lambda] * \lambda^2 / c ] - 48.590

copy/tab {P4}/{P2} work
N = m$existc("work","wave")
if N .eq. -1 comp/table work :wave = :lambda
N = m$existc("work","mag")
if N .eq. -1 comp/table work :mag  = -2.5 * log10(:f_ny) + 16.410
!! Account for binning only in one (y) dimension since the other dimension (x) is collapsed anyway...
comp/table work :phan = {texp}*(({EffArea}*5.5*10**10)/:WAVE)*10**(-0.4*:mag)*{ybin}


N = 0
do order = {ORDER_FIRST} {ORDER_LAST} 1
  N = N+1
  start = {ext{P3}{N},START}
  step = {ext{P3}{N},STEP} 
  npoints = {ext{P3}{N},NPIX(1)}
  nstand = {work.tbl,tblcontr(4)}
  do j = 1 {nstand}
    wave(1) = {work,:wave,@{j}}
    wave(2) = wave(1)-1.
    wave(3) = wave(1)+1.
    end = {START}+{NPOINTS}*{STEP}
    if {wave(2)} .gt. {START} then
      if wave(3) .lt. {end} then
        phot = {work,:phan,@{j}}
        stat/ima ext{P3}{N}[{wave(2)}:{wave(3)}] >Null
        numcoun = {outputr(3)}
        if {numcoun} .gt. 1 then
          effi = ({CCD_GAIN}*{numcoun}) / ({phot} * {wave} * 1.e-5 )
          write/out "WL/Ord={wave(1)}/{N}/[{order}] {numcoun} Phot={phot} Eff={effi}"
          num = num+1
          efficiency,:wave,@{num} = {wave(1)}
          efficiency,:eff,@{num} = {effi}
          efficiency,:order,@{num} = {n}
        endif
      endif
    endif
  enddo
enddo
!
sort/table efficiency :wave
crea/column efficiency :wave{P3}
crea/column efficiency :eff{P3}
nstand = {efficiency.tbl,tblcontr(4)}
do j = 1 {nstand}
 n = j+1
 wave(1) = {efficiency.tbl,:wave,@{j}}
 wave(2) = {efficiency.tbl,:wave,@{n}}
 eff1 = {efficiency.tbl,:eff,@{j}}
 eff2 = {efficiency.tbl,:eff,@{n}}
 diff = wave(1)-wave(2)
 if M$abs({diff}) .lt. 0.5  then
  effi = {eff1} + {eff2}
  efficiency,:wave{P3},@{j} = {wave(1)}
  efficiency,:eff{P3},@{j} = {effi}
 else
  efficiency,:wave{P3},@{j} = {wave(1)}
  efficiency,:eff{P3},@{j} = {eff1} 
 endif
enddo
!
sele/tab efficiency all
copy/tab efficiency efficiency{p1}{P3}
write/out "-------------------------------------------------------------------"
write/out "Computed efficiencies for fiber {P3} stored in table efficiency{p1}{P3}"
write/out "Use columns :wave{P3} :eff{P3}"
write/out "-------------------------------------------------------------------"
!
@@ efficiency_plots {p1} {p2} {p3}
!
outd/fits efficiency{p1}{P3}.tbl efficiency{OFLABEL}.{mtfits}