Hi! everybody. In the last topic, I asked for help to find the roots of a nonlinear function.
Today, I get into trouble with interpolation. I don’t know why it doesn’t work.
I use subroutines spline and splint from Fortran 95 recipe. but the compiler checking said that it has errors: "type mismatch in Fortran
Here are my code:
program test_28_08
implicit none
integer, parameter :: nc=1001
integer, parameter :: ni=1000
double precision xmin, xmax
double precision xi(ni), yi(ni)
double precision y2
double precision step, deriv3
double precision yp1,ypn,r,y
integer i, step1
common A,Z,R0,am
real A,Z,R0,am
double precision:: Vn,Vc,Vl
double precision rmin,rmax
c--------------------------------------------------------- input array
! input parameter
A=102.0d+00
Z=50.0d+00
R0=1.25*A**0.3333333
write(*,*) 'nuclear radius' ,R0
am=931.5*A*4.0/(A+4.0) !average mass
write(*,*) 'average mass', am
xmin = 0.1d+00
xmax = 100.0d+00
! generate xi,yi
step =(xmax-xmin)/float(nc-1)
do i=1,ni
r=xmin+step*dfloat(i)
call compute_Vn(r,Vn)
call compute_Vc(r,Vc)
call compute_Vl(r,Vl)
xi(i)=step*dfloat(i)
yi(i)=Vl+Vc+Vn
end do
! output data into a file
open(1,file='out_table.dat',status='unknown')
do i=1,ni
write(1,*) xi(i),yi(i)
end do
close(1)
write(*,*) 'array is created'
c--------------------------------------------------------first integrate
! step 2: calculate derivatives
rmin=xi(1)
rmax=xi(ni)
yp1 = deriv3(xmin, xi, yi, ni, 1)
ypn = deriv3(xmax, xi, yi, ni, 1)
write (*,201) rmin,yp1,rmax,ypn
201 format (7f12.5)
c------------------------------------------------------------spline
call spline(xi,yi,ni,yp1,ypn,y2)
c-----------------------------------------------------splint + integrate
step1=nint((4.587679-0.022480)/ni)
do i =1,(step1-1)
r=step*dfloat(i)
call splint(xi,yi,y2,ni,r,y)
write (*,*) y
end do
end program test_28_08
c-----------------------subroutine------------------------
! Function f(x) to get derivation
function deriv3(xx, xi, yi, ni, m)
implicit none
integer, parameter :: n=3
double precision deriv3, xx
integer ni, m
double precision xi(ni), yi(ni)
double precision x(n), f(n)
integer i, j, k, ix
! exit if too high-order derivative was needed,
if (m > 2) then
deriv3 = 0.0
return
end if
! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
if (xx < xi(1) .or. xx > xi(ni)) then
deriv3 = 0.0
return
end if
! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
i = 1
j = ni
do while (j > i+1)
k = (i+j)/2
if (xx < xi(k)) then
j = k
else
i = k
end if
end do
! shift i that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
i = i + 1 - n/2
! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1
! just wanted to use index i
ix = i
! initialization of f(n) and x(n)
do i=1,n
f(i) = yi(ix+i-1)
x(i) = xi(ix+i-1)
end do
! calculate the first-order derivative using Lagrange interpolation
if (m == 1) then
deriv3 =(2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-
1 x(1))*(x(2)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-
2 x(1))*(x(3)-x(2)))
! calculate the second-order derivative using Lagrange interpolation
else
deriv3 = 2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
end if
end function deriv3
c-------------------------------------COMPUTE POTENTIAL WITH CONDITIONS
!--------------------------------------------------------------------
subroutine compute_Vn(r,Vn)
implicit none
double precision r,Vn
real::A,Z,R0,am
common A,Z,R0,am
Vn=-50.0/(1.0+exp((r-R0)/0.5))
end subroutine compute_Vn
!----------------------------------------------------------------------
subroutine compute_Vc(r,Vc)
implicit none
double precision:: r,phi,up,Vc
real::A,Z,R0,am
common A,Z,R0,am
if (r.lt.R0) then
Vc=(Z*1.44/R0)*(3.0-(r/R0)**2.0)
else if (r.gt.R0) then
Vc=2.88*Z/r
else if (r.eq.R0) then
phi=1.5d0/R0
up=phi*r
Vc=(Z*2.88/R0)*(1-exp(-up-0.5*(up**2.0)-0.35*(up**3.0)))
end if
end subroutine compute_Vc
!----------------------------------------------------------------------
subroutine compute_Vl(rr,Vl)
implicit none
double precision rr,Vl
real::A,Z,R0,am
common A,Z,R0,am
real,parameter:: hbarc =197.33
Vl=((hbarc)**2.0)/(8.0*am*rr**2.0)
end subroutine compute_Vl
c---------------------------------
c----------------------------------------------------------------spline from F95. recipe
SUBROUTINE spline(x,y,n,yp1,ypn,y2)
implicit none
INTEGER n,NMAX
REAL yp1,ypn,x(n),y(n),y2(n)
PARAMETER (NMAX=1000)
INTEGER i,k
REAL p,qn,sig,un,u(NMAX)
if (yp1.gt..99e30) then
y2(1)=0.
u(1)=0.
else
y2(1)=-0.5
u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
endif
do i=2,n-1
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
p=sig*y2(i-1)+2.
y2(i)=(sig-1.)/p
u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-
1 y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
end do
if (ypn.gt..99e30) then
qn=0.
un=0.
else
qn=0.5
un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
endif
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
do k=n-1,1,-1
y2(k)=y2(k)*y2(k+1)+u(k)
enddo
return
END SUBROUTINE spline
C--------------------------------------------------------------splint from F95. recipe
SUBROUTINE splint(xa,ya,y2a,n,x,y)
implicit none
INTEGER n
REAL x,y,xa(n),y2a(n),ya(n)
INTEGER k,khi,klo
REAL a,b,h
klo=1
khi=n
1 if (khi-klo.gt.1) then
k=(khi+klo)/2
if(xa(k).gt.x)then
khi=k
else
klo=k
endif
goto 1
endif
h=xa(khi)-xa(klo)
if (h.eq.0.) stop
a=(xa(khi)-x)/h
b=(x-xa(klo))/h
y=a*ya(klo)+b*ya(khi)+
2 ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
return
END subroutine splint
