125 lines
3.0 KiB
FortranFixed
125 lines
3.0 KiB
FortranFixed
!-----------4 parameters---------------------------------
|
|
! Non-rectangular hyperbola function and derivatives
|
|
! y=a+b(x-sqrt(c+dx+x2)) or y=(ax+b-sqrt((ax+b)^2-4abcx))/2c-d
|
|
!
|
|
subroutine fnonrecthypb(a,b,c,d,x,y,iwrong)
|
|
implicit none
|
|
integer iwrong
|
|
double precision a,b,c,d,x,y,p
|
|
iwrong=0
|
|
goto 10
|
|
!----------------------------------------
|
|
!a=alpha, b=Amax, c=theta, d=rd
|
|
p=(a*x+b)*(a*x+b)-4.0d0*a*b*c*x
|
|
if(p.lt.0.0d0)then
|
|
iwrong=1
|
|
return
|
|
endif
|
|
y=(a*x+b-dsqrt(p))/(2.0d0*c)-d
|
|
return
|
|
!----------------------------------------
|
|
10 p=c+d*x+x*x
|
|
if(p.lt.0.0d0)then
|
|
iwrong=1
|
|
return
|
|
endif
|
|
y=a+b*(x-dsqrt(p))
|
|
return
|
|
end
|
|
|
|
subroutine indices_fnonrecthypb(a,b,c,d,root,
|
|
& der_root,fmax,iwrong)
|
|
implicit none
|
|
double precision a,b,c,d,root,der_root,fmax,p
|
|
integer iwrong
|
|
iwrong=0
|
|
goto 10
|
|
!---------------------------
|
|
root=(c*d*d-b*d)/(a*d-a*b)
|
|
fmax=b
|
|
p=(a*root+b)*(a*root+b)-4.0d0*a*b*c*root
|
|
if(p.lt.0.0d0)then
|
|
iwrong=1
|
|
return
|
|
endif
|
|
der_root=(a-a*(a*root+b-2.0d0*b*c)/dsqrt(p))
|
|
& /(2.0d0*c)
|
|
return
|
|
!------------------------------
|
|
10 root=(b*b*c-a*a)/(2.0d0*a*b-b*b*d)
|
|
fmax=a-0.5d0*b*d
|
|
p=root*root+d*root+c
|
|
if(p.lt.0.0d0)then
|
|
iwrong=1
|
|
return
|
|
endif
|
|
der_root=b*(1.0d0-(2.0d0*root+d)/
|
|
& (2.0d0*dsqrt(p)))
|
|
return
|
|
end
|
|
|
|
subroutine der_fnonrecthypb(a,b,c,d,x,da,db,dc,dd,dx,
|
|
& iwrong)
|
|
implicit none
|
|
integer iwrong
|
|
double precision a,b,c,d,x,da,db,dc,dd,dx
|
|
double precision p
|
|
iwrong=0
|
|
goto 10
|
|
!-----------------------------------------
|
|
p=(a*x+b)*(a*x+b)-4.0d0*a*b*c*x
|
|
if(p.lt.0.0d0)then
|
|
iwrong=1
|
|
return
|
|
endif
|
|
p=dsqrt(p)
|
|
da=(x-x*(a*x+b-2.0d0*b*c)/p)/(2.0d0*c)
|
|
db=(1.0d0-(a*x+b-2.0d0*a*c*x)/p)/(2.0d0*c)
|
|
dd=-1.0d0
|
|
dc=a*b*x/(c*p)-(a*x+b-p)/(2.0d0*c*c)
|
|
dx=(a-a*(a*x+b-2.0d0*b*c)/p)/(2.0d0*c)
|
|
return
|
|
!------------------------------------------
|
|
10 p=c+d*x+x*x
|
|
if(p.lt.0.0d0)then
|
|
iwrong=1
|
|
return
|
|
endif
|
|
p=dsqrt(p)
|
|
da=1.0d0
|
|
db=x-p
|
|
dc=-b/(2.0d0*p)
|
|
dd=dc*x
|
|
dx=b*(1.0d0-(d+2.0d0*x)/(2.0d0*p))
|
|
return
|
|
end
|
|
|
|
!-------3 parameters----------------------
|
|
subroutine recthypb(a,b,c,x,y)
|
|
implicit none
|
|
double precision a,b,c,x,y
|
|
y=(a*x+b)/(x+c)
|
|
return
|
|
end
|
|
|
|
subroutine indices_frecthypb(a,b,c,root,
|
|
& der_root,fmax)
|
|
implicit none
|
|
double precision a,b,c,root,
|
|
& der_root,fmax
|
|
root=-b/a
|
|
der_root=a*a/(a*c-b)
|
|
fmax=a
|
|
return
|
|
end
|
|
|
|
subroutine der_recthypb(a,b,c,x,da,db,dc,dx)
|
|
implicit none
|
|
double precision a,b,c,x,da,db,dc,dx
|
|
da=x/(x+c)
|
|
db=1.0d0/(x+c)
|
|
dc=-(a*x+b)/((x+c)*(x+c))
|
|
dx=a/(x+c)-(a*x+b)/((x+c)*(x+c))
|
|
return
|
|
end
|