program rbf
! rbf.f90
! This file creates a Radial basis function matrix.
! 
!
! We store this matrix in matrixfile.dat
! 

integer :: m,n,d,i,j,init(4)
real(kind=8) :: l,r,coef
real(kind=8), allocatable, dimension(:) :: xx,yy,value,value2
real(kind=8), allocatable, dimension(:,:) :: A
real        ::rand

open(unit=20,file='inputfile.dat',status='old')
read(20,*)init
n=init(1)
close(20)
m=n-3
allocate(value(2*m),value2(m),xx(m),yy(m),A(n,n))
open(unit=30,file='vecrand.dat',status='old')
read(30,*)value
do i=1,m
    xx(i)=value(2*i-1)
    yy(i)=value(2*i)
end do
close(30)
do i=1,m
    d=0
    do j=1,m
        r=dsqrt((xx(i)-xx(j))**2+(yy(j)-yy(i))**2)
        if (r>0) then
            A(i,j)=r*r*log(r)
            if (d==0)then
                d=j
                coef=0.5*d*rand(0)
                A(i,i)=coef*coef*log(coef)
            end if
        end if
    end do
end do
do j=1,m
    A(m+1,j)=xx(j)
    A(m+2,j)=yy(j)
    A(m+3,j)=1.d0
    A(j,m+1)=xx(j)
    A(j,m+2)=yy(j)
    A(j,m+3)=1.d0
end do
do i=m+1,n
    do j=m+1,n
        A(i,j)=0.d0
    end do
end do
deallocate(xx,yy,value,value2)
open(unit=10,file='matrixfile.dat',status='unknown')
do i=1,n
    write(10,*) (A(i,j),j=1,n)
end do
close(10)
deallocate(A)
end program rbf