c23456789012345678901234567890123456789012345678901234567890123456789012
c
c program to test if a rotation is the identity
c
c input data file: file code 11
c   line 1  : lat, lon, angle (best rotation)
c   line 2  : 3 numbers in free format
c   line 3  : rotation matrix
c   line 4  : 3 numbers in free format
c   line 5  : 3 numbers in free format
c   line 6  : 3 numbers in free format
c   line 7  : inverse of covariance matrix       
c   line 8  : 3 numbers in free format
c   line 9  : 3 numbers in free format
c   line 10 : 3 numbers in free format
c
c this program is to be compiled with the following libraries
c   dtrans  : utility routines
c   numlib  : numerical libraries
c
c***********************************************************************
c
      implicit double precision(a-h,o-z)
      real xchi,plev,df,d(2)
      dimension qhat(4),h(3),sigma(3,3)
      character filnam*20
1000  format(a)
c
  1   write (*,*) ' Enter rotation file name'
      read(*,1000) filnam
      open (unit=11,file=filnam,status='old',err=1)
c
c read 2 first comment lines
      read(11,1000) filnam
      read(11,1000) filnam
c read bestfit rotation: alat,alon and angle
      read(11,1000) filnam
      read(11,*) alat,along,rho
c read confidence level
      read(11,1000) filnam
      read(11,*) dum
c read hatkap and degree of freedom
      read(11,1000) filnam
      read(11,*) hatkap,df
      if (df.eq.0.) df=10000.
c read nbr of points, sections
      read(11,1000) filnam
      read(11,*) dum
c read rotation matrix
      read(11,1000) filnam
      do 5 i=1,3
    5 read(11,*) dum
c read covariance matrix
      read(11,1000) filnam
      do 10 i=1,3
   10 read(11,*) dum
c read inverse of covariance matrix: H11.2
      read(11,1000) filnam
      do 15 i=1,3
   15 read(11,*) (sigma(i,j),j=1,3)
      close(11)
c
      if(hatkap.gt.1.) then
	 write(*,*)  ' kappahat = ',hatkap
         write (*,*) ' Do you want kappahat to be set to 1?'
         write (*,*) '          1=no    2=yes'
         read (*,*) ichoice
      endif
      call trans3(alat,along,rho,qhat)
      call trans8(qhat,h)
      if (df.ge.1000.or.ichoice.eq.2) hatkap=1
      dxchi=0.
      do 20 i=1,3
      do 20 j=1,3
   20 dxchi=dxchi+h(i)*sigma(i,j)*h(j)*hatkap
      xchi=dxchi
      if (df.ge.1000.or.ichoice.eq.2) then
         call xdch(plev,3.,xchi,ier)
         if (ier.ne.0) write(6,*) 'subroutine xdch--ier ',ier
         write(*,*) 'chi-sq statistic, df:',xchi,3.
      else
         xchi=xchi/3.
         d(1)=3.
         d(2)=df
         call xdf(plev,d,xchi,ier)
         if (ier.ne.0) write(6,*) 'subroutine xdf--ier ',ier
         write(*,*) 'f statistic, df:',xchi,d
      endif
      plev=1-plev
      write(*,*) 'Significance level: ',plev
      stop
      end