c***************trikap.f************************************************
c23456789012345678901234567890123456789012345678901234567890123456789012
c June, 1993
c
c modified Sept. '95 declared plev,df,xchi as real*4
c modified Dec, 1994 and June '95 moved subroutines to teclib
c
c triple junction estimation program, an adaptation of hellinger3.f to be
c used when kappas are unequal. Assumes independent estimates of kappas
c have been obtained. Sections along any boundary must be numbered
c consecutively.
c
c
c
c input data file: file code 10
c first line: for boundary of plates 1 and 2, # of sections, first section
c number, kappahat(1)
c second line: for boundary of plates 1 and 3, # of sections, first sect. #,
c kappahat(2)
c third line: for boundary of plates 2 and 3, # of secs, first sect. #,
c kappahat(3)
c fifth 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 numlib--numerical libraries
c numlib2--LAPACK routines
c teclib--design matrix routines, function r2, etc.
c CAUTION: changes in the parameters ndatmx and msect MUST also be made
c in all subroutines in teclib
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.
c Output from this program can be used by 'conreg' to generate files suitable
c for mapping the confidence regions
c
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 trikap
implicit double precision(a-h,o-z)
parameter (msect=70,msig=2*msect*msect+13*msect+21,msig2=2*msect+6
&, mwork=2*msig2,ndatmx=400)
c msect must be at least 2 or workspace will be too small!
real*4 xchi,plev,df,plev1,xchi1
dimension 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),
& etai(msect,3,2),sigma3(msig2,msig2),ind(2,3),sigma5(9,9)
dimension npts(3,msect),num(6),u(3,msect,ndatmx,3),b(3,3),
& wt(3,msect,ndatmx),rkap(3),hatkappa(3),
& x(ndatmx,ndatmx),xt(ndatmx,ndatmx),xtx(ndatmx,ndatmx)
dimension tempo(ndatmx,ndatmx),pzero(ndatmx,ndatmx),
& xk(ndatmx,ndatmx),p1(ndatmx,ndatmx),a(3,3),cov(3,3),covinv(3,3)
character filnam*20,ans*1
external r2
common nsect,sigma,qhati,eta,etai,wt,npts,u,fr1,fr2,fr3,n1,n2,n3
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./,rfact2/6366.1977/
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
npts(i1,i2)=0
do 50 i3=1,3
do 50 i4=1,3
50 sigma(i1,i2,i3,i4)=0.
c
write(6,*) 'data file name? '
read(5,1000) filnam
open(10,file=filnam,status='unknown')
read(10,*) nsect1, nfsect1,hatkappa(1)
read(10,*) nsect2, nfsect2,hatkappa(2)
read(10,*) nsect3, nfsect3,hatkappa(3)
read(10,*) nsect
do 80 i = 1,3
80 rkap(i)=hatkappa(i)**.5
ndat=0
100 read(10,*,end=250,err=100) iside,isect,alat,along,sd
if (iside .gt. 3) goto 100
ndat=ndat+1
npts(iside,isect)=npts(iside,isect)+1
if ((isect.ge.nfsect1).and.(isect.lt.nfsect1+nsect1)) then
wt(iside,isect,npts(iside,isect))=rkap(1)/sd
else
if ((isect.ge.nfsect2) .and. (isect.lt.nfsect2+nsect2)) then
wt(iside,isect,npts(iside,isect))=rkap(2)/sd
else
wt(iside,isect,npts(iside,isect))=rkap(3)/sd
endif
endif
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)*wt(iside,isect,npts(iside,isect))**2
go to 100
c
c
250 do 115 i = 1,6
115 num(i)=0
c
c num(1) is the number of points on plate 2 used to est. rot. of P1 to P2
c num(2) is the number of points on plate 2 used to est. rot. of P1 to P2
c num(3) is the number of points on plate 3 used to est. rot. of P1 to P3
c num(4) is the number of points on plate 1 used to est. rot. of P1 to P3
c num(5) is the number of points on plate 2 used to est. rot. of P2 to P3
c num(6) is the number of points on plate 3 used to est. rot. of P2 to P3
c
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)
c
c minimization section
c
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,*)' '
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,*)' '
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
c order of columns (rows) of sigma3 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
c Construction of the design matrix X
c
ncolx = 2*nsect + 6
if (ndat .gt. ndatmx) write(6,*)'STOP. ndatmx must be increased'
if (ncolx .gt. msig2) write(6,*)'STOP. msect must be increased'
call trimat(ndat,ncolx,ahat,num,nsect1,nfsect1,nsect2,nfsect2,
& nsect3,nfsect3,x)
c
do 350 i = 1,ndat
do 350 j = 1,ncolx
350 xt(j,i)=x(i,j)
call matmul(xt,x,xtx,ncolx,ndat,ncolx)
c
c open(13,file='x.new',status='unknown')
c do 356 i = 1,ndat
c write(13,*)(x(i,j),j=1,ncolx)
c 356 write(13,*)' '
c close(13)
c
c calculation of sigma3 =(xtx)^-1. this is full variance
c covariance matrix (INCLUDING factors of kappahat and rfact) for
c rotation and normal parameters
c
call invert(xtx,ncolx)
c xtx is now replaced by xtx-inverse
do 305 i = 1,ncolx
do 305 j = 1,ncolx
305 sigma3(i,j)=xtx(i,j)
c
c calculation of the projection matrix P-zero
c
call matmul(x,xtx,tempo,ndat,ncolx,ncolx)
call matmul(tempo,xt,pzero,ndat,ndat,ndat)
c
c
c
c calculation of covariance matrix sigma5
c 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)
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
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
c
write(6,*)
write(6,*) 'enter confidence 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
c output to files to be used by conreg to plot confidence regions
c
c calculation of the projection matrix P1 for rotation A1
n1=num(1)+num(2)
n2=num(3)+num(4)
n3=num(5)+num(6)
fr1=n1-3-2*nsect1
fr2=n2-3-2*nsect2
fr3=n3-3-2*nsect3
If (fr1 .lt .fr2 )then
fmin=fr1
fmax=fr2
else
fmin=fr2
fmax = fr1
endif
if (fr3 .lt. fmin) fmin=fr3
if (fr3 .gt. fmax) fmax=fr3
do 360 i = 1,ndat
do 360 j = 1,ncolx-3
360 xk(i,j)=x(i,j+3)
ncolxk=ncolx-3
rank=ncolx
call projection(xk,ndat,ncolxk,p1)
call abc(pzero,p1,a,b,c)
C
call critcalc(a,b,c,plev,cfact,dfr,cxf)
call trans5(qhat(1),axis(1),axis(2),axis(3))
hatkapc = 3./cfact
do 500 i = 1,3
do 500 j = 1,3
cov(i,j) = sigma5(i,j)
500 a(i,j)=ahat(i,j,1)
call inv(cov,covinv)
call resultsout(axis,a,cov,covinv,plev,hatkapc,df,ndat,1,2,
& cxf,fmin,fmax,rank)
c
c calculate matrix P1 and confidence critical value for rotation A2
c
do 400 i = 1,ndat
do 402 j = 1,3
402 xk(i,j) = x(i,j)
do 400 j = 4,ncolxk
400 xk(i,j) = x(i,j+3)
call projection(xk,ndat,ncolxk,p1)
call abc(pzero,p1,a,b,c)
call critcalc(a,b,c,plev,cfact,dfr,cxf)
call trans5(qhat(5),axis(1),axis(2),axis(3))
hatkapc = 3./cfact
do 405 i = 1,3
do 405 j = 1,3
cov(i,j) = sigma5(i+3,j+3)
405 a(i,j) = ahat(i,j,2)
call inv(cov,covinv)
call resultsout(axis,a,cov,covinv,plev,hatkapc,df,ndat,
& 1,3,cxf,fmin,fmax,rank)
c
c calculate confidence critical value for rotation A3
c
do 420 i = 1,n1
do 422 j = 1,3
422 xk(i,j)=x(i,j)
do 420 j = 4,ncolxk
420 xk(i,j)=x(i,j+3)
do 425 i = n1+1,n1+n2
do 425 j = 1,ncolxk
425 xk(i,j)=x(i,j+3)
do 440 i = n1+n2+1,n1+n2+num(5)
do 442 j = 1,3
442 xk(i,j) = x(i,j)
do 440 j = 4,ncolxk
440 xk(i,j) = x(i,j+3)
do 450 i = n1+n2+num(5)+1,ndat
do 450 j = 1,ncolxk
450 xk(i,j) = x(i,j+3)
c
c open(13,file='xk.dat',status='unknown')
c do 454 i = 1,ndat
c write(13,*)(xk(i,j),j=1,ncolxk)
c 454 write(13,*)' '
c close(13)
c
call projection(xk,ndat,ncolxk,p1)
call abc(pzero,p1,a,b,c)
call critcalc(a,b,c,plev,cfact,dfr,cxf)
call trans5(qhat(9),axis(1),axis(2),axis(3))
hatkapc=3./cfact
do 455 i = 1,3
do 455 j = 1,3
cov(i,j) = sigma5(i+6,j+6)
455 a(i,j) = ahat(i,j,3)
call inv(cov,covinv)
call resultsout(axis,a,cov,covinv,plev,hatkapc,df,ndat,
& 2,3,cxf,fmin,fmax,rank)
c
c
stop
end
c
c