! @(#)ccdrongain.prg	17.1.1.1 (ESO-IPG) 01/25/02 17:52:07
! @(#)ccdrongain.prg
! *************************************************************************
!.IDENTIFICATION: ccdrongain.prg
!
!.PURPOSE: Computes the conversion factor and read-out-noise using two bias
!          framed and two flat field frames
! 
!.AUTHOR: T. Augusteijn
! 
!.USE: @@ccdrongain p1 p2 p3 p4
!             where p1 = first bias frame
!                   p2 = second bias frame
!                   p3 = first flat field frame
!                   p4 = second flat field frame
! 
!.ASSUMPTIONS: All frames have the same size. For a reliable result the 
!              flat fields should be taken through the same filter.
! 
!.METHOD: Of each frame the outer 51 pixels are excluded. The remainder is
!         divided in 64 equal sized rectangles. To calculate the RON the
!         bias frames are subtracted, where the o
! 
!.RESULT:
!                       gain = computed gain (e-/ADU)
!                       egain = 
!                       OUTPUTR(2) = ron (e-)
! *************************************************************************
define/par p1 ? ima "Bias 1: "
define/par p2 ? ima "Bias 2: "
define/par p3 ? ima "Flat field 1: "
define/par p4 ? ima "Flat field 2: "

define/local sizex/i/1/1 0
define/local sizey/i/1/1 0
define/local xbegin/i/1/1 0
define/local xstart/i/1/1 0
define/local ybegin/i/1/1 0
define/local ystart/i/1/1 0
define/local xend/i/1/1 0
define/local xstep/i/1/1 0
define/local yend/i/1/1 0
define/local ystep/i/1/1 0
define/local xlow/r/1/1 0.0
define/local xhigh/r/1/1 0.0
define/local avebias/r/1/2 0.0 all
define/local aveff/d/1/2 0.0 all
define/local diffaveb/r/1/1 0.0
define/local rataveff/d/1/1 0.0
define/local n/i/1/1 0
define/local m/i/1/1 0
define/local nm/i/1/1 0
define/local gain/d/1/1 0.0
define/local egain/d/1/1 0.0
define/local ron/d/1/1 0.0
define/local eron/d/1/1 0.0

write/out " "
write/out "************************************************************************"
write/out "!                                                                      !"
write/out "!  This programme calculates the Gain and RON using 2 bias frames and  !"
write/out "!  2 flat fields. All frames should have the same size, and the flat   !"
write/out "!  fields should be taken through the same filter.                     !"
write/out "!                                                                      !"
write/out "************************************************************************"
write/out " "
write/out " "

! Create tables

create/tab bias 1 64 null
create/col bias :avebias r*4
create/col bias :stdbias r*4

create/tab ff 4 64 null
create/col ff :averatff r*4
create/col ff :stdratff r*4
create/col ff :aveff1 r*4
create/col ff :aveff2 r*4

! determine 100x100 area in the center of the image

sizex = '{p1},NPIX(1)' 
sizey = '{p1},NPIX(2)'

xbegin = '{p1},start(1)' + sizex/2 - 50
ybegin = '{p1},start(2)' + sizey/2 - 50

xend = xbegin + 100
yend = ybegin + 100

! Determine average of the bias frames using one sided 2 sigma clipping
! on the central 100x100 pixels and compare the values

write/out " The average values of the bias frames in the central 100x100 pixels:"
write/out " "

set/midas output=no

stat/ima {p1} [{xbegin},{ybegin}:{xend},{yend}]

xhigh = {outputr(8)} + 3*{outputr(4)}

stat/ima {p1} [{xbegin},{ybegin}:{xend},{yend}] ? 0,{xhigh}

avebias(1) = {outputr(3)}

set/midas output=yes

set/format f7.3 f7.3

write/out " Bias frame" {p1} " Ave:" {outputr(3)} "ADU  Sig:" {outputr(4)}

set/midas output=no

stat/ima {p2} [{xbegin},{ybegin}:{xend},{yend}]

xhigh = {outputr(8)} + 3*{outputr(4)}

stat/ima {p2} [{xbegin},{ybegin}:{xend},{yend}] ? 0,{xhigh}

avebias(2) = {outputr(3)}

set/midas output=yes

write/out " Bias frame" {p2} " Ave:" {outputr(3)} "ADU  Sig:" {outputr(4)}
write/out " "

diffaveb = m$abs(avebias(1)-avebias(2))

if {diffaveb} .gt. {outputr(4)} then

write/out "                 !! Error!! "
write/out " Difference average of bias frames ({diffaveb}) is larger than the "
write/out " standard deviation ({outputr(3)}) in one frame. "
write/out " "

return

endif

! Determine the average of the flat field frames in the central 100x100 pixels
! using two sided 2 sigma clipping, and see which is higher

set/midas output=no

stat/ima {p3} [{xbegin},{ybegin}:{xend},{yend}]

xlow = {outputr(8)} - 5*{outputr(4)}

xhigh = {outputr(8)} + 5*{outputr(4)}

stat/ima {p3} [{xbegin},{ybegin}:{xend},{yend}] ? {xlow},{xhigh}

aveff(1) = {outputr(3)} - {avebias(1)}

stat/ima {p4} [{xbegin},{ybegin}:{xend},{yend}]

xlow = {outputr(8)} - 5*{outputr(4)}

xhigh = {outputr(8)} + 5*{outputr(4)}

stat/ima {p4} [{xbegin},{ybegin}:{xend},{yend}] ? {xlow},{xhigh}

aveff(2) = {outputr(3)} - {avebias(1)}

rataveff = {aveff(1)} / {aveff(2)}

set/midas output=yes

if {rataveff} .gt. 100.0 then

  write/out "                 !! Error!! "
  write/out " The ratio of intensity in the two flat fields ({rataveff}) is larger "
  write/out " than 100. No reliable Gain can be derived. "
  write/out " "

  return

 else

  rataveff = 1.0 / rataveff

  if {rataveff} .gt. 100.0 then

   write/out "                 !! Error!!"
   write/out " The ratio of intensity in the two flat fields ({rataveff}) is larger "
   write/out " than 100. No reliable Gain can be derived. "
   write/out " "

   return

  endif

endif

! Define 64 areas in the image excluding the outer 51 pixels

xbegin = '{p1},start(1)' + 51
ybegin = '{p1},start(2)' + 51

xstep = (sizex - 102)/8
ystep = (sizey - 102)/8

! Compute difference of bias frame and determine average std in the 16 areas

write/out " "
write/out " The RON and Gain will be calculated now. This will take some time.... "

comp/ima &a = {p1} - {p2}

set/midas output=no

do m = 1 8 1

 xstart = xbegin + (m-1)*xstep
 xend   = xbegin + m*xstep

  do n = 1 8 1

   nm = (m-1)*8 + n

   ystart = ybegin + (n-1)*ystep
   yend = ybegin + n*ystep

   stat/ima &a [{xstart},{ystart}:{xend},{yend}]

   xhigh = {outputr(3)} + 3*{outputr(4)}

   xlow = {outputr(3)} - 3*{outputr(4)}

   stat/ima &a [{xstart},{ystart}:{xend},{yend}] ? {xlow},{xhigh}

   bias,:avebias,@{nm} = {outputr(3)}
   bias,:stdbias,@{nm} = {outputr(4)}

  enddo

enddo

stat/tab bias :stdbias

xhigh = {outputr(3)} + 2*{outputr(4)}

xlow = {outputr(3)} - 2*{outputr(4)}

sel/tab bias :stdbias.lt.{xhigh}.and.:stdbias.gt.{xlow}

stat/tab bias :stdbias

ron = outputr(3)/1.414214
eron = outputr(4)/1.414214

! Substract bias from flatfields

comp/ima &b = {p3} - {p1}
comp/ima &c = {p4} - {p1}

! Compute ratio of the flat fields such that the average is ~>1

if {aveff(1)} .gt. {aveff(2)} then

  comp/ima &a = 1000.0 * &b / &c

 else

  comp/ima &a = 1000.0 * &c / &b

endif

do m = 1 8 1

 xstart = xbegin + (m-1)*xstep
 xend   = xbegin + m*xstep

  do n = 1 8 1

   nm = (m-1)*8 + n

   ystart = ybegin + (n-1)*ystep
   yend = ybegin + n*ystep

   stat/ima &a [{xstart},{ystart}:{xend},{yend}]

   xhigh = {outputr(8)} + 10*{outputr(4)}

   xlow = {outputr(8)} - 10*{outputr(4)}

   stat/ima &a [{xstart},{ystart}:{xend},{yend}] ? {xlow},{xhigh}

   ff,:averatff,@{nm} = {outputr(3)}/1000.0
   ff,:stdratff,@{nm} = {outputr(4)}/1000.0

   stat/ima &b [{xstart},{ystart}:{xend},{yend}]

   xhigh = {outputr(8)} + 5*{outputr(4)}

   xlow = {outputr(8)} - 5*{outputr(4)}

   stat/ima &b [{xstart},{ystart}:{xend},{yend}] ? {xlow},{xhigh}

   ff,:aveff1,@{nm} = {outputr(3)}

   stat/ima &c [{xstart},{ystart}:{xend},{yend}]

   xhigh = {outputr(8)} + 5*{outputr(4)}

   xlow = {outputr(8)} - 5*{outputr(4)}

   stat/ima &c [{xstart},{ystart}:{xend},{yend}] ? {xlow},{xhigh}

   ff,:aveff2,@{nm} = {outputr(3)}

  enddo

enddo

! Compute gain with two way 2 sigma clipping

comp/tab ff :gain = (:averatff/:stdratff)*(:averatff/:stdratff)*(1/:aveff1+1/:aveff2)

stat/tab ff :gain

xlow = {outputr(3)} - 2*{outputr(4)}
xhigh = {outputr(3)} + 2*{outputr(4)}

sel/tab ff :gain.lt.{xhigh}.and.:gain.gt.{xlow}

stat/tab ff :gain

gain = outputr(3)
egain = outputr(4)

ron = ron * gain
eron = ron * m$sqrt(eron*eron*gain*gain/(ron*ron)+egain*egain/(gain*gain))

set/midas output=yes

set/format ,f7.3 ,f7.3

write/out " "
write/out "   ---------------------------------------------------------"
write/out "                                                            "
write/out "      The results are: Gain =  {gain} +/- {egain}  e-/ADU   "
write/out "                       RON  =  {ron} +/- {eron}  e-         "
write/out "                                                            "
write/out "   ---------------------------------------------------------"
write/out " "

delete/imag &a no
delete/imag &b no
delete/imag &c no

return