1 files changed, 408 insertions, 0 deletions
diff --git a/noao/rv/rvvfit.x b/noao/rv/rvvfit.x
new file mode 100644
index 00000000..d900d9b3
--- /dev/null
+++ b/noao/rv/rvvfit.x
@@ -0,0 +1,408 @@
+include <gset.h>
+include "rvpackage.h"
+include "rvflags.h"
+include "rvcont.h"
+
+# RV_VERBOSE_FIT - Write a verbose description of the fit and correlation.
+
+procedure rv_verbose_fit (rv, ofd)
+
+pointer	rv				#I RV struct pointer
+int	ofd				#I Output file descriptor
+
+pointer	sp, fname
+int	fd
+int	open()
+
+errchk	mktemp,	open
+
+begin
+	if (ofd != STDOUT) {
+	    if (RV_VERBOSE(rv) == OF_TXTONLY || 
+		RV_VERBOSE(rv) == OF_NOLOG ||
+	 	RV_VERBOSE(rv) == OF_STXTONLY) {
+	            return
+	    }
+	}
+
+	# Allocate some space
+	call smark (sp)
+	call salloc (fname, SZ_PATHNAME, TY_CHAR)
+
+	# Make a temporary file name and open it up
+	if (ofd == STDOUT) {
+	    call mktemp ("uparm$tmp", Memc[fname], SZ_PATHNAME)
+	    iferr (fd = open (Memc[fname], NEW_FILE, TEXT_FILE))
+	        call error (0, "Error opening temp file.")
+
+	    call wrt_fit (rv, fd) 		# Start printing it out
+	    call close (fd) 			# Close the file
+
+	    # Page the file
+	    call gpagefile (RV_GP(rv), Memc[fname], "")
+	    call delete (Memc[fname]) 		# Clean up
+
+	} else if (ofd != NULL) {
+	    call wrt_fit (rv, ofd)
+	    call fprintf (ofd, "\f")
+	}
+
+	call sfree (sp)
+end
+
+
+# RV_MEAN_RESID - Compute the median and sigma of the residuals of the fit.
+
+procedure rv_mean_resid (rv, mean, sigma)
+
+pointer	rv				#I RV struct pointer
+real	mean				#O Mean of residuals
+real	sigma				#O Sigma of residuals
+
+pointer	sp, resx, resy
+real	x, y
+int	npts, i
+
+real	model()
+
+begin
+	if (RV_FITDONE(rv) == NO) {
+	    if (RV_INTERACTIVE(rv) == YES)
+	        call rv_errmsg ("Error: No fit yet done to the data.")
+	    return
+	}
+	if (RV_FITFUNC(rv) == SINC) {
+	    mean = 0.0
+	    sigma = 0.0
+	    return
+	}
+	npts = RV_IEND(rv) - RV_ISTART(rv) + 1
+
+	call smark (sp)
+	call salloc (resx, npts, TY_REAL)
+	call salloc (resy, npts, TY_REAL)
+
+	# Compute the residuals or ratio of the fit
+	x =  WRKPIXX(rv,RV_ISTART(rv))
+	do i = 1, npts {
+	    Memr[resx+i-1] =  x
+            if (IS_DBLSTAR(rv) == NO) {
+                switch (RV_FITFUNC(rv)) {
+                case GAUSSIAN:
+                    call cgauss1d (x, 1, COEFF(rv,1), 4, y)
+                case LORENTZIAN:
+                    call lorentz (x, 1, COEFF(rv,1), 4, y)
+                case PARABOLA:
+                    call polyfit (x, 1, COEFF(rv,1), 3, y)
+                 }
+            } else {
+                y = model (x, DBL_COEFFS(rv,1), 3*DBL_NSHIFTS(rv)+2)
+                y = DBL_SCALE(rv) * y +
+                    (DBL_Y1(rv)+DBL_SLOPE(rv)*(x-DBL_X1(rv)))
+            }
+
+	    Memr[resy+i-1] = WRKPIXY(rv,i+RV_ISTART(rv)-1) - y
+	    x = x + 1.
+	}
+	call aavgr (Memr[resy], npts, mean, sigma)
+		
+	call sfree (sp)
+end
+
+
+# WRT_FIT - Write a verbose description of the fit and correlation
+
+procedure wrt_fit (rv, fd)
+
+pointer	rv				#I RV struct pointer
+int	fd				#I Tmp file descriptor
+
+pointer	sp, ffunc, orange, rrange, system_id, title
+bool	itob()
+
+include "fitcom.com"
+
+begin
+	# Allocate some space
+	call smark (sp)
+	call salloc (ffunc, SZ_FNAME, TY_CHAR)
+	call salloc (orange, SZ_LINE, TY_CHAR)
+	call salloc (rrange, SZ_LINE, TY_CHAR)
+	call salloc (system_id, SZ_LINE, TY_CHAR)
+	call salloc (title, 2*SZ_LINE, TY_CHAR)
+
+	# Get those string valued parameters
+	call rv_make_range_string (RV_OSAMPLE(rv), Memc[orange])
+	call rv_make_range_string (RV_RSAMPLE(rv), Memc[rrange])
+	if (IS_DBLSTAR(rv) == NO)
+	    call nam_fitfunc (rv, Memc[ffunc])
+	else
+	    call strcpy ("deblend", Memc[ffunc], SZ_FNAME)
+	call sysid (Memc[system_id], SZ_LINE)
+
+	call sprintf (Memc[title], 2*SZ_LINE, "\n%14t%s\n\t   %s\n\n")
+	        call pargstr (
+		    "Description of Fit to CCF Peak and Cross-Correlation")
+		call pargstr (Memc[system_id])
+	call fprintf (fd, Memc[title])
+
+	# Write out the image stuff
+	#call fprintf (fd, "Obj  = `%.24s[%.4d]'%40tstar = `%.24s'\n")
+	call fprintf (fd, "Obj  = `%24s[%.4d]'%40tstar = `%.24s'\n")
+	    call pargstr (IMAGE(rv))
+	    call pargi (RV_APNUM(rv))
+	    call pargstr (OBJNAME(rv))
+	#call fprintf (fd, "Temp = `%.24s[%.4d]'%40tstar = `%.24s'\n")
+	call fprintf (fd, "Temp = `%24s[%.4d]'%40tstar = `%.24s'\n")
+	    call pargstr (RIMAGE(rv))
+	    call pargi (RV_APNUM(rv))
+	    call pargstr (TEMPNAME(rv))
+	call fprintf (fd, "Deltav = %.3f Km/sec%40tTempvel = %.3f Km/sec\n\n")
+	    call pargr (RV_DELTAV(rv))
+	    call pargr (TEMPVEL(rv,RV_TEMPNUM(rv)))
+
+	# Now do the fitting parameters
+	call fprintf (fd, "Fit Parameters:\n")
+	call fprintf (fd, "%10tFunction = `%s'%46tWidth = %g\n")
+	    call pargstr (Memc[ffunc])
+	    call pargr (RV_FITWIDTH(rv))
+	call fprintf (fd, "%12tHeight = %g%43tMinwidth = %g\n")
+	    call pargr (RV_FITHGHT(rv))
+	    call pargr (RV_MINWIDTH(rv))
+	call fprintf (fd, "%14tPeak = %g%43tMaxwidth = %g\n")
+	    call pargb (itob(RV_PEAK(rv)))
+	    call pargr (RV_MAXWIDTH(rv))
+	call fprintf (fd, "%11tWeights = %b%41tBackground = %g\n")
+	    call pargr (RV_WEIGHTS(rv))
+	    call pargr (RV_BACKGROUND(rv))
+	call fprintf (fd, "%9tWincenter = %g%45tWindow = %d\n")
+	    call pargr (RV_WINCENPAR(rv))
+	    call pargr (RV_WINPAR(rv))
+
+	# Write out some more fitting information
+	call fprintf (fd, "\n\tNumber of points fit = %d\n")
+	    call pargi (nfit)
+	if (RV_FITFUNC(rv) != SINC && RV_FITFUNC(rv) != CENTER1D) {
+	    call fprintf (fd, "\tNumber of iterations = %d\n")
+	        call pargi (niter)
+	    call fprintf (fd, "\tNumber of coeffs fit = 1 - %d\n")
+	        call pargi (nfitpars)
+	    call fprintf (fd, "\tChi Squared of fit = %.4g\n")
+	        call pargr (chisqr)
+	    call fprintf (fd, "\tFit Coefficients:\n")
+	    call wrt_coeffs (rv, fd)
+	}
+
+	call rv_mean_resid (rv, mresid, sresid)
+	call fprintf (fd, "\n\tMean Residual = %f\n\tSigma of Residuals = %f\n")
+	    call pargr (mresid)
+	    call pargr (sresid)
+	call fprintf (fd, "\tMaximum of cross-correlation is in bin = %d.\n")
+	    call pargi (binshift)
+	call fprintf (fd, "\tVariance of cross-correlation = %g\n")
+	    call pargr (ccfvar)
+	call fprintf (fd, "\tHJD of observation = %.5f %50tMJD = %.5f\n")
+	    if (RV_DCFLAG(rv) == -1) {
+	        call pargd (INDEFD)
+	        call pargd (INDEFD)
+	    } else {
+	        call pargd (RV_HJD(rv))
+	        call pargd (RV_MJD_OBS(rv))
+	    }
+	call fprintf (fd, "\tObject sample used in correlation = `%s'\n")
+	    call pargstr (Memc[orange])
+	call fprintf (fd, "\tTemplate sample used in correlation = `%s'\n")
+	    call pargstr (Memc[rrange])
+	call fprintf (fd, "\tTonry&Davis R value = %g\n\n")
+	    call pargr (RV_R(rv))
+
+	# Now print out some velocity information
+	call fprintf (fd, "Velocity Results:\n")
+	call wrt_velocity (rv, fd)
+
+	# Lastly print out any error comments
+	call fprintf (fd, "\nComments:\n")
+	if (RV_ERRCOMMENTS(rv) != NULL) {
+	    call fprintf (fd, "%s\n")
+	        call pargstr (ERRCOMMENTS(rv))
+	}
+	call fprintf (fd, "\n")
+
+	# Clean up
+	call flush (fd)
+	call sfree (sp)
+end
+
+
+# WRT_VELOCITY - Write out the velocity information.
+
+procedure wrt_velocity (rv, fd)
+
+pointer	rv					#I RV struct pointer
+pointer	fd					#I Output file descriptor
+
+int	i
+double	rv_shift2vel()
+
+begin
+	if (IS_DBLSTAR(rv) == NO) {
+	    call fprintf (fd, "\tShift of peak = %.4f pixels\n")
+	        call pargr (RV_SHIFT(rv))
+	    call fprintf (fd, "\tCorrelation height = %.3f\n")
+	        call pargr (RV_HEIGHT(rv))
+	    call fprintf (fd, "\tFWHM of peak = %g Km/sec\t(=%g pixels)\n\n")
+	        if (RV_DCFLAG(rv) == -1 || IS_INDEF(RV_FWHM(rv)))
+	            call pargr (INDEF)
+	        else
+	            call pargr (RV_FWHM(rv)*RV_DELTAV(rv))
+	        call pargr (RV_FWHM(rv))
+	    call fprintf (fd, "\tVelocity computed from shift = %.4f Km/sec\n")
+	        call pargr (RV_VREL(rv))
+	    call fprintf (fd, "\tObserved velocity = %.4f Km/sec\n")
+	        call pargd (RV_VOBS(rv))
+	    call fprintf (fd,"\tHeliocentric velocity = %.4f +/- %.3f Km/sec\n")
+	        call pargd (RV_VCOR(rv))
+		if (RV_DCFLAG(rv) != -1)
+	            call pargd (RV_ERROR(rv))
+		else
+	            call pargd (INDEFD)
+	} else {
+	    call fprintf (fd, "\tShift of peak[1] = %.4f pixels\n")
+	        call pargr (DBL_SHIFT(rv,1))
+	    do i = 2, DBL_NSHIFTS(rv) {
+		call fprintf (fd, "%22t[%d] = %.4f pixels\n")
+		    call pargi (i)
+		    call pargr (DBL_SHIFT(rv,i))
+	    }
+	    call fprintf (fd, "\tCorrelation height[1] = %.3f\n")
+	        call pargr (DBL_HEIGHT(rv,1))
+	    do i = 2, DBL_NSHIFTS(rv) {
+		call fprintf (fd, "%27t[%d] = %.3f\n")
+		    call pargi (i)
+		    call pargr (DBL_HEIGHT(rv,i))
+	    }
+	    call fprintf (fd, "\tFWHM of peak[1] = %f Km/s\n")
+	        call pargr (DBL_FWHM(rv,1))
+	    do i = 2, DBL_NSHIFTS(rv) {
+		call fprintf (fd, "%21t[%d] = %f Km/s\n")
+		    call pargi (i)
+		    call pargr (DBL_FWHM(rv,i))
+	    }
+
+	    call fprintf (fd, 
+		"\n\tVelocity computed from shift[1] = %.4f Km/s\n")
+	            call pargd (rv_shift2vel(rv,DBL_SHIFT(rv,1)))
+	    do i = 2, DBL_NSHIFTS(rv) {
+		call fprintf (fd, "%37t[%d] = %.4f Km/s\n")
+		    call pargi (i)
+	            call pargd (rv_shift2vel(rv,DBL_SHIFT(rv,i)))
+	    }
+	    call fprintf (fd, "\tObserved velocity[1] = %.4f Km/s\n")
+	        call pargr (DBL_VOBS(rv,1))
+	    do i = 2, DBL_NSHIFTS(rv) {
+		call fprintf (fd, "%26t[%d] = %.4f Km/s\n")
+		    call pargi (i)
+		    call pargr (DBL_VOBS(rv,i))
+	    }
+	    call fprintf(fd,"\tHeliocentric velocity[1] = %.4f +/- %.3f Km/s\n")
+	        call pargr (DBL_VHELIO(rv,1))
+	        call pargr (DBL_VERR(rv,1))
+	    do i = 2, DBL_NSHIFTS(rv) {
+		call fprintf (fd, "%30t[%d] = %.4f +/- %.3f Km/s\n")
+		    call pargi (i)
+		    call pargr (DBL_VHELIO(rv,i))
+		    call pargr (DBL_VERR(rv,i))
+	    }
+	}
+	call flush (fd)
+end
+
+
+# WRT_COEFFS - Write the fit coefficients and errors
+
+procedure wrt_coeffs (rv, fd)
+
+pointer	rv					#I RV struct pointer
+int	fd					#I File descriptor
+
+begin
+	if (fd == NULL)
+	    return
+
+	if (IS_DBLSTAR(rv) == YES) {
+	    call wrt_debl_coeffs (rv, fd)
+
+	} else if (RV_FITFUNC(rv) == GAUSSIAN || RV_FITFUNC(rv) == LORENTZIAN) {
+	    call fprintf (fd, "\t\tc[1] = %8.4f  +/- %6.4f%65t# Amplitude\n")
+		call pargr (COEFF(rv,1))
+		call pargr (ECOEFF(rv,1))
+	    call fprintf (fd, "\t\tc[2] = %8.4f  +/- %6.4f%65t# Center\n")
+		call pargr (COEFF(rv,2))
+		call pargr (ECOEFF(rv,2))
+	    if (RV_FITFUNC(rv) == GAUSSIAN)
+	        call fprintf (fd, "\t\tc[3] = %8.4f  +/- %6.4f%65t# Sigma^2\n")
+	    else
+	        call fprintf (fd, "\t\tc[3] = %8.4f  +/- %6.4f%65t# FWHM\n")
+		    call pargr (COEFF(rv,3))
+		    call pargr (ECOEFF(rv,3))
+	    if (IS_INDEF(RV_BACKGROUND(rv))) {
+	        call fprintf (fd, 
+		    "\t\tc[4] = %8.4f  +/- %6.4f%65t# Background\n")
+		        call pargr (COEFF(rv,4))
+		        call pargr (ECOEFF(rv,4))
+	    } else {
+	        call fprintf (fd, 
+		    "\t\tc[4] = %8.4f  (fixed)%65t# Background\n")
+		        call pargr (RV_BACKGROUND(rv))
+	    }
+
+	} else if (RV_FITFUNC(rv) == PARABOLA) {
+	    call fprintf (fd, "\t\tc[1] = %8.4f  +/- %6.4f\n")
+		call pargr (COEFF(rv,1))
+		call pargr (ECOEFF(rv,1))
+	    call fprintf (fd, "\t\tc[2] = %8.4f  +/- %6.4f\n")
+		call pargr (COEFF(rv,2))
+		call pargr (ECOEFF(rv,2))
+	    call fprintf (fd, "\t\tc[3] = %8.4f  +/- %6.4f\n")
+		call pargr (COEFF(rv,3))
+		call pargr (ECOEFF(rv,3))
+	}
+	call flush (fd)
+end
+
+
+# WRT_DEBL_COEFFS - Write the fit coefficients and errors for a deblended fit.
+
+procedure wrt_debl_coeffs (rv, fd)
+
+pointer	rv					#I RV struct pointer
+int	fd					#I File descriptor
+
+int	i
+
+begin
+	if (fd == NULL)
+	    return
+
+	call fprintf (fd, "\t\tc[1] = %8.4f %45t# First line position\n")
+	    call pargr (DBL_COEFFS(rv,1))
+	call fprintf (fd, "\t\tc[2] = %8.4f %45t# First line sigma\n")
+	    call pargr (DBL_COEFFS(rv,2))
+	do i = 1, DBL_NSHIFTS(rv) {
+	    call fprintf (fd, "\t\tc[%d] = %8.4f %45t# Line #%d Amplitude\n")
+		call pargi (3*i)
+	        call pargr (DBL_COEFFS(rv,3*i))
+		call pargi (i)
+	    call fprintf (fd, 
+		"\t\tc[%d] = %8.4f %45t# Line #%d Center (relative)\n")
+		    call pargi (3*i+1)
+	            call pargr (DBL_COEFFS(rv,3*i+1))
+		    call pargi (i)
+	    call fprintf (fd, 
+		"\t\tc[%d] = %8.4f %45t# Line #%d Sigma (relative)\n")
+		    call pargi (3*i+2)
+	            call pargr (DBL_COEFFS(rv,3*i+2))
+		    call pargi (i)
+	}
+	call flush (fd)
+end