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 (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

      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------------------------------------------------------------------------
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			real angle(npoints) !angle between the tangent line in the input points and
35			! the independent variable axis
36
37			real residual(npoints) !residuals
38			real chi !sum of squared residuals
39			real xc,zc,radius !circle parameters
40
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			parameter (big=1.e4) !just a number greater than C(npoints,3)
55			double precision xxc(big),zzc(big),rrr(big) !centres and radii to be averaged
56
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	mocchiut	1.1	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			endif
184			if(jfail.eq.-1) then
185	pam-fi	1.2	if(DEBUG)then
186	mocchiut	1.1	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			endif
216			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			xc=xc+xxc(i)
300			zc=zc+zzc(i)
301			radius=radius+rrr(i)
302			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			residual(i)=tmp2 - tmp(i)
339			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			angle(i)=(zc-indep(i)) / tmp(i)
351			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------------------------------------------------------------------------