PROGRAM CORDM C C FROM SEPARATE CAMPAIGN ESTIMATIONS OF COORDINATES C UNDER MINIMAL CONSTRAINTS (SINGLE STATION FIXED + PSEUDO- C AZIMUTH + PSEUDO-DISTANCE) C C DISPLACEMENTS VERSUS 'RELATIVE' ELLIPSES AT SPECIFIED ALPHA LEVEL C - FROM CHOICE OF MINIMAL CONSTRAINTS C C MODELLING OF DEFORMATION WITH PARAMETERS C A0 B0 EX EY EXY OMEGA, K, A1 A2 A3 A4 A5 B1 B2 B3 B4 B5 C C C INPUT : C C TITLE FIRST LINE OF DATA TAKEN AS DESCRIPTION (CHARACTER) C ITNZ TOTAL NUMBER OF ZONES IN DEFORMATION MODEL C IDEFM NUMBER OF PARAMETERS IN DEFORMATION MODEL C IALPH " 1 " TESTING AND REGIONS TO LEVEL OF ALPHA C ALPHA LEVEL OF SIGNIFICANCE FOR TESTING C IPREA " 1 " PREANALYSIS FOR ACCURACY IN SELECTED PARAMETERS C C C FOR THE FIRST OF THE TWO CAMPAIGNS BEING COMPARED : C C EPOCH1 CAMPAIGN IDENTIFIER CHARACTER STRING ; C SIG1 THE ESTIMATED VARIANCE FACTOR OF THE C INDIVIDUAL CAMPAIGN ADJUSTMENT; C IDF1 THE DEGREES OF FREEDOM; C INS1 TOTAL NUMBER OF STATIONS; C INB1 NUMBER OF BLAHA STATIONS; C INP1 NUMBER OF WEIGHTED STATIONS; C INF1 NUMBER OF FIXED STATIONS; C C FORM1 5A5 FORMAT FOR STN1,E1,N1,IZ1 C C VECTORS OF : C STN1 STATION IDENTIFIER NAMES; C E1,N1 EASTING AND NORTHING OF STATION ; C IZ1 ZONE CODES . C CX1 VARIANCE-COVARIANCE MATRIX FOR THE E1,N1 FROM C INDIVIDUAL EPOCH ADJUSTMENT (UPPER TRIANGULAR BY ROWS) C C AND SIMILARLY FOR THE SECOND EPOCH C C C USES LIBRARY SUBROUTINES COPYD GDATE MADDD MMULD MOUTD MSUBD C TRNSD C USES IMSL FUNCTION SUBROUTINES C C C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) REAL*8 IAT,IFR,ITO,ISTNG,ITOA REAL*4 DF1,DF2,SCRV90,SCRV95,SCRV99,PP,X95,X05,X99,X01 @,SXT95,SXT99,SXTM95,SXTM99,SN95,SN99,SNM95,SNM99,SCRV LOGICAL*1 DATE(18),TIME(6),CODE(6) CHARACTER*5 FORM1(5),FORM2(5) CHARACTER*8 BLK/' '/ DIMENSION STN1(50),STN2(50),E1(50),E2(50),N1(50),N2(50) @,DX(90),CX1(90,90),CX2(90,90),CDX(90,90),B(90,50),E(50) @,IZ1(50),IZ2(50),DXT(1,90),DXTP(1,90),DISP(50),XG(50) @,AZ(50),CV(90,90),IAT(200),IFR(200),ITO(200),ISTNG(50),YG(50) DIMENSION PHI(50),CHPARA(50),IPC(50) @,IZZC(50),CE(50,50),SE(50),XY1(90),XY2(90),STNG(50) DIMENSION AE(50),BE(50),AZE(50),DE(50),DP(50),S(90,90),ST(90,90), @ TITLE(10),ICOLA(50),IZG1(50),IZG2(50),IZCG(50),IZGC(50) @,XY(90),SCX(90,90),AP(200,50),STNG1(50),STNG2(50),CXG1(90,90) @,CXG2(90,90),X(50),Y(50),EG(50),NG(50),IZPC(30,50),P(90,90) DIMENSION CEE(3,3),VDX(90),CVDX(90,90),ANAT(90,90),AN(90,50) @,AT(50,90),AX(90),ATP(50,90),N(50,50),U(50),VT(1,90),VTP(1,90) @,D(3,3),DT(3,3),C1(3,3),CIDT(3,3) IDRA=60 IDCB=50 IDCA=90 IDRC=IDCA IRV=50 ICB=50 IDB=90 IDAA=200 RTD=180.0D0/(4.0D0*DATAN(1.0D0)) IALPH=0 C CALL GDATE (DATE,TIME,CODE) WRITE(6,4) CODE(5),CODE(6),(CODE(I),I=1,4),(TIME(I),I=1,6) 4 FORMAT('1','19',3(2A1,1X),50X,2A1,':',2A1,':',2A1) READ(5,2) (TITLE(I),I=1,10) 2 FORMAT(10A8) WRITE(6,3) (TITLE(I),I=1,10) 3 FORMAT('-',10A8,//) C READ(5,*) ITNZ,IDEFM, IALPH,ALPHA,IPREA IF(IDEFM.GT.0) WRITE(6,53) ITNZ,IDEFM 53 FORMAT(' ',' DEFORMATION MODEL WITH ',I5,' ZONES AND ',I5, @ ' PARAMETERS'/) IF(IALPH.EQ.1) WRITE(6,54) ALPHA 54 FORMAT(' ',' TESTING AT ALPHA LEVEL OF ',F7.5/) ALPH1=0.95D0 IF(IALPH.EQ.1) ALPH1=1.0D0-ALPHA C IF(IDEFM.LE.0) GO TO 100 IDFM=IDEFM IF(IDEFM.GT.40) IDFM=40 READ(5,204) (IPC(I),I=1,IDFM) WRITE(6,205) (IPC(I),I=1,IDFM) DO 60 I=1,ITNZ READ(5,204) (IZPC(I,J),J=1,IDFM) WRITE(6,206) I,(IZPC(I,J),J=1,IDFM) 60 CONTINUE IF(IDEFM.LE.40) GO TO 75 WRITE(6,69) 69 FORMAT('-') READ(5,204) (IPC(I),I=41,IDEFM) WRITE(6,205) (IPC(I),I=41,IDEFM) DO 65 I=1,ITNZ READ(5,204) (IZPC(I,J),J=41,IDEFM) WRITE(6,206) I,(IZPC(I,J),J=41,IDEFM) 65 CONTINUE 75 CONTINUE 204 FORMAT(40I2) 205 FORMAT(6X,40I3) 206 FORMAT(2X,I2,2X,40I3) C 100 READ(5,1) EPOCH1,SIG1,IDF1,INS1,INB1,INP1,INF1 IECM1=0 1 FORMAT(A8,2X,F15.12,5I5) READ(5,9) (FORM1(I),I=1,5) 9 FORMAT(5A5) WRITE(6,5) EPOCH1 5 FORMAT('1',T6,'ADJUSTED COORDINATES FOR EPOCH ',A8//' ',T11, @'STATION',T31,'EASTING',T51,'NORTHING'/) INS=INS1-INF1 DO 10 I=1,INS READ(5,FORM1) STN1(I),E1(I),N1(I),IZ1(I) WRITE(6,15) STN1(I),E1(I),N1(I),IZ1(I) 10 CONTINUE 15 FORMAT('0',T11,A8,T31,F17.8,T51,F17.8,I4) WRITE(6,16) 16 FORMAT('0',T6,'FIXED STATION(S)') J=INS+1 JJ=INS+INF1 DO 11 I=J,JJ READ(5,FORM1) STN1(I),E1(I),N1(I),IZ1(I) WRITE(6,15) STN1(I),E1(I),N1(I),IZ1(I) 11 CONTINUE IRB=2*INS DO 31 I=1,IRB READ(5,*) (CX1(I,J),J=I,IRB) DO 9031 JJ=I,IRB CX1(JJ,I)=CX1(I,JJ) 9031 CONTINUE C 90 02 02 J=1 TO J=I PLUS DO 9031 LOOP 31 CONTINUE IF(IECM1.EQ.1) THEN DO 12 I=1,IRB DO 13 I2=1,IRB CX1(I,I2)=CX1(I,I2)/SIG1 13 CONTINUE 12 CONTINUE ENDIF WRITE(6,6) IDF1 6 FORMAT('-',T6,'DEGREES OF FREEDOM IN ADJUSTMENT :',I5/) WRITE(6,7) SIG1 7 FORMAT('-',T6,'ESTIMATED VARIANCE FACTOR :',F15.12///) IJJ=2*INS+1 IJ=2*INS1 DO 70 I=IJJ,IJ DO 70 II=1,IJ CX1(I,II)=0.0D0 CX1(II,I)=0.0D0 70 CONTINUE IGC=0 CALL SEGREG (STN1,E1,N1,IZ1,CX1,IDRC,IDRC/2,INS1,IGC,STNG1,IZG1 @,XY1,CXG1,INPG1,XY,S,ST,SCX) C WRITE(6,520) 520 FORMAT('1',/) C C READ(5,1) EPOCH2,SIG2,IDF2,INS2,INB2,INP2,INF2 IECM2=0 READ(5,9) (FORM2(I),I=1,5) WRITE(6,5) EPOCH2 INS=INS2-INF2 DO 30 I=1,INS READ(5,FORM2) STN2(I),E2(I),N2(I),IZ2(I) WRITE(6,15) STN2(I),E2(I),N2(I),IZ2(I) 30 CONTINUE WRITE(6,16) J=INS+1 JJ=INS+INF2 DO 32 I=J,JJ READ(5,FORM2) STN2(I),E2(I),N2(I),IZ2(I) WRITE(6,15) STN2(I),E2(I),N2(I),IZ2(I) 32 CONTINUE IRB=2*INS DO 21 I=1,IRB READ(5,*) (CX2(I,J),J=I,IRB) DO 9021 JJ=I,IRB CX2(JJ,I)=CX2(I,JJ) 9021 CONTINUE C 90 02 02 J=1 TO J=I, PLUS DO 9021 LOOP 21 CONTINUE IF(IECM2.EQ.1) THEN DO 22 I=1,IRB DO 23 I2=1,IRB CX2(I,I2)=CX2(I,I2)/SIG2 23 CONTINUE 22 CONTINUE ENDIF WRITE(6,6) IDF2 WRITE(6,7) SIG2 IJJ=2*INS+1 IJ=2*INS2 DO 90 I=IJJ,IJ DO 90 II=1,IJ CX2(I,II)=0.0D0 CX2(II,I)=0.0D0 90 CONTINUE IGC=0 CALL SEGREG (STN2,E2,N2,IZ2,CX2,IDRC,IDRC/2,INS2,IGC,STNG2,IZG2 @,XY2,CXG2,INPG2,XY,S,ST,SCX) C C ID=0 IF(INF2.EQ.INF1) ID=INF1+2 IF(INF2.EQ.INF1) INF=INF1 IRB=0 IF(INPG1.NE.INPG2) GO TO 9000 IRB=2*INPG1 INS=INPG1 DO 40 I=1,INPG1 IF(STNG1(I).NE.STNG2(I)) GO TO 41 STNG(I)=STNG1(I) IZCG(I)=IZG1(I) JJ=2*I J=JJ-1 EG(I)=XY1(J) NG(I)=XY1(JJ) GO TO 40 41 WRITE(6,45) GO TO 9000 40 CONTINUE 45 FORMAT('1 STATIONS NOT MATCHED IN SEGREGATION'//) IDF12=IDF1+IDF2 SIG12=(IDF1*SIG1+IDF2*SIG2)/(IDF12) WRITE(6,8) SIG12,IDF12 8 FORMAT('1',T6,'WEIGHTED VARIANCE FACTOR :',F15.12/' ', @ ' DEGREES OF FREEDOM :',I5) DO 33 II=1,IRB DX(II)=XY2(II)-XY1(II) DO 34 I=1,IRB CDX(I,II)=CXG1(I,II)+CXG2(I,II) 34 CONTINUE 33 CONTINUE C CALL MOUTD (CDX,IDB,IRB,IRB) C C DISPLACEMENTS REFERRED TO CHOICE OF MINIMAL CONSTRAINTS C CALL DAZ1 (DX,IDCA,IRB,DISP,IDCA/2,INS,AZ) CALL ELLAN (DX,IDCA,IRB,CDX,0,SIG12,IDF12,SNGL(ALPHA),0,0,IDCA/2 @ ,IJ,AZ,AE,BE,AZE,PHI,DE,DP) WRITE(6,610) ALPH1 610 FORMAT('1',T7,'STATION',T17,'DX (M)',T27,'DY (M)',T42,'A (M)',T52, @'B (M)',T62,'AZ A (DEG) AT 0.39 AND AT ',F6.3//) II=2*INS FAK=1.0D0 IF(IALPH.EQ.1) SCRV = CHIIN(SNGL(1.0D0-ALPHA),2.0) IF(IALPH.EQ.1) FAK=DSQRT(DBLE(SCRV)) DO 630 I=1,II,2 J=(I+1)/2 WRITE(6,620) STNG(J),DX(I),DX(I+1),AE(J),BE(J),AZE(J) WRITE(6,621) FAK*AE(J),FAK*BE(J) 630 CONTINUE 620 FORMAT('0',T6,A8,T17,F7.4,T27,F7.4,T42,F6.4,T52,F6.4,T62,F10.4) 621 FORMAT('0',T42,F6.4,T52,F6.4) WRITE(6,640) ALPH1 640 FORMAT('1',T6,'STATION',T17,'TOTAL',T47,'TO PERIPHERY OF'/' ',T17, @'DISPLACEMENT AZ (DEG)',T47,'STANDARD ELLIPSE ( ',F6.4, @' ) ELLIPSE'//) DO 650 I=1,INS WRITE(6,660) STNG(I),DISP(I),AZ(I),DE(I),FAK*DE(I) 650 CONTINUE 660 FORMAT('0',T6,A8,T17,F7.4,T35,F10.4,T49,F7.4,T72,F7.4) WRITE(6,850) 850 FORMAT('1'//) WRITE(6,863) INS,INS 863 FORMAT(' ',I3,' 0',' 0',I3,' 200.00',' 1.000',' 10.0', @ ' 1.000',' 20.0'/' ','(A8,2F15.4,I3)') DO 729 I=1,INS WRITE(6,862) STNG(I),EG(I),NG(I) 729 CONTINUE 862 FORMAT(' ',A8,2F15.4,' 4') WRITE(6,864) 864 FORMAT(' ','(A8,4F10.4,I4,2I3)') DO 730 I=1,II,2 J=(I+1)/2 CALL DEGDMS (AZE(J),IDE,IME,ISE) WRITE(6,865) STNG(J),DX(I),DX(I+1),FAK*AE(J),FAK*BE(J) @ ,IDE,IME,ISE 730 CONTINUE 865 FORMAT(' ',A8,4F10.4,I4,2I3) C C C IDEFM NE ZERO DEFORMATION PARAMETERS FROM COORDINATE DIFFERENCES C 900 IF(IDEFM.EQ.0) GO TO 9000 IF(IDEFM.LT.0) GO TO 933 CALL BMAT (EG,NG,IZCG,IDCB,INS,IZPC,IPC,B,IDCA,IDCB,IDEFM,CHPARA) C WRITE(6,9904) C9904 FORMAT('1','B') C CALL MOUTD (B,IDB,2*INS,IDEFM) IJ=INS-INF C WRITE (6,9907) C9907 FORMAT('1','CDX') C CALL MOUTD (CDX,IDB,2*INS,2*INS) CALL COPYD (P,IDB,CDX,IDB,2*INS,2*INS) CALL SPIN (P,2*IJ,IDB,DET,IDEXP) INOBS=2*IJ C WRITE(6,9905) DET,IDEXP C9905 FORMAT('1','P = INVERSE OF CDX',5X,F7.4,'D',I4) C CALL MOUTD (P,IDB,INOBS,INOBS) CALL PARADJ (IPREA,B,IDB,IDCB,INOBS,IDEFM,DX,P,E,CE,VDX,CVDX, @SIGE,ANAT,AN,AT,AX,ATP,N,U,VT,VTP) C WRITE(6,9901) C9901 FORMAT('1','A*N*AT') C CALL MOUTD (ANAT,IDB,INOBS,INOBS) C WRITE(6,9902) C9902 FORMAT('1','A*N') C CALL MOUTD (AN,IDB,INOBS,IDEFM) C WRITE(6,9903) C9903 FORMAT('1','AT*P') C CALL MOUTD (ATP,IDCB,IDEFM,INOBS) SIG=SIGE*(INOBS-IDEFM) IDFE=2*INS-(INF+2)-IDEFM SIGE=SIG/IDFE IF(IPREA.EQ.1) SIGE=1.0D0 WRITE(6,906) 906 FORMAT('1',T6,'ESTIMATED DEFORMATION COMPONENTS, UNSCALED STD. ', @'DEV. AND "(1-ALPHA)" LEVELS'/) DO 910 I=1,IDEFM PP=0.0000 IF(IPREA.NE.1) SE(I)=E(I)**2/(CE(I,I)*SIGE) IF(IPREA.NE.1) PP = FDF(SNGL(SE(I)),1.0,FLOAT(IDFE)) SDE=DSQRT(CE(I,I)) DO 908 J=1,ITNZ IF(IZPC(J,I).EQ.1) WRITE(6,907) J,CHPARA(I),E(I),SDE,PP 908 CONTINUE 910 CONTINUE 907 FORMAT('0',5X,I4,5X,A8,E17.10,2X,'+/-',E17.10,F15.4/) WRITE(6,915) SIGE 915 FORMAT('1',5X,'VARIANCE-COVARIANCE MATRIX FOR ESTIMATED DEFORMA', @'TION COMPONENTS ( NOT MULTIPLIED BY ',F15.10,' )') CALL MOUTD (CE,IDCB,IDEFM,IDEFM) WRITE(6,918) 918 FORMAT('1',5X,'DERIVED DEFORMATIONS, UNSCALED STD. DEV. AND', @' "(1-ALPHA)" LEVELS'/) DO 920 I=1,IDEFM IF(IPC(I).NE.1) GO TO 923 Q=(E(I)**2*CE(I+1,I+1)+E(I+1)**2*CE(I,I)-2.0D0*E(I)*E(I+1)*CE(I, @I+1))/((CE(I,I)*CE(I+1,I+1)-CE(I,I+1)**2)*2.0D0*SIGE) PP = FDF(SNGL(Q),2.0,FLOAT(IDFE)) DISS=DSQRT(E(I)**2+E(I+1)**2) IF(DISS.LE.0.0D0) THEN AZZ=0.0D0 GO TO 922 ENDIF A11=E(I)/DISS A12=E(I+1)/DISS A21=E(I+1)/DISS**2 A22=-E(I)/DISS**2 SDD=DSQRT(A11*A11*CE(I,I)+2.0D0*A11*A12*CE(I,I+1)+A12*A12*CE(I+1, @ I+1)) SDA=DSQRT(A21*A21*CE(I,I)+2.0D0*A21*A22*CE(I,I+1)+A22*A22*CE(I+1, @ I+1))*RTD AZZ=DATAN2(E(I),E(I+1)) IF(AZZ.LT.0.0D0) AZZ=AZZ+8.0D0*DATAN(1.0D0) AZZ=AZZ*RTD 922 IF(IPREA.NE.1) GOTO 919 CALL ELLIP (CE(I,I),CE(I,I+1),CE(I+1,I+1),AR,BR,AZR,PHIR) WRITE(6,9865) BLK,E(I),E(I+1),AR,BR,IFIX(SNGL(AZR)) 9865 FORMAT(' ',A8,4F10.4,I4) GOTO 920 919 DO 921 J=1,ITNZ IF(IZPC(J,I).EQ.1) WRITE(6,925) J,CHPARA(I),CHPARA(I+1), @ DISS,SDD,AZZ,SDA,PP 921 CONTINUE 923 IF(IPC(I).NE.3) GO TO 920 IF(IPC(I+1).NE.4) GO TO 920 IF(IPC(I+2).NE.5) GO TO 920 WRITE(6,929) 929 FORMAT('-',T42,'E MAX',T92,'AZ E MAX',T111,'(1-ALPHA)'/ @ '-',T42,'E MIN',T111,'(1-ALPHA)'/) CALL STRPR (E(I),E(I+1),E(I+2),CE(I,I),CE(I+1,I+1),CE(I+2,I+2), @ CE(I,I+1),CE(I,I+2),CE(I+1,I+2),EMAX,EMIN,PHIE,CEE,C1,D,DT,CIDT) Q=EMAX**2/(CEE(1,1)*SIGE) AZZE=90.0D0-PHIE*RTD PP = FDF(SNGL(Q),3.0,FLOAT(IDFE)) SDMAX=DSQRT(CEE(1,1)) SDMIN=DSQRT(CEE(2,2)) SDAZ=DSQRT(CEE(3,3))*RTD DO 924 J=1,ITNZ IF(IZPC(J,I).EQ.1) WRITE(6,926) J,CHPARA(I),CHPARA(I+1), @ CHPARA(I+2),EMAX,SDMAX,AZZE,SDAZ,PP 924 CONTINUE Q=EMIN**2/(CEE(2,2)*SIGE) PP = FDF(SNGL(Q),3.0,FLOAT(IDFE)) DO 928 J=1,ITNZ IF(IZPC(J,I).EQ.1) WRITE(6,927) J,CHPARA(I),CHPARA(I+1), @ CHPARA(I+2),EMIN,SDMIN,PP 928 CONTINUE 920 CONTINUE 925 FORMAT('0',5X,I4,5X,2A8,F15.10,' +/-',F15.10,F15.5,' +/-' @,F10.4,F15.4/) 926 FORMAT('0',5X,I4,5X,3A8,F15.10,' +/-',F15.10,F15.5,' +/-' @,F10.4,F15.4/) 927 FORMAT('0',5X,I4,5X,3A8,F15.10,' +/-',F15.10,30X,F15.4/) WRITE(6,6918) 6918 FORMAT('1',5X,'A0, B0, A(0.95), B(0.95), AZ(D,M,S)'/) SCRV95 = FIN(0.95,2.0,FLOAT(IDFE)) EFAK=DSQRT(SIGE)*DBLE(SQRT(SCRV95)) DO 6920 I=1,IDEFM IF(IPC(I).NE.1) GOTO 6920 CALL ELLIP (CE(I,I),CE(I,I+1),CE(I+1,I+1),AR,BR,AZR,PHIR) CALL DEGDMS (AZR,ID,IM,IS) WRITE(6,865) BLK,E(I),E(I+1),AR*EFAK,BR*EFAK,ID,IM,IS 6920 CONTINUE IF(IPREA.EQ.1) GOTO 9000 WRITE(6,691) SIGE,SIG12 C 690 FORMAT('1'/) 691 FORMAT('1 STATION, DISPLACEMENT COMPONENT RESIDUAL,', @ ' SCALED STANDARDIZED DCR'/ @ ' ','ESTIMATED VARIANCE FACTOR ',F15.10/ @ ' ','VARIANCE FACTOR A PRIORI ',F15.10/) DO 930 I=1,INOBS J=I/2 II=I IF(I.NE.2*(II/2)) J=J+1 VD=VDX(I)/DSQRT(CVDX(I,I)*SIGE) WRITE(6,710) STNG(J),VDX(I),VD 930 CONTINUE 710 FORMAT(' ',10X,A8,5X,F10.6,5X,F10.6) AL05=0.05D0 AL01=0.01D0 CALL TAURE (1,IDFE,AL05,XT95) CALL TAURE (1,IDFE,AL01,XT99) CALL TAURE (INOBS,IDFE,AL05,XTM95) CALL TAURE (INOBS,IDFE,AL01,XTM99) WRITE(6,932) XT95,XT99,XTM95,XTM99 932 FORMAT('-',T6,'TAU (0.05) : ',F6.4,10X,'TAU (0.01) : ',F6.4/'0', @T6,'TAU MAX (0.05) : ',F6.4,10X,'TAU MAX (0.01) : ',F6.4//) AL05N=1.0D0-AL05/2.0D0 AL01N=1.0D0-AL01/2.0D0 AL05M=1.0D0-AL05/(2.0D0*FLOAT(INOBS)) AL01M=1.0D0-AL01/(2.0D0*FLOAT(INOBS)) SXT95 = TIN(SNGL(AL05N),FLOAT(INOBS)) SXT99 = TIN(SNGL(AL01N),FLOAT(INOBS)) SXTM95 = TIN(SNGL(AL05M),FLOAT(INOBS)) SXTM99 = TIN(SNGL(AL01M),FLOAT(INOBS)) WRITE(6,942) SXT95,SXT99,SXTM95,SXTM99 942 FORMAT('-',T6,'T (0.05) : ',F6.4,10X,'T (0.01) : ',F6.4/'0', @T6,'T MAX (0.05) : ',F6.4,10X,'T MAX (0.01) : ',F6.4//) SN95 = ANORIN(SNGL(AL05N)) SN99 = ANORIN(SNGL(AL01N)) SNM95 = ANORIN(SNGL(AL05M)) SNM99 = ANORIN(SNGL(AL01M)) WRITE(6,952) SN95,SN99,SNM95,SNM99 952 FORMAT('-',T6,'NORM(0.05) : ',F6.4,10X,'NORM(0.01) : ',F6.4/'0', @T6,'NORM MAX(0.05) : ',F6.4,10X,'NORM MAX(0.01) : ',F6.4//) GO TO 934 933 CALL TRNSD (VT,1,DX,IDB,IRB,1) IRB=IRB-2*INF CALL SPIN (CDX,IRB,IDB,DET,IDEXP) CALL MMULD (VTP,1,VT,1,CDX,IDB,1,IRB,IRB) CALL MMULD (SIG,1,VTP,1,DX,IDB,1,IRB,1) ICON=1 IF(INF.EQ.2) ICON=0 IDFE=IRB-ICON SIGE=SIG/ IDFE 934 WRITE(6,935) SIG 935 FORMAT('- THE QUADRATIC FORM : ',F15.10/) WRITE(6,720) SIGE 720 FORMAT('- ESTIMATED VARIANCE FACTOR : ',F15.10) WRITE(6,944) IDFE 944 FORMAT('- DEGREES OF FREEDOM : ',I4) DF1=FLOAT(IDFE) DF2=FLOAT(IDF12) SCRV95 = FIN(0.95,DF1,DF2) SCRV99 = FIN(0.99,DF1,DF2) WRITE(6,945) SCRV95 WRITE(6,946) SCRV99 945 FORMAT('-',T6,'F CRITICAL VALUE AT 0.05 : ',F10.4) 946 FORMAT('-',T6,'F CRITICAL VALUE AT 0.01 : ',F10.4) 9000 CALL GDATE (DATE,TIME,CODE) WRITE(6,9005) (TIME(I),I=1,6) 9005 FORMAT(' EXECUTION FINISHED AT ',2A1,':',2A1,':',2A1) STOP END C SUBROUTINE BMAT (E,N,IZC,IDV,INS,IZPC,IPC,B,IDRB,IDCB,IWCB,CHPARA) C C C INPUT : C E VECTOR OF X OR EASTING COORDINATES C N VECTOR OF Y OR NORTHING COORDINATES C IZC VECTOR OF ZONE CODES WITH SAME SUBSCRIPTS AS E AND N C IDV DECLARED DIMENSION OF E, N, IZC C INS TOTAL NUMBER OF STATIONS, WORKING DIMENSION OF E, N, IZC C LAST STATION IN LIST FIXED IN ADJUSTMENT FOR E,N OR DX C IZPC MATRIX OF ZONE CODES FOR PARAMETERS (COLUMNS OF B MATRIX) C IPC VECTOR OF PARAMETER CODES : C '1' RIGID BODY DISPLACEMENT PARALLEL TO E OR X " A O " C '2' RIGID BODY DISPLACEMENT PARALLEL TO N OR Y " B O " C '3' STRAIN IN E OR X DIRECTION " E X " C '4' STRAIN IN N OR Y DIRECTION " E Y " C '5' SHEARING STRAIN " E XY " C '6' DIFFERENTIAL ROTATION " OMEGA " C '7' SCALE " K " C DX = A1*X + A2*Y + A3*X*Y + A4*X*X + A5*Y*Y + A6*X*X*X + A7*X*X*Y C + A8*X**4 + A9*X**5 C C DY = B1*X + B2*Y + B3*X*Y + B4*X*X + B5*Y*Y + B6*X*X*X + B7*X*X*Y C C '11' "A1" C '12' "A2" C '13' "A3" C '14' "A4" C '15' "A5" C '16' "A6" C '17' "A7" C '18' "A8" C '19' "A9" C '21' "B1" C '22' "B2" C '23' "B3" C '24' "B4" C '25' "B5" C '26' "B6" C '27' "B7" C IZPC, IPC HAVE SAME SUBSCRIPTS ,THE COLUMNS OF THE B MATRIX C IWCB WORKING NUMBER OF COLUMNS (TOTAL NUMBER OF PARAMETERS) C C OUTPUT : C B POPULATED DESIGN MATRIX DECLARED AS IDRB BY IDCB, C WORKING AS 2*(INS-1) BY IWCB C CHPARA VECTOR OF PARAMETER CHARACTERS (A8) WITH SAME SUBSCRIPTS C AS IZZC, IPC, IE. COLUMNS OF B MATRIX C C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) REAL*8 AO/' A O '/,BO/' B O '/,EX/' E X '/,EY/' E Y '/, @EXY/' E XY '/,OM/' OMEGA '/,K/' K '/,A1/' A1 '/,A3/' A3 '/, @A4/' A4 '/,A5/' A5 '/,B2/' B2 '/,B3/' B3 '/,B4/' B4 '/, @B5/' B5 '/,A2/' A2 '/,B1/' B1 '/,A6/' A6 '/,A7/' A7 '/, @B6/' B6 '/,B7/' B7 '/,A8/' A8 '/,A9/' A9 '/ DIMENSION E(IDV),N(IDV),IZC(IDV),IZPC(30,IDCB),IPC(IDCB),B(IDRB, @IDCB),CHPARA(IDCB) DO 10 I=1,IDRB DO 20 J=1,IDCB B(I,J)=0.0D0 20 CONTINUE 10 CONTINUE INS1=INS-1 DO 30 I=1,INS1 JJ=2*I J=JJ-1 X=E(I)-E(INS) Y=N(I)-N(INS) DO 40 II=1,IWCB IF (IZPC(IZC(I),II).EQ.1) THEN IF (IPC(II).EQ.1) THEN B(J,II)=1.0D0 CHPARA(II)=AO ENDIF IF (IPC(II).EQ.2) THEN B(JJ,II)=1.0D0 CHPARA(II)=BO ENDIF IF (IPC(II).EQ.3) THEN B(J,II)=X CHPARA(II)=EX ENDIF IF (IPC(II).EQ.4) THEN B(JJ,II)=Y CHPARA(II)=EY ENDIF IF (IPC(II).EQ.5) THEN B(J,II)=Y B(JJ,II)=X CHPARA(II)=EXY ENDIF IF (IPC(II).EQ.6) THEN B(J,II)=-Y B(JJ,II)=X CHPARA(II)=OM ENDIF IF (IPC(II).EQ.7) THEN B(J,II)=X B(JJ,II)=Y CHPARA(II)=K ENDIF IF (IPC(II).EQ.11) THEN B(J,II)=X CHPARA(II)=A1 ENDIF IF (IPC(II).EQ.12) THEN B(J,II)=Y CHPARA(II)=A2 ENDIF IF (IPC(II).EQ.13) THEN B(J,II)=X*Y CHPARA(II)=A3 ENDIF IF (IPC(II).EQ.14) THEN B(J,II)=X**2 CHPARA(II)=A4 ENDIF IF (IPC(II).EQ.15) THEN B(J,II)=Y**2 CHPARA(II)=A5 ENDIF IF (IPC(II).EQ.16) THEN B(J,II)=X*X*X CHPARA(II)=A6 ENDIF IF (IPC(II).EQ.17) THEN B(J,II)=X*X*Y CHPARA(II)=A7 ENDIF IF (IPC(II).EQ.18) THEN B(J,II)=X*X*X*X CHPARA(II)=A8 ENDIF IF (IPC(II).EQ.19) THEN B(J,II)=X*X*X*X*X CHPARA(II)=A9 ENDIF IF (IPC(II).EQ.21) THEN B(JJ,II)=X CHPARA(II)=B1 ENDIF IF (IPC(II).EQ.22) THEN B(JJ,II)=Y CHPARA(II)=B2 ENDIF IF (IPC(II).EQ.23) THEN B(JJ,II)=X*Y CHPARA(II)=B3 ENDIF IF (IPC(II).EQ.24) THEN B(JJ,II)=X**2 CHPARA(II)=B4 ENDIF IF (IPC(II).EQ.25) THEN B(JJ,II)=Y**2 CHPARA(II)=B5 ENDIF IF (IPC(II).EQ.26) THEN B(JJ,II)=X*X*X CHPARA(II)=B6 ENDIF IF (IPC(II).EQ.27) THEN B(JJ,II)=X*X*Y CHPARA(II)=B7 ENDIF ENDIF 40 CONTINUE 30 CONTINUE RETURN END C SUBROUTINE DEGDMS (DEG,ID,IM,IS) C C DEG ANGLE IN DECIMAL SEXAGESIMAL DEGREES C C ID,IM,IS RETURNED AS DEGREES, MINUTES, SECONDS C C IMPLICIT REAL*8 (A-H,O-Z) ID=DABS(DEG) F=(DEG-ID)*60.0D0 IM=F SEC=(F-IM)*60.0D0 IF(DABS(SEC-60.0D0).GT.1.0D-6) GOTO 10 SEC=0.0D0 IM=IM+1 IF(IM.NE.60) GOTO 10 IM=0 ID=ID+1 10 CONTINUE ID=ID*DSIGN(1.0D0,DEG) IF(ID.EQ.0) IM=IM*DSIGN(1.0D0,DEG) IF(IM.EQ.0.AND.ID.EQ.0) SEC=SEC*DSIGN(1.0D0,DEG) IS=SEC RETURN END C SUBROUTINE PARADJ (IPREA,A,IDRA,IDCA,IRA,ICA,L,P,X,CX,V,CV,SIG, @ ANAT,AN,AT,AX,ATP,N,U,VT,VTP) C C LINEAR PARAMETRIC ADJUSTMENT OF L+V=A*X C C INPUT : C IPREA 0 ADJUSTMENT C 1 PREANALYSIS C A DESIGN MATRIX DECLARED AS IDRA BY IDCA C WORKING AS IRA BY ICA C L VECTOR OF OBSERVATIONS C P WEIGHT MATRIX FOR L C C OUTPUT : C X VECTOR OF ESTIMATED PARAMETERS C CX VARIANCE-COVARIANCE MATRIX FOR X C V VECTOR OF RESIDUALS C CV VARIANCE-COVARIANCE MATRIX FOR V C SIG ESTIMATED VARIANCE FACTOR C C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) DIMENSION A(IDRA,IDCA),L(IDRA),P(IDRA,IDRA),X(IDCA),CX(IDCA,IDCA) @,V(IDRA),AT(IDCA,IDRA),AX(IDRA),ATP(IDCA,IDRA),N(IDCA,IDCA), @U(IDCA),VT(1,IDRA),VTP(1,IDRA),AN(IDRA,IDCA),ANAT(IDRA,IDRA), @CV(IDRA,IDRA) CALL TRNSD (AT,IDCA,A,IDRA,IRA,ICA) CALL MMULD (ATP,IDCA,AT,IDCA,P,IDRA,ICA,IRA,IRA) CALL MMULD (N,IDCA,ATP,IDCA,A,IDRA,ICA,IRA,ICA) C WRITE(6,9906) C9906 FORMAT('1','N = BT * P * B') C CALL MOUTD (N,IDCA,ICA,ICA) CALL SPIN (N,ICA,IDCA,DET,IDEXP) C WRITE(6,9908) DET,IDEXP C9908 FORMAT(' ','INVERSE OF N HAS DETERMINANT ',F7.4,'D',I4) CALL MMULD (AT,IDCA,N,IDCA,ATP,IDCA,ICA,ICA,IRA) C WRITE(6,9907) C9907 FORMAT('1','CE * BTP') C CALL MOUTD (AT,IDCA,ICA,IRA) CALL TRNSD (AT,IDCA,A,IDRA,IRA,ICA) CALL COPYD (CX,IDCA,N,IDCA,ICA,ICA) IF(IPREA.EQ.1) GOTO 100 CALL MMULD (U,IDCA,ATP,IDCA,L,IDRA,ICA,IRA,1) CALL MMULD (X,IDCA,N,IDCA,U,IDCA,ICA,ICA,1) CALL MMULD (AX,IDRA,A,IDRA,X,IDCA,IRA,ICA,1) CALL MSUBD (V,IDRA,AX,IDRA,L,IDRA,IRA,1) CALL TRNSD (VT,1,V,IDRA,IRA,1) CALL MMULD (VTP,1,VT,1,P,IDRA,1,IRA,IRA) CALL MMULD (SIG,1,VTP,1,V,IDRA,1,IRA,1) SIG=SIG/(IRA-ICA) CALL MMULD (AN,IDRA,A,IDRA,N,IDCA,IRA,ICA,ICA) CALL MMULD (ANAT,IDRA,AN,IDRA,AT,IDCA,IRA,ICA,IRA) CALL SPIN (P,IRA,IDRA,DET,IDEXP) CALL MADDD (CV,IDRA,P,IDRA,ANAT,IDRA,IRA,IRA) RETURN 100 SIG=1.0D0 DO 110 I=1,ICA X(I)=0.0D0 110 CONTINUE RETURN END C C SUBROUTINE STRPR (EX,EY,EXY,CEX,CEY,CEXY,CEXEY,CEXEXY,CEYEXY, @ EMAX,EMIN,PHIE,C2,C1,D,DT,CIDT) C C C INPUT : C EX,EY,EXY ELEMENTS OF STRAIN TENSOR C CEX,CEXEY,CEXEXY C CEY ,CEYEXY C SYM CEXY ELEMENTS OF ASSOCIATED VAR-COV MATRIX C C C OUTPUT : C EMAX,EMIN MAXIMUM AND MINIMUM SHEAR C PHIE CCW ANGLE (RADIANS) FROM X AXIS TO DIRECTION OF EMAX C C2 ASSOCIATED VAR-COV MATRIX C C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) DIMENSION C1(3,3),D(3,3),DT(3,3),CIDT(3,3),C2(3,3) P1=(EX+EY)/2.0D0 P2=(EX-EY) P3=P2**2/4.0D0+EXY**2 P4=DSQRT(P3) EMAX=P1+P4 EMIN=P1-P4 PHIE=DATAN2(2.0D0*EXY,P2)/2.0D0 C1(1,1)=CEX C1(2,2)=CEY C1(3,3)=CEXY C1(1,2)=CEXEY C1(1,3)=CEXEXY C1(2,1)=CEXEY C1(2,3)=CEYEXY C1(3,1)=CEXEXY C1(3,2)=CEYEXY D(1,1)=0.5D0+P2/(4.0D0*P4) D(1,2)=0.5D0-P2/(4.0D0*P4) D(1,3)=EXY/P4 D(2,1)=D(1,2) D(2,2)=D(1,1) D(2,3)=-D(1,3) D(3,1)=-EXY/(4.0D0*P3) D(3,2)=-D(3,1) D(3,3)=P2/(4.0D0*P3) CALL TRNSD (DT,3,D,3,3,3) CALL MMULD (CIDT,3,C1,3,DT,3,3,3,3) CALL MMULD (C2,3,D,3,CIDT,3,3,3,3) RETURN END C C SUBROUTINE SEGREG (STN,X,Y,IZC,CX,IDRC,IDRC2,INS,IGC,STNG,IZG,XYG @,CXG,INPG,XY,S,ST,SCX) IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) DIMENSION STN(IDRC2),X(IDRC2),Y(IDRC2),CX(IDRC,IDRC) @,STNG(IDRC2),XY(IDRC),CXG(IDRC,IDRC),XYG(IDRC),S(IDRC,IDRC) @,ST(IDRC,IDRC),SCX(IDRC,IDRC),IZG(IDRC2),IZC(IDRC2) DO 10 I=1,IDRC DO 10 J=1,IDRC S(I,J)=0.0D0 ST(J,I)=0.0D0 10 CONTINUE DO 20 I=1,INS J=2*I-1 XY(J)=X(I) XY(J+1)=Y(I) 20 CONTINUE INPG=0 J=1 DO 30 I=1,INS IF(IZC(I).EQ.IGC) GO TO 30 INPG=INPG+1 STNG(INPG)=STN(I) IZG(INPG)=IZC(I) J=J+2 J2=2*I J1=J2-1 I2=2*INPG I1=I2-1 S(I1,J1)=1.0D0 S(I2,J2)=1.0D0 30 CONTINUE IWR=2*INPG IWC=2*INS CALL TRNSD (ST,IDRC,S,IDRC,IWR,IWC) CALL MMULD (XYG,IDRC,S,IDRC,XY,IDRC,IWR,IWC,1) CALL MMULD (SCX,IDRC,S,IDRC,CX,IDRC,IWR,IWC,IWC) CALL MMULD (CXG,IDRC,SCX,IDRC,ST,IDRC,IWR,IWC,IWR) RETURN END C C C C SUBROUTINE IMAXX (VEC,IDIM,IRV,ISUB,MAX) C C FINDS THE ELEMENT WITH MAXIMUM VALUE IN GIVEN VECTOR AND C RECOGNIZES ITS POSITION IN THE VECTOR C INPUT : C VEC VECTOR OF ELEMENTS TO BE EXAMINED (INTEGER) C IDIM DECLARED DIMENSION OF VEC C IRV WORKING DIMENSION OF VEC C OUTPUT : C ISUB SUBSCRIPT OF THE ELEMENT WITH MAXIMUM VALUE C MAX THE VALUE OF THAT ELEMENT C IMPLICIT INTEGER*4 (A-Z) DIMENSION VEC(IDIM) MAX=VEC(1) ISUB=1 IF(IRV.EQ.1) GO TO 20 DO 10 I=2,IRV IF(VEC(I).LE.MAX) GO TO 10 MAX=VEC(I) ISUB=I 10 CONTINUE 20 RETURN END C C SUBROUTINE ISORT (VEC,ISUBV,IDIM,IRV) C C SORTS VECTOR ELEMENTS AND RETURNS VECTOR WITH ELEMENTS IN C ASCENDING ORDER C INPUT : C VEC VECTOR OF ELEMENTS TO BE SORTED (INTEGER) C ISUBV VECTOR OF ORIGINAL SUBSCRIPTS TO ELEMENTS IN VEC C IDIM DECLARED DIMENSION OF VEC, ISUBV C IRV WORKING DIMENSION OF THE VECTORS C OUTPUT : C VEC, ISUBV RETURNED WITH SORTED ELEMENTS C ISUBV IS THEN KEY TO ORIGINAL PLACEMENT OF THE C SORTED ELEMENTS C USES SUBROUTINE IMAXX C IMPLICIT INTEGER*4(A-Z) DIMENSION VEC(IDIM),ISUBV(IDIM) DO 10 I=1,IRV CALL IMAXX (VEC,IDIM,IRV-I+1,ISUB,MAX) VEC(ISUB)=VEC(IRV-I+1) II=ISUBV(ISUB) ISUBV(ISUB)=ISUBV(IRV-I+1) ISUBV(IRV-I+1)=II VEC(IRV-I+1)=MAX 10 CONTINUE 20 RETURN END C C C SUBROUTINE DMAXX (VEC,IDIM,IRV,ISUB,MAX) C C FINDS THE ELEMENT WITH MAXIMUM VALUE IN GIVEN VECTOR AND C RECOGNIZES ITS POSITION IN THE VECTOR C INPUT : C VEC VECTOR OF ELEMENTS TO BE EXAMINED C IDIM DECLARED DIMENSION OF VEC C IRV WORKING DIMENSION OF VEC C OUTPUT : C ISUB SUBSCRIPT OF THE ELEMENT WITH MAXIMUM VALUE C MAX THE VALUE OF THAT ELEMENT C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) DIMENSION VEC(IDIM) MAX=VEC(1) ISUB=1 IF(IRV.EQ.1) GO TO 20 DO 10 I=2,IRV IF(VEC(I).LE.MAX) GO TO 10 MAX=VEC(I) ISUB=I 10 CONTINUE 20 RETURN END C C SUBROUTINE DSORT (VEC,ISUBV,IDIM,IRV) C C SORTS VECTOR ELEMENTS AND RETURNS VECTOR WITH ELEMENTS IN C ASCENDING ORDER C INPUT : C VEC VECTOR OF ELEMENTS TO BE SORTED C ISUBV VECTOR OF ORIGINAL SUBSCRIPTS TO ELEMENTS IN VEC C IDIM DECLARED DIMENSION OF VEC, ISUBV C IRV WORKING DIMENSION OF THE VECTORS C OUTPUT : C VEC, ISUBV RETURNED WITH SORTED ELEMENTS C ISUBV IS THEN KEY TO ORIGINAL PLACEMENT OF THE C SORTED ELEMENTS C USES SUBROUTINE DMAXX C IMPLICIT REAL*8(A-H,K-Z), INTEGER*4(I,J) DIMENSION VEC(IDIM),ISUBV(IDIM) DO 10 I=1,IRV CALL DMAXX (VEC,IDIM,IRV-I+1,ISUB,MAX) VEC(ISUB)=VEC(IRV-I+1) II=ISUBV(ISUB) ISUBV(ISUB)=ISUBV(IRV-I+1) ISUBV(IRV-I+1)=II VEC(IRV-I+1)=MAX 10 CONTINUE 20 RETURN END C C C C SUBROUTINE DAZ1 (DX,IDR1,IR1,DISP,IDR2,IR2,AZ) C C INPUT: C DX VECTOR OF DISPLACEMENT COMPONENTS BY DIRECT C SOLUTION FROM OBSERVATION DIFFERENCES C IDR1 DECLARED DIMENSION OF DX C IR1 WORKING DIMENSION OF DX C IDR2 DECLARED DIMENSION OF DISP, AZ C IR2 WORKING DIMENSION OF DISP, AZ C OUTPUT: C DISP VECTOR OF DISPLACEMENTS FROM FIRST TO SECOND EPOCH C IN THE SENSE OF DX C AZ VECTOR OF AZIMUTHS (DEGREES) OF THE DISPLACEMENTS C C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) DIMENSION DX(IDR1),DISP(IDR2),AZ(IDR2) PI=4.0D0*DATAN(1.0D0) IRF=2*(IR1/2) IF(IRF.NE.IR1) IR1=IR1-1 DO 100 I=1,IR1,2 J=(I+1)/2 DISP(J)=DSQRT(DX(I)**2+DX(I+1)**2) IF(DISP(J).LE.0.0D0) THEN AZ(J)=0.0D0 GO TO 100 ENDIF AZ(J)=DATAN2(DX(I),DX(I+1)) IF(AZ(J).LT.0.0D0) AZ(J)=AZ(J)+2.0D0*PI AZ(J)=AZ(J)*180.0D0/PI 100 CONTINUE RETURN END C C C C SUBROUTINE ELLIP (QXX,QXY,QYY,A,B,PHI,PHII) C C C INPUT : C QXX,QXY,QYY ELEMENTS OF VARIANCE-COVARIANCE 2 BY 2 MATRIX C C OUTPUT : C A,B LENGTHS OF MAJOR AND MINOR SEMIAXES C PHI AZIMIUTH (DEGREES) OF MAJOR SEMIAXIS C PHII CCW ANGLE (RADIANS) FROM X-AXIS TO MAJOR SEMIAXIS C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) PI=DARCOS(-1.0D0) P1=(QXX+QYY)/2.0D0 P2=DSQRT((QXX-QYY)**2/4.0D0+QXY**2) A=DSQRT(P1+P2) B=DSQRT(P1-P2) IF(QXX.LT.1.0D-20.AND.QYY.LT.1.0D-20) THEN PHI=0.0D0 PHII=PI/2.0D0 GO TO 100 ENDIF PHI=-DATAN2(-2.0D0*QXY,(QYY-QXX))/2.0D0 IF(PHI.LT.0.0D0)PHI=PHI+2.0D0*PI PHI=180.0D0*PHI/PI PHII=DATAN2(2.0D0*QXY,(QXX-QYY))/2.0D0 100 RETURN END C C C C SUBROUTINE ELLAN (X,IDX,IWX,CX,IECX,SIGX,IDF,ALPHA,IALPH,ISCALE, @ IDY,IWY,AZ,AE,BE,AZE,PHI,DE,DP) C C INPUT : C X VECTOR OF 'PARAMETERS' DECLARED AS IDX, WORKING AS IWX C CX VARIANCE-COVARIANCE MATRIX FOR X C IECX '1' CX WAS ESTIMATED ( MULTIPLIED BY SIGX ) C '0' CX WAS NOT ESTIMATED ( NOT MULTIPLIED BY SIGX ) C SIGX ESTIMATED VARIANCE FACTOR C IDF DEGREES OF FREEDOM IN ESTIMATION OF SIGX C ALPHA SIGNIFICANCE LEVEL, E.G. 0.05, FOR TESTING C IALPHA '1' BRING TO LEVEL OF ALPHA C ISCALE '1' SCALE BY SIGX C IDY,IWY DECLARED AND WORKING DIMENSIONS OF AZ,AE,BE,AZE,DE,DP C C USES SUBROUTINE ELLIP C C IMPLICIT REAL*8 (A-H,K-Z), INTEGER*4 (I,J) REAL*4 ALPHA,SCRV,ALPH,XALPH DIMENSION X(IDX),CX(IDX,IDX),AE(IDY),BE(IDY), @ AZE(IDY),DE(IDY),DP(IDY),AZ(IDY),PHI(IDY) PI=DARCOS(-1.0D0) IF(IECX.EQ.1.AND.ISCALE.NE.1) THEN DO 10 I=1,IWX DO 20 J=1,IWX CX(I,J)=CX(I,J)/SIGX 20 CONTINUE 10 CONTINUE ENDIF IF(IECX.NE.1.AND.ISCALE.EQ.1) THEN DO 30 I=1,IWX DO 40 J=1,IWX CX(I,J)=CX(I,J)*SIGX 40 CONTINUE 30 CONTINUE ENDIF FAK=1.0D0 IF(IALPH.EQ.1) SCRV = CHIIN(ALPHA,2.0) IF(IALPH.EQ.1) FAK=DSQRT(DBLE(SCRV)) DO 80 I=1,IWX,2 J=(I+1)/2 IF(CX(I,I).LT.1.0D-12) GO TO 85 IF(CX(I+1,I+1).LT.1.0D-12) GO TO 85 CALL ELLIP (CX(I,I),CX(I,I+1),CX(I+1,I+1),AE(J),BE(J),AZE(J), @ PHI(J)) DAZE=(AZ(J)-AZE(J))*PI/180.0D0 DE(J)=FAK/DSQRT((DCOS(DAZE)/AE(J))**2+(DSIN(DAZE)/BE(J))**2) DP(J)=FAK*DSQRT((AE(J)*DCOS(DAZE))**2+(BE(J)*DSIN(DAZE))**2) GO TO 80 85 AE(J)=0.0D0 BE(J)=0.0D0 AZE(J)=0.0D0 PHI(J)=0.0D0 DE(J)=0.0D0 DP(J)=0.0D0 80 CONTINUE RETURN END C C SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) C C SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC C POSITIVE-DEFINITE MATRICES. THE MATRIX INVERTED IS THE UPPER C N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM C IN THE CALLING ROUTINE. C C INPUT: C Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE C MATRIX TO BE INVERTED. C C N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER) C OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST C NOT BE LARGER THAN MM) . C MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. C C C OUTPUT: C C Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF C THE INPUT UPPER LEFT N BY N PORTION. C C DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE C INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE C IDEXP BELOW. C C IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED C UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN C TWO PARTS CORRESPONDING TO C DETERMINANT = DET * 10 ** IDEXP . C THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE C COMPUTATION OF THE DETERMINANT. TO PRINT THE DETERM- C INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS; C (FOR EXAMPLE) C PRINT 10,DET,IDEXP C 10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4) C C R.R. STEEVES C SEPT., 1979 C C IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N) DIMENSION Q(MM,MM) DET=0.D0 DO 4 J=1,N DO 4 I=1,J IF(I.EQ.1) GO TO 2 20 M=I-1 SUM=0.0D0 DO 1 K=1,M 1 SUM=SUM+Q(K,I)*Q(K,J) Q(I,J)=Q(I,J)-SUM 2 IF(I.EQ.J) GO TO 3 Q(I,J)=Q(I,J)/Q(I,I) GO TO 4 3 CONTINUE DET=DET+DLOG10(Q(I,I)) Q(I,I)=DSQRT(Q(I,I)) 4 CONTINUE IDEXP=DET RPART=DET-IDEXP APART=DABS(RPART) IF(APART.LT.1.D-20) DET=1.D0 IF(APART.LT.1.D-20) GO TO 10 DET=10**RPART 10 CONTINUE DO 7 J=1,N DO 7 I=1,J IF(I.LT.J) GO TO 5 Q(J,J)=1.0D0/Q(J,J) GO TO 7 5 SUM=0.0D0 M=J-1 DO 6 K=I,M 6 SUM=SUM-Q(I,K)*Q(K,J) Q(I,J)=SUM/Q(J,J) 7 CONTINUE DO 9 J=1,N DO 9 I=1,J SUM=0.0D0 DO 8 K=J,N 8 SUM=SUM+Q(I,K)*Q(J,K) Q(I,J)=SUM IF(I.EQ.J) GO TO 9 Q(J,I)=SUM 9 CONTINUE RETURN END C C SUBROUTINE TAURE( NT,NU,ALPH,CRTAU ) C COMPUTES THE REJECTION LEVEL FOR NORMALISED RESIDUALS FOR A GIVEN NU C OBSERVATIONS , DEGREES OF FREEDOM AND DESIRED LEVEL OF TYPE I ERROR C PARAMETERS C NT - NUMBER OF OBSERVATIONS C NU - DEGREES OF FREEDOM C ALPH - DESIRED PROBABILITY OF TYPE I ERROR C CRTAU - CRITICAL VALUE ( MAX-TAU ) PRODUCED BY THE SUBROUTINE C REFERENCE : A.J. POPE (1976) - "THE STATISTICS OF RESIDUALS AND THE C DETECTION OF OUTLIERS" , U.S. DEPT. OF COMMERCE , NOAA T C REPORT NOS 65 NGS 1. IMPLICIT REAL*8(A-H,O-Z) DATA PI/ 3.1415926535898 / PD = 2. /PI S = 1. WNU = NU U = WNU -1. IF( U.EQ.0. ) GO TO 72 IF ( ALPH.EQ.0. ) GO TO 72 IF ( ALPH.LT.1. ) GO TO 10 CRTAU = 0. 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