# Assumptions:
# badpix.dat object.fits thar.fits flat.fits bias.fits all exist
# We know the spatial binning, in this case 1x1
# We know the trimsection for the data, in this case [1:2048,1:4608]
# We know whether this is single- and dual-aperture data
# We know the number of orders, and which ones, we are going to extract
# We know which slice we are going to trace on
# We know how wide the orders are in pixels

observatory set gemini-south
# Setting package information to echelle defautls; exectue ":q", ":q"
setinstrument echelle
# Set observ=gemini-south dispaxis=2 databas=database
epar echelle

# Make sure CCDSEC, DETSEC and DATASEC all agree
hcopy input1=object.fits[1] input2=object.fits[1] keywor1=DATASEC keywor2=CCDSEC keytype=string
hcopy input1=object.fits[1] input2=object.fits[1] keywor1=DATASEC keywor2=DETSEC keytype=string
hcopy input1=flat.fits[1] input2=flat.fits[1] keywor1=DATASEC keywor2=CCDSEC keytype=string
hcopy input1=flat.fits[1] input2=flat.fits[1] keywor1=DATASEC keywor2=DETSEC keytype=string
hcopy input1=thar.fits[1] input2=thar.fits[1] keywor1=DATASEC keywor2=CCDSEC keytype=string
hcopy input1=thar.fits[1] input2=thar.fits[1] keywor1=DATASEC keywor2=DETSEC keytype=string
hcopy input1=bias.fits[0] input2=bias.fits[0] keywor1=DATASEC keywor2=CCDSEC keytype=string
hcopy input1=bias.fits[0] input2=bias.fits[0] keywor1=DATASEC keywor2=DETSEC keytype=string

# Bias subtraction
!/bin/ls object.fits thar.fits flat.fits | sed 's/.fits/.fits[1]/' > in1
!/bin/ls object.fits thar.fits flat.fits | sed 's/.fits/_b.fits/' > in2
ccdproc images=@in1 output=@in2 ccdtype=' ' fixpix+ oversca+ trim+ zerocor+ darkcor- flatcor- readaxis=line fixfile=badpix.dat biassec=image trimsec=[1:2048,1:4608] zero=bias.fits interac+

# Copy over relevant header keywords
hcopy input1=object.fits[0] input2=object_b.fits keywor1=OBSEPOCH keywor2=' ' keytype=real
hcopy input1=object.fits[0] input2=object_b.fits keywor1=OBSEPOCH keywor2=' ' keytype=real
hcopy input1=object.fits[0] input2=object_b.fits keywor1=DATE-OBS keywor2=' ' keytype=string
hcopy input1=object.fits[0] input2=object_b.fits keywor1=UTSTART keywor2=' ' keytype=string
hcopy input1=object.fits[0] input2=object_b.fits keywor1=UTSTART keywor2=' ' keytype=string
hedit images=object_b.fits fields=LAT_OBS value=-30.228333 add+ addonly- delete- verify- show+ update+
hedit images=object_b.fits fields=LONG_OBS value=70.723333 add+ addonly- delete- verify- show+ update+
hedit images=object_b.fits fields=ALT_OBS value=2737.0 add+ addonly- delete- verify- show+ update+
hcopy input1=object.fits[0] input2=object_b.fits keywor1=RA keywor2=' ' keytype=string
hcopy input1=object.fits[0] input2=object_b.fits keywor1=RA keywor2=' ' keytype=string
hcopy input1=object.fits[0] input2=object_b.fits keywor1=DEC keywor2=' ' keytype=string
hcopy input1=object.fits[0] input2=object_b.fits keywor1=DEC keywor2=' ' keytype=string
hcopy input1=object.fits[0] input2=object_b.fits keywor1=EQUINOX keywor2=' ' keytype=real
hcopy input1=object.fits[0] input2=object_b.fits keywor1=EXPTIME keywor2=' ' keytype=real

# Define flats for each sub-aperture
!ln -s flat_b.fits trace.fits
!ln -s flat_b.fits order_trace.fits

# Define an aperture covering the slice with the most coverage for
#  tracing each order and do the trace
apfind trace.fits apertur=9 referen=' ' interac+ find+ recente+ resize+ edit+ nsum=30 nfind=9 minsep=200.0 maxsep=400.0 order=increasing
aptrace trace.fits interac+ find- recente- resize- edit- trace+ fittrac+ nsum=30 step=30 nlost=10 functio=chebyshev order=4 niterat=4

# Redefine the apertures for tracing the entire orders and apply the
#  above trace. Note that you must try to include an extra two or three
#  pixels either side of the order.
apfind order_trace.fits referen=trace.fits apertur=9 interac+ find- recente- resize- edit+ nsum=30 nfind=9 order=increasing

# Flatten the flatfield.
apflatten flat_b.fits flat_bn.fits referen=order_trace.fits interac+ find- recente- resize- edit- trace- fittrac- nsum=30 thresho=20.0 readnoi=rdnoise gain=gain functio=chebyshev order=11 niterat=3

!ln -s flat_b.fits flat_b_01.fits
!ln -s flat_b.fits flat_b_02.fits
!ln -s flat_b.fits flat_b_03.fits
!ln -s flat_b.fits flat_b_04.fits
!ln -s flat_b.fits flat_b_05.fits
!ln -s flat_b.fits flat_b_06.fits
!ln -s flat_b.fits flat_b_07.fits
!ln -s flat_b.fits flat_b_08.fits
!ln -s flat_b.fits flat_b_09.fits
!ln -s flat_b.fits flat_b_10.fits
!ln -s flat_b.fits flat_b_11.fits
!ln -s flat_b.fits flat_b_12.fits
!ln -s flat_b.fits flat_b_13.fits
!ln -s flat_b.fits flat_b_14.fits

# Create shifted traces for sub-apertures: TO BE EXECUTED IN SHELL, NOT CL
cd database; ../trace_create.pl trace flat_b 5 7 201 1 ; cd ..

# Flatfield the object and thar files
!echo "object_b.fits\nthar_b.fits" > in1
!sed 's/_b/_bn/' in1 > in2
ccdproc images=@in1 output=@in2 fixpix+ oversca- trim- zerocor- flatcor+ fixfile=badpix.dat flat=flat_bn.fits

# Divide flattened flat into original flat to obtain blaze function
!ln -s flat_b.fits blaze_b.fits
ccdproc images=blaze_b.fits output=blaze_bn.fits fixpix+ oversca- trim- zerocor- flatcor+ fixfile=badpix.dat flat=flat_bn.fits

# Fit scattered light distribution in object and blaze exposure
# Set functio=chebyshev order=2
epar apscat1
# Set functio=chebyshev order=5
epar apscat2
apscatter input=object_bn.fits output=object_bns.fits referen=order_trace.fits interac+ find- recente- resize- edit- trace- fittrac- subtrac+ smooth+ fitscat+ fitsmoo+ nsum=50 buffer=2
apscatter input=blaze_bn.fits output=blaze_bns.fits referen=order_trace.fits interac+ find- recente- resize- edit- trace- fittrac- subtrac+ smooth+ fitscat+ fitsmoo+ nsum=50 buffer=2

!ln -s object_bns.fits object_bns_01.fits
!ln -s object_bns.fits object_bns_02.fits
!ln -s object_bns.fits object_bns_03.fits
!ln -s object_bns.fits object_bns_04.fits
!ln -s object_bns.fits object_bns_05.fits
!ln -s object_bns.fits object_bns_06.fits
!ln -s object_bns.fits object_bns_07.fits
!ln -s object_bns.fits object_bns_08.fits
!ln -s object_bns.fits object_bns_09.fits
!ln -s object_bns.fits object_bns_10.fits
!ln -s object_bns.fits object_bns_11.fits
!ln -s object_bns.fits object_bns_12.fits
!ln -s object_bns.fits object_bns_13.fits
!ln -s object_bns.fits object_bns_14.fits
!ln -s blaze_bns.fits blaze_bns_01.fits
!ln -s blaze_bns.fits blaze_bns_02.fits
!ln -s blaze_bns.fits blaze_bns_03.fits
!ln -s blaze_bns.fits blaze_bns_04.fits
!ln -s blaze_bns.fits blaze_bns_05.fits
!ln -s blaze_bns.fits blaze_bns_06.fits
!ln -s blaze_bns.fits blaze_bns_07.fits
!ln -s blaze_bns.fits blaze_bns_08.fits
!ln -s blaze_bns.fits blaze_bns_09.fits
!ln -s blaze_bns.fits blaze_bns_10.fits
!ln -s blaze_bns.fits blaze_bns_11.fits
!ln -s blaze_bns.fits blaze_bns_12.fits
!ln -s blaze_bns.fits blaze_bns_13.fits
!ln -s blaze_bns.fits blaze_bns_14.fits

# Extract object spectra from each subaperture
!/bin/ls object_bns_*.fits > in1
!sed 's/_bns/_bnse/' in1 > in2
!/bin/ls flat_b_*.fits > in3
apall input=@in1 output=@in2 format=echelle referen=@in3 interac- find- recente- resize- edit- trace- fittrac- extract+ extras+ review- nsum=5 backgro=none weights=variance clean+ saturat=30000.0 readnoi=rdnoise gain=gain lsigma=4.0 usigma=4.0 nsubaps=1

# Extract blaze spectra from each subaperture
!/bin/ls blaze_bns_*.fits > in1
!sed 's/_bns/_bnse/' in1 > in2
!/bin/ls flat_b_*.fits > in3
apall input=@in1 output=@in2 format=echelle referen=@in3 interac- find- recente- resize- edit- trace- fittrac- extract+ extras- review- nsum=5 backgro=none weights=variance clean+ saturat=30000.0 readnoi=rdnoise gain=gain lsigma=4.0 usigma=4.0 nsubaps=1

# Normalize the blaze spectra by a fitting a constant value
!/bin/ls blaze_bnse_*.fits > in1
continuum input=@in1 output=@in1 type=ratio replace- wavesca- logscal- overrid- listonl- interac- naverag=1 function=chebyshev order=1 low_rej=3.0 high_re=3.0 niterat=5 grow=0.0 markrej-

# Split the flux and error arrays of object spectra
!/bin/ls object_bnse_*.fits > in1
!sed 's/object_bnse/flx/' in1 > in2
!sed 's/object_bnse/sig/' in1 > in3
scopy input=@in1 output=@in2 bands=1 renumbe- offset=0 merge- rebin-
scopy input=@in1 output=@in3 bands=3 renumbe- offset=0 merge- rebin-

# Blaze-correct flux and sigma spectra
!/bin/ls flx_*.fits > in1
!/bin/ls sig_*.fits > in2
!/bin/ls blaze_bnse_*.fits > in3
sarith input1=@in1 op=/ input2=@in3 output=@in1 rebin- clobber+
sarith input1=@in2 op=/ input2=@in3 output=@in2 rebin- clobber+

noao.imred.echelle
# Set nsum=9 extras- t_step=30 t_funct=chebyshev t_order=4 t_niter=5
#  weights=variance thresho=40.0 coordli=linelists$thar.dat match=0.02
#  fwidth=5.0 cradius=9.0 i_funct=chebyshev i_xorde=4 i_yorde=4
#  i_niter=5 select=match sort=none group=none refit+ lineari- log-
#  flux+
epar sparams

!ln -s thar_bn.fits thar_bn_01.fits
!ln -s thar_bn.fits thar_bn_02.fits
!ln -s thar_bn.fits thar_bn_03.fits
!ln -s thar_bn.fits thar_bn_04.fits
!ln -s thar_bn.fits thar_bn_05.fits
!ln -s thar_bn.fits thar_bn_06.fits
!ln -s thar_bn.fits thar_bn_07.fits
!ln -s thar_bn.fits thar_bn_08.fits
!ln -s thar_bn.fits thar_bn_09.fits
!ln -s thar_bn.fits thar_bn_10.fits
!ln -s thar_bn.fits thar_bn_11.fits
!ln -s thar_bn.fits thar_bn_12.fits
!ln -s thar_bn.fits thar_bn_13.fits
!ln -s thar_bn.fits thar_bn_14.fits

# Extract and identify ThAr spectra of subaperture with best coverage
apall input=thar_bn_09.fits output=wav_09.fits format=echelle referen=flat_b_09.fits interac- find- recente- resize- edit- trace- fittrac- extract+ extras+ review- nsum=5 backgro=none weights=variance clean+ saturat=30000.0 readnoi=rdnoise gain=gain lsigma=3.0 usigma=3.0 nsubaps=1
ecidentify images=wav_09.fits databas=database coordli=linelists$thar.dat match=0.02 maxfeat=100 zwidth=10.0 ftype=emission fwidth=5.0 cradius=9.0 thresho=80.0 minsep=4.0 functio=chebyshev xorder=3 yorder=3 niterat=5 lowreje=3.0 highrej=3.0 autowri-

# Extract and identify ThAr spectra for other apertures using previous
#  solution as first guess.
!sed 's/wav_09/wav_01/' database/ecwav_09 > database/ecwav_01
!sed 's/wav_09/wav_02/' database/ecwav_09 > database/ecwav_02
!sed 's/wav_09/wav_03/' database/ecwav_09 > database/ecwav_03
!sed 's/wav_09/wav_04/' database/ecwav_09 > database/ecwav_04
!sed 's/wav_09/wav_05/' database/ecwav_09 > database/ecwav_05
!sed 's/wav_09/wav_06/' database/ecwav_09 > database/ecwav_06
!sed 's/wav_09/wav_07/' database/ecwav_09 > database/ecwav_07
!sed 's/wav_09/wav_08/' database/ecwav_09 > database/ecwav_08
!sed 's/wav_09/wav_10/' database/ecwav_09 > database/ecwav_10
!sed 's/wav_09/wav_11/' database/ecwav_09 > database/ecwav_11
!sed 's/wav_09/wav_12/' database/ecwav_09 > database/ecwav_12
!sed 's/wav_09/wav_13/' database/ecwav_09 > database/ecwav_13
!sed 's/wav_09/wav_14/' database/ecwav_09 > database/ecwav_14
!/bin/ls thar_bn_*.fits | grep -v "09" > in1
!sed 's/thar_bn/wav/' in1 > in2
!/bin/ls flat_b_*.fits | grep -v "09" > in3
apall input=@in1 output=@in2 format=echelle referen=@in3 interac- find- recente- resize- edit- trace- fittrac- extract+ extras+ review- nsum=5 backgro=none weights=variance clean+ saturat=30000.0 readnoi=rdnoise gain=gain lsigma=3.0 usigma=3.0 nsubaps=1
ecidentify images=@in2 databas=database coordli=linelists$thar.dat match=0.02 maxfeat=100 zwidth=10.0 ftype=emission fwidth=5.0 cradius=9.0 thresho=80.0 minsep=4.0 functio=chebyshev xorder=5 yorder=5 niterat=5 lowreje=3.0 highrej=3.0 autowri-

# Attach wavelength solutions to corresponding flux and sigma spectra
!/bin/ls flx_*.fits > in1
!/bin/ls sig_*.fits > in2
!/bin/ls wav_*.fits > in3
refspectra input=@in1 referen=@in3 refaps=' ' ignorea+ select=match sort=none group=none time=no overrid+ confirm- assign+
dispcor input=@in1 output=@in1 lineari=no databas=database log=no flux=yes blank=0.0 samedis- global- ignorea- confirm- listonl- verbose-
refspectra input=@in2 referen=@in3 refaps=' ' ignorea+ select=match sort=none group=none time=no overrid+ confirm- assign+
dispcor input=@in2 output=@in2 lineari=no databas=database log=no flux=yes blank=0.0 samedis- global- ignorea- confirm- listonl- verbose-

# UVES_popler hd140283_red.list -filetype 1 -disp 0.7 -ordsig 13.0 -ordsignbr 8.0 -ordsigzero 11.0 -clipsig 9.0 -scalclip 9.0 -pctlred 0.05 -rsigredl 7.0 -rsigredu 9.0 -vcred 3500.0