include	<smw.h>


# Evaluate SMW coordinate transform.  These procedures call the
# MWCS procedures unless the WCS is a split MULTISPEC WCS.  In that
# case the appropriate piece needs to be determined and the physical
# line numbers manipulated.  For log sampled spectra conversions
# must be made for EQUISPEC/NDSPEC.  The convention is that coordinates
# are always input and output and linear.  Note that the MULTISPEC
# function driver already takes care of this.
#
# SMW_CTRANR -- N dimensional real coordinate transformation.
# SMW_CTRAND -- N dimensional double coordinate transformation.
# SMW_C1TRANR -- One dimensional real coordinate transformation.
# SMW_C1TRAND -- One dimensional double coordinate transformation.
# SMW_C2TRANR -- Two dimensional real coordinate transformation.
# SMW_C2TRAND -- Two dimensional double coordinate transformation.


$for (rd)
# SMW_CTRAN -- N dimensional coordinate transformation.

procedure smw_ctran$t (ct, p1, p2, ndim)

pointer	ct		#I SMW CT pointer
PIXEL	p1[ndim]	#I Input coordinate
PIXEL	p2[ndim]	#O Output coordinate
int	ndim		#I Dimensionality

int	i, j, format, daxis, aaxis, dtype, naps
pointer	smw, aps
errchk	mw_ctran$t

begin
	if (SMW_NCT(ct) != 1)
	    call error (1, "SMW_CTRAN: Wrong WCS type")

	call amov$t (p1, p2, ndim)

	smw = SMW_SMW(ct)
	format = SMW_FORMAT(smw)
	daxis = SMW_DAXIS(ct)
	aaxis = SMW_AAXIS(ct)
	dtype = SMW_DTYPE(smw)
	naps = SMW_NSPEC(smw)
	aps = SMW_APS(smw)
	switch (format) {
	case SMW_ND, SMW_ES:
	    switch (SMW_CTTYPE(ct)) {
	    case SMW_LW, SMW_PW:
		call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim)
		if (daxis > 0 && dtype == DCLOG)
		    p2[daxis] = 10. ** max (-20$F, min (20$F, p2[daxis]))
		if (aaxis > 0 && format == SMW_ES) {
		    i = max (1, min (naps, nint (p2[aaxis])))
		    p2[aaxis] = Memi[aps+i-1]
		}
	    case SMW_WL, SMW_WP:
		if (daxis > 0 && dtype == DCLOG)
		    p2[daxis] = log10 (p2[daxis])
		if (aaxis > 0 && format == SMW_ES) {
		    j = nint (p2[aaxis])
		    p2[aaxis] = 1
		    do i = 1, naps {
			if (j == Memi[aps+i-1]) {
			    p2[aaxis] = i
			    break
			}
		    }
		}
		call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim)
	    default:
		call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim)
	    }
	case SMW_MS:
	    call mw_ctran$t (SMW_CT(ct,0), p1, p2, ndim)
	}
end


# SMW_C1TRAN -- One dimensional coordinate transformation.

PIXEL procedure smw_c1tran$t (ct, x1)

pointer	ct		#I SMW CT pointer
PIXEL	x1		#I Input coordinate
PIXEL	x2		#O Output coordinate

errchk	mw_ctran$t

begin
	call smw_ctran$t (ct, x1, x2, 1)
	return (x2)
end


# SMW_C2TRAN -- Two dimensional coordinate transformation.

procedure smw_c2tran$t (ct, x1, y1, x2, y2)

pointer	ct		#I SMW CT pointer
PIXEL	x1, y1		#I Input coordinates
PIXEL	x2, y2		#O Output coordinates

PIXEL	p1[2], p2[2]
int	i, j, naps
PIXEL	xp, yp
pointer	aps, smw_ct()
errchk	smw_ct, mw_c2tran$t

begin
	# Unsplit WCS.
	if (SMW_NCT(ct) == 1) {
	    p1[1] = x1
	    p1[2] = y1
	    call smw_ctran$t (ct, p1, p2, 2)
	    x2 = p2[1]
	    y2 = p2[2]
	    return
	}

	# If we get here then we are dealing with a split MULTISPEC WCS.
	# Depending on the systems being transformed there may need to
	# transformation made on the physical coordinate system.

	switch (SMW_CTTYPE(ct)) {
	case SMW_LW:
	    call mw_c2tran$t (SMW_CTL(ct), x1, y1, xp, yp)
	    i = nint (yp)
	    yp = mod (i-1, SMW_NSPLIT) + 1
	    call mw_c2tran$t (smw_ct(ct,i), xp, yp, x2, y2)
	case SMW_PW:
	    i = nint (y1)
	    yp = mod (i-1, SMW_NSPLIT) + 1
	    call mw_c2tran$t (smw_ct(ct,i), x1, yp, x2, y2)
	case SMW_WL:
	    aps = SMW_APS(SMW_SMW(ct))
	    naps = SMW_NSPEC(SMW_SMW(ct))
	    j = nint (y1)
	    do i = 1, naps {
		if (Memi[aps+i-1] == j) {
		    call mw_c2tran$t (smw_ct(ct,i), x1, y1, xp, yp)
		    yp = i
		    call mw_c2tran$t (SMW_CTL(ct), xp, yp, x2, y2)
		    return
		}
	    }
	    call error (1, "Aperture not found")
	case SMW_WP:
	    aps = SMW_APS(SMW_SMW(ct))
	    naps = SMW_NSPEC(SMW_SMW(ct))
	    j = nint (y1)
	    do i = 1, naps {
		if (Memi[aps+i-1] == j) {
		    call mw_c2tran$t (smw_ct(ct,i), x1, y1, x2, y2)
		    y2 = i
		    return
		}
	    }
	    call error (1, "Aperture not found")
	default:
	    x2 = x1
	    y2 = y1
	}
end
$endfor