c Detecting aperiodic photometric variability
c with the Friends of Friends (FoF) algorithm
c Author: M. Rozyczka, version April 2018
c
****************************************************************************
*   This program is free software; you can redistribute it and/or modify   *
*   it under the terms of the  GNU General  Public  License as published   *
*   by the Free  Software  Foundation;  either version 3 of the License,   *
*   or (at your option) any later version.                                 *
*   This  program  is  distributed  in the hope that it will be  useful,   *
*   but  WITHOUT  ANY  WARRANTY;   without  even  the  implied  warranty   *
*   of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU   *
*   General Public License for more details.                               *
****************************************************************************
c
c
c Details of the approach are given in Rozyczka et al. 2018, Acta Astr. 68,63
c Briefly, each examined light curve is split into clusters of points;
c the brightest cluster is identfied, and, depending on its location,
c on the (normalized time - normalized brightness) plane a variability 
c index vi, where [0 <= vi <= 1], is assigned to the curve. Curves with 
c highest indices are then inspected individually by eye. 
c                           CAUTION
c It is advisable to calibrate the code as described in the paper 
c before applying it to a full-scale search for variables. 
c
c
c parameters of the code:
c maxobs - maximum number of points in a light curve
c maxnei - maximum size of a cluster
c maxcl  - maximum number of clusters (per light curve)
c
c input:
c fof.control  - file with control parameters
c                npmin = minimum number of points per cluster
c                      i.e. aggregatess with fewer than npmin points 
c                      are not counted as clusters [default: 5]
c                eps   = clustering length (CL) parameter, where  
c                      CL = eps*(mean distance between points
c                      of the normalized light curve) [default: 0.03]
c                ntops = maximum number of brightest aggregates with 
c                      fewer than npmin points which will be stripped 
c                      from the light curve before the variability 
c                      index is calculated (ntops > 0 allows to get 
c                      rid of unphysically brigth points, originating
c                      e.g. from hot pixels)  [default: 1]
c fof.dat      - list of light curves to be examined
c light curves - record structure: HJD, magnitudo
c
c output: file fof.log 
c         header: list of control parameters from file fof.control
c         body  : sequential number, name of the light curve ,
c         mean distance between points of the normalized light curve,
c         clustering length, initial and final number of points,
c         number of seed clusters, total number of identified clusters,
c         final number of identified clusters, variability index,
c         number of brightest aggregates with fewer than npmin points 
c         which had been stripped from the light curve.
c
c   
      implicit double precision (a-h,o-z)
      parameter (maxobs = 3000, maxnei = 3000, maxcl =1500)
      dimension t(maxobs), to(maxobs), v(maxobs), vo(maxobs)
      dimension d(maxobs,maxobs)
      dimension clmean(maxcl), clmin(maxcl), clmax(maxcl)
      dimension ids(maxobs), nnei(maxobs), nei(maxobs, maxnei)
      dimension nncl(maxcl), ncl(maxobs,maxcl), iclmax(maxcl)
      dimension nncl1(maxcl), ncl1(maxobs,maxcl)
      logical, dimension (1:maxobs):: visited
      logical, dimension (1:maxcl,1:maxcl):: mergecl
      logical link_ij
      character(len=30) fname, fnam, line
c
c open file with names of light curves, output file, and file with control parameters
c
      open(10,file="fof.dat")
      open(20,file="fof.log")
      open(30,file="fof.control") 
c
      read(30,*) npmin
      read(30,*) eps0
      read(30,*) ntops
c      
      write(20,"(14x,a17,i3,a9,f6.3,a14,i3)") 
     &"Run with  npmin =",npmin,",  eps0 =",eps0,",  and ntops =",ntops
      write(20,*)" "
      write(20,"(14x,a4,4x,2(5x,a2),4x,a6,4x,a5,2x,a8,2x,a8)")
     &"file","md","cl","points","seeds","clusters","var_indx"
      istar = 0
c
  100 continue             !start loop over lightcurves
c
c clean arrays
c
      ids     = 0
      nei     = 0
      nnei    = 0
      ncl     = 0
      ncl1    = 0
      nncl    = 0
      nncl1   = 0
      clmean  = 0.d0
      clmin   = 0.d0
      d       = 0.d0
      t       = 0.d0
      v       = 0.d0
      mergecl = .false.
      visited = .false.
c
      read(10,"(a30)",end=500) fname
      open(11,file=trim(fname),err=501)
      istar = istar+1
c
c skip header if present
c
      do i = 1, 20
        read(11,*,end=501) line
        if (line(1:1) .ne. "#") goto 90 
      end do
   90 backspace(11)
c
c read in the next light curve
c
      do i = 1,maxobs
        read(11,*,end=101) to(i),v(i)
        vo(i) = v(i)
        if (to(i) .gt. 2d6) to(i) = to(i) - 2.45d6
      end do
      write(20,*)istar,trim(fname)," maxobs too small"
      goto 100
  101 nobs = i-1 
c
c sort time
c
  102 ll = 0
      do i = 1,nobs-1
        if (to(i+1).lt.to(i)) then
          aux = t(i)
          t(i) = t(i+1)
          t(i+1) = aux
          aux = v(i)
          v(i) = v(i+1)
          v(i+1) = aux
          aux = to(i)
          to(i) = to(i+1)
          to(i+1) = aux
          aux = vo(i)
          vo(i) = vo(i+1)
          vo(i+1) = aux
          ll = 1
        end if
      end do
      if (ll .eq. 1) goto 102
c
c compress time
c
      hjdn = 0
      hjd  = int(to(1))
      do i = 1,nobs
        if (int(to(i)).eq.hjd) then
          t(i) = to(i)-hjd+hjdn
        else
          hjd  = int(to(i))
          hjdn = hjdn+1
          t(i) = to(i)-hjd+hjdn
        end if
      end do
c
c sort V and remove points fainter than median V
c
  103 ll = 0
      do i = 1,nobs-1
        if (v(i+1).gt.v(i)) then
          aux = v(i)
          v(i) = v(i+1)
          v(i+1) = aux
          aux = t(i)
          t(i) = t(i+1)
          t(i+1) = aux
          aux = to(i)
          to(i) = to(i+1)
          to(i+1) = aux
          aux = vo(i)
          vo(i) = vo(i+1)
          vo(i+1) = aux
          ll = 1
        end if
      end do
      if (ll .eq. 1) goto 103
      ll = 0
      do i = nobs/2+2,nobs
        ll = ll+1
        t(ll) = t(i)
        v(ll) = v(i)
        ids(ll) = i
        t(i)  = 0
        v(i)  = 0
      end do
      nobs = ll
c
c normalize t and V 1vectors
c
      tmin = 1.d20
      vmin = 1.d20
      do i = 1,nobs
        if (t(i) .lt. tmin) tmin = t(i)
        if (v(i) .lt. vmin) vmin = v(i)
      end do
      tmax = -1.d20
      vmax = -1.d20
      do i = 1,nobs
        if (t(i) .gt. tmax) tmax = t(i)
        if (v(i) .gt. vmax) vmax = v(i)
      end do
      dt   = tmax - tmin
      dv   = vmax - vmin
      do i = 1,nobs
        t(i) = (t(i)-tmin)/dt
        v(i) = (vmax-v(i))/dv
      end do
c 
c find distances between points
c
      dd = 0.d0
      do i = 1,nobs
        tc = t(i)
        vc = v(i)
        do j = 1,nobs
          d(i,j) = sqrt((tc-t(j))**2+(vc-v(j))**2)
          dd = dd + d(i,j)
        end do
      end do
c
c find mean distance and clustering length
c
      dd = dd/(dble(nobs))**2
      eps = eps0*dd
c
c find neighbours for each (t,v) point 
c
      do i = 1,nobs
        ii = 0
        do j = 1,nobs
          if (d(i,j) .le. eps) then
            ii = ii+1
            if (ii .gt. maxnei) then 
              write(20,*)istar,trim(fname)," maxnei too small"
              goto 100
            end if
            nei(i,ii) = j  ! ids of neighbors to point i
          end if
        end do
        nnei(i) = ii       ! number of neighbors to point i
      end do
c
c define disjunct cluster seeds
c
      do i = 1, maxobs
        visited(i) = .false.
      end do
      ntotcl = 0
      do i = 1,nobs
        if (.not.visited(i)) then
          ntotcl = ntotcl+1
          icl    = 0
          do j = 1,nnei(i)
            if (.not.visited(nei(i,j))) then  
              icl = icl+1
              if (icl .gt. maxcl) then
                write(20,*)istar,trim(fname)," maxcl too small"
                goto 100
              end if
              ncl(ntotcl,icl) = nei(i,j)  !ids of unique members of ntotcl-th seed
              visited(nei(i,j)) = .true.
            end if
          end do   
          visited(i) = .true.
          nncl(ntotcl) = icl              !number of unique members of ntotcl-th seed
        end if
      end do
      mm1 = 0 
      do i = 1,ntotcl
        mm1 = mm1 + nncl(i)               !total number of seed members, must be equal to nobs
      end do
      if (mm1 .ne. nobs) stop "n(seed members) .ne. nobs"
c
c mark cluster seeds for merging
c
      do i = 1,ntotcl
        do j = 1,ntotcl 
          mergecl(i,j) = .false.
        end do
      end do
      do i = 1,ntotcl
        do j = i,ntotcl
          do i1 = 1,nncl(i)
            do j1 = 1,nncl(j)
              if (d(ncl(i,i1),ncl(j,j1)) .le. eps) mergecl(i,j) = .true.
            end do
          end do
        end do
      end do
      do i = ntotcl-1,1,-1
        do j = i+1, ntotcl
          link_ij = .false.
          do k = 1, ntotcl
            if (mergecl(i,k) .and. mergecl(j,k)) link_ij = .true.
          end do
          if (link_ij) then  
            do k = 1,ntotcl
              mergecl(i,k) = mergecl(i,k) .or. mergecl(j,k)
              mergecl(j,k) = .false.
            end do
          end if
        end do
      end do
c
      nm1 = 0
      nm2 = 0
      nm3 = 0 
      do i = 1,ntotcl
        kk = 0
        do j = 1,ntotcl
          if (mergecl(i,j)) then
            nm1 = nm1+nncl(j)
            nm2 = nm2+1
            kk = 1
          end if        
        end do
        if (kk.eq.1) nm3 = nm3+1
      end do
c
c merge cluster seeds
c   
      ntotcl1 = 0
      do i = 1,ntotcl
        if (mergecl(i,i)) then
          ntotcl1 = ntotcl1+1
          ll = 0
          do j = 1,ntotcl
            if (mergecl(i,j)) then 
              do k=1,nncl(j) 
                ll = ll+1
                if (ll .gt. maxcl) then
                  write(20,*)
     &            istar,trim(fname)," maxcl too small for merged cl"
                  goto 100
                end if
                ncl1(ntotcl1,ll) = ncl(j,k)
              end do
            end if
          end do
          nncl1(ntotcl1) = ll
        end if
      end do
c
c if there are top clusters with n < npmin remove them
c and rescale the v - array
c
      do i = 1,ntotcl1
        iclmax(i) = i
        clmax(i) = 0.d0
        do j = 1, nncl1(i)
          ij = ncl1(i,j)
          if (clmax(i) .lt. v(ij)) clmax(i) = v(ij) 
        end do
      end do
  202 ll = 0
      do i = 1,ntotcl1-1
        if (clmax(i+1).lt.clmax(i)) then
          aux = clmax(i)
          clmax(i) = clmax(i+1)
          clmax(i+1) = aux
          naux = iclmax(i)
          iclmax(i) = iclmax(i+1)
          iclmax(i+1) = naux
          ll = 1
        end if
      end do
      if (ll .eq. 1) goto 202
      nred = 0
      ntotcl2 = ntotcl1
      do i = ntotcl1, ntotcl1-ntops+1, -1
        if (nncl1(iclmax(i)) .ge. npmin) goto 203 
        ntotcl2 = ntotcl2 -1
        nred = nred + 1
      end do 
  203 ntotcl1 = ntotcl2
      teps = 1.d-20
      do i = 1,maxobs
        v(i) = v(i)/(clmax(ntotcl1)+teps)
      end do
c
c reject clusters of .le. npmin points
c
      ntotcl = 0
      do i = 1,ntotcl1
        if (nncl1(i) .ge. npmin) then 
          ntotcl = ntotcl+1
          nncl(ntotcl) = nncl1(i)
          do j = 1, nncl1(i)
            ncl(ntotcl,j) = ncl1(i,j)
          end do
        end if
      end do     
      iscl = 0
      do i = 1,ntotcl
        iscl = iscl + nncl(i)
      end do
c
c for all clusters find mean normalized V and minimum normalized V
c
      do i = 1,ntotcl
        clmean(i) = 0.d0
        clmin(i)  = 1.d1
        do j = 1,nncl(i)
          ij = ncl(i,j)
          clmean(i) = clmean(i) + v(ij)
          if(clmin(i) .gt. v(ij)) clmin(i) = v(ij) 
        end do
        clmean(i) = clmean(i)/nncl(i)
      end do
c
c find maximum mean normalized V and minimum normalized V among clusters
c
      clmaxmean = 0.d0
      clmaxmin  = 0.d0
      if (ntotcl.gt.0) then
        clmaxmean = maxval(clmean)
        clmaxmin  = maxval(clmin)
      end if
c
c printout: 
c nm1    = total number of merged points
c iscl   = total number of of final cluster members
c nm2    = total number of cluster seeds
c nm3    = number of all clusters
c ntotcl = number of final clusters (each with with more than npmin members)
c 
      var_index = clmaxmin
      if (clmaxmin .gt. 0.4) var_index = clmaxmean
      write(20
     &,"(i6,2x,a15,2f7.3,2i5,2x,i5,2x,i4,i5,2x,f7.3,i3)")
     &istar,trim(fname),dd,eps,nm1,iscl,nm2,nm3,ntotcl,var_index,nred
c
      goto 100
c
  500 stop "last lc processed"
  501 write(20,*)
     &istar,trim(fname)," error opening file"
      goto 100
      end
