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