-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathclose_match_radec.pro
158 lines (150 loc) · 4.39 KB
/
close_match_radec.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
pro close_match_radec,t1,s1,t2,s2,m1,m2,ep,allow,miss1,silent=silent
;+
; NAME:
; close_match_radec
;
; PURPOSE:
; this will find close matches between 2 sets of points (t1,s1)
; and (t2,s2) (note t1,t2,s1,s2 are all arrays) in ra dec space.
;
; CALLING SEQUENCE:
; close_match_radec,t1,s1,t2,s2,m1,m2,ep,allow,miss1
;
; INPUTS:
; t1,s1: the ra dec of the first set
; t2,s2: the ra dec of the second set
; ep: this is the error which defines a close match. A pair is considered
; a match if |t1-t2|/cos(dec) AND |s1-s2| are both less than ep. This is faster
; than doing a euclidean measure on the innermost loop of the program
; and just as good.
; allow: how many matches in the (t2,s2) space will you allow for
; each (t1,s1)
;
; OUTPUTS:
; m1,m2: the indices of the matches in each space. That is
; (t1(m1),s1(m1)) matches (t2(m2),s2(m2))
; miss1: this gives the index of the things in x1 NOT found to match (optional)
;
; OPTIONAL KEYWORD PARAMETERS:
; none
;
; NOTES:
; It sorts the t2 list so that it can do a binary search for each t1.
; It then carves out the sublist of t2 where it is off by less than ep.
; It then only searches that sublist for points where s1 and s2 differ
; by less than ep.
; PROCEDURES CALLED:
; binary_search, rem_dup
; REVISION HISTORY:
; written by David Johnston -University of Michigan June 97
;
; Tim McKay August 97
; bug fixed, matcharr extended to "long" to accomodate ROTSE I images
; Tim McKay 6/23/98
; Altered to be an ra dec match, with appropriate scaling of ra range...
; Tim McKay 7/8/99
; Altered to order matches by distance, rather than just ra dec distance
;-
On_error,2 ;Return to caller
if N_params() LT 8 then begin
print,'Syntax - close_match,ra1,dec1,ra2,dec2,m1,m2,ep,allow,miss1,silent=silent'
return
endif
; first sort out the allowed errors in ra and dec.....
epdec=ep
n1=n_elements(t1)
n2=n_elements(t2)
matcharr=lonarr(n1,allow) ;the main book-keeping device for
matcharr(*,*)=-1 ;matches -initialized to -1
ind=lindgen(n2)
sor=sort(t2) ;sort t2
t2s=t2[sor]
s2s=s2[sor]
ind=ind[sor] ;keep track of index
runi=0
endt=t2s[n2-1]
for i=0l , n1-1l do begin ;the main top level loop over t1
t=t1[i]
dec=s1[i]
epra=ep/cos(dec*0.01745329)
tm=t-epra ;sets the range of good ones
tp=t+epra
binary_search,t2s,tm,in1
;searched for the first good one
if in1 eq -1 then if tm lt endt then in1=0
;in case tm smaller than all t2 but still may be some matches
if in1 ne -1 then begin
in1=in1+1
in2=in1-1
jj=in2+1
while jj lt n2 do begin
if t2s[in2+1] lt tp then begin
in2=in2+1 & jj=jj+1
endif else jj=n2
endwhile
if (n2 eq 1) then in2 = 1
;while loop carved out sublist to check
;a little tricky,be careful
if in1 le in2 then begin
if (n2 ne 1) then begin
check=s2s[in1:in2] ;the sublist to check
tcheck=t2s[in1:in2]
endif else begin
check=s2s[0]
tcheck=t2s[0]
endelse
s=s1[i]
offby=abs(check-s)
toffby=abs(tcheck-t1[i])
good=where(offby lt epdec and toffby lt epra,ngood)+in1
;the selection made here
if ngood ne 0 then begin
if ngood gt allow then begin
;now calculate real distances
;Dave's old way....
;offby=offby(good-in1)
;The new way....
offby=sphdist(t1[i],s1[i],$
t2s[good],s2s[good],/degrees)
good=good[sort(offby)];sorts by closeness
ngood=allow
;not more than you are allowed by 'allow'
endif
good=good[0:ngood-1]
matcharr[i,0:ngood-1]=good
;finally the matches are entered in
runi=runi+ngood ;a running total
endif
endif
endif
endfor
if not keyword_set(silent) then print,'total put in bytarr',runi
matches=where(matcharr ne -1,this)
if this eq 0 then begin
if not keyword_set(silent) then $
print,'no matches found'
m1=-1 & m2=-1
return
endif
m1=matches mod n1 ;a neat trick to extract them correctly
m2=matcharr[matches] ;from the matcharr matrix
if not keyword_set(silent) then $
print,n_elements(m1),' matches'
m2=ind[m2] ;remember, must unsort
dif=m1[uniq(m1,sort(m1))]
if not keyword_set(silent) then $
print,n_elements(dif),' different matches'
if n_params() eq 9 then begin
if n_elements(m1) lt n1 then begin
miss1=lindgen(n1)
remove,dif,miss1
if not keyword_set(silent) then $
print,n_elements(miss1),' misses'
endif else begin
miss1=-1
if not keyword_set(silent) then $
print,'no misses'
endelse
endif
return
end