import sys
import os
import casac
import string
import inspect
import pdb

from parameter_check import *

casapath=os.environ['AIPSPATH']
import asap as sd
os.environ['AIPSPATH']=casapath

def sdcal(infile=None,telescope=None,fluxunit=None,specunit=None,frame=None,doppler=None,calmode=None,scanlist=None,field=None,iflist=None,scanaverage=None,timeaverage=None,polaverage=None,kernel=None,kwidth=None,tau=None,blmode=None,blpoly=None,interactive=None,masklist=None,sdfile=None,outform=None,plotlevel=None,thresh=None,avg_limit=None,edge=None):
	"""ASAP SD calibration task (prototype):

	Task Version 2007-03-06 STM
	Keyword arguments:
        infile -- name of input SD dataset
	        options: (str) file name
		default: '' (none set) REQUIRED
		example: 'mopra-2005-05-08_0350.rpf'
		         Supported formats: ASAP,MS,RPFITS,SDFITS
	telescope -- the telescope name or characteristics
	        options: (str) name or (list) list of gain info
		default: '' (none set)
		example: telescope='GBT' for the GBT
		         telescope='AT' for one of the default scopes
			 known to ASAP.
			 telescope=[104.9,0.43] diameter(m), ap.eff.
			 telescope=[0.743] gain in Jy/K
			 telescope='FIX' to change default fluxunit
			 see description below
	fluxunit -- units for line flux
	        options: 'K','Jy',''
		default: '' (keep current fluxunit)
		WARNING: For GBT data, see description below.
	specunit -- units for spectral axis
	        options: (str) 'channel','km/s','GHz','MHz','kHz','Hz'
	        default: 'channel'
		example: this will be the units for masklist
	frame -- frequency frame for spectral axis
	        options: (str) 'LSRK','REST','TOPO','LSRD','BARY',
	                 'GEO','GALACTO','LGROUP','CMB'
		default: currently set frame in scantable
		WARNING: frame='REST' not yet implemented
	doppler -- doppler mode
	        options: (str) 'RADIO','OPTICAL','Z','BETA','GAMMA'
	        default: currently set doppler in scantable
	calmode -- calibration mode
	        options: 'ps','nod','fs','fsotf','quotient','none'
	        default: 'none'
		example: choose mode 'none' if you have
		         already calibrated and want to
	                 try baselines or averaging
	scanlist -- list of scan numbers to process
	        default: [] (use all scans)
		example: [21,22,23,24]
		         this selection is in addition to field
			 and iflist
	field -- selection string for selecting scans by name
	        default: '' (no name selection)
		example: 'FLS3a*'
		         this selection is in addition to scanlist
			 and iflist
	iflist -- list of IF id numbers to select
	        default: [] (use all IFs)
		example: [15]
		         this selection is in addition to scanlist
			 and field
	scanaverage -- average integrations within scans
	        options: (bool) True,False
		default: False
		example: if True, this happens in read-in
	                 For GBT, set False!
	timeaverage -- average times for multiple scan cycles
	        options: (bool) True,False
		default: False
		example: if True, this happens after calibration
	polaverage -- average polarizations
	        options: (bool) True,False
		default: False
	kernel -- type of spectral smoothing
	        options: 'hanning','gaussian','boxcar','none'
	        default: 'none'
	kwidth -- width of spectral smoothing kernel
	        options: (int) in channels 
	        default: 5
		example: 5 or 10 seem to be popular for boxcar
		         ignored for hanning (fixed at 5 chans)
                         (0 will turn off gaussian or boxcar)
	tau -- atmospheric optical depth
	        default: 0.0 (no correction)
	blmode -- mode for baseline fitting
	        options: (str) 'auto','list','none'
		default: 'none'
		example: blmode='none' turns off baseline fitting
		         blmode='auto' uses AUTOPARS (see below)
		         in addition to blpoly to run linefinder
			 to determine line-free regions
			 USE WITH CARE! May need to tweak AUTOPARS.
	blpoly -- order of baseline polynomial
	        options: (int) (<0 turns off baseline fitting)
	        default: 5
		example: typically in range 2-9 (higher values
		         seem to be needed for GBT)
	interactive -- interactive mode for baseline fitting
	        options: (bool) True,False
		default: False
		WARNING: Currently this just asks whether you accept
		         the displayed fit and if not, continues
	                 without doing any baseline fit.
	masklist -- list of mask regions to INCLUDE in BASELINE fit
	        default: [] (entire spectrum)
		example: [[1000,3000],[5000,7000]]
		         if blmode='auto' then this mask will be applied
			 before fitting
	sdfile -- Name of output file
	        default: '' (<infile>_cal)
		example: note that sdfile is the OUTPUT of sdcal and
	                 is the INPUT filename param for the other
			 sdtasks to steamline processing
	        WARNING: output file will be overwritten
	outform -- format of output file
	        options: 'ASCII','SDFITS','MS','ASAP'
		default: 'ASAP'
                example: the ASAP format is easiest for further sd
	                 processing; use MS for CASA imaging.
			 If ASCII, then will append some stuff to
			 the sdfile name
	plotlevel -- control for plotting of results
	        options: (int) 0=none, 1=some, 2=more, <0=hardcopy
		default: 0 (no plotting)
		example: plotlevel<0 as abs(plotlevel), e.g.
		         -1 => hardcopy of final plot (will be named
		         <sdfile>_calspec.eps)
		WARNING: be careful plotting in fsotf mode!

	AUTOPARS: the following parameters are used for blmode='auto' ONLY
	-------------------------------------------------------------------
	thresh -- S/N threshold for linefinder
	        default: 5
		example: a single channel S/N ratio above which the channel is
	                 considered to be a detection
	avg_limit -- channel averaging for broad lines
	        default: 4
		example: a number of consequtive channels not greater than
                         this parameter can be averaged to search for broad lines
	edge -- channels to drop at beginning and end of spectrum
	        default: 0
		example: [1000] drops 1000 channels at beginning AND end
		         [1000,500] drops 1000 from beginning and 500 from end

                Note: For bad baselines threshold should be increased,
                and avg_limit decreased (or even switched off completely by
                setting this parameter to 1) to avoid detecting baseline
                undulations instead of real lines.


        DESCRIPTION:

	Task sdcal performs data selection, calibration, and/or spectral
	baseline fitting for single-dish spectra.  By setting calmode='none'
	one can run sdcal on already calibrated data, for further selection
	or baseline fitting.  Likewise, one can set blmode='none' to bypass
	baseline fitting.

	If you give multiple IFs in iflist, then your scantable will have
	multiple IFs.  This can be handled, but there can be funny interactions
	later on.  We recommend you split each IF out into separate files
	by re-running sdcal with each IF in turn.

	ASAP recognizes the data of the "AT" telescopes, but currently
	does not know about the GBT or any other telescope. This task
	does know about GBT.  Therefore, if you wish to change the
	fluxunit (see below), then you need to tell it what to do.  If
	you set telescope='AT' it will use its internal defaults.  If
	you set telescope 'GBT', it will use an approximate aperture
	efficiency conversion.  If you give it a list instead of a
	string, then if the list has a single float it is assumed to
	be the gain in Jy/K, if two or more elements they are assumed
	to be telescope diameter (m) and aperture efficiency
	respectively.

	Note that sdcal assumes that the fluxunit is set correctly in
	the data already.  If not, then set telescope='FIX' and it
	will set the default units to fluxunit without conversion.
	WARNING: If the data in infile is an ms from GBT, it will
	currently say its in 'Jy' but its really 'K', so set
	telescope='FIX' and fluxunit='K' to fix this.
	
	"""
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['__last_taskname']='sdcal'
        tb=myf['tb']
        im = myf['im']
	qa = myf['qa']
        me = myf['me']
	sm = myf['sm']
        ia = myf['ia']
        cl = myf['cl']
        pl = myf['pl']
	rg = myf['rg']
	default=myf['default']
	plotxy=myf['plotxy']
        ###
        ###
        #Handle globals or user over-ride of arguments
        #
        function_signature_defaults=dict(zip(sdcal.func_code.co_varnames,sdcal.func_defaults))
        for item in function_signature_defaults.iteritems():
		key,val = item
                keyVal = eval(key)
                if (keyVal == None):
                        #user hasn't set it - use global/default
                        pass
                else:
                        #user has set it - use over-ride
			#        sys._getframe(1).f_globals[key]=keyVal
			myf[key]=keyVal
			
	###
	#
        #Add type/menu/range error checking here
#	if type(mask)==str: mask=[mask]
        #
	# LIST PARAMETER NAMES HERE
        arg_names=['infile','telescope','fluxunit','specunit','frame','doppler',\
		   'calmode','scanlist','field','iflist',\
		   'scanaverage','timeaverage','polaverage',\
		   'kernel','kwidth','tau','blmode','blpoly','interactive',\
		   'masklist','sdfile','outform','plotlevel',\
		   'thresh','avg_limit','edge']
	#
	# Assign parameter values
	arg_values = []
	for ptst in arg_names:
		stst = myf[ptst]
		if ( type(stst) == str ):
			# Strip off leading and trailing whitespace
			stst = stst.strip()
			myf[ptst] = stst
			sex = ptst+"='"+stst+"'"
		else:
			sex = ptst+"="+str(stst)
		exec(sex)
		arg_values = arg_values + [stst]
	#
	# LIST ALLOWED PARAMETER TYPES HERE
        arg_types=[str,(str,list),str,str,str,str,str,(list,int),str,(list,int),bool,bool,bool,str,(int,float),(float,int),str,int,bool,list,str,str,int,(float,int),(int,float),(list,int,float)]
	#parameter_printvalues(arg_names,arg_values,arg_types)
	#
	# Check allowed types and values
        try:
                parameter_checktype(arg_names,arg_values,arg_types)
		parameter_checkmenu('outform',outform,['ASCII','SDFITS','MS','MS2','ASAP','ascii','sdfits','ms','ms2','asap'])
		parameter_checkmenu('calmode',calmode,['ps','nod','fs','fsotf','quotient','none'])
		parameter_checkmenu('fluxunit',fluxunit,['K','k','Jy','JY','jy',''])
		parameter_checkmenu('specunit',specunit,['channel','km/s','GHz','MHz','kHz','Hz'])
		parameter_checkmenu('frame',frame,['REST','TOPO','LSRD','LSRK','BARY','GEO','GALACTO','LGROUP','CMB',''])
		parameter_checkmenu('doppler',doppler,['RADIO','OPTICAL','Z','BETA','GAMMA','radio','optical','z',''])
		parameter_checkmenu('blmode',blmode,['auto','list','none'])
		parameter_checkmenu('kernel',kernel,['hanning','gaussian','boxcar','none'])
        except TypeError, e:
                print "sdcal -- TypeError: ", e
		return
        except ValueError, e:
                print "sdcal -- OptionError: ", e
                return
        ###

        ###
        ### Now the actual task code
        ###

	try:
#	    pdb.set_trace()
            # strip off leading and trailing whitespace
	    datfile = infile.strip()

            # load the data with or without averaging
            s=sd.scantable(datfile,scanaverage)
	    if ( abs(plotlevel) > 1 ):
		    # print summary of input data
		    print "Initial Raw Scantable:"
		    print s
	    
	    # Default file name
	    if ( sdfile == '' ):
		    project = datfile + '_cal'
	    else:
		    sdfile = sdfile.strip()
		    project = sdfile

	    # determine current fluxunit
	    fluxunit_now = s.get_fluxunit()
	    print "Current fluxunit = "+fluxunit_now
	    
	    # set flux unit string (be more permissive than ASAP)
	    if ( fluxunit == 'k' ):
		    fluxunit = 'K'
	    elif ( fluxunit == 'JY' or fluxunit == 'jy' ):
		    fluxunit = 'Jy'
	    
	    # fix the fluxunit if necessary
	    # No way to query scantable to find what telescope is
	    # so rely on user input
	    if ( telescope == 'FIX' or telescope == 'fix' ):
			    if ( fluxunit != '' ):
				    if ( fluxunit == fluxunit_now ):
					    print "No need to change default fluxunits"
				    else:
					    s.set_fluxunit(fluxunit)
					    print "Reset default fluxunit to "+fluxunit
					    fluxunit_now = s.get_fluxunit()
			    else:
				    print "Warning - no fluxunit for set_fluxunit"
		    
	    elif ( fluxunit=='' or fluxunit==fluxunit_now ):
		    print "No need to convert fluxunits"
	    elif ( telescope == 'GBT' or telescope == 'gbt' ):
		    print "Convert fluxunit to "+fluxunit
		    # THIS IS THE CHEESY PART
		    # Calculate ap.eff eta at rest freq
		    # Use Ruze law
		    #   eta=eta_0*exp(-(4pi*eps/lambda)**2)
		    # with
		    print "Using GBT parameters"
		    eps = 0.390  # mm
		    eta_0 = 0.71 # at infinite wavelength
		    # Ideally would use a freq in center of
		    # band, but rest freq is what I have
		    rf = s.get_restfreqs()[0]*1.0e-9 # GHz
		    eta = eta_0*pl.exp(-0.001757*(eps*rf)**2)
		    print "Calculated ap.eff. eta = %5.3f " % (eta)
		    print "At rest frequency %5.3f GHz" % (rf)
		    D = 104.9 # 100m x 110m
		    print "Assume phys.diam D = %5.1f m" % (D)
		    s.convert_flux(eta=eta,d=D)
		    
		    print "Successfully converted fluxunit to "+fluxunit
	    elif ( telescope == 'AT' or telescope == 'at' ):
		    s.convert_flux()
	    elif ( type(telescope) == list ):
		    # User input telescope params
		    if ( len(telescope) > 1 ):
			    D = telescope[0]
			    eta = telescope[1]
			    print "Use phys.diam D = %5.1f m" % (D)
			    print "Use ap.eff. eta = %5.3f " % (eta)
			    s.convert_flux(eta=eta,d=D)
		    elif ( len(telescope) > 0 ):
			    jypk = telescope[0]
			    print "Use gain = %6.4f Jy/K " % (jypk)
			    s.convert_flux(jyperk=jypk)
		    else:
			    print "Empty telescope list"
	    else:
		    # Unknown telescope type
		    print "Unknown telescope - cannot convert"

	    # set default spectral axis unit
	    if ( specunit != '' ):
		    s.set_unit(specunit)

	    # reset frame and doppler if needed
	    if ( frame != '' ):
		    s.set_freqframe(frame)
	    else:
		    print 'Using current frequency frame'

	    if ( doppler != '' ):
		    if ( doppler == 'radio' ):
			    ddoppler = 'RADIO'
		    elif ( doppler == 'optical' ):
			    ddoppler = 'OPTICAL'
		    elif ( doppler == 'z' ):
			    ddoppler = 'Z'
		    else:
			    ddoppler = doppler
		    
		    s.set_doppler(ddoppler)
	    else:
		    print 'Using current doppler convention'

	    # Select scan and field
	    sel = sd.selector()
	    
	    # Set up scanlist
	    if ( type(scanlist) == list ):
		    # is a list
		    scans = scanlist
	    else:
		    # is a single int, make into list
		    scans = [ scanlist ]
	    # Now select them
	    if ( len(scans) > 0 ):
		    sel.set_scans(scans)
	    
	    # Select source names
	    if ( field != '' ):
		    sel.set_name(field)
		    # NOTE: currently can only select one
		    # set of names this way, will probably
		    # need to do a set_query eventually

	    # Select IFs
	    if ( type(iflist) == list ):
		    # is a list
		    ifs = iflist
	    else:
		    # is a single int, make into list
		    ifs = [ iflist ]
	    if ( len(ifs) > 0 ):
		    # Do any IF selection
		    sel.set_ifs(ifs)

	    # Apply the selection (if any)
	    s.set_selection(sel)

	    del sel

	    scanns = s.getscannos()
	    sn=list(scanns)
	    print "Number of scans to be processed:", len(sn)

	    if ( calmode == 'nod' ):
		    # nod data
		    print "Calibrating nod data"
		    scal=sd.calnod(s,sn)
	    elif ( (calmode=='fs') or (calmode=='fsotf') ):
		    # frequency-switched data    
		    print "Calibrating frequency-switched data"
		    scal=sd.calfs(s,sn)
	    elif ( calmode == 'ps' ):
		    # position switch (default)
		    print "Calibrating position-switched data"
		    scal=sd.calps(s,sn)
	    elif ( calmode == 'quotient' ):
	            # quotient (for Mopra etc.)
		    print "Calibrating using quotient"
		    scal=s.auto_quotient()
	    else:
		    # no calibration
		    print "No calibration"
		    scal=s.copy()

	    # Done with scantable s - clean up to free memory
	    del s
	    
	    # do opacity (atmospheric optical depth) correction
	    if ( tau > 0.0 ):
		    # recalculate az/el (needed for GBT data)
		    scal.recalc_azel()
		    scal.opacity(tau)

	    # Average in time if desired
	    if ( timeaverage ):
		    stave=sd.average_time(scal,weight='tintsys')
		    del scal
		    # Now average over polarizations;Tsys-weighted (1/Tsys**2) average
		    if ( polaverage ):
			    np = stave.npol()
			    if ( np > 1 ):
				    spave=stave.average_pol(weight='tsys')
			    else:
				    # only single polarization
				    print "Single polarization data - no need to average"
				    spave=stave.copy()
		    else:
			    spave=stave.copy()
		    del stave
	    else:
		    if ( polaverage ):
			    np = scal.npol()
			    if ( np > 1 ):
				    spave=scal.average_pol(weight='tsys')
			    else:
				    # only single polarization
				    print "Single polarization data - no need to average"
				    spave=scal.copy()
		    else:
			    spave=scal.copy()
		    del scal
	    
	    if ( abs(plotlevel) > 1 ):
		    # print summary of calibrated data
		    print "Final Calibrated Scantable:"
		    print spave

	    # Smooth the spectrum (if desired)
	    if ( kernel != 'none' and (not (kwidth<=0 and kernel!='hanning'))):
		    if ( abs(plotlevel) > 1 ):
			    # plot spectrum before smoothing
			    # each IF is separate panel, pols stacked
			    sd.plotter.plot(spave)
			    sd.plotter.set_mode(stacking='p',panelling='i')
			    if ( plotlevel < -1 ):
				    # Hardcopy - currently no way w/o screen display first
				    pltfile=project+'_rawspec.eps'
				    sd.plotter.save(pltfile)
		    print "Smoothing spectrum with kernel "+kernel
		    spave.smooth(kernel,kwidth)

	    # Polynomial baseline (if requested, e.g. blpoly=5)
	    if ( blmode == 'auto' ):
		    # Auto-fit baseline
		    if( blpoly >= 0 ):
			    if (len(masklist) > 0):
				    # Use baseline mask for regions to INCLUDE in baseline fit
				    # Create mask using list, e.g. masklist=[[500,3500],[5000,7500]]
				    basemask=spave.create_mask(masklist)
				    spave.auto_poly_baseline(mask=basemask,order=blpoly,edge=edge,
							     threshold=thresh,chan_avg_limit=avg_limit,
							     plot=interactive)
				    del basemask
			    else:
				    # Automatic polynomial baseline without mask
				    spave.auto_poly_baseline(order=blpoly,edge=edge,threshold=thresh,
							     chan_avg_limit=avg_limit,plot=interactive)
	    elif ( blmode == 'list' ):
		    # Fit in specified region
		    if( blpoly >= 0 ):
			    if (len(masklist) > 0):
				    # Use baseline mask for regions to INCLUDE in baseline fit
				    # Create mask using list, e.g. masklist=[[500,3500],[5000,7500]]
				    basemask=spave.create_mask(masklist)
				    spave.poly_baseline(mask=basemask,order=blpoly,plot=interactive)
				    del basemask
			    else:
				    # Without mask
				    spave.poly_baseline(order=blpoly,plot=interactive)
	    
	    # Plot final spectrum
	    if ( abs(plotlevel) > 0 ):
		    # each IF is separate panel, pols stacked
		    sd.plotter.plot(spave)
		    sd.plotter.set_mode(stacking='p',panelling='i')
		    sd.plotter.set_histogram(hist=True)
		    sd.plotter.axhline(color='r',linewidth=2)
		    if ( plotlevel < 0 ):
			    # Hardcopy - currently no way w/o screen display first
			    pltfile=project+'_calspec.eps'
			    sd.plotter.save(pltfile)
			    
	    
	    # Now save the spectrum and write out final ms
	    if ( (outform == 'ASCII') or (outform == 'ascii') ):
		    outform = 'ASCII'
		    spefile = project + '_'
	    elif ( (outform == 'ASAP') or (outform == 'asap') ):
		    outform = 'ASAP'
		    spefile = project
	    elif ( (outform == 'SDFITS') or (outform == 'sdfits') ):
		    outform = 'SDFITS'
		    spefile = project
	    elif ( (outform == 'MS') or (outform == 'ms') or (outform == 'MS2') or (outform == 'ms2') ):
		    outform = 'MS2'
		    spefile = project
	    else:
		    outform = 'ASAP'
		    spefile = project
		    
		    
	    spave.save(spefile,outform,True)
	    print "Wrote output "+outform+" file "+spefile

	    # Clean up scantable
	    del spave
	    
	    # DONE
	
	except TypeError, e:
                print "sdcal -- TypeError: ", e
		return
        except ValueError, e:
                print "sdcal -- OptionError: ", e
                return
	except Exception, instance:
		print '***Error***',instance
		return

        saveinputs=myf['saveinputs']
        saveinputs('sdcal','sdcal.last')

###End of sdcal

def sdcal_defaults():
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['infile']=''
	myf['telescope']=''
	myf['fluxunit']=''
	myf['specunit']='channel'
	myf['frame']=''
	myf['doppler']=''
	myf['calmode']='none'
        myf['scanlist']=[]
        myf['field']=''
        myf['iflist']=[]
	myf['scanaverage']=False
	myf['timeaverage']=False
	myf['polaverage']=False
        myf['kernel']='none'
        myf['kwidth']=5
        myf['tau']=0.0
        myf['blmode']='none'
        myf['blpoly']=5
        myf['interactive']=False
        myf['masklist']=[]
        myf['sdfile']=''
        myf['outform']='ASAP'
	myf['plotlevel']=0
	myf['thresh']=5
	myf['avg_limit']=4
	myf['edge']=[0]
	
def sdcal_description(key='sdcal'):
	desc={'sdcal': 'ASAP SD calibration task\n',
	'infile': 'name of input SD dataset',
	'telescope': 'name or param of telescope for flux conversion',
	'fluxunit': 'units for line flux (K,Jy) (''=current)',
	'specunit': 'units for spectral axis (channel,km/s,GHz)',
	'frame': 'frequency reference frame, e.g. LSRK (''=current)',
	'doppler': 'doppler convention, e.g. RADIO (''=current)',
	'calmode': 'SD calibration mode (ps,nod,fs,fsotf,none)',
        'scanlist': 'list of scans to use (e.g. [1,2,3,4])',
        'field': 'string for selection by source name (e.g. FLS3a*)',
        'iflist': 'list of IF ids to select (e.g. [0,1])',
	'scanaverage': 'average integs within scans (True,False)',
	'timeaverage': 'average scans over time (True,False)',
	'polaverage': 'average over polarizations (True,False)',
        'kernel': 'type of spectral smoothing (e.g. hanning)',
	'kwidth': 'width of spectral kernel in channels',
        'tau': 'atmospheric optical depth for correction',
	'blmode': 'baseline fitting mode (auto,list,none)',
        'blpoly': 'order of polynomial for baseline fitting (<0=none)',
        'interactive': 'interactive baseline fitting (True/False)',
        'masklist': 'list of channel regions for baseline fitting ([]=auto)',
        'sdfile': 'output file name',
        'outform': 'output file format (ASCII,MS,SDFITS,ASAP)',
	'plotlevel': 'plot results (0=none,1+=some,<0=hardcopy)',
	'thresh': '(Auto) S/N ratio limit for linefinder',
	'avg_limit': '(Auto) max number of chans to average for broad lines',
	'edge': '(Auto) channels to drop from beginning and end',
	}
	return desc[key]

###End of task

def sdfit(sdfile=None,telescope=None,fluxunit=None,specunit=None,frame=None,doppler=None,scanlist=None,field=None,iflist=None,fitmode=None,maskline=None,invertmask=None,nfit=None,fitfile=None,plotlevel=None,thresh=None,min_nchan=None,avg_limit=None,box_size=None,edge=None,xstat=None):
	"""ASAP SD line-fitting task (prototype):

	Task Version 2007-03-08 STM
	Keyword arguments:
        sdfile -- name of input SD dataset
	        default: none - must input file name
		example: 'mysd.asap'
		         See sdcal for allowed formats.
	telescope -- the telescope name or characteristics
	        options: (str) name or (list) list of gain info
		default: '' (none set)
		example: telescope='GBT' for the GBT
		         telescope='AT' for one of the default scopes
			 known to ASAP.
			 telescope=[104.9,0.43] diameter(m), ap.eff.
			 telescope=[0.743] gain in Jy/K
			 telescope='FIX' to change default fluxunit
			 see description below
	fluxunit -- units for line flux
		options: (str) 'K','Jy',''
	        default: '' (keep current fluxunit)
		WARNING: For GBT data, see description below.
	specunit -- units for spectral axis
	        options: (str) 'channel','km/s','GHz','MHz','kHz','Hz',''
	        default: '' (keep current specunit)
	frame -- frequency frame for spectral axis
	        options: (str) 'LSRK','REST','TOPO','LSRD','BARY',
	                 'GEO','GALACTO','LGROUP','CMB'
		default: currently set frame in scantable
		WARNING: frame='REST' not yet implemented
	doppler -- doppler mode
	        options: (str) 'RADIO','OPTICAL','Z','BETA','GAMMA'
	        default: currently set doppler in scantable
	scanlist -- list of scan numbers to process
	        default: [] (use all scans)
		example: [21,22,23,24]
	field -- selection string for selecting scans by name
	        default: '' (no name selection)
		example: 'FLS3a*'
		         this selection is in addition to scanlist
			 and iflist
	iflist -- list of IF id numbers to select
	        default: [] (use all IFs)
		example: [15]
	fitmode -- mode for fitting
	        options: (str) 'list','auto'
		default: 'auto'
		example: 'list' will use maskline to define regions to
	                        fit for lines with nfit in each
	                 'auto' will use the linefinder to fit for lines
			        using autopars (see below)
	maskline -- list of mask regions to INCLUDE in LINE fitting
	        default: all
		example: maskline=[[3900,4300]] for a single region, or
		         maskline=[[3900,4300],[5000,5400]] for two, etc.
	invertmask -- invert mask (EXCLUDE masklist instead)
	        options: (bool) True, False
	        default: False
		example: invertmask=True, then will make one region that is
		         the exclusion of the maskline regions
	nfit -- list of number of gaussian lines to fit in in maskline region
	        default: 0 (no fitting)
		example: nfit=[1] for single line in single region,
		         nfit=[2] for two lines in single region,
	                 nfit=[1,1] for single lines in each of two regions, etc.
	fitfile -- name of output file for fit results
	        default: no output fit file
		example: 'mysd.fit'
	plotlevel -- control for plotting of results
	        options: (int) 0=none, 1=some, 2=more
		default: 0 (no plotting)
		example: plotlevel=1 plots fit and residual
		         no hardcopy available for fitter
		WARNING: be careful plotting OTF data with lots of fields

	AUTOPARS: the following parameters are used for fitmode='auto' ONLY
	-------------------------------------------------------------------
	thresh -- S/N threshold for linefinder
	        default: 5
		example: a single channel S/N ratio above which the channel is
	                 considered to be a detection
        min_nchan -- minimum number of consecutive channels for linefinder
	        default: 3
                example: minimum number of consecutive channels required to pass threshold
	avg_limit -- channel averaging for broad lines
	        default: 4
		example: a number of consequtive channels not greater than
                         this parameter can be averaged to search for broad lines
	box_size -- running mean box size
	        default: 0.2
		example: a running mean box size specified as a fraction
                         of the total spectrum length
	edge -- channels to drop at beginning and end of spectrum
	        default: 0
		example: [1000] drops 1000 channels at beginning AND end
		         [1000,500] drops 1000 from beginning and 500 from end

                Note: For bad baselines threshold should be increased,
                and avg_limit decreased (or even switched off completely by
                setting this parameter to 1) to avoid detecting baseline
                undulations instead of real lines.

        -------------------------------------------------------------------
	xstat -- RETURN ONLY: a Python dictionary of line statistics
	        keys:    'peak','cent','fwhm'
		example: each value is a list of lists with one list of
		         2 entries [fitvalue,error] per component.
			 e.g. xstat['peak']=[[234.9, 4.8],[234.2, 5.3]]
			 for 2 components.


        DESCRIPTION:

	Task sdfit is a basic line-fitter for single-dish spectra.
	It assumes that the spectra have been calibrated in sdcal.
	Furthermore, it assumes that any selection of scans, IFs,
	polarizations, and time and channel averaging/smoothing has
	also already been done (in sdcal) as there are no controls
	for these.  Note that you can run sdcal with calmode='none'
	and do selection, writing out a new scantable.

	Note that multiple scans and IFs can in principle be handled, but
	we recommend that you use scanlist, field, and iflist to give a
	single selection for each fit.

	For complicated spectra, sdfit does not do a good job of
	"auto-guessing" the starting model for the fit.  We recommend
	you use sd.fitter in the toolkit which has more options, such
	as fixing components in the fit and supplying starting guesses
	by hand.

	WARNING: sdfit will currently return the fit for the first
	row in the scantable.  Does not handle multiple polarizations.

	ASAP recognizes the data of the "AT" telescopes, but currently
	does not know about the GBT or any other telescope. This task
	does know about GBT.  Therefore, if you wish to change the
	fluxunit (see below), then you need to tell it what to do.  If
	you set telescope='AT' it will use its internal defaults.  If
	you set telescope 'GBT', it will use an approximate aperture
	efficiency conversion.  See sdcal description for information
	on fluxunit conversion.
	
	"""
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['__last_taskname']='sdfit'
        tb=myf['tb']
        im = myf['im']
	qa = myf['qa']
        me = myf['me']
	sm = myf['sm']
        ia = myf['ia']
        cl = myf['cl']
        pl = myf['pl']
	rg = myf['rg']
	default=myf['default']
	plotxy=myf['plotxy']
        ###
        ###
        #Handle globals or user over-ride of arguments
        #
        function_signature_defaults=dict(zip(sdfit.func_code.co_varnames,sdfit.func_defaults))
        for item in function_signature_defaults.iteritems():
		key,val = item
                keyVal = eval(key)
                if (keyVal == None):
                        #user hasn't set it - use global/default
                        pass
                else:
                        #user has set it - use over-ride
			#        sys._getframe(1).f_globals[key]=keyVal
			myf[key]=keyVal
			
	###
	#
        #Add type/menu/range error checking here
#	if type(mask)==str: mask=[mask]
        #
	# LIST PARAMETER NAMES HERE
        arg_names=['sdfile','telescope','fluxunit','specunit','frame','doppler',\
		   'scanlist','field','iflist','fitmode','maskline','invertmask',\
		   'nfit','fitfile','plotlevel',\
		   'thresh','min_nchan','avg_limit','box_size','edge','xstat']
	#
	# Assign parameter values
	arg_values = []
	for ptst in arg_names:
		stst = myf[ptst]
		if ( type(stst) == str ):
			# Strip off leading and trailing whitespace
			stst = stst.strip()
			myf[ptst] = stst
			sex = ptst+"='"+stst+"'"
		else:
			sex = ptst+"="+str(stst)
		exec(sex)
		arg_values = arg_values + [stst]
	#
	# LIST ALLOWED PARAMETER TYPES HERE
        arg_types=[str,(str,list),str,str,str,str,(list,int),str,(int,list),str,\
		   list,bool,(int,list),str,int,(float,int),(int,float),\
		   (int,float),(float,int),(list,int),dict]
	#parameter_printvalues(arg_names,arg_values,arg_types)
	#
	# Check allowed types and values
        try:
                parameter_checktype(arg_names,arg_values,arg_types)
		parameter_checkmenu('fitmode',fitmode,['list','auto'])
		parameter_checkmenu('fluxunit',fluxunit,['K','k','Jy','JY','jy',''])
		parameter_checkmenu('specunit',specunit,['channel','km/s','GHz','MHz','kHz','Hz',''])
		parameter_checkmenu('frame',frame,['REST','TOPO','LSRD','LSRK','BARY','GEO','GALACTO','LGROUP','CMB',''])
		parameter_checkmenu('doppler',doppler,['RADIO','OPTICAL','Z','BETA','GAMMA','radio','optical','z',''])
        except TypeError, e:
                print "sdfit -- TypeError: ", e
		return
        except ValueError, e:
                print "sdfit -- OptionError: ", e
                return
        ###

        ###
        ### Now the actual task code
        ###

	try:
#	    pdb.set_trace()
            #load the data without averaging
            s=sd.scantable(sdfile,False)
	    
	    # determine current fluxunit
	    fluxunit_now = s.get_fluxunit()
	    print "Current fluxunit = "+fluxunit_now
	    
	    # set flux unit string (be more permissive than ASAP)
	    if ( fluxunit == 'k' ):
		    fluxunit = 'K'
	    elif ( fluxunit == 'JY' or fluxunit == 'jy' ):
		    fluxunit = 'Jy'
	    
	    # fix the fluxunit if necessary
	    # No way to query scantable to find what telescope is
	    # so rely on user input
	    if ( telescope == 'FIX' or telescope == 'fix' ):
			    if ( fluxunit != '' ):
				    if ( fluxunit == fluxunit_now ):
					    print "No need to change default fluxunits"
				    else:
					    s.set_fluxunit(fluxunit)
					    print "Reset default fluxunit to "+fluxunit
					    fluxunit_now = s.get_fluxunit()
			    else:
				    print "Warning - no fluxunit for set_fluxunit"
		    
	    elif ( fluxunit=='' or fluxunit==fluxunit_now ):
		    print "No need to convert fluxunits"
	    elif ( telescope == 'GBT' or telescope == 'gbt' ):
		    print "Convert fluxunit to "+fluxunit
		    # THIS IS THE CHEESY PART
		    # Calculate ap.eff eta at rest freq
		    # Use Ruze law
		    #   eta=eta_0*exp(-(4pi*eps/lambda)**2)
		    # with
		    print "Using GBT parameters"
		    eps = 0.390  # mm
		    eta_0 = 0.71 # at infinite wavelength
		    # Ideally would use a freq in center of
		    # band, but rest freq is what I have
		    rf = s.get_restfreqs()[0]*1.0e-9 # GHz
		    eta = eta_0*pl.exp(-0.001757*(eps*rf)**2)
		    print "Calculated ap.eff. eta = %5.3f " % (eta)
		    print "At rest frequency %5.3f GHz" % (rf)
		    D = 104.9 # 100m x 110m
		    print "Assume phys.diam D = %5.1f m" % (D)
		    s.convert_flux(eta=eta,d=D)
		    
		    print "Successfully converted fluxunit to "+fluxunit
	    elif ( telescope == 'AT' or telescope == 'at' ):
		    s.convert_flux()
	    elif ( type(telescope) == list ):
		    # User input telescope params
		    if ( len(telescope) > 1 ):
			    D = telescope[0]
			    eta = telescope[1]
			    print "Use phys.diam D = %5.1f m" % (D)
			    print "Use ap.eff. eta = %5.3f " % (eta)
			    s.convert_flux(eta=eta,d=D)
		    elif ( len(telescope) > 0 ):
			    jypk = telescope[0]
			    print "Use gain = %6.4f Jy/K " % (jypk)
			    s.convert_flux(jyperk=jypk)
		    else:
			    print "Empty telescope list"
	    else:
		    # Unknown telescope type
		    print "Unknown telescope - cannot convert"

	    # set spectral axis unit
	    if ( specunit != '' ):
		    s.set_unit(specunit)
		    specunit_now = specunit
	    else:
		    specunit_now = s.get_unit()

	    # reset frame and doppler if needed
	    if ( frame != '' ):
		    s.set_freqframe(frame)
	    else:
		    print 'Using current frequency frame'

	    if ( doppler != '' ):
		    if ( doppler == 'radio' ):
			    ddoppler = 'RADIO'
		    elif ( doppler == 'optical' ):
			    ddoppler = 'OPTICAL'
		    elif ( doppler == 'z' ):
			    ddoppler = 'Z'
		    else:
			    ddoppler = doppler
		    
		    s.set_doppler(ddoppler)
	    else:
		    print 'Using current doppler convention'

	    # Prepare a selection
	    sel=sd.selector()

	    # Scan selection
	    if ( type(scanlist) == list ):
		    # is a list
		    scans = scanlist
	    else:
		    # is a single int, make into list
		    scans = [ scanlist ]
	    if ( len(scans) > 0 ):
		    sel.set_scans(scans)
	    
	    # Select source names
	    if ( field != '' ):
		    sel.set_name(field)
		    # NOTE: currently can only select one
		    # set of names this way, will probably
		    # need to do a set_query eventually

	    # IF selection
	    if ( type(iflist) == list ):
		    # is a list
		    ifs = iflist
	    else:
		    # is a single int, make into list
		    ifs = [ iflist ]
	    if ( len(ifs) > 0 ):
		    # Do any IF selection
		    sel.set_ifs(ifs)

	    # Apply the selection (if any)
	    s.set_selection(sel)

	    del sel

	    # Make line region masks and list of line regions
	    if ( fitmode == 'list' ):
		    # Assume the user has given a list of lines
		    # e.g. maskline=[[3900,4300]] for a single line
		    if ( len(maskline) > 0 ):
			    # There is a user-supplied channel mask for lines
			    if ( not invertmask ):
				    linemask=s.create_mask(maskline)
				    domask = True
				    doguess = True
				    # Make sure this is a list-of-lists
				    # (e.g. [[1,10],[20,30]])
				    if ( type(maskline[0]) == list ):
					    ll = len(maskline)
					    nlines = len(maskline)
					    linelist = maskline
				    else:
					    # Not a list, turn into one
					    nlines = len(maskline)/2
					    linelist = []
					    for i in range(nlines):
						    lo = maskline[2*i]
						    up = maskline[2*i+1]
						    linelist = linelist + [[lo,up]]
			    else:
				    # invert regions
				    linemask=s.create_mask(maskline,invert=True)
				    domask = True
				    nlines = 1
				    linelist=[]
				    doguess = False
		    else:
			    # Use whole region
			    nlines = 1
			    linelist=[]
			    doguess = True
			    domask = False

		    print "Identified ",nlines," regions for fitting"
		    if ( doguess ):
			    "Will use these as starting guesses"
		    else:
			    "No starting guesses available"
		    
	    else:
		    # Fit mode AUTO
		    print "Trying AUTO mode"
		    if ( len(maskline) > 0 ):
			    # There is a user-supplied channel mask for lines
			    linemask=s.create_mask(maskline,invert=invertmask)
			    domask = True
		    else:
			    domask = False
		    # A single region
		    nlines = 1
		    doguess = False
	    
	    # If we have line regions, get starting guesses
	    linemax=[]
	    lineeqw=[]
	    linecen=[]
	    if ( doguess ):
		    # For each line get guess of max, cen and estimated
		    # equivalent width (sum/max)
		    if( len(linelist) > 0):
			    for x in linelist:
				    msk = s.create_mask(x)
				    maxl = s.stats('max',msk)
				    suml = s.stats('sum',msk)
				    if ( maxl[0] != 0.0 ):
					    eqw = suml[0]/maxl[0]
				    else:
					    eqw = 0.0
				    linemax = linemax + [maxl[0]]
				    lineeqw = lineeqw + [eqw]
				    cen = 0.5*(x[0] + x[1])
				    linecen = linecen + [cen]
				    del msk
		    else:
			    # For what its worth, do stats on unmasked region
			    if ( domask ):
				    maxl = s.stats('max',linemask)
				    suml = s.stats('sum',linemask)
			    else:
				    maxl = s.stats('max')
				    suml = s.stats('sum')
				    
			    if ( maxl[0] != 0.0 ):
				    eqw = suml[0]/maxl[0]
			    else:
				    eqw = 0.0
			    linemax = linemax + [maxl[0]]
			    lineeqw = lineeqw + [eqw]
			    # I dont know a better way to specify the center of the spectrum
			    cen = s.nchan/2
			    linecen = linecen + [cen]
						    
	    # Now the line fitting
	    if (nlines > 0):
		    if ( fitmode == 'list' ):
			    # Get number of things to fit from nfit list
			    if ( type(nfit) == list ):
				    numfit = len(nfit)
				    nlist = nfit
			    else:
				    numfit = 1
				    nlist = [nfit]
			    # Drop extra over nlines
			    numfit = min(numfit,nlines)
		    
			    ncomps = 0
			    comps = []
			    for i in range(numfit):
				    ncomps = ncomps + nlist[i]
				    comps += [nlist[i]]
		    else:
			    # Auto mode - one comp per line region
			    # Overwriting user-supplied nfit
			    numfit = nlines
			    ncomps = nlines
			    comps = []
			    for i in range(nlines):
				    comps += [1]

		    print "Will fit ",ncomps," components in ",numfit," regions"
		    if (numfit > 0):
			    # Fit the line using numfit gaussians
			    # Assume the nfit list matches maskline
			    f=sd.fitter()
			    f.set_function(gauss=ncomps)
			    if ( domask ):
				    f.set_scan(s,linemask)
			    else:
				    f.set_scan(s)
			    if ( doguess ):
				    # Guesses using max, cen, and eqw
				    # NOTE: should there be user options here?
				    n = 0
				    for i in range(numfit):
					    if ( comps[i] == 1 ):
						    # use guess
						    maxl = linemax[i]
						    fwhm = 0.7*abs( lineeqw[i] )
						    cenl = linecen[i]
						    f.set_gauss_parameters(maxl, cenl, fwhm,\
									   component=n)
						    n = n + 1
					    else:
						    # cannot guess for multiple comps yet
						    n = n + comps[i]
			    else:
				    # No guesses
				    print "Fitting lines without starting guess"
			    
			    # Now fit
			    f.fit()
			    fstat = f.get_parameters()

			    # Check for convergence
			    goodfit = ( len(fstat['errors']) > 0 )
			    if ( goodfit ):
				    # Plot residuals
				    if ( abs(plotlevel) > 0):
					    # each IF is separate panel, pols stacked
					    f.plot(residual=True)
				    
				    # Retrieve fit parameters
				    xstat['nfit'] = ncomps
				    xstat['peak'] = []
				    xstat['cent'] = []
				    xstat['fwhm'] = []
				    for i in range(ncomps):
				       	    fstat = f.get_parameters(i)
					    xstat['peak'] = xstat['peak'] \
							    + [[fstat['params'][0],\
							       fstat['errors'][0]]]
					    xstat['cent'] = xstat['cent'] \
							    + [[fstat['params'][1],\
							       fstat['errors'][1]]]
					    xstat['fwhm'] = xstat['fwhm'] \
							    + [[fstat['params'][2],\
							       fstat['errors'][2]]]
				    # Store fit
				    if ( fitfile != '' ):
					    f.store_fit(fitfile)
			    else:
				    # Did not converge
				    print "ERROR: Fit failed to converge"
				    xstat['nfit'] = -ncomps
			    
			    # Clean up
			    del f

	    # Assign variables for return
	    myf['xstat']=xstat
	    
	    # Final clean up
	    del s
	    if (domask): del linemask
	    
	    # DONE
	
	except TypeError, e:
                print "sdfit -- TypeError: ", e
		return
        except ValueError, e:
                print "sdfit -- OptionError: ", e
                return
	except Exception, instance:
		print '***Error***',instance
		return

        saveinputs=myf['saveinputs']
        saveinputs('sdfit','sdfit.last')

###End of sdfit

def sdfit_defaults():
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['sdfile']=''
	myf['telescope']=''
	myf['fluxunit']=''
	myf['specunit']=''
	myf['frame']=''
	myf['doppler']=''
        myf['scanlist']=[]
        myf['field']=''
        myf['iflist']=[]
	myf['fitmode']='auto'
        myf['maskline']=[]
	myf['invertmask']=False
        myf['nfit']=0
        myf['fitfile']=''
	myf['plotlevel']=0
	myf['thresh']=5
	myf['min_nchan']=3
	myf['avg_limit']=4
	myf['box_size']=0.2
	myf['edge']=[0]
	myf['xstat']={}
	
def sdfit_description(key='sdfit'):
	desc={'sdfit': 'ASAP SD line fitting task\n',
	'sdfile': 'name of input SD dataset',
	'telescope': 'name or param of telescope for flux conversion',
	'fluxunit': 'units for line flux (K,Jy,def=current)',
	'specunit': 'units for spectral axis (def=current)',
	'frame': 'frequency reference frame, e.g. LSRK (def=current)',
	'doppler': 'doppler convention, e.g. RADIO (def=current)',
        'scanlist': 'list of scans to use (e.g. [1,2,3,4])',
        'field': 'string for selection by source name (e.g. FLS3a*)',
        'iflist': 'list of IF ids to select (e.g. [0,1])',
	'fitmode': 'SD line fitting mode (auto,list)',
        'maskline': 'list of regions for line fitting (e.g. [[2500,3000],[3500,4500]])',
	'invertmask': 'invert maskline regions (True,False)',
        'nfit': 'list of number of Gaussians to fit in each line region (e.g. [2,1])',
        'fitfile': 'output file name for fit results',
	'plotlevel': 'control for plotting of results (0=none,1+=some)',
	'thresh': '(Auto) S/N ratio limit for linefinder',
	'min_nchan': '(Auto) min number of consecutive chans above threshold',
	'avg_limit': '(Auto) max number of chans to average for broad lines',
	'box_size': '(Auto) running mean box size as fraction of spectrum',
	'edge': '(Auto) channels to drop from beginning and end',
	'xstat': '(RETURN ONLY) Python dictionary with line fitting output',
	}
	return desc[key]

###End of task

def sdstat(sdfile=None,telescope=None,fluxunit=None,specunit=None,frame=None,doppler=None,scanlist=None,field=None,iflist=None,masklist=None,invertmask=None,xstat=None):
	"""ASAP SD spectral region statistics task (prototype):

	Task Version 2007-03-06 STM
	Keyword arguments:
        sdfile -- name of input SD dataset
	        default: none - must input file name
		example: 'mysd.asap'
		         See sdcal for allowed formats.
	telescope -- the telescope name or characteristics
	        options: (str) name or (list) list of gain info
		default: '' (none set)
		example: telescope='GBT' for the GBT
		         telescope='AT' for one of the default scopes
			 known to ASAP.
			 telescope=[104.9,0.43] diameter(m), ap.eff.
			 telescope=[0.743] gain in Jy/K
			 telescope='FIX' to change default fluxunit
			 see description below
	fluxunit -- units for line flux
	        options: 'K','Jy',''
		default: '' (keep current fluxunit)
		WARNING: For GBT data, see description below.
	specunit -- units for spectral axis
	        options: 'channel','km/s','GHz','MHz','kHz','Hz',''
		default: '' (keep current specunit)
		example: make sure this is the same units of masklist
	frame -- frequency frame for spectral axis
	        options: (str) 'LSRK','REST','TOPO','LSRD','BARY',
	                 'GEO','GALACTO','LGROUP','CMB'
		default: currently set frame in scantable
	doppler -- doppler mode
	        options: (str) 'RADIO','OPTICAL','Z','BETA','GAMMA'
	        default: currently set doppler in scantable
	scanlist -- list of scan numbers to process
	        default: [] (use all scans)
		example: [21,22,23,24]
	field -- selection string for selecting scans by name
	        default: '' (no name selection)
		example: 'FLS3a*'
		         this selection is in addition to scanlist
			 and iflist
	iflist -- list of IF id numbers to select
	        default: [] (use all IFs)
		example: [15]
	masklist -- list of mask regions to INCLUDE in stats
	        default: [] (whole spectrum)
		example: [4000,4500] for one region
		         [[1000,3000],[5000,7000]]
			 these must be pairs of [lo,hi] boundaries
	invertmask -- invert mask (EXCLUDE masklist instead)
	        options: (bool) True,False
		default: false

        -------------------------------------------------------------------
	xstat -- RETURN ONLY: a Python dictionary of line statistics
	           keys: 'rms','stddev','max','min','sum','median','mean',
		         'eqw'
		example: print "rms = ",xstat['rms']
		         these can be used for testing in scripts or
			 for regression
			 
			 'eqw' is equivalent width (sum/mag) where mag
			 is either max or min depending on which has
			 greater magnitude.


        DESCRIPTION:

	Task sdstat computes basic statistics (rms,mean,median,sum)
	for single-dish spectra.  It assumes that the spectra have
	been calibrated in sdcal.  Furthermore, it assumes that any
	selection of scans, IFs, polarizations, and time and channel
	averaging/smoothing has also already been done (in sdcal) as
	there are no controls for these.  Note that you can run sdcal
	with calmode='none' and do selection, writing out a new
	scantable.
	
	Note that multiple scans and IFs can in principle be handled, but
	we recommend that you use scanlist, field, and iflist to give a
	single selection for each run.

	WARNING: If you do have multiple scantable rows, then xstat
	values will be lists.

	ASAP recognizes the data of the "AT" telescopes, but currently
	does not know about the GBT or any other telescope. This task
	does know about GBT.  Therefore, if you wish to change the
	fluxunit (see below), then you need to tell it what to do.  If
	you set telescope='AT' it will use its internal defaults.  If
	you set telescope 'GBT', it will use an approximate aperture
	efficiency conversion.  See sdcal description for information
	on fluxunit conversion.
	
	"""
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['__last_taskname']='sdstat'
        tb=myf['tb']
        im = myf['im']
	qa = myf['qa']
        me = myf['me']
	sm = myf['sm']
        ia = myf['ia']
        cl = myf['cl']
        pl = myf['pl']
	rg = myf['rg']
	default=myf['default']
	plotxy=myf['plotxy']
        ###
        ###
        #Handle globals or user over-ride of arguments
        #
        function_signature_defaults=dict(zip(sdstat.func_code.co_varnames,sdstat.func_defaults))
        for item in function_signature_defaults.iteritems():
		key,val = item
                keyVal = eval(key)
                if (keyVal == None):
                        #user hasn't set it - use global/default
                        pass
                else:
                        #user has set it - use over-ride
			#        sys._getframe(1).f_globals[key]=keyVal
			myf[key]=keyVal
	
	###
	#
        #Add type/menu/range error checking here
        #
	# LIST PARAMETER NAMES HERE
        arg_names=['sdfile','telescope','fluxunit','specunit','frame','doppler','scanlist','field','iflist','masklist','invertmask','xstat']
	#
	# Assign parameter values
	arg_values = []
	for ptst in arg_names:
		stst = myf[ptst]
		if ( type(stst) == str ):
			# Strip off leading and trailing whitespace
			stst = stst.strip()
			myf[ptst] = stst
			sex = ptst+"='"+stst+"'"
		else:
			sex = ptst+"="+str(stst)
		exec(sex)
		arg_values = arg_values + [stst]
	#
	# LIST ALLOWED PARAMETER TYPES HERE
        arg_types=[str,(str,list),str,str,str,str,(list,int),str,(int,list),list,bool,dict]
	#parameter_printvalues(arg_names,arg_values,arg_types)
        try:
                parameter_checktype(arg_names,arg_values,arg_types)
		parameter_checkmenu('fluxunit',fluxunit,['K','k','Jy','JY','jy',''])
		parameter_checkmenu('specunit',specunit,['channel','km/s','GHz','MHz','kHz','Hz',''])
		parameter_checkmenu('frame',frame,['REST','TOPO','LSRD','LSRK','BARY','GEO','GALACTO','LGROUP','CMB',''])
		parameter_checkmenu('doppler',doppler,['RADIO','OPTICAL','Z','BETA','GAMMA','radio','optical','z',''])
        except TypeError, e:
                print "sdstat -- TypeError: ", e
		return
        except ValueError, e:
                print "sdstat -- OptionError: ", e
                return
        ###

        ###
        ### Now the actual task code
        ###

	try:
#	    pdb.set_trace()
            #load the data without averaging
            s=sd.scantable(sdfile,False)
	    
	    # determine current fluxunit
	    fluxunit_now = s.get_fluxunit()
	    print "Current fluxunit = "+fluxunit_now
	    
	    # set flux unit string (be more permissive than ASAP)
	    if ( fluxunit == 'k' ):
		    fluxunit = 'K'
	    elif ( fluxunit == 'JY' or fluxunit == 'jy' ):
		    fluxunit = 'Jy'
	    
	    # fix the fluxunit if necessary
	    # No way to query scantable to find what telescope is
	    # so rely on user input
	    if ( telescope == 'FIX' or telescope == 'fix' ):
			    if ( fluxunit != '' ):
				    if ( fluxunit == fluxunit_now ):
					    print "No need to change default fluxunits"
				    else:
					    s.set_fluxunit(fluxunit)
					    print "Reset default fluxunit to "+fluxunit
					    fluxunit_now = s.get_fluxunit()
			    else:
				    print "Warning - no fluxunit for set_fluxunit"
		    
	    elif ( fluxunit=='' or fluxunit==fluxunit_now ):
		    print "No need to convert fluxunits"
	    elif ( telescope == 'GBT' or telescope == 'gbt' ):
		    print "Convert fluxunit to "+fluxunit
		    # THIS IS THE CHEESY PART
		    # Calculate ap.eff eta at rest freq
		    # Use Ruze law
		    #   eta=eta_0*exp(-(4pi*eps/lambda)**2)
		    # with
		    print "Using GBT parameters"
		    eps = 0.390  # mm
		    eta_0 = 0.71 # at infinite wavelength
		    # Ideally would use a freq in center of
		    # band, but rest freq is what I have
		    rf = s.get_restfreqs()[0]*1.0e-9 # GHz
		    eta = eta_0*pl.exp(-0.001757*(eps*rf)**2)
		    print "Calculated ap.eff. eta = %5.3f " % (eta)
		    print "At rest frequency %5.3f GHz" % (rf)
		    D = 104.9 # 100m x 110m
		    print "Assume phys.diam D = %5.1f m" % (D)
		    s.convert_flux(eta=eta,d=D)
		    
		    print "Successfully converted fluxunit to "+fluxunit
	    elif ( telescope == 'AT' or telescope == 'at' ):
		    s.convert_flux()
	    elif ( type(telescope) == list ):
		    # User input telescope params
		    if ( len(telescope) > 1 ):
			    D = telescope[0]
			    eta = telescope[1]
			    print "Use phys.diam D = %5.1f m" % (D)
			    print "Use ap.eff. eta = %5.3f " % (eta)
			    s.convert_flux(eta=eta,d=D)
		    elif ( len(telescope) > 0 ):
			    jypk = telescope[0]
			    print "Use gain = %6.4f Jy/K " % (jypk)
			    s.convert_flux(jyperk=jypk)
		    else:
			    print "Empty telescope list"
	    else:
		    # Unknown telescope type
		    print "Unknown telescope - cannot convert"

	    # set default spectral axis unit
	    if ( specunit != '' ):
		    s.set_unit(specunit)
	    
	    # reset frame and doppler if needed
	    if ( frame != '' ):
		    s.set_freqframe(frame)
	    else:
		    print 'Using current frequency frame'

	    if ( doppler != '' ):
		    if ( doppler == 'radio' ):
			    ddoppler = 'RADIO'
		    elif ( doppler == 'optical' ):
			    ddoppler = 'OPTICAL'
		    elif ( doppler == 'z' ):
			    ddoppler = 'Z'
		    else:
			    ddoppler = doppler
		    
		    s.set_doppler(ddoppler)
	    else:
		    print 'Using current doppler convention'

	    # Prepare a selection
	    sel=sd.selector()

	    # Scan selection
	    if ( type(scanlist) == list ):
		    # is a list
		    scans = scanlist
	    else:
		    # is a single int, make into list
		    scans = [ scanlist ]
	    if ( len(scans) > 0 ):
		    sel.set_scans(scans)
	    
	    # Select source names
	    if ( field != '' ):
		    sel.set_name(field)
		    # NOTE: currently can only select one
		    # set of names this way, will probably
		    # need to do a set_query eventually

	    # IF selection
	    if ( type(iflist) == list ):
		    # is a list
		    ifs = iflist
	    else:
		    # is a single int, make into list
		    ifs = [ iflist ]
	    if ( len(ifs) > 0 ):
		    # Do any IF selection
		    sel.set_ifs(ifs)

	    # Apply the selection (if any)
	    s.set_selection(sel)

	    # set the mask region
	    if ( len(masklist) > 0):
		    msk=s.create_mask(masklist,invert=invertmask)
		    
		    maxl=s.stats('max',msk)
		    minl=s.stats('min',msk)
		    suml=s.stats('sum',msk)
		    meanl=s.stats('mean',msk)
		    medianl=s.stats('median',msk)
		    rmsl=s.stats('rms',msk)
		    stdvl=s.stats('stddev',msk)
		    del msk
	    else:
		    # Full region
		    print 'Using full region'
		    
		    maxl=s.stats('max')
		    minl=s.stats('min')
		    suml=s.stats('sum')
		    meanl=s.stats('mean')
		    medianl=s.stats('median')
		    rmsl=s.stats('rms')
		    stdvl=s.stats('stddev')

	    # Put into output dictionary
	    ns = len(maxl)
	    if ( ns > 1 ):
		    # User selected multiple scans,ifs
		    # put into lists
		    xstat['rms']=list(rmsl)
		    xstat['stddev']=list(stdvl)
		    xstat['max']=list(maxl)
		    xstat['min']=list(minl)
		    xstat['sum']=list(suml)
		    xstat['median']=list(medianl)
		    xstat['mean']=list(meanl)

		    eqw=[]
		    for i in range(ns):
			    if ( maxl[i] != 0.0 or minl[i] != 0.0 ):
				    if ( abs(maxl[i]) >= abs(minl[i]) ):
					    eqw = eqw + [ suml[i]/maxl[i] ]
				    else:
					    eqw = eqw + [ suml[i]/minl[i] ]
			    else:
				    eqw = eqw + [0.0]
					    
		    xstat['eqw']=eqw
	    else:
		    # Single scantable only
		    # put into scalars
		    xstat['rms']=rmsl[0]
		    xstat['stddev']=stdvl[0]
		    xstat['max']=maxl[0]
		    xstat['min']=minl[0]
		    xstat['sum']=suml[0]
		    xstat['median']=medianl[0]
		    xstat['mean']=meanl[0]

		    # Construct equivalent width (=sum/max)
		    if ( maxl[0] != 0.0 or minl[0] != 0.0 ):
			    if ( abs(maxl[0]) >= abs(minl[0]) ):
				    eqw = suml[0]/maxl[0]
			    else:
				    eqw = suml[0]/minl[0]
		    else:
			    eqw = 0.0
		    xstat['eqw']=eqw

	    # Assign variables for return
	    myf['xstat']=xstat
	    
	    # Final clean up
	    del s
	    
	    # DONE
	
	except TypeError, e:
                print "sdstat -- TypeError: ", e
		return
        except ValueError, e:
                print "sdstat -- OptionError: ", e
                return
	except Exception, instance:
		print '***Error***',instance
		return

        saveinputs=myf['saveinputs']
        saveinputs('sdstat','sdstat.last')

###End of sdstat

def sdstat_defaults():
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['sdfile']=''
	myf['telescope']=''
	myf['fluxunit']=''
	myf['specunit']=''
	myf['frame']=''
	myf['doppler']=''
        myf['scanlist']=[]
        myf['field']=''
        myf['iflist']=[]
        myf['masklist']=[]
	myf['invertmask']=False
	myf['xstat']={}
	
def sdstat_description(key='sdstat'):
	desc={'sdstat': 'ASAP SD spectrum statistics task\n',
	'sdfile': 'name of input SD dataset',
	'telescope': 'name or param of telescope for flux conversion',
	'fluxunit': 'units for line flux (K,Jy) (def=current)',
	'specunit': 'units for spectral axis (def=current)',
	'frame': 'frequency reference frame, e.g. LSRK (def=current)',
	'doppler': 'doppler convention, e.g. RADIO (def=current)',
        'scanlist': 'list of scans to use (e.g. [1,2,3,4])',
        'field': 'string for selection by source name (e.g. FLS3a*)',
        'iflist': 'list of IF ids to select (e.g. [0,1])',
        'masklist': 'list of channel [lo,hi] regions for line fitting',
	'invertmask': 'invert masklist (exclude regions) (True,False)',
	'xstat': '(RETURN ONLY) Python dictionary with line fitting output',
	}
	return desc[key]

###End of task

def sdplot(sdfile=None,telescope=None,fluxunit=None,specunit=None,frame=None,doppler=None,scanlist=None,field=None,iflist=None,scanaverage=None,timeaverage=None,polaverage=None,kernel=None,kwidth=None,stack=None,panel=None,flrange=None,sprange=None,linecat=None,linedop=None,plotfile=None):
	"""ASAP SD plotting task (prototype):

	Task Version 2007-03-06 STM
	Keyword arguments:
        sdfile -- name of input SD dataset
	        default: none - must input file name
		example: 'mysd.asap'
		         See sdcal for allowed formats.
	telescope -- the telescope name or characteristics
	        options: (str) name or (list) list of gain info
		default: '' (none set)
		example: telescope='GBT' for the GBT
		         telescope='AT' for one of the default scopes
			 known to ASAP.
			 telescope=[104.9,0.43] diameter(m), ap.eff.
			 telescope=[0.743] gain in Jy/K
			 telescope='FIX' to change default fluxunit
			 see description below
	fluxunit -- units for line flux
	        options: 'K','Jy',''
		default: '' (keep current unit)
		WARNING: For GBT data, see description below.
	specunit -- units for spectral axis
	        options: (str) 'channel','km/s','GHz','MHz','kHz','Hz',''
	        default: '' (keep current specunit)
		example: this will be the units for masklist
	frame -- frequency frame for spectral axis
	        options: (str) 'LSRK','REST','TOPO','LSRD','BARY',
	                 'GEO','GALACTO','LGROUP','CMB'
		default: currently set frame in scantable
		WARNING: frame='REST' not yet implemented
	doppler -- doppler mode
	        options: (str) 'RADIO','OPTICAL','Z','BETA','GAMMA'
	        default: currently set doppler in scantable
	scanlist -- list of scan numbers to process
	        default: [] (use all scans)
		example: [21,22,23,24]
	field -- selection string for selecting scans by name
	        default: '' (no name selection)
		example: 'FLS3a*'
		         this selection is in addition to scanlist
			 and iflist
	iflist -- list of IF id numbers to select
	        default: [] (use all IFs)
		example: [15]
	scanaverage -- average integrations within scans
	        options: (bool) True,False
		default: False
		example: if True, this happens in read-in
	timeaverage -- average times for multiple scans
	        options: (bool) True,False
		default: True
	polaverage -- average polarizations
	        options: (bool) True,False
		default: True
	kernel -- type of spectral smoothing
	        options: 'hanning','gaussian','boxcar','none'
	        default: 'none'
	kwidth -- width of spectral smoothing kernel
	        options: (int) in channels 
	        default: 5
		example: 5 or 10 seem to be popular for boxcar
		         ignored for hanning (fixed at 5 chans)
                         (0 will turn off gaussian or boxcar)
	stack -- code for stacking on single plot
	        options: 'p','b','i','t','s' or
		         'pol', 'beam', 'if', 'time', 'scan'
		default: 'p'
		example: maximum of 25 stacked spectra
	                 stack by pol, beam, if, time, scan
	panel -- code for splitting into multiple panels
	        options: 'p','b','i','t','s' or
		         'pol', 'beam', 'if', 'time', 'scan'
		default: 'i'
		example: maximum of 25 panels
	                 panel by pol, beam, if, time, scan
	flrange -- range for flux axis of plot
	        options: (list) [min,max]
		default: [] (full range)
		example: flrange=[-0.1,2.0] if 'K'
		         assumes current fluxunit
	sprange -- range for spectral axis of plot
	        options: (list) [min,max]
		default: [] (full range)
		example: sprange=[42.1,42.5] if 'GHz'
		         assumes current specunit
	linecat -- control for line catalog plotting
	        options: (str) 'all','none' or by molecule
		default: 'none' (no lines plotted)
		example: linecat='SiO' for SiO lines
		         linecat='*OH' for alcohols
	                 uses sprange to limit catalog
		WARNING: specunit must be in frequency (*Hz)
			 to plot from the line catalog!
			 and must be 'GHz' or 'MHz' to use 
			 sprange to limit catalog
	linedop -- doppler offset for line catalog plotting
	        options: (float) doppler velocity (km/s)
		default: 0.0
		example: linedop=-30.0
	plotfile -- file name for hardcopy output
	        options: (str) filename.eps,.ps,.png
		default: '' (no hardcopy)
		example: 'specplot.eps','specplot.png'
	                 Note this autodetects the format from
	                 the suffix (.eps,.ps,.png).
	                 

        DESCRIPTION:

	Task sdplot displays single-dish spectra.  
	It assumes that the spectra have been calibrated in sdcal.
	It does allow selection of scans, IFs,polarizations, and
	some time and channel averaging/smoothing options also,
	but does not write out this data.

	Some plot options, like annotation and changing titles,
	legends, colors, fonts, and the like are not supported
	in this task.  You should use sd.plotter from the ASAP
	toolkit directly for this.

	This task uses the JPL line catalog as supplied by ASAP.
	If you wish to use a different catalog, or have it plot
	the line IDs from top or bottom (rather than alternating),
	then you will need to explore the sd toolkit also.	
	
	Note that multiple scans and IFs can in principle be handled
	through stacking and panelling, but this is fairly rudimentary
	at present and you have little control of what happens in
	individual panels.  We recommend that you use scanlist, field,
	and iflist to give a single selection for each run.

	Currently, setting specunit='GHz' fixes the x-axis span of
	each IF panel to be the same (an example of the limitations
	of ASAP plotting at present).

	ASAP recognizes the data of the "AT" telescopes, but currently
	does not know about the GBT or any other telescope. This task
	does know about GBT.  Therefore, if you wish to change the
	fluxunit (see below), then you need to tell it what to do.  If
	you set telescope='AT' it will use its internal defaults.  If
	you set telescope 'GBT', it will use an approximate aperture
	efficiency conversion.  See sdcal description for information
	on fluxunit conversion.
	
	WARNING: be careful plotting otf data with lots of fields!

	"""
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['__last_taskname']='sdplot'
        tb=myf['tb']
        im = myf['im']
	qa = myf['qa']
        me = myf['me']
	sm = myf['sm']
        ia = myf['ia']
        cl = myf['cl']
        pl = myf['pl']
	rg = myf['rg']
	default=myf['default']
	plotxy=myf['plotxy']
        ###
        ###
        #Handle globals or user over-ride of arguments
        #
        function_signature_defaults=dict(zip(sdplot.func_code.co_varnames,sdplot.func_defaults))
        for item in function_signature_defaults.iteritems():
		key,val = item
                keyVal = eval(key)
                if (keyVal == None):
                        #user hasn't set it - use global/default
                        pass
                else:
                        #user has set it - use over-ride
			#        sys._getframe(1).f_globals[key]=keyVal
			myf[key]=keyVal
			
	###
	#
        #Add type/menu/range error checking here
#	if type(mask)==str: mask=[mask]
        #
	# LIST PARAMETER NAMES HERE
        arg_names=['sdfile','telescope','fluxunit','specunit','frame','doppler',\
		   'scanlist','field','iflist',\
		   'scanaverage','timeaverage','polaverage',\
		   'kernel','kwidth','stack','panel','flprange','sprange',\
		   'linecat','linedop','plotfile']
	#
	# Assign parameter values
	arg_values = []
	for ptst in arg_names:
		stst = myf[ptst]
		if ( type(stst) == str ):
			# Strip off leading and trailing whitespace
			stst = stst.strip()
			myf[ptst] = stst
			sex = ptst+"='"+stst+"'"
		else:
			sex = ptst+"="+str(stst)
		exec(sex)
		arg_values = arg_values + [stst]
	#
	# LIST ALLOWED PARAMETER TYPES HERE
        arg_types=[str,(str,list),str,str,str,str,(list,int),str,(int,list),\
		   bool,bool,bool,str,(int,float),str,str,list,list,str,\
		   (float,int),str]
	#parameter_printvalues(arg_names,arg_values,arg_types)
	#
	# Check allowed types and values
        try:
                parameter_checktype(arg_names,arg_values,arg_types)
		parameter_checkmenu('fluxunit',fluxunit,['K','k','Jy','JY','jy',''])
		parameter_checkmenu('specunit',specunit,['channel','km/s','GHz','MHz','kHz','Hz',''])
		parameter_checkmenu('frame',frame,['REST','TOPO','LSRD','LSRK','BARY','GEO','GALACTO','LGROUP','CMB',''])
		parameter_checkmenu('doppler',doppler,['RADIO','OPTICAL','Z','BETA','GAMMA','radio','optical','z',''])
		parameter_checkmenu('kernel',kernel,['hanning','gaussian','boxcar','none'])
		parameter_checkmenu('stack',stack,['p','i','t','b','s', 'pol', 'beam', 'if', 'time', 'scan'])
		parameter_checkmenu('panel',panel,['p','i','t','b','s', 'pol', 'beam', 'if', 'time', 'scan'])
        except TypeError, e:
                print "sdplot -- TypeError: ", e
		return
        except ValueError, e:
                print "sdplot -- OptionError: ", e
                return
        ###

        ###
        ### Now the actual task code
        ###

	try:
#	    pdb.set_trace()
            #load the data with or without averaging
            s=sd.scantable(sdfile,scanaverage)
	    
	    # determine current fluxunit
	    fluxunit_now = s.get_fluxunit()
	    print "Current fluxunit = "+fluxunit_now
	    
	    # set flux unit string (be more permissive than ASAP)
	    if ( fluxunit == 'k' ):
		    fluxunit = 'K'
	    elif ( fluxunit == 'JY' or fluxunit == 'jy' ):
		    fluxunit = 'Jy'
	    
	    # fix the fluxunit if necessary
	    # No way to query scantable to find what telescope is
	    # so rely on user input
	    if ( telescope == 'FIX' or telescope == 'fix' ):
			    if ( fluxunit != '' ):
				    if ( fluxunit == fluxunit_now ):
					    print "No need to change default fluxunits"
				    else:
					    s.set_fluxunit(fluxunit)
					    print "Reset default fluxunit to "+fluxunit
					    fluxunit_now = s.get_fluxunit()
			    else:
				    print "Warning - no fluxunit for set_fluxunit"
		    
	    elif ( fluxunit=='' or fluxunit==fluxunit_now ):
		    print "No need to convert fluxunits"
	    elif ( telescope == 'GBT' or telescope == 'gbt' ):
		    print "Convert fluxunit to "+fluxunit
		    # THIS IS THE CHEESY PART
		    # Calculate ap.eff eta at rest freq
		    # Use Ruze law
		    #   eta=eta_0*exp(-(4pi*eps/lambda)**2)
		    # with
		    print "Using GBT parameters"
		    eps = 0.390  # mm
		    eta_0 = 0.71 # at infinite wavelength
		    # Ideally would use a freq in center of
		    # band, but rest freq is what I have
		    rf = s.get_restfreqs()[0]*1.0e-9 # GHz
		    eta = eta_0*pl.exp(-0.001757*(eps*rf)**2)
		    print "Calculated ap.eff. eta = %5.3f " % (eta)
		    print "At rest frequency %5.3f GHz" % (rf)
		    D = 104.9 # 100m x 110m
		    print "Assume phys.diam D = %5.1f m" % (D)
		    s.convert_flux(eta=eta,d=D)
		    
		    print "Successfully converted fluxunit to "+fluxunit
	    elif ( telescope == 'AT' or telescope == 'at' ):
		    s.convert_flux()
	    elif ( type(telescope) == list ):
		    # User input telescope params
		    if ( len(telescope) > 1 ):
			    D = telescope[0]
			    eta = telescope[1]
			    print "Use phys.diam D = %5.1f m" % (D)
			    print "Use ap.eff. eta = %5.3f " % (eta)
			    s.convert_flux(eta=eta,d=D)
		    elif ( len(telescope) > 0 ):
			    jypk = telescope[0]
			    print "Use gain = %6.4f Jy/K " % (jypk)
			    s.convert_flux(jyperk=jypk)
		    else:
			    print "Empty telescope list"
	    else:
		    # Unknown telescope type
		    print "Unknown telescope - cannot convert"

	    # set spectral axis unit
	    if ( specunit != '' ):
		    print "Changing spectral axis to "+specunit
		    s.set_unit(specunit)

	    # reset frame and doppler if needed
	    if ( frame != '' ):
		    print "Changing frequency frame to "+frame
		    s.set_freqframe(frame)
	    else:
		    print 'Using current frequency frame'

	    if ( doppler != '' ):
		    if ( doppler == 'radio' ):
			    ddoppler = 'RADIO'
		    elif ( doppler == 'optical' ):
			    ddoppler = 'OPTICAL'
		    elif ( doppler == 'z' ):
			    ddoppler = 'Z'
		    else:
			    ddoppler = doppler
		    
		    s.set_doppler(ddoppler)
	    else:
		    print 'Using current doppler convention'

	    # Prepare a selection
	    sel=sd.selector()

	    # Scan selection
	    if ( type(scanlist) == list ):
		    # is a list
		    scans = scanlist
	    else:
		    # is a single int, make into list
		    scans = [ scanlist ]
	    if ( len(scans) > 0 ):
		    sel.set_scans(scans)
	    
	    # Select source names
	    if ( field != '' ):
		    sel.set_name(field)
		    # NOTE: currently can only select one
		    # set of names this way, will probably
		    # need to do a set_query eventually

	    # IF selection
	    if ( type(iflist) == list ):
		    # is a list
		    ifs = iflist
	    else:
		    # is a single int, make into list
		    ifs = [ iflist ]
	    if ( len(ifs) > 0 ):
		    # Do any IF selection
		    sel.set_ifs(ifs)

	    # Apply the selection (if any)
	    s.set_selection(sel)

	    # Average in time if desired
	    if ( timeaverage ):
		    stave=sd.average_time(s,weight='tintsys')
		    del s
		    # Now average over polarizations;Tsys-weighted (1/Tsys**2) average
		    if ( polaverage ):
			    np = stave.npol()
			    if ( np > 1 ):
				    spave=stave.average_pol(weight='tsys')
			    else:
				    # only single polarization
				    print "Single polarization data - no need to average"
				    spave=stave.copy()
		    else:
			    spave=stave.copy()
		    del stave
	    else:
		    if ( polaverage ):
			    np = s.npol()
			    if ( np > 1 ):
				    spave=s.average_pol(weight='tsys')
			    else:
				    # only single polarization
				    print "Single polarization data - no need to average"
				    spave=s.copy()
		    else:
			    spave=s.copy()
		    del s
	    
	    # Smooth the spectrum (if desired)
	    if ( kernel != 'none' and (not (kwidth<=0 and kernel!='hanning'))):
		    print "Smoothing spectrum with kernel "+kernel
		    spave.smooth(kernel,kwidth)

	    # Plot final spectrum
	    # each IF is separate panel, pols stacked
	    sd.plotter.plot(spave)
	    sd.plotter.set_mode(stacking=stack,panelling=panel)

	    # Plot red x-axis at y=0 (currently disabled)
	    # sd.plotter.set_histogram(hist=True)
	    # sd.plotter.axhline(color='r',linewidth=2)

	    # Set axis ranges (if requested)
	    if len(flrange)==1: print "flrange needs 2 limits - ignoring"
	    if len(sprange)==1: print "sprange needs 2 limits - ignoring"
	    if ( len(sprange) > 1 ):
		    if ( len(flrange) > 1 ):
			    sd.plotter.set_range(sprange[0],sprange[1],flrange[0],flrange[1])
		    else:
			    sd.plotter.set_range(sprange[0],sprange[1])
	    elif ( len(flrange) > 1 ):
		    sd.plotter.set_range(ystart=flrange[0],yend=flrange[1])
	    else:
		    # Set default range explicitly (in case range was ever set)
		    sd.plotter.set_range()

	    # Line catalog
	    dolinc=False
	    if ( linecat != 'none' and linecat != '' ):
		    # Use jpl catalog for now (later allow custom catalogs)
		    catpath=casapath+'/data/catalogs/lines'
		    catname=catpath+'/jpl_asap.tbl'
		    # TEMPORARY: hard-wire to my version
		    # catname='/home/sandrock/smyers/NAUG/Tasks/jpl.tbl'
		    # FOR LOCAL CATALOGS:
		    # catname='jpl.tbl'
		    try:
			    linc=sd.linecatalog(catname)
			    dolinc=True
		    except:
			    print "Could not find catalog at "+catname
			    dolinc=False
		    if ( dolinc ):
			    if ( len(sprange)>1 ):
				    if ( specunit=='GHz' or specunit=='MHz' ):
					    linc.set_frequency_limits(sprange[0],sprange[1],specunit)
				    else:
					    print "ERROR: sd.linecatalog.set_frequency_limits accepts only GHz and MHz"
					    print "continuing without sprange selection on catalog"
			    if ( linecat != 'all' and linecat != 'ALL' ):
				    # do some molecule selection
				    linc.set_name(linecat)
			    # Plot up the selected part of the line catalog
			    # use doppler offset
			    sd.plotter.plot_lines(linc,doppler=linedop)

	    # Hardcopy
	    if ( plotfile != '' ):
		    # currently no way w/o screen display first
		    sd.plotter.save(plotfile)

	    # Do some clean up
	    del spave
	    if dolinc: del linc
	    
	    # DONE
	
	except TypeError, e:
                print "sdplot -- TypeError: ", e
		return
        except ValueError, e:
                print "sdplot -- OptionError: ", e
                return
	except Exception, instance:
		print '***Error***',instance
		return

        saveinputs=myf['saveinputs']
        saveinputs('sdplot','sdplot.last')

###End of sdplot

def sdplot_defaults():
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['sdfile']=''
	myf['telescope']=''
	myf['fluxunit']=''
	myf['specunit']=''
	myf['frame']=''
	myf['doppler']=''
        myf['scanlist']=[]
        myf['field']=''
        myf['iflist']=[]
	myf['scanaverage']=False
	myf['timeaverage']=False
	myf['polaverage']=False
	myf['kernel']='none'
        myf['kwidth']=0
	myf['stack']='p'
	myf['panel']='i'
	myf['flrange']=[]
	myf['sprange']=[]
	myf['linecat']='none'
	myf['linedop']=0.0
	myf['plotfile']=''
	
def sdplot_description(key='sdplot'):
	desc={'sdplot': 'ASAP SD calibration task\n',
	'sdfile': 'name of input SD dataset',
	'telescope': 'name or param of telescope for flux conversion',
	'fluxunit': 'units for line flux (K,Jy)',
	'specunit': 'units for spectral axis (channel,km/s,GHz)',
	'frame': 'frequency reference frame, e.g. LSRK (''=current)',
	'doppler': 'doppler convention, e.g. RADIO (''=current)',
        'scanlist': 'list of scans to use (e.g. [1,2,3,4])',
        'field': 'string for selection by source name (e.g. FLS3a*)',
        'iflist': 'list of IF ids to select (e.g. [0,1])',
	'scanaverage': 'average integs within scans (True,False)',
	'timeaverage': 'average scans over time (True,False)',
	'polaverage': 'average over polarizations (True,False)',
        'kernel': 'type of spectral smoothing (e.g. hanning)',
	'kwidth': 'width of spectral kernel in channels',
	'stack': 'code for plot stacking (p,i,b,t,s)',
	'panel': 'code for plot panelling (p,i,b,t,s)',
	'flrange': 'flux axis range [min,max] ([]=all)',
	'sprange': 'spectral axis range [min,max] ([]=all)',
	'linecat': 'line catalog control (none,all,or by molecule)',
	'linedop': 'line catalog doppler offset (km/s)',
	'plotfile': 'file name for hardcopy (''=none,.eps,.ps,.png)',
	}
	return desc[key]

###End of task

def sdlist(infile=None,scanaverage=None,listfile=None):
	"""ASAP SD listing task (prototype):

	Task Version 2007-03-06 STM
	Keyword arguments:
        infile -- name of input SD dataset
	scanaverage -- average integrations within scans
	        options: (bool) True,False
		default: False
		example: if True, this happens in read-in
	                 For GBT, set False!
	listfile -- Name of output file for summary list
	        default: '' (no output file)
		example: 'mysd_summary.txt'
	        WARNING: output file will be overwritten


        DESCRIPTION:

	Task sdlist lists the scan summary of the dataset after importing
	as a scantable into ASAP.  It will optionally output this summary
	as file.

	Note that if your PAGER environment variable is set to 'less' and
	you have set the 'verbose' ASAP environment variable to True
	(the default), then the screen version of the summary will page.
	You can disable this for sdlist by setting
	sd.rcParams['verbose']=False
	before running sdlist.  Set it back afterwards if you want lots
	of information.

	"""
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['__last_taskname']='sdlist'
        tb=myf['tb']
        im = myf['im']
	qa = myf['qa']
        me = myf['me']
	sm = myf['sm']
        ia = myf['ia']
        cl = myf['cl']
        pl = myf['pl']
	rg = myf['rg']
	default=myf['default']
	plotxy=myf['plotxy']
        ###
        ###
        #Handle globals or user over-ride of arguments
        #
        function_signature_defaults=dict(zip(sdlist.func_code.co_varnames,sdlist.func_defaults))
        for item in function_signature_defaults.iteritems():
		key,val = item
                keyVal = eval(key)
                if (keyVal == None):
                        #user hasn't set it - use global/default
                        pass
                else:
                        #user has set it - use over-ride
			#        sys._getframe(1).f_globals[key]=keyVal
			myf[key]=keyVal
			
	###
	#
        #Add type/menu/range error checking here
#	if type(mask)==str: mask=[mask]
        #
	# LIST PARAMETER NAMES HERE
        arg_names=['infile','scanaverage','listfile']
	#
	# Assign parameter values
	arg_values = []
	for ptst in arg_names:
		stst = myf[ptst]
		if ( type(stst) == str ):
			# Strip off leading and trailing whitespace
			stst = stst.strip()
			myf[ptst] = stst
			sex = ptst+"='"+stst+"'"
		else:
			sex = ptst+"="+str(stst)
		exec(sex)
		arg_values = arg_values + [stst]
	#
	# LIST ALLOWED PARAMETER TYPES HERE
        arg_types=[str,bool,str]
	#parameter_printvalues(arg_names,arg_values,arg_types)
	#
	# Check allowed types and values
        try:
                parameter_checktype(arg_names,arg_values,arg_types)
        except TypeError, e:
                print "sdlist -- TypeError: ", e
		return
        except ValueError, e:
                print "sdlist -- OptionError: ", e
                return
        ###

        ###
        ### Now the actual task code
        ###

	try:
#	    pdb.set_trace()
            #load the data with or without averaging
            s=sd.scantable(infile,scanaverage)

	    if ( listfile == '' ):
		    sum = s.summary()
	    else:
		    sum = s.summary(listfile)

	    if ( sd.rcParams['verbose'] == 'False'):
		    # print the summary to the screen manually
		    print sum
	    
	    # Clean up scantable
	    del s
	    
	    # DONE
	
	except TypeError, e:
                print "sdlist -- TypeError: ", e
		return
        except ValueError, e:
                print "sdlist -- OptionError: ", e
                return
	except Exception, instance:
		print '***Error***',instance
		return

        saveinputs=myf['saveinputs']
        saveinputs('sdlist','sdlist.last')

###End of sdlist

def sdlist_defaults():
	a=inspect.stack()
	stacklevel=0
	for k in range(len(a)):
   	  if (string.find(a[k][1], 'ipython console') > 0):
       		stacklevel=k
       		break
	myf=sys._getframe(stacklevel).f_globals
	myf['infile']=''
	myf['scanaverage']=False
        myf['listfile']=''
	
def sdlist_description(key='sdlist'):
	desc={'sdlist': 'ASAP SD listing task\n',
	'infile': 'name of input SD dataset',
	'scanaverage': 'average integs within scans (True,False)',
        'listfile': 'output file name',
	}
	return desc[key]

###End of task

def sdinp(taskname=None):
	try:
		a=inspect.stack()
		stacklevel=0
		for k in range(len(a)):
			if (string.find(a[k][1], 'ipython console') > 0):
				stacklevel=k
				break
		myf=sys._getframe(stacklevel).f_globals
		if taskname==None: taskname=myf['taskname']
		myf['taskname']=taskname
		myf['taskname']=taskname
		if type(taskname)!=str:
		       	taskname=taskname.__name__
		       	myf['taskname']=taskname
		       	myf['taskname']=taskname
                try:
                        parameter_checktype(['taskname'],taskname,str)
                        parameter_checkmenu('taskname',taskname,['sdcal','sdfit','sdplot','sdstat','sdlist'])
                except TypeError, e:
                        print "inp -- TypeError: ", e
                        return
                except ValueError, e:
                        print "inp -- OptionError: ", e
                        return
			
                ###Check if task exists by checking if task_defaults is defined
                if (myf.has_key(taskname+'_defaults') == False):
                        raise TypeError, "task %s is not defined " %taskname
		if(myf.has_key('__last_taskname')==False):
			myf[taskname+'_defaults']()
			# eval(taskname+'_defaults()')
		       	myf['__last_taskname']=taskname
		###Comment  this section if when running inp you don't want to
		###reset task parameter when last task run is different
		#if(myf['__last_taskname'] != taskname):
		#	myf[taskname+'_defaults']()
		#	myf['__last_taskname']=taskname
                ##############################################################
		f=zip(myf[taskname].func_code.co_varnames,myf[taskname].func_defaults)
		print '# %-12s       --   %s '%(taskname,myf[taskname+'_description']())
                for j in range(len(f)) :
                        k=f[j][0]
                        #keyVal=eval(k)
                        #if (f[j][1]==None): 
                        #        myf[k]=eval(k)
                        print '%-15s      = \'%10s\'\t#  %s '%(str(k),myf[k],myf[taskname+'_description'](k))
        except TypeError, e:
                print "sdinp --error: ", e
        except Exception, e:
                print "---",e


def sdsave(taskname=None, outfile=''):
	""" Save current input values to file on disk for a specified task: 

	taskname -- Name of task
		default: <unset>; example: taskname='bandpass'
		<Options: type tasklist() for the complete list>
	outfile -- Output file for the task inputs
		default: taskname.saved; example: outfile=taskname.orion

	"""
        try:
		a=inspect.stack()
		stacklevel=0
		for k in range(len(a)):
			if (string.find(a[k][1], 'ipython console') > 0):
				stacklevel=k
				break
		myf=sys._getframe(stacklevel).f_globals
		if taskname==None: taskname=myf['taskname']
		myf['taskname']=taskname
                if type(taskname)!=str:
                        taskname=taskname.__name__
                        myf['taskname']=taskname
                parameter_checktype(['taskname','outfile'],[taskname,outfile],[str,str])
                parameter_checkmenu('taskname',taskname,['sdcal','sdfit','sdplot','sdstat','sdlist'])

                ###Check if task exists by checking if task_defaults is defined
                if (myf.has_key(taskname+'_defaults') == False):
                        raise TypeError, "task %s is not defined " %taskname
                if taskname==None: taskname=myf['taskname']
                myf['taskname']=taskname
		if outfile=='': outfile=taskname+'.saved'
		taskparameterfile=open(outfile,'w')
		print >>taskparameterfile, '%-15s    = "%s"'%('taskname', taskname)
                f=zip(myf[taskname].func_code.co_varnames,myf[taskname].func_defaults)
                scriptstring='#'+str(taskname)+'('
                for j in range(len(f)):
                        k=f[j][0]
                        if(type(myf[k])==str):
                                print >>taskparameterfile, '%-15s    =  "%s"'%(k, myf[k])
                                scriptstring=scriptstring+k+'="'+myf[k]+'",'
                        else :
                                print >>taskparameterfile, '%-15s    =  %s'%(k, myf[k])
                                scriptstring=scriptstring+k+'='+str(myf[k])+','

                scriptstring=scriptstring.rstrip(',')
                scriptstring=scriptstring+')'           
                print >>taskparameterfile,scriptstring
		taskparameterfile.close()
        except TypeError, e:
                print "sdsave --error: ", e

def sdget(taskname=None, outfile=''):
	""" Get last input values from file on disk for a specified task: 

	taskname -- Name of task
		default: <unset>; example: taskname='bandpass'
		<Options: type tasklist() for the complete list>
	outfile -- Output file for the task inputs
		default: taskname.last then taskname.saved
		example: outfile=taskname.orion

	"""
        try:
		a=inspect.stack()
		stacklevel=0
		for k in range(len(a)):
			if (string.find(a[k][1], 'ipython console') > 0):
				stacklevel=k
				break
		myf=sys._getframe(stacklevel).f_globals
		if taskname==None: taskname=myf['taskname']
		myf['taskname']=taskname
                if type(taskname)!=str:
                        taskname=taskname.__name__
                        myf['taskname']=taskname
                parameter_checktype(['taskname','outfile'],[taskname,outfile],[str,str])
                parameter_checkmenu('taskname',taskname,['sdcal','sdfit','sdplot','sdstat','sdlist'])

                ###Check if task exists by checking if task_defaults is defined
                if (myf.has_key(taskname+'_defaults') == False):
                        raise TypeError, "task %s is not defined " %taskname
                if taskname==None: taskname=myf['taskname']
                myf['taskname']=taskname
                f=zip(myf[taskname].func_code.co_varnames,myf[taskname].func_defaults)
                for j in range(len(f)):
                        k=f[j][0]
			stst = myf[k]
			if ( type(stst) == str ):
				sex = k+"='"+stst+"'"
			else:
				sex = k+"="+str(stst)
			exec(sex)
		if outfile=='':
			outfile=taskname+'.last'
			try:
				taskparameterfile=open(outfile,'r')
			except:
				outfile=taskname+'.saved'
				try:
					taskparameterfile=open(outfile,'r')
				except:
					print "ERROR - no taskname.last or .saved"
					return
		
			taskparameterfile.close()
		execfile(outfile)
		# Put the task parameters back into the global namespace
                f=zip(myf[taskname].func_code.co_varnames,myf[taskname].func_defaults)
                for j in range(len(f)):
                        k=f[j][0]
                        myf[k] = eval(k)
		print "Restored parameters from file "+outfile
        except TypeError, e:
                print "sdget --error: ", e