*************************************************************************
*     
*     Subroutine tricircle.f
*     
*     - find the best circle passing through npoints: compute the circle
*     passing through every combination of 3 of them and average the 
*     resulting centres and radii
*     
*     output variables:
*     - angle(npoints)                      angle of the tangent in the input points, in degrees, range -90..+90
*     - residual(npoints)                   residuals
*     - chi                                 sum of squared residuals
*     - xc,zc,radius                        circle parameters
*     - eflag                               error flag
*     
*     to be called inside ./fitxy.f
*     
*************************************************************************


      subroutine tricircle(npoints,dep,indep,angle,residual,chi
     $     ,xc,zc,radius,eflag)


c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      integer npoints           !fit number of points
      real dep(npoints),indep(npoints) !dependent and independent variables

c      real angle(npoints)       !angle between the tangent line in the input points and 
c                                ! the independent variable axis
c      real residual(npoints)    !residuals
c      real chi                  !sum of squared residuals
c      real xc,zc,radius         !circle parameters   
      double precision angle(npoints)       !angle between the tangent line in the input points and 
                                ! the independent variable axis
      double precision residual(npoints)    !residuals EM GCC4.7
      double precision chi              !sum of squared residuals EM GCC4.7
      double precision xc,zc,radius         !circle parameters   EM GCC4.7
      
      integer eflag             !error flag =1 if the procedure fails

c     integer nloops            !number of combinations of 3 points out of npoints
      integer index(3)          !indexes of the 3 chosen points

      parameter(scale=1000.)
      double precision z(3),x(3),unit(3),zzxx(3) !temp variables
      double precision a(3,3),d(3,3),e(3,3),f(3,3)
      double precision ir(3)
      double precision deta,detd,dete,detf !determinants
      integer ifail             !=-1 if singular matrix error, =0 if not singular
      integer jfail             !=0 if determinant can be evaluated, =-1 if determinat is probably too small, =+1 if too large

      integer ibig ! EM GCC4.7
c      parameter (ibig=1.e4)      !just a number greater than C(npoints,3)
      parameter (ibig=10000)      !just a number greater than C(npoints,3)   EM GCC 4.7
      double precision xxc(ibig),zzc(ibig),rrr(ibig) !centres and radii to be averaged

      double precision tmp1,tmp2,tmp(npoints)    !temp variables

      logical DEBUG

      real pigr                 !3.1415...
      pigr=ACOS(-1.)

      
      DEBUG = .false.
      if(eflag.eq.1)DEBUG = .true.

      eflag = 0


c------------------------------------------------------------------------
c     choose 3 points out of npoints
c------------------------------------------------------------------------
c     nloops = fact(npoints) / (fact(3) * fact(npoints-3))
      
      k=0
      do i1=1,npoints-2
        index(1)=i1
        do i2=i1+1,npoints-1
          index(2)=i2
          do i3=i2+1,npoints
            index(3)=i3

            k=k+1               !number of combinations
c$$$  print*,' '          !???
c$$$  print*,'k =',k,' index =',index

c------------------------------------------------------------------------
c     build temp vectors
c------------------------------------------------------------------------
            do i=1,3
              z(i)=indep(index(i))/scale !to avoid too big numbers in matrix computation
              x(i)=dep(index(i))/scale
              unit(i)=1.
              zzxx(i)=z(i)**2.+x(i)**2.
            enddo
c$$$  print*,'z =',z,' x =',x !???
c$$$  print*,'unit =',unit,' zzxx =',zzxx

c------------------------------------------------------------------------
c     build the matrixes
c------------------------------------------------------------------------
            do i=1,3
              a(i,1)=z(i)       !A has (z x 1) as columns
              a(i,2)=x(i)
              a(i,3)=unit(i)
              d(i,1)=zzxx(i)    !D has (zzxx x 1) as columns
              d(i,2)=x(i)
              d(i,3)=unit(i)
              e(i,1)=zzxx(i)    !E has (zzxx z 1) as columns
              e(i,2)=z(i)
              e(i,3)=unit(i)
              f(i,1)=zzxx(i)    !F has (zzxx z x) as columns
              f(i,2)=z(i)
              f(i,3)=x(i)
            enddo

c$$$  print*,'matrix A:'  !???
c$$$  do i=1,3
c$$$  print*,(a(i,j),j=1,3)
c$$$  enddo
c$$$  print*,'matrix D:'  !???
c$$$  do i=1,3
c$$$  print*,(d(i,j),j=1,3)
c$$$  enddo
c$$$  print*,'matrix E:'  !???
c$$$  do i=1,3
c$$$  print*,(e(i,j),j=1,3)
c$$$  enddo
c$$$  print*,'matrix F:'  !???
c$$$  do i=1,3
c$$$  print*,(f(i,j),j=1,3)
c$$$  enddo

c------------------------------------------------------------------------
c     compute the determinants of A, D, E and F matrixes
c     using DFACT (http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/f011/top.html)
c------------------------------------------------------------------------

            ifail=0
            jfail=0
            call DFACT(3,a,3,ir,ifail,deta,jfail)
            if(ifail.eq.-1) then
              if(DEBUG)then
                 print*,'tricircle: ERROR: singular matrix A:'
                 do i=1,3
                    print*,(a(i,j),j=1,3)
                 enddo
              endif
              eflag=1
            endif
            if(jfail.eq.-1) then
              if(DEBUG)then
                 print*
     $           ,'tricircle: ERROR: matrix A: determinant too small?'
                 do i=1,3
                    print*,(d(i,j),j=1,3)
                 enddo
              endif
              eflag=1
            elseif(jfail.eq.1) then
              if(DEBUG)then
                 print*
     $             ,'tricircle: ERROR: matrix A: determinant too large?'
                 do i=1,3
                    print*,(d(i,j),j=1,3)
                 enddo
              endif
              eflag=1
            endif
            
            ifail=0
            jfail=0
            call DFACT(3,d,3,ir,ifail,detd,jfail)
            if(ifail.eq.-1) then
              if(DEBUG)then
                 print*,'tricircle: ERROR: singular matrix D:'
                 do i=1,3
                    print*,(d(i,j),j=1,3)
                 enddo
              endif
              eflag=1
           endif
           if(jfail.eq.-1) then
              if(DEBUG)then
                 print*
     $      ,'tricircle: ERROR: matrix D: determinant too small?'
                 do i=1,3
                    print*,(d(i,j),j=1,3)
                 enddo
              endif
              eflag=1
            elseif(jfail.eq.1) then
              if(DEBUG)then
              print*
     $             ,'tricircle: ERROR: matrix D: determinant too large?'
              do i=1,3
                print*,(d(i,j),j=1,3)
              enddo
              endif
              eflag=1
            endif

            ifail=0
            jfail=0
            call DFACT(3,e,3,ir,ifail,dete,jfail)
            if(ifail.eq.-1) then
              if(DEBUG)then
              print*,'tricircle: ERROR: singular matrix E:'
              do i=1,3
                print*,(e(i,j),j=1,3)
              enddo
              endif
              eflag=1
           endif
            if(jfail.eq.-1) then
              if(DEBUG)then
              print*
     $             ,'tricircle: ERROR: matrix E: determinant too small?'
              do i=1,3
                print*,(e(i,j),j=1,3)
              enddo
              endif
              eflag=1
            elseif(jfail.eq.1) then
              if(DEBUG)then
              print*
     $             ,'tricircle: ERROR: matrix E: determinant too large?'
              do i=1,3
                print*,(e(i,j),j=1,3)
              enddo
              endif
              eflag=1
            endif

            ifail=0
            jfail=0
            call DFACT(3,f,3,ir,ifail,detf,jfail)
c     ifail=-1            !???
            if(ifail.eq.-1) then
              if(DEBUG)then
              print*,'tricircle: ERROR: singular matrix F:'
              do i=1,3
                print*,(f(i,j),j=1,3)
              enddo
              endif
              eflag=1
            endif
            if(jfail.eq.-1) then
              if(DEBUG)then
              print*
     $             ,'tricircle: ERROR: matrix F: determinant too small?'
              do i=1,3
                print*,(f(i,j),j=1,3)
              enddo
              endif
              eflag=1
            elseif(jfail.eq.1) then
              if(DEBUG)then
              print*
     $             ,'tricircle: ERROR: matrix F: determinant too large?'
              do i=1,3
                print*,(f(i,j),j=1,3)
              enddo
              endif
              eflag=1
            endif

c------------------------------------------------------------------------
c     compute the centre and radius
c------------------------------------------------------------------------
            detd=-detd
            detf=-detf

c            xxc(k)=-detd/(2.*deta)
c            zzc(k)=-dete/(2.*deta)
            xxc(k)=-dete/(2.*deta)
            zzc(k)=-detd/(2.*deta)
            rrr(k)=SQRT((detd**2+dete**2)/(4.*deta**2.)-detf/deta)

c$$$  write(30,*) xxc(k)*scale !???
c$$$  write(31,*) zzc(k)*scale !???
c$$$  write(32,*) rrr(k)*scale !???
c$$$  print*,'xxc =',xxc(k)*scale,' zzc =',zzc(k)*scale
c$$$  $           ,' rrr =',rrr(k)*scale !???

          enddo                 !index loops
        enddo
      enddo


c------------------------------------------------------------------------
c     averages the centres and the radii
c------------------------------------------------------------------------
      xc=0.
      zc=0.
      radius=0.
      do i=1,k
        xc=xc+xxc(i)
        zc=zc+zzc(i)
        radius=radius+rrr(i)
      enddo
      xc=xc/k * scale           !back to micrometers
      zc=zc/k * scale
      radius=radius/k * scale

c$$$  c------------------------------------------------------------------------
c$$$  c     check for small radius...!???
c$$$  c------------------------------------------------------------------------
c$$$  num=200
c$$$  height=ABS(indep(1)-indep(6))
c$$$  if(radius.lt.(num*height)) then
c$$$  xc=0.
c$$$  zc=0.
c$$$  radius=0.
c$$$  print*,'tricircle: ERROR: bad circle'
c$$$  print*,'radius' ,radius,' < ', num,' x',height
c$$$  c$$$  print*,dep              !???
c$$$  c$$$  print*,indep             !???
c$$$  eflag=1
c$$$  endif


c------------------------------------------------------------------------
c     computes residuals and chi-squared
c------------------------------------------------------------------------
      chi=0.
      
c     print*,xc,zc,radius !???
      do i=1,npoints
        tmp1 = SQRT(radius**2.-(indep(i)-zc)**2.)
        tmp2 = dep(i)-xc
        if(ABS(tmp2-tmp1).le.ABS(tmp2+tmp1)) then !it chooses the right sign
          tmp(i)=tmp1           !according to residuals
        else
          tmp(i)=-tmp1
        endif
        residual(i)=tmp2 - tmp(i)
        chi=chi + residual(i)**2.
c     print*,dep(i)                 !???
c     print*,indep(i)              !???
c     print*,tmp1,tmp2,tmp(i),residual(i) !???
      enddo

c------------------------------------------------------------------------
c     it computes the angle between the tangent to the circumference and the
c     independent variable axis
c------------------------------------------------------------------------
      do i=1,npoints
        angle(i)=(zc-indep(i)) / tmp(i)
        angle(i)=ATAN(angle(i)) !-pi/2 <= angle <= pi/2
        angle(i)=angle(i)/pigr*180.
      enddo
      
      return
      end










c------------------------------------------------------------------------
c------------------------------------------------------------------------
c------------------------------------------------------------------------


c$$$c------------------------------------------------------------------------
c$$$c     Function to find the factorial value
c$$$c------------------------------------------------------------------------
c$$$
c$$$c     from http://www.digitalcoding.com/programming/fortran/tutorial/ftute10.htm
c$$$
c$$$
c$$$      FUNCTION FACT(N)
c$$$      FACT=1
c$$$      DO 10 J=2,N
c$$$         FACT=FACT*J
c$$$10    CONTINUE
c$$$      RETURN
c$$$      END
c$$$
c$$$c------------------------------------------------------------------------