c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c program invrot ====== November 3 1988 J-Y R =======
c
c modified june,1993 to work with addplus.f
c
c calculate inverse rotation and its covariance matrix
c a=b**t
c cova=b*covb*(b**t)
c
c for more details on notations and meanings of the purpose
c of this program, see:
c Chang T., Estimating the relative rotation of two tectonic
c plates from boundary crossings.
c J. Amer. Stat. Assn., scheduled to appear Dec. 1988
c
program invrot
implicit double precision (a-h,o-z)
character filnam*70,name*15
dimension a(3,3),cova(3,3),b(3,3),covb(3,3),
&w(3,3),xa(3),xb(3)
c
1 write (*,1000)
1000 format ('Input file: ',$)
read (*,'(a)') name
open (unit=10,file=name,status='old',err=1)
write (*,'(a)') name
2 write(*,1002)
1002 format('Output file: ',$)
read(*,'(a)') filnam
open(unit=15,file=filnam,status='new',err=2)
write(*,'(a)') filnam
c
c read two first comment lines
read (10,'(a)') filnam
write (*,'(a)') filnam
read (10,'(a)') filnam
write (*,'(a)') filnam
c read alat,alon,angle
read(10,'(a)') filnam
write (*,'(a)') filnam
read(10,*) xb
write (*,'(3x,3f8.2)') xb
c read conf. level
read(10,'(a)') filnam
write (*,'(a)') filnam
read (10,*) plev,fkap2,fkap1
write(*,'(2x,3f6.2)') plev,fkap2,fkap1
c read kappa estimated, degrees of freedom
read(10,'(a)') filnam
write (*,'(a)') filnam
read(10,*) hatkap,df,xchi,flag
write(*,'(3x,f5.2,f5.0,f6.2)') hatkap,df,xchi,flag
c read number of points, sections ...
read(10,'(a)') dum
if (flag .ne. 1.) then
read(10,*) ndat,nsect
else
read(10,*) npts,fmin,fmax,rank
endif
c read rotation matrix
read(10,'(a)') filnam
write (*,'(a)') filnam
do 3 i=1,3
read(10,*) (b(i,j),j=1,3)
3 write (*,*) (b(i,j),j=1,3)
c read covariance matrix
read(10,'(a)') filnam
write (*,'(a)') filnam
do 5 i=1,3
read (10,*) (covb(i,j),j=1,3)
5 write (*,*) (covb(i,j),j=1,3)
write(*,*) ' '
close(10)
c
c calculate covariance matrix for a=b**t
35 call transp(b,a)
call mul(covb,a,w)
call mul(b,w,cova)
xa(1)=xb(1)
xa(2)=xb(2)
xa(3)=-xb(3)
call inv(cova,w)
c
50 write (*,*) ' Inverse rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') xa
write(*,*) 'kappa and degrees of freedom'
write(*,'(3x,f5.2,f5.0)') hatkap,df
write (*,*) 'Rotation matrix'
do 55 i=1,3
55 write (*,*) (a(i,j),j=1,3)
write (*,*) 'Covariance matrix'
do 60 i=1,3
60 write (*,*) (cova(i,j),j=1,3)
c
write(15,1060) name
1060 format ('Inverse rotation from ',a,/)
write(15,*) 'Fitted rotation--alat,along,rho: '
write(15,*) xa
write(15,*) 'conf. level, conf. interval for kappa: '
write(15,*) plev,fkap2,fkap1
write(15,*) 'kappahat, degrees of freedom,xchi,flag'
1061 format(1x,f8.5,2x,f4.0,3x,f8.4,3x,f2.0)
write(15,1061) hatkap,df,xchi,flag
if (flag .ne. 1.) then
write(15,*) 'Number of points, sections'
write(15,*) ndat,nsect
else
write(15,*)'number of points, fmin,fmax,rank'
write(15,'(7x,i4,3f6.0)')npts,fmin,fmax,rank
endif
write(15,*) 'ahat: '
do 70 i=1,3
70 write(15,*) (a(i,j),j=1,3)
write(15,*) 'covariance matrix'
do 75 i=1,3
75 write(15,*) (cova(i,j),j=1,3)
write(15,*) 'H11.2 matrix: '
do 80 i=1,3
80 write(15,*) (w(i,j),j=1,3)
close (15)
c
stop
end
c****************************************************************
c
subroutine inv(a,b)
c invert 3x3 matrix
double precision a(3,3), b(3,3),det
det=a(1,1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3))-
& a(1,2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3))+
& a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
if (det.eq.0.) goto 5
b(1,1)=a(2,2)*a(3,3)-a(3,2)*a(2,3)
b(1,2)=a(3,2)*a(1,3)-a(1,2)*a(3,3)
b(1,3)=a(1,2)*a(2,3)-a(2,2)*a(1,3)
b(2,1)=a(3,1)*a(2,3)-a(2,1)*a(3,3)
b(2,2)=a(1,1)*a(3,3)-a(3,1)*a(1,3)
b(2,3)=a(2,1)*a(1,3)-a(1,1)*a(2,3)
b(3,1)=a(2,1)*a(3,2)-a(3,1)*a(2,2)
b(3,2)=a(3,1)*a(1,2)-a(1,1)*a(3,2)
b(3,3)=a(1,1)*a(2,2)-a(2,1)*a(1,2)
do 1 i=1,3
do 1 j=1,3
1 b(i,j)=b(i,j)/det
return
5 write (*,*) 'Trouble: det=0'
stop
end
c****************************************************************
c
subroutine mul(a,b,c)
c multiply 3x3 matrices: c=b.a
double precision a(3,3),b(3,3),c(3,3)
do 5 i=1,3
do 5 j=1,3
c(i,j)=0.
do 5 k=1,3
5 c(i,j)=c(i,j)+a(i,k)*b(k,j)
return
end
c****************************************************************
c
subroutine transp(a,b)
c calculate transpose of a: b=a**t
double precision a(3,3),b(3,3)
do 5 i=1,3
do 5 j=1,3
5 b(i,j)=a(j,i)
return
end
c****************************************************************