*************************************************************************
*
* 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
real angle(npoints) !angle between the tangent line in the input points and
! the independent variable axis
real residual(npoints) !residuals
real chi !sum of squared residuals
real xc,zc,radius !circle parameters
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
parameter (big=1.e4) !just a number greater than C(npoints,3)
double precision xxc(big),zzc(big),rrr(big) !centres and radii to be averaged
double precision tmp1,tmp2,tmp(npoints) !temp variables
real pigr !3.1415...
pigr=ACOS(-1.)
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
print*,'tricircle: ERROR: singular matrix A:'
do i=1,3
print*,(a(i,j),j=1,3)
enddo
eflag=1
endif
if(jfail.eq.-1) then
print*
$ ,'tricircle: ERROR: matrix A: determinant too small?'
do i=1,3
print*,(d(i,j),j=1,3)
enddo
eflag=1
elseif(jfail.eq.1) then
print*
$ ,'tricircle: ERROR: matrix A: determinant too large?'
do i=1,3
print*,(d(i,j),j=1,3)
enddo
eflag=1
endif
ifail=0
jfail=0
call DFACT(3,d,3,ir,ifail,detd,jfail)
if(ifail.eq.-1) then
print*,'tricircle: ERROR: singular matrix D:'
do i=1,3
print*,(d(i,j),j=1,3)
enddo
eflag=1
endif
if(jfail.eq.-1) then
print*
$ ,'tricircle: ERROR: matrix D: determinant too small?'
do i=1,3
print*,(d(i,j),j=1,3)
enddo
eflag=1
elseif(jfail.eq.1) then
print*
$ ,'tricircle: ERROR: matrix D: determinant too large?'
do i=1,3
print*,(d(i,j),j=1,3)
enddo
eflag=1
endif
ifail=0
jfail=0
call DFACT(3,e,3,ir,ifail,dete,jfail)
if(ifail.eq.-1) then
print*,'tricircle: ERROR: singular matrix E:'
do i=1,3
print*,(e(i,j),j=1,3)
enddo
eflag=1
endif
if(jfail.eq.-1) then
print*
$ ,'tricircle: ERROR: matrix E: determinant too small?'
do i=1,3
print*,(e(i,j),j=1,3)
enddo
eflag=1
elseif(jfail.eq.1) then
print*
$ ,'tricircle: ERROR: matrix E: determinant too large?'
do i=1,3
print*,(e(i,j),j=1,3)
enddo
eflag=1
endif
ifail=0
jfail=0
call DFACT(3,f,3,ir,ifail,detf,jfail)
c ifail=-1 !???
if(ifail.eq.-1) then
print*,'tricircle: ERROR: singular matrix F:'
do i=1,3
print*,(f(i,j),j=1,3)
enddo
eflag=1
endif
if(jfail.eq.-1) then
print*
$ ,'tricircle: ERROR: matrix F: determinant too small?'
do i=1,3
print*,(f(i,j),j=1,3)
enddo
eflag=1
elseif(jfail.eq.1) then
print*
$ ,'tricircle: ERROR: matrix F: determinant too large?'
do i=1,3
print*,(f(i,j),j=1,3)
enddo
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------------------------------------------------------------------------