include include include # RG_GETREFTIE -- Get the reference pixel coordinate tie points by reading # the image cursor or a file. int procedure rg_getreftie (fd, xreftie, yreftie, ntie, file_type, interactive) int fd #I the input file descriptor real xreftie[ARB] #O the output x coordinates of the tie points real yreftie[ARB] #O the output y coordinates of the tie points int ntie #I the number of tie points int file_type #I the input file type bool interactive #I the int nref, wcs, key pointer sp, str int clgcur(), fscan(), nscan() begin # Allocate temporary space. call smark (sp) call salloc (str, SZ_LINE, TY_CHAR) # Print the prompt string. if (interactive) { # Issue prompt. if (file_type == RG_REFFILE) { call printf ( "\nMark 1-%d reference objects on the display\n") } else { call printf ( "\nMark the same %d input objects on the display\n") } call pargi (ntie) # Mark the points nref = 0 while (clgcur ("icommands", xreftie[nref+1], yreftie[nref+1], wcs, key, Memc[str], SZ_LINE) != EOF) { nref = nref + 1 if (file_type == RG_REFFILE) { call printf (" Reference coordinate %d %0.3f %0.3f\n") call pargi (nref) call pargr (xreftie[nref]) call pargr (yreftie[nref]) } else { call printf (" Input coordinate %d %0.3f %0.3f\n") call pargi (nref) call pargr (xreftie[nref]) call pargr (yreftie[nref]) } if (nref >= ntie) break } } else { # Issue prompt. if (fd == STDIN) { if (file_type == RG_REFFILE) { call printf ( "\nEnter coordinates of 1-%d reference objects\n") } else { call printf ( "Enter coordinates of %d corresponding input objects\n") } call pargi (ntie) } nref = 0 while (fscan (fd) != EOF) { call gargr (xreftie[nref+1]) call gargr (yreftie[nref+1]) if (nscan() < 2) break nref = nref + 1 if (nref >= ntie) break call gargr (xreftie[nref+1]) call gargr (yreftie[nref+1]) if (nscan() < 4) break nref = nref + 1 if (nref >= ntie) break call gargr (xreftie[nref+1]) call gargr (yreftie[nref+1]) if (nscan() < 6) break nref = nref + 1 break } } call sfree (sp) return (nref) end # RG_GETREFCEL -- Get the reference pixel coordinate tie points by reading # the image cursor or a file. int procedure rg_getrefcel (fd, xreftie, yreftie, ntie, projection, reflng, reflat, lngunits, latunits, file_type) int fd #I the input file descriptor real xreftie[ARB] #O the output x coordinates of the tie points real yreftie[ARB] #O the output y coordinates of the tie points int ntie #I the number of tie points char projection[ARB] #I the sky projection geometry double reflng #I the ra / longitude of the reference point double reflat #I the dec / latitude of the reference point int lngunits #I the ra / longitude units int latunits #I the dec / latitude units int file_type #I the input file type int nref pointer sp, dxref, dyref, str int fscan(), nscan() begin # Allocate temporary space. call smark (sp) call salloc (str, SZ_LINE, TY_CHAR) call salloc (dxref, ntie, TY_DOUBLE) call salloc (dyref, ntie, TY_DOUBLE) # Issue prompt. if (fd == STDIN) { if (file_type == RG_REFFILE) { call printf ( "\nEnter coordinates of 1-%d reference objects\n") } else { call printf ( "Enter coordinates of %d corresponding input objects\n") } call pargi (ntie) } # Read in the tie point. nref = 0 while (fscan (fd) != EOF) { call gargd (Memd[dxref+nref]) call gargd (Memd[dyref+nref]) if (nscan() < 2) break nref = nref + 1 if (nref >= ntie) break call gargd (Memd[dxref+nref]) call gargd (Memd[dyref+nref]) if (nscan() < 4) break nref = nref + 1 if (nref >= ntie) break call gargd (Memd[dxref+nref]) call gargd (Memd[dyref+nref]) if (nscan() < 6) break nref = nref + 1 break } # Convert to standard coordinates. if (nref > 0) { call rg_celtostd (projection, Memd[dxref], Memd[dyref], Memd[dxref], Memd[dyref], nref, reflng, reflat, lngunits, latunits) call amulkd (Memd[dxref], 3600.0d0, Memd[dxref], nref) call amulkd (Memd[dyref], 3600.0d0, Memd[dyref], nref) call achtdr (Memd[dxref], xreftie, nref) call achtdr (Memd[dyref], yreftie, nref) } call sfree (sp) return (nref) end # RG_PLINCOEFF -- Print the computed transformation on the standard output. procedure rg_plincoeff (xlabel, ylabel, xref, yref, xlist, ylist, ntie, coeff, ncoeff) char xlabel[ARB] #I the x equation label char ylabel[ARB] #I the x equation label real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates real xlist[ARB] #I the input x input coordinates real ylist[ARB] #I the input y input coordinates int ntie #I number of tie points real coeff[ARB] #I the output coefficient array int ncoeff #I the number of coefficients int i real xmag, ymag, xrot, yrot begin # List the tie points on the standard output. if (ntie > 0) { do i = 1, ntie { call printf ( " tie point: %3d ref: %9.3f %9.3f input: %9.3f %9.3f\n") call pargi (i) call pargr (xref[i]) call pargr (yref[i]) call pargr (xlist[i]) call pargr (ylist[i]) } call printf ("\n") } # Print the transformation coefficients to the standard output. call printf ("Initial linear transformation\n") call printf (" %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (xlabel) call pargr (coeff[3]) call pargr (coeff[1]) call pargr (coeff[2]) call printf (" %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (ylabel) call pargr (coeff[6]) call pargr (coeff[4]) call pargr (coeff[5]) call rg_ctogeo (coeff[1], -coeff[2], -coeff[4], coeff[5], xmag, ymag, xrot, yrot) call printf ( " dx: %0.2f dy: %0.2f xmag: %0.3f ymag: %0.3f xrot: %0.1f yrot: %0.1f\n") call pargr (coeff[3]) call pargr (coeff[6]) call pargr (xmag) call pargr (ymag) call pargr (xrot) call pargr (yrot) call printf ("\n") end # RG_PMLINCOEFF -- Print the computed transformation on the standard output. procedure rg_pmlincoeff (xlabel, ylabel, coeff, ncoeff) char xlabel[ARB] #I the x equation label char ylabel[ARB] #I the x equation label real coeff[ARB] #I the output coefficient array int ncoeff #I the number of coefficients real xmag, ymag, xrot, yrot begin # Write the matched transformation coefficients to the standard output. call printf ("Matched triangles transformation\n") call printf (" %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (xlabel) call pargr (coeff[3]) call pargr (coeff[1]) call pargr (coeff[2]) call printf (" %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (ylabel) call pargr (coeff[6]) call pargr (coeff[4]) call pargr (coeff[5]) call rg_ctogeo (coeff[1], -coeff[2], -coeff[4], coeff[5], xmag, ymag, xrot, yrot) call printf ( " dx: %0.2f dy: %0.2f xmag: %0.3f ymag: %0.3f xrot: %0.1f yrot: %0.1f\n") call pargr (coeff[3]) call pargr (coeff[6]) call pargr (xmag) call pargr (ymag) call pargr (xrot) call pargr (yrot) call printf ("\n") end # RG_WLINCOEFF -- Write the computed transformation to the output file. procedure rg_wlincoeff (fd, xlabel, ylabel, xref, yref, xlist, ylist, ntie, coeff, ncoeff) int fd #I pointer to the output file char xlabel[ARB] #I the x equation label char ylabel[ARB] #I the x equation label real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates real xlist[ARB] #I the input x input coordinates real ylist[ARB] #I the input y input coordinates int ntie #I number of tie points real coeff[ARB] #I the output coefficient array int ncoeff #I the number of coefficients int i real xmag, ymag, xrot, yrot begin # List the tie points. if (ntie > 0) { do i = 1, ntie { call fprintf (fd, "# tie point: %3d ref: %9.3f %9.3f input: %9.3f %9.3f\n") call pargi (i) call pargr (xref[i]) call pargr (yref[i]) call pargr (xlist[i]) call pargr (ylist[i]) } call fprintf (fd, "#\n") } # Write the transformation coefficients to the output file. call fprintf (fd, "# Initial linear transformation\n") call fprintf (fd, "# %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (xlabel) call pargr (coeff[3]) call pargr (coeff[1]) call pargr (coeff[2]) call fprintf (fd, "# %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (ylabel) call pargr (coeff[6]) call pargr (coeff[4]) call pargr (coeff[5]) call rg_ctogeo (coeff[1], -coeff[2], -coeff[4], coeff[5], xmag, ymag, xrot, yrot) call fprintf (fd, "# dx: %0.2f dy: %0.2f xmag: %0.3f ymag: %0.3f xrot: %0.1f yrot: %0.1f\n") call pargr (coeff[3]) call pargr (coeff[6]) call pargr (xmag) call pargr (ymag) call pargr (xrot) call pargr (yrot) call fprintf (fd, "#\n") end # RG_WMLINCOEFF -- Print the computed transformation on the standard output. procedure rg_wmlincoeff (ofd, xlabel, ylabel, coeff, ncoeff) int ofd #I the output file descriptor char xlabel[ARB] #I the x equation label char ylabel[ARB] #I the x equation label real coeff[ARB] #I the output coefficient array int ncoeff #I the number of coefficients real xmag, ymag, xrot, yrot begin # Write the matched transformation coefficients to the output file. call fprintf (ofd, "# Matched triangles transformation\n") call fprintf (ofd, "# %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (xlabel) call pargr (coeff[3]) call pargr (coeff[1]) call pargr (coeff[2]) call fprintf (ofd, "# %4.4s[tie] = %10g + %10g * x[tie] + %10g * y[tie]\n") call pargstr (ylabel) call pargr (coeff[6]) call pargr (coeff[4]) call pargr (coeff[5]) call rg_ctogeo (coeff[1], -coeff[2], -coeff[4], coeff[5], xmag, ymag, xrot, yrot) call fprintf (ofd, "# dx: %0.2f dy: %0.2f xmag: %0.3f ymag: %0.3f xrot: %0.1f yrot: %0.1f\n") call pargr (coeff[3]) call pargr (coeff[6]) call pargr (xmag) call pargr (ymag) call pargr (xrot) call pargr (yrot) call fprintf (ofd, "#\n") end # RG_LINCOEFF -- Compute the transformation given one to three tie points. int procedure rg_lincoeff (xref, yref, xlist, ylist, ntie, coeff, ncoeff) real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates real xlist[ARB] #I the input x input coordinates real ylist[ARB] #I the input y input coordinates int ntie #I number of tie points real coeff[ARB] #O the output coefficient array int ncoeff #I the number of coefficients int ier, xier, yier, nfcoeff pointer sp, wts, fcoeff, sx, sy real xmin, xmax, ymin, ymax int rg_onestar(), rg_twostar(), rg_threestar() begin switch (ntie) { case 0: ier = ERR case 1: ier = rg_onestar (xref, yref, xlist, ylist, ntie, coeff, ncoeff) case 2: ier = rg_twostar (xref, yref, xlist, ylist, ntie, coeff, ncoeff) case 3: ier = rg_threestar (xref, yref, xlist, ylist, ntie, coeff, ncoeff) default: call smark (sp) call salloc (fcoeff, 3, TY_REAL) call salloc (wts, ntie, TY_REAL) call alimr (xlist, ntie, xmin, xmax) call alimr (ylist, ntie, ymin, ymax) call gsinit (sx, GS_POLYNOMIAL, 2, 2, GS_XNONE, xmin, xmax, ymin, ymax) call gsinit (sy, GS_POLYNOMIAL, 2, 2, GS_XNONE, xmin, xmax, ymin, ymax) call amovkr (1.0, Memr[wts], ntie) call gsfit (sx, xlist, ylist, xref, Memr[wts], ntie, WTS_UNIFORM, xier) call gsfit (sy, xlist, ylist, yref, Memr[wts], ntie, WTS_UNIFORM, yier) if (xier == OK && xier == OK) { call gscoeff (sx, Memr[fcoeff], nfcoeff) coeff[3] = Memr[fcoeff] coeff[1] = Memr[fcoeff+1] coeff[2] = Memr[fcoeff+2] call gscoeff (sy, Memr[fcoeff], nfcoeff) coeff[6] = Memr[fcoeff] coeff[4] = Memr[fcoeff+1] coeff[5] = Memr[fcoeff+2] ier = OK } else ier = ERR call gsfree (sx) call gsfree (sy) call sfree (sp) } return (ier) end # RG_COMPUTE -- Transform the input list coordinates. The transformation # may be done in place. procedure rg_compute (xlist, ylist, xtrans, ytrans, nstars, coeff, ncoeff) real xlist[ARB] #I the input x coordinates real ylist[ARB] #I the input y coordinates real xtrans[ARB] #O the output x transformed coordinates real ytrans[ARB] #O the output y transformed coordinates int nstars #I the number of points real coeff[ARB] #I the input coefficient array int ncoeff #I the number of coefficients int i real xval, yval begin do i = 1, nstars { xval = xlist[i] yval = ylist[i] xtrans[i] = coeff[1] * xval + coeff[2] * yval + coeff[3] ytrans[i] = coeff[4] * xval + coeff[5] * yval + coeff[6] } end # RG_INTERSECT -- Compute the intersection of two sorted lists given a # matching tolerance. int procedure rg_intersect (ofd, xref, yref, refindex, rlineno, nrefstars, xlist, ylist, xtrans, ytrans, listindex, ilineno, nliststars, tolerance, xformat, yformat) int ofd #I the output file descriptor real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates int refindex[ARB] #I the input reference coordinates sort index int rlineno[ARB] #I the input reference coordinate line numbers int nrefstars #I the number of reference stars real xlist[ARB] #I the input x list coordinates real ylist[ARB] #I the input y list coordinates real xtrans[ARB] #I the input x transformed list coordinates real ytrans[ARB] #I the input y transformed list coordinates int listindex[ARB] #I the input list sort index int ilineno[ARB] #I the input input line numbers int nliststars #I the number of input stars real tolerance #I the matching tolerance char xformat[ARB] #I the output x coordinate format char yformat[ARB] #I the output y coordinate format int blp, rp, rindex, lp, lindex, rmatch, lmatch, ninter pointer sp, fmtstr real dx, dy, tol2, rmax2, r2 begin call smark (sp) call salloc (fmtstr, SZ_LINE, TY_CHAR) # Construct the fromat string call sprintf (Memc[fmtstr], SZ_LINE, "%s %s %s %s %%5d %%5d\n") if (xformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (xformat) if (yformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (yformat) if (xformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (xformat) if (yformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (yformat) # Initialize the intersection routine. tol2 = tolerance ** 2 blp = 1 ninter = 0 # Loop over the reference list stars. for (rp = 1; rp <= nrefstars; rp = rp + 1) { # Get the index of the reference star in question. rindex = refindex[rp] # Compute the start of the search range. for (; blp <= nliststars; blp = blp + 1) { lindex = listindex[blp] dy = yref[rindex] - ytrans[lindex] if (dy < tolerance) break } # Break if the end of the input list is reached. if (blp > nliststars) break # If one is outside the tolerance limits skip to next reference # object. if (dy < -tolerance) next # Find the closest match to the reference object. rmax2 = tol2 rmatch = 0 lmatch = 0 for (lp = blp; lp <= nliststars; lp = lp + 1) { # Compute the distance between the two points. lindex = listindex[lp] dy = yref[rindex] - ytrans[lindex] if (dy < -tolerance) break dx = xref[rindex] - xtrans[lindex] r2 = dx ** 2 + dy ** 2 # A match has been found. if (r2 <= rmax2) { rmax2 = r2 rmatch = rindex lmatch = lindex } } # A match was found so write the results to the output file. if (rmatch > 0 && lmatch > 0) { ninter = ninter + 1 call fprintf (ofd, Memc[fmtstr]) call pargr (xref[rmatch]) call pargr (yref[rmatch]) call pargr (xlist[lmatch]) call pargr (ylist[lmatch]) call pargi (rlineno[rmatch]) call pargi (ilineno[lmatch]) } } call sfree (sp) return (ninter) end # RG_LLINTERSECT -- Compute the intersection of two sorted lists given a # matching tolerance. int procedure rg_llintersect (ofd, lngref, latref, xref, yref, refindex, rlineno, nrefstars, xlist, ylist, xtrans, ytrans, listindex, ilineno, nliststars, tolerance, lngformat, latformat, xformat, yformat) int ofd #I the output file descriptor double lngref[ARB] #I the input ra/longitude reference coordinates double latref[ARB] #I the input dec/latitude reference coordinates real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates int refindex[ARB] #I the input reference coordinates sort index int rlineno[ARB] #I the input reference coordinate line numbers int nrefstars #I the number of reference stars real xlist[ARB] #I the input x list coordinates real ylist[ARB] #I the input y list coordinates real xtrans[ARB] #I the input x transformed list coordinates real ytrans[ARB] #I the input y transformed list coordinates int listindex[ARB] #I the input list sort index int ilineno[ARB] #I the input input line numbers int nliststars #I the number of input stars real tolerance #I the matching tolerance char lngformat[ARB] #I the output ra / longitude coordinate format char latformat[ARB] #I the output dec / latitude coordinate format char xformat[ARB] #I the output x coordinate format char yformat[ARB] #I the output y coordinate format int blp, rp, rindex, lp, lindex, rmatch, lmatch, ninter pointer sp, fmtstr real dx, dy, tol2, rmax2, r2 begin call smark (sp) call salloc (fmtstr, SZ_LINE, TY_CHAR) # Construct the fromat string call sprintf (Memc[fmtstr], SZ_LINE, "%s %s %s %s %%5d %%5d\n") if (lngformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (lngformat) if (latformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (latformat) if (xformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (xformat) if (yformat[1] == EOS) call pargstr ("%13.7g") else call pargstr (yformat) # Initialize the intersection routine. tol2 = tolerance ** 2 blp = 1 ninter = 0 # Loop over the reference list stars. for (rp = 1; rp <= nrefstars; rp = rp + 1) { # Get the index of the reference star in question. rindex = refindex[rp] # Compute the start of the search range. for (; blp <= nliststars; blp = blp + 1) { lindex = listindex[blp] dy = yref[rindex] - ytrans[lindex] if (dy < tolerance) break } # Break if the end of the input list is reached. if (blp > nliststars) break # If one is outside the tolerance limits skip to next reference # object. if (dy < -tolerance) next # Find the closest match to the reference object. rmax2 = tol2 rmatch = 0 lmatch = 0 for (lp = blp; lp <= nliststars; lp = lp + 1) { # Compute the distance between the two points. lindex = listindex[lp] dy = yref[rindex] - ytrans[lindex] if (dy < -tolerance) break dx = xref[rindex] - xtrans[lindex] r2 = dx ** 2 + dy ** 2 # A match has been found. if (r2 <= rmax2) { rmax2 = r2 rmatch = rindex lmatch = lindex } } # A match was found so write the results to the output file. if (rmatch > 0 && lmatch > 0) { ninter = ninter + 1 call fprintf (ofd, Memc[fmtstr]) call pargd (lngref[rmatch]) call pargd (latref[rmatch]) call pargr (xlist[lmatch]) call pargr (ylist[lmatch]) call pargi (rlineno[rmatch]) call pargi (ilineno[lmatch]) } } call sfree (sp) return (ninter) end # RG_LMKCOEFF -- Given the geometry of a linear transformation compute # the coefficients required to tranform from the input to the reference # system. procedure rg_lmkcoeff (xin, yin, xmag, ymag, xrot, yrot, xout, yout, coeff, ncoeff) real xin, yin #I the origin of the input coordinates real xmag, ymag #I the input x and y scale factors real xrot, yrot #I the iput x and y rotation factors real xout, yout #I the origin of the reference coordinates real coeff[ARB] #O the output coefficient array int ncoeff #I the number of coefficients begin # Compute the x fit coefficients. coeff[1] = xmag * cos (DEGTORAD(xrot)) coeff[2] = -ymag * sin (DEGTORAD(yrot)) coeff[3] = xout - coeff[1] * xin - coeff[2] * yin # Compute the y fit coefficients. coeff[4] = xmag * sin (DEGTORAD(xrot)) coeff[5] = ymag * cos (DEGTORAD(yrot)) coeff[6] = yout - coeff[4] * xin - coeff[5] * yin end # RG_ONESTAR -- Compute the transformation coefficients for a simple # shift operation. int procedure rg_onestar (xref, yref, xlist, ylist, ntie, coeff, ncoeff) real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates real xlist[ARB] #I the input x list coordinates real ylist[ARB] #I the input y list coordinates int ntie #I the number of tie points real coeff[ARB] #O the output coefficient array int ncoeff #I the number of coefficients begin # Compute the x transformation. coeff[1] = 1.0 coeff[2] = 0.0 coeff[3] = xref[1] - xlist[1] # Compute the y transformation. coeff[4] = 0.0 coeff[5] = 1.0 coeff[6] = yref[1] - ylist[1] return (OK) end # RG_TWOSTAR -- Compute the transformation coefficients of a simple shift, # magnification, and rotation. int procedure rg_twostar (xref, yref, xlist, ylist, ntie, coeff, ncoeff) real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates real xlist[ARB] #I the input x list coordinates real ylist[ARB] #I the input y list coordinates int ntie #I the number of tie points real coeff[ARB] #O the output coefficient array int ncoeff #I the number of coefficients real rot, mag, dxlis, dylis, dxref, dyref, cosrot, sinrot real rg_posangle() begin # Compute the deltas. dxlis = xlist[2] - xlist[1] dylis = ylist[2] - ylist[1] dxref = xref[2] - xref[1] dyref = yref[2] - yref[1] # Compute the rotation angle. rot = rg_posangle (dxref, dyref) - rg_posangle (dxlis, dylis) cosrot = cos (rot) sinrot = sin (rot) # Compute the magnification factor. mag = dxlis ** 2 + dylis ** 2 if (mag <= 0.0) mag = 0.0 else mag = sqrt ((dxref ** 2 + dyref ** 2) / mag) # Compute the transformation coefficicents. coeff[1] = mag * cosrot coeff[2] = - mag * sinrot coeff[3] = xref[1] - mag * cosrot * xlist[1] + mag * sinrot * ylist[1] coeff[4] = mag * sinrot coeff[5] = mag * cosrot coeff[6] = yref[1] - mag * sinrot * xlist[1] - mag * cosrot * ylist[1] return (OK) end # RG_THREESTAR -- Compute the transformation coefficients using a simple # shift, magnification in x and y, rotation, and skew. int procedure rg_threestar (xref, yref, xlist, ylist, ntie, coeff, ncoeff) real xref[ARB] #I the input x reference coordinates real yref[ARB] #I the input y reference coordinates real xlist[ARB] #I the input x list coordinates real ylist[ARB] #I the input y list coordinates int ntie #I the number of tie points real coeff[ARB] #O the output coefficient array int ncoeff #I the number of coefficients real dx23, dx13, dx12, dy23, dy13, dy12, det bool fp_equalr() int rg_twostar() begin # Compute the deltas. dx23 = xlist[2] - xlist[3] dx13 = xlist[1] - xlist[3] dx12 = xlist[1] - xlist[2] dy23 = ylist[2] - ylist[3] dy13 = ylist[1] - ylist[3] dy12 = ylist[1] - ylist[2] # Compute the determinant. det = xlist[1] * dy23 - xlist[2] * dy13 + xlist[3] * dy12 if (fp_equalr (det, 0.0)) return (rg_twostar (xref, yref, xlist, ylist, ntie, coeff, ncoeff)) # Compute the x transformation. coeff[1] = (xref[1] * dy23 - xref[2] * dy13 + xref[3] * dy12) / det coeff[2] = (-xref[1] * dx23 + xref[2] * dx13 - xref[3] * dx12) / det coeff[3] = (xref[1] * (xlist[2] * ylist[3] - xlist[3] * ylist[2]) + xref[2] * (ylist[1] * xlist[3] - xlist[1] * ylist[3]) + xref[3] * (xlist[1] * ylist[2] - ylist[1] * xlist[2])) / det # Compute the y transformation. coeff[4] = (yref[1] * dy23 - yref[2] * dy13 + yref[3] * dy12) / det coeff[5] = (-yref[1] * dx23 + yref[2] * dx13 - yref[3] * dx12) / det coeff[6] = (yref[1] * (xlist[2] * ylist[3] - xlist[3] * ylist[2]) + yref[2] * (ylist[1] * xlist[3] - xlist[1] * ylist[3]) + yref[3] * (xlist[1] * ylist[2] - ylist[1] * xlist[2])) / det return (OK) end # RG_POSANGLE -- Compute the position angle of a 2D vector. The angle is # measured counter-clockwise from the positive x axis. real procedure rg_posangle (x, y) real x #I the x vector component real y #I the y vector component real theta bool fp_equalr() begin if (fp_equalr (y, 0.0)) { # 0-valued y component if (x > 0.0) theta = 0.0 else if (x < 0.0) theta = PI else theta = 0.0 } else if (fp_equalr (x, 0.0)) { # 0-valued x component if (y > 0.0) theta = PI / 2.0 else if (y < 0.0) theta = 3.0 * PI / 2.0 else theta = 0.0 } else if (x > 0.0 && y > 0.0) { # 1st quadrant theta = atan (y / x) } else if (x > 0.0 && y < 0.0) { # 4th quadrant theta = 2.0 * PI + atan (y / x) } else if (x < 0.0 && y > 0.0) { # 2nd quadrant theta = PI + atan (y / x) } else if (x < 0.0 && y < 0.0) { # 3rd quadrant theta = PI + atan (y / x) } return (theta) end # RG_CTOGEO -- Transform the linear transformation coefficients to useful # geometric parameters. procedure rg_ctogeo (a, b, c, d, xscale, yscale, xrot, yrot) real a #I the x coefficient of the x coordinate fit real b #I the y coefficient of the x coordinate fit real c #I the x coefficient of the y coordinate fit real d #I the y coefficient of the y coordinate fit real xscale #I output x scale real yscale #I output y scale real xrot #I rotation of point on x axis real yrot #I rotation of point on y axis bool fp_equalr() begin xscale = sqrt (a * a + c * c) yscale = sqrt (b * b + d * d) # Get the x and y axes rotation factors. if (fp_equalr (a, 0.0) && fp_equalr (c, 0.0)) xrot = 0.0 else xrot = RADTODEG (atan2 (-c, a)) if (xrot < 0.0) xrot = xrot + 360.0 if (fp_equalr (b, 0.0) && fp_equalr (d, 0.0)) yrot = 0.0 else yrot = RADTODEG (atan2 (b, d)) if (yrot < 0.0) yrot = yrot + 360.0 end