*-------------------------------------------------------------------------
* SP - various general purpose 2d spline interpolation routines
* based on STS_A model.
*-------------------------------------------------------------------------
subroutine sp1(a,x,y,n,ro,s,t,r,b,aa,c1,c2,c3,max)
implicit real(8) (a-g,o-z)
real(4) x(n),y(n),a(n)
real(8) s(n),t(3),ro
real(8) r(max*(max-1)/2),b(max),aa((max-3)*(max-2)/2)
real(8) c1(max),c2(max),c3(max)
*----------- computation of kernel matrix K (I,J) I<J -------------------
n3=n-3
do j=2,n
jm1=j-1
ij=(jm1*(j-2))/2
do i=1,jm1
dx=dble(x(i)-x(j))
dy=dble(y(i)-y(j))
a3=dx*dx+dy*dy
r(ij+i)=a3*dlog(a3)
c write(*,*)i,j,ij+i
end do
end do
*------------------- computation of C1,C2,C3(I) --------------------------
v1=dble(y(3)-y(1))
v4=dble(x(2)-x(1))
v2=dble(x(1)-x(3))
v3=dble(y(1)-y(2))
a1=v1*v4-v2*v3
v1=v1/a1
v2=v2/a1
v3=v3/a1
v4=v4/a1
do i=1,n3
ip3=i+3
a1=dble(x(1)-x(ip3))
a2=dble(y(1)-y(ip3))
c3(i)=v3*a1+v4*a2
c2(i)=v1*a1+v2*a2
c1(i)=-(1.d0+c2(i)+c3(i))
end do
*-------------------- computation of AA(I,J) I <= J ----------------------
do i=1,n3
i1=(i*(i+1))/2
ij=i1+i+2
a1=c2(i)*r(1)+c3(i)*r(2)+r(ij)
a2=c3(i)*r(3)+r(ij+1)
a3=ro*(c1(i)*c1(i)+c2(i)*c2(i)+c3(i)*c3(i)+1.d0)
aa(i1)=2.d0*(c1(i)*a1+c2(i)*a2+c3(i)*r(ij+2))+a3
end do
do j=2,n3
jm1=j-1
ij=(j*(j+3)+4)/2
it=(j*jm1)/2
do i=1,jm1
i1=it+i
j1=(i*(i+3)+4)/2
a1=c2(j)*r(1)+c3(j)*r(2)+r(ij)
a2=c1(j)*r(1)+c3(j)*r(3)+r(ij+1)
a3=c1(j)*r(2)+c2(j)*r(3)+r(ij+2)
a4=c1(j)*r(j1)+c2(j)*r(j1+1)+c3(j)*r(j1+2)+r(ij+i+2)
a5=ro*(c1(i)*c1(j)+c2(i)*c2(j)+c3(i)*c3(j))
aa(i1)=c1(i)*a1+c2(i)*a2+c3(i)*a3+a4+a5
end do
end do
*------------------ Cholesky decomposition -------------------------------
! print*,'sp11',(max-3)*(max-2)/2
call ludecp(aa,aa,n3,d1,d2,ier,(max-3)*(max-2)/2)
*---------- forward/back substitution for A ------------------------------
! print*,"sp115"
do i=1,n3
b(i)=c1(i)*a(1)+c2(i)*a(2)+c3(i)*a(3)+a(i+3)
end do
call luelmp(aa,b,n3,b,(max-3)*(max-2)/2)
! print*,"sp12"
do i=4,n
s(i)=b(i-3)
end do
*-------------------- compute A, spline coefficients ---------------------
a1=0.d0
a2=0.d0
a3=0.d0
do i=1,n3
ip3=i+3
a1=a1+c1(i)*s(ip3)
a2=a2+c2(i)*s(ip3)
a3=a3+c3(i)*s(ip3)
end do
s(1)=a1
s(2)=a2
s(3)=a3
a1=s(2)*r(1)+s(3)*r(2)+ro*s(1)
a2=s(1)*r(1)+s(3)*r(3)+ro*s(2)
a3=s(1)*r(2)+s(2)*r(3)+ro*s(3)
do i=4,n
i1=(i*(i-3)+4)/2
a1=a1+s(i)*r(i1)
a2=a2+s(i)*r(i1+1)
a3=a3+s(i)*r(i1+2)
end do
a1=a(1)-a1
a2=a(2)-(a2+a1)
a3=a(3)-(a3+a1)
t(1)=v1*a2+v3*a3
t(2)=v2*a2+v4*a3
t(3)=a1-(x(1)*t(1)+y(1)*t(2))
* ----------- return with spline coefficients in S & T ------------------
return
end
*-------------------------------------------------------------------------
* SP2 - solve for U at XI,YI given spline coefficients in S & T
* from SP1. X and Y are coordinate arrays same as those sent to SP1.
*-------------------------------------------------------------------------
subroutine sp2(xi,yi,x,y,n,s,t,u)
implicit real(8) (a-g,o-z)
real(4) x(n),y(n),xi,yi,u
real(8) s(n),t(3)
u=0.d0
*------------ interpolate at Xi,Yi from S(N) coefficients ---------------
do k=1,n
dx=dble(xi-x(k))
dy=dble(yi-y(k))
tt=dx*dx+dy*dy
if (tt.gt.1.d-30) then
dltt=dlog(tt)
!dxtt=dx*(dltt+1.d0)
!dytt=dy*(dltt+1.d0)
u=u+s(k)*tt*dltt
endif
end do
*----------------- add the contribution from the Ts ----------------------
u=u+(t(1)*xi+t(2)*yi+t(3))
! print*,u,t(1),t(2),t(3),xi,yi
return
end
subroutine sp2dx(xi,yi,x,y,n,s,t,u)
implicit real(8) (a-g,o-z)
real(4) x(n),y(n),xi,yi,u
real(8) s(n),t(3)
u=0.
*------------ interpolate at Xi,Yi from S(N) coefficients ---------------
do k=1,n
dx=dble(xi-x(k))
dy=dble(yi-y(k))
tt=dx*dx+dy*dy
if (tt.gt.1.d-30) then
dltt=dlog(tt)
dxtt=dx*(dltt+1.d0)
dytt=dy*(dltt+1.d0)
u=u+sngl(s(k)*(dx*(dltt+1.d0)))
endif
end do
*----------------- add the contribution from the Ts ----------------------
u=2.d0*u+t(1)
! u=u+sngl(t(1)*xi+t(2)*yi+t(3))
return
end
subroutine sp2dy(xi,yi,x,y,n,s,t,u)
implicit real(8) (a-g,o-z)
real(4) x(n),y(n),xi,yi,u
real(8) s(n),t(3)
u=0.
*------------ interpolate at Xi,Yi from S(N) coefficients ---------------
do k=1,n
dx=dble(xi-x(k))
dy=dble(yi-y(k))
tt=dx*dx+dy*dy
if (tt.gt.1.d-30) then
dltt=dlog(tt)
dxtt=dx*(dltt+1.d0)
dytt=dy*(dltt+1.d0)
u=u+sngl(s(k)*dy*(dltt+1.d0))
endif
end do
*----------------- add the contribution from the Ts ----------------------
u=2.d0*u+t(2)
! u=u+sngl(t(1)*xi+t(2)*yi+t(3))
return
end
*---------------------------------------------------------------------
* EXCHG - swap individual elements in a grid field data set accor-
* ding to NLOC markers in LOC1,2.
*---------------------------------------------------------------------
subroutine exchg(x,y,u,n,loc1,loc2,nloc)
real(4) x(n),y(n),u(n)
integer(4) loc1(nloc),loc2(nloc)
do i=1,nloc
tmp=x(loc1(i))
x(loc1(i))=x(loc2(i))
x(loc2(i))=tmp
tmp=y(loc1(i))
y(loc1(i))=y(loc2(i))
y(loc2(i))=tmp
tmp=u(loc1(i))
u(loc1(i))=u(loc2(i))
u(loc2(i))=tmp
end do
return
end