interpo.f

来自「计算能量损失,重离子经过介质后的能量沉积LET,射程」· F 代码 · 共 82 行

*$debug
*$declare

	subroutine interpo(n,a,f1,x,y1,ierr,m,eps)
	real	v1(20,20),b(20)
	real	a(*),f1(*)

	real	s,bv,vv1,el,ec
	integer	m2,is,ib,im,il,ir,i,j,k

	m2=2*m
	is=1
	ib=n
1	continue
		im=(is+ib)/2
		s=sign(1.,a(n)-a(1))
********		write(*,*) is,ib,im,a(im),x
		if((x-a(im))*s.lt.0)then
		   ib=im
		else if((x-a(im))*s.gt.0)then
		   is=im
		else
		   y1=f1(im)
		   ierr=0
		   return
		end if
	if(is.lt.ib-1) go to 1

	il=is-m+1
	ir=ib+m-1
********	write(*,*) il,ir
	if(il.lt.1) then
	   il=1
	else if (ir.gt.n) then
	   il=n-m2+1
	end if
********	write(*,*) il,ir

	do i=1,m2
		b(i)=a(i+il-1)-x
		v1(1,i)=f1(i+il-1)
	enddo

	do i=2,m2
	    k=i-1
	    do j=i,m2
		if(abs(b(j)).lt.abs(b(k))) k=j
	    enddo
	    if(k.ne.i-1) then
		bv=b(k)
		b(k)=b(i-1)
		b(i-1)=bv
		vv1=v1(1,k)
		v1(1,k)=v1(1,i-1)
		v1(1,i-1)=vv1
	    endif
	enddo

	ierr=-9999
	el=1
	do i=2,m2
	    do j=2,i
	      v1(j,i)=(v1(j-1,j-1)*b(i)-v1(j-1,i)*b(j-1))/(b(i)-b(j-1))
	    enddo
	    ec=abs(1.-v1(i-1,i-1)/v1(i,i))
	    if(ec.le.eps) then
		y1=v1(i,i)
		ierr=i
		return
	    elseif(ec.gt.el)then
	    	y1=v1(i-1,i-1)
		ierr=1-i
		return
	    endif
	    el=ec
	enddo

	if(ierr.ne.0) y1=v1(m2,m2)

	return
	end