! @(#)irac_focus.prg	17.1.1.1 (ES0-DMD) 01/25/02 17:52:42
! @(#)irac_focus.prg	17.1.1.1 (ESO-Chile) 01/25/02 17:52:42
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!.COPYRIGHT   (C) 1995 European Southern Observatory
!.IDENT       irac_focus.prg
!.AUTHOR      C. Lidman,  ESO-Chile
!.KEYWORDS    IRAC package
!.PURPOSE     procedure to determine the correct focus
!.NOTE        Adapted from comb.prg 
!.VERSION     950826 C. Lidman Creation
!             951010 CEL
!             951031 CEL Selection by sequence name corrected
!-------------------------------------------------------
crossref select seqname focout create

define/par p1 null C "Enter sequence name:"
define/par p2 focus C "Name of output/input image:"
define/par p3 Y C "Create Image:"

define/maxpar 3

write/out {p1},{p2},{p3}

!Check Input parameter bounds.

def/local  last/I/1/1 100
def/local  n/I/1/1 1
def/local  upper/I/1/1 1
def/local  readout/I/1/1 2
def/local  gain/R/1/1 6.6
def/local  ndit/I/1/1/ 10
def/local  ron/r/1/1 49.
def/local  med/r/1/1 0.0
def/local  scale/r/1/1 0.49
def/local  lens/c/1/2 LC
def/local  filter/c/1/4/ K
def/local  s_lens/c/1/2 LC
def/local  s_filter/c/1/4/ K
def/local  s_readou/I/1/1 2
def/local  error/I/1/1 0
def/local  cont/C/1/1 Y
def/local  exist/I/1/1 1
def/local  medran/I/1/1 0
def/local  focus/R/1/1 3000.0
def/local  a2/d/1/1 0.0
def/local  a3/d/1/1 0.0
def/local  xfocus/d/1/1 3000.0
def/local  yfocus/d/1/1 3000.0
def/local  average/d/1/1 3000.0
define/loc focname/C/1/20 "focfunction"
IF M$EXIST("{focname}.fit") .EQ. 0 THEN  
   IF AUX_MODE(1) .LE. 1 THEN                                    ! VMS
      define/local fmtdir/c/1/60 -
         "MID_DISK:[&MIDASHOME.&MIDVERS.STDRED.IRAC2.LIB]"
   ELSE                                                          ! UNIX
      define/local fmtdir/c/1/60 -
         "$MIDASHOME/$MIDVERS/stdred/irac2/lib/"
   ENDIF
   -COPY {fmtdir}{focname}.fit {focname}.fit                       ! copy
ENDIF

copy/table irac2b_ost work

delete/col work :MJD :DATE :AZIMUTH :ALTITUDE

select/table work :SEQID .eq. {p1}

copy/table work work2

set/midas output=no
show/tab work2
last = {outputi(2)}
set/midas output=yes

if p3 .eq. "Y" then

!Determine the readout mode

copy/dk {work2,FILENAME,@1} _ED_NCORRS readout
copy/dk {work2,FILENAME,@1} _ED_NDIT ndit
copy/dk {work2,FILENAME,@1} _EIO3_NAME lens
copy/dk {work2,FILENAME,@1} _EIO2_NAME filter

If lens .eq. "LA" then
  scale = 0.151
elseif lens .eq. "LB" then
  scale = 0.278
elseif lens .eq. "LC" then
  scale = 0.507
elseif lens .eq. "LD" then
  scale = 0.708
elseif lens .eq. "LE" then
  scale = 1.061
endif

!Check the compatability of filters, lens, readout mode etc.
!This could be extended to other optical elements

do n = 2 {last} 1

copy/dk {work2,FILENAME,@1} _ED_NCORRS s_readou
copy/dk {work2,FILENAME,@1} _EIO2_NAME s_filter
copy/dk {work2,FILENAME,@1} _EIO3_NAME s_lens

if {s_readou} .ne. {readout} then
  write/out Incompatable readout modes.
  error = 2 
endif

if "{s_filter}" .ne. "{filter}" then
  write/out Incompatable filters.
  error = 1
endif

if "{s_lens}" .ne. "{lens}" then
  write/out Incompatible lenses.
  error = 2
endif

enddo

if {error} .eq. 2 then
   return
elseif {error} .eq. 1 then
   inquire/key cont "Do wish to continue (Y/N):"
   if cont .eq. "N" then
       return
   endif
endif

create/icat obj.cat work2,:FILENAME

!The average command will need to be replaced with the combine/ccd command,
!as the combine/ccd command is likely to create the better skies.
!The Combine/ccd command requires, at the very least, an association table to
!work. There is a strange problem with the calculation of the mode.
!For now we shall stick with the average command.

write/out skynum={last}

medran = {last}-3

if {medran} .le. 0 then
  medran = 1
endif

average/images sky = obj.cat ? ? min,{medran}

average/images obj = obj.cat
comp/ima {p2} = obj - sky

!Display the image and ask the user to select the object to determine the focus

endif

@@ irac_acuts {p2}

center/gauss CURSOR focus ? ? 10,10

create/col focus :focus r

do n = 1 {last} 1
  copy/dk {work2,FILENAME,@{n}} _ET_FOC_POS focus
  write/table focus :focus @{n} {focus}
enddo

set/graph binmode=off

plot/table focus :focus :xsig ? 2 0
set/graph colour=2
overplot/table focus :focus :ysig ? 3 0
set/graph colour=1

set/midas output=no
fit/table 100,0.000001,0.2 focus :xsig :focus focfunction
set/midas output=yes
compute/fit focus :xout = focfunction(:focus)
overplot/table focus :focus :xout ? 0 2

write/out

copy/dk focfunction.fit fitparam/d/2/1 a2
copy/dk focfunction.fit fitparam/d/3/1 a3
xfocus = {a2}/(-2.*a3)
write/out x-focus is {xfocus}

set/midas output=no
fit/table 100,0.000001,0.2 focus :ysig :focus focfunction
set/midas output=yes
compute/fit focus :yout = focfunction(:focus)
set/graph colour=2
overplot/table focus :focus :yout ? 0 2
set/graph colour=1

copy/dk focfunction.fit fitparam/d/2/1 a2
copy/dk focfunction.fit fitparam/d/3/1 a3
yfocus = {a2}/(-2.*a3)
write/out y-focus is {yfocus}

average = ({xfocus} + {yfocus})/2.

write/out

write/out Average = {average}

!Clean up loose files

if p3 .eq. "Y" then

del/image sky NO
del/image obj NO

endif

del/tab work NO
del/tab work2 NO
del/tab focus NO