	program mkinvert
C rewritten 12 April 2002
	implicit double precision (a-h,o-z)
	parameter(maxstars=250000)
	parameter(maxfiles=5)
	integer irvec1(4),idvec1(4)
	character*1 sign
	character*9 fileout
	character*1 filts(maxfiles)
	character*13 filen
	dimension ra(maxstars,maxfiles),dec(maxstars,maxfiles),
     &  am(maxstars,maxfiles),ame(maxstars,maxfiles),
     &  niter(maxstars,maxfiles), 
     &  istat(maxstars,maxfiles),nstars(maxfiles) 
	dimension cterm(maxfiles),cterm2(maxfiles)
*  SECONDS TO RADIANS
      PARAMETER (S2R=0.72722052166430398D-4)
*  ARC SECONDS TO RADIANS
      PARAMETER (AS2R=0.48481368110953599D-5)
C
*  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)

C Distance to be considered the same star (in arcsec)
	data hit/0.75/
	data filts/'U','B','V','R','I'/
C Now read EVERYTHING into memory
C for every chip
	do 1000 ichip=1,8
	write(6,1010) ichip
1010	format('Working with chip:',i3)
C For this chip do a lot of things:
	do 50 icolor=1,maxfiles
	write(6,1011) icolor
1011	format('And color:',i5)
	write(filen,3) filts(icolor),ichip
3	format(a1,'im',i1,'merge.dat')
	write(6,2) filen
2	format('Now looking at:',a13)
C let's see if the file exists!
	istar=0
	nstars(icolor)=0
	open(unit=1,file=filen,status='old',err=6)
	go to 5
6	write(6,7)
7	format('File does not exist.')
C allows us to be missing a filter...
	go to 50
5	continue
C get the color terms:
	write(6,1005) filen(1:4)//'colorterms'
1005	format('Attempting to open file:',a14)
	open(unit=2,file=filen(1:4)//'colorterms',status='old')
 	read(2,*) c,c2
	write(6,1006) c,c2
1006	format('Extracted color terms',2f9.3)
	close(unit=2)
	cterm(icolor)=c
	cterm2(icolor)=c2
C OK, now read everything into memmory
11	read(1,12,end=99) ihr,imin,sec,sign,ideg,iamin,asec,ave,err,nit
12	format(i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2,2x,f7.3,2x,f7.3,
     &  2x,i2)
	istar=istar+1
	if(istar.gt.maxstars) stop "too many stars in mkinvert!"
	ra(istar,icolor)=stor(ihr,imin,sec)
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	dec(istar,icolor)=astor(ideg,iamin,asec)
	am(istar,icolor)=ave
	ame(istar,icolor)=err
	niter(istar,icolor)=nit
	go to 11
99	close(unit=1)
	nstars(icolor)=istar
	write(6,51) icolor,filts(icolor),nstars(icolor)
51	format('Color',i3,1x,a1,' had',i7,' stars')
C go get another color!
50	continue
	if(nstars(3).eq.0) go to 1000
	if(nstars(2).eq.0) go to 1000
	if(nstars(4).eq.0) go to 1000
	write(fileout,1003) ichip
1003	format('IM',i1,'.ubvri')
	write(6,1004) fileout
1004	format('Fileout:',1x,a9)
	open(unit=3,file=fileout,status='new')
C Ok, we now have everything read in!
C now star with V and go looking for matches!
	do 100 iv=1,nstars(3)
Cif((iv/100)*100.eq.iv) write(6,3001) iv
3001	format('Looking at v photometry for star',i5)
C does the star have a match in the B filter?
	  do 150 ib=1,nstars(2)
C         if((ib/100)*100.eq.ib) write(6,3002) ib
3002	format('Looking at b star',i5)
	  if(abs(dec(iv,3)-dec(ib,2)).lt.1.0e-5) then
	  call apa(ra(iv,3),dec(iv,3),ra(ib,2),dec(ib,2),sep,pa)
	  else
	  sep=1.
	  end if
	  dist=sep/as2r
	    if(dist.lt.hit) then
C We have a winner!
Cwrite(6,3003) iv,ib
3003	format('V star',i7,' matches B star',i7)
C Oh, boy!
	    be=ame(ib,2)
	    ve=ame(iv,3)
	    bmv=am(ib,2)-am(iv,3)
	    bmve=sqrt(be**2+ve**2)
	    iflag=0
C	write(6,3004) bmv,bmve,ve,be
3004		format('Beginning b-v, b-ve,ve,be:',4f8.2)
52	    v=am(iv,3)-cterm(3)*bmv
	    ve=sqrt(ame(iv,3)**2+(cterm(3)*bmve)**2)
	    b=am(ib,2)-cterm(2)*bmv
	    be=sqrt(ame(ib,2)**2+(cterm(2)*bmve)**2)
	    bmv=b-v
	    bmve=sqrt(ve**2+be**2)
C  write(6,3005) iflag,bmv,bmve,v,ve,b,be
3005	format('iflag,bmv,bmve,v,ve,b,be',i3,6f8.3)
	    iflag=iflag+1
	    if(iflag.lt.10) go to 52
C We have v, we have b
C Let's see if there is an R as well
	    do 155 ir=1,nstars(4)
	    if(abs(dec(iv,3)-dec(ir,4)).lt.1.0e-5) then
	    call apa(ra(iv,3),dec(iv,3),ra(ir,4),dec(ir,4),sep,pa)
	    else
	    sep=1.
	    end if
	    dist=sep/as2r
	    if(dist.lt.hit) then
C    write(6,3006) iv,ib,ir
3006	    format('We have a 3-way hit! iv,ib,ir:',3i8)
C We have another winner!
	      iflag=0
	      vmr=v-am(ir,4)
	      re=ame(ir,4)
              vmre=sqrt(ve**2+re**2)
C      write(6,3007) vmr,vmre
3007	      format('Begining vmr,vmre:',2f8.2)
C      write(6,4001) cterm(4)
4001	      format('cterm-4',f8.2)
53	      r=am(ir,4)-cterm(4)*vmr
              re=sqrt(ame(ir,4)**2+(cterm(4)*vmre)**2)
	      vmr=v-r
	      vmre=sqrt(ve**2+re**2)
	      iflag=iflag+1
C      write(6,3008) iflag,vmr,vmre
3008	      format('iflag,vmr,vmre',i5,2f8.2)
	      if(iflag.lt.10) go to 53
C Now we should see if there is a U or an I
              umb=99.999
	      umbe=99.999
	      if(nstars(1).eq.0) go to 260
	      do 160 iu=1,nstars(1)
	      if(abs(dec(iv,3)-dec(iu,1)).lt.1.0e-5) then
	      call apa(ra(iv,3),dec(iv,3),ra(iu,1),dec(iu,1),sep,pa)
	      else
	      sep=1.
              end if
              dist=sep/as2r
              if(dist.lt.hit) then
C      write(6,3009) iv,ib,ir,iu
3009	      format('Wow, a U match!',4i5)
              iflag=0
	      umb=am(iu,1)-b
              ue=ame(iu,1)
	      umbe=sqrt(be**2+ue**2)
C      write(6,3010) umb,umbe
3010	      format('Beginning U-B,U-Be',2f8.2)
54            continue
              if(umb.gt.0.) then
	      u=am(iu,1)-cterm(1)*umb
              ue=sqrt(ame(iu,1)**2+(cterm(1)*umbe)**2)
              else
              u=am(iu,1)-cterm2(1)*umb
              ue=sqrt(ame(iu,1)**2+(cterm2(1)*umbe)**2)
              end if
              umb=u-b
              umbe=sqrt(ue**2+be**2)
              iflag=iflag+1
C      write(6,3011) iflag,umb,umbe
3011	      format('iflag,umb,umbe',i5,2f8.2)
              if(iflag.lt.10) go to 54
	      go to 260
              end if
160           continue
C      write(6,3012) iv,ib,ir
3012	      format('There was no U match for iv,ib,ir:',3i5)
C See if there's an I
260           rmi=99.999
	      rmie=99.999
	      if(nstars(5).eq.0) go to 2000
	      do 360 ii=1,nstars(5)
	      if(abs(dec(iv,3)-dec(ii,5)).lt.1.0e-5) then
	      call apa(ra(iv,3),dec(iv,3),ra(ii,5),dec(ii,5),sep,pa)
	      else
              sep=1.
              end if
              dist=sep/as2r
              if(dist.lt.hit) then
C             write(6,30122) iv,ib,ir,ii
30122	      format('I match!:',4i7)
              iflag=0
              rmi=r-am(ii,5)
              aie=ame(ii,5)
	      rmie=sqrt(re**2+aie**2)
C             write(6,3013) rmi,rmie
3013	      format('Begining r-i,r-ie:',2f8.2)
55            continue
              ai=am(ii,5)-cterm(5)*rmi
              aie=(ame(ii,5)**2+(cterm(5)*rmie)**2)
              rmi=r-ai
              rmie=sqrt(re**2+aie**2)
              iflag=iflag+1
C      write(6,3014) iflag,rmi,rmie
3014	      format('iflag,r-i,r-ie',i5,2f8.2)
              if(iflag.lt.10) go to 55
C write it out!
              go to 2000
              end if
360           continue
C      write(6,3015) iv,ib,ir
3015	      format('Failed to find an I match',3i7)
C even if no I, go write it out
              go to 2000
C ends the r loop:
		end if
155		continue
C ends the b loop:
              end if
150	      continue
C forget it if there's no B
              go to 100
2000          call dr2tf(3,ra(iv,3),sign,irvec1)
              call dr2af(2,dec(iv,3),sign,idvec1)
              write(3,2001) irvec1,sign,idvec1,v,ve,
     &        bmv,bmve,umb,umbe,vmr,vmre,rmi,rmie
2001	format(i2.2,':',i2.2,':',i2.2,'.',i3.3,1x,a1,
     &         i2.2,':',i2.2,':',i2.2,'.',i2.2,
     &        10f7.3)
100	continue
	close(unit=3)
1000	continue
	stop
	end
              
 

	subroutine apa(rac,dcc,rao,dco,dist,pa)
c computes the distance (delt) and position (pa) along the
c sky FROM (rac,dec) to (rao,dco)
c written aug 14, 1987
C rewritten August 21, 1989, to do using "standard coords", rather
C than the old and peculiar (but possibly equivalent!) way
c plm
      implicit double precision (a-h,o-z)
        parameter (PI=3.141592653589793238462643)
C CONVERT TO PROPER COORDS
        CALL DS2TP(RAO,DCO,RAC,DCC,XI,ETA)
        pa=0.
        if(xi*eta.ne.0.)PA=DATAN2(XI,ETA)
        if(pa.le.0.) pa=pa+2.*pi
        DIST=DSQRT(XI**2+ETA**2)
        return
        end
	include '/tofu/home/massey/observing/astromsub.f'
