1 subroutine fitsq(rms,x,y,nn,t,b)
2 c x and y are the vectors of coordinates (dimensioned (3,n)) of the two
3 c structures to be superimposed. nn is 3*n, where n is the number of
4 c points. t and b are respectively the translation vector and the
5 c rotation matrix that transforms the second set of coordinates to the
6 c frame of the first set.
7 c eta = machine-specific variable
9 c dimension x(1),y(1),t(3)
10 dimension x(3*nn),y(3*nn),t(3)
11 dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3)
13 c small=25.0*rmdcon(3)
16 c the following is a very lenient value for 'small'
28 c print*,'x = ',x(nc+i),' y = ',y(nc+i)
29 xav(i)=xav(i)+x(nc+i)/fn
30 20 yav(i)=yav(i)+y(nc+i)/fn
34 * print*,'xav = ',(xav(j),j=1,3)
35 * print*,'yav = ',(yav(j),j=1,3)
41 rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn
43 b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn
50 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3)
51 call mvvad(b,xav,yav,t)
54 rms=rms+2.0*c(j,i)*b(j,i)
56 if (abs(rms).gt.small) go to 140
59 140 if (rms.gt.0.0) go to 150
63 * 150 write (6,302) sqrt(rms)
66 301 format (5x,'rms deviation negligible')
67 302 format (5x,'rms deviation ',f14.6)
68 303 format (//,5x,'negative ms deviation - ',f14.6)
70 subroutine sivade(x,q,r,dt)
71 c computes q,e and r such that q(t)xr = diag(e)
72 dimension x(3,3),q(3,3),r(3,3),e(3)
73 dimension h(3,3),p(3,3),u(3,3),d(3)
81 xnrm=xnrm+x(j,i)*x(j,i)
91 30 if (abs(x(j,n)).gt.xmax) xmax=abs(x(j,n))
99 h(n,n)=h(n,n)+sign(a,h(n,n))
106 60 x(j,i)=x(j,i)-s*h(j,n)
108 if (n.gt.1) go to 110
109 xmax=amax1(abs(x(1,2)),abs(x(1,3)))
112 a=sqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3))
113 den=a*(a+abs(h(2,3)))
115 h(2,3)=h(2,3)+sign(a,h(2,3))
122 90 x(i,j)=x(i,j)-s*h(j,3)
127 120 p(j,i)=-d(1)*h(j,1)*h(i,1)
128 130 p(i,i)=1.0+p(i,i)
131 u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2)
132 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3)
136 if (abs(x(2,3)).gt.small*(abs(x(2,2))+abs(x(3,3)))) go to 160
139 160 if (abs(x(1,2)).gt.small*(abs(x(1,1))+abs(x(2,2)))) go to 180
141 if (x(2,3).ne.0.0) go to 170
145 180 if (nq.eq.3) go to 310
147 if (np.gt.npq) go to 230
151 if (abs(x(nn,nn)).gt.small*xnrm) go to 220
153 if (x(nn,nn+1).eq.0.0) go to 220
155 go to (190,210,220),nn
157 200 call givns(x,q,1,j)
159 210 call givns(x,q,2,3)
161 if (n0.ne.0) go to 150
164 if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn)
165 b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1)
166 c=x(nn,nn)*x(nn,nn+1)
169 rt=b-xn2/(dd+sign(sqrt(dd*dd+xn2),dd))
170 y=x(np,np)*x(np,np)-rt
171 z=x(np,np)*x(np,np+1)
173 if (abs(y).lt.abs(z)) go to 240
189 260 r(j,n+1)=-s*a+c*b
192 if (abs(y).lt.abs(z)) go to 270
208 290 q(j,n+1)=-s*a+c*b
209 if (n.ge.nn) go to 300
218 if (e(i).ge.0.0) go to 350
222 350 if (i.eq.1) go to 360
223 if (abs(e(i)).lt.abs(e(i-1))) go to 360
224 call switch(i,1,q,r,e)
227 if (n0.ne.0) go to 330
228 if (abs(e(3)).gt.small*xnrm) go to 370
230 if (abs(e(2)).gt.small*xnrm) go to 370
232 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3))
233 * write (1,501) (e(i),i=1,3)
235 501 format (/,5x,'singular values - ',3e15.5)
237 subroutine givns(a,b,m,n)
238 dimension a(3,3),b(3,3)
239 if (abs(a(m,n)).lt.abs(a(n,n))) go to 10
258 subroutine switch(n,m,u,v,d)
259 dimension u(3,3),v(3,3),d(3)
274 c subroutine mvvad(a,b,c,d)
275 subroutine mvvad(b,xav,yav,t)
276 dimension b(3,3),xav(3),yav(3),t(3)
277 c dimension a(3,3),b(3),c(3),d(3)
281 c 10 d(j)=d(j)+a(j,i)*b(i)
285 10 t(j)=t(j)+b(j,i)*xav(i)
288 real function det (a,b,c)
289 dimension a(3),b(3),c(3)
290 det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3))
291 1 +a(3)*(b(1)*c(2)-b(2)*c(1))
294 subroutine mmmul(a,b,c)
295 dimension a(3,3),b(3,3),c(3,3)
300 10 c(i,j)=c(i,j)+a(i,k)*b(k,j)
303 subroutine matvec(uvec,tmat,pvec,nback)
304 real*4 tmat(3,3),uvec(3,nback), pvec(3,nback)
310 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j)