diff options
Diffstat (limited to 'sys/mwcs/wfzpn.x')
| -rw-r--r-- | sys/mwcs/wfzpn.x | 600 |
1 files changed, 600 insertions, 0 deletions
diff --git a/sys/mwcs/wfzpn.x b/sys/mwcs/wfzpn.x new file mode 100644 index 00000000..6c8db38a --- /dev/null +++ b/sys/mwcs/wfzpn.x @@ -0,0 +1,600 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math.h> +include <ctype.h> +include "imwcs.h" +include "mwcs.h" + +.help WFZPN +.nf ------------------------------------------------------------------------- +# WFZPN -- WCS function driver for the zenithal / azimuthal polynomial +# projection. + +Driver routines: + + FN_INIT wf_zpn_init (fc, dir) + FN_DESTROY wf_zpn_destroy (fc) + FN_FWD wf_zpn_fwd (fc, v1, v2) + FN_INV wf_zpn_inv (fc, v1, v2) + +.endhelp -------------------------------------------------------------------- + +define MAX_NITER 20 + +# Driver specific fields of function call (FC) descriptor. +define FC_LNGCOR Memi[$1+FCU] # RA axis correction +define FC_LATCOR Memi[$1+FCU+1] # DEC axis correction +define FC_IRA Memi[$1+FCU+2] # RA axis (1 or 2) +define FC_IDEC Memi[$1+FCU+3] # DEC axis (1 or 2) +define FC_NP Memd[P2D($1+FCU+4)] # poly order (0-9) +define FC_LONGP Memd[P2D($1+FCU+6)] # LONGPOLE (rads) +define FC_COLATP Memd[P2D($1+FCU+8)] # (90 - DEC) (rads) +define FC_COSLATP Memd[P2D($1+FCU+10)] # cosine (90 - DEC) +define FC_SINLATP Memd[P2D($1+FCU+12)] # sine (90 - DEC) +define FC_SPHTOL Memd[P2D($1+FCU+14)] # trig tolerance +define FC_PC Memd[P2D($1+FCU+16)+($2)] # poly coefficients (9) +define FC_RODEG Memd[P2D($1+FCU+36)] # RO (degs) +define FC_ZD Memd[P2D($1+FCU+38)] # colat of FIP (degs) +define FC_R Memd[P2D($1+FCU+40)] # radius of FIP (degs) +define FC_BADCVAL Memd[P2D($1+FCU+42)] # Bad coordinate value +define FC_W Memd[P2D($1+FCU+44)+($2)-1] # CRVAL (axis 1 and 2) + + +# WF_ZPN_INIT -- Initialize the zenithal/azimuthal polynomial forward or inverse +# transform. Initialization for this transformation consists of, determining +# which axis is RA / LON and which is DEC / LAT, computing the celestial +# longitude and colatitude of the native pole, reading in the the native +# longitude of the pole of the celestial coordinate system LONGPOLE from the +# attribute list, precomputing the Euler angles and various intermediary +# functions of the reference coordinates, reading in the projection parameter +# RO from the attribute list, reading in up to ten polynomial coefficients, +# and, for polynomial orders greater than 2 computing the colatitude and radius +# of the first point of inflection. If LONGPOLE is undefined then a value of +# 180.0 degrees is assumed. If RO is undefined a value of 180.0 / PI is +# assumed. If the polynomial coefficients are all zero then an error condition +# is posted. If the order of the polynomial is 2 or greater and there is no +# point of inflection an error condition is posted. The ZPN projection with +# an order of 1 and 0th and 1st coefficients of 0.0 and 1.0 respectively is +# equivalent to the ARC projtection. 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 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_zpn_init (fc, dir) + +pointer fc #I pointer to FC descriptor +int dir #I direction of transform + +int i, j, np, szatstr, maxorder, ualen, index, ip +double dec, zd1, d1, zd2, d2, zd, d, r, tol, dval +pointer sp, atname, atvalue, ct, mw, wp, wv, im, idb, rp +char compare[4] +bool match +int ctod(), strlen(), idb_nextcard(), itoc() +pointer wf_gsopen(), idb_open() +data tol/1.0d-13/ +errchk wf_decaxis(), mw_gwattrs() + +begin + # Allocate space for the attribute string. + call smark (sp) + call salloc (atname, SZ_ATNAME, TY_CHAR) + call salloc (atvalue, SZ_LINE, TY_CHAR) + + # Get the required mwcs pointers. + ct = FC_CT(fc) + mw = CT_MW(ct) + wp = FC_WCS(fc) + im = MI_REFIM(mw) + + # 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] + + # Get the celestial coordinates of the native pole which are in + # this case the ra and 90 - dec of the reference point. + + dec = DDEGTORAD(90.0d0 - FC_W(fc,FC_IDEC(fc))) + + # Determine the native longitude of the pole of the celestial + # coordinate system corresponding to the FITS keyword LONGPOLE. + # This number has no default and should normally be set to 180 + # degrees. Search both axes for this quantity. + + FC_LONGP(fc) = DDEGTORAD(180.0d0) + + # Precompute the trigomometric functions used by the spherical geometry + # code to improve efficiency. + + FC_COLATP(fc) = dec + FC_COSLATP(fc) = cos(dec) + FC_SINLATP(fc) = sin(dec) + + # 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. + + FC_RODEG(fc) = 180.0d0/DPI + szatstr = SZ_LINE + + # Fetch the longitude correction surface. Note that the attribute + # string may be of any length so the length of atvalue may have + # to be adjusted. + + FC_LNGCOR(fc) = NULL + + # Fetch the latitude correction surface. Note that the attribute + # string may be of any length so the length of atvalue may have + # to be adjusted. + + FC_LATCOR(fc) = NULL + + # Read through the fits header once more and pick up the PV matrix + # cards. Read the values and store them, keeping track of what is + # the highest order coefficient. With this projection only the dec + # axis coefficients matter. Technically we can have up to 99 + # coefficients. But we restrict this to 10 for the moment. + + maxorder = -1 + idb = idb_open(im,ualen) + compare[1] = 'P' + compare[2] = 'V' + i = itoc(FC_IDEC(fc),compare[3],1) + compare[4] = '_' + while (idb_nextcard(idb,rp) != EOF) { + match = true + do i = 0,3 { + if (Memc[rp+i] != compare[i+1]) { + match = false + break; + } + } + if (! match) + next + if (! IS_DIGIT(Memc[rp+4])) + next + index = TO_INTEG(Memc[rp+4]) + do i = 5,7 { + if (! IS_DIGIT(Memc[rp+i])) + break + else + index = 10*index + TO_INTEG(Memc[rp+i]) + } + if (index > 9) + next + ip = IDB_STARTVALUE + if (ctod(Memc[rp],ip,dval) <= 0) + dval = 0.0d0 + if (index > maxorder) + maxorder = index + FC_PC(fc,index) = dval + } + call idb_close(idb) + + # If all the coefficients are 0.0 the polynomial is undefined. + if (maxorder < 0) { + call sfree (sp) + call error (0, "WFT_ZPN_INIT: The polynomial is undefined") + } + + # Determine the number of coefficients. + FC_NP(fc) = double(maxorder) + np = maxorder + + if (np >= 3) { + # Find the point of inflection closest to the pole. + zd1 = 0.0d0 + d1 = FC_PC(fc,1) + if (d1 <= 0.0d0) { + call sfree (sp) + call error (0, + "WFT_ZPN_INIT: The point of inflection does not exist") + } + + # Find the point where the derivative first goes negative. + do i = 1, 180 { + zd2 = DPI * double (i) / 180.0d0 + d2 = 0.0d0 + do j = np, 1, -1 + d2 = d2 * zd2 + j * FC_PC(fc,j) + if (d2 <= 0.0d0) + break + zd1 = zd2 + d1 = d2 + } + + # Find where the derivative is 0. + if (d2 <= 0.0d0) { + do i = 1, 10 { + zd = zd1 - d1 * (zd2 - zd1) / (d2 - d1) + d = 0.0d0 + do j = np, 1, -1 + d = d * zd + j * FC_PC(fc,j) + if (abs(d) < tol) + break + if (d < 0.0d0) { + zd2 = zd + d2 = d + } else { + zd1 = zd + d1 = d + } + } + + # No negative derivative. + } else + zd = DPI + + r = 0.0d0 + do j = np, 0, -1 + r = r * zd + FC_PC(fc,j) + FC_ZD(fc) = zd + FC_R(fc) = r + } + + # Set the spherical trigonometric tolerance. + FC_SPHTOL(fc) = 1.0d-5 + + # Set the bad coordinate value. + FC_BADCVAL(fc) = INDEFD + + # Free working space. + call sfree (sp) + +end + + +# WF_ZPN_FWD -- Forward transform (physical to world) for the zenithal / +# azimuthal polynomial projection. + +procedure wf_zpn_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, i, j, k +double x, y, r, zd, a, b, c, d, zd1, zd2, r1, r2, lambda, rt, tol +double phi, theta, costhe, sinthe, dphi, cosphi, sinphi, ra, dec, dlng, z +double wf_gseval() +data tol/1.0d-13/ + +define phitheta_ 11 + +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. + + k = nint (FC_NP(fc)) + if (FC_LNGCOR(fc) == NULL) + x = p[ira] + else + x = p[ira] + wf_gseval (FC_LNGCOR(fc), p[ira], p[idec]) + if (FC_LATCOR(fc) == NULL) + y = p[idec] + else + y = p[idec] + wf_gseval (FC_LATCOR(fc), p[ira], p[idec]) + r = sqrt (x * x + y * y) / FC_RODEG(fc) + + # Solve. + + # Constant no solution + if (k < 1) { + w[ira] = FC_BADCVAL(fc) + w[idec] = FC_BADCVAL(fc) + return + + # Linear. + } else if (k == 1) { + zd = (r - FC_PC(fc,0)) / FC_PC(fc,1) + + # Quadratic. + } else if (k == 2) { + + a = FC_PC(fc,2) + b = FC_PC(fc,1) + c = FC_PC(fc,0) - r + d = b * b - 4.0d0 * a * c + if (d < 0.0d0) { + w[ira] = FC_BADCVAL(fc) + w[idec] = FC_BADCVAL(fc) + return + } + d = sqrt (d) + + # Choose solution closet to the pole. + zd1 = (-b + d) / (2.0d0 * a) + zd2 = (-b - d) / (2.0d0 * a) + zd = min (zd1, zd2) + if (zd < -tol) + zd = max (zd1, zd2) + if (zd < 0.0d0) { + if (zd < -tol) { + w[ira] = FC_BADCVAL(fc) + w[idec] = FC_BADCVAL(fc) + return + } + zd = 0.0d0 + } else if (zd > DPI) { + if (zd > (DPI + tol)) { + w[ira] = FC_BADCVAL(fc) + w[idec] = FC_BADCVAL(fc) + return + } + zd = DPI + } + + # Higher order solve iteratively. + } else { + + zd1 = 0.0d0 + r1 = FC_PC(fc,0) + zd2 = FC_ZD(fc) + r2 = FC_R(fc) + + if (r < r1) { + if (r < (r1 - tol)) { + w[ira] = FC_BADCVAL(fc) + w[idec] = FC_BADCVAL(fc) + return + } + zd = zd1 + goto phitheta_ + } else if (r > r2) { + if (r > (r2 + tol)) { + w[ira] = FC_BADCVAL(fc) + w[idec] = FC_BADCVAL(fc) + return + } + zd = zd2 + goto phitheta_ + } else { + do j = 1, 100 { + lambda = (r2 - r) / (r2 - r1) + if (lambda < 0.1d0) + lambda = 0.1d0 + else if (lambda > 0.9d0) + lambda = 0.9d0 + zd = zd2 - lambda * (zd2 - zd1) + rt = 0.0d0 + do i = k, 0, -1 + rt = (rt * zd) + FC_PC(fc,i) + if (rt < r) { + if ((r - rt) < tol) + goto phitheta_ + r1 = rt + zd1 = zd + } else { + if ((rt - r) < tol) + goto phitheta_ + r2 = rt + zd2 = zd + } + if (abs(zd2 - zd1) < tol) + goto phitheta_ + } + } + + } + +phitheta_ + + # Compute PHI. + if (r == 0.0d0) + phi = 0.0d0 + else + phi = atan2 (x, -y) + + # Compute THETA. + theta = DHALFPI - zd + + # 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_SINLATP(fc) - costhe * FC_COSLATP(fc) * cosphi + if (abs (x) < FC_SPHTOL(fc)) + x = -cos (theta + FC_COLATP(fc)) + costhe * FC_COSLATP(fc) * + (1.0d0 - cosphi) + y = -costhe * sinphi + if (x != 0.0d0 || y != 0.0d0) { + dlng = atan2 (y, x) + } else { + dlng = dphi + DPI + } + ra = FC_W(fc,ira) + DRADTODEG(dlng) + + # Normalize the RA. + if (FC_W(fc,ira) >= 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_COLATP(fc)) + if (dec > 90.0d0) + dec = 180.0d0 - dec + if (dec < -90.0d0) + dec = -180.0d0 - dec + } else { + z = sinthe * FC_COSLATP(fc) + costhe * FC_SINLATP(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_ZPN_INV -- Inverse transform (world to physical) for the zenithal / +# azimuthal polynomial projection. + +procedure wf_zpn_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, i, niter +double ra, dec, cosdec, sindec, cosra, sinra, x, y, phi, theta, s, r, dphi, z +double xm, ym, f, fx, fy, g, gx, gy, denom, dx, dy +double wf_gseval(), wf_gsder() + +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_W(fc,ira)) + dec = DDEGTORAD (w[idec]) + cosra = cos (ra) + sinra = sin (ra) + cosdec = cos (dec) + sindec = sin (dec) + + # Compute PHI. + x = sindec * FC_SINLATP(fc) - cosdec * FC_COSLATP(fc) * cosra + if (abs(x) < FC_SPHTOL(fc)) + x = -cos (dec + FC_COLATP(fc)) + cosdec * FC_COSLATP(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_COLATP(fc) + if (theta > DHALFPI) + theta = DPI - theta + if (theta < -DHALFPI) + theta = -DPI - theta + } else { + z = sindec * FC_COSLATP(fc) + cosdec * FC_SINLATP(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. + + s = DHALFPI - theta + r = 0.0d0 + do i = 9, 0, -1 + r = r * s + FC_PC(fc,i) + r = FC_RODEG(fc) * r + + if (FC_LNGCOR(fc) == NULL && FC_LATCOR(fc) == NULL) { + p[ira] = r * sin (phi) + p[idec] = -r * cos (phi) + + } else { + xm = r * sin (phi) + ym = -r * cos (phi) + x = xm + y = ym + niter = 0 + + repeat { + if (FC_LNGCOR(fc) != NULL) { + f = x + wf_gseval (FC_LNGCOR(fc), x, y) - xm + fx = wf_gsder (FC_LNGCOR(fc), x, y, 1, 0) + fx = 1.0 + fx + fy = wf_gsder (FC_LNGCOR(fc), x, y, 0, 1) + } else { + f = x - xm + fx = 1.0 + fy = 0.0 + } + if (FC_LATCOR(fc) != NULL) { + g = y + wf_gseval (FC_LATCOR(fc), x, y) - ym + gx = wf_gsder (FC_LATCOR(fc), x, y, 1, 0) + gy = wf_gsder (FC_LATCOR(fc), x, y, 0, 1) + gy = 1.0 + gy + } else { + g = y - ym + gx = 0.0 + gy = 1.0 + } + denom = fx * gy - fy * gx + if (denom == 0.0d0) + break + dx = (-f * gy + g * fy) / denom + dy = (-g * fx + f * gx) / denom + x = x + max (-1.0D0, min (1.0D0, dx)) + y = y + max (-1.0D0, min (1.0D0, dy)) + if (max (abs (dx), abs (dy), abs(f), abs(g)) < 2.80d-7) + break + + niter = niter + 1 + + } until (niter >= MAX_NITER) + + p[ira] = x + p[idec] = y + } +end + + +# WF_ZPN_DESTROY -- Free up the distortion surface pointers. + +procedure wf_zpn_destroy (fc) + +pointer fc #I pointer to the FC descriptor + +begin + if (FC_LNGCOR(fc) != NULL) + call wf_gsclose (FC_LNGCOR(fc)) + if (FC_LATCOR(fc) != NULL) + call wf_gsclose (FC_LATCOR(fc)) +end |