anchor.diff:  c     t    = array {t(i)}, containing the time values of the time series
   c     x    = array {x(i)}, containing the data values of the time series
   c     u    = temporal/work/dummy array, must be dimensioned in the
  -c            calling program in the same way as  {t(i)}
  +c     calling program in the same way as  {t(i)}
   c     v    = the same as  {u(i)}
   c     nf   = number of frequency points in which the spectrum is computed
   c     fmin = minimum frequency (MUST be > 0)
  
positive.code: c******************************
  c
        do 100 jf=1,nf
        f0=fmin+df*dfloat(jf-1)
        p0=1.0d0/f0
  c
  c======================================================
  c     Compute folded time series with  *p0*  period
  c======================================================
  c
        do 101 j=1,nb
          y(j) = 0.0d0
        ibi(j) = 0
   101  continue
  c
        do 102 i=1,n
        ph     = u(i)*f0
        ph     = ph-idint(ph)
        j      = 1 + idint(nb*ph)
        ibi(j) = ibi(j) + 1
          y(j) =   y(j) + v(i)
   102  continue
  c
  c-----------------------------------------------
  c     Extend the arrays  ibi()  and  y() beyond
  c     nb   by  wrapping
  c
        do 104 j=nb1,nbkma
        jnb    = j-nb
        ibi(j) = ibi(jnb)
          y(j) =   y(jnb)
   104  continue
  c-----------------------------------------------
  c
  c===============================================
  c     Compute BLS statistics for this period
  c===============================================
  c
        power=0.0d0
  c
        do 1 i=1,nb
        s     = 0.0d0
        k     = 0
        kk    = 0
        nb2   = i+kma
        do 2 j=i,nb2
        k     = k+1
        kk    = kk+ibi(j)
        s     = s+y(j)
        if(k.lt.kmi) go to 2
        if(kk.lt.kkmi) go to 2
        rn1   = dfloat(kk)
        pow   = s*s/(rn1*(rn-rn1))
        if(pow.lt.power) go to 2
        power = pow
        jn1   = i
        jn2   = j
        rn3   = rn1
        s3    = s
   2    continue
   1    continue
  c
        power = dsqrt(power)
        p(jf) = power
  c
        if(power.lt.bpow) go to 100
        bpow  =  power
        in1   =  jn1
        in2   =  jn2
        qtran =  rn3/rn
        depth = -s3*rn/(rn3*(rn-rn3))
        bper  =  p0
  c
   100  continue
  c
  
negative.code: c
  c
        subroutine eebls(n,t,x,u,v,nf,fmin,df,nb,qmi,qma,
       &     p,bper,bpow,depth,qtran,in1,in2)
  c
  c------------------------------------------------------------------------
  c     >>>>>>>>>>>> This routine computes BLS spectrum <<<<<<<<<<<<<<
  c
  c     [ see Kovacs, Zucker & Mazeh 2002, A&A, Vol. 391, 369 ]
  c
  c     This is the slightly modified version of the original BLS routine
  c     by considering Edge Effect (EE) as suggested by
  