C
C Program GRDCKN to compute the values of the Von Karman correlation C function according to the formula (K.4) of the paper Klimes, L.: C Correlation functions of random media. In: Seismic waves in 3-D C structures, Report 6, Department of Geophysics, Charles C University, Prague (1997). C C Version: 5.40 C Date: 2000, February 21 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: bulant@seis.karlov.mff.cuni.cz 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: C ACORG .LT. 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 Default: ACOR=999999. (infinity) C CKNMAX=real ... Maximum value of the correlation function. C Default: CKNMAX=999999. (infinity) 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=0. C D2=real... Grid interval along the X2 axis. C Default: D2=0. C D3=real... Grid interval along the X3 axis. C Default: D3=0. 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,RMAT,WMAT,GAMMLN,BESSIK REAL GAMMLN C ERROR ... File error.for. C RSEP1,RSEP3T,RSEP3R,RSEP3I ... C File sep.for. C RMAT,WMAT ... File forms.for. C GAMMLN ... File gammln.for. C BESSIK ... File bessik.for. 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 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.) 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 C Computing the value of the Gamma function: GAMMA=EXP(GAMMLN(DIM/2.+VN)) C Computing the x-independent part of the correlation function: CKN0=KAPPA*KAPPA*2.**(1.-DIM-VN)*PI**(-DIM/2.)/GAMMA*ACOR**VN CKN=0. C OPEN(LU2,FILE=FILOUT) 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 XX=SQRT((COOR(1)-X01)**2+(COOR(2)-X02)**2+(COOR(3)-X03)**2) IF (XX.NE.0.) THEN 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 CKN=CKNMAX ENDIF IF (CKN.GT.CKNMAX) CKN=CKNMAX WRITE(LU2,*) CKN 21 CONTINUE 22 CONTINUE 23 CONTINUE 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 INCLUDE 'bessik.for' C bessik.for INCLUDE 'beschb.for' C beschb.for INCLUDE 'chebev.for' C chebev.for C C======================================================================= C