c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c
c program addplus June 1, 1993
c
c Sept. 1995 declared plev,df,xf as real*4
c
c
c This program composes two independently estimated rotations and
c calculates the covariance matrix of the composition.
c The rotations are assumed to have unequal kappas,
c or to be composites of rotations with unequal kappas. A conservative
c estimate of the critical value for construction of a 95% confidence region
c is produced. The output file is formatted to be used by conreg.f to
c produce plots of the confidence regions, or to be used as input to this
c program for further compositions.
c Output from invrot, hellinger1 or hellinger3 can be used provided
c that a flag value of 0. is included at the end of the line that gives values
c for hatkap, degrees of freedom and xchi.
c
c COMPILE WITH DTRANS AND NUMLIB
c
program addplus
implicit double precision (a-h,o-z)
character filnam*70,name(2)*15
dimension a(3,3),cova(3,3),b(3,3),covb(3,3),c(3,3),covc(3,3),
&w(3,3),xa(3),xb(3),xc(3),h(3)
REAL*4 PLEV,XF,XCHI
DATA NDAT/0/,NSECT/0/
c
itour=1
1 write (*,1000) itour
1000 format ('Input file ',i1,': ',$)
read (*,'(a)') name(itour)
open (unit=10,file=name(itour),status='old',err=1)
write (*,' (a)') name(itour)
c read two first comment lines
read (10,'(a)') filnam
write (*,' (a)') filnam
read (10,'(a)') filnam
write (*,' (a)') filnam
c read alat,alon,angle
read(10,'(a)') filnam
write (*,' (a)') filnam
read(10,*) xb
write (*,'(3x,3f8.2)') xb
c read conf. level
read(10,'(a)') filnam
write (*,*) ' confidence level'
read (10,*) plev
write(*,'(f6.2)') plev
c read kappa estimated, degrees of freedom
c flagb = 1 if the rotation is a composition of rotations with unequal
c kappas.
read(10,'(a)') filnam
write (*,' (a)') filnam
read(10,*) hatkapb,dfb,xchi,flagb
write(*,'(3x,f5.2,f5.0,f8.3,f3.0)') hatkapb,dfb,xchi,flagb
c read number of points
read(10,'(a)') filnam
if (flagb .ne. 1.) then
read(10,*) nptsb
fmaxb=dfb
fminb=dfb
rankb=nptsb-dfb
else
read (10,*) nptsb, fminb,fmaxb,rankb
endif
write(6,*)'npts, fmin, fmax, rank'
write(6,7) nptsb,fminb,fmaxb,rankb
7 format(1x,i5,3x,f6.1,3x,f6.1,3x,f6.1)
c read rotation matrix
read(10,'(a)') filnam
write (*,*)' ahat:'
do 3 i=1,3
read(10,*) (b(i,j),j=1,3)
3 write (*,*) (b(i,j),j=1,3)
c read covariance matrix
read(10,'(a)') filnam
write (*,'(a)') filnam
do 5 i=1,3
read (10,*) (covb(i,j),j=1,3)
5 write (*,*) (covb(i,j),j=1,3)
write(*,*) ' '
close(10)
c
itour=itour+1
goto (20,35) itour-1
20 do 30 i=1,3
xa(i)=xb(i)
do 30 j=1,3
a(i,j)=b(i,j)
30 cova(i,j)=covb(i,j)
hatkapa=hatkapb
dfa=dfb
nptsa=nptsb
fmina=fminb
fmaxa=fmaxb
ranka=rankb
flaga = flagb
goto 1
c
35 if (fmina .lt. fminb) then
fmin=fmina
else
fmin=fminb
endif
if (fmaxa .gt. fmaxb) then
fmax = fmaxa
else
fmax = fmaxb
endif
n = nptsa + nptsb
r= ranka + rankb
c
z = (r-3.)**2/n + 6 - r
if (z .lt. 0.) then
al=z/fmin
else
al = z/fmax
endif
if (al .gt. -1.5 ) then
bigk = 6*(1 + 9/fmin)/(3+2*al)**2
df = (12*bigk + 2)/(3*bigk - 2)
else
df = 4.0
endif
cfact = (3 + 6/fmin)*(1-2/df)
call xidf(plev,3.0,df,xf,ier)
cxf = cfact*xf
if (flagb .ne. 1.) then
DO 36 I=1,3
DO 36 J=1,3
36 COVB(I,J)=COVB(I,J)/HATKAPB
endif
if (flaga .ne. 1.) then
do 37 i = 1,3
do 37 j = 1,3
37 cova(i,j)=cova(i,j)/hatkapa
endif
C
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
call mul(b,a,c)
call trans9(c,h,rho,errchk)
call trans6(h,xc(1),xc(2))
xc(3)=rho*57.29577951
c
50 write (*,*) ' Product rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') xc
write(*,*) 'adjusted kappa and degrees of freedom'
write(*,'(3x,f9.4,2x,f9.4)') hatkapc,df
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,*) 'conf. level '
write(15,*) plev,1.0,1.0
write(15,*) '3./cfact, degrees of freedom, crit.value, flag'
write(15,'(f12.6,2x,f9.4,2x,f9.4,2x,f3.1)') hatkapc,df,cxf,1.0
write(15,*) 'Number of points, fmin,fmax,rank'
write(15,'(7x,i4,3f5.0)') n,fmin,fmax,r
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 inv(a,b)
c invert 3x3 matrix
double precision a(3,3), b(3,3),det
det=a(1,1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3))-
& a(1,2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3))+
& a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
if (det.eq.0.) goto 5
b(1,1)=a(2,2)*a(3,3)-a(3,2)*a(2,3)
b(1,2)=a(3,2)*a(1,3)-a(1,2)*a(3,3)
b(1,3)=a(1,2)*a(2,3)-a(2,2)*a(1,3)
b(2,1)=a(3,1)*a(2,3)-a(2,1)*a(3,3)
b(2,2)=a(1,1)*a(3,3)-a(3,1)*a(1,3)
b(2,3)=a(2,1)*a(1,3)-a(1,1)*a(2,3)
b(3,1)=a(2,1)*a(3,2)-a(3,1)*a(2,2)
b(3,2)=a(3,1)*a(1,2)-a(1,1)*a(3,2)
b(3,3)=a(1,1)*a(2,2)-a(2,1)*a(1,2)
do 1 i=1,3
do 1 j=1,3
1 b(i,j)=b(i,j)/det
return
5 write (*,*) 'Trouble: det=0'
stop
end
c****************************************************************
c
subroutine mul(a,b,c)
c multiply 3x3 matrices: c=ab
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****************************************************************