*** splin2NEW.c	Thu Jan 28 00:57:11 1993
--- splin2.c	Wed Jan 27 00:46:47 1993
***************
*** 1,37 ****
! void splin2(x1a,x2a,ya,y2a,m,n,x1,x2,y)
! float x1a[],x2a[],**ya,**y2a,x1,x2,*y;
! int m,n;
  {
! 	int j;
! 	float *ytmp,*yytmp,*vector();
! 	void spline(),splint(),free_vector();
  
! 	ytmp=vector(1,n);
! 	yytmp=vector(1,n);
  	for (j=1;j<=m;j++)
! 		splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);
! 	spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp);
! 	splint(x1a,yytmp,ytmp,m,x1,y);
! 	free_vector(yytmp,1,n);
! 	free_vector(ytmp,1,n);
  }
- void splie2(x1a,x2a,ya,m,n,y2a)
- float x1a[],x2a[],**ya,**y2a;
- int m,n;
- {
- 	int j;
- 	void spline();
  
  	for (j=1;j<=m;j++)
  		spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);
  }
! void splint(xa,ya,y2a,n,x,y)
! float xa[],ya[],y2a[],x,*y;
! int n;
  {
! 	int klo,khi,k;
! 	float h,b,a;
! 	void nrerror();
  
  	klo=1;
  	khi=n;
--- 1,60 ----
! /* 2-D bicubic spline from Numerical Recipes in C by Press. W.H et al. */
! /* ###### BUG NOTE ######
!    NOTE that there is a bug in Numerical Recipes (in both 1-st and 2-nd editions)
!    in file splin2.c which is used here with the bug fixed.
! */
! 
! #include <stdio.h>
! #include <stdlib.h>
! #include "splin2.h"
! /* #include "version.h" */
! 
! #ifndef APPNAME
! #define APPNAME "ContourView"
! #endif
! 
! /* float *ytmp,*yytmp;  for debuggin */
! 
! /* Given x1a, x2a, ya (Data), and y2a (computed by ssplie2),
!  * this routine computes interpolated value y at point x1, x2 by
!  * bicubic spline.  From Press et al., Num Rec for C.
!  * Original in the book is in error.  See the code below.
!  */
! void splin2(float *x1a, float *x2a, float **ya, float **y2a,
! 		int m, int n, float x1, float x2, float *y)
  {
! int j;
! float *ytmp,*yytmp;
  
! 	ytmp=fvector(1, m);	/* was fvector(1, n) which is incorrect */
! 	yytmp=fvector(1, m);	/* was fvector(1, n) which is incorrect */
! 
  	for (j=1;j<=m;j++)
! 		splint(x2a, ya[j], y2a[j], n, x2, &yytmp[j]);
! 	spline(x1a, yytmp, m, 1.0e30, 1.0e30, ytmp);
! 	splint(x1a, yytmp, ytmp, m, x1, y);
! 	free_fvector(yytmp, 1, m);	/* was free_fvector(yytmp, 1, n) which is incorrect */
! 	free_fvector(ytmp, 1, m);	/* was free_fvector(ytmp, 1, n) which is incorrect */
  }
  
+ 
+ /* From Press, et al.  Num Rec for C.
+  *	Returns second derivatives in array y2a[1..m][1..n]
+  */
+ void splie2(float *x1a, float *x2a, float **ya, int m, int n, float **y2a)
+ {
+ int j;
  	for (j=1;j<=m;j++)
  		spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);
  }
! 
! 
! 
! 
! void splint(float *xa, float *ya, float *y2a, int n, float x, float *y)
  {
! int klo,khi,k;
! float h, b, a;
  
  	klo=1;
  	khi=n;
***************
*** 41,60 ****
  		else klo=k;
  	}
  	h=xa[khi]-xa[klo];
! 	if (h == 0.0) nrerror("Bad XA input to routine SPLINT");
  	a=(xa[khi]-x)/h;
  	b=(x-xa[klo])/h;
! 	*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
  }
! void spline(x,y,n,yp1,ypn,y2)
! float x[],y[],yp1,ypn,y2[];
! int n;
  {
! 	int i,k;
! 	float p,qn,sig,un,*u,*vector();
! 	void free_vector();
  
! 	u=vector(1,n-1);
  	if (yp1 > 0.99e30)
  		y2[1]=u[1]=0.0;
  	else {
--- 64,86 ----
  		else klo=k;
  	}
  	h=xa[khi]-xa[klo];
! 	if (h == 0.0) {
! 		fprintf(stderr,"%s: Bad XA input to routine splint()\n", APPNAME);
! 		exit(1);
! 	}
  	a=(xa[khi]-x)/h;
  	b=(x-xa[klo])/h;
! 	*y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
  }
! 
! 
! 
! void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
  {
! int i,k;
! float p,qn,sig,un,*u;
  
! 	u=fvector(1,n-1);
  	if (yp1 > 0.99e30)
  		y2[1]=u[1]=0.0;
  	else {
***************
*** 77,290 ****
  	y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  	for (k=n-1;k>=1;k--)
  		y2[k]=y2[k]*y2[k+1]+u[k];
! 	free_vector(u,1,n-1);
  }
  
- #include <malloc.h>
- #include <stdio.h>
  
! void nrerror(error_text)
! char error_text[];
! {
! 	void exit();
  
! 	fprintf(stderr,"Numerical Recipes run-time error...\n");
! 	fprintf(stderr,"%s\n",error_text);
! 	fprintf(stderr,"...now exiting to system...\n");
! 	exit(1);
! }
  
! 
! 
! float *vector(nl,nh)
! int nl,nh;
  {
  	float *v;
! 
! 	v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
! 	if (!v) nrerror("allocation failure in vector()");
! 	return v-nl;
  }
  
- int *ivector(nl,nh)
- int nl,nh;
- {
- 	int *v;
  
! 	v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
! 	if (!v) nrerror("allocation failure in ivector()");
! 	return v-nl;
! }
! 
! double *dvector(nl,nh)
! int nl,nh;
  {
! 	double *v;
! 
! 	v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
! 	if (!v) nrerror("allocation failure in dvector()");
! 	return v-nl;
  }
  
  
! 
! float **matrix(nrl,nrh,ncl,nch)
! int nrl,nrh,ncl,nch;
  {
! 	int i;
! 	float **m;
! 
! 	m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*));
! 	if (!m) nrerror("allocation failure 1 in matrix()");
! 	m -= nrl;
! 
! 	for(i=nrl;i<=nrh;i++) {
  		m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
! 		if (!m[i]) nrerror("allocation failure 2 in matrix()");
! 		m[i] -= ncl;
  	}
  	return m;
  }
  
- double **dmatrix(nrl,nrh,ncl,nch)
- int nrl,nrh,ncl,nch;
- {
- 	int i;
- 	double **m;
  
! 	m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
! 	if (!m) nrerror("allocation failure 1 in dmatrix()");
! 	m -= nrl;
! 
! 	for(i=nrl;i<=nrh;i++) {
! 		m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
! 		if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
! 		m[i] -= ncl;
! 	}
! 	return m;
! }
! 
! int **imatrix(nrl,nrh,ncl,nch)
! int nrl,nrh,ncl,nch;
  {
- 	int i,**m;
- 
- 	m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
- 	if (!m) nrerror("allocation failure 1 in imatrix()");
- 	m -= nrl;
- 
- 	for(i=nrl;i<=nrh;i++) {
- 		m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
- 		if (!m[i]) nrerror("allocation failure 2 in imatrix()");
- 		m[i] -= ncl;
- 	}
- 	return m;
- }
- 
- 
- 
- float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
- float **a;
- int oldrl,oldrh,oldcl,oldch,newrl,newcl;
- {
- 	int i,j;
- 	float **m;
- 
- 	m=(float **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(float*));
- 	if (!m) nrerror("allocation failure in submatrix()");
- 	m -= newrl;
- 
- 	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;
- 
- 	return m;
- }
- 
- 
- 
- void free_vector(v,nl,nh)
- float *v;
- int nl,nh;
- {
- 	free((char*) (v+nl));
- }
- 
- void free_ivector(v,nl,nh)
- int *v,nl,nh;
- {
- 	free((char*) (v+nl));
- }
- 
- void free_dvector(v,nl,nh)
- double *v;
- int nl,nh;
- {
- 	free((char*) (v+nl));
- }
- 
- 
- 
- void free_matrix(m,nrl,nrh,ncl,nch)
- float **m;
- int nrl,nrh,ncl,nch;
- {
  	int i;
  
- 	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
- 	free((char*) (m+nrl));
  }
  
- void free_dmatrix(m,nrl,nrh,ncl,nch)
- double **m;
- int nrl,nrh,ncl,nch;
- {
- 	int i;
  
- 	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
- 	free((char*) (m+nrl));
- }
  
- void free_imatrix(m,nrl,nrh,ncl,nch)
- int **m;
- int nrl,nrh,ncl,nch;
- {
- 	int i;
- 
- 	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
- 	free((char*) (m+nrl));
- }
- 
- 
- 
- void free_submatrix(b,nrl,nrh,ncl,nch)
- float **b;
- int nrl,nrh,ncl,nch;
- {
- 	free((char*) (b+nrl));
- }
- 
- 
- 
- float **convert_matrix(a,nrl,nrh,ncl,nch)
- float *a;
- int nrl,nrh,ncl,nch;
- {
- 	int i,j,nrow,ncol;
- 	float **m;
- 
- 	nrow=nrh-nrl+1;
- 	ncol=nch-ncl+1;
- 	m = (float **) malloc((unsigned) (nrow)*sizeof(float*));
- 	if (!m) nrerror("allocation failure in convert_matrix()");
- 	m -= nrl;
- 	for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
- 	return m;
- }
- 
- 
- 
- void free_convert_matrix(b,nrl,nrh,ncl,nch)
- float **b;
- int nrl,nrh,ncl,nch;
- {
- 	free((char*) (b+nrl));
- }
--- 103,171 ----
  	y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  	for (k=n-1;k>=1;k--)
  		y2[k]=y2[k]*y2[k+1]+u[k];
! 	free_fvector(u,1,n-1);
  }
  
  
! /* memory allocation routines */
  
! /* ###### 1-d float */
  
! float *fvector(int nl, int nh)
  {
  	float *v;
! 	v=(float *)malloc((unsigned)(nh-nl+1)*sizeof(float));
! 	if(!v) fprintf(stderr, "allocation failure in fvector()\n");
! 	return(v-nl);
  }
  
  
! void free_fvector(float *v, int nl, int nh)
  {
! 	free((void *)(v+nl));
  }
  
  
! /* ##### 2-d float */
!  
! float **fmatrix(int nrl, int nrh, int ncl, int nch)
  {
!    int i=2,j;
!    int error=0;
!    float **m;
! 	m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float *));
! 	if (!m) fprintf(stderr, "allocation failure 1 in fmatrix()\n");
! 	else{
! 	   m -= nrl;
! 	   i=nrl;
! 	   do{
  		m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
! 		if (!m[i]) {
! 			fprintf(stderr, "allocation failure 2 in fmatrix()\n");
! 			error=1;
! 		}else
! 			m[i] -= ncl;
! 		i++;
!            }while( (i<=nrh) && (!error));
  	}
+ 	if(error){
+ 		fprintf(stderr, "Freeing allocated submatrices...\n");
+ 		for(j=i-2;j>=nrl;j--)	       /* free those allocated */
+ 			free((void *)(m[j]+ncl));
+ 		free((void *)(m+nrl));
+ 		m=NULL;
+  	}
  	return m;
  }
  
  
! void free_fmatrix(float **m, int nrl, int nrh, int ncl, int nch)
  {
  	int i;
+ 	for(i=nrh;i>=nrl;i--)free((void *)(m[i]+ncl));
+ 	free((void *) (m+nrl));
  
  }