C rewrote to deal with n6822 data from Chuck
C dec 16 1987
C modified June 1991 to use photcal output files, e.g.,
C (variable length name), V, VErr, BMV, BMVErr, UMB, UMBErr
	Parameter(maxfiles=10,maxstars=15000)
	implicit double precision (a-h,o-z)
	character*80 filein(maxfiles)
	character*111 inputline
	character*1 sign
	integer irvec1(4),idvec1(4)
	double precision raf(maxstars),decf(maxstars),
     &                   ran(maxstars),decn(maxstars),
     &                   as2r,s2r,
     &                   stor,astor
	dimension xf(maxstars),yf(maxstars),amagf(maxstars,3,maxfiles),
     .amagef(maxstars,3,maxfiles),xn(maxstars),yn(maxstars),
     .amagn(maxstars,3),amagen(maxstars,3),df(maxstars),dn(maxstars),
     .nums(maxfiles),numhere(maxstars),vvalue(maxfiles),
     .vsigma(maxfiles),
     .umbvalue(maxfiles),umbsigma(maxfiles),bmvvalue(maxfiles),
     .bmvsigma(maxfiles)
*  SECONDS TO RADIANS
      PARAMETER (S2R=0.72722052166430398D-4)
*  ARC SECONDS TO RADIANS
      PARAMETER (AS2R=0.48481368110953599D-5)
*  SECONDS TO RADIANS
      STOR(RH,RM,RS)=S2R*(60D0*(60D0*RH+RM)+RS)

*  ARC SECONDS TO RADIANS
      ASTOR(DD,DM,DS)=AS2R*(60D0*(60D0*DD+DM)+DS)

	write(6,1)
1	format(1x,'number of files?')
	read(5,2) nfiles
2	format(i3)
C so far, so good!
	if(nfiles.gt.maxfiles) then
	write(6,33) 
33	format(1x,'too many files!')
	stop
	end if
Copen(unit=2,file='total.out',status='new')
	do 5 i=1,nfiles
	write(6,3) i 
3	format('input file name',i3,':',$)
	read (5,4) filein(i)
4	format(a80)
	open(unit=1,file=filein(i),status='old')
	n=0
7	read(1,9,end=8) inputline
9	format(a111)
	if((inputline(1:1).eq.' ').or.(inputline(1:1).eq.'#')) go to 7
	n=n+1
	go to 7
8	close(unit=1)
	nums(i)=n
	write(6,11) n
11	format(1x,i4,' stars this file...')	
5	continue
C OK now get the data for each file; treat the 1st file a little bit
C special
	n=0
	open(unit=1,file=filein(1),status='old')
21	n=n+1
23	read(1,9,end=25) inputline
	if(inputline(1:1).eq.'#') go to 23
	if(inputline(1:1).eq.' ') go to 23
        read(inputline(1:35),24) ihr,imin,sec,sign,ideg,iamin,asec
24      format(i2,3x,i2,3x,f6.3,3x,a1,i2,3x,i2,3x,f5.2)
        if(sign.eq.'-') then
        ideg=-ideg
        iamin=-iamin
        asec=-asec
        end if
        raf(n)=stor(ihr,imin,sec)
        decf(n)=astor(ideg,iamin,asec)
Cwrite(6,*) ideg,iamin,asec,decf(n)
        read(inputline(51:111),*) xf(n),yf(n),
     &        (amagf(n,j,1),amagef(n,j,1), j=1,3)
	numhere(n)=1
	go to 21
25	numf=n-1
c OK, now go get the next bunch of files....
	close(unit=1)
	do 101 i=2,nfiles
	open(unit=1,file=filein(i),status='old')
	n=0
100	n=n+1
123	read(1,9,end=105) inputline
	if(inputline(1:1).eq.'#') go to 123
	if(inputline(1:1).eq.' ') go to 123
        read(inputline(1:35),24) ihr,imin,sec,sign,ideg,iamin,asec
        if(sign.eq.'-') then
        ideg=-ideg
        iamin=-iamin
        asec=-asec
        end if
        ran(n)=stor(ihr,imin,sec)
        decn(n)=astor(ideg,iamin,asec)
        read(inputline(51:111),*) xn(n),yn(n),
     &        (amagn(n,j),amagen(n,j), j=1,3)
	dn(n)=amagn(n,1)
	go to 100
105	close(unit=1)
	n=n-1
	do 106 j=1,numf
106	df(j)=amagf(j,1,1)
	call moff(xf,yf,df,numf,xn,yn,dn,n,x1m2,y1m2)
c go through each of the 'new' stars and see if they match any of
c the old.  If they do, add them in, if they don't, tack them on the
c end
	numnew=0
	do 108 k=1,n
	small=5e6
	do 107 j=1,numf
c match?
	dif2=(xf(j)-xn(k)-x1m2)**2+
     .       (yf(j)-yn(k)-y1m2)**2
	dif=sqrt(dif2)
	if(dif.lt.small) then
	ical=j
	small=dif
	end if
107	continue
c do we have a match?
	if(small.lt.2.0) then
c yes!
		numhere(ical)=numhere(ical)+1
		numh=numhere(ical)
		if(numh.gt.i) then
		write(6,1112) ical,k,i,xf(ical),yf(ical),xn(k),yn(k)
1112	format('Problem with star',i4,' and star',i4,' of file',
     &          i3,/'Star in first file:',2f8.2,
     &             /'Star in 2nd   file:',2f8.2)
	end if
		do 109 j=1,3
		amagf(ical,j,numh)=amagn(k,j)
		amagef(ical,j,numh)=amagen(k,j)
109	continue
	else
	numnew=numnew+1
	do 110 j=1,3
	amagf(numnew+numf,j,1)=amagn(k,j)
	amagef(numnew+numf,j,1)=amagen(k,j)
110	continue
	xf(numnew+numf)=xn(k)+x1m2
	yf(numnew+numf)=yn(k)+y1m2
	raf(numnew+numf)=ran(k)
	decf(numnew+numf)=decn(k)
	numhere(numnew+numf)=1
	end if
108	continue
	numf=numf+numnew
101	continue
	do 202 i=1,numf
	k1=numhere(i)
	   do 203 j=1,k1
	   vvalue(j)=amagf(i,1,j)
	   vsigma(j)=amagef(i,1,j)
	   bmvvalue(j)=amagf(i,2,j)
	   bmvsigma(j)=amagef(i,2,j)
	   umbvalue(j)=amagf(i,3,j)
203	   umbsigma(j)=amagef(i,3,j)
	   call xfit (vvalue,vsigma,k1,1,vmean,vsigm,vsig)
	   call xfit (bmvvalue,bmvsigma,k1,1,bmvmean,bmvsigm,bmvsig)
	   call xfit (umbvalue,umbsigma,k1,1,umbmean,umbsigm,umbsig)
	call dr2tf(3,raf(i),sign,irvec1)
	call dr2af(2,decf(i),sign,idvec1)
	   write(2,205) irvec1,sign,idvec1,
     .                  xf(i),yf(i),vmean,bmvmean,umbmean,
     .k1,vsigm,bmvsigm,umbsigm
205	format(i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,
     &      a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2,
     &         2f7.1,3f6.2,i2,3f6.2)
202	continue
	stop
	end
c
c
	subroutine xfit(x,sigmax,npts,mode,xmean,sigmam,sigma)
	implicit double precision (a-h,o-z)
	double precision sum,sumx,weight,free
	dimension x(1),sigmax(1)
11	sum=0.
	sumx=0.
	sigma=0.
	sigmam=0.
	if(npts.eq.0) return
	if(npts.eq.1) then
	xmean=x(1)
	sigma=sigmax(1)
	sigmam=sigmax(1)
	return
	end if
20	do 32 i=1,npts
21	if(mode) 22,22,24
22	weight=1.
	go to 31
24	if(sigmax(i).lt.0.007) sigmax(i)=0.007
	weight=1./sigmax(i)**2
C small little cheat removed here, damn it!
31	sum=sum+weight
32	sumx=sumx+weight*x(i)
c find mean and s.d.
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)
	go to 70
64	sigmam=sigma/dsqrt(sum)
	go to 70
66	sigmam=dsqrt(1./sum)
70	return
	end
	subroutine moff(x1,y1,v1,n1,x2,y2,v2,n2,offx1m2,offy1m2)
	implicit double precision (a-h,o-z)
	dimension x1(1),y1(1),v1(1),x2(1),y2(1),v2(1)
c input arrays x1,y1 containing positions of n1 stars
c              x2,y2 containing positions of n2 stars
c              v1,v2 contain mags (must agree within .2 to be considered)
c  returns offx1m2,offx1m2 offsets
	dimension ixbin(4096),iybin(4096)
	do 1 i=1,4096
	ixbin(i)=0
	iybin(i)=0
1	continue
c zero the bins
	do 100 i=1,n1
	do 101 j=1,n2
	xdif=x1(i)-x2(j)
	ydif=y1(i)-y2(j)
c get rid of things with very different magnitudes
	vdif=v1(i)-v2(j)
	if(abs(vdif).gt.0.2) go to 101
	ix=((xdif+512.)/.25)+1
	iy=((ydif+512.)/.25)+1
c figure out which if any bin it will go into
	if((ix.lt.1).or.(ix.gt.4096)) go to 101
	if((iy.lt.1).or.(iy.gt.4096)) go to 101
c just checking.
	ixbin(ix)=ixbin(ix)+1
	iybin(iy)=iybin(iy)+1
101	continue
100	continue
	ixm=0
	iym=0
	do 200 i=1,4096
c OK, so which bins have the most in them?
	write(17,5) i,ixbin(i),iybin(i)
5	format(1x,i5,2i6)
	if(ixbin(i).ge.ixm)then
		ix=i
		ixm=ixbin(i)
	end if
	if(iybin(i).ge.iym) then
		iy=i
		iym=iybin(i)
	end if
200	continue
		xave=0.
		ixave=0
		yave=0.
		iyave=0
		do 50 i=1,3
c average over the 5 bins on either side
		i1=ix-2+i
		j1=iy-2+i
		ixave=ixave+ixbin(i1)
		xave=xave+ixbin(i1)*i1
		iyave=iyave+iybin(j1)
		yave=yave+iybin(j1)*j1
50		continue
	xave=xave/ixave
	yave=yave/iyave
	offx1m2=(xave-1)*.25-512.
	offy1m2=(yave-1)*.25-512.
	if((ixm.lt.6).or.(iym.lt.6)) then
	write(6,1110) x1(5),y1(5),v1(5),x2(5),y2(5),v2(5),n1,n2
1110	format(1x,'5th stars',6f8.3,2i5)
	write(6,1102) ixm,offx1m2,offy1m2
1102	format(1x,'Not many stars in common! Offsets are', i3,2f10.2)
	write(6,1103) 
1103	format(1x,'Enter new offsets or CR to keep the old')
	read(5,1104) xa,ya
1104	format(2f8.3)
C	if((xa.eq.0.000).and.(ya.eq.0.000)) return
	offx1m2=xa
	offy1m2=ya
	end if
	write(6,1101) offx1m2,offy1m2
1101	format(1x,'average offsets:',2f10.2)
	return
	end


      SUBROUTINE DR2AF(NDP,ANGLE,SIGN,IDMSF)
*  
*  
*  CONVERT ANGLE FOR OUTPUT
*  ------------------------
*  
*  
*  GIVEN:  
*     NDP           NUMBER OF DECIMAL PLACES OF ARC SECONDS
*     ANGLE         ANGLE IN RADIANS (DOUBLE PRECISION)
*  
*  RETURNED:
*     SIGN          '+' OR '-'
*     IDMSF         DEGREES, ARC MINUTES, ARC SECONDS, FRACTION
*  
*  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*  
      INTEGER NDP  
      DOUBLE PRECISION ANGLE
      CHARACTER*1 SIGN
      INTEGER IDMSF(4)
   
*  RADIANS TO ARC SECONDS  
      DOUBLE PRECISION R2AS
      PARAMETER (R2AS=206264.80624709636D0)

*  HANDLE SIGN 
      SIGN='+' 
      IF(ANGLE.LT.0D0) SIGN='-'

*  FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE
      RAS=REAL(10**MAX(NDP,0))
      RAM=RAS*60D0 
      RD=RAM*60D0  

*  ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED
      A=NINT(RAS*R2AS*ABS(ANGLE))

*  SEPARATE INTO FIELDS
      AD=INT(A/RD) 
      A=A-AD*RD
      AM=INT(A/RAM)
      A=A-AM*RAM
      AS=INT(A/RAS)
      AF=A-AS*RAS


*  RETURN RESULTS  
      IDMSF(1)=INT(AD) 
      IDMSF(2)=INT(AM) 
      IDMSF(3)=INT(AS) 
      IDMSF(4)=INT(AF) 

      END  
      SUBROUTINE DR2TF(NDP,ANGLE,SIGN,IHMSF)
*  
*  
*  CONVERT ANGLE FOR OUTPUT, EXPRESSED AS TIME
*  -------------------------------------------
*  
*  
*  GIVEN:  
*     NDP           NUMBER OF DECIMAL PLACES OF SECONDS
*     ANGLE         ANGLE IN RADIANS (DOUBLE PRECISION)
*  
*  RETURNED:
*     SIGN          '+' OR '-'
*     IHMSF         HOURS, MINUTES, SECONDS, FRACTION
*  

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      INTEGER NDP  
      DOUBLE PRECISION ANGLE
      CHARACTER*1 SIGN 
      INTEGER IHMSF(4) 

*  RADIANS TO SECONDS  
      DOUBLE PRECISION R2S
      PARAMETER (R2S=13750.987083139757D0)

*  HANDLE SIGN 
      SIGN='+' 
      IF(ANGLE.LT.0D0) SIGN='-'

*  FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE
      RS=REAL(10**MAX(NDP,0))
      RM=RS*60D0
      RH=RM*60D0

*  ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED
      A=NINT(RS*R2S*ABS(ANGLE))

*  SEPARATE INTO FIELDS
      AH=INT(A/RH) 
      A=A-AH*RH
      AM=INT(A/RM) 
      A=A-AM*RM
      AS=INT(A/RS) 
      AF=A-AS*RS

*  RETURN RESULTS  
      IHMSF(1)=INT(AH) 
      IHMSF(2)=INT(AM) 
      IHMSF(3)=INT(AS) 
      IHMSF(4)=INT(AF) 

      END  

