c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c program add2.f April 4, 1993
c
c modified Sept. 1995 declaring plev,df,xf as real*4
c modified June 1994 improving input from data files
c
c this version uses LAPACK for matrix inversion
c
c can be shortened by constructing XK directly, not using X
c
c This is an adaptation of addrot.f to be used when kappas are unequal.
c Critical value for the confidence region construction is calculated
c based on estimates developed by S. Johansen
c
c
C COMPILE WITH NUMLIB.F , numlib2.f and dtrans.f
C
c add two rotations and calculate the covariance matrix of the
c product rotation
c chat=bhat.ahat
c =b*phi(tbhat)*a*phi(tahat)
c =b*a*(a**t)*phi(tbhat)*a*phi(tahat) ..
c =b*a*phi((a**t)*tbhat)*phi(tahat)
c ~b*a*phi((a**t)*tbhat +tahat) up to linear approximations
c
c tchat~(a**t)*tbhat +tahat
c cov(tchat)~(a**t)*cov(tbhat)*a + cov(tahat)
C
c
c
c
program add2
implicit double precision (a-h,o-z)
parameter (ndatmx=500,ndtmx2=2*ndatmx)
parameter (msect = 100,msig=2*msect*msect*msect+7*msect+6,msig2=
& 2*msect+3,mwork=2*msig2)
character filnam*70,name(2)*15,dum*1,datdum*30,datnam(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)
dimension axis(3),h(3),etai(msect,3,2),qhati(4),eta(msect,3,3),
& num(2),wt(2,msect,20)
dimension u(2,msect,20,3),x2(ndatmx,ndatmx),
& x(ndatmx,ndatmx),xt(ndatmx,ndatmx),npts(2,msect),axm(3,3),
& tempo(ndatmx,ndatmx),xtsx(ndatmx,ndatmx)
dimension xtsxinv(ndatmx,ndatmx),pone(ndatmx,ndatmx),
& xk(ndatmx,ndatmx),p11(ndatmx,ndatmx),p22(ndatmx,ndatmx)
dimension p11sq(ndatmx,ndatmx),p22sq(ndatmx,ndatmx),
& sigma(2,msect,3,3)
REAL*4 PLEV,df,XF
external R1
Common /vectors/nsect,sigma,qhati,eta, etai
DATA NDAT/0/,rfact/4.0528473E07/,rfact2/6366.1977/,
& nsig/4/,maxfn/1000/,fmin/5000./,fmax/0./
1010 format(A)
1012 format(A29,A15)
1013 format(1x,A29,A15)
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,1012) datdum,datnam(itour)
write (*,1013) datdum,datnam(itour)
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 (*,'(a)') filnam
read (10,*) plev,d1,d2
write(*,'(3f6.2)') plev,d1,d2
c read kappa estimated, degrees of freedom
read(10,'(a)') filnam
write (*,*) 'kappahat, degrees of freedom'
read(10,*) hatkapb,dfb
write(*,'(3x,f5.2,f5.0)') hatkapb,dfb
if (dfb .lt. fmin) fmin = dfb
if (dfb .gt. fmax) fmax = dfb
c read number of points, sections ...
read(10,'(a)') dum
read(10,*) nrowx1,LS
c read rotation matrix
read(10,'(a)') filnam
write (*,'(a)') filnam
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
c
17 write(*,1001)
1001 format(' enter name of data file for this rotation ' ,$)
read(*,'(a)',err=17)filnam
OPEN(10,FILE=FILNAM,STATUS='UNKNOWN')
REWIND 10
READ(10,*) nsect
NDAT=0
num(iside)=0
NUM(1)=0
NUM(2)=0
do 51 I1 = 1,2
do 51 I2 = 1,nsect
npts(i1,i2) = 0
do 51 i3 = 1,3
do 51 i4 = 1,3
51 sigma(i1,i2,i3,i4) = 0.
do 200 ktimes = 1,nrowx1
100 READ(10,*, err=100) ISIDE,ISECT,ALAT,ALONG,SD
if (iside.gt.2) goto 100
npts(iside,isect) = npts(iside,isect) + 1
CALL TRANS1(ALAT,ALONG,AXIS)
n = npts(iside,isect)
wt(iside,isect,n) = 1.0/sd
do 105 i = 1,3
105 u(iside,isect,n ,i) = axis(i)
NUM(ISIDE)=Num(iside) + 1
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)
200 continue
201 close (10)
c
c setting up the eta vectors
c
250 alat = xb(1)
along = xb(2)
rho = xb(3)
CALL TRANS3(ALAT,ALONG,RHO,QHATI)
DO 207 I=1,3
207 H(I)=0.
RMIN=R1(H)*RFACT
c
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
c
c construct estimate of separate fit design matrix X, beginning
c with first 3 columns.
c
nrowx = nrowx1
ncolx = 3 + 2*nsect
do 505 k = 1, nrowx
do 505 j = 1, ncolx
505 x(k,j) = 0.
nrow = 1
do 610 isect = 1, nsect
do 570 k = 1,3
do 570 j = 1,3
570 axm(k,j) = 0.
do 580 k = 1,3
do 580 j = 1,3
do 580 m = 1,3
580 axm(k,j) = axm(k,j) + b(k,m)*eta(isect,j,m)
do 600 np = 1,npts(2,isect)
do 595 k = 1,3
do 590 j = 1,3
590 x(nrow,k) = x(nrow,k) + u(2,isect,np,j)*axm(j,k)
595 x(nrow,k) = x(nrow,k)*wt(2,isect,np)*rfact2
600 nrow = nrow + 1
610 continue
nr = 1
nc = 4
do 650 isect = 1,nsect
do 620 k = 1,3
do 620 j = 1,2
620 axm(k,j) = 0.
do 630 k = 1,3
do 630 j = 1,2
do 630 m = 1,3
630 axm(k,j) = axm(k,j) + b(k,m)*etai(isect,m,j)
do 645 np = 1,npts(2,isect)
do 640 k = 1,3
x(nr,nc) = x(nr,nc) + u(2,isect,np,k)*axm(k,1)
640 x(nr,nc+1) = x(nr,nc+1) + u(2,isect,np,k)*axm(k,2)
x(nr,nc) = wt(2,isect,np)*x(nr,nc)*rfact2
x(nr,nc+1) = wt(2,isect,np)*x(nr,nc+1)*rfact2
645 nr = nr + 1
650 nc = nc + 2
nc = 4
do 680 isect = 1,nsect
do 670 np = 1,npts(1,isect)
do 660 k = 1,3
x(nr,nc) = x(nr,nc) + u(1,isect,np,k)*etai(isect,k,1)
660 x(nr,nc+1) = x(nr,nc+1) +
& u(1,isect,np,k)*etai(isect,k,2)
x(nr,nc) = x(nr,nc)*wt(1,isect,np)*rfact2
x(nr,nc+1) = x(nr,nc+1)*wt(1,isect,np)*rfact2
670 nr = nr + 1
680 nc = nc + 2
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
ks = ls
nrowx2 = nrowx
ncolx2=ncolx
do 31 i = 1,nrowx2
do 31 j = 1,ncolx2
31 x2(i,j) = x(i,j)
goto 1
c
c Construct estimated design matrix X from the separate fit matrices
c
35 nrowx1=nrowx
ncolx1=ncolx
nrowx = nrowx1 + nrowx2
if (nrowx .gt. ndatmx) write(*,*)'STOP. number of data points',
& 'exceeds ndatmx. Parameter statement must be changed.'
ncolx = ncolx1 + ncolx2
if(ncolx.gt.3+2*msect) write(*,*) 'STOP! msect exceeded!'
do 32 i = 1,nrowx2
do 32 j = 1,ncolx2
32 x(i+nrowx1,j+ncolx1)=x2(i,j)
do 33 i = 1,nrowx1
do 33 j = ncolx1+1,ncolx
33 x(i,j) = 0.0
do 34 i = nrowx1+1,nrowx
do 34 j = 1,ncolx1
34 x(i,j) = 0.0
c
c transform X to S^(-1/2)*X
c
do 303 i = 1,nrowx1
do 303 j = 1,ncolx
303 x(i,j) = hatkapb**.5*x(i,j)
do 305 i = nrowx1 + 1,nrowx
do 305 j = 1,ncolx
305 x(i,j) = hatkapa**.5*x(i,j)
c
c construct XK
c
ncolxk = ncolx - 3
do 320 i = 1,nrowx1
do 320 j = 1,ncolxk
320 xk(i,j) = x(i,j)
do 360 i = nrowx1+1,nrowx
do 330 j = 1,3
330 xk(i,j)=-x(i,ncolx1+j)
do 340 j = 4,ncolx1
340 xk(i,j) = 0.0
do 350 j = ncolx1 + 1,ncolxk
350 xk(i,j)=x(i,j+3)
360 continue
c
c multiply columns 1-3, rows nrowx1+1 and down, by A-transpose
c
do 368 i = nrowx1 + 1, nrowx
do 364 j = 1,3
tempo(i,j) = 0.0
do 364 k = 1,3
364 tempo(i,j) = tempo(i,j) + xk(i,k)*a(j,k)
do 366 j = 1,3
366 xk(i,j) = tempo(i,j)
368 continue
c
c open(13,file='x.out',status='unknown')
c do 370 i = 1,nrowx
c write(13,*) (x(i,j),j = 1,ncolx)
c 370 write(13,*) ' '
c close (13)
c open(13,file='xk.out',status='unknown')
c do 380 i = 1,nrowx
c write(13,*)(xk(i,j),j = 1,ncolxk)
c 380 write (13,*)' '
c close(13)
c
c
c calculate P-one
c
do 400 i = 1,nrowx
do 400 j = 1,ncolxk
400 xt(j,i) = xk(i,j)
call matmul(xt,xk,xtsx,ncolxk,nrowx,ncolxk)
do 401 i = 1,ncolxk
do 401 j = 1,ncolxk
401 xtsxinv(i,j)=xtsx(i,j)
c
c open(13,file='xtsx',status='unknown')
c do 405 i = 1,ncolxk
c 405 write(13,*)(xtsx(i,j),j=1,ncolxk)
c close(13)
write(6,*) 'inverting matrices'
c
call dpotrf('U',ncolxk,xtsxinv,ndatmx,info)
if (info.ne.0) write(6,*)'failure in LAPACK routine. info = ',info
call dpotri('U',ncolxk,xtsxinv,ndatmx,info)
if (info.ne.0) write(6,*)'failure in LAPACK routine. info = ',info
do 407 i = 2,ncolxk
do 407 j = 1,i-1
407 xtsxinv(i,j)=xtsxinv(j,i)
c
c
call matmul(xtsxinv, xtsx, tempo,ncolxk,ncolxk,ncolxk)
c open(13,file='Lapackproduct',status='unknown')
c do 408 i = 1,ncolxk
c write(13,*)(tempo(i,k),k=1,ncolxk)
c 408 write(13,*)' '
c close(13)
c
call matmul(xtsxinv,xt,tempo,ncolxk,ncolxk,nrowx)
call matmul(xk,tempo,pone,nrowx,ncolxk,nrowx)
c
c open(13,file='pone.out',status='unknown')
c do 780 i = 1,nrowx
c write(13,*)(pone(i,j),j = 1,nrowx)
c 780 write(13,*)' '
c close(13)
c
c Calculation of correction terms A and C
c
do 800 i = 1,nrowx1
do 800 j = 1,nrowx1
800 p11(i,j)=pone(i,j)
do 805 i = 1,nrowx2
do 805 j = 1,nrowx2
805 p22(i,j)=pone(i+nrowx1,j+nrowx1)
call matmul(p11,p11,p11sq,nrowx1,nrowx1,nrowx1)
call matmul(p22,p22,p22sq,nrowx2,nrowx2,nrowx2)
call trace(p11,nrowx1,trp11)
call trace(p11sq,nrowx1,trp11sq)
call trace(p22,nrowx2,trp22)
call trace(p22sq,nrowx2,trp22sq)
c dfb=f1,dfa=f2
capa=(ncolx1+trp11sq-2*trp11)/dfb + (ncolx2+trp22sq-2*trp22)/dfa
capc = (ncolx1-trp11)**2/dfb + (ncolx2-trp22)**2/dfa - capa
write(6,*)'A, C',capa,capc
c
c calculate c and f, where Q is cF(3,f)
c
exq = 3 + 2*capa
vq = 6 + 14*capa + 2*capc
write (6,*) 'mean of Q = ', exq
write (6,*) 'variance of Q = ', vq
ex2 = exq*exq
write (6,*)
df = (12*vq + 2*ex2)/(3*vq - 2*ex2)
cfact = exq*(1 - 2/df)
write (6,*) 'deg. freedom ', df
write (6,*) ' c = ',cfact
write(*,*)'enter confidence level'
read(*,*) plev
call xidf(plev,3.0,df,xf,ier)
cxf = cfact*xf
write(6,*)'xf, critical value= c*f(3,df) :',xf,cxf
write(6,*)' '
write(6,*)' '
c
DO 36 I=1,3
DO 36 J=1,3
COVA(I,J)=COVA(I,J)/HATKAPA
36 COVB(I,J)=COVB(I,J)/HATKAPB
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 sumrot(xa,xb,xc)
call mul(b,a,c)
c
50 write (*,*) ' Product rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') xc
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,*)'plev'
write(15,*) plev,1.0,1.0
c flag = 1 is output for use by the prgram addplus.f
write(15,*) '3/cfact, degrees of freedom, crit. value, flag'
write(15,'(f12.6,2x,f8.4,2x,f12.4,2x,f3.1)') hatkapc,df,cxf,1.0
write(15,*) 'Number of points, fmin, fmax, rank'
write(15,*) nrowx,fmin,fmax,ncolx
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=b.a
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****************************************************************
c
subroutine sumrot (a,b,c)
c
c finds the total rotation given the two successive rotations
c a(1)=latitude, a(2)=longitude, a(3)=angle
c
implicit double precision (a-h,o-z)
dimension a(3),b(3),c(3)
data rad/.0174532925199432958/
call cdtrn (a(1),a(2),a(3),w1,x1,y1,z1)
call cdtrn (b(1),b(2),b(3),w2,x2,y2,z2)
wt=w1*w2-x1*x2-y1*y2-z1*z2
xt=w1*x2+x1*w2-y1*z2+z1*y2
yt=w1*y2+x1*z2+y1*w2-z1*x2
zt=w1*z2-x1*y2+y1*x2+z1*w2
if(wt) 5,1,1
1 tt=acos(wt)
alt=acos(zt/( sin(tt)))
go to 10
5 tt=acos(-wt)
alt=acos(zt/(-sin(tt)))
10 pt=atan2(yt,xt)
c(1)=90.-alt/rad
gs=pt/rad
if(gs.gt.180.) gs=gs-360.
c(2)=gs
c(3)=tt*2./rad
return
end
c*****************************************************************
c
subroutine cdtrn (f,g,o,w,x,y,z)
implicit double precision (a-h,o-z)
data rad/.0174532925199432958/
pi=rad*180.
al=(90.-f)*rad
p=g*rad
t=o*rad
st=sin(t/2.)
w=cos(t/2.)
x=st*sin(al)*cos(p)
y=st*sin(al)*sin(p)
z=st*cos(al)
if(t-pi) 5,5,1
1 w=-w
x=-x
y=-y
z=-z
5 return
end
c
c
C
C
FUNCTION R1(H)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER (MSECT=100)
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 /vectors/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
c
c
subroutine matmul(a,b,c,nra,nrb,ncb)
c
c this subroutine only works for matrices dimensioned ndatmx by ndatmx
c
implicit double precision (a-h,o-z)
parameter (ndatmx = 500)
dimension a(ndatmx, ndatmx),b(ndatmx,ndatmx),c(ndatmx,ndatmx)
do 10 i = 1,nra
do 10 j = 1, ncb
10 c(i,j) = 0.
do 21 i = 1, nra
do 21 j = 1,ncb
do 20 k = 1, nrb
20 c(i,j) = c(i,j) + a(i,k)*b(k,j)
21 continue
return
end
c
c
c
subroutine trace(a,ndima,t)
c
implicit double precision (a-h,o-z)
parameter (ndatmx = 500)
dimension a(ndatmx,ndatmx)
t = 0.
do 10 k = 1,ndima
10 t = t + a(k,k)
return
end