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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show 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 *************************************************************************
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