require "numru/ssl2"
include NumRu

narray = SSL2.with_narray?

def gen(x,y,fxy,narray)
  nx = x.length
  ny = y.length
  if narray then 
    fxy[0,0] = 9*x[0]*x[0]*y[0]*y[0]
    fxy[nx+1,0] = 9*x[nx-1]*x[nx-1]*y[0]*y[0]
    fxy[0,ny+1] = 9*x[0]*x[0]*y[ny-1]*y[ny-1]
    fxy[nx+1,ny+1] = 9*x[nx-1]*x[nx-1]*y[ny-1]*y[ny-1]
  else
    fxy[0][0] = 9*x[0]*x[0]*y[0]*y[0]
    fxy[0][nx+1] = 9*x[nx-1]*x[nx-1]*y[0]*y[0]
    fxy[ny+1][0] = 9*x[0]*x[0]*y[ny-1]*y[ny-1]
    fxy[ny+1][nx+1] = 9*x[nx-1]*x[nx-1]*y[ny-1]*y[ny-1]
  end

  x1 = x[0]*x[0]*3
  xn = x[nx-1]*x[nx-1]*3
  y1 = y[0]*y[0]*3
  yn = y[ny-1]*y[ny-1]*3

  for i in 0..nx-1
    fx = x[i]*x[i]*x[i]
    if narray then
      fxy[i+1,0] = y1*fx
      fxy[i+1,ny+1] = yn*fx
    else
      fxy[0][i+1] = y1*fx
      fxy[ny+1][i+1] = yn*fx
    end
  end
  for j in 0..ny-1
    fy = y[j]*y[j]*y[j]
    if narray then
      fxy[0,j+1] = x1*fy
      fxy[nx+1,j+1] = xn*fy
    else
      fxy[j+1][0] = x1*fy
      fxy[j+1][nx+1] = xn*fy
    end
  end

  for i in 0..nx-1
    fx = x[i]*x[i]*x[i]
    for j in 0..ny-1
      fy = y[j]*y[j]*y[j]
      if narray then
	fxy[i+1,j+1] = fx*fy
      else
	fxy[j+1][i+1] = fx*fy
      end
    end
  end

end

nx=10
ny=10
m=3
if narray then
  x = NArray.sfloat(nx)
  y = NArray.sfloat(nx)
  fxy = NArray.sfloat(nx+m-1,ny+m-1)
else
  x = Array.new(nx)
  y = Array.new(nx)
  fxy = Array.new(ny+m-1)
  for i in 0..fxy.length-1
    fxy[i] = Array.new(nx+m-1)
  end
end

px=0
hx=1.0/(nx-1)
for i in 0..nx-1
  x[i] = px+i*hx
end
py=0
hy=1.0/(ny-1)
for j in 0..ny-1
  y[j] = py+j*hy
end

gen(x,y,fxy,narray)

c = SSL2.bicd1(x,y,fxy,m)

ix=4
vx=0.5
iy=4
vy=0.5
for iswx in -1..m
  iswy=iswx
  ix,iy,f = SSL2.bifd1(x,y,m,c,iswx,vx,ix,iswy,vy,iy)  # ix and iy can be nil
  case iswx
  when -1
    printf("integral = %12.6e\n",f)
  when 0
    printf("value = %12.6e\n",f)
  else
    printf("derivative %i = %12.6e\n",iswx,f)
  end
end

vx = vy = NArray.sfloat(3,2).indgen!.mul!(0.1)
iswy= iswx = 0
ix,iy,f = SSL2.bifd1(x,y,m,c,iswx,vx,nil,iswy,vy,nil)
puts "multiple interpolation:"
p f