c23456789012345678901234567890123456789012345678901234567890123456789012
c
c program constructing a 95% confidence region around a bestfit rotation
c [this program is based on Hellinger1n.f]
c
c Aug. 95 added J-Y.R's computation of volume to output
c
c input data file: file code 11
c line 1 : lat, lon, angle (best rotation)
c line 2 : 3 numbers in free format
c line 3 : rotation matrix
c line 4 : 3 numbers in free format
c line 5 : 3 numbers in free format
c line 6 : 3 numbers in free format
c line 7 : inverse of covariance matrix
c line 8 : 3 numbers in free format
c line 9 : 3 numbers in free format
c line 10 : 3 numbers in free format
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 output files:
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***********************************************************************
c
implicit double precision(a-h,o-z)
real xchi,plev,df,qbing(10),axisr(3)
dimension ahatm(3,4),axis(3),qhat(4),h(3),
& work(6),ahat(3,3),d(3),z(3,3),sigma(3,3),sigma2(3,3)
character filnam*20
data rfact/4.0528473e07/,nsig/4/,maxfn/1000/
1000 format(a)
c
1 write (*,*) ' Enter rotation file name'
read(*,1000) filnam
write (*,*) ' Do you want kappahat to be set to 1?'
write (*,*) ' 1=no 2=yes'
read (*,*) ichoice
open (unit=11,file=filnam,status='old',err=1)
c
c read 2 first comment lines
read(11,1000) filnam
read(11,1000) filnam
c read bestfit rotation: alat,alon and angle
read(11,1000) filnam
read(11,*) alat,along,rho
c=====translate best lat & lon into x,y,z vector
call trans1(alat,along,axis)
c=====translate best lat & lon & angle into a quaternion
call trans3(alat,along,rho,qhat)
c read confidence level
read(11,1000) filnam
read(11,*) plev
c read hatkap and degree of freedom
read(11,1000) filnam
read(11,*) hatkap,df
if (df.eq.0.) df=10000.
c read nbr of points, sections
read(11,1000) filnam
read(11,*) dum
c read rotation matrix
read(11,1000) filnam
do 5 i=1,3
5 read(11,*) (ahat(i,j),j=1,3)
c read covariance matrix
read(11,1000) filnam
do 10 i=1,3
10 read(11,*) dum
c read inverse of covariance matrix: H11.2
read(11,1000) filnam
do 15 i=1,3
15 read(11,*) (sigma(i,j),j=1,3)
do 16 i=1,3
do 16 j=1,3
16 sigma2(i,j)=sigma(i,j)
close(11)
c
if (df.ge.1000.or.ichoice.eq.2) then
call xidch(plev,3.,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidch--ier ',ier
hatkap=1
else
call xidf(plev,3.,df,xchi,ier)
if (ier.ne.0) write(6,*) 'subroutine xidf--ier ',ier
xchi=3.*xchi
endif
xchi=xchi/hatkap
write(6,*) 'Significance level: ',plev
write(*,*) 'Kappahat :',hatkap,' Deg. of freedom: ',df
write(6,*) ' xchi: ',xchi
c
write(6,*)
write(6,*) 'h11.2: ',((sigma(i,j),j=1,i),i=1,3)
c call jacobi(sigma,3,3,d,z,3,work,nrot)
call jacobi(sigma2,3,3,d,z,3,work,nrot)
write(6,*) 'subroutine jacobi--nrot: ',nrot
write(6,*) 'eigenvalues: ',(d(j),j=1,3)
do 12 j=1,3
12 write(6,*) 'eigenvector: ',(z(i,j),i=1,3)
vol=1.
write(6,*)
write(6,*) '"confidence region endmembers"'
do 13 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.))
vol=vol*angled
do 14 i=1,3
h(i)=0.
do 14 k=1,3
14 h(i)=h(i)+ahat(i,k)*z(k,j)
write(6,*) 'v-side eigenvector: ',(sngl(h(i)),i=1,3)
call trans6(h,alat,along)
write(6,*) 'latitude, longitude: ',sngl(alat),sngl(along)
13 write(6,*) 'angle of rotation: radians, degrees ',
& sngl(angle),sngl(angled)
c
c compute volume of confidence region: vol=4/3*pi*a*b*c
c
vol=4*3.14159*vol/3
vol=vol*111.111*111.111*111.111
write(6,*) 'Vol= ',vol,' km3'
c
ahatm(1,1)=-qhat(2)
ahatm(1,2)=qhat(1)
ahatm(1,3)=qhat(4)
ahatm(1,4)=-qhat(3)
ahatm(2,1)=-qhat(3)
ahatm(2,2)=-qhat(4)
ahatm(2,3)=qhat(1)
ahatm(2,4)=qhat(2)
ahatm(3,1)=-qhat(4)
ahatm(3,2)=qhat(3)
ahatm(3,3)=-qhat(2)
ahatm(3,4)=qhat(1)
k=0
do 20 i=1,4
do 20 j=1,i
k=k+1
qbing(k)=0.
do 20 k1=1,3
axisr(k1)=axis(k1)
do 20 k2=1,3
20 qbing(k)=qbing(k)+4.*ahatm(k1,i)*sigma(k1,k2)*ahatm(k2,j)
call bingham(qbing,0.,xchi,axisr,jer)
stop
end