# File rvsao/Xcsao/xcinit.x
# September 19, 2000
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics

# Copyright(c) 2000 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.
 
# XCINIT sets up cross-correlation vectors and template parameters
 
include	<imhdr.h>
include	<imio.h>
include	<fset.h>
include <smw.h>
include "rvsao.h"
include	"emv.h"
include	"xcv.h"
include	"contin.h"

procedure xcinit ()

char	vel_init[SZ_LINE]	# type of velocity for initial value:
				# |zero|guess|zguess|correlation|emission|
				# combination|"
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)
char	chopstr[SZ_LINE]	# Chop out emission lines:
				# (yes | no | tempfile | specfile)
char	zero_pad[SZ_LINE]	# Zero-pad before correlating or from file:
				# (no | yes | tempfile)
char	vel_plot[SZ_LINE]	# type of velocity for redshifting plot:
				# correlation|emission|combination|search

int	ntspec0			# Number of template multspec numbers
int	ntcor			# ntemp * ncor2
int	itemp, ncor2
double	tmpvel

bool	clgetb()
int	clgeti(), clscan()
real	clgetr()
double	clgetd()
int	strdic()
int	strlen(), ldir
int	decode_ranges()
pointer	imtopen()

include	"rvsao.com"
include "results.com"
include "contin.com"
include "template.com"
include "emv.com"

begin
	c0 = 299792.5
	dindef = INDEFD

# 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)

# 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

# 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)

	emrej1 = clgetr ("s_em_reject")
	abrej1 = clgetr ("s_abs_reject")
	emrej2 = clgetr ("t_em_reject")
	abrej2 = clgetr ("t_abs_reject")

# Renormalize spectrum if requested
	renorm = clgetb ("renormalize")

# Optional wavelength limits for cross-correlation
	minwav0 = clgetd ("st_lambda")
	if (minwav0 == 0.d0)
	    minwav0 = dindef
	maxwav0 = clgetd ("end_lambda")
	if (maxwav0 == 0.d0)
	    maxwav0 = dindef

# Templates against which to correlate spectra
	call clgstr ("templates",tempfiles,SZ_PATHNAME)
	call clgstr ("tempdir",tempdir,SZ_PATHNAME)
	ldir = strlen (tempdir)
	if (ldir > 0 && tempdir[ldir] != '/' && tempdir[ldir] != '$') {
	    tempdir[ldir+1] = '/'
	    tempdir[ldir+2] = EOS
	    }

# Multispec template numbers (use only first if multiple files)
	tband = clgeti ("tempband")
	call clgstr ("tempnum",tempnums0,SZ_LINE)

# Echelle switch (template spectrum tracks object spectrum number if true)
	echelle = clgetb ("echelle")

# Maximum number of templates to cache
	nctmax = clgeti ("tempcache")
	if (nctmax > MAXTEMPS) nctmax = MAXTEMPS

	call countemp (tempfiles,tempdir,tempnums0,echelle,debug,itemp)
	if (debug) {
	    call printf ("XCFIT: %d templates being cross-correlated\n")
		call pargi (itemp)
	    call flush (STDOUT)
	    }
	templist = imtopen (tempfiles)
	if (decode_ranges (tempnums0, tspec_range, MAX_RANGES, ntspec0) == ERR)
	    call error (1, "XCINIT: Illegal multispec list")

# Zero padding on, off, or read from template file
	call clgstr ("zeropad",zero_pad,SZ_LINE)
	tzpad = strdic (zero_pad,zero_pad,SZ_LINE, XC_ZTYPES)

# Emission line chopping flag
	call clgstr ("s_emchop",chopstr,SZ_LINE)
        sfchop = strdic (chopstr, chopstr, SZ_LINE, XC_CTYPES)
	call clgstr ("t_emchop",chopstr,SZ_LINE)
        tfchop = strdic (chopstr, chopstr, SZ_LINE, XC_CTYPES)

# Fraction of spectrum to apodize at each end
	han = clgetr ("bell_window")

# Filter type and filter constants [always cosine-bell]
	filter = COSBELL
	lo     = clgeti ("low_bin")
	toplo  = clgeti ("top_low")
	nrun   = clgeti ("nrun")
	topnrn = clgeti ("top_nrun")
	if (debug) {
	    call printf ("XCFIT: filter is %d %d %d %d\n")
		call pargi (lo)
		call pargi (toplo)
		call pargi (topnrn)
		call pargi (nrun)
	    }

# Optional tshift to be used for all spectra and templates
# tshift can also be read from the TSHIFT parameter in a template header
	tshift = clgetd ("tshift")

# If npts not a power of two, determine next largest power of two:
	npts = clgeti ("ncols")
	if (npts != 256 && npts != 512 && npts != 1024 && npts != 2048 &&
	    npts != 4096 && npts != 8192 && npts != 16384) {
	    npts = alog (real (npts)) / LN2 +  1.
	    npts = 2 ** npts
	    }

# Check for unreasonable filtering
	if (lo < 1) lo = 1
	if (lo > npts) lo = 1
	if (toplo  < 1) toplo = 1
	if (toplo > npts) toplo = 1
	if (nrun < 1) nrun = npts
	if (nrun > npts) nrun = npts
	if (topnrn < 1) topnrn = npts
	if (topnrn > npts) topnrn = npts

# Allocate correlation and velocity vectors
	ncor2 = npts * 2
	call countemp (tempfiles,tempdir,tempnums0,echelle,debug,itemp)
	if (debug) {
	    call printf ("XCFIT: %d points in spectrum, up to %d points in transforms\n")
		call pargi (npts)
		call pargi (ncor2)
	    call printf ("XCFIT: %d templates being cross-correlated\n")
		call pargi (itemp)
	    call flush (STDOUT)
	    }
	ntcor = ncor2 * itemp
	call malloc (xcor, ntcor, TY_REAL)
	call malloc (xvel, ntcor, TY_REAL)

# Optional object velocity to be used to shift templates
	czguess0 = 0.d0
	czguess0 = clgetd ("czguess")

# Number of times to shift template
	nzpass = 0
	nzpass = clgeti ("nzpass")
	nzpass = nzpass + 1

# If npts not a power of two, determine next largest power of two:
	npts = clgeti ("ncols")
	if (npts != 256 && npts != 512 && npts != 1024 && npts != 2048 &&
	    npts != 4096 && npts != 8192 && npts != 16384) {
	    npts = alog (real (npts)) / LN2 +  1.
	    npts = 2 ** npts
	    }

# Allocate rebinned spectra and wavelength vectors
	call malloc (shspec, npts, TY_REAL)
	call malloc (shtemp, npts, TY_REAL)
	call malloc (wltemp, npts, TY_REAL)

# Set report mode
	rmode = 2
	rmode = clgeti ("report_mode")

# Set flag to remove regions from spectrum
	if (clscan("fixbad") != EOF)
	    fixbad = clgetb ("fixbad")
	else
	    fixbad = FALSE

# Set fraction of correlation peak to fit
	pkfrac0 = clgetd ("pkfrac")


# Set flag to save correlation information to archive file
	arcflag = FALSE
	arcflag = clgetb ("archive")

# Set redshift velocity limits
	minvel = clgetd ("minvel")
	maxvel = clgetd ("maxvel")
	if (maxvel != dindef && minvel != dindef && maxvel < minvel) {
	    tmpvel = minvel
	    minvel = maxvel
	    maxvel = tmpvel
	    }

# Print the filtered cross-correlation function if requested
	filexcor = FALSE
	filexcor = clgetb ("xcor_file")
 
	return

end


procedure countemp (tempfiles, tempdir, tempnums0, echelle, debug, ntemp)

char	tempfiles[ARB]		# List of template spectra
char	tempdir[ARB]		# Directory for template spectra
char	tempnums0[ARB]		# List of multispec spectra to read
bool	echelle			# If true, template apertures track object
bool	debug
int	ntemp			# Number of templates (returned)

char	tempfile[SZ_PATHNAME]	# Template spectrum file name
char	temppath[SZ_PATHNAME]	# Template spectrum path name
char	tempnums[SZ_LINE]	# List of multispec spectra to read
pointer templist		# List of template files
int     ntspec,ntspec0		# Number of template multispec spectra
int     tspec_range[3,MAX_RANGES]
int	tspec			# Template spectrum to read from multispec file
int	ip, jp, lfile
char	lbracket[3]	# left bracket, parenthesis, brace
char	rbracket[3]	# right bracket, parenthesis, brace
char	cdot

int	stridx(), stridxs()
int	decode_ranges(),get_next_number()
int	imtgetim(), imaccess(), strlen()
pointer	imtopen()

define  newtemp_ 40
define  newtap_  50
define  endspec_ 60

begin
	call sprintf (lbracket,3,"[({")
	call sprintf (rbracket,3,"])}")
	cdot = '.'

	templist = imtopen (tempfiles)
	if (decode_ranges (tempnums0, tspec_range, MAX_RANGES, ntspec0) == ERR)
	    call error (1, "COUNTEMP: Illegal multispec list")
	ntemp = 0

# Initialize new template file
newtemp_
	tspec = -1

# Get next template spectrum file name from the list
	if (imtgetim (templist, tempfile, SZ_PATHNAME) == EOF)
	    go to endspec_

# Check for specified apertures in multispec template file
	ip = stridxs (lbracket,tempfile)
	if (echelle)
	    ntspec = 1
	else if (ip > 0) {
	    lfile = strlen (tempfile)
	    tempfile[ip] = EOS
	    ip = ip + 1
	    jp = 1
	    while (stridx (tempfile[ip],rbracket) == 0 && ip <= lfile) {
		tempnums[jp] = tempfile[ip]
		tempfile[ip] = EOS
		ip = ip + 1
		jp = jp + 1
		}
	    tempnums[jp] = EOS
	    if (decode_ranges (tempnums,tspec_range,MAX_RANGES,ntspec) == ERR)
		call error (1, "COUNTEMP: Illegal template multispec list")
	    }
	else {
	    call strcpy (tempnums0,tempnums,SZ_LINE)
	    if (decode_ranges (tempnums,tspec_range,MAX_RANGES,ntspec) == ERR)
		call error (1, "COUNTEMP: Illegal template multispec list")
	    }
	if (debug) {
	    call printf ("COUNTEMP: template file %d is %s [%s] = %d aps\n")
		call pargi (ntemp+1)
		call pargstr (tempfile)
		call pargstr (tempnums)
		call pargi (ntspec)
	    call flush (STDOUT)
	    }

# Check for readability of template spectrum
	if (stridx ("/",tempfile) == 1 || stridx ("$",tempfile) > 0 ||
	    stridx (cdot,tempfile) == 1)
	    call strcpy (tempfile,temppath,SZ_PATHNAME)
	else {
	    call strcpy (tempdir,temppath,SZ_PATHNAME)
	    call strcat (tempfile,temppath,SZ_PATHNAME)
	    }
	if (imaccess (temppath, READ_ONLY) == NO) {
	    call eprintf ("COUNTEMP: cannot read template file %s \n")
		call pargstr (temppath)
	    go to newtemp_
	    }
	if (debug) {
	    call printf ("COUNTEMP: template path %d is %s\n")
		call pargi (ntemp+1)
		call pargstr (temppath)
	    call flush (STDOUT)
	    }

# Get next multispec aperture number from list
newtap_
	if (echelle)
	    tspec = 1
	else if (get_next_number (tspec_range, tspec) == EOF)
	    go to newtemp_

# Add this template to counter
	ntemp = ntemp + 1
	if (debug) {
	    call printf ("COUNTEMP: %d templates\n")
		call pargi (ntemp)
	    call flush (STDOUT)
	    }

	ntspec = ntspec - 1
	if (ntspec > 0)
	    go to newtap_
	go to newtemp_

endspec_
	call imtclose (templist)
	if (debug) {
	    call printf ("COUNTEMP: counted %d templates\n")
		call pargi (ntemp)
	    call flush (STDOUT)
	    }
	return

end


# Free rebinned spectra, wavelength, and correlation vectors

procedure xcfree ()

include "template.com"

begin
	call mfree (xcor, TY_REAL)
	call mfree (xvel, TY_REAL)
	call mfree (shspec, TY_REAL)
	call mfree (shtemp, TY_REAL)
	call mfree (wltemp, TY_REAL)
	return
end

# Nov 18 1999	New subroutines
# Dec  8 1999	Read maximum number of templates to cache