C * ******************************************************************10 C * *11 C * *13 C * *14 C * PROGRAM STRAIN *15 C * -------------- *16 C * *17 C * *18 C * *19 C * THIS PROGRAM COMPUTES THE ELEMENTS OF SYMMETRIC STRAIN TENSOR *20 C * VALUES OF AVERAGE DIFFERENTIAL ROTATION AND OTHER COMPONENTS OF *21 C * *23 C * *24 C * NOTE THAT THE PROGRAM 'STRAIN' IS WRITTEN TO HANDLE UPTO *25 C * 800 STATIONS BUT IT CAN BE MADE TO HANDLE MORE THAN 800 *26 C * STATIONS BY SIMPLY CHANGING THE DIMENSIONS OF MATRICES WITH *27 C * 800 ARRAYS. IN ADDITION,IT IS ASSUMED THAT A POINT AT WHICH *28 C * STRAIN IS BEING COMPUTED IS NOT DIRECTLY CONNECTED BY *29 C * OBSERVATIONS TO MORE THAN 20 STATIONS. *30 C * *31 C * JOB CONTROL LANGUAGE *33 C * -------------------- *34 C * *35 C * *36 C * 1.0 INPUT DATA SETS *37 C * *38 C * STATION NAMES ARRANGED IN TWO COLUMNS AS WILL BE EXPLAINED *39 C * LATER *40 C * *41 C * //GO.FT03F001 DD DSN=A.M6068.TDIRU1.DATA, *42 C * //UNIT=P3350,VOL=SER=USER11,DISP=SHR *43 C * *44 C * STATION NAMES AND CO-ORDINATES OF FIRST ADJUSTMENT *45 C * *46 C * //GO.FT01F001 DD DSN=A.M6068.ADJCRD31.DATA, *47 C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR *48 C * *49 C * STATION NAMES AND CO-ORDINATES OF SECOND ADJUSTMENT *50 C * *51 C * //GO.FT15F001 DD DSN=A.M6068.TCRAA1D.DATA, *52 C * // UNIT P3350,VOL=SER=USER11,DISP=SHR *53 C * *54 C * *55 C * 2.0 OUTPUT DATA SETS *56 C * *57 C * ELEMENTS OF THE SYMMETRIC STRAIN TENSOR (AS INPUT FOR *58 C * PROGRAMS 'EVALUE AND 'NETPLCT') *59 C * *60 C * //GO.FT17F001 DD DSN=A.M6068.TMATRIX2.DATA, *61 C * // UNIT=P3350,VOL=SER=USER11,DISP=(NEW,CATIG,DELETE) *62 C * // DCB=(RECFM=FB,LRECL=60,BLKSIZE=6000),SPACE(TRK,(10,15),RLSE) *63 C * *64 C * VALUES OF AVERAGE DIFFERENTIAL ROTATION AND DISPLACEMENTS *65 C * (AS INPUT FOR PROGRAM 'NETPLOT') *66 C * //GO.FT09F001 DD DSN=A.M6068.STNDXDY2.OMEGA, *67 C * // UNIT=DASD,VOL=SER=USER11,DISP=(NEW,CATLG,DELETE), *68 C * // DCB=(RECFM=FB,LRECL=80,BLKSIZE=8000),SPACE=(TRK,(10,15),RLSE) *69 C * //GO.SYSIN DD * *70 C * // *71 C * *72 C * *73 C * IMPUT= *74 C * *75 C * ALL INPUT FOR THIS PROGRAM IS FROM FILES AND NO DATA CARDS ARE *76 C * REQUIRED. THE CONTROL LANGUAGE REQUIREMENTS ARE GIVEN *77 C * IN THE PROGRAM LISTING *78 C * *79 C * 79.1 C * N,X,Y STATION NUMBER AND CO-ORDINATES OF FIRST ADJUSTMENT *80 C * N3,X2,Y2 STATION NUMBER AND CO-ORDINATES OF SECOND ADJUSTMENT *81 C * FORMAT FOR ABOVE INPUT CO-ORDINATES IS:FORMAT(I8,2F15.4) *82 C * *83 C * N1,N2 STATION NUMBERS ARRANGED IN TWO COLUMNS SUCH THAT THE *84 C * FIRST COLUMN CONTAINS STATION NUMBERS TO WHICH ALL THE POINTS *85 C * IN THE SECOND COLUMN ARE CONNECTED BY OBSERVATIONS *86 C * FORMAT FOR N1 AND N2 IS FORMAT(2X,I8,2X,I8) *87 C * *88 C * NOTE THAT THE ABOVE CO-ORDINATES CAN BE THE DIRECT OUTPUT *89 C * FROM ADJUSTMENT PROGRAM OF THE NETWORK *90 C * *91 C * *92 C * OUTPUT: *93 C * *94 C * E(1) AND F(2) ARE THE DIAGONAL ELEMENTS OF THE SYMMETRIC *95 C * 2BY2 STRAIN TENSOR MATRIX AND EP12 IS ITS OFF-DIAGONAL ELEMENT *96 C * THESE THREE ELEMENTS ARE OUTPUT IN THE FORMAT FORMAT(1X,3F18,10 *97 C * ) AND THESE ARE USED BY BOTH PROGRAMS 'EVALUE'AND 'NETPLOT' *98 C * *99 C * *100 C * OMEGA IS THE AVERAGE DIFFERENTIAL ROTATION *101 C * DX1 AND DY1 ARE THE DISPLACEMENTS IN X AND Y-DIRECTIONS *102 C * RESPECTIVELY *103 C * THESE ARE OUTPUT IN THE FORMAT FORMAT(I8,2F10.4,F16.10) *104 C * AND ARE USED BY THE PROGRAM 'NETPLOT'. *105 C * *106 C * *107 C * GM1=GAMMA ONE *108 C * GM2=GAMMA TWO *109 C * GM=GAMMA *110 C * DILAT=DILATATION *111 C * SHR=SHEAR *112 C * THESE ARE VARIOUS OTHER COMPONENTS OF STRAIN USEFUL TO STUDY *113 C * THE CRUSTAL MOVEMENT. THESE ARE HOWEVER NOT USEFUL FOR OUR *114 C * INVESTIGATIONS *115 C * *116 C * THE READER IS ADVISED TO REFER TO THE REFERENCES TO FIND OUT *117 C * AS TO WHAT THE ABOVE QUANTITIES MEAN *118 C * *119 C * REFERENCES: *120 C * *121 C * THAPA,K.(1980).STRAIN AS A DIAGNOSTIC TOOL TO SPOT INCONSIS- *122 C * TENT OBSERVATIONS AND CONSTRAINTS.M.SC.E. THESIS IN THE *123 C * DEPT. OF SURVEYING ENGINEERING UNB,FREDERICTON,N.B. CANADA *124 C * *125 C * RIKITAKE,T.(1976). EARTHQUAKE PREDICTION,ELSEVIER *126 C * *127 C * *128 C * WRITTEN BY KHAGENDRA THAPA NOVEMBER, 1979 *129 C * *130 C * *131 C * ******************************************************************132 C 133 C 134 PROGRAM STRAIN(INPUT,OUTPUT,TAPE1,TAPE15,TAPE3,TAPE2,TAPE11, 134.1 +TAPE17,TAPE9,TAPE4,TAPE5,TAPE7=INPUT,TAPE6=OUTPUT,TAPE8) 134.2 LOGICAL LSOBS,LAMB 134.3 DIMENSION A(60,3),E(3),F(3),T(3,60),P(3,3) 136 DIMENSION N(800),N3(800),X(800),Y(800),X2(800),Y2(800) 137 DIMENSION DX(60),DY(60),ST(60) 138 DIMENSION G(60),H(60),NETWORK(8),C(7) 139 DATA ICON,ORLAT,ORLON,CRLAT/60,42.,90.,80./ 139.05 C 139.1 LSOBS=.FALSE. 139.2 LAMB=.FALSE. 139.3 PI=4.*ATAN(1.) 140 C 140.001 C READ/WRITE TITEL 140.002 C 140.003 READ(7,97)NETWORK 140.01 97 FORMAT(8A10) 140.02 WRITE(6,88)NETWORK 140.03 WRITE(11,97) NETWORK 140.035 88 FORMAT("1",4X,8A10//) 140.04 C 140.05 C READ PARAMETER CARD 140.06 C 140.07 READ(7,98)CM,RCRIT,BLAT,CLAT 140.08 98 FORMAT(F10.4,F10.0,2F10.4) 140.081 IF(BLAT.NE.0)LAMB=.TRUE. 140.0815 IF(RCRIT.NE.0.)LSOBS=.TRUE. 140.082 WRITE(6,99) 140.083 99 FORMAT(1H ,3X,"OPTION SELECTED:"//) 140.084 IF(LSOBS)WRITE(6,101)RCRIT 140.085 IF(.NOT.LSOBS)WRITE(6,102) 140.086 IF(LAMB)WRITE(6,103)BLAT,CLAT 140.0865 101 FORMAT(3X,"SIMULATED CONNECTIONS , CRITICAL RADIUS ›M!=", 140.087 +F11.0/) 140.088 102 FORMAT(3X,"CONNECTIONS READ FROM OBSERVATION FILE"/) 140.089 103 FORMAT(3X,"OPTION LAMB: LAMBERT CONIC PROJECTION:"/ 140.09 +3X,"LAT1 ›DEG!=",F10.4,3X,"LAT2 ›DEG!=",F10.4/) 140.095 IF(.NOT.LAMB)WRITE(6,193)CM 140.096 193 FORMAT(1H ,3X,"CENTRAL MERIDIAN USED IN UTM-PROJECTION ›DEG!:", 140.097 +F10.4/) 140.098 WRITE(6,91) 140.1 91 FORMAT(13X,"FIRST ADJUSTMENT",35X,"SECOND ADJUSTMENT"//) 140.2 WRITE(6,92) 140.3 92 FORMAT(2X,"NUMBER",13X,"X1",18X,"Y1",14X,"NUMBER",13X,"X2",16X, 140.4 +"Y2",15X,"IZONE"/) 140.5 N(1)=0 141 C 141.001 C THESE SUBROUTINES ARE USED TO SORT THE STATION NUMBERS AND 141.01 C COORDINATE OF THE FIRST ADJUSTMENT IN A CERTAIN ORDER 141.02 C 141.03 CALL SMSORT(80) 141.1 CALL SMFILE("SORT","CODED",1,"REWIND") 141.2 CALL SMFILE("OUTPUT","CODED",2,"REWIND") 141.3 CALL SMKEY(8,1,8,0,"DISPLAY") 141.4 CALL SMEND 141.5 C 141.51 C THESE SUBROUTINES ARE CALLED TO SORT THE STATION NUMBERS AND 141.52 C COORDINATES OF THE SECOND ADJUSTMENT IN A CERTAIN ORDER 141.53 C 141.54 CALL SMSORT(80) 141.6 CALL SMFILE("SORT","CODED",15,"REWIND") 141.7 CALL SMFILE("OUTPUT","CODED",4,"REWIND") 141.8 CALL SMKEY(8,1,8,0,"DISPLAY") 141.9 CALL SMEND 141.91 DO 1 I=2,800 142 C 143 C READ CO-ORDINATES OF FIRST ADJUSTMENT AND STATION NUMBERS 144 C 145 C 146.09 READ(2,200)MC,N(I),M1,M2,SDLA,SMLA,SLA,SDLO,SMLO,SLO 146.1 C THIS STATEMENT IS INSERTATED TO INDICATE THE END OF DATA IN THE FILE 146.11 C 146.12 IF(EOF(2).NE.0) GOTO 2 146.15 C CONVERT LATITUDES AND LONGITUDES INTO DECIMAL DEGREES 146.16 C 146.17 SLAT1=(SDLA+SMLA/60.0+SLA/3600.0) 146.2 SLON1=(SDLO+SMLO/60.0+SLO/3600.0) 146.3 IF(I.GT.2.OR.CM.GT.0.) GOTO 39 146.31 CM=SLON1 146.32 39 CONTINUE 146.33 C 146.34 C OPTION UTM: COMPUTE UTM-COORDINATES 146.35 C 146.37 81 IF(LAMB)GOTO 82 146.38 CALL TAS002(SLAT1,SLON1,SNORTH1,SEAST1,IDUM,CM) 146.4 X(I)=SEAST1 146.5 Y(I)=SNORTH1 146.6 82 CONTINUE 146.7 C 147 C OPTION LAMB: COMPUTE COORDINATES OF THE LAMBERT CONIC PROJECTION 147.1 83 IF(.NOT.LAMB)GOTO 84 147.2 C DEFINE CONSTANTS 147.3 C(1)=BLAT 147.4 C(2)=CLAT 147.5 C(3)=ORLAT 147.6 C(4)=ORLON 147.7 C(5)=C(6)=C(7)=0. 147.8 CALL LAMCON(C) 147.9 CALL LAM(C,SLAT1,SLON1,X(I),Y(I),AM1) 147.91 84 CONTINUE 147.92 C 147.93 C READ CO-ORDINATES OF SECOND ADJUSTMENT AND STATION NUMBERS 148 C 149 READ(4,200)MC2,N3(I),M3,M4,SDLA1,SMLA1,SLA1,SDLO1,SMLO1,SLO1 150 200 FORMAT(1X,I2,4X,3A8,8X,2F3.0,F10.5,2F3.0,F10.5) 150.01 SLAT2=(SDLA1+SMLA1/60.0+SLA1/3600.0) 150.1 SLON2=(SDLO1+SMLO1/60.0+SLO1/3600.0) 150.2 CALL TAS002(SLAT2,SLON2,SNORTH2,SEAST2,IDUM,CM) 150.3 X2(I)=SEAST2 150.4 Y2(I)=SNORTH2 150.5 WRITE(6,66)N(I),X(I),Y(I),N3(I),X2(I),Y2(I),IDUM 150.7 66 FORMAT(3X,A8,2F19.5,10X,A8,2F19.5,5X,I4) 150.8 1 CONTINUE 150.9 2 ICOUNT=I-1 151 I2=I-2 151.1 WRITE(6,191)I2 151.5 191 FORMAT(1H0,3X,"NUMBER OF STATIONS:",I4/) 151.6 WRITE(11,192)I2 151.7 192 FORMAT(I4) 151.8 51 IF(LSOBS)GOTO 52 151.871 C 151.8710 C THESE SUBROUTINES ARE CALLED TO SORT THE OBSERVED DIRECTIONS AND 151.8710 C DISTANCES INTO BLOCKS SO THAT THEY ARE ARRANGED IN SUCH A WAY 151.8711 C THAT ALL THE STATIONS TO WHICH OBSERVATIONS ARE MADE ARE IN 151.8712 C BLOCKS 151.8712 C 151.8713 CALL SMSORT(80) 151.8714 CALL SMFILE("SORT","CODED",3,"REWIND") 151.8714 CALL SMFILE("OUTPUT","CODED",5,"REWIND") 151.8715 CALL SMKEY(8,1,8,0,"DISPLAY") 151.8716 CALL SMKEY(2,1,2,0,"DISPLAY") 151.8717 CALL SMEND 151.8717 IFILE=5 151.88 GOTO 78 151.902 C 151.903 C OPTION SOBS: SIMULATED CONNECTIONS BETWEEN NEIGHBOURING 151.904 C POINTS FOR THE STRAIN COMPUTATION 151.905 C 151.906 52 WRITE(6,195) 151.9064 195 FORMAT(1H1,3X,"CONNECTIONS BETWEEN NEIGHBOURING POINTS:"/ 151.9068 +1X,"TYP",3X,"FROM",5X,"TO"/) 151.9072 RCRIT2=RCRIT**2 151.908 M5=0 151.9085 41 DO 42 KSIM=2,ICOUNT 151.909 JCON=0 151.91 43 DO 44 JSIM=2,ICOUNT 151.911 IF(KSIM.EQ.JSIM)GOTO 44 151.9115 IF((X(JSIM)-X(KSIM))**2+(Y(JSIM)-Y(KSIM))**2.GT.RCRIT2)GOTO 44 151.912 53 IF(JCON.LT.ICON)GOTO 54 151.9121 WRITE(6,196)ICON 151.9122 196 FORMAT(3X,"NUMBER OF CONNECTIONS EXCEEDS:",I3,3X,"EXCEEDING CON- 151.9123 +NECTIONS NEGLECTED"/) 151.9124 GOTO 42 151.9125 54 JCON=JCON+1 151.9126 N1=N(KSIM) 151.913 N2=N3(JSIM) 151.914 WRITE(6,194)M5,N1,N2 151.915 WRITE(8,100)M5,N1,N2 151.916 194 FORMAT(1H ,3X,I2,2X,A8,2X,A8) 151.917 44 CONTINUE 151.918 42 CONTINUE 151.919 C 153 IFILE=8 154 78 REWIND 8 154.5 C 155 WRITE(6,86) 155.1 86 FORMAT("1",2X,"STATION NAME",2X,"1ST DIAGONAL ELEMENT",2X, 155.2 +"OFF DIAGONAL ELEMENT",2X,"2ND DIAGONAL ELEMENT",5X,"OMEGA", 155.3 +10X,"CONSTANTS",11X,"DISPLACEMENTS"//) 155.4 C 156 K=1 156.1 C 156.2 C READ NAMES OF CONNECTED STATIONS FROM FILE 5 OR 8 156.3 C 156.4 30 READ(IFILE,100)M5,N1,N2 157 100 FORMAT(1X,I2,4X,A8,16X,A8) 157.0001 IF(EOF(IFILE).NE.0)GO TO 50 157.001 IF(M5.NE.0.AND.M5.NE.1.AND.M5.NE.2) GOTO 30 157.01 IF(N1.EQ.N(K))GOTO 4 159 IF(K.EQ.1)GOTO 6 160 C 162 C BY CALLING THESE SUBROUTINES SOLUTION FOR DISPLACEMENT 163 C GRADIENTS ISE OBTAINED 164 C 165 IF(J.LE.2) GOTO 61 165.1 50 CALL MATTNP(A,ICON,T,3,J,3) 166 CALL MATMUL(T,3,A,ICON,P,3,3,J,3) 167 CALL CHOLD(P,3,3,DETA) 168 CALL MTVMUL(T,3,DX,ICON,G,ICON,3,J) 169 CALL MTVMUL(T,3,DY,ICON,H,ICON,3,J) 170 CALL MTVMUL(P,3,G,ICON,E,3,3,3) 171 CALL MTVMUL(P,3,H,ICON,F,3,3,3) 172 GOTO 62 172.001 61 CONTINUE 172.01 E(1)=0.0 172.02 E(2)=0.0 172.03 E(3)=0.0 172.04 F(1)=0.0 172.05 F(2)=0.0 172.06 F(3)=0.0 172.07 62 CONTINUE 172.08 DX1=X2(K)-X(K) 172.1 DY1=Y2(K)-Y(K) 172.2 EP12=0.5*(E(2)+F(1)) 173 OMEGA=0.5*(E(2)-F(1)) 173.01 WRITE(6,76)N(K),E(1),EP12,F(2),OMEGA,E(3),F(3),DX1,DY1 173.1 76 FORMAT(3X,A8,8X,E14.6,9X,E14.6,5X,E14.6,4X,E14.6, 173.2 +1X,2F10.4,1X,2F10.4/) 173.4 C 176 C COMPUTE OFF-DIAGONAL ELEMENT OF THE SYMMETRIC STRAIN TENSOR 177 C 178 C 183 C WRITE THE ELEMENTS OF SYMMETRIC STRAIN TENSOR IN A FILE 184 C WITH UNIT 17 185 C 186 WRITE(17,455)E(1),EP12,F(2) 187 455 FORMAT(1X,3F18.10) 187.1 C 189 C 191 WRITE(11,122)N(K),X(K),Y(K) 192 122 FORMAT(A8,2F15.4) 193 C 194 C 208 C WRITE STATION NUMBERS DISPLACEMENTS AND VALUES OF 209 C AVERAGE DIFFERENTIAL ROTATION ON A FILE I N UNIT 9 210 C 211 WRITE(9,135)N(K),DX1,DY1,OMEGA 212 135 FORMAT(A8,2F10.4,F16.10) 213 IF(K.EQ.ICOUNT)STOP 227 6 CONTINUE 228 K=K+1 229 JJ=1 230 ST(1)=N1 231 SX=X(K) 232 SY=Y(K) 233 C 234 C COMPUTE THE DISPLACEMENTS AND ELEMENTS OF "A" MATRIX 235 C FOR THE STATION AT WHICH THE STRAIN ELLIPSE IS BEING COMPUTED 236 C 237 DX(1)=X2(K)-X(K) 238 DY(1)=Y2(K)-Y(K) 239 A(1,1)=0.0 240 A(1,2)=0.0 241 A(1,3)=1.0 242 4 CONTINUE 243 C CHECK IF REPEATED CONNECTION 243.5 DO 10 I=1,JJ 244 10 IF(N2.EQ.ST(I))GO TO 30 245 C WRITE CONNECTION IN FILE 8 245.1 47 IF(LSOBS)GOTO 48 245.2 WRITE(8,100)M5,N1,N2 245.3 48 CONTINUE 245.4 JJ=I 246 J=JJ 247 ST(I)=N2 248 DO 5 I=1,ICOUNT 249 IF(N2.NE.N3(I))GOTO 5 250 C 251 C COMPUTE THE ELEMENTS OF 'A'MATRIX AND THE DISPLACEMENTS 252 C FOR STATIONS WHICH ARE CONNECTED BY OBSERVATIONS TO THE POINT 253 C AT WHICH VALUE OF AVE. DIFF. ROTATION AND STRAIN ELLIPSE 254 C ARE BEING COMPUTED 255 C 256 DX(J)=X2(I)-X(I) 257 DY(J)=Y2(I)-Y(I) 258 A(J,1)=X(I)-SX 259 A(J,2)=Y(I)-SY 260 A(J,3)=1.0 261 5 CONTINUE 262 GOTO 30 264 END 266 C * ******************************************************************267 C * *268 C * *269 C * SUBROUTINE MTVMUL *270 C * ----------------- *271 C * *272 C * *273 C * THIS SUBROUTINE MULTIPLIES A MATRIX OF DIMENSION M BY N WITH *274 C * A COLUMN VECTOR OF N BY 1 *275 C * A AND B ARE THE MATRIX AND A COLUMN VECTOR TO BE MULTI- *276 C * PLIED RESPECTIVELY *277 C * C IS THE RESULTANT VECTOR. *278 C * IA,IB,IC ARE THE ROW DIMENSIONS OF A AND COLUMN DIMENSIONS *279 C * OF B AND C RESPECTIVELY *280 C * M BY N IS THE SIZE OF A *281 C * *282 C * *283 C * INPUT: *284 C * *285 C * MATRIX A AND COLUMN VECTOR B *286 C * *287 C * *288 C * OUTPUT: *289 C * *290 C * COLUMN VECTOR C *291 C * *292 C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. *293 C * *294 C * *295 C * ******************************************************************296 C 297 C 298 SUBROUTINE MTVMUL(A,IA,B,IB,C,IC,M,N) 299 DIMENSION A(IA,N),B(IB),C(IC) 300 DO 10 I=1,M 301 C(I)=0. 302 DO 10 J=1,N 303 10 C(I)=A(I,J)*B(J)+C(I) 304 RETURN 305 END 306 C * ******************************************************************307 C * *308 C * *309 C * SUBROUTINE MATMUL *310 C * ----------------- *311 C * *312 C * *313 C * THIS SUBROUTINE MULTIPLES TWO MATRICES A AND B AND C IS THE *314 C * RESULTING MATRIX *315 C * *316 C * IA,IB,IC ARE THE ROW DIMENSIONS OF A,B AND C RESPECTIVELY *317 C * L BY M IS THE SIZE OF A *318 C * M BY N IS THE SIZE OF B *319 C * L BY N IS THE SIZE OF C *320 C * *321 C * *322 C * INPUT: *323 C * *324 C * MATRICES A AND B *325 C * *326 C * *327 C * OUTPUT: *328 C * *329 C * MATRIX C *330 C * *331 C * WRITTEN BY KHANGENDRA THAPA NOVEMBER, 1979. *332 C * *333 C * *334 C * ******************************************************************335 C 336 C 337 SUBROUTINE MATMUL(A,IA,B,IB,C,IC,L,M,N) 338 DIMENSION A(IA,M),B(IB,N),C(IC,N) 340 DO 2 I=1,L 341 DO 2 J=1,N 342 C(I,J)=0.0E0 343 DO 1 K=1,M 344 C(I,J)=C(I,J)+A(I,K)*B(K,J) 345 1 CONTINUE 346 2 CONTINUE 347 RETURN 348 END 349 C * ******************************************************************350 C * *351 C * *352 C * *353 C * SUBROUTINE MATTNP *354 C * ----------------- *355 C * *356 C * *357 C * THIS SUBROUTINE TRANSPOSES MATRIX A TO MATRIX B *358 C * *359 C * IDA AND IDB ARE THE ROW DIMENSIONS OF A AND B RESPECTIVELY *360 C * M BY N IS THE SIZE OF A *361 C * N BY M IS THE SIZE OF B *362 C * *363 C * *364 C * INPUT: *365 C * *366 C * MATRIX A *367 C * *368 C * *369 C * *370 C * OUTPUT: *371 C * *372 C * MATRIX B *373 C * *374 C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. *375 C * *376 C * *377 C * ******************************************************************378 C 379 C 380 SUBROUTINE MATTNP (A,IDA,B,IDB,M,N) 381 DIMENSION A(IDA,N),B(IDB,M) 383 DO 1 I=1,M 384 DO 2 J=1,N 385 B(J,I)=A(I,J) 386 2 CONTINUE 387 1 CONTINUE 388 RETURN 389 END 390 C * ******************************************************************391 C * *392 C * *393 C * SUBROUTINE CHOLD *394 C * ---------------- *395 C * *396 C * *397 C * THIS SUBROUTINE DOES MATRIX INVERSION USING CHCLESKI *398 C * DECOMPOSITION *399 C * *400 C * *401 C * INPUT: *402 C * *403 C * A=POSITIVE DEFINITE SYMMETRIC MATRIX *404 C * IRDA=ROW DIMENSION OF MATRIX A *405 C * NA=SIZE OF MATRIX A *406 C * *407 C * *408 C * OUTPUT: *409 C * *410 C * A=INVERSE OF INPUT MATRIX(INPUT DESTROYED *411 C * DETA=DETERMINANT OF MATRIX A *412 C * *413 C * *414 C * ******************************************************************415 C 416 C 417 SUBROUTINE CHOLD(A,IRDA,NA,DETA) 418 C 419 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION 420 C 421 C INPUT ARGUMENTS 422 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX 423 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX 424 C NA = SIZE OF INPUT MATRIX 425 C OUTPUT ARGUMENTS 426 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) 427 C DETA = DETERMINANT OF INPUT MATRIX 428 C * = ERROR RETURN (TAKEN IF NA .T. 1 OR IF DETA .LT. SING) 429 C 430 DIMENSION A(IRDA,NA) 432 DATA SING/1E-10/ 435 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 436 IF(NA .LT. 1) GO TO 18 437 DETA = A(1,1) 438 A(1,1)=SQRT(DETA) 439 IF(NA .EQ. 1) GO TO 6 440 DO 1 I = 2,NA 441 IF(A(1,1).LE.SING)STOP "DIVISION BY ZERO" 441.1 1 A(I,1) = A(I,1) / A(1,1) 442 DO 5 J = 2,NA 443 SUM = 0. 444 J1 = J - 1 445 DO 2 K = 1,J1 446 2 SUM = SUM + A(J,K) ** 2 447 DETA = DETA * (A(J,J) - SUM) 448 A(J,J) = SQRT(A(J,J) - SUM) 449 IF(J .EQ. NA) GO TO 5 450 J2 = J + 1 451 DO 4 I = J2,NA 452 SUM = 0. 453 DO 3 K = 1,J1 454 3 SUM = SUM + A(I,K) * A(J,K) 455 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 456 5 CONTINUE 457 6 IF(ABS(DETA) .LT. SING) GO TO 16 458 C INVERSION OF LOWER TRIANGULAR MATRIX 459 DO 7 I = 1,NA 460 7 A(I,I) = 1. / A(I,I) 461 IF(NA .EQ. 1) GO TO 10 462 N1 = NA - 1 463 DO 9 J = 1,N1 464 J2 = J + 1 465 DO 9 I = J2,NA 466 SUM = 0. 467 I1 = I - 1 468 DO 8 K = J,I1 469 8 SUM = SUM + A(I,K) * A(K,J) 470 9 A(I,J) = - A(I,I) * SUM 471 C CONSTRUCTION OF INVERSE OF INPUT MATRIX 472 10 DO 15 J = 1,NA 473 IF(J .EQ. 1) GO TO 12 474 J1 = J - 1 475 DO 11 I = 1,J1 476 11 A(I,J) = A(J,I) 477 12 DO 14 I = J,NA 478 SUM = 0. 479 DO 13 K = I,NA 480 13 SUM = SUM + A(K,I) * A(K,J) 481 14 A(I,J) = SUM 482 15 CONTINUE 483 RETURN 484 16 WRITE(6,17) DETA 485 17 FORMAT(10X, "SINGULAR MATRIX IN CHOILD. DET =",E20.5) 486 REWIND 1 487 18 WRITE(6,19) 488 19 FORMAT(10X,"MATRIX OF DIMENSION ZERO IN CHOLD") 489 REWIND 1 490 END 491 C * ******************************************************************496 C * *497 C * *498 C * *499 C * PROGRAM EVALUE *500 C * -------------- *501 C * *502 C * *503 C * THIS PROGRAM COMPUTES SEMI-MAJOR EXIS,SEMI-MINOR AXIS AND *504 C * ORIENTATION OF SEMI-MAJOR AXIS OF THE STRAIN ELLIPSE. *505 C * *506 C * *507 C * JOB CONTROL LANGUAGE *508 C * -------------------- *509 C * *510 C * *511 C * 1.0 INPUT DATA SETS *512 C * *513 C * STATION NUMBERS AND CO-ORDINATES *514 C * *515 C * //GO.FT01F001 DD DSN=A.M6068.ADJCRD31.DTAA, *516 C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR *517 C * *518 C * ELEMENTS OF SYMMETRIC STRAIN TENSOR MATRIX *519 C * //GC.FT15F001 DD DSN=A.M6068.TMATRIX2.DATA, *520 C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR *521 C * *522 C * *523 C * 2.0 OUTPUT DATA SETS *524 C * *525 C * AXES AND ORIENTATION OF STRAIN ELLIPSES *526 C * *527 C * //GO.FT31F001 DD DSN=A.M6068.TEVALUE2.DATA, *528 C * // UNIT=DASD,VOL=SER=USER11,DISP=(NREW,CATLG,DELETE) *529 C * // DCB=(RECFM=FB,LRECL=33,BLKSIZE=3300),SPACE=(TRK,(10,10)RLSE) *530 C * //GO.SYSIN DD * *531 C * // *532 C * *533 C * *534 C * INPUT: *535 C * *536 C * CARD 1 FORMAT(I4) *537 C * NS=NUMBER OF STATIONS IN THE NETWORK *538 C * *539 C * ALL OTHER INPUT IS DONE FROM FILES *540 C * *541 C * N1,X,Y STATION NAMES AND CO-ORDINATES OF POINTS FROM EITHER *542 C * OF THE TWO ADJUSTMENTS. *543 C * FORMAT(I8,2F15.4) *544 C * *545 C * A= SYMMETRIC STRAIN TENSOR MATRIX IN THE FORM 1ST DIAGONAL *546 C * ELEMENT,OFF-DIAGONAL ELEMENT AND SECOND DIAGONAL ELEMENT *547 C * (DIRECT OUTPUT OF PROGRAM 'STRAIN') *548 C * FORMAT(1X,3F18.10) *549 C * *550 C * *551 C * OUTPUT: *552 C * D(1) AND D(2) ARE THE TWO AXES OF OF THE STRAIN ELLIPSES, *553 C * THAT IS, THE EIGEN VALUES OF STRAIN TENSOR MATRIX. *554 C * PHI=ORIENTATION ANGLE OF THE AXES OF STRAIN ELLIPSE *555 C * FORMAT(1X,2F16.10,F10.6) *556 C * (USED IN PROGRAM 'NETPLOT') *557 C * *558 C * WK=PERFORMANCE INDEX OF THE SUBROUTINE EIGRS. *559 C * IER=GIVES INDICATION OF ANY ERROR. *560 C * *561 C * *562 C * REFERENCES: *563 C * *564 C * IMS MANUAL *565 C * *566 C * *567 C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. *568 C * *569 C * *570 C * ******************************************************************571 C 572 C 573 PROGRAM EVALUE(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,TAPE11, 575 +TAPE17,TAPE3) 575.1 DIMENSION A(3),D(2),NETWORK(8),P(2) 576 C 576.1 C READ TITEL 577 READ(11,131)NETWORK 578 131 FORMAT(8A10) 579 C 580 C READ NUMBER OF STATION NUMBERS 581 C 582 READ(11,126)NS 583 126 FORMAT(I4) 584 WRITE(6,132)NS 584.01 132 FORMAT(1H1,3X,"NUMBER OF STATIONS:",I4//) 584.02 WRITE(6,25) 584.1 25 FORMAT(3X,"NUMBER",7X,"MAJOR AXIS",6X,"MINOR AXIS",5X, 584.2 +"ORIENTATION ANGLE"/16X,"MICROSTRAIN",5X,"MICROSTRAIN", 584.3 +10X,"RAD"/) 584.4 DO 77 I=1,NS 585 C 586 C READ STATION NUMBER AND CO-ORDINATES OBTAINED FROM 587 C EITHER OF THE TWO ADJUSTMENTS 588 C 589 READ(11,100)N1,X,Y 590 C 591 C READ ELEMENTS OF STRAIN TENSOR MATRIX IN UNIT 17(THIS IS 592 C AN OUTPUT OF PROGRAM STRAIN) 593 C 594 READ(17,455)A 595 455 FORMAT(1X,3F18.10) 596 C 597 C COMPUTE THE EIGEN VALUES AND 598 C EIGENVECTORS OF THE STRAIN TENSOR 599 C 600 AM=A(1)-A(3) 601 AM2=AM*AM 601.1 AA13=A(1)+A(3) 601.2 A2=A(2)*A(2) 601.3 DIETER=SQRT(AM2+4*A2) 601.4 D(1)=(AA13+DIETER)/2. 601.5 D(2)=(AA13-DIETER)/2. 601.6 100 FORMAT(A8,2F15.4) 602 P(1)=ABS(D(1)) 603 P(2)=ABS(D(2)) 604 IF(P(1) .EQ.0.0.OR.P(2) .EQ.0.0)GOTO 99 609 C 610 C DECIDE WHICH OF THE EIGEN VALUES IS GREATER 611 C 612 IF(A(2).EQ.0.0) GOTO 99 612.1 IF(P(1).GT.P(2)) GOTO 29 613 C 614 C COMPUTE ORIENTATION ANGLE FOR THE SEMI-MAJOR AXIS OF 615 C THE STRAIN ELLIPSE 616 C 617 ANGL=-(A(1)-D(2))/A(2) 618 PHI=ATAN(ANGL) 619 GOTO 69 620 29 ANGL=-(A(1)-D(1))/A(2) 621 PHI=ATAN(ANGL) 622 GOTO 69 623 99 PHI=0.0 624 69 CONTINUE 625 D(1)=D(1)*1.0E6 626 D(2)=D(2)*1.0E6 627 IF(P(2).GT.P(1)) GOTO 177 628 C 629 C WRITE THE VALUES OF AXES OF STRAIN ELLIPSE AND ORIENTATION ANGLE 630 C ON FILE IN UNIT 31 AS IN INPUT FOR PROGRAM NETPLOT 631 C 632 WRITE(3,140)(D(J),J=1,2),PHI 633 WRITE(6,141)N1,(D(J),J=1,2),PHI 634 GOTO 48 635 177 WRITE(3,140)D(2),D(1),PHI 636 WRITE(6,141)N1,D(2),D(1),PHI 637 48 CONTINUE 638 140 FORMAT(1X,2F16.10,F10.6) 639 141 FORMAT(2X,A8,2(3X,2F15.8),5X,F10.6) 641 77 CONTINUE 650 STOP 651 END 652 C * ******************************************************************655 C * *656 C * *657 C * PROGRAM NETPLOT *658 C * --------------- *659 C * *660 C * *661 C * THIS PROGRAM PLOTS THE STRAIN ELLIPSES AS WELL AS THE VALUES *662 C * OF AVERAGE DIFFERENTIAL ROTATION FOR NETWORK STATIONS *663 C * NETWORK POINTS ARE IDENTIFIED BY PLOTTING STATION NUMBERS *664 C * *665 C * *666 C * JOB CONTROL LANGUAGE *667 C * -------------------- *668 C * *669 C * *670 C * 1.0 INPUT DATA SETS *671 C * *672 C * STATION NUMBERS,VALUES OF DISPLACEMENTS AND VALUES OF AVERAGE *673 C * DIFFERENTIAL ROTATION(DIRECT OUTPUT OF PROGRAM 'STRAIN') *674 C * *675 C * //GO.FT08F001 DD DSN=A.M60689STNDXDY2.OMEGA, *676 C * // UNIT=DASD,VOL=SER=USER11,DISP=SHR *677 C * *678 C * ELEMENTS OF THE SYMMETRIC STRAIN TENSOR MATRIX *679 C * (DIRECT OUTPUT OF PROGRAM 'STRAIN') *680 C * *681 C * //GO.FT04F001 DD DSN=A.M6068.TMATRIX2.DATA, *682 C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR *683 C * *684 C * VALUES OF THE AXES OF STRAIN ELLIPSES AND THE ORIENTATION *685 C * ANGLE OF THE SEMI-MAJOR AXIS OF THE STRAIN ELLIPSE *686 C * (DIRECT OUTPUT OF PROGRAM 'EVALUE') *687 C * *688 C * //GO.FT03F001 DD DSN=A.M6068.TEVALUE2.DATA, *689 C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR *690 C * *691 C * STATION NUMBERS AND CO-ORDINATES *692 C * *693 C * //GO.FT02F001 DD DSN=A.M6068.ADJCRD31.DATA, *694 C * // UNIT=P3350,VOL=SER=USER11,DISP=SHR *695 C * //GO.SYSIN DD * *696 C * // *697 C * INPUT: *698 C * *699 C * CARD 1 FORMAT(I4) *700 C * NS=NUMBER OF STATIONS IN THE NETWORK *701 C * *702 C * ALL OTHER INPUT IS DONE FROM FILES *703 C * *704 C * ISTN,W,Z STATION NUMBER AND CO-ORDINATES OF THE NETWORK POINTS *705 C * FORMAT(1X,7A1,2F15.4) *706 C * *707 C * RMAJ,RMIN,PHI ARE SEMI-MAJOR, SEMI-MINOR AXES AND ORIENTATION *708 C * ANGLE OF SEMI-MAJOR AXIS OF STRAIN ELLIPSE RESPECTIVELY *709 C * FORMAT(1X,2F16.10,F10.6) *710 C * *711 C * ITON,DX1,DY1 AND OMEGA ARE STATION NUMBER DISPLACEMENTS *712 C * IN X AND Y DIRECTION AND THE VALUES OF AVERAGE DIFFERENTIAL *713 C * ROTATION RESPECTIVELY *714 C * FORMAT(1X,A1,6A1,2F10.4,F16.10) *715 C * *716 C * A CONTAINS THE ELEMENTS OF THE SYMMETRIC STRAIN TENSOR MATRIX *717 C * ORIENTATION ANGLE PHI AND VALUES OF AVERAGE DIFF. ROTATION *718 C * FORMAT(1X,3F18.10) *719 C * *720 C * *721 C * OUTPUT: *722 C * *723 C * ORIENTATION ANGLE PHI OF SEMI MAJOR AXIS OF STRAIN ELLIPSE *724 C * AND VALUES OF AVERAGE DIFFERENTIAL ROTATION *725 C * FORMAT(7X,7A1,2F16.8/) *726 C * *727 C * *728 C * REFERENCES: *729 C * *730 C * GUJAR,U.G.(1975).USERS GUIDE.COMPUTER PLOTTING VOLUME 11, *731 C * UNIVERSITY OF NEW BRUNSWICK COMPUTER CENTRE. *732 C * *733 C * *734 C * WRITTEN BY KHAGENDRA THAPA NOVEMBER,1979. *735 C * *736 C * *737 C * ******************************************************************738 C 739 C 740 PROGRAM NETPLOT(INPUT,OUTPUT,TAPE7=INPUT,TAPE11,TAPE9,TAPE3, 740.1 +TAPE12,TAPE6=OUTPUT,TAPE8) 740.2 LOGICAL NOELLPS,LDXY,LCON,NODXY 741 DIMENSION A(3),W(800),Z(800),N(800),TITLE(8) 744 C MARGINS=MAX.ELLIPSE DIM.›CM!,X-,Y-EXTEND ›IN! 745 DATA RMARG,XMAX,YMAX/5.,300.,33./ 745.1 C 745.2 PI=4.*ATAN(1.) 745.3 LDXY=.FALSE. 745.31 LCON=.FALSE. 745.32 C 745.4 READ(11,105)TITLE 746 105 FORMAT(8A10) 747 WRITE(6,106)TITLE 747.1 106 FORMAT("1",3X,8A10) 747.2 C 748 C READ THE NUMBER OF STATION NUMBERS 749 READ(11,111)NS 749.2 111 FORMAT(I4) 749.3 C 749.4 C READ PLOTPARAMETERS 749.5 C 750 READ(7,201)RADIUS,SCALE,SCALEL,SCALROT,SCALDXY,ICON 751 201 FORMAT(5F10.3,I2) 752 WRITE(6,202)NS,RADIUS,SCALE,SCALEL,SCALROT 752.1 202 FORMAT(3X,"STATION NUMBER=",I4,3X,"RADIUS=",F5.2,3X,"SCALE=1:", 752.2 +F8.1/3X,"SCALE OF ELLIPSES=",F7.3,"CM PER MICROSTRAIN",3X, 752.3 +"SCALE OF ROTATIONS",F7.3,"DEG. PER 0.1SEC"/) 752.4 C 753 23 IF(SCALDXY.EQ.0.)GOTO 24 754 LDXY=.TRUE. 754.1 WRITE(6,203)SCALDXY 754.2 203 FORMAT(3X,"SCALE OF DISPLACEMENTS=1:",F7.3/) 754.3 24 CONTINUE 754.4 25 IF(ICON.EQ.0)GOTO 26 754.5 LCON=.TRUE. 754.6 WRITE(6,204) 754.7 204 FORMAT(3X,"OPTION CON: CONNECTIONS USED FOR STRAIN COMPUTATION 754.8 +TO BE PLOTTED"/) 754.9 26 CONTINUE 755 C 755.1 CALL PLOTS(BUFF,1) 756 C 757 C SUBROUTINE PLOT IS USED TO SHIFT THE PAPER TO PROPER SETTING 759.1 CALL PLOT(0.,0.,-3) 760 WRITE(12,222) 760.1 222 FORMAT(9X,"X-COORDINATE",8X,"Y-COORDINATE",6X,"X IN INCHES", 760.2 +,3X,"Y IN INCHES"/) 760.3 C 760.4 C READ STATION NUMBERS AND CO-ORDINATES OF STATIONS 761.4 C COMPUTE PERIMETER 761.41 C 761.42 DO 22 I=1,NS 761.5 READ(11,450)N(I),W(I),Z(I) 762 450 FORMAT(A8,2F15.4) 763 IF(I.GT.1) GOTO 21 763.1 WMIN=W(1) 764.1 ZMIN=Z(1) 764.2 WMAX=W(1) 764.3 ZMAX=Z(1) 764.4 21 CONTINUE 765 IF(WMIN.GT.W(I)) WMIN=W(I) 766 IF(ZMIN.GT.Z(I)) ZMIN=Z(I) 767 IF(WMAX.LT.W(I)) WMAX=W(I) 767.1 IF(ZMAX.LT.Z(I)) ZMAX=Z(I) 767.2 22 CONTINUE 768.1 RMARG=RMARG/2.54 768.2 ZM1=0. 768.3 WM1=5./2.54 768.31 W1=WM1+RMARG 768.4 Z1=ZM1+RMARG 768.5 3 DW=(WMAX-WMIN)*100./2.54 768.6 WM2=W1+DW/SCALE+RMARG 768.65 DZ=(ZMAX-ZMIN)*100./2.54 768.7 ZM2=Z1+DZ/SCALE+RMARG 768.75 C 768.771 C CHECK PLOT SIZE 768.772 C 768.773 1 IF(WM2.LT.XMAX.OR.ZM2.LT.YMAX)GOTO2 768.774 SCALE1=DW/(XMAX-2.*RMARG-WM1) 768.775 SCALE2=DZ/(YMAX-2.*RMARG-ZM1) 768.776 SCALE=SCALE1 768.777 IF(SCALE2.GT.SCALE1)SCALE=SCALE2 768.778 WRITE(6,451)SCALE 768.779 451 FORMAT(3X,"COMPUTED SCALE=1:",F8.1/) 768.78 GOTO 3 768.781 2 CONTINUE 768.782 WRITE(6,312) 768.7820 312 FORMAT(2X,"NUMBER",8X,"MAJOR",8X,"MINOR",8X,"OMEGA",8X,"PHI"/ 768.7820 +18X,"CM",11X,"CM",11X,"DEG",8X,"DEG"/) 768.7820 C 768.7820 C WRITE TITLE AND LEGEND 768.7821 CALL SYMBOL(0.21,0.,0.2,"NETWORK: ",90.,9) 768.7822 CALL SYMBOL(0.21,2.,0.2,TITLE,90.,80) 768.7822 CALL SYMBOL(0.5,0.,0.15,"MAP SCALE 1:",90.,12) 768.7823 CALL NUMBER(0.5,2.,0.15,SCALE,90.,-1) 768.7824 CALL SYMBOL(0.8,0.,0.15,"ELLIPSE SCALE: CM PER MICROSTRAIN", 768.7824 +90.,36) 768.7825 CALL NUMBER(0.8,2.2,0.15,SCALEL,90.,1) 768.7826 CALL SYMBOL(1.1,0.,0.15,"SCALE OF ROTATIONS: DEG. PER 0.1 SEC."768.7827 +,90.,40) 768.7827 CALL NUMBER(1.1,2.9,0.15,SCALROT,90.,1) 768.7828 C 768.783 31 IF(.NOT.LDXY)GOTO 32 768.7831 CALL SYMBOL(0.8,8.,0.15,"SCALE OF DISPLACEMENTS 1:",90.,25) 768.7832 CALL NUMBER(0.8,12.,0.15,SCALDXY,90.,1) 768.7833 32 CONTINUE 768.7834 C 768.7835 C PLOT PERIMETER 768.784 C 768.785 CALL PLOT(WM1,ZM1,3) 768.8 CALL PLOT(WM2,ZM1,2) 768.81 CALL PLOT(WM2,ZM2,2) 768.82 CALL PLOT(WM1,ZM2,2) 768.83 CALL PLOT(WM1,ZM1,2) 768.84 CALL PLOT(WM1,ZM1,3) 768.85 C 784 C READ MAJOR AND MINOR AXES OF STRAIN ELLIPSES AND THE 785 C ORIENTATION ANGLE OF THE SEMI-MAJOR AXES 786 C 787 98 DO 99 J=1,NS 787.1 NOELLPS=.FALSE. 787.2 NODXY=.FALSE. 787.3 READ(3,150)RMAJOR,RMINOR,PHI 788 150 FORMAT(1X,2F16.10,F10.6) 789 C 794 C READ STATION NUMBER DISPLACEMENTS OF POINTS AND 795 C THE VALUES OF AVERAGE DIFFERENTIAL ROTATION 796 C 797 READ(9,415)N(J),DX,DY,OMEGA 798 415 FORMAT(A8,2F10.4,F16.10) 799 RMAJ=RMAJOR*SCALEL 799.1 RMIN=RMINOR*SCALEL 799.2 RMAJOR=RMAJ/2.54 799.3 RMINOR=RMIN/2.54 799.4 C 799.41 DX=DX*100./2.54/SCALDXY 799.42 DY=DY*100./2.54/SCALDXY 799.43 C 799.5 C CHECK ELLIPSE SIZE 799.6 IF(ABS(RMAJOR).GT.RMARG)NOELLPS=.TRUE. 799.7 C 799.8 C CHECK DISPLACEMENTS 799.9 IF(ABS(DX).GT.RMARG.OR.ABS(DY).GT.RMARG)NODXY=.TRUE. 799.91 C 805 C COMPUTE THE COORDINATES AT THE ENDS OF EACH AXIS 806 C OF THE STRAIN ELLIPSES 807 C 808 COSPHI=COS(PHI) 810 SINPHI=SIN(PHI) 811 RMAJC=RMAJOR*COSPHI 811.1 RMAJS=RMAJOR*SINPHI 811.2 RMINC=RMINOR*COSPHI 811.3 RMINS=RMINOR*SINPHI 811.4 X=(W(J)-WMIN) 815.1 Y=(Z(J)-ZMIN) 815.2 X=X*100./2.54/SCALE+W1 815.3 Y=Y*100./2.54/SCALE+Z1 815.4 X1=X+RMAJC 816 Y1=Y+RMAJS 817 X2=X-RMAJC 818 Y2=Y-RMAJS 819 X3=X-RMINS 820 Y3=Y+RMINC 821 X4=X+RMINS 822 Y4=Y-RMINC 823 RMAJABS=ABS(RMAJOR) 824 RMINABS=ABS(RMINOR) 825 C 835 C COMPUTE THE CO-ORDINATES OF THE CORNERS OF THE SEGMENT 836 C * OF THE CIRCLE *837 C 838 OMEGA=OMEGA*SCALROT*36000. 838.5 X5=X+RADIUS 839 Y5=Y 840 X6=X+RADIUS*COS(OMEGA) 841 Y6=Y+RADIUS*SIN(OMEGA) 842 X7=X6+.1 843 DOMEGA=OMEGA*180.0/PI 843.1 DPHI=PHI*180.0/PI 843.15 WRITE(6,155)N(J),RMAJ,RMIN,DOMEGA,DPHI 843.2 155 FORMAT(/2X,A8,4(3X,F10.5)) 843.3 IF(NOELLPS)WRITE(6,156) 843.31 156 FORMAT(1H+,64X,"ELLIPSE PLOTTING SUPPRESSED"/) 843.32 WRITE(12,444)W(J),Z(J),X,Y 843.4 444 FORMAT(3X,2(3X,F15.4),2(3X,F12.3)) 843.5 C 843.6 C COMPUTE COOORDINATES OF THE ENDS OF THE DISPLACEMENT VECTORS 843.61 C 843.62 11 IF(.NOT.LDXY.OR.NODXY)GOTO 12 843.621 X8=X+DX 843.63 Y8=Y+DY 843.64 12 CONTINUE 843.65 C 844.8 C SUBROUTINE PLOT ENABLES ONE TO DRAW POINTS AND ST. LINES 845 C 846 IF(NOELLPS)GOTO 83 846.05 CALL NEWPEN(1) 846.1 CALL PLOT(X1,Y1,3) 847 C 848 C SUBROUTINE DASHP ENABLES ONE TO PLOT DASHLINES 849 C 850 C PLOT DASH LINES FOR NEGATIVE VALUES OF STRAIN ELLIPSES 851 C 852 IF(RMAJOR.LT.0.) CALL DASHP(X2,Y2,0.07) 853 C 854 C PLOT SOLID LINES FOR POSITIVE VALUES OF AXES OF STRAIN ELLIPSES 855 C 856 IF(RMAJOR.GE.0.)CALL PLOT(X2,Y2,2) 857 CALL PLOT(X4,Y4,3) 858 IF(RMINOR.LT.0.0)CALL DASHP(X3,Y3,0.07) 859 IF(RMINOR.GT.0.0)CALL PLOT(X3,Y3,2) 860 CALL PLOT(X1,Y1,3) 860.1 CALL PLOT(X1,Y1,2) 860.2 C 861 C PLOT THE STRAIN ELLIPSE BY CALLING THE SUBROUTINE'ELIPS' 862 C 863 IF(RMAJOR.LT.0.0)GO TO 68 863.1 CALL ELIPS(X1,Y1,RMAJABS,RMINABS,DPHI,0.0,360.,2) 864 GOTO 83 864.1 68 CONTINUE 864.2 CALL PLOT(X2,Y2,3) 864.25 CALL PLOT(X2,Y2,2) 864.26 CALL ELIPS(X2,Y2,RMAJABS,RMINABS,DPHI,0.0,360.0,2) 864.3 83 CONTINUE 864.5 IF(OMEGA.EQ.0.0)GOTO 88 865 IF(OMEGA.LT.0.0)GOTO 73 866 C 867 C * PLOT THE ARC OF THE CIRCLE BY CALLING THE SUBROUTINE 'CIRCL' 868 C 869 CALL NEWPEN(2) 869.01 CALL PLOT(X5,Y5,3) 869.1 CALL PLOT(X5,Y5,2) 869.2 CALL CIRCL(X5,Y5,0.0,DOMEGA,RADIUS,RADIUS,0.0) 870 C 871 C PLOT SOLID LINES FOR THE RADII OF THE SEGMENT OF THE CIRCLE 872 C FOR POSITIVE VALUES OF OMEGA AND DRAW THE SEGMENT OF THE CIRCLE 873 C 874 CALL PLOT(X,Y,3) 875 CALL PLOT(X5,Y5,2) 876 CALL PLOT(X,Y,3) 877 CALL PLOT(X6,Y6,2) 878 GOTO 88 879 73 CONTINUE 879.1 CALL NEWPEN(2) 879.15 CALL PLOT(X5,Y5,3) 879.2 CALL PLOT(X5,Y5,2) 879.3 CALL CIRCL(X5,Y5,0.0,DOMEGA,RADIUS,RADIUS,0.0) 880 C 881 C PLOT DASH LINES FOR THE RADII OF THE SEGMENT OF THE CIRCLE 882 C FOR NEGATIVE VALUES OF OMEGA AND DRAW THE SEGMENT OF THE CIRCLE 883 C 884 CALL PLOT(X5,Y5,3) 884.1 CALL DASHP(X,Y,0.07) 885.1 CALL PLOT(X6,Y6,3) 885.2 CALL DASHP(X,Y,0.07) 887 88 CONTINUE 889 C 889.1 C OPTION DISP: PLOT DISPLACEMENTS 889.2 C 889.3 27 IF(.NOT.LDXY.OR.NODXY)GOTO 28 889.4 CALL NEWPEN(3) 889.5 CALL PLOT(X,Y,3) 889.6 CALL PLOT(X8,Y8,2) 889.7 28 CONTINUE 889.8 C 890 C PLOT STATION NUMBERS BY CALLING SUBROUTINE 'SYMBOL' 891 C 892 CALL NEWPEN(1) 892.5 CALL SYMBOL(X7,Y,0.08,N(J),0.0,8) 893 C 894 99 CONTINUE 898 C 898.1 C OPTION CON: PLOT CONNECTIONS 898.2 C 898.3 29 IF(.NOT.LCON)GOTO 30 898.4 CALL NEWPEN(4) 898.405 47 DO 48 K=1,2000 898.41 READ(8,501)M5,N1,N2 898.5 501 FORMAT(1X,I2,4X,A8,16X,A8) 898.6 IF(EOF(8).NE.0)GOTO 30 898.63 IF(M5.NE.0.AND.M5.NE.1.AND.M5.NE.2)GOTO 48 898.66 41 DO 42 I=1,NS 898.7 43 DO 44 J=1,NS 898.8 45 IF(I.EQ.J)GOTO 46 898.9 IF(N(I).NE.N1)GOTO 46 898.91 IF(N(J).NE.N2)GOTO 46 898.916 X=(W(I)-WMIN)/SCALE*100./2.54+W1 898.922 Y=(Z(I)-ZMIN)/SCALE*100./2.54+Z1 898.928 X9=(W(J)-WMIN)/SCALE*100./2.54+W1 898.934 Y9=(Z(J)-ZMIN)/SCALE*100./2.54+Z1 898.94 X9=(X+X9)/2. 898.946 Y9=(Y+Y9)/2. 898.952 CALL PLOT(X,Y,3) 898.958 CALL PLOT(X9,Y9,2) 898.964 46 CONTINUE 898.97 44 CONTINUE 898.976 42 CONTINUE 898.982 48 CONTINUE 898.983 30 CONTINUE 898.988 C 899 C SUBROUTINE PLOT IS CALLED TO INDICATE THE END OF ALL PLOTTING 900 C 901 CALL PLOT(0.,0.,999) 902 STOP 903 END 904