/[PAMELA software]/tracker/ground/source/analysis/tricircle.f
ViewVC logotype

Annotation of /tracker/ground/source/analysis/tricircle.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Wed Mar 8 15:00:39 2006 UTC (18 years, 9 months ago) by pam-fi
Branch: MAIN, trk-ground
CVS Tags: R3v02, HEAD
Changes since 1.1: +0 -0 lines
First CVS release of tracker ground software (R3v02) 

1 pam-fi 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     real pigr !3.1415...
60     pigr=ACOS(-1.)
61    
62     eflag = 0
63    
64    
65     c------------------------------------------------------------------------
66     c choose 3 points out of npoints
67     c------------------------------------------------------------------------
68     c nloops = fact(npoints) / (fact(3) * fact(npoints-3))
69    
70     k=0
71     do i1=1,npoints-2
72     index(1)=i1
73     do i2=i1+1,npoints-1
74     index(2)=i2
75     do i3=i2+1,npoints
76     index(3)=i3
77    
78     k=k+1 !number of combinations
79     c$$$ print*,' ' !???
80     c$$$ print*,'k =',k,' index =',index
81    
82     c------------------------------------------------------------------------
83     c build temp vectors
84     c------------------------------------------------------------------------
85     do i=1,3
86     z(i)=indep(index(i))/scale !to avoid too big numbers in matrix computation
87     x(i)=dep(index(i))/scale
88     unit(i)=1.
89     zzxx(i)=z(i)**2.+x(i)**2.
90     enddo
91     c$$$ print*,'z =',z,' x =',x !???
92     c$$$ print*,'unit =',unit,' zzxx =',zzxx
93    
94     c------------------------------------------------------------------------
95     c build the matrixes
96     c------------------------------------------------------------------------
97     do i=1,3
98     a(i,1)=z(i) !A has (z x 1) as columns
99     a(i,2)=x(i)
100     a(i,3)=unit(i)
101     d(i,1)=zzxx(i) !D has (zzxx x 1) as columns
102     d(i,2)=x(i)
103     d(i,3)=unit(i)
104     e(i,1)=zzxx(i) !E has (zzxx z 1) as columns
105     e(i,2)=z(i)
106     e(i,3)=unit(i)
107     f(i,1)=zzxx(i) !F has (zzxx z x) as columns
108     f(i,2)=z(i)
109     f(i,3)=x(i)
110     enddo
111    
112     c$$$ print*,'matrix A:' !???
113     c$$$ do i=1,3
114     c$$$ print*,(a(i,j),j=1,3)
115     c$$$ enddo
116     c$$$ print*,'matrix D:' !???
117     c$$$ do i=1,3
118     c$$$ print*,(d(i,j),j=1,3)
119     c$$$ enddo
120     c$$$ print*,'matrix E:' !???
121     c$$$ do i=1,3
122     c$$$ print*,(e(i,j),j=1,3)
123     c$$$ enddo
124     c$$$ print*,'matrix F:' !???
125     c$$$ do i=1,3
126     c$$$ print*,(f(i,j),j=1,3)
127     c$$$ enddo
128    
129     c------------------------------------------------------------------------
130     c compute the determinants of A, D, E and F matrixes
131     c using DFACT (http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/f011/top.html)
132     c------------------------------------------------------------------------
133    
134     ifail=0
135     jfail=0
136     call DFACT(3,a,3,ir,ifail,deta,jfail)
137     if(ifail.eq.-1) then
138     print*,'tricircle: ERROR: singular matrix A:'
139     do i=1,3
140     print*,(a(i,j),j=1,3)
141     enddo
142     eflag=1
143     endif
144     if(jfail.eq.-1) then
145     print*
146     $ ,'tricircle: ERROR: matrix A: determinant too small?'
147     do i=1,3
148     print*,(d(i,j),j=1,3)
149     enddo
150     eflag=1
151     elseif(jfail.eq.1) then
152     print*
153     $ ,'tricircle: ERROR: matrix A: determinant too large?'
154     do i=1,3
155     print*,(d(i,j),j=1,3)
156     enddo
157     eflag=1
158     endif
159    
160     ifail=0
161     jfail=0
162     call DFACT(3,d,3,ir,ifail,detd,jfail)
163     if(ifail.eq.-1) then
164     print*,'tricircle: ERROR: singular matrix D:'
165     do i=1,3
166     print*,(d(i,j),j=1,3)
167     enddo
168     eflag=1
169     endif
170     if(jfail.eq.-1) then
171     print*
172     $ ,'tricircle: ERROR: matrix D: determinant too small?'
173     do i=1,3
174     print*,(d(i,j),j=1,3)
175     enddo
176     eflag=1
177     elseif(jfail.eq.1) then
178     print*
179     $ ,'tricircle: ERROR: matrix D: determinant too large?'
180     do i=1,3
181     print*,(d(i,j),j=1,3)
182     enddo
183     eflag=1
184     endif
185    
186     ifail=0
187     jfail=0
188     call DFACT(3,e,3,ir,ifail,dete,jfail)
189     if(ifail.eq.-1) then
190     print*,'tricircle: ERROR: singular matrix E:'
191     do i=1,3
192     print*,(e(i,j),j=1,3)
193     enddo
194     eflag=1
195     endif
196     if(jfail.eq.-1) then
197     print*
198     $ ,'tricircle: ERROR: matrix E: determinant too small?'
199     do i=1,3
200     print*,(e(i,j),j=1,3)
201     enddo
202     eflag=1
203     elseif(jfail.eq.1) then
204     print*
205     $ ,'tricircle: ERROR: matrix E: determinant too large?'
206     do i=1,3
207     print*,(e(i,j),j=1,3)
208     enddo
209     eflag=1
210     endif
211    
212     ifail=0
213     jfail=0
214     call DFACT(3,f,3,ir,ifail,detf,jfail)
215     c ifail=-1 !???
216     if(ifail.eq.-1) then
217     print*,'tricircle: ERROR: singular matrix F:'
218     do i=1,3
219     print*,(f(i,j),j=1,3)
220     enddo
221     eflag=1
222     endif
223     if(jfail.eq.-1) then
224     print*
225     $ ,'tricircle: ERROR: matrix F: determinant too small?'
226     do i=1,3
227     print*,(f(i,j),j=1,3)
228     enddo
229     eflag=1
230     elseif(jfail.eq.1) then
231     print*
232     $ ,'tricircle: ERROR: matrix F: determinant too large?'
233     do i=1,3
234     print*,(f(i,j),j=1,3)
235     enddo
236     eflag=1
237     endif
238    
239     c------------------------------------------------------------------------
240     c compute the centre and radius
241     c------------------------------------------------------------------------
242     detd=-detd
243     detf=-detf
244    
245     c xxc(k)=-detd/(2.*deta)
246     c zzc(k)=-dete/(2.*deta)
247     xxc(k)=-dete/(2.*deta)
248     zzc(k)=-detd/(2.*deta)
249     rrr(k)=SQRT((detd**2+dete**2)/(4.*deta**2.)-detf/deta)
250    
251     c$$$ write(30,*) xxc(k)*scale !???
252     c$$$ write(31,*) zzc(k)*scale !???
253     c$$$ write(32,*) rrr(k)*scale !???
254     c$$$ print*,'xxc =',xxc(k)*scale,' zzc =',zzc(k)*scale
255     c$$$ $ ,' rrr =',rrr(k)*scale !???
256    
257     enddo !index loops
258     enddo
259     enddo
260    
261    
262     c------------------------------------------------------------------------
263     c averages the centres and the radii
264     c------------------------------------------------------------------------
265     xc=0.
266     zc=0.
267     radius=0.
268     do i=1,k
269     xc=xc+xxc(i)
270     zc=zc+zzc(i)
271     radius=radius+rrr(i)
272     enddo
273     xc=xc/k * scale !back to micrometers
274     zc=zc/k * scale
275     radius=radius/k * scale
276    
277     c$$$ c------------------------------------------------------------------------
278     c$$$ c check for small radius...!???
279     c$$$ c------------------------------------------------------------------------
280     c$$$ num=200
281     c$$$ height=ABS(indep(1)-indep(6))
282     c$$$ if(radius.lt.(num*height)) then
283     c$$$ xc=0.
284     c$$$ zc=0.
285     c$$$ radius=0.
286     c$$$ print*,'tricircle: ERROR: bad circle'
287     c$$$ print*,'radius' ,radius,' < ', num,' x',height
288     c$$$ c$$$ print*,dep !???
289     c$$$ c$$$ print*,indep !???
290     c$$$ eflag=1
291     c$$$ endif
292    
293    
294     c------------------------------------------------------------------------
295     c computes residuals and chi-squared
296     c------------------------------------------------------------------------
297     chi=0.
298    
299     c print*,xc,zc,radius !???
300     do i=1,npoints
301     tmp1 = SQRT(radius**2.-(indep(i)-zc)**2.)
302     tmp2 = dep(i)-xc
303     if(ABS(tmp2-tmp1).le.ABS(tmp2+tmp1)) then !it chooses the right sign
304     tmp(i)=tmp1 !according to residuals
305     else
306     tmp(i)=-tmp1
307     endif
308     residual(i)=tmp2 - tmp(i)
309     chi=chi + residual(i)**2.
310     c print*,dep(i) !???
311     c print*,indep(i) !???
312     c print*,tmp1,tmp2,tmp(i),residual(i) !???
313     enddo
314    
315     c------------------------------------------------------------------------
316     c it computes the angle between the tangent to the circumference and the
317     c independent variable axis
318     c------------------------------------------------------------------------
319     do i=1,npoints
320     angle(i)=(zc-indep(i)) / tmp(i)
321     angle(i)=ATAN(angle(i)) !-pi/2 <= angle <= pi/2
322     angle(i)=angle(i)/pigr*180.
323     enddo
324    
325     return
326     end
327    
328    
329    
330    
331    
332    
333    
334    
335    
336    
337     c------------------------------------------------------------------------
338     c------------------------------------------------------------------------
339     c------------------------------------------------------------------------
340    
341    
342     c$$$c------------------------------------------------------------------------
343     c$$$c Function to find the factorial value
344     c$$$c------------------------------------------------------------------------
345     c$$$
346     c$$$c from http://www.digitalcoding.com/programming/fortran/tutorial/ftute10.htm
347     c$$$
348     c$$$
349     c$$$ FUNCTION FACT(N)
350     c$$$ FACT=1
351     c$$$ DO 10 J=2,N
352     c$$$ FACT=FACT*J
353     c$$$10 CONTINUE
354     c$$$ RETURN
355     c$$$ END
356     c$$$
357     c$$$c------------------------------------------------------------------------

  ViewVC Help
Powered by ViewVC 1.1.23