polint.c
No OneTemporary
Actions

Subscribers

None

File Metadata

	#include <math.h>
	#include "lt.h"

	void polint(double xa,double ya,int n,double x,double y,double dy)
	{
	int i,m,ns=1;
	double den,dif,dift,ho,hp,w;
	double c,d;

	dif=fabs(x-xa[1]);
	c=vector(1,n);
	d=vector(1,n);
	for (i=1;i<=n;i++) {
	if ( (dift=fabs(x-xa[i])) < dif) {
	ns=i;
	dif=dift;
	}
	c[i]=ya[i];
	d[i]=ya[i];
	}
	*y=ya[ns--];
	for (m=1;m<n;m++) {
	for (i=1;i<=n-m;i++) {
	ho=xa[i]-x;
	hp=xa[i+m]-x;
	w=c[i+1]-d[i];
	// Change the error limit from == 0.0 to < 1e-10 EJ030806
	if ( (den=ho-hp) == 0 ) nrerror("Error in routine POLINT");
	den=w/den;
	d[i]=hp*den;
	c[i]=ho*den;
	}
	y += (dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
	}
	free_vector(d,1);
	free_vector(c,1);
	}