Files
piscal/dataassim/math/othersupmath/AsymGaussians.f
2016-02-03 18:52:05 +00:00

115 lines
3.1 KiB
FortranFixed

!This file contains gaussian functions even though the subroutine names are sigmoid
double precision function sigmoidfunc(y0,a,b,c,x0,x)
implicit none
!y=y0+a*exp(-((x-x0)/b)**c)) for x > x0 or y=y0+a*exp(-((x0-x)/b)**c)) for x < x0
double precision y0,a,b,c,x0,x
if(x.gt.x0)then
sigmoidfunc=y0+a*dexp(-((x-x0)/b)**c)
else
sigmoidfunc=y0+a*dexp(-((x0-x)/b)**c)
endif
return
end
!-------------------------------------------------------------------
subroutine gradsigmoidfunc(y0,a,b,c,x0,x,grad)
implicit none
double precision y0,a,b,c,x0,x,grad(6),term,term1
! a<->grad(1)
! b<->grad(2)
! c<->grad(3)
! x0<->grad(4)
! y0<->grad(5)
! x<->grad(6)
grad(5)=1.0d0
if(x.gt.x0)then
term=dexp(-((x-x0)/b)**c)
term1=-((x-x0)/b)**(c-1.0d0)
grad(1)=term
grad(6)=term*c*term1/b
grad(4)=-term*c*term1/b
grad(2)=-term*c*term1*(x-x0)/(b*b)
grad(3)=term*(-((x-x0)/b)**c)*dlog((x-x0)/b)
else
term=dexp(-((x0-x)/b)**c)
term1=-((x0-x)/b)**(c-1.0d0)
grad(1)=term
grad(6)=-term*c*term1/b
grad(4)=term*c*term1/b
grad(2)=-term*c*term1*(x0-x)/(b*b)
grad(3)=term*(-((x0-x)/b)**c)*dlog((x0-x)/b)
endif
return
end
!--------------------------------------------------------------------
double precision function twoexpfunc(y0,a1,b1,c1,x01,
&a2,b2,c2,x02,x)
implicit none
double precision y0,a1,b1,c1,x01,a2,b2,c2,x02,x,sigmoidfunc,
&x0,a,b,c
!In Asymmetrical Gaussians, c1 and c2 have no use.
!y=y0+a1*exp(-((x-x01)/b1)**a2)) for x > x01 or y=y0+a1*exp(-((x01-x)/b2)**x02)) for x < x01
x0=x01
a=a1
if(x.gt.x01)then
b=b1
c=a2
else
b=b2
c=x02
endif
twoexpfunc=sigmoidfunc(y0,a,b,c,x0,x)
return
end
!---------------------------------------------------------------------
subroutine gradtwoexp(y0,a1,b1,c1,x01,a2,b2,c2,x02,x,grad)
implicit none
double precision y0,a1,b1,c1,x01,a2,b2,c2,x02,x,grad(10),grad6(6),
&x0,a,b,c
integer i
! a1<->grad(1)
! b1<->grad(2)
! c1<->grad(3)
! x01<->grad(4)
! y0<->grad(5)
! x<->grad(6)
! a2<->grad(7)
! b2<->grad(8)
! c2<->grad(9)
! x02<->grad(10)
do i=1,10
grad(i)=0.0d0
enddo
x0=x01
a=a1
if(x.gt.x01)then
b=b1
c=a2
call gradsigmoidfunc(y0,a,b,c,x0,x,grad6)
! a<->grad6(1)
! b<->grad6(2)
! c<->grad6(3)
! x0<->grad6(4)
! y0<->grad6(5)
! x<->grad6(6)
grad(1)=grad6(1)
grad(2)=grad6(2)
grad(4)=grad6(4)
grad(5)=grad6(5)
grad(6)=grad6(6)
grad(7)=grad6(3)
else
b=b2
c=x02
call gradsigmoidfunc(y0,a,b,c,x0,x,grad6)
grad(1)=grad6(1)
grad(4)=grad6(4)
grad(5)=grad6(5)
grad(6)=grad6(6)
grad(8)=grad6(2)
grad(10)=grad6(3)
endif
return
end
!------------------------------------------------------------------------------