# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. include include "mwcs.h" .help WFTSC .nf ------------------------------------------------------------------------- WFTSC -- WCS function driver for the tangentil spherical cube projection. Driver routines: FN_INIT wf_tsc_init (fc, dir) FN_DESTROY (none) FN_FWD wf_tsc_fwd (fc, v1, v2) FN_INV wf_tsc_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_C1 Memd[P2D($1+FCU+16)] # RO * (PI / 4) define FC_C2 Memd[P2D($1+FCU+18)] # (4 / PI) / RO define FC_BADCVAL Memd[P2D($1+FCU+20)] # bad coordinate value define FC_W Memd[P2D($1+FCU+22)+($2)-1] # CRVAL axis (1 and 2) # WF_TSC_INIT -- Initialize the forward or inverse tangential spherical cube # projection transform. Initialization for this transformation consists of, # determining which axis is RA / LON and which is DEC / LAT, reading in 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 angles and various # intermediary functions of the reference point, reading in the projection # parameter RO from the attribute list, and precomputing the various required # intermediate quantities. If LONGPOLE is undefined then a value of 180.0 # degrees is assumed if the celestial latitude is less than 0, otherwise 0 # is assumed. If LATPOLE is undefined then the most northerly of the two # possible solutions 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_tsc_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_TSC_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_TSC_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_TSC_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_TSC_INIT: Invalid projection parameters") # Compute the required intermediate quantities. FC_C1(fc) = FC_RODEG(fc) * (DPI / 4.0d0) FC_C2(fc) = 1.0d0 / FC_C1(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_TSC_FWD -- Forward transform (physical to world) for the tangential # spherical projection. procedure wf_tsc_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, face double l, m, n, phi, theta, costhe, sinthe, dphi, cosphi, sinphi, x, y, z double xf, yf, ra, dec, dlng 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. xf = p[ira] * FC_C2(fc) yf = p[idec] * FC_C2(fc) if (xf > 5.0d0) { face = 4 xf = xf - 6.0d0 l = -1.0d0 / sqrt (1.0d0 + xf * xf + yf * yf) m = -l * xf n = -l * yf } else if (xf > 3.0d0) { face = 3 xf = xf - 4.0d0 m = -1.0d0 / sqrt (1.0d0 + xf * xf + yf * yf) l = m * xf n = -m * yf } else if (xf > 1.0d0) { face = 2 xf = xf - 2.0d0 l = 1.0d0 / sqrt (1.0d0 + xf * xf + yf * yf) m = -l * xf n = l * yf } else if (yf > 1.0d0) { face = 0 yf = yf - 2.0d0 n = 1.0d0 / sqrt (1.0d0 + xf * xf + yf * yf) l = n * xf m = -n * yf } else if (yf < -1.0d0) { face = 5 yf = yf + 2.0d0 n = -1.0d0 / sqrt (1.0d0 + xf * xf + yf * yf) l = -n * xf m = -n * yf } else { face = 1 m = 1.0d0 / sqrt (1.0d0 + xf * xf + yf * yf) l = m * xf n = m * yf } # Compute PHI. if (l == 0.0d0 && m == 0.0d0) phi = 0.0d0 else phi = atan2 (l, m) # Compute THETA. theta = asin(n) # 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_TSC_INV -- Inverse transform (world to physical) for the tangential # spherical projection. procedure wf_tsc_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, face double ra, dec, cosdec, sindec, cosra, sinra, x, y, z, dphi, phi, theta double costhe, l, m, n, rho, tol, x0, y0, xf, yf data tol /1.0d-12/ 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) l = costhe * sin (phi) m = costhe * cos (phi) n = sin (theta) face = 0 rho = n if (m > rho) { face = 1 rho = m } if (l > rho) { face = 2 rho = l } if (-m > rho) { face = 3 rho = -m } if (-l > rho) { face = 4 rho = -l } if (-n > rho) { face = 5 rho = -n } switch (face) { case 0: xf = l / rho yf = -m / rho x0 = 0.0d0 y0 = 2.0d0 case 1: xf = l / rho yf = n / rho x0 = 0.0d0 y0 = 0.0d0 case 2: xf = -m / rho yf = n / rho x0 = 2.0d0 y0 = 0.0d0 case 3: xf = -l / rho yf = n / rho x0 = 4.0d0 y0 = 0.0d0 case 4: xf = m / rho yf = n / rho x0 = 6.0d0 y0 = 0.0d0 case 5: xf = l / rho yf = m / rho x0 = 0.0d0 y0 = -2.0d0 } if (abs(xf) > 1.0d0) { if (abs(xf) > (1.0d0 + tol)) { p[ira] = FC_BADCVAL(fc) p[idec] = FC_BADCVAL(fc) return } if (xf >= 0.0d0) xf = 1.0d0 else xf = -1.0d0 } if (abs(yf) > 1.0d0) { if (abs(yf) > (1.0d0 + tol)) { p[ira] = FC_BADCVAL(fc) p[idec] = FC_BADCVAL(fc) return } if (yf >= 0.0d0) yf = 1.0d0 else yf = -1.0d0 } p[ira] = FC_C1(fc) * (x0 + xf) p[idec] = FC_C1(fc) * (y0 + yf) end