c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c program addrot ====== November 3 1988 J-Y R =======
c
C MODIFIED OCTOBER 14, 1991 -- MODIFICATIONS IN CAPS
C COMPILE WITH NUMLIB.F
C
c add two rotations and calculate the covariance matrix of the
c product rotation
c chat=bhat.ahat
c =b*phi(tbhat)*a*phi(tahat)
c =b*a*(a**t)*phi(tbhat)*a*phi(tahat)
c =b*a*phi((a**t)*tbhat)*phi(tahat)
c ~b*a*phi((a**t)*tbhat +tahat) up to linear approximations
c
c tchat~(a**t)*tbhat +tahat
c cov(tchat)~(a**t)*cov(tbhat)*a + cov(tahat)
C
C IF THE DEGREES OF FREEDOM IN KAPPAHAT FOR A AND FOR B ARE BOTH
C GREATER THAN 50, THEY ARE EACH ASSUMED TO THE THE TRUE KAPPAS.
C OTHERWISE A POOLED KAPPAHAT IS USED. IN THIS CASE, AN ERROR
C MESSAGE IS PRODUCED IF AN F-TEST WOULD CONTRADICT THE EQUALITY
C OF THE KAPPAS.
C
C IF KAPPA IS KNOWN FOR ONE ROTATION, BUT ESTIMATED FOR THE OTHER,
C AND IF THE ESTIMATED KAPPA HAS FEWER THAN 50 DEGREES OF FREEDOM,
C THE ESTIMATED KAPPA IS ASSUMED TO BE CORRECT, AND AN ERROR MESSAGE
C IS PRINTED.
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 addrot
implicit double precision (a-h,o-z)
character filnam*70,name(2)*15,dum*1
dimension a(3,3),cova(3,3),b(3,3),covb(3,3),c(3,3),covc(3,3),
&w(3,3),xa(3),xb(3),xc(3)
REAL PLEV,PLEVF,ALPHA,ALPHA1,DFC,D(2),XF,XCHI
DATA NDAT/0/,NSECT/0/
c
itour=1
1 write (*,1000) itour
1000 format ('Input file ',i1,': ',$)
read (*,'(a)') name(itour)
open (unit=10,file=name(itour),status='old',err=1)
write (*,'(a)') name(itour)
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,d1,d2
write(*,'(3f6.2)') plev,d1,d2
c read kappa estimated, degrees of freedom
read(10,'(a)') filnam
write (*,'(a)') filnam
read(10,*) hatkapb,dfb
write(*,'(3x,f5.2,f5.0)') hatkapb,dfb
c read number of points, sections ...
read(10,'(a)') dum
read(10,*) d1
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
itour=itour+1
goto (20,35) itour-1
20 do 30 i=1,3
xa(i)=xb(i)
do 30 j=1,3
a(i,j)=b(i,j)
30 cova(i,j)=covb(i,j)
hatkapa=hatkapb
dfa=dfb
goto 1
c
C ADDITIONS OF OCTOBER 14, 1991
C CALCULATION OF COMBINED HATKAPC
C
35 IFLAG1=0
IFLAG2=0
IF ((DFA.GT.50.).AND.(DFB.GT.50.)) THEN
DO 36 I=1,3
DO 36 J=1,3
COVA(I,J)=COVA(I,J)/HATKAPA
36 COVB(I,J)=COVB(I,J)/HATKAPB
DFC = 10000.
HATKAPC = 1.
PLEV = -1.
FKAP1=0.
FKAP2=0.
ELSE IF ((DFA.GT.9999.).OR.(DFB.GT.9999)) THEN
WRITE(*,9999)
9999 FORMAT('***A KNOWN KAPPA HAS BEEN COMBINED WITH AN ESTIMATED'/
&'KAPPA WITH FEWER THAN 50 DEGREES OF FREEDOM***')
IFLAG1=1
DO 37 I=1,3
DO 37 J=1,3
COVA(I,J)=COVA(I,J)/HATKAPA
37 COVB(I,J)=COVB(I,J)/HATKAPB
DFC = 10000.
HATKAPC = 1.
PLEV = -1.
FKAP1=0.
FKAP2=0.
ELSE
PLEV = .95
ALPHA=(1.-PLEV)/2.
ALPHA1=(1.+PLEV)/2.
D(1)=DFA
D(2)=DFB
XF=HATKAPB/HATKAPA
CALL XDF(PLEVF,D,XF,IER)
IF (IER.NE.0) WRITE(*,*) 'ERROR IN XDF'
IF ((PLEVF.LT.ALPHA).OR.(PLEVF.GT.ALPHA1)) THEN
IFLAG2=1
WRITE(*,9998)
9998 FORMAT('***POOLED KAPPAC USED WHEN TWO KAPPAS ARE SIGNICANTLY',
&' DIFFERENT AT ALPHA=.05***')
ENDIF
DFC=DFA+DFB
HATKAPC=DFC/(DFA/HATKAPA+DFB/HATKAPB)
CALL XIDCH(ALPHA,DFC,XCHI,IER)
IF (IER.NE.0) WRITE(*,*) 'ERROR IN XIDCH'
FKAP1=HATKAPC*XCHI/DFC
CALL XIDCH(ALPHA1,DFC,XCHI,IER)
IF (IER.NE.0) WRITE(*,*) 'ERROR IN XIDCH'
FKAP2=HATKAPC*XCHI/DFC
ENDIF
C
c calculate covariance matrix for c
call transp(a,c)
call mul(covb,a,w)
call mul(c,w,covc)
do 40 i=1,3
do 40 j=1,3
40 covc(i,j)=covc(i,j)+cova(i,j)
c calculate the inverse of covariance matrix H11.2
call inv(covc,w)
c calculate c=b.a
call sumrot(xa,xb,xc)
call mul(b,a,c)
c calculate kappac and dfc if relevant
C hatkapc=1.
C dfc=10000.
C if (hatkapa.eq.1..or.hatkapb.eq.1.) goto 50
C dfc=dfa+dfb
C hatkapc=dfc/(dfa/hatkapa+dfb/hatkapb)c
50 write (*,*) ' Product rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') xc
write(*,*) 'kappa and degrees of freedom'
write(*,'(3x,f5.2,f5.0)') hatkapc,dfc
write (*,*) 'Rotation matrix'
do 55 i=1,3
55 write (*,*) (c(i,j),j=1,3)
write (*,*) 'Covariance matrix'
do 60 i=1,3
60 write (*,*) (covc(i,j),j=1,3)
c
65 write(*,1065)
1065 format(/,' Enter name of the output file: ',$)
read(*,'(a)') filnam
open(unit=15,file=filnam,status='new',err=65)
write(*,'(a)') filnam
write(15,1060) name(1),name(2)
1060 format ('Results from adding ',a,' to ',a,/)
write(15,*) 'Fitted rotation--alat,along,rho: '
write(15,*) xc
write(15,*) 'conf. level, conf. interval for kappa: '
write(15,*) plev,fkap2,fkap1
write(15,*) 'kappahat, degrees of freedom,xchi'
write(15,*) hatkapc,dfc,xchi
write(15,*) 'Number of points, sections'
write(15,*) ndat,nsect
write(15,*) 'ahat: '
do 70 i=1,3
70 write(15,*) (c(i,j),j=1,3)
write(15,*) 'covariance matrix'
do 75 i=1,3
75 write(15,*) (covc(i,j),j=1,3)
write(15,*) 'H11.2 matrix: '
do 80 i=1,3
80 write(15,*) (w(i,j),j=1,3)
IF (IFLAG1.NE.0) WRITE(15,9999)
IF (IFLAG2.NE.0) WRITE(15,9998)
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****************************************************************
c
subroutine sumrot (a,b,c)
c
c finds the total rotation given the two successive rotations
c a(1)=latitude, a(2)=longitude, a(3)=angle
c
implicit double precision (a-h,o-z)
dimension a(3),b(3),c(3)
data rad/.0174532925199432958/
call cdtrn (a(1),a(2),a(3),w1,x1,y1,z1)
call cdtrn (b(1),b(2),b(3),w2,x2,y2,z2)
wt=w1*w2-x1*x2-y1*y2-z1*z2
xt=w1*x2+x1*w2-y1*z2+z1*y2
yt=w1*y2+x1*z2+y1*w2-z1*x2
zt=w1*z2-x1*y2+y1*x2+z1*w2
if(wt) 5,1,1
1 tt=acos(wt)
alt=acos(zt/( sin(tt)))
go to 10
5 tt=acos(-wt)
alt=acos(zt/(-sin(tt)))
10 pt=atan2(yt,xt)
c(1)=90.-alt/rad
gs=pt/rad
if(gs.gt.180.) gs=gs-360.
c(2)=gs
c(3)=tt*2./rad
return
end
c*****************************************************************
c
subroutine cdtrn (f,g,o,w,x,y,z)
implicit double precision (a-h,o-z)
data rad/.0174532925199432958/
pi=rad*180.
al=(90.-f)*rad
p=g*rad
t=o*rad
st=sin(t/2.)
w=cos(t/2.)
x=st*sin(al)*cos(p)
y=st*sin(al)*sin(p)
z=st*cos(al)
if(t-pi) 5,5,1
1 w=-w
x=-x
y=-y
z=-z
5 return
end
c*****************************************************************