1 subroutine fitsq(rms,x,y,nn,t,b,non_conv)
2 implicit real*8 (a-h,o-z)
3 include 'DIMENSIONS.ZSCOPT'
4 include 'COMMON.IOUNITS'
5 c x and y are the vectors of coordinates (dimensioned (3,n)) of the two
6 c structures to be superimposed. nn is 3*n, where n is the number of
7 c points. t and b are respectively the translation vector and the
8 c rotation matrix that transforms the second set of coordinates to the
9 c frame of the first set.
10 c eta = machine-specific variable
12 dimension x(3*nn),y(3*nn),t(3)
13 dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3)
16 c small=25.0*rmdcon(3)
19 c the following is a very lenient value for 'small'
32 crc write(iout,*)'x = ',x(nc+i),' y = ',y(nc+i)
33 xav(i)=xav(i)+x(nc+i)/fn
34 20 yav(i)=yav(i)+y(nc+i)/fn
44 rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2
49 c write(iout,*)'xav = ',(xav(j),j=1,3)
50 c write(iout,*)'yav = ',(yav(j),j=1,3)
51 c write(iout,*)'t = ',(t(j),j=1,3)
52 c write(iout,*)'rms=',rms
53 if (rms.lt.small) return
60 rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn
62 b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn
65 call sivade(b,q,r,d,non_conv)
69 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3)
70 call mvvad(b,xav,yav,t)
73 rms=rms+2.0*c(j,i)*b(j,i)
75 if (dabs(rms).gt.small) go to 140
78 140 if (rms.gt.0.0d0) go to 150
79 c write (iout,303) rms
82 c 150 write (iout,302) dsqrt(rms)
85 301 format (5x,'rms deviation negligible')
86 302 format (5x,'rms deviation ',f14.6)
87 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)
103 xnrm=xnrm+x(j,i)*x(j,i)
113 30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n))
119 den=a*(a+dabs(h(n,n)))
121 h(n,n)=h(n,n)+dsign(a,h(n,n))
128 60 x(j,i)=x(j,i)-s*h(j,n)
130 if (n.gt.1) go to 110
131 xmax=dmax1(dabs(x(1,2)),dabs(x(1,3)))
134 a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3))
135 den=a*(a+dabs(h(2,3)))
137 h(2,3)=h(2,3)+sign(a,h(2,3))
144 90 x(i,j)=x(i,j)-s*h(j,3)
149 120 p(j,i)=-d(1)*h(j,1)*h(i,1)
150 130 p(i,i)=1.0+p(i,i)
153 u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2)
154 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3)
159 if (nit.gt.10000) then
160 print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
164 if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160
167 160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180
169 if (x(2,3).ne.0.0d0) go to 170
173 180 if (nq.eq.3) go to 310
175 if (np.gt.npq) go to 230
179 if (dabs(x(nn,nn)).gt.small*xnrm) go to 220
181 if (x(nn,nn+1).eq.0.0d0) go to 220
183 go to (190,210,220),nn
185 200 call givns(x,q,1,j)
187 210 call givns(x,q,2,3)
189 if (n0.ne.0) go to 150
192 if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn)
193 b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1)
194 c=x(nn,nn)*x(nn,nn+1)
197 rt=b-xn2/(dd+sign(dsqrt(dd*dd+xn2),dd))
198 y=x(np,np)*x(np,np)-rt
199 z=x(np,np)*x(np,np+1)
201 if (dabs(y).lt.dabs(z)) go to 240
203 c=1.0/dsqrt(1.0d0+t*t)
207 s=1.0/dsqrt(1.0d0+t*t)
217 260 r(j,n+1)=-s*a+c*b
220 if (dabs(y).lt.dabs(z)) go to 270
236 290 q(j,n+1)=-s*a+c*b
237 if (n.ge.nn) go to 300
247 if (nit.gt.10000) then
248 print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
253 if (e(i).ge.0.0d0) go to 350
257 350 if (i.eq.1) go to 360
258 if (dabs(e(i)).lt.dabs(e(i-1))) go to 360
259 call switch(i,1,q,r,e)
262 if (n0.ne.0) go to 330
263 if (dabs(e(3)).gt.small*xnrm) go to 370
265 if (dabs(e(2)).gt.small*xnrm) go to 370
267 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3))
268 * write (1,501) (e(i),i=1,3)
270 501 format (/,5x,'singular values - ',3e15.5)
272 subroutine givns(a,b,m,n)
273 implicit real*8 (a-h,o-z)
274 dimension a(3,3),b(3,3)
275 if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10
294 subroutine switch(n,m,u,v,d)
295 implicit real*8 (a-h,o-z)
296 dimension u(3,3),v(3,3),d(3)
311 subroutine mvvad(b,xav,yav,t)
312 implicit real*8 (a-h,o-z)
313 dimension b(3,3),xav(3),yav(3),t(3)
314 c dimension a(3,3),b(3),c(3),d(3)
318 c 10 d(j)=d(j)+a(j,i)*b(i)
322 10 t(j)=t(j)+b(j,i)*xav(i)
325 double precision function det (a,b,c)
326 implicit real*8 (a-h,o-z)
327 dimension a(3),b(3),c(3)
328 det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3))
329 1 +a(3)*(b(1)*c(2)-b(2)*c(1))
332 subroutine mmmul(a,b,c)
333 implicit real*8 (a-h,o-z)
334 dimension a(3,3),b(3,3),c(3,3)
339 10 c(i,j)=c(i,j)+a(i,k)*b(k,j)
342 subroutine matvec(uvec,tmat,pvec,nback)
343 implicit real*8 (a-h,o-z)
344 real*8 tmat(3,3),uvec(3,nback), pvec(3,nback)
350 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j)