diff options
Diffstat (limited to 'noao/onedspec/t_fitprofs.x')
| -rw-r--r-- | noao/onedspec/t_fitprofs.x | 1151 |
1 files changed, 1151 insertions, 0 deletions
diff --git a/noao/onedspec/t_fitprofs.x b/noao/onedspec/t_fitprofs.x new file mode 100644 index 00000000..9aa389bc --- /dev/null +++ b/noao/onedspec/t_fitprofs.x @@ -0,0 +1,1151 @@ +include <error.h> +include <imhdr.h> +include <smw.h> +include <gset.h> +include <ctotok.h> + + +# Profile types. +define PTYPES "|gaussian|lorentzian|voigt|" +define GAUSS 1 # Gaussian profile +define LORENTZ 2 # Lorentzian profile +define VOIGT 3 # Voigt profile + +# Type of constraints. +define FITTYPES "|fixed|single|all|" +define FIXED 1 # Fixed parameter +define SINGLE 2 # Fit a single value for all lines +define INDEP 3 # Fit independent values for all lines + +# Elements of fit array. +define BKG 1 # Background +define POS 2 # Position +define INT 3 # Intensity +define GAU 4 # Gaussian FWHM +define LOR 5 # Lorentzian FWHM + +# Output image options. +define OPTIONS "|difference|fit|" +define DIFF 1 +define FIT 2 + +# Monte-Carlo errors +define MC_N 50 # Monte-Carlo samples (overridden by users) +define MC_P 10 # Percent done interval (percent) +define MC_SIG 68 # Sigma sample point (percent) + +define NSUB 3 # Number of pixel subsamples + + +# T_FITPROFS -- Fit image profiles. + +procedure t_fitprofs() + +int inlist # List of input spectra +pointer aps # Aperture list +pointer bands # Band list + +int ptype # Profile type +pointer pg, xg, yg, sg, lg # Fitting region and initial components +real gfwhm # Default gfwhm +real lfwhm # Default lfwhm +int fit[5] # Fit flags: background, position, gfwhm, lfwhm + +int nerrsample # Number of error samples to use +real sigma0 # Constant noise +real invgain # Inverse gain + +pointer components # List of components +bool verbose # Verbose? +int log # Log file +int plot # Plot file +int outlist # List of output spectra +int option # Output image option +bool clobber # Clobber existing images? +bool merge # Merge with existing images? + +real x, y, g, l +bool complement +int i, p, ng, nalloc +pointer sp, input, output, ptr + +real clgetr() +bool clgetb() +int clgeti(), clgwrd(), clscan() +int imtopenp(), imtgetim(), imtlen() +int open(), fscan(), nscan(), strdic(), nowhite() +pointer rng_open() +errchk open + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + + # Get parameters. + inlist = imtopenp ("input") + outlist = imtopenp ("output") + if (imtlen (outlist) > 1 && imtlen (outlist) != imtlen (inlist)) + call error (1, "Input and output image lists do not make sense") + + verbose = clgetb ("verbose") + call clgstr ("logfile", Memc[output], SZ_FNAME) + if (nowhite (Memc[output], Memc[output], SZ_FNAME) == 0) + log = NULL + else + log = open (Memc[output], APPEND, TEXT_FILE) + call clgstr ("plotfile", Memc[output], SZ_FNAME) + if (nowhite (Memc[output], Memc[output], SZ_FNAME) == 0) + plot = NULL + else + plot = open (Memc[output], APPEND, BINARY_FILE) + + ptype = clgwrd ("profile", Memc[output], SZ_FNAME, PTYPES) + gfwhm = clgetr ("gfwhm") + lfwhm = clgetr ("lfwhm") + + if (clgetb ("fitbackground")) + fit[BKG] = SINGLE + else + fit[BKG] = FIXED + fit[POS] = clgwrd ("fitpositions", Memc[output], SZ_FNAME, FITTYPES) + fit[INT] = INDEP + fit[GAU] = clgwrd ("fitgfwhm", Memc[output], SZ_FNAME, FITTYPES) + fit[LOR] = clgwrd ("fitlfwhm", Memc[output], SZ_FNAME, FITTYPES) + option = clgwrd ("option", Memc[output], SZ_FNAME, OPTIONS) + clobber = clgetb ("clobber") + merge = clgetb ("merge") + nerrsample = clgeti ("nerrsample") + sigma0 = clgetr ("sigma0") + invgain = clgetr ("invgain") + if (IS_INDEF(sigma0) || IS_INDEF(invgain) || sigma0<0. || invgain<0.) { + sigma0 = INDEF + invgain = INDEF + } + + # Get the initial positions/peak/ptype/gfwhm/lfwhm. + call clgstr ("positions", Memc[input], SZ_FNAME) + if (nowhite (Memc[input], Memc[input], SZ_FNAME) == 0) { + call sfree (sp) + call error (1, "A 'positions' file must be specified") + } + i = open (Memc[input], READ_ONLY, TEXT_FILE) + ng = 0 + while (fscan (i) != EOF) { + call gargr (x) + call gargr (y) + call gargwrd (Memc[output], SZ_FNAME) + call gargr (g) + call gargr (l) + p = strdic (Memc[output], Memc[output], SZ_FNAME, PTYPES) + if (p == 0) + p = ptype + switch (nscan()) { + case 0: + next + case 1: + y = INDEF + p = ptype + g = gfwhm + l = lfwhm + case 2: + p = ptype + g = gfwhm + l = lfwhm + case 3: + g = gfwhm + l = lfwhm + case 4: + switch (p) { + case GAUSS: + l = lfwhm + case LORENTZ: + l = g + g = gfwhm + case VOIGT: + l = lfwhm + } + } + + if (ng == 0) { + nalloc = 10 + call malloc (pg, nalloc, TY_INT) + call malloc (xg, nalloc, TY_REAL) + call malloc (yg, nalloc, TY_REAL) + call malloc (sg, nalloc, TY_REAL) + call malloc (lg, nalloc, TY_REAL) + } else if (ng == nalloc) { + nalloc = nalloc + 10 + call realloc (pg, nalloc, TY_INT) + call realloc (xg, nalloc, TY_REAL) + call realloc (yg, nalloc, TY_REAL) + call realloc (sg, nalloc, TY_REAL) + call realloc (lg, nalloc, TY_REAL) + } + switch (p) { + case GAUSS: + Memi[pg+ng] = p + Memr[xg+ng] = x + Memr[yg+ng] = y + Memr[sg+ng] = g + Memr[lg+ng] = 0. + case LORENTZ: + Memi[pg+ng] = p + Memr[xg+ng] = x + Memr[yg+ng] = y + Memr[sg+ng] = 0. + Memr[lg+ng] = g + case VOIGT: + Memi[pg+ng] = p + Memr[xg+ng] = x + Memr[yg+ng] = y + Memr[sg+ng] = g + Memr[lg+ng] = l + } + ng = ng + 1 + } + call close (i) + if (ng == 0) + call error (1, "No profiles defined") + + call realloc (xg, ng+2, TY_REAL) + call realloc (yg, ng+2, TY_REAL) + call realloc (sg, ng+2, TY_REAL) + call realloc (lg, ng+2, TY_REAL) + + # Get fitting region and add to end of xg array. + i = clscan ("region") + call gargr (Memr[xg+ng]) + call gargr (Memr[xg+ng+1]) + if (i == EOF || nscan() < 1) + + # Decode range strings and set complement if needed. + complement = false + call clgstr ("lines", Memc[input], SZ_FNAME) + ptr = input + if (Memc[ptr] == '!') { + complement = true + ptr = ptr + 1 + } + iferr (aps = rng_open (Memc[ptr], INDEF, INDEF, INDEF)) + call error (1, "Bad lines/column/aperture list") + + call clgstr ("bands", Memc[input], SZ_FNAME) + ptr = input + if (Memc[ptr] == '!') { + complement = true + ptr = ptr + 1 + } + iferr (bands = rng_open (Memc[ptr], INDEF, INDEF, INDEF)) + call error (1, "Bad band list") + + # Decode components. + call clgstr ("components", Memc[input], SZ_FNAME) + iferr (components = rng_open (Memc[input], INDEF, INDEF, INDEF)) + call error (1, "Bad component list") + + while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) { + if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF) + Memc[output] = EOS + + call fp_ms (Memc[input], aps, bands, complement, Memi[pg], Memr[xg], + Memr[yg], Memr[sg], Memr[lg], ng, fit, nerrsample, + sigma0, invgain, components, verbose, log, plot, Memc[output], + option, clobber, merge) + } + + if (log != NULL) + call close (log) + if (plot != NULL) + call close (plot) + call rng_close (aps) + call rng_close (bands) + call rng_close (components) + call imtclose (inlist) + call imtclose (outlist) + call mfree (pg, TY_INT) + call mfree (xg, TY_REAL) + call mfree (yg, TY_REAL) + call mfree (sg, TY_REAL) + call mfree (lg, TY_REAL) + call sfree (sp) +end + + +# FP_MS -- Handle I/O and call fitting procedure. + +procedure fp_ms (input, aps, bands, complement, pg, xg, yg, sg, lg, ng, fit, + nerrsample, sigma0, invgain, components, verbose, log, plot, output, + option, clobber, merge) + +char input[ARB] # Input image +pointer aps # Apertures +pointer bands # Bands +bool complement # Complement aperture selection + +int pg[ng] # Profile type +real xg[ng] # Positions +real yg[ng] # Peaks +real sg[ng] # Gaussian FWHM +real lg[ng] # Lorentzian FWHM +int ng # Number of profiles +int fit[5] # Fit flags + +int nerrsample # Number of error samples +real sigma0 # Constant noise +real invgain # Inverse gain + +pointer components # Output Component list +bool verbose # Verbose output? +int log # Log file descriptor +int plot # Plot file descriptor +char output[ARB] # Output image +int option # Output image option +bool clobber # Clobber existing image? +bool merge # Merge with existing image? + +real aplow[2], aphigh[2] +double a, b, w1, wb, dw, z, p1, p2, p3 +bool select +int i, j, k, l, ap, beam, dtype, nw, ninaps, noutaps, nbands, naps, last +int mwoutdim, axis[3] +pointer ptr, in, out, tmp, mwin, mwout, sh, shout +pointer sp, str, key, temp, ltm1, ltv1, ltm2, ltv2, coeff, outaps +pointer model + +double shdr_lw() +int imaccess(), imgnfn() +bool streq(), strne(), rng_elementi(), fp_equald() +pointer smw_openim(), mw_open() +pointer immap(), imgl3r(), impl3r(), imofnlu() +errchk immap, smw_openim, mw_open, shdr_open, imunmap, imgstr, imgl3r, impl3r +errchk imdelete +data axis/1,2,3/ + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (key, SZ_LINE, TY_CHAR) + call salloc (temp, SZ_FNAME, TY_CHAR) + call salloc (ltm1, 3*3, TY_DOUBLE) + call salloc (ltv1, 3, TY_DOUBLE) + call salloc (ltm2, 3*3, TY_DOUBLE) + call salloc (ltv2, 3, TY_DOUBLE) + coeff = NULL + + # Initialize. + in = NULL; out = NULL; tmp = NULL + mwin = NULL; mwout = NULL + sh = NULL; shout = NULL + ninaps = 0; noutaps = 0; nbands = 0 + + iferr { + # Check for existing output image and abort if clobber is not set. + if (output[1] != EOS && imaccess (output, READ_ONLY) == YES) { + if (!clobber) { + call sprintf (Memc[str], SZ_LINE, + "Output spectrum %s already exists") + call pargstr (output) + call error (1, Memc[str]) + } else if (merge) { + # Merging when the input and output are the same is a nop. + if (streq (input, output)) { + call sfree (sp) + return + } + + # Open the output and check the type. + ptr = immap (output, READ_ONLY, 0); out = ptr + ptr = smw_openim (out); mwout = ptr + if (SMW_FORMAT(mwout) == SMW_ND) { + call sprintf (Memc[str], SZ_LINE, "%s - Wrong format") + call pargstr (output) + call error (1, Memc[str]) + } + + # Determine existing apertures. + noutaps = SMW_NSPEC(mwout) + nbands = SMW_NBANDS(mwout) + call salloc (outaps, noutaps, TY_INT) + do i = 1, noutaps { + call shdr_open (out, mwout, i, 1, INDEFI, SHHDR, sh) + Memi[outaps+i-1] = AP(sh) + } + } + call mktemp ("temp", Memc[temp], SZ_FNAME) + } else + call strcpy (output, Memc[temp], SZ_FNAME) + + # Open the input and determine the number of final output + # apertures in order to set the output dimensions. + + ptr = immap (input, READ_ONLY, 0); in = ptr + ptr = smw_openim (in); mwin = ptr + + naps = noutaps + + j = 1 + if (SMW_FORMAT(mwin) != SMW_ND) { + j = 0 + do i = 1, SMW_NBANDS(mwin) { + select = rng_elementi (bands, i) + if (!select) + next + j = j + 1 + } + if (j == 0) + call error (1, "No bands selected in image") + } + nbands = max (j, nbands) + + do i = 1, SMW_NSPEC(mwin) { + call shdr_open (in, mwin, i, 1, INDEFI, SHHDR, sh) + ap = AP(sh) + if (SMW_FORMAT(mwin) == SMW_ND) { + call smw_mw (mwin, i, 1, ptr, j, k) + select = rng_elementi (aps, j) && rng_elementi (bands, k) + } else + select = rng_elementi (aps, ap) + + if ((complement && select) || (!complement && !select)) + next + for (j=0; j<noutaps && Memi[outaps+j]!=ap; j=j+1) + ; + if (j == noutaps) + naps = naps + 1 + ninaps = ninaps + 1 + } + if (ninaps == 0) { + call sprintf (Memc[str], SZ_LINE, "No apertures selected in %s") + call pargstr (input) + call error (1, Memc[str]) + } + + # Set the output spectrum. For merging with an existing output + # copy to a temporary spectrum with size set appropriately. + # For a new output setup copy the input header, reset the + # physical line mapping, and clear all dispersion parameters. + + if (out != NULL) { + ptr = immap (Memc[temp], NEW_COPY, out); tmp = ptr + if (IM_PIXTYPE(tmp) != TY_DOUBLE) + IM_PIXTYPE(tmp) = TY_REAL + + IM_LEN(tmp,1) = max (SMW_LLEN(mwin,1), IM_LEN(out,1)) + IM_LEN(tmp,2) = naps + IM_LEN(tmp,3) = max (nbands, IM_LEN(out,3)) + if (nbands > 1) + IM_NDIM(tmp) = 3 + else if (naps > 1) + IM_NDIM(tmp) = 2 + else + IM_NDIM(tmp) = 1 + + do j = 1, IM_LEN(out,3) + do i = 1, IM_LEN(out,2) { + ptr = impl3r (tmp, i, j) + call aclrr (Memr[ptr], IM_LEN(tmp,1)) + call amovr (Memr[imgl3r(out,i,j)], Memr[ptr], IM_LEN(out,1)) + } + do j = 1, IM_LEN(out,3) + do i = IM_LEN(out,2)+1, IM_LEN(tmp,2) { + ptr = impl3r (tmp, i, j) + call aclrr (Memr[ptr], IM_LEN(tmp,1)) + } + do j = IM_LEN(out,3)+1, nbands + do i = 1, IM_LEN(tmp,2) { + ptr = impl3r (tmp, i, j) + call aclrr (Memr[ptr], IM_LEN(tmp,1)) + } + call imunmap (out) + out = tmp + tmp = NULL + } else if (Memc[temp] != EOS) { + ptr = immap (Memc[temp], NEW_COPY, in); out = ptr + if (IM_PIXTYPE(out) != TY_DOUBLE) + IM_PIXTYPE(out) = TY_REAL + + # Set header + IM_LEN(out,1) = SMW_LLEN(mwin,1) + IM_LEN(out,2) = naps + IM_LEN(out,3) = nbands + if (nbands > 1) + IM_NDIM(out) = 3 + else if (naps > 1) + IM_NDIM(out) = 2 + else + IM_NDIM(out) = 1 + mwoutdim = IM_NDIM(out) + + j = imofnlu (out, "DISPAXIS,APID*,BANDID*") + while (imgnfn (j, Memc[key], SZ_LINE) != EOF) + call imdelf (out, Memc[key]) + call imcfnl (j) + + i = SMW_PDIM(mwin) + j = SMW_PAXIS(mwin,1) + + ptr = mw_open (NULL, mwoutdim); mwout = ptr + call mw_newsystem (mwout, "equispec", mwoutdim) + call mw_swtype (mwout, axis, mwoutdim, "linear", "") + if (LABEL(sh) != EOS) + call mw_swattrs (mwout, 1, "label", LABEL(sh)) + if (UNITS(sh) != EOS) + call mw_swattrs (mwout, 1, "units", UNITS(sh)) + + call mw_gltermd (SMW_MW(mwin,0), Memd[ltm1], Memd[ltv1], i) + call mw_gltermd (mwout, Memd[ltm2], Memd[ltv2], mwoutdim) + Memd[ltv2] = Memd[ltv1+(j-1)] + Memd[ltm2] = Memd[ltm1+(i+1)*(j-1)] + call mw_sltermd (mwout, Memd[ltm2], Memd[ltv2], mwoutdim) + call smw_open (mwout, NULL, out) + } + + if (out != NULL) { + # Check dispersion function compatibility + # Nonlinear functions can be copied to different physical + # coordinate system though the linear dispersion can be + # modified. + + call mw_gltermd (SMW_MW(mwout,0), Memd[ltm2], Memd[ltv2], mwoutdim) + a = Memd[ltv2] + b = Memd[ltm2] + if (DC(sh) == DCFUNC) { + i = SMW_PDIM(mwin) + j = SMW_PAXIS(mwin,1) + + call mw_gltermd (SMW_MW(mwin,0), Memd[ltm1], Memd[ltv1], i) + Memd[ltv1] = Memd[ltv1+(j-1)] + Memd[ltm1] = Memd[ltm1+(i+1)*(j-1)] + if (!fp_equald (a,Memd[ltv1]) || !fp_equald (b,Memd[ltm1])) { + call error (1, + "Physical basis for nonlinear dispersion functions don't match") + } + } + } + + # Now do the actual fitting + call salloc (model, SMW_LLEN(mwin,1), TY_REAL) + last = noutaps + do i = 1, SMW_NSPEC(mwin) { + call shdr_open (in, mwin, i, 1, INDEFI, SHHDR, sh) + + # Check apertures. + ap = AP(sh) + if (SMW_FORMAT(mwin) == SMW_ND) { + call smw_mw (mwin, i, 1, ptr, j, k) + select = rng_elementi (aps, j) && rng_elementi (bands, k) + } else + select = rng_elementi (aps, ap) + + if ((complement && select) || (!complement && !select)) + next + + call fp_title (sh, Memc[str], verbose, log) + + call shdr_open (in, mwin, i, 1, INDEFI, SHDATA, sh) + if (SN(sh) < SMW_LLEN(mwin,1)) + call aclrr (Memr[model], SMW_LLEN(mwin,1)) + iferr (call fp_fit (sh, Memr[SX(sh)], Memr[SY(sh)], SN(sh), pg, + xg, yg, sg, lg, ng, fit, nerrsample, sigma0, invgain, + components, verbose, log, plot, Memc[str], Memr[model])) { + call erract (EA_WARN) + } + + if (out != NULL) { + for (j=0; j<noutaps && Memi[outaps+j]!=ap; j=j+1) + ; + + # Set output logical and physical lines + if (j < noutaps) + l = j + 1 + else { + l = last + 1 + last = l + } + + # Copy and adjust dispersion info + call smw_gwattrs (mwin, i, 1, AP(sh), beam, + dtype, w1, dw, nw, z, aplow, aphigh, coeff) + + w1 = shdr_lw (sh, 1D0) + wb = shdr_lw (sh, double (SN(sh))) + p1 = (NP1(sh) - a) / b + p2 = (NP2(sh) - a) / b + p3 = (IM_LEN(out,1) - a) / b + nw = nint (min (max (p1 ,p3), max (p1 ,p2))) + NP1(sh) - 1 + if (p1 != p2) + dw = (wb - w1) / (p2 - p1) * (1 + z) + w1 = w1 * (1 + z) - (p1 - 1) * dw + + call smw_swattrs (mwout, l, 1, ap, beam, dtype, + w1, dw, nw, z, aplow, aphigh, Memc[coeff]) + + # Copy titles + call smw_sapid (mwout, l, 1, TITLE(sh)) + if (Memc[SID(sh,1)] != EOS) + call imastr (out, "BANDID1", Memc[SID(sh,1)]) + + # Copy the data + switch (option) { + case DIFF: + call asubr (Memr[SY(sh)], Memr[model], + Memr[impl3r(out,l,1)+NP1(sh)-1], SN(sh)) + case FIT: + call amovr (Memr[model], Memr[impl3r(out,l,1)+NP1(sh)-1], + SN(sh)) + } + + # Verify copy + if (verbose) { + call shdr_open (out, mwout, l, 1, INDEFI, SHHDR, shout) + call printf ("%s%s(%d) --> %s%s(%s)\n") + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargi (AP(sh)) + call pargstr (IMNAME(shout)) + call pargstr (IMSEC(shout)) + call pargi (AP(shout)) + call flush (STDOUT) + } + } + } + + call smw_close (MW(sh)) + if (out != NULL) { + call smw_saveim (mwout, out) + if (shout != NULL) + call smw_close (MW(shout)) + call imunmap (out) + if (strne (Memc[temp], output)) { + call imdelete (output) + call imrename (Memc[temp], output) + } + } + call imunmap (in) + } then { + if (shout != NULL) + call smw_close (MW(shout)) + else if (mwout != NULL) + call smw_close (mwout) + if (sh != NULL) + call smw_close (MW(sh)) + else if (mwin != NULL) + call smw_close (mwin) + if (tmp != NULL) + call imunmap (tmp) + if (out != NULL) + call imunmap (out) + if (in != NULL) + call imunmap (in) + call erract (EA_WARN) + } + + call shdr_close (shout) + call shdr_close (sh) + call mfree (coeff, TY_CHAR) + call sfree (sp) +end + + +define SQ2PI 2.5066283 + +# FP_FIT -- Fit profile functions + +procedure fp_fit (sh, x, y, n, ptypes, pos, peaks, gfwhms, lfwhms, ng, fit, + nerrsample, sigma0, invgain, components, verbose, log, plot, title, mod) + +pointer sh # Spectrum data structure +real x[n] # Coordinates +real y[n] # Data +int n # Number of data points + +int ptypes[ARB] # Profile types +real pos[ARB] # Fitting region and initial positions +real peaks[ARB] # Peak values +real gfwhms[ARB] # Background levels and initial gfwhm +real lfwhms[ARB] # Initial lfwhm +int ng # Number of gaussian components + +int fit[5] # Fit flags + +int nerrsample # Number of error samples +real sigma0 # Constant noise +real invgain # Inverse gain + +pointer components # Component list +bool verbose # Output to STDOUT? +int log # Log file descriptor +int plot # Plot file descriptor +char title[ARB] # Plot title +real mod[n] # Model + +int i, j, k, i1, i2, nfit, nsub, mc_n, mc_p, mc_sig +long seed +real xc, x1, x2, dx, y1, dy, z1, dz, w, z, scale, sscale +real peak, flux, cont, gfwhm, lfwhm, eqw, chisq +real flux1, cont1, eqw1, wyc1, slope1, v, u +bool doerr +pointer sp, str, xd, yd, sd, xg, yg, sg, lg, pg, yd1, xg1, yg1, sg1, lg1 +pointer ym, conte, xge, yge, sge, lge, fluxe, eqwe +pointer gp, gopen() +bool rng_elementi() +real model(), gasdev(), asumr() +double shdr_lw(), shdr_wl +errchk fp_background, dofit, dorefit + +begin + # Determine fitting region. + x1 = pos[ng+1] + x2 = pos[ng+2] + i1 = nint (shdr_wl (sh, double(x1))) + i2 = nint (shdr_wl (sh, double(x2))) + i = min (n, max (i1, i2)) + i1 = max (1, min (i1, i2)) + i2 = i + nfit = i2 - i1 + 1 + if (nfit < 3) { + call aclrr (mod, n) + call error (1, "Too few data points in fitting region") + } + x1 = shdr_lw (sh, double(i1)) + x2 = shdr_lw (sh, double(i2)) + + # Allocate memory. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (xd, nfit, TY_REAL) + call salloc (yd, nfit, TY_REAL) + call salloc (sd, nfit, TY_REAL) + call salloc (xg, ng, TY_REAL) + call salloc (yg, ng, TY_REAL) + call salloc (sg, ng, TY_REAL) + call salloc (lg, ng, TY_REAL) + call salloc (pg, ng, TY_INT) + + # Subtract the continuum and scale the data. + call fp_background (sh, x, y, n, x1, x2, y1, dy) + scale = 0. + doerr = !IS_INDEF(sigma0) + do i = i1, i2 { + Memr[xd+i-i1] = x[i] + Memr[yd+i-i1] = y[i] - (y1 + dy * (x[i]-x1)) + if (y[i] <= 0.) + doerr = false + scale = max (scale, abs (Memr[yd+i-i1])) + } + if (doerr) { + do i = i1, i2 + Memr[sd+i-i1] = sqrt (sigma0 ** 2 + invgain * y[i]) + sscale = asumr (Memr[sd], nfit) / nfit + } else { + call amovkr (1., Memr[sd], nfit) + sscale = 1. + } + call adivkr (Memr[yd], scale, Memr[yd], nfit) + call adivkr (Memr[sd], sscale, Memr[sd], nfit) + y1 = y1 / scale + dy = dy / scale + + # Setup initial estimates. + do i = 1, ng { + Memr[xg+i-1] = pos[i] + Memr[sg+i-1] = gfwhms[i] + Memr[lg+i-1] = lfwhms[i] + Memi[pg+i-1] = ptypes[i] + if (IS_INDEF(peaks[i])) { + j = max (1, min (nfit, nint (shdr_wl(sh,double(pos[i])))-i1+1)) + Memr[yg+i-1] = Memr[yd+j-1] + } else + Memr[yg+i-1] = peaks[i] / scale + } + z1 = 0. + dz = 0. + dx = (x[n] - x[1]) / (n - 1) + nsub = NSUB + call dofit (fit, Memr[xd], Memr[yd], Memr[sd], + nfit, dx, nsub, z1, dz, Memr[xg], Memr[yg], Memr[sg], + Memr[lg], Memi[pg], ng, chisq) + + # Compute Monte-Carlo errors. + mc_n = nerrsample + mc_p = nint (mc_n * MC_P / 100.) + mc_sig = nint (mc_n * MC_SIG / 100.) + if (doerr && mc_sig > 9) { + call salloc (yd1, nfit, TY_REAL) + call salloc (ym, nfit, TY_REAL) + call salloc (xg1, ng, TY_REAL) + call salloc (yg1, ng, TY_REAL) + call salloc (sg1, ng, TY_REAL) + call salloc (lg1, ng, TY_REAL) + call salloc (conte, mc_n*ng, TY_REAL) + call salloc (xge, mc_n*ng, TY_REAL) + call salloc (yge, mc_n*ng, TY_REAL) + call salloc (sge, mc_n*ng, TY_REAL) + call salloc (lge, mc_n*ng, TY_REAL) + call salloc (fluxe, mc_n*ng, TY_REAL) + call salloc (eqwe, mc_n*ng, TY_REAL) + do i = 1, nfit { + w = Memr[xd+i-1] + Memr[ym+i-1] = model (w, dx, nsub, Memr[xg], Memr[yg], + Memr[sg], Memr[lg], Memi[pg], ng) + } + seed = 1 + do i = 0, mc_n-1 { + do j = 1, nfit + Memr[yd1+j-1] = Memr[ym+j-1] + + sscale / scale * Memr[sd+j-1] * gasdev (seed) + wyc1 = z1 + slope1 = dz + call amovr (Memr[xg], Memr[xg1], ng) + call amovr (Memr[yg], Memr[yg1], ng) + call amovr (Memr[sg], Memr[sg1], ng) + call amovr (Memr[lg], Memr[lg1], ng) + call dorefit (fit, Memr[xd], Memr[yd1], Memr[sd], + nfit, dx, nsub, wyc1, slope1, + Memr[xg1], Memr[yg1], Memr[sg1], + Memr[lg1], Memi[pg], ng, chisq) + + do j = 0, ng-1 { + cont = y1 + z1 + (dy + dz) * Memr[xg+j] - dy * x1 + cont1 = y1 + wyc1 + (dy + slope1) * Memr[xg+j] - dy * x1 + switch (Memi[pg+j]) { + case GAUSS: + flux = 1.064467 * Memr[yg+j] * Memr[sg+j] + flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] + case LORENTZ: + flux = 1.570795 * Memr[yg+j] * Memr[lg+j] + flux1 = 1.570795 * Memr[yg1+j] * Memr[lg1+j] + case VOIGT: + call voigt (0., 0.832555*Memr[lg+j]/Memr[sg+j], v, u) + flux = 1.064467 * Memr[yg+j] * Memr[sg+j] / v + call voigt (0., 0.832555*Memr[lg1+j]/Memr[sg1+j], v, u) + flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] / v + } + if (cont > 0. && cont1 > 0.) { + eqw = -flux / cont + eqw1 = -flux1 / cont1 + } else { + eqw = 0. + eqw1 = 0. + } + Memr[conte+j*mc_n+i] = abs (cont1 - cont) + Memr[xge+j*mc_n+i] = abs (Memr[xg1+j] - Memr[xg+j]) + Memr[yge+j*mc_n+i] = abs (Memr[yg1+j] - Memr[yg+j]) + Memr[sge+j*mc_n+i] = abs (Memr[sg1+j] - Memr[sg+j]) + Memr[lge+j*mc_n+i] = abs (Memr[lg1+j] - Memr[lg+j]) + Memr[fluxe+j*mc_n+i] = abs (flux1 - flux) + Memr[eqwe+j*mc_n+i] = abs (eqw1 - eqw) + } + } + do j = 0, ng-1 { + call asrtr (Memr[conte+j*mc_n], Memr[conte+j*mc_n], mc_n) + call asrtr (Memr[xge+j*mc_n], Memr[xge+j*mc_n], mc_n) + call asrtr (Memr[yge+j*mc_n], Memr[yge+j*mc_n], mc_n) + call asrtr (Memr[sge+j*mc_n], Memr[sge+j*mc_n], mc_n) + call asrtr (Memr[lge+j*mc_n], Memr[lge+j*mc_n], mc_n) + call asrtr (Memr[fluxe+j*mc_n], Memr[fluxe+j*mc_n], mc_n) + call asrtr (Memr[eqwe+j*mc_n], Memr[eqwe+j*mc_n], mc_n) + } + call amulkr (Memr[conte], scale, Memr[conte], mc_n*ng) + call amulkr (Memr[yge], scale, Memr[yge], mc_n*ng) + call amulkr (Memr[fluxe], scale, Memr[fluxe], mc_n*ng) + } + + call amulkr (Memr[yg], scale, Memr[yg], ng) + y1 = (y1 + z1 + dz * x1) * scale + dy = (dy + dz) * scale + + # Log computed values + call sprintf (Memc[str], SZ_LINE, + "# Nfit=%d, background=%b, positions=%s, gfwhm=%s, lfwhm=%s\n") + call pargi (ng) + call pargb (fit[BKG] == SINGLE) + if (fit[POS] == FIXED) + call pargstr ("fixed") + else if (fit[POS] == SINGLE) + call pargstr ("single") + else + call pargstr ("all") + if (fit[GAU] == FIXED) + call pargstr ("fixed") + else if (fit[GAU] == SINGLE) + call pargstr ("single") + else + call pargstr ("all") + if (fit[LOR] == FIXED) + call pargstr ("fixed") + else if (fit[LOR] == SINGLE) + call pargstr ("single") + else + call pargstr ("all") + if (log != NULL) + call fprintf (log, Memc[str]) + if (verbose) + call printf (Memc[str]) + call sprintf (Memc[str], SZ_LINE, "# %8s%10s%10s%10s%10s%10s%10s\n") + call pargstr ("center") + call pargstr ("cont") + call pargstr ("flux") + call pargstr ("eqw") + call pargstr ("core") + call pargstr ("gfwhm") + call pargstr ("lfwhm") + if (log != NULL) + call fprintf (log, Memc[str]) + if (verbose) + call printf (Memc[str]) + do i = 1, ng { + if (!rng_elementi (components, i)) + next + xc = Memr[xg+i-1] + cont = y1 + dy * (xc - x1) + peak = Memr[yg+i-1] + gfwhm = Memr[sg+i-1] + lfwhm = Memr[lg+i-1] + switch (Memi[pg+i-1]) { + case 1: + flux = 1.064467 * peak * gfwhm + case 2: + flux = 1.570795 * peak * lfwhm + case 3: + call voigt (0., 0.832555*lfwhm/gfwhm, v, u) + flux = 1.064467 * peak * gfwhm / v + } + + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + + call sprintf (Memc[str], SZ_LINE, + " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargr (xc) + call pargr (cont) + call pargr (flux) + call pargr (eqw) + call pargr (peak) + call pargr (gfwhm) + call pargr (lfwhm) + if (log != NULL) + call fprintf (log, Memc[str]) + if (verbose) + call printf (Memc[str]) + if (doerr && mc_sig > 9) { + call sprintf (Memc[str], SZ_LINE, + " (%7.7g) (%7.7g) (%7.6g) (%7.4g) (%7.6g) (%7.4g) (%7.4g)\n") + call pargr (Memr[xge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[conte+(i-1)*mc_n+mc_sig]) + call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[yge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[sge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[lge+(i-1)*mc_n+mc_sig]) + if (log != NULL) + call fprintf (log, Memc[str]) + if (verbose) + call printf (Memc[str]) + } + } + + # Compute model. + call aclrr (mod, n) + do i = 0, ng-1 { + if (!rng_elementi (components, i+1)) + next + do j = 1, n + #mod[j] = model (x[j], dx, nsub, Memr[xg+i], Memr[yg+i], + # Memr[sg+i], Memr[lg+i], Memi[pg+i], ng) + mod[j] = mod[j] + model (x[j], dx, nsub, Memr[xg+i], Memr[yg+i], + Memr[sg+i], Memr[lg+i], Memi[pg+i], 1) + } + + # Draw graphs + if (plot != NULL) { + gp = gopen ("stdvdm", NEW_FILE, plot) + call gascale (gp, y[i1], nfit, 2) + call asubr (y[i1], mod[i1], Memr[yd], nfit) + call grscale (gp, Memr[yd], nfit, 2) + do i = i1, i2 + Memr[yd+i-i1] = mod[i] + y1 + dy * (x[i] - x1) + call grscale (gp, Memr[yd], nfit, 2) + call gswind (gp, x1, x2, INDEF, INDEF) + call glabax (gp, title, "", "") + call gseti (gp, G_PLTYPE, 1) + call gpline (gp, Memr[xd], y[i1], nfit) + call gseti (gp, G_PLTYPE, 2) + call gpline (gp, Memr[xd], Memr[yd], nfit) + call gline (gp, x1, y1, x2, y1+dy*(x2-x1)) + call gseti (gp, G_PLTYPE, 3) + call asubr (y[i1], mod[i1], Memr[yd], nfit) + call gpline (gp, Memr[xd], Memr[yd], nfit) + call gseti (gp, G_PLTYPE, 4) + do i = 0, ng-1 { + if (!rng_elementi (components, i+1)) + next + k = 0 + do j = i1, i2 { + w = x[j] + z = model (w, dx, nsub, Memr[xg+i], Memr[yg+i], + Memr[sg+i], Memr[lg+i], Memi[pg+i], 1) + z = z + y1 + dy * (w - x1) + if (k == 0) { + call gamove (gp, w, z) + k = 1 + } else + call gadraw (gp, w, z) + } + } + call gclose (gp) + } + + call sfree (sp) +end + + +# FP_BACKGROUND -- Iniital background. + +procedure fp_background (sh, x, y, n, x1, x2, y1, dy) + +pointer sh #I Spectrum pointer +real x[n] #I Coordinate values +real y[n] #I Data +int n #I Number of data points +real x1, x2 #I Fit endpoints +real y1, dy #O Background + +int i, j, k, m, func +real xval[2], yval[2] +double z1, z2, z3 +pointer sp, bkg, str + +int ctotok(), ctor(), ctod(), strdic(), nscan() +real asumr(), amedr() +double shdr_wl(), shdr_lw() + +define err_ 10 + +begin + call smark (sp) + call salloc (bkg, SZ_LINE, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + xval[1] = x1 + xval[2] = x2 + + call clgstr ("background", Memc[bkg], SZ_LINE) + call sscan (Memc[bkg]) + do j = 1, 2 { + call gargwrd (Memc[bkg], SZ_LINE) + if (nscan() != j) { + i = max (1, min (n, nint (shdr_wl (sh, double(xval[j]))))) + xval[j] = shdr_lw (sh, double(i)) + yval[j] = y[i] + next + } + + k = 1 + if (ctor (Memc[bkg], k, yval[j]) == 0) { + if (ctotok (Memc[bkg], k, Memc[str], SZ_LINE) != TOK_IDENTIFIER) + goto err_ + func = strdic (Memc[str], Memc[str], SZ_LINE, "|avg|med|") + if (func == 0) + goto err_ + k = k + 1 + if (ctod (Memc[bkg], k, z1) == 0) + goto err_ + k = k + 1 + if (ctod (Memc[bkg], k, z2) == 0) + goto err_ + k = k + 1 + if (ctod (Memc[bkg], k, z3) == 0) + z3 = 1 + + z1 = shdr_wl (sh, z1) + z2 = shdr_wl (sh, z2) + i = max (1, nint(min(z1,z2))) + m = min (n, nint(max(z1,z2))) - i + 1 + if (m < 1) + goto err_ + + # This is included to eliminate an optimizer bug on solaris. + call sprintf (Memc[bkg], SZ_LINE, "%g %g %g %d %d\n") + call pargd (z1) + call pargd (z2) + call pargd (z3) + call pargi (i) + call pargi (m) + + switch (func) { + case 1: + xval[j] = z3 * asumr (x[i], m) / m + yval[j] = z3 * asumr (y[i], m) / m + case 2: + xval[j] = z3 * asumr (x[i], m) / m + yval[j] = z3 * amedr (y[i], m) + } + } + } + + if (xval[1] == xval[2]) { + dy = 0. + y1 = (yval[1] + yval[2]) / 2. + } else { + dy = (yval[2] - yval[1]) / (xval[2] - xval[1]) + y1 = yval[1] + dy * (x1 - xval[1]) + } + return + +err_ + call sfree (sp) + call error (1, "Syntax error in background specification") +end + + +include <time.h> + +# FP_TITLE -- Set title string and print. + +procedure fp_title (sh, str, verbose, log) + +pointer sh # Spectrum header structure +char str[SZ_LINE] # Title string +bool verbose # Verbose? +int log # Log file descriptor + +pointer sp, time, smw +long clktime() + +begin + # Select title format. + smw = MW(sh) + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call sprintf (str, SZ_LINE, "%s%s: %s") + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargstr (TITLE(sh)) + case SMW_ES, SMW_MS: + call sprintf (str, SZ_LINE, "%s - Ap %d: %s") + call pargstr (IMNAME(sh)) + call pargi (AP(sh)) + call pargstr (TITLE(sh)) + } + + # Set time and log header. + call smark (sp) + call salloc (time, SZ_DATE, TY_CHAR) + call cnvdate (clktime(0), Memc[time], SZ_DATE) + if (log != NULL) { + call fprintf (log, "# %s %s\n") + call pargstr (Memc[time]) + call pargstr (str) + } + if (verbose) { + call printf ("# %s %s\n") + call pargstr (Memc[time]) + call pargstr (str) + } + + call sfree (sp) +end |