diff options
Diffstat (limited to 'noao/nproto/skygroup.cl')
| -rw-r--r-- | noao/nproto/skygroup.cl | 195 |
1 files changed, 195 insertions, 0 deletions
diff --git a/noao/nproto/skygroup.cl b/noao/nproto/skygroup.cl new file mode 100644 index 00000000..24fa05aa --- /dev/null +++ b/noao/nproto/skygroup.cl @@ -0,0 +1,195 @@ +# SKYGROUP -- Group coordinate list on the sky. +# +# The input is list of ra (0h-24h or 0d-360d) and dec (-90d to 90d) in the +# first two columns followed by arbitrary data (usually a filename). +# This is complicated by the periodicities at 0h. + +procedure skygroup (input, output) + +file input {prompt="Input list"} +string output {prompt="Output rootname"} +string extn = "" {prompt="Optional output extension"} +real sep = 60 {prompt="Separation between groups (arcsec)"} +string raunit = "hr" {prompt="RA unit (hr|deg)"} +bool keepcoords = yes {prompt="Keep coordinates in output?"} +string raformat = "%.2h" {prompt="Output RA format"} +string decformat = "%.1h" {prompt="Output DEC format"} + +struct *fd1, *fd2 + +begin + file in, out, fname, temp1, temp2, temp3 + int i, j, n, n1 + real dra, r1, d1, r2, d2, r3, d3, r4 + struct fmtstr, data1, data2, data3 + + # Temporary files. + fname = mktemp ("tmp") + in = fname // "in" + temp1 = fname // "1" + temp2 = fname // "2" + temp3 = fname // "3" + + # Set parameters. + fname = input + out = output + + # Check for existing output files. + files (out//"_[0-9][0-9][0-9][0-9]*"//extn, > temp1) + count (temp1) | scan (n); delete (temp1, verify-) + if (access(out) || n > 0) + error (1, "Output files already exist") + + if (raunit == "hr") + dra = 24 + else + dra = 360 + + if (keepcoords) + fmtstr = "%d " // raformat // " " // decformat // " %s" + else + fmtstr = "%d %s" + fmtstr += "\n" + + # We start by sorting in dec. + sort (fname, col=2, num+, > in) + + # Find jumps in dec bigger than the separation and then sort + # in ra and find jumps in ra bigger than separation. Handle + # the wrap around at 0h by duplicating to extend beyond 24h. + # The duplicates will be eliminated during the merging process. + + n = 1 + fd1 = in + if (fscan (fd1, r1, d1, data1) == EOF) + error (1, "No data or badly formated data") + while (fscan (fd1, r2, d2, data2) != EOF) { + print (r1, d1, data1, >> temp1) + if (r1 / max (0.001, dcos(d1)) * 3600 <= sep) { + r4 = r1 + dra + print (r4, d1, data1, >> temp1) + } + if (abs(d2-d1) <= sep) { + r1 = r2; d1 = d2; data1 = data2 + next + } + + r3 = r2; d3 = d2; data3 = data2 + + sort (temp1, col=1, num+, > temp2) + delete (temp1, verify-) + + fd2 = temp2 + if (fscan (fd2, r1, d1, data1) == EOF); + while (fscan (fd2, r2, d2, data2) != EOF) { + if (keepcoords) + printf (fmtstr, n, r1, d1, data1, >> temp3) + else + printf (fmtstr, n, data1, >> temp3) + skysep (r1, d1, r2, d2, raunit=raunit, verb-) + if (skysep.sep <= sep) { + r1 = r2; d1 = d2; data1 = data2 + next + } + + n += 1 + r1 = r2; d1 = d2; data1 = data2 + } + fd2 = ""; delete (temp2, verify-) + + if (keepcoords) + printf (fmtstr, n, r1, d1, data1, >> temp3) + else + printf (fmtstr, n, data1, >> temp3) + n += 1 + + r1 = r3; d1 = d3; data1 = data3 + } + fd1 = ""; delete (in, verify-) + + print (r1, d1, data1, >> temp1) + if (r1 / max (0.001, dcos(d1)) * 3600 <= sep) { + r4 = r1 + dra + print (r4, d1, data1, >> temp1) + } + + sort (temp1, col=1, num+, > temp2) + delete (temp1, verify-) + + fd2 = temp2 + if (fscan (fd2, r1, d1, data1) == EOF); + while (fscan (fd2, r2, d2, data2) != EOF) { + if (keepcoords) + printf (fmtstr, n, r1, d1, data1, >> temp3) + else + printf (fmtstr, n, data1, >> temp3) + skysep (r1, d1, r2, d2, raunit=raunit, verb-) + if (skysep.sep <= sep) { + r1 = r2; d1 = d2; data1 = data2 + next + } + + n += 1 + r1 = r2; d1 = d2; data1 = data2 + } + fd2 = ""; delete (temp2, verify-) + + if (keepcoords) + printf (fmtstr, n, r1, d1, data1, >> temp3) + else + printf (fmtstr, n, data1, >> temp3) + + # Now write out the lists and check for duplicate which must be + # merged. + + sort (temp3, col=1, num+, > temp1); delete (temp3, verify-) + touch (temp2) + fd1 = temp1 + if (fscan (fd1, i, data1) == EOF); + while (fscan (fd1, j, data2) != EOF) { + if (data1 == data2) { + print (j, i, >> temp2) + next + } + printf ("%s_%03d%s\n", out, i, extn) | scan (fname) + print (data1, >> fname) + i = j; data1 = data2 + } + fd1 = ""; delete (temp1, verify-) + printf ("%s_%03d%s\n", out, i, extn) | scan (fname) + print (data1, >> fname) + sort (temp2, col=1, num+, rev+) | unique (> temp1) + delete (temp2, verify-) + + # Merge the lists. + n1 = n + fd1 = temp1 + while (fscan (fd1, j, i) != EOF) { + printf ("%s_%03d%s\n", out, j, extn) | scan (fname) + if (access (fname)) { + printf ("%s_%03d%s\n", out, i, extn) | scan (in) + concat (in, fname, append+) + delete (in, verify-) + n1 -= 1 + } + } + fd1 = ""; delete (temp1, verify-) + + # Renumber if needed. + if (n1 != n) { + i = 1 + for (j=1; j<=n; j+=1) { + printf ("%s_%03d%s\n", out, j, extn) | scan (fname) + if (access(fname)) { + if (i != j) { + printf ("%s_%03d%s\n", out, i, extn) | scan (in) + rename (fname, in) + } + i += 1 + } + } + } + + # Create the final output list of lists. + files (out//"_[0-9]*", > out//extn) +end |