! @(#)necdispres.prg	17.1.1.1  01/25/02 17:51:08
!+++ NECDISPRES +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.COPYRIGHT     (C) 1998  European Southern Observatory
!
!.IDENT         necdispres.prg
!
!.KEYWORDS      spectroscopy, echelle, resolution
!
!.USAGE         @s (IDENT) lintab lincat wlc param yrange disp title subt
!
!.INPUT         lintab  (P1): name of a line table
!				product of IDENT/ECHELLE
!               lincat  (P2): line reference catalog (thar.tbl)
!               wlc     (P3): wavelength calibration frame
!               param   (P4): parameters used for extraction, resampling
!				and measuring the line width:
!				slit,sample,detect,ordsta,ordend
!				slit	window size for the order extraction
!				sample	step size for the resampling
!				detect	detection window size to find line
!					centers (CENTER/GAUSS)
!				ordsta	first order to use
!				ordend	last order to use
!               yrange  (P5): plot ranges in y (4 values)
!				min1,max1,min2,max2
!				for the deltaX and resolution, respectively
!               disp    (P6): display flag,free spectral range flag:
!				YES: show result on graphics window
!				NO:  creates no plot.
!				YES: show result only in free spectral range
!				NO:  take the entire CCD.
!
!.OUTPUT        postscript.ps, resolution.dat, free_sp_rg.dat.
!
!.RETURN        
!
!.DEFAULT	P1	{LINTAB}
!		P2	{LINCAT}
!		P3	{WLC}
!		P4	{SLIT},{SAMPLE},30,{ECHORD(4)},{ECHORD(5)}
!		P5	+
!		P6	Y,N
!		P7	+
!		P8	+
!
!.PURPOSE       MIDAS procedure to display Echelle Spectra line's FWHMs
!		as a function of pixel number in X and Y as well as the
!		resolution (lambda/dlambda).
!
!		Creates 6 graphs and the postscript file with these 6 plots :
!		. plots column :DX (FWHM) vs column :X_1 and :YNEW_1,
!		. plots column :RESOL vs column :X_1 and :YNEW_1,
!		. plots column :RESOL vs column :IDENT_1 and :X_1, 
!
!		Computes the Median, Average and Standard Deviation of DX
!		(FWHM) and RESOL.
!
!		Finally, it creates an ASCII file with columns:
!		  :X_1,:YNEW_1,:IDENT_1,:ORDER_1,:DX,:DLAMBDA,:RESOL
!
!		necdispres,FREE_SP_RG computes the central wavelength
!		corresponding to the central column of the CCD, the wavelength
!		range of the free spectral range, selects the values from table
!		contained within the the FREE SPECTRAL RANGE before ploting
!		and computation of statiscal values. 
!		It also creates an ASCII file with columns: 
!		  :X :YNEW_1 :Y_1 :ORDER_1 :lamb_ct :lamb_st :lamb_en .
!
!
!.ALGORITHM     Measures the line width of the lines stored in the line table
!		(P1) by means of the command CENTER/GAUSS. These measurements
!		will be made in pixel-space as well as in the wavelength
!		calibrated space of the arlamp exposure (P3).
!		The result will be stored in a ASCII table and as a postscript
!		file.
!		The title may be configured using the parameters 
!
!		At the end the result may be displayed controled via the
!		parameter P4.
!
!.ENVIRON       ECHELLE context
!   
!.LANGUAGE      MIDAS 98NOV
!
!.AUTHOR        P.Ballester, O.Boitquin, S.Wolf ESO/DMD
!
!.VERSION       0.5 1999/05/10  Created (PB,OB).
!		    1999/05/20	User interface improved (SW).
!		    1999/05/21	Possibility to pass title values.
!		0.6 1999/05/31  necdispres,FREE_SP_RG subroutine created (OB). 
!		    1999/06/01  option p6 extended to get the free spectral
!				range.
!
!.COMMENT
!------------------------------------------------------------------------------

define/param  P1 {LINTAB}      TAB  "Line Table:"
define/param  P2 {LINCAT}     CHAR  "Line catalog" 
define/param  P3 {WLC}         IMA  "WLC frame:" 
define/param  P4 {SLIT},{SAMPLE},30,{ECHORD(4)},{ECHORD(5)}  NUMB "Slit, sample, line width, orders to merge (start, stop)"
define/param  P5 +             NUMB  "Y-ranges for the plots:"
define/param  P6 Y,N           CHAR  "Display result and Free Spectral range: Yes or No"
define/param  P7 + c "TITLE: central-wavel.,slit(arcsec),CCD-temp,CCD-name:"
define/param  p8 + c "TITLE: Enter file name,Observing date:"
!define/param  p9 N             CHAR  "Display only the free spectral range: Yes or No"
!
! local variables
!
define/local CNT/I/1/1  1
define/local K/I/1/1  0
define/local catrord/C/1/25 catrord
define/local disp/C/1/3 " " ALL
define/local sprg/C/1/3 " " ALL
define/local SPEC/C/1/25 " " ALL
define/local KORD/I/1/1 0
define/local od1/i/1/1 0
define/local od2/i/1/1 0
define/local slt/i/1/1 0 
define/local sam/r/1/1 0
define/local wid/r/1/1 0
!
!
!
! Get the parameters for the order extraction
!
write/key outputr {SLIT},{SAMPLE},30,{ECHORD(4)},{ECHORD(5)}
write/key outputr {p4}				!Insert user values
slt = outputr(1)				!Remember as local variables
sam = outputr(2)				!sampling
wid = outputr(3)				!line width
od1 = outputr(4)				!order number: start
od2 = outputr(5)				!order number: end
!
! Display info
!
set/format I1 F7.5
write/out "Determine the resolution using:"
write/out "  Line table:            {p1}"
write/out "  Line catalog:          {p2}"
write/out "  Wavel.calib. frame:    {p3}"
write/out "  Slit width (extract.): {slt}"
write/out "  Sampling:              {sam}"
write/out "  Use orders:            [{od1},{od2}]"
write/out "  Line detection window: {wid}"
if p5(1:1) .ne. "+" write/out "  Y plot ranges:         {p5}"
!
! Do the order extraction and resampling of the arc lamp exposure
!
save/eche dispr >Null
extract/echelle   {p3}  &w {slt},0 AVERAGE
rebin/echelle &w &r {sam}
set/format I1
!
! Merge the orders while storing orders separately
!
merge/echelle &r rord {od1},{od2} NOAPPEND >Null
creat/icat &cr rord*.bdf >Null
!
! Joining tables: measurments and line catalog
!
select/table {P2} all >Null
select/table {P1} all >Null
join/table {P1} :IDENT {P2} :WAVE  &l 0.001
!
! Measure the line width in pixel space using CENTER/GAUSS
!
write/out "Measuring FWHM in pixel space" 
compute/table     &l :XSTART = :X_1 - {wid}/2.
compute/table     &l :XEND   = :X_1 + {wid}/2.
compute/table     &l :YSTART = :Y_1
compute/table     &l :YEND   = :Y_1

center/gauss      middummw.bdf,middumml.tbl middumml.tbl >Null
select/table      &l :IDENT_1.gt.0. >Null
compute/table     &l :DX = :XFWHM
statistic/table   &l :DX

compute/table     &l :X = :X_1
regress/poly      &l :AUX_1 :X,:YNEW_1  {DC},{DC} KEYLONG >Null
save/regress      &l WAVMAP KEYLONG >Null
compute/table     &l :X = :XSTART
compute/regress   &l :WAV_ST = WAVMAP >Null
compute/table     &l :XSTART = :WAV_ST/:ORDER_1
compute/table     &l :X = :XEND
compute/regress   &l :WAV_ED = WAVMAP >Null
compute/table     &l :XEND = :WAV_ED/:ORDER_1
!
!        
!
$$ ls rord*.bdf > lsrord.lst
$$ echo "DEFINE/FIELD 1 56 C :NAME" > fmtrord.fmt
CREATE/TABL &ctor * * lsrord.lst fmtrord.fmt
create/colu &ctor :SEQ I
comput/tabl &ctor :SEQ = seq
create/colu &ctor :ABSORD I
comput/tabl &ctor :ABSORD = {echord(6)} - seq

create/colu &ctor :WLSTART R
execut/tabl middummctor &ctor,:WLSTART,@[:SEQ] = m$value([:name],start)

write/out "Measuring FWHM in wavelength space" 
do k = {od1} {od2} 1 
  kord = m$value(&ctor,:ABSORD,@{k})
  select/table &l :ORDER_1.eq.{kord} >Null
  spec = m$value(&ctor,:NAME,@{k})
  center/gauss {spec},middumml.tbl middumml.tbl >Null
enddo

select/table      &l :IDENT_1.gt.0. >Null
compute/table     &l :DLAMBDA = :XFWHM  
compute/table     &l :RESOL   = :IDENT_1/:DLAMBDA
!statistic/table   &l :RESOL

select/table      &l :STATUS.eq.0..AND.:RESOL.gt.0. >Null

k = m$parse(p6,"flag")
if k .gt. 1 then
  disp = "{flag01}"
  sprg = "{flag02}"
else if k .eq. 1 then
  disp = "{flag01}"
  sprg = "NO"
endif
disp = m$upper(disp)
sprg = m$upper(sprg)

if sprg(1:1) .eq. "N" then
  copy/table   &l  &resol             !No selection of the free spectral range
else
  @s necdispres,FREE_SP_RG &l &resol  !Selection of the free spectral range
endif

!2. Plotting

set/graph

define/local STEP/R/1/1 1.
define/local LABEL/C/1/20  " " 
define/local offset_graph/R/1/1 1.
define/local y_min_graph/R/1/1 1.
define/local y_max_graph/R/1/1 1.

define/local y_scale_graph/R/1/1 1.
define/local plot_check/I/1/1 0

set/format I1

set/graph 
offset_graph = 15
set/graph font=1 yscale=-30 xoffset=20 yoffset={offset_graph}

!PLOT/AXES 0,4096 2,10 10,.3,1,-1


! :FWHM is the x-FWHM of spectral lines. It can be measured during
! SEARCH/ECHELLE which finds the lines.
! :RESOL is the spectral resolution (Dlambda/lambda = :FWHM*Disp/:IDENT)
! The dispersion can only be known after wavelength calibration 
!

select/table &resol :IDENT_1.GT.0..AND.:RESOL.LT.500000. 0. >Null
set/graph
if disp(1:1) .eq. "Y" @s necdispres,PLOT &resol {p3} {p5} {p6} {p7} {p8}
!
! store the data in an ASCII table
!
assign/print FILE resolution.dat
print/table &resol :X_1 :YNEW_1 :IDENT_1 :ORDER_1 :DX :DLAMBDA :RESOL
!
! reset the previous echelle settings
!
initial/eche dispr
-delete dispr*.tbl
-delete rord*.bdf


!=== SUBROUTINES =============================================================

ENTRY PLOT

define/param p1 &resol c "Enter name of resolution table:"
define/param p2 {wlc} c "Enter name of wavelength calibration frame WLC:"
define/param p3 + c "Enter y-ranges for the plots:"
define/param p4 Y c "Display result on graphics window: Yes or No"
define/param p5 + c "TITLE: central-wavel.,slit(arcsec),CCD-temp,CCD-name:"
define/param p6 + c "TITLE: Enter file name,Observing date:"

define/local K/I/1/1  0
!define/local colour/I/1/1  0
define/local Ncolour/I/1/1 4
define/local colsel/I/1/4 6,4,3,2
define/local typsel/I/1/4 5,3,2,8
define/local NX/I/1/1 0
define/local NY/I/1/1 0
define/local UPY/I/1/1 0
define/local LWY/I/1/1 0
define/local MINident/R/1/1  0.
define/local MAXident/R/1/1  0.
!define/local CENident/R/1/1  0.
define/local MNDX/I/1/1 1
define/local MXDX/I/1/1 6
define/local MNRESOL/R/1/1 30000.0
define/local MXRESOL/R/1/1 150000.0
!
! statistic
!
define/local MeanDX/R/1/1  0.
define/local StdDX/R/1/1  0.
define/local N2DX/I/1/1  0.
define/local MeanRES/R/1/1  0.
define/local StdRES/R/1/1  0.
define/local N2RES/I/1/1  0.
define/local MedDX/R/1/1  0.
define/local MedRES/R/1/1  0.
!
! title
!
define/local hfname/c/1/100 "Filename: "
define/local fname/c/1/100 " " ALL
define/local cname/c/1/100 " " ALL
define/local datobs/c/1/100 " " ALL
define/local swid/c/1/10 " " ALL
define/local temp/c/1/10 " " ALL
define/local wlen/c/1/10 " " ALL
define/local maintit/i/1/1 0
define/local subtit/i/1/1 0
!
! configure title
!
if p5(1:1) .ne. "+" then
  wlen = "?"
  swid = "?"
  temp = "?"
  cname = "???"
  k = m$parse(p5,"tit")
  set/format I2
  if k .ge. 1 wlen = "{tit01}"
  if k .ge. 2 slit = "{tit02}"
  if k .ge. 3 temp = "{tit03}"
  if k .ge. 4 cname = "{tit04}"
  maintit = 1
endif
!
! configure subtitle
!
if p6(1:1) .ne. "+" then
  fname = "?"
  datobs = "?"
  k = m$parse(p6,"sub")
  set/format I2
  if k .ge. 1 fname = "{sub01}"
  if k .ge. 2 datobs = "{sub02}"
  subtit = 1
endif


if p3(1:1) .ne. "+" then
  define/local yr/r/1/4 {MNDX},{MXDX},{MNRESOL},{MXRESOL}	!set defaults
  write/out "Y-ranges [default]: [{yr(1)},{yr(2)}]; [{yr(3)},{yr(4)}]"
  write/key yr/r/1/4 {p3}					!overwrite by
  write/out "Y-ranges [user]:    [{yr(1)},{yr(2)}]; [{yr(3)},{yr(4)}]"
  MNDX = yr(1)							!the user sett.
  MXDX = yr(2)
  MNRESOL = yr(3)
  MXRESOL = yr(4)
endif
!
! get the dimension of the WLC-image
!
NX = m$value({p2},NPIX(1))
NY = m$value({p2},NPIX(2))
!
! Determine the statistic of the table
!
select/table {P1} all >Null
statis/table {P1} :DX >Null
MeanDX = {outputr(3)}
StdDX = {outputr(4)}
N2DX = {outputi(2)}/2
statis/table {P1} :RESOL >Null
MeanRES = {outputr(3)}
StdRES = {outputr(4)}
N2RES = {outputi(2)}/2
sort/tabl {P1} :DX
MedDX = m$value({P1},:DX,@{N2DX})
sort/tabl {P1} :RESOL
MedRES = m$value({P1},:RESOL,@{N2RES})
!
!
! Configure the plot
!
set/grap
assign/grap POSTSCRIPT.l
set/grap font=1 !yaxis=auto !pmode=0
!
! 1.plot
!
select/table {P1} all >Null !:Y_1.EQ.{ECHORD(4)}
statis/table {P1} :IDENT_1 >Null
MINident = {outputr(1)} - 50.
MAXident = {outputr(2)} + 50.
!CENident = {outputr(3)}
plot/axes  0,{NX} {MNDX},{MXDX} -65,-65,13,100 "\+ \+ X  (pixels)" "\+ \+ \Delta X (line FWHM in pixels)"

UPY = 0
!DO K = {ECHORD(4)} {ECHORD(5)} 1  	!order selection
do K = 1 {Ncolour} 1
  LWY = {UPY}
  UPY = {NY}/{Ncolour}*{K}
! SELECT/TABLE {P1} :Y_1.EQ.{K}   	!order selection
  select/table {P1} :YNEW_1.GE.{LWY}.AND.:YNEW_1.LT.{UPY} >Null
! colour = {K}/{Ncolour}  		!order selection
  set/grap colour={colsel({K})}
  overpl/table {P1}   :X_1 :DX -65,-65,13,100 {typsel({K})}
enddo
!
! 2.plot
!
set/grap colour=1
select/table {P1} all >Null
overplot/axes  0,{NY} {MNDX},{MXDX} -65,-65,13,15 "\+ \+ Y  (pixels)" "\+ \+ \Delta X (line FWHM in pixels)"
UPY = 0
do K = 1 {Ncolour} 1
  LWY = {UPY}
  UPY = {NY}/{Ncolour}*{K} 
  select/table {P1} :YNEW_1.GE.{LWY}.AND.:YNEW_1.LT.{UPY} >Null
  set/grap colour={colsel({K})}
  overpl/table {P1}   :YNEW_1 :DX -65,-65,13,15 {typsel({K})}
enddo
!
! 3.plot
!
set/grap colour=1
select/table {P1} all >Null
overplot/axes  0,{NX} {MNRESOL},{MXRESOL} -65,-65,104,100 "\+ \+ X  (pixels)" "\+ \+ Resolution : \lambda / \Delta \lambda"
UPY = 0
do K = 1 {Ncolour} 1
  LWY = {UPY}
  UPY = {NY}/{Ncolour}*{K} 
  select/table {P1} :YNEW_1.GE.{LWY}.AND.:YNEW_1.LT.{UPY} >Null
  set/grap colour={colsel({K})}
  overpl/table {P1}   :X_1 :RESOL -65,-65,104,100 {typsel({K})}
enddo
!
! 4.plot
!
set/grap colour=1
select/table {P1} all >Null
overplot/axes 0,{NY} {MNRESOL},{MXRESOL} -65,-65,104,15 "\+ \+ Y  (pixels)" "\+ \+ Resolution : \lambda / \Delta \lambda"

UPY = 0
do K = 1 {Ncolour} 1
  LWY = {UPY}
  UPY = {NY}/{Ncolour}*{K} 
  select/table {P1} :YNEW_1.GE.{LWY}.AND.:YNEW_1.LT.{UPY} >Null
  set/grap colour={colsel({K})}
  overpl/table {P1}   :YNEW_1 :RESOL -65,-65,104,15 {typsel({K})}
enddo
!
! 5.plot
!
set/grap colour=1
select/table {P1} all >Null
overplot/axes {MINident},{MAXident} {MNRESOL},{MXRESOL} -65,-65,194,100 "\+ \+ \lambda  (\AA )" "\+ \+ Resolution : \lambda / \Delta \lambda"

UPY = 0
do K = 1 {Ncolour} 1
  LWY = {UPY}
  UPY = {NY}/{Ncolour}*{K} 
  select/table {P1} :YNEW_1.GE.{LWY}.AND.:YNEW_1.LT.{UPY} >Null
  set/grap colour={colsel({K})}
  overpl/table {P1}   :IDENT_1  :RESOL -65,-65,194,100 {typsel({K})}
enddo
!
! 6.plot
!
set/grap colour=1
select/table {P1} all >Null
overplot/axes 0,{NX} 0,{NY} -65,-65,194,15 "\+ \+ X  (pixels)" "\+ \+ Y  (pixels)"

UPY = 0
do K = 1 {Ncolour} 1
  LWY = {UPY}
  UPY = {NY}/{Ncolour}*{K} 
  select/table {P1} :YNEW_1.GE.{LWY}.AND.:YNEW_1.LT.{UPY} >Null
  set/grap colour={colsel({K})}
  overpl/table {P1}   :X_1 :YNEW_1 -65,-65,194,15 {typsel({K})}
enddo
!
!
! Plot the title
!
set/gra pmode=-1 colour=1
overpl/axes >Null

set/grap tsize=1.3					!main title
if maintit .eq. 1 then
  label/grap "Central wavelength: {wlen} nm, slit: {swid} arcsec, CCD: {cname}, {temp} \circ C" 90,167,mm
else
  label/grap "Resolution plot" 90,167,mm
endif

set/grap tsize=1					!sub title
if subtit .eq. 1 then
  outputi(1) = m$len(datobs)
  label/grap "{hfname} {fname}, Observ. date: {datobs(1:{outputi(1)})}" 90,162,mm
else
  label/grap "Wavelength calibration file used: {wlc}" 90,162,mm
endif

set/grap tsize=1.0					!sub title: statistic
SET/FORM f5.3
label/grap "Median(DX) = {MedDX} pix , Mean(DX) = {MeanDX} pix , \sigma (DX) = {StdDX} pix" 90,157,mm
label/grap "Median(RES) = {MedRES} , Mean(RES) = {MeanRES} , \sigma (RES) = {StdRES}" 90,152,mm
!
!
! Check if plot has to be displayed on a graphics window:
!
if p4(1:1) .eq. "Y" then
  set/format I1
  show/disp >Null			!check if graph.window already exists
  if outputi(12) .lt. 0 then
    create/grap 1 1000,707
    outputi(12) = 1
  endif
  copy/grap  g,{outputi(12)} >Null
endif
!..............................................................................

ENTRY FREE_SP_RG

define/param p1 &l c "Enter name of input table:"
define/param p2 &resol c "Enter name of output resolution table:"
define/param p3 {wlc} c "Enter name of wavelength calibration frame WLC:"

!define/local od1/i/1/1 0
!define/local od2/i/1/1 0
define/local K/I/1/1  0
define/local Kini/I/1/1  0
define/local NX2/I/1/1 0
define/local wave_ct/R/1/1  0.
define/local wave_st/R/1/1  0.
define/local wave_en/R/1/1  0.

! Central column of the CCD
NX2 = m$value({p3},NPIX(1)) /2

! Compute the Yfit from the order table 
copy/table {ORDTAB} &yfit
select/table &yfit all
DELETE/ROW  &yfit @1..{outputi(1)} 
create/row &yfit @0 {echord(1)}
comput/table &yfit :ORDER = seq
comput/table &yfit :X = {NX2}
compute/regress &yfit :yfit = COEFF

! Compute the central wavelength and the free spectral range
copy/table   {p1}  &free
select/table &free all
DELETE/ROW  &free @1..{outputi(1)}
create/row &free @0 {echord(1)}
compute/table &free :Y_1 = seq
compute/table &free :X = {NX2}
copy/tt &yfit :yfit &free :ynew_1
compute/regress &free :lambda = WAVMAP
compute/table &free :ORDER_1 = {echord(6)} - :Y_1
compute/table &free :lamb_ct = :lambda/:ORDER_1
compute/table &free :lamb_st = :lamb_ct - :lamb_ct/(2*:ORDER_1)
compute/table &free :lamb_en = :lamb_ct + :lamb_ct/(2*:ORDER_1)
 
! Select context for the selection of the free spectral range
!wave_ct = m$value(&free,:lamb_ct,@{ECHORD(4)})
wave_st = m$value(&free,:lamb_st,@{ECHORD(4)})
wave_en = m$value(&free,:lamb_en,@{ECHORD(4)})
SELECT/TABLE {P1} SELECT.AND.:Y_1.EQ.{ECHORD(4)}.AND.:IDENT_1.gt.{wave_st}.AND.:IDENT_1.lt.{wave_en} >Null

! Selection of the free spectral range 
Kini = {ECHORD(4)} +1
do k = {Kini} {ECHORD(5)} 1 
!  wave_ct = m$value(&free,:lamb_ct,@{k})
  wave_st = m$value(&free,:lamb_st,@{k})
  wave_en = m$value(&free,:lamb_en,@{k})
  SELECT/TABLE {P1} SELECT.OR.(:Y_1.EQ.{K}.AND.:IDENT_1.gt.{wave_st}.AND.:IDENT_1.lt.{wave_en}) >Null
enddo   

copy/table {p1} {p2}
!copy/table {p1} &TEMP1
!copy/table &TEMP1 {P1}

!
! store the data in an ASCII table
!
assign/print FILE free_sp_rg.dat
print/table &free :X :YNEW_1 :Y_1 :ORDER_1 :lamb_ct :lamb_st :lamb_en
!

!..............................................................................