! %W% (ESO-IPG) %G% %U%
! @@ xcorall.prg	  
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! MIDAS procedure xcorall
! A. Kaufer	981130
! 
! cross-correlation for the FEROS Object-Calibration mode 
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 
!  version 8 : order by order Xcorr, order by order corrections, 
!              histogramm fitting
!
define/para p1 ? TAB "Enter table name: "
define/para p2 ? N "Enter frame number: "
define/para p3 ? N "Enter reference frame number: "
define/para p4 ? N "Enter ThAr frame number: "
define/para p5 ? N "Enter reference ThAr frame number: "
define/para p6 0.15 N "Enter histogram binsize : "
define/para p7 ? C "Enter action [create/enter]: "
!
!
define/local i/i/1/1      0
define/local j/i/1/1      0
define/local first/i/1/1  1
define/local last/i/1/1  39
define/local num/i/1/1    0
define/local temp/r/1/1 0.0
define/local lam/r/1/1  0.0
define/local dis/r/1/1  0.0
define/local sig/r/1/1  0.0
define/local dv/r/1/1   0.0
define/local lb/r/1/1   0.0
define/local lbmm/r/1/1 0.0
define/local bc/r/1/1   0.0
define/local bcmm/r/1/1 0.0
define/local c/r/1/1    299792.458
define/local win/i/1/1  50
define/local hpara/r/1/3 {p6},0,0

goto {p7}

create:
dele/tab {p1} NO
crea/tab {p1} 8 500
crea/col {p1} :num               I*2
crea/col {p1} :JD24        F14.6 R*8
crea/col {p1} :bc   "km/s" F8.3  R*4 
crea/col {p1} :dbc  "km/s" F8.3  R*4 
crea/col {p1} :ddbc "km/s" F8.3  R*4 
goto finish

enter:
copy/dk {p1}.tbl TBLCONTR/I/4/1 line/i/1/1 
create/row {p1} @{line} 1

set/format F14.6
line = line+1
{p1},:num,@{line}  = {p2}
{p1},:JD24,@{line} = {f{p2}ext1,O_TIME(4)}
{p1},:bc,@{line}   = {f{p2}ext1,BARY_CORR}

crea/tab median 8 100
crea/col median :num
crea/col median :order
crea/col median :obj
crea/col median :thar
crea/col median :calib
crea/col median :thar2
crea/col median :dotcb

num = 0
do i = {first} {last}
 num = num+1
 median,:order,@{num} = i
 median,:num,@{num} = {p2}
!
 write/out "---------------------"
 write/out "Processing order {i} "
 write/out "---------------------"
! get dispersion of order for fiber 1
 j = 1+(i-1)*5
 lam = {ThAr{p4}lines1.tbl,dcoef({j})}
 j = j+1
 dis = {ThAr{p4}lines1.tbl,dcoef({j})}
 dv = c*dis*0.015/lam
! write/out "dv/pixel (fiber 1) = {dv} km/s "

! create gaussian acc. to blaze for order shaping
  j = i+1
! length of the order
  if {i} .lt. 38 then 
   lb   = {BLAZE.tbl,:wave,{j}}-{BLAZE.tbl,:wave,{i}}
   lbmm  = lb/dis
   sig = lbmm/5
   crea/ima gauss 1,4102 -30.765,0.015 GAUSS 0,{sig}
  endif

! obj-obj
 extr/ima &frame = f{p2}ext1[<,@{i}:>,@{i}]
 extr/ima &templ = f{p3}ext1[<,@{i}:>,@{i}] 
 write/out "object*objectref fiber 1 : f{p2}ext1 * f{p3}ext1"
! ***
! comp/ima &frame = &frame*gauss
! stat/ima &frame >null
! comp/ima &frame = &frame/{outputr(2)}
! comp/ima &templ = &templ*gauss
! stat/ima &templ >null
! comp/ima &templ = &templ/{outputr(2)}
! ***
 stat/ima &frame [-1:1] >null
 comp/ima &frame = {outputr(2)}-&frame
 comp/ima &frame = &frame*gauss
 stat/ima &templ [-1:1] >null
 comp/ima &templ = {outputr(2)}-&templ
 comp/ima &templ = &templ*gauss
! ***
! plo &frame
! plo &templ
 xcorr  &templ &frame xcorr {win} 
 center/gauss xcorr >null
 median,:obj,@{num} = outputr(5)*dis/lam*c
! plo xcorr

! thar - thar
 extr/ima &frame = f{p4}ext1[<,@{i}:>,@{i}]
 extr/ima &templ = f{p5}ext1[<,@{i}:>,@{i}] 
 write/out "thar  *tharref   fiber 1 : f{p4}ext1 * f{p5}ext1"
 comp/ima &frame = &frame*gauss
 stat/ima &frame >null
 comp/ima &frame = &frame/{outputr(2)}
 comp/ima &templ = &templ*gauss
 stat/ima &templ >null
 comp/ima &templ = &templ/{outputr(2)}
 xcorr  &templ &frame xcorr {win} 
 center/gauss xcorr >null
 median,:thar,@{num} = outputr(5)*dis/lam*c
! plo xcorr

! get dispersion of order for fiber 2
! j = 1+(i-1)*5
! lam = {ThAr{p4}lines2.tbl,dcoef({j})}
! j = j+1
! dis = {ThAr{p4}lines2.tbl,dcoef({j})}
! dv = c*dis*0.015/lam
! write/out "dv/pixel (fiber 2) = {dv} km/s "

! obj-thar  calib fiber
 extr/ima &frame = f{p2}ext2[<,@{i}:>,@{i}]
 extr/ima &templ = f{p4}ext2[<,@{i}:>,@{i}] 
 write/out "object*thar      fiber 2 : f{p2}ext2 * f{p4}ext2"
 comp/ima &frame = &frame*gauss
 stat/ima &frame >null
 comp/ima &frame = &frame/{outputr(2)}
 comp/ima &templ = &templ*gauss
 stat/ima &templ >null
 comp/ima &templ = &templ/{outputr(2)}
 xcorr  &templ &frame xcorr {win} 
 center/gauss xcorr >null
 median,:calib,@{num} = outputr(5)*dis/lam*c
! plo xcorr

! thar - thar fib 2
! extr/ima &frame = f{p4}ext2[<,@{i}:>,@{i}]
! extr/ima &templ = f{p5}ext2[<,@{i}:>,@{i}] 
! write/out "thar  *tharref   fiber 2 : f{p4}ext2 * f{p5}ext2"
! comp/ima &frame = &frame*gauss
! stat/ima &frame >null
! comp/ima &frame = &frame/{outputr(2)}
! comp/ima &templ = &templ*gauss
! stat/ima &templ >null
! comp/ima &templ = &templ/{outputr(2)}
! xcorr  &templ &frame xcorr {win} 
! center/gauss xcorr >null
! median,:thar2,@{num} = outputr(5)*dis/lam*c
enddo 
!
! correction order by order:  object - zeropoint - drift  - barycentric
!
comp/tab median :dotcb = :obj-:thar-:calib-{{p1},:bc,@{line}}
!
! find center of histogram
!
comp/histo hist = median :dotcb {hpara(1)} 
stat/ima   hist 
hpara(2) = {hist,START(1)}+outputi(3)*{hist,STEP(1)}-5.
hpara(3) = {hist,START(1)}+outputi(3)*{hist,STEP(1)}+5.
!
! center histogram with gaussian
!
comp/histo hist = median :dotcb {hpara(1)} {hpara(2)} {hpara(3)} 
set/gra xa ya
plo hist
center/gauss hist
if  outputr(14) .lt. 5 then
 write/out "----------------------------------"
 write/out " {p2}, {p4}  rejected, sorry ! " 
 write/out "----------------------------------"
 goto finish
endif
!
! write/results to output table
!
{p1},:dbc,@{line}  = outputr(5)
{p1},:ddbc,@{line} = outputr(12)/outputr(14)
!
sele/tab median all
copy/tab median median{p2}
dele/tab median NO
!
read/tab {p1} @{line}
goto finish

finish:
set/gra color=1
sele/tab {p1} all

end: