##################################################################
# hdsql : Subaru HDS Quick Look Script 
#  Originaly developed by Wako Aoki
#    revised by Akito Tajitsu <tajitsu@subaru.naoj.org>
#              2011.10.28 ver.2.00
###################################################################
procedure hdsql(inid)
### Input parameters
 string inid {prompt = 'Input frame ID'}
 string indirec {prompt = 'directory of Input data\n'}

 bool   overscan=yes {prompt = 'Overscan?'}
 bool   biassub=no  {prompt = 'BIAS / Dark Subtraction?'}
 bool   maskbad=yes {prompt = 'Mask Bad Pixels?'}
 bool   linear=no   {prompt = 'Linearity Correction?'}
 bool   cosmicra=no {prompt = 'Cosmicray-event rejection?'}
 bool   scatter=no {prompt = 'Scattered light subtraction?'}
 bool   flat=no {prompt = 'Flat fielding?'}
 bool   apall=no {prompt = 'Extract spectra with apall?'}
 bool   wavecal=no  {prompt = 'Wavelength calibration?'}
 bool   remask=no {prompt = 'Re-Mask wavlength calibrated spectrum?'}
 bool   rvcorrect=no  {prompt = 'Heliocentric wavelength correction?'}
 bool   splot=no    {prompt = 'Splot Spectrum?\n\n### Overscan ###'}

# Parameters for overscan
 bool   os_save=yes {prompt = 'Save overscaned data?\n\n### Bias/Dark Subtraction ###'}

# Parameters for BIAS Subtraction
 bool   bs_save=yes {prompt = 'Save BIAS / Dark subtracted data?'}
 string bs_in="" {prompt = 'Input frame for BIAS / Dark subtraction (if necessary)'}
 string bs_style="bias" {prompt = 'Subtraction style (bias|dark)', enum="bias|dark|both"}
 string bs_refer {prompt = 'BIAS frame'}
 string bs_dark  {prompt = 'Dark + BIAS frame\n\n### Masking Bad Pixels ###'}

# Parameters for maskbad
 bool   mb_save=yes {prompt = 'Save masked data?'}
 string mb_in="" {prompt = 'Input frame for Masking Bad Pixels (if necessary)'}
 string mb_refer {prompt = 'Bad Pix Mask frame'}
 bool   mb_auto=no {prompt = 'Auto mask creation from BIAS(bsrefer)?'}
 int   mb_upper=300 {prompt = 'Upper limit for mb_auto'}
 int   mb_lower=-100 {prompt = 'Lower limit for mb_auto'}
 bool  mb_clean=yes  {prompt = 'Cleaning by wacosm11?'}
 int   mb_base=1     {prompt = 'Baseline for wacosm11\n\n### Linearity Correction ###'}

# Parameters for Linearity Correction
 bool   ln_save=yes {prompt = 'Save Linearity Corrected data?'}
 string ln_in="" {prompt = 'Input frame for Linearity Correction (if necessary)\n\n### Cosmic-ray Rejection ###'}

# Parameters for cosmicray-event rejection
 bool   cr_save=yes {prompt = 'Save cosmicray processed data?'}
 string cr_in="" {prompt = 'Input frame for wacosm1 (if necessary)'}
 real   cr_base=2000  {prompt = 'Baseline for wacosm1\n\n### Scattered-light Subtraction ###'}
# real   cr_lower  {prompt = 'Lower limit of replacement window'}

# scattered light subtraction
 bool   sc_save=yes {prompt = 'Save scattered light subtracted data?'}
 string sc_in="" {prompt = 'Input frame for scattered light subtraction (if necessary)'}
 string sc_refer {prompt = 'Reference for aperture finding'}
 bool   sc_inter=yes {prompt = 'Run apscatter interactively?'}
 bool   sc_recen=yes {prompt = 'Recenter apertures for apscatter?'}
 bool   sc_resiz=yes {prompt = 'Resize apertures for apscatter?'}
 bool   sc_edit=yes  {prompt = 'Edit apertures for apscatter?'}
 bool   sc_trace=no {prompt = 'Trace apertures for apscatter?'}
 bool   sc_fittr=no {prompt = 'Fit the traced points interactively for apscatter?\n\n### Flat Fielding ###'}

# Flat fielding
 bool   fl_save=yes {prompt = 'Save flat-fielded data?'}
 string fl_in="" {prompt = 'Input frame for flat fielding (if necessary)'}
 string fl_refer {prompt = 'Flat frame\n\n### Aperture Extraction ###'}

# Parameters for apall
 bool   ap_save=yes {prompt = 'Save apalled data?'}
 string   ap_in="" {prompt = 'Input frame for apall (if necessary)?'}
 string ap_refer {prompt = 'Reference frame for apall'}
 bool   ap_inter=yes {prompt = 'Run apall interactively?'}
 bool   ap_recen=yes {prompt = 'Recenter apertures?'}
 bool   ap_resiz=yes {prompt = 'Resize apertures?'}
 bool   ap_edit=yes  {prompt = 'Edit apertures?'}
 bool   ap_trace=no {prompt = 'Trace apertures?'}
 bool   ap_fittr=no {prompt = 'Fit the traced points interactively?'}
 real   ap_llimi=-30 {prompt = 'Lower aperture limit relative to center'}
 real   ap_ulimi=30 {prompt = 'Upper aperture limit relative to center'}
 real   ap_yleve=0.05 {prompt = 'Fraction of peak for automatic width determination?'} 
 bool   ap_peak=yes  {prompt = 'Is ylevel a fraction of the peak?'}
 string   ap_bg="none"  {prompt = 'Background to subtract\n\n### Wavelength Calibration ###', enum="none|average|median|minimum|fit"}

# Parameters for wavelength calibration
 bool   wv_save=yes  {prompt = 'Save wavelength-calibrated data?'}
 string wv_in=""    {prompt = 'Input frame for wavelength calibration (if necessary)'}
 string wv_refer {prompt = 'Reference frame for refspectra'}
 bool wv_log=no   {prompt = 'Logarithmic wavelength scale?\n\n### Re-Mask after Wavelength Calibration###'}

# Parameters for remask
 bool  zm_save=yes  {prompt = 'Save re-masked data?'}
 string zm_in=""    {prompt = 'Input frame for re-mask afre wave calib. (if necessary)'}
 real   zm_val=1.0 {prompt = 'Pixel Value replaced to All Bad Pixels'}
 real   zm_thresh=0.1 {prompt = 'Threshold pixel value for bad column [0-1]\n\n### Heliocentric Wavelength Correction ###'}

# Parameters for Heliocentric wavelength correction
 string rv_in=""    {prompt = 'Input frame for radial velocity correction (if necessary)'}
 string rv_obs="subaru"    {prompt = 'Observatory\n\n### Splot ###'}

 int    sp_line=1 {prompt = 'Splot image line/aperture to plot\n'}

begin
string version="2.00 (10-28-2011)"
string input_id
string input0,input,output
string flag0,flag
string osfile
string mbinfile,mbfile
string bsinfile,bsfile
string irinfile,irfile
string lninfile,lnfile
string crinfile,crfile
string flinfile,flfile
string scinfile,scfile
string apinfile,apfile
string wvinfile,wvfile
string zminfile,zmfile
string rvinfile,rvfile
string apbg
bool wvlog
string rvobs
string mask_ap, mask_wv,temp1, temp2, temp3, maskimg, masktemp, mask0, mskap
int x_in, y_in, x_msk, y_msk, ccd_in, ccd_msk
string dktemp1, dktemp2
int et_in, et_dk
real et_ratio
bool os_done, bs_done, mb_done, ln_done, cr_done,sc_done,fl_done,ap_done
bool wv_done, zm_done,rv_done
string hq_tmp
bool d_ans
string nextin


os_done=no
bs_done=no
mb_done=no
ln_done=no
cr_done=no
sc_done=no
fl_done=no
ap_done=no
wv_done=no
zm_done=no
rv_done=no

wvlog=wv_log
apbg=ap_bg
rvobs=rv_obs

flag=""
nextin=""

input_id=inid

 output  = ("H"+input_id)
 print("output ID" + output)

 input0=(indirec+input_id+".fits[0]")
 input=(indirec+input_id+".fits")
 print("input data= "+input)

# overscan
if (overscan){
printf("\n")
printf("##################################\n")
printf("# [ 1/12] Overscan\n")
printf("##################################\n")

   flag0=flag
   flag="o"
   if (os_save){
#      osfile=os_out
      osfile=(output+flag)
   } else
      osfile=mktemp("tmp_os")
   if((access(osfile))||access(osfile//".fits")){
     printf("*** OverScanned file \"%s\" already exsits!!\n",osfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(osfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
   print("output overscaned data= "+osfile)
   print("# Overscan is now processing...")
   overscan(inimage=input0,outimage=osfile)
   hedit(osfile,'HQ_OS',"done",add+,del-, ver-,show-,update+)
   os_done=yes
   nextin=osfile
} else {
   print("overscan not processing")
   os_done=no
}


# bias
if (biassub){
printf("\n")
printf("##################################\n")
printf("# [ 2/12] Bias / Dark Subtraction\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"b"	
   if(flag0==""){
      if(bs_in==""){
        bsinfile=input
      }
      else{
        bsinfile=bs_in
        output =bs_in
      }
   }
   else{
     bsinfile=nextin
   }
 #
   if (bs_save){
      bsfile=(output+flag)
      print("output Bias/Dark Subtracted data= "+bsfile)
   } else
      bsfile=mktemp("tmp_bs")
   if((access(bsfile))||access(bsfile//".fits")){
     printf("*** Bias/Dark subtracted file \"%s\" already exsits!!\n",bsfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(bsfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
   if (bs_style == "bias"){
     print("# bias subtraction is now processing...")
     imarith(bsinfile,'-', bs_refer, bsfile, ver-,noact-)
   }
   else {
     dktemp1=mktemp("tmp_dk")
     imarith(bs_dark,'-', bs_refer, dktemp1, ver-,noact-)

     imgets(bsinfile,'EXPTIME')
     et_in=int(imgets.value)
     imgets(dktemp1,'EXPTIME')
     et_dk=int(imgets.value)

     if(et_in==et_dk){
       print("# dark subtraction is now processing...")
       imarith(bsinfile,'-', dktemp1, bsfile, ver-,noact-)
     }
     else{
       dktemp2=mktemp("tmp_dk")
       printf("# scaling dark frame x %d/%d...\n",et_in,et_dk)
       et_ratio=real(et_in)/real(et_dk)
       imarith(dktemp1,'*', et_ratio, dktemp2, ver-,noact-)
       print("# dark subtraction is now processing...")
       imarith(bsinfile,'-', dktemp2, bsfile, ver-,noact-)
       
       imdelete(dktemp2)
     }
     imdelete(dktemp1)

     if (bs_style == "both"){
       dktemp1=mktemp("tmp_dk")
       imrename(bsfile,dktemp1)

       print("# bias subtraction is now processing...")
       imarith(dktemp1,'-', bs_refer, bsfile, ver-,noact-)

       imdelete(dktemp1)
     }
   }
   hedit(bsfile,'HQ_BS',"done",add+,del-, ver-,show-,update+)
   bs_done=yes
   nextin=bsfile
} else  {
   print("bias/dark subtraction not processing")     
   bs_done=no
}


# mask bad column
if (maskbad){
printf("\n")
printf("##################################\n")
printf("# [ 3/12] Masking Bad Pixels\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"m"	
   if(flag0==""){
      if(mb_in==""){
        mbinfile=input
      }
      else{
        mbinfile=mb_in
        output =mb_in
      }
   }
   else{
     mbinfile=nextin
   }
#
   if (mb_save){
      mbfile=(output+flag)
      print("output bad pix masked data= "+mbfile)
   } else
      mbfile=mktemp("tmp_mb")
   if((access(mbfile))||access(mbfile//".fits")){
     printf("*** Masked file \"%s\" already exsits!!\n",mbfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(mbfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
   print("# masking bad pixels...")

## Auto Creation
   if(mb_auto){
     masktemp=mktemp("tmp_msk")
     printf("# Creating mask from BIAS \"%s\"...\n",bs_refer)
     mkbadmask(bs_refer, masktemp, lower=mb_lower, upper=mb_upper, clean=mb_clean, base=mb_base)
     maskimg=masktemp
   }
   else{
     maskimg=mb_refer
   }

## Get Header Information
    imgets(mbinfile,'i_naxis1')
    x_in=int(imgets.value)
    imgets(mbinfile,'i_naxis2')
    y_in=int(imgets.value)

    imgets(maskimg,'i_naxis1')
    x_msk=int(imgets.value)
    imgets(maskimg,'i_naxis2')
    y_msk=int(imgets.value)

    if( (x_msk!=x_in) || (y_msk!=y_in) ){
      printf("####### Error !! MASK dimenstion mismatch!! ABORT!! #######\n")
      bye
    }

    imgets(mbinfile,'DET-ID')
    ccd_in=int(imgets.value)
    imgets(maskimg,'DET-ID')
    ccd_msk=int(imgets.value)

    if(ccd_msk!=ccd_in){
      printf("####### Error !! MASK CCD Color mismatch!! ABORT!! #######\n")
      bye
    }

   imcopy(mbinfile, mbfile)
   fixpix(mbfile, maskimg,linterp=INDEF, cinterp=INDEF, ver+, pixels-)
   hedit(mbfile,'H_MASK0',maskimg,add+,del-, ver-,show-,update+)
   hedit(mbfile,'HQ_MB',"done",add+,del-, ver-,show-,update+)
   mb_done=yes
   nextin=mbfile
   mask0=maskimg
} else {
   print("bad column mask not processing")     
   mb_done=no
}

# linearity
if (linear){
printf("\n")
printf("##################################\n")
printf("# [ 4/12] Linearity Correction\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"l"	
   if(flag0==""){
      if(ln_in==""){
        lninfile=input
      }
      else{
        lninfile=ln_in
        output=ln_in
      }

      imgets(lninfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(lninfile,'H_MASK0')
         mask0=imgets.value
      }
   }
   else{
     lninfile=nextin
   }
#
   if (ln_save){
      lnfile=(output+flag)
      print("output linearity corrected data= "+lnfile)
   } else
      lnfile=mktemp("tmp_ln")
   if((access(lnfile))||access(lnfile//".fits")){
     printf("*** Linearity Corrected file \"%s\" already exsits!!\n",lnfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(lnfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
   print("# linearity correction is now processing...")
   hdslinear(lninfile,lnfile,outtype="real", auto_b+)
   hedit(lnfile,'HQ_LN',"done",add+,del-, ver-,show-,update+)
   ln_done=yes
   nextin=lnfile
} else {
   print("linearity correction not processing")     
   ln_done=no
}



# wacosm11
if (cosmicra){
printf("\n")
printf("##################################\n")
printf("# [ 5/12] Cosmic Ray Removal\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"c"	
   if(flag0==""){
      if(cr_in==""){
        crinfile=input
      }
      else{
        crinfile=cr_in
        output=cr_in
      }

      imgets(crinfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(crinfile,'H_MASK0')
         mask0=imgets.value
      }
   }
   else{
     crinfile=nextin
   }
#
   if (cr_save){
      crfile=(output+flag)
      print("output cosmic-ray removed data= "+crfile)
   } else
      crfile=mktemp("tmp_cr")
   if((access(crfile))||access(crfile//".fits")){
     printf("*** Cosmic Ray Rejected file \"%s\" already exsits!!\n",crfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(crfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
   print("# wacosm11 is now processing...")
   wacosm11 (in_f=crinfile,out_f=crfile,base=cr_base)
   hedit(crfile,'HQ_CR',"done",add+,del-, ver-,show-,update+)
   cr_done=yes
   nextin=crfile
} else {
   print("wacosm11 not processing")     
   cr_done=no
}

# scattered light subtraction
if (scatter){
printf("\n")
printf("##################################\n")
printf("# [ 6/12] Scattered Light Subtraction\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"s"	
   if(flag0==""){
      if(sc_in==""){
        scinfile=input
      }
      else{
        scinfile=sc_in
        output=sc_in
      }

      imgets(scinfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(scinfile,'H_MASK0')
         mask0=imgets.value
      }
   }
   else{
     scinfile=nextin
   }
#
   if (sc_save){
#      scfile=sc_out
      scfile=output+flag
   } else
      scfile=mktemp("tmp_sc")
   if((access(scfile))||access(scfile//".fits")){
     printf("*** Scattered Light Subtracted file \"%s\" already exsits!!\n",scfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(scfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
    imgets(scinfile,'DET-ID')
    ccd_in=int(imgets.value)
    imgets(sc_refer,'DET-ID')
    ccd_msk=int(imgets.value)

    if(ccd_msk!=ccd_in){
      printf("####### Error !! Aperture Reference CCD Color mismatch!! ABORT!! #######\n")
      bye
    }

   print("# Scattered light subtracting is now processing...")
   apscatter(scinfile,scfile,interac=sc_inter,referen=sc_refer,recente=sc_recen,resize=sc_resiz,edit=sc_edit,trace=sc_trace,fittrac=sc_fittr)
   hedit(scfile,'HQ_SC',"done",add+,del-, ver-,show-,update+)
   sc_done=yes
   nextin=scfile
} else {
   print("Scattered light subtraction not processing")     
   sc_done=no
}

# flat fielding
if (flat){
printf("\n")
printf("##################################\n")
printf("# [ 7/12] Flat Fielding\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"f"	
   if(flag0==""){
      if(fl_in==""){
        flinfile=input
      }
      else{
        flinfile=fl_in
        output=fl_in
      }

      imgets(flinfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(flinfile,'H_MASK0')
         mask0=imgets.value
      }
   }
   else{
     flinfile=nextin
   }
#
   if (fl_save){
#      flfile=fl_out
      flfile=(output+flag)
      print("output flatted data= "+flfile)
   } else
      flfile=mktemp("tmp_fl")
   if((access(flfile))||access(flfile//".fits")){
     printf("*** Flat Fielded file \"%s\" already exsits!!\n",flfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(flfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
    imgets(flinfile,'DET-ID')
    ccd_in=int(imgets.value)
    imgets(fl_refer,'DET-ID')
    ccd_msk=int(imgets.value)

    if(ccd_msk!=ccd_in){
      printf("####### Error !! Flat CCD Color mismatch!! ABORT!! #######\n")
      bye
    }

   print("# Flat fielding is now processing...")
   print(flinfile,"/",fl_refer,"=>",flfile)     
   imarith(flinfile,"/",fl_refer,flfile)
   hedit(flfile,'HQ_FL',"done",add+,del-, ver-,show-,update+)
   fl_done=yes
   nextin=flfile
} else {
   print("flat-fielding not processing")     
   fl_done=no
}

# apall
if (apall){
  printf("\n")
  printf("##################################\n")
  printf("# [ 8/12] Aperture Extraction\n")
  printf("##################################\n")

   flag0=flag
   flag=flag+"_ec"	
   if(flag0==""){
      if(ap_in==""){
        apinfile=input
      }
      else{
        apinfile=ap_in
        output=ap_in
      }

      imgets(apinfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(apinfile,'H_MASK0')
         mask0=imgets.value
      }
   }
   else{
     apinfile=nextin
   }
#
   if (ap_save){
#      apfile=ap_out
      apfile=output+flag
   } else
      apfile=mktemp("tmp_ap")
   if((access(apfile))||access(apfile//".fits")){
     printf("*** Aperture Extracted file \"%s\" already exsits!!\n",apfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(apfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#
    imgets(apinfile,'DET-ID')
    ccd_in=int(imgets.value)
    imgets(ap_refer,'DET-ID')
    ccd_msk=int(imgets.value)

    if(ccd_msk!=ccd_in){
      printf("####### Error !! Aperture Reference CCD Color mismatch!! ABORT!! #######\n")
      bye
    }

   print("# apall is now processing...")

  if(apbg=='none'){
     apall (input=apinfile,output=apfile, apertures="", format="echelle",
       references=ap_refer,profiles=" ", interactive=ap_inter, find=no, 
       recenter=ap_recen, resize=ap_resiz, edit=ap_edit,trace=ap_trace,
       fittrace=ap_fittr, extract=yes, extras=yes, review=no, line=INDEF,
       nsum=100, lower=-10., upper=10., apidtable="",
       width=10.,radius=5.,threshold=0., 
       minsep=5., maxsep=1000.,order="increasing",aprecenter="",
       npeaks=INDEF, shift=no,llimit=ap_llimi,ulimit=ap_ulimi, 
       ylevel=ap_yleve, peak=yes, 
       bkg=no,r_grow=0.1, avglimits=yes, t_nsum=10, t_step=3, t_nlost=10,
       t_function="legendre",t_order=3,t_sample="*",t_naverage=1,t_niterate=2,
       t_low_reject=3., t_high_rejec=3., t_grow=0., background=apbg, skybox=1,
       weights="none", pfit="fit1d", clean=no, saturation=INDEF, readnoise="0.",
       gain="1.", lsigma=4., usigma=4., nsubaps=1)
   }else{
     apall (input=apinfile,output=apfile, apertures="", format="echelle",
       references=ap_refer,profiles=" ", interactive=ap_inter, find=no, 
       recenter=ap_recen, resize=ap_resiz, edit=ap_edit,trace=ap_trace,
       fittrace=ap_fittr, extract=yes, extras=no, review=no, line=INDEF,
       nsum=100, lower=-10., upper=10., apidtable="",
       width=10.,radius=5.,threshold=0., 
       minsep=5., maxsep=1000.,order="increasing",aprecenter="",
       npeaks=INDEF, shift=no,llimit=ap_llimi,ulimit=ap_ulimi, 
       ylevel=ap_yleve, peak=yes, 
       bkg=no,r_grow=0.1, avglimits=yes, t_nsum=10, t_step=3, t_nlost=10,
       t_function="legendre",t_order=3,t_sample="*",t_naverage=1,t_niterate=2,
       t_low_reject=3., t_high_rejec=3., t_grow=0., background=apbg, skybox=1,
       weights="none", pfit="fit1d", clean=no, saturation=INDEF, readnoise="0.",
       gain="1.", lsigma=4., usigma=4., nsubaps=1)
  }

   hedit(apfile,'HQ_AP',"done",add+,del-, ver-,show-,update+)
   hedit(apfile,'H_MSKAP',apinfile,add+,del-, ver-,show-,update+)
   mskap=apinfile
   ap_done=yes
   nextin=apfile
} else {
    print("apall not processing")
   ap_done=no
}



# Wavelength calibration (refspectra + dispcor)
  if (wavecal){
printf("\n")
printf("##################################\n")
printf("# [ 9/12] Wavelength Calibration\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"w"	
   if(flag0==""){
      if(wv_in==""){
        wvinfile=input
      }
      else{
        wvinfile=wv_in
        output=wv_in
      }

      imgets(wvinfile,'HQ_AP')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         ap_done=yes
         hq_tmp=""
      }
      else {
         printf("####### Error !! Input file for Wavelength Calibration must be aperture extracyed!! ABORT!! #######\n")
         bye
      }

      imgets(wvinfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(wvinfile,'H_MASK0')
         mask0=imgets.value

         imgets(wvinfile,'H_MSKAP')
         mskap=imgets.value
      }
   }
   else{
      if(!ap_done){
         printf("####### Error !! Input file for Wavelength Calibration must be aperture extracted!! ABORT!! #######\n")
         bye
      }
      wvinfile=nextin
   }

   if (wv_save){
      wvfile=output+flag
   } else
      wvfile=mktemp("tmp_wv")
   if((access(wvfile))||access(wvfile//".fits")){
     printf("*** Wavelength Calibrated file \"%s\" already exsits!!\n",wvfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(wvfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }
#

   imgets(wvinfile,'DET-ID')
   ccd_in=int(imgets.value)
   imgets(wv_refer,'DET-ID')
   ccd_msk=int(imgets.value)

   if(ccd_msk!=ccd_in){
     printf("####### Error !! Aperture Reference CCD Color mismatch!! ABORT!! #######\n")
     bye
   }

   print("# Wavelength calibration is now processing...")
   refspectra(input=wvinfile,answer=yes,referen=wv_refer,sort="MJD", group=" ")
   if(wvlog){
      dispcor(input=wvinfile,output=wvfile, log+)
   }
   else{
      dispcor(input=wvinfile,output=wvfile, log-)
    }
    hedit(wvfile,'HQ_WV',"done",add+,del-, ver-,show-,update+)
    wv_done=yes
    nextin=wvfile
} else {
     print("wavelength-calibration not processing")
     wv_done=no
}

if(remask){
printf("\n")
printf("##################################\n")
printf("# [10/12] Re-Masking Wavelength Calibrated Spectrum\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"z"	
   if(flag0==""){
      if(zm_in==""){
        zminfile=input
      }
      else{
        zminfile=zm_in
        output=zm_in
      }

      imgets(zminfile,'HQ_AP')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         ap_done=yes
         hq_tmp=""
      }
      else {
         printf("####### Error !! Input file for Re-Masking must be aperture extracted!! ABORT!! #######\n")
         bye
      }

      imgets(zminfile,'HQ_WV')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         wv_done=yes
         hq_tmp=""
      }
      else {
         printf("####### Warning! Re-Masking only for Aperture Extracted data\n")
      }


      imgets(zminfile,'HQ_MB')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         mb_done=yes
         hq_tmp=""
          
         imgets(zminfile,'H_MASK0')
         mask0=imgets.value

         imgets(zminfile,'H_MSKAP')
         mskap=imgets.value
      }
      else{
         printf("####### Error !! Input file has not been masked!! ABORT!! #######\n")
         bye
      }
   }
   else{
      if(!mb_done){
         printf("####### Error !! Input file for Re-masking has not been masked!! ABORT!! #######\n")
         bye
      }
      if(!ap_done){
         printf("####### Error !! Input file for Re-masking must be aperture extracted!! ABORT!! #######\n")
         bye
      }
      if(!wv_done){
         printf("####### Warning! Re-Masking only for Aperture Extracted data\n")
      }
      zminfile=nextin
    }
  
    printf("# Extracting mask referring \"%s\"...\n", mskap)
    mask_ap=mktemp("tmp.mask")
    
    if(access(mask_ap//".fits")) imdelete(mask_ap)

    apall (input=mask0 ,output=mask_ap, apertures=" ", format="echelle",
       references=mskap,profiles=" ", interactive=no, find=no, 
       recenter=no, resize=no, edit=no,trace=no,
       fittrace=no, extract=yes, extras=no, review=no)

    imreplace(mask_ap, 1, imagina=0.,upper=INDEF, lower=1.0, radius=0.)

    printf(">>> Mask for fixpix \"%s\" has been created!\n",mask_ap)

    if (zm_save){
       zmfile=output+flag
     } else 
       zmfile=mktemp("tmp_zm")
    
   if((access(zmfile))||access(zmfile//".fits")){
     printf("*** Re-masked file \"%s\" already exsits!!\n",zmfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(zmfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }

    temp1=mktemp("tmp_mask")
    temp2=mktemp("tmp_mask")

   if(!wv_done){
        imreplace(mask_ap, 0, imagina=0.,upper=zm_thresh, lower=INDEF, radius=0.)
        imreplace(mask_ap, 1, imagina=0.,upper=INDEF, lower=zm_thresh, radius=0.)

        imarith(zminfile,"*", mask_ap, temp1)

        temp3=mktemp("tmp_mask")
        imarith(mask_ap,"*", zm_val, temp3)
        imarith(temp1,"-", temp3, temp2)
        imdelete(temp3)

        imarith(zminfile,"-", temp2, zmfile)

        mask_wv="Mask_"//nextin
        if(access(mask_wv//".fits")) imdelete(mask_wv)
        imrename (mask_ap, mask_wv)
   }else{
      print("# Wavelength calibration for MASK is now processing...")
      if(access(mask_ap//".fits")){
        mask_wv="Mask_"//wvfile
        if(access(mask_wv//".fits")) imdelete(mask_wv)
 
        refspectra(input=mask_ap,answer=yes,referen=wv_refer,sort="MJD", group=" ")
        if(wvlog){
           dispcor(input=mask_ap,output=mask_wv, log+)
        }
        else{
          dispcor(input=mask_ap,output=mask_wv, log-)
        }
        imreplace(mask_wv, 0, imagina=0.,upper=zm_thresh, lower=INDEF, radius=0.)
        imreplace(mask_wv, 1, imagina=0.,upper=INDEF, lower=zm_thresh, radius=0.)
      }

      imarith(zminfile,"*", mask_wv, temp1)

      temp3=mktemp("tmp_mask")
      imarith(mask_wv,"*", zm_val, temp3)
      imarith(temp1,"-", temp3, temp2)
      imdelete(temp3)

      imarith(zminfile,"-", temp2, zmfile)
      imdelete(mask_ap)

    }
    imdelete(temp2)
    imdelete(temp1)

    hedit(zmfile,'H_MASK',mask_wv,add+,del-, ver-,show-,update+)
    hedit(zmfile,'HQ_ZM',"done",add+,del-, ver-,show-,update+)
    zm_done=yes
    nextin=zmfile
}
else{
     print("re-masking not processing")
     zm_done=no
}


# Heliocentric wavelength  (rvcorrect + rvhds)
if (rvcorrect){
printf("\n")
printf("##################################\n")
printf("# [11/12] Heliocentric Wavelength Correction\n")
printf("##################################\n")

   flag0=flag
   flag=flag+"r"	
   if(flag0==""){
      if(rv_in==""){
        rvinfile=input
      }
      else{
        rvinfile=rv_in
        output=rv_in
      }

      imgets(rvinfile,'HQ_WV')
      hq_tmp=imgets.value
      if(hq_tmp=="done"){
         wv_done=yes
         hq_tmp=""
      }
      else {
         printf("###### ERROR !! Input for RV correction must be Wavelength Calibrated !!\n")
	 bye
      }
   }
   else{
     if(!wv_done){
         printf("###### ERROR !! Input for RV correction must be Wavelength Calibrated !!\n")
	 bye
     }
     rvinfile=nextin
   }

   rvfile=output+flag

   if((access(rvfile))||access(rvfile//".fits")){
     printf("*** Heliocentric RV Corrected file \"%s\" already exsits!!\n",rvfile)
     printf(">>> Do you want to overwrite this file? (y/n) : ")
     while(scan(d_ans)!=1) {}
     if(d_ans){
        imdelete(rvfile)
     }
     else{
        printf("!!! Please remove exsiting file !!! ABORT !!!\n")
        bye
     }
   }

   rvhds(inimage=rvinfile,outimage=rvfile, observa=rvobs)
   hedit(rvfile,'HQ_RV',"done",add+,del-, ver-,show-,update+)
   rv_done=yes
   nextin=rvfile
   printf(">>> Heliocentric wavelength corrected \"%s\" has been created.\n", rvfile)
}
else{
     print("heliocentric wavelength correction not processing")
     rv_done=no
}
 


#
if (overscan){
  if (!os_save){
    imdel(images=osfile)
  }
}
if (biassub){
  if (!bs_save){
    imdel(images=bsfile)
  }
}
if (maskbad){
  if (!mb_save){
    imdel(images=mbfile)
  }
  if (mb_auto){
    imdel(images=masktemp)
  }
}
if (cosmicra){
  if (!cr_save){
    imdel(images=crfile)
  }
}
if (flat){
  if (!fl_save){
    imdel(images=flfile)
  }
}
if (scatter){
  if (!sc_save){
    imdel(images=scfile)
  }
}
if (apall){
  if (!ap_save){
    imdel(images=apfile)
  }
}
if (remask){
  if (!zm_save){
    imdel(images=zmfile)
  }
}
if (wavecal){
  if (!wv_save){
    imdel(images=wvfile)
  }
}

if (splot && ap_done){
printf("\n")
printf("##################################\n")
printf("# [12/12] Plotting Resultant Spectrum\n")
printf("##################################\n")

      splot (images=nextin,line=sp_line,band=1)
}


printf("\n")
printf("##############################################################\n")
printf("# hdsql : FINISH\n")
printf("#   ver %s developped by W.Aoki and A.Tajitsu\n",version)
printf("#\n")
printf("#  Resultant File :   %s%s.fits\n",output,flag)
if ((zm_done)&&(mb_done)){
  printf("#       Mask File :   %s.fits\n",mask_wv)
}
printf("##############################################################\n")

#endofp:


bye
end