C
C Program GRDCKN to compute correlation functions C C Version: 6.10 C Date: 2007, June 8 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C C The values of the Von Karman correlation function or of the C power-law correlatiopn function are calculated according to the C equations (48) or (66) of the paper: Klimes, L. (2002): Correlation C functions of random media. Pure and Applied Geophysics, 159, C 1811-1831. C C If the value of (NDIM/2+POWERN) equals zero, the Gamma function in C equations (48) or (66) can not be calculated, and the correlation C function becomes to be Dirac distribution realized by the product of C values of (1-|Xi-X0i|/|Di|)/|Di| in the points closer to the point C X0 than one grid interval, and by zeros in other points. C C Otherwise, if the value of ACOR equals to 999999., the power-law C correlation function according to the equation (66) is calculated. C C Otherwise, the Von Karman correlation function according to the C equation (66) is calculated. C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Data for the correlation function: C NDIM=positive integer ... Spatial dimension. C Default: NDIM equal to the dimension of the input grid. C KAPPA=real ... Multiplicative factor. C Default: KAPPA=1. C POWERN=real... Exponent or index related to fractal dimension: C Medium is self-affine at distances L: L .LT. ACOR C Reasonable values for geology: -0.5 .LT. POWERN .LT. 0.0 C Default: POWERN=0.0 C ACOR=positive real... Von Karman (large-scale) correlation length: C Suppresses large heterogeneities (larger than ACOR). C For infinit ACOR (ACOR=999999.), the power-law C correlation function is calculated. C For other values of ACOR, the Von Karman correlation C function is calculated. C Default: ACOR=999999. (infinity) C CKNMAX=real ... Maximum value of the correlation function. C Default: CKNMAX=999999. C Names of input and output formatted files: C PTS='string' ... Name of the file with coordinates of point X0 C in the form PTS. C Default: PTS=' ' means that the coordinates are [0.,0.,0]. C CKNOUT='string'... Name of the output file with the values C of the Von Karman correlation function in gridpoints. C Default: CKNOUT='ckn.out' C For general description of the files with gridded data refer to C file forms.htm. C Data specifying the parameters of the grid: C O1=real, O2=real, O3=real ... Coordinates of the origin of the C grid. C Default: O1=0. O2=0. O3=0. C D1=real... Grid interval along the X1 axis. C Default: D1=1. C D2=real... Grid interval along the X2 axis. C Default: D2=1. C D3=real... Grid interval along the X3 axis. C Default: D3=1. C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1 C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,GAMMLN,BESSIK REAL GAMMLN C ERROR ... File error.for. C RSEP1,RSEP3T,RSEP3R,RSEP3I ... C File sep.for. C GAMMLN ... File C gammln.for of Numerical Recipes. C BESSIK ... File C bessik.for of Numerical Recipes. C INTEGER NDIM REAL DIM,KAPPA,VN,ACOR,CKNMAX CHARACTER*80 FILSEP,FILEX0,FILOUT INTEGER LU1,LU2 PARAMETER (LU1=1,LU2=2) CHARACTER*3 TEXT INTEGER IGROUP,K1,K2,K3,N1,N2,N3,I1,I2,I3 INTEGER MX PARAMETER (MX=300) REAL X(MX,3),COOR(3) REAL PI PARAMETER (PI=3.1415926) REAL XX,GAMMA,CKN0,CKN,RI,RK,RIP,RKP,DIM2VN,ABSD1,ABSD2,ABSD3 REAL O1,O2,O3,D1,D2,D3 REAL X01,X02,X03 C----------------------------------------------------------------------- C C Reading a name of the file with the input data: FILSEP=' ' WRITE(*,'(A)') '+GRDCKN: Enter input filename: ' READ(*,*) FILSEP C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU1,FILSEP) ELSE C GRDCKN-01 CALL ERROR('GRDCKN-01: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C WRITE(*,'(A)') '+GRDCKN: Working ... ' C C Reading the file with coordinates of point X0: CALL RSEP3T('PTS',FILEX0,' ') IF (FILEX0.NE.' ') THEN OPEN(LU1,FILE=FILEX0,STATUS='OLD') READ(LU1,*) (TEXT,I1=1,20) READ(LU1,*) TEXT,X01,X02,X03 CLOSE(LU1) ELSE X01=0. X02=0. X03=0. ENDIF C Reading the values describing the grid: CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) ABSD1=ABS(D1) ABSD2=ABS(D2) ABSD3=ABS(D3) C Reading filename of the output file: CALL RSEP3T('CKNOUT',FILOUT,'ckn.out') C Reading numerical constants: NDIM=3 IF (N1.EQ.1) NDIM=NDIM-1 IF (N2.EQ.1) NDIM=NDIM-1 IF (N3.EQ.1) NDIM=NDIM-1 IF (NDIM.EQ.0) NDIM=NDIM+1 I1=NDIM CALL RSEP3I('NDIM',NDIM,I1) DIM=FLOAT(NDIM) CALL RSEP3R('KAPPA',KAPPA,1.) CALL RSEP3R('POWERN',VN,0.) CALL RSEP3R('ACOR',ACOR,999999.) CALL RSEP3R('CKNMAX',CKNMAX,999999.) IF (ACOR.LE.0.) THEN C GRDCKN-02 CALL ERROR('GRDCKN-02: ACOR less or equal zero.') ENDIF C DIM2VN=DIM/2.+VN IF (DIM2VN.NE.0.) THEN C Computing the value of the Gamma function for d/2+N: GAMMA=EXP(GAMMLN(DIM2VN)) C Computing the x-independent part of the correlation function: IF (ACOR.NE.999999.) THEN CKN0=KAPPA*KAPPA*2.**(1.-DIM-VN)*PI**(-DIM/2.)/GAMMA*ACOR**VN ELSE CKN0=KAPPA*KAPPA*2.**(-DIM-2*VN)*PI**(-DIM/2.)* * EXP(GAMMLN(ABS(VN)))/GAMMA ENDIF ENDIF C OPEN(LU2,FILE=FILOUT) C C Loop over points x: DO 23 I3=1,N3 COOR(3)=O3+FLOAT(I3-1)*D3 DO 22 I2=1,N2 COOR(2)=O2+FLOAT(I2-1)*D2 DO 21 I1=1,N1 COOR(1)=O1+FLOAT(I1-1)*D1 IF (DIM2VN.EQ.0.) THEN C Correlation function is the Dirac distribution C realized by several nonzero values around X0: IF (((N1.EQ.1).OR.(ABS(COOR1-X01).LT.ABSD1)).AND. * ((N2.EQ.1).OR.(ABS(COOR2-X02).LT.ABSD2)).AND. * ((N3.EQ.1).OR.(ABS(COOR3-X03).LT.ABSD3))) THEN CKN=1. IF (N1.NE.1) CKN=CKN*(1-ABS(COOR1-X01)/ABSD1)/ABSD1 IF (N2.NE.1) CKN=CKN*(1-ABS(COOR2-X02)/ABSD2)/ABSD2 IF (N3.NE.1) CKN=CKN*(1-ABS(COOR3-X03)/ABSD3)/ABSD3 ELSE CKN=0. ENDIF ELSE XX=SQRT((COOR(1)-X01)**2+ * (COOR(2)-X02)**2+(COOR(3)-X03)**2) IF (XX.EQ.0.) THEN C The value of correlation function is infinite, C realized by CKNMAX at the point X0: CKN=CKNMAX ELSE IF (ACOR.NE.999999.) THEN C Correlation function is computed according to the C equation 48: C Computing the value of the MacDonald function: CALL BESSIK(XX/ACOR,ABS(VN),RI,RK,RIP,RKP) C Computing the correlation function: CKN=CKN0*XX**VN*RK ELSE C Correlation function is computed according to the C equation 66: CKN=CKN0*XX**(2.*VN) ENDIF ENDIF IF (CKN.GT.CKNMAX) CKN=CKNMAX ENDIF WRITE(LU2,*) CKN 21 CONTINUE 22 CONTINUE 23 CONTINUE C WRITE(LU2,*) '/' CLOSE(LU2) WRITE(*,'(A)') '+GRDCKN: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'gammln.for' C gammln.for of Numerical Recipes INCLUDE 'bessik.for' C bessik.for of Numerical Recipes INCLUDE 'beschb.for' C beschb.for of Numerical Recipes INCLUDE 'chebev.for' C chebev.for of Numerical Recipes C C======================================================================= C