# File rvsao/Xcor/t_xcsao.x
# August 19, 1999
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After Gerard Kriss, Johns Hopkins University and others

# Copyright(c) 1999 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.
 
# XCSAO is an IRAF task for obtaining redshifts and velocity dispersions
# using cross correlation methods.  For arguments, see parameter file xcsao.par.
# Information is shared in common blocks defined in "rvsao.com".
 
include	<imhdr.h>
include	<imio.h>
include	<fset.h>
include	<smw.h>
include "rvsao.h"
include	"xcv.h"
include	"emv.h"

procedure t_xcsao ()

int	i
char	specfile[SZ_PATHNAME]	# Object spectrum file name
char	specpath[SZ_PATHNAME]	# Object spectrum path name
char	tempfiles[SZ_PATHNAME]	# List of template spectra
char	specdir[SZ_PATHNAME]	# Directory for object spectra
#char	filt_type[SZ_LINE]	# Filter for transform
				# (ramp | welch | hanning | cos-bell)
char	svel_corr[SZ_LINE]	# Type of velocity correction for spectrum
				# (none | file | heliocentric | barycentric)
char	tvel_corr[SZ_LINE]	# Type of velocity correction for template
				# (none | file | heliocentric | barycentric)
bool	savevel0		# Save velocity, error, and R in data file header
int	rmode		# Report format (1=normal,2=one-line)
int	logfiles	# List of log files
char	logfile[SZ_PATHNAME] # Log file name
char	wtitle[20]	# Title for wavelength plots of spectrum

int	mspec		# Object aperture to read from multispec file
int	mband		# Object band to read from multispec file

pointer	speclist	# List of spectrum files
char	str[SZ_LINE]
int	fd
char	vel_plot[SZ_LINE]	# type of velocity for redshifting plot
				# correlation|emission|combination|search

int	nmspec0		# Number of object multispec spectra
int	mspec_range[3,MAX_RANGES]
int	ip,jp,lfile	# Limits for multispec aperture decoding
char	lbracket[3]	# "[({"
char	rbracket[3]	# "])}"
 
bool	clgetb()
int	clpopnu(), clgeti(), clgfil(), open()
int	strdic(), stridx(), stridxs()
real	clgetr()
int	decode_ranges(),get_next_number()
char	vel_init[SZ_LINE]	# type of velocity for initial value
int	imtgetim(), imaccess(), strlen(), ldir, clscan()
pointer	imtopenp()

define	newspec_ 10
define	newap_	 20
define	endxc_	 90

include	"rvsao.com"
include "results.com"
include	"emv.com"
 
begin
	c0 = 299792.5
	qplot = FALSE
	nfound = 0
	nctemp = 0
	call sprintf (lbracket,3,"[({")
	call sprintf (rbracket,3,"])}")
	call sprintf (wtitle,20,"Wavelength")

# Get task parameters.

# Spectra to cross-correlate
	speclist = imtopenp ("spectra")

# Multispec spectrum numbers (use only first if multiple files)
	call clgstr ("specnum",specnums,SZ_LINE)
	if (decode_ranges (specnums, mspec_range, MAX_RANGES, nmspec0) == ERR){
	    call sprintf (str, SZ_LINE, "T_XCSAO: Illegal multispec list <%s>")
		call pargstr (specnums)
	    call error (1, str)
	    }
	call clgstr ("specdir",specdir,SZ_PATHNAME)
	ldir = strlen (specdir)
	if (specdir[1] != EOS && specdir[ldir] != '/') {
	    specdir[ldir+1] = '/'
	    specdir[ldir+2] = EOS
	    }
	mband = clgeti ("specband")

# Templates against which to correlate spectra
	call clgstr ("templates",tempfiles,SZ_PATHNAME)

# Optional correlation plot, where peak may be selected by cursor
	pltcor  = clgetb ("xcor_plot")

# Optional intermediate data plot switches
	pltspec = clgetb ("obj_plot")
	plttemp = clgetb ("temp_plot")
	pltcon  = clgetb ("contsub_plot")
	pltapo  = clgetb ("apodize_plot")
	pltfft  = clgetb ("fft_plot")
	pltuc   = clgetb ("uxcor_plot")
	plttft  = clgetb ("tfft_plot")

# Print processing information
	debug  = clgetb ("debug")

# Continuum fit parameter pset
	call cont_get_pars()

# Number of times to smooth (1-2-1) final data plot
	nsmooth = clgeti ("nsmooth")

# Velocity center and width of summary page cross-correlation plot
	xcr0 = clgetr ("cvel")
	xcrdif = clgetr ("dvel")

# Type of fit for correlation peak and fraction of peak to fit
	pkmode0 = clgeti ("pkmode")

# Type of heliocentric velocity correction to be used
	call clgstr ("svel_corr",svel_corr,SZ_LINE)
	svcor = strdic (svel_corr,svel_corr,SZ_LINE, HC_VTYPES)
	call clgstr ("tvel_corr",tvel_corr,SZ_LINE)
	tvcor = strdic (tvel_corr,tvel_corr,SZ_LINE, HC_VTYPES)

# Type of velocity for initial redshift
        call clgstr ("vel_init",vel_init,SZ_LINE)
        vinit = strdic (vel_init,vel_init,SZ_LINE,XC_VTYPES)

# Image header result flag
	savevel0 = FALSE
	savevel0 = clgetb ("save_vel")

# Report mode for log file
	rmode = 1
	rmode = clgeti ("report_mode")

# Initialize emission and absorption lines for labelling
	call eminit (FALSE)
 
# Open log files and write a header.

	logfiles = clpopnu ("logfiles")
	call fseti (STDOUT, F_FLUSHNL, YES)
	i = 0
	call strcpy ("rvsao.xcsao",taskname,SZ_LINE)
	while (clgfil (logfiles, logfile, SZ_PATHNAME) != EOF) {
	    fd = open (logfile, APPEND, TEXT_FILE)
	    if (fd == ERR) break
	    if (rmode == 1) {
		call loghead (taskname,str)
		call fprintf (fd, "%s\n")
		    call pargstr (str)
		}
	    i = i + 1
	    logfd[i] = fd
	    }
	nlogfd = i
	call clpcls (logfiles)

# Type of velocity for plotting emission and absorption lines
	if (clscan("vel_plot") != EOF) {
	    call clgstr ("vel_plot",vel_plot,SZ_LINE)
	    vplot = strdic (vel_plot,vel_plot,SZ_LINE,PL_VTYPES)
	    }
	else
	    vplot = VCORREL

# Print tab table headers
	if (rmode < 0)
	    call xcrshead (rmode)

# Get next object spectrum file name from the list
newspec_
	if (imtgetim (speclist, specfile, SZ_PATHNAME) == EOF)
	   go to endxc_

# Check for specified apertures in multispec spectrum file
	ip = stridxs (lbracket,specfile)
	if (ip > 0) {
	    lfile = strlen (specfile)
	    specfile[ip] = EOS
	    jp = 0
	    ip = ip + 1
	    while (stridx (specfile[ip],rbracket) == 0 && ip <= lfile) {
		jp = jp + 1
		specnums[jp] = specfile[ip]
		specfile[ip] = EOS
		ip = ip + 1
		}
	    if (jp > 0)
		specnums[jp+1] = EOS
	    else
		call strcpy ("0",specnums,SZ_LINE)
	    if (decode_ranges (specnums,mspec_range,MAX_RANGES,nmspec) == ERR){
		call sprintf (str, SZ_LINE, "T_XCSAO: Illegal multispec list <%s>")
		    call pargstr (specnums)
		call error (1, str)
		}
	    }
	else
	    nmspec = nmspec0
	if (debug) {
	    call printf ("XCSAO: next file is %s [%s] = %d aps\n")
		call pargstr (specfile)
		call pargstr (specnums)
		call pargi (nmspec)
	    }

# Check for readability of object spectrum
	call strcpy (specdir,specpath,SZ_PATHNAME)
	call strcat (specfile,specpath,SZ_PATHNAME)
	if (imaccess (specpath, READ_ONLY) == NO) {
	    call eprintf ("XCSAO: cannot read spectrum file %s \n")
		call pargstr (specpath)
	    go to newspec_
	    }

# Get next multispec number from list
	mspec = -1
newap_
	savevel = savevel0
	if (nmspec <= 0)
	    go to newspec_
	if (get_next_number (mspec_range, mspec) == EOF)
	    go to newspec_

	call xcfit (specfile, specdir, mspec, mband)

# Move on to next aperture or next image
	nmspec = nmspec - 1
	if (nmspec > 0)
	    go to newap_
	go to newspec_
 
# Close the log files
endxc_	do i = 1, nlogfd {
	    call close (logfd[i])
	    }                                              

#  Close spectrum list
	call imtclose (speclist)

end
 
# Mid 1970's	Paul Schechter -- Data windowing routines.  Some functions
#				  from the IBM SSP library.
# June	1987	Gerard Kriss -- Wrote major portions in SPP.
# June	1988	Stephen Levine -- major revisions
#  { {Sept	1989	Doug Mink--add parameters for 2-d echelle
# April	1990	add more fft filters and bcv
#		add option to resize velocity graph
# June	1990	add emission line chopping
# July	1990	rebin spectrum for each template
#		add John Tonry's fft
#		add Guillermo Torres' peak fit
# Sept	1990	add optional archive record output
# Oct	1990	add pixel limits for redshift
#		add normalization of spectra
#		add smoothing to spectrum display
# Dec	1990	add archive flag
# Jan	1991	write velocity to image header
#		add velocities or multiply 1+v/c's
# Mar	1991	plot final velocity correlation
# June	1991	make all velocities double
#		add IRAF continuum fitting
# Sept	1991	add second pass to align template
#		add plot options at end
# Nov 14 1991	Change label strings to char variables instead of literals
#		move vcombine and write combined velocities to header
# Nov 15 1991	Clean up debug formats
# Nov 18 1991	Change lsrplot calls to spplot and xcplot
# Nov 20 1991	Free str and logfile if necessary
# Dec  5 1991	Set velxc, vxerr, and vr for absorbtion line plot
# Dec 12 1991	Free object and template spectra when done
# Dec 16 1991	Set wavelength vector for region being cross-correlated

# Feb 18 1992	If multispec, write results to aperture-dependent strings
# Mar 27 1992	Get shift limits in velocity; get velocity vector from xcorfit
# Apr 20 1992	Pass mspec as argument
# Apr 22 1992	Put all cursor interaction and plotting into xcplot
# May 22 1992	Initialize mspec and tspec to -1
# Aug 12 1992	Add HISTORY line about XCSAO; drop polynomial continuum option
# Oct  9 1992	Read "echelle" parameter
# Nov 30 1992	Move spectrum smoothing to xcplot

# Feb  1 1993	Allow INDEF for wavelength limit parameters
# Feb  2 1993	Exit if a zero spectrum is found when renormalizing
# May  4 1993	Add option to save correlation to a file
# May 10 1993	Fix correlation writing
# May 18 1993	Move correlation file output to subroutine xcprint
# Jun  2 1993	Move wavelength <-> pixel conversion to subroutines
# Jun 14 1993	Use run-time allocation for all spectrum and wavelength vectors
# Jun 16 1993	Set version to 1.1; use shdr for wavelength conversions
# Jul  1 1993	Pass spectrum header structure to rebinl
# Jul  7 1993	Add spectrum header to getspec and gettemp
# Jul  8 1993	Close sh properly
# Aug  4 1993	Print error message and exit if there is no overlap
# Aug 20 1993	Deal with double wavelength vector
# Dec  2 1993	Only write history line once for multispec files
# Dec  2 1993	Pass mspec to xc_rslts
# Dec  3 1993	Update to version 1.2

# Feb  3 1994	Proceed gracefully if files are not useable
# Feb 11 1994	Add correlation file name list.
# Mar 23 1994	Add multispec aperture specification for each image
# Mar 23 1994	Update to version 1.3
# Apr  6 1994	Change smooth parameter to nsmooth for consistency
# Apr 12 1994	Return MWCS header pointer from getspec and gettemp
# Apr 13 1994	Drop unused variable imax
# Apr 15 1994	Initialize world coordinate system for spectrum
# Apr 19 1994	Pass object spectrum label as variable, not literal
# Apr 20 1994	Drop 2nd argument from IMTOPEN call
# Apr 26 1994	Use spectrum header wavelength vector directly
# Apr 26 1994	Use arrays rather than pointers when possible
# May  9 1994	Add number of times to smooth spectrum argument to PLOTSPEC call
# May 23 1994	Fix error message when not able to write to file
# May 23 1994	Keep smoothed spectrum local to XCPLOT
# Jun 15 1994	Move filtered-template flag from parameter list to image file
# Jun 15 1994	Ignore template directory for template names with / in them
# Jun 23 1994	Keep MWCS pointer in SHDR structure
# Jun 24 1994	Set spectrum velocity in common
# Jun 24 1994	Update to version 1.4
# Jul 29 1994	Read version from fquot.h
# Aug  3 1994	Change emission line chopping flags
# Aug  3 1994	Change common and header from fquot to rvsao
# Aug  4 1994	Keep nmspec in labelled common
# Aug  8 1994	Add specdir as parameter
# Aug 10 1994	Write log heading in one line instead of two
# Aug 24 1994	Fix bug so SPECDIR can be null
# Dec 19 1994	Move renorm and filter parameters to XCFIT

# Jan 31 1995	Change lengths of file and directory names to sz_pathname
# Mar  9 1995	Move zero-padding to XCFIT
# Apr  5 1995	Change solar system velocity correction flags
# May 15 1995	Change all sz_fname to sz_line, which is 100 chars longer
# May 17 1995	Fix bug in bracket -> mspec code
# Aug 11 1995	Move PKFRAC to XCFIT
# Oct  4 1995	Default SPECNUMS to 0 if null string is input

# Aug  7 1996	Use smw.h

# Jan 24 1997	Drop phase plot; add transformed transform plot
# Apr  9 1997	Deal with null aperture list in DECODE_RANGES, not here
# Apr 14 1997	Add option to plot template spectra
# Aug 27 1997	Add parameter for spectrum multspec band

# Apr  6 1999	Initialize number of cached templates to zero
# Jul 16 1999	Initialize qplot to FALSE
# Jul 28 1999	Print tab table headers, if appropriate
# Aug 19 1999	Set vplot from vel_plot parameter