You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
simulacion-permeabilidad/tools/connec/connec-src/simpConnec/CONNEC2D.FOR

426 lines
15 KiB
Fortran

C CONNEC2D.FOR VER.21-FEB-2002
C
C PROGRAM FOR CONNECTIVITY ANALYSIS OF A 2D INDICATOR MAP
C
C GIVEN AN INDICATOR FIELDS (VALUES 0 AND 1 ONLY), THE RANDOM SET
C WITH VALUES 1 IS ANALYSED FOR CONNECTIVITY.
C
C INPUT PARAMETER FILE WITH
C
C IPHA : 0 OR 1 FOR CONNECTIVITY ANALYSIS OF PHASE 0 OR 1.
C ICON : 4 FOR 4-CONNECTIVITY AND 8 FOR 8-CONNECTIVITY
C CINP : INPUT FILE WITH INDICATOR VARIABLE (VALUS 0/1 ONLY)
C NX NY : NUMBER OF POINTS IN X AND Y
C DX DY : GRID DIMENSIONS IN X AND Y
C COUT : OUTPUT FILE WITH STATISTICS
C COU2 : OUTPUT FILE WITH CONNECTED COMPONENTS
C COU3 : OUTPUT FILE WITH CONNECTIVITY FUNCITON
C
C DX,DY ARE ONLY USED FOR CONVERTING GRID UNITS TO REAL UNITS.
C NX : NUMBER OF POINTS ALONG THE X DIRECTION
C NY : NUMBER OF POINTS ALONG THE Y DIRECTION
C
C IND(I,J) : 0/1 INDICATOR VALUE
C AT THE TERMINATION OF THE PROGRAM THE MATRIX IND(I,J) CONTAINS
C THE CONNECTED COMPONENTS WITH VALUES 1,2,3,4,5,... FOR FIRST,
C SECOND, THIRD, ... ETC CONNECTED COMMPONENTS.
C THE VALUE 0 REMAINS 0.
C COUT : OUTPUT FILE WITH STATISTICS AND CONNECTIVITY FUNCTION.
C COU2 : OUTPUT FILE WITH THE CONNECTED COMPONENTS.
C COU3 : OUTPUT FILE WITH CONNECTIVITY FUNCTION
C
C FOR THE CONNECTIVITY FUNCTION:
C NPX(K): NUMBER OF PAIRS OF VALUES SEPARATED A DISTANCE K ALONG
C THE X DIRECTION AND THAT BELONG TO FACIES 1.
C NPCX(K): NUMBER OF THE PREVIOUS VALUES THAT ARE CONNECTED.
C NPY(K): NUMBER OF PAIRS OF VALUES SEPARATED A DISTANCE K ALONG
C THE Y DIRECTION AND THAT BELONG TO FACIES 1.
C NPCY(K): NUMBER OF THE PREVIOUS VALUES THAT ARE CONNECTED.
C
C
C PPHA : PROPORTION OF FACIES 1.
C NCC : NUMBER OF CONNECTED COMPONENTS.
C RISM : MEAN CONNECTED COMPONENT SIZE IN PIXELS.
C RSME : MRAN SIZE IN REAL UNITS.
C RTOT : MEAN SIZE RELATIVE TO SIZE OF FACIES 1.
C RXME : MEAN LENGTH ALONG THE X DIRECTION.
C RYME : MEAN LENGTH ALONG THE Y DIRECTION.
C ICOM : NUMBER OF THE LARGEST COMPONENT.
C ISMA : SIZE IN PIXELS (OR LARGEST COMPONENT).
C RSIZ : SIZE RELATIVE TO SIZE OF FACIES 1.
C IXMA : MAXIMUM LENGTH ALONG X.
C IYMA : MAXIMUM LENGTH ALONG Y.
C ISMI : SIZE OF SMALLEST COMPONENT.
C IXMI : MINIMUM LENGTH ALONG X.
C IYMI : MINIMUM LENGTH ALONG Y.
C IPX : NUMBER OF PERCOLATING COMPONENTS ALONG X.
C IPY : NUMBER OF PERCOLATING COMPONENTS ALONG Y.
C CCFU(.,.) : CONNECTIVITY FUNCTION.
C
C CCFU(.,1): CONNECTIVITY FUNCTION ALONG X (E-W).
C CCFU(.,2): CONNECTIVITY FUNCTION ALONG Y (N-S).
C CCFU(.,3): CONNECTIVITY FUNCTION ALONG NE-SW DIRECTION.
C CCFU(.,4): CONNECTIVITY FUNCTION ALONG NW-SE DIRECTION.
C CCFU(.,5): MEAN CONNECTIVITY FUNCTION ALONG X AND Y.
C CCFU(.,6): MEAN CONNECTIVITY FUNCTION ALONG THE DIAGONALS.
C
C
C ----------------------------------------------------------------
C
DIMENSION IND(4096,4096)
C,NPX(256),NPY(256),NPCX(256),NPCY(256)
C DIMENSION NPE(256),NPW(256),NPCE(256),NPCW(256),CCFU(256,6)
C
CHARACTER*12 CPAR,CINP,COU2
C
C MATRIX DIMENSION LIMITATION
C
NXM=4096
NYM=4096
C
CPAR="coninput.txt"
OPEN (1,FILE=CPAR)
READ (1,*) IPHA
IF (IPHA.NE.0.AND.IPHA.NE.1) THEN
WRITE (6,*)'THE FIRST LINE IN THE PARAMETER FILE'
WRITE (6,*)'MUST BE A 0 OR A 1'
WRITE (6,*)'0 FOR CONNECTIVITY ANALYSIS OF PHASE 0'
WRITE (6,*)'1 FOR CONNECTIVITY ANALYSIS OF PHASE 1'
WRITE (6,*)'THE ACTUAL VALUE IN THE PARAMETER FILE IS: ',IPHA
STOP
END IF
READ (1,*) ICON
IF (ICON.NE.4.AND.ICON.NE.8) THEN
WRITE (6,*)'THE SECOND LINE IN THE PARAMETER FILE'
WRITE (6,*)'MUST BE A 4 OR A 8'
WRITE (6,*)'4 MEANS 4-CONNECTIVITY'
WRITE (6,*)'8 MEANS 8-CONNECTIVITY'
WRITE (6,*)'THE ACTUAL VALUE IN THE PARAMETER FILE IS: ',ICON
STOP
END IF
READ (1,100) CINP
READ (1,*) NX,NY
READ (1,*) DX,DY
READ (1,100) COU2
CLOSE (1)
IF (NX.GT.NXM) THEN
WRITE (6,*)'MAXIMUM NX IS ',NXM
WRITE (6,*)'ACTUAL VALUE IS ',NX
STOP
END IF
IF (NY.GT.NYM) THEN
WRITE (6,*)'MAXIMUM NY IS ',NYM
WRITE (6,*)'ACTUAL VALUE IS ',NY
STOP
END IF
C
C READING EXPERIMENTAL INDICATOR DATA 2D FIELD
C
OPEN (1,FILE=CINP)
JPHA=0
DO 1 I=1,NX
DO 2 J=1,NY
READ (1,*) VAL
IND(I,J)=INT(VAL)
IF (IPHA.EQ.0) THEN
IF (INT(VAL).EQ.0) THEN
IND(I,J)=1
ELSE
IND(I,J)=0
END IF
END IF
IF (IND(I,J).NE.1.AND.IND(I,J).NE.0) THEN
WRITE (6,*)'EXPERIMENTAL DATA MUST BE INDICATOR DATA 0/1'
WRITE (6,*)'ACTUAL VALUE : ',IND(I,J)
STOP
END IF
IF (IPHA.EQ.0) THEN
IF (INT(VAL).EQ.0) THEN
IND(I,J)=1
ELSE
IND(I,J)=0
END IF
END IF
IF (IND(I,J).EQ.1) JPHA=JPHA+1
2 CONTINUE
1 CONTINUE
CLOSE (1)
IARE=NX*NY
PPHA=FLOAT(JPHA)/IARE
WRITE (6,*)
IF (ICON.EQ.4) WRITE (6,*)'4-CONNECTIVITY ANALYSIS'
IF (ICON.EQ.8) WRITE (6,*)'8-CONNECTIVITY ANALYSIS'
WRITE (6,*)
WRITE (6,*)'PROPORTION OF PHASE 1 IS: ',PPHA
C
C
C LOOKING FOR CONNECTED COMPONENTS
C
NCC=0
C
5 DO 6 I=1,NX
DO 7 J=1,NY
IF (IND(I,J).EQ.1) THEN
IAI=I
IAJ=J
GOTO 8
END IF
7 CONTINUE
6 CONTINUE
C NO MORE CONNECTED COMPONENTS SENDING CONTROL TO LABEL 20
GOTO 20
C
C A NEW COMPONEND HAS BEEN FOUND, INCREASE NUMBER OF
C CONNECTED COMPONENTS BY 1
C
C
8 NCC=NCC+1
C
C LOOKING FOR COMPONENT NCC WITH LABEL NNC+1
C
IND(IAI,IAJ)=NCC+1
DO 9 I=IAI,NX
DO 10 J=IAJ,NY
IF ((I-1).GT.0) THEN
IF (IND(I-1,J).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I-1,J)=NCC+1
END IF
END IF
IF (I+1.LE.NX) THEN
IF (IND(I+1,J).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I+1,J)=NCC+1
END IF
END IF
IF (J-1.GT.0) THEN
IF (IND(I,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I,J-1)=NCC+1
END IF
END IF
IF (J+1.LE.NY) THEN
IF (IND(I,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I,J+1)=NCC+1
END IF
END IF
IF (ICON.EQ.8) THEN
IF ((I-1).GT.0.AND.(J-1).GT.0) THEN
IF (IND(I-1,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I-1,J-1)=NCC+1
END IF
END IF
IF ((I+1).LE.NX.AND.(J-1).GT.0) THEN
IF (IND(I+1,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I+1,J-1)=NCC+1
END IF
END IF
IF ((I-1).GT.0.AND.(J+1).LE.NY) THEN
IF (IND(I-1,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I-1,J+1)=NCC+1
END IF
END IF
IF ((I+1).LE.NX.AND.(J+1).LE.NY) THEN
IF (IND(I+1,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I+1,J+1)=NCC+1
END IF
END IF
END IF
10 CONTINUE
9 CONTINUE
DO 40 I=IAI,1,-1
DO 41 J=IAJ,1,-1
IF ((I-1).GT.0) THEN
IF (IND(I-1,J).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I-1,J)=NCC+1
END IF
END IF
IF (I+1.LE.NX) THEN
IF (IND(I+1,J).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I+1,J)=NCC+1
END IF
END IF
IF (J-1.GT.0) THEN
IF (IND(I,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I,J-1)=NCC+1
END IF
END IF
IF (J+1.LE.NY) THEN
IF (IND(I,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I,J+1)=NCC+1
END IF
END IF
IF (ICON.EQ.8) THEN
IF ((I-1).GT.0.AND.(J-1).GT.0) THEN
IF (IND(I-1,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I-1,J-1)=NCC+1
END IF
END IF
IF ((I+1).LE.NX.AND.(J-1).GT.0) THEN
IF (IND(I+1,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I+1,J-1)=NCC+1
END IF
END IF
IF ((I-1).GT.0.AND.(J+1).LE.NY) THEN
IF (IND(I-1,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I-1,J+1)=NCC+1
END IF
END IF
IF ((I+1).LE.NX.AND.(J+1).LE.NY) THEN
IF (IND(I+1,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) IND(I+1,J+1)=NCC+1
END IF
END IF
END IF
41 CONTINUE
40 CONTINUE
C
62 NNN=0
DO 60 I=1,NX
DO 61 J=1,NY
IF ((I-1).GT.0) THEN
IF (IND(I-1,J).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I-1,J)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF (I+1.LE.NX) THEN
IF (IND(I+1,J).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I+1,J)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF (J-1.GT.0) THEN
IF (IND(I,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I,J-1)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF (J+1.LE.NY) THEN
IF (IND(I,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I,J+1)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF (ICON.EQ.8) THEN
IF ((I-1).GT.0.AND.(J-1).GT.0) THEN
IF (IND(I-1,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I-1,J-1)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF ((I+1).LE.NX.AND.(J-1).GT.0) THEN
IF (IND(I+1,J-1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I+1,J-1)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF ((I-1).GT.0.AND.(J+1).LE.NY) THEN
IF (IND(I-1,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I-1,J+1)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
IF ((I+1).LE.NX.AND.(J+1).LE.NY) THEN
IF (IND(I+1,J+1).EQ.1) THEN
IF (IND(I,J).EQ.NCC+1) THEN
IND(I+1,J+1)=NCC+1
NNN=NNN+1
END IF
END IF
END IF
END IF
61 CONTINUE
60 CONTINUE
IF (NNN.GT.0) GOTO 62
GOTO 5
C
C THE CONNECTED COMPONENT IS THE COMPONENT WITH VALUES
C IA(I,J)=NCC, WE SUBSTRACT 1 FROM THE LEBEL
C
20 DO 50 I=1,NX
DO 51 J=1,NY
IF (IND(I,J).NE.0) IND(I,J)=IND(I,J)-1
51 CONTINUE
50 CONTINUE
C
C
C OUTPUT FILE WITH CONNECTED COMPONENTS
C
OPEN (1,FILE=COU2)
C OPEN (2,FILE=CINP)
C
DO 15 I=1,NX
DO 16 J=1,NY
C READ (2,*) ZVA
WRITE (1,*) IND(I,J)
16 CONTINUE
15 CONTINUE
C
CLOSE (1)
C
100 FORMAT (A)
110 FORMAT (1X,'INPUT INDICATOR FILE : ',A12)
120 FORMAT (1X,'NX : ',I6)
130 FORMAT (1X,'NY : ',I6)
140 FORMAT (1X,'DX : ',F10.3)
150 FORMAT (1X,'DY : ',F10.3)
160 FORMAT (1X,'PROPORTION OF FACIES ',I1,' : ',F10.4)
170 FORMAT (/1X,'CONNECTED COMPONENTS STATISTICS'//
* /1X,'NUMBER OF CONNECTED COMPONENTS : ',I6)
175 FORMAT (/1X,'AVERAGES'/1X,8('='))
176 FORMAT (/1X,'MINIMA'/1X,6('='))
177 FORMAT (/1X,'MAXIMA'/1X,6('='))
178 FORMAT (/1X,'PERCOLATION'/1X,11('='))
180 FORMAT (1X,'MEAN SIZE IN PIXELS : ',F10.4)
190 FORMAT (1X,'MEAN SIZE REAL UNITS : ',F12.4)
200 FORMAT (1X,'MEAN SIZE RELATIVE TO TOTAL AREA OF FACIES ',I1,
* ' : ',F10.4)
210 FORMAT (1X,'MEAN LENGTH ALONG X (IN PIXELS) : ',F10.4)
220 FORMAT (1X,'MINIMUM LENGTH ALONG X (IN PIXELS) : ',I6)
230 FORMAT (1X,'MAXIMUM LENGTH ALONG X (IN PIXELS) : ',I6)
235 FORMAT (1X,'THE LARGEST COMPONENT IS NUMBER : ',I6)
236 FORMAT (1X,'WITH MAXIMUM SIZE IN PIXELS : ',I6)
237 FORMAT (1X,'AND RELATIVE TO TOTAL AREA OF FACIES '
* ,I1,' : ',F10.4)
238 FORMAT (1X,'SIZE IN PIXELS OF SMALLEST COMPONENT : ',I6)
240 FORMAT (1X,'MEAN LENGTH ALONG Y (IN PIXELS) : ',F10.4)
250 FORMAT (1X,'MIMIMUN LENGTH ALONG Y (IN PIXELS) : ',I6)
260 FORMAT (1X,'MAXIMUM LENGTH ALONG Y (IN PIXELS) : ',I6)
270 FORMAT (1X,'NUMBER OF PERCOLATING COMPONENTS IN X : ',I6)
280 FORMAT (1X,'NUMBER OF PERCOLATING COMPONENTS IN Y : ',I6)
290 FORMAT (/1X,'OUTPUT FILE WITH CONNECTED COMPONENTS : ',A12)
295 FORMAT (/1X,'OUTPUT FILE WITH CONNECTIVIY FUNCTION : ',A12)
300 FORMAT (/1X,'CONNECTIVITY FUNCTION'/1X,21('=')
* //1X,'ALONG THE X DIRECTION')
310 FORMAT (F10.3,F12.6,2I12)
320 FORMAT (/1X,'ALONG THE Y DIRECTION')
321 FORMAT (/1X,'ALONG THE NE-SW DIRECTION')
322 FORMAT (/1X,'ALONG THE NW-SE DIRECTION')
323 FORMAT (/1X,'AVERAGE ALONG THE DIAGONALS')
330 FORMAT (/1X,'AVERAGE ALONG X AND Y')
400 FORMAT (25(/),1X,'CONNEC2D PROGRAM VER. 1.0'/1X,27('=')///
* 1X,'INPUT PARAMETER FILE ---> ',$)
410 FORMAT (/1X,'OUTPUT FILE WITH STATISTICS AND CONNECTIVITY '
* 'FUNCTION: ',A12)
420 FORMAT (/1X,'OUTPUT FILE WITH CONNECTED COMPONENTS: ',A12)
425 FORMAT (/1X,'OUTPUT FILE WITH CONNECTIVITY FUNCTION: ',A12)
430 FORMAT (1X,'OUTPUT RESULTS OF CONNEC2D'//1X'4-CONNECTIVITY'
* ' ANALYSIS'/1X,'PARAMETER FILE WAS : ',A12//)
440 FORMAT (1X,'OUTPUT RESULTS OF CONNEC2D'//1X'4-CONNECTIVITY'
* ' ANALYSIS'/1X,'PARAMETER FILE WAS : ',A12//)
END