C ******************************************************************* C * * C * * C * * C * * C * PROGRAM STRAIN * C * -------------- * C * * C * * C * * C * THIS PROGRAM COMPUTES THE ELEMENTS OF SYMMETRIC STRAIN TENSOR * C * VALUES OF AVERAGE DIFFERENTIAL ROTATION AND OTHER COMPONENTS OF* C * STRAIN AS WELL AS DISPLACEMENTS OF THE NETWORK POINTS * C * * C * * C * NOTE THAT THE PROGRAM 'STRAIN' IS WRITTEN TO HANDLE UPTO * C * 800 STATIONS BUT IT CAN BE MADE TO HANDLE MORE THAN 800 * C * STATIONS BY SIMPLY CHANGING THE DIMENSIONS OF MATRICES WITH * C * 800 ARRAYS. IN ADDITION,IT IS ASSUMED THAT A POINT AT WHICH * C * STRAIN IS BEING COMPUTED IS NOT DIRECTLY CONNECTED BY * C * OBSERVATIONS TO MORE THAN 20 STATIONS. * C * * C * * C * JOB CONTROL LANGUAGE * C * -------------------- * C * * C * * C * 1.0 INPUT DATA SETS * C * * C * STATION NAMES ARRANGED IN TWO COLUMNS AS WILL BE EXPLAINED * C * LATER * C * * C * //GO.FT03F001 DD DSN=A.M6068.TDIRU1.DATA, * C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR * C * * C * STATION NAMES AND CO-ORDINATES OF FIRST ADJUSTMENT * C * * C * //GO.FT01F001 DD DSN=A.M6068.ADJCRD31.DATA, * C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR * C * * C * STATION NAMES AND CO-ORDINATES OF SECOND ADJUSTMENT * C * * C * //GO.FT15F001 DD DSN=A.M6068.TCRAA1D.DATA, * C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR * C * * C * * C * 2.0 OUTPUT DATA SETS * C * * C * ELEMENTS OF THE SYMMETRIC STRAIN TENSOR (AS INPUT FOR * C * PROGRAMS 'EVALUE AND 'NETPLOT') * C * * C * //GO.FT17F001 DD DSN=A.M6068.TMATRIX2.DATA, * C * // UNIT=P3350,VOL=SER=USER11,DISP=(NEW,CATLG,DELETE), * C * // DCB=(RECFM=FB,LRECL=60,BLKSIZE=6000),SPACE(TRK,(10,15),RLSE)* C * * C * VALUES OF AVERAGE DIFFERENTIAL ROTATION AND DISPLACEMENTS * C * (AS INPUT FOR PROGRAM 'NETPLOT') * C * //GO.FT09F001 DD DSN=A.M6068.STNDXDY2.OMEGA, * C * // UNIT=DASD,VOL=SER=USER11,DISP=(NEW,CATLG,DELETE), * C * // DCB=(RECFM=FB,LRECL=80,BLKSIZE=8000),SPACE=(TRK,(10,15),RLSE)* C * //GO.SYSIN DD * * C * // * C * * C * * C * INPUT: * C * * C * ALL INPUT FOR THIS PROGRAM IS FROM FILES AND NO DATA CARDS ARE * C * REQUIRED. THE CONTROL LANGUAGE REQUIREMENTS ARE GIVEN * C * IN THE PROGRAM LISTING * C * * C * N,X,Y STATION NUMBER AND CO-ORDINATES OF FIRST ADJUSTMENT * C * N3,X2,Y2 STATION NUMBER AND CO-ORDINATES OF SECOND ADJUSTMENT * C * FORMAT FOR ABOVE INPUT CO-ORDINATES IS:FORMAT(I8,2F15.4) * C * * C * N1,N2 STATION NUMBERS ARRANGED IN TWO COLUMNS SUCH THAT THE * C * FIRST COLUMN CONTAINS STATION NUMBERS TO WHICH ALL THE POINTS * C * IN THE SECOND COLUMN ARE CONNECTED BY OBSERVATIONS * C * FORMAT FOR N1 AND N2 IS FORMAT(2X,I8,2X,I8) * C * * C * NOTE THAT THE ABOVE CO-ORDINATES CAN BE THE DIRECT OUTPUT * C * FROM ADJUSTMENT PROGRAM OF THE NETWORK * C * * C * * C * OUTPUT: * C * * C * E(1) AND F(2) ARE THE DIAGONAL ELEMENTS OF THE SYMMETRIC * C * 2BY2 STRAIN TENSOR MATRIX AND EP12 IS ITS OFF-DIAGONAL ELEMENT * C * THESE THREE ELEMENTS ARE OUTPUT IN THE FORMAT FORMAT(1X,3F18.10* C * ) AND THESE ARE USED BY BOTH PROGRAMS 'EVALUE'AND 'NETPLOT' * C * * C * * C * OMEGA IS THE AVERAGE DIFFERENTIAL ROTATION * C * DX1 AND DY1 ARE THE DISPLACEMENTS IN X AND Y-DIRECTIONS * C * RESPECTIVELY * C * THESE ARE OUTPUT IN THE FORMAT FORMAT(I8,2F10.4,F16.10) * C * AND ARE USED BY THE PROGRAM 'NETPLOT'. * C * * C * * C * GM1=GAMMA ONE * C * GM2=GAMMA TWO * C * GM=GAMMA * C * DILAT=DILATATION * C * SHR=SHEAR * C * THESE ARE VARIOUS OTHER COMPONENTS OF STRAIN USEFUL TO STUDY * C * THE CRUSTAL MOVEMENT. THESE ARE HOWEVER NOT USEFUL FOR OUR * C * INVESTIGATIONS * C * * C * THE READER IS ADVISED TO REFER TO THE REFERENCES TO FIND OUT * C * AS TO WHAT THE ABOVE QUANTITIES MEAN * C * * C * REFERENCES: * C * * C * THAPA,K.(1980).STRAIN AS A DIAGNOSTIC TOOL TO SPOT INCONSIS- * C * TENT OBSERVATIONS AND CONSTRAINTS.M.SC.E. THESIS IN THE * C * DEPT.OF SURVEYING ENGINEERING UNB,FREDERICTON,N.B. CANADA * C * * C * RIKITAKE,T.(1976). EARTHQUAKE PREDICTION,ELSEVIER * C * * C * * C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979 * C * * C * * C ******************************************************************* C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(20,3),B(20),C(3,3),D(20),E(3),F(3),T(3,20),P(3,3) DIMENSION N(800),N3(800),X(800),Y(800),X2(800),Y2(800) DIMENSION DX(20),DY(20),ST(20) DIMENSION G(20),H(20) PI=DARCOS(-1.D0) N(1)=0 DO 1 I=2,800 C C READ CO-ORDINATES OF FIRST ADJUSTMENT AND STATION NUMBERS C READ(1,200,END=2)N(I),X(I),Y(I) C C READ CO-ORDINATES OF SECOND ADJUSTMENT AND STATION NUMBERS C 1 READ(15,200)N3(I),X2(I),Y2(I) 2 ICOUNT=I-1 K=1 C C READ THE STATION NUMBERS ARRANGED IN TWO COLUMNS AS C EXPLAINED EARLIER C 30 READ(3,100,END=50)N1,N2 100 FORMAT(2X,I8,2X,I8) IF(N1.EQ.N(K))GOTO 4 IF(K.EQ.1)GOTO 6 IF(N1.EQ.6910170)GOTO 30 C C BY CALLING THESE SUBROUTINES SOLUTION FOR DISPLACEMENT C GRADIENTS ARE OBTAINED C 50 CALL MATTNP(A,20,T,3,J,3) CALL MATMUL(T,3,A,20,P,3,3,J,3) CALL CHOLD(P,3,3,DETA) CALL MTVMUL(T,3,DX,20,G,20,3,J) CALL MTVMUL(T,3,DY,20,H,20,3,J) CALL MTVMUL(P,3,G,20,E,3,3,3,3) CALL MTVMUL(P,3,H,20,F,3,3,3,3) WRITE(6,128)N(K) 128 FORMAT(7X,'THE VALUES OF DISPLACEMENT GRADIENTS'/, @ ' FOR STATION NUMBER',I8/) C C COMPUTE OFF-DIAGONAL ELEMENT OF THE SYMMETRIC STRAIN TENSOR C EP12=0.5*(E(2)+F(1)) WRITE(6,451)E(1),EP12,F(2) 451 FORMAT(4X,3(3X,D14.8)/) 455 FORMAT(1X,3F18.10) C C WRITE THE ELEMENTS OF SYMMETRIC STRAIN TENSOR IN A FILE C WITH UNIT 17 C WRITE(17,455)E(1),EP12,F(2) WRITE(6,410)E(3),F(3) C C PRINT THE VALUES OF CONSTANTS C 410 FORMAT(7X,'THE VALUES OF C AND D ARE RESPECTIVELY'/,3X, @2(4X,D14.8)/) C C COMPUTE AND PRINT VARIOUS OTHER COMPONENTS OF STRAIN TENSOR C THESE ARE ONLY USEFUL FOR CRUSTAL ANALYSIS C GM1=E(1)-F(2) GM2=2*E(2) GM=DSQRT(GM1*GM1+GM2*GM2) OMEGA=0.5*(E(2)-F(1)) WRITE(6,360)GM1,GM2,GM,OMEGA 360 FORMAT(7X,'THE VALUES OF GM1,GM2,GM AND OMEGA ARE RESPECTIVELY'/ @4X,2(3X,D14.8)/,4X,2(3X,D14.8)/) DX1=X2(K)-X(K) DY1=Y2(K)-Y(K) WRITE(6,138)N(K),DX1,DY1,OMEGA C C WRITE STATION NUMBERS DISPLACEMENTS AND VALUES OF C AVERAGE DIFFERENTIAL ROTATION ON A FILE I N UNIT 9 C WRITE(9,135)N(K),DX1,DY1,OMEGA 135 FORMAT(I8,2F10.4,F16.10) 138 FORMAT(7X,'THE STATION NAME DISPLACEMENTS & OMEGA ARE'/ @,3X,I8,2F10.4,D14.8/) SHR=2*EP12 DILAT=E(1)+F(2) SIGMA=DSQRT(GM1*GM1+SHR*SHR) EPSLN=SHR+SIGMA IF(DILAT.EQ.0.0) GOTO 23 TNTH=(-GM1+SIGMA)/DILAT DIRETN=DATAN(TNTH)*180.D0/PI WRITE(6,500)SHR,DILAT,SIGMA,EPSLN,DIRETN 23 CONTINUE 500 FORMAT(7X,'THE VALUES OF SHEAR DILATION,MAX.SHEAR STRAIN'/, @7X,'AND ITS DIRECTION ARE RESPECTIVELY'/,7X,3F14.9/,7X,2F14.9/) IF(K.EQ.ICOUNT)GO TO 20 6 CONTINUE K=K+1 JJ=1 ST(1)=N1 SX=X(K) SY=Y(K) C C COMPUTE THE DIS PLACEMENTS AND ELEMENTS OF 'A' MATRIX C FOR THE STATION AT WHICH THE STRAIN ELLIPSE IS BEING COMPUTED C DX(1)=X2(K)-X(K) DY(1)=Y2(K)-Y(K) A(1,1)=0.0 A(1,2)=0.0 A(1,3)=1.0 4 CONTINUE DO 10 I=1,JJ 10 IF(N2.EQ.ST(I))GO TO 30 JJ=I J=JJ ST(I)=N2 DO 5 I=1,ICOUNT IF(N2.NE.N3(I))GOTO 5 C C COMPUTE THE ELEMENTS OF 'A' MATRIX AND THE DISPLACEMENTS C FOR STATIONS WHICH ARE CONNECTED BY OBSERVATIONS TO THE POINT C AT WHICH VALUE OF AVE. DIFF. ROTATION AND STRAIN ELLIPSE C ARE BEING COMPUTED C DX(J)=X2(I)-X(I) DY(J)=Y2(I)-Y(I) A(J,1)=X(I)-SX A(J,2)=Y(I)-SY A(J,3)=1.0 5 CONTINUE 200 FORMAT(I8,2F15.4) GOTO 30 20 STOP END C ******************************************************************* C * * C * * C * SUBROUTINE MTVMUL * C * ----------------- * C * * C * * C * THIS SUBROUTINE MULTIPLIES A MATRIX OF DIMENSION M BY N WITH * C * A COLUMN VECTOR OF N BY 1 C * A AND B ARE THE MATRIX AND A COLUMN VECTOR TO BE MULTI- * C * PLIED RESPECTIVELY * C * C IS THE RESULTANT VECTOR. * C * IA,IB,IC ARE THE ROW DIMENSIONS OF A AND COLUMN DIMENSIONS * C * OF B AND C RESPECTIVELY * C * M BY N IS THE SIZE OF A * C * * C * * C * INPUT: * C * * C * MATRIX A AND COLUMN VECTOR B * C * * C * * C * OUTPUT: * C * * C * COLUMN VECTOR C * C * * C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. * C * * C * * C ******************************************************************* C C SUBROUTINE MTVMUL(A,IA,B,IB,C,IC,M,N) REAL*8 A(IA,N),B(IB),C(IC) DO 10 I=1,M C(I)=0 DO 10 J=1,N 10 C(I)=A(I,J)*B(J)+C(I) RETURN END C ******************************************************************* C * * C * * C * SUBROUTINE MATMUL * C * ----------------- * C * * C * * C * THIS SUBROUTINE MULTIPLIES TWO MATRICES A AND B AND C IS THE * C * RESULTING MATRIX * C * * C * IA,IB,IC ARE THE ROW DIMENSIONS OF A,B AND C RESPECTIVELY. * C * L BY M IS THE SIZE OF A * C * M BY N IS THE SIZE OF B * C * L BY N IS THE SIZE OF C * C * * C * * C * INPUT: * C * * C * MATRICES A AND B * C * * C * * C * OUTPUT: * C * * C * MATRIX C * C * * C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. * C * * C * * C ******************************************************************* C C SUBROUTINE MATMUL(A,IA,B,IB,C,IC,L,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(IA,M),B(IB,N),C(IC,N) DO 2 I=1,L DO 2 J=1,N C(I,J)=0.0D0 DO 1 K=1,M C(I,J)=C(I,J)+A(I,K)*B(K,J) 1 CONTINUE 2 CONTINUE RETURN END C ******************************************************************* C * * C * * C * * C * SUBROUTINE MATTNP * C * ----------------- * C * * C * * C * THIS SUBROUTINE TRANSPOSES MATRIX A TO MATRIX B * C * * C * IDA AND IDB ARE THE ROW DIMENSIONS OF A AND B RESPECTIVELY * C * M BY N IS THE SIZE OF A * C * N BY M IS THENSIZE OF B * C * * C * * C * INPUT: * C * * C * MATRIX A * C * * C * * C * OUTPUT: * C * * C * MATRIX B * C * * C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. * C * * C * * C ******************************************************************* C C SUBROUTINE MATTNP(A,IDA,B,IDB,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(IDA,N),B(IDB,M) DO 1 I=1,M DO 2 J=1,N B(J,I)=A(I,J) 2 CONTINUE 1 CONTINUE RETURN END C ******************************************************************* C * * C * * C * SUBROUTINE CHOLD * C * ---------------- * C * * C * * C * THIS SUBROUTINE DOES MATRIX INVERSION USING CHOLESKI * C * DECOMPOSITION * C * * C * * C * INPUT: * C * * C * A=POSITIVE DEFINITE SYMMETRIC MATRIX * C * IRDA =ROW DIMENSION OF MATRIX A * C * NA=SIZE OF MATRIX A * C * * C * * C * OUTPUT: * C * * C * A=INVERSE OF INPUT MATRIX(INPUT DESTROYED * C * DETA=DETERMINANT OF MATRIX A * C * * C * * C ******************************************************************* C SUBROUTINE CHOLD(A,IRDA,NA,DETA) C 0430 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION 0440 C 0450 C INPUT ARGUMENTS 0460 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX 0470 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX 0480 C NA = SIZE OF INPUT MATRIX 0490 C OUTPUT ARGUMENTS 0500 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) 0510 C DETA = DETERMINANT OF INPUT MATRIX 0520 C * = ERROR RETURN (TAKEN IF NA .LT. 1 OR IF DETA .LT. SING) 0530 C 0540 DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING 0550 DIMENSION A(IRDA,NA) 0560 SQRT(SUM) = DSQRT(SUM) 0570 ABS(DETA) = DABS(DETA) 0580 DATA SING/1D-10/ 0590 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 0600 IF(NA .LT. 1) GO TO 18 0610 DETA = A(1,1) 0620 A(1,1) = SQRT(A(1,1)) 0630 IF(NA .EQ. 1) GO TO 6 0640 DO 1 I = 2,NA 0650 1 A(I,1) = A(I,1) / A(1,1) 0660 DO 5 J = 2,NA 0670 SUM = 0. 0680 J1 = J - 1 0690 DO 2 K = 1,J1 0700 2 SUM = SUM + A(J,K) ** 2 0710 DETA = DETA * (A(J,J) - SUM) 0720 A(J,J) = SQRT(A(J,J) - SUM) 0730 IF(J .EQ. NA) GO TO 5 0740 J2 = J + 1 0750 DO 4 I = J2,NA 0760 SUM = 0. 0770 DO 3 K = 1,J1 0780 3 SUM = SUM + A(I,K) * A(J,K) 0790 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 0800 5 CONTINUE 0810 6 IF(ABS(DETA) .LT. SING) GO TO 16 0820 C INVERSION OF LOWER TRIANGULAR MATRIX 0830 DO 7 I = 1,NA 0840 7 A(I,I) = 1. / A(I,I) 0850 IF(NA .EQ. 1) GO TO 10 0860 N1 = NA - 1 0870 DO 9 J = 1,N1 0880 J2 = J + 1 0890 DO 9 I = J2,NA 0900 SUM = 0. 0910 I1 = I - 1 0920 DO 8 K = J,I1 0930 8 SUM = SUM + A(I,K) * A(K,J) 0940 9 A(I,J) = - A(I,I) * SUM 0950 C CONSTRUCTION OF INVERSE OF INPUT MATRIX 0960 10 DO 15 J = 1,NA 0970 IF(J .EQ. 1) GO TO 12 0980 J1 = J - 1 0990 DO 11 I = 1,J1 1000 11 A(I,J) = A(J,I) 1010 12 DO 14 I = J,NA 1020 SUM = 0. 1030 DO 13 K = I,NA 1040 13 SUM = SUM + A(K,I) * A(K,J) 1050 14 A(I,J) = SUM 1060 15 CONTINUE 1070 RETURN 1080 16 WRITE(6,17) DETA 1090 17 FORMAT(10X, 'SINGULAR MATRIX IN CHOLD. DET =',E20.5) 1100 RETURN 1 1110 18 WRITE(6,19) 1120 19 FORMAT(10X,'MATRIX OF DIMENSION ZERO IN CHOLD') 1130 RETURN 1 1140 END 1150