86 lines
2.1 KiB
Fortran
86 lines
2.1 KiB
Fortran
|
|
SUBROUTINE postmedcauchy(m,x,weight,maxiter,muhat)
|
|
implicit none
|
|
|
|
integer :: m,maxiter
|
|
real*8 :: x(m),weight,muhat(m)
|
|
|
|
real*8 :: magdata(m),magdatatmp(m),weight2(m),hi(m),Tol,conTol,lo(m),midpoint(m),y(m),pi,Xvec(m),normconstants,fx(m),delta(maxiter)
|
|
real*8 :: z(m),y2(m),y3(m),fx2(m),fmidpoint(m),zeromd(m),mu,sigma,yr(m),FUNC_sign
|
|
integer :: i,idx(m),numiter
|
|
|
|
muhat = 0
|
|
magdata = abs(x)
|
|
magdatatmp = magdata
|
|
idx = 1
|
|
delta = 0
|
|
|
|
do i=1,m
|
|
if(magdata(i)>=20) then
|
|
idx(i)=0
|
|
magdata(i) = sqrt(-1.0)
|
|
end if
|
|
end do
|
|
|
|
weight2 = weight
|
|
lo = 0
|
|
hi = maxval(magdata)
|
|
Tol = 1e-9
|
|
numiter = 0
|
|
conTol = 1e30
|
|
mu = 0
|
|
sigma = 1
|
|
pi = 3.1415926
|
|
zeromd = 0
|
|
|
|
do while(conTol>Tol)
|
|
numiter = numiter +1
|
|
midpoint = (lo+hi)/2.0
|
|
|
|
y = magdata-midpoint
|
|
Xvec = abs((y-mu)**2)
|
|
normconstants = 1.0/(sqrt(2.0*pi)*sigma)
|
|
fx = normconstants*exp(-Xvec/(2.0*sigma))
|
|
|
|
z = (y-mu)/sigma
|
|
y2 = 1.0/2.0*erfc(-z/sqrt(2.0))
|
|
|
|
z = (-midpoint-mu)/sigma
|
|
y3 = 1.0/2.0*erfc(-z/sqrt(2.0))
|
|
|
|
Xvec = abs((midpoint-mu)**2)
|
|
fx2 = normconstants*exp(-Xvec/(2.0*sigma))
|
|
|
|
yr = y2 - magdata*fx + ((magdata*midpoint-1)*fx*y3/fx2)
|
|
|
|
y3 = magdata**2.0*(1.0/weight2-1.0)-1.0
|
|
|
|
y2 = 1+exp(-magdata**2/2.0)*y3
|
|
|
|
fmidpoint = y2/2.0-yr
|
|
|
|
do i=1,m
|
|
if(fmidpoint(i) <= zeromd(i)) then
|
|
lo(i) = midpoint(i)
|
|
else
|
|
hi(i) = midpoint(i)
|
|
end if
|
|
end do
|
|
delta(numiter) = maxval(abs(hi-lo))
|
|
conTol = delta(numiter)
|
|
if(numiter > maxiter) exit
|
|
end do
|
|
|
|
muhat = (lo+hi)/2.0
|
|
|
|
do i=1,m
|
|
if(idx(i)==0) then
|
|
muhat(i) = magdatatmp(i)-2.0/magdatatmp(i)
|
|
end if
|
|
if(muhat(i)<1e-7) muhat(i)=0
|
|
muhat(i) = FUNC_sign(x(i))*muhat(i)
|
|
if(abs(muhat(i)) > abs(x(i))) muhat(i) = x(i)
|
|
end do
|
|
|
|
end subroutine
|
|
|