c23456789 averotkap
c
c estimates a single rotation by combining several different independent
c estimates from data sets with different kappas. Reads fit files and
c data sets.
c
c Compile with dtrans,numlib,numlib2,teclib
c CAUTION: changes in the parameters ndatmx and msect MUST also be made
c in all routines in teclib
c
program averotkap
c
implicit double precision (a-h,o-z)
parameter (ndatmx=400,msect=70)
character*40 filename
real*4 plev,df,xf
dimension work(6),evd(3),t(3),ttot(3),uc(3),
&qhat(10,3),arot(3),ahat(4),bhat(4),chat(4),
&a(3,3),b(3,3),h(3,3),evz(3,3),v(3,3),p(10,3,3)
dimension xk(ndatmx,ndatmx),numsect(10),np(10),hkap(10),
& qhati(4),wt(2,msect,30),sigma(2,msect,3,3),
& eta(msect,3,3),axis(3),h2(3),numpts(2,msect),axm(3,3)
dimension u(2,msect,ndatmx,3),pone(ndatmx,ndatmx),
& pii(10,ndatmx,ndatmx),trpii(10),trpiisq(10),dfr(10),
& etai(msect,3,2),w(3,3)
common nsect,sigma,qhati,eta,etai,wt,npts,u
data rad/.0174532925199432958/,rfact/4.0528473E07/,
& rfact2/6366.1977/, fmin/5000./,fmax/0./
c
if=0
idet=0
nrowdone=0
ncoldone=3
nsectdone=0
do 17 i =1,ndatmx
do 17 j = 1,ndatmx
17 xk(i,j)=0.0
detmax=-1.d+10
1 write (*,1000)
1000 format ('Enter rotation file name: (=end to end) ',$)
read (*,'(a)') filename
if (filename.eq.'end') goto 10
call readrot(filename,ior,arot(1),arot(2),arot(3),a,b,h,plev,
& hatkap,df,ndum1)
if (ior.ne.0) goto 1
if=if+1
write (*,'(3f7.2)') (arot(i),i=1,3)
do 5 i=1,3
qhat(if,i)=arot(i)
do 5 j=1,3
5 p(if,i,j)=h(i,j)*hatkap
call detr(h,d)
d=d*hatkap
write (*,'(2i3,d12.3)') if,idet,d
c
c read data
c
1009 write(*,*)' '
write(*,1010)
1010 format(' Enter data filename ',$)
read(*,'(a)')filename
write(*,*)' '
open(10,file=filename,status='old',err=1009)
read(10,*)numsect(if)
do 20 i1=1,2
do 20 i2=nsectdone+1,nsectdone+numsect(if)
numpts(i1,i2)=0
do 20 i3=1,3
do 20 i4=1,3
20 sigma(i1,i2,i3,i4)=0.0
hkap(if)=hatkap
np(if)=ndum1
dfr(if)=np(if)-2*numsect(if)-3
if (fmin .gt. dfr(if)) fmin=dfr(if)
if (fmax .lt. dfr(if)) fmax=dfr(if)
do 100 nk = 1,np(if)
101 read(10,*,err=101)iside,isect,alat,along,sd
if (iside.gt.2) goto 101
isect=nsectdone+isect
numpts(iside,isect)=numpts(iside,isect)+1
call trans1(alat,along,axis)
wt(iside,isect,numpts(iside,isect))=1/sd
do 105 i = 1,3
105 u(iside,isect,numpts(iside,isect),i)=axis(i)
do 120 j1=1,3
do 120 j2=1,3
120 sigma(iside,isect,j1,j2)=sigma(iside,isect,j1,j2)+
& axis(j1)*axis(j2)/sd**2
100 continue
nsectdone=nsectdone+numsect(if)
close(10)
if (detmax.gt.d) goto 1
idet=if
detmax=d
goto 1
c------- calculate covariance matrix v
c
10 nsect=nsectdone
nrowx=0
do 13 i = 1,if
13 nrowx=nrowx+np(i)
ncolxk=3+2*nsect
if (nsect.gt.msect) write(*,*)'STOP! msect exceeded'
if (nrowx.gt.ndatmx)write(*,*)'STOP! ndatmx exceeded'
if ((nsect.gt.msect).or.(nrowx.gt.ndatmx)) stop
do 15 k=1,if
do 15 i=1,3
do 15 j=1,3
15 v(i,j)=v(i,j)+p(k,i,j)
call detr(v,d)
write (*,*) 'Det(v) ',d
c
call jacobi(v,3,3,evd,evz,3,work,nrot)
do 220 i=1,3
do 220 j=1,3
v(i,j)=0.
do 220 k = 1,3
220 v(i,j)=v(i,j)+evz(i,k)*evz(j,k)/evd(k)
c
c------- start main loop
c
21 do 22 i=1,3
ttot(i)=0.
22 arot(i)=qhat(idet,i)
c
25 continue
c
call trans3(arot(1),arot(2),-arot(3),ahat)
c
do 35 k=1,if
if (k.eq.idet) goto 35
c
c -- convert lat,lon,angle into quaternion bhat
c
call trans3(qhat(k,1),qhat(k,2),qhat(k,3),bhat)
c
c -- add the 2 quaternions
c
call qmult(ahat,bhat,chat)
c
c -- convert quaternion chat into lat,lon,angle
c
call trans5(chat,alat,alon,rho)
c write (*,'(3f9.3)') alat,alon,rho
call trans1(alat,alon,t)
do 27 i=1,3
27 t(i)=t(i)*rho*rad
c write(*,'(3E16.8)') t
do 29 i=1,3
do 29 j=1,3
29 h(i,j)=p(k,i,j)
call rotvec(h,t,uc)
do 30 i=1,3
30 ttot(i)=ttot(i)+uc(i)
c write(*,'(3E16.8)') ttot
c
35 continue
c
c------------ end of main loop
c
call rotvec(v,ttot,t)
write(*,*) 't: ',t
c t converges to 0 as a-hat converges to A
rho=0.
do 40 i=1,3
40 rho=rho+t(i)*t(i)
rho=sqrt(rho)
do 42 i=1,3
42 t(i)=t(i)/rho
call trans6(t,alat,alon)
call trans3(alat,alon,rho/rad,bhat)
do 45 i=2,4
45 ahat(i)=-ahat(i)
call qmult(ahat,bhat,chat)
call trans5(chat,alat,alon,rho)
write (*,*) ' '
write (*,'(3f9.3)') alat,alon,rho
write (*,*) 'Happy with this? 1=yes 2=no'
read (*,*) ian
if (ian.eq.1) goto 50
c
idet=if+1
qhat(idet,1)=alat
qhat(idet,2)=alon
qhat(idet,3)=rho
c
goto 21
c
c set up eta vectors
c
c chat is the quaternion form of new estimate
c
50 call trans4(chat,a)
c
write(*,*)'rotation matrix'
do 51 i = 1,3
51 write(*,*)(a(i,j),j=1,3)
c
do 55 i = 1,4
55 qhati(i)=chat(i)
do 130 i = 1,3
130 h2(i)=0.
rmin=r1(h2)*rfact
c
C
nrowdone=0
nsectdone=0
ncoldone=3
do 615 nrot=1,if
nrow = 1+nrowdone
do 610 isect = 1+nsectdone, numsect(nrot)+nsectdone
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) + a(k,m)*eta(isect,j,m)
do 600 kp = 1,numpts(2,isect)
do 595 k = 1,3
do 590 j = 1,3
590 xk(nrow,k) = xk(nrow,k) + u(2,isect,kp,j)*axm(j,k)
595 xk(nrow,k) = xk(nrow,k)*wt(2,isect,kp)*rfact2
600 nrow = nrow + 1
610 continue
nr=1+nrowdone
nc = 1+ncoldone
do 650 isect = 1+nsectdone,numsect(nrot)+nsectdone
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) + a(k,m)*etai(isect,m,j)
do 645 kp = 1,numpts(2,isect)
do 640 k = 1,3
xk(nr,nc) = xk(nr,nc) + u(2,isect,kp,k)*axm(k,1)
640 xk(nr,nc+1) = xk(nr,nc+1) + u(2,isect,kp,k)*axm(k,2)
xk(nr,nc) = wt(2,isect,kp)*xk(nr,nc)*rfact2
xk(nr,nc+1) = wt(2,isect,kp)*xk(nr,nc+1)*rfact2
645 nr=nr+1
650 nc = nc + 2
nc = 1+ncoldone
do 680 isect = 1+nsectdone,numsect(nrot)+nsectdone
do 670 kp = 1,numpts(1,isect)
do 660 k = 1,3
xk(nr,nc) = xk(nr,nc) + u(1,isect,kp,k)*etai(isect,k,1)
660 xk(nr,nc+1) = xk(nr,nc+1) +
& u(1,isect,kp,k)*etai(isect,k,2)
xk(nr,nc) = xk(nr,nc)*wt(1,isect,kp)*rfact2
xk(nr,nc+1) = xk(nr,nc+1)*wt(1,isect,kp)*rfact2
670 nr=nr+1
680 nc = nc + 2
nrowdone=nrowdone+np(nrot)
nsectdone=nsectdone+numsect(nrot)
615 ncoldone=ncoldone+2*numsect(nrot)
c
c
c transform XK to S^(-1/2)*XK
c
nrowdone=0
ncoldone=3
do 510 k = 1,if
rkap=hkap(k)**.5
do 530 i = 1,np(k)
do 530 j = 1,3
530 xk(nrowdone+i,j)=rkap*xk(nrowdone+i,j)
do 520 j = ncoldone+1,ncoldone+2*numsect(k)
520 xk(nrowdone+i,j)=rkap*xk(nrowdone+i,j)
nrowdone=nrowdone+np(k)
510 ncoldone=ncoldone+2*numsect(k)
npts=nrowx
c
c open(unit=15,file='sx.dat',status='unknown')
c do 511 i = 1,nrowx
c 511 write(15,*)(xk(i,j),j=1,ncolxk)
c close(15)
c
c calculate P-one
c
call projection(xk,nrowx,ncolxk,pone)
c
c
c open(unit=15,file='pone.dat',status='unknown')
c do 513 i = 1,nrowx
c 513 write(15,*)(pone(i,j),j=1,nrowx)
c close(15)
c
c
c calculation of correction terms A and C
c
nrow=0
do 710 k=1,if
do 700 i=1,np(k)
do 700 j = 1,np(k)
700 pii(k,i,j)=pone(i+nrow,j+nrow)
710 nrow=nrow+np(k)
do 730 k=1,if
trpii(k)=0.
trpiisq(k)=0.
do 730 i = 1,np(k)
trpii(k)=trpii(k)+pii(k,i,i)
do 730 j = 1,np(k)
730 trpiisq(k)=trpiisq(k) + pii(k,i,j)*pii(k,j,i)
capa=0.
do 740 k = 1,if
740 capa=capa+(3+2*numsect(k)+trpiisq(k)-2*trpii(k))/dfr(k)
capc=0.
do 750 k = 1,if
750 capc=capc+(3+2*numsect(k)-trpii(k))**2/dfr(k)
capc=capc-capa
c
c calculate c and f, where Q is cF(3,f)
c
exq = 3 + 2*capa
vq = 6 + 14*capa + 2*capc
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
c write(6,*)'enter plev'
c read(6,*) plev
plev = .95
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
c
c adjust hatkapc for use in conreg.f
HATKAPC = 3./cfact
C
c calculate the inverse of covariance matrix H11.2
call inv(v,w)
write (*,*) ' rotation:'
write (*,*) ' '
write (*,*) 'Latitude, longitude, angle'
write (*,'(3x,3f8.2)') alat,alon,rho
write (*,*) 'Rotation matrix'
do 56 i=1,3
56 write (*,*) (a(i,j),j=1,3)
write (*,*) 'Covariance matrix'
do 60 i=1,3
60 write (*,*) (v(i,j),j=1,3)
c
65 write(*,1065)
1065 format(/,' Enter name of the output file: ',$)
read(*,'(a)') filename
open(unit=15,file=filename,status='unknown',err=65)
write(*,'(a)') filename
write(15,*)' results from averotkap'
write(15,*)' '
write(15,*) 'Fitted rotation--alat,along,rho: '
write(15,*) alat,alon,rho
write(15,*)'plev'
write(15,*) plev,1.0,1.0
c flag = 1 is output for use by the program addplus.f
write(15,*) '3/cfact, degrees of freedom, crit. value, flag'
write(15,1066) hatkapc,df,cxf,1.0
1066 format(f12.6,2x,f10.4,2x,f12.4,2x,f4.1)
write(15,*) 'Number of points, fmin, fmax, rank'
write(15,*) nrowx,fmin,fmax,ncolxk
write(15,*) 'ahat: '
do 70 i=1,3
70 write(15,*) (a(i,j),j=1,3)
write(15,*) 'covariance matrix'
do 75 i=1,3
75 write(15,*) (v(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)
stop
end
c************** readrot.f **********************************************
c subroutine to read a rotation file, with covariance matrix, kappahat ...
c
c input data file: file code 11
c
c 1. best fit rotation (ahat) in form: alat, along, rho
c
c 2. covariance matrix (except for division by kappahat) of ahat
c
c 3. kappahat (or kappa) and degrees of freedom
c
c 4. matrix h11.2
c
c234567
subroutine readrot(filname,ior,alat,alon,rho,ahat,b,h,plev,hatkap,
& df,ndum1)
c
implicit double precision(a-h,o-z)
real plev,df
character filname*40,dumline*40
dimension ahat(3,3),b(3,3),h(3,3)
1000 format (a)
c
ior=0
open(11,file=filname,status='old',err=20)
c
c read 2 first comment lines
c
read(11,1000) dumline
read(11,1000) dumline
c
c read bestfit rotation: alat,alon and angle
c
read(11,1000) dumline
read(11,*) alat,alon,rho
call trans7(alat,alon,rho,ahat)
c
c read confidence level
c
read(11,1000) dumline
read(11,*) plev
c
c read hatkap and degrees of freedom
c
read(11,1000) dumline
read(11,*) hatkap,df
c
c read nbr of points, sections
c
read(11,1000) dumline
write(*,*) dumline
read(11,*) ndum1
c
c read rotation matrix
c
read(11,1000) dumline
write(*,*) dumline
do 5 i=1,3
5 read (11,*) dum1
c
c read covariance matrix
c
read(11,1000) dumline
do 10 i=1,3
10 read(11,*) (b(i,j),j=1,3)
c
c read inverse of covariance matrix: H11.2
c
read(11,1000) dumline
do 15 i=1,3
15 read(11,*) (h(i,j),j=1,3)
c
close(11)
c
return
20 write (*,'(a)') 'Trouble opening filename: ',filname
ior=1
return
end
c=================================================================
c234567
subroutine rotvec(a,x,y)
c
c a*x=y
c
double precision a,x,y
dimension a(3,3),x(3),y(3)
do 1 i=1,3
y(i)=0.
do 1 j=1,3
1 y(i)=y(i)+a(i,j)*x(j)
return
end
c=================================================================
c234567
subroutine detr(a,d)
c
c calculate determinant of a
c
double precision a,d
dimension a(3,3)
d=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))
return
end
C
C
FUNCTION R1(HR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER (MSECT=70)
integer nsect
DIMENSION HR(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(HR,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