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