C Written Phil Massey, Lowell Obs.; made available 10/29/01
C Takes "kfiles" worth of stellar photometry, and looks for variables,
C assumes that the same stars are in the same order in each file, and allows
C finds the shifts in the photometry and the sky
C file names are "allmags"
C testing
	implicit double precision (a-h,o-z)
	parameter(kfiles=15)
	parameter(maxstars=100000)
	character*1 sign
CC
CC
CC
        dimension amag(maxstars,kfiles),sigma(maxstars,kfiles),
     &  sky(maxstars,kfiles),istat(maxstars),number(kfiles),
     &  ra(maxstars,kfiles),dec(maxstars,kfiles),
     &  avem(kfiles),asig(maxstars),skym(kfiles),
     &  dra(maxstars),ddec(maxstars),dmag(maxstars),
     &  dsig(maxstars),ratio(maxstars),xx(kfiles),
     &  x(maxstars),s(maxstars),corr(kfiles),skies(maxstars),
     &  aves(kfiles),skyma(kfiles),skysig(kfiles)
CC
CC
CC
	character*25 newfile
	character*25 namelist
C
*  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 set ifixum to 0 in order to IGNORE any photometric offsets
	ifixum=1
	open(unit=1,file='temp',status='old')
	read(1,2) namelist
	close(unit=1)
	open(unit=1,file=namelist,status='old')
	nfiles=0
C read everything into memory....again....
1	read(1,2,end=20) newfile
Cwrite(6,2) newfile
2	format(a25)
	alast=lnblnk(newfile)
	nfiles=nfiles+1
	if(nfiles.gt.kfiles) stop 'too many files'
	open(unit=2,file=newfile(1:alast)//'com',status='old')
	istar=0
9	read(2,11,end=15) ihr,imin,sec,sign,ideg,iamin,asec,
     &   am,sig,sk
11      format(i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2,
     &   1x,f6.3,1x,f6.3,1x,f11.5,1x,a10)         
Cwrite(6,11) ihr,imin,sec,sign,ideg,iamin,asec,
C&   am,sig,sk
       	istar=istar+1
	ra(istar,nfiles)=stor(ihr,imin,sec)
	if(sign.ne.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	dec(istar,nfiles)=astor(ideg,iamin,asec)
	amag(istar,nfiles)=am
	sigma(istar,nfiles)=sig
	sky(istar,nfiles)=sk
        go to 9                                  
15	number(nfiles)=istar
	close(unit=2)
	write(6,16) nfiles,istar
16	format('File=',i3,' Number=',i5)
	go to 1
CCCCCCCCCCCCCCCCCCCCCCCCC
20	close(unit=1)
	write(6,21) nfiles,(number(j),j=1,nfiles)
21	format(i5,30i5)
	do 22 j=2,nfiles
	if(number(j).ne.number(1)) stop 'not all the same number of stars'
22	continue
	num=number(1)
C OK, now check if we've marked the same stars in the same order:
	corr(1)=0.
	do 25 j=2,nfiles
	do 26 i=1,num
	dra(i)=(ra(i,j)-ra(i,1))/s2r
	ddec(i)=(dec(i,j)-dec(i,1))/as2r
	dmag(i)=amag(i,j)-amag(i,1)
	dsig(i)=sigma(i,j)+sigma(i,1)
	ratio(i)=dmag(i)/dsig(i)
Cwrite(6,23) j,i,dra(i),ddec(i),dmag(i),dsig(i),ratio(i)
23	format('F:',i3,' S:',i3,' Dra(s)=',f6.3,' Ddec(a)=',f5.2,
     &  ' Dmag=',f6.3,' Sig=',f5.3,' rat:',f5.1)
26	continue
	call median (dra,num,ara)
        call median (ddec,num,adec)
        call median (dmag,num,amagh)
        call xfit(dmag,dsig,num,1,aamag,aasigm,aasig)
	write(6,24) j,ara,adec,amagh,aamag,aasigm
24	format(i3,'Median: ra', f7.3,' dec', f7.3,' mag=',f7.3,
     &  ' Wt ave ',f7.3,' sigm=',f7.3)
	corr(j)=amagh
25	continue
C Now try to do things "right"
C for each star, average the errors:
	do 105 i=1,num
	do 100 j=1,nfiles
	xx(j)=sigma(i,j)
100	continue
	call xfit(xx,xx,nfiles,0,avesig,sigmam,sx)
	asig(i)=avesig
105	continue
C set up the status buffer
C start with the assumption that all stars are constant
	do 201 i=1,num
	istat(i)=0
201	continue
C Now test for variables, adopting a correction based upon the
C median mag difference with file 1:
	do 501 i=1,num
	do 501 j=1,nfiles-1
	do 501 k=j+1,nfiles
	fixum=corr(j)-corr(k)
	if(ifixum.eq.0) fixum=0.
	dif=abs(amag(i,j)-amag(i,k)-fixum)
	error=sigma(i,j)+sigma(i,k)
	if((dif.gt.3.0*error).and.(dif.gt.0.05)) then
Cif(dif.gt.0.1) then
	istat(i)=1
Cwrite(6,601) i,j,k,amag(i,j),amag(i,k),dif,error,sky(i,k)
590	format(3i5,'amags(2),corr,dif,error,sky',6f8.3)
	end if
501	continue
	ivar=0
	do 502 i=1,number(1)
	if(istat(i).eq.1) ivar=ivar+1
502	continue
	write(6,503) ivar
503	format('Pass 0: number of variables is',i5)
C
C Now find the average for each file
C 
	iflag=0
C
310	continue
	do 202 j=1,nfiles
	icount=0
	do 203 i=1,num
	if (istat(i).eq.0) then
	icount=icount+1
	x(icount)=amag(i,j)
	s(icount)=asig(i)
	skies(icount)=sky(i,j)
	end if
203	continue
	call xfit(x,s,icount,1,ave,sigmam,sx)
	call median(skies,icount,sk)
	call xfit(skies,skies,icount,0,ska,sigmms,sxs)
	avem(j)=ave
	aves(j)=sigmam
	skym(j)=sk
	skyma(j)=ska
	skysig(j)=sigmms
	write(6,204) icount,j,ave
204	format(2i5,f8.3)
202	continue
C First pass! let's see who's variable...
	ivar=0
	do 550 i=1,num
	istat(i)=0
550	continue
	do 215 i=1,num
	do 215 j=1,nfiles-1
	do 215 k=j+1,nfiles
	corr(k)=avem(k)-avem(1)
	corr(j)=avem(j)-avem(1)
	fixum=corr(j)-corr(k)
	if(ifixum.eq.0) fixum=0. 
	dif=abs(amag(i,j)-amag(i,k)-fixum)
	error=sigma(i,j)+sigma(i,k)
	if((dif.gt.3.0*error).and.(dif.gt.0.05)) then
Cif(dif.gt.0.1) then
	istat(i)=1
Cwrite(6,601) i,j,k,amag(i,j),amag(i,k),dif,error
	end if
601	format('Declaring a variable',3i5,5f8.3)
215	continue
	do 220 i=1,number(1)
	if(istat(i).eq.1) ivar=ivar+1
220	continue
	iflag=iflag+1
	write(6,221) iflag,ivar
221	format('Pass',i5,' Number of variables=',i5)
	if(iflag.lt.5) go to 310
	open(unit=2,file='temp2',status='unknown')
	open(unit=1,file=namelist,status='old')
	do 229 j=1,nfiles
	read(1,2) newfile
C The way imcombine works, is that it will multiply by the number
C (unless it's "mode"), so we really want the recipical.  I.e.,
C if a star has 1000 counts in frame 1, and 1500 in frame 2, factor=1000/1500.
C
	factor=10.**((avem(j)-avem(1))/2.5)
	factsig=(10.**-(aves(j)/2.5)-1.)
C Similarly, imcombine will ADD the factor before scaling and combining.
C So, if frame 1 had 10 counts and frame 2 has 100 counts, msczero should
C be 0 and -90.
	write(6,226) newfile,factor,factsig,-skym(j),-skyma(j),skysig(j)
	write(2,225) newfile,factor,-skym(j)
226	format(a25,f10.3,f6.3,2f10.2,f6.2)
225	format(a25,f10.3,f10.2)
229	continue
	close(unit=2)
	close(unit=1)
	stop
	end
	subroutine xfit(x,sigmax,npts,mode,xmean,sigmam,sigma)
	implicit double precision (a-h,o-z)
	double precision sum,sumx,weight,free
	dimension x(*),sigmax(*)
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 median(data,n,amed)
	implicit double precision(a-h,o-z)
	dimension data(*),indx(100000)
	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)
	implicit double precision(a-h,o-z)
        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
        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'

