C
C *** HRCL8 ***
C (6/16/98)
C
C THIS PROGRAM FITS A SMOOTH HISTORY TO ROTATION DATA EXTENDING INTO THE
C PAST. IT ALLOWS FOR ELLIPTICAL ERROR DISTRIBUTIONS AROUND THE DATA
C POINTS, CONSIDERED AS POINTS ON S**3. THE INPUT DATA FILE SHOULD
C HAVE THE FOLLOWING FORMAT:
C FIRST RECORD: NDATA (NO. OF DATA POINTS)
C FOR EACH DATA POINT, THERE ARE FIVE RECORDS.
C THE FIRST OF THESE CONTAINS:
C TAU = TIME
C THETA = AXIS LATITUDE (DEGREES)
C PHI = AXIS LONGITUDE (DEGREES)
C RHO = ANGLE OF ROTATION (DEGREES)
C THE REMAINING FOUR RECORDS CONTAIN THE ROWS OF THE 4X4 SYMMETRIC
C MATRIX QTILDE CORRESP. TO THE DATA POINT. (SEE APPENDIX 3 IN
C HANNA AND CHANG (1990) FOR THE DEFINITION OF THIS MATRIX.) WE
C ASSUME THAT THE MATRIX QTILDE IN THE INPUT CORRESPONDS TO A
C CONFIDENCE REGION OF 95% LEVEL; THE 95% CONFIDENCE REGION
C CONSISTS OF ALL VECTORS Q IN S**3 SUCH THAT
C (Q**T)*QTILDE*Q
C IS LESS OR EQUAL TO 1.0. NOTE THAT THE MATRIX QTILDE HAS ONE ZERO
C EIGENVALUE AND THREE POSITIVE ONES. THE EIGENVECTOR OF QTILDE
C CORRESP. TO THE ZERO EIGENVALUE IS THE POINT ON S**3 CORRESP. TO
C THE ROTATION DETERMINED BY THE TRIPLE (THETA,PHI,RHO).
C NOTE: NDATA SHOULD BE AT MOST 50. TIMES SHOULD BE STRICTLY INCREASING.
C
C NOTE: UNROLLING IS ACCOMPLISHED BY LEFT-TRANSLATION. CORRESPONDING
C TO A PATH F(T) ON S**3, THE ROTATION PATH R(T) IN SO(4) IS
C DEFINED BY:
C
C R(T)*X = CONJ(F(T))*X .
C
C HERE X IS A VECTOR IN R**4. THE MULTIPLICATION ON THE LEFT IS
C MATRIX MULTIPLICATION AND THE MULTIPLICATION ON THE RIGHT IS
C QUATERNION MULTIPLICATION. R(T) MAPS F(T) INTO 1. IF X IS IN
C THE TANGENT SPACE TO S**3 AT F(T), THEN R(T)*X IS IN THE TANGENT
C SPACE TO S**3 AT 1.
C
C FROM THE TERMINAL, HRCL8 REQUESTS THE FOLLOWING INPUT:
C ANGDIV - MAX. LENGTH OF SEGMENT (IN DEGREES) FOR A PIECEWISE-GEODESIC
C REPRESENTATION OF A SMOOTH PATH ON S**3.
C ANGACC - ACCURACY (IN DEGREES) TO WHICH EACH POINT OF THE FITTED SPLINE
C PATH IS TO BE DETERMINED; MEASURED ON S**3.
C ISMTH -CODE FOR CHOICE OF THE SMOOTHING PARAMETER ALAMBDA:
C 0 - SPECIFY ALAMBDA
C 1 - DETERMINE ALAMBDA BY CROSS-VALIDATION
C 2 - DETERMINE ALAMBDA BY REQUIRING THE FITTED SPLINE
C POINTS TO LIE IN THE KNOWN CONFIDENCE REGIONS
C 3 - DETERMINE ALAMBDA BY THE DISCREPANCY METHOD
C (NOTE: ALAMBDA IS A POSITIVE NO. WHICH CONTROLS THE AMOUNT
C OF SMOOTHING. THE LARGER ALAMBDA, THE MORE SMOOTHING.)
C IF THE USER CHOOSES ISMTH = 0, THEN HE/SHE MUST SPECIFY
C THE VALUE OF ALAMBDA. THIS IS DONE BY FIRST TYPING A CODE
C ISLAMB:
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C IF ISLAMB = 1, THE USER THEN SPECIFIES THE FINITE VALUE
C DESIRED.
C
C AS OUTPUT, HRCL8 PRODUCES:
C (1) A TABLE OF THE FITTED SPLINE PATH. THIS TABLE CONSISTS OF NDIV
C LINES, WHERE NDIV = THE NO. OF POINTS (ON S**3) AT WHICH THE FITTED
C SPLINE IS CALCULATED. EACH LINE CONTAINS: TIME, AXIS LATITUDE (IN
C DEGREES), AXIS LONGITUDE (IN DEGREES), AND ANGLE OF ROTATION (IN
C DEGREES).
C (2) A TABLE OF INFORMATION ON THE ITERATIONS. EACH LINE OF THIS TABLE
C CONTAINS: ITERATION NO., APPROX. ACCURACY (IN DEGREES ON S**3) OF
C THE FITTED PATH, MAX. DIFFERENCE BETWEEN DATA AND FITTED SPLINE (IN
C DEGREES ON S**3), AND THE VALUE OF THE SMOOTHNESS OBJECTIVE FUNCTION
C (FOR THIS ITERATION). (NOTE: THE APPROX. ACCURACY OF THE FITTED
C PATH IS THE MAX. DIFFERENCE BETWEEN SUCCESSIVE FITTED SPLINES.)
C (3) A TABLE COMPARING THE DATA POINTS AND THE LAST FITTED SPLINE.
C EACH LINE CONTAINS: TIME, AXIS LATITUDE OF DATA POINT (DEGREES),
C AXIS LONGITUDE OF DATA POINT (DEGREES), ANGLE OF ROTATION OF DATA
C POINT (DEGREES), AND DISTANCE (IN DEGREES MEASURED ON S**3) BETWEEN
C THE FITTED SPLINE AND THE DATA POINT.
C
PARAMETER (NDATAT=50)
PARAMETER (NDIVT=40000)
PARAMETER (NUTOT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),THETA(NDATAT),PHI(NDATAT),RHO(NDATAT)
DIMENSION CAPQ(0:3,0:3,NDATAT), W(4,NDATAT)
DIMENSION T(NDIVT),Q(4,NDIVT)
DIMENSION ACC(50),ACCD(50),FOBJ(50)
DIMENSION ALAMB(50),XVFVAL(50)
DIMENSION ILAMB(50)
DIMENSION DISTW(NDATAT)
DIMENSION QVECT(4)
DIMENSION ELAMB(3,NDATAT),EQ(4,3,NDATAT)
DIMENSION WORKQ(4,4),DB(4),EB(4),E2(4),ZB(4,4)
DIMENSION QFIT(4,NDATAT),QTQQ(NDATAT)
DIMENSION SD(3,NDATAT)
DIMENSION KINDEX(NDATAT)
DIMENSION TSTAR(NUTOT)
LOGICAL MATZ
CHARACTER SW,ANS,SWCF,SWGPH
CHARACTER FLNAM1*20
CHARACTER FLNAM3*20
CHARACTER FNAME*4
CHARACTER FLNAM4*10,SUF4*6
CHARACTER FLNAM5*10,SUF5*6
C
C NOTE:
C W(4,NDATAT) - DATA POINTS ON S**3 (EUCLIDEAN COORDINATES)
C T(NDIVT) - TIMES CORRESP. TO THE NODES ON PIECEWISE-GEODESIC
C REPRESENTATION OF A SMOOTH CURVE ON S**3. (NDIV = NO. OF
C NODES; NDIVT = MAX. NO. OF NODES.)
C Q(4,NDIVT) - NODES OF THE FITTED SPLINE CURVE ON S**3.
C (IF G(T) IS THE FITTED CURVE ON S**3, THEN
C Q(K) = G(T(K)).)
C ACC(50) - VECTOR OF APPROX. ACCURACIES (DEGREES ON S**3) OF FITTED
C PATH; ONE ENTRY PER ITERATION.
C ACCD(50) - VECTOR OF MAX. DIFFERENCES (DEGREES ON S**3) BETWEEN DATA
C AND FITTED SPLINE; ONE ENTRY PER ITERATION.
C FOBJ(50) - VECTOR OF VALUES OF SMOOTHNESS OBJECTIVE FUNCTION; ONE
C ENTRY PER ITERATION.
C DISTW(NDATAT) - DISTANCES (DEGREES ON S**3) BETWEEN FITTED SPLINE AND
C DATA POINTS.
C
PRINT *,' '
PRINT *,'FITTING A SMOOTH HISTORY TO ROTATION DATA EXTENDING INTO'
PRINT *,'THE PAST (ELLIPTICAL ERRORS CASE, UNROLLING DEFINED'
PRINT *,'BY LEFT-TRANSLATION)'
C
PRINT *,' '
PRINT *,'TYPE SLONG = STARTING VALUE FOR REPORTING AXIS'
PRINT *,'LONGITUDES (AXIS LONGITUDES WILL BE REPORTED IN'
PRINT *,'THE RANGE [SLONG,SLONG + 360).)'
READ *,SLONG
C
C READ INPUT DATA.
C
PRINT *,' '
PRINT *,'TYPE DATA FILE NAME (AT MOST 20 CHARACTERS)'
READ 999,FLNAM1
OPEN(10,FILE=FLNAM1,STATUS='UNKNOWN')
REWIND (10)
READ(10,*) NDATA
DO 5 I=1,NDATA
READ(10,*) TAU(I),DUMMY,DUMMY,DUMMY
DO 8 J=0,3
8 READ(10,*) (CAPQ(J,M,I),M=0,3)
5 CONTINUE
CLOSE(10)
C
C NOTE THAT THE MATRICES QTILDE(I) ARE STORED IN THE ARRAY CAPQ.
C
C WE IGNORE THE VALUES OF THETA(I), PHI(I), RHO(I) FROM THE INPUT FILE.
C WE TAKE THE DATA POINT W(I) ON S**3 AS THE EIGENVECTOR OF QTILDE(I)
C CORRESP. TO THE ZERO EIGENVALUE. THE VALUES OF THETA(I), PHI(I), AND
C RHO(I) ARE THEN CALCULATED FROM W(I). (THE REASON FOR PROCEEDING IN
C THIS WAY IS TO BE ABLE TO CALCULATE W(I) MORE ACCURATELY.)
C
C ASSURE THAT THE MATRICES CAPQ(I) ARE SYMMETRIC.
C
DO 15 I=1,NDATA
DO 12 J=1,3
DO 12 M=0,J - 1
TEMP=(CAPQ(J,M,I) + CAPQ(M,J,I))/2.0D0
CAPQ(J,M,I)=TEMP
12 CAPQ(M,J,I)=TEMP
15 CONTINUE
C
C CALCULATE THE EIGENVALUES AND EIGENVECTORS OF THE MATRICES CAPQ(I).
C
C TO CALCULATE THE EIGENVALUES AND EIGENVECTORS OF THE MATRICES CAPQ(I),
C WE USE THE EISPACK ROUTINES BANDR AND IMTQL2.
C NM = ROW DIMENSION OF WORKQ = 4
C NEIG = ORDER OF CAPQ(I) = 4
C MB = NO. OF SUB-DIAGONALS OF CAPQ(I) + 1 = 4
C MATZ = .TRUE. (CALCULATE ORTH. TRANSF. MATRIX ZB)
C
NM=4
NEIG=4
MB=4
MATZ=.TRUE.
C
DO 80 I=1,NDATA
C
C SET UP THE ARRAY WORKQ. WORKQ CONTAINS THE MATRIX CAPQ(I) IN BAND
C SYMMETRIC STORAGE FORM (EISPACK VERSION). SPECIFICALLY, THE MAIN
C DIAGONAL OF CAPQ(I) GOES INTO THE LAST COLUMN OF WORKQ. THE FIRST
C SUB-DIAGONAL OF CAPQ(I) GOES INTO THE NEXT-TO-LAST COLUMN OF WORKQ
C (PUSHED TO THE BOTTOM OF THE COLUMN), ETC.
C
DO 50 J=1,4
DO 50 M=1,4
50 WORKQ(J,M)=0.0D0
DO 52 J=1,4
52 WORKQ(J,4)=CAPQ(J-1,J-1,I)
DO 56 K=1,3
DO 54 J=K+1,4
54 WORKQ(J,4-K)=CAPQ(J-1,J-K-1,I)
56 CONTINUE
C
CALL BANDR(NM,NEIG,MB,WORKQ,DB,EB,E2,MATZ,ZB)
CALL IMTQL2(NM,NEIG,DB,EB,ZB,IERR)
IF (IERR.NE.0) THEN
PRINT *,' '
PRINT *,'ERROR IN FINDING THE EIGENVALUES AND EIGENVECTORS'
PRINT *,'OF CAPQ(I), USING IMTQL2'
PRINT *,'I = ',I
PRINT *,'IERR = ',IERR
STOP
END IF
C
C THE VECTOR DB CONTAINS THE EIGENVALUES OF CAPQ(I) IN ASCENDING ORDER.
C THE COLUMNS OF THE ARRAY ZB CONTAIN THE CORRESP. ORTHONORMAL EIGEN-
C VECTORS.
C
C WE ASSURE THAT THE SMALLEST EIGENVALUE OF CAPQ(I) IS ZERO BY SUBTRACTING
C DB(1)*IDENTITY FROM CAPQ(I).
C
DO 60 J=0,3
60 CAPQ(J,J,I)=CAPQ(J,J,I) - DB(1)
C
DO 62 J=2,4
62 DB(J)=DB(J) - DB(1)
DB(1)=0.0D0
C
C WE MULTIPLY THE EIGENVALUES OF CAPQ(I) BY 7.815 (THE UPPER 5% VALUE
C FOR THE CHI-SQUARE DISTRIBUTION WITH 3 DEGREES OF FREEDOM) IN ORDER
C TO MAKE THE UNROLLED ERROR MATRICES SIGINV(I) (SEE THE ROUTINE ESSPLC
C BELOW) ACTUALLY INVERSES OF COVARIANCE MATRICES.
C
DO 70 L=1,3
70 ELAMB(L,I)=DB(L+1)*7.815D0
DO 72 L=1,3
DO 72 J=1,4
72 EQ(J,L,I)=ZB(J,L+1)
C
C TAKE THE DATA POINT W(I) AS THE EIGENVECTOR OF CAPQ(I) CORRESP. TO THE
C EIGENVALUE ZERO.
C
DO 75 J=1,4
75 W(J,I)=ZB(J,1)
C
C MAKE SURE THAT W(1,I) IS NON-NEGATIVE, SO THAT W(I) CORRESPONDS TO A
C ROTATION OF R**3 OF ANGLE BETWEEN ZERO AND 180 DEGREES.
C
IF (W(1,I).LT.0.0D0) THEN
DO 78 J=1,4
78 W(J,I)= -W(J,I)
END IF
80 CONTINUE
C
C NOTE: EQ(J,L,I) IS THE J-TH COMPONENT OF THE L-TH EIGENVECTOR OF THE
C ERROR MATRIX CORRESP. TO W(I).
C
C CALCULATE THE ARRAYS THETA(I),PHI(I),RHO(I), BASED ON THE DATA POINTS
C W(I).
C
PI=4.0D0*DATAN(1.0D0)
DO 90 I=1,NDATA
DO 95 J=1,4
95 QVECT(J)=W(J,I)
CALL CONV(QVECT,SLONG,PI,THETA(I),PHI(I),RHO(I))
90 CONTINUE
C
PRINT *,' '
PRINT *,'NO. OF DATA POINTS = ',NDATA
PRINT *,'TABLE OF DATA POINTS:'
PRINT 900
PRINT 901
DO 110 I=1,NDATA
110 PRINT 902,I,TAU(I),THETA(I),PHI(I),RHO(I)
900 FORMAT(' ',2X,'I',6X,'TIME',6X,' AXIS LAT.',2X,'AXIS LONG.',2X,
&' ROT. ANG.')
901 FORMAT(' ',21X,'(DEG.)',6X,'(DEG.)',6X,'(DEG.)')
902 FORMAT(' ',2X,I3,2X,G12.5,3(2X,F10.2))
C
C CALCULATE THE STANDARD DEVIATIONS OF THE ERROR DISTRIBUTIONS IN THE
C PRINCIPAL ERROR DIRECTIONS.
C
DO 112 I=1,NDATA
DO 112 L=1,3
112 SD(L,I)=(1.0D0/DSQRT(ELAMB(L,I)))*(180.0D0/PI)
C
PRINT *,' '
PRINT *,'STANDARD DEVIATIONS OF THE ERROR DISTRIBUTIONS, IN'
PRINT *,'THE PRINCIPAL ERROR DIRECTIONS (IN DEGREES ON S**3):'
PRINT *,' '
PRINT 925
DO 115 I=1,NDATA
115 PRINT 926,I,TAU(I),(SD(L,I),L=1,3)
925 FORMAT(' ',2X,'I',5X,'TAU(I)',15X,'STANDARD DEVIATIONS')
926 FORMAT(' ',I3,2X,G12.5,3(2X,G12.5))
C
125 PRINT *,' '
PRINT *,'TYPE ANGDIV = MAX. LENGTH OF SEGMENT (IN DEGREES'
PRINT *,'ON S**3) FOR A PIECEWISE-GEODESIC REPRESENTATION'
PRINT *,'OF A SMOOTH PATH. (SHOULD BE AT LEAST AS SMALL AS'
PRINT *,'THE STANDARD DEVIATIONS LISTED ABOVE.)'
READ *,ANGDIV
PRINT *,'TYPE ANGACC = ACCURACY (IN DEGREES ON S**3) TO WHICH'
PRINT *,'EACH POINT ON THE FITTED SPLINE PATH IS TO BE DETER-'
PRINT *,'MINED. (SHOULD BE AT LEAST AS SMALL AS THE STANDARD'
PRINT *,'DEVIATIONS LISTED ABOVE.)'
READ *,ANGACC
PRINT *,'TYPE CODE FOR METHOD OF DETERMINING THE SMOOTHING'
PRINT *,'PARAMETER ALAMBDA:'
PRINT *,' 0 - SPECIFY ALAMBDA'
PRINT *,' 1 - DETERMINE ALAMBDA BY CROSS-VALIDATION'
PRINT *,' 2 - DETERMINE ALAMBDA BY REQUIRING THE FITTED SPLINE'
PRINT *,' POINTS TO LIE IN THE KNOWN CONFIDENCE REGIONS'
PRINT *,' 3 - DETERMINE ALAMBDA BY THE DISCREPANCY METHOD'
READ *,ISMTH
IF (ISMTH.EQ.0) THEN
PRINT *,' '
PRINT *,'TYPE CODE FOR VALUE OF ALAMBDA:'
PRINT *,' 1 - FINITE VALUE'
PRINT *,' 2 - ALAMBDA = INFINITY'
READ *,ISLAMB
IF (ISLAMB.EQ.1) THEN
PRINT *,'TYPE VALUE OF ALAMBDA'
READ *,ASLAMB
END IF
ELSE IF (ISMTH.EQ.1) THEN
PRINT *,' '
PRINT *,'TO DETERMINE THE OPTIMAL VALUE OF ALAMBDA, THE'
PRINT *,'PROGRAM CALCULATES CROSS-VALIDATION SCORES ON A GRID'
PRINT *,'OF MMAX+1 ALAMBDA-VALUES. YOU WILL BE ASKED TO'
PRINT *,'SPECIFY MMAX AT EACH ITERATION. (MMAX = 10 IS A'
PRINT *,'REASONABLE VALUE.)'
ELSE IF (ISMTH.EQ.2) THEN
PRINT *,' '
PRINT *,'UNDER THIS OPTION, THE PROGRAM USES THE LARGEST'
PRINT *,'VALUE OF ALAMBDA SUCH THAT ALL BUT K OF THE FITTED'
PRINT *,'SPLINE POINTS LIE IN THE CORRESP. CONFIDENCE'
PRINT *,'REGIONS. YOU WILL BE ASKED TO SPECIFY K LATER.'
PRINT *,'(K MUST BE BETWEEN 0 AND NDATA-1, INCLUSIVE.)'
ELSE IF (ISMTH.EQ.3) THEN
PRINT *,' '
PRINT *,'UNDER THIS OPTION, THE PROGRAM USES THE LARGEST'
PRINT *,'VALUE OF ALAMBDA SUCH THAT:'
PRINT *,' '
PRINT *,' (SUM OF'
PRINT *,' (D(ALAMBDA,I)**T)*SIGINV(I)*D(ALAMBDA,I))'
PRINT *,' /(3*NDATA)'
PRINT *,' '
PRINT *,'IS LESS OR EQUAL 1.0. HERE SIGINV(I) IS THE'
PRINT *,'INVERSE OF THE COVARIANCE MATRIX AT THE I-TH DATA'
PRINT *,'POINT AND D(ALAMBDA,I) IS THE DISCREPANCY AT THE'
PRINT *,'I-TH DATA POINT, I.E.,'
PRINT *,' D(ALAMBDA,I) = WSTAR(I) - G(ALAMBDA,TAU(I)) ,'
PRINT *,'WHERE G(ALAMBDA,T) IS THE FITTED SPLINE.'
END IF
C
PRINT *,' '
PRINT *,'TYPE MAX. NO. OF ITERATIONS (5 IS USUALLY SUFFICIENT)'
READ *,ITERMX
C
PRINT *,' '
PRINT *,'TYPE CODE FOR PRINTING OPTION:'
PRINT *,' 0 - STANDARD'
PRINT *,' 105 - THE MATRICES QTILDE(I)'
PRINT *,' 110 - DATA POINTS IN EUCLIDEAN COORDINATES (ON S**3)'
PRINT *,' 115 - POS. EIGENVALUES OF THE MATRICES QTILDE(I)'
PRINT *,' AND THEIR CORRESP. EIGENVECTORS'
PRINT *,' 116 - THE MATRIX OMEGA'
PRINT *,' 120 - TIMES AND NODES OF PIECEWISE-GEODESIC APPROX.'
PRINT *,' (INTIAL CURVE)'
PRINT *,' 140 - THE ARRAYS VSTAR, WSTAR, EQSTAR FROM UNRLLC'
PRINT *,' 145 - THE UNWRAPPED ERROR MATRICES SIGINV'
PRINT *,' 160 - NODES OF FITTED SPLINE IN R**3'
PRINT *,' 170 - THE ARRAY Q FROM WRAPC: FITTED POINTS WRAPPED'
PRINT *,' ONTO S**3'
PRINT *,' 201 - THE MATRIX OMEGA (FROM SEFIT)'
PRINT *,' 202 - THE ARRAYS BMAT, ETAU (FROM SEFIT)'
PRINT *,' 203 - THE ARRAY GMAT (FROM SEFIT)'
PRINT *,' 204 - THE ARRAYS Z, H (FROM SEFIT)'
PRINT *,' 205 - THE ARRAY SMAT (FROM SEFIT)'
PRINT *,' 206 - THE EIGENVALUES AND EIGENVECTORS'
PRINT *,' OF THE MATRIX OMEGA-S (FROM SEFIT)'
PRINT *,' 207 - THE VECTOR DELTA (FROM SEFIT)'
PRINT *,' 208 - THE MATRIX ALPHA (FROM SEFIT)'
PRINT *,' 209 - THE SOL. VECTOR EPS (FROM SEFIT)'
PRINT *,' 210 - THE ARRAY GHAT (FROM SEFIT)'
PRINT *,' 211 - THE ARRAY YFIT (FROM SEFIT)'
PRINT *,' 221 - DETAILS OF THE CALC. OF XVSC (FROM EVXV)'
PRINT *,' 231 - THE MATRIX OMEGA (FROM CLCEIG)'
PRINT *,' 232 - THE MATRIX B (FROM CLCEIG)'
PRINT *,' 233 - THE MATRIX G (FROM CLCEIG)'
PRINT *,' 234 - THE EIGENVALUES OF OMEGA (FROM CLCEIG)'
PRINT *,' 235 - THE EIGENVALUES OF G (FROM CLCEIG)'
PRINT *,' 251 - THE RATIO RHO (FROM FOLAM)'
PRINT *,' 253 - TABLE OF EST. RECIPS. OF COND. NOS. (FROM'
PRINT *,' FOLAM)'
PRINT *,' 261 - THE ARRAYS GHAT AND WFIT (FROM ESMTH)'
PRINT *,' 262 - THE OBJ. FN. AND TERMS MAKING IT UP (FROM'
PRINT *,' ESMTH)'
PRINT *,' 270 - DETAILS OF ITERATION IN FOLAM2'
PRINT *,' 275 - DETAILS OF ITERATION IN FOLAM3'
PRINT *,' 280 - THE ARRAYS KSTAR(NUTOT) AND GFHAT(4,NUTOT),'
PRINT *,' IN CASE TIMES TSTAR(NU) ARE SPECIFIED FOR'
PRINT *,' POSTERIOR PROBABILITY REGIONS'
READ *,IPRINT
C
IF (IPRINT.EQ.105) THEN
PRINT *,' '
PRINT *,'THE MATRICES QTILDE(I):'
DO 118 I=1,NDATA
PRINT *,' '
PRINT 910,I
DO 120 J=0,3
120 PRINT 911,(CAPQ(J,M,I),M=0,3)
118 CONTINUE
910 FORMAT(' ','DATA POINT NO. ',I3)
911 FORMAT(' ',4(2X,G12.5))
END IF
C
IF (IPRINT.EQ.110) THEN
PRINT *,' '
PRINT *,'DATA POINTS IN EUCLIDEAN COORDINATES (ON S**3):'
PRINT 920
DO 140 I=1,NDATA
140 PRINT 921,I,(W(J,I),J=1,4)
920 FORMAT(' ',2X,'I',5X,'W(1,I)',8X,'W(2,I)',8X,'W(3,I)',8X,
&'W(4,I)')
921 FORMAT(' ',I3,4(2X,G12.5))
END IF
C
IF (IPRINT.EQ.115) THEN
PRINT *,' '
PRINT *,'THE POSITIVE EIGENVALUES OF THE MATRICES QTILDE(I),'
PRINT *,'AND THEIR CORRESP. EIGENVECTORS:'
DO 190 I=1,NDATA
PRINT 930,I
PRINT *,' EIGENVALUES'
PRINT 931,(ELAMB(L,I)/7.815D0,L=1,3)
PRINT *,' EIGENVECTORS (THE COLUMNS BELOW)'
DO 165 J=1,4
165 PRINT 931,(EQ(J,L,I),L=1,3)
190 CONTINUE
930 FORMAT(' ','DATA POINT NO. ',I3)
931 FORMAT(' ',3(2X,G12.5))
END IF
C
C THE USER HAS THE OPTION OF WRITING INFORMATION TO A FILE, TO PREPARE
C FOR CALCULATING POSTERIOR PROBABILITY REGIONS.
C
PRINT *,' '
PRINT *,'DO YOU WANT TO WRITE INFORMATION TO A FILE, TO PREPARE'
PRINT *,'FOR CALCULATING POSTERIOR PROBABILITY REGIONS? (Y OR N)'
READ 999,SWCF
IF ((SWCF.EQ.'Y').OR.(SWCF.EQ.'y')) THEN
ICF=1
PRINT *,' '
PRINT *,'TYPE NAME OF FILE TO RECEIVE INFORMATION (AT MOST 20'
PRINT *,'CHARACTERS)'
READ 999,FLNAM3
OPEN(30,FILE=FLNAM3,STATUS='UNKNOWN')
REWIND (30)
WRITE(30,932) FLNAM1
WRITE(30,933) SLONG,ANGDIV,ANGACC
932 FORMAT(' ','FITTED SPLINE TO DATA FROM FILE: ',A20)
933 FORMAT(' ','SLONG = ',G12.5,'; ANGDIV = ',G12.5,'; ANGACC = ',
%G12.5)
ELSE
ICF=0
END IF
C
IF (ICF.EQ.1) THEN
PRINT *,' '
PRINT *,'YOU MUST SPECIFY A SEQUENCE OF TIMES TSTAR(NU)'
PRINT *,'IN THE INTERVAL [TAU(1),TAU(NDATA)], FOR WHICH'
PRINT *,'YOU WANT POSTERIOR PROBABILITY REGIONS FOR THE'
PRINT *,'POINTS G(TSTAR(NU)).'
PRINT *,' '
PRINT *,'TYPE THE NUMBER OF TIMES TSTAR(NU) YOU WANT TO'
PRINT *,'SPECIFY (AT LEAST 1, AT MOST 50):'
READ *,NUMAX
PRINT *,'TYPE THE TIMES TSTAR(NU) WHICH YOU WANT, IN'
PRINT *,'INCREASING ORDER, ONE TIME PER LINE:'
PRINT 1010,TAU(1),TAU(NDATA)
1010 FORMAT(' ',' (NOTE: THE TIMES TSTAR(NU) MUST BE BETWEEN ',
&'TAU(1) = ',G12.5,/' ',' AND TAU(NDATA) = ',G12.5,
&', INCLUSIVE.)')
DO 310 NU=1,NUMAX
310 READ *,TSTAR(NU)
C
PRINT *,' '
PRINT *,'THE TIMES FOR POSTERIOR PROBABILITY REGIONS:'
PRINT *,' '
PRINT 1015
DO 315 NU=1,NUMAX
315 PRINT 1016,NU,TSTAR(NU)
1015 FORMAT(' ',1X,'NU',4X,'TSTAR(NU)')
1016 FORMAT(' ',I3,3X,G12.5)
END IF
C
PRINT *,' '
PRINT *,'DO YOU WANT TO CREATE FILES SUITABLE FOR GRAPHING THE'
PRINT *,'FITTED CURVE ALONG WITH THE INPUT ROTATIONS AND'
PRINT *,'CONFIDENCE REGIONS? (Y OR N)'
READ 999,SWGPH
IF ((SWGPH.EQ.'Y').OR.(SWGPH.EQ.'y')) THEN
IGPH=1
PRINT *,'FOR PURPOSES OF GRAPHING, THE PROGRAM HRCL8 WRITES'
PRINT *,'TWO FILES:'
PRINT *,' XXXX.INPUT - CONTAINS INPUT ROTATIONS AND THEIR'
PRINT *,' Q-MATRICES'
PRINT *,' XXXX.CURVE - CONTAINS THE FITTED CURVE'
PRINT *,'TYPE THE NAME YOU WANT TO REPLACE "XXXX" (USE 4'
PRINT *,'NON-BLANK CHARACTERS)'
READ 999,FNAME
SUF4='.input'
FLNAM4=FNAME//SUF4
SUF5='.curve'
FLNAM5=FNAME//SUF5
ELSE
IGPH=0
END IF
C
C CALL ESSPLC.
C
IFAIL=0
CALL ESSPLC(NDATA,TAU,W,ELAMB,EQ,ANGDIV,ANGACC,ISMTH,ISLAMB,
&ASLAMB,ITERMX,ICF,NUMAX,TSTAR,SLONG,IGPH,FNAME,IPRINT,
&NDIV,T,Q,NITER,ACC,ACCD,FOBJ,ILAMB,ALAMB,XVFVAL,DISTW,ACC1,
&QFIT,KINDEX,IFAIL)
IF (IFAIL.NE.0) THEN
GO TO 290
END IF
C
IF (ICF.EQ.1) THEN
CLOSE(30)
PRINT *,' '
PRINT *,'INFORMATION FOR CALCULATING POSTERIOR PROBABILITY'
PRINT *,'REGIONS HAS BEEN WRITTEN TO THE FILE: '
PRINT 934,FLNAM3
934 FORMAT(' ',5X,A20)
END IF
C
PRINT *,' '
PRINT *,'NO. OF POINTS AT WHICH FITTED SPLINE IS CALC. = '
PRINT *,NDIV
C
PRINT *,' '
198 PRINT *,'TYPE C TO CONTINUE.'
READ 999,ANS
IF ((ANS.EQ.'C').OR.(ANS.EQ.'c')) THEN
GO TO 199
ELSE
GO TO 198
END IF
199 CONTINUE
C
PRINT *,' '
PRINT *,'TABLE OF FITTED SPLINE:'
PRINT 935,NDIV
935 FORMAT(' ',2X,'THERE ARE ',I5,' POINTS IN THE FITTED SPLINE.')
PRINT *,' YOU MAY PRINT ONLY EVERY M-TH POINT BETWEEN DATA'
PRINT *,' TIMES. PLEASE TYPE M.'
200 READ *,M
IF (M.LT.1) THEN
PRINT *,' '
PRINT *,'M MUST BE AT LEAST 1. PLEASE RETYPE.'
GO TO 200
END IF
PRINT 938
PRINT 940
DO 202 I=1,NDATA
IF (I.LT.NDATA) THEN
KSTART=KINDEX(I)
KSTOP=KINDEX(I+1) - 1
ELSE
KSTART=NDIV
KSTOP=NDIV
END IF
DO 202 K=KSTART,KSTOP
IF (MOD(K-KSTART,M).EQ.0) THEN
DO 204 J=1,4
204 QVECT(J)=Q(J,K)
CALL CONV(QVECT,SLONG,PI,ALAT,ALONG,RHOS)
PRINT 945,K,T(K),ALAT,ALONG,RHOS
END IF
202 CONTINUE
938 FORMAT(' ',3X,'K',6X,'TIME',6X,' AXIS LAT.',2X,'AXIS LONG.',
&2X,' ROT. ANG.')
940 FORMAT(' ',22X,'(DEG.)',6X,'(DEG.)',6X,'(DEG.)')
945 FORMAT(' ',1X,I5,2X,G12.5,3(2X,F10.2))
C
PRINT *,' '
PRINT *,'NO. OF ITERATIONS NEEDED FOR CONVERGENCE = ',NITER
PRINT *,' '
PRINT *,'INFORMATION ON ITERATIONS:'
PRINT *,'NOTE:'
PRINT *,' APPROX.ACC. = APPROX. ACCURACY (IN DEGREES ON S**3)'
PRINT *,' OF THE FITTED PATH; MEASURED BY TAKING THE MAX.'
PRINT *,' DIFFERENCE BETWEEN SUCCESSIVE FITTED SPLINES.'
PRINT *,' MAX.DIFF. = MAX. DIFFERENCE BETWEEN DATA AND THE'
PRINT *,' FITTED SPLINE (IN DEGREES ON S**3).'
PRINT 950
DO 210 ITER=1,NITER
210 PRINT 955,ITER,ACC(ITER),ACCD(ITER),FOBJ(ITER)
950 FORMAT(' ','ITER',3X,'APPROX.ACC.',9X,'MAX.DIFF.',17X,'OBJ.FN.')
955 FORMAT(' ',2X,I2,2X,G12.5,2(2X,G22.15))
C
PRINT *,' '
PRINT *,'VALUES OF ALAMBDA USED, TOGETHER WITH CORRESP.'
PRINT *,'CROSS-VALIDATION SCORES'
PRINT 957
DO 215 ITER=1,NITER
IF (ILAMB(ITER).EQ.1) THEN
PRINT 958,ITER,ALAMB(ITER),XVFVAL(ITER)
ELSE
PRINT 959,ITER,XVFVAL(ITER)
END IF
215 CONTINUE
957 FORMAT(' ','ITER',5X,'ALAMBDA',7X,'XVFVAL')
958 FORMAT(' ',I4,2X,G12.5,2X,G22.15)
959 FORMAT(' ',I4,2X,' INFINITY ',2X,G22.15)
C
PRINT *,' '
PRINT *,'COMPARISON OF DATA POINTS AND LAST FITTED SPLINE'
PRINT *,'NOTE: AXIS-LATITUDES, AXIS-LONGITUDES, AND ROTATION-'
PRINT *,'ANGLES ARE FOR THE DATA POINTS. DIST. = DISTANCE'
PRINT *,'(IN DEGREES ON S**3) BETWEEN FITTED SPLINE AND DATA'
PRINT *,'POINT.'
PRINT 960
PRINT 965
DO 220 I=1,NDATA
PHIM=CLONG(PHI(I),SLONG)
220 PRINT 970,I,TAU(I),THETA(I),PHIM,RHO(I),DISTW(I)
960 FORMAT(' ',2X,'I',6X,'TIME',6X,' AXIS LAT.',2X,'AXIS LONG.',
&2X,' ROT. ANG.',3X,'DIST.')
965 FORMAT(' ',21X,'(DEG.)',6X,'(DEG.)',6X,'(DEG.)',5X,'(S**3)')
970 FORMAT(' ',I3,2X,G12.5,3(2X,F10.2),2X,F7.2)
C
C FOR EACH DATA TIME TAU(I), CALCULATE THE QUANTITY:
C (QHAT**T)*QTILDE*QHAT ,
C WHERE QHAT IS THE FITTED SPLINE POINT AND QTILDE IS THE 4X4 MATRIX
C DEFINING THE 95% CONFIDENCE REGION.
C
DO 230 I=1,NDATA
SUM=0.0D0
DO 225 K=1,4
DO 225 L=1,4
225 SUM=SUM + QFIT(K,I)*CAPQ(K-1,L-1,I)*QFIT(L,I)
230 QTQQ(I)=SUM
C
PRINT *,' '
PRINT *,'MORE INFORMATION ON THE FITTED SPLINE POINTS IN RELATION'
PRINT *,'TO THE DATA POINTS: FOR EACH DATA TIME TAU(I), WE GIVE'
PRINT *,'THE QUANTITY:'
PRINT *,' (QHAT**T)*QTILDE*QHAT .'
PRINT *,'THE FITTED SPLINE POINT QHAT LIES WITHIN THE 95% CONFI-'
PRINT *,'DENCE REGION AT TAU(I) IF AND ONLY IF THIS QUANTITY IS'
PRINT *,'LESS OR EQUAL 1.0.'
PRINT 972
DO 235 I=1,NDATA
235 PRINT 973,I,TAU(I),QTQQ(I)
972 FORMAT(' ',2X,'I',5X,'TAU(I)',7X,'QTQQ(I)')
973 FORMAT(' ',I3,2X,G12.5,2X,G12.5)
C
PRINT *,' '
PRINT 975,ACC1
975 FORMAT(' ','MAX. DIFFERENCE BETWEEN FIRST AND FINAL FITTED',
&' SPLINES'/' ','(IN DEGREES ON S**3) = ',G12.5)
C
C If files suitable for graphing have been requested, they are wriiten
C here.
C
IF (IGPH.EQ.1) THEN
C
C Writing the file XXXX.input. This file contains the input rotations
C and their Q-matrices.
C
OPEN(40,FILE=FLNAM4,STATUS='UNKNOWN')
REWIND (40)
DO 240 I=1,NDATA
WRITE(40,*) TAU(I),THETA(I),PHI(I),RHO(I)
DO 242 L=0,3
242 WRITE(40,*) (7.815D0*CAPQ(L,M,I),M=0,3)
240 CONTINUE
CLOSE (40)
C
C Writing the file XXXX.curve. This file contains the fitted curve.
C
OPEN(50,FILE=FLNAM5,STATUS='UNKNOWN')
REWIND (50)
PRINT *,' '
PRINT 978,FLNAM5
978 FORMAT(' ','PREPARING TO WRITE THE FITTED CURVE TO THE',
&' FILE: ',A10)
PRINT 979,NDIV
979 FORMAT(' ','THERE ARE ',I5,' POINTS IN THE FITTED SPLINE.')
PRINT *,' YOU MAY WRITE ONLY EVERY M-TH POINT BETWEEN'
PRINT *,' DATA TIMES. PLEASE TYPE M.'
243 READ *,M
IF (M.LT.1) THEN
PRINT *,' '
PRINT *,'M MUST BE AT LEAST 1. PLEASE RETYPE.'
GO TO 243
END IF
DO 250 I=1,NDATA
IF (I.LT.NDATA) THEN
KSTART=KINDEX(I)
KSTOP=KINDEX(I+1) - 1
ELSE
KSTART=NDIV
KSTOP=NDIV
END IF
DO 248 K=KSTART,KSTOP
IF (MOD(K-KSTART,M).EQ.0) THEN
DO 245 J=1,4
245 QVECT(J)=Q(J,K)
CALL CONV(QVECT,SLONG,PI,FLAT,FLONG,FRHO)
WRITE(50,*) T(K),FLAT,FLONG,FRHO
END IF
248 CONTINUE
250 CONTINUE
CLOSE (50)
C
PRINT *,' '
PRINT *,'INFORMATION FOR GRAPHING HAS BEEN WRITTEN TO THE'
PRINT *,'FILES:'
PRINT 980,FLNAM4,FLNAM5
980 FORMAT(' ',5X,A10,' AND ',A10)
END IF
C
290 CONTINUE
PRINT *,'DO YOU WANT TO REPEAT THE FITTING PROCESS USING'
PRINT *,'DIFFERENT PARAMETERS? (Y OR N)'
READ 999,SW
999 FORMAT (A)
IF ((SW.EQ.'Y').OR.(SW.EQ.'y')) GO TO 125
C
STOP
END
C
C *** ESSPLC ***
C
C THIS SUBROUTINE DOES THE ACTUAL CALCULATION FOR THE DRIVER PROGRAM.
C (THE MAIN FUNCTION OF HRCL8 IS TO READ IN THE INPUT DATA AND TO PRINT
C THE RESULTS.)
C
C INPUTS:
C NDATA - NO. OF DATA POINTS
C TAU(NDATAT) - TIMES OF THE DATA POINTS
C W(4,NDATAT) - THE DATA POINTS, EXPRESSED AS POINTS ON S**3
C ELAMB(3,NDATAT) - THE POS. EIGENVALUES OF THE ERROR MATRICES ASSOC.
C WITH THE DATA POINTS
C EQ(4,3,NDATAT) - THE EIGENVECTORS OF THE ERROR MATRICES CORRESP. TO
C THE POSITIVE EIGENVALUES
C ANGDIV - MAX. LENGTH OF SEGMENT (IN DEGREES) OF A PIECEWISE-
C GEODESIC REPRESENTATION OF A SMOOTH PATH ON S**3
C ANGACC - ACCURACY (IN DEGREES) TO WHICH EACH POINT OF THE
C FITTED SPLINE PATH IS TO BE DETERMINED
C ISMTH - CODE FOR CHOICE OF THE SMOOTHING PARAMETER ALAMBDA:
C 0 - ALAMBDA IS SPECIFIED
C 1 - DETERMINE ALAMBDA BY CROSS-VALIDATION
C 2 - DETERMINE ALAMBDA BY REQUIRING THE FITTED
C SPLINE POINTS TO LIE IN THE KNOWN CONFIDENCE
C REGIONS
C 3 - DETERMINE ALAMBDA BY THE DISCREPANCY METHOD
C ISLAMB - IF ISMTH = 0, ISLAMB CONTAINS CODE FOR VALUE OF
C ALAMBDA:
C 1 - FINITE VALUE
C 2 - ALAMBDA = INFINITY
C ASLAMB - IF ISMTH = 0 AND ISLAMB = 1, ASLAMB CONTAINS THE
C VALUE OF ALAMBDA DESIRED
C ITERMX - MAX. NO. OF ITERATIONS
C ICF - CODE FOR OPTION OF WRITING INFORMATION FOR POSTERIOR
C PROBABILITY REGIONS TO FILE 30:
C 0 - DO NOT WRITE
C 1 - DO WRITE
C NUMAX - NUMBER OF TIMES TSTAR(NU) AT WHICH TO CALCULATE
C POSTERIOR PROBABILTY REGIONS FOR G(TSTAR(NU)),
C ASSUMING ICF=1
C TSTAR(NUTOT) - TABLE OF TIMES TSTAR(NU), IN CASE ICF=1
C SLONG - STARTING VALUE FOR REPORTING AXIS LONGITUDES (AXIS
C LONGITUDES WILL BE REPORTED IN THE RANGE [SLONG,
C SLONG + 360)
C IGPH - CODE FOR OPTION OF WRITING INFORMATION FOR GRAPHING
C TO FILES
C 0 - DO NOT WRITE
C 1 - DO WRITE
C FNAME - NAME TO BE USED IN THE PREFIX OF FILE NAMES FOR FILES
C TO BE USED FOR GRAPHING (4 NON-BLANK CHARACTERS)
C IPRINT - CODE FOR PRINTING OPTION
C
C OUTPUTS:
C NDIV - NO. OF NODES IN PIECEWISE-GEODESIC REPRESENTATION
C OF SMOOTH PATH ON S**3
C T(NDIVT) - TIMES OF THE NODES IN PIECEWISE-GEODESIC REPRESENTATION
C OF SMOOTH PATH ON S**3
C Q(4,NDIVT) - NODES OF FITTED SPLINE PATH ON S**3
C (IF G(T) IS THE FITTED CURVE ON S**3, THEN
C Q(K) = G(T(K)).)
C NITER - NO. OF ITERATIONS NEEDED FOR CONVERGENCE
C ACC(50) - VECTOR OF APPROX. ACCURACIES (IN DEGREES ON S**3)
C OF FITTED PATHS; ONE ENTRY PER ITERATION.
C ACCD(50) - VECTOR OF MAX. DIFFERENCES (IN DEGREES ON S**3)
C BETWEEN DATA AND FITTED SPLINE; ONE ENTRY PER ITERATION.
C FOBJ(50) - VECTOR OF VALUES OF THE SMOOTHNESS OBJECTIVE FUNCTION;
C ONE ENTRY PER ITERATION
C ILAMB(50) - VECTOR OF INDICATORS FOR THE VALUE OF ALAMBDA USED,
C ONE ENTRY PER ITERATION. (ILAMB(ITER) = 1 MEANS
C FINITE VALUE; ILAMB(ITER) = 2 MEANS ALAMBDA = INFINITY.)
C ALAMB(50) - VECTOR OF SMOOTHING PARAMETER VALUES; ONE ENTRY PER
C ITERATION. (IF ILAMB(ITER) = 1, THEN ALAMB(ITER)
C CONTAINS THE FINITE VALUE OF ALAMBDA USED.)
C XVFVAL(50) - VECTOR OF CROSS-VALIDATION SCORES; ONE ENTRY PER
C ITERATION
C DISTW(NDATAT) - DISTANCES (IN DEGREES ON S**3) BETWEEN FINAL FITTED
C SPLINE AND DATA POINTS.
C ACC1 - MAX. DIFFERENCE BETWEEN FIRST AND FINAL FITTED SPLINE
C (IN DEGREES ON S**3).
C QFIT(4,NDATAT) - FITTED SPLINE POINTS ON S**3 CORRESP. TO DATA TIMES
C KINDEX(NDATAT) - INTEGER ARRAY: FOR EACH I, KINDEX(I) IS THE INDEX
C OF THE NODE ON THE FITTED SPLINE WHICH CORRESP. IN
C TIME TO THE I-TH DATA POINT
C IFAIL - ERROR INDICATOR.
C
SUBROUTINE ESSPLC(NDATA,TAU,W,ELAMB,EQ,ANGDIV,ANGACC,ISMTH,ISLAMB,
&ASLAMB,ITERMX,ICF,NUMAX,TSTAR,SLONG,IGPH,FNAME,IPRINT,
&NDIV,T,Q,NITER,ACC,ACCD,FOBJ,ILAMB,ALAMB,XVFVAL,DISTW,ACC1,
&QFIT,KINDEX,IFAIL)
PARAMETER (NDATAT=50)
PARAMETER (NDIVT=40000)
PARAMETER (NUTOT=50)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION TAU(NDATAT),W(4,NDATAT)
DIMENSION ELAMB(3,NDATAT),EQ(4,3,NDATAT)
DIMENSION T(NDIVT),Q(4,NDIVT)
DIMENSION ACC(50),ACCD(50),FOBJ(50)
DIMENSION ALAMB(50),XVFVAL(50)
DIMENSION ILAMB(50)
DIMENSION DISTW(NDATAT)
DIMENSION OM(0:3,NDATAT)
DIMENSION KINDEX(NDATAT)
DIMENSION V(4,NDIVT),VSTAR(3,NDIVT)
DIMENSION WSTAR(3,NDATAT)
DIMENSION EQSTAR(4,3,NDATAT),SIGINV(3,3,NDATAT)
DIMENSION GHAT(3,NDATAT),WFIT(3,NDATAT)
DIMENSION QSTAR(3,NDIVT),QFIRST(4,NDIVT)
DIMENSION ETAU(-2:NDATAT+3)
DIMENSION QFIT(4,NDATAT)
DIMENSION TSTAR(NUTOT),GFHAT(4,NUTOT)
DIMENSION KSTAR(NUTOT)
CHARACTER FNAME*4
C
PI=4.0D0*DATAN(1.0D0)
C
C CALCULATE THE MATRIX OMEGA.
C
CALL CALCOM(NDATA,TAU,OM)
IF (IPRINT.EQ.116) THEN
PRINT *,' '
PRINT *,'THE DIAGONALS OF THE MATRIX OMEGA'
PRINT 900
DO 300 I=1,NDATA-3
300 PRINT 901,I,OM(0,I),OM(1,I),OM(2,I),OM(3,I)
DO 310 I=NDATA-2,NDATA
310 PRINT 901,I,(OM(K,I),K=0,NDATA-I)
900 FORMAT(' ',2X,'I',4X,'OMEGA(I,I)',2X,'OMEGA(I,I+1)',
& 2X,'OMEGA(I,I+2)',2X,'OMEGA(I,I+3)')
901 FORMAT(' ',I3,4(2X,G12.5))
END IF
C
C DIVIDE THE TIME INTERVALS COMING FROM THE DATA INTO SHORT SEGMENTS
C AND CONSTRUCT AN INITIAL BASE CURVE BY CONNECTING THE NODES BY
C GEODESIC ARCS.
C
AD=ANGDIV*(PI/180.D0)
NDIV=1
T(1)=TAU(1)
DO 10 J=1,4
10 V(J,1)=W(J,1)
KINDEX(1)=1
C
DO 40 I=1,NDATA-1
C
C CALCULATE DIST = ANGULAR DISTANCE BETWEEN W(I) AND W(I+1), MEASURED
C IN RADIANS ON S**3.
C
PROD=0.0D0
DO 20 J=1,4
20 PROD=PROD + W(J,I)*W(J,I+1)
DIST=DACOS(PROD)
C
C NSPLIT = THE NO. OF SUB-INTERVALS INTO WHICH WE WILL DIVIDE THE
C INTERVAL (TAU(I),TAU(I+1)).
C
NSPLIT=DIST/AD + 1
IF ((NDIV+NSPLIT).GT.NDIVT) THEN
IFAIL=1
PRINT *,'NO. OF NODES ON PIECEWISE-GEODESIC PATH EXCEEDS'
PRINT *,'NDIVT. EITHER ANGDIV OR NDIVT MUST BE INCREASED.'
RETURN
END IF
C
C CALCULATE THE TIMES T(NDIV+J) AND NODES V(NDIV+J) OF THAT PART
C OF THE PIECEWISE-GEODESIC PATH WHICH CORRESPONDS TO THE INTERVAL
C (TAU(I),TAU(I+1)).
C
SD=DSIN(DIST)
DELT=(TAU(I+1) - TAU(I))/NSPLIT
DO 30 J=1,NSPLIT
T(NDIV+J)=TAU(I+1) - (NSPLIT - J)*DELT
F1=DBLE(NSPLIT-J)/DBLE(NSPLIT)
F2=DBLE(J)/DBLE(NSPLIT)
ST1=DSIN(F1*DIST)/SD
ST2=DSIN(F2*DIST)/SD
DO 35 L=1,4
35 V(L,NDIV+J)=ST1*W(L,I) + ST2*W(L,I+1)
30 CONTINUE
C
C SET KINDEX(I+1) = THE INDEX OF THE NODE ON THE PIECEWISE-GEODESIC
C PATH WHICH CORRESPONDS IN TIME TO W(I+1).
C
KINDEX(I+1)=NDIV + NSPLIT
C
C UPDATE NDIV.
C
NDIV=NDIV + NSPLIT
40 CONTINUE
C
C PRINT NDIV.
C
PRINT *,' '
PRINT *,'NO. OF NODES IN PIECEWISE-GEODESIC REPRESENTATION OF'
PRINT *,'SMOOTH PATH ON S**3 (NDIV) = ',NDIV
C
IF (IPRINT.EQ.120) THEN
PRINT *,' '
PRINT *,'TIMES AND NODES ON PIECEWISE-GEODESIC PATH'
PRINT *,' (STARTING PATH)'
PRINT 910
DO 320 K=1,NDIV
320 PRINT 911,K,T(K),(V(L,K),L=1,4)
910 FORMAT(' ',2X,'K',6X,'T(K)',9X,'V(1,K)',7X,'V(2,K)',7X,
&'V(3,K)',7X,'V(4,K)')
911 FORMAT(' ',I3,2X,G12.5,1X,4(1X,G12.5))
PRINT *,' '
PRINT *,'THE TABLE KINDEX:'
PRINT *,' (KINDEX(I) = INDEX OF THE NODE WHICH CORRESP. IN'
PRINT *,' TIME TO I-TH DATA POINT.)'
PRINT 912
DO 325 I=1,NDATA
325 PRINT 913,I,KINDEX(I)
912 FORMAT(' ',4X,'I',2X,'KINDEX(I)')
913 FORMAT(' ',2X,I3,5X,I4)
END IF
C
C MAIN LOOP. THIS LOOP IS ON ITER = ITERATION NUMBER.
C
DO 210 ITER=1,ITERMX
C
PRINT *,' '
PRINT *,'START OF ITERATION NO. ',ITER
C
C UNROLL THE POINTS V(K), K = 1, ... , NDIV, ONTO R**3, GETTING
C IMAGE POINTS VSTAR(K). ALSO, UNWRAP THE DATA POINTS W(I), I = 1, ... ,
C NDATA, RELATIVE TO THE BASE PATH REPRESENTED BY THE V(K)'S, TO
C GET IMAGE POINTS WSTAR(I) IN R**3. IN ADDITION, UNWRAP THE EIGENVECTORS
C OF THE ERROR MATRICES CORRESP. TO THE DATA POINTS.
C
PRINT *,' UNWRAPPING TO R**3'
C
CALL UNRLLC(NDIV,V,NDATA,W,EQ,KINDEX,ITER,IPRINT,VSTAR,
&WSTAR,EQSTAR)
C
C CALCULATE THE UNWRAPPED ERROR MATRICES SIGINV(IDATA) CORRESP. TO THE
C UNWRAPPED DATA POINTS WSTAR(IDATA).
C
DO 60 I=1,NDATA
DO 55 J=1,3
DO 55 L=1,3
SIGINV(J,L,I)=0.0D0
DO 50 M=1,3
TERM=ELAMB(M,I)*EQSTAR(J+1,M,I)*EQSTAR(L+1,M,I)
50 SIGINV(J,L,I)=SIGINV(J,L,I) + TERM
55 CONTINUE
60 CONTINUE
C
IF ((IPRINT.EQ.145).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'THE UNWRAPPED ERROR MATRICES SIGINV(IDATA) CORRESP. TO'
PRINT *,'THE UNWRAPPED DATA POINTS WSTAR(IDATA)'
DO 360 I=1,NDATA
PRINT *,'DATA POINT NO. ',I
DO 365 J=1,3
365 PRINT 940,(SIGINV(J,L,I),L=1,3)
360 CONTINUE
940 FORMAT(' ',3(2X,G12.5))
END IF
C
C CALCULATE THE ARRAY GHAT, WHICH CONTAINS THE COEFFICIENTS (WITH
C RESPECT TO THE B-SPLINE BASIS) OF THE OPTIMAL SMOOTHING SPLINE TO
C THE POINTS WSTAR(I). (NOTE: IF ISMTH = 1, 2, OR 3, THE SUBROUTINE
C ESMTH ALSO FINDS THE OPTIMAL SMOOTHING PARAMETER ALAMBDA.)
C
IF (ISMTH.EQ.0) THEN
ILAMB(ITER)=ISLAMB
ALAMB(ITER)=ASLAMB
END IF
C
PRINT *,' FITTING IN R**3'
C
CALL ESMTH(NDATA,TAU,WSTAR,SIGINV,ISMTH,ILAMB(ITER),
&ALAMB(ITER),IPRINT,ITER,XVFVAL(ITER),GHAT,WFIT,FOBJ(ITER))
C
C CALCULATE THE FITTED POINTS QSTAR(K) IN R**3, BASED ON THE OPTIMAL
C COEFFICIENTS GHAT(L,M). (IF GSTAR(T) IS THE FITTED SPLINE IN R**3,
C THEN QSTAR(K) = GSTAR(T(K)).)
C
CALL CETAU(NDATA,TAU,ETAU)
C
C THE CASE I = 1.
C
I=1
KEND=KINDEX(2) - 1
DO 75 K=1,KEND
DO 70 L=1,3
70 QSTAR(L,K)=0.0D0
DO 75 M=1,3
CALL CBETA(NDATA,ETAU,M,T(K),I,BETAV)
DO 72 L=1,3
72 QSTAR(L,K)=QSTAR(L,K) + GHAT(L,M)*BETAV
75 CONTINUE
C
C THE CASES I = 2, 3, ... , NDATA - 2.
C
DO 90 I=2,NDATA-2
KSTART=KINDEX(I)
KEND=KINDEX(I+1) - 1
DO 85 K=KSTART,KEND
DO 80 L=1,3
80 QSTAR(L,K)=0.0D0
DO 85 M=I-1,I+2
CALL CBETA(NDATA,ETAU,M,T(K),I,BETAV)
DO 82 L=1,3
82 QSTAR(L,K)=QSTAR(L,K) + GHAT(L,M)*BETAV
85 CONTINUE
90 CONTINUE
C
C THE CASE I = NDATA - 1.
C
I=NDATA - 1
KSTART=KINDEX(NDATA-1)
KEND=KINDEX(NDATA)
DO 105 K=KSTART,KEND
DO 100 L=1,3
100 QSTAR(L,K)=0.0D0
DO 105 M=NDATA-2,NDATA
CALL CBETA(NDATA,ETAU,M,T(K),I,BETAV)
DO 102 L=1,3
102 QSTAR(L,K)=QSTAR(L,K) + GHAT(L,M)*BETAV
105 CONTINUE
C
IF ((IPRINT.EQ.160).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'NODES OF FITTED SPLINE IN R**3'
PRINT 950
DO 370 K=1,NDIV
370 PRINT 951,K,(QSTAR(J,K),J=1,3)
950 FORMAT(' ',2X,'K',10X,'QSTAR(K)')
951 FORMAT(' ',I3,3(2X,G12.5))
END IF
C
C WRAP THE FITTED POINTS QSTAR(K) ONTO S**3, RELATIVE TO THE BASE PATH
C REPRESENTED BY THE V(K)'S. THE IMAGE POINTS OBTAINED ARE: Q(1),
C Q(2), ... , Q(NDIV).
C
PRINT *,' WRAPPING BACK TO S**3'
C
CALL WRAPC(NDIV,V,QSTAR,VSTAR,ITER,IPRINT,Q)
C
C IF ITER.EQ.1, WE SAVE THE FITTED SPLINE PATH FOR COMPARISON WITH
C FUTURE ITERATIONS.
C
IF (ITER.EQ.1) THEN
DO 172 K=1,NDIV
DO 172 J=1,4
172 QFIRST(J,K)=Q(J,K)
END IF
C
C CALCULATE ACC(ITER) = MAX. DIFFERENCE BETWEEN SUCCESSIVE FITTED
C SPLINES (MEASURED IN DEGREES ON S**3).
C
PMIN=1.0D0
DO 175 K=1,NDIV
PROD=0.0D0
DO 180 J=1,4
180 PROD=PROD + Q(J,K)*V(J,K)
PMIN=DMIN1(PMIN,PROD)
175 CONTINUE
ACC(ITER)=DACOS(PMIN)*(180.0D0/PI)
C
C CALCULATE THE DISTW(I)'S, I.E., THE DISTANCES (IN DEGREES ON S**3)
C BETWEEN THE DATA POINTS AND THE FITTED SPLINE. ALSO, CALCULATE
C ACCD(ITER) = THE MAX. OF THESE DISTANCES.
C
DO 185 I=1,NDATA
PROD=0.0D0
K=KINDEX(I)
DO 190 J=1,4
190 PROD=PROD + W(J,I)*Q(J,K)
PROD=DMIN1(1.0D0,PROD)
PROD=DMAX1(-1.0D0,PROD)
DISTW(I)=DACOS(PROD)*(180.0D0/PI)
185 CONTINUE
ACCD(ITER)=DISTW(1)
DO 195 I=2,NDATA
195 ACCD(ITER)=DMAX1(DISTW(I),ACCD(ITER))
C
C TEST FOR CONVERGENCE.
C IF CONVERGENCE HAS BEEN ACHIEVED, CALCULATE MAX. DIFFERENCE BETWEEN
C FIRST AND FINAL FITTED SPLINES.
C IF NOT, PREPARE FOR NEXT ITERATION.
C
IF (ACC(ITER).LE.ANGACC) THEN
NITER=ITER
PRINT *,'CONVERGENCE AT ITERATION NO. ',NITER
GO TO 212
ELSE
DO 200 K=1,NDIV
DO 205 J=1,4
205 V(J,K)=Q(J,K)
200 CONTINUE
END IF
C
210 CONTINUE
C
C END OF MAIN LOOP.
C
NITER=ITERMX
PRINT *,'NO CONVERGENCE AFTER MAX. NO. OF ITERATIONS'
C
C CALCULATE ACC1 = MAX. DIFFERENCE BETWEEN FIRST AND FINAL FITTED
C SPLINES (MEASURED IN DEGREES ON S**3).
C
212 CONTINUE
PMIN=1.0D0
DO 215 K=1,NDIV
PROD=0.0D0
DO 218 J=1,4
218 PROD=PROD + QFIRST(J,K)*Q(J,K)
PMIN=DMIN1(PMIN,PROD)
215 CONTINUE
ACC1=DACOS(PMIN)*180.0D0/PI
C
C IF ICF=1, WE CALCULATE THE FITTED POINTS GFHAT(TSTAR(NU)).
C
C FIRST, CALCULATE THE ARRAY KSTAR(NUTOT). THE VALUE K = KSTAR(NU) IS
C SUCH THAT TSTAR(NU) LIES IN THE INTERVAL [T(K),T(K+1)). THE CASE
C TSTAR(NU).EQ.T(NDIV) IS A SPECIAL CASE; IN THIS CASE WE SET
C KSTAR(NU) = NDIV.
C
C SECOND, CALCULATE THE ARRAY GFHAT(4,NUTOT).
C
IF (ICF.EQ.1) THEN
DO 220 NU=1,NUMAX
IF (NU.EQ.1) THEN
KSTART=1
ELSE
KSTART=KSTAR(NU-1)
END IF
DO 225 K=KSTART,NDIV
IF (T(K).GT.TSTAR(NU)) THEN
KSTAR(NU)=K-1
GO TO 223
END IF
225 CONTINUE
KSTAR(NU)=NDIV
223 CONTINUE
220 CONTINUE
C
DO 230 NU=1,NUMAX
IF (KSTAR(NU).LT.NDIV) THEN
K=KSTAR(NU)
PROD=0.0D0
DO 232 J=1,4
232 PROD=PROD + Q(J,K)*Q(J,K+1)
DIST=DACOS(PROD)
SD=DSIN(DIST)
F1=(T(K+1) - TSTAR(NU))/(T(K+1) - T(K))
F2=(TSTAR(NU) - T(K))/(T(K+1) - T(K))
ST1=DSIN(F1*DIST)/SD
ST2=DSIN(F2*DIST)/SD
DO 235 J=1,4
235 GFHAT(J,NU)=ST1*Q(J,K) + ST2*Q(J,K+1)
ELSE
DO 238 J=1,4
238 GFHAT(J,NU)=Q(J,NDIV)
END IF
230 CONTINUE
END IF
C
IF (IPRINT.EQ.280) THEN
IF (ICF.EQ.1) THEN
PRINT *,' '
PRINT *,'THE ARRAY KSTAR(NUTOT):'
PRINT 955
DO 375 NU=1,NUMAX
375 PRINT 956,NU,KSTAR(NU)
955 FORMAT(' ',1X,'NU',3X,'KSTAR(NU)')
956 FORMAT(' ',I3,5X,I5)
PRINT *,' '
PRINT *,'THE ARRAY GFHAT(4,NUTOT):'
PRINT 960
DO 380 NU=1,NUMAX
380 PRINT 961,NU,(GFHAT(J,NU),J=1,4)
960 FORMAT(' ',1X,'NU',10X,'GFHAT(TSTAR(NU))')
961 FORMAT(' ',I3,4(2X,G12.5))
ELSE
PRINT *,' '
PRINT *,'SINCE NO TIMES TSTAR(NU) WERE SPECIFIED, THERE ARE'
PRINT *,'NO ARRAYS KSTAR(NUTOT) AND GFHAT(4,NUTOT).'
END IF
END IF
C
IF (ICF.EQ.1) THEN
WRITE(30,*) ILAMB(NITER)
IF (ILAMB(NITER).EQ.1) THEN
WRITE(30,*) ALAMB(NITER)
END IF
WRITE(30,*) SLONG,IGPH
IF (IGPH.EQ.1) THEN
WRITE(30,999) FNAME
999 FORMAT (A)
END IF
WRITE(30,*) NDATA
DO 285 I=1,NDATA
WRITE(30,*) TAU(I),KINDEX(I)
WRITE(30,*) (WSTAR(L,I),L=1,3)
DO 290 J=1,3
290 WRITE(30,*) (SIGINV(J,L,I),L=1,3)
WRITE(30,*) (GHAT(L,I),L=1,3)
K=KINDEX(I)
WRITE(30,*) (QSTAR(L,K),L=1,3)
285 CONTINUE
WRITE(30,*) NDIV
DO 293 K=1,NDIV
293 WRITE(30,*) T(K)
WRITE(30,*) NUMAX
DO 294 NU=1,NUMAX
WRITE(30,*) TSTAR(NU),KSTAR(NU)
294 WRITE(30,*) (GFHAT(J,NU),J=1,4)
END IF
C
C SET UP THE ARRAY QFIT, SO THAT IT CONTAINS THE FITTED SPLINE POINTS
C CORRESPONDING TO THE DATA TIMES.
C
DO 295 I=1,NDATA
K=KINDEX(I)
DO 298 J=1,4
298 QFIT(J,I)=Q(J,K)
295 CONTINUE
C
RETURN
END
C
C *** UNRLLC ***
C
C THIS SUBROUTINE UNROLLS THE POINTS V(1), ... , V(NDIV) ONTO R**3,
C GETTING IMAGE POINTS VSTAR(1), ... , VSTAR(NDIV). ALSO, IN THE PROCESS,
C THIS SUBROUTINE UNWRAPS THE DATA POINTS W(1), ... , W(NDATA), RELATIVE
C TO THE BASE PATH REPRESENTED BY THE V(K)'S, TO GET IMAGE POINTS
C WSTAR(1), ... , WSTAR(NDATA) IN R**3. FINALLY, THIS SUBROUTINE UNWRAPS
C THE EIGENVECTORS EQ(L,IDATA) OF THE ERROR MATRICES ASSOC. WITH THE DATA
C POINTS, RELATIVE TO THE BASE PATH REPRESENTED BY THE V(K)'S.
C
C NOTE: UNROLLING IS ACCOMPLISHED BY LEFT-TRANSLATION. CORRESPONDING
C TO A PATH F(T) ON S**3, THE ROTATION PATH R(T) IN SO(4) IS
C DEFINED BY:
C
C R(T)*X = CONJ(F(T))*X .
C
C HERE X IS A VECTOR IN R**4; THE MULTIPLICATION ON THE LEFT IS
C MATRIX MULTIPLICATION, WHILE THE MULTIPLICATION ON THE RIGHT IS
C QUATERNION MULTIPLICATION. THE BASE PATH F(T) IS REPRESENTED IN
C DISCRETE FORM BY THE V(K)'S, I.E., V(K) = F(T(K)).
C
C INPUTS:
C NDIV - NO. OF NODES IN PIECEWISE-GEODESIC REPRESENTATION OF
C BASE PATH ON S**3
C V(4,NDIVT) - NODES OF BASE PATH ON S**3
C NDATA - NO. OF DATA POINTS
C W(4,NDATAT) - DATA POINTS ON S**3
C EQ(4,3,NDATAT) - EIGENVECTORS OF THE ERROR MATRICES ASSOC. WITH THE
C DATA POINTS. EQ(J,L,I) IS THE J-TH COMPONENT OF
C THE L-TH EIGENVECTOR OF THE ERROR MATRIX CORRESP.
C TO W(I).
C KINDEX(NDATAT) - INTEGER ARRAY: FOR EACH I, KINDEX(I) IS THE INDEX
C OF THE NODE ON THE BASE PATH WHICH CORRESP. IN TIME
C TO THE I-TH DATA POINT
C ITER - ITERATION NO.
C IPRINT - CODE FOR PRINTING OPTION
C 140 - THE ARRAYS VSTAR, WSTAR, EQSTAR
C
C OUTPUTS:
C VSTAR(3,NDIVT) - NODES OF UNROLLED BASE PATH IN R**3
C WSTAR(3,NDATAT)- UNWRAPPED DATA POINTS IN R**3
C EQSTAR(4,3,NDATAT) - UNWRAPPED EIGENVECTORS FOR THE ERROR MATRICES
C
SUBROUTINE UNRLLC(NDIV,V,NDATA,W,EQ,KINDEX,ITER,IPRINT,
&VSTAR,WSTAR,EQSTAR)
PARAMETER (NDATAT=50)
PARAMETER (NDIVT=40000)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION V(4,NDIVT),W(4,NDATAT),EQ(4,3,NDATAT)
DIMENSION KINDEX(NDATAT)
DIMENSION VSTAR(3,NDIVT),WSTAR(3,NDATAT),EQSTAR(4,3,NDATAT)
DIMENSION CV(4),V1(4),EX(4),X(3)
DIMENSION VZ(4),WZ(4),EH(4),H(3)
DIMENSION CW(4),EQZ(4),EQSTRZ(4)
C
C PREPARE FOR LOOP. SET VSTAR(1) = 0.
C
DO 10 J=1,3
10 VSTAR(J,1)=0.0D0
IDATA=1
C
C MAIN LOOP. THIS LOOP IS ON THE NODES IN THE BASE PATH.
C
DO 250 K=1,NDIV
C
C IF K.LT.NDIV, CALCULATE VSTAR(K+1):
C
C VSTAR(K+1) = VSTAR(K) + X(K) ,
C
C WHERE: EXP(X(K)) = CONJ(V(K))*V(K+1) .
C
CV(1)=V(1,K)
DO 20 J=2,4
20 CV(J)= -V(J,K)
C
IF (K.LT.NDIV) THEN
DO 25 J=1,4
25 V1(J)=V(J,K+1)
CALL QMULT(CV,V1,EX)
CALL IVEXP2(EX,X,IERR)
IF (IERR.EQ.1) THEN
PRINT *,' '
PRINT *,'ERROR IN USE OF INVEXP, IN UNRLLC, TO CALCU-'
PRINT *,'LATE X.'
PRINT *,'K = ',K
STOP
END IF
DO 30 J=1,3
30 VSTAR(J,K+1)=VSTAR(J,K) + X(J)
END IF
C
C UNWRAP THE DATA POINT CORRESP. TO THE TIME T(K), IF THERE IS ONE.
C
120 CONTINUE
IF (IDATA.GT.NDATA) GO TO 240
IF (KINDEX(IDATA).NE.K) GO TO 240
C
C AT THIS JUNCTURE, KINDEX(IDATA) = K. HENCE, THE IDATA-TH DATA POINT
C CORRESPONDS IN TIME TO THE K-TH NODE ON THE BASE PATH. THE UNWRAPPED
C DATA POINT WSTAR(IDATA) IS GIVEN BY:
C
C WSTAR(IDATA) = VSTAR(K) + H(IDATA) ,
C
C WHERE: EXP(H(IDATA)) = CONJ(V(K))*W(IDATA) .
C
C THE UNWRAPPED EIGENVECTORS EQSTAR(L,IDATA) ARE GIVEN BY:
C
C EQSTAR(L,IDATA) = CONJ(W(IDATA))*EQ(L,IDATA) , L = 1, 2, 3 .
C
DO 125 J=1,4
VZ(J)=V(J,K)
125 WZ(J)=W(J,IDATA)
CALL QMULT(CV,WZ,EH)
CALL IVEXP2(EH,H,IERR)
IF (IERR.EQ.1) THEN
PRINT *,' '
PRINT *,'ERROR IN USE OF INVEXP, IN UNRLLC, TO CALCULATE H.'
PRINT *,'K = ',K
STOP
END IF
DO 130 J=1,3
130 WSTAR(J,IDATA)=VSTAR(J,K) + H(J)
C
C CALCULATE EQSTAR(L,IDATA), L = 1, 2, 3.
C
CW(1)=W(1,IDATA)
DO 135 J=2,4
135 CW(J)= -W(J,IDATA)
DO 150 L=1,3
DO 140 J=1,4
140 EQZ(J)=EQ(J,L,IDATA)
CALL QMULT(CW,EQZ,EQSTRZ)
DO 145 J=1,4
145 EQSTAR(J,L,IDATA)=EQSTRZ(J)
150 CONTINUE
C
IDATA=IDATA + 1
GO TO 120
C
240 CONTINUE
C
250 CONTINUE
C
C END OF MAIN LOOP.
C
IF ((IPRINT.EQ.140).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'THE ARRAY VSTAR (UNROLLED BASE PATH)'
PRINT 900
DO 300 K=1,NDIV
300 PRINT 901,K,(VSTAR(J,K),J=1,3)
900 FORMAT(' ',2X,'K',10X,'VSTAR(K)')
901 FORMAT(' ',I3,4(2X,G12.5))
PRINT *,' '
PRINT *,'THE ARRAY WSTAR (UNWRAPPED DATA POINTS)'
PRINT 902
DO 305 I=1,NDATA
305 PRINT 901,I,(WSTAR(J,I),J=1,3)
902 FORMAT(' ',2X,'I',10X,'WSTAR(I)')
PRINT *,' '
PRINT *,'THE UNWRAPPED EIGENVECTORS OF THE ERROR MATRICES'
PRINT *,'CORRESP. TO THE DATA POINTS'
DO 400 I=1,NDATA
PRINT *,'DATA POINT NO. ',I
PRINT *,' EIGENVECTORS CORRESP. TO POS. EIGENVALUES'
PRINT *,' (COLUMNS BELOW)'
DO 410 J=1,4
410 PRINT 930,(EQSTAR(J,L,I),L=1,3)
930 FORMAT(' ',3(2X,G12.5))
400 CONTINUE
END IF
C
RETURN
END
C
C *** WRAPC ***
C
C THIS SUBROUTINE WRAPS THE FITTED POINTS QSTAR(1), QSTAR(2), ... ,
C QSTAR(NDIV) ONTO S**3, RELATIVE TO THE BASE PATH REPRESENTED BY THE
C V(K)'S. THE IMAGE POINTS OBTAINED ARE: Q(1), Q(2), ... , Q(NDIV).
C
C NOTE: IF F(T) IS THE BASE PATH ON S**3, THE ROTATION PATH R(T) IN
C SO(4) IS DEFINED BY:
C
C R(T)*X = CONJ(F(T))*X .
C
C HERE X IS A VECTOR IN R**4; THE MULTIPLICATION ON THE LEFT IS
C MATRIX MULTIPLICATION, WHILE THE MULTIPLICATION ON THE RIGHT IS
C QUATERNION MULTIPLICATION.
C
C INPUTS:
C NDIV - NO. OF NODES IN DISCRETE REPRESENTATION OF BASE
C PATH ON S**3
C V(4,NDIVT) - NODES OF BASE PATH ON S**3
C QSTAR(3,NDIVT)- FITTED POINTS IN R**3
C VSTAR(3,NDIVT)- NODES OF UNROLLED BASE PATH IN R**3
C ITER - ITERATION NO.
C IPRINT - CODE FOR PRINTING OPTION
C 170 - THE ARRAY Q (FITTED POINTS WRAPPED ONTO
C S**3)
C
C OUTPUTS:
C Q(4,NDIVT) - IMAGES OF THE FITTED POINTS QSTAR(K), ON S**3
C
SUBROUTINE WRAPC(NDIV,V,QSTAR,VSTAR,ITER,IPRINT,Q)
PARAMETER (NDIVT=40000)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION V(4,NDIVT),QSTAR(3,NDIVT),VSTAR(3,NDIVT)
DIMENSION Q(4,NDIVT)
DIMENSION DIFF(3),U(0:3),VZ(4),QZ(4)
C
C MAIN LOOP. THIS LOOP IS ON THE NODES IN THE BASE PATH.
C
DO 200 K=1,NDIV
C
C CALCULATE DIFF = QSTAR(K) - VSTAR(K) .
C
DO 25 J=1,3
25 DIFF(J)=QSTAR(J,K) - VSTAR(J,K)
C
C CALCULATE Q(K) = V(K)*(EXP(DIFF)) .
C
CALL CLCEXP(DIFF,U)
DO 35 J=1,4
35 VZ(J)=V(J,K)
CALL QMULT(VZ,U,QZ)
DO 45 J=1,4
45 Q(J,K)=QZ(J)
C
200 CONTINUE
C
C END OF MAIN LOOP.
C
IF ((IPRINT.EQ.170).AND.(ITER.LE.5)) THEN
PRINT *,' '
PRINT *,'ITER = ',ITER
PRINT *,'THE ARRAY Q (FITTED POINTS WRAPPED ONTO S**3)'
PRINT 900
DO 300 K=1,NDIV
300 PRINT 901,K,(Q(J,K),J=1,4)
900 FORMAT(' ',2X,'K',10X,'Q(K)')
901 FORMAT(' ',I3,4(2X,G12.5))
END IF
C
RETURN
END