! @(#)irac_dcomb.prg	17.1.1.1 (ES0-DMD) 01/25/02 17:52:42
! @(#)irac_dcomb.prg	17.1.1.1 (ESO-Chile) 01/25/02 17:52:42
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!.COPYRIGHT   (C) 1995 European Southern Observatory
!.IDENT       irac_dcomb.prg
!.AUTHOR      C. Lidman,  ESO-Chile
!.KEYWORDS    IRAC package
!.PURPOSE     procedure to combine dithered IR images
!.VERSION     950803 C. Lidman Creation
!             951010 CEL       Included automatic alignment and better sky 
!                              subtraction.
!             951105           Improved the sky subtraction.                  
!             951105           Added checks for automatic alignment.
!-------------------------------------------------------
crossref select seqname accsky align output trim tag

define/par p1 null C "Enter selection string:"
define/par p2 null C "Enter sequence name:"
define/par p3 N C "Accurate Sky subtraction before combining"
define/par p4 Y C "Align object frames (N/Y/A/B):"
define/par p5 ? C "Name of output image:"
define/par p6 Y C "Trim Image:" 
define/par p7 Y C "Tag all frames as skies:"

define/maxpar 7

write/out {p1},{p2},{p3},{p4},{p5},{p6},{p7}

def/local  last/I/1/1 100
def/local  xc/r/1/1 0.0
def/local  yc/r/1/1 0.0
def/local  n/I/1/1 1
def/local  upper/I/1/1 1
def/local  xshift/R/1/1 0.0
def/local  yshift/R/1/1 0.0
def/local  xstart/R/1/1 0.0
def/local  ystart/R/1/1 0.0
def/local  yend/R/1/1 0.0
def/local  xend/R/1/1 0.0
def/local  xlow/R/1/1 1.0
def/local  xhigh/R/1/1 256.0
def/local  ylow/R/1/1 1.0
def/local  yhigh/R/1/1 256.0
def/local  yshmin/R/1/1 0.0
def/local  yshmax/R/1/1 0.0
def/local  xshmin/R/1/1 0.0
def/local  xshmax/R/1/1 0.0
def/local  medran/I/1/1 0
def/local  kappa/R/1/1 1.5
def/local  thresh/R/1/1 0.0
def/local  avemed/R/1/1 1.0
def/local  skynum/I/1/1 1
def/local  skymed/R/1/1 0.0
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  med2/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  xapprox/r/1/1 0.0
def/local  yapprox/r/1/1 0.0
def/local  bias/r/1/1 1050.
def/local  shortcut/r/1/1 0.0
def/local  er/r/1/1 0.0
def/local  reload/C/1/1 F

if p3 .eq. "Y" then
   inquire/key kappa "Enter Kappa Threshold(Typically between 0.5 and 5.0):"
endif

copy/table irac2b_ost work
create/col work :COMBINE I*1
delete/col work :MJD :DATE :AZIMUTH :ALTITUDE


if p1 .ne. "null" then
   write/tab work :COMBINE {p1} 1
   select/table work :COMBINE .eq. 1
else
   if p2 .ne. "null" then
      select/table work :SEQID .eq. {p2}
   else
      write/out You must select either a sequence or a list of files
      return
   endif
endif

copy/table work work2

!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

if {readout} .eq. 4 then
  gain = 6.2
endif

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

!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,@{n}} _ED_NCORRS s_readou
copy/dk {work2,FILENAME,@{n}} _EIO2_NAME s_filter
copy/dk {work2,FILENAME,@{n}} _EIO3_NAME s_lens

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

if "{s_filter}" .ne. "{filter}" then
  write/out Incompatible 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

!Compute the median of each frame and write it to the table work2.

create/col work2.tbl :MEDIAN r

do n = 1 {last} 1

stat/image {work2.tbl,FILENAME,@{n}} ? ? ? FN
write/table work2 :MEDIAN @{n} {OUTPUTR(8)} 

enddo 

!Tag skies - this could be more sophisticated

create/col work2.tbl :EXPTYPE c*8

!If you do not want to tag all object frames as skies

if p7 .eq. "N" then
  write/descr work2.tbl SKY/C/1/50 ":IDENT.EQ."~*sky*" .OR. :TYPE .EQ. "SKY""
  classify/image work2 SKY :EXPTYPE sky
else
  write/table work2.tbl :EXPTYPE @1..{last} sky
endif

!One could add a procedure so that the user has the option of changing the 
!allocation made by the programme.

select/tab work2 :EXPTYPE .eq. "sky"

!From the median of each sky add a zero point offset so all medians
!are equal. Only then will pixel rejection work well.
!This should be updated with combine/ccd in the next version.

copy/tab work2 sky

stat/tab sky :median
create/col sky.tbl :ZEROIMAGE c*13
avemed = {OUTPUTR(3)}

set/midas output=no
show/tab sky
skynum = {outputi(2)}
set/midas output=yes

write/out skynum={skynum}

do n = 1 {skynum} 1
    write/table sky :ZEROIMAGE @{n} sky{n}.bdf
    comp/ima  sky{n} = {sky.tbl,FILENAME,@{n}}+{avemed}-{sky.tbl,median,@{n}}
enddo

create/icat sky.cat sky,:ZEROIMAGE

medran = {skynum}-3

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

!The following command is very crude. It should be replaced with the 
!COMBINE/CCD command.

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

stat/image sky ? ? ? FN
skymed = {OUTPUTR(8)}

select/table work2 ALL

!If tag all frames as skies is false

if p7 .eq. "N" then
   select/tab work2 :EXPTYPE .ne. "*sky*"
else
   select/tab work2 ALL
endif

!If alignment is not required

if p4 .eq. "N" then 
  create/icat obj.cat work2,:FILENAME
  average/images obj = obj.cat
  comp/ima obj = obj - sky
else

!Give the user the choice of centering all frames individually or
!allow the computer to rely on the approximate offsets from the 
!telesope. The telescope offsetting is not accurate enough for 
!blind offsetting. Nevertheless the option is available.

set/midas output=no
show/tab work2
upper = {OUTPUTI(8)}
! This gives me the number of selected rows
set/midas output=yes

copy/tab work2 work3
create/col work3 :SHIFILENAME ? ? C*15
write/out {upper}
create/tab center 1 10 null
create/col center :XSTART r
create/col center :XEND r
create/col center :YSTART r
create/col center :YEND r

do n = 1 {upper} 1

!First, subtract a scaled sky.

stat/image {work3.tbl,FILENAME,@{n}} ? ? ? FN

comp/ima skysub{n} = {work3.tbl,FILENAME,@{n}} - sky * {OUTPUTR(8)} / {skymed}

!Second, subtract the remaining median.

stat/ima skysub{n} ? ? ? FN

comp/ima skysub{n} = skysub{n} - {OUTPUTR(8)}

if p3 .eq. "Y" then 
!Using the median to determine the statistics.
   stat/image skysub{n} ? ? ? FN
   med = OUTPUTR(8)
!Calculate the expected sky noise from the count average.
!For DCR this requires the an approximate value for the bias. There
!could be an error if DCR was used with low background conditions.
   if {readout} .eq. 2 then
      med2 = {work3,MEDIAN,@{n}}+{bias}
      thresh = {med}+ 2.*{kappa}/{gain}*M$SQRT(({med2}+{ron})/{ndit}*{gain})
   else
      shortcut = M$SQRT((({work3,MEDIAN,@{n}})*{gain}+{gain}*{ron})/{ndit})
      thresh = {med}+ 2.*{kappa}/{gain}*{shortcut}
   endif
   write/out {thresh},{med}
   repl/image skysub{n} &a {thresh},>={med}
   repl/image &a &b <,-100={med}
   filter/sm &b &c 5,5,0.0 YA
!There are many filtering options; this one is one of the quickest.
   comp/ima skysub{n} = skysub{n} - &c
   del/image &a NO
   del/image &b NO
   del/image &c NO
endif

if {n} .eq. 1 .and. p4 .ne. "B" then  

!Should put in a loop to check the centering

  load/image skysub{n}
  write/out "Select an object common to all frames"
  center/gauss CURSOR ? ? ? 20,20  !box size is 20 x 20
  xc = {OUTPUTR(5)}
  yc = {OUTPUTR(6)}

elseif p4 .eq. "A" then

!Compute image shifts

  xapprox = {xc} - {xshift} + {work3,OFFSET_DEC,@{n}}/{scale}
  yapprox = {yc} - {yshift} + {work3,OFFSET_RA,@{n}}/{scale}
  xstart = {xapprox} - 10
  xend   = {xapprox} + 10
  ystart = {yapprox} - 10
  yend   = {yapprox} + 10
  write/table center :XSTART @1 {xstart}
  write/table center :XEND @1 {xend}
  write/table center :YSTART @1 {ystart}
  write/table center :YEND @1 {yend}
  center/gauss skysub{n},center ? ? ? 

  er = 1./{scale}
  
  if {OUTPUTR(7)} .ge. {er} .or. {OUTPUTR(8)} .ge. {er} then
	reload = "T"
  endif
  if {OUTPUTR(7)} .eq. 0.0 then
	reload = "T"
  endif
  if {OUTPUTR(12)}/{OUTPUTR(13)} .lt. 0.33 .or. {OUTPUTR(12)}/{OUTPUTR(13)} .gt. 3.o then
	reload = "T"
  endif
 
  if "{reload}" .eq. "T" then
     write/out WARNING: This field may not be accurately centered. 
     load/image skysub{n}
     write/out Please select the object with the mouse.
     center/gauss CURSOR ? ? ? 20,20  !box size is 20 x 20
  endif
     reload = "F"
     xshift = {xc}-{OUTPUTR(5)}
     yshift = {yc}-{OUTPUTR(6)}

elseif p4 .eq. "B" then

! Blind shifts - not recommended

  xshift = {xshift} + {work3,OFFSET_DEC,@{n}}/{scale}
  yshift = {yshift} + {work3,OFFSET_RA,@{n}}/{scale}

elseif p4 .eq. "Y" then

  load/image skysub{n}
  write/out "Select the same bright bright object in each frame"
  center/gauss CURSOR ? ? ? 20,20  !box size is 10 x 10
  xshift = {xc}-{OUTPUTR(5)}
  yshift = {yc}-{OUTPUTR(6)}

endif

!Shift images - except the first 

if {n} .ge. 2 then
 write/out "{n} xshift = {xshift}  yshift = {yshift}"

 if {xshift} .lt. {xshmin} then
  xshmin = {xshift}
 endif
 if {xshift} .gt. {xshmax} then
  xshmax = {xshift}
 endif
 if {yshift} .lt. {yshmin} then
  yshmin = {yshift}
 endif
 if {yshift} .gt. {yshmax} then
  yshmax = {yshift}
 endif

 create/image &a 2,512,512 1,1,1,1 ? -10000
 insert/image skysub{n} &a @129,@129

!There is strange feature in the MIDAS shift command. Zero shifts do not
!produce a resultant image.


 if M$ABS({xshift}) .lt. 0.5 .and. M$ABS({yshift}) .lt. 0.5 then
   copy/ii &a imsh{n}
 else
   shift/ima  &a imsh{n} {xshift},{yshift}
 endif

 del/image &a NO

 write/tab work3 :SHIFILENAME @{n} imsh{n}.bdf
 del/ima &a NO

else

create/image imsh{n} 2,512,512 1,1,1,1 ? -10000
insert/image skysub{n} imsh{n} @129,@129
write/tab work3 :SHIFILENAME @1 imsh{n}.bdf

endif

del/ima skysub{n} NO

enddo

   create/icat obj.cat work3,:SHIFILENAME
   medran = M$NINT({upper}/2)-2
   if {medran} .le. 0 then
     medran = 1
   endif
   write/out Median Combining with {medran}

   average/images {p5} = obj.cat ? ? median,{medran},{medran} -1000,20000

!Again the combine/ccd procedure may give better results.

endif

if p6 .eq. "Y" then 
 xshmin = 384 + M$NINT({xshmin})
 yshmin = 384 + M$NINT({yshmin}) 
 xshmax = 129 + M$NINT({xshmax})
 yshmax = 129 + M$NINT({yshmax})
 extract/image {p5} = {p5}[{xshmax},{yshmax}:{xshmin},{yshmin}]
else
 xshmax = 384 + M$NINT({xshmax})
 yshmax = 384 + M$NINT({yshmax}) 
 xshmin = 129 + M$NINT({xshmin})
 yshmin = 129 + M$NINT({yshmin})
 extract/image {p5} = {p5}[{xshmin},{yshmin}:{xshmax},{yshmax}]
endif

@s irac_acuts {p5}

!Write output file to a results table that should not be deleted

!You should create the table here if it does not exist.

exist = M$EXIST("reductions.tbl") 

if {exist} .eq. 0 then
   create/tab reductions 2 100 null
   create/col reductions :FILENAME c*20
   create/col reductions :COMBINED i*2
endif

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

last = {last} + 1

write/table reductions :FILENAME @{last} {p5}
write/table reductions :COMBINED @{last} 1
write/descr {p5} _EIO3_NAME/C/1/2 {lens}

!CLEAN UP ALL THE LOOSE FILES

del/icat obj NOCONF
del/icat sky NOCONF
del/tab work NO
del/tab work2 NO
del/tab work3 NO
del/tab sky NO
del/tab obj NO
del/tab center NO