aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imcoords/src/ccfunc.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imcoords/src/ccfunc.x')
-rw-r--r--pkg/images/imcoords/src/ccfunc.x639
1 files changed, 639 insertions, 0 deletions
diff --git a/pkg/images/imcoords/src/ccfunc.x b/pkg/images/imcoords/src/ccfunc.x
new file mode 100644
index 00000000..9f60498a
--- /dev/null
+++ b/pkg/images/imcoords/src/ccfunc.x
@@ -0,0 +1,639 @@
+include <imhdr.h>
+include <math.h>
+include <math/gsurfit.h>
+include <mwset.h>
+include <pkg/skywcs.h>
+
+
+# CC_RPROJ -- Read the projection parameters from a file into an IRAF string
+# containing the projection type followed by an MWCS WAT string, e.g
+# "zpn projp1=value projp2=value" .
+
+int procedure cc_rdproj (fd, projstr, maxch)
+
+int fd #I the input file containing the projection parameters
+char projstr[ARB] #O the output projection parameters string
+int maxch #I the maximum size of the output projection string
+
+int projection, op
+pointer sp, keyword, value, param
+int fscan(), nscan(), strdic(), gstrcpy()
+
+begin
+ projstr[1] = EOS
+ if (fscan (fd) == EOF)
+ return (0)
+
+ call smark (sp)
+ call salloc (keyword, SZ_FNAME, TY_CHAR)
+ call salloc (value, SZ_FNAME, TY_CHAR)
+ call salloc (param, SZ_FNAME, TY_CHAR)
+
+ call gargwrd (Memc[keyword], SZ_FNAME)
+ projection = strdic (Memc[keyword], Memc[keyword], SZ_FNAME,
+ WTYPE_LIST)
+ if (projection <= 0 || projection == WTYPE_LIN || nscan() == 0) {
+ call sfree (sp)
+ return (0)
+ }
+
+ # Copy the projection function into the projection string.
+ op = 1
+ op = op + gstrcpy (Memc[keyword], projstr[op], maxch)
+
+ # Copy the keyword value pairs into the projection string.
+ while (fscan(fd) != EOF) {
+ call gargwrd (Memc[keyword], SZ_FNAME)
+ call gargwrd (Memc[value], SZ_FNAME)
+ if (nscan() != 2)
+ next
+ call sprintf (Memc[param], SZ_FNAME, " %s = %s")
+ call pargstr (Memc[keyword])
+ call pargstr (Memc[value])
+ op = op + gstrcpy (Memc[param], projstr[op], maxch - op + 1)
+ }
+
+ call sfree (sp)
+
+ return (projection)
+end
+
+
+define NEWCD Memd[ncd+(($2)-1)*ndim+($1)-1]
+
+# CC_WCSIM -- Update the image world coordinate system.
+
+procedure cc_wcsim (im, coo, projection, lngref, latref, sx1, sy1, transpose)
+
+pointer im #I the pointer to the input image
+pointer coo #I the pointer to the coordinate structure
+char projection[ARB] #I the sky projection geometry
+double lngref, latref #I the position of the reference point.
+pointer sx1, sy1 #I pointer to linear surfaces
+bool transpose #I transpose the wcs
+
+int ndim, naxes, ax1, ax2, axmap, wtype
+double xshift, yshift, a, b, c, d, denom, xpix, ypix, tlngref, tlatref
+pointer mw, sp, str, r, w, cd, ltm, ltv, iltm, nr, ncd, axes, axno, axval
+int mw_stati(), sk_stati(), strdic()
+pointer mw_openim()
+
+begin
+ mw = mw_openim (im)
+ ndim = mw_stati (mw, MW_NPHYSDIM)
+
+ # Allocate working memory for the vectors and matrices.
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call salloc (axno, IM_MAXDIM, TY_INT)
+ call salloc (axval, IM_MAXDIM, TY_INT)
+ call salloc (axes, IM_MAXDIM, TY_INT)
+ call salloc (r, ndim, TY_DOUBLE)
+ call salloc (w, ndim, TY_DOUBLE)
+ call salloc (cd, ndim * ndim, TY_DOUBLE)
+ call salloc (ltm, ndim * ndim, TY_DOUBLE)
+ call salloc (ltv, ndim, TY_DOUBLE)
+ call salloc (iltm, ndim * ndim, TY_DOUBLE)
+ call salloc (nr, ndim, TY_DOUBLE)
+ call salloc (ncd, ndim * ndim, TY_DOUBLE)
+
+ # Compute the original logical to world transformation.
+ call mw_gaxmap (mw, Memi[axno], Memi[axval], ndim)
+ call mw_gltermd (mw, Memd[ltm], Memd[ltv], ndim)
+
+ # Get the axis map.
+ call mw_gaxlist (mw, 03B, Memi[axes], naxes)
+ axmap = mw_stati (mw, MW_USEAXMAP)
+ ax1 = Memi[axes]
+ ax2 = Memi[axes+1]
+
+ # Set the system.
+ iferr (call mw_newsystem (mw, "image", ndim))
+ ;
+
+ # Set the axes and projection type.
+ if (projection[1] == EOS) {
+ call mw_swtype (mw, Memi[axes], ndim, "linear", "")
+ } else {
+ call mw_swtype (mw, Memi[axes], ndim, projection,
+ "axis 1: axtype=ra axis 2: axtype=dec")
+ }
+
+ # Compute the new referemce point.
+ switch (sk_stati(coo, S_NLNGUNITS)) {
+ case SKY_DEGREES:
+ tlngref = lngref
+ case SKY_RADIANS:
+ tlngref = RADTODEG(lngref)
+ case SKY_HOURS:
+ tlngref = 15.0d0 * lngref
+ default:
+ tlngref = lngref
+ }
+ switch (sk_stati(coo, S_NLATUNITS)) {
+ case SKY_DEGREES:
+ tlatref = latref
+ case SKY_RADIANS:
+ tlatref = RADTODEG(latref)
+ case SKY_HOURS:
+ tlatref = 15.0d0 * latref
+ default:
+ tlatref = latref
+ }
+ if (! transpose) {
+ Memd[w+ax1-1] = tlngref
+ Memd[w+ax2-1] = tlatref
+ } else {
+ Memd[w+ax1-1] = tlatref
+ Memd[w+ax2-1] = tlngref
+ }
+
+
+ # Fetch the linear coefficients of the fit.
+ call geo_gcoeffd (sx1, sy1, xshift, yshift, a, b, c, d)
+
+ # Compute the new reference pixel.
+ denom = a * d - c * b
+ if (denom == 0.0d0)
+ xpix = INDEFD
+ else
+ xpix = (b * yshift - d * xshift) / denom
+ if (denom == 0.0d0)
+ ypix = INDEFD
+ else
+ ypix = (c * xshift - a * yshift) / denom
+ Memd[nr+ax1-1] = xpix
+ Memd[nr+ax2-1] = ypix
+
+ # Compute the new CD matrix.
+ if (! transpose) {
+ NEWCD(ax1,ax1) = a / 3600.0d0
+ NEWCD(ax1,ax2) = c / 3600.0d0
+ NEWCD(ax2,ax1) = b / 3600.0d0
+ NEWCD(ax2,ax2) = d / 3600.0d0
+ } else {
+ NEWCD(ax1,ax1) = c / 3600.0d0
+ NEWCD(ax1,ax2) = a / 3600.0d0
+ NEWCD(ax2,ax1) = d / 3600.0d0
+ NEWCD(ax2,ax2) = b / 3600.0d0
+ }
+
+ # Reset the axis map.
+ call mw_seti (mw, MW_USEAXMAP, axmap)
+
+ # Recompute and store the new wcs if update is enabled.
+ call mw_saxmap (mw, Memi[axno], Memi[axval], ndim)
+ if (sk_stati (coo, S_PIXTYPE) == PIXTYPE_PHYSICAL) {
+ call mw_swtermd (mw, Memd[nr], Memd[w], Memd[ncd], ndim)
+ } else {
+ call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], ndim)
+ call mwinvertd (Memd[ltm], Memd[iltm], ndim)
+ call asubd (Memd[nr], Memd[ltv], Memd[r], ndim)
+ call mwvmuld (Memd[iltm], Memd[r], Memd[nr], ndim)
+ call mw_swtermd (mw, Memd[nr], Memd[w], Memd[cd], ndim)
+ }
+
+ # Save the fit.
+ if (! transpose) {
+ call sk_seti (coo, S_PLNGAX, ax1)
+ call sk_seti (coo, S_PLATAX, ax2)
+ } else {
+ call sk_seti (coo, S_PLNGAX, ax2)
+ call sk_seti (coo, S_PLATAX, ax1)
+ }
+ call sk_saveim (coo, mw, im)
+ call mw_saveim (mw, im)
+ call mw_close (mw)
+
+ # Force the CDELT keywords to update. This will be unecessary when
+ # mwcs is updated to deal with non-quoted and / or non left-justified
+ # CTYPE keywords..
+ wtype = strdic (projection, Memc[str], SZ_FNAME, WTYPE_LIST)
+ if (wtype > 0)
+ call sk_seti (coo, S_WTYPE, wtype)
+ call sk_ctypeim (coo, im)
+
+ # Reset the fit. This will be unecessary when wcs is updated to deal
+ # with non-quoted and / or non left-justified CTYPE keywords.
+ call sk_seti (coo, S_WTYPE, 0)
+ call sk_seti (coo, S_PLNGAX, 0)
+ call sk_seti (coo, S_PLATAX, 0)
+
+ call sfree (sp)
+end
+
+
+# CC_NWCSIM -- Update the image world coordinate system.
+
+procedure cc_nwcsim (im, coo, projection, lngref, latref, sx1, sy1, sx2, sy2,
+ transpose)
+
+pointer im #I the pointer to the input image
+pointer coo #I the pointer to the coordinate structure
+char projection[ARB] #I the sky projection geometry
+double lngref, latref #I the position of the reference point.
+pointer sx1, sy1 #I pointer to linear surfaces
+pointer sx2, sy2 #I pointer to distortion surfaces
+bool transpose #I transpose the wcs
+
+int l, i, ndim, naxes, ax1, ax2, axmap, wtype, szatstr
+double xshift, yshift, a, b, c, d, denom, xpix, ypix, tlngref, tlatref
+pointer mw, sp, r, w, cd, ltm, ltv, iltm, nr, ncd, axes, axno, axval
+pointer projstr, projpars, wpars, mwnew, atstr
+bool streq()
+int mw_stati(), sk_stati(), strdic(), strlen(), itoc()
+pointer mw_openim(), mw_open()
+errchk mw_gwattrs(), mw_newsystem()
+
+begin
+ # Open the image wcs and determine its size.
+ mw = mw_openim (im)
+ ndim = mw_stati (mw, MW_NPHYSDIM)
+
+ # Allocate working memory for the wcs attributes, vectors, and
+ # matrices.
+ call smark (sp)
+ call salloc (projstr, SZ_FNAME, TY_CHAR)
+ call salloc (projpars, SZ_LINE, TY_CHAR)
+ call salloc (wpars, SZ_LINE, TY_CHAR)
+ call salloc (axno, IM_MAXDIM, TY_INT)
+ call salloc (axval, IM_MAXDIM, TY_INT)
+ call salloc (axes, IM_MAXDIM, TY_INT)
+ call salloc (r, ndim, TY_DOUBLE)
+ call salloc (w, ndim, TY_DOUBLE)
+ call salloc (cd, ndim * ndim, TY_DOUBLE)
+ call salloc (ltm, ndim * ndim, TY_DOUBLE)
+ call salloc (ltv, ndim, TY_DOUBLE)
+ call salloc (iltm, ndim * ndim, TY_DOUBLE)
+ call salloc (nr, ndim, TY_DOUBLE)
+ call salloc (ncd, ndim * ndim, TY_DOUBLE)
+
+ # Open the new wcs and set the system type.
+ mwnew = mw_open (NULL, ndim)
+ call mw_gsystem (mw, Memc[projstr], SZ_FNAME)
+ iferr {
+ call mw_newsystem (mw, "image", ndim)
+ } then {
+ call mw_newsystem (mwnew, Memc[projstr], ndim)
+ } else {
+ call mw_newsystem (mwnew, "image", ndim)
+ }
+
+ # Set the LTERM.
+ call mw_gltermd (mw, Memd[ltm], Memd[ltv], ndim)
+ call mw_sltermd (mwnew, Memd[ltm], Memd[ltv], ndim)
+
+ # Store the old axis map for later use.
+ call mw_gaxmap (mw, Memi[axno], Memi[axval], ndim)
+
+ # Get the celestial coordinate axes list.
+ call mw_gaxlist (mw, 03B, Memi[axes], naxes)
+ axmap = mw_stati (mw, MW_USEAXMAP)
+ ax1 = Memi[axes]
+ ax2 = Memi[axes+1]
+
+ # Set the axes and projection type for the celestial coordinate
+ # axes. Don't worry about the fact that the axes may in fact be
+ # glon and glat, elon and elat, or slon and slat, instead of
+ # ra and dec. This will be fixed up later.
+ if (projection[1] == EOS) {
+ call mw_swtype (mwnew, Memi[axes], ndim, "linear", "")
+ } else {
+ call sscan (projection)
+ call gargwrd (Memc[projstr], SZ_FNAME)
+ call gargstr (Memc[projpars], SZ_LINE)
+ call sprintf (Memc[wpars], SZ_LINE,
+ "axis 1: axtype = ra %s axis 2: axtype = dec %s")
+ call pargstr (Memc[projpars])
+ call pargstr (Memc[projpars])
+ if (streq (Memc[projstr], "tnx") && sx2 == NULL && sy2 == NULL)
+ call strcpy ("tan", Memc[projstr], SZ_FNAME)
+ call mw_swtype (mwnew, Memi[axes], ndim, Memc[projstr], Memc[wpars])
+ }
+
+ # Copy the attributes of the remaining axes to the new wcs.
+ szatstr = SZ_LINE
+ call malloc (atstr, szatstr, TY_CHAR)
+ do l = 1, ndim {
+ if (l == ax1 || l == ax2)
+ next
+ iferr {
+ call mw_gwattrs (mw, l, "wtype", Memc[projpars], SZ_LINE)
+ } then {
+ call mw_swtype (mwnew, l, 1, "linear", "")
+ } else {
+ call mw_swtype (mwnew, l, 1, Memc[projpars], "")
+ }
+ for (i = 1; ; i = i + 1) {
+ if (itoc (i, Memc[projpars], SZ_LINE) <= 0)
+ Memc[projpars] = EOS
+ repeat {
+ iferr (call mw_gwattrs (mw, l, Memc[projpars],
+ Memc[atstr], szatstr))
+ Memc[atstr] = EOS
+ if (strlen(Memc[atstr]) < szatstr)
+ break
+ szatstr = szatstr + SZ_LINE
+ call realloc (atstr, szatstr, TY_CHAR)
+ }
+ if (Memc[atstr] == EOS)
+ break
+ call mw_swattrs (mwnew, l, Memc[projpars], Memc[atstr])
+ }
+ }
+ call mfree (atstr, TY_CHAR)
+
+ # Compute the new referemce point.
+ switch (sk_stati(coo, S_NLNGUNITS)) {
+ case SKY_DEGREES:
+ tlngref = lngref
+ case SKY_RADIANS:
+ tlngref = RADTODEG(lngref)
+ case SKY_HOURS:
+ tlngref = 15.0d0 * lngref
+ default:
+ tlngref = lngref
+ }
+ switch (sk_stati(coo, S_NLATUNITS)) {
+ case SKY_DEGREES:
+ tlatref = latref
+ case SKY_RADIANS:
+ tlatref = RADTODEG(latref)
+ case SKY_HOURS:
+ tlatref = 15.0d0 * latref
+ default:
+ tlatref = latref
+ }
+ if (! transpose) {
+ Memd[w+ax1-1] = tlngref
+ Memd[w+ax2-1] = tlatref
+ } else {
+ Memd[w+ax1-1] = tlatref
+ Memd[w+ax2-1] = tlngref
+ }
+ # Fetch the linear coefficients of the fit.
+ call geo_gcoeffd (sx1, sy1, xshift, yshift, a, b, c, d)
+
+ # Compute the new reference pixel.
+ denom = a * d - c * b
+ if (denom == 0.0d0)
+ xpix = INDEFD
+ else
+ xpix = (b * yshift - d * xshift) / denom
+ if (denom == 0.0d0)
+ ypix = INDEFD
+ else
+ ypix = (c * xshift - a * yshift) / denom
+ Memd[nr+ax1-1] = xpix
+ Memd[nr+ax2-1] = ypix
+
+ # Compute the new CD matrix.
+ if (! transpose) {
+ NEWCD(ax1,ax1) = a / 3600.0d0
+ NEWCD(ax1,ax2) = c / 3600.0d0
+ NEWCD(ax2,ax1) = b / 3600.0d0
+ NEWCD(ax2,ax2) = d / 3600.0d0
+ } else {
+ NEWCD(ax1,ax1) = c / 3600.0d0
+ NEWCD(ax1,ax2) = a / 3600.0d0
+ NEWCD(ax2,ax1) = d / 3600.0d0
+ NEWCD(ax2,ax2) = b / 3600.0d0
+ }
+
+ # Recompute and store the new wcs.
+ call mw_saxmap (mwnew, Memi[axno], Memi[axval], ndim)
+ if (sk_stati (coo, S_PIXTYPE) == PIXTYPE_PHYSICAL) {
+ call mw_swtermd (mwnew, Memd[nr], Memd[w], Memd[ncd], ndim)
+ } else {
+ call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], ndim)
+ call mwinvertd (Memd[ltm], Memd[iltm], ndim)
+ call asubd (Memd[nr], Memd[ltv], Memd[r], ndim)
+ call mwvmuld (Memd[iltm], Memd[r], Memd[nr], ndim)
+ call mw_swtermd (mwnew, Memd[nr], Memd[w], Memd[cd], ndim)
+ }
+
+ # Add the second order terms in the form of the wcs attributes
+ # lngcor and latcor. These are not FITS standard and can currently
+ # be understood only by IRAF.
+ if ((streq(Memc[projstr], "zpx") || streq (Memc[projstr], "tnx")) &&
+ (sx2 != NULL || sy2 != NULL)) {
+ if (! transpose)
+ call cc_wcscor (im, mwnew, sx1, sx2, sy1, sy2, "lngcor",
+ "latcor", ax1, ax2)
+ else
+ call cc_wcscor (im, mwnew, sx1, sx2, sy1, sy2, "lngcor",
+ "latcor", ax2, ax1)
+ }
+
+ # Save the fit.
+ if (! transpose) {
+ call sk_seti (coo, S_PLNGAX, ax1)
+ call sk_seti (coo, S_PLATAX, ax2)
+ } else {
+ call sk_seti (coo, S_PLNGAX, ax2)
+ call sk_seti (coo, S_PLATAX, ax1)
+ }
+ call sk_saveim (coo, mwnew, im)
+ call mw_saveim (mwnew, im)
+ call mw_close (mwnew)
+ call mw_close (mw)
+
+ # Force the CTYPE keywords to update. This will be unecessary when
+ # mwcs is updated to deal with non-quoted and / or non left-justified
+ # CTYPE keywords..
+ wtype = strdic (Memc[projstr], Memc[projstr], SZ_FNAME, WTYPE_LIST)
+ if (wtype > 0)
+ call sk_seti (coo, S_WTYPE, wtype)
+ call sk_ctypeim (coo, im)
+
+ # Reset the fit.
+ call sk_seti (coo, S_WTYPE, 0)
+ call sk_seti (coo, S_PLNGAX, 0)
+ call sk_seti (coo, S_PLATAX, 0)
+
+ call sfree (sp)
+end
+
+
+# CC_WCSCOR -- Reformulate the higher order surface fit into a correction
+# term in degrees that can be written into the header as a wcs attribute.
+# This attribute will be written as string containing the surface definition.
+
+procedure cc_wcscor (im, mw, sx1, sx2, sy1, sy2, xiname, etaname, xiaxis,
+ etaaxis)
+
+pointer im #I pointer to the input image
+pointer mw #I pointer to the wcs structure
+pointer sx1, sx2 #I pointer to the linear and distortion xi surfaces
+pointer sy1, sy2 #I pointer to the linear and distortion eta surfaces
+char xiname[ARB] #I the wcs xi correction attribute name
+char etaname[ARB] #I the wcs eta correction attribute name
+int xiaxis #I the xi axis number
+int etaaxis #I the eta axis number
+
+int i, j, function, xxorder, xyorder, xxterms, yxorder, yyorder, yxterms
+int nx, ny, npix, ier
+double sxmin, sxmax, symin, symax, ratio, x, y, xstep, ystep, ximin, ximax
+double etamin, etamax
+pointer sp, xpix, ypix, xilin, etalin, dxi, deta, wgt, nsx2, nsy2
+int dgsgeti()
+double dgsgetd()
+begin
+ if (sx2 == NULL && sy2 == NULL)
+ return
+ if (dgsgeti (sx1, GSTYPE) != dgsgeti (sy1, GSTYPE))
+ return
+
+ # Get the function, xmin, xmax, ymin, and ymax parameters for the
+ # surfaces.
+ function = min (dgsgeti (sx1, GSTYPE), dgsgeti (sy1, GSTYPE))
+ sxmin = max (dgsgetd (sx1, GSXMIN), dgsgetd (sy1, GSXMIN))
+ sxmax = min (dgsgetd (sx1, GSXMAX), dgsgetd (sy1, GSXMAX))
+ symin = max (dgsgetd (sx1, GSYMIN), dgsgetd (sy1, GSYMIN))
+ symax = min (dgsgetd (sx1, GSYMAX), dgsgetd (sy1, GSYMAX))
+
+ # Get the order and cross-terms parameters from the higher order
+ # functions.
+ if (sx2 != NULL) {
+ xxorder = dgsgeti (sx2, GSXORDER)
+ xyorder = dgsgeti (sx2, GSYORDER)
+ xxterms = dgsgeti (sx2, GSXTERMS)
+ } else {
+ xxorder = dgsgeti (sx1, GSXORDER)
+ xyorder = dgsgeti (sx1, GSYORDER)
+ xxterms = dgsgeti (sx1, GSXTERMS)
+ }
+ if (sy2 != NULL) {
+ yxorder = dgsgeti (sy2, GSXORDER)
+ yyorder = dgsgeti (sy2, GSYORDER)
+ yxterms = dgsgeti (sy2, GSXTERMS)
+ } else {
+ yxorder = dgsgeti (sy1, GSXORDER)
+ yyorder = dgsgeti (sy1, GSYORDER)
+ yxterms = dgsgeti (sy1, GSXTERMS)
+ }
+
+ # Choose a reasonable coordinate grid size based on the x and y order
+ # of the fit and the number of rows and columns in the image.
+ ratio = double (IM_LEN(im,2)) / double (IM_LEN(im,1))
+ nx = max (xxorder + 3, yxorder + 3, 10)
+ ny = max (yyorder + 3, xyorder + 3, nint (ratio * 10))
+ npix = nx * ny
+
+ # Allocate some working space.
+ call smark (sp)
+ call salloc (xpix, npix, TY_DOUBLE)
+ call salloc (ypix, npix, TY_DOUBLE)
+ call salloc (xilin, npix, TY_DOUBLE)
+ call salloc (etalin, npix, TY_DOUBLE)
+ call salloc (dxi, npix, TY_DOUBLE)
+ call salloc (deta, npix, TY_DOUBLE)
+ call salloc (wgt, npix, TY_DOUBLE)
+
+ # Compute the grid of x and y points.
+ xstep = (sxmax - sxmin) / (nx - 1)
+ ystep = (symax - symin) / (ny - 1)
+ y = symin
+ npix = 0
+ do j = 1, ny {
+ x = sxmin
+ do i = 1, nx {
+ Memd[xpix+npix] = x
+ Memd[ypix+npix] = y
+ x = x + xstep
+ npix = npix + 1
+ }
+ y = y + ystep
+ }
+
+
+ # Compute the weights
+ call amovkd (1.0d0, Memd[wgt], npix)
+
+ # Evalute the linear surfaces and convert the results from arcseconds
+ # to degrees.
+ call dgsvector (sx1, Memd[xpix], Memd[ypix], Memd[xilin], npix)
+ call adivkd (Memd[xilin], 3600.0d0, Memd[xilin], npix)
+ call alimd (Memd[xilin], npix, ximin, ximax)
+ call dgsvector (sy1, Memd[xpix], Memd[ypix], Memd[etalin], npix)
+ call adivkd (Memd[etalin], 3600.0d0, Memd[etalin], npix)
+ call alimd (Memd[etalin], npix, etamin, etamax)
+
+ # Evalute the distortion surfaces, convert the results from arcseconds
+ # to degrees, and compute new distortion surfaces.
+ if (sx2 != NULL) {
+ call dgsvector (sx2, Memd[xpix], Memd[ypix], Memd[dxi], npix)
+ call adivkd (Memd[dxi], 3600.0d0, Memd[dxi], npix)
+ call dgsinit (nsx2, function, xxorder, xyorder, xxterms,
+ ximin, ximax, etamin, etamax)
+ call dgsfit (nsx2, Memd[xilin], Memd[etalin], Memd[dxi],
+ Memd[wgt], npix, WTS_UNIFORM, ier)
+ call cc_gsencode (mw, nsx2, xiname, xiaxis)
+ } else
+ nsx2 = NULL
+ if (sy2 != NULL) {
+ call dgsvector (sy2, Memd[xpix], Memd[ypix], Memd[deta], npix)
+ call adivkd (Memd[deta], 3600.0d0, Memd[deta], npix)
+ call dgsinit (nsy2, function, yxorder, yyorder, yxterms,
+ ximin, ximax, etamin, etamax)
+ call dgsfit (nsy2, Memd[xilin], Memd[etalin], Memd[deta],
+ Memd[wgt], npix, WTS_UNIFORM, ier)
+ call cc_gsencode (mw, nsy2, etaname, etaaxis)
+ } else
+ nsy2 = NULL
+
+ # Store the string in the mcs structure in the format of a wcs
+ # attribute.
+
+ # Free the new surfaces.
+ if (nsx2 != NULL)
+ call dgsfree (nsx2)
+ if (nsy2 != NULL)
+ call dgsfree (nsy2)
+
+ call sfree (sp)
+end
+
+
+# CC_GSENCODE -- Encode the surface in an mwcs attribute.
+
+procedure cc_gsencode (mw, gs, atname, axis)
+
+pointer mw #I pointer to the mwcs structure
+pointer gs #I pointer to the surface to be encoded
+char atname[ARB] #I attribute name for the encoded surface
+int axis #I axis for which the encode surface is encoded
+
+int i, op, nsave, szatstr, szpar
+pointer sp, coeff, par, atstr
+int dgsgeti(), strlen(), gstrcpy()
+
+begin
+ nsave = dgsgeti (gs, GSNSAVE)
+ call smark (sp)
+ call salloc (coeff, nsave, TY_DOUBLE)
+ call salloc (par, SZ_LINE, TY_CHAR)
+ call dgssave (gs, Memd[coeff])
+
+ szatstr = SZ_LINE
+ call malloc (atstr, szatstr, TY_CHAR)
+ op = 0
+ do i = 1, nsave {
+ call sprintf (Memc[par], SZ_LINE, "%g ")
+ call pargd (Memd[coeff+i-1])
+ szpar = strlen (Memc[par])
+ if (szpar > (szatstr - op)) {
+ szatstr = szatstr + SZ_LINE
+ call realloc (atstr, szatstr, TY_CHAR)
+ }
+ op = op + gstrcpy (Memc[par], Memc[atstr+op], SZ_LINE)
+
+ }
+
+ call mw_swattrs (mw, axis, atname, Memc[atstr])
+ call mfree (atstr, TY_CHAR)
+ call sfree (sp)
+end
+
+
+
generated by cgit (git 2.34.1) at 2025-09-20 10:45:35 -0400