C Takes "kfiles" worth of stellar photometry, and looks for the difference C in pixels... C assumes that the same stars are in the same order in each file, and allows C finds the shifts in the xy positions implicit double precision (a-h,o-z) parameter(kfiles=12) parameter(maxstars=200000) character*1 sign CC CC CC dimension amag(maxstars,kfiles),sigma(maxstars,kfiles), & sky(maxstars,kfiles),istat(maxstars),number(kfiles), & xk(maxstars,kfiles),yk(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 C set ifixum to 0 in order to IGNORE any photometric offsets ifixum=0 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,*,end=15) xa,ya,am,sig,sk istar=istar+1 xk(istar,nfiles)=xa yk(istar,nfiles)=ya 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)=xk(i,j)-xk(i,1) ddec(i)=(yk(i,j)-yk(i,1)) dmag(i)=amag(i,j)-amag(i,1) dsig(i)=sigma(i,j)+sigma(i,1) ratio(i)=dmag(i)/dsig(i) 26 continue call median (dra,num,ara) call median (ddec,num,adec) write(6,24) j,ara,adec 24 format(i3,'Median: x', f7.3,' y') 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(20000) 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'