	dimension Aabip(2000,20),Aabipr(2000,20),
     &            Aphilp(2000,20),Aphilpr(2000,20)
	dimension arraya(20),arrayar(20),arrayp(20),arraypr(20)
	dimension difa(20),difar(20),difp(20),difpr(20),
     &            indxa(20),indxar(20),indxp(20),indxpr(20)
	dimension periods(2000)
C Consider 2000 different periods, 1.2 to 80 days.
C But optimized for P=8-15 days
C periods=(20-1.0)/500.=0.038
C start with one period
C let fobs be the shortest interval: 0.25 days in the original
C let zobs be the longest period we care about; 12 in the original
C let nobs be the number of observations=15 in the original
	nobs=15
	fobs=0.50
	zobs=30.
	nobs1=nobs-1
C Attempt to find a good factor:
	fact=goodfact(nobs,fobs,zobs)
Cfact=1.175
Cfact=(zobs/fobs)**(1./nobs1)
	write(6,901) fobs,fact,nobs,zobs
901	format('Gaps=',f6.2,' times',f6.3,' **n,n=1,',i5,' to cover',f6.2)
	do 100 i=1,2000
	P=1+(i-1)*0.038
	periods(i)=P
C OK, now for this period....pick a random starting time:
	a=rand(0)
Cwrite(6,931) a
931	format('Random value=',f6.3)
C I'm going to assume now that a is between 0 and 1.0
C Now let's observe the star using two different methods, over nobs points:
	abi=0.
	do 105 j=1,nobs
	if(j.eq.1.) then
	val=0.0	
	else
	val=fobs*fact**(j-2)
	end if
	abi=abi+val
Cwrite(6,106) abi,val
106	format('Next abi time:',f6.2,' interval=',f6.2)
	abip=a+(abi/P)
	Aabip(i,j)=abip-int(abip)
	b=rand(0)
	abir=abi+(b-0.5)*0.2
	abipr=a+(abir/P)
	Aabipr(i,j)=abipr-int(abipr)
C OK, now figure out the Phil-phase
	phil=j-1
	philp=a+phil/P
	Aphilp(i,j)=philp-int(philp)
911	format('The Phil phase is',5f6.2)
	philr=phil+(b-0.5)*0.2
	philpr=a+philr/P
	Aphilpr(i,j)=philpr-int(philpr)
105	continue
100	continue
C OK, time for some statistics:
	do 200 i=1,2000
	suma=0.
	suma2=0.
	sumar=0.
	sumar2=0.
	sump=0.
	sump2=0.
	sumpr=0.
	sumpr2=0.
	do 205 j=1,nobs
	arraya(j)=Aabip(i,j)
	arrayar(j)=Aabipr(i,j)
	arrayp(j)=Aphilp(i,j)
	arraypr(j)=Aphilpr(i,j)
205	continue
C OK, at this point we want to sort the various arrays....
	call sort(nobs,arraya,indxa)
	call sort(nobs,arrayar,indxar)
	call sort(nobs,arrayp,indxp)
	call sort(nobs,arraypr,indxpr)
C Now figure out the differences in phases between successive points:
	do 206 j=1,nobs1
5	format(i3,8f8.2)
	difa(j)=arraya(indxa(j+1))-arraya(indxa(j))
	difar(j)=arrayar(indxar(j+1))-arrayar(indxar(j))
	difp(j)=arrayp(indxp(j+1))-arrayp(indxp(j))
	difpr(j)=arraypr(indxpr(j+1))-arraypr(indxpr(j))
206	continue
	difa(nobs)=arraya(indxa(1))+1.0-arraya(indxa(nobs))
	difar(nobs)=arrayar(indxar(1))+1.0-arrayar(indxar(nobs))
	difp(nobs)=arrayp(indxp(1))+1.0-arrayp(indxp(nobs))
	difpr(nobs)=arraypr(indxpr(1))+1.0-arraypr(indxpr(nobs))
	sa=1.
	sar=1.
	sp=1.
	spr=1.
	ba=-1
	bar=-1
	bp=-1
	bpr=-1
	do 210 j=1,nobs
	suma=suma+difa(j)
	if(difa(j).gt.ba) ba=difa(j)
	if(difa(j).lt.sa) sa=difa(j)
	sumar=sumar+difar(j)
	if(difar(j).gt.bar) bar=difar(j)
	if(difar(j).lt.sar) sar=difar(j)
	sump=sump+difp(j)
	if(difp(j).gt.bp) bp=difp(j)
	if(difp(j).lt.sp) sp=difp(j)
	sumpr=sumpr+difpr(j) 
	if(difpr(j).gt.bpr) bpr=difpr(j)
	if(difpr(j).lt.spr) spr=difpr(j)
210	continue
	avea=suma/float(nobs)
	avear=sumar/float(nobs)
	avep=sump/float(nobs)
	avepr=sumpr/float(nobs)
	do 211 j=1,nobs
	suma2=suma2+(difa(j)-avea)**2
	sumar2=sumar2+(difar(j)-avear)**2
	sump2=sump2+(difp(j)-avep)**2
	sumpr2=sumpr2+(difpr(j)-avepr)**2
211	continue
	siga=sqrt(suma2/nobs)
	sigar=sqrt(sumar2/nobs)
	sigp=sqrt(sump2/nobs)
	sigpr=sqrt(sumpr2/nobs)
	write(3,301) periods(i),
     &  avea,siga,sa,ba,
     &  avear,sigar,sar,bar,
     &  avep,sigp,sp,bp,
     &  avepr,sigpr,spr,bpr
301     format('P=',f7.3,' Abi=',4f6.3,' Abir=',4f6.3,' Phil=',4f6.3,
     &         ' Philr=',4f6.3)

Cwrite(6,301) periods(i),
C    &  avea,siga,sa,ba,
C    &  avear,sigar,sar,bar,
C    &  avep,sigp,sp,bp,
C    &  avepr,sigpr,spr,bpr
200	continue
	stop
	end
        subroutine sort(n,arrin,indx)
        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
	function goodfact(nobs,fobs,zobs)
        big=(zobs/fobs)**(1./(nobs-1.))
        small=fobs
        db=resid(fobs,zobs,nobs,big)
        if(db.lt.0.) stop 'I am confused big'
        ds=resid(fobs,zobs,nobs,small)
        if(ds.gt.0.) stop 'I am confused small'
        factnew=small
        value=.0001
1       factnew=factnew+value
        dnew=resid(fobs,zobs,nobs,factnew)
C   write(6,*) factnew,dnew
        if(dnew.ge.0.) then
	goodfact=factnew
	return
	end if
        go to 1
        end
        function resid(fobs,zobs,nobs,fact)
        sum=0.
        do 1 i=1,nobs
        sum=sum+fobs*fact**(i-1)
1       continue
        resid=sum-zobs
        return
        end
