aboutsummaryrefslogtreecommitdiff
path: root/sys/mwcs/wfpco.x
diff options
context:
space:
mode:
Diffstat (limited to 'sys/mwcs/wfpco.x')
-rw-r--r--sys/mwcs/wfpco.x518
1 files changed, 518 insertions, 0 deletions
diff --git a/sys/mwcs/wfpco.x b/sys/mwcs/wfpco.x
new file mode 100644
index 00000000..a9cd8e12
--- /dev/null
+++ b/sys/mwcs/wfpco.x
@@ -0,0 +1,518 @@
+# 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
generated by cgit (git 2.34.1) at 2025-09-20 05:11:28 -0400