c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c program addtj4k.f June, 1995
c
c modified Sept. 1995 declared plev,xf,df to be real*4
c
c This is an adaptation of addrot.f to be used for combining a triple
c junction estimate with another rotation when kappas on the branches
c of the triple junction are unequal. It estimates BA where A is
c estimated at a triple junction. Assumes A rotates
c plate 1 to plate 2.
c
c Critical value for the confidence region construction is calculated
c based on estimates developed by S. Johansen
c
c
C COMPILE WITH NUMLIB.F , numlib2.f , dtrans.f, teclib.f
c
C reads output files from trikap fit of triple junction and output
c from hellingerfit for rotation B as well as data files for both A and B
c
c data files must be formatted as follows:
c
c triple junction data file, file code 10
c first line: for boundary of plates 1 and 2, # of sections, first sec.#,
c estimated kappa
c second line: for boundary of plates 1 and 3, # of sections, first sec.#,
c estimated kappa
c third line: for boundary of plates 2 and 3, # of sections, first sec.#,
c estimated kappa
c fourth line: total number of sections
c data lines follow, each line to represent one crossing and contain
c iside--side number (1,2,or3)
c isect--section number
c crossing latitude
c crossing longitude
c estimated error in crossing location (expressed in kilometers)
c
c
c for second rotation, file code 10
c first line: number of sections
c data lines, one for each crossing
c
c
program addtj4k
implicit double precision (a-h,o-z)
parameter (ndatmx=400,ndtmx2=2*ndatmx)
parameter (msect = 70,msig=2*msect*msect+13*msect+21,msig2=
& 2*msect+6,mwork=2*msig2)
character filnam*70,name(2)*15,datdum*50,datnam*15
dimension sigma(3,msect,3,3),axis(3),w(3,3),xa(3),xb(3),xc(3)
dimension h(3),etai(msect,3,2),ind(2,3),yp(7),num(6),
& qhati(12),eta(msect,3,3),ahat(3,3,3),wt(3,msect,ndatmx)
dimension u(3,msect,ndatmx,3),x2(ndatmx,ndatmx),
& x(ndatmx,ndatmx),npts(3,msect),axm(3,3),
& tempo(ndatmx,ndatmx),b(3,3),c(3,3)
dimension ub(3,msect,ndatmx,3),nptsb(2,msect),numb(2),
& wtb(2,msect,ndatmx),covb(3,3),pone(ndatmx,ndatmx)
dimension np(4),fr(4),cova(3,3),
& xk(ndatmx,ndatmx),p1(4,ndatmx,ndatmx),p0(4,ndatmx,ndatmx)
dimension trp1(5),trp0(4),trp1sq(4),trp0sq(4),trp1p0(4),
& pzero(ndatmx,ndatmx),a(3,3),covc(3,3)
REAL*4 PLEV,XF,df
external R1
Common nsect,sigma,qhati,eta,etai,wt,npts,u
data ind/1,2,1,3,2,3/
DATA rfact/4.0528473E07/,rfact2/6366.1977/,
& nsig/4/,maxfn/1000/,fmin/5000./,fmax/0./
1010 format(A)
1011 format(1x,A)
1012 format(A29,A15)
1013 format(1x,A29,A15)
c
c read in estimates and data for the triple junction rotation
c
ndat=0
nsect=0
do 50 i1 = 1,3
do 50 i2 = 1,msect
npts(i1,i2)=0
do 50 i3 = 1,3
do 50 i4 = 1,3
50 sigma(i1,i2,i3,i4) = 0.
c
1 write(*,1011)'filename for fit of rotation A'
write(*,1011)' from plate 1 to plate 2?'
read (*,1010) name(1)
open (unit=10,file=name(1),status='old',err=1)
write (*,1011) name(1)
read(10,1010) datdum
write(*,1011)datdum
read(10,1010)datdum
read(10,1010) datdum
write(*,1011)datdum
read(10,*)alat,along,rho
write(*,*)alat,along,rho
xa(1)=alat
xa(2)=along
xa(3)=rho
call trans3(alat,along,rho,qhati)
do 15 i = 1,3
15 read(10,1010) datdum
read (10,*) hatkapa,dfa
read (10,1010) datdum
read(10,*) nrowx1,nsect
read(10,1010)datdum
write(*,1011)'rotation matrix'
do 5 i = 1,3
read(10,*)(ahat(i,j,1),j=1,3)
5 write(*,*)(ahat(i,j,1),j=1,3)
write(*,1011)' '
read(10,'(a)') filnam
write (*,1011) filnam
do 6 i = 1,3
read(10,*)(cova(i,j),j=1,3)
6 write(*,*)(cova(i,j),j=1,3)
write(*,*) ' '
close(10)
c
16 write(*,1011)'filename for fit of rotation: plate 1 to plate 3?'
read(*,'(a)')datdum
open(unit=10,file=datdum,status='old',err=16)
write(*,1011)datdum
c read only the rotation and rotation matrix
read(10,'(a)')datdum
write(*,1011)datdum
read(10,1010)datdum
read(10,1010)datdum
write(*,1011) datdum
read(10,*)alat,along,rho
write(*,*)alat, along,rho
call trans3(alat,along,rho,qhati(5))
call qmulti(qhati)
do 7 i = 1,7
7 read (10,'(a)')datdum
write(*,1011)'rotation matrix'
do 8 i= 1,3
read(10,*)(ahat(i,j,2),j=1,3)
8 write(*,*) (ahat(i,j,2),j=1,3)
write(*,*)' '
close(10)
9 write(*,1011)'triple junction data file name?'
read(*,'(a)')filnam
write(*,1011)filnam
open(unit=10,file=filnam,status='old',err=9)
rewind 10
read(10,*)nsect1,nfsect1,hatkap1
read(10,*)nsect2,nfsect2,hatkap2
read(10,*)nsect3,nfsect3,hatkap3
read(10,*)nsect
do 110 k = 1, nrowx1
13 read(10,*,end=250,err=13) iside,isect,alat,along,sd
if (iside .gt. 3) goto 13
npts(iside,isect)=npts(iside,isect)+1
wt(iside,isect,npts(iside,isect))=1/sd
call trans1(alat,along,axis)
n=npts(iside,isect)
do 105 i = 1,3
105 u(iside,isect,n,i)=axis(i)
do 110 j1 = 1,3
do 110 j2 = 1,3
110 sigma(iside,isect,j1,j2)=sigma(iside,isect,j1,j2)
& +axis(j1)*axis(j2)*1/sd**2
250 do 115 i = 1,6
115 num(i)=0
do 120 j = nfsect1,nfsect1+nsect1 - 1
num(1) = num(1) + npts(2,j)
120 num(2) =num(2) + npts(1,j)
do 125 j = nfsect2,nfsect2+nsect2-1
num(3) = num(3) + npts(3,j)
125 num(4) = num(4) + npts(1,j)
do 130 j = nfsect3,nfsect3+nsect3-1
num(5)=num(5) + npts(2,j)
130 num(6)=num(6) + npts(3,j)
close(10)
c
c
n1=num(1)+num(2)
n2=num(3)+num(4)
n3=num(5)+num(6)
fr(1)=n1-3-2*nsect1
fr(2)=n2-3-2*nsect2
fr(3)=n3-3-2*nsect3
c
c set up eta vectors
c
do 210 i = 1,6
210 yp(i)=0.
rmin=r2(yp)*rfact
c
c set up first part of design matrix
c
ncolx1 = 6 + 2*nsect
call trimat(nrowx1,ncolx1,ahat,num,nsect1,nfsect1,nsect2,nfsect2,
& nsect3,nfsect3,x)
c
c read fit file for rotation B
c
2 write(*,*)' '
write(*,1011)'enter file name for fit of second rotation, B'
read(*,1010)name(2)
open(unit=10,file=name(2),status='old',err=2)
rewind 10
write (*,1011)name(2)
read(10,1012)datdum, datnam
write(*,1013)datdum,datnam
read(10,1010)filnam
write(*,1010)filnam
c read alat alon,angle
read(10,*)xb
write(*,*)xb
read(10,1010)filnam
read(10,1010)filnam
read(10,1010)filnam
read(10,*)hatkapb,dfb
read(10,1010)filnam
read(10,*)nrowx2,ls
read(10,1010)filnam
c read rotation matrix
do 305 i = 1,3
305 read(10,*)(b(i,j),j=1,3)
read(10,1010)filnam
c read covariance matrix
do 310 i = 1,3
310 read(10,*)(covb(i,j),j=1,3)
write(*,*)' '
close(10)
c
c read data file for rotation B
c
OPEN(10,FILE=datNAM,STATUS='UNKNOWN')
REWIND 10
write(*,*)' '
READ(10,*) nsectb
NUMb(1)=0
NUMb(2)=0
do 51 I1 = 1,2
do 51 I2 = 1,nsectb
nptsb(i1,i2) = 0
do 51 i3 = 1,3
do 51 i4 = 1,3
51 sigma(i1,i2,i3,i4) = 0.
do 200 ktimes = 1,nrowx2
100 READ(10,*,err=100) ISIDE,ISECT,ALAT,ALONG,SD
nptsb(iside,isect) = nptsb(iside,isect) + 1
CALL TRANS1(ALAT,ALONG,AXIS)
n = nptsb(iside,isect)
wtb(iside,isect,n) = 1/sd
do 151 i = 1,3
151 ub(iside,isect,n ,i) = axis(i)
NUMb(ISIDE)=Numb(iside) + 1
do 114 j1 = 1,3
do 114 j2=1,3
114 sigma(iside,isect,j1,j2)=sigma(iside,isect,j1,j2)
& + axis(j1)*axis(j2)*1/(sd**2)
200 continue
201 close (10)
c
c setting up the eta vectors
c
alat = xb(1)
along = xb(2)
rho = xb(3)
CALL TRANS3(ALAT,ALONG,RHO,QHATI)
DO 207 I=1,3
207 H(I)=0.
RMIN=R1(H)*RFACT
c
C
C ETA IS THE MATRIX M(ETA) OF PAPER; ETAI IS THE MATRIX O-SUB I.
C ETA AND ETAI ARE SET BY R1.
C
c
c construct estimate of design matrix X2 (for rotation B), beginning
c with first 3 columns.
c
nrowx = nrowx2
ncolx = 3 + 2*nsectb
do 505 k = 1, nrowx
do 505 j = 1, ncolx
505 x2(k,j) = 0.
nrow = 1
do 610 isect = 1, nsectb
do 570 k = 1,3
do 570 j = 1,3
570 axm(k,j) = 0.
do 580 k = 1,3
do 580 j = 1,3
do 580 m = 1,3
580 axm(k,j) = axm(k,j) + b(k,m)*eta(isect,j,m)
do 600 n = 1,nptsb(2,isect)
do 595 k = 1,3
do 590 j = 1,3
590 x2(nrow,k) = x2(nrow,k) + ub(2,isect,n,j)*axm(j,k)
595 x2(nrow,k) = x2(nrow,k)*wtb(2,isect,n)*rfact2
600 nrow = nrow + 1
610 continue
nr = 1
nc = 4
do 650 isect = 1,nsectb
do 620 k = 1,3
do 620 j = 1,2
620 axm(k,j) = 0.
do 630 k = 1,3
do 630 j = 1,2
do 630 m = 1,3
630 axm(k,j) = axm(k,j) + b(k,m)*etai(isect,m,j)
do 645 n = 1,nptsb(2,isect)
do 640 k = 1,3
x2(nr,nc) = x2(nr,nc) + ub(2,isect,n,k)*axm(k,1)
640 x2(nr,nc+1) = x2(nr,nc+1) + ub(2,isect,n,k)*axm(k,2)
x2(nr,nc) = wtb(2,isect,n)*x2(nr,nc)*rfact2
x2(nr,nc+1) = wtb(2,isect,n)*x2(nr,nc+1)*rfact2
645 nr = nr + 1
650 nc = nc + 2
nc = 4
do 680 isect = 1,nsectb
do 670 n = 1,nptsb(1,isect)
do 660 k = 1,3
x2(nr,nc) = x2(nr,nc) + ub(1,isect,n,k)*etai(isect,k,1)
660 x2(nr,nc+1) = x2(nr,nc+1) +
& ub(1,isect,n,k)*etai(isect,k,2)
x2(nr,nc) = x2(nr,nc)*wtb(1,isect,n)*rfact2
x2(nr,nc+1) = x2(nr,nc+1)*wtb(1,isect,n)*rfact2
670 nr = nr + 1
680 nc = nc + 2
close(10)
c
c
c Construct estimated design matrix X from the separate fit matrices
c
ncolx2=3+2*nsectb
nrowx = nrowx1 + nrowx2
if (nrowx .gt. ndatmx) write(6,*)'STOP. number of data points',
& 'exceeds ndatmx. Parameter statement must be increased.'
ncolx = ncolx1 + ncolx2
if (ncolx .gt. msig2) write(6,*)'STOP. msect must be increased.'
do 32 i = 1,nrowx2
do 32 j = 1,ncolx2
32 x(i+nrowx1,j+ncolx1)=x2(i,j)
do 33 i = 1,nrowx1
do 33 j = ncolx1+1,ncolx
33 x(i,j) = 0.0
do 34 i = nrowx1+1,nrowx
do 34 j = 1,ncolx1
34 x(i,j) = 0.0
c
c transform X to S^(-1/2)*X
c
rkap1=hatkap1**.5
rkap2=hatkap2**.5
rkap3=hatkap3**.5
rkapb=hatkapb**.5
do 303 i = 1,n1
do 303 j = 1,ncolx1
303 x(i,j)=rkap1*x(i,j)
do 306 i = n1+1,n1+n2
do 306 j = 1, ncolx1
306 x(i,j)=rkap2*x(i,j)
do 307 i = n1 + n2 + 1,nrowx1
do 307 j = 1,ncolx1
307 x(i,j) = rkap3*x(i,j)
do 309 i = nrowx1+1,nrowx
do 309 j = 1,ncolx1+1,ncolx
309 x(i,j)=rkapb*x(i,j)
c
c
c construct pzero
c
call projection(x,nrowx,ncolx,pzero)
c
c construct XK
c
ncolxk = ncolx - 3
do 320 i = 1,nrowx1
do 320 j = 1,ncolxk
320 xk(i,j) = x(i,j)
do 360 i = nrowx1+1,nrowx
do 330 j = 1,3
330 xk(i,j)=-x(i,ncolx1+j)
do 340 j = 4,ncolx1
340 xk(i,j) = 0.0
do 350 j = ncolx1 + 1,ncolxk
350 xk(i,j)=x(i,j+3)
360 continue
c
c multiply columns 1-3, rows nrowx1+1 and down,on the right by A
c
do 368 i = nrowx1 + 1, nrowx
do 364 j=1,3
364 tempo(i,j)=0.0
do 366 j = 1,3
do 366 k = 1,3
366 tempo(i,j)=tempo(i,j) + xk(i,k)*ahat(k,j,1)
do 367 j = 1,3
367 xk(i,j)=tempo(i,j)
368 continue
c open(13,file='xk.out',status='unknown')
c do 380 i = 1,nrowx
c write(13,*)(xk(i,j),j = 1,ncolxk)
c 380 write (13,*)' '
c close(13)
c
c
c calculate P-one
c
call projection(xk,nrowx,ncolxk,pone)
c
c open(unit=15,file='pone.dat',status='unknown')
c do 777 i = 1,nrowx
c 777 write(15,*)(pone(i,j),j=1,nrowx)
c close(15)
c
c Calculation of correction terms A,B and C
c
fr(4)=nrowx2-3-2*nsectb
do 802 i = 1,4
if (fr(i) .lt. fmin) fmin=fr(i)
802 if (fr(i) .gt. fmax) fmax=fr(i)
np(1)=n1
np(2)=n2
np(3)=n3
np(4)=nrowx2
c
nrowdone=0
do 800 k = 1,4
do 805 i = 1,np(k)
do 805 j = 1,np(k)
p1(k,i,j)=pone(i+nrowdone,j+nrowdone)
805 p0(k,i,j)=pzero(i+nrowdone,j+nrowdone)
800 nrowdone=nrowdone+np(k)
c
do 820 k = 1,4
trp1(k)=0.0
trp0(k)=0.0
trp1sq(k)=0.0
trp0sq(k)=0.0
trp1p0(k)=0.0
do 820 i = 1,np(k)
trp1(k)=trp1(k)+p1(k,i,i)
trp0(k)=trp0(k)+p0(k,i,i)
do 825 m = 1,np(k)
trp1sq(k)=trp1sq(k)+p1(k,i,m)*p1(k,m,i)
trp0sq(k)=trp0sq(k)+p0(k,i,m)*p0(k,m,i)
825 trp1p0(k)=trp1p0(k)+p1(k,i,m)*p0(k,m,i)
820 continue
c
c
capa=0.
capb=0.
capc=0.
do 830 k=1,4
capa=capa+(trp0(k)-trp1(k)-trp1p0(k)+trp1sq(k))/fr(k)
capb=capb+(trp0(k)-trp1(k)-trp0sq(k)+trp1p0(k))/fr(k)
830 capc=capc+((trp0(k)-trp1(k))**2-trp0sq(k)+2*trp1p0(k)-
& trp1sq(k))/fr(k)
write(*,1011)'results for combined rotation'
write(*,*)' '
write(6,*)'A,B,C',capa,capb,capc
c
c calculate c and f, where Q is cF(3,f)
c
exq = 3 + 2*capa + 2*capb
vq = 6 + 14*capa + 2*capc + 2*capb
write (6,*) 'mean of Q = ', exq
write (6,*) 'variance of Q = ', vq
ex2 = exq*exq
write (6,*)
df = (12*vq + 2*ex2)/(3*vq - 2*ex2)
cfact = exq*(1 - 2/df)
write (6,*) 'deg. freedom ', df
write (6,*) ' c = ',cfact
c write(6,*)'enter plev'
c read(6,*) plev
plev = .95
call xidf(plev,3.0,df,xf,ier)
cxf = cfact*xf
write(6,*)'xf, critical value= c*f(3,df) :',xf,cxf
write(6,*)' '
c
DO 36 I=1,3
DO 36 J=1,3
a(i,j)=ahat(i,j,1)
COVA(I,J)=COVA(I,J)/HATKAPA
36 COVB(I,J)=COVB(I,J)/hatkapb
c
c adjust hatkapc for use in conreg.f
HATKAPC = 3./cfact
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
c
call sumrot(xa,xb,xc)
call mul(b,a,c)
c
c
write (*,*) ' Product rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') xc
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,*)'plev'
write(15,*) plev,1.0,1.0
c flag = 1 is output for use by the prgram addplus.f
write(15,*) '3/cfact, degrees of freedom, crit. value, flag'
write(15,'(f12.6,2x,f8.4,2x,f8.4,2x,f3.1)') hatkapc,df,cxf,1.0
write(15,*) 'Number of points, fmin, fmax, rank'
write(15,*) nrowx,fmin,fmax,ncolx
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)
close (15)
c
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
c
C
C
FUNCTION R1(H)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER (MSECT=70)
DIMENSION H(3),SIGMA(2,MSECT,3,3),QHATI(4),ETA(MSECT,3,3),
& ETAI(MSECT,3,2),AHAT(3,3),SIG(3,3),D(3),Z(3,3),WK(6),QHAT(4)
COMMON NSECT,SIGMA,QHATI,ETA,ETAI
C
CALL TRANS2(H,QHATI,QHAT)
CALL TRANS4(QHAT,AHAT)
R1=0.
DO 100 I=1,NSECT
DO 110 J=1,3
DO 110 K=1,3
SIG(J,K)=SIGMA(1,I,J,K)
DO 110 K1=1,3
DO 110 K2=1,3
110 SIG(J,K)=SIG(J,K)+AHAT(K1,J)*SIGMA(2,I,K1,K2)*AHAT(K2,K)
CALL JACOBI(SIG,3,3,D,Z,3,WK,NROT)
IF (NROT.LT.0) WRITE(6,*) 'SUBROUTINE JACOBI(3)--NROT: ',NROT
ETA(I,1,1)=0.
ETA(I,2,1)=Z(3,1)
ETA(I,3,1)=-Z(2,1)
ETA(I,1,2)=-Z(3,1)
ETA(I,2,2)=0.
ETA(I,3,2)=Z(1,1)
ETA(I,1,3)=Z(2,1)
ETA(I,2,3)=-Z(1,1)
ETA(I,3,3)=0.
IF (ABS(Z(3,1)).GT.(.2)) THEN
DO 120 J=1,3
ETAI(I,J,1)=ETA(I,J,1)
120 ETAI(I,J,2)=ETA(I,J,2)
ELSE IF (ABS(Z(1,1)).GT.(.2)) THEN
DO 125 J=1,3
ETAI(I,J,1)=ETA(I,J,2)
125 ETAI(I,J,2)=ETA(I,J,3)
ELSE
DO 130 J=1,3
ETAI(I,J,1)=ETA(I,J,1)
130 ETAI(I,J,2)=ETA(I,J,3)
ENDIF
100 R1=R1+D(1)
RETURN
END
C