/* bicubic spline interpolating for evenly distributed grids */

#include "../iolib/message.h"
#include "typedef.h"

float bicubinterp();
float *farray1();

float bicubinterpe(dim,coord,f,x,y,iffirst)
Milmt2 *dim;  /* dimension of the matrix */
Mflmt2 *coord; /* corresponding coordinates of x and y */
float **f; 
float x,y;
int iffirst;	/* =1: preparing xx, yy */
{
static int mm=0, nn=0;
static float *xx,*yy;
int dm,dn,m1,m2,n1,n2;
float xmin,xmax,ymin,ymax;
int i;
	
	dm=1+dim->m2-dim->m1;
	dn=1+dim->n2-dim->n1;
	m1=dim->m1;
	m2=dim->m2;
	n1=dim->n1;
	n2=dim->n2;
	xmin=coord->x1;
	xmax=coord->x2;
	ymin=coord->y1;
	ymax=coord->y2;
	
	if(dm<2||dn<2) STOP(bicubinterpe: illegal dimension);
	if(x<xmin || x>xmax) STOP(bicubinterpe: x out off range);
	if(y<ymin || y>ymax) STOP(bicubinterpe: y out off range);

	if(dm>mm)
	{ if(mm>0) free_farray1(xx,m1); 
	  xx=farray1(m1,m2);
	  mm=dm;
	}
	if(dn>nn)
	{ if(nn>0) free_farray1(yy,n1); 
	  yy=farray1(n1,n2);
	  nn=dn;
	}
	
	if(iffirst)
	{ for(i=m1;i<=m2;i++) 
	    xx[i]=xmin+(xmax-xmin)*(float)(i-m1)/(float)(m2-m1);
	  for(i=n1;i<=n2;i++) 
	    yy[i]=ymin+(ymax-ymin)*(float)(i-n1)/(float)(n2-n1);
	}
	
	return bicubinterp(xx,yy,f,dim,x,y,iffirst);
	
}