C Assumes that the files have been sorted by y!
C Finished 6/21/00 for Duffy SMC project
C Rewritten Feb 15, 2002 for finding WRs in IC10. Again.
	character*40 wcfile,wnfile,ctfile
	parameter(nmax=10000)
	dimension wcx(nmax),wcy(nmax),
     &            wnx(nmax),wny(nmax),
     &            ctx(nmax),cty(nmax),
     &            wcm(nmax),wnm(nmax),ctm(nmax),
     &            wce(nmax),wne(nmax),cte(nmax)
	dimension wcwnm(nmax),wcctm(nmax),wnctm(nmax),
     &            wcwne(nmax),wccte(nmax),wncte(nmax),
     &            wcchi(nmax),wnchi(nmax),ctchi(nmax)
	dimension matwc(nmax),matwn(nmax),matct(nmax)

	data hit/2.0/
	data cut/5.0/
	data zero/8.05/
C in pixel values.
	hit2=hit**2
5	write(6,1) 
1	format('Enter wc filename; e.g., wc1s...:',$)
	read(5,22) wcfile
22	format(a40)
	open(unit=1,file=wcfile,status='old',err=2)
	go to 4
2	write(6,3) 
3	format('PROBLEM WITH THIS FILE. TRY AGAIN!')
	go to 5
4	nwc=0
6	read(1,*,end=11) x,y,amagm,sigexp,chi
	nwc=nwc+1
	wcx(nwc)=x
	wcy(nwc)=y
	wcm(nwc)=amagm-zero
	wce(nwc)=sigexp
	wcchi(nwc)=chi
Cwrite(6,991) nwc,wcx(nwc),wcy(nwc),wcm(nwc),wce(nwc),wcchi(nwc)
991	format(i5,2f8.2,3f8.3)
	go to 6
11	close(unit=1)
	write(6,12) wcfile,nwc
	write(7,12) wcfile,nwc
12	format('File:',a40,' contains',i5,'stars.')
16	write(6,13)
13	format('Enter wn filename; e.g., wn1s...:',$)
	read(5,22) wnfile
	open(unit=1,file=wnfile,status='old',err=15)
	go to 14
15	write(6,3)
	go to 16
14	nwn=0
18	read(1,*,end=21) x,y,amagm,sigexp,chi
	nwn=nwn+1
	wnx(nwn)=x
	wny(nwn)=y
	wnm(nwn)=amagm-zero
	wne(nwn)=sigexp
	wnchi(nwn)=chi
Cwrite(6,991) nwn,wnx(nwn),wny(nwn),wnm(nwn),wne(nwn),wnchi(nwn)
	go to 18
21	close(unit=1)
	write(6,12) wnfile,nwn
	write(7,12) wnfile,nwn
26	write(6,23)
23	format('Enter ct filename; e.g., ct1s...:',$)
	read(5,22) ctfile
	open(unit=1,file=ctfile,status='old',err=25)
	go to 24
25	write(6,3)
	go to 26
24	nct=0
28	read(1,*,end=31) x,y,amagm,sigexp,chi
	nct=nct+1
	ctx(nct)=x
	cty(nct)=y
	ctm(nct)=amagm-zero
	cte(nct)=sigexp
	ctchi(nct)=chi
Cwrite(6,991) nct,ctx(nct),cty(nct),ctm(nct),cte(nct),ctchi(nct)
	go to 28
31	close(unit=1)
	write(6,12) ctfile,nct
	write(7,12) ctfile,nct
C whew!
C Now get offsets
	write(6,32) 
32	format('Enter X(WC-WN), Y(WC-WN):',$)
	read(5,*) dxwcwn,dywcwn
	write(6,33)
33	format('Enter X(WC-CT), Y(WC-CT):',$)
	read(5,*) dxwcct,dywcct
C therefore:
	dxwnct=dxwcct-dxwcwn
	dywnct=dywnct-dywnct
C Now go looking for hits:
	do 100 iwc=1,nwc
	ihitwcwn=0
	smawcwn=5*hit2
	do 200 j=1,nwn
	wnxp=wnx(j)+dxwcwn
	wnyp=wny(j)+dywcwn
C find the clostest match
C fist, give up if the star is too far down on the list
	if((wnyp-wcy(iwc).gt.10.*hit)) go to 202
	dist2=(wnxp-wcx(iwc))**2+(wnyp-wcy(iwc))**2
		if(dist2.lt.smawcwn) then
		smawcwn=dist2
		iwn=j
		end if
200	continue
202     if(smawcwn.gt.hit2) go to 100
C give up if there wasn't a match
C OK, now go looking for a match in the ct filter!
	smawcct=5*hit2
	do 300 k=1,nct
	ctxp=ctx(k)+dxwcct
	ctyp=cty(k)+dywcct
C find the cloest match with CT
C first, give up if the star is too far down on the list
	if((ctyp-wcy(iwc).gt.10.*hit)) go to 302
	dist2=(ctxp-wcx(iwc))**2+(ctyp-wcy(iwc))**2
		if(dist2.lt.smawcct) then
		smawcct=dist2
		ict=k
		end if
300	continue
302	if(smawcct.gt.hit2) go to 100
Cwrite(6,993) iwc,iwn,ict,smawcwn,smawcct
993	format('Found match WC to WN to CT:', i5,i5,i5,2f6.2)
C matched in all 3!
	imatch=imatch+1
	matwc(imatch)=iwc
	matwn(imatch)=iwn
	matct(imatch)=ict
	wcwnm(imatch)=wcm(iwc)-wnm(iwn)
	wcwne(imatch)=sqrt( (wce(iwc))**2 + (wne(iwn))**2)
	wcctm(imatch)=wcm(iwc)-ctm(ict)
	wccte(imatch)=sqrt( (wce(iwc))**2 + (cte(ict))**2)
	wnctm(imatch)=wnm(iwn)-ctm(ict)
	wncte(imatch)=sqrt( (wne(iwn))**2 + (cte(ict))**2)
Cwrite(6,994) imatch,wcctm(imatch),wccte(imatch),
C    &                      wnctm(imatch),wncte(imatch),
C    &                      wcwnm(imatch),wcwne(imatch)
994	format(i5,6f8.3)
100	continue
	nmatch=imatch
	write(6,101) nmatch
	write(7,101) nmatch
101	format(i5, 'matched in all 3 frames....')
	call median(wcctm,nmatch,wcctc)
	call median(wnctm,nmatch,wnctc)
	call median(wcwnm,nmatch,wcwnc)
	wcwncexp=wcctc - wnctc
	write(6,102) wcctc,wnctc, wcwnc,wcwncexp
	write(7,102) wcctc,wnctc, wcwnc,wcwncexp
102	format('WC-CT correction:',f8.3,
     &  /'WN-CT correction:',f8.3, 
     &  /'WC-WN correction:',f8.3,' expected:',f8.3)
C GREAT!  It's working, Bullwinkle!  now look for significant guys!
C
C
C
	write(6,502)
	write(7,502)
502	format('  X(WC) ','  Y(WC) ','   WC  ','  WN  ','  CT  ',
     &  ' ---WC-CT----',2x,' ---WN-CT---',2x,' ---WC-WN---')
C        123456789012
C               12345678   12345678   123456   123456   123456
	do 500 i=1,nmatch
C for every star that's matched, find the mag differences with the corrections
	dwcct=wcctm(i)-wcctc
	sigma=wccte(i)
	if(sigma.lt.0.005) sigma=0.005
	sigwcct=-dwcct/sigma
	dwnct=wnctm(i)-wnctc
	sigma=wncte(i)
	if(sigma.lt.0.005) sigma=0.005
	sigwnct=-dwnct/sigma
	dwcwn=wcwnm(i)-wcwnc
	sigma=wcwne(i)
	if(sigma.lt.0.005) sigma=0.005
	sigwcwn=abs(dwcwn)/sigma
C Whew!
C Now see if anyone is good.
	if((sigwcct.gt.cut).and.(dwcct.le.-.10)) then
	write (6,501) wcx(matwc(i)),wcy(matwc(i)),
     &                wcm(matwc(i)),wnm(matwn(i)),ctm(matct(i)),
     &                dwcct,sigwcct,dwnct,sigwnct,dwcwn,sigwcwn
	write (7,501) wcx(matwc(i)),wcy(matwc(i)),
     &                wcm(matwc(i)),wnm(matwn(i)),ctm(matct(i)),
     &                dwcct,sigwcct,dwnct,sigwnct,dwcwn,sigwcwn
501	format(2f8.2,1x,3f6.2,1x,f6.3,f6.2,2x,f6.3,f6.2,2x,f6.3,f6.2)
	else if ((sigwnct.gt.cut).and.(dwwnct.le.-.10)) then
	write (6,501) wcx(matwc(i)),wcy(matwc(i)),
     &                wcm(matwc(i)),wnm(matwn(i)),ctm(matct(i)),
     &                dwcct,sigwcct,dwnct,sigwnct,dwcwn,sigwcwn
	write (7,501) wcx(matwc(i)),wcy(matwc(i)),
     &                wcm(matwc(i)),wnm(matwn(i)),ctm(matct(i)),
     &                dwcct,sigwcct,dwnct,sigwnct,dwcwn,sigwcwn
	else if ((sigwcwn.gt.2.*cut).and.(abs(dwcwn).ge.0.2)) then
	write (6,501) wcx(matwc(i)),wcy(matwc(i)),
     &                wcm(matwc(i)),wnm(matwn(i)),ctm(matct(i)),
     &                dwcct,sigwcct,dwnct,sigwnct,dwcwn,sigwcwn
	write (7,501) wcx(matwc(i)),wcy(matwc(i)),
     &                wcm(matwc(i)),wnm(matwn(i)),ctm(matct(i)),
     &                dwcct,sigwcct,dwnct,sigwnct,dwcwn,sigwcwn
	else
	go to 500
	end if
500	continue
	stop
	end
	subroutine momment(data,n,ave,adev,sdev,var,skew,curt)
	dimension data(*)
C Given an array of data (1:n) this routine returns a bunch of junk
C from Numerical Recipies p 607
	if(n.lt.1) return
	s=0.
	do 11 j=i1,n
	s=s+data(j)
11	continue
	ave=s/n
	adev=0.
	var=0.
	skew=0.
	curt=0.
	ep=0.
	do 12 j=1,n
	s=data(j)-ave
	ep=ep+s
	adev=adev+abs(s)
	p=s*s
	var=var+p
	p=p*s
	skew=skew+p
	p=p*s
	curt=curt+p
12	continue
	adev=adev/n
	var=(var-ep**2/n)/(n-1)
	sdev=sqrt(var)
	if(var.ne.0.) then
	skew=skew/(n*sdev**3)
	curt=curt/(n*var**2)-3.
	else
	write(6,13)
13	format('No skew or kurtosis when zero variance in moment')
	end if
	return
	end
	subroutine median(data,n,amed)
	parameter (nmax=10000)
	dimension data(*),indx(nmax)
	call sort(n,data,indx)
	if( ((n/2)*2).eq.n) then
	amed=(data(indx(n/2))+data(indx(n/2+1)))/2.
	else
	amed=data(indx(n/2+1))
	end if
	return
	end
        subroutine sort(n,arrin,indx)
        dimension arrin(*),indx(*)
        do 11 j=1,n
        indx(j)=j
11      continue
        if(n.eq.1) return
        l=n/2+1
        ir=n
10      continue
                if(l.gt.1) then
                        l=l-1
                        indxt=indx(l)
                        q=arrin(indxt)
                ELSE
                        indxt=indx(ir)
                        q=arrin(indxt)
                        indx(ir)=indx(1)
                        ir=ir-1
                        if(ir.eq.1) then
                                indx(1)=indxt
                                return
                        endif
                endif
        i=l
        j=l+l
20      if(j.le.ir) then
                if(j.lt.ir) then
                        if(arrin(indx(j)).lt.arrin(indx(j+1)))j=j+1
                endif
                if(q.lt.arrin(indx(j)))then
                        indx(i)=indx(j)
                        i=j
                        j=j+j
                else
                        j=ir+1
                endif
        go to 20
        endif
        indx(i)=indxt
        go to 10
        end
c subroutine xfit.f
c
c source
c   Bevington, page 76.
c
c purpose
c   calculate the mean and estimated errors for a set of data points
c
c usage
c   call xfit (x, sigmax, npts, mode, xmean, sigmam, sigma)
c
c description of parameters
c   x      - array of data points
c   sigmax - array of standard deviations for data points
c   npts   - number of data points
c   mode   - determines method of weighting
c            +1 (instrumental) weight(i) = 1./sigmax(i)**2
c             0 (no weighting) weight(i) = 1.
c            -1 (statistical)  weight(i) = 1.
c   xmean  - weighted mean
c   sigmam - standard deviation of mean
c   sigma  - standard deviation of data
c
c subroutines and function subprograms required
c   none
c
	subroutine xfit (x,sigmax,npts,mode,xmean,sigmam,sigma)
	double precision sum,sumx,weight,free
	dimension x(1),sigmax(1)
c
c accumulate weighted sums
c
11	sum=0.
	sumx=0.
	sigma=0.
	sigmam=0.
20	do 32 i=1,npts
21	if (mode) 22,22,24
22	weight=1.
	goto 31
24	weight=1./sigmax(i)**2
31	sum=sum+weight
32	sumx=sumx+weight*x(i)
c
c evaluate mean and standard deviations
c
41	xmean=sumx/sum
51	do 52 i=1,npts
52	sigma=sigma+(x(i)-xmean)**2
	free=npts-1
54	sigma=dsqrt(sigma/free)
61	if (mode) 62,64,66
62	sigmam=dsqrt(xmean/sum)
	goto 70
64	sigmam=sigma/dsqrt(sum)
	goto 70
66	sigmam=dsqrt(1./sum)
70	return
	end
