### overscan region subtraction
ls *.fits > list
!sed -e 's/FCSA00//' -e 's/fits/bs.fits/' list > list.bs
ovsub_hpk @list @list.bs 
# for spectroscopic data of FOCAS, DO NOT FORGET to reset WCS information
wcsreset @list.bs world


### create and subtract standard-bias
hsel *.bs.fits $I 'OBJECT=="BIAS"&&@"DET-ID"==1' > list.bias1
hsel *.bs.fits $I 'OBJECT=="BIAS"&&@"DET-ID"==2' > list.bias2
imcomb @list.bias1 bias1
imcomb @list.bias2 bias2
hsel *.bs.fits $I '@"DATA-TYP"=="OBJECT"&&@"DET-ID"==1' > list.obj1
hsel *.bs.fits $I '@"DATA-TYP"=="OBJECT"&&@"DET-ID"==2' > list.obj2
!sed 's/bs/bs2/' list.obj1 > list.obj1.bs2
!sed 's/bs/bs2/' list.obj2 > list.obj2.bs2
imarith @list.obj1 - bias1 @list.obj1.bs2
imarith @list.obj2 - bias2 @list.obj2.bs2


### flat fielding
hsel @list.bs $I 'OBJECT=="DOMEFLAT"&&EXPTIME==2.5&&@"DET-ID"==1' > list.flat1
hsel @list.bs $I 'OBJECT=="DOMEFLAT"&&EXPTIME==2.5&&@"DET-ID"==2' > list.flat2
imcomb @list.flat1 flat1
imcomb @list.flat2 flat2
imarith flat1 / 20000 flat1
imarith flat2 / 20000 flat2

hsel *.bs2.fits $I 'OBJECT=="Mrk 6"&&@"DET-ID"==1' > list.mos1
hsel *.bs2.fits $I 'OBJECT=="Mrk 6"&&@"DET-ID"==2' > list.mos2
!sed 's/bs2/ff/' list.mos1 > list.mos1.ff
!sed 's/bs2/ff/' list.mos2 > list.mos2.ff
imarith @list.mos1 / flat1 @list.mos1.ff
imarith @list.mos2 / flat2 @list.mos2.ff


### distortion correction and combine chip1/2
bigimage 129973.ff 129974.ff mos1 obsmode=spec
bigimage 129975.ff 129976.ff mos2 obsmode=spec
bigimage 129977.ff 129978.ff mos3 obsmode=spec
bigimage 129981.ff 129982.ff mos4 obsmode=spec
bigimage 129983.ff 129984.ff mos5 obsmode=spec
bigimage 129985.ff 129986.ff mos6 obsmode=spec
bigimage 129987.ff 129988.ff mos7 obsmode=spec

bigimage 129991.bs 129992.bs compblue obsmode=spec
bigimage 129993.bs 129994.bs compred obsmode=spec


### registration ###
## example of using xregister task
# for measuring shifts, use data before flat-fielding
# 129974.bs2.fits, 129976.bs2.fits ... 129988.bs2.fits
hsel @list.obj2.bs2 $I 'EXPTIME==900' > list.mrk6.chip2
# the reference image is mos7.fits (129988.fits) which is closest to the comparison frames
# output fits files are used just to check the calculation by xregister
# usually lorejec and hirejec are necessary to reject cosmicrays for calculation
# a region with an isolated slit and bright sky line(s) is used for the cross-correlation
xregister @list.mrk6.chip2 129988.bs2 [155:265,1400:1550] xreg.log output=@list.mrk6.chip2.sh lorejec=3 hirejec=3 > xreg.dat
# apply the calculated shifts to mos*.fits (results of bigimage)
!grep pixels xreg.dat | awk '{print "imshift mos"NR".fits mos"NR".sh.fits",$8+0.07255363,$9-0.1795726}' > shift.cl
cl < shift.cl


### extraction
# because the task produces many files, a new directory is created for working place
mkdir extract
cd extract
cp ../compblue.fits ./
cp ../compred.fits ./
cp ../mos?.sh.fits ./

cp ../mrk6obs.sbr ./
cp ../../mask/mrk6obs.mdp ./

# a small correction dx=-1 is applied
mosextract compred.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract compblue.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp

mosextract mos1.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract mos2.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract mos3.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract mos4.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract mos5.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract mos6.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp
mosextract mos7.sh.fits mrk6obs.sbr dx=-1 mdpfile=mrk6obs.mdp


### create comparison frames
# e.g., compblue.001.fits and compred.001.fits -> comp.001.fits
!ls compred.*.fits |sed 's/.fits$//' > list.red
!ls compblue.*.fits |sed 's/.fits$//' > list.blue
!sed 's/red//' list.red > list.comp
mkcomp @list.blue @list.red @list.comp wavelen=6900


### identify
mosident @list.comp linedata=thar.300.dat order=5
# redo unsuccessful one
mosident comp.012 linedata=thar.300.dat order=5

# fitcoords for all spectra
# usually, xorder=2 is enough for 10"-long slit
# yorder should be same as for mosident
fitcoords @list.comp xorder=2 yorder=5


### transform (wavelength calibration)
ls mos1.sh.0*.fits > list.mos1
!sed 's/sh/wc/' list.mos1 > list.mos1.wc
# task mostrans is used to set the dispersion function of all spectra identical
mostrans @list.mos1 @list.mos1.wc @list.comp y1=3700 y2=8300 dy=1.37

ls mos2.sh.0*.fits > list.mos2
!sed 's/sh/wc/' list.mos2 > list.mos2.wc
mostrans @list.mos2 @list.mos2.wc @list.comp y1=3700 y2=8300 dy=1.37

# and so on