! ******************************************************************* ! * * ! * Function dlafs3 * ! * * ! ******************************************************************* ! Double Precision Version 1.4 ! Written by Gordon A. Fenton, TUNS, Mar. 15, 1995 ! Latest Update: Jan 11, 2005 ! ! PURPOSE returns the covariance between two points in a partially ! separable 3-D fractional Gaussian noise field. ! ! Returns the covariance between points in a 3-D fractional Gaussian noise ! process which is isotropic in the X-Y plane. The covariance function is ! partially separable, and is given by ! ! B(X,Y,Z) = var*r(X,Y)*r(Z) ! ! where r(Z) is the correlation function in the z-direction, ! ! 1 2G 2G 2G ! r(Z) = -------- [ |Z + pb| - 2|Z| + |Z - pb| ] (1) ! 2G ! 2*pb ! ! while, in the X-Y plane, the isotropic correlation function r(X,Y) is ! given by ! ! 1 2H 2H 2H ! r(X,Y) = r(s) = -------- [ |s + pb| - 2|s| + |s - pb| ] (2) ! 2H ! 2*pb ! ! where s = sqrt(X*X + Y*Y), and `pb' is the length over which the fractional ! Brownian motion is averaged in order to make this, the derivative process, ! exist. Normally `pb' is selected to be quite small (of the order of the ! size of the discretization interval). The parameters `G' and `H' are the ! self-similar parameters, or Hurst exponents, (0.5 < G,H < 1) with ! G = H = 1 giving a total correlated field and G = H = 0.5 giving a white ! noise field (for arbitrarily small pb). ! ! The parameters var, pb, H, and G are brought in through the common block ! /dparam/. X, Y, and Z are the function arguments and are the elements ! of the lag vector separating the two points in the field for which the ! covariance is desired. ! ! If var < 0, then this function computes the variance of a local average ! of the random field, averaged over a domain having side dimensions {X,Y,Z}. ! In this case, the variance function (which is also separable) ! ! V(X,Y,Z) = V(X,Y)*V(Z) ! ! is used, and the function returns |var|*V(X,Y,Z). The 1-D ! variance function corresponding to Eq. (1) is given by ! ! r r r r ! |Z + pb| - 2( |Z| + |pb| ) + |Z - pb| ! V(Z) = ----------------------------------------- (3) ! 2 2H ! Z * (2H+1)*(2H+2)*pb ! ! where r = (2*G+2). Notice that as Z goes to zero the above becomes 0/0, but ! its limiting value is 1.0. ! ! The variance function which yields a close approximation to Eq. (2)c is given by ! ! r r r r ! |D + pb| - 2( |D| + |pb| ) + |D - pb| ! V(X,Y) = -------------------------------------------------- (4) ! 2 2G ! D * (2G+1)*(2G+2)*pb ! ! where r = (2H+2) and D = |X| + |Y| (the 1-norm of {X,Y}). The actual ! covariance function corresponding to (4) is equal to (2) along the ! X,Y coordinate axes, and remains close to (2) for G > 0.8. For G < 0.7, ! the covariance in the diagonal direction is too high, but since the ! covariance in such fields dies off very rapidly, the error is ! negligible. ! ! Parameters to this process are brought in through the common block ! /dparam/ and are described as follows; ! ! var the point variance of the fGn process. ! ! pb the averaging length discussed above. ! ! H the Hurst exponent governing the process in the X,Y plane. ! In this implementation, 0.5 < H < 1; values of H near 1 yield ! a correlation structure which remains very high (and thus a ! covariance matrix which may be nearly singular). Values of H ! near 0.5 yield a band- limited white noise process. ! ! G the Hurst exponent governing the process in the Z direction. ! In this implementation, 0.5 < G < 1; values of G near 1 yield ! a correlation structure which remains very high (and thus a ! covariance matrix which may be nearly singular). Values of G ! near 0.5 yield a band- limited white noise process. ! ! da dummy placeholder provided so that the common block /dparam/ ! has the same form for a variety of variance functions. ! ! Arguments to this function are just the X, Y, and Z dimensions of the ! lag distance (or the physical averaging region). ! ! REVISION HISTORY: ! 1.1 added lvarfn flag to end of /dparam/ common block. Set it to ! true if the variance function is to be returned. Otherwise this ! function returns the covariance. (Apr 29/00) ! 1.11 minor changes to documentation above (Jul 14/00) ! 1.2 eliminated lvarfn -- now return covariances only if var < 0 (Apr 9/01) ! 1.3 reversed default - now return covariances if var > 0 (Apr 11/01) ! 1.31 revised above documentation to reflect revision 1.3 (May 9/01) ! 1.4 made X-Y the isotropic plane (Z is usually vertical) (Jan 11/05) !---------------------------------------------------------------------- real*8 function dlafs3( X, Y, Z ) implicit real*8 (a-h,o-z) common/dparam/ var, pb, H, G, da data zero/0.d0/, one/1.d0/, two/2.d0/ abs(q) = dabs(q) sqrt(q) = dsqrt(q) if( var .lt. zero ) then ! return variance function ! Z-direction variance function if( Z .eq. zero ) then v1 = one else e0 = two*G e1 = e0 + one e2 = e0 + two t1 = abs(Z+pb)**e2 - two*(abs(Z)**e2+pb**e2) + abs(Z-pb)**e2 v1 = t1/(Z*Z*e1*e2*(pb**e0)) endif ! X-Y plane variance function if( (X .eq. zero) .and. (Y .eq. zero) ) then v2 = one else D = abs(Y) + abs(X) e0 = two*H e1 = e0 + one e2 = e0 + two t1 = (D+pb)**e2 - two*(D**e2 + pb**e2) + abs(D-pb)**e2 v2 = t1/(D*D*e1*e2*pb**e0) endif dlafs3 = -var*v1*v2 else ! return covariance ! Z-direction covariance if( Z .eq. zero ) then v1 = one else e0 = two*G t1 = abs(Z+pb)**e0 - two*(abs(Z)**e0) + abs(Z-pb)**e0 v1 = t1/(two*(pb**e0)) endif ! Y-direction covariance if( (Y .eq. zero) .and. (X .eq. zero) ) then v2 = one else e0 = two*H D = sqrt(Y*Y + X*X) t1 = (D+pb)**e0 - two*(D**e0) + abs(D-pb)**e0 v2 = t1/(two*(pb**e0)) endif dlafs3 = var*v1*v2 endif return end