1 subroutine fitsq(rms,x,y,nn,t,b,non_conv)
2 implicit real*8 (a-h,o-z)
3 include 'COMMON.IOUNITS'
4 c x and y are the vectors of coordinates (dimensioned (3,n)) of the two
5 c structures to be superimposed. nn is 3*n, where n is the number of
6 c points. t and b are respectively the translation vector and the
7 c rotation matrix that transforms the second set of coordinates to the
8 c frame of the first set.
9 c eta = machine-specific variable
11 dimension x(3*nn),y(3*nn),t(3)
12 dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3)
15 c small=25.0*rmdcon(3)
18 c the following is a very lenient value for 'small'
31 c write(iout,*)'x = ',x(nc+i),' y = ',y(nc+i)
32 xav(i)=xav(i)+x(nc+i)/fn
33 20 yav(i)=yav(i)+y(nc+i)/fn
43 rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2
48 c write(iout,*)'xav = ',(xav(j),j=1,3)
49 c write(iout,*)'yav = ',(yav(j),j=1,3)
50 c write(iout,*)'t = ',(t(j),j=1,3)
51 c write(iout,*)'rms=',rms
52 if (rms.lt.small) return
59 rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn
61 b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn
64 call sivade(b,q,r,d,non_conv)
68 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3)
69 call mvvad(b,xav,yav,t)
72 rms=rms+2.0*c(j,i)*b(j,i)
74 if (dabs(rms).gt.small) go to 140
77 140 if (rms.gt.0.0d0) go to 150
78 c write (iout,303) rms
81 c 150 write (iout,302) dsqrt(rms)
84 301 format (5x,'rms deviation negligible')
85 302 format (5x,'rms deviation ',f14.6)
86 303 format (//,5x,'negative ms deviation - ',f14.6)
89 subroutine sivade(x,q,r,dt,non_conv)
90 implicit real*8(a-h,o-z)
91 c computes q,e and r such that q(t)xr = diag(e)
92 dimension x(3,3),q(3,3),r(3,3),e(3)
93 dimension h(3,3),p(3,3),u(3,3),d(3)
96 c write (2,*) "SIVADE"
100 c small=2.0*rmdcon(3)
104 xnrm=xnrm+x(j,i)*x(j,i)
114 30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n))
120 den=a*(a+dabs(h(n,n)))
122 h(n,n)=h(n,n)+dsign(a,h(n,n))
129 60 x(j,i)=x(j,i)-s*h(j,n)
131 if (n.gt.1) go to 110
132 xmax=dmax1(dabs(x(1,2)),dabs(x(1,3)))
135 a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3))
136 den=a*(a+dabs(h(2,3)))
138 h(2,3)=h(2,3)+sign(a,h(2,3))
145 90 x(i,j)=x(i,j)-s*h(j,3)
150 120 p(j,i)=-d(1)*h(j,1)*h(i,1)
151 130 p(i,i)=1.0+p(i,i)
154 u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2)
155 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3)
160 c write (2,*) "nit",nit," e",(x(i,i),i=1,3)
161 if (nit.gt.10000) then
162 print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
166 if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160
169 160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180
171 if (x(2,3).ne.0.0d0) go to 170
175 180 if (nq.eq.3) go to 310
177 c write (2,*) "np",np," npq",npq
178 if (np.gt.npq) go to 230
182 c write (2,*) "nn",nn
183 if (dabs(x(nn,nn)).gt.small*xnrm) go to 220
185 if (x(nn,nn+1).eq.0.0d0) go to 220
187 c write (2,*) "nn",nn
188 go to (190,210,220),nn
190 200 call givns(x,q,1,j)
192 210 call givns(x,q,2,3)
194 c write (2,*) "nn",nn," np",np," nq",nq," n0",n0
195 c write (2,*) "x",(x(i,i),i=1,3)
196 if (n0.ne.0) go to 150
199 if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn)
200 b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1)
201 c=x(nn,nn)*x(nn,nn+1)
204 rt=b-xn2/(dd+sign(dsqrt(dd*dd+xn2),dd))
205 y=x(np,np)*x(np,np)-rt
206 z=x(np,np)*x(np,np+1)
208 c write (2,*) "n",n," a",a," b",b," c",c," y",y," z",z
209 if (dabs(y).lt.dabs(z)) go to 240
211 c=1.0/dsqrt(1.0d0+t*t)
215 s=1.0/dsqrt(1.0d0+t*t)
225 260 r(j,n+1)=-s*a+c*b
228 if (dabs(y).lt.dabs(z)) go to 270
244 290 q(j,n+1)=-s*a+c*b
245 if (n.ge.nn) go to 300
255 if (nit.gt.10000) then
256 print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
260 c write (2,*) "e",(e(i),i=1,3)
262 if (e(i).ge.0.0d0) go to 350
266 350 if (i.eq.1) go to 360
267 if (dabs(e(i)).lt.dabs(e(i-1))) go to 360
268 call switch(i,1,q,r,e)
271 if (n0.ne.0) go to 330
272 c write (2,*) "e",(e(i),i=1,3)
273 if (dabs(e(3)).gt.small*xnrm) go to 370
275 if (dabs(e(2)).gt.small*xnrm) go to 370
277 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3))
278 c write (2,*) "nit",nit
279 c write (2,501) (e(i),i=1,3)
281 501 format (/,5x,'singular values - ',3e15.5)
283 subroutine givns(a,b,m,n)
284 implicit real*8 (a-h,o-z)
285 dimension a(3,3),b(3,3)
286 if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10
305 subroutine switch(n,m,u,v,d)
306 implicit real*8 (a-h,o-z)
307 dimension u(3,3),v(3,3),d(3)
322 subroutine mvvad(b,xav,yav,t)
323 implicit real*8 (a-h,o-z)
324 dimension b(3,3),xav(3),yav(3),t(3)
325 c dimension a(3,3),b(3),c(3),d(3)
329 c 10 d(j)=d(j)+a(j,i)*b(i)
333 10 t(j)=t(j)+b(j,i)*xav(i)
336 double precision function det (a,b,c)
337 implicit real*8 (a-h,o-z)
338 dimension a(3),b(3),c(3)
339 det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3))
340 1 +a(3)*(b(1)*c(2)-b(2)*c(1))
343 subroutine mmmul(a,b,c)
344 implicit real*8 (a-h,o-z)
345 dimension a(3,3),b(3,3),c(3,3)
350 10 c(i,j)=c(i,j)+a(i,k)*b(k,j)
353 subroutine matvec(uvec,tmat,pvec,nback)
354 implicit real*8 (a-h,o-z)
355 real*8 tmat(3,3),uvec(3,nback), pvec(3,nback)
361 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j)