c****************hellinger1.f******************************************
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c simplified September 9, 2000, July 2001
c
c modified by BHK, 6 jan 93, to suppress some output
c 4 March 93 moved output to file 15
c 15 april 93 deleted output
c 3 June 93 added flag to output for use by addplus.f
c Sept 95 declared as real*4: df,plev,xchi,plev1,xchi1
c
c last modified october 31, 1988
C MINOR CHANGE: OCTOBER 14, 1991 -- CHANGE IN CAPS
c
c program implementing section 1 of 'on reconstructing tectonic plate
c motion from ship track crossings'
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--moving or 2--fixed)
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
c this program estimates the rotation to rotate side 1 into side 2.
c
c starting with an initial guess of the optimal rotation, the program
c does a grid search of all rotations within .2 radians of the initial
c guess to refine the guess.
c it then 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***********************************************************************
c
program Hellinger1
implicit double precision(a-h,o-z)
parameter (msect=50,msig=2*msect*msect+7*msect+6,msig2=2*msect+3,
& mwork=2*msig2)
real*4 xchi,plev,df,plev1,xchi1
dimension sigma(2,msect,3,3),axis(3),qhat(4),h(3),cov(3,3),
& hp(3,4),yp(4),work(mwork),qhati(4),eta(msect,3,3),ahat(3,3),
& temp(3,3),d(msig2),z(msig2,msig2),etai(msect,3,2),sigma2(msig),
& sigma3(msig2,msig2),sigma4(3,3)
character filnam*20,filnam2*20,ans*1
external r1
common nsect,sigma,qhati,eta,etai
data rfact2/6366.197724/,nsig/4/,maxfn/1000/
DATA HATKAP/0./,DF/0./,FKAP2/0./,FKAP1/0./
1000 format(a)
c
do 50 i1=1,2
do 50 i2=1,msect
do 50 i3=1,3
do 50 i4=1,3
50 sigma(i1,i2,i3,i4)=0.
c
c calculating matrices sigma and sigma-tilde
c
55 write(6,*) 'Data file name? '
read(5,1000) filnam
open(10,file=filnam,status='old',err=55)
write(6,*) 'Output file name? '
read(5,1000) filnam2
open(15,file=filnam2,status='unknown')
read(10,*) nsect
ndat=0
100 read(10,*,end=250,err=100) iside,isect,alat,along,sd
if (iside.gt.2) goto 100
sd = sd/rfact2
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)
go to 100
c
c minimization section--grid search
c
250 write(6,*) 'Initial guess: alat, along, rho? '
read(5,*) alat,along,rho
write(6,*) 'Search radius (radians)? '
read(5,*) eps
call trans3(alat,along,rho,qhati)
write(6,*) 'qhat: ',qhati
do 207 i=1,3
207 h(i)=0.
rmin=r1(h)
251 write(6,*) 'do you want a grid search, eps = ',eps
read(5,'(a)') ans
if ((ans.ne.'Y').and.(ans.ne.'y')) go to 260
write(6,*) 'calling grid search--r = ',rmin
call grds(qhati,eps,rmin)
if (rmin.gt.1.e30) then
write(6,*) 'error return from grid--try another initial guess'
go to 250
endif
write(6,*) 'return from grid search--r = ',rmin
write(6,*) 'qhat: ',qhati
go to 251
c
c minimization section--calling amoeba
c
260 do 215 j=1,4
do 210 i=1,3
210 hp(i,j)=0.
if (j.ne.1) hp(j-1,j)=eps
215 yp(j)=r1(hp(1,j))
rmin=yp(1)
ftol=(.1)**(2*nsig+1)
iter=maxfn
write(6,*)
write(6,*) 'calling amoeba--r = ',rmin
rprev=rmin
call amoeba(hp,yp,3,3,ftol,r1,iter,work)
rmin=r1(hp(1,1))
call trans2(hp(1,1),qhati,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
call trans5(qhat,alat,along,rho)
write(6,*) 'Fitted rotation: ',alat,along,rho
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,4
212 qhati(i)=qhat(i)
go to 260
endif
call trans1(alat,along,axis)
call trans4(qhat,ahat)
write(6,*) 'Fitted rotation--alat,along,rho: '
write(6,*) alat,along,rho
write(6,*) 'ahat: '
do 211 i=1,3
211 write(6,*) (ahat(i,j),j=1,3)
c
c
c replace sigma-tilde with (ahat**t)*sigma-tilde*ahat
c
219 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)*sigma(2,isect,k1,k2)*ahat(k2,j)
do 230 i=1,3
do 230 j=1,3
230 sigma(2,isect,i,j)=temp(i,j)
c
c calculation of (x**t)*x matrix--stored in sigma2 in symmetric mode
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
do 300 i=1,msig
300 sigma2(i)=0.
k=0
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)
do 302 isect=1,nsect
do 302 i=1,2
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(2,isect,k1,k2)*eta(isect,j,k2)
k=k+2*(isect-1)
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))*etai(isect,k2,j)
c
ndim=2*nsect+3
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 sigma3 is the full variance-covariance matrix of the rotation and
c normal section parameters (except possibly for division by hatkap).
c rotation parameters are stored in upper 3 by 3 corner.
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)
write(6,*)
write(6,*) 'sigma3: '
write(6,*) ((sigma3(i,j),j=1,i),i=1,3)
c
c calculation of critical point
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,3.,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidch--ier ',ier
else
df=ndat-2*nsect-3
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,3.,df,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidf--ier: ',ier
xchi=xchi*rmin*3./df
endif
c
c calculation of h11.2 matrix--stored sigma4
c
do 317 i=1,3
do 317 j=1,3
317 cov(i,j)=sigma3(i,j)
call jacobi(sigma3,3,msig2,d,z,msig2,work,nrot)
if (nrot.lt.0) write(6,*) 'subroutine jacobi(2)--nrot ',nrot
do 311 i=1,3
do 311 j=1,3
sigma4(i,j)=0.
do 311 k=1,3
311 sigma4(i,j)=sigma4(i,j)+z(i,k)*z(j,k)/d(k)
do 316 i=1,3
316 d(i)=1./d(i)
write(6,*)
write(6,*) 'h11.2: ',((sigma4(i,j),j=1,i),i=1,3)
write(6,*) 'eigenvalues: ',(d(j),j=1,3)
do 312 j=1,3
312 write(6,*) 'eigenvector: ',(z(i,j),i=1,3)
write(6,*)
write(6,*) '"confidence region endmembers"'
do 313 j=1,3
write(6,*) 'u-side eigenvector: ',(sngl(z(i,j)),i=1,3)
call trans6(z(1,j),alat,along)
write(6,*) 'latitude, longitude: ',sngl(alat),sngl(along)
angle=(d(j)/xchi)**(-.5)
angled=angle*45./atan(dble(1.))
j1=2*j
j2=j1+1
do 314 i=1,3
h(i)=0.
do 314 k=1,3
h(i)=h(i)+ahat(i,k)*z(k,j)
314 continue
write(6,*) 'v-side eigenvector: ',(sngl(h(i)),i=1,3)
call trans6(h,alat,along)
write(6,*) 'latitude, longitude: ',sngl(alat),sngl(along)
313 write(6,*) 'angle of rotation: radians, degrees ',
& sngl(angle),sngl(angled)
c
write(15,1460) filnam
1460 format ('Results from Hellinger1 using ',a,/)
write(15,*) 'Fitted rotation--alat,along,rho: '
call trans5(qhat,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 used by addplus.f to indicate origin of estimates
write(15,*) 'kappahat, degrees of freedom,xchi,flag'
write(15,*) hatkap,df,xchi,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),j=1,3)
write(15,*) 'covariance matrix'
do 465 i=1,3
465 write(15,*) (cov(i,j),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
c
stop
end
c
c
function r1(h)
implicit double precision(a-h,o-z)
parameter (msect=50)
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
subroutine grds(qhati,eps,rmin)
implicit double precision (a-h,o-z)
dimension h(3),qhati(4),qhat(4)
rmin=1.1e30
do 100 i=-5,5
do 100 j=-5,5
do 100 k=-5,5
h(1)=i*eps/5.
h(2)=j*eps/5.
h(3)=k*eps/5.
r=r1(h)
if (r.lt.rmin) then
imin=i
jmin=j
kmin=k
rmin=r
endif
100 continue
h(1)=imin*eps/5.
h(2)=jmin*eps/5.
h(3)=kmin*eps/5.
call trans2(h,qhati,qhat)
do 20 i=1,4
20 qhati(i)=qhat(i)
eps=eps/5.
return
end