/* bicubic spline interpolating f(x1,x2)
 * using functions in recipes which is in FORTRAN convention
 */

#include <stdlib.h>
#include "message.h"
#include "typedef.h"

float **farray2();
void splie2();
void splin2();

float bicubinterp(x1,x2,f,dim,xx1,xx2,iffirst)
float *x1,*x2,**f;
float xx1,xx2;
Milmt2 *dim;  /* dimension of the matrix */
int iffirst;   /* =1: recalculating the spline */
{
static float **f2;
static int mm=0, nn=0;
float ff;
static float *x1a,*x2a,**fa;  /* C-- FORTRAN */
int dm,dn,i,m1,m2,n1;
	
	dm=1+dim->m2-dim->m1;
	dn=1+dim->n2-dim->n1;
	m1=dim->m1;
	m2=dim->m2;
	n1=dim->n1;
	
	if(dm<2||dn<2) stop("bicubinterp: illegal dimension");
	
	if(dm>mm)
	{ if(mm>0) free(++fa);
	  fa=((float **) malloc(dm*sizeof(float*)))-1;
	}
	
	x1a=x1-1-m1;
	x2a=x2-1-n1;
	for(i=m1;i<=m2;i++) 
	  fa[i+1-m1]=f[i]-1-m1;

	if(iffirst)
	{ if(dm>mm || dn>nn)
	  { if(mm>0 ||nn>0) free_farray2(f2,1,dm,1);
	    f2=farray2(1,dm,1,dn);
	    mm=dm;
	    nn=dn;
	  }
	  splie2(x1a,x2a,fa,dm,dn,f2);
	}
	
	splin2(x1a,x2a,fa,f2,dm,dn,xx1,xx2,&ff);
	
	return(ff);
		
}