1 files changed, 184 insertions, 0 deletions
diff --git a/noao/imred/dtoi/hdshift.x b/noao/imred/dtoi/hdshift.x
new file mode 100644
index 00000000..cde846ad
--- /dev/null
+++ b/noao/imred/dtoi/hdshift.x
@@ -0,0 +1,184 @@
+include	<math/curfit.h>
+
+# T_HDSHIFT -- Task hdshift in the dtoi package.  This task is provided
+# to support Kormendy's method of combining related characteristic curves.
+# A zero point shift in log exp unique to each set of spots is calculated and
+# subtracted.  A single curve is fit to the combined, shifted data in
+# a separate task (hdfit).
+
+procedure t_hdshift ()
+
+pointer	sp, de, fun, save, db, cv, ps_coeff, den, exp, cvf
+int	fd, rec, nsave, ncoeff, nvalues, i, nfile, junk
+real	a0, ref_a0
+
+pointer	ddb_map()
+int	clpopni(), ddb_locate(), ddb_geti(), cvstati(), strncmp(), clgfil()
+
+begin
+	# Allocate space on stack for string buffers
+	call smark (sp)
+	call salloc (de, SZ_FNAME, TY_CHAR)
+	call salloc (cvf, SZ_FNAME, TY_CHAR)
+	call salloc (fun, SZ_FNAME, TY_CHAR)
+
+	# Get list of the database names.  The curfit information is retrieved
+	# from the first file in the list, the list is then rewound.
+
+	fd = clpopni ("database")
+	junk = clgfil (fd, Memc[cvf], SZ_FNAME)
+	call clprew (fd)
+
+	# Get coefficients of common fit from cv_file
+	db = ddb_map (Memc[cvf], READ_ONLY)
+	rec = ddb_locate (db, "cv")
+	nsave = ddb_geti (db, rec, "save")
+	call salloc (save, nsave, TY_REAL)
+	call ddb_gar (db, rec, "save", Memr[save], nsave, nsave)
+	call ddb_gstr (db, rec, "function", Memc[fun], SZ_LINE)
+
+	call cvrestore (cv, Memr[save])
+	ncoeff = cvstati (cv, CVNCOEFF)
+	call salloc (ps_coeff, ncoeff, TY_REAL)
+
+	if (strncmp (Memc[fun], "power", 5) == 0)
+	    call cvcoeff (cv, Memr[ps_coeff], ncoeff)
+	else
+	    call cvpower (cv, Memr[ps_coeff], ncoeff)
+
+	do i = 1, ncoeff {
+	    call eprintf ("%d   %.7g\n")
+	        call pargi (i)
+	        call pargr (Memr[ps_coeff+i-1])
+	}
+
+	call ddb_unmap (db)
+
+	nfile = 0
+	while (clgfil (fd, Memc[de], SZ_FNAME) != EOF) {
+	    db = ddb_map (Memc[de], READ_ONLY)
+	    call hds_read (db, den, exp, nvalues)
+	    call hds_calc (den, exp, nvalues, Memr[ps_coeff], ncoeff, a0)
+	    nfile = nfile + 1
+	    if (nfile == 1)
+		ref_a0 = a0
+	    a0 = a0 - ref_a0
+
+	    call printf ("file %s: subtracting zero point a0 = %.7g\n")
+		call pargstr (Memc[de])
+		call pargr (a0)
+
+	    # Write new log exposure information to database
+	    db = ddb_map (Memc[de], APPEND)
+	    call hds_wdb  (db, exp, nvalues, a0)
+	    call mfree (den, TY_REAL)
+	    call mfree (exp, TY_REAL)
+	    call ddb_unmap (db)
+	}
+
+	call clpcls (fd)
+	call sfree (sp)
+end
+
+
+# HDS_READ -- Read the density and exposure values from the database file.
+# The density above fog and log exposure values are returned, as well as
+# the number of data pairs read.
+
+procedure hds_read (db, den, exp, nvalues)
+
+pointer	db			# Pointer to input database file
+pointer	den			# Pointer to density array - returned
+pointer	exp			# Pointer to exposure array - returned
+int	nvalues			# Number of data pairs read - returned
+
+real	fog
+int	nden, nexp, rec
+int	ddb_locate(), ddb_geti()
+real	ddb_getr()
+
+begin
+	# Get fog value to be subtracted from density
+	rec = ddb_locate (db, "fog")
+	fog = ddb_getr (db, rec, "density")
+
+	# Get density array
+	rec = ddb_locate (db, "density")
+	nden = ddb_geti (db, rec, "den_val")
+	call malloc (den, nden, TY_REAL)
+	call ddb_gar (db, rec, "den_val", Memr[den], nden, nden)
+	call asubkr (Memr[den], fog, Memr[den], nden)
+
+	# Get exposure array
+	rec = ddb_locate (db, "exposure")
+	nexp = ddb_geti (db, rec, "log_exp")
+	call malloc (exp, nexp, TY_REAL)
+	call ddb_gar (db, rec, "log_exp", Memr[exp], nexp, nexp)
+
+	nvalues = min (nden, nexp)
+end
+
+
+# HDS_CALC -- Calculate the individual shift, a0.
+
+procedure hds_calc (den, exp, nvalues, ps_coeff, ncoeff, a0)
+
+pointer	den
+pointer	exp
+int	nvalues
+real	ps_coeff[ARB]
+int	ncoeff
+real	a0
+
+int	i
+real	yavg, ycalc, xavg
+
+begin
+	# Calculate average density and log exposure values
+	xavg = 0.0
+	yavg = 0.0
+
+	do i = 1, nvalues {
+	    xavg = xavg + Memr[den+i-1]
+	    yavg = yavg + Memr[exp+i-1]
+	}
+
+	xavg = xavg / real (nvalues)
+	yavg = yavg / real (nvalues)
+
+	ycalc = 0.0
+	do i = 2, ncoeff
+	    ycalc = ycalc + ps_coeff[i] * (xavg ** real (i-1))
+
+	# Subtraction yields the zero point shift in question
+	a0 = yavg - ycalc
+end
+
+
+# HDS_WDB -- Write shifted log exposure values to database.
+
+procedure hds_wdb  (db, exp, nvalues, a0)
+
+pointer	db		# Pointer to database
+pointer	exp		# Pointer to array of exposure values
+int	nvalues		# Number of exposure values in sample
+real	a0		# Shift to be subtracted
+
+pointer	sp, expsub
+
+begin
+	call smark (sp)
+	call salloc (expsub, nvalues, TY_REAL)
+
+	call ddb_ptime (db)
+	call ddb_prec (db, "exposure")
+
+	call eprintf ("a0 = %g\n")
+	    call pargr (a0)
+
+	call asubkr (Memr[exp], a0, Memr[expsub], nvalues)
+	call ddb_par (db, "log_exp", Memr[expsub], nvalues)
+
+	call ddb_putr (db, "A0 shift", a0)
+	call sfree (sp)
+end