Re: Welch-Stetson technique, Mark IV plans, and software

[Doug Welch <welch@physics.mcmaster.ca>]

Hi Mike,

That is a pretty good summary of how the W-S technique works.
It can also be extended to any number of "observations per
epoch" fairly easily. (Here, an "epoch" is the interval over 
which you expect the photometric variation of a star to be 
correlated.) Below is some code which calculates the statistic
for multiple observations per epoch. (It is FORTRAN - don't
hate me. I love C, I really do!!) The file enum contains the
groupings, but this, too, could be programmed.

Another thing you can do with the variability index is to
pair residuals from similar brightness but *different* stars.
Since the resulting index cannot be do to real variations, it
is a check on systematic errors which mimic stellar variability
which might be a function of flat-fielding, seeing, crowding,
time since last election, etc.

Cheers,
Doug

	implicit real*4 (a-h,o-z)

	parameter(iepmax=8, nepoch=13)

	character*30  filvres
	character*30  filstat
	integer       nf(iepmax)
	real          dv(nepoch)
	real          del(7)
	real          wt(7)

	wt(1)=0.0
	wt(2)=1.0	
	wt(3)=2.0/3.0
	wt(4)=3.0/6.0
	wt(5)=4.0/10.0
	wt(6)=5.0/15.0
	wt(7)=6.0/21.0

	print 100
 100	format('> Input residual file : ',$)
	read(5,*) filvres

	print 120
 120	format('>Output stats    file : ',$)
	read(5,*) filstat

	open(unit=3, file='enum', status='old')
	do i=1,iepmax
		read(3,*) ii,nf(i)
	end do
	close(unit=3)

	open(unit=1, file=filvres, status='old')
	open(unit=2, file=filstat, status='new')

c --- read in all the residuals

	ii=0
 200	ii=ii+1
        read(1,*,end=899) id,xc,yc,vavg,evavg,(dv(i), i=1,nepoch)

c	print 210, id,xc,yc,vavg,(dv(i), i=1,nepoch)
c210	format(i5,2f9.3,1f9.4,<nepoch>f10.3)

c --- for each epoch

	statvar=0.
	nstat0=0
	index=0
	do i=1,iepmax

c ... figure out how many useful residuals at this epoch
c ... and store in an auxiliary array del
		
	  nuse=0
	  do iii=1,nf(i)
	    index=index+1		
   	    if (dv(index).ne.100.000) then
		nuse=nuse+1
		del(nuse)=dv(index)
	    end if
	  end do

c ... for each useful epoch calculate product of all independent
c ... pairs of residuals, sum, then weight appropriately and
c ... add to final variability index statistic

	  statepo=0.0
	  if (nuse.eq.0.or.nuse.eq.1) then
	    nstat0=nstat0+1
	  else
	    na=nuse
	    nb=nuse-1
	    do j=1,nuse-1
	      do k=j+1,nuse
c		statepo=statepo+del(j)*del(k)
	        if (del(j).ne.0.0.and.del(k).ne.0.0) then
		  sign=del(j)*del(k)/abs(del(j)*del(k))
		  statepo=statepo+sign*sqrt(abs(del(j)*del(k)))
	        end if
	      end do
            end do
	    statvar=statvar+wt(nuse)*statepo
	  end if

	end do

c ... once statistic has been calculated, normalize for number
c ... of distinguishable epochs, and save as -1000.0 if no useful
c ,., statistic can be formed

	ncor=iepmax-nstat0
	if (ncor.ne.0) then
	  statvar=statvar/sqrt(real(ncor*(ncor-1)))
	else
	  statvar=-1000.0
	end if

	write(2,300) id,xc,yc,vavg,evavg,statvar,ncor
 300	format(i5,2f9.3,2f9.4,f9.3,i3)

c ... write output and then go to next star

	goto 200

 899	print *, ii-1,' stars processed'

	close(unit=2)
	close(unit=1)

	end