c****************hellinger3.f******************************************
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c modified Sept. 1995 real declaration changed to real*4
c minor modification june 4,1993 flag added to output on device 15 to
c be used by program addplus.f
c
c last modified november 1, 1988
C LAST MODIFIED MARCH 11, 1991: CORRECT ERROR IN CALCULATION OF
C COVARIANCE BETWEEN ROTATION 2-3 AND THE ROTATIONS 1-2 AND 1-3.
C NOTE: OLD CALCULATION OF COVARIANCE MATRIX FOR ROTATION 2-3 WITH
C ITSELF, ROTATION 1-2 WITH ITSELF, ROTATION 1-3 WITH ITSELF, AND
C ROTATION 1-2 WITH ROTATION 1-3 ARE CORRECT. CHANGED LINES IN CAPS.
C
C MINOR MODIFICATION OCTOBER 14, 1991: MODIFICATION IN CAPS
C
C NOTE: VERSIONS OF THIS PROGRAM DATED OCT. 31, 1988 AND EARLIER HAVE A
C SERIOUS BUG. THIS BUG WAS CORRECTED IN VERSION OF NOV. 1, 1988.
c
c triple junction estimation program
c reference: chang, t. estimating the relative rotation of two tectonic
c plates from boundary crossings, j. amer. stat. assn., scheduled to
c appear dec. 1988.
c
c input data file: file code 10
c first line is number of sections
c data lines follow: each line to represent one crossing and contain
c iside--side number (either 1, 2, or 3)
c isect--section number
c crossing latitude
c crossing longitude
c estimated error in crossing location (expressed in kilometers)
c
c this program is to be compiled with the following libraries
c dtrans--utility routines
c bingham--program to produce grid matrices for use in contouring
c confidence region
c numlib--numerical libraries
c
c this program simultaneously estimates the rotations of side 1 into
c side 2, side 1 into side 3, and side 2 into side 3. a
c confidence region for each rotation is produced in a form suitable for
c mapping using surface2. the output files (for each rotation) are
c 30--file of points on the boundary of the region of possible axes
c 31--upper surface grid matrix
c 32--lower surface grid matrix
c
c starting with an initial guess of the optimal rotations, the program
c uses subroutine amoeba, a downhill simplex method, for
c minimization. it is recommended that amoeba be called at least twice
c and at least enough times to produce a minimum which is different from
c the initial guess.
c
c amoeba works with rotations using a quaternion representation. qhati
c is the initial guess to amoeba, with the rotation side 1 to side 2
c stored in the first 4 entries, side 1 to side 3 in entries 5 to 8, and
c side 2 to side 3 in entries 9 to 12. amoeba then parameterizes the
c candidate rotation qhat using a vector h of length 6 defined by
c qhat(1:4)=qhati(1:4)*exp(h(1:3))
c qhat(5:8)=qhati(5:8)*exp(h(4:6))
c qhat(9:12)=qhat(5:8)*qhatbar(1:4)
c here * is quaternion multiplication, qhatbar is the quaternionic
c conjugate of qhat, and exp is the map which takes a 3 vector h into
c right hand rule rotation of |h| radians around the axis h/|h|. amoeba
c has not changed qhati exactly when h(1:6)=0.
c ftol is the estimated relative difference between the value of r at
c qhat and its value at the true minimum. the estimated relative error
c in qhat is thus sqrt(ftol).
c
c***********************************************************************
c
program Hellinger3
implicit double precision(a-h,o-z)
parameter (msect=50,msig=2*msect*msect+13*msect+21,msig2=2*msect+6
&, mwork=2*msig2)
c msect must be at least 2 or workspace will be too small!
real*4 xchi,plev,df,plev1,qbing(10),axisr(3),xchi1
dimension ahatm(3,4),sigma(3,msect,3,3),axis(3),qhat(12),
& hp(6,7),yp(7),work(mwork),qhati(12),eta(msect,3,3),ahat(3,3,3),
& temp(3,3),d(msig2),z(msig2,msig2),etai(msect,3,2),sigma2(msig),
& sigma3(msig2,msig2),ind(2,3),sigma4(3,3),sigma5(9,9)
character filnam*20,filnam2*20,f1*1,f2*1,ans*1
external r2
common nsect,sigma,qhati,eta,etai
data ind/1,2,1,3,2,3/
data rfact/4.0528473e07/,nsig/4/,maxfn/1000/
DATA HATKAP/0./,DF/0./,FKAP1/0./,FKAP2/0./
1000 format(a)
c
c calculation of the matrices sigma and sigma-tilde of reference
c
do 50 i1=1,3
do 50 i2=1,msect
do 50 i3=1,3
do 50 i4=1,3
50 sigma(i1,i2,i3,i4)=0.
c
c
write(6,*) 'data file name? '
read(5,1000) filnam
open(10,file=filnam,status='unknown')
read(10,*) nsect
ndat=0
100 read(10,*,end=250,err=250) iside,isect,alat,along,sd
ndat=ndat+1
call trans1(alat,along,axis)
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)/(sd**2)
c & +axis(j1)*axis(j2)/(sd**2)
go to 100
c
c minimization section
c
250 write(6,*) 'initial guess, side 1 to side 2: alat, along, rho? '
read(5,*) alat,along,rho
call trans3(alat,along,rho,qhati)
write(6,*) 'initial guess, side 1 to side 3: alat, along, rho? '
read(5,*) alat,along,rho
call trans3(alat,along,rho,qhati(5))
call qmulti(qhati)
write(6,*) 'qhat: ',qhati
eps=.2
c
260 do 215 j=1,7
do 210 i=1,6
210 hp(i,j)=0.
if (j.ne.1) hp(j-1,j)=eps
215 yp(j)=r2(hp(1,j))
rmin=yp(1)*rfact
ftol=(.1)**(2*nsig+1)
iter=maxfn
write(6,*)
write(6,*) 'calling amoeba--r = ',rmin
rprev=rmin
call amoeba(hp,yp,6,6,ftol,r2,iter,work)
rmin=r2(hp(1,1))*rfact
call trans2(hp(1,1),qhati,qhat)
call trans2(hp(4,1),qhati(5),qhat(5))
call qmulti(qhat)
write(6,*) 'return from amoeba--r = ',rmin
write(6,*) 'r improvement: ',rprev-rmin
write(6,*) 'h: ',(hp(i,1),i=1,3)
write(6,*) 'qhat: ',qhat
write(6,*) 'ftol, iter: ',ftol,iter
write(6,*) 'do you want to call amoeba again? '
read(5,'(a)') ans
if ((ans.eq.'Y').or.(ans.eq.'y')) then
eps=eps/20.
do 212 i=1,12
212 qhati(i)=qhat(i)
go to 260
endif
do 211 k=1,3
write(6,*) 'fitted rotation side ',ind(1,k),
& ' to side ',ind(2,k),'--alat,along,rho: '
k1=4*k-3
call trans5(qhat(k1),alat,along,rho)
call trans4(qhat(k1),ahat(1,1,k))
write(6,*) alat,along,rho
write(6,*) 'ahat: '
do 211 i=1,3
211 write(6,*) (ahat(i,j,k),j=1,3)
c
c calculation of the matrix (x**t)*x. sigma2 stores this matrix in
c symmetric mode.
c order of columns (rows) of sigma2 is: rotation side 1 to 2
c (3 indices), rotation side 1 to 3 (3 indices), normal vectors
c (2*nsect indices)
c at this point sigma is no longer needed. convert sigma to
c (ahat**t)*sigma*ahat.
c eta is the matrix m(eta) of paper; etai is the matrix o-sub i.
c eta and etai are set by r2.
c
219 do 230 ks=2,3
ks1=ks-1
do 230 isect=1,nsect
do 220 i=1,3
do 220 j=1,3
temp(i,j)=0.
do 220 k1=1,3
do 220 k2=1,3
220 temp(i,j)=temp(i,j)+
& ahat(k1,i,ks1)*sigma(ks,isect,k1,k2)*ahat(k2,j,ks1)
do 230 i=1,3
do 230 j=1,3
230 sigma(ks,isect,i,j)=temp(i,j)
c
do 300 i=1,msig
300 sigma2(i)=0.
k=0
c lower triangular part of 3 by 3 block for rotation side 1 to 2
do 301 i=1,3
do 301 j=1,i
k=k+1
do 301 isect=1,nsect
do 301 k1=1,3
do 301 k2=1,3
301 sigma2(k)=sigma2(k)+
& eta(isect,i,k1)*sigma(2,isect,k1,k2)*eta(isect,j,k2)
c 3 by 3 block of zeros followed by lower triangular part of
c 3 by 3 block for rotation side 1 to 3
do 305 i=1,3
k=k+3
do 305 j=1,i
k=k+1
do 305 isect=1,nsect
do 305 k1=1,3
do 305 k2=1,3
305 sigma2(k)=sigma2(k)+
& eta(isect,i,k1)*sigma(3,isect,k1,k2)*eta(isect,j,k2)
c for each section insert 2 by 3 block for the normal by
c rotation from side 1 to 2, followed by a 2 by 3 block for the
c normal by rotation from side 1 to 3
do 302 isect=1,nsect
do 302 i=1,2
do 303 ks=2,3
do 303 j=1,3
k=k+1
do 303 k1=1,3
do 303 k2=1,3
303 sigma2(k)=sigma2(k)+
& etai(isect,k1,i)*sigma(ks,isect,k1,k2)*eta(isect,j,k2)
c skip over (isect-1) 2 by 2 blocks of zeros
k=k+2*(isect-1)
c 2 by 2 block for normal by normal
do 302 j=1,i
k=k+1
do 302 k1=1,3
do 302 k2=1,3
302 sigma2(k)=sigma2(k)+etai(isect,k1,i)*(sigma(1,isect,k1,k2)+
& sigma(2,isect,k1,k2)+sigma(3,isect,k1,k2))*etai(isect,k2,j)
c
ndim=2*nsect+6
k=0
do 304 i=1,ndim
do 304 j=1,i
k=k+1
sigma3(i,j)=sigma2(k)
304 sigma3(j,i)=sigma2(k)
c
c
c calculation of inverse sigma3 to sigma2. this is full variance
c covariance matrix (except for factors of kappahat and rfact) for
c rotation and normal parameters
c
call jacobi(sigma3,ndim,msig2,d,z,msig2,work,nrot)
if (nrot.lt.0) then
write(6,*) 'subroutine jacobi(1)--nrot ',nrot
write(6,*) 'ndim,k: ',ndim,k
write(6,*) 'sigma3: ',(sigma2(i),i=1,k)
do 309 j=1,ndim
write(6,*) 'eigenvalue: ',d(j)
309 write(6,*) 'eigenvector: ',(z(i,j),i=1,ndim)
endif
do 310 i=1,ndim
do 310 j=1,ndim
sigma3(i,j)=0.
do 310 k=1,ndim
310 sigma3(i,j)=sigma3(i,j)+z(i,k)*z(j,k)/d(k)
c
c calculation of covariance matrix sigma5 (except possibly for division
c by hatkap) of the 9 rotation parameters. this matrix is necessarily
c singular.
c
do 320 i=1,6
do 320 j=1,6
320 sigma5(i,j)=sigma3(i,j)/rfact
do 321 i=1,3
i3=i+3
i6=i+6
do 321 j=1,3
j3=j+3
j6=j+6
sigma5(i6,j)=0.
sigma5(i6,j3)=0.
sigma5(i6,j6)=0.
do 322 k1=1,3
k13=k1+3
SIGMA5(I6,J)=SIGMA5(I6,J)+AHAT(I,K1,1)*(SIGMA5(K13,J)-
& SIGMA5(K1,J))
SIGMA5(I6,J3)=SIGMA5(I6,J3)+AHAT(I,K1,1)*(SIGMA5(K13,J3)-
& SIGMA5(K1,J3))
do 322 k2=1,3
k23=k2+3
C (MARCH 11, 1991): NEXT TWO LINES HAVE BEEN CORRECTED AND PLACED ABOVE
C sigma5(i6,j)=sigma5(i6,j)+ahat(i,k1,1)*(sigma5(k13,k2)-
C & sigma5(k1,k2))*ahat(j,k2,1)
C sigma5(i6,j3)=sigma5(i6,j3)+ahat(i,k1,1)*(sigma5(k13,k23)-
C & sigma5(k1,k23))*ahat(j,k2,1)
322 sigma5(i6,j6)=sigma5(i6,j6)+ahat(i,k1,1)*(sigma5(k1,k2)+
& sigma5(k13,k23)-sigma5(k1,k23)-sigma5(k13,k2))*ahat(j,k2,1)
sigma5(i,j6)=sigma5(j6,i)
321 sigma5(i3,j6)=sigma5(j6,i3)
write(6,*)
write(6,*) 'sigma5: '
write(6,*) ((sigma5(i,j),j=1,i),i=1,9)
c
c calculation of critical points: xchi for simultaneous regions;
c xchi1 for individual regions.
c
write(6,*)
write(6,*) 'enter significance level: '
read(5,*) plev
write(6,*)
write(6,*) 'do you want an estimated kappa? '
read(5,1000) ans
if ((ans.eq.'n').or.(ans.eq.'n')) then
call xidch(plev,6.,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidch--ier ',ier
call xidch(plev,3.,xchi1,ier)
if (ier.ne.0) write(6,*) 'subroutine xidch--ier ',ier
else
df=ndat-2*nsect-6
plev1=(1.-plev)/2
call xidch(plev1,df,xchi1,ier)
if (ier.ne.0) write(6,*) 'subroutine xidch--ier: ',ier
plev1=1-plev1
call xidch(plev1,df,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidch--ier: ',ier
fkap1=xchi/rmin
fkap2=xchi1/rmin
write(6,*) plev,' confidence interval for kappa'
write(6,*) fkap2,fkap1
hatkap=df/rmin
write(6,*) 'kappahat, degrees of freedom: ',hatkap,df
call xidf(plev,6.,df,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidf--ier: ',ier
xchi=xchi*rmin*6./df
call xidf(plev,3.,df,xchi1,ier)
if (ier.ne.0) write(6,*) 'subroutine xidf--ier: ',ier
xchi1=xchi1*rmin*3./df
endif
c
c plotting of confidence regions using bingham. approach is outlined in
c section iv of hanna and chang: on graphically representing the
c confidence region for an unknown rotation in three dimensions.
c
c sigma4 will store the matrix formerly called h11.2.
c
do 480 kr=1,3
write(6,*)
write(6,'(72(1h=))')
write(6,*) 'Plotting simultaneous confidence region for rotation',
& 'side ',ind(1,kr),' to ',ind(2,kr)
kk=3*(kr-1)
do 461 i=1,3
do 461 j=1,3
461 temp(i,j)=sigma5(i+kk,j+kk)
call jacobi(temp,3,3,d,z,msig2,work,nrot)
if (nrot.lt.0) write(6,*) 'subroutine jacobi(4)--nrot,kr ',nrot,kr
do 462 i=1,3
do 462 j=1,3
sigma4(i,j)=0.
do 462 k=1,3
462 sigma4(i,j)=sigma4(i,j)+z(i,k)*z(j,k)/d(k)
do 464 i=1,3
464 d(i)=1./d(i)
write(6,*) 'matrix h11.2: '
write(6,*) ((sigma4(i,j),j=1,i),i=1,3)
write(6,*) 'eigenvalues: '
write(6,*) (d(i),i=1,3)
write(6,*) 'eigenvectors: '
do 463 j=1,3
463 write(6,*) (z(i,j),i=1,3)
c
kr2=(kr-1)*4
ahatm(1,1)=-qhat(2+kr2)
ahatm(1,2)=qhat(1+kr2)
ahatm(1,3)=qhat(4+kr2)
ahatm(1,4)=-qhat(3+kr2)
ahatm(2,1)=-qhat(3+kr2)
ahatm(2,2)=-qhat(4+kr2)
ahatm(2,3)=qhat(1+kr2)
ahatm(2,4)=qhat(2+kr2)
ahatm(3,1)=-qhat(4+kr2)
ahatm(3,2)=qhat(3+kr2)
ahatm(3,3)=-qhat(2+kr2)
ahatm(3,4)=qhat(1+kr2)
k=0
do 451 i=1,4
do 451 j=1,i
k=k+1
qbing(k)=0.
do 451 k1=1,3
do 451 k2=1,3
451 qbing(k)=qbing(k)+4.*ahatm(k1,i)*sigma4(k1,k2)*ahatm(k2,j)
kr1=4*kr-3
call trans5(qhat(kr1),alat,along,rho)
call trans1(alat,along,axis)
do 454 i=1,3
454 axisr(i)=axis(i)
call bingham(qbing,0.,xchi,axisr,jer)
write(6,*)
write(6,*) 'Plotting individual confidence region for rotation',
& 'side ',ind(1,kr),' to ',ind(2,kr)
call bingham(qbing,0.,xchi1,axisr,jer)
c
write(f1,'(i1)') ind(1,kr)
write(f2,'(i1)') ind(2,kr)
filnam2='tj_par'//f1//f2
open (unit=15,file=filnam2,status='unknown')
write(15,1460) filnam,ind(1,kr),ind(2,kr)
1460 format ('Results from Hellinger3 using ',a,/,
&'===== Rotation from side ',i1,' to ',i1)
write(15,*) 'Fitted rotation--alat,along,rho: '
write(15,*) alat,along,rho
write(15,*) 'conf. level, conf. interval for kappa: '
write(15,*) plev,fkap2,fkap1
if (hatkap.eq.0.) hatkap=1.
if (df.eq.0.) df=10000.
c flag is for use by addplus.f
write(15,*) 'kappahat, degrees of freedom,xchi,flag'
write(15,'(f10.5,2x,f4.0,2x,f10.5,2x,f3.1)') hatkap,df,xchi1,0.0
write(15,*) 'Number of points, sections'
write(15,*) ndat,nsect
write(15,*) 'ahat: '
do 460 i=1,3
460 write(15,*) (ahat(i,j,kr),j=1,3)
write(15,*) 'covariance matrix'
do 465 i=1,3
465 write(15,*) (sigma5(i+kk,j+kk),j=1,3)
write(15,*) 'H11.2 matrix: '
do 470 i=1,3
470 write(15,*) (sigma4(i,j),j=1,3)
close (15)
c
480 continue
c
stop
end
c
c
function r2(h)
c calculation of hellinger criterion r for a triple junction.
implicit double precision(a-h,o-z)
parameter (msect=50)
dimension h(6),sigma(3,msect,3,3),qhati(12),eta(msect,3,3),
& etai(msect,3,2),ahat(3,3,2),sig(3,3),d(3),z(3,3),wk(6),qhat(8)
common nsect,sigma,qhati,eta,etai
c
call trans2(h,qhati,qhat)
call trans4(qhat,ahat(1,1,1))
call trans2(h(4),qhati(5),qhat(5))
call trans4(qhat(5),ahat(1,1,2))
r2=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 ks=2,3
ks1=ks-1
do 110 k1=1,3
do 110 k2=1,3
110 sig(j,k)=sig(j,k)+ahat(k1,j,ks1)*sigma(ks,i,k1,k2)*ahat(k2,k,ks1)
call jacobi(sig,3,3,d,z,3,wk,nrot)
if (nrot.lt.0) then
write(6,*) 'subroutine jacobi(3)--nrot: ',nrot
write(6,*) 'vector h: ',h
endif
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 r2=r2+d(1)
return
end
c
c
subroutine qmulti(q)
implicit double precision(a-h,o-z)
dimension q(12),qbar(4)
do 10 i=2,4
10 qbar(i)=-q(i)
qbar(1)=q(1)
call qmult(q(5),qbar,q(9))
return
end