# File rvsao/Emsao/emfitx.x
# October 3, 1995
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics

# Copyright(c) 1995 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.

#  Compute radial velocity from emission line shift

include "rvsao.h"
include "emv.h"

procedure emfitx (spectrum,wlspec,specfile,mspec,specim,pblue,pred,
                 wblue,wred)

real	spectrum[ARB]	# Unsmoothed object spectrum (wavelength-binned)
real	wlspec[ARB]	# Wavelengths for object spectrum
char	specfile[ARB]	# Name of spectrum file
int	mspec		# Spectrum number if multispec file, else 0
pointer	specim		# Header structure for spectrum
int	pblue		# Blue limit of spectrum in pixels
int	pred		# Red limit of spectrum in pixels
double	wblue		# Blue limit of spectrum in wavelength
double	wred		# Red limit of spectrum in wavelength

int	rmode		# Report format (1=normal,2=one-line)
double	tcz		# Cz for object (returned)
double	tczrms		# 1+Z rms error for object (returned)
double	tczchi		# 1+Z  chi^2 for object (returned)
int	i
double	drms, meanwidth
double	czvel		# Initial velocity guess in km/sec
int	nz
int	npix
pointer	emset		# Pointer to EMSAO parameters
pointer	clopset()

pointer	sp
pointer	smplot		# Smoothed object spectrum (wavelength-binned)
pointer	smspec		# Spectrum smoothed for searching (continuum removed)
pointer	smcont		# Continuum smoothed for searching and plotting
pointer	sfspec		# Spectrum smoothed for fitting
pointer	sfcont		# Continuum smoothed for fitting

pointer	smused		# Smoothed object spectrum within wavelength limits
pointer	smcusd		# Smoothed continuum spectrum within wavelength limits
pointer	sfused		# Object spectrum for fitting within wavelength limits
pointer	sfcusd		# Continuum smoothed for fitting within limits

bool	linefit		# True to fit line profiles (centers)
bool	velfit		# True to fit velocity for multiple lines
int	spectype	# =1 for object spectrum continuum substraction
pointer	emset		# Pset pointer for emission line information
char	wtitle[20]
bool	clpgetb()
double	clpgetd()
real	clpgetr()
int	clpgeti()

include	"emv.com"
include	"rvsao.com"
include	"contin.com"

define	restart_  10
define	fitlines_ 20
define	fitvel_   30

begin

# Get emission line fitting parameters
	iferr (pp = clopset("emsao"))
	    call error (0, "XCSAO: Error opening `emsao' pset")
	newresults = FALSE
	gvel = 0.d0
	drms = 0.d0
	drms = clpgetd (emset, "disperr")
	meanwidth = 0.d0
	spectype = 1
	do i = 1, MAXREF {
	    override[i] = 0
	    }
	call sprintf (wtitle,20,"Wavelength")

# Number of points to fit at line peak
	npfit = 2
	npfit = clpgeti (emset,"npfit")

# Number of coefficients to fit for line continuum
	nlcont=1
	nlcont = clpgeti (emset,"nlcont")

# Line search wavelength half-width in angstroms
	wspan = 10.d0
	wspan = clpgeti (emset,"wspan")

	zsig = 2.d0
	zsig = clgetd (emset, "linesig")

	tczrms = 0.d0
	tczchi = 0.d0

	debug = clpgetb (emset, "debug")

#  Spectrum smoothing for line fitting (Use with care!)
	esmooth = 0
	esmooth = clpgeti (emset,"esmooth")

#  Continuum parameters
	call cont_get_pars()
	abrej[1] = 0.
	emrej[1] = 0.

# Type of velocity for plotting emission and absorption lines
	call clgstr ("vel_plot",vel_plot,SZ_LINE)
	vplot = strdic (vel_plot,vel_plot,SZ_LINE,PL_VTYPES)

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

# Initialize wavelengths from files
	call eminitx (emset, TRUE)

	sp = NULL
	sfspec = NULL
	sfcont = NULL
	smcont = NULL
	smspec = NULL

#  Allocate unsmoothed spectrum to be fit
	call smark (sp)
	call salloc (sfspec,specpix,TY_REAL)
	if (sfspec == NULL)
	    call printf (" EMFIT:  cannot allocate sfspec\n")

#  Allocate unsmoothed continuum for fit
	call salloc (sfcont, specpix, TY_REAL)
	if (sfcont == NULL)
	    call printf (" EMFIT:  cannot allocate sfcont\n")

#  Allocate smoothed spectrum for plotting
	call salloc (smplot,specpix,TY_REAL)
	if (smplot == NULL)
	    call printf (" EMFIT:  cannot allocate smplot\n")

#  Allocate smoothed spectrum for line search
	call salloc (smspec,specpix,TY_REAL)
	if (smspec == NULL)
	    call printf (" EMFIT:  cannot allocate smspec\n")

#  Allocate smoothed continuum for line search
	call salloc (smcont, specpix, TY_REAL)
	if (smcont == NULL)
	    call printf (" EMFIT:  cannot allocate smcont\n")

#  Set pointers for wavelength region actually used
	npix = pred - pblue + 1
	smused = smspec + pblue - 1
	smcusd = smcont + pblue - 1
	sfused = sfspec + pblue - 1
	sfcusd = sfcont + pblue - 1

	nsmooth = 10
	nsmooth = clpgeti (emset,"nsmooth")
	skyspec = FALSE

# If plot enabled, show the object and spectra.     
	if (pltspec) {
	    call plotspec (npix,spectrum[pblue],specname,wlspec[pblue],wtitle,nsmooth)
	    }

#  Refit lines or velocity if switches set while in interactive cursor mode
	if (sfit) {
	    linefit = TRUE
	    velfit = TRUE
	    goto restart_
	    }
	if (lfit) {
	    linefit = TRUE
	    velfit = TRUE
	#  Fit and subtract continuum from smoothed spectrum for line search
	    call amovr (Memr[smplot],Memr[smspec],specpix)
	    call icsubcon (npix,Memr[smused],wlspec[pblue],specname,spectype,Memr[smcusd])
	    goto fitlines_
	    }
	if (vfit) {
	    velfit = TRUE
	    goto fitvel_
	    }

# Smooth spectrum vector for plotting
	call amovr (spectrum,Memr[smplot],specpix)
        call smooth (Memr[smplot],specpix,nsmooth)

	linefit = TRUE
	linefit = clpgetb (emset, "linefit")
	if (linefit)
	    velfit = TRUE
	else {
	    velfit = FALSE
	    savevel = FALSE
	    }

#  Make copy of spectrum from which to subtract continuum for line finding
restart_

#  Fit and subtract continuum for search
	call amovr (Memr[smplot],Memr[smspec],specpix)
	call icsubcon (npix,Memr[smused],wlspec[pblue],specname,spectype,Memr[smcusd])
#	if (debug) call printf ("EMFIT: Continuum subtracted from spectrum\n")

	gvel = clpgetd (emset, "czguess")
	if (gvel == INDEFD)
	    gvel = 0.d0
	if (gvel < 1.d0 && gvel > -1.d0 && gvel != 0.d0)
	     gvel = gvel * c0

#  Set initial value of Cz before fitting lines
	if (linefit) {
	    switch (vinit) {
		case VSEARCH:
		    call emguess (Memr[smspec],Memr[smcont],specpix,
				  wblue,wred,czvel,debug)
		case VGUESS:
		    czvel = gvel
		case VCORREL:
		    czvel = spxvel
		case VEMISS:
		    czvel = spevel
		case VCOMB:
		    czvel = spvel
		default:
		}
	    }
	cvel = 1.d0 + (czvel / c0)
	if (debug) {
	    call printf ("Initial velocity guess is %8.2f -> z= %6.4f\n")
		call pargd (czvel)
		call pargd (cvel)
	    }
	call flush (STDOUT)

#  Find emission lines
fitlines_
	if (linefit) {
	    call emsrch (Memr[smspec],Memr[smcont],specpix,
			 wblue,wred,cvel,debug)

#	Smooth spectrum for line fitting
	    call amovr (spectrum,Memr[sfspec],specpix)
	    call smooth (Memr[sfused],npix,esmooth)

#	Set up continuum to be subtracted before fitting lines
	    call icsubcon (npix,Memr[sfused],wlspec[pblue],specname,spectype,Memr[sfcusd])

#	Fit profiles to lines
	    call emlfit (Memr[sfspec],Memr[sfcont],sky,specpix,debug)
	    }

#  Compute velocity for all lines
fitvel_
	if (velfit) {
	    newresults = TRUE
	    if (nfound > 0)
	        call emvfit (wblue,wred,drms,meanwidth,debug,
			     tcz,tczrms,tczchi,nz)
	    else {
		tcz = 0.d0
		tczrms = 0.d0
		}
	    if (tcz != 0.d0)
		spevel = (tcz - 1.d0) * c0 + spechcv
	    else
		spevel = spechcv
	    speerr = tczrms * c0
	    nfit = nz
	    if (nfit == 1) {
		speerr = 0.d0
		speerr = clpgetr (emset, "sigline")
		if (speerr <= 0) {
		    do i = 1, nfound {
			if (wtobs[i] > 0 || override[i] > 0)
			    speerr = c0 * emparams[9,2,i]
			}
		    }
		}

#  Combine emission with cross-correlation velocity
	    spnlf = nfit
	    call vcombine (spxvel,spxerr,spxr,spevel,speerr,spnlf,spvel,sperr,debug)
	    }

#  Print table of results
	rmode = 1
	rmode = clpgeti (emset,"report_mode")
	call emrslts (specfile,mspec,specim,rmode)

	call clcpset (emset)
	call sfree (sp)

	return

end