diff options
Diffstat (limited to 'noao/twodspec/apextract/approfile.x')
| -rw-r--r-- | noao/twodspec/apextract/approfile.x | 765 |
1 files changed, 765 insertions, 0 deletions
diff --git a/noao/twodspec/apextract/approfile.x b/noao/twodspec/apextract/approfile.x new file mode 100644 index 00000000..eeb31a6d --- /dev/null +++ b/noao/twodspec/apextract/approfile.x @@ -0,0 +1,765 @@ +include <mach.h> +include <gset.h> +include <math/curfit.h> +include "apertures.h" + + +# AP_PROFILE -- Determine spectrum profile with pixel rejection. +# +# The profile is determined by dividing each dispersion point by an estimate +# of the spectrum and then smoothing and normalizing to unit integral. +# This routine has two algorithms (procedures) for smoothing, one for nearly +# aligned spectra and one for tilted spectra. The selection is determined +# by the calling program and signaled by whether there is a variation in +# the profile offsets. For both smoothing algorithms the same iterative +# rejection algorithm may be used to eliminate deviant points from affecting +# the profile. This rejection is selected by the "clean" parameter. +# A plot of the final profile along the dispersion may be made for the +# special plotfile "debugfits" or "debugall". +# +# Dispersion points with saturated pixels are ignored as well a when the +# total sky subtracted flux is negative. + +procedure ap_profile (im, ap, dbuf, nc, nl, c1, l1, sbuf, svar, profile, nx, ny, + xs, ys, asi) + +pointer im # IMIO pointer +pointer ap # Aperture structure +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +pointer sbuf # Sky values (NULL if none) +pointer svar # Sky variances +real profile[ny,nx] # Profile (returned) +int nx, ny # Size of profile array +int xs[ny], ys # Origin of profile array +pointer asi # Image interpolator for edge pixel weighting + +real gain # Gain +real rdnoise # Readout noise +real saturation # Maximum value for an unsaturated pixel +bool clean # Clean cosmic rays? +real lsigma, usigma # Rejection sigmas. + +int fd, ix, iy, ix1, ix2, xs1, xs2, nsum +int i, niterate, ixrej, iyrej, nrej, nreject +real p, s, chisq, tfac, rrej, predict, var0, var, vmin, resid, wt1, wt2, dat +pointer sp, str, spec, x1, x2, y, reject, xreject, data, sky, cv, gp + +int apgeti() +real apgetr(), ap_cveval(), apgimr() +bool apgetb() +errchk salloc, ap_horne, ap_marsh, apgimr, ap_asifit + +begin + # Allocate memory. Adjust pointers to be one indexed. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (spec, ny, TY_REAL) + call salloc (x1, ny, TY_REAL) + call salloc (x2, ny, TY_REAL) + call salloc (y, ny, TY_REAL) + call salloc (reject, nx*ny, TY_BOOL) + if (sbuf == NULL) { + call salloc (sky, nx, TY_REAL) + sky = sky - 1 + } + spec=spec-1; x1=x1-1; x2=x2-1; y=y-1 + + # Get task parameters. + gain = apgimr ("gain", im) + rdnoise = apgimr ("readnoise", im) ** 2 + saturation = apgetr ("saturation") + if (!IS_INDEF(saturation)) + saturation = saturation * gain + lsigma = apgetr ("lsigma") + usigma = apgetr ("usigma") + clean = apgetb ("clean") + if (clean) + niterate = apgeti ("niterate") + else + niterate = 0 + + # Initialize. + if (rdnoise == 0.) + vmin = 1. + else + vmin = rdnoise + if (sbuf == NULL) { + call aclrr (Memr[sky+1], nx) + var0 = rdnoise + } + cv = AP_CV(ap) + + # Set aperture limits and initialize rejection flags. + call alimi (xs, ny, xs1, xs2) + i = AP_AXIS(ap) + p = AP_CEN(ap,i) + AP_LOW(ap,i) + s = AP_CEN(ap,i) + AP_HIGH(ap,i) + xreject = reject + do iy = 1, ny { + dat = ap_cveval (cv, real (iy + ys - 1)) - c1 + 1 + Memr[x1+iy] = p + dat + Memr[x2+iy] = s + dat + Memr[x1+iy] = max (0.5, Memr[x1+iy]) + c1 - xs[iy] + Memr[x2+iy] = min (nc + 0.49, Memr[x2+iy]) + c1 - xs[iy] + ix1 = nint (Memr[x1+iy]) + ix2 = nint (Memr[x2+iy]) + Memr[y+iy] = iy + do ix = 1, nx { + if (ix < ix1 || ix > ix2) + Memb[xreject] = false + else + Memb[xreject] = true + xreject = xreject + 1 + } + } + + # Estimate spectrum by summing across the aperture with partial + # pixel estimates at the aperture edges. The initial profile + # estimates are obtained by normalizing by the spectrum estimate. + # Profiles where the spectrum is below sky are set to zero. + + call aclrr (profile, nx * ny) + nrej = 0 + do iy = 1, ny { + if (Memr[x1+iy] >= Memr[x2+iy]) { + Memr[spec+iy] = 0. + do ix = 1, nx + profile[iy,ix] = 0. + next + } + + call ap_asifit (dbuf+(iy+ys-1-l1)*nc, nc, xs[iy]-c1+1, + Memr[x1+iy]-c1+xs[iy], Memr[x2+iy]-c1+xs[iy], data, asi) + if (sbuf != NULL) + sky = sbuf + (iy - 1) * nx - 1 + call ap_edge (asi, Memr[x1+iy]+1, Memr[x2+iy]+1, wt1, wt2) + ix1 = nint (Memr[x1+iy]) + ix2 = nint (Memr[x2+iy]) + s = 0. + do ix = ix1, ix2 { + if (!IS_INDEF(saturation)) + if (Memr[data+ix] > saturation) { + s = 0. + nrej = nrej + 1 + break; + } + dat = Memr[data+ix] - Memr[sky+ix] + if (ix1 == ix2) + dat = wt1 * dat + else if (ix == ix1) + dat = wt1 * dat + else if (ix == ix2) + dat = wt2 * dat + s = s + dat + } + + if (s > 0.) { + do ix = ix1, ix2 + profile[iy,ix] = max (0., (Memr[data+ix]-Memr[sky+ix])/s) + } else { + do ix = ix1, ix2 + profile[iy,ix] = 0. + } + Memr[spec+iy] = s + } + + if (nrej == ny) + call error (1, "All profiles contain saturated pixels") + else if (nrej > 0) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT: %d profiles with saturated pixels in aperture %d") + call pargi (nrej) + call pargi (AP_ID(ap)) + if (nrej < ny / 3) + call ap_log (Memc[str], YES, NO, NO) + else + call ap_log (Memc[str], YES, NO, YES) + } + + # Smooth the profile and possibly reject deviant pixels. + nreject = 0 + tfac = 2. + do i = 0, niterate { + + # Estimate profile. + if (xs1 == xs2) + call ap_horne (im, cv, dbuf, nc, nl, c1, l1, Memr[spec+1], sbuf, + svar, Memb[reject], profile, nx, ny, xs, ys, + Memr[x1+1], Memr[x2+1]) + else + call ap_marsh (im, dbuf, nc, nl, c1, l1, Memr[spec+1], sbuf, + svar, Memb[reject], profile, nx, ny, xs, ys, + Memr[x1+1], Memr[x2+1]) + + if (i == niterate) + break + + # Reject pixels. The rejection threshold is based on the overall + # chi square. Pixels are rejected on the basis of the current + # chi square and the largest residual not rejected is compared + # against the final chi square to possibly trigger another round + # of rejections. + + chisq = 0.; nsum = 0; ixrej = 0; iyrej = 0; rrej = 0.; nrej = 0 + do iy = 1, ny { + s = Memr[spec+iy] + if (s <= 0.) + next + call ap_asifit (dbuf+(iy+ys-1-l1)*nc, nc, xs[iy]-c1+1, + Memr[x1+iy]-c1+xs[iy], Memr[x2+iy]-c1+xs[iy], data, asi) + if (sbuf != NULL) { + sky = sbuf + (iy - 1) * nx - 1 + var0 = rdnoise + Memr[svar+iy-1] + } + call ap_edge (asi, Memr[x1+iy]+1, Memr[x2+iy]+1, wt1, wt2) + xreject = reject + (iy - 1) * nx - 1 + ix1 = nint (Memr[x1+iy]) + ix2 = nint (Memr[x2+iy]) + do ix = ix1, ix2 { + if (Memb[xreject+ix]) { + nsum = nsum + 1 + predict = max (0., s * profile[iy,ix] + Memr[sky+ix]) + var = max (vmin, var0 + predict) + resid = (Memr[data+ix] - predict) / sqrt (var) + chisq = chisq + resid**2 + if (resid < -tfac*lsigma || resid > tfac*usigma) { + if (ix < ix1 || ix > ix2) + p = 0. + else if (ix1 == ix2) + p = wt1 + else if (ix == ix1) + p = wt1 + else if (ix == ix2) + p = wt2 + else + p = 1 + Memr[spec+iy] = Memr[spec+iy] - + p * (Memr[data+ix] - predict) + nrej = nrej + 1 + Memb[xreject+ix] = false + } else if (abs (resid) > abs (rrej)) { + rrej = resid + if (ix < ix1 || ix > ix2) + p = 0. + else if (ix1 == ix2) + p = wt1 + else if (ix == ix1) + p = wt1 + else if (ix == ix2) + p = wt2 + else + p = 1 + dat = p * (Memr[data+ix] - predict) + ixrej = ix + iyrej = iy + } + } + } + } + + if (nsum == 0) + call error (1, "All pixels rejected") + tfac = sqrt (chisq / nsum) + if (rrej < -tfac * lsigma || rrej > tfac * usigma) { + Memr[spec+iyrej] = Memr[spec+iyrej] - dat + xreject = reject + (iyrej - 1) * nx - 1 + Memb[xreject+ixrej] = false + nrej = nrej + 1 + } + + nreject = nreject + nrej + if (nrej == 0) + break + } + + # These plots are too big for production work but can be turned on + # for debugging. + + call ap_popen (gp, fd, "fits") + if (gp != NULL) { + ix1 = xs1 + ix2 = xs2 + nx - 1 + if (xs1 != xs2) { + ix1 = ix1 + 1 + ix2 = ix2 - 1 + } + do ix = ix1, ix2 { + nrej = 0 + do iy = 1, ny { + i = ix - xs[iy] + 1 + if (i < 1 || i > nx) + next + if (Memr[spec+iy] <= 0.) + next + data = dbuf + (iy + ys - 1 - l1) * nc + ix - c1 - 1 + if (sbuf != NULL) + s = Memr[sbuf+(iy-1)*nx+i-1] + else + s = Memr[sky+i] + nrej = nrej + 1 + Memr[y+nrej] = iy + ys - 1 + Memr[x1+nrej] = max (-.1, min (1.1, + (Memr[data+1] - s) / Memr[spec+iy])) + Memr[x2+nrej] = profile[iy,i] + } + call gclear (gp) + call gascale (gp, Memr[x1+1], nrej, 2) + call grscale (gp, Memr[x2+1], nrej, 2) + call gswind (gp, Memr[y+1], Memr[y+nrej], INDEF, INDEF) + if (AP_AXIS(ap) == 1) { + call sprintf (Memc[str], SZ_LINE, "Column %d") + call pargi (ix) + call glabax (gp, Memc[str], "Line", "Profile") + } else { + call sprintf (Memc[str], SZ_LINE, "Line %d") + call pargi (ix) + call glabax (gp, Memc[str], "Column", "Profile") + } + call gpmark (gp, Memr[y+1], Memr[x1+1], nrej, GM_POINT, 1., 1.) + call gpline (gp, Memr[y+1], Memr[x2+1], nrej) + } + } + call ap_pclose (gp, fd) + + # Log the number of rejected pixels. + if (clean) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT: %d pixels rejected for profile from aperture %d") + call pargi (nreject) + call pargi (AP_ID(ap)) + call ap_log (Memc[str], YES, NO, NO) + } + + call sfree (sp) +end + + +# AP_HORNE -- Determine profile by fitting a low order function parallel to +# dispersion along image lines or columns after dividing by a spectrum +# estimate. An initial profile estimate and a rejection array are +# required for setting the weights. This is a straightforward algorithm +# similar to images.fit1d except that it is noninteractive. The fitting +# function is fixed at a cubic spline and the number of pieces is set by +# the amount of tilt such that there is one cubic spline piece per +# passage across the tilted spectrum plus an amount based on the order +# of the tracing function. It is named after Keith Horne +# since this is what is outlined in his paper. The profile array is used +# cleverly to minimize memory requirements. The storage order of the +# profile array, which is transposed relative to the data, is determined +# by this procedure. + +procedure ap_horne (im, cvtrace, dbuf, nc, nl, c1, l1, spec, sbuf, svar, reject, + profile, nx, ny, xs, ys, x1, x2) + +pointer im # IMIO pointer +pointer cvtrace # Trace pointer +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Spectrum estimate +pointer sbuf # Sky values (NULL if none) +pointer svar # Sky variances +bool reject[nx,ny] # Rejection flags +real profile[ny,nx] # Initial profile in, improved profile out +int nx, ny # Size of profile array +int xs[ny], ys # Origin of profile array +real x1[ny], x2[ny] # Aperture limits in profile array + +int cvtype # Curfit type +int order # Order of curfit function. +real rdnoise # Readout noise in RMS data numbers. + +int ix, iy, ierr +real p, s, sk, var, vmin, var0, wmin +pointer sp, y, w, cv, dbuf1, data, sky + +#int apgeti() +int cvstati() +real apgimr() +errchk salloc, apgimr + +begin + call smark (sp) + call salloc (y, ny, TY_REAL) + call salloc (w, ny, TY_REAL) + + # Get CL parameters + #cvtype = apgeti ("e_function") + #order = apgeti ("e_order") + rdnoise = apgimr ("readnoise", im) ** 2 + + # Initialize. + call alimr (x1, ny, p, s) + cvtype = SPLINE3 + order = int (s - p + 1) + max (0, cvstati (cvtrace, CVNCOEFF) - 2) + #order = min (20, order) + order = 2 * order + call cvinit (cv, cvtype, order, 1., real (ny)) + do iy = 1, ny + Memr[y+iy-1] = iy + if (rdnoise == 0.) + vmin = 1. + else + vmin = rdnoise + dbuf1 = dbuf + (ys - l1 - 1) * nc - c1 - 1 + if (sbuf == NULL) { + sk = 0. + var0 = rdnoise + } + + # For each line parallel to the dispersion divide by a spectrum + # estimate and then fit the smoothing function. Use the input + # profile and rejection array to set the weights. + + do ix = 1, nx { + data = dbuf1 + ix + if (sbuf != NULL) + sky = sbuf - nx - 1 + ix + wmin = MAX_REAL + do iy = 1, ny { + s = spec[iy] + if (s > 0. && reject[ix,iy]) { + if (sbuf != NULL) { + sk = Memr[sky+iy*nx] + var0 = rdnoise + Memr[svar+iy-1] + } + p = profile[iy,ix] + var = max (vmin, var0 + max (0., s * p + sk)) + var = (s ** 2) / var + wmin = min (wmin, var) + Memr[w+iy-1] = var + profile[iy,ix] = (Memr[data+iy*nc+xs[iy]] - sk) / s + } else + Memr[w+iy-1] = 0. + } + if (wmin == MAX_REAL) + call amovkr (1., Memr[w], ny) + else + call amaxkr (Memr[w], wmin / 10., Memr[w], ny) + call cvfit (cv, Memr[y], profile[1,ix], Memr[w], ny, WTS_USER, ierr) + call cvvector (cv, Memr[y], profile[1,ix], ny) + call amaxkr (profile[1,ix], 0., profile[1,ix], ny) + } + + call cvfree (cv) + call sfree (sp) +end + + +# AP_MARSH -- Determine profile by Marsh algorithm (PASP V101, P1032, 1989). +# This algorithm fits low order polynomials to weighted points sampled +# at uniform intervals parallel to the aperture trace. The polynomials +# are coupled through the weights and so requires a 2D matrix inversion. +# This is a relatively slow algorithm but does provide low order smoothing +# for arbitrary profile shapes in highly tilted spectra. An estimate +# of the profile, a rejection array, sky and sky variance, and aperture +# limit arrays are required. + +procedure ap_marsh (im, dbuf, nc, nl, c1, l1, spec, sbuf, svar, reject, + profile, nx, ny, xs, ys, x1, x2) + +pointer im # IMIO pointer +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Spectrum estimate +pointer sbuf # Sky values (NULL if none) +pointer svar # Sky variances +bool reject[nx,ny] # Rejection flags +real profile[ny,nx] # Initial profile in, improved profile out +int nx, ny # Size of profile array +int xs[ny], ys # Origin of profile array +real x1[ny], x2[ny] # Aperture limits in profile array + +real spix # Polynomial pixel separation +int npols # Number of polynomials +int order # Order of function. +real rdnoise # Readout noise in RMS data numbers. + +int il, jl, kl, ll, ix, iy, ix1, ix2, nside, nadd +int ip, ip1, ip2, index1, index2, index3 +real p, s, s2, dat, sk, var, vmin, var0 +real dx0, dx1, dx2, dx3, dx4, xj, xk, xt, xz, qj, qk, xadd +double sum1, sum2 +pointer sp, work, work1, work2, work3, work4, ysum, data, sky + +int apgeti() +real apgetr(), apgimr() +errchk salloc, apgimr + +begin + # Get CL parameters + #npols = apgeti ("npols") + spix = apgetr ("polysep") + order = apgeti ("polyorder") + rdnoise = apgimr ("readnoise", im) ** 2 + + # Set dimensions. + npols = (x2[1] - x1[1] + 2) / spix + spix = (x2[1] - x1[1] + 2) / real (npols) + nside = npols * order + nadd = nside * nside + if (spix > 1.) + call error (4, "Polynomial separation too large") + + # Allocate memory. One index pointers. + call smark (sp) + call salloc (work, nadd+3*nside, TY_REAL) + call salloc (work4, nside, TY_INT) + call salloc (ysum, ny, TY_REAL) + work = work - 1 + work1 = work + nadd + work2 = work1 + nside + work3 = work2 + nside + work4 = work4 - 1 + ysum=ysum-1 + if (sbuf == NULL) { + call salloc (sky, nx, TY_REAL) + sky = sky - 1 + } + + # Initialize. + call aclrr (Memr[work+1], nadd+3*nside) + call aclri (Memi[work4+1], nside) + if (rdnoise == 0.) + vmin = 1. + else + vmin = rdnoise + if (sbuf == NULL) { + call aclrr (Memr[sky+1], nx) + var0 = rdnoise + } + + # Factors for weights. + dx0 = 0.5 - spix + dx1 = abs (dx0) + dx2 = 1. - (dx0 / spix) ** 2 + dx3 = 0.5 + spix + dx4 = sqrt (2.) * spix + + # Accumulate least terms for least squares matrix equation AX = B. + + # First accumulate B. + do jl = 0, npols-1 { + do iy = 1, ny { + if (spec[iy] <= 0.) + next + + xj = x1[iy] - 1 + spix * (real (jl) + 0.5) + ix1 = nint (xj - spix) + ix2 = nint (xj + spix) + if (ix1 < 1 || ix2 > nx) { + Memr[ysum+iy] = 0. + next + } + + data = dbuf + (iy + ys - 1 - l1) * nc + xs[iy] - c1 - 1 + if (sbuf != NULL) { + sky = sbuf + (iy - 1) * nx - 1 + var0 = rdnoise + Memr[svar+iy-1] + } + + # Evaluate qj, the contribution of polynomial number jl+1 + # for the pixel ix1,jj. Four cases are considered. The + # first two account for the triangular interpolation + # function partially overlapping a pixel, on one side + # only. The third is for the function wholly inside a + # pixel, and finally for the pixel wholly covered by the + # interpolation function. + + s = spec[iy] + sum1 = 0. + do ix = ix1, ix2 { + if (!reject[ix,iy]) + next + p = profile[iy,ix] + sk = Memr[sky+ix] + dat = Memr[data+ix] - sk + var = max (vmin, var0 + max (0., s * p + sk)) + + xz = xj - real (ix) + xt = abs (xz) + if (xt >= dx1) { + if (xt >= 0.5) + qj = ((xt - dx3) / dx4) ** 2 + else + qj = 1.- ((xt - dx0) / dx4) ** 2 + + } else if (xt <= dx0) + qj = 1. + else + qj = dx2 - (xz / spix) ** 2 + sum1 = sum1 + qj * s * dat / var + } + Memr[ysum+iy] = sum1 + } + + index1 = order * jl + do il = 1, order { + sum1 = 0. + ip = il - 1 + do iy = 1, ny + if (spec[iy] > 0.) + sum1 = sum1 + Memr[ysum+iy] * ((real (iy) / ny) ** ip) + Memr[work1+index1+il] = sum1 + } + } + + # Now accumulate matrix A. Since it is symmetric we only need to + # evaluate half of it. Since it is banded we only need to evaluate + # contribution if two polynomial terms can be affected by the same + # pixel. + + ip1 = nside - 1 + ip2 = order * ip1 + do jl = 0, npols-1 { + do kl = 0, jl { + if (spix * (jl - kl - 2) > 0.) + next + do iy = 1, ny { + if (spec[iy] <= 0.) + next + if (sbuf != NULL) { + sky = sbuf + (iy - 1) * nx - 1 + var0 = rdnoise + Memr[svar+iy-1] + } + + # Compute left and right limits of polynomials jl+1 + # and kl+1 for this value of y Evaluate sum over row + # of qj[jl+1] times qj[kl+1] where qj[i] is fraction + # of polynomial i which contributes to to pixel ix,jj. + + xj = x1[iy] - 1 + spix * (real (jl) + 0.5) + xk = x1[iy] - 1 + spix * (real (kl) + 0.5) + ix1 = nint (xj - spix) + ix2 = nint (xk + spix) + + if (ix2 < ix1 || ix1 < 1 || ix2 > nx) { + Memr[ysum+iy] = 0. + next + } + + s = spec[iy] + s2 = s * s + sum1 = 0. + do ix = ix1, ix2 { + if (reject[ix,iy]) { + p = profile[iy,ix] + sk = Memr[sky+ix] + var = max (vmin, var0 + max (0., s * p + sk)) + + xz = xj - real (ix) + xt = abs (xz) + if (xt >= dx1) { + if (xt >= 0.5) + qj = ((xt-dx3)/dx4)**2 + else + qj = 1.- ((xt-dx0)/dx4)**2 + } else if (xt <= dx0) + qj = 1. + else + qj = dx2 - (xz / spix) ** 2 + if (kl != jl) { + xz = xk - real (ix) + xt = abs (xz) + if (xt >= dx1) { + if (xt >= 0.5) + qk = ((xt-dx3)/dx4)**2 + else + qk = 1.-((xt-dx0)/dx4)**2 + } else if (xt <= dx0) + qk = 1. + else + qk = dx2 - (xz / spix) ** 2 + } else + qk = qj + sum1 = sum1 + qj * qk * s2 / var + } + } + Memr[ysum+iy] = sum1 + } + + do il = 1, order { + do ll = 1, il { + sum1 = 0. + ip = il + ll - 2 + do iy = 1, ny + if (spec[iy] > 0.) + sum1 = sum1 + + Memr[ysum+iy] * ((real (iy) / ny)**ip) + index1 = nside * (order*jl+il-1) + order * kl + ll + Memr[work+index1] = sum1 + if (ll != il) { + ip = ip1 * (ll - il) + index2 = index1 + ip + Memr[work+index2] = sum1 + } else + index2 = index1 + if (kl != jl) { + index3 = index2 + ip2 * (kl - jl) + Memr[work+index3] = sum1 + if (ll != il) + Memr[work+index3-ip] = sum1 + } + } + } + } + } + + # Solve matrix equation AX = B for X. A is a real symmetric, + # positive definite matrix, dimension (order*npols)**2. X is + # the vector representing the coefficients fitted to the + # normalized profile. Coefficients are reordered for later speed. + + call hfti (Memr[work+1], nside, nside, nside, Memr[work1+1], 1, 1, + 0.01, ip, p, Memr[work2+1], Memr[work3+1], Memi[work4+1]) + + do jl = 1, order { + do il = 1, npols { + index1 = order * (il - 1) + jl + index2 = npols * (jl - 1) + il + Memr[work+index2] = Memr[work1+index1] + } + } + + # Evaluate fit and make profile positive only. + do iy = 1, ny { + ix1 = nint (x1[iy]) + ix2 = nint (x2[iy]) + xadd = x1[iy] - 1 + s = 0. + do ix = 1, nx { + xj = real (ix) - xadd - 0.5 + xk = real (ix) - xadd + 0.5 + ip1 = int (xj / spix + 0.5) + ip2 = int (xk / spix + 1.5) + ip1 = max (1, min (ip1, npols)) + ip2 = max (1, min (ip2, npols)) + sum1 = 0. + do jl = 0, order-1 { + index1 = npols * jl + sum2 = 0. + do il = ip1, ip2 { + xz = xadd + spix * (real (il-1) + 0.5) - real (ix) + xt = abs (xz) + if (xt >= dx1) { + if (xt >= 0.5) + qj = ((xt - dx3) / dx4) ** 2 + else + qj = 1. - ((xt - dx0) / dx4) ** 2 + } else if (xt <= dx0) + qj = 1. + else + qj = dx2 - (xz / spix) ** 2 + sum2 = sum2 + qj * Memr[work+index1+il] + } + sum1 = sum1 + sum2 * ((real (iy)/ ny) ** jl) + } + profile[iy,ix] = max (0.d0, sum1) + } + } + + call sfree (sp) +end |