📄 spline3.cs
字号:
{
dxjp1 = x[j]-x[j-1];
dyjp1 = y[j]-y[j-1];
dxp = dxj+dxjp1;
ctbl[1,j-1] = dxjp1/dxp;
ctbl[2,j-1] = 1-ctbl[1,j-1];
ctbl[3,j-1] = 6*(dyjp1/dxjp1-dyj/dxj)/dxp;
dxj = dxjp1;
dyj = dyjp1;
j = j+1;
}
}
ctbl[1,0] = -(b1/2);
ctbl[2,0] = b2/2;
if( n!=2 )
{
j = 2;
while( j<=nxm1 )
{
pj = ctbl[2,j-1]*ctbl[1,j-2]+2;
ctbl[1,j-1] = -(ctbl[1,j-1]/pj);
ctbl[2,j-1] = (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,j-2])/pj;
j = j+1;
}
}
yppb = (b4-b3*ctbl[2,nxm1-1])/(b3*ctbl[1,nxm1-1]+2);
i = 1;
while( i<=nxm1 )
{
j = n-i;
yppa = ctbl[1,j-1]*yppb+ctbl[2,j-1];
dx = x[j]-x[j-1];
ctbl[3,j-1] = (yppb-yppa)/dx/6;
ctbl[2,j-1] = yppa/2;
ctbl[1,j-1] = (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]*dx)*dx;
yppb = yppa;
i = i+1;
}
for(i=1; i<=n; i++)
{
ctbl[0,i-1] = y[i-1];
ctbl[4,i-1] = x[i-1];
}
}
}
/*************************************************************************
Obsolete subroutine, left for backward compatibility.
*************************************************************************/
public static double spline3interpolate(int n,
ref double[,] c,
double x)
{
double result = 0;
int i = 0;
int l = 0;
int half = 0;
int first = 0;
int middle = 0;
n = n-1;
l = n;
first = 0;
while( l>0 )
{
half = l/2;
middle = first+half;
if( c[4,middle]<x )
{
first = middle+1;
l = l-half-1;
}
else
{
l = half;
}
}
i = first-1;
if( i<0 )
{
i = 0;
}
result = c[0,i]+(x-c[4,i])*(c[1,i]+(x-c[4,i])*(c[2,i]+c[3,i]*(x-c[4,i])));
return result;
}
/*************************************************************************
Internal subroutine. Heap sort.
*************************************************************************/
private static void heapsortpoints(ref double[] x,
ref double[] y,
int n)
{
int i = 0;
int j = 0;
int k = 0;
int t = 0;
double tmp = 0;
bool isascending = new bool();
bool isdescending = new bool();
//
// Test for already sorted set
//
isascending = true;
isdescending = true;
for(i=1; i<=n-1; i++)
{
isascending = isascending & x[i]>x[i-1];
isdescending = isdescending & x[i]<x[i-1];
}
if( isascending )
{
return;
}
if( isdescending )
{
for(i=0; i<=n-1; i++)
{
j = n-1-i;
if( j<=i )
{
break;
}
tmp = x[i];
x[i] = x[j];
x[j] = tmp;
tmp = y[i];
y[i] = y[j];
y[j] = tmp;
}
return;
}
//
// Special case: N=1
//
if( n==1 )
{
return;
}
//
// General case
//
i = 2;
do
{
t = i;
while( t!=1 )
{
k = t/2;
if( x[k-1]>=x[t-1] )
{
t = 1;
}
else
{
tmp = x[k-1];
x[k-1] = x[t-1];
x[t-1] = tmp;
tmp = y[k-1];
y[k-1] = y[t-1];
y[t-1] = tmp;
t = k;
}
}
i = i+1;
}
while( i<=n );
i = n-1;
do
{
tmp = x[i];
x[i] = x[0];
x[0] = tmp;
tmp = y[i];
y[i] = y[0];
y[0] = tmp;
t = 1;
while( t!=0 )
{
k = 2*t;
if( k>i )
{
t = 0;
}
else
{
if( k<i )
{
if( x[k]>x[k-1] )
{
k = k+1;
}
}
if( x[t-1]>=x[k-1] )
{
t = 0;
}
else
{
tmp = x[k-1];
x[k-1] = x[t-1];
x[t-1] = tmp;
tmp = y[k-1];
y[k-1] = y[t-1];
y[t-1] = tmp;
t = k;
}
}
}
i = i-1;
}
while( i>=1 );
}
/*************************************************************************
Internal subroutine. Heap sort.
*************************************************************************/
private static void heapsortdpoints(ref double[] x,
ref double[] y,
ref double[] d,
int n)
{
int i = 0;
int j = 0;
int k = 0;
int t = 0;
double tmp = 0;
bool isascending = new bool();
bool isdescending = new bool();
//
// Test for already sorted set
//
isascending = true;
isdescending = true;
for(i=1; i<=n-1; i++)
{
isascending = isascending & x[i]>x[i-1];
isdescending = isdescending & x[i]<x[i-1];
}
if( isascending )
{
return;
}
if( isdescending )
{
for(i=0; i<=n-1; i++)
{
j = n-1-i;
if( j<=i )
{
break;
}
tmp = x[i];
x[i] = x[j];
x[j] = tmp;
tmp = y[i];
y[i] = y[j];
y[j] = tmp;
tmp = d[i];
d[i] = d[j];
d[j] = tmp;
}
return;
}
//
// Special case: N=1
//
if( n==1 )
{
return;
}
//
// General case
//
i = 2;
do
{
t = i;
while( t!=1 )
{
k = t/2;
if( x[k-1]>=x[t-1] )
{
t = 1;
}
else
{
tmp = x[k-1];
x[k-1] = x[t-1];
x[t-1] = tmp;
tmp = y[k-1];
y[k-1] = y[t-1];
y[t-1] = tmp;
tmp = d[k-1];
d[k-1] = d[t-1];
d[t-1] = tmp;
t = k;
}
}
i = i+1;
}
while( i<=n );
i = n-1;
do
{
tmp = x[i];
x[i] = x[0];
x[0] = tmp;
tmp = y[i];
y[i] = y[0];
y[0] = tmp;
tmp = d[i];
d[i] = d[0];
d[0] = tmp;
t = 1;
while( t!=0 )
{
k = 2*t;
if( k>i )
{
t = 0;
}
else
{
if( k<i )
{
if( x[k]>x[k-1] )
{
k = k+1;
}
}
if( x[t-1]>=x[k-1] )
{
t = 0;
}
else
{
tmp = x[k-1];
x[k-1] = x[t-1];
x[t-1] = tmp;
tmp = y[k-1];
y[k-1] = y[t-1];
y[t-1] = tmp;
tmp = d[k-1];
d[k-1] = d[t-1];
d[t-1] = tmp;
t = k;
}
}
}
i = i-1;
}
while( i>=1 );
}
/*************************************************************************
Internal subroutine. Tridiagonal solver.
*************************************************************************/
private static void solvetridiagonal(double[] a,
double[] b,
double[] c,
double[] d,
int n,
ref double[] x)
{
int k = 0;
double t = 0;
a = (double[])a.Clone();
b = (double[])b.Clone();
c = (double[])c.Clone();
d = (double[])d.Clone();
x = new double[n-1+1];
a[0] = 0;
c[n-1] = 0;
for(k=1; k<=n-1; k++)
{
t = a[k]/b[k-1];
b[k] = b[k]-t*c[k-1];
d[k] = d[k]-t*d[k-1];
}
x[n-1] = d[n-1]/b[n-1];
for(k=n-2; k>=0; k--)
{
x[k] = (d[k]-c[k]*x[k+1])/b[k];
}
}
/*************************************************************************
Internal subroutine. Three-point differentiation
*************************************************************************/
private static double diffthreepoint(double t,
double x0,
double f0,
double x1,
double f1,
double x2,
double f2)
{
double result = 0;
double a = 0;
double b = 0;
t = t-x0;
x1 = x1-x0;
x2 = x2-x0;
a = (f2-f0-x2/x1*(f1-f0))/(AP.Math.Sqr(x2)-x1*x2);
b = (f1-f0-a*AP.Math.Sqr(x1))/x1;
result = 2*a*t+b;
return result;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -