SUBROUTINE
LOGGAUSSPDF_DV_DV
(n
, x
, y
, pars
, parsd0
, parsd
, parsdd
, res
, &
&
resd0
, resd
, resdd
, nbdirs
, nbdirs0
)
USE
DIFFSIZES_DV
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(nbdirsmax0, 5
), INTENT(IN)
:: parsd0
REAL*8
, DIMENSION(
nbdirsmax, 5
)
, INTENT(IN)
:: parsd
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax, 5
), INTENT(IN)
:: parsdd
REAL*8
, DIMENSION(
n)
, INTENT(OUT)
:: res
REAL*8
, DIMENSION(nbdirsmax0, n), INTENT(OUT)
:: resd0
REAL*8
, DIMENSION(
nbdirsmax, n)
, INTENT(OUT)
:: resd
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax, n), INTENT(OUT)
:: resdd
REAL*8
, DIMENSION(
n)
:: cen1
, cen2
REAL*8
, DIMENSION(nbdirsmax0, n) :: cen1d0
, cen2d0
REAL*8
, DIMENSION(
nbdirsmax, n)
:: cen1d
, cen2d
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax, n) :: cen1dd
, cen2dd
REAL*8
:: t1
, f1
, f2
, f12
REAL*8
, DIMENSION(nbdirsmax0) :: t1d0
, f1d0
, f2d0
, f12d0
REAL*8
, DIMENSION(
nbdirsmax)
:: t1d
, f1d
, f2d
, f12d
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax) :: t1dd
, f1dd
, f2dd
, f12dd
REAL*8
:: arg1
REAL*8
, DIMENSION(nbdirsmax0) :: arg1d0
REAL*8
, DIMENSION(
nbdirsmax)
:: arg1d
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax) :: arg1dd
REAL*8
:: result1
REAL*8
, DIMENSION(nbdirsmax0) :: result1d0
REAL*8
, DIMENSION(
nbdirsmax)
:: result1d
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax) :: result1dd
REAL*8
:: arg2
REAL*8
, DIMENSION(nbdirsmax0) :: arg2d0
REAL*8
, DIMENSION(
nbdirsmax)
:: arg2d
REAL*8
, DIMENSION(nbdirsmax0, nbdirsmax) :: arg2dd
INTEGER
:: nd
INTEGER
:: nbdirs
INTRINSIC
LOG
INTRINSIC
SQRT
REAL*8
:: result10
REAL*8
, DIMENSION(nbdirsmax0) :: result10d
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