src/F77/tricircle.f

*************************************************************************
*     
*     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 (ibig=10000)      !just a number greater than C(npoints,3) !EM GCC4.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+REAL(xxc(i)) !EM GCC4.7
        zc=zc+REAL(zzc(i)) !EM GCC4.7
        radius=radius+REAL(rrr(i)) !EM GCC4.7
      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)=REAL(tmp2 - tmp(i)) !EM GCC4.7
        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)=REAL((zc-indep(i)) / tmp(i)) !EM GCC4.7
        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------------------------------------------------------------------------
1	mocchiut	1.1	*************************************************************************
2			*
3			* Subroutine tricircle.f
4			*
5			* - find the best circle passing through npoints: compute the circle
6			* passing through every combination of 3 of them and average the
7			* resulting centres and radii
8			*
9			* output variables:
10			* - angle(npoints) angle of the tangent in the input points, in degrees, range -90..+90
11			* - residual(npoints) residuals
12			* - chi sum of squared residuals
13			* - xc,zc,radius circle parameters
14			* - eflag error flag
15			*
16			* to be called inside ./fitxy.f
17			*
18			*************************************************************************
19
20
21			subroutine tricircle(npoints,dep,indep,angle,residual,chi
22			$ ,xc,zc,radius,eflag)
23
24
25			c------------------------------------------------------------------------
26			c
27			c local variables
28			c
29			c------------------------------------------------------------------------
30
31			integer npoints !fit number of points
32			real dep(npoints),indep(npoints) !dependent and independent variables
33
34	mocchiut	1.6	real angle(npoints) !angle between the tangent line in the input points and
35	mocchiut	1.1	! the independent variable axis
36	mocchiut	1.6
37			real residual(npoints) !residuals
38			real chi !sum of squared residuals
39			real xc,zc,radius !circle parameters
40	mocchiut	1.1
41			integer eflag !error flag =1 if the procedure fails
42
43			c integer nloops !number of combinations of 3 points out of npoints
44			integer index(3) !indexes of the 3 chosen points
45
46			parameter(scale=1000.)
47			double precision z(3),x(3),unit(3),zzxx(3) !temp variables
48			double precision a(3,3),d(3,3),e(3,3),f(3,3)
49			double precision ir(3)
50			double precision deta,detd,dete,detf !determinants
51			integer ifail !=-1 if singular matrix error, =0 if not singular
52			integer jfail !=0 if determinant can be evaluated, =-1 if determinat is probably too small, =+1 if too large
53
54	mocchiut	1.6	parameter (ibig=10000) !just a number greater than C(npoints,3) !EM GCC4.7
55	mocchiut	1.4	double precision xxc(ibig),zzc(ibig),rrr(ibig) !centres and radii to be averaged
56	mocchiut	1.1
57			double precision tmp1,tmp2,tmp(npoints) !temp variables
58
59	pam-fi	1.2	logical DEBUG
60
61	mocchiut	1.1	real pigr !3.1415...
62			pigr=ACOS(-1.)
63
64	pam-fi	1.2
65			DEBUG = .false.
66			if(eflag.eq.1)DEBUG = .true.
67
68	mocchiut	1.1	eflag = 0
69
70
71			c------------------------------------------------------------------------
72			c choose 3 points out of npoints
73			c------------------------------------------------------------------------
74			c nloops = fact(npoints) / (fact(3) * fact(npoints-3))
75
76			k=0
77			do i1=1,npoints-2
78			index(1)=i1
79			do i2=i1+1,npoints-1
80			index(2)=i2
81			do i3=i2+1,npoints
82			index(3)=i3
83
84			k=k+1 !number of combinations
85			c$$$ print*,' ' !???
86			c$$$ print*,'k =',k,' index =',index
87
88			c------------------------------------------------------------------------
89			c build temp vectors
90			c------------------------------------------------------------------------
91			do i=1,3
92			z(i)=indep(index(i))/scale !to avoid too big numbers in matrix computation
93			x(i)=dep(index(i))/scale
94			unit(i)=1.
95			zzxx(i)=z(i)2.+x(i)2.
96			enddo
97			c$$$ print*,'z =',z,' x =',x !???
98			c$$$ print*,'unit =',unit,' zzxx =',zzxx
99
100			c------------------------------------------------------------------------
101			c build the matrixes
102			c------------------------------------------------------------------------
103			do i=1,3
104			a(i,1)=z(i) !A has (z x 1) as columns
105			a(i,2)=x(i)
106			a(i,3)=unit(i)
107			d(i,1)=zzxx(i) !D has (zzxx x 1) as columns
108			d(i,2)=x(i)
109			d(i,3)=unit(i)
110			e(i,1)=zzxx(i) !E has (zzxx z 1) as columns
111			e(i,2)=z(i)
112			e(i,3)=unit(i)
113			f(i,1)=zzxx(i) !F has (zzxx z x) as columns
114			f(i,2)=z(i)
115			f(i,3)=x(i)
116			enddo
117
118			c$$$ print*,'matrix A:' !???
119			c$$$ do i=1,3
120			c$$$ print*,(a(i,j),j=1,3)
121			c$$$ enddo
122			c$$$ print*,'matrix D:' !???
123			c$$$ do i=1,3
124			c$$$ print*,(d(i,j),j=1,3)
125			c$$$ enddo
126			c$$$ print*,'matrix E:' !???
127			c$$$ do i=1,3
128			c$$$ print*,(e(i,j),j=1,3)
129			c$$$ enddo
130			c$$$ print*,'matrix F:' !???
131			c$$$ do i=1,3
132			c$$$ print*,(f(i,j),j=1,3)
133			c$$$ enddo
134
135			c------------------------------------------------------------------------
136			c compute the determinants of A, D, E and F matrixes
137			c using DFACT (http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/f011/top.html)
138			c------------------------------------------------------------------------
139
140			ifail=0
141			jfail=0
142			call DFACT(3,a,3,ir,ifail,deta,jfail)
143			if(ifail.eq.-1) then
144	pam-fi	1.2	if(DEBUG)then
145			print*,'tricircle: ERROR: singular matrix A:'
146			do i=1,3
147			print*,(a(i,j),j=1,3)
148			enddo
149			endif
150	mocchiut	1.1	eflag=1
151			endif
152			if(jfail.eq.-1) then
153	pam-fi	1.2	if(DEBUG)then
154			print*
155			$ ,'tricircle: ERROR: matrix A: determinant too small?'
156			do i=1,3
157			print*,(d(i,j),j=1,3)
158			enddo
159			endif
160	mocchiut	1.1	eflag=1
161			elseif(jfail.eq.1) then
162	pam-fi	1.2	if(DEBUG)then
163			print*
164	mocchiut	1.1	$ ,'tricircle: ERROR: matrix A: determinant too large?'
165	pam-fi	1.2	do i=1,3
166			print*,(d(i,j),j=1,3)
167			enddo
168			endif
169	mocchiut	1.1	eflag=1
170			endif
171
172			ifail=0
173			jfail=0
174			call DFACT(3,d,3,ir,ifail,detd,jfail)
175			if(ifail.eq.-1) then
176	pam-fi	1.2	if(DEBUG)then
177	pam-fi	1.3	print*,'tricircle: ERROR: singular matrix D:'
178			do i=1,3
179			print*,(d(i,j),j=1,3)
180			enddo
181	pam-fi	1.2	endif
182	mocchiut	1.1	eflag=1
183	pam-fi	1.3	endif
184			if(jfail.eq.-1) then
185			if(DEBUG)then
186			print*
187			$ ,'tricircle: ERROR: matrix D: determinant too small?'
188			do i=1,3
189			print*,(d(i,j),j=1,3)
190			enddo
191	pam-fi	1.2	endif
192	mocchiut	1.1	eflag=1
193			elseif(jfail.eq.1) then
194	pam-fi	1.2	if(DEBUG)then
195	mocchiut	1.1	print*
196			$ ,'tricircle: ERROR: matrix D: determinant too large?'
197			do i=1,3
198			print*,(d(i,j),j=1,3)
199			enddo
200	pam-fi	1.2	endif
201	mocchiut	1.1	eflag=1
202			endif
203
204			ifail=0
205			jfail=0
206			call DFACT(3,e,3,ir,ifail,dete,jfail)
207			if(ifail.eq.-1) then
208	pam-fi	1.2	if(DEBUG)then
209	mocchiut	1.1	print*,'tricircle: ERROR: singular matrix E:'
210			do i=1,3
211			print*,(e(i,j),j=1,3)
212			enddo
213	pam-fi	1.2	endif
214	mocchiut	1.1	eflag=1
215	pam-fi	1.3	endif
216	mocchiut	1.1	if(jfail.eq.-1) then
217	pam-fi	1.2	if(DEBUG)then
218	mocchiut	1.1	print*
219			$ ,'tricircle: ERROR: matrix E: determinant too small?'
220			do i=1,3
221			print*,(e(i,j),j=1,3)
222			enddo
223	pam-fi	1.2	endif
224	mocchiut	1.1	eflag=1
225			elseif(jfail.eq.1) then
226	pam-fi	1.2	if(DEBUG)then
227	mocchiut	1.1	print*
228			$ ,'tricircle: ERROR: matrix E: determinant too large?'
229			do i=1,3
230			print*,(e(i,j),j=1,3)
231			enddo
232	pam-fi	1.2	endif
233	mocchiut	1.1	eflag=1
234			endif
235
236			ifail=0
237			jfail=0
238			call DFACT(3,f,3,ir,ifail,detf,jfail)
239			c ifail=-1 !???
240			if(ifail.eq.-1) then
241	pam-fi	1.2	if(DEBUG)then
242	mocchiut	1.1	print*,'tricircle: ERROR: singular matrix F:'
243			do i=1,3
244			print*,(f(i,j),j=1,3)
245			enddo
246	pam-fi	1.2	endif
247	mocchiut	1.1	eflag=1
248			endif
249			if(jfail.eq.-1) then
250	pam-fi	1.2	if(DEBUG)then
251	mocchiut	1.1	print*
252			$ ,'tricircle: ERROR: matrix F: determinant too small?'
253			do i=1,3
254			print*,(f(i,j),j=1,3)
255			enddo
256	pam-fi	1.2	endif
257	mocchiut	1.1	eflag=1
258			elseif(jfail.eq.1) then
259	pam-fi	1.2	if(DEBUG)then
260	mocchiut	1.1	print*
261			$ ,'tricircle: ERROR: matrix F: determinant too large?'
262			do i=1,3
263			print*,(f(i,j),j=1,3)
264			enddo
265	pam-fi	1.2	endif
266	mocchiut	1.1	eflag=1
267			endif
268
269			c------------------------------------------------------------------------
270			c compute the centre and radius
271			c------------------------------------------------------------------------
272			detd=-detd
273			detf=-detf
274
275			c xxc(k)=-detd/(2.*deta)
276			c zzc(k)=-dete/(2.*deta)
277			xxc(k)=-dete/(2.*deta)
278			zzc(k)=-detd/(2.*deta)
279			rrr(k)=SQRT((detd2+dete2)/(4.deta*2.)-detf/deta)
280
281			c$$$ write(30,) xxc(k)scale !???
282			c$$$ write(31,) zzc(k)scale !???
283			c$$$ write(32,) rrr(k)scale !???
284			c$$$ print,'xxc =',xxc(k)scale,' zzc =',zzc(k)*scale
285			c$$$ $ ,' rrr =',rrr(k)*scale !???
286
287			enddo !index loops
288			enddo
289			enddo
290
291
292			c------------------------------------------------------------------------
293			c averages the centres and the radii
294			c------------------------------------------------------------------------
295			xc=0.
296			zc=0.
297			radius=0.
298			do i=1,k
299	mocchiut	1.6	xc=xc+REAL(xxc(i)) !EM GCC4.7
300			zc=zc+REAL(zzc(i)) !EM GCC4.7
301			radius=radius+REAL(rrr(i)) !EM GCC4.7
302	mocchiut	1.1	enddo
303			xc=xc/k * scale !back to micrometers
304			zc=zc/k * scale
305			radius=radius/k * scale
306
307			c$$$ c------------------------------------------------------------------------
308			c$$$ c check for small radius...!???
309			c$$$ c------------------------------------------------------------------------
310			c$$$ num=200
311			c$$$ height=ABS(indep(1)-indep(6))
312			c$$$ if(radius.lt.(num*height)) then
313			c$$$ xc=0.
314			c$$$ zc=0.
315			c$$$ radius=0.
316			c$$$ print*,'tricircle: ERROR: bad circle'
317			c$$$ print*,'radius' ,radius,' < ', num,' x',height
318			c$$$ c$$$ print*,dep !???
319			c$$$ c$$$ print*,indep !???
320			c$$$ eflag=1
321			c$$$ endif
322
323
324			c------------------------------------------------------------------------
325			c computes residuals and chi-squared
326			c------------------------------------------------------------------------
327			chi=0.
328
329			c print*,xc,zc,radius !???
330			do i=1,npoints
331			tmp1 = SQRT(radius2.-(indep(i)-zc)2.)
332			tmp2 = dep(i)-xc
333			if(ABS(tmp2-tmp1).le.ABS(tmp2+tmp1)) then !it chooses the right sign
334			tmp(i)=tmp1 !according to residuals
335			else
336			tmp(i)=-tmp1
337			endif
338	mocchiut	1.6	residual(i)=REAL(tmp2 - tmp(i)) !EM GCC4.7
339	mocchiut	1.1	chi=chi + residual(i)**2.
340			c print*,dep(i) !???
341			c print*,indep(i) !???
342			c print*,tmp1,tmp2,tmp(i),residual(i) !???
343			enddo
344
345			c------------------------------------------------------------------------
346			c it computes the angle between the tangent to the circumference and the
347			c independent variable axis
348			c------------------------------------------------------------------------
349			do i=1,npoints
350	mocchiut	1.6	angle(i)=REAL((zc-indep(i)) / tmp(i)) !EM GCC4.7
351	mocchiut	1.1	angle(i)=ATAN(angle(i)) !-pi/2 <= angle <= pi/2
352			angle(i)=angle(i)/pigr*180.
353			enddo
354
355			return
356			end
357
358
359
360
361
362
363
364
365
366
367			c------------------------------------------------------------------------
368			c------------------------------------------------------------------------
369			c------------------------------------------------------------------------
370
371
372			c$$$c------------------------------------------------------------------------
373			c$$$c Function to find the factorial value
374			c$$$c------------------------------------------------------------------------
375			c$$$
376			c$$$c from http://www.digitalcoding.com/programming/fortran/tutorial/ftute10.htm
377			c$$$
378			c$$$
379			c$$$ FUNCTION FACT(N)
380			c$$$ FACT=1
381			c$$$ DO 10 J=2,N
382			c$$$ FACT=FACT*J
383			c$$$10 CONTINUE
384			c$$$ RETURN
385			c$$$ END
386			c$$$
387			c$$$c------------------------------------------------------------------------