# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<math.h>
include	"mwcs.h"

.help WFPCO
.nf -------------------------------------------------------------------------
WFPCO -- WCS function driver for the polyconic projection.

Driver routines:

	FN_INIT		    wf_pco_init (fc, dir)
	FN_DESTROY			(none)
	FN_FWD		     wf_pco_fwd (fc, v1, v2)
	FN_INV		     wf_pco_inv (fc, v1, v2)

.endhelp --------------------------------------------------------------------

# Driver specific fields of function call (FC) descriptor.
define	FC_IRA		Memi[$1+FCU]		     # RA axis (1 or 2)
define	FC_IDEC		Memi[$1+FCU+1]		     # DEC axis (1 or 2)
define  FC_NATRA        Memd[P2D($1+FCU+2)]          # RA of native pole (rads)
define  FC_NATDEC       Memd[P2D($1+FCU+4)]          # DEC of native pole (rads)
define  FC_LONGP        Memd[P2D($1+FCU+6)]          # LONGPOLE (rads)
define  FC_COSDEC       Memd[P2D($1+FCU+8)]          # cosine (NATDEC)
define  FC_SINDEC       Memd[P2D($1+FCU+10)]         # sine (NATDEC)
define  FC_SPHTOL       Memd[P2D($1+FCU+12)]         # trig tolerance
define	FC_RODEG	Memd[P2D($1+FCU+14)]	     # RO (degs)
define	FC_RORAD	Memd[P2D($1+FCU+16)]	     # RO (rads)
define	FC_RECRORAD	Memd[P2D($1+FCU+18)]	     # 1 / RO (rads)
define	FC_2RODEG	Memd[P2D($1+FCU+20)]	     # 2 * RO 
define	FC_BADCVAL	Memd[P2D($1+FCU+22)]	     # bad coordinate value
define	FC_W		Memd[P2D($1+FCU+24)+($2)-1]  # CRVAL axis (1 and 2)


# WF_PCO_INIT -- Initialize the polyconic forward or inverse
# transform. Initialization for this transformation consists of, determining
# which axis is RA / LON and which is DEC / LAT, reading in the the native
# longitude and latitude  of the pole in celestial coordinates LONGPOLE and
# LATPOLE from the attribute list, computing the celestial longitude and
# colatitude of the native pole, precomputing the Euler anges and various
# intermediary functions of the reference point, and reading in the
# projection parameter RO from the attribute list. If LONGPOLE is undefined
# then a value of 180.0 degrees is assumed if the native latitude of the
# reference point is less than 0, otherwise 0 is assumed. If LATPOLE is
# undefined then the most northerly of the two possible solutions for the
# latitude of the native pole is chosen, otherwise the solution closest to
# LATPOLE is chosen. If RO is undefined a value of 180.0 / PI is assumed.
# In order to determine the axis order, the parameter "axtype={ra|dec}
# {xlon|xlat}" must have been set in the attribute list for the function.
# The LONGPOLE, LATPOLE and RO parameters may be set in either or both of the
# axes attribute lists, but the value in the RA axis attribute list takes
# precedence. 

procedure wf_pco_init (fc, dir)

pointer	fc			#I pointer to FC descriptor
int	dir			#I direction of transform

int	i
double  dec, latpole, theta0, clat0, slat0, cphip, sphip, cthe0, sthe0, x, y, z
double  u, v, latp1, latp2, latp, maxlat, tol
pointer	sp, atvalue, ct, mw, wp, wv
int	ctod()
data    tol/1.0d-10/
errchk	wf_decaxis(), mw_gwattrs()

begin
	# Allocate space for the attribute string.
	call smark (sp)
	call salloc (atvalue, SZ_LINE, TY_CHAR)

	# Get the required mwcs pointers.
	ct = FC_CT(fc)
	mw = CT_MW(ct)
	wp = FC_WCS(fc)

	# Determine which is the DEC axis, and hence the axis order.
	call wf_decaxis (fc, FC_IRA(fc), FC_IDEC(fc))

	# Get the value of W for each axis, i.e. the world coordinates at
	# the reference point.
	wv = MI_DBUF(mw) + WCS_W(wp) - 1
	do i = 1, 2
	    FC_W(fc,i) = Memd[wv+CT_AXIS(ct,FC_AXIS(fc,i))-1]

        # Determine the native longitude and latitude of the pole of the
        # celestial coordinate system corresponding to the FITS keywords
        # LONGPOLE and LATPOLE. LONGPOLE has no default but will be set
        # to 180 or 0 depending on the value of the declination of the
        # reference point. LATPOLE has no default but will be set depending
        # on the values of LONGPOLE and the reference declination.
        iferr {
            call mw_gwattrs (mw, FC_IRA(fc), "longpole", Memc[atvalue], SZ_LINE)
        } then {
            iferr {
                call mw_gwattrs (mw, FC_IDEC(fc), "longpole", Memc[atvalue],
                    SZ_LINE)
            } then {
                FC_LONGP(fc) = INDEFD
            } else {
                i = 1
                if (ctod (Memc[atvalue], i, FC_LONGP(fc)) <= 0)
                    FC_LONGP(fc) = INDEFD
            }
        } else {
            i = 1
            if (ctod (Memc[atvalue], i, FC_LONGP(fc)) <= 0)
                FC_LONGP(fc) = INDEFD
        }

       iferr {
            call mw_gwattrs (mw, FC_IRA(fc), "latpole", Memc[atvalue], SZ_LINE)
        } then {
            iferr {
                call mw_gwattrs (mw, FC_IDEC(fc), "latpole", Memc[atvalue],
                    SZ_LINE)
            } then {
                latpole = INDEFD
            } else {
                i = 1
                if (ctod (Memc[atvalue], i, latpole) <= 0)
                    latpole = INDEFD
            }
        } else {
            i = 1
            if (ctod (Memc[atvalue], i, latpole) <= 0)
                latpole = INDEFD
        }

       # Fetch the RO projection parameter which is the radius of the
        # generating sphere for the projection. If RO is absent which
        # is the usual case set it to 180 / PI. Search both axes for
        # this quantity.
        iferr {
            call mw_gwattrs (mw, FC_IRA(fc), "ro", Memc[atvalue], SZ_LINE)
        } then {
            iferr {
                call mw_gwattrs (mw, FC_IDEC(fc), "ro", Memc[atvalue],
                    SZ_LINE)
            } then {
                FC_RODEG(fc) = 180.0d0 / DPI
            } else {
                i = 1
                if (ctod (Memc[atvalue], i, FC_RODEG(fc)) <= 0)
                    FC_RODEG(fc) = 180.0d0 / DPI
            }
        } else {
            i = 1
            if (ctod (Memc[atvalue], i, FC_RODEG(fc)) <= 0)
                FC_RODEG(fc) = 180.0d0 / DPI
        }

        # Compute the native longitude of the celestial pole.
        dec = DDEGTORAD(FC_W(fc,FC_IDEC(fc)))
        theta0 = 0.0d0
        if (IS_INDEFD(FC_LONGP(fc))) {
            if (dec < theta0)
                FC_LONGP(fc) = DPI
            else
                FC_LONGP(fc) = 0.0d0
        } else
            FC_LONGP(fc) = DDEGTORAD(FC_LONGP(fc))

        # Compute the celestial longitude and latitude of the native pole.
        clat0 = cos (dec)
        slat0 = sin (dec)
        cphip = cos (FC_LONGP(fc))
        sphip = sin (FC_LONGP(fc))
        cthe0 = cos (theta0)
        sthe0 = sin (theta0)

        x = cthe0 * cphip
        y = sthe0
        z = sqrt (x * x + y * y)

        # The latitude of the native pole is determined by LATPOLE in this
        # case.
        if (z == 0.0d0) {

            if (slat0 != 0.0d0)
                call error (0, "WF_PCO_INIT: Invalid projection parameters")
            if (IS_INDEFD(latpole))
                latp = 999.0d0
            else
                latp = DDEGTORAD(latpole)

        } else {
            if (abs (slat0 / z) > 1.0d0)
                call error (0, "WF_PCO_INIT: Invalid projection parameters")


            u = atan2 (y, x)
            v = acos (slat0 / z)
            latp1 = u + v
            if (latp1 > DPI)
                latp1 = latp1 - DTWOPI
            else if (latp1 < -DPI)
                latp1 = latp1 + DTWOPI

            latp2 = u - v
            if (latp2 > DPI)
                latp2 = latp2 - DTWOPI
            else if (latp2 < -DPI)
                latp2 = latp2 + DTWOPI

            if (IS_INDEFD(latpole))
                maxlat = 999.0d0
            else
                maxlat = DDEGTORAD(latpole)
            if (abs (maxlat - latp1) < abs (maxlat - latp2)) {
                if (abs (latp1) < (DHALFPI + tol))
                    latp = latp1
                else
                    latp = latp2
            } else {
                if (abs (latp2) < (DHALFPI + tol))
                    latp = latp2
                else
                    latp = latp1
            }
        }

        FC_NATDEC(fc) = DHALFPI - latp

        z = cos (latp) * clat0
        if (abs(z) < tol) {

            # Celestial pole at the reference point.
            if (abs(clat0) < tol) {
                FC_NATRA(fc) = DDEGTORAD(FC_W(fc,FC_IRA(fc)))
                FC_NATDEC(fc) = DHALFPI - theta0
            # Celestial pole at the native north pole.
            } else if (latp > 0.0d0) {
                FC_NATRA(fc) = DDEGTORAD(FC_W(fc,FC_IRA(fc))) + FC_LONGP(fc) -
                    DPI
                FC_NATDEC(fc) = 0.0d0
            # Celestial pole at the native south pole.
            } else if (latp < 0.0d0) {
                FC_NATRA(fc) = DDEGTORAD(FC_W(fc,FC_IRA(fc))) - FC_LONGP(fc)
                FC_NATDEC(fc) = DPI
            }

        } else {
            x = (sthe0 - sin (latp) * slat0) / z
            y = sphip * cthe0 / clat0
            if (x == 0.0d0 && y == 0.0d0)
                call error (0, "WF_PCO_INIT: Invalid projection parameters")
            FC_NATRA(fc) = DDEGTORAD(FC_W(fc,FC_IRA(fc))) - atan2 (y,x)
        }

        if (FC_W(fc,FC_IRA(fc)) >= 0.0d0) {
            if (FC_NATRA(fc) < 0.0d0)
                FC_NATRA(fc) = FC_NATRA(fc) + DTWOPI
        } else {
            if (FC_NATRA(fc) > 0.0d0)
                FC_NATRA(fc) = FC_NATRA(fc) - DTWOPI
        }
        FC_COSDEC(fc) = cos (FC_NATDEC(fc))
        FC_SINDEC(fc) = sin (FC_NATDEC(fc))

        # Check for ill-conditioned parameters.
        if (abs(latp) > (DHALFPI+tol))
            call error (0, "WF_PCO_INIT: Invalid projection parameters")

	# Compute the required intermediate quantities.
	FC_RORAD(fc) = DDEGTORAD(FC_RODEG(fc))
	FC_RECRORAD(fc) = 1.0d0 / FC_RORAD(fc)
	FC_2RODEG(fc) = 2.0d0 * FC_RODEG(fc)

        # Set the bad coordinate value.
        FC_SPHTOL(fc) = 1.0d-5

	# Set the bad coordinate value.
	FC_BADCVAL(fc) = INDEFD

	# Free working space.
	call sfree (sp)
end


# WF_PCO_FWD -- Forward transform (physical to world) for the polyconic
# projection.

procedure wf_pco_fwd (fc, p, w)

pointer	fc			#I pointer to FC descriptor
double	p[2]			#I physical coordinates (x, y)
double	w[2]			#O world coordinates (ra, dec)

int	ira, idec, j
double	x, y, z, phi, theta, costhe, sinthe, dphi, cosphi, sinphi
double	ra, dec, wconst, tol, thepos, theneg, xx, ymthe, fpos, fneg, lambda
double	tanthe, f, dlng
double	xp, yp
data	tol / 1.0d-12/

begin
	# Get the axis numbers.
	ira = FC_IRA(fc)
	idec = FC_IDEC(fc)

	# Compute native spherical coordinates PHI and THETA in degrees from
	# the projected coordinates. This is the projection part of the
	# computation.

	x = p[ira]
	y = p[idec]

	# Compute PHI and THETA.
	wconst = abs (y * FC_RECRORAD(fc))
	if (wconst < tol) {
	    phi = x * FC_RECRORAD(fc)
	    theta = 0.0d0
	} else if (abs (wconst - 90.0d0) < tol) {
	    phi = 0.0d0
	    if (y >= 0.0d0)
		theta = DHALFPI
	    else
		theta = -DHALFPI
	} else {
	    if (y > 0.0d0)
		thepos = 90.0d0
	    else
		thepos = -90.0d0
	    theneg = 0.0d0

	    xx = x * x
	    ymthe = y - FC_RORAD(fc) * thepos
	    fpos = xx + ymthe * ymthe
	    fneg = -999.0d0

	    do j = 1, 64 {

		# Compute the required interval.
		if (fneg < -100.0d0)
		    theta = (thepos + theneg) / 2.0d0
		else {
		    lambda = fpos / (fpos - fneg)
		    if (lambda < 0.1d0)
			lambda = 0.1d0
		    else if (lambda > 0.9d0)
			lambda = 0.9d0
		    theta = thepos - lambda * (thepos - theneg)
		}

		# Compute the residue.
		ymthe = y - FC_RORAD(fc) * theta
		tanthe = tan (DDEGTORAD(theta))
		f = xx + ymthe * (ymthe - FC_2RODEG(fc) / tanthe)

		# Check for convergence.
		if (abs(f) < tol)
		    break
		if (abs (thepos - theneg) < tol)
		    break

		# Redefine the interval
		if (f > 0.0d0) {
		    thepos = theta
		    fpos = f
		} else {
		    theneg = theta
		    fneg = f
		}
	    }

	    theta = DDEGTORAD(theta)
	    xp = FC_RODEG(fc) - ymthe * tanthe
	    yp = x * tanthe
	    if (xp == 0.0d0 && yp == 0.0d0)
		phi = 0.0d0
	    else
		phi = atan2 (yp, xp) / sin (theta)
	}

	# Compute the celestial coordinates RA and DEC from the native
	# coordinates PHI and THETA. This is the spherical geometry part
	# of the computation.

	costhe = cos (theta)
	sinthe = sin (theta)
	dphi = phi - FC_LONGP(fc)
	cosphi = cos (dphi)
	sinphi = sin (dphi)

	# Compute the RA.
	x = sinthe * FC_SINDEC(fc) - costhe * FC_COSDEC(fc) * cosphi
        if (abs (x) < FC_SPHTOL(fc))
            x = -cos (theta + FC_NATDEC(fc)) + costhe * FC_COSDEC(fc) *
                (1.0d0 - cosphi)
	y = -costhe * sinphi
	if (x != 0.0d0 || y != 0.0d0) {
	    dlng = atan2 (y, x)
	} else {
	    dlng = dphi + DPI
	}
	ra = DRADTODEG(FC_NATRA(fc) + dlng)

	# Normalize the RA.
        if (FC_NATRA(fc) >= 0.0d0) {
            if (ra < 0.0d0)
                ra = ra + 360.0d0
        } else {
            if (ra > 0.0d0)
                ra = ra - 360.0d0
        }
	if (ra > 360.0d0)
	    ra = ra - 360.0d0
	else if (ra < -360.0d0)
	    ra = ra + 360.0d0

	# Compute the DEC.
	if (mod (dphi, DPI) == 0.0d0) {
            dec = DRADTODEG(theta + cosphi * FC_NATDEC(fc))
            if (dec > 90.0d0)
                dec = 180.0d0 - dec
            if (dec < -90.0d0)
                dec = -180.0d0 - dec
        } else {
            z = sinthe * FC_COSDEC(fc) + costhe * FC_SINDEC(fc) * cosphi
            if (abs(z) > 0.99d0) {
                if (z >= 0.0d0)
                    dec = DRADTODEG(acos (sqrt(x * x + y * y)))
                else
                    dec = DRADTODEG(-acos (sqrt(x * x + y * y)))
            } else
                dec = DRADTODEG(asin (z))
        }

	# Store the results.
	w[ira]  = ra
	w[idec] = dec
end


# WF_PCO_INV -- Inverse transform (world to physical) for the polyconic
# projection.

procedure wf_pco_inv (fc, w, p)

pointer	fc			#I pointer to FC descriptor
double	w[2]			#I input world (RA, DEC) coordinates
double	p[2]			#I output physical coordinates

int	ira, idec
double	ra, dec, cosdec, sindec, cosra, sinra, x, y, phi, theta, costhe
double	a, sinthe, cotthe, dphi, z

begin
	# Get the axes numbers.
	ira = FC_IRA(fc)
	idec = FC_IDEC(fc)

	# Compute the transformation from celestial coordinates RA and
	# DEC to native coordinates PHI and THETA. This is the spherical
	# geometry part of the transformation.

	ra  = DDEGTORAD (w[ira]) - FC_NATRA(fc)
	dec = DDEGTORAD (w[idec])
	cosra = cos (ra)
	sinra = sin (ra)
	cosdec = cos (dec)
	sindec = sin (dec)

	# Compute PHI.
	x = sindec * FC_SINDEC(fc) - cosdec * FC_COSDEC(fc) * cosra
        if (abs(x) < FC_SPHTOL(fc))
            x = -cos (dec + FC_NATDEC(fc)) + cosdec * FC_COSDEC(fc) *
                (1.0d0 - cosra)
	y = -cosdec * sinra
	if (x != 0.0d0 || y != 0.0d0)
	    dphi = atan2 (y, x) 
	else
	    dphi = ra - DPI
	phi = FC_LONGP(fc) + dphi
	if (phi > DPI)
	    phi = phi - DTWOPI
	else if (phi < -DPI)
	    phi = phi + DTWOPI

	# Compute THETA.
	if (mod (ra, DPI) == 0.0) {
            theta = dec + cosra * FC_NATDEC(fc)
            if (theta > DHALFPI)
                theta = DPI - theta
            if (theta < -DHALFPI)
                theta = -DPI - theta
        } else {
            z = sindec * FC_COSDEC(fc) + cosdec * FC_SINDEC(fc) * cosra
            if (abs (z) > 0.99d0) {
                if (z >= 0.0)
                    theta = acos (sqrt(x * x + y * y))
                else
                    theta = -acos (sqrt(x * x + y * y))
            } else
                theta = asin (z)
        }

	# Compute the transformation from native coordinates PHI and THETA
	# to projected coordinates X and Y.

	costhe = cos (theta)
	sinthe = sin (theta)
	a = phi * sinthe
	if (sinthe == 0.0d0) {
	    p[ira] = FC_RODEG(fc) * phi 
	    p[idec] = 0.0d0
	} else {
	    cotthe = costhe / sinthe
	    p[ira] = FC_RODEG(fc) * cotthe * sin (a)
	    p[idec] = FC_RODEG(fc) * (cotthe * (1.0d0 - cos(a)) + theta) 
	}

end