SUBROUTINE
LGOBFUN_DV_DV
(n
, x
, y
, wts
, x0
, y0
, pp
, ppd0
, ppd
, hx
, hy
, ll
&
&
, lld
, lldd
, cv
, fixrho
, nbdirs
, nbdirs0
)
USE
DIFFSIZES
IMPLICIT NONE
INTEGER, INTENT(IN)
:: n
REAL*8, DIMENSION(n), INTENT(IN)
:: x
REAL*8, DIMENSION(n), INTENT(IN)
:: y
REAL*8, DIMENSION(n), INTENT(IN)
:: wts
REAL*8, INTENT(IN)
:: x0
REAL*8, INTENT(IN)
:: y0
REAL*8, DIMENSION(5), INTENT(IN)
:: pp
REAL*8, DIMENSION(nbdirsmax, 5), INTENT(IN)
:: ppd0
REAL*8, DIMENSION(5, 5), INTENT(IN)
:: ppd
REAL*8, INTENT(IN)
:: hx
REAL*8, INTENT(IN)
:: hy
LOGICAL, INTENT(IN)
:: cv
REAL*8, INTENT(IN)
:: fixrho
REAL*8, INTENT(OUT)
:: ll
REAL*8, INTENT(OUT)
:: lld
(5
)
REAL*8, INTENT(OUT)
:: lldd
(nbdirsmax
, 5
)
REAL*8, DIMENSION(n)
:: lgauss
REAL*8, DIMENSION(nbdirsmax, n)
:: lgaussd0
REAL*8, DIMENSION(5, n)
:: lgaussd
REAL*8, DIMENSION(nbdirsmax, 5, n)
:: lgaussdd
REAL*8, DIMENSION(5)
:: pars2
REAL*8, DIMENSION(nbdirsmax, 5)
:: pars2d0
REAL*8, DIMENSION(5, 5)
:: pars2d
REAL*8, DIMENSION(nbdirsmax, 5, 5)
:: pars2dd
REAL*8, DIMENSION(1)
:: xtmp
, ytmp
, restmp
REAL*8, DIMENSION(nbdirsmax, 1)
:: xtmpd0
, ytmpd0
, restmpd0
REAL*8, DIMENSION(5, 1)
:: xtmpd
, ytmpd
, restmpd
REAL*8, DIMENSION(nbdirsmax, 5, 1)
:: restmpdd
REAL*8, DIMENSION(5)
:: pars
REAL*8, DIMENSION(nbdirsmax, 5)
:: parsd0
REAL*8, DIMENSION(5, 5)
:: parsd
REAL*8, DIMENSION(nbdirsmax, 5, 5)
:: parsdd
REAL*8, DIMENSION(n)
:: arg1
REAL*8, DIMENSION(5, n)
:: arg1d
REAL*8, DIMENSION(nbdirsmax, 5, n)
:: arg1dd
REAL*8
:: arg10
REAL*8
:: arg10d0
(nbdirsmax
)
REAL*8
:: arg10d
(5
)
REAL*8
:: arg10dd
(nbdirsmax
, 5
)
INTEGER
:: nd
INTEGER
:: nbdirs
INTRINSIC
EXP
INTRINSIC
ABS
INTRINSIC
SUM
REAL*8
:: abs2
REAL*8
:: abs1
INTRINSIC
SQRT
REAL*8
:: result1
REAL*8
:: result1d
(nbdirsmax
)
INTEGER
:: nd0
INTEGER
:: nbdirs0
ll = 0.0_8
IF
(cv) THEN
DO
nd0=1
,nbdirs0
parsdd(nd0, :, :) = 0.0_8
END DO
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
parsdd(nd0, nd, 1
:2
) = 0.0_8
parsdd(nd0, nd, 3
:4
) = ppd(nd, 3
:4
)*ppd0(nd0, 3
:4
)*EXP
(pp(3
:4
))
END DO
parsd(nd, 1
:2
) = ppd(nd, 1
:2
)
parsd(nd, 3
:4
) = ppd(nd, 3
:4
)*EXP
(pp(3
:4
))
END DO
DO
nd0=1
,nbdirs0
parsd0(nd0, 1
:2
) = ppd0(nd0, 1
:2
)
parsd0(nd0, 3
:4
) = ppd0(nd0, 3
:4
)*EXP
(pp(3
:4
))
END DO
pars(1
:2
) = pp(1
:2
)
pars(3
:4
) = EXP
(pp(3
:4
))
IF
(fixrho .GE. 0.
) THEN
abs1 = fixrho
ELSE
abs1 = -fixrho
END IF
IF
(abs1 .LT. 1.0_8
) THEN
DO
nd0=1
,nbdirs0
lldd(nd0, :) = 0.0_8
END DO
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
parsdd(nd0, nd, 5
) = 0.0_8
lldd(nd0, nd) = -(0.5_8
*2
*ppd(nd, 5
)*ppd0(nd0, 5
))
END DO
parsd(nd, 5
) = 0.0_8
lld(nd) = -(0.5_8
*2
*pp(5
)*ppd(nd, 5
))
END DO
DO
nd0=1
,nbdirs0
parsd0(nd0, 5
) = 0.0_8
END DO
pars(5
) = fixrho
ll = -(0.5_8
*pp(5
)**2
)
ELSE
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
parsdd(nd0, nd, 5
) = ((2.0_8
*ppd(nd, 5
)*(ppd0(nd0, 5
)*EXP
(pp(5
&
&
))*(1.0_8
+EXP
(pp(5
)))+EXP
(pp(5
))**2
*ppd0(nd0, 5
))-2.0_8
*ppd(&
&
nd, 5
)*2
*EXP
(pp(5
))**2
*ppd0(nd0, 5
))*(1.0_8
+EXP
(pp(5
)))**2
-(&
&
2.0_8
*ppd(nd, 5
)*EXP
(pp(5
))*(1.0_8
+EXP
(pp(5
)))-2.0_8
*EXP
(pp(&
&
5
))**2
*ppd(nd, 5
))*2
*(1.0_8
+EXP
(pp(5
)))*ppd0(nd0, 5
)*EXP
(pp(&
&
5
)))/((1.0_8
+EXP
(pp(5
)))**2
)**2
END DO
parsd(nd, 5
) = (2.0_8
*ppd(nd, 5
)*EXP
(pp(5
))*(1.0_8
+EXP
(pp(5
)))-&
&
2.0_8
*EXP
(pp(5
))**2
*ppd(nd, 5
))/(1.0_8
+EXP
(pp(5
)))**2
END DO
DO
nd0=1
,nbdirs0
parsd0(nd0, 5
) = (2.0_8
*ppd0(nd0, 5
)*EXP
(pp(5
))*(1.0_8
+EXP
(pp(5
)&
&
))-2.0_8
*EXP
(pp(5
))**2
*ppd0(nd0, 5
))/(1.0_8
+EXP
(pp(5
)))**2
END DO
pars(5
) = -1.0_8
+ 2.0_8
*EXP
(pp(5
))/(1.0_8
+EXP
(pp(5
)))
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
lldd(nd0, nd) = 0.0_8
END DO
lld(nd) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
lldd(nd0, :) = 0.0_8
END DO
END IF
ELSE
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
parsdd(nd0, nd, :) = 0.0_8
END DO
parsd(nd, :) = ppd(nd, :)
END DO
DO
nd0=1
,nbdirs0
parsd0(nd0, :) = ppd0(nd0, :)
END DO
pars = pp
IF
(fixrho .GE. 0.
) THEN
abs2 = fixrho
ELSE
abs2 = -fixrho
END IF
IF
(abs2 .LT. 1.0_8
) THEN
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
parsdd(nd0, nd, 5
) = 0.0_8
END DO
parsd(nd, 5
) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
parsd0(nd0, 5
) = 0.0_8
END DO
pars(5
) = fixrho
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
lldd(nd0, nd) = 0.0_8
END DO
lld(nd) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
lldd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
parsdd(nd0, :, :) = 0.0_8
END DO
ELSE
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
lldd(nd0, nd) = 0.0_8
END DO
lld(nd) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
lldd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
parsdd(nd0, :, :) = 0.0_8
END DO
END IF
END IF
CALL
LOGGAUSSPDF_DV_DV
(n, x, y, pars, parsd0, parsd, parsdd, lgauss, &
&
lgaussd0, lgaussd, lgaussdd, nbdirs, nbdirs0)
DO
nd0=1
,nbdirs0
arg10d0(nd0) = 2
*pars(3
)*parsd0(nd0, 3
)
END DO
arg10 = pars(3
)**2
+ hx**2
DO
nd0=1
,nbdirs0
arg1dd(nd0, :, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
pars2dd(nd0, :, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
arg10dd(nd0, :) = 0.0_8
END DO
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
arg1dd(nd0, nd, :) = wts*lgaussdd(nd0, nd, :)
lldd(nd0, nd) = lldd(nd0, nd) + SUM
(arg1dd(nd0, nd, :))/(1.0_8
*n)
pars2dd(nd0, nd, 1
:2
) = parsdd(nd0, nd, 1
:2
)
arg10dd(nd0, nd) = 2
*(parsd0(nd0, 3
)*parsd(nd, 3
)+pars(3
)*parsdd(&
&
nd0, nd, 3
))
END DO
arg1d(nd, :) = wts*lgaussd(nd, :)
lld(nd) = lld(nd) + SUM
(arg1d(nd, :))/(1.0_8
*n)
pars2d(nd, 1
:2
) = parsd(nd, 1
:2
)
arg10d(nd) = 2
*pars(3
)*parsd(nd, 3
)
IF
(arg10 .EQ. 0.0
) THEN
DO
nd0=1
,nbdirs0
pars2dd(nd0, nd, 3
) = 0.0_8
END DO
pars2d(nd, 3
) = 0.0_8
ELSE
result1 = SQRT
(arg10)
DO
nd0=1
,nbdirs0
IF
(arg10 .EQ. 0.0
) THEN
result1d(nd0) = 0.0_8
ELSE
result1d(nd0) = arg10d0(nd0)/(2.0
*SQRT
(arg10))
END IF
pars2dd(nd0, nd, 3
) = (arg10dd(nd0, nd)*2.0
*result1-arg10d(nd)*&
&
2.0
*result1d(nd0))/(2.0
*result1)**2
END DO
pars2d(nd, 3
) = arg10d(nd)/(2.0
*result1)
END IF
DO
nd0=1
,nbdirs0
arg10dd(nd0, nd) = 2
*(parsd0(nd0, 4
)*parsd(nd, 4
)+pars(4
)*parsdd(&
&
nd0, nd, 4
))
END DO
arg10d(nd) = 2
*pars(4
)*parsd(nd, 4
)
xtmpd(nd, 1
) = 0.0_8
ytmpd(nd, 1
) = 0.0_8
END DO
arg1(:) = wts*lgauss
ll = ll + SUM
(arg1(:))/(1.0_8
*n)
DO
nd0=1
,nbdirs0
pars2d0(nd0, 1
:2
) = parsd0(nd0, 1
:2
)
IF
(arg10 .EQ. 0.0
) THEN
pars2d0(nd0, 3
) = 0.0_8
ELSE
pars2d0(nd0, 3
) = arg10d0(nd0)/(2.0
*SQRT
(arg10))
END IF
arg10d0(nd0) = 2
*pars(4
)*parsd0(nd0, 4
)
END DO
pars2(1
:2
) = pars(1
:2
)
pars2(3
) = SQRT
(arg10)
arg10 = pars(4
)**2
+ hy**2
DO
nd0=1
,nbdirs0
IF
(arg10 .EQ. 0.0
) THEN
pars2d0(nd0, 4
) = 0.0_8
ELSE
pars2d0(nd0, 4
) = arg10d0(nd0)/(2.0
*SQRT
(arg10))
END IF
END DO
pars2(4
) = SQRT
(arg10)
DO
nd=1
,nbdirs
IF
(arg10 .EQ. 0.0
) THEN
DO
nd0=1
,nbdirs0
pars2dd(nd0, nd, 4
) = 0.0_8
END DO
pars2d(nd, 4
) = 0.0_8
ELSE
result1 = SQRT
(arg10)
DO
nd0=1
,nbdirs0
IF
(arg10 .EQ. 0.0
) THEN
result1d(nd0) = 0.0_8
ELSE
result1d(nd0) = arg10d0(nd0)/(2.0
*SQRT
(arg10))
END IF
pars2dd(nd0, nd, 4
) = (arg10dd(nd0, nd)*2.0
*result1-arg10d(nd)*&
&
2.0
*result1d(nd0))/(2.0
*result1)**2
END DO
pars2d(nd, 4
) = arg10d(nd)/(2.0
*result1)
END IF
DO
nd0=1
,nbdirs0
pars2dd(nd0, nd, 5
) = ((((parsdd(nd0, nd, 5
)*pars(3
)+parsd(nd, 5
)*&
&
parsd0(nd0, 3
)+parsd0(nd0, 5
)*parsd(nd, 3
)+pars(5
)*parsdd(nd0, &
&
nd, 3
))*pars(4
)+(parsd(nd, 5
)*pars(3
)+pars(5
)*parsd(nd, 3
))*&
&
parsd0(nd0, 4
)+(parsd0(nd0, 5
)*pars(3
)+pars(5
)*parsd0(nd0, 3
))*&
&
parsd(nd, 4
)+pars(5
)*pars(3
)*parsdd(nd0, nd, 4
))*pars2(3
)*pars2(&
&
4
)+((parsd(nd, 5
)*pars(3
)+pars(5
)*parsd(nd, 3
))*pars(4
)+pars(5
)*&
&
pars(3
)*parsd(nd, 4
))*(pars2d0(nd0, 3
)*pars2(4
)+pars2(3
)*pars2d0&
&
(nd0, 4
))-((parsd0(nd0, 5
)*pars(3
)+pars(5
)*parsd0(nd0, 3
))*pars(&
&
4
)+pars(5
)*pars(3
)*parsd0(nd0, 4
))*(pars2d(nd, 3
)*pars2(4
)+pars2&
&
(3
)*pars2d(nd, 4
))-pars(5
)*pars(3
)*pars(4
)*(pars2dd(nd0, nd, 3
)*&
&
pars2(4
)+pars2d(nd, 3
)*pars2d0(nd0, 4
)+pars2d0(nd0, 3
)*pars2d(nd&
&
, 4
)+pars2(3
)*pars2dd(nd0, nd, 4
)))*pars2(3
)**2
*pars2(4
)**2
-(((&
&
parsd(nd, 5
)*pars(3
)+pars(5
)*parsd(nd, 3
))*pars(4
)+pars(5
)*pars(&
&
3
)*parsd(nd, 4
))*pars2(3
)*pars2(4
)-pars(5
)*pars(3
)*pars(4
)*(&
&
pars2d(nd, 3
)*pars2(4
)+pars2(3
)*pars2d(nd, 4
)))*2
*pars2(3
)*pars2&
&
(4
)*(pars2d0(nd0, 3
)*pars2(4
)+pars2(3
)*pars2d0(nd0, 4
)))/((pars2&
&
(3
)*pars2(4
))**2
)**2
END DO
pars2d(nd, 5
) = (((parsd(nd, 5
)*pars(3
)+pars(5
)*parsd(nd, 3
))*pars(4
&
&
)+pars(5
)*pars(3
)*parsd(nd, 4
))*pars2(3
)*pars2(4
)-pars(5
)*pars(3
)*&
&
pars(4
)*(pars2d(nd, 3
)*pars2(4
)+pars2(3
)*pars2d(nd, 4
)))/(pars2(3
)&
&
*pars2(4
))**2
END DO
DO
nd0=1
,nbdirs0
pars2d0(nd0, 5
) = (((parsd0(nd0, 5
)*pars(3
)+pars(5
)*parsd0(nd0, 3
))*&
&
pars(4
)+pars(5
)*pars(3
)*parsd0(nd0, 4
))*pars2(3
)*pars2(4
)-pars(5
)*&
&
pars(3
)*pars(4
)*(pars2d0(nd0, 3
)*pars2(4
)+pars2(3
)*pars2d0(nd0, 4
)&
&
))/(pars2(3
)*pars2(4
))**2
xtmpd0(nd0, 1
) = 0.0_8
ytmpd0(nd0, 1
) = 0.0_8
END DO
pars2(5
) = pars(5
)*pars(3
)*pars(4
)/(pars2(3
)*pars2(4
))
xtmp(1
) = x0
ytmp(1
) = y0
CALL
LOGGAUSSPDF_DV_DV
(1
, xtmp, ytmp, pars2, pars2d0, pars2d, pars2dd&
&
, restmp, restmpd0, restmpd, restmpdd, nbdirs, &
&
nbdirs0)
DO
nd=1
,nbdirs
DO
nd0=1
,nbdirs0
lldd(nd0, nd) = lldd(nd0, nd) - restmpdd(nd0, nd, 1
)*EXP
(restmp(1
)&
&
) - restmpd(nd, 1
)*restmpd0(nd0, 1
)*EXP
(restmp(1
))
END DO
lld(nd) = lld(nd) - restmpd(nd, 1
)*EXP
(restmp(1
))
END DO
ll = ll - EXP
(restmp(1
))
END
SUBROUTINE
LGOBFUN_DV_DV
SUBROUTINE
LOGGAUSSPDF_DV_DV
(n
, x
, y
, pars
, parsd0
, parsd
, parsdd
, res
, &
&
resd0
, resd
, resdd
, nbdirs
, nbdirs0
)
USE
DIFFSIZES
IMPLICIT NONE
REAL*8, PARAMETER
:: twopi
=6.283185307179586e+00_8
INTEGER, INTENT(IN)
:: n
REAL*8, DIMENSION(n), INTENT(IN)
:: x
REAL*8, DIMENSION(n), INTENT(IN)
:: y
REAL*8, DIMENSION(5), INTENT(IN)
:: pars
REAL*8, DIMENSION(nbdirsmax, 5), INTENT(IN)
:: parsd0
REAL*8, DIMENSION(5, 5), INTENT(IN)
:: parsd
REAL*8, DIMENSION(nbdirsmax, 5, 5), INTENT(IN)
:: parsdd
REAL*8, DIMENSION(n), INTENT(OUT)
:: res
REAL*8, DIMENSION(nbdirsmax, n), INTENT(OUT)
:: resd0
REAL*8, DIMENSION(5, n), INTENT(OUT)
:: resd
REAL*8, DIMENSION(nbdirsmax, 5, n), INTENT(OUT)
:: resdd
REAL*8, DIMENSION(n)
:: cen1
, cen2
REAL*8, DIMENSION(nbdirsmax, n)
:: cen1d0
, cen2d0
REAL*8, DIMENSION(5, n)
:: cen1d
, cen2d
REAL*8, DIMENSION(nbdirsmax, 5, n)
:: cen1dd
, cen2dd
REAL*8
:: t1
, f1
, f2
, f12
REAL*8
:: t1d0
(nbdirsmax
), f1d0
(nbdirsmax
), f2d0
(nbdirsmax
), f12d0
(&
&
nbdirsmax
)
REAL*8
:: t1d
(5
), f1d
(5
), f2d
(5
), f12d
(5
)
REAL*8
:: t1dd
(nbdirsmax
, 5
), f1dd
(nbdirsmax
, 5
), f2dd
(nbdirsmax
, 5
), &
&
f12dd
(nbdirsmax
, 5
)
REAL*8
:: arg1
REAL*8
:: arg1d0
(nbdirsmax
)
REAL*8
:: arg1d
(5
)
REAL*8
:: arg1dd
(nbdirsmax
, 5
)
REAL*8
:: result1
REAL*8
:: result1d0
(nbdirsmax
)
REAL*8
:: result1d
(5
)
REAL*8
:: result1dd
(nbdirsmax
, 5
)
REAL*8
:: arg2
REAL*8
:: arg2d0
(nbdirsmax
)
REAL*8
:: arg2d
(5
)
REAL*8
:: arg2dd
(nbdirsmax
, 5
)
INTEGER
:: nd
INTEGER
:: nbdirs
INTRINSIC
LOG
INTRINSIC
SQRT
REAL*8
:: result10
REAL*8
:: result10d
(nbdirsmax
)
INTEGER
:: nd0
INTEGER
:: nbdirs0
t1 = -(0.5_8
/(1.0_8
-pars(5
)**2
))
arg1 = 1.0_8
- pars(5
)**2
result1 = SQRT
(arg1)
DO
nd0=1
,nbdirs0
t1d0(nd0) = (-(0.5_8
*2
*pars(5
)*parsd0(nd0, 5
)))/(1.0_8
-pars(5
)**2
)**&
&
2
f1d0(nd0) = (t1d0(nd0)*pars(3
)**2
-t1*2
*pars(3
)*parsd0(nd0, 3
))/(pars&
&
(3
)**2
)**2
f2d0(nd0) = (t1d0(nd0)*pars(4
)**2
-t1*2
*pars(4
)*parsd0(nd0, 4
))/(pars&
&
(4
)**2
)**2
f12d0(nd0) = -((2.0
*(parsd0(nd0, 5
)*t1+pars(5
)*t1d0(nd0))*pars(3
)*&
&
pars(4
)-2.0
*pars(5
)*t1*(parsd0(nd0, 3
)*pars(4
)+pars(3
)*parsd0(nd0&
&
, 4
)))/(pars(3
)*pars(4
))**2
)
cen1d0(nd0, :) = -parsd0(nd0, 1
)
cen2d0(nd0, :) = -parsd0(nd0, 2
)
arg1d0(nd0) = -(2
*pars(5
)*parsd0(nd0, 5
))
IF
(arg1 .EQ. 0.0
) THEN
result1d0(nd0) = 0.0_8
ELSE
result1d0(nd0) = arg1d0(nd0)/(2.0
*SQRT
(arg1))
END IF
arg2d0(nd0) = twopi*((parsd0(nd0, 3
)*result1+pars(3
)*result1d0(nd0))&
&
*pars(4
)+pars(3
)*result1*parsd0(nd0, 4
))
END DO
f1 = t1/pars(3
)**2
f2 = t1/pars(4
)**2
f12 = -(2.0
*pars(5
)*t1/(pars(3
)*pars(4
)))
cen1 = x - pars(1
)
cen2 = y - pars(2
)
arg2 = twopi*pars(3
)*pars(4
)*result1
DO
nd0=1
,nbdirs0
resdd(nd0, :, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
arg1dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
f12dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
cen2dd(nd0, :, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
f1dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
arg2dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
f2dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
result1dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
t1dd(nd0, :) = 0.0_8
END DO
DO
nd0=1
,nbdirs0
cen1dd(nd0, :, :) = 0.0_8
END DO
DO
nd=1
,nbdirs
t1d(nd) = (-(0.5_8
*2
*pars(5
)*parsd(nd, 5
)))/(1.0_8
-pars(5
)**2
)**2
DO
nd0=1
,nbdirs0
t1dd(nd0, nd) = (-(0.5_8
*2
*(parsd0(nd0, 5
)*parsd(nd, 5
)+pars(5
)*&
&
parsdd(nd0, nd, 5
))*(1.0_8
-pars(5
)**2
)**2
)-0.5_8
*2
**3
*pars(5
)**2
&
&
*parsd(nd, 5
)*(1.0_8
-pars(5
)**2
)*parsd0(nd0, 5
))/((1.0_8
-pars(5
)&
&
**2
)**2
)**2
f1dd(nd0, nd) = ((t1dd(nd0, nd)*pars(3
)**2
+t1d(nd)*2
*pars(3
)*&
&
parsd0(nd0, 3
)-2
*((t1d0(nd0)*pars(3
)+t1*parsd0(nd0, 3
))*parsd(nd&
&
, 3
))-2
*(t1*pars(3
)*parsdd(nd0, nd, 3
)))*pars(3
)**4
-(t1d(nd)*&
&
pars(3
)**2
-t1*2
*pars(3
)*parsd(nd, 3
))*2
**2
*pars(3
)**3
*parsd0(nd0&
&
, 3
))/((pars(3
)**2
)**2
)**2
f2dd(nd0, nd) = ((t1dd(nd0, nd)*pars(4
)**2
+t1d(nd)*2
*pars(4
)*&
&
parsd0(nd0, 4
)-2
*((t1d0(nd0)*pars(4
)+t1*parsd0(nd0, 4
))*parsd(nd&
&
, 4
))-2
*(t1*pars(4
)*parsdd(nd0, nd, 4
)))*pars(4
)**4
-(t1d(nd)*&
&
pars(4
)**2
-t1*2
*pars(4
)*parsd(nd, 4
))*2
**2
*pars(4
)**3
*parsd0(nd0&
&
, 4
))/((pars(4
)**2
)**2
)**2
f12dd(nd0, nd) = -(((2.0
*((parsdd(nd0, nd, 5
)*t1+parsd(nd, 5
)*t1d0&
&
(nd0)+parsd0(nd0, 5
)*t1d(nd)+pars(5
)*t1dd(nd0, nd))*pars(3
)*pars&
&
(4
)+(parsd(nd, 5
)*t1+pars(5
)*t1d(nd))*(parsd0(nd0, 3
)*pars(4
)+&
&
pars(3
)*parsd0(nd0, 4
)))-2.0
*((parsd0(nd0, 5
)*t1+pars(5
)*t1d0(&
&
nd0))*(parsd(nd, 3
)*pars(4
)+pars(3
)*parsd(nd, 4
))+pars(5
)*t1*(&
&
parsdd(nd0, nd, 3
)*pars(4
)+parsd(nd, 3
)*parsd0(nd0, 4
)+parsd0(&
&
nd0, 3
)*parsd(nd, 4
)+pars(3
)*parsdd(nd0, nd, 4
))))*pars(3
)**2
*&
&
pars(4
)**2
-(2.0
*(parsd(nd, 5
)*t1+pars(5
)*t1d(nd))*pars(3
)*pars(4
&
&
)-2.0
*pars(5
)*t1*(parsd(nd, 3
)*pars(4
)+pars(3
)*parsd(nd, 4
)))*2
*&
&
pars(3
)*pars(4
)*(parsd0(nd0, 3
)*pars(4
)+pars(3
)*parsd0(nd0, 4
)))&
&
/((pars(3
)*pars(4
))**2
)**2
)
cen1dd(nd0, nd, :) = -parsdd(nd0, nd, 1
)
cen2dd(nd0, nd, :) = -parsdd(nd0, nd, 2
)
arg1dd(nd0, nd) = -(2
*(parsd0(nd0, 5
)*parsd(nd, 5
)+pars(5
)*parsdd(&
&
nd0, nd, 5
)))
END DO
f1d(nd) = (t1d(nd)*pars(3
)**2
-t1*2
*pars(3
)*parsd(nd, 3
))/(pars(3
)**2
&
&
)**2
f2d(nd) = (t1d(nd)*pars(4
)**2
-t1*2
*pars(4
)*parsd(nd, 4
))/(pars(4
)**2
&
&
)**2
f12d(nd) = -((2.0
*(parsd(nd, 5
)*t1+pars(5
)*t1d(nd))*pars(3
)*pars(4
)-&
&
2.0
*pars(5
)*t1*(parsd(nd, 3
)*pars(4
)+pars(3
)*parsd(nd, 4
)))/(pars(&
&
3
)*pars(4
))**2
)
cen1d(nd, :) = -parsd(nd, 1
)
cen2d(nd, :) = -parsd(nd, 2
)
arg1d(nd) = -(2
*pars(5
)*parsd(nd, 5
))
IF
(arg1 .EQ. 0.0
) THEN
DO
nd0=1
,nbdirs0
result1dd(nd0, nd) = 0.0_8
END DO
result1d(nd) = 0.0
ELSE
result10 = SQRT
(arg1)
DO
nd0=1
,nbdirs0
IF
(arg1 .EQ. 0.0
) THEN
result10d(nd0) = 0.0_8
ELSE
result10d(nd0) = arg1d0(nd0)/(2.0
*SQRT
(arg1))
END IF
result1dd(nd0, nd) = (arg1dd(nd0, nd)*2.0
*result10-arg1d(nd)*2.0
&
&
*result10d(nd0))/(2.0
*result10)**2
END DO
result1d(nd) = arg1d(nd)/(2.0
*result10)
END IF
arg2d(nd) = twopi*((parsd(nd, 3
)*result1+pars(3
)*result1d(nd))*pars(&
&
4
)+pars(3
)*result1*parsd(nd, 4
))
DO
nd0=1
,nbdirs0
arg2dd(nd0, nd) = twopi*((parsdd(nd0, nd, 3
)*result1+parsd(nd, 3
)*&
&
result1d0(nd0)+parsd0(nd0, 3
)*result1d(nd)+pars(3
)*result1dd(nd0&
&
, nd))*pars(4
)+(parsd(nd, 3
)*result1+pars(3
)*result1d(nd))*&
&
parsd0(nd0, 4
)+(parsd0(nd0, 3
)*result1+pars(3
)*result1d0(nd0))*&
&
parsd(nd, 4
)+pars(3
)*result1*parsdd(nd0, nd, 4
))
resdd(nd0, nd, :) = f1dd(nd0, nd)*cen1**2
+ f1d(nd)*2
*cen1*cen1d0(&
&
nd0, :) - (arg2dd(nd0, nd)*arg2-arg2d(nd)*arg2d0(nd0))/arg2**2
+&
&
2
*((f1d0(nd0)*cen1+f1*cen1d0(nd0, :))*cen1d(nd, :)) + 2
*(f1*cen1&
&
*cen1dd(nd0, nd, :)) + f2dd(nd0, nd)*cen2**2
+ f2d(nd)*2
*cen2*&
&
cen2d0(nd0, :) + 2
*((f2d0(nd0)*cen2+f2*cen2d0(nd0, :))*cen2d(nd&
&
, :)) + 2
*(f2*cen2*cen2dd(nd0, nd, :)) + (f12dd(nd0, nd)*cen1+&
&
f12d(nd)*cen1d0(nd0, :)+f12d0(nd0)*cen1d(nd, :)+f12*cen1dd(nd0, &
&
nd, :))*cen2 + (f12d(nd)*cen1+f12*cen1d(nd, :))*cen2d0(nd0, :) +&
&
(f12d0(nd0)*cen1+f12*cen1d0(nd0, :))*cen2d(nd, :) + f12*cen1*&
&
cen2dd(nd0, nd, :)
END DO
resd(nd, :) = f1d(nd)*cen1**2
- arg2d(nd)/arg2 + f1*2
*cen1*cen1d(nd&
&
, :) + f2d(nd)*cen2**2
+ f2*2
*cen2*cen2d(nd, :) + (f12d(nd)*cen1+&
&
f12*cen1d(nd, :))*cen2 + f12*cen1*cen2d(nd, :)
END DO
DO
nd0=1
,nbdirs0
resd0(nd0, :) = f1d0(nd0)*cen1**2
- arg2d0(nd0)/arg2 + f1*2
*cen1*&
&
cen1d0(nd0, :) + f2d0(nd0)*cen2**2
+ f2*2
*cen2*cen2d0(nd0, :) + (&
&
f12d0(nd0)*cen1+f12*cen1d0(nd0, :))*cen2 + f12*cen1*cen2d0(nd0, :)
END DO
res = -LOG
(arg2) + f1*cen1**2
+ f2*cen2**2
+ f12*cen1*cen2
END
SUBROUTINE
LOGGAUSSPDF_DV_DV