program tbmvmt
      parameter(ntbmax=120)
      real tbnoma(ntbmax,3),tbnomb(ntbmax,3),tbmsa(ntbmax,3)
      real tbmsb(ntbmax,3),tbsrta(ntbmax,3),tbsrtb(ntbmax,3)
      real trans(3,3),delzb(ntbmax),delxb(ntbmax),delyb(ntbmax)
      real delsb(ntbmax)
      character*1 aaa(72)
c
      open(unit=11,file='tbnoms.dat',form='formatted',status='old')
      open(unit=12,file='tbmeas.dat',form='formatted',status='old')
      open(unit=15,file='tbdxyz.dat',form='formatted',status='unknown')
c
c-----define coil angle, 20deg for a "B" coil, 40deg for a "C"
c
      phi=40.
c
c-----read in the nominal tooling ball points for the top coil
c
      do 120 i=1,ntbmax
      read(11,*,end=121)(tbnoma(i,j),j=1,3)
c     read(11,1201,end=121)(tbnoma(i,j),j=1,3)
c     write(15,1004)(tbnoma(i,j),j=1,3)
120   continue
121   continue
c-----number of nominal points
      numnom=i-1
c
c-----read in the tooling ball points measured on the top coil
c
      do 130 i=1,ntbmax
      read(12,*,end=131)(tbmsa(i,j),j=1,3)
130   continue
131   continue
c-----number of measured points
      nmeas=i-1
c
c-----sort the array of measured points so that each entry is aligned with
c-----the corresponding nominal point
c
      do 200 imeas=1,nmeas
c-----compute the distance squared between this measured point and each nominal
      do 150 inom=1,numnom
      dsq=0.
      do 145 j=1,3
      dsq=dsq+(tbmsa(imeas,j)-tbnoma(inom,j))**2
145   continue
c-----if (dsq .lt. 1.) we have a match-->exit the loop
      if(dsq .lt. 1.)go to 155
150   continue
      write(*,1101)(tbmsa(imeas,j),j=1,3)
155   continue
c-----put the "inom" entry of "tbnoma" into the "imeas" entry of "tbsrta"
      do 160 j=1,3
      tbsrta(imeas,j)=tbnoma(inom,j)
160   continue
c
200   continue
c
c-----apply a coordinate trausformation to "tbnoma" and "tbsrta"
c-----rotate through minus phi
c-----this will place the z coordinates normal to the flange
c
c-----convert phi to radians
c
      pi=acos(-1.)
      phi=phi*pi/180.
c
c-----define the transformation matrix
c
      do 210 i=1,3
      do 210 j=1,3
      trans(i,j)=0.
210   continue
c
      trans(1,1)=cos(-phi)
      trans(2,2)=1.
      trans(3,3)=cos(-phi)
      trans(1,3)=sin(-phi)
      trans(3,1)=sin(phi)
c
      do 230 imeas=1,nmeas
c
      do 225 j=1,3
      tbsrtb(imeas,j)=0.
      tbnomb(imeas,j)=0.
c
      do 220 k=1,3
      tbsrtb(imeas,j)=tbsrtb(imeas,j)+trans(j,k)*tbsrta(imeas,k)
      tbmsb(imeas,j)=tbmsb(imeas,j)+trans(j,k)*tbmsa(imeas,k)
220   continue
c
225   continue
c
c-----calculate the x,y,and z differences, now that z is normal to the flange
c
c     delzb(imeas)=tbsrtb(imeas,3)-tbnomb(imeas,3)
      delzb(imeas)=tbmsb(imeas,3)-tbsrtb(imeas,3)
      delxb(imeas)=tbmsb(imeas,1)-tbsrtb(imeas,1)
      delyb(imeas)=tbmsb(imeas,2)-tbsrtb(imeas,2)
      delsb(imeas)=delxb(imeas)**2+delyb(imeas)**2+delzb(imeas)**2
      delsb(imeas)=sqrt(delsb(imeas))
230   continue
c
c-----compute averages and standard deviations
c
      xavg=0.
      yavg=0.
      zavg=0.
      savg=0.
      xstd=0.
      ystd=0.
      zstd=0.
      sstd=0.
      do 235 i=1,nmeas
      xavg=xavg+delxb(i)/nmeas
      yavg=yavg+delyb(i)/nmeas
      zavg=zavg+delzb(i)/nmeas
      savg=savg+delsb(i)/nmeas
235   continue
c
      do 240 i=1,nmeas
      xstd=xstd+((delxb(i)-xavg)**2)/nmeas
      ystd=ystd+((delyb(i)-yavg)**2)/nmeas
      zstd=zstd+((delzb(i)-zavg)**2)/nmeas
      sstd=sstd+((delsb(i)-savg)**2)/nmeas
240   continue
      xstd=sqrt(xstd)
      ystd=sqrt(ystd)
      zstd=sqrt(zstd)
      sstd=sqrt(sstd)
c
c-----write the results
c
      write(15,1003)numnom,nmeas
      write(15,1102)
      write(15,1103)
      write(15,1104)
      do 250 i=1,nmeas
      write(15,1001)(tbmsa(i,j),j=1,3),(tbsrta(i,j),j=1,3)
250   continue
c
      write(15,1102)
c
      write(15,1105)
      write(15,1106)
      do 260 i=1,nmeas
      write(15,1002)(tbsrtb(i,j),j=1,3),(tbmsb(i,j),j=1,3),delxb(i)
     1,delyb(i),delzb(i),delsb(i)
260   continue
c
c-----write the averages
c
      write(15,1005)xavg,yavg,zavg,savg
      write(15,1006)xstd,ystd,zstd,sstd
      stop
c
1001  format(9f8.3)
1002  format(10f8.3)
1003  format(' number of nominals = ',i3,5x,'number measured = ',i3)
1004  format(2f20.2)
1005  format(1x,'xavg = ',f8.4,4x,'yavg = ',f8.4,4x,'zavg = ',f8.4,4x
     1,'savg = ',f8.4)
1006  format(1x,'xstd = ',f8.4,4x,'ystd = ',f8.4,4x,'zstd = ',f8.4,4x
     1,'sstd = ',f8.4)
c
1101  format(' error..failure to match',3f15.5)
1102  format(' ')
1103  format(' locations in half period coordinate system')
1104  format(' xmeas',t12,'ymeas',t20,'zmeas',t28,'xnsrt',t36,'ynsrt'
     1,t44,'znsrt')
1105  format(' locations in coordinate system normal to interface')
1106  format('  xnomb',t11,'ynomb',t19,'znomb',t27,'xmeasb',t35
     1,'ymeasb',t43,'zmeasb',t52,'delx',t60,'dely',t68,'delz'
     2,t76,'dels')
c
1201  format(f7.0,f17.0,f16.0)
1202  format(72a1)
1203  format(3f10.0)
      end