###############################
#
# ORION-S HC3N Reduction Script
# tool version
# Position-Switched data
#
###############################
# clean up files from the previous run
import os
os.system('rm -rf orions_pscal2123.asap orions_hc3n_cal_tpave.asap SCAN0_CYCLE0_BEAM0_IF0.txt orions_hc3n_reduced.eps orions_hc3n_fit.txt')

# if you haven't loaded ASAP module yet..
asap_init()

s=sd.scantable('OrionS_rawACSmod',False)#load the data without averaging

# dislay summary of the data set (currently all ASAP messages go to the terminal)
s.summary()
#flux unit is missing, so set it to the GBT default unit, Kelvin
s.set_fluxunit('K')

# calibrate HC3N position switched data
scal=sd.calps(s,[20,21,22,23])
# save the intermidiate data (to be used in plotting_demo.py)
scal.save('orions_pscal2123.asap')
# clean up
del s
### opacity correction
scal.opacity(0.09)

### data selection
sel=sd.selector()                       # Prepare a selection
sel.set_ifs(0)                          # select HC3N IF
scal.set_selection(sel)                 # get this IF

### averaging (time and polarization)
stave=sd.average_time(scal,weight='tintsys')    # average in time
spave=stave.average_pol(weight='tsys')  # average polarizations;Tsys-weighted (1/Tsys**2) average
sd.plotter.plot(spave)                  # plot

### boxcar smoothing with width 5
spave.smooth('boxcar',5)
spave.save('orions_hc3n_cal_tpave.asap')

### baselining
spave.auto_poly_baseline(order=2, edge=[1000])
sd.plotter.plot(spave)

spave.set_unit('GHz')
sd.plotter.plot(spave)
sd.plotter.set_histogram(hist=True)
sd.plotter.axhline(color='r',linewidth=2)
sd.plotter.save('orions_hc3n_reduced.eps')

###statistics
spave.set_unit('channel')
# rms on line-free region
rmsmask=spave.create_mask([5000,7000])
rms=spave.stats(stat='rms',mask=rmsmask)
# line statitics
linemask=spave.create_mask([3900,4200])
max=spave.stats('max',linemask)         #  IF[0] = 0.918
sum=spave.stats('sum',linemask)         #  IF[0] = 64.994
median=spave.stats('median',linemask)   #  IF[0] = 0.091
mean=spave.stats('mean',linemask)       #  IF[0] = 0.210

###fitting
spave.set_unit('channel')
sd.plotter.plot(spave)
f=sd.fitter()
msk=spave.create_mask([3928,4255])
# fit a Gaussian
f.set_function(gauss=1)
f.set_scan(spave,msk)
f.fit()
# plot fit result
f.plot(residual=True)
fitparam=f.get_parameters()
#   0: peak = 0.807 K , centre = 4091.256 channel, FWHM = 70.811 channel
#      area = 60.860 K channel
f.store_fit('orions_hc3n_fit.txt')

# try fitting with 2 Gaussians
f.set_function(gauss=2)
f.fit()
f.plot(components=[0,1,-1],residual=True,plotparms=True)

### Save the spectrum
# to ASCII file
spave.save('orions_hc3n_reduced','ASCII',True)
# to MS2  (SDFITS does not work for pol averaged data)
spave.save('orions_hc3n_reduced.ms','MS2',True)

###############################
#
# ORION-S HC3N Reduction Script
# task version
# Position-Switched data
#
###############################
# if you haven't loaded ASAP module yet..
#asap_init()

#clean up for previous run
import os
os.system('rm -rf sdtaskdemo_ori_hc3n.fit')
#sdlist() -summary info
default('sdlist')
infile='OrionS_rawACSmod'
go

#sdcal() -import, calibrate, smooth, average, baseline subtraction
inp('sdcal')

#fix fluxunit
telescope='FIX'
fluxunit='K'

#set calibration mode to 'ps'(=position switched)
calmode='ps'

#select scan numbers
scanlist=[20,21,22,23]

#select IF 0 for HC3N line
iflist=[0]

#averaging
timeaverage=True
polaverage=True

#kernel for smoothing
kernel='boxcar'

#opacity correction
tau=0.09

#baseline subtraction
blmode='auto'
blpoly=2
edge=[1000] # drop 1000 channels from ecah edge for line search

plotlevel=2

#output
sdfile='sdtaskdemo_ori_hc3n.asap'
inp
#save the input parameters
saveinputs('sdcal', 'sdcal.orions.save')
sdcal()

#plotting with sdplot()
default(sdplot)
sdfile='sdtaskdemo_ori_hc3n.asap'
inp
sdplot()

specunit='GHz'
sdplot()

specunit='km/s'
sdplot()

### sdstat() - get some statistics
inp(sdstat)
telescope=''
fluxunit=''
specunit=''
scanlist=[]
sdfile='sdtaskdemo_ori_hc3n.asap'
# rms on line-free region
rmsmask=[5000,7000]
inp
sdstat()

### sdfit() - line fitting
inp(sdstat)
sdfile='sdtaskdemo_ori_hc3n.asap'
telescope=''
fluxunit=''
specunit=''
fitmode='list'  # use list of regions
maskline=[3928,4255] # line fitting region
nfit=[1]
fitfile='sdtaskdemo_ori_hc3n.fit'
inp
sdfit()

#plotting demo
# Plot each scan of calibrated position switched observations of OrionS
#######################################################################
print "plotting demo..."
#create a scantable object (from saved calibrated data after running
#ori_hc3n_demo.py)
s=sd.scantable('orions_pscal2123.asap')
# time average each scan
scanav = s.average_time(scanav=True)
# create a selector object and make selections
sel=sd.selector()
sel.set_polarisations(0)
# apply selections to the scantable
scanav.set_selection(sel)
# set stacking = ifs and panelling = time
sd.plotter.set_mode('if','time')
sd.plotter.set_range(1000,7000,-1,5)
print "Plotting spectra in channel, stacking IF,multiple pannels for different scans"
sd.plotter.plot(scanav)

#reset the selection
sel.reset()
scanav.set_selection(sel)
#new selection
sel.set_ifs(0)
scanav.set_selection(sel)
scanav.set_unit('GHz')
sd.plotter.set_range()
sd.plotter.set_mode('p','time')
print "Plotting in frequency, stacking=pol, panelling=time"
sd.plotter.plot(scanav)

# use of plot_lines?
# get the reduced data
s=sd.scantable('orions_hc3n_reduced.ms')
s.set_unit('GHz')
sd.plotter.plot(s)
jpl=sd.linecatalog('asap_jpl.txt')
jpl.set_frequency_limits(44,46,'GHz')
jpl.set_name('HCCCN')
sd.plotter.plot_lines(jpl, doppler=-27, location='top')

# use of linefinder
# start from calibrated, time and polarization averaged data
s=sd.scantable('orions_hc3n_cal_tpave.asap')
sd.plotter.plot(s)
lf=sd.linefinder()
lf.set_scan(s)
lf.set_options(threshold=3)
lf.find_lines(edge=1000)
baselined_s=s.poly_baseline(mask=lf.get_mask(),order=3,insitu=False)
sd.plotter.plot(baselined_s)
sd.plotter.set_title('orions hc3n baselined spectra')
## auto_poly_baseline

#line-free region rms
rms=baselined_s.stats('rms',lf.get_mask())

#line statistics
msk=baselined_s.create_mask(lf.get_ranges())
max=baselined_s.stats('max',mask=msk)
min=baselined_s.stats('min',mask=msk)
mean=baselined_s.stats('mean',mask=msk)
max=baselined_s.stats('max',mask=msk)