From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- sys/mwcs/wftpv.x | 556 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 556 insertions(+) create mode 100644 sys/mwcs/wftpv.x (limited to 'sys/mwcs/wftpv.x') diff --git a/sys/mwcs/wftpv.x b/sys/mwcs/wftpv.x new file mode 100644 index 00000000..812362be --- /dev/null +++ b/sys/mwcs/wftpv.x @@ -0,0 +1,556 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include +include "imwcs.h" +include "mwcs.h" + +.help WFTPV +.nf ------------------------------------------------------------------------- +# WFTPV -- WCS function driver for the TPV polynomial projection. + +Driver routines: + + FN_INIT wf_tpv_init (fc, dir) + FN_DESTROY wf_tpv_destroy (fc) + FN_FWD wf_tpv_fwd (fc, v1, v2) + FN_INV wf_tpv_inv (fc, v1, v2) + +.endhelp -------------------------------------------------------------------- + +define MAX_NITER 20 + +# 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_COSDEC Memd[P2D($1+FCU+2)] # cosine(dec) +define FC_SINDEC Memd[P2D($1+FCU+4)] # sine(dec) +define FC_W Memd[P2D($1+FCU+6)+($2)-1] # W (CRVAL) for each axis +define FC_NPRA Memi[$1+FCU+10] # poly order (0-39) +define FC_NPDEC Memi[$1+FCU+11] # poly order (0-39) +define FC_PV Memi[$1+FCU+12] # pointer to PV data (double) + +define FC_A Memd[FC_PV($1)+($2)] # RA coefficient +define FC_B Memd[FC_PV($1)+40+($2)] # DEC coefficient + +# WF_TPV_INIT -- Initialize the tan 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 TPV 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_tpv_init (fc, dir) + +pointer fc #I pointer to FC descriptor +int dir #I direction of transform + +int i, ualen, index, ip +double dec, dval +pointer ct, mw, wp, wv, im, idb, rp + +int idb_nextcard(), itoc(), ctod() +pointer idb_open() +errchk wf_decaxis() + +begin + # Allocate PV storage. This is freed in wf_tpv_destroy. + call calloc (FC_PV(fc), 80, TY_DOUBLE) + + # Set non-zero defaults. + FC_NPRA(fc) = 1 + FC_NPDEC(fc) = 1 + FC_A(fc,1) = 1D0 + FC_B(fc,1) = 1D0 + + # Get the required mwcs pointers. + + # 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. + + ct = FC_CT(fc) + mw = CT_MW(ct) + wp = FC_WCS(fc) + 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] + + # Precompute the sin and cos of the declination at the reference pixel. + dec = DEGTORAD(FC_W(fc,FC_IDEC(fc))) + FC_COSDEC(fc) = cos(dec) + FC_SINDEC(fc) = sin(dec) + + # Read through the fits header once more and pick up the PV cards. + # Read the values and store them, keeping track of what is + # the highest order coefficient. + + im = MI_REFIM(mw) + idb = idb_open(im,ualen) + while (idb_nextcard(idb,rp) != EOF) { + if (Memc[rp] != 'P' || Memc[rp+1] != 'V' || Memc[rp+3] != '_') + next + if (Memc[rp+2] != '1' && Memc[rp+2] != '2') + 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 > 39) + next + ip = IDB_STARTVALUE + if (ctod(Memc[rp],ip,dval) <= 0) + dval = 0.0d0 + i = TO_INTEG(Memc[rp+2]) + if (i == FC_IRA(fc)) { + FC_A(fc,index) = dval + if (index > FC_NPRA(fc)) + FC_NPRA(fc) = double(index) + } else { + FC_B(fc,index) = dval + if (index > FC_NPDEC(fc)) + FC_NPDEC(fc) = double(index) + } + } + call idb_close(idb) + +end + + +# WF_TPV_FWD -- Forward transform (physical to world) for the tangent plane +# with polynomial distortion. + +procedure wf_tpv_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) + +double x, y, z, a, b, ra, dec + +begin + # Compute the standard coordinates. + + x = p[1] + y = p[2] + call tpv_poly (fc, x, y, a, b) + + # Rotate the rectangular coordinate system of the vector [1,xi,eta] + # by the declination so that the X axis will pass through the equator. + + a = DEGTORAD(a) + b = DEGTORAD(b) + + x = FC_COSDEC(fc) - b * FC_SINDEC(fc) + y = a + z = FC_SINDEC(fc) + b * FC_COSDEC(fc) + + # Compute RA and DEC in radians. + if (x == 0.0D0 && y == 0.0D0) + ra = 0.0D0 + else + ra = atan2 (y, x) + dec = atan2 (z, sqrt (x*x + y*y)) + + # Return RA and DEC in degrees. + dec = RADTODEG(dec) + ra = RADTODEG(ra) + FC_W(fc,FC_IRA(fc)) + + if (ra < 0D0) + ra = ra + 360D0 + else if (ra > 360D0) + ra = ra - 360D0 + + w[FC_IRA(fc)] = ra + w[FC_IDEC(fc)] = dec + +end + + +# WF_TPV_INV -- Inverse transform (world to physical) for the tangent plane +# projection with polynomials. + +procedure wf_tpv_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, niter +double ra, dec +double cosra, cosdec, sinra, sindec, cosdist +double a, b, x, y, f, g, fx, gx, fy, gy, denom, dx, dy, dmax + +begin + ira = FC_IRA(fc) + idec = FC_IDEC(fc) + + ra = DEGTORAD (w[ira] - FC_W(fc,ira)) + dec = DEGTORAD (w[idec]) + + cosra = cos (ra) + sinra = sin (ra) + cosdec = cos (dec) + sindec = sin (dec) + cosdist = sindec * FC_SINDEC(fc) + cosdec * FC_COSDEC(fc) * cosra + + a = RADTODEG(cosdec * sinra / cosdist) + b = RADTODEG((sindec * FC_COSDEC(fc) - cosdec * FC_SINDEC(fc) * cosra) / + cosdist) + + x = a + y = b + dmax = 30. / 3600. + niter = 0 + + repeat { + call tpv_poly (fc, x, y, f, g) + call tpv_der (fc, x, y, fx, gx, fy, gy) + + f = f - a + g = g - b + + 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 (-dmax, min (dmax, dx)) + y = y + max (-dmax, min (dmax, dy)) + + if (max (abs (dx), abs (dy), abs(f), abs(g)) < 2.0d-7) + break + + niter = niter + 1 + + } until (niter >= MAX_NITER) + + p[ira] = x + p[idec] = y + +end + + +# WF_TPV_DESTROY -- Free up the distortion surface pointers. + +procedure wf_tpv_destroy (fc) + +pointer fc #I pointer to the FC descriptor + +begin + call mfree (FC_PV(fc), TY_DOUBLE) +end + + +# TPV_POLY -- Evaluate TPV polynomial (x,y -> xi,eta) + +procedure tpv_poly (fc, x, y, a, b) + +pointer fc #I pointer to FC descriptor +double x, y #I physical coordinates +double a, b #O standard coordinates (xi, eta) in deg + +int n +double r, r2, r3, r5, r7, x2, x3, x4, x5, x6, x7, y2, y3, y4, y5, y6, y7 + +begin + # Compute the standard coordinates. + # This depends on undefined coefficients being zero. + + x2 = x * x + y2 = y * y + r2 = x2 + y2 + r = sqrt (r2) + n = max (FC_NPRA(fc), FC_NPDEC(fc)) + + a = FC_A(fc,0) + FC_A(fc,1) * x + FC_A(fc,2) * y + FC_A(fc,3) * r + b = FC_B(fc,0) + FC_B(fc,1) * y + FC_B(fc,2) * x + FC_B(fc,3) * r + if (n <= 3) + return + a = a + FC_A(fc,4) * x2 + FC_A(fc,5) * x*y + FC_A(fc,6) * y2 + b = b + FC_B(fc,4) * y2 + FC_B(fc,5) * x*y + FC_B(fc,6) * x2 + if (n <= 6) + return + x3 = x * x2 + y3 = y * y2 + r3 = r * r2 + a = a + FC_A(fc,7) * x3 + b = b + FC_B(fc,7) * y3 + a = a + FC_A(fc,8) * x2*y + b = b + FC_B(fc,8) * y2*x + a = a + FC_A(fc,9) * x*y2 + b = b + FC_B(fc,9) * y*x2 + a = a + FC_A(fc,10) * y3 + b = b + FC_B(fc,10) * x3 + a = a + FC_A(fc,11) * r3 + b = b + FC_B(fc,11) * r3 + if (n <= 11) + return + x4 = x * x3 + y4 = y * y3 + a = a + FC_A(fc,12) * x4 + b = b + FC_B(fc,12) * y4 + a = a + FC_A(fc,13) * x3*y + b = b + FC_B(fc,13) * y3*x + a = a + FC_A(fc,14) * x2*y2 + b = b + FC_B(fc,14) * y2*x2 + a = a + FC_A(fc,15) * x*y3 + b = b + FC_B(fc,15) * y*x3 + a = a + FC_A(fc,16) * y4 + b = b + FC_B(fc,16) * x4 + if (n <= 16) + return + x5 = x * x4 + y5 = y * y4 + r5 = r3 * r2 + a = a + FC_A(fc,17) * x5 + b = b + FC_B(fc,17) * y5 + a = a + FC_A(fc,18) * x4*y + b = b + FC_B(fc,18) * y4*x + a = a + FC_A(fc,19) * x3*y2 + b = b + FC_B(fc,19) * y3*x2 + a = a + FC_A(fc,20) * x2*y3 + b = b + FC_B(fc,20) * y2*x3 + a = a + FC_A(fc,21) * x*y4 + b = b + FC_B(fc,21) * y*x4 + a = a + FC_A(fc,22) * y5 + b = b + FC_B(fc,22) * x5 + a = a + FC_A(fc,23) * r5 + b = b + FC_B(fc,23) * r5 + if (n <= 23) + return + x6 = x * x5 + y6 = y * y5 + a = a + FC_A(fc,14) * x6 + b = b + FC_B(fc,24) * y6 + a = a + FC_A(fc,25) * x5*y + b = b + FC_B(fc,25) * y5*x + a = a + FC_A(fc,26) * x4*y2 + b = b + FC_B(fc,26) * y4*x2 + a = a + FC_A(fc,27) * x3*y3 + b = b + FC_B(fc,27) * y3*x3 + a = a + FC_A(fc,28) * x2*y4 + b = b + FC_B(fc,28) * y2*x4 + a = a + FC_A(fc,29) * x*y5 + b = b + FC_B(fc,29) * y*x5 + a = a + FC_A(fc,30) * y6 + b = b + FC_B(fc,30) * x6 + if (n <= 30) + return + x7 = x * x6 + y7 = y * y6 + r7 = r5 * r2 + a = a + FC_A(fc,31) * x7 + b = b + FC_B(fc,31) * y7 + a = a + FC_A(fc,32) * x6*y + b = b + FC_B(fc,32) * y6*x + a = a + FC_A(fc,33) * x5*y2 + b = b + FC_B(fc,33) * y5*x2 + a = a + FC_A(fc,34) * x4*y3 + b = b + FC_B(fc,34) * y4*x3 + a = a + FC_A(fc,35) * x3*y4 + b = b + FC_B(fc,35) * y3*x4 + a = a + FC_A(fc,36) * x2*y5 + b = b + FC_B(fc,36) * y2*x5 + a = a + FC_A(fc,37) * x*y6 + b = b + FC_B(fc,37) * y*x6 + a = a + FC_A(fc,38) * y7 + b = b + FC_B(fc,38) * x7 + a = a + FC_A(fc,39) * r7 + b = b + FC_B(fc,39) * r7 + +end + + +# TPV_DER -- Evaluate TPV polynomial (x,y -> xi,eta) + +procedure tpv_der (fc, x, y, ax, bx, ay, by) + +pointer fc #I pointer to FC descriptor +double x, y #I physical coordinates +double ax, bx #O standard coordinates (xi, eta) in deg +double ay, by #O standard coordinates (xi, eta) in deg + +int n +double r, r2, r4, r6, x2, x3, x4, x5, x6, y2, y3, y4, y5, y6, rx, ry + +begin + x2 = x * x + y2 = y * y + r2 = x2 + y2 + r = sqrt (r2) + if (r < 2.0d-7) { + rx = 1D0 + ry = 1D0 + } else { + rx = x / r + ry = y / r + } + n = max (FC_NPRA(fc), FC_NPDEC(fc)) + + ax = FC_A(fc,1) + FC_A(fc,3) * rx + by = FC_B(fc,1) + FC_B(fc,3) * ry + ay = FC_A(fc,2) + FC_A(fc,3) * ry + bx = FC_B(fc,2) + FC_B(fc,3) * rx + if (n <= 3) + return + ax = ax + 2 * FC_A(fc,4) * x + FC_A(fc,5) * y + by = by + 2 * FC_B(fc,4) * y + FC_B(fc,5) * x + ay = ay + FC_A(fc,5) * x + 2 * FC_A(fc,6) * y + bx = bx + FC_B(fc,5) * y + 2 * FC_B(fc,6) * x + if (n <= 6) + return + ax = ax + 3 * FC_A(fc,7) * x2 + by = by + 3 * FC_B(fc,7) * y2 + ax = ax + 2 * FC_A(fc,8) * x*y + by = by + 2 * FC_B(fc,8) * y*x + ax = ax + FC_A(fc,9) * y2 + by = by + FC_B(fc,9) * x2 + ax = ax + 3 * FC_A(fc,11) * r2 * rx + by = by + 3 * FC_B(fc,11) * r2 * ry + ay = ay + FC_A(fc,8) * x2 + bx = bx + FC_B(fc,8) * y2 + ay = ay + 2 * FC_A(fc,9) * x*y + bx = bx + 2 * FC_B(fc,9) * y*x + ay = ay + 3 * FC_A(fc,10) * y2 + bx = bx + 3 * FC_B(fc,10) * x2 + ay = ay + 3 * FC_A(fc,11) * r2 * ry + bx = bx + 3 * FC_B(fc,11) * r2 * rx + if (n <= 11) + return + x3 = x * x2 + y3 = y * y2 + ax = ax + 4 * FC_A(fc,12) * x3 + by = by + 4 * FC_B(fc,12) * y3 + ax = ax + 3 * FC_A(fc,13) * x2*y + by = by + 3 * FC_B(fc,13) * y2*x + ax = ax + 2 * FC_A(fc,14) * x*y2 + by = by + 2 * FC_B(fc,14) * y*x2 + ax = ax + FC_A(fc,15) * y3 + by = by + FC_B(fc,15) * x3 + ay = ay + FC_A(fc,13) * x3 + bx = bx + FC_B(fc,13) * y3 + ay = ay + 2 * FC_A(fc,14) * x2*y + bx = bx + 2 * FC_B(fc,14) * y2*x + ay = ay + 3 * FC_A(fc,15) * x*y2 + bx = bx + 3 * FC_B(fc,15) * y*x2 + ay = ay + 4 * FC_A(fc,16) * y3 + bx = bx + 4 * FC_B(fc,16) * x3 + if (n <= 16) + return + x4 = x * x3 + y4 = y * y3 + r4 = r2 * r2 + ax = ax + 5 * FC_A(fc,17) * x4 + by = by + 5 * FC_B(fc,17) * y4 + ax = ax + 4 * FC_A(fc,18) * x3*y + by = by + 4 * FC_B(fc,18) * y3*x + ax = ax + 3 * FC_A(fc,19) * x2*y2 + by = by + 3 * FC_B(fc,19) * y2*x2 + ax = ax + 2 * FC_A(fc,20) * x*y3 + by = by + 2 * FC_B(fc,20) * y*x3 + ax = ax + FC_A(fc,21) * y4 + by = by + FC_B(fc,21) * x4 + ax = ax + 5 * FC_A(fc,23) * r4 * rx + by = by + 5 * FC_B(fc,23) * r4 * ry + ay = ay + FC_A(fc,18) * x4 + bx = bx + FC_B(fc,18) * y4 + ay = ay + 2 * FC_A(fc,19) * x3*y + bx = bx + 2 * FC_B(fc,19) * y3*x + ay = ay + 3 * FC_A(fc,20) * x2*y2 + bx = bx + 3 * FC_B(fc,20) * y2*x2 + ay = ay + 4 * FC_A(fc,21) * x*y3 + bx = bx + 4 * FC_B(fc,21) * y*x3 + ay = ay + 5 * FC_A(fc,22) * y4 + bx = bx + 5 * FC_B(fc,22) * x4 + ay = ay + 5 * FC_A(fc,23) * r4 * ry + bx = bx + 5 * FC_B(fc,23) * r4 * rx + if (n <= 23) + return + x5 = x * x4 + y5 = y * y4 + ax = ax + 6 * FC_A(fc,14) * x5 + by = by + 6 * FC_B(fc,24) * y5 + ax = ax + 5 * FC_A(fc,25) * x4*y + by = by + 5 * FC_B(fc,25) * y4*x + ax = ax + 4 * FC_A(fc,26) * x3*y2 + by = by + 4 * FC_B(fc,26) * y3*x2 + ax = ax + 3 * FC_A(fc,27) * x2*y3 + by = by + 3 * FC_B(fc,27) * y2*x3 + ax = ax + 2 * FC_A(fc,28) * x*y4 + by = by + 2 * FC_B(fc,28) * y*x4 + ax = ax + FC_A(fc,29) * y5 + by = by + FC_B(fc,29) * x5 + ay = ay + FC_A(fc,25) * x5 + bx = bx + FC_B(fc,25) * y5 + ay = ay + 2 * FC_A(fc,26) * x4*y + bx = bx + 2 * FC_B(fc,26) * y4*x + ay = ay + 3 * FC_A(fc,27) * x3*y2 + bx = bx + 3 * FC_B(fc,27) * y3*x2 + ay = ay + 4 * FC_A(fc,28) * x2*y3 + bx = bx + 4 * FC_B(fc,28) * y2*x3 + ay = ay + 5 * FC_A(fc,29) * x*y4 + bx = bx + 5 * FC_B(fc,29) * y*x4 + ay = ay + 6 * FC_A(fc,30) * y5 + bx = bx + 6 * FC_B(fc,30) * x5 + if (n <= 30) + return + x6 = x * x5 + y6 = y * y5 + r6 = r4 * r2 + ax = ax + 7 * FC_A(fc,31) * x6 + by = by + 7 * FC_B(fc,31) * y6 + ax = ax + 6 * FC_A(fc,32) * x5*y + by = by + 6 * FC_B(fc,32) * y5*x + ax = ax + 5 * FC_A(fc,33) * x4*y2 + by = by + 5 * FC_B(fc,33) * y4*x2 + ax = ax + 4 * FC_A(fc,34) * x3*y3 + by = by + 4 * FC_B(fc,34) * y3*x3 + ax = ax + 3 * FC_A(fc,35) * x2*y4 + by = by + 3 * FC_B(fc,35) * y2*x4 + ax = ax + 2 * FC_A(fc,36) * x*y5 + by = by + 2 * FC_B(fc,36) * y*x5 + ax = ax + FC_A(fc,37) * y6 + by = by + FC_B(fc,37) * x6 + ax = ax + 7 * FC_A(fc,39) * r6 * rx + by = by + 7 * FC_B(fc,39) * r6 * ry + ay = ay + FC_A(fc,32) * x6 + bx = bx + FC_B(fc,32) * y6 + ay = ay + 2 * FC_A(fc,33) * x5*y + bx = bx + 2 * FC_B(fc,33) * y5*x + ay = ay + 3 * FC_A(fc,34) * x4*y2 + bx = bx + 3 * FC_B(fc,34) * y4*x2 + ay = ay + 4 * FC_A(fc,35) * x3*y3 + bx = bx + 4 * FC_B(fc,35) * y3*x3 + ay = ay + 5 * FC_A(fc,36) * x2*y4 + bx = bx + 5 * FC_B(fc,36) * y2*x4 + ay = ay + 6 * FC_A(fc,37) * x*y5 + bx = bx + 6 * FC_B(fc,37) * y*x5 + ay = ay + 7 * FC_A(fc,38) * y6 + bx = bx + 7 * FC_B(fc,38) * x6 + ay = ay + 7 * FC_A(fc,39) * r6 * ry + bx = bx + 7 * FC_B(fc,39) * r6 * rx + +end -- cgit