Context Navigation

source: trunk/sources/ALGLIB/spline3.cs @ 2430

Visit:

Last change on this file since 2430 was 2430, checked in by gkronber, 15 years ago
Moved ALGLIB code into a separate plugin. #783
File size: 40.2 KB

Line
1	/*************************************************************************
2	Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
3
4	>>> SOURCE LICENSE >>>
5	This program is free software; you can redistribute it and/or modify
6	it under the terms of the GNU General Public License as published by
7	the Free Software Foundation (www.fsf.org); either version 2 of the
8	License, or (at your option) any later version.
9
10	This program is distributed in the hope that it will be useful,
11	but WITHOUT ANY WARRANTY; without even the implied warranty of
12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13	GNU General Public License for more details.
14
15	A copy of the GNU General Public License is available at
16	http://www.fsf.org/licensing/licenses
17
18	>>> END OF LICENSE >>>
19	*************************************************************************/
20
21	using System;
22
23	namespace alglib
24	{
25	public class spline3
26	{
27	/*************************************************************************
28	This subroutine builds linear spline coefficients table.
29
30	Input parameters:
31	X - spline nodes, array[0..N-1]
32	Y - function values, array[0..N-1]
33	N - points count, N>=2
34
35	Output parameters:
36	C - coefficients table. Used by SplineInterpolation and other
37	subroutines from this file.
38
39	-- ALGLIB PROJECT --
40	Copyright 24.06.2007 by Bochkanov Sergey
41	*************************************************************************/
42	public static void buildlinearspline(double[] x,
43	double[] y,
44	int n,
45	ref double[] c)
46	{
47	int i = 0;
48	int tblsize = 0;
49
50	x = (double[])x.Clone();
51	y = (double[])y.Clone();
52
53	System.Diagnostics.Debug.Assert(n>=2, "BuildLinearSpline: N<2!");
54
55	//
56	// Sort points
57	//
58	heapsortpoints(ref x, ref y, n);
59
60	//
61	// Fill C:
62	// C[0] - length(C)
63	// C[1] - type(C):
64	// 3 - general cubic spline
65	// C[2] - N
66	// C[3]...C[3+N-1] - x[i], i = 0...N-1
67	// C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
68	//
69	tblsize = 3+n+(n-1)*4;
70	c = new double[tblsize-1+1];
71	c[0] = tblsize;
72	c[1] = 3;
73	c[2] = n;
74	for(i=0; i<=n-1; i++)
75	{
76	c[3+i] = x[i];
77	}
78	for(i=0; i<=n-2; i++)
79	{
80	c[3+n+4*i+0] = y[i];
81	c[3+n+4*i+1] = (y[i+1]-y[i])/(x[i+1]-x[i]);
82	c[3+n+4*i+2] = 0;
83	c[3+n+4*i+3] = 0;
84	}
85	}
86
87
88	/*************************************************************************
89	This subroutine builds cubic spline coefficients table.
90
91	Input parameters:
92	X - spline nodes, array[0..N-1]
93	Y - function values, array[0..N-1]
94	N - points count, N>=2
95	BoundLType - boundary condition type for the left boundary
96	BoundL - left boundary condition (first or second derivative,
97	depending on the BoundLType)
98	BoundRType - boundary condition type for the right boundary
99	BoundR - right boundary condition (first or second derivative,
100	depending on the BoundRType)
101
102	Output parameters:
103	C - coefficients table. Used by SplineInterpolation and
104	other subroutines from this file.
105
106	The BoundLType/BoundRType parameters can have the following values:
107	* 0, which corresponds to the parabolically terminated spline
108	(BoundL/BoundR are ignored).
109	* 1, which corresponds to the first derivative boundary condition
110	* 2, which corresponds to the second derivative boundary condition
111
112	-- ALGLIB PROJECT --
113	Copyright 23.06.2007 by Bochkanov Sergey
114	*************************************************************************/
115	public static void buildcubicspline(double[] x,
116	double[] y,
117	int n,
118	int boundltype,
119	double boundl,
120	int boundrtype,
121	double boundr,
122	ref double[] c)
123	{
124	double[] a1 = new double[0];
125	double[] a2 = new double[0];
126	double[] a3 = new double[0];
127	double[] b = new double[0];
128	double[] d = new double[0];
129	int i = 0;
130	int tblsize = 0;
131	double delta = 0;
132	double delta2 = 0;
133	double delta3 = 0;
134
135	x = (double[])x.Clone();
136	y = (double[])y.Clone();
137
138	System.Diagnostics.Debug.Assert(n>=2, "BuildCubicSpline: N<2!");
139	System.Diagnostics.Debug.Assert(boundltype==0 \| boundltype==1 \| boundltype==2, "BuildCubicSpline: incorrect BoundLType!");
140	System.Diagnostics.Debug.Assert(boundrtype==0 \| boundrtype==1 \| boundrtype==2, "BuildCubicSpline: incorrect BoundRType!");
141	a1 = new double[n-1+1];
142	a2 = new double[n-1+1];
143	a3 = new double[n-1+1];
144	b = new double[n-1+1];
145
146	//
147	// Special case:
148	// * N=2
149	// * parabolic terminated boundary condition on both ends
150	//
151	if( n==2 & boundltype==0 & boundrtype==0 )
152	{
153
154	//
155	// Change task type
156	//
157	boundltype = 2;
158	boundl = 0;
159	boundrtype = 2;
160	boundr = 0;
161	}
162
163	//
164	//
165	// Sort points
166	//
167	heapsortpoints(ref x, ref y, n);
168
169	//
170	// Left boundary conditions
171	//
172	if( boundltype==0 )
173	{
174	a1[0] = 0;
175	a2[0] = 1;
176	a3[0] = 1;
177	b[0] = 2*(y[1]-y[0])/(x[1]-x[0]);
178	}
179	if( boundltype==1 )
180	{
181	a1[0] = 0;
182	a2[0] = 1;
183	a3[0] = 0;
184	b[0] = boundl;
185	}
186	if( boundltype==2 )
187	{
188	a1[0] = 0;
189	a2[0] = 2;
190	a3[0] = 1;
191	b[0] = 3(y[1]-y[0])/(x[1]-x[0])-0.5boundl*(x[1]-x[0]);
192	}
193
194	//
195	// Central conditions
196	//
197	for(i=1; i<=n-2; i++)
198	{
199	a1[i] = x[i+1]-x[i];
200	a2[i] = 2*(x[i+1]-x[i-1]);
201	a3[i] = x[i]-x[i-1];
202	b[i] = 3(y[i]-y[i-1])/(x[i]-x[i-1])(x[i+1]-x[i])+3(y[i+1]-y[i])/(x[i+1]-x[i])(x[i]-x[i-1]);
203	}
204
205	//
206	// Right boundary conditions
207	//
208	if( boundrtype==0 )
209	{
210	a1[n-1] = 1;
211	a2[n-1] = 1;
212	a3[n-1] = 0;
213	b[n-1] = 2*(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
214	}
215	if( boundrtype==1 )
216	{
217	a1[n-1] = 0;
218	a2[n-1] = 1;
219	a3[n-1] = 0;
220	b[n-1] = boundr;
221	}
222	if( boundrtype==2 )
223	{
224	a1[n-1] = 1;
225	a2[n-1] = 2;
226	a3[n-1] = 0;
227	b[n-1] = 3(y[n-1]-y[n-2])/(x[n-1]-x[n-2])+0.5boundr*(x[n-1]-x[n-2]);
228	}
229
230	//
231	// Solve
232	//
233	solvetridiagonal(a1, a2, a3, b, n, ref d);
234
235	//
236	// Now problem is reduced to the cubic Hermite spline
237	//
238	buildhermitespline(x, y, d, n, ref c);
239	}
240
241
242	/*************************************************************************
243	This subroutine builds cubic Hermite spline coefficients table.
244
245	Input parameters:
246	X - spline nodes, array[0..N-1]
247	Y - function values, array[0..N-1]
248	D - derivatives, array[0..N-1]
249	N - points count, N>=2
250
251	Output parameters:
252	C - coefficients table. Used by SplineInterpolation and
253	other subroutines from this file.
254
255	-- ALGLIB PROJECT --
256	Copyright 23.06.2007 by Bochkanov Sergey
257	*************************************************************************/
258	public static void buildhermitespline(double[] x,
259	double[] y,
260	double[] d,
261	int n,
262	ref double[] c)
263	{
264	int i = 0;
265	int tblsize = 0;
266	double delta = 0;
267	double delta2 = 0;
268	double delta3 = 0;
269
270	x = (double[])x.Clone();
271	y = (double[])y.Clone();
272	d = (double[])d.Clone();
273
274	System.Diagnostics.Debug.Assert(n>=2, "BuildHermiteSpline: N<2!");
275
276	//
277	// Sort points
278	//
279	heapsortdpoints(ref x, ref y, ref d, n);
280
281	//
282	// Fill C:
283	// C[0] - length(C)
284	// C[1] - type(C):
285	// 3 - general cubic spline
286	// C[2] - N
287	// C[3]...C[3+N-1] - x[i], i = 0...N-1
288	// C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
289	//
290	tblsize = 3+n+(n-1)*4;
291	c = new double[tblsize-1+1];
292	c[0] = tblsize;
293	c[1] = 3;
294	c[2] = n;
295	for(i=0; i<=n-1; i++)
296	{
297	c[3+i] = x[i];
298	}
299	for(i=0; i<=n-2; i++)
300	{
301	delta = x[i+1]-x[i];
302	delta2 = AP.Math.Sqr(delta);
303	delta3 = delta*delta2;
304	c[3+n+4*i+0] = y[i];
305	c[3+n+4*i+1] = d[i];
306	c[3+n+4i+2] = (3(y[i+1]-y[i])-2d[i]delta-d[i+1]*delta)/delta2;
307	c[3+n+4i+3] = (2(y[i]-y[i+1])+d[i]delta+d[i+1]delta)/delta3;
308	}
309	}
310
311
312	/*************************************************************************
313	This subroutine builds Akima spline coefficients table.
314
315	Input parameters:
316	X - spline nodes, array[0..N-1]
317	Y - function values, array[0..N-1]
318	N - points count, N>=5
319
320	Output parameters:
321	C - coefficients table. Used by SplineInterpolation and
322	other subroutines from this file.
323
324	-- ALGLIB PROJECT --
325	Copyright 24.06.2007 by Bochkanov Sergey
326	*************************************************************************/
327	public static void buildakimaspline(double[] x,
328	double[] y,
329	int n,
330	ref double[] c)
331	{
332	int i = 0;
333	double[] d = new double[0];
334	double[] w = new double[0];
335	double[] diff = new double[0];
336
337	x = (double[])x.Clone();
338	y = (double[])y.Clone();
339
340	System.Diagnostics.Debug.Assert(n>=5, "BuildAkimaSpline: N<5!");
341
342	//
343	// Sort points
344	//
345	heapsortpoints(ref x, ref y, n);
346
347	//
348	// Prepare W (weights), Diff (divided differences)
349	//
350	w = new double[n-2+1];
351	diff = new double[n-2+1];
352	for(i=0; i<=n-2; i++)
353	{
354	diff[i] = (y[i+1]-y[i])/(x[i+1]-x[i]);
355	}
356	for(i=1; i<=n-2; i++)
357	{
358	w[i] = Math.Abs(diff[i]-diff[i-1]);
359	}
360
361	//
362	// Prepare Hermite interpolation scheme
363	//
364	d = new double[n-1+1];
365	for(i=2; i<=n-3; i++)
366	{
367	if( Math.Abs(w[i-1])+Math.Abs(w[i+1])!=0 )
368	{
369	d[i] = (w[i+1]diff[i-1]+w[i-1]diff[i])/(w[i+1]+w[i-1]);
370	}
371	else
372	{
373	d[i] = ((x[i+1]-x[i])diff[i-1]+(x[i]-x[i-1])diff[i])/(x[i+1]-x[i-1]);
374	}
375	}
376	d[0] = diffthreepoint(x[0], x[0], y[0], x[1], y[1], x[2], y[2]);
377	d[1] = diffthreepoint(x[1], x[0], y[0], x[1], y[1], x[2], y[2]);
378	d[n-2] = diffthreepoint(x[n-2], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
379	d[n-1] = diffthreepoint(x[n-1], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
380
381	//
382	// Build Akima spline using Hermite interpolation scheme
383	//
384	buildhermitespline(x, y, d, n, ref c);
385	}
386
387
388	/*************************************************************************
389	This subroutine calculates the value of the spline at the given point X.
390
391	Input parameters:
392	C - coefficients table. Built by BuildLinearSpline,
393	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
394	X - point
395
396	Result:
397	S(x)
398
399	-- ALGLIB PROJECT --
400	Copyright 23.06.2007 by Bochkanov Sergey
401	*************************************************************************/
402	public static double splineinterpolation(ref double[] c,
403	double x)
404	{
405	double result = 0;
406	int n = 0;
407	int l = 0;
408	int r = 0;
409	int m = 0;
410
411	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
412	n = (int)Math.Round(c[2]);
413
414	//
415	// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
416	//
417	l = 3;
418	r = 3+n-2+1;
419	while( l!=r-1 )
420	{
421	m = (l+r)/2;
422	if( c[m]>=x )
423	{
424	r = m;
425	}
426	else
427	{
428	l = m;
429	}
430	}
431
432	//
433	// Interpolation
434	//
435	x = x-c[l];
436	m = 3+n+4*(l-3);
437	result = c[m]+x(c[m+1]+x(c[m+2]+x*c[m+3]));
438	return result;
439	}
440
441
442	/*************************************************************************
443	This subroutine differentiates the spline.
444
445	Input parameters:
446	C - coefficients table. Built by BuildLinearSpline,
447	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
448	X - point
449
450	Result:
451	S - S(x)
452	DS - S'(x)
453	D2S - S''(x)
454
455	-- ALGLIB PROJECT --
456	Copyright 24.06.2007 by Bochkanov Sergey
457	*************************************************************************/
458	public static void splinedifferentiation(ref double[] c,
459	double x,
460	ref double s,
461	ref double ds,
462	ref double d2s)
463	{
464	int n = 0;
465	int l = 0;
466	int r = 0;
467	int m = 0;
468
469	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
470	n = (int)Math.Round(c[2]);
471
472	//
473	// Binary search
474	//
475	l = 3;
476	r = 3+n-2+1;
477	while( l!=r-1 )
478	{
479	m = (l+r)/2;
480	if( c[m]>=x )
481	{
482	r = m;
483	}
484	else
485	{
486	l = m;
487	}
488	}
489
490	//
491	// Differentiation
492	//
493	x = x-c[l];
494	m = 3+n+4*(l-3);
495	s = c[m]+x(c[m+1]+x(c[m+2]+x*c[m+3]));
496	ds = c[m+1]+2xc[m+2]+3AP.Math.Sqr(x)c[m+3];
497	d2s = 2c[m+2]+6x*c[m+3];
498	}
499
500
501	/*************************************************************************
502	This subroutine makes the copy of the spline.
503
504	Input parameters:
505	C - coefficients table. Built by BuildLinearSpline,
506	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
507
508	Result:
509	CC - spline copy
510
511	-- ALGLIB PROJECT --
512	Copyright 29.06.2007 by Bochkanov Sergey
513	*************************************************************************/
514	public static void splinecopy(ref double[] c,
515	ref double[] cc)
516	{
517	int s = 0;
518	int i_ = 0;
519
520	s = (int)Math.Round(c[0]);
521	cc = new double[s-1+1];
522	for(i_=0; i_<=s-1;i_++)
523	{
524	cc[i_] = c[i_];
525	}
526	}
527
528
529	/*************************************************************************
530	This subroutine unpacks the spline into the coefficients table.
531
532	Input parameters:
533	C - coefficients table. Built by BuildLinearSpline,
534	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
535	X - point
536
537	Result:
538	Tbl - coefficients table, unpacked format, array[0..N-2, 0..5].
539	For I = 0...N-2:
540	Tbl[I,0] = X[i]
541	Tbl[I,1] = X[i+1]
542	Tbl[I,2] = C0
543	Tbl[I,3] = C1
544	Tbl[I,4] = C2
545	Tbl[I,5] = C3
546	On [x[i], x[i+1]] spline is equals to:
547	S(x) = C0 + C1t + C2t^2 + C3*t^3
548	t = x-x[i]
549
550	-- ALGLIB PROJECT --
551	Copyright 29.06.2007 by Bochkanov Sergey
552	*************************************************************************/
553	public static void splineunpack(ref double[] c,
554	ref int n,
555	ref double[,] tbl)
556	{
557	int i = 0;
558
559	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineUnpack: incorrect C!");
560	n = (int)Math.Round(c[2]);
561	tbl = new double[n-2+1, 5+1];
562
563	//
564	// Fill
565	//
566	for(i=0; i<=n-2; i++)
567	{
568	tbl[i,0] = c[3+i];
569	tbl[i,1] = c[3+i+1];
570	tbl[i,2] = c[3+n+4*i];
571	tbl[i,3] = c[3+n+4*i+1];
572	tbl[i,4] = c[3+n+4*i+2];
573	tbl[i,5] = c[3+n+4*i+3];
574	}
575	}
576
577
578	/*************************************************************************
579	This subroutine performs linear transformation of the spline argument.
580
581	Input parameters:
582	C - coefficients table. Built by BuildLinearSpline,
583	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
584	A, B- transformation coefficients: x = A*t + B
585	Result:
586	C - transformed spline
587
588	-- ALGLIB PROJECT --
589	Copyright 30.06.2007 by Bochkanov Sergey
590	*************************************************************************/
591	public static void splinelintransx(ref double[] c,
592	double a,
593	double b)
594	{
595	int i = 0;
596	int n = 0;
597	double v = 0;
598	double dv = 0;
599	double d2v = 0;
600	double[] x = new double[0];
601	double[] y = new double[0];
602	double[] d = new double[0];
603
604	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
605	n = (int)Math.Round(c[2]);
606
607	//
608	// Special case: A=0
609	//
610	if( a==0 )
611	{
612	v = splineinterpolation(ref c, b);
613	for(i=0; i<=n-2; i++)
614	{
615	c[3+n+4*i] = v;
616	c[3+n+4*i+1] = 0;
617	c[3+n+4*i+2] = 0;
618	c[3+n+4*i+3] = 0;
619	}
620	return;
621	}
622
623	//
624	// General case: A<>0.
625	// Unpack, X, Y, dY/dX.
626	// Scale and pack again.
627	//
628	x = new double[n-1+1];
629	y = new double[n-1+1];
630	d = new double[n-1+1];
631	for(i=0; i<=n-1; i++)
632	{
633	x[i] = c[3+i];
634	splinedifferentiation(ref c, x[i], ref v, ref dv, ref d2v);
635	x[i] = (x[i]-b)/a;
636	y[i] = v;
637	d[i] = a*dv;
638	}
639	buildhermitespline(x, y, d, n, ref c);
640	}
641
642
643	/*************************************************************************
644	This subroutine performs linear transformation of the spline.
645
646	Input parameters:
647	C - coefficients table. Built by BuildLinearSpline,
648	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
649	A, B- transformation coefficients: S2(x) = A*S(x) + B
650	Result:
651	C - transformed spline
652
653	-- ALGLIB PROJECT --
654	Copyright 30.06.2007 by Bochkanov Sergey
655	*************************************************************************/
656	public static void splinelintransy(ref double[] c,
657	double a,
658	double b)
659	{
660	int i = 0;
661	int n = 0;
662	double v = 0;
663	double dv = 0;
664	double d2v = 0;
665	double[] x = new double[0];
666	double[] y = new double[0];
667	double[] d = new double[0];
668
669	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
670	n = (int)Math.Round(c[2]);
671
672	//
673	// Special case: A=0
674	//
675	for(i=0; i<=n-2; i++)
676	{
677	c[3+n+4i] = ac[3+n+4*i]+b;
678	c[3+n+4i+1] = ac[3+n+4*i+1];
679	c[3+n+4i+2] = ac[3+n+4*i+2];
680	c[3+n+4i+3] = ac[3+n+4*i+3];
681	}
682	}
683
684
685	/*************************************************************************
686	This subroutine integrates the spline.
687
688	Input parameters:
689	C - coefficients table. Built by BuildLinearSpline,
690	BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
691	X - right bound of the integration interval [a, x]
692	Result:
693	integral(S(t)dt,a,x)
694
695	-- ALGLIB PROJECT --
696	Copyright 23.06.2007 by Bochkanov Sergey
697	*************************************************************************/
698	public static double splineintegration(ref double[] c,
699	double x)
700	{
701	double result = 0;
702	int n = 0;
703	int i = 0;
704	int l = 0;
705	int r = 0;
706	int m = 0;
707	double w = 0;
708
709	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineIntegration: incorrect C!");
710	n = (int)Math.Round(c[2]);
711
712	//
713	// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
714	//
715	l = 3;
716	r = 3+n-2+1;
717	while( l!=r-1 )
718	{
719	m = (l+r)/2;
720	if( c[m]>=x )
721	{
722	r = m;
723	}
724	else
725	{
726	l = m;
727	}
728	}
729
730	//
731	// Integration
732	//
733	result = 0;
734	for(i=3; i<=l-1; i++)
735	{
736	w = c[i+1]-c[i];
737	m = 3+n+4*(i-3);
738	result = result+c[m]*w;
739	result = result+c[m+1]*AP.Math.Sqr(w)/2;
740	result = result+c[m+2]AP.Math.Sqr(w)w/3;
741	result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
742	}
743	w = x-c[l];
744	m = 3+n+4*(l-3);
745	result = result+c[m]*w;
746	result = result+c[m+1]*AP.Math.Sqr(w)/2;
747	result = result+c[m+2]AP.Math.Sqr(w)w/3;
748	result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
749	return result;
750	}
751
752
753	/*************************************************************************
754	Obsolete subroutine, left for backward compatibility.
755	*************************************************************************/
756	public static void spline3buildtable(int n,
757	int diffn,
758	double[] x,
759	double[] y,
760	double boundl,
761	double boundr,
762	ref double[,] ctbl)
763	{
764	bool c = new bool();
765	int e = 0;
766	int g = 0;
767	double tmp = 0;
768	int nxm1 = 0;
769	int i = 0;
770	int j = 0;
771	double dx = 0;
772	double dxj = 0;
773	double dyj = 0;
774	double dxjp1 = 0;
775	double dyjp1 = 0;
776	double dxp = 0;
777	double dyp = 0;
778	double yppa = 0;
779	double yppb = 0;
780	double pj = 0;
781	double b1 = 0;
782	double b2 = 0;
783	double b3 = 0;
784	double b4 = 0;
785
786	x = (double[])x.Clone();
787	y = (double[])y.Clone();
788
789	n = n-1;
790	g = (n+1)/2;
791	do
792	{
793	i = g;
794	do
795	{
796	j = i-g;
797	c = true;
798	do
799	{
800	if( x[j]<=x[j+g] )
801	{
802	c = false;
803	}
804	else
805	{
806	tmp = x[j];
807	x[j] = x[j+g];
808	x[j+g] = tmp;
809	tmp = y[j];
810	y[j] = y[j+g];
811	y[j+g] = tmp;
812	}
813	j = j-1;
814	}
815	while( j>=0 & c );
816	i = i+1;
817	}
818	while( i<=n );
819	g = g/2;
820	}
821	while( g>0 );
822	ctbl = new double[4+1, n+1];
823	n = n+1;
824	if( diffn==1 )
825	{
826	b1 = 1;
827	b2 = 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-boundl);
828	b3 = 1;
829	b4 = 6/(x[n-1]-x[n-2])*(boundr-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
830	}
831	else
832	{
833	b1 = 0;
834	b2 = 2*boundl;
835	b3 = 0;
836	b4 = 2*boundr;
837	}
838	nxm1 = n-1;
839	if( n>=2 )
840	{
841	if( n>2 )
842	{
843	dxj = x[1]-x[0];
844	dyj = y[1]-y[0];
845	j = 2;
846	while( j<=nxm1 )
847	{
848	dxjp1 = x[j]-x[j-1];
849	dyjp1 = y[j]-y[j-1];
850	dxp = dxj+dxjp1;
851	ctbl[1,j-1] = dxjp1/dxp;
852	ctbl[2,j-1] = 1-ctbl[1,j-1];
853	ctbl[3,j-1] = 6*(dyjp1/dxjp1-dyj/dxj)/dxp;
854	dxj = dxjp1;
855	dyj = dyjp1;
856	j = j+1;
857	}
858	}
859	ctbl[1,0] = -(b1/2);
860	ctbl[2,0] = b2/2;
861	if( n!=2 )
862	{
863	j = 2;
864	while( j<=nxm1 )
865	{
866	pj = ctbl[2,j-1]*ctbl[1,j-2]+2;
867	ctbl[1,j-1] = -(ctbl[1,j-1]/pj);
868	ctbl[2,j-1] = (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,j-2])/pj;
869	j = j+1;
870	}
871	}
872	yppb = (b4-b3ctbl[2,nxm1-1])/(b3ctbl[1,nxm1-1]+2);
873	i = 1;
874	while( i<=nxm1 )
875	{
876	j = n-i;
877	yppa = ctbl[1,j-1]*yppb+ctbl[2,j-1];
878	dx = x[j]-x[j-1];
879	ctbl[3,j-1] = (yppb-yppa)/dx/6;
880	ctbl[2,j-1] = yppa/2;
881	ctbl[1,j-1] = (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]dx)dx;
882	yppb = yppa;
883	i = i+1;
884	}
885	for(i=1; i<=n; i++)
886	{
887	ctbl[0,i-1] = y[i-1];
888	ctbl[4,i-1] = x[i-1];
889	}
890	}
891	}
892
893
894	/*************************************************************************
895	Obsolete subroutine, left for backward compatibility.
896	*************************************************************************/
897	public static double spline3interpolate(int n,
898	ref double[,] c,
899	double x)
900	{
901	double result = 0;
902	int i = 0;
903	int l = 0;
904	int half = 0;
905	int first = 0;
906	int middle = 0;
907
908	n = n-1;
909	l = n;
910	first = 0;
911	while( l>0 )
912	{
913	half = l/2;
914	middle = first+half;
915	if( c[4,middle]<x )
916	{
917	first = middle+1;
918	l = l-half-1;
919	}
920	else
921	{
922	l = half;
923	}
924	}
925	i = first-1;
926	if( i<0 )
927	{
928	i = 0;
929	}
930	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])));
931	return result;
932	}
933
934
935	/*************************************************************************
936	Internal subroutine. Heap sort.
937	*************************************************************************/
938	private static void heapsortpoints(ref double[] x,
939	ref double[] y,
940	int n)
941	{
942	int i = 0;
943	int j = 0;
944	int k = 0;
945	int t = 0;
946	double tmp = 0;
947	bool isascending = new bool();
948	bool isdescending = new bool();
949
950
951	//
952	// Test for already sorted set
953	//
954	isascending = true;
955	isdescending = true;
956	for(i=1; i<=n-1; i++)
957	{
958	isascending = isascending & x[i]>x[i-1];
959	isdescending = isdescending & x[i]<x[i-1];
960	}
961	if( isascending )
962	{
963	return;
964	}
965	if( isdescending )
966	{
967	for(i=0; i<=n-1; i++)
968	{
969	j = n-1-i;
970	if( j<=i )
971	{
972	break;
973	}
974	tmp = x[i];
975	x[i] = x[j];
976	x[j] = tmp;
977	tmp = y[i];
978	y[i] = y[j];
979	y[j] = tmp;
980	}
981	return;
982	}
983
984	//
985	// Special case: N=1
986	//
987	if( n==1 )
988	{
989	return;
990	}
991
992	//
993	// General case
994	//
995	i = 2;
996	do
997	{
998	t = i;
999	while( t!=1 )
1000	{
1001	k = t/2;
1002	if( x[k-1]>=x[t-1] )
1003	{
1004	t = 1;
1005	}
1006	else
1007	{
1008	tmp = x[k-1];
1009	x[k-1] = x[t-1];
1010	x[t-1] = tmp;
1011	tmp = y[k-1];
1012	y[k-1] = y[t-1];
1013	y[t-1] = tmp;
1014	t = k;
1015	}
1016	}
1017	i = i+1;
1018	}
1019	while( i<=n );
1020	i = n-1;
1021	do
1022	{
1023	tmp = x[i];
1024	x[i] = x[0];
1025	x[0] = tmp;
1026	tmp = y[i];
1027	y[i] = y[0];
1028	y[0] = tmp;
1029	t = 1;
1030	while( t!=0 )
1031	{
1032	k = 2*t;
1033	if( k>i )
1034	{
1035	t = 0;
1036	}
1037	else
1038	{
1039	if( k<i )
1040	{
1041	if( x[k]>x[k-1] )
1042	{
1043	k = k+1;
1044	}
1045	}
1046	if( x[t-1]>=x[k-1] )
1047	{
1048	t = 0;
1049	}
1050	else
1051	{
1052	tmp = x[k-1];
1053	x[k-1] = x[t-1];
1054	x[t-1] = tmp;
1055	tmp = y[k-1];
1056	y[k-1] = y[t-1];
1057	y[t-1] = tmp;
1058	t = k;
1059	}
1060	}
1061	}
1062	i = i-1;
1063	}
1064	while( i>=1 );
1065	}
1066
1067
1068	/*************************************************************************
1069	Internal subroutine. Heap sort.
1070	*************************************************************************/
1071	private static void heapsortdpoints(ref double[] x,
1072	ref double[] y,
1073	ref double[] d,
1074	int n)
1075	{
1076	int i = 0;
1077	int j = 0;
1078	int k = 0;
1079	int t = 0;
1080	double tmp = 0;
1081	bool isascending = new bool();
1082	bool isdescending = new bool();
1083
1084
1085	//
1086	// Test for already sorted set
1087	//
1088	isascending = true;
1089	isdescending = true;
1090	for(i=1; i<=n-1; i++)
1091	{
1092	isascending = isascending & x[i]>x[i-1];
1093	isdescending = isdescending & x[i]<x[i-1];
1094	}
1095	if( isascending )
1096	{
1097	return;
1098	}
1099	if( isdescending )
1100	{
1101	for(i=0; i<=n-1; i++)
1102	{
1103	j = n-1-i;
1104	if( j<=i )
1105	{
1106	break;
1107	}
1108	tmp = x[i];
1109	x[i] = x[j];
1110	x[j] = tmp;
1111	tmp = y[i];
1112	y[i] = y[j];
1113	y[j] = tmp;
1114	tmp = d[i];
1115	d[i] = d[j];
1116	d[j] = tmp;
1117	}
1118	return;
1119	}
1120
1121	//
1122	// Special case: N=1
1123	//
1124	if( n==1 )
1125	{
1126	return;
1127	}
1128
1129	//
1130	// General case
1131	//
1132	i = 2;
1133	do
1134	{
1135	t = i;
1136	while( t!=1 )
1137	{
1138	k = t/2;
1139	if( x[k-1]>=x[t-1] )
1140	{
1141	t = 1;
1142	}
1143	else
1144	{
1145	tmp = x[k-1];
1146	x[k-1] = x[t-1];
1147	x[t-1] = tmp;
1148	tmp = y[k-1];
1149	y[k-1] = y[t-1];
1150	y[t-1] = tmp;
1151	tmp = d[k-1];
1152	d[k-1] = d[t-1];
1153	d[t-1] = tmp;
1154	t = k;
1155	}
1156	}
1157	i = i+1;
1158	}
1159	while( i<=n );
1160	i = n-1;
1161	do
1162	{
1163	tmp = x[i];
1164	x[i] = x[0];
1165	x[0] = tmp;
1166	tmp = y[i];
1167	y[i] = y[0];
1168	y[0] = tmp;
1169	tmp = d[i];
1170	d[i] = d[0];
1171	d[0] = tmp;
1172	t = 1;
1173	while( t!=0 )
1174	{
1175	k = 2*t;
1176	if( k>i )
1177	{
1178	t = 0;
1179	}
1180	else
1181	{
1182	if( k<i )
1183	{
1184	if( x[k]>x[k-1] )
1185	{
1186	k = k+1;
1187	}
1188	}
1189	if( x[t-1]>=x[k-1] )
1190	{
1191	t = 0;
1192	}
1193	else
1194	{
1195	tmp = x[k-1];
1196	x[k-1] = x[t-1];
1197	x[t-1] = tmp;
1198	tmp = y[k-1];
1199	y[k-1] = y[t-1];
1200	y[t-1] = tmp;
1201	tmp = d[k-1];
1202	d[k-1] = d[t-1];
1203	d[t-1] = tmp;
1204	t = k;
1205	}
1206	}
1207	}
1208	i = i-1;
1209	}
1210	while( i>=1 );
1211	}
1212
1213
1214	/*************************************************************************
1215	Internal subroutine. Tridiagonal solver.
1216	*************************************************************************/
1217	private static void solvetridiagonal(double[] a,
1218	double[] b,
1219	double[] c,
1220	double[] d,
1221	int n,
1222	ref double[] x)
1223	{
1224	int k = 0;
1225	double t = 0;
1226
1227	a = (double[])a.Clone();
1228	b = (double[])b.Clone();
1229	c = (double[])c.Clone();
1230	d = (double[])d.Clone();
1231
1232	x = new double[n-1+1];
1233	a[0] = 0;
1234	c[n-1] = 0;
1235	for(k=1; k<=n-1; k++)
1236	{
1237	t = a[k]/b[k-1];
1238	b[k] = b[k]-t*c[k-1];
1239	d[k] = d[k]-t*d[k-1];
1240	}
1241	x[n-1] = d[n-1]/b[n-1];
1242	for(k=n-2; k>=0; k--)
1243	{
1244	x[k] = (d[k]-c[k]*x[k+1])/b[k];
1245	}
1246	}
1247
1248
1249	/*************************************************************************
1250	Internal subroutine. Three-point differentiation
1251	*************************************************************************/
1252	private static double diffthreepoint(double t,
1253	double x0,
1254	double f0,
1255	double x1,
1256	double f1,
1257	double x2,
1258	double f2)
1259	{
1260	double result = 0;
1261	double a = 0;
1262	double b = 0;
1263
1264	t = t-x0;
1265	x1 = x1-x0;
1266	x2 = x2-x0;
1267	a = (f2-f0-x2/x1(f1-f0))/(AP.Math.Sqr(x2)-x1x2);
1268	b = (f1-f0-a*AP.Math.Sqr(x1))/x1;
1269	result = 2at+b;
1270	return result;
1271	}
1272	}
1273	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences