subroutine fitsq(rms,x,y,nn,t,b) c x and y are the vectors of coordinates (dimensioned (3,n)) of the two c structures to be superimposed. nn is 3*n, where n is the number of c points. t and b are respectively the translation vector and the c rotation matrix that transforms the second set of coordinates to the c frame of the first set. c eta = machine-specific variable c dimension x(1),y(1),t(3) dimension x(3*nn),y(3*nn),t(3) dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3) eta = z00100000 c small=25.0*rmdcon(3) c small=25.0*eta c small=25.0*10.e-10 c the following is a very lenient value for 'small' small = 0.0001 fn=nn do 10 i=1,3 xav(i)=0.0 yav(i)=0.0 do 10 j=1,3 10 b(j,i)=0.0 nc=0 c do 30 n=1,nn do 20 i=1,3 c print*,'x = ',x(nc+i),' y = ',y(nc+i) xav(i)=xav(i)+x(nc+i)/fn 20 yav(i)=yav(i)+y(nc+i)/fn 30 nc=nc+3 c * print*,'xav = ',(xav(j),j=1,3) * print*,'yav = ',(yav(j),j=1,3) nc=0 rms=0.0 do 50 n=1,nn do 40 i=1,3 rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn do 40 j=1,3 b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn 40 c(j,i)=b(j,i) 50 nc=nc+3 call sivade(b,q,r,d) sn3=sign(1.0,d) do 120 i=1,3 do 120 j=1,3 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3) call mvvad(b,xav,yav,t) do 130 i=1,3 do 130 j=1,3 rms=rms+2.0*c(j,i)*b(j,i) 130 b(j,i)=-b(j,i) if (abs(rms).gt.small) go to 140 * write (6,301) return 140 if (rms.gt.0.0) go to 150 * write (6,303) rms rms=0.0 * stop * 150 write (6,302) sqrt(rms) 150 rms=sqrt(rms) return 301 format (5x,'rms deviation negligible') 302 format (5x,'rms deviation ',f14.6) 303 format (//,5x,'negative ms deviation - ',f14.6) end subroutine sivade(x,q,r,dt) c computes q,e and r such that q(t)xr = diag(e) dimension x(3,3),q(3,3),r(3,3),e(3) dimension h(3,3),p(3,3),u(3,3),d(3) eta = z00100000 small=25.0*10.e-10 c small=25.0*eta c small=2.0*rmdcon(3) xnrm=0.0 do 20 i=1,3 do 10 j=1,3 xnrm=xnrm+x(j,i)*x(j,i) u(j,i)=0.0 r(j,i)=0.0 10 h(j,i)=0.0 u(i,i)=1.0 20 r(i,i)=1.0 xnrm=sqrt(xnrm) do 110 n=1,2 xmax=0.0 do 30 j=n,3 30 if (abs(x(j,n)).gt.xmax) xmax=abs(x(j,n)) a=0.0 do 40 j=n,3 h(j,n)=x(j,n)/xmax 40 a=a+h(j,n)*h(j,n) a=sqrt(a) den=a*(a+abs(h(n,n))) d(n)=1.0/den h(n,n)=h(n,n)+sign(a,h(n,n)) do 70 i=n,3 s=0.0 do 50 j=n,3 50 s=s+h(j,n)*x(j,i) s=d(n)*s do 60 j=n,3 60 x(j,i)=x(j,i)-s*h(j,n) 70 continue if (n.gt.1) go to 110 xmax=amax1(abs(x(1,2)),abs(x(1,3))) h(2,3)=x(1,2)/xmax h(3,3)=x(1,3)/xmax a=sqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3)) den=a*(a+abs(h(2,3))) d(3)=1.0/den h(2,3)=h(2,3)+sign(a,h(2,3)) do 100 i=1,3 s=0.0 do 80 j=2,3 80 s=s+h(j,3)*x(i,j) s=d(3)*s do 90 j=2,3 90 x(i,j)=x(i,j)-s*h(j,3) 100 continue 110 continue do 130 i=1,3 do 120 j=1,3 120 p(j,i)=-d(1)*h(j,1)*h(i,1) 130 p(i,i)=1.0+p(i,i) do 140 i=2,3 do 140 j=2,3 u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2) 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3) call mmmul(p,u,q) 150 np=1 nq=1 if (abs(x(2,3)).gt.small*(abs(x(2,2))+abs(x(3,3)))) go to 160 x(2,3)=0.0 nq=nq+1 160 if (abs(x(1,2)).gt.small*(abs(x(1,1))+abs(x(2,2)))) go to 180 x(1,2)=0.0 if (x(2,3).ne.0.0) go to 170 nq=nq+1 go to 180 170 np=np+1 180 if (nq.eq.3) go to 310 npq=4-np-nq if (np.gt.npq) go to 230 n0=0 do 220 n=np,npq nn=n+np-1 if (abs(x(nn,nn)).gt.small*xnrm) go to 220 x(nn,nn)=0.0 if (x(nn,nn+1).eq.0.0) go to 220 n0=n0+1 go to (190,210,220),nn 190 do 200 j=2,3 200 call givns(x,q,1,j) go to 220 210 call givns(x,q,2,3) 220 continue if (n0.ne.0) go to 150 230 nn=3-nq a=x(nn,nn)*x(nn,nn) if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn) b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1) c=x(nn,nn)*x(nn,nn+1) dd=0.5*(a-b) xn2=c*c rt=b-xn2/(dd+sign(sqrt(dd*dd+xn2),dd)) y=x(np,np)*x(np,np)-rt z=x(np,np)*x(np,np+1) do 300 n=np,nn if (abs(y).lt.abs(z)) go to 240 t=z/y c=1.0/sqrt(1.0+t*t) s=c*t go to 250 240 t=y/z s=1.0/sqrt(1.0+t*t) c=s*t 250 do 260 j=1,3 v=x(j,n) w=x(j,n+1) x(j,n)=c*v+s*w x(j,n+1)=-s*v+c*w a=r(j,n) b=r(j,n+1) r(j,n)=c*a+s*b 260 r(j,n+1)=-s*a+c*b y=x(n,n) z=x(n+1,n) if (abs(y).lt.abs(z)) go to 270 t=z/y c=1.0/sqrt(1.0+t*t) s=c*t go to 280 270 t=y/z s=1.0/sqrt(1.0+t*t) c=s*t 280 do 290 j=1,3 v=x(n,j) w=x(n+1,j) a=q(j,n) b=q(j,n+1) x(n,j)=c*v+s*w x(n+1,j)=-s*v+c*w q(j,n)=c*a+s*b 290 q(j,n+1)=-s*a+c*b if (n.ge.nn) go to 300 y=x(n,n+1) z=x(n,n+2) 300 continue go to 150 310 do 320 i=1,3 320 e(i)=x(i,i) 330 n0=0 do 360 i=1,3 if (e(i).ge.0.0) go to 350 e(i)=-e(i) do 340 j=1,3 340 q(j,i)=-q(j,i) 350 if (i.eq.1) go to 360 if (abs(e(i)).lt.abs(e(i-1))) go to 360 call switch(i,1,q,r,e) n0=n0+1 360 continue if (n0.ne.0) go to 330 if (abs(e(3)).gt.small*xnrm) go to 370 e(3)=0.0 if (abs(e(2)).gt.small*xnrm) go to 370 e(2)=0.0 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3)) * write (1,501) (e(i),i=1,3) return 501 format (/,5x,'singular values - ',3e15.5) end subroutine givns(a,b,m,n) dimension a(3,3),b(3,3) if (abs(a(m,n)).lt.abs(a(n,n))) go to 10 t=a(n,n)/a(m,n) s=1.0/sqrt(1.0+t*t) c=s*t go to 20 10 t=a(m,n)/a(n,n) c=1.0/sqrt(1.0+t*t) s=c*t 20 do 30 j=1,3 v=a(m,j) w=a(n,j) x=b(j,m) y=b(j,n) a(m,j)=c*v-s*w a(n,j)=s*v+c*w b(j,m)=c*x-s*y 30 b(j,n)=s*x+c*y return end subroutine switch(n,m,u,v,d) dimension u(3,3),v(3,3),d(3) do 10 i=1,3 tem=u(i,n) u(i,n)=u(i,n-1) u(i,n-1)=tem if (m.eq.0) go to 10 tem=v(i,n) v(i,n)=v(i,n-1) v(i,n-1)=tem 10 continue tem=d(n) d(n)=d(n-1) d(n-1)=tem return end c subroutine mvvad(a,b,c,d) subroutine mvvad(b,xav,yav,t) dimension b(3,3),xav(3),yav(3),t(3) c dimension a(3,3),b(3),c(3),d(3) c do 10 j=1,3 c d(j)=c(j) c do 10 i=1,3 c 10 d(j)=d(j)+a(j,i)*b(i) do 10 j=1,3 t(j)=yav(j) do 10 i=1,3 10 t(j)=t(j)+b(j,i)*xav(i) return end real function det (a,b,c) dimension a(3),b(3),c(3) det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3)) 1 +a(3)*(b(1)*c(2)-b(2)*c(1)) return end subroutine mmmul(a,b,c) dimension a(3,3),b(3,3),c(3,3) do 10 i=1,3 do 10 j=1,3 c(i,j)=0.0 do 10 k=1,3 10 c(i,j)=c(i,j)+a(i,k)*b(k,j) return end subroutine matvec(uvec,tmat,pvec,nback) real*4 tmat(3,3),uvec(3,nback), pvec(3,nback) c do 2 j=1,nback do 1 i=1,3 uvec(i,j) = 0.0 do 1 k=1,3 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j) 2 continue return end