diff options
author | Joe Hunkeler <jhunkeler@gmail.com> | 2015-08-11 16:51:37 -0400 |
---|---|---|
committer | Joe Hunkeler <jhunkeler@gmail.com> | 2015-08-11 16:51:37 -0400 |
commit | 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 (patch) | |
tree | 4464880c571602d54f6ae114729bf62a89518057 /noao/digiphot/daophot/allstar/dpregroup.x | |
download | iraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz |
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'noao/digiphot/daophot/allstar/dpregroup.x')
-rw-r--r-- | noao/digiphot/daophot/allstar/dpregroup.x | 219 |
1 files changed, 219 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/allstar/dpregroup.x b/noao/digiphot/daophot/allstar/dpregroup.x new file mode 100644 index 00000000..a7561fb8 --- /dev/null +++ b/noao/digiphot/daophot/allstar/dpregroup.x @@ -0,0 +1,219 @@ +# DP_ALSORT -- Sort the stars on the y coordinate. + +procedure dp_alsort (id, x, y, mag, sky, sumwt, dxold, dyold, xclamp, yclamp, + nstar) + +int id[ARB] # array if star ids +real x[ARB] # array of star x values +real y[ARB] # array of y values +real mag[ARB] # array of star magnitudes +real sky[ARB] # array of star sky values +real sumwt[ARB] # array of weighted sums +real dxold[ARB] # the previous star x array +real dyold[ARB] # the previous star y array +real xclamp[ARB] # the x array clamps +real yclamp[ARB] # the y array clamps +int nstar # the number of stars + +int ier +pointer sp, index + +begin + # Allocate some working memory. + call smark (sp) + call salloc (index, nstar, TY_INT) + + # Sort the star list on y. + call quick (y, nstar, Memi[index], ier) + + # Recitfy the remaining arrays. + call dp_irectify (id, Memi[index], nstar) + call dp_rectify (x, Memi[index], nstar) + call dp_rectify (mag, Memi[index], nstar) + call dp_rectify (sky, Memi[index], nstar) + call dp_rectify (sumwt, Memi[index], nstar) + call dp_rectify (dxold, Memi[index], nstar) + call dp_rectify (dyold, Memi[index], nstar) + call dp_rectify (xclamp, Memi[index], nstar) + call dp_rectify (yclamp, Memi[index], nstar) + + # Release memory. + call sfree (sp) +end + + +# DP_REGROUP -- Group stars into physical associations based on proximity. + +procedure dp_regroup (id, x, y, mag, sky, chi, dxold, dyold, xclamp, yclamp, + nstar, radius, last) + +int id[ARB] # array if star ids +real x[ARB] # array of star x values +real y[ARB] # array of y values +real mag[ARB] # array of star magnitudes +real sky[ARB] # array of star sky values +real chi[ARB] # array of chis +real dxold[ARB] # the previous star x array +real dyold[ARB] # the previous star y array +real xclamp[ARB] # the x array clamps +real yclamp[ARB] # the y array clamps +int nstar # the number of stars +real radius # the fitting radius +int last[ARB] # the last array (NO or YES for last star) + +int itop, itest, j, k, i, ier +pointer sp, index +real critrad, critradsq, xtest, ytest, dx, dy + +begin + # If there is only one stars its value of last is YES. + if (nstar <= 1) { + last[1] = YES + return + } + + # Allocate some working memory. + call smark (sp) + call salloc (index, nstar, TY_INT) + + # Compute the critical radius. + critrad = 2.0 * radius + critradsq = critrad * critrad + + # Sort the star list on y and rectify the accompanying x array. + call quick (y, nstar, Memi[index], ier) + call dp_rectify (x, Memi[index], nstar) + + # At this point the stars are in a stack NSTAR stars long. The + # variable ITEST will point to the position in the stack + # occupied by the star which is currently the center of a + # circle of the CRITRAD pixels, within which we are + # looking for other stars. ITEST starts out with a value of one. + # ITOP points to the top position in the stack of the + # stars which have not yet been assigned to groups. ITOP starts + # out with the value of 2. Each time through, the program goes + # down through the stack from ITOP and looks for stars within + # the critrad pixels from the star at stack position ITEST. + # When such a star is found, it changes places in the + # stack with the star at ITOP and ITOP is incremented by one. + # When the search reaches a star J such that Y(J) - Y(ITEST) + # > CRITRAD we know that no further stars will be found + # within the CRITRAD pixels, the pointer ITEST is incremented + # by one, and the search proceeds again from the new + # value of ITOP. If the pointer ITEST catches up with the + # pointer ITOP, then the current group is complete + # with the star at the position ITEST-1, and + # ITOP is incremented by one. If ITOP catches up with NSTAR + # the grouping process is complete. + + itop = 2 + do itest = 1, nstar { + + # Initialize the last array. + last[itest] = NO + + # ITEST has reached ITOP; so no other unassigned stars are + # within CRITRADIUS pixels of any member of the current + # group. The group is therefore complete. Signify this by + # setting LAST[ITEST-1] = YES, and incrementing the value of + # ITOP by one. If ITOP is greater than NSTAR at this point + # list is finished. + + if (itest == itop) { + j = itest - 1 + if (j > 0) + last[j] = YES + itop = itop + 1 + if (itop > nstar) { + last[itest] = YES + break + } + } + + # Now go through the list of unassigned stars, occupying + # positions ITOP through NSTAR in the stack, to look for + # stars within CRITRADIUS pixels of the star at + # position ITEST inthe stack. If one is found, move it + # up to stack position ITOP and increment ITOP by one. + + xtest = x[itest] + ytest = y[itest] + j = itop + + do i = j, nstar { + + # Check the distance between the stars. + dy = y[i] - ytest + if (dy > critrad) + break + dx = x[i] - xtest + if (abs (dx) > critrad) + next + if ((dx * dx + dy * dy) > critradsq) + next + + # This star is within CRITRAD pixels of the + # star at stack position ITEST. Therefore it is + # added to the current group by moving it up to position + # ITOP in the stack, where the pointer ITEST may reach it. + + call dp_rgswap (itop, i, x, y, Memi[index]) + + # Now increment ITOP by 1 to point at the top most + # unassigned star in the stack. + itop = itop + 1 + if (itop > nstar) { + do k = itest, nstar - 1 + last[k] = NO + last[nstar] = YES + break + } + } + + if (itop > nstar) + break + } + + # Reorder the remaining arrays to match the x and y arrays. + call dp_irectify (id, Memi[index], nstar) + call dp_rectify (mag, Memi[index], nstar) + call dp_rectify (sky, Memi[index], nstar) + call dp_rectify (chi, Memi[index], nstar) + call dp_rectify (dxold, Memi[index], nstar) + call dp_rectify (dyold, Memi[index], nstar) + call dp_rectify (xclamp, Memi[index], nstar) + call dp_rectify (yclamp, Memi[index], nstar) + + call sfree (sp) +end + + +# DP_RGSWAP -- Swap the i-th and j-th stars in the stack without otherwise +# altering the order of the stars.. + +procedure dp_rgswap (i, j, x, y, index) + +int i, j # indices of the two stars to be swapped +real x[ARB] # the array of x values +real y[ARB] # the array of y values +int index[ARB] # the index array + +int ihold, k,l +real xhold, yhold + +begin + xhold = x[j] + yhold = y[j] + ihold = index[j] + + do k = j, i + 1, -1 { + l = k - 1 + x[k] = x[l] + y[k] = y[l] + index[k] = index[l] + } + + x[i] = xhold + y[i] = yhold + index[i] = ihold +end |