@PROCESS DC(RNBLK, @PROCESS RUBLK, @PROCESS AXBLK, @PROCESS RN1BLK) C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: TDAP - THREE DIMENSIONAL ADJUSTMENT PROGRAM * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: MAIN * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS PROGRAM IS TO PERFORM A THREE * C * DIMENSIONAL PRE-ANALYSIS OR ADJUSTMENT OF A GEODETIC NETWORK. * C * * C * * C * * C * SUBROUTINES REQUIRED: CONF3D, CPUTIM (UNBCC), GDATE(UNBCC), * C * NORVEC, READOB, SPIN (UNBCC-SELIB). * C * * C * FUNCTION SUBPROGRAMS REQUIRED: N/A * C * * C * ALGORITHM/METHOD: N/A * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C * * C * COPYRIGHT (C) 1988 MARK WILLIAM ROHDE * C * * C * >>>>>>>> WARNING <<<<<<<< * C * * C * ALL RIGHTS RESERVED. NO PART OF THIS WORK MAY BE REPRODUCED IN ANY * C * FORM OR BY ANY MEANS OR USED TO MAKE A DERIVATIVE WORK (SUCH AS A * C * TRANSLATION, TRANSFORMATION, OR ADAPTATION) WITHOUT PERMISSION IN * C * WRITING FROM THE AUTHOR. * C * * C ********************************************************************** C * * C * VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * IND I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (DEGREES) * C * INM I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (MINUTES) * C * INS R*8 ARRAY OF ORTHOGONAL ZENITH ANGLES (SECONDS) * C * ISET I*2 CURRENT OBSERVATION SET NUMBER BEING PROCESSED * C * * C * * C * * C * I0 R*8 INCLINATION COMPONENT IN THE DIRECTION 0 DEGREES * C * I1 R*8 INCLINATION COMPONENT IN THE DIRECTION 90 DEGREES * C * * C * * C * * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FILE(2),FIXSTA(100),FRESTA(500),STA(3500,3) CHARACTER*40 PROJ(3) CHARACTER*78 TITLE INTEGER*2 APOPT,NIO(5) INTEGER*4 OBCODE(3500) LOGICAL FAIL,FINAL,REFVAR C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) COMMON /RUBLK/ RU(1000),IPVT(1000) COMMON /AXBLK/ AUX(40000) COMMON /RN1BLK/ RN1(1000,1000) DIMENSION A(3500,6),FIXCOR(100,3),FRECOR(500,3),ICA(3500,6) DIMENSION OBS(3500,4),W(3500),DELTA(1000) DIMENSION RHAT(3500),RHATS(3500),RSTD(3500) C C INITIALIZE DATA CONSTANTS DATA (NIO(I), I = 1, 5) /0, 0, 0, 0, 0/ C C ASSIGN MATRIX SIZE CONTROL VARIABLES NOB = 3500 NFI = 100 NFR = 500 NUK = 1000 C C SET ELAPSED CPU TIME COUNTER TO ZERO CALL CLOCKX(ELAPS1) C C INITIALIZE PROGRAM CONTROL PARAMETERS FAIL = .FALSE. FINAL = .FALSE. REFVAR = .FALSE. ITER = -1 C C READ IN OBSERVATION DATA INPUT FILE CALL READOB(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,NFIX,NFRE,NOBS,OBCOD &E,OBS,STA,NFI,NFR,NOB,TITLE,FILE,PROJ,NIO) NUKN = NFRE * 3 C C PRINT OUT TITLE PAGE AND INPUT DATA CALL PRTDAT(APOPT,FILE,FIXCOR,FIXSTA,FRECOR,FRESTA,PROJ,TITLE,NFI, &NFIX,NFR,NFRE,FINAL) C C COMPUTE THE FIRST DESIGN MATRIX (A) AND THE WEIGHT MATRIX (P) 5000 ITER = ITER + 1 CALL NORVEC(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,ICA,NFIX,NFRE,NIDB,N &OBS,NUKN,OBCODE,OBS,STA,A,W,NFI,NFR,NOB,NUK,ITER) DO 2 I = 1,NUKN DO 3 J = 1,NUKN RN1(I,J) = RN(I,J) 3 CONTINUE 2 CONTINUE C C SOLVE SYSTEM OF NORMAL EQUATIONS IF (APOPT .EQ. 0) THEN CALL DGEF(RN1,NUK,NUKN,IPVT) CALL DGES(RN1,NUK,NUKN,IPVT,RU) ELSE GOTO 5001 ENDIF C C INVERSE THE COEFFICIENT MATRIX OF THE NORMAL EQUATIONS C5001 CALL CLOCKX(TIME1) C CALL SPIN(NUKN,1000,DETRN,IEXPRN) C CALL CLOCKX(TIME2) C TIME = (TIME2 - TIME1) / 1.0E06 C WRITE(6,2000) TIME C2000 FORMAT('0','MATRIX INVERSION ELAPSED TIME = ',F12.2,' SECONDS') C C SOLVE SYSTEM OF NORMAL EQUATIONS C CALL MMULD(DELTA,NUK,RN,NUK,RU,NUK,NUKN,NUKN,1) C DO 201 I = 1,NUKN C RU(I) = DELTA(I) C 201 CONTINUE C C UPDATE APPROXIMATE COORDINATES CALL UPDATE(FRECOR,NFRE,ITER,NFR,NUK,0,FRESTA) C C CHECK FOR SOLUTION DIVERGENCE OR CONVERGENCE C CALL CHKDIV( IF (ITER .EQ. 10) THEN WRITE(6,*) ITER GOTO 6000 ENDIF DO 1 I = 1,NUKN IF (DABS(RU(I)) .GE. 1.D-06) GOTO 5000 1 CONTINUE C C INVERSE THE COEFFICIENT MATRIX OF THE NORMAL EQUATIONS 5001 NAUX = 40000 CALL CLOCKX(TIME1) CALL DGEICD(RN,NUK,NUKN,1,RCOND,DET,AUX,NAUX) CALL CLOCKX(TIME2) TIME = (TIME2 - TIME1) / 1.0E06 WRITE(6,2000) TIME 2000 FORMAT('0'//1X,'MATRIX INVERSION ELAPSED TIME = ',F9.5,' SECONDS') RCOND = 1.D0 / RCOND WRITE(6,1004) RCOND IDF = NOBS - NUKN - NIDB IF (APOPT .EQ. 1) GOTO 5003 C C COMPUTE ADJUSTED RESIDUALS CALL RESID(A,ICA,NOBS,OBCODE,OBS,RHAT,W,NOB,NUK,IDF,APOST,RSTD,RHA &TS,STA) C C PRINT STATISTICAL SUMMARY ON RESULTS FROM ADJUSTMENT CALL STATSUM(APOST,FAC,FAIL,IDF,ITER,NIO,NIDB,NOBS,NUKN) C C DETERMINE IF VARIANCE-COVARIANCE MATRIX OF THE PARAMETERS MUST BE C SCALED BY THE A POSTERIORI VARIANCE FACTOR IF (FAIL) CALL SCALAR(NUKN,NUKN,APOST) C C PRINT OUT ADJUSTED COORDINATES AND COVARIANCE MATRIX TO AUXILIARY FILE FINAL = .TRUE. C WRITE(14,*) NFIX, NFRE, NUKN CALL PRTDAT(APOPT,FILE,FIXCOR,FIXSTA,FRECOR,FRESTA,PROJ,TITLE,NFI, &NFIX,NFR,NFRE,FINAL) C C PRINT FINAL ADJUSTED COORDINATES WRITE(14,1010) IDF WRITE(14,1011) APOST WRITE(14,1007) DO 100 I = 1,NFRE WRITE (14,1005) FRESTA(I),FRECOR(I,1),FRECOR(I,2),FRECOR(I,3) 100 CONTINUE WRITE(14,1008) DO 101 I = 1,NFIX WRITE (14,1005) FIXSTA(I),FIXCOR(I,1),FIXCOR(I,2),FIXCOR(I,3) 101 CONTINUE C C PRINT OUT VARIANCE-COVARIANCE MATRIX AS A VECTOR, COLUMN-BY-COLUMN WRITE (14,1009) DO 102 J = 1,NUKN DO 103 I = 1,J WRITE(14,1006) RN(I,J) 103 CONTINUE 102 CONTINUE C C COMPUTE RELATIVE VARIANCES AND COVARIANCES FOR SELECTED ELEMENTS C CALL DIPEL(FRESTA,NFR,NFRE) C C CHECK FOR RESIDUAL OUTLIERS IF (.NOT. FAIL) REFVAR = .TRUE. CALL OUTLIER(RHATS,STA,5.D0,1.D0,IDF,REFVAR,NOB,NOBS,OBCODE,OBS,RH &AT) C C COMPUTE AND PRINT ABSOLUTE 95% CONFIDENCE ELLIPSOIDS 5003 CALL CONF3D(0.05D0,IDF,FAIL,FRESTA,NUKN,NFR,NUK,APOST) C C COMPUTE AND PRINT ABSOLUTE 95% CONFIDENCE ELLIPSES (CROSS-SECTIONS OF C 95% CONFIDENCE ELLIPSOIDS) CALL CONF2D(0.05D0,IDF,FAIL,FRESTA,NUKN,NFR,NUK,APOST) C C DETERMINE PROGRAM EXECUTION TIME 6000 CALL CLOCKX(ELAPS2) ELAPS = (ELAPS2 - ELAPS1) / 1.0E06 WRITE(6,1003) ELAPS C C I/O FORMAT STATEMENTS 1000 FORMAT(A78) 1001 FORMAT(I1,2(2X,I3)) 1002 FORMAT('1',80('*')/1X,'*',78X,'*'/1X,'*',A78,'*'/1X,'*',78X,'*'/1X &,80('*')/1X,'DATE: ',18A1,42X,'TIME: ',2A1,':',2A1,':',2A1///) 1003 FORMAT('0','PROGRAM EXECUTION TIME = ',F12.2,' (SEC)'//) 1004 FORMAT('0','CONDITION NUMBER FOR THE SYSTEM OF NORMAL EQUATIONS = &',F20.10//) 1005 FORMAT(1X,A8,2X,3(E22.16,1X)) 1006 FORMAT(1X,E22.16) 1007 FORMAT(' ','FREE STATIONS') 1008 FORMAT(' ','FIXED STATIONS') 1009 FORMAT(' ','VARIANCE-COVARIANCE MATRIX') 1010 FORMAT(' ','NUMBER OF DEGREES OF FREEDOM = ',I5) 1011 FORMAT(' ','A POSTERIORI VARIANCE FACTOR = ',E22.16) C STOP END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: READOB * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO READ IN THE * C * OBSERVATION DATA INPUT FILE. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL READOB(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,NFIX,NFRE, * C * NOBS,OBCODE,OBS,STA) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * APOP I*2 ADJUSTMENT/PRE-ANALYSIS OPTION CODE * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * OBS R*8 MATRIX CONTAINING OBSERVATIONS * C * STA C*8 MATRIX OF AT/FROM/TO STATION NAMES * C * NOBS I*2 NUMBER OF OBSERVATIONS READ IN * C * OBCODE I*2 ARRAY OF OBSERVATION CODE NUMBERS * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C SUBROUTINE READOB(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,NFIX,NFRE,NOBS &,OBCODE,OBS,STA,NFI,NFR,NOB,TITLE,FILE,PROJ,NIO) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FILE(2),FIXSTA(NFI),FRESTA(NFR),STA(NOB,3),TEMP CHARACTER*40 PROJ(3) CHARACTER*78 TITLE INTEGER*2 APOPT,NIO(5) INTEGER*4 OBCODE(NOB) C C ASSIGN MATRIX DIMENSIONS DIMENSION FIXCOR(NFI,3),FRECOR(NFR,3),OBS(NOB,4) C C INITIALIZE OBSERVATION COUNTERS NOBS = 0 C C READ IN INPUT/OUTPUT FILE NAMES READ (5,1000) FILE(1),FILE(2) C C READ IN PROJECT DESCRIPTORS READ (5,1001) PROJ(1) READ (5,1001) PROJ(2) READ (5,1001) PROJ(3) C C READ IN TITLE READ (5,1002) TITLE C C READ IN PROGRAM CONTROL PARAMETERS READ (5,1003) APOPT,NFIX,NFRE C C READ IN FIXED STATION NAMES AND COORDINATES DO 5 I = 1,NFIX READ(5,1004) FIXSTA(I),FIXCOR(I,1),FIXCOR(I,2),FIXCOR(I,3) C WRITE(6,*) FIXSTA(I),FIXCOR(I,1),FIXCOR(I,2),FIXCOR(I,3) 5 CONTINUE C C READ IN FREE STATION NAMES AND APPROXIMATE COORDINATES DO 10 I = 1,NFRE READ(5,1004) FRESTA(I),FRECOR(I,1),FRECOR(I,2),FRECOR(I,3) C WRITE(6,*) FRESTA(I),FRECOR(I,1),FRECOR(I,2),FRECOR(I,3) 10 CONTINUE C C READ IN OBSERVATION INPUT DATA DO 15 I = 1,NOB READ(5,1005) OBCODE(I),STA(I,1),STA(I,2),STA(I,3),OBS(I,1),OBS( &I,2),OBS(I,3),OBS(I,4) IF (OBS(I,1) .LE. 0.D0) OBS(I,1) = 0.001D0 NOBS = NOBS + 1 IF (OBCODE(I) .EQ. 1) NIO(1) = NIO(1) + 1 IF (IABS(OBCODE(I)) .EQ. 2) NIO(2) = NIO(2) + 1 IF (OBCODE(I) .EQ. 3) NIO(3) = NIO(3) + 1 IF (OBCODE(I) .EQ. 4) NIO(4) = NIO(4) + 1 IF (OBCODE(I) .EQ. 5) NIO(5) = NIO(5) + 1 IF (OBCODE(I) .EQ. 1) OBS(I,1) = OBS(I,1) / 1000.D0 IF (OBCODE(I) .EQ. -9) GOTO 5000 C IF (OBCODE(I) .NE. 1 .AND. OBCODE(I) .NE. 4) OBS(I,1) = 0.5D0 15 CONTINUE C C I/O FORMAT STATEMENTS 1000 FORMAT(2(A8,2X)) 1001 FORMAT(A40) 1002 FORMAT(A78) 1003 FORMAT(I1,2X,I3,2X,I3) 1004 FORMAT(A8,2X,3(F16.4,4X)) 1005 FORMAT(I2,2X,3(A8,2X),F7.5,3(1X,F9.5)) C1006 FORMAT(I2,8X,3(A8,2X),F12.6,3(F3.1)) 5000 NOBS = NOBS - 1 C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: PRTDAT * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO READ IN THE * C * OBSERVATION DATA INPUT FILE. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL PRTDAT(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,NFIX,NFRE, * C * NOBS,OBCODE,OBS,STA) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * APOP I*2 ADJUSTMENT/PRE-ANALYSIS OPTION CODE * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * OBS R*8 MATRIX CONTAINING OBSERVATIONS * C * STA C*8 MATRIX OF AT/FROM/TO STATION NAMES * C * NOBS I*2 NUMBER OF OBSERVATIONS READ IN * C * OBCODE I*2 ARRAY OF OBSERVATION CODE NUMBERS * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK) SUBROUTINE PRTDAT(APOPT,FILE,FIXCOR,FIXSTA,FRECOR,FRESTA,PROJ,TITL &E,NFI,NFIX,NFR,NFRE,FINAL) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FILE(2),FIXSTA(NFI),FRESTA(NFR) CHARACTER*40 PROJ(3) CHARACTER*78 TITLE INTEGER*2 APOPT LOGICAL FINAL C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) DIMENSION FIXCOR(NFI,3),FRECOR(NFR,3),INOW(14) C C DETERMINE CURRENT DATE AND TIME IF (.NOT. FINAL) THEN CALL DATIMX(INOW) C C PRINT TITLE PAGE AND OPTION SUMMARY WRITE (6,1000) TITLE IF (INOW(1) .NE. -1) THEN WRITE (6,1013) INOW(8),INOW(7),INOW(6),INOW(5),INOW(4),INOW(3) ENDIF WRITE (6,1001) FILE(1),FILE(2) WRITE (6,1002) PROJ(1),PROJ(2),PROJ(3) C WRITE (6,1003) C C PRINT INITIAL APPROXIMATE COORDINATES WRITE (6,1004) WRITE (6,1005) WRITE (6,1006) DO 1 I = 1,NFRE WRITE (6,1007) FRESTA(I),FRECOR(I,1),FRECOR(I,2),FRECOR(I,3) 1 CONTINUE WRITE (6,1008) WRITE (6,1006) DO 2 I = 1,NFIX WRITE (6,1007) FIXSTA(I),FIXCOR(I,1),FIXCOR(I,2),FIXCOR(I,3) 2 CONTINUE C C PRINT INITIAL APPROXIMATE COORDINATES AND CORRESPONDING STANDARD C DEVIATIONS TO A FILE C WRITE(14,1004) C WRITE(14,1011) C J = 1 C DO 6 I = 1,NFRE C XCOR = FRECOR(I,1) * 1000.D0 C YCOR = FRECOR(I,2) * 1000.D0 C ZCOR = FRECOR(I,3) * 1000.D0 C XDEV = DSQRT(RN(J,J)) * 1000.D0 C YDEV = DSQRT(RN(J+1,J+1)) * 1000.D0 C ZDEV = DSQRT(RN(J+2,J+2)) * 1000.D0 C WRITE(14,1012) FRESTA(I),XCOR,XDEV,YCOR,YDEV,ZCOR,ZDEV C J = J + 3 C 6 CONTINUE ELSE C C PRINT FINAL ADJUSTED COORDINATES WRITE (6,1009) WRITE (6,1005) WRITE (6,1006) DO 3 I = 1,NFRE WRITE (6,1007) FRESTA(I),FRECOR(I,1),FRECOR(I,2),FRECOR(I,3) C WRITE (16,1010) FRESTA(I),FRECOR(I,1),FRECOR(I,2),FRECOR(I,3) 3 CONTINUE WRITE (6,1008) WRITE (6,1006) DO 4 I = 1,NFIX WRITE (6,1007) FIXSTA(I),FIXCOR(I,1),FIXCOR(I,2),FIXCOR(I,3) C WRITE (14,1010) FIXSTA(I),FIXCOR(I,1),FIXCOR(I,2),FIXCOR(I,3) 4 CONTINUE C C PRINT FINAL ADJUSTED COORDINATES AND CORRESPONDING STANDARD C DEVIATIONS TO A FILE WRITE(14,1009) WRITE(14,1011) J = 1 DO 5 I = 1,NFRE XCOR = FRECOR(I,1) * 1000.D0 YCOR = FRECOR(I,2) * 1000.D0 ZCOR = FRECOR(I,3) * 1000.D0 XDEV = DSQRT(RN(J,J)) * 1000.D0 YDEV = DSQRT(RN(J+1,J+1)) * 1000.D0 ZDEV = DSQRT(RN(J+2,J+2)) * 1000.D0 WRITE(14,1012) FRESTA(I),XCOR,XDEV,YCOR,YDEV,ZCOR,ZDEV J = J + 3 5 CONTINUE ENDIF C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('1',80('*')/1X,'*',78X,'*'/1X,'*',A78,'*'/1X,'*',78X,'*'/1X &,80('*')) 1001 FORMAT('0',T27,'F I L E M A N A G E M E N T'/T27,28('-')//T28,'IN &PUT FILE NAME : ',A8//T28,'OUTPUT FILE NAME: ',A8///) 1002 FORMAT('0',T24,'P R O J E C T M A N A G E M E N T'/T24,34('-')//T &12,'PROJECT NAME : ',A40//T12,'PROJECT NUMBER : ',A40//T12,'PR &OJECT LOCATION: ',A40///) 1003 FORMAT('0',T16,'S U M M A R Y O F S E L E C T E D O P T I O N S &'/T16,50('-')//) 1004 FORMAT('1',T11,'I N I T I A L A P P R O X I M A T E C O O R D I &N A T E S'/T11,59('-')///) 1005 FORMAT(' ',T34,'FREE STATIONS'/T34,13('-')) 1006 FORMAT('0',T13,'STATION',T25,'X COORDINATE',T41,'Y COORDINATE',T57 &,'Z COORDINATE') 1007 FORMAT('0',T13,A8,T25,F12.7,T41,F12.7,T57,F12.7) 1008 FORMAT('0',T33,'FIXED STATIONS'/T33,14('-')) 1009 FORMAT('1',T17,'F I N A L A D J U S T E D C O O R D I N A T E S' &/T17,49('-')///) 1010 FORMAT(1X,A8,2X,3(F16.4,4X)) 1011 FORMAT('0',T5,'STATION',T17,'X (MM)',T27,'+/- (MM)',T38,'Y (MM)',T &47,'+/- (MM)',T59,'Z (MM)',T68,'+/- (MM)') 1012 FORMAT('0',T5,A8,T15,F9.2,T29,F5.3,T36,F9.2,T49,F5.3,T57,F9.2,T70, &F5.3) 1013 FORMAT(' ','DATE: ',I4,1X,I2,1X,I2,50X,'TIME: ',I2,':',I2,':',I2// &/) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: NORVEC * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO CONTROL THE * C * COMPUTATION AND FORMATION OF THE NORMAL EQUATIONS MATRIX (N) AND * C * THE CONSTANT VECTOR (W). * C * * C * * C * SUBROUTINES REQUIRED: AZIM, DIRN, DIST, ZENI. * C * * C * FUNCTION SUBPROGRAMS REQUIRED: N/A * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL NORVEC(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,ICA,NFIX, * C * NFRE,NOBS,NUKN,OBCODE,OBS,STA,A,RN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * IND I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (DEGREES) * C * INM I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (MINUTES) * C * INS R*8 ARRAY OF ORTHOGONAL ZENITH ANGLES (SECONDS) * C * ISET I*2 CURRENT OBSERVATION SET NUMBER BEING PROCESSED * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * I0 R*8 INCLINATION COMPONENT IN THE DIRECTION 0 DEGREES * C * I1 R*8 INCLINATION COMPONENT IN THE DIRECTION 90 DEGREES * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK, @PROCESS RUBLK) SUBROUTINE NORVEC(APOPT,FIXCOR,FIXSTA,FRECOR,FRESTA,ICA,NFIX,NFRE, &NIDB,NOBS,NUKN,OBCODE,OBS,STA,A,W,NFI,NFR,NOB,NUK,ITER) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FIXSTA(NFI),FRESTA(NFR),STA(NOB,3) INTEGER*2 APOPT INTEGER*4 OBCODE(NOB) C C ASSIGN MATRIX AND VECTOR DIMENSIONS COMMON /RNBLK/ RN(1000,1000) COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),FIXCOR(NFI,3),FRECOR(NFR,3),ICA(NOB,6) DIMENSION OBS(NOB,4),W(NOB) C C INITIALIZE THE FIRST DESIGN MATRIX (A) DO 5 I = 1,NOBS DO 10 J = 1,6 A(I,J) = 0.D0 ICA(I,J) = 0 10 CONTINUE 5 CONTINUE C C INITIALIZE MATRIX OF NORMAL EQUATIONS AND CONSTANT VECTOR DO 15 I = 1,NUKN RU(I) = 0.D0 DO 20 J = 1,NUKN RN(I,J) = 0.D0 20 CONTINUE 15 CONTINUE C C FORM THE MATRIX OF NORMAL EQUATIONS I = 1 NIDB = 0 C C C ADD CONTRIBUTION OF DISTANCES TO NORMAL EQUATIONS & CONSTANT VECTOR IF (APOPT .EQ. 0 .AND. ITER .EQ. 0) WRITE(6,1000) 1 IOB = OBCODE(I) GOTO (2,3,4,7,8),IOB 2 CALL DIST(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE,NOBS,NUKN, &OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) GOTO 6 C C ADD CONTRIBUTION OF DIRECTIONS TO NORMAL EQUATIONS & CONSTANT VECTOR 3 CALL DIRN(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE,NOBS,NUKN, &OBCODE,OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) NIDB = NIDB + 1 GOTO 6 4 CONTINUE C C ADD CONTRIBUTION OF AZIMUTHS TO NORMAL EQUATIONS & CONSTANT VECTOR 7 CALL AZIM(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE,NOBS,NUKN, &OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) GOTO 6 C C ADD CONTRIBUTION OF ZENITH DIST. TO NORMAL EQUATIONS & CONSTANT VECTOR 8 CALL ZENI(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE,NOBS,NUKN, &OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) 6 IF (I .LE. NOBS) GOTO 1 C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('1',T13,'SUMMARY OF INPUT OBSERVATIONS AND INITIAL MISCLOSU &RES'/T13,53('-')///1X,'OBSERVATION',T18,'AT',T27,'FROM',T38,'TO',T &47,'OBSERVED',T59,'STD.DEV.',T69,'MISCLOSURE') C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: NORADD * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * COEFFICIENTS FOR THE FIRST DESIGN MATRIX (A) FOR DISTANCES. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL NORADD(A,I,ICA,NOBS,NUKN,P,RN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 12 23 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 FIRST DESIGN MATRIX (A) * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * FIXCOR R*8 MATRIX OF FIXED STATION COORDINATES * C * FIXSTA C*8 ARRAY CONTAINING FIXED STATION NAMES * C * FRECOR R*8 MATRIX OF FREE STATION COORDINATES * C * FRESTA C*8 ARRAY CONTAINING FREE STATION NAMES * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 UPDATED FIRST DESIGN MATRIX (A) * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C C @PROCESS DC(RNBLK) SUBROUTINE NORADD(A,I,ICA,NOBS,NUKN,P,NOB,NUK) C C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) DIMENSION A(NOB,6),ICA(NOB,6) C C ADD CONTRIBUTION TO MATRIX OF NORMAL EQUATIONS DO 5 J = 1,6 DO 5 K = 1,6 IF (ICA(I,J) .EQ. 0 .OR. ICA(I,K) .EQ. 0) GOTO 5 RN(ICA(I,J),ICA(I,K)) = RN(ICA(I,J),ICA(I,K)) + A(I,J) * A(I &,K) * P 5 CONTINUE C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: WVADD * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * COEFFICIENTS FOR THE FIRST DESIGN MATRIX (A) FOR DISTANCES. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL WVADD(A,I,ICA,P,RU,W) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 12 23 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 FIRST DESIGN MATRIX (A) * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * FIXCOR R*8 MATRIX OF FIXED STATION COORDINATES * C * FIXSTA C*8 ARRAY CONTAINING FIXED STATION NAMES * C * FRECOR R*8 MATRIX OF FREE STATION COORDINATES * C * FRESTA C*8 ARRAY CONTAINING FREE STATION NAMES * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 UPDATED FIRST DESIGN MATRIX (A) * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C C @PROCESS DC(RUBLK) SUBROUTINE WVADD(A,I,ICA,P,W,NOB,NUK) C C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) C C ASSIGN MATRIX DIMENSIONS COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),ICA(NOB,6) C C ADD CONTRIBUTION TO CONSTANT VECTOR (RU) DO 1 J = 1,6 IF (ICA(I,J) .EQ. 0) GOTO 1 RU(ICA(I,J)) = RU(ICA(I,J)) + A(I,J) * P * W 1 CONTINUE C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: UPDATE * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO ADD COMPUTED * C * ITERATIVE CORRECTIONS TO THE INITIAL APPROXIMATE COORDINATES OF * C * THE FREE STATION COORDINATE VARIATES. IF THE OPTION, NPRT, IS * C * ACTIVE, A SUMMARY OF THE OLD COORDINATES, ITERATIVE CORRECTIONS, * C * AND THE NEW COORDINATES WILL BE PRINTED WITH EACH ITERATION DURING * C * THE LEAST SQUARES ADJUSTMENT. * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL UPDATE(DELTA,FRECOR,NUKN,ITER) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1988 01 02 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * ITER I*2 ITERATION COUNTER * C * NUKN I*2 NUMBER OF UNKNOWN PARAMETERS * C * DELTA R*8 SOLUTION VECTOR * C * FRECOR R*8 MATRIX OF FREE STATION COORDINATES * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * FRECOR R*8 UPDATED FREE STATION COORDINATES * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * XC R*8 X COORDINATE VARIATE ITERATIVE CORRECTION * C * YC R*8 Y COORDINATE VARIATE ITERATIVE CORRECTION * C * ZC R*8 Z COORDINATE VARIATE ITERATIVE CORRECTION * C * NEWX R*8 UPDATED X COORDINATE VARIATE * C * NEWY R*8 UPDATED Y COORDINATE VARIATE * C * NEWZ R*8 UPDATED Z COORDINATE VARIATE * C * OLDX R*8 PREVIOUS X COORDINATE VARIATE * C * OLDY R*8 PREVIOUS Y COORDINATE VARIATE * C * OLDZ R*8 PREVIOUS Z COORDINATE VARIATE * C * * C ********************************************************************** C ********************************************************************** C C C @PROCESS DC(RUBLK) SUBROUTINE UPDATE(FRECOR,NFRE,ITER,NFR,NUK,NPRT,CNAM) C C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 CNAM(NFR) REAL*8 NEWX,NEWY,NEWZ C C ASSIGN MATRIX AND VECTOR DIMENSIONS COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION FRECOR(NFR,3) C C PRINT SUMMARY PAGE HEADING IF (ITER .EQ. 0) WRITE (6,1000) IF (ITER .GT. 0 .AND. NPRT .EQ. 1) GOTO 2 WRITE (6,1001) ITER C C UPDATE APPROXIMATE COORDINATES 2 J = 1 DO 1 I = 1,NFRE OLDX = FRECOR(I,1) OLDY = FRECOR(I,2) OLDZ = FRECOR(I,3) DX = -RU(J) DY = -RU(J+1) DZ = -RU(J+2) NEWX = OLDX + DX NEWY = OLDY + DY NEWZ = OLDZ + DZ FRECOR(I,1) = NEWX FRECOR(I,2) = NEWY FRECOR(I,3) = NEWZ IF (ITER .GT. 0 .AND. NPRT .EQ. 1) GOTO 1 IF (ITER .GE. 0) THEN WRITE (6,1002) CNAM(I),OLDX,DX,NEWX,OLDY,DY,NEWY,OLDZ,DZ,NEWZ ENDIF J = J + 3 1 CONTINUE C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('1',T8,'SUMMARY OF ITERATIVE CORRECTIONS TO INITIAL APPROXI &MATE COORDINATES'/T8,67('-')///) 1001 FORMAT(' ',T34,'ITERATION # ',I2//1X,'STATION',T18,'OLD COORDINATE &',T44,'CORRECTION',T66,'NEW COORDINATE'/) 1002 FORMAT(' ',A8,3(6X,'X = ',F14.6)/T10,3(6X,'Y = ',F14.6)/T10,3(6X,' &Z = ',F14.6)//) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: DIST * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * COEFFICIENTS FOR THE FIRST DESIGN MATRIX (A) FOR DISTANCES. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL DIST(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE, * C * NOBS,NUKN,OBS,STA,RN) * C * * C * COMMENTS: N/A * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 FIRST DESIGN MATRIX (A) * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * FIXCOR R*8 MATRIX OF FIXED STATION COORDINATES * C * FIXSTA C*8 ARRAY CONTAINING FIXED STATION NAMES * C * FRECOR R*8 MATRIX OF FREE STATION COORDINATES * C * FRESTA C*8 ARRAY CONTAINING FREE STATION NAMES * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 UPDATED FIRST DESIGN MATRIX (A) * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK, @PROCESS RUBLK) SUBROUTINE DIST(A,FIXCOR,FIXSTA,FRECOR,FRESTA,IAR,ICA,NFIX,NFRE,NO &BS,NUKN,OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FIXSTA(NFI),FRESTA(NFR),FRSTA,STA(NOB,3),TOSTA INTEGER*2 APOPT C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000) COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),FIXCOR(NFI,3),FRECOR(NFR,3),ICA(NOB,6) DIMENSION OBS(NOB,4),W(NOB) C C DETERMINE THE FIXED STATION COORDINATES FRSTA = STA(IAR,2) TOSTA = STA(IAR,3) DO 5 I = 1,NFIX IF (FRSTA .EQ. FIXSTA(I)) THEN XFR = FIXCOR(I,1) YFR = FIXCOR(I,2) ZFR = FIXCOR(I,3) ENDIF IF (TOSTA .EQ. FIXSTA(I)) THEN XTO = FIXCOR(I,1) YTO = FIXCOR(I,2) ZTO = FIXCOR(I,3) ENDIF 5 CONTINUE C C DETERMINE THE FREE STATION COORDINATES ICFR = 0 ICTO = 0 DO 10 I = 1,NFRE IF (FRSTA .EQ. FRESTA(I)) THEN XFR = FRECOR(I,1) YFR = FRECOR(I,2) ZFR = FRECOR(I,3) ICFR = I * 3 ENDIF IF (TOSTA .EQ. FRESTA(I)) THEN XTO = FRECOR(I,1) YTO = FRECOR(I,2) ZTO = FRECOR(I,3) ICTO = I * 3 ENDIF 10 CONTINUE C C COMPUTE SPATIAL DISTANCE BETWEEN STATIONS SFRTO = DSQRT((XTO - XFR)**2 + (YTO - YFR)**2 + (ZTO - ZFR)**2) C C COMPUTE THE COEFFICIENTS OF THE FIRST DESIGN MATRIX IF (ICFR .GE. 1) THEN A(IAR,1) = (XFR - XTO) / SFRTO A(IAR,2) = (YFR - YTO) / SFRTO A(IAR,3) = (ZFR - ZTO) / SFRTO ICA(IAR,1) = ICFR - 2 ICA(IAR,2) = ICFR - 1 ICA(IAR,3) = ICFR ENDIF IF (ICTO .GE. 1) THEN A(IAR,4) = (XTO - XFR) / SFRTO A(IAR,5) = (YTO - YFR) / SFRTO A(IAR,6) = (ZTO - ZFR) / SFRTO ICA(IAR,4) = ICTO - 2 ICA(IAR,5) = ICTO - 1 ICA(IAR,6) = ICTO ENDIF C P = 1.D0 / (OBS(IAR,1)**2) CALL NORADD(A,IAR,ICA,NOBS,NUKN,P,NOB,NUK) C C COMPUTE CONTRIBUTION TO MISCLOSURE VECTOR IF ADJUSTMENT IF (APOPT .EQ. 1) GOTO 5000 W(IAR) = SFRTO - OBS(IAR,3) IF (ITER .GT. 0) GOTO 5001 STD = OBS(IAR,1) * 1000.D0 IF (APOPT .EQ. 0) WRITE(6,1000) STA(IAR,1),STA(IAR,2),STA(IAR,3),O &BS(IAR,3),STD,W(IAR) 5001 CALL WVADD(A,IAR,ICA,P,W(IAR),NOB,NUK) 5000 IAR = IAR + 1 C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('0','SLOPE DIST.',T15,A8,T25,A8,T35,A8,T45,F12.5,T60,F6.4,T &68,F11.3) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: DIRN * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * COEFFICIENTS FOR THE FIRST DESIGN MATRIX (A) FOR DIRECTIONS. * C * * C * * C * * C * SUBROUTINES REQUIRED: NORADD * C * * C * FUNCTION SUBPREGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL DIRN(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE, * C * NOBS,NUKN,OBS,STA,RN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * IND I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (DEGREES) * C * INM I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (MINUTES) * C * INS R*8 ARRAY OF ORTHOGONAL ZENITH ANGLES (SECONDS) * C * ISET I*2 CURRENT OBSERVATION SET NUMBER BEING PROCESSED * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * I0 R*8 INCLINATION COMPONENT IN THE DIRECTION 0 DEGREES * C * I1 R*8 INCLINATION COMPONENT IN THE DIRECTION 90 DEGREES * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK, @PROCESS RUBLK) SUBROUTINE DIRN(A,FIXCOR,FIXSTA,FRECOR,FRESTA,IAR,ICA,NFIX,NFRE,NO &BS,NUKN,OBCODE,OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FIXSTA(NFI),FRESTA(NFR),FRSTA,STA(NOB,3),TOSTA INTEGER*2 APOPT INTEGER*4 OBCODE(NOB) C C ASSIGN MATRIX AND VECTOR DIMENSIONS COMMON /RNBLK/ RN(1000,1000) COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),FIXCOR(NFI,3),FRECOR(NFR,3),ICA(NOB,6) DIMENSION OBS(NOB,4),P(40),U(40,40) DIMENSION W(NOB) C C COMPUTE DATA CONSTANTS AND CONVERSION FACTORS PI = DARCOS(-1.D0) TWOPI = 2.D0 * PI RO = 3600.D0 * 180.D0 / PI C C DETERMINE NUMBER OF DIRECTIONS WITHIN THE DIRECTION BUNDLE II = IAR + 40 DO 1 J = IAR,II M = J IF (OBCODE(J) .EQ. -2 .OR. IABS(OBCODE(J)) .NE. 2) GOTO 100 1 CONTINUE 100 NUM = M - IAR + 1 C C COMPUTE AND FORM FIRST DESIGN MATRIX (A) C WRITE(6,*) IAR,M,NUM DO 2 J = IAR,M FRSTA = STA(J,2) TOSTA = STA(J,3) C C DETERMINE THE FIXED STATION COORDINATES DO 5 I = 1,NFIX IF (FRSTA .EQ. FIXSTA(I)) THEN XFR = FIXCOR(I,1) YFR = FIXCOR(I,2) ENDIF IF (TOSTA .EQ. FIXSTA(I)) THEN XTO = FIXCOR(I,1) YTO = FIXCOR(I,2) ENDIF 5 CONTINUE C C DETERMINE THE FREE STATION COORDINATES ICFR = 0 ICTO = 0 DO 10 I = 1,NFRE IF (FRSTA .EQ. FRESTA(I)) THEN XFR = FRECOR(I,1) YFR = FRECOR(I,2) ICFR = I * 3 ENDIF IF (TOSTA .EQ. FRESTA(I)) THEN XTO = FRECOR(I,1) YTO = FRECOR(I,2) ICTO = I * 3 ENDIF 10 CONTINUE C C COMPUTE HORIZONTAL DISTANCE BETWEEN STATIONS (SQUARED) SFRTO = (XTO - XFR)**2 + (YTO - YFR)**2 C C COMPUTE THE COEFFICIENTS OF THE FIRST DESIGN MATRIX IF (ICFR .GE. 1) THEN A(J,1) = (YFR - YTO) / SFRTO * RO A(J,2) = (XTO - XFR) / SFRTO * RO ICA(J,1) = ICFR - 2 ICA(J,2) = ICFR - 1 ENDIF IF (ICTO .GE. 1) THEN A(J,3) = (YTO - YFR) / SFRTO * RO A(J,4) = (XFR - XTO) / SFRTO * RO ICA(J,3) = ICTO - 2 ICA(J,4) = ICTO - 1 ENDIF C C COMPUTE CORRESPONDING OBSERVATION WEIGHT K = J - IAR + 1 P(K) = 1.D0 / (OBS(J,1)**2) C C DETERMINE IF ADJUSTMENT OR PRE-ANALYSIS IF (APOPT .EQ. 1) GOTO 2 C C COMPUTE ORIENTATION UNKNOWN OF DIRECTION BUNDLE AND ADD CONTRIBUTION C TO MISCLOSURE VECTOR IF ADJUSTMENT IF (J .GT. IAR) GOTO 3 D1 = (OBS(J,2) + OBS(J,3)/60.D0 + OBS(J,4)/3600.D0) * PI/180.D0 Z = DATAN2((XTO - XFR),(YTO - YFR)) IF (Z .LT. 0.D0) Z = Z + TWOPI 3 AL = DATAN2((XTO - XFR),(YTO - YFR)) IF (AL .LT. 0.D0) AL = AL + TWOPI IF (AL .LT. Z) AL = AL + TWOPI W(J) = AL - Z - (OBS(J,2) + OBS(J,3)/60.D0 + OBS(J,4)/3600.D0) * &PI / 180.D0 + D1 W(J) = W(J) * RO IF (ITER .GT. 0) GOTO 2 ID = OBS(J,2) IM = OBS(J,3) IF (APOPT .EQ. 0) WRITE(6,1000) STA(J,1),STA(J,2),STA(J,3),ID,IM,O &BS(J,4),OBS(J,1),W(J) 2 CONTINUE C C COMPUTE TRANSFORMATION MATRIX (U) FOR WEIGHT MATRIX (P) SUM = 0.D0 DO 6 J = 1,NUM SUM = SUM + P(J) 6 CONTINUE DO 7 J = 1,NUM DO 7 K = 1,NUM U(J,K) = -P(J) * P(K) / SUM IF (J .EQ. K) U(J,K) = U(J,K) + P(K) 7 CONTINUE C C ADD CONTRIBUTION OF DIRECTION BUNDLE TO MATRIX OF NORMAL EQUATIONS C CALL MOUTD(U,30,NUM,NUM) DO 8 J = IAR,M DO 8 K = IAR,M DO 8 L1 = 1,4 DO 8 L2 = 1,4 C IF(ICA(J,L1) .GT. ICA(K,L2)) GOTO 8 IF(ICA(J,L1) .EQ. 0 .OR. ICA(K,L2) .EQ. 0) GOTO 8 RN(ICA(J,L1),ICA(K,L2)) = RN(ICA(J,L1),ICA(K,L2)) + A(J,L1) * A(K, &L2) * U(J-IAR+1,K-IAR+1) 8 CONTINUE C C ADD CONTRIBUTION TO CONSTANT VECTOR(RU) IF ADJUSTMENT IF (APOPT .EQ. 1) GOTO 5000 DO 9 J = IAR,M DO 9 K = 1,NUM DO 9 L = 1,4 IF (ICA(J,L) .EQ. 0) GOTO 9 RU(ICA(J,L)) = RU(ICA(J,L)) + A(J,L)*U(J-IAR+1,K)*W(K+IAR-1) 9 CONTINUE 5000 IAR = IAR + NUM C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('0','DIRECTION',T15,A8,T25,A8,T35,A8,T45,I3,T49,I2,T52,F5.2 &,T60,F5.2,T69,F10.2) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: AZIM * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * COEFFICIENTS OF THE FIRST DESIGN MATRIX (A) FOR AZIMUTHS. * C * * C * * C * * C * SUBROUTINES REQUIRED: N/A * C * * C * FUNCTION SUBPROGRAMS REQUIRED: N/A * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL AZIM(A,FIXCOR,FIXSTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE, * C * NOBS,NUKN,OBS,STA,RN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * IND I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (DEGREES) * C * INM I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (MINUTES) * C * INS R*8 ARRAY OF ORTHOGONAL ZENITH ANGLES (SECONDS) * C * ISET I*2 CURRENT OBSERVATION SET NUMBER BEING PROCESSED * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * I0 R*8 INCLINATION COMPONENT IN THE DIRECTION 0 DEGREES * C * I1 R*8 INCLINATION COMPONENT IN THE DIRECTION 90 DEGREES * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK, @PROCESS RUBLK) SUBROUTINE AZIM(A,FIXCOR,FIXSTA,FRECOR,FRESTA,IAR,ICA,NFIX,NFRE,NO &BS,NUKN,OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FIXSTA(NFI),FRESTA(NFR),FRSTA,STA(NOB,3),TOSTA INTEGER*2 APOPT C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),FIXCOR(NFI,3),FRECOR(NFR,3),ICA(NOB,6) DIMENSION OBS(NOB,4),W(NOB) C C COMPUTE DATA CONSTANTS AND CONVERSION FACTORS PI = DARCOS(-1.D0) TWOPI = 2.D0 * PI RO = 3600.D0 * 180.D0 / PI C C DETERMINE THE FIXED STATION COORDINATES FRSTA = STA(IAR,2) TOSTA = STA(IAR,3) DO 5 I = 1,NFIX IF (FRSTA .EQ. FIXSTA(I)) THEN XFR = FIXCOR(I,1) YFR = FIXCOR(I,2) ENDIF IF (TOSTA .EQ. FIXSTA(I)) THEN XTO = FIXCOR(I,1) YTO = FIXCOR(I,2) ENDIF 5 CONTINUE C C DETERMINE THE FREE STATION COORDINATES ICFR = 0 ICTO = 0 DO 10 I = 1,NFRE IF (FRSTA .EQ. FRESTA(I)) THEN XFR = FRECOR(I,1) YFR = FRECOR(I,2) ICFR = I * 3 ENDIF IF (TOSTA .EQ. FRESTA(I)) THEN XTO = FRECOR(I,1) YTO = FRECOR(I,2) ICTO = I * 3 ENDIF 10 CONTINUE C C COMPUTE HORIZONTAL DISTANCE BETWEEN STATIONS (SQUARED) SFRTO = (XTO - XFR)**2 + (YTO - YFR)**2 C C COMPUTE THE COEFFICIENTS OF THE FIRST DESIGN MATRIX IF (ICFR .GE. 1) THEN A(IAR,1) = (YFR - YTO) / SFRTO * RO A(IAR,2) = (XTO - XFR) / SFRTO * RO ICA(IAR,1) = ICFR - 2 ICA(IAR,2) = ICFR - 1 ENDIF IF (ICTO .GE. 1) THEN A(IAR,3) = (YTO - YFR) / SFRTO * RO A(IAR,4) = (XFR - XTO) / SFRTO * RO ICA(IAR,3) = ICTO - 2 ICA(IAR,4) = ICTO - 1 ENDIF C P = 1.D0 / (OBS(IAR,1)**2) CALL NORADD(A,IAR,ICA,NOBS,NUKN,P,NOB,NUK) C C COMPUTE CONTRIBUTION TO MISCLOSURE VECTOR (W) IF ADJUSTMENT IF (APOPT .EQ. 1) GOTO 5000 AL = DATAN2((XTO - XFR),(YTO - YFR)) IF (AL .LT. 0.D0) AL = AL + TWOPI W(IAR) = AL - (OBS(IAR,2) + OBS(IAR,3) / 60.D0 + OBS(IAR,4) / 3600 &.D0) * PI / 180.D0 W(IAR) = W(IAR) * RO IF (ITER .GT. 0) GOTO 5001 ID = OBS(IAR,2) IM = OBS(IAR,3) IF (APOPT .EQ. 0) WRITE(6,1000) STA(IAR,1),STA(IAR,2),STA(IAR,3),I &D,IM,OBS(IAR,4),OBS(IAR,1),W(IAR) 5001 CALL WVADD(A,IAR,ICA,P,W(IAR),NOB,NUK) 5000 IAR = IAR + 1 C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('0','AZIMUTH',T15,A8,T25,A8,T35,A8,T45,I3,T49,I2,T52,F5.2,T &60,F5.2,T69,F10.2) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: ZENI * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * COEFFICIENTS OF THE FIRST DESIGN MATRIX (A) FOR ZENITH ANGLES. * C * * C * * C * * C * SUBROUTINES REQUIRED: N/A * C * * C * FUNCTION SUBPROGRAMS REQUIRED: N/A * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL ZENI(A,FIXCOR,FXISTA,FRECOR,FRESTA,I,ICA,NFIX,NFRE, * C * NOBS,NUKN,OBS,STA,RN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * IND I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (DEGREES) * C * INM I*2 ARRAY OF ORTHOGONAL ZENITH ANGLES (MINUTES) * C * INS R*8 ARRAY OF ORTHOGONAL ZENITH ANGLES (SECONDS) * C * ISET I*2 CURRENT OBSERVATION SET NUMBER BEING PROCESSED * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * I0 R*8 INCLINATION COMPONENT IN THE DIRECTION 0 DEGREES * C * I1 R*8 INCLINATION COMPONENT IN THE DIRECTION 90 DEGREES * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK, @PROCESS RUBLK) SUBROUTINE ZENI(A,FIXCOR,FIXSTA,FRECOR,FRESTA,IAR,ICA,NFIX,NFRE,NO &BS,NUKN,OBS,STA,APOPT,W,NFI,NFR,NOB,NUK,ITER) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) INTEGER*2 APOPT CHARACTER*8 FIXSTA(NFI),FRESTA(NFR),FRSTA,STA(NOB,3),TOSTA C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),FIXCOR(NFI,3),FRECOR(NFR,3),ICA(NOB,6) DIMENSION OBS(NOB,4),W(NOB) C C COMPUTE DATA CONSTANTS AND CONVERSION FACTORS PI = DARCOS(-1.D0) RO = 3600.D0 * 180.D0 / PI C C DETERMINE THE FIXED STATION COORDINATES FRSTA = STA(IAR,2) TOSTA = STA(IAR,3) DO 5 I = 1,NFIX IF (FRSTA .EQ. FIXSTA(I)) THEN XFR = FIXCOR(I,1) YFR = FIXCOR(I,2) ZFR = FIXCOR(I,3) ENDIF IF (TOSTA .EQ. FIXSTA(I)) THEN XTO = FIXCOR(I,1) YTO = FIXCOR(I,2) ZTO = FIXCOR(I,3) ENDIF 5 CONTINUE C C DETERMINE THE FREE STATION COORDINATES ICFR = 0 ICTO = 0 DO 10 I = 1,NFRE IF (FRSTA .EQ. FRESTA(I)) THEN XFR = FRECOR(I,1) YFR = FRECOR(I,2) ZFR = FRECOR(I,3) ICFR = I * 3 ENDIF IF (TOSTA .EQ. FRESTA(I)) THEN XTO = FRECOR(I,1) YTO = FRECOR(I,2) ZTO = FRECOR(I,3) ICTO = I * 3 ENDIF 10 CONTINUE C C COMPUTE HORIZONTAL/SPATIAL DISTANCES BETWEEN STATIONS (SQUARED) SFRTO = (XTO - XFR)**2 + (YTO - YFR)**2 SFRTO1 = DSQRT(SFRTO) SFRTO2 = SFRTO + (ZTO - ZFR)**2 C SFRTO3 = DSQRT(SFRTO2) C WRITE(6,1010) C1010 FORMAT(' ','HORZ. & SLOPE DISTANCES') C WRITE(6,*) SFRTO1,SFRTO3 C C COMPUTE THE COEFFICIENTS OF THE FIRST DESIGN MATRIX IF (ICFR .GE. 1) THEN A(IAR,1) =-((ZTO - ZFR) * (XTO - XFR)) / (SFRTO2 * SFRTO1) * RO A(IAR,2) =-((ZTO - ZFR) * (YTO - YFR)) / (SFRTO2 * SFRTO1) * RO A(IAR,3) =-(-SFRTO1 / SFRTO2) * RO ICA(IAR,1) = ICFR - 2 ICA(IAR,2) = ICFR - 1 ICA(IAR,3) = ICFR ENDIF IF (ICTO .GE. 1) THEN A(IAR,4) =-((ZFR - ZTO) * (XTO - XFR)) / (SFRTO2 * SFRTO1) * RO A(IAR,5) =-((ZFR - ZTO) * (YTO - YFR)) / (SFRTO2 * SFRTO1) * RO A(IAR,6) =-(SFRTO1 / SFRTO2) * RO ICA(IAR,4) = ICTO - 2 ICA(IAR,5) = ICTO - 1 ICA(IAR,6) = ICTO ENDIF C P = 1.D0 / (OBS(IAR,1)**2) CALL NORADD(A,IAR,ICA,NOBS,NUKN,P,NOB,NUK) C C ADD CONTRIBUTION TO MISCLOSURE VECTOR(W) IF ADJUSTMENT IF (APOPT .EQ. 1) GOTO 5000 ARG = (ZTO - ZFR) / SFRTO1 VA = DATAN(ARG) ZA = PI / 2.D0 - VA W(IAR) = ZA - (OBS(IAR,2) + OBS(IAR,3) / 60.D0 + OBS(IAR,4) / 3600 &.D0) * PI / 180.D0 W(IAR) = W(IAR) * RO IF (ITER .GT. 0) GOTO 5001 ID = OBS(IAR,2) IM = OBS(IAR,3) IF (APOPT .EQ. 0) WRITE(6,1000) STA(IAR,1),STA(IAR,2),STA(IAR,3),I &D,IM,OBS(IAR,4),OBS(IAR,1),W(IAR) 5001 CALL WVADD(A,IAR,ICA,P,W(IAR),NOB,NUK) 5000 IAR = IAR + 1 C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('0','ZENITH ANGLE',T15,A8,T25,A8,T35,A8,T45,I3,T49,I2,T52,F &5.2,T60,F5.2,T69,F10.2) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: STATSUM * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO READ IN THE * C * OBSERVATION DATA INPUT FILE. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL STATSUM( * C * * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * APOP I*2 ADJUSTMENT/PRE-ANALYSIS OPTION CODE * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * OBS R*8 MATRIX CONTAINING OBSERVATIONS * C * STA C*8 MATRIX OF AT/FROM/TO STATION NAMES * C * NOBS I*2 NUMBER OF OBSERVATIONS READ IN * C * OBCODE I*2 ARRAY OF OBSERVATION CODE NUMBERS * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C SUBROUTINE STATSUM(APOST,FAC,FAIL,IDF,ITER,NIO,NIDB,NOBS,NUKN) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) INTEGER*2 NIO(5) LOGICAL FAIL C C PRINT PAGE HEADER AND STATISTICS SUMMARY NUKNS = NUKN + NIDB WRITE (6,1000) WRITE (6,1001) ITER,5 WRITE (6,1002) NIO(1),NIDB,NIO(2),NIO(3),NIO(4),NIO(5),NUKN,NOBS,N &UKNS WRITE (6,1003) IDF C C PERFORM CHI-SQUARE TEST ON THE A PRIORI VARIANCE FACTOR CALL CHITST(0.05D0,APOST,1.D0,IDF,FAC,FAIL) C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('1',T23,'S T A T I S T I C A L S U M M A R Y'/T23,36('-')/ &//) 1001 FORMAT('0',T14,'NUMBER OF ITERATIONS REQUIRED FOR CONVERGENCE --> & ',I3//T14,'MAXIMUM NUMBER OF ITERATIONS ALLOWED ',11('-'),'> ',I &3/) 1002 FORMAT('0',T11,'NUMBER OF OBSERVATIONS',T38,'| NUMBER OF UNKNOWN P &ARAMETERS'/T11,27('-'),'|',32('-')//T11,'SLOPE DISTANCES',T31,I5,2 &X,'|',' ORIENTATION',T66,I5/T11,'HORZ. DIRECTIONS',T31,I5,' |'/T1 &1,'HORZ. ANGLES',T31,I5,' |'/T11,'AZIMUTHS',T31,I5,' |'/T11,'ZEN &ITH ANGLES',T31,I5,' | COORDINATES',T66,I5/T31,5('-'),T66,5('-')/ &T11,'TOTALS',T31,I5,T66,I5/) 1003 FORMAT('0',T20,'TOTAL NUMBER OF DEGREES OF FREEDOM = ',I4/) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: RESID * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS PROGRAM IS TO COMPUTE THE VECTOR * C * OF ADJUSTED RESIDUALS. * C * * C * * C * * C * SUBROUTINES REQUIRED: NONE * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: * C * * C * * C * COMMENTS: N/A * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 12 31 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 FIRST DESIGN MATRIX (A) * C * NFIX I*2 NUMBER OF FIXED STATIONS * C * NFRE I*2 NUMBER OF FREE STATIONS * C * FIXCOR R*8 MATRIX OF FIXED STATION COORDINATES * C * FIXSTA C*8 ARRAY CONTAINING FIXED STATION NAMES * C * FRECOR R*8 MATRIX OF FREE STATION COORDINATES * C * FRESTA C*8 ARRAY CONTAINING FREE STATION NAMES * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * A R*8 UPDATED FIRST DESIGN MATRIX (A) * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * SUM R*8 ARRAY OF ORTHOGONAL MEASUREMENT CHECK SUMS * C * INCL R*8 ARRAY OF ORTHOGONAL MEASUREMENTS (DECIMAL DEGREES) * C * IPTR I*2 POINTER FOR INPUT ARRAYS CORRESPONDING TO ISET * C * CHECK L*0 LOGICAL ARGUMENT FOR PASS/FAIL OF CHECK TOLERANCE * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RUBLK) SUBROUTINE RESID(A,ICA,NOBS,OBCODE,OBS,RHAT,W,NOB,NUK,IDF,APOST,RS &TD,RHATS,STA) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 STA(NOB,3) INTEGER*4 OBCODE(NOB) C C ASSIGN MATRIX AND VECTOR DIMENSIONS COMMON /RUBLK/ RU(1000),IPVT(1000) DIMENSION A(NOB,6),ICA(NOB,6),OBS(NOB,4),RHAT(NOB),RHATS(NOB) DIMENSION RSTD(NOB),W(NOB) C C DETERMINE OBSERVATION TYPE QSUM = 0.D0 I = 1 10 IOB = OBCODE(I) GOTO(1,3,1,1,1) IOB C C COMPUTE DISTANCE, AZIMUTH AND ZENITH ANGLE RESIDUALS 1 W1 = 0.D0 DO 2 J = 1,6 IF (ICA(I,J) .EQ. 0) GOTO 2 W1 = W1 - A(I,J) * RU(ICA(I,J)) 2 CONTINUE RHAT(I) = W(I) + W1 QSUM = QSUM + RHAT(I)**2 / OBS(I,1)**2 I = I + 1 GOTO 24 C C COMPUTE DIRECTION RESIDUALS 3 II = I + 40 DO 4 J = I,II M = J IF (OBCODE(J) .EQ. -2) GOTO 5 4 CONTINUE 5 NUM = M - I + 1 SUM = 0.D0 DO 7 J = I,M SUM = SUM + 1.D0 / OBS(J,1)**2 W1 = 0.D0 DO 6 K = 1,4 IF(ICA(J,K) .EQ. 0) GOTO 6 W1 = W1 - A(J,K) * RU(ICA(J,K)) 6 CONTINUE RHAT(J) = W(J) + W1 7 CONTINUE SUM1 = 0.D0 DO 8 J = I,M W1 = 0.D0 DO 15 K = 1,4 IF(ICA(J,K) .EQ. 0) GOTO 15 W1 = W1 + A(J,K) * RU(ICA(J,K)) 15 CONTINUE W1 = (W1 - W(J)) / OBS(J,1)**2 SUM1 = SUM1 + W1 8 CONTINUE SUM2 = SUM1 / SUM DO 9 J = I,M 9 RHAT(J) = RHAT(J) + SUM2 DO 13 J = I,M 13 QSUM = QSUM + RHAT(J)**2 / OBS(J,1)**2 I = I + NUM 24 IF (I .LE. NOBS) GOTO 10 C C COMPUTE STANDARDIZED RESIDUALS USING POPE'S APPROXIMATION DF = DFLOAT(IDF) APOST = QSUM / DF ARG1 = DSQRT(DF / DFLOAT(NOBS)) ARG2 = (DSQRT(APOST) / DSQRT(1.D0)) * ARG1 WRITE(6,1000) DO 25 I = 1,NOBS RSTD(I) = ARG2 * OBS(I,1) RHATS(I) = RHAT(I) / RSTD(I) IF (OBCODE(I) .NE. 1) THEN ID = OBS(I,2) IM = OBS(I,3) ENDIF IF (OBCODE(I) .EQ. 1) RHAT(I) = RHAT(I) * 1000.D0 IF (OBCODE(I) .EQ. 1) OBS(I,1) = OBS(I,1) * 1000.D0 IF (OBCODE(I) .EQ. 1) WRITE(6,1001) I,STA(I,1),STA(I,2),STA(I,3 &),OBS(I,3),OBS(I,1),RHAT(I) IF (IABS(OBCODE(I)) .EQ. 2) WRITE(6,1002) I,STA(I,1),STA(I,2),S &TA(I,3),ID,IM,OBS(I,4),OBS(I,1),RHAT(I) IF (OBCODE(I) .EQ. 4) WRITE(6,1003) I,STA(I,1),STA(I,2),STA(I,3 &),ID,IM,OBS(I,4),OBS(I,1),RHAT(I) IF (OBCODE(I) .EQ. 5) WRITE(6,1004) I,STA(I,1),STA(I,2),STA(I,3 &),ID,IM,OBS(I,4),OBS(I,1),RHAT(I) 25 CONTINUE C C INPUT/OUTPUT FORMAT STATEMENTS 1000 FORMAT('1',T22,'SUMMARY OF OBSERVATIONS AND RESIDUALS'/T22,37('-') &///T7,'OBSERVATION',T20,'AT',T29,'FROM',T38,'TO',T50,'OBSERVED',T6 &2,'STD.DEV.',T72,'RESIDUAL') 1001 FORMAT('0',I4,' SLOPE DIST.',T20,A8,T29,A8,T38,A8,T48,F12.5,T62,F7 &.4,T71,F10.5) 1002 FORMAT('0',I4,' DIRECTION',T20,A8,T29,A8,T38,A8,T48,I3,T52,I2,T55, &F5.2,T62,F5.2,T72,F9.5) 1003 FORMAT('0',I4,' AZIMUTH',T20,A8,T29,A8,T38,A8,T48,I3,T52,I2,T55,F5 &.2,T62,F5.2,T72,F9.5) 1004 FORMAT('0',I4,' ZENITH ANGLE',T20,A8,T29,A8,T38,A8,T48,I3,T52,I2,T &55,F5.2,T62,F5.2,T72,F9.5) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: CONF3D * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * THREE DIMENSIONAL CONFIDENCE ELLIPSOID AT A STATED SIGNIFICANCE * C * LEVEL. * C * * C * * C * SUBROUTINES REQUIRED: EIGEN3 (UNBCC-SELIB) * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL CONF3D(ALPHA,CX,DF,FAIL,NUKN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 28 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * CX R*8 VARIANCE-COVARIANCE MATRIX OF ADJUSTED PARAMETERS * C * DF R*8 NUMBER OF DEGREES OF FREEDOM * C * NUKN I*2 NUMBER OF UNKNOWN PARAMETERS * C * ALPHA R*8 SIGNIFICANCE LEVEL * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * * C * * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * S R*8 VECTOR CONTAINING EIGENVALUES (SEMIMAJOR AXES) * C * CXEX R*8 EXTRACTED PORTION OF VARIANCE-COVARIANCE MATRIX OF * C * ADJUSTED PARAMETERS * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK) SUBROUTINE CONF3D(ALPHA,IDF,FAIL,FRESTA,NUKN,NFR,NUK,APOST) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FRESTA(NFR) LOGICAL FAIL REAL*8 DCHIIN,DFIN C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) DIMENSION CXEX(3,3),S(3,1) C C EXTRACT APPROPRIATE PORTION OF VARIANCE-COVARIANCE MATRIX IFRE = 1 C C COMPUTE UNIVARIATE (ONE DIMENSIONAL) EXPANSION FACTOR TO 95.0% PROB = 1.D0 - ALPHA D1 = 1.D0 D2 = DFLOAT(IDF) IF (.NOT. FAIL) THEN TVAL = DCHIIN(PROB,D1) XFAC = DSQRT(TVAL) ELSE TVAL = DFIN(PROB,D1,D2) XFAC = DSQRT(TVAL) ENDIF C C CALCULATE MULTIVARIATE (THREE DIMENSIONAL) EXPANSION FACTOR TO 95.0% D1 = 3.D0 D2 = DFLOAT(IDF) IF (.NOT. FAIL) THEN TVAL = DCHIIN(PROB,D1) CBAR = DSQRT(TVAL) ELSE TVAL = DFIN(PROB,D1,D2) CBAR = DSQRT(3.D0 * TVAL) ENDIF WRITE(6,1000) IF (.NOT. FAIL) THEN WRITE (6,1001) CBAR WRITE (6,1002) XFAC IF (APOST .NE. 0.D0) WRITE (6,1003) APOST ELSE WRITE (6,1004) CBAR WRITE (6,1005) XFAC WRITE (6,1006) APOST ENDIF WRITE(6,1007) DO 5 I = 1,NUKN,3 I1 = 1 J1 = 1 K1 = I + 2 DO 10 J = I,K1 DO 15 K = I,K1 CXEX(I1,J1) = RN(J,K) J1 = J1 + 1 15 CONTINUE J1 = 1 I1 = I1 + 1 10 CONTINUE C CALL MOUTD(CXEX,3,3,3) C C COMPUTE EIGENVALUES AND EIGENVECTORS OF EXTRACTED COVARIANCE MATRIX CALL EIGEN3(CXEX,PROB,3.D0,IDF,FAIL,FRESTA,IFRE,S,CBAR,NFR,XFAC &) IFRE = IFRE + 1 5 CONTINUE C C INPUT/FORMAT STATEMENTS 1000 FORMAT('1',80('*')/1X,'*',78X,'*'/1X,'*',7X,'SUMMARY OF ABSOLUTE 9 &5.0% CONFIDENCE ELLIPSOIDS (OUT-OF-CONTEXT)',7X,'*'/1X,'*',78X,'*' &/1X,80('*')//) 1001 FORMAT('0','FACTOR USED FOR OBTAINING THESE ELLIPSOIDS FROM STANDA &RD ELLIPSOIDS: (VARIANCE FACTOR KNOWN) = ',F6.4/) 1002 FORMAT('0','FACTOR USED FOR OBTAINING XYZ CONFIDENCE INTERVALS: (V &ARIANCE FACTOR KNOWN) = ',F6.4/) 1003 FORMAT('0','COVARIANCE MATRIX OF PARAMETERS NOT MULTIPLIED BY A PO &STERIORI VARIANCE FACTOR (',F12.6,')'/T34,3('-')/) 1004 FORMAT('0','FACTOR USED FOR OBTAINING THESE ELLIPSOIDS FROM STANDA &RD ELLIPSOIDS: (VARIANCE FACTOR UNKNOWN) = ',F6.4/) 1005 FORMAT('0','FACTOR USED FOR OBTAINING XYZ CONFIDENCE INTERVALS: (V &ARIANCE FACTOR UNKNOWN) = ',F6.4/) 1006 FORMAT('0','COVARIANCE MATRIX OF PARAMETERS IS MULTIPLIED BY A POS &TERIORI VARIANCE FACTOR (',F12.6,')'/T34,2('-')/) 1007 FORMAT('0',40X,'SPATIAL ANGLES (IN DEGREES) BETWEEN AXES'//1X,'POI &NT',7X,'SEMI-MAJOR AXES (METRES)',16X,'X(MODEL)',2X,'Y(MODEL)',2X, &'Z(MODEL)',10X,'+/- XYZ VALUES (AT STD. AND 95%) (MM)') C RETURN END C ****************************************************************** C * * C * SUBROUTINE 'EIGEN3' COMPUTES THE EIGEN VALUES AND * C * EIGEN VECTORS FOR A 3-D CASE. * C * NOTE: THE IMSL SUBROUTINE MDFI USES REAL*4 PARAMETERS * C * * C * INPUT: * C * A - THE VARIANCE CO-VARIANCE MATRIX * C * PROB - IS THE PROBABILITY REQUESTED IN THE * C * EXCLUSIVE RANGE (0,1), * C * U & DF - ARE THE FIRST & SECOND DEGREES OF FREEDOM * C * RESPECTIVELY, * C * * C * OUTPUT: * C * A - AN ARRAY CONTAINING THE EIGEN VECTORS, * C * NOTE: THE VARIANCE CO-VARIANCE MATRIX IS * C * DESTROYED, * C * S - THE SEMI MAJOR AXES. * C * * C * IFLAG - THOUGH NOT AN INPUT/OUTPUT PARAMETER OF * C * EIGEN3 IT IS USED AS AN INPUT/OUTPUT ERROR * C * PARAMETER IN MDFI. EXPLANITION- * C * TERMINAL ERROR = 128 + N * C * WARNING ERROR = 32 + N * C * N=1 MEANS THAT AN ERROR OCCURED IN * C * MDBETI * C * N=2 MEANS PROB WAS NOT IN THE EXCLUSIVE * C * RANGE (0,1) * C * N=3 MEANS THE COMPUTED VALUE OF TVAL * C * WOULD PRODUCE AN OVERFLOW. TVAL IS * C * SET TO MACHINE INFINITY. * C * * C * * C * BY: H. MONIWA (1973) * C * MODIFIED BY: R. CALDWELL (1976) * C * * C ****************************************************************** C C SUBROUTINE EIGEN3 (A,PROB,U,IDF,FAIL,FRESTA,IFRE,S,CBAR,NFR,XFAC) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N) CHARACTER*1 XYZ(3)/'X','Y','Z'/ CHARACTER*8 FRESTA(NFR) LOGICAL FAIL REAL*8 SDEV(3),SDEVE(3) C C ASSIGN MATRIX DIMENSIONS DIMENSION A(3,3),R(9),EV(3),WS(6),S(3),B(3,3) EQUIVALENCE (R(1),B(1,1)) C C CONSTRUCT PARTIAL VARIANCE-COVARIANCE MATRIX A(2,1)=A(1,2) A(3,1)=A(1,3) A(3,2)=A(2,3) C C COMPUTE STANDARD DEVIATIONS IN XYZ COORDINATE AXIS DIRECTIONS SDEV(1) = DSQRT(A(1,1)) * 1000.D0 SDEV(2) = DSQRT(A(2,2)) * 1000.D0 SDEV(3) = DSQRT(A(3,3)) * 1000.D0 SDEVE(1) = SDEV(1) * XFAC SDEVE(2) = SDEV(2) * XFAC SDEVE(3) = SDEV(3) * XFAC C C CALCULATING THE EXPANDING FACTOR C C IF( PROB .EQ. 0.0E0 ) GOTO 40 C D1 = U C D2 = DF C TVAL = DFIN(PROB,D1,D2) C C CHECKING THE SUBROUTINE FOR ERRORS C C IF (.NOT. FAIL) THEN C CBAR = 2.8D0 C GOTO 40 C ENDIF C TTVAL = TVAL C CBAR = DSQRT( 3.0D0 * TTVAL ) C 40 CONTINUE C SUMA=0.0 DO 50 I=1,3 50 SUMA=SUMA+A(I,I) C C WRITE(6,11) C WRITE(6,12) ((A(I,J),J=1,3),I=1,3) C C EIGENVALUES AND EIGENVECTORS C CALL EIGENZ(A,R,EV,WS,3,3,0) C SUMS=0.0 NEG=1 DO 60 I=1,3 IF(EV(I).LT.0.0) NEG=-1 SUMS=SUMS+EV(I) S(I)=DSQRT(DABS(EV(4-I))) DO 60 J=1,3 60 A(I,J)=R(9-I*3+J) C C WRITE(6,13) (EV(4-I),(A(I,J),J=1,3),I=1,3) C DO 65 I=1,3 DO 65 J=1,3 B(I,J)=0.0 DO 65 K=1,3 65 B(I,J)=B(I,J)+A(K,I)*A(K,J) C PAI=3.1415926535898 DO 70 I=1,3 DO 70 J=1,3 70 A(I,J)=DARCOS(A(I,J))*180.0D0/PAI C IF(NEG.EQ.-1) GO TO 100 C C WRITE(6,14) C WRITE(6,17) C WRITE(6,21) C IF (.NOT. FAIL) THEN C WRITE(6,22) CBAR C ELSE C WRITE(6,23) CBAR C ENDIF C WRITE(6,15) (XYZ(I),S(I),XYZ(I),(A(I,J),J=1,3),I=1,3) IF( PROB .EQ. 0.0E0 ) GOTO 115 C C CALCULATING THE INCREASED CONFIDENCE REGION C DO 90 I = 1,3 90 S(I) = CBAR * S(I) C WRITE(6,16) ((B(I,J),J=1,3),I=1,3) C GO TO 110 C C 110 WRITE(6,19) PROB WRITE(6,25) FRESTA(IFRE),(XYZ(I),S(I),XYZ(I),(A(I,J),J=1,3),XYZ(I) &,SDEV(I),XYZ(I),SDEVE(I),I=1,3) GOTO 115 100 WRITE(6,18) C 115 RETURN C 11 FORMAT(1H1,//, T35,'*********************************',/, C * T35,'* *',/, C * T35,'* ANALYSIS OF ERROR-ELLIPSOID *',/, C * T35,'* *',/, C * T35,'*********************************',///) C 12 FORMAT(1H0, 5X, '(1) READ-IN COVARIANCE MATRIX:',//,(T10,3D16.6)) C 13 FORMAT(1H0,/,6X,'(2) EIGENVALUES AND EIGENVECTORS OF COVARIANCE ', C * 'MATRIX:',// ,T14,'EIGENVALUES',T41,'EIGENVECTORS (DIRECTI', C * 'ON COSINES)',//,(T9, F15.6,T34,F15.6,2F12.6)) C 14 FORMAT(1H0,/, 6X,'(3) PARAMETERS OF STANDARD ERROR-ELLIPSOID (19', C * '.9% CONFIDENCE REGION):') C 15 FORMAT(1H0,T41,'SPATIAL ANGLES (IN DEGREES) BETWEEN AXES',//, C * T14,'SEMI-MAJOR AXES', T53,'X(MODEL)', 2X,'Y(MODEL)', 2X, C * 'Z(MODEL)',//,(T14,A1,' =',F12.6,T41,A1,'(ELLIPSOID)', C * F7.2,2F10.2)) C 16 FORMAT(1H0,/, 6X,'(4) ORTHOGONALITY TEST',11H (M'M = I):,//, C * (T10,3F16.6)) C 17 FORMAT(1H+,T25,'________') 18 FORMAT(1H0,//,6X,'***(ERROR)*** NEGATIVE EIGENVALUE(S) DETECTED', * //,20X,'CHECK INPUT DATA FOR ROUND-OFF ERRORS') C 19 FORMAT (1H0,/,6X,'(5) PARAMETERS OF THE ERROR-ELLIPSOID', C * ' INCREASED TO A PROBABILITY OF',2X,F5.3,'% CONFIDENCE', C * ' REGION:') C 20 FORMAT ( 5(/),T30,'ERROR IN MDFI INPUT PARAMETERS. IFLAG IS',I5) 25 FORMAT(' '/1X,A8,(T14,A1,' =',F12.6,13X,A1,'(ELLIPSOID)',F7.2,2(F1 &0.2),11X,'(SIGMA)',A1,' = ',F6.4,3X,'(@ 95%)',A1,' = ',F6.4)) END SUBROUTINE DEGREE(RAD,IDEG,IMIN,SEC) C C SUBROUTINE TO CONVERT AN ANGLE EXPRESSED IN RADIANS TO AN ANGLE C EXPRESSED IN DEGREES, MINUTES, AND SECONDS. C C ARGUMENT DEFINITIONS -- C INPUT ARGUMENTS C RAD - ANGLE EXPRESSED IN RADIANS C OUTPUT ARGUMENTS C IDEG - DEGREES VALUE OF ANGLE C IMIN - MINUTES VALUE OF ANGLE C SEC - SECONDS VALUE OF ANGLE C IMPLICIT REAL*8 (A-H,M-Z) C C CONVERT THE ANGLE TO DECIMAL DEGREES DEG = (RAD * 180.0D0) / DARCOS(-1.0D0) C C CONVERT DECIMAL DEGREES TO DEGREES, MINUTES, AND SECONDS IDEG = IDINT(DEG) D = (DEG - DFLOAT(IDEG)) MIN = (DABS(D * 60.0D0)) IMIN = IDINT(MIN) SEC = (MIN - DFLOAT(IMIN)) * 60.0D0 C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: CONF2D * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * TWO DIMENSIONAL ABSOLUTE AND RELATIVE CONFIDENCE ELLIPSES FOR A * C * STATED SIGNIFICANCE LEVEL. * C * * C * * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL CONF3D(ALPHA,CX,DF,FAIL,NUKN) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 12 30 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * CX R*8 VARIANCE-COVARIANCE MATRIX OF ADJUSTED PARAMETERS * C * DF R*8 NUMBER OF DEGREES OF FREEDOM * C * NUKN I*2 NUMBER OF UNKNOWN PARAMETERS * C * ALPHA R*8 SIGNIFICANCE LEVEL * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * * C * * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * S R*8 VECTOR CONTAINING EIGENVALUES (SEMIMAJOR AXES) * C * CXEX R*8 EXTRACTED PORTION OF VARIANCE-COVARIANCE MATRIX OF * C * ADJUSTED PARAMETERS * C * * C ********************************************************************** C ********************************************************************** C C @PROCESS DC(RNBLK) SUBROUTINE CONF2D(ALPHA,IDF,FAIL,FRESTA,NUKN,NFR,NUK,APOST) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 FRESTA(NFR) CHARACTER*9 XYZP(3)/'X-Y PLANE','X-Z PLANE','Y-Z PLANE'/ REAL*8 DCHIIN,DFIN LOGICAL FAIL C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) DIMENSION CXEX(3,3) C C EXTRACT APPROPRIATE PORTION OF VARIANCE-COVARIANCE MATRIX IFRE = 1 PI = DARCOS(-1.D0) C C CALCULATING THE EXPANDING FACTOR PROB = 1.D0 - ALPHA D1 = 2.D0 D2 = DFLOAT(IDF) IF (.NOT. FAIL) THEN TVAL = DCHIIN(PROB,D1) XFAC = DSQRT(TVAL) ELSE TVAL = DFIN(PROB,D1,D2) XFAC = DSQRT(2.D0 * TVAL) ENDIF WRITE (6,1000) IF (.NOT. FAIL) THEN WRITE (6,1001) XFAC IF (APOST .NE. 0.D0) WRITE (6,1002) APOST ELSE WRITE (6,1003) XFAC WRITE (6,1004) APOST ENDIF WRITE (6,1005) WRITE (6,1006) DO 5 I = 1,NUKN,3 I1 = 1 J1 = 1 K1 = I + 2 DO 10 J = I,K1 DO 15 K = I,K1 CXEX(I1,J1) = RN(J,K) J1 = J1 + 1 15 CONTINUE J1 = 1 I1 = I1 + 1 10 CONTINUE C C COMPUTE CONFIDENCE ELLIPSES IN X-Y PLANE CXX = CXEX(1,1) CXY = CXEX(1,2) CYY = CXEX(2,2) CALL ELIPS(CXX,CXY,CYY,XFAC,AXY,BXY,AZIXY) C AXY1 = AXY * 1000.D0 C BXY1 = BXY * 1000.D0 THETA = PI - AZIXY C WRITE(14,1009) AXY,BXY,THETA CALL DEGREE(AZIXY,IDEG,IMIN,SEC) ISEC = SEC WRITE (6,1007) XYZP(1),FRESTA(IFRE),AXY,BXY,IDEG,IMIN,ISEC C C COMPUTE CONFIDENCE ELLIPSES IN X-Z PLANE CXZ = CXEX(1,3) CZZ = CXEX(3,3) CALL ELIPS(CXX,CXZ,CZZ,XFAC,AXZ,BXZ,AZIXZ) CALL DEGREE(AZIXZ,IDEG,IMIN,SEC) ISEC = SEC WRITE (6,1007) XYZP(2),FRESTA(IFRE),AXZ,BXZ,IDEG,IMIN,ISEC C C COMPUTE CONFIDENCE ELLIPSES IN Y-Z PLANE CYZ = CXEX(2,3) CALL ELIPS(CYY,CYZ,CZZ,XFAC,AYZ,BYZ,AZIYZ) CALL DEGREE(AZIYZ,IDEG,IMIN,SEC) ISEC = SEC WRITE (6,1007) XYZP(3),FRESTA(IFRE),AYZ,BYZ,IDEG,IMIN,ISEC WRITE (6,1008) IFRE = IFRE + 1 5 CONTINUE C C INPUT/FORMAT STATEMENTS 1000 FORMAT('1',80('*')/1X,'*',78X,'*'/1X,'*',8X,'SUMMARY OF ABSOLUTE 9 &5.0% CONFIDENCE ELLIPSES (OUT-OF-CONTEXT)',8X,'*'/1X,'*',78X,'*'/1 &X,80('*')//) 1001 FORMAT('0','FACTOR USED FOR OBTAINING THESE ELLIPSES FROM STANDARD & ELLIPSES: (VARIANCE FACTOR KNOWN) = ',F6.4/) 1002 FORMAT('0','COVARIANCE MATRIX OF PARAMETERS NOT MULTIPLIED BY A PO &STERIORI VARIANCE FACTOR (',F12.6,')'/T34,3('-')/) 1003 FORMAT('0','FACTOR USED FOR OBTAINING THESE ELLIPSES FROM STANDARD & ELLIPSES: (VARIANCE FACTOR UNKNOWN) = ',F6.4/) 1004 FORMAT('0','COVARIANCE MATRIX OF PARAMETERS IS MULTIPLIED BY A POS &TERIORI VARIANCE FACTOR (',F12.6,')'/T34,2('-')/) 1005 FORMAT('0','(SEMI-MAJOR AND SEMI-MINOR AXES IN MILLIMETRES)'/) 1006 FORMAT('0','CROSS-SECTION',T21,'POINT',T33,'SEMI-MAJOR AXIS',T53,' &SEMI-MINOR AXIS',T73,'AZIMUTH'/) 1007 FORMAT('0',T4,A9,T20,A8,T37,F7.4,T57,F7.4,T72,I3,1X,I2,1X,I2) 1008 FORMAT(' ',80('-')) 1009 FORMAT(1X,F9.4,2X,F9.4,2X,F11.6) C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: ELIPS * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * PARAMETERS OF THE TWO DIMENSIONAL CONFIDENCE ELLIPSE. * C * * C * * C * * C * SUBROUTINES REQUIRED: EIGEN3 (UNBCC-SELIB) * C * * C * FUNCTION SUBPROGRAMS REQUIRED: NONE * C * * C * ALGORITHM/METHOD: N/A * C * * C * USAGE: CALL ELIPS(CXX,CXY,CYY,XFAC,A,B,AZIM) * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 12 30 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * CX R*8 VARIANCE-COVARIANCE MATRIX OF ADJUSTED PARAMETERS * C * DF R*8 NUMBER OF DEGREES OF FREEDOM * C * NUKN I*2 NUMBER OF UNKNOWN PARAMETERS * C * ALPHA R*8 SIGNIFICANCE LEVEL * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * * C * * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * S R*8 VECTOR CONTAINING EIGENVALUES (SEMIMAJOR AXES) * C * CXEX R*8 EXTRACTED PORTION OF VARIANCE-COVARIANCE MATRIX OF * C * ADJUSTED PARAMETERS * C * * C ********************************************************************** C ********************************************************************** C C SUBROUTINE ELIPS(CXX,CXY,CYY,XFAC,A,B,THETA) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) C C COMPUTE PROGRAM DATA CONSTANTS PI = DARCOS(-1.D0) C C COMPUTE SEMI-MAJOR AND SEMI-MINOR AXES OF CONFIDENCE ELLIPSE ARG1 = (CXX + CYY) / 2.D0 ARG2 = DSQRT((CXX - CYY)**2 / 4.D0 + CXY**2) AMAX = (ARG1 + ARG2) A = DSQRT(AMAX) * XFAC * 1000.D0 B = DSQRT(ARG1 - ARG2) * XFAC * 1000.D0 C C COMPUTE AZIMUTH OF SEMI-MAJOR AXIS OF ELLIPSE ARG1 = AMAX - CXX ARG2 = DSQRT(CXX**2 + ARG1**2) THETA = DARCOS(ARG1 / ARG2) THETA = DSIGN(THETA,CXY) IF (THETA .LT. 0.D0) THETA = THETA + 2.D0 * PI C RETURN END @PROCESS DC(RNBLK) SUBROUTINE SCALAR(K,L,S) C C SCALAR MULTIPLICATION OF AN M X N MATRIX C C ARGUMENT DEFINITIONS -- C INPUT ARGUMENTS C X - M X N MATRIX TO BE MULTIPLIED C K - NUMBER OF ROWS OF MATRIX C L - NUMBER OF COLUMNS OF MATRIX C S - SCALAR C OUTPUT ARGUMENTS C Z - RESULTING MATRIX AFTER SCALAR MULTIPLICATION C IMPLICIT REAL*8(A-H,O-Z) C C ASSIGN MATRIX DIMENSIONS COMMON /RNBLK/ RN(1000,1000) C C PERFORM MATRIX SCALAR MULTIPLICATION DO 10 I = 1, K DO 20 J = 1, L RN(I,J) = S * RN(I,J) 20 CONTINUE 10 CONTINUE C RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: CHITST - CHI SQUARE TEST ON APRIOR VARIANCE FACTOR * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL - FORTRAN 77) * C * * C * DESCRIPTION: THE PURPOSE OF THIS PROGRAM IS TO TEST THE APRIORI * C * AND APOSTERIORI VARIANCE FACTORS FOR COMPATIBILITY AT THE GIVEN * C * CONFIDENCE LEVEL. * C * * C * FUNCTION CALLS: CHISQ(UNBCC-SELIB). * C * * C * ALGORITHM: N/A * C * * C * COMMENTS: NONE * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1987 11 15 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * IDF I*4 NUMBER OF DEGREES OF FREEDOM * C * ALPHA R*8 CHOSEN STATISTICAL SIGNIFICANCE LEVEL * C * APOST R*8 COMPUTED APOSTERIORI VARIANCE FACTOR * C * APRIOR R*8 APRIORI VARIANCE FACTOR * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * FAC R*8 FACTOR USED TO SCALE VARIANCE-COVARIANCE MATRIX OF * C * ADJUSTED PARAMETERS DEPENDING ON OUTCOME OF CHI * C * SQUARE TEST ON APRIORI VARIANCE FACTOR * C * FAIL L*0 PROGRAM CONTROL PARAMETER DEPENDING ON OUTCOME OF * C * CHI SQUARE TEST ON APRIORI VARIANCE FACTOR * C * * C ********************************************************************** C ********************************************************************** C C SUBROUTINE CHITST(ALPHA,APOST,APRIOR,IDF,FAC,FAIL) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z), INTEGER*4(I-N) LOGICAL FAIL REAL*8 DCHIIN C C PRINT OUT A POSTERIORI VARIANCE FACTOR WRITE(6,1000) APOST C C COMPUTE CONFIDENCE LEVEL CONF = 100.D0 * (1.D0 - ALPHA) C C COMPUTE ABSCISSA VALUES D1 = DFLOAT(IDF) ALPH1 = 1.D0 - (ALPHA / 2.D0) XCHI1 = DCHIIN(ALPH1,D1) ALPH2 = ALPHA / 2.D0 XCHI2 = DCHIIN(ALPH2,D1) XL = (DFLOAT(IDF) * APOST) / XCHI1 XR = (DFLOAT(IDF) * APOST) / XCHI2 C C DETERMINE IF APRIORI AND APOSTERIORI VARIANCE FACTORS ARE COMPATIBLE WRITE(6,1001) WRITE(6,1002) XL,APRIOR,XR IF (XL .LT. APRIOR .AND. XR .GT. APRIOR) THEN FAC = APRIOR FAIL = .FALSE. WRITE(6,1003) CONF WRITE(6,1004) ELSE FAC = APOST FAIL = .TRUE. WRITE(6,1005) CONF WRITE(6,1006) ENDIF C C I/O FORMAT STATEMENTS 1000 FORMAT('0',T19,'A POSTERIORI VARIANCE FACTOR = ',F12.6) 1001 FORMAT('0',T17,'CHI-SQUARE TEST ON THE A PRIORI VARIANCE FACTOR'/T &17,47('-')) 1002 FORMAT('0',T17,F17.8,' < ',F8.6,' < ',F17.8/) 1003 FORMAT('0','CHI-SQUARE TEST ON A PRIORI VARIANCE FACTOR PASSED AT &',F6.3,' % CONFIDENCE LEVEL'/T46,6('-')) 1004 FORMAT('0',T25,'(A PRIORI VARIANCE FACTOR KNOWN)') 1005 FORMAT('0','CHI-SQUARE TEST ON A PRIORI VARIANCE FACTOR FAILED AT &',F6.3,' % CONFIDENCE LEVEL'/T46,6('-')) 1006 FORMAT('0',T24,'(A PRIORI VARIANCE FACTOR UNKNOWN)') C RETURN END SUBROUTINE EIGENZ(T,R,VA,A,N,NX,NV) C IMPLICIT REAL *8(A-H,O-Z) DIMENSION T(NX,NX),VA(NX),A(6),R(9) C K=0 DO 10 J=1,N DO 10 I=1,J K=K+1 A(K)=T(I,J) 10 CONTINUE NN=N*(N+1)/2 20 IF (NV-1) 30,60,30 30 IQ=-N DO 50 J=1,N IQ=IQ+N DO 50 I=1,N IJ=IQ+I R(IJ)=0.0D0 IF (I-J) 50,40,50 40 R(IJ)=1.0D0 50 CONTINUE C C COMPUTE INITIAL AND FINAL NORMS ANORM AND ANORMX C 60 ANORM=0.0D0 DO 80 I=1,N DO 80 J=I,N IF (I-J) 70,80,70 70 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 80 CONTINUE IF (ANORM) 370,370,90 90 ANORM=DSQRT(2.0D2*ANORM) ANRMX=ANORM*1.0D-06/DFLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOULD, THR C IND=0 THR=ANORM 100 THR=THR/DFLOAT(N) 110 L=1 120 M=L+1 C C COMPUTE SIN AND COS C 130 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 140 IF (DABS(A(LM))-THR) 300,150,150 150 IND=1 LL=L+LQ MM=M+MQ X=0.5D0*(A(LL)-A(MM)) 160 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) IF (X) 170,180,180 170 Y=-Y 180 SINX=Y/DSQRT(2.0D0*(1.0D0+(DSQRT(1.0D0-Y*Y)))) SINX2=SINX*SINX 190 COSX=1.0D0-SINX2 COSX=DSQRT(COSX) COSX2=COSX*COSX SINCS=SINX*COSX C C ROTATE L AND M COLUMNS C ILQ=N*(L-1) IMQ=N*(M-1) DO 290 I=1,N IQ=(I*I-I)/2 IF (I-L) 200,270,200 200 IF (I-M) 210,270,220 210 IM=I+MQ GO TO 230 220 IM=M+IQ 230 IF (I-L) 240,250,250 240 IL=I+LQ GO TO 260 250 IL=L+IQ 260 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 270 IF (NV-1) 280,290,280 280 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 290 CONTINUE X=2.0D0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X C C TESTS FOR COMPLITION C C TEST FOR M LAST COLUMN C 300 IF (M-N) 310,320,310 310 M=M+1 GO TO 130 C C TEST FOR L SECOND FROM LAST COLUMN C 320 IF (L-(N-1)) 330,340,330 330 L=L+1 GO TO 120 340 IF (IND-1) 360,350,360 350 IND=0 GO TO 110 C C COMPARE THRESHOULD WITH FINAL NORM C 360 IF (THR-ANRMX) 370,370,100 C C SORT EIGENVALUES AND EIGENVECTORS C 370 IQ=-N DO 410 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 410 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF (A(LL)-A(MM)) 380,410,410 380 X=A(LL) A(LL)=A(MM) A(MM)=X IF (NV-1) 390,410,390 390 DO 400 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 400 R(IMR)=X 410 CONTINUE K=0 DO 420 I=1,N K=K+I 420 VA(I)=A(K) RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: TAURE * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C * CRITICAL VALUE FOR REJECTION OF STANDARDIZED RESIDUALS WITH * C * CONTROL OF TYPE I ERROR. * C * * C * SUBROUTINES REQUIRED: N/A * C * * C * FUNCTION SUBPROGRAMS REQUIRED: N/A * C * * C * ALGORITHM/METHOD: POPE'S APPROXIMATION. * C * REFERENCE - THE STATISTICS OF RESIDUALS AND THE DETECTION OF * C * OUTLIERS, A.J. POPE (1976). * C * U.S. DEPARTMENT OF COMMERCE * C * NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION * C * NATIONAL OCEAN SURVEY * C * NATIONAL GEODETIC SURVEY * C * NOAA TECHNICAL REPORT NO. 65 NGS1 * C * * C * USAGE: CALL TAURE(NT,NU,ALPH,CRTAU) * C * * C * COMMENTS: N/A * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1988 03 11 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * INPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * NT I*2 NUMBER OF OBSERVATIONS * C * NU I*2 NUMBER OF DEGREES OF FREEDOM * C * ALPH R*8 DESIRED PROBABILITY OF TYPE I ERROR * C * * C ********************************************************************** C * * C * OUTPUT VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * CRTAU R*8 CRITICAL VALUE (TAU-MAX) * C * * C ********************************************************************** C * * C * LOCAL VARIABLE NAMES : DEFINITION * C * * C * NAME TYPE DESCRIPTION * C * ------ ---- -------------------------------------------------- * C * * C ********************************************************************** C ********************************************************************** C C SUBROUTINE TAURE(NT,NU,ALPH,CRTAU) C C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) DATA PI/ 3.1415926535898 / PD = 2.D0 / PI S = 1.D0 WNU = NU U = WNU - 1.D0 IF (U .EQ. 0.D0) GO TO 72 IF (ALPH .EQ. 0.D0) GO TO 72 IF (ALPH .LT. 1.D0) GO TO 10 CRTAU = 0.D0 C RETURN C 10 Q = NT IF ( ALPH.GT.0.5 ) GO TO 19 AA = ALPH / Q DELT = AA DO 18 I = 1,100 XI = I DELT = DELT * ALPH * (( XI*Q - 1.)/(( XI+1.)*Q)) IF ( DELT.LT.1.D-14 ) GO TO 20 18 AA = AA + DELT 19 AA = 1. - (1.-ALPH)**(1./Q) 20 P = 1. - AA IF(U.EQ.1. ) GO TO 71 F = 1.3862943611199 - 2.*DLOG(AA) G = DSQRT(F) X = G - (2.515517 + 0.802853*G + 0.010328*F) $ / (1. + 1.432788*G + F*(0.189269 + 0.001308*G)) Y = X*X A = X*(1. + Y)/4. B = X*(3. + Y*(16. + 5.*Y))/96. C = X*(-15. + Y*(17. + Y*(19. + 3.*Y)))/384. E = X*(-945. + Y*(-1920. + Y*(1482. + Y*(776. + 79.*Y))))/92160. V = 1./U T = X + V*(A + V*(B + V*(C + E*V))) S = T/DSQRT(U + T*T) UM = U - 1. DELL = 1. DO 70 M = 1,50 H = 1. - S*S R = 0.0 IF ( DMOD(U,2.D0).EQ.0.0 ) GO TO 49 DD = DSQRT(H) N = 0.5*UM DO 45 I = 1,N Z = 2*I R = R + DD D = DD 45 DD = DD * H * (Z/(Z+1.)) R = PD*(R*S + DARSIN(S)) D = PD*D*UM GO TO 61 49 DD = 1. N = 0.5*U DO 55 I = 1,N Z = 2*I R = R + DD D = DD 55 DD = DD*H*((Z-1.)/Z) R = R*S D = D*UM 61 CONTINUE DEL = (P-R)/D IF( DABS( DEL/DELL ) .GT.0.5) GO TO 72 DELL = DEL S = S + DEL IF( DABS(DEL) .LT. 1.D-8 ) GO TO 72 70 CONTINUE GO TO 72 71 S =DSIN(P/PD) 72 CRTAU = S*DSQRT(WNU) RETURN END C ********************************************************************** C ********************************************************************** C * * C * PROGRAM NAME: OUTLIER * C * * C * VERSION: 1.0.0 * C * * C * CATEGORY: SUBROUTINE * C * * C * COMPUTER: IBM 3090 PROCESSOR COMPLEX * C * * C * LANGUAGE: IBM VS FORTRAN VERSION 2 (1987 06) (LEVEL 2.0.0) * C * * C * DESCRIPTION: THE PURPOSE OF THIS SUBROUTINE IS TO PERFORM THE * C * OUT-OF-CONTEXT AND IN-CONTEXT TESTS FOR RESIDUAL OUTLIERS. THE * C * APPROPRIATE TESTS ARE SELECTED DEPENDING ON IF THE A PRIORI * C * REFERENCE VARIANCE IS KNOWN OR UNKNOWN. * C * * C * * C * * C * SUBROUTINES REQUIRED: N/A * C * * C * FUNCTION SUBPROGRAMS REQUIRED: DNORIN, TAURE * C * * C * ALGORITHM/METHOD: * C * * C * USAGE: CALL OUTLIER( * C * * C * COMMENTS: * C * * C * HISTORY: * C * * C * DATE PROGRAMMER MODIFICATION * C * ---------- -------------------- -------------------------------- * C * 1988 06 23 MARK W. ROHDE TYPED IN PROGRAM SOURCE CODE * C * * C * * C * * C * * C ********************************************************************** C ********************************************************************** C * * C * VARIABLES : DEFINITION * C * * C * * C * INPUT ARGUMENTS * C * * C * DMS - REAL*8 VALUE OF ANGLE EXPRESSED IN DDD.MMSSS FORMAT * C * * C * * C * OUTPUT ARGUMENTS * C * * C * DEC - REAL*8 VALUE OF ANGLE EXPRESSE IN DECIMAL DEGREE FORMAT * C * * C * * C * LOCAL VARIABLES * C * * C * TEMP - REAL*8 VALUE USED TEMPORARILY IN CONVERSION COMPUTATION * C * TEMP1 - REAL*8 VALUE USED TEMPORARILY IN CONVERSION COMPUTATION * C * * C ********************************************************************** C ********************************************************************** C C SUBROUTINE OUTLIER(RHATS,STA,ALPHA,APRIOR,IDF,REFVAR,NOB,NOBS,OBCO &DE,OBS,RHAT) C C ASSIGN VARIABLE TYPE DECLARATIONS IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 STA(NOB,3) INTEGER*4 OBCODE(NOB) LOGICAL CTEXT,FNDOUT,REFVAR REAL*8 DNORIN C C ASSIGN MATRIX DIMENSIONS DIMENSION OBS(NOB,4),RHAT(NOB),RHATS(NOB) C C PERFORM OUT-OF-CONTEXT TESTING FOR RESIDUAL OUTLIERS CTEXT = .FALSE. CONF = 100.D0 - ALPHA WRITE(6,2000) CONF C C DETERMINE IF A PRIORI REFERENCE VARIANCE IS KNOWN OR UNKNOWN ALPH1 = 1.D0 - (ALPHA / 200.D0) 5000 FNDOUT = .FALSE. WRITE(6,2001) IF (REFVAR) THEN J = 1 NREJ = 0 DO 10 I = 1,NOBS CRNOR = DNORIN(ALPH1) IF (DABS(RHATS(I)) .GE. CRNOR) THEN FNDOUT = .TRUE. NREJ = NREJ + 1 IF (OBCODE(I) .EQ. 1) WRITE(6,2002) I IF (OBCODE(I) .EQ. 2) WRITE(6,2003) I,J IF (OBCODE(I) .EQ. 4) WRITE(6,2004) I IF (OBCODE(I) .EQ. 5) WRITE(6,2005) I WRITE(6,2006) STA(I,1),STA(I,2),STA(I,3),RHATS(I), & OBS(I,1),CRNOR ENDIF IF (OBCODE(I) .NE. -2) THEN J = J + 1 ELSE J = 1 ENDIF 10 CONTINUE IF (FNDOUT) THEN IPCT = (NREJ * 100) / NOBS WRITE(6,2007) NREJ,IPCT WRITE(6,2008) ELSE WRITE(6,2009) ENDIF ELSE J = 1 NREJ = 0 DO 20 I = 1,NOBS IF (CTEXT) THEN NT = NOBS ELSE NT = 1 ENDIF ALPH = ALPHA / 100.D0 CALL TAURE(NT,IDF,ALPH,CRTAU) IF (DABS(RHATS(I)) .GE. CRTAU) THEN FNDOUT = .TRUE. NREJ = NREJ + 1 IF (OBCODE(I) .EQ. 1) WRITE(6,2002) I IF (OBCODE(I) .EQ. 2) WRITE(6,2003) I,J IF (OBCODE(I) .EQ. 4) WRITE(6,2004) I IF (OBCODE(I) .EQ. 5) WRITE(6,2005) I WRITE(6,2006) STA(I,1),STA(I,2),STA(I,3),RHATS(I), & OBS(I,1),CRTAU ENDIF IF (OBCODE(I) .NE. -2) THEN J = J + 1 ELSE J = 1 ENDIF 20 CONTINUE IF (FNDOUT) THEN IPCT = (NREJ * 100) / NOBS WRITE(6,2007) NREJ,IPCT WRITE(6,2008) ELSE WRITE(6,2009) ENDIF ENDIF IF (CTEXT) GOTO 5001 C C PERFORM IN-CONTEXT TESTING FOR RESIDUAL OUTLIERS WRITE(6,2010) CONF ALPH1 = 1.D0 - (ALPHA / (200.D0 * DFLOAT(NOBS))) CTEXT = .TRUE. GOTO 5000 C C INPUT/OUTPUT FORMAT STATEMENTS 2000 FORMAT('1',18X,'OUT-OF-CONTEXT TESTING FOR RESIDUAL OUTLIERS'/19X, &44('-')//29X,F6.3,'% CONFIDENCE LEVEL'//) 2001 FORMAT('0',T58,'STD.DEV.'/T7,'OBSERVATION',T22,'AT',T31,'FROM',T40 &,'TO',T49,'RESIDUAL',T58,'RESIDUAL',T67,'CRITICAL VALUE') 2002 FORMAT('0',I4,T7,'SLOPE DIST.') 2003 FORMAT('0',I4,T7,'DIRECTION',T19,I2) 2004 FORMAT('0',I4,T7,'AZIMUTH') 2005 FORMAT('0',I4,T7,'ZENITH ANGLE',T19,I2) 2006 FORMAT('+',T22,A8,T31,A8,T40,A8,T49,F8.4,T59,F6.4,T70,F7.4) 2007 FORMAT('0',T6,I4,' RESIDUAL(S) (',I3,'% OF THE OBSERVATIONS) WERE &FLAGGED FOR REJECTION'//) 2008 FORMAT('0',35('*'),' WARNING ',36('*')/1X,'*',78X,'*'/1X,'* OBSERV &ATIONS CORRESPONDING TO REJECTED RESIDUALS WERE USED IN THE ADJUST &MENT *'/1X,'*',78X,'*'/1X,80('*')//) 2009 FORMAT('0',T19,'******** NO RESIDUAL OUTLIERS FOUND ********'//) 2010 FORMAT('0',20X,'IN-CONTEXT TESTING FOR RESIDUAL OUTLIERS'/21X,40(' &-')//29X,F6.3,'% CONFIDENCE LEVEL'//) C 5001 RETURN END