Menu
▾
▴
[Hugin-cvs] SF.net SVN: hugin:[4327] hugin/trunk/src/lens_calibrate
|
From: <bl...@us...> - 2009-09-05 11:14:14
|
Revision: 4327
http://hugin.svn.sourceforge.net/hugin/?rev=4327&view=rev
Author: blimbo
Date: 2009-09-05 11:14:06 +0000 (Sat, 05 Sep 2009)
Log Message:
-----------
Removed old files
Modified Paths:
--------------
hugin/trunk/src/lens_calibrate/CMakeLists.txt
Removed Paths:
-------------
hugin/trunk/src/lens_calibrate/Spline.cpp
hugin/trunk/src/lens_calibrate/Spline.h
Modified: hugin/trunk/src/lens_calibrate/CMakeLists.txt
===================================================================
--- hugin/trunk/src/lens_calibrate/CMakeLists.txt 2009-09-05 10:45:09 UTC (rev 4326)
+++ hugin/trunk/src/lens_calibrate/CMakeLists.txt 2009-09-05 11:14:06 UTC (rev 4327)
@@ -1,23 +1,3 @@
-/***************************************************************************
- * Copyright (C) 2009 by Tim Nugent *
- * tim...@gm... *
- * *
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation; either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program; if not, write to the *
- * Free Software Foundation, Inc., *
- * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
- ***************************************************************************/
-
ADD_EXECUTABLE(calibrate_lens
Main.cpp
Globals.cpp
Deleted: hugin/trunk/src/lens_calibrate/Spline.cpp
===================================================================
--- hugin/trunk/src/lens_calibrate/Spline.cpp 2009-09-05 10:45:09 UTC (rev 4326)
+++ hugin/trunk/src/lens_calibrate/Spline.cpp 2009-09-05 11:14:06 UTC (rev 4327)
@@ -1,6646 +0,0 @@
-/***************************************************************************
- * Copyright (C) John Burkardt *
- * The computer code and data files described are distributed under the *
- * GNU LGPL license. *
- ***************************************************************************/
-
-# include <cstdlib>
-# include <iostream>
-# include <iomanip>
-# include <cmath>
-# include <ctime>
-#include <string.h>
-
-using namespace std;
-
-# include "Spline.h"
-
-//****************************************************************************80
-
-double basis_function_b_val ( double tdata[], double tval )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_FUNCTION_B_VAL evaluates the B spline basis function.
-//
-// Discussion:
-//
-// The B spline basis function is a piecewise cubic which
-// has the properties that:
-//
-// * it equals 2/3 at TDATA(3), 1/6 at TDATA(2) and TDATA(4);
-// * it is 0 for TVAL <= TDATA(1) or TDATA(5) <= TVAL;
-// * it is strictly increasing from TDATA(1) to TDATA(3),
-// and strictly decreasing from TDATA(3) to TDATA(5);
-// * the function and its first two derivatives are continuous
-// at each node TDATA(I).
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 24 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// Alan Davies, Philip Samuels,
-// An Introduction to Computational Geometry for Curves and Surfaces,
-// Clarendon Press, 1996,
-// ISBN: 0-19-851478-6,
-// LC: QA448.D38.
-//
-// Parameters:
-//
-// Input, double TDATA(5), the nodes associated with the basis function.
-// The entries of TDATA are assumed to be distinct and increasing.
-//
-// Input, double TVAL, a point at which the B spline basis function is
-// to be evaluated.
-//
-// Output, double BASIS_FUNCTION_B_VAL, the value of the function at TVAL.
-//
-{
-# define NDATA 5
-
- int left;
- int right;
- double u;
- double yval;
-
- if ( tval <= tdata[0] || tdata[NDATA-1] <= tval )
- {
- yval = 0.0;
- return yval;
- }
-//
-// Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL.
-//
- r8vec_bracket ( NDATA, tdata, tval, &left, &right );
-//
-// U is the normalized coordinate of TVAL in this interval.
-//
- u = ( tval - tdata[left-1] ) / ( tdata[right-1] - tdata[left-1] );
-//
-// Now evaluate the function.
-//
- if ( tval < tdata[1] )
- {
- yval = pow ( u, 3 ) / 6.0;
- }
- else if ( tval < tdata[2] )
- {
- yval = ( ( ( - 3.0
- * u + 3.0 )
- * u + 3.0 )
- * u + 1.0 ) / 6.0;
- }
- else if ( tval < tdata[3])
- {
- yval = ( ( ( + 3.0
- * u - 6.0 )
- * u + 0.0 )
- * u + 4.0 ) / 6.0;
- }
- else if ( tval < tdata[4] )
- {
- yval = pow ( ( 1.0 - u ), 3 ) / 6.0;
- }
-
- return yval;
-
-# undef NDATA
-}
-//****************************************************************************80
-
-double basis_function_beta_val ( double beta1, double beta2, double tdata[],
- double tval )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_FUNCTION_BETA_VAL evaluates the beta spline basis function.
-//
-// Discussion:
-//
-// With BETA1 = 1 and BETA2 = 0, the beta spline basis function
-// equals the B spline basis function.
-//
-// With BETA1 large, and BETA2 = 0, the beta spline basis function
-// skews to the right, that is, its maximum increases, and occurs
-// to the right of the center.
-//
-// With BETA1 = 1 and BETA2 large, the beta spline becomes more like
-// a linear basis function; that is, its value in the outer two intervals
-// goes to zero, and its behavior in the inner two intervals approaches
-// a piecewise linear function that is 0 at the second node, 1 at the
-// third, and 0 at the fourth.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 24 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// Alan Davies, Philip Samuels,
-// An Introduction to Computational Geometry for Curves and Surfaces,
-// Clarendon Press, 1996,
-// ISBN: 0-19-851478-6,
-// LC: QA448.D38.
-//
-// Parameters:
-//
-// Input, double BETA1, the skew or bias parameter.
-// BETA1 = 1 for no skew or bias.
-//
-// Input, double BETA2, the tension parameter.
-// BETA2 = 0 for no tension.
-//
-// Input, double TDATA[5], the nodes associated with the basis function.
-// The entries of TDATA are assumed to be distinct and increasing.
-//
-// Input, double TVAL, a point at which the B spline basis function is
-// to be evaluated.
-//
-// Output, double BASIS_FUNCTION_BETA_VAL, the value of the function at TVAL.
-//
-{
-# define NDATA 5
-
- double a;
- double b;
- double c;
- double d;
- int left;
- int right;
- double u;
- double yval;
-
- if ( tval <= tdata[0] || tdata[NDATA-1] <= tval )
- {
- yval = 0.0;
- return yval;
- }
-//
-// Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL.
-//
- r8vec_bracket ( NDATA, tdata, tval, &left, &right );
-//
-// U is the normalized coordinate of TVAL in this interval.
-//
- u = ( tval - tdata[left-1] ) / ( tdata[right-1] - tdata[left-1] );
-//
-// Now evaluate the function.
-//
- if ( tval < tdata[1] )
- {
- yval = 2.0 * u * u * u;
- }
- else if ( tval < tdata[2] )
- {
- a = beta2 + 4.0 * beta1 + 4.0 * beta1 * beta1
- + 6.0 * ( 1.0 - beta1 * beta1 )
- - 3.0 * ( 2.0 + beta2 + 2.0 * beta1 )
- + 2.0 * ( 1.0 + beta2 + beta1 + beta1 * beta1 );
-
- b = - 6.0 * ( 1.0 - beta1 * beta1 )
- + 6.0 * ( 2.0 + beta2 + 2.0 * beta1 )
- - 6.0 * ( 1.0 + beta2 + beta1 + beta1 * beta1 );
-
- c = - 3.0 * ( 2.0 + beta2 + 2.0 * beta1 )
- + 6.0 * ( 1.0 + beta2 + beta1 + beta1 * beta1 );
-
- d = - 2.0 * ( 1.0 + beta2 + beta1 + beta1 * beta1 );
-
- yval = a + b * u + c * u * u + d * u * u * u;
- }
- else if ( tval < tdata[3] )
- {
- a = beta2 + 4.0 * beta1 + 4.0 * beta1 * beta1;
-
- b = - 6.0 * beta1 * ( 1.0 - beta1 * beta1 );
-
- c = - 3.0 * ( beta2 + 2.0 * beta1 * beta1
- + 2.0 * beta1 * beta1 * beta1 );
-
- d = 2.0 * ( beta2 + beta1 + beta1 * beta1 + beta1 * beta1 * beta1 );
-
- yval = a + b * u + c * u * u + d * u * u * u;
- }
- else if ( tval < tdata[4] )
- {
- yval = 2.0 * pow ( beta1 * ( 1.0 - u ), 3 );
- }
-
- yval = yval / ( 2.0 + beta2 + 4.0 * beta1 + 4.0 * beta1 * beta1
- + 2.0 * beta1 * beta1 * beta1 );
-
- return yval;
-# undef NDATA
-}
-//****************************************************************************80
-
-double *basis_matrix_b_uni ( void )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_B_UNI sets up the uniform B spline basis matrix.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 11 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// James Foley, Andries vanDam, Steven Feiner, John Hughes,
-// Computer Graphics, Principles and Practice,
-// Second Edition,
-// Addison Wesley, 1995,
-// ISBN: 0201848406,
-// LC: T385.C5735.
-//
-// Parameters:
-//
-// Output, double BASIS_MATRIX_B_UNI[4*4], the basis matrix.
-//
-{
- int i;
- int j;
- double *mbasis;
- double mbasis_save[4*4] = {
- -1.0 / 6.0,
- 3.0 / 6.0,
- -3.0 / 6.0,
- 1.0 / 6.0,
- 3.0 / 6.0,
- -6.0 / 6.0,
- 0.0,
- 4.0 / 6.0,
- -3.0 / 6.0,
- 3.0 / 6.0,
- 3.0 / 6.0,
- 1.0 / 6.0,
- 1.0 / 6.0,
- 0.0,
- 0.0,
- 0.0 };
-
- mbasis = new double[4*4];
-
- for ( j = 0; j < 4; j++ )
- {
- for ( i = 0; i < 4; i++ )
- {
- mbasis[i+j*4] = mbasis_save[i+j*4];
- }
- }
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_beta_uni ( double beta1, double beta2 )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_BETA_UNI sets up the uniform beta spline basis matrix.
-//
-// Discussion:
-//
-// If BETA1 = 1 and BETA2 = 0, then the beta spline reduces to
-// the B spline.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// James Foley, Andries vanDam, Steven Feiner, John Hughes,
-// Computer Graphics, Principles and Practice,
-// Second Edition,
-// Addison Wesley, 1995,
-// ISBN: 0201848406,
-// LC: T385.C5735.
-//
-// Parameters:
-//
-// Input, double BETA1, the skew or bias parameter.
-// BETA1 = 1 for no skew or bias.
-//
-// Input, double BETA2, the tension parameter.
-// BETA2 = 0 for no tension.
-//
-// Output, double BASIS_MATRIX_BETA_UNI[4*4], the basis matrix.
-//
-{
- double delta;
- int i;
- int j;
- double *mbasis;
-
- mbasis = new double[4*4];
-
- mbasis[0+0*4] = - 2.0 * beta1 * beta1 * beta1;
- mbasis[0+1*4] = 2.0 * beta2
- + 2.0 * beta1 * ( beta1 * beta1 + beta1 + 1.0 );
- mbasis[0+2*4] = - 2.0 * ( beta2 + beta1 * beta1 + beta1 + 1.0 );
- mbasis[0+3*4] = 2.0;
-
- mbasis[1+0*4] = 6.0 * beta1 * beta1 * beta1;
- mbasis[1+1*4] = - 3.0 * beta2
- - 6.0 * beta1 * beta1 * ( beta1 + 1.0 );
- mbasis[1+2*4] = 3.0 * beta2 + 6.0 * beta1 * beta1;
- mbasis[1+3*4] = 0.0;
-
- mbasis[2+0*4] = - 6.0 * beta1 * beta1 * beta1;
- mbasis[2+1*4] = 6.0 * beta1 * ( beta1 - 1.0 ) * ( beta1 + 1.0 );
- mbasis[2+2*4] = 6.0 * beta1;
- mbasis[2+3*4] = 0.0;
-
- mbasis[3+0*4] = 2.0 * beta1 * beta1 * beta1;
- mbasis[3+1*4] = 4.0 * beta1 * ( beta1 + 1.0 ) + beta2;
- mbasis[3+2*4] = 2.0;
- mbasis[3+3*4] = 0.0;
-
- delta = ( ( 2.0
- * beta1 + 4.0 )
- * beta1 + 4.0 )
- * beta1 + 2.0 + beta2;
-
- for ( j = 0; j < 4; j++ )
- {
- for ( i = 0; i < 4; i++ )
- {
- mbasis[i+j*4] = mbasis[i+j*4] / delta;
- }
- }
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_bezier ( void )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_BEZIER_UNI sets up the cubic Bezier spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points are stored as
-// ( P1, P2, P3, P4 ). P1 is the function value at T = 0, while
-// P2 is used to approximate the derivative at T = 0 by
-// dP/dt = 3 * ( P2 - P1 ). Similarly, P4 is the function value
-// at T = 1, and P3 is used to approximate the derivative at T = 1
-// by dP/dT = 3 * ( P4 - P3 ).
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// James Foley, Andries vanDam, Steven Feiner, John Hughes,
-// Computer Graphics, Principles and Practice,
-// Second Edition,
-// Addison Wesley, 1995,
-// ISBN: 0201848406,
-// LC: T385.C5735.
-//
-// Parameters:
-//
-// Output, double BASIS_MATRIX_BEZIER[4*4], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[4*4];
-
- mbasis[0+0*4] = -1.0;
- mbasis[0+1*4] = 3.0;
- mbasis[0+2*4] = -3.0;
- mbasis[0+3*4] = 1.0;
-
- mbasis[1+0*4] = 3.0;
- mbasis[1+1*4] = -6.0;
- mbasis[1+2*4] = 3.0;
- mbasis[1+3*4] = 0.0;
-
- mbasis[2+0*4] = -3.0;
- mbasis[2+1*4] = 3.0;
- mbasis[2+2*4] = 0.0;
- mbasis[2+3*4] = 0.0;
-
- mbasis[3+0*4] = 1.0;
- mbasis[3+1*4] = 0.0;
- mbasis[3+2*4] = 0.0;
- mbasis[3+3*4] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_hermite ( void )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_HERMITE sets up the Hermite spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points are stored as
-// ( P1, P2, P1', P2' ), with P1 and P1' being the data value and
-// the derivative dP/dT at T = 0, while P2 and P2' apply at T = 1.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// James Foley, Andries vanDam, Steven Feiner, John Hughes,
-// Computer Graphics, Principles and Practice,
-// Second Edition,
-// Addison Wesley, 1995,
-// ISBN: 0201848406,
-// LC: T385.C5735.
-//
-// Parameters:
-//
-// Output, double BASIS_MATRIX_HERMITE[4*4], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[4*4];
-
- mbasis[0+0*4] = 2.0;
- mbasis[0+1*4] = -2.0;
- mbasis[0+2*4] = 1.0;
- mbasis[0+3*4] = 1.0;
-
- mbasis[1+0*4] = -3.0;
- mbasis[1+1*4] = 3.0;
- mbasis[1+2*4] = -2.0;
- mbasis[1+3*4] = -1.0;
-
- mbasis[2+0*4] = 0.0;
- mbasis[2+1*4] = 0.0;
- mbasis[2+2*4] = 1.0;
- mbasis[2+3*4] = 0.0;
-
- mbasis[3+0*4] = 1.0;
- mbasis[3+1*4] = 0.0;
- mbasis[3+2*4] = 0.0;
- mbasis[3+3*4] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_overhauser_nonuni ( double alpha, double beta )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_OVERHAUSER_NONUNI sets the nonuniform Overhauser spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points P1, P2, P3 and
-// P4 are not uniformly spaced in T, and that P2 corresponds to T = 0,
-// and P3 to T = 1.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, double ALPHA, BETA.
-// ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || )
-// BETA = || P3 - P2 || / ( || P4 - P3 || + || P3 - P2 || ).
-//
-// Output, double BASIS_MATRIX_OVERHAUSER_NONUNI[4*4], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[4*4];
-
- mbasis[0+0*4] = - ( 1.0 - alpha ) * ( 1.0 - alpha ) / alpha;
- mbasis[0+1*4] = beta + ( 1.0 - alpha ) / alpha;
- mbasis[0+2*4] = alpha - 1.0 / ( 1.0 - beta );
- mbasis[0+3*4] = beta * beta / ( 1.0 - beta );
-
- mbasis[1+0*4] = 2.0 * ( 1.0 - alpha ) * ( 1.0 - alpha ) / alpha;
- mbasis[1+1*4] = ( - 2.0 * ( 1.0 - alpha ) - alpha * beta ) / alpha;
- mbasis[1+2*4] = ( 2.0 * ( 1.0 - alpha )
- - beta * ( 1.0 - 2.0 * alpha ) ) / ( 1.0 - beta );
- mbasis[1+3*4] = - beta * beta / ( 1.0 - beta );
-
- mbasis[2+0*4] = - ( 1.0 - alpha ) * ( 1.0 - alpha ) / alpha;
- mbasis[2+1*4] = ( 1.0 - 2.0 * alpha ) / alpha;
- mbasis[2+2*4] = alpha;
- mbasis[2+3*4] = 0.0;
-
- mbasis[3+0*4] = 0.0;
- mbasis[3+1*4] = 1.0;
- mbasis[3+2*4] = 0.0;
- mbasis[3+3*4] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_overhauser_nul ( double alpha )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_OVERHAUSER_NUL sets the nonuniform left Overhauser spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points P1, P2, and
-// P3 are not uniformly spaced in T, and that P1 corresponds to T = 0,
-// and P2 to T = 1. (???)
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, double ALPHA.
-// ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || )
-//
-// Output, double BASIS_MATRIX_OVERHAUSER_NUL[3*3], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[3*3];
-
- mbasis[0+0*3] = 1.0 / alpha;
- mbasis[0+1*3] = - 1.0 / ( alpha * ( 1.0 - alpha ) );
- mbasis[0+2*3] = 1.0 / ( 1.0 - alpha );
-
- mbasis[1+0*3] = - ( 1.0 + alpha ) / alpha;
- mbasis[1+1*3] = 1.0 / ( alpha * ( 1.0 - alpha ) );
- mbasis[1+2*3] = - alpha / ( 1.0 - alpha );
-
- mbasis[2+0*3] = 1.0;
- mbasis[2+1*3] = 0.0;
- mbasis[2+2*3] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_overhauser_nur ( double beta )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_OVERHAUSER_NUR sets the nonuniform right Overhauser spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points PN-2, PN-1, and
-// PN are not uniformly spaced in T, and that PN-1 corresponds to T = 0,
-// and PN to T = 1. (???)
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 14 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, double BETA.
-// BETA = || P(N) - P(N-1) || / ( || P(N) - P(N-1) || + || P(N-1) - P(N-2) || )
-//
-// Output, double BASIS_MATRIX_OVERHAUSER_NUR[3*3], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[3*3];
-
- mbasis[0+0*3] = 1.0 / beta;
- mbasis[0+1*3] = - 1.0 / ( beta * ( 1.0 - beta ) );
- mbasis[0+2*3] = 1.0 / ( 1.0 - beta );
-
- mbasis[1+0*3] = - ( 1.0 + beta ) / beta;
- mbasis[1+1*3] = 1.0 / ( beta * ( 1.0 - beta ) );
- mbasis[1+2*3] = - beta / ( 1.0 - beta );
-
- mbasis[2+0*3] = 1.0;
- mbasis[2+1*3] = 0.0;
- mbasis[2+2*3] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_overhauser_uni ( void)
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_OVERHAUSER_UNI sets the uniform Overhauser spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points P1, P2, P3 and
-// P4 are uniformly spaced in T, and that P2 corresponds to T = 0,
-// and P3 to T = 1.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 14 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// James Foley, Andries vanDam, Steven Feiner, John Hughes,
-// Computer Graphics, Principles and Practice,
-// Second Edition,
-// Addison Wesley, 1995,
-// ISBN: 0201848406,
-// LC: T385.C5735.
-//
-// Parameters:
-//
-// Output, double BASIS_MATRIX_OVERHASUER_UNI[4*4], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[4*4];
-
- mbasis[0+0*4] = - 1.0 / 2.0;
- mbasis[0+1*4] = 3.0 / 2.0;
- mbasis[0+2*4] = - 3.0 / 2.0;
- mbasis[0+3*4] = 1.0 / 2.0;
-
- mbasis[1+0*4] = 2.0 / 2.0;
- mbasis[1+1*4] = - 5.0 / 2.0;
- mbasis[1+2*4] = 4.0 / 2.0;
- mbasis[1+3*4] = - 1.0 / 2.0;
-
- mbasis[2+0*4] = - 1.0 / 2.0;
- mbasis[2+1*4] = 0.0;
- mbasis[2+2*4] = 1.0 / 2.0;
- mbasis[2+3*4] = 0.0;
-
- mbasis[3+0*4] = 0.0;
- mbasis[3+1*4] = 2.0 / 2.0;
- mbasis[3+2*4] = 0.0;
- mbasis[3+3*4] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_overhauser_uni_l ( void )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_OVERHAUSER_UNI_L sets the left uniform Overhauser spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points P1, P2, and P3
-// are not uniformly spaced in T, and that P1 corresponds to T = 0,
-// and P2 to T = 1.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 14 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Output, double BASIS_MATRIX_OVERHASUER_UNI_L[3*3], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[3*3];
-
- mbasis[0+0*3] = 2.0;
- mbasis[0+1*3] = - 4.0;
- mbasis[0+2*3] = 2.0;
-
- mbasis[1+0*3] = - 3.0;
- mbasis[1+1*3] = 4.0;
- mbasis[1+2*3] = - 1.0;
-
- mbasis[2+0*3] = 1.0;
- mbasis[2+1*3] = 0.0;
- mbasis[2+2*3] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double *basis_matrix_overhauser_uni_r ( void )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_OVERHAUSER_UNI_R sets the right uniform Overhauser spline basis matrix.
-//
-// Discussion:
-//
-// This basis matrix assumes that the data points P(N-2), P(N-1),
-// and P(N) are uniformly spaced in T, and that P(N-1) corresponds to
-// T = 0, and P(N) to T = 1.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 14 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Output, double BASIS_MATRIX_OVERHASUER_UNI_R[3*3], the basis matrix.
-//
-{
- double *mbasis;
-
- mbasis = new double[3*3];
-
- mbasis[0+0*3] = 2.0;
- mbasis[0+1*3] = - 4.0;
- mbasis[0+2*3] = 2.0;
-
- mbasis[1+0*3] = - 3.0;
- mbasis[1+1*3] = 4.0;
- mbasis[1+2*3] = - 1.0;
-
- mbasis[2+0*3] = 1.0;
- mbasis[2+1*3] = 0.0;
- mbasis[2+2*3] = 0.0;
-
- return mbasis;
-}
-//****************************************************************************80
-
-double basis_matrix_tmp ( int left, int n, double mbasis[], int ndata,
- double tdata[], double ydata[], double tval )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BASIS_MATRIX_TMP computes Q = T * MBASIS * P
-//
-// Discussion:
-//
-// YDATA is a vector of data values, most frequently the values of some
-// function sampled at uniformly spaced points. MBASIS is the basis
-// matrix for a particular kind of spline. T is a vector of the
-// powers of the normalized difference between TVAL and the left
-// endpoint of the interval.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 14 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int LEFT, indicats that TVAL is in the interval
-// [ TDATA(LEFT), TDATA(LEFT+1) ], or that this is the "nearest"
-// interval to TVAL.
-// For TVAL < TDATA(1), use LEFT = 1.
-// For TDATA(NDATA) < TVAL, use LEFT = NDATA - 1.
-//
-// Input, int N, the order of the basis matrix.
-//
-// Input, double MBASIS[N*N], the basis matrix.
-//
-// Input, int NDATA, the dimension of the vectors TDATA and YDATA.
-//
-// Input, double TDATA[NDATA], the abscissa values. This routine
-// assumes that the TDATA values are uniformly spaced, with an
-// increment of 1.0.
-//
-// Input, double YDATA[NDATA], the data values to be interpolated or
-// approximated.
-//
-// Input, double TVAL, the value of T at which the spline is to be
-// evaluated.
-//
-// Output, double BASIS_MATRIX_TMP, the value of the spline at TVAL.
-//
-{
- double arg;
- int first;
- int i;
- int j;
- double temp;
- double tm;
- double *tvec;
- double yval;
-
- tvec = new double[n];
-
- if ( left == 1 )
- {
- arg = 0.5 * ( tval - tdata[left-1] );
- first = left;
- }
- else if ( left < ndata - 1 )
- {
- arg = tval - tdata[left-1];
- first = left - 1;
- }
- else if ( left == ndata - 1 )
- {
- arg = 0.5 * ( 1.0 + tval - tdata[left-1] );
- first = left - 1;
- }
-//
-// TVEC(I) = ARG**(N-I).
-//
- tvec[n-1] = 1.0;
- for ( i = n-2; 0 <= i; i-- )
- {
- tvec[i] = arg * tvec[i+1];
- }
-
- yval = 0.0;
- for ( j = 0; j < n; j++ )
- {
- tm = 0.0;
- for ( i = 0; i < n; i++ )
- {
- tm = tm + tvec[i] * mbasis[i+j*n];
- }
- yval = yval + tm * ydata[first - 1 + j];
- }
-
- delete [] tvec;
-
- return yval;
-}
-//****************************************************************************80
-
-void bc_val ( int n, double t, double xcon[], double ycon[], double *xval,
- double *yval )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BC_VAL evaluates a parameterized Bezier curve.
-//
-// Discussion:
-//
-// BC_VAL(T) is the value of a vector function of the form
-//
-// BC_VAL(T) = ( X(T), Y(T) )
-//
-// where
-//
-// X(T) = Sum ( 0 <= I <= N ) XCON(I) * BERN(I,N)(T)
-// Y(T) = Sum ( 0 <= I <= N ) YCON(I) * BERN(I,N)(T)
-//
-// BERN(I,N)(T) is the I-th Bernstein polynomial of order N
-// defined on the interval [0,1],
-//
-// XCON(0:N) and YCON(0:N) are the coordinates of N+1 "control points".
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// David Kahaner, Cleve Moler, Steven Nash,
-// Numerical Methods and Software,
-// Prentice Hall, 1989,
-// ISBN: 0-13-627258-4,
-// LC: TA345.K34.
-//
-// Parameters:
-//
-// Input, int N, the order of the Bezier curve, which
-// must be at least 0.
-//
-// Input, double T, the point at which the Bezier curve should
-// be evaluated. The best results are obtained within the interval
-// [0,1] but T may be anywhere.
-//
-// Input, double XCON[0:N], YCON[0:N], the X and Y coordinates
-// of the control points. The Bezier curve will pass through
-// the points ( XCON(0), YCON(0) ) and ( XCON(N), YCON(N) ), but
-// generally NOT through the other control points.
-//
-// Output, double *XVAL, *YVAL, the X and Y coordinates of the point
-// on the Bezier curve corresponding to the given T value.
-//
-{
- double *bval;
- int i;
-
- bval = bp01 ( n, t );
-
- *xval = 0.0;
- for ( i = 0; i <= n; i++ )
- {
- *xval = *xval + xcon[i] * bval[i];
- }
-
- *yval = 0.0;
- for ( i = 0; i <= n; i++ )
- {
- *yval = *yval + ycon[i] * bval[i];
- }
-
- delete [] bval;
-
- return;
-}
-//****************************************************************************80
-
-double bez_val ( int n, double x, double a, double b, double y[] )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BEZ_VAL evaluates a Bezier function at a point.
-//
-// Discussion:
-//
-// The Bezier function has the form:
-//
-// BEZ(X) = Sum ( 0 <= I <= N ) Y(I) * BERN(N,I)( (X-A)/(B-A) )
-//
-// BERN(N,I)(X) is the I-th Bernstein polynomial of order N
-// defined on the interval [0,1],
-//
-// Y(0:N) is a set of coefficients,
-//
-// and if, for I = 0 to N, we define the N+1 points
-//
-// X(I) = ( (N-I)*A + I*B) / N,
-//
-// equally spaced in [A,B], the pairs ( X(I), Y(I) ) can be regarded as
-// "control points".
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// David Kahaner, Cleve Moler, Steven Nash,
-// Numerical Methods and Software,
-// Prentice Hall, 1989,
-// ISBN: 0-13-627258-4,
-// LC: TA345.K34.
-//
-// Parameters:
-//
-// Input, int N, the order of the Bezier function, which
-// must be at least 0.
-//
-// Input, double X, the point at which the Bezier function should
-// be evaluated. The best results are obtained within the interval
-// [A,B] but X may be anywhere.
-//
-// Input, double A, B, the interval over which the Bezier function
-// has been defined. This is the interval in which the control
-// points have been set up. Note BEZ(A) = Y(0) and BEZ(B) = Y(N),
-// although BEZ will not, in general pass through the other
-// control points. A and B must not be equal.
-//
-// Input, double Y[0:N], a set of data defining the Y coordinates
-// of the control points.
-//
-// Output, double BEZ_VAL, the value of the Bezier function at X.
-//
-{
- double *bval;
- int i;
- double value;
- double x01;
-
- if ( b - a == 0.0 )
- {
- cout << "\n";
- cout << "BEZ_VAL - Fatal error!\n";
- cout << " Null interval, A = B = " << a << "\n";
- exit ( 1 );
- }
-//
-// X01 lies in [0,1], in the same relative position as X in [A,B].
-//
- x01 = ( x - a ) / ( b - a );
-
- bval = bp01 ( n, x01 );
-
- value = 0.0;
- for ( i = 0; i <= n; i++ )
- {
- value = value + y[i] * bval[i];
- }
-
- delete [] bval;
-
- return value;
-}
-//****************************************************************************80
-
-double bp_approx ( int n, double a, double b, double ydata[], double xval )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BP_APPROX evaluates the Bernstein polynomial for F(X) on [A,B].
-//
-// Formula:
-//
-// BERN(F)(X) = sum ( 0 <= I <= N ) F(X(I)) * B_BASE(I,X)
-//
-// where
-//
-// X(I) = ( ( N - I ) * A + I * B ) / N
-// B_BASE(I,X) is the value of the I-th Bernstein basis polynomial at X.
-//
-// Discussion:
-//
-// The Bernstein polynomial BERN(F) for F(X) is an approximant, not an
-// interpolant; in other words, its value is not guaranteed to equal
-// that of F at any particular point. However, for a fixed interval
-// [A,B], if we let N increase, the Bernstein polynomial converges
-// uniformly to F everywhere in [A,B], provided only that F is continuous.
-// Even if F is not continuous, but is bounded, the polynomial converges
-// pointwise to F(X) at all points of continuity. On the other hand,
-// the convergence is quite slow compared to other interpolation
-// and approximation schemes.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// David Kahaner, Cleve Moler, Steven Nash,
-// Numerical Methods and Software,
-// Prentice Hall, 1989,
-// ISBN: 0-13-627258-4,
-// LC: TA345.K34.
-//
-// Parameters:
-//
-// Input, int N, the degree of the Bernstein polynomial to be used.
-//
-// Input, double A, B, the endpoints of the interval on which the
-// approximant is based. A and B should not be equal.
-//
-// Input, double YDATA[0:N], the data values at N+1 equally spaced points
-// in [A,B]. If N = 0, then the evaluation point should be 0.5 * ( A + B).
-// Otherwise, evaluation point I should be ( (N-I)*A + I*B ) / N ).
-//
-// Input, double XVAL, the point at which the Bernstein polynomial
-// approximant is to be evaluated. XVAL does not have to lie in the
-// interval [A,B].
-//
-// Output, double BP_APPROX, the value of the Bernstein polynomial approximant
-// for F, based in [A,B], evaluated at XVAL.
-//
-{
- double *bvec;
- int i;
- double yval;
-//
-// Evaluate the Bernstein basis polynomials at XVAL.
-//
- bvec = bpab ( n, a, b, xval );
-//
-// Now compute the sum of YDATA(I) * BVEC(I).
-//
- yval = 0.0;
-
- for ( i = 0; i <= n; i++ )
- {
- yval = yval + ydata[i] * bvec[i];
- }
- delete [] bvec;
-
- return yval;
-}
-//****************************************************************************80
-
-double *bp01 ( int n, double x )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BP01 evaluates the Bernstein basis polynomials for [0,1] at a point.
-//
-// Discussion:
-//
-// For any N greater than or equal to 0, there is a set of N+1 Bernstein
-// basis polynomials, each of degree N, which form a basis for
-// all polynomials of degree N on [0,1].
-//
-// Formula:
-//
-// BERN(N,I,X) = [N!/(I!*(N-I)!)] * (1-X)**(N-I) * X**I
-//
-// N is the degree;
-//
-// 0 <= I <= N indicates which of the N+1 basis polynomials
-// of degree N to choose;
-//
-// X is a point in [0,1] at which to evaluate the basis polynomial.
-//
-// First values:
-//
-// B(0,0,X) = 1
-//
-// B(1,0,X) = 1-X
-// B(1,1,X) = X
-//
-// B(2,0,X) = (1-X)**2
-// B(2,1,X) = 2 * (1-X) * X
-// B(2,2,X) = X**2
-//
-// B(3,0,X) = (1-X)**3
-// B(3,1,X) = 3 * (1-X)**2 * X
-// B(3,2,X) = 3 * (1-X) * X**2
-// B(3,3,X) = X**3
-//
-// B(4,0,X) = (1-X)**4
-// B(4,1,X) = 4 * (1-X)**3 * X
-// B(4,2,X) = 6 * (1-X)**2 * X**2
-// B(4,3,X) = 4 * (1-X) * X**3
-// B(4,4,X) = X**4
-//
-// Special values:
-//
-// B(N,I,1/2) = C(N,K) / 2**N
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// David Kahaner, Cleve Moler, Steven Nash,
-// Numerical Methods and Software,
-// Prentice Hall, 1989,
-// ISBN: 0-13-627258-4,
-// LC: TA345.K34.
-//
-// Parameters:
-//
-// Input, int N, the degree of the Bernstein basis polynomials.
-//
-// Input, double X, the evaluation point.
-//
-// Output, double BP01[0:N], the values of the N+1 Bernstein basis
-// polynomials at X.
-//
-{
- double *bern;
- int i;
- int j;
-
- bern = new double[n+1];
-
- if ( n == 0 )
- {
- bern[0] = 1.0;
- }
- else if ( 0 < n )
- {
- bern[0] = 1.0 - x;
- bern[1] = x;
-
- for ( i = 2; i <= n; i++ )
- {
- bern[i] = x * bern[i-1];
- for ( j = i-1; 1 <= j; j-- )
- {
- bern[j] = x * bern[j-1] + ( 1.0 - x ) * bern[j];
- }
- bern[0] = ( 1.0 - x ) * bern[0];
- }
-
- }
-
- return bern;
-}
-//****************************************************************************80
-
-double *bpab ( int n, double a, double b, double x )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// BPAB evaluates the Bernstein basis polynomials for [A,B] at a point.
-//
-// Formula:
-//
-// BERN(N,I,X) = [N!/(I!*(N-I)!)] * (B-X)**(N-I) * (X-A)**I / (B-A)**N
-//
-// First values:
-//
-// B(0,0,X) = 1
-//
-// B(1,0,X) = ( B-X ) / (B-A)
-// B(1,1,X) = ( X-A ) / (B-A)
-//
-// B(2,0,X) = ( (B-X)**2 ) / (B-A)**2
-// B(2,1,X) = ( 2 * (B-X) * (X-A) ) / (B-A)**2
-// B(2,2,X) = ( (X-A)**2 ) / (B-A)**2
-//
-// B(3,0,X) = ( (B-X)**3 ) / (B-A)**3
-// B(3,1,X) = ( 3 * (B-X)**2 * (X-A) ) / (B-A)**3
-// B(3,2,X) = ( 3 * (B-X) * (X-A)**2 ) / (B-A)**3
-// B(3,3,X) = ( (X-A)**3 ) / (B-A)**3
-//
-// B(4,0,X) = ( (B-X)**4 ) / (B-A)**4
-// B(4,1,X) = ( 4 * (B-X)**3 * (X-A) ) / (B-A)**4
-// B(4,2,X) = ( 6 * (B-X)**2 * (X-A)**2 ) / (B-A)**4
-// B(4,3,X) = ( 4 * (B-X) * (X-A)**3 ) / (B-A)**4
-// B(4,4,X) = ( (X-A)**4 ) / (B-A)**4
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 February 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Reference:
-//
-// David Kahaner, Cleve Moler, Steven Nash,
-// Numerical Methods and Software,
-// Prentice Hall, 1989,
-// ISBN: 0-13-627258-4,
-// LC: TA345.K34.
-//
-// Parameters:
-//
-// Input, integer N, the degree of the Bernstein basis polynomials.
-// For any N greater than or equal to 0, there is a set of N+1
-// Bernstein basis polynomials, each of degree N, which form a basis
-// for polynomials on [A,B].
-//
-// Input, double A, B, the endpoints of the interval on which the
-// polynomials are to be based. A and B should not be equal.
-//
-// Input, double X, the point at which the polynomials are to be
-// evaluated. X need not lie in the interval [A,B].
-//
-// Output, double BERN[0:N], the values of the N+1 Bernstein basis
-// polynomials at X.
-//
-{
- double *bern;
- int i;
- int j;
-
- if ( b == a )
- {
- cout << "\n";
- cout << "BPAB - Fatal error!\n";
- cout << " A = B = " << a << "\n";
- exit ( 1 );
- }
-
- bern = new double[n+1];
-
- if ( n == 0 )
- {
- bern[0] = 1.0;
- }
- else if ( 0 < n )
- {
- bern[0] = ( b - x ) / ( b - a );
- bern[1] = ( x - a ) / ( b - a );
-
- for ( i = 2; i <= n; i++ )
- {
- bern[i] = ( x - a ) * bern[i-1] / ( b - a );
- for ( j = i-1; 1 <= j; j-- )
- {
- bern[j] = ( ( b - x ) * bern[j] + ( x - a ) * bern[j-1] ) / ( b - a );
- }
- bern[0] = ( b - x ) * bern[0] / ( b - a );
- }
- }
-
- return bern;
-}
-//****************************************************************************80
-
-int chfev ( double x1, double x2, double f1, double f2, double d1, double d2,
- int ne, double xe[], double fe[], int next[] )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// CHFEV evaluates a cubic polynomial given in Hermite form.
-//
-// Discussion:
-//
-// This routine evaluates a cubic polynomial given in Hermite form at an
-// array of points. While designed for use by SPLINE_PCHIP_VAL, it may
-// be useful directly as an evaluator for a piecewise cubic
-// Hermite function in applications, such as graphing, where
-// the interval is known in advance.
-//
-// The cubic polynomial is determined by function values
-// F1, F2 and derivatives D1, D2 on the interval [X1,X2].
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 12 August 2005
-//
-// Author:
-//
-// Original FORTRAN77 version by Fred Fritsch, Lawrence Livermore National Laboratory.
-// C++ version by John Burkardt.
-//
-// Reference:
-//
-// Fred Fritsch, Ralph Carlson,
-// Monotone Piecewise Cubic Interpolation,
-// SIAM Journal on Numerical Analysis,
-// Volume 17, Number 2, April 1980, pages 238-246.
-//
-// David Kahaner, Cleve Moler, Steven Nash,
-// Numerical Methods and Software,
-// Prentice Hall, 1989,
-// ISBN: 0-13-627258-4,
-// LC: TA345.K34.
-//
-// Parameters:
-//
-// Input, double X1, X2, the endpoints of the interval of
-// definition of the cubic. X1 and X2 must be distinct.
-//
-// Input, double F1, F2, the values of the function at X1 and
-// X2, respectively.
-//
-// Input, double D1, D2, the derivative values at X1 and
-// X2, respectively.
-//
-// Input, int NE, the number of evaluation points.
-//
-// Input, double XE[NE], the points at which the function is to
-// be evaluated. If any of the XE are outside the interval
-// [X1,X2], a warning error is returned in NEXT.
-//
-// Output, double FE[NE], the value of the cubic function
-// at the points XE.
-//
-// Output, int NEXT[2], indicates the number of extrapolation points:
-// NEXT[0] = number of evaluation points to the left of interval.
-// NEXT[1] = number of evaluation points to the right of interval.
-//
-// Output, int CHFEV, error flag.
-// 0, no errors.
-// -1, NE < 1.
-// -2, X1 == X2.
-//
-{
- double c2;
- double c3;
- double del1;
- double del2;
- double delta;
- double h;
- int i;
- int ierr;
- double x;
- double xma;
- double xmi;
-
- if ( ne < 1 )
- {
- ierr = -1;
- cout << "\n";
- cout << "CHFEV - Fatal error!\n";
- cout << " Number of evaluation points is less than 1.\n";
- cout << " NE = " << ne << "\n";
- return ierr;
- }
-
- h = x2 - x1;
-
- if ( h == 0.0 )
- {
- ierr = -2;
- cout << "\n";
- cout << "CHFEV - Fatal error!\n";
- cout << " The interval [X1,X2] is of zero length.\n";
- return ierr;
- }
-//
-// Initialize.
-//
- ierr = 0;
- next[0] = 0;
- next[1] = 0;
- xmi = r8_min ( 0.0, h );
- xma = r8_max ( 0.0, h );
-//
-// Compute cubic coefficients expanded about X1.
-//
- delta = ( f2 - f1 ) / h;
- del1 = ( d1 - delta ) / h;
- del2 = ( d2 - delta ) / h;
- c2 = -( del1 + del1 + del2 );
- c3 = ( del1 + del2 ) / h;
-//
-// Evaluation loop.
-//
- for ( i = 0; i < ne; i++ )
- {
- x = xe[i] - x1;
- fe[i] = f1 + x * ( d1 + x * ( c2 + x * c3 ) );
-//
-// Count the extrapolation points.
-//
- if ( x < xmi )
- {
- next[0] = next[0] + 1;
- }
-
- if ( xma < x )
- {
- next[1] = next[1] + 1;
- }
-
- }
-
- return 0;
-}
-//****************************************************************************80
-
-double *d3_mxv ( int n, double a[], double x[] )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// D3_MXV multiplies a D3 matrix times a vector.
-//
-// Discussion:
-//
-// The D3 storage format is used for a tridiagonal matrix.
-// The superdiagonal is stored in entries (1,2:N), the diagonal in
-// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
-// original matrix is "collapsed" vertically into the array.
-//
-// Example:
-//
-// Here is how a D3 matrix of order 5 would be stored:
-//
-// * A12 A23 A34 A45
-// A11 A22 A33 A44 A55
-// A21 A32 A43 A54 *
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 15 November 2003
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int N, the order of the linear system.
-//
-// Input, double A[3*N], the D3 matrix.
-//
-// Input, double X[N], the vector to be multiplied by A.
-//
-// Output, double D3_MXV[N], the product A * x.
-//
-{
- double *b;
- int i;
-
- b = new double[n];
-
- for ( i = 0; i < n; i++ )
- {
- b[i] = a[1+i*3] * x[i];
- }
- for ( i = 0; i < n-1; i++ )
- {
- b[i] = b[i] + a[0+(i+1)*3] * x[i+1];
- }
- for ( i = 1; i < n; i++ )
- {
- b[i] = b[i] + a[2+(i-1)*3] * x[i-1];
- }
-
- return b;
-}
-//****************************************************************************80
-
-double *d3_np_fs ( int n, double a[], double b[] )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// D3_NP_FS factors and solves a D3 system.
-//
-// Discussion:
-//
-// The D3 storage format is used for a tridiagonal matrix.
-// The superdiagonal is stored in entries (1,2:N), the diagonal in
-// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
-// original matrix is "collapsed" vertically into the array.
-//
-// This algorithm requires that each diagonal entry be nonzero.
-// It does not use pivoting, and so can fail on systems that
-// are actually nonsingular.
-//
-// Example:
-//
-// Here is how a D3 matrix of order 5 would be stored:
-//
-// * A12 A23 A34 A45
-// A11 A22 A33 A44 A55
-// A21 A32 A43 A54 *
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 15 November 2003
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int N, the order of the linear system.
-//
-// Input/output, double A[3*N].
-// On input, the nonzero diagonals of the linear system.
-// On output, the data in these vectors has been overwritten
-// by factorization information.
-//
-// Input, double B[N], the right hand side.
-//
-// Output, double D3_NP_FS[N], the solution of the linear system.
-// This is NULL if there was an error because one of the diagonal
-// entries was zero.
-//
-{
- int i;
- double *x;
- double xmult;
-//
-// Check.
-//
- for ( i = 0; i < n; i++ )
- {
- if ( a[1+i*3] == 0.0 )
- {
- return NULL;
- }
- }
- x = new double[n];
-
- for ( i = 0; i < n; i++ )
- {
- x[i] = b[i];
- }
-
- for ( i = 1; i < n; i++ )
- {
- xmult = a[2+(i-1)*3] / a[1+(i-1)*3];
- a[1+i*3] = a[1+i*3] - xmult * a[0+i*3];
- x[i] = x[i] - xmult * x[i-1];
- }
-
- x[n-1] = x[n-1] / a[1+(n-1)*3];
- for ( i = n-2; 0 <= i; i-- )
- {
- x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3];
- }
-
- return x;
-}
-//****************************************************************************80
-
-void d3_print ( int n, double a[], char *title )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// D3_PRINT prints a D3 matrix.
-//
-// Discussion:
-//
-// The D3 storage format is used for a tridiagonal matrix.
-// The superdiagonal is stored in entries (1,2:N), the diagonal in
-// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
-// original matrix is "collapsed" vertically into the array.
-//
-// Example:
-//
-// Here is how a D3 matrix of order 5 would be stored:
-//
-// * A12 A23 A34 A45
-// A11 A22 A33 A44 A55
-// A21 A32 A43 A54 *
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 20 September 2003
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int N, the order of the matrix.
-// N must be positive.
-//
-// Input, double A[3*N], the D3 matrix.
-//
-// Input, char *TITLE, a title to print.
-//
-{
- if ( 0 < s_len_trim ( title ) )
- {
- cout << "\n";
- cout << title << "\n";
- }
-
- cout << "\n";
-
- d3_print_some ( n, a, 1, 1, n, n );
-
- return;
-}
-//****************************************************************************80
-
-void d3_print_some ( int n, double a[], int ilo, int jlo, int ihi, int jhi )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// D3_PRINT_SOME prints some of a D3 matrix.
-//
-// Discussion:
-//
-// The D3 storage format is used for a tridiagonal matrix.
-// The superdiagonal is stored in entries (1,2:N), the diagonal in
-// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
-// original matrix is "collapsed" vertically into the array.
-//
-// Example:
-//
-// Here is how a D3 matrix of order 5 would be stored:
-//
-// * A12 A23 A34 A45
-// A11 A22 A33 A44 A55
-// A21 A32 A43 A54 *
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 05 January 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int N, the order of the matrix.
-// N must be positive.
-//
-// Input, double A[3*N], the D3 matrix.
-//
-// Input, int ILO, JLO, IHI, JHI, designate the first row and
-// column, and the last row and column, to be printed.
-//
-{
-# define INCX 5
-
- int i;
- int i2hi;
- int i2lo;
- int inc;
- int j;
- int j2;
- int j2hi;
- int j2lo;
-//
-// Print the columns of the matrix, in strips of 5.
-//
- for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
- {
- j2hi = j2lo + INCX - 1;
- j2hi = i4_min ( j2hi, n );
- j2hi = i4_min ( j2hi, jhi );
-
- inc = j2hi + 1 - j2lo;
-
- cout << "\n";
- cout << " Col: ";
- for ( j = j2lo; j <= j2hi; j++ )
- {
- j2 = j + 1 - j2lo;
- cout << setw(7) << j << " ";
- }
-
- cout << "\n";
- cout << " Row\n";
- cout << " ---\n";
-//
-// Determine the range of the rows in this strip.
-//
- i2lo = i4_max ( ilo, 1 );
- i2lo = i4_max ( i2lo, j2lo - 1 );
-
- i2hi = i4_min ( ihi, n );
- i2hi = i4_min ( i2hi, j2hi + 1 );
-
- for ( i = i2lo; i <= i2hi; i++ )
- {
-//
-// Print out (up to) 5 entries in row I, that lie in the current strip.
-//
- cout << setw(6) << i << " ";
-
- for ( j2 = 1; j2 <= inc; j2++ )
- {
- j = j2lo - 1 + j2;
-
- if ( 1 < i-j || 1 < j-i )
- {
- cout << " ";
- }
- else if ( j == i+1 )
- {
- cout << setw(12) << a[0+(j-1)*3] << " ";
- }
- else if ( j == i )
- {
- cout << setw(12) << a[1+(j-1)*3] << " ";
- }
- else if ( j == i-1 )
- {
- cout << setw(12) << a[2+(j-1)*3] << " ";
- }
-
- }
- cout << "\n";
- }
- }
-
- cout << "\n";
-
- return;
-# undef INCX
-}
-//****************************************************************************80
-
-double *d3_uniform ( int n, int *seed )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// D3_UNIFORM randomizes a D3 matrix.
-//
-// Discussion:
-//
-// The D3 storage format is used for a tridiagonal matrix.
-// The superdiagonal is stored in entries (1,2:N), the diagonal in
-// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
-// original matrix is "collapsed" vertically into the array.
-//
-// Example:
-//
-// Here is how a D3 matrix of order 5 would be stored:
-//
-// * A12 A23 A34 A45
-// A11 A22 A33 A44 A55
-// A21 A32 A43 A54 *
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 January 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int N, the order of the linear system.
-//
-// Input/output, int *SEED, a seed for the random number generator.
-//
-// Output, double D3_UNIFORM[3*N], the D3 matrix.
-//
-{
- double *a;
- int i;
- double *u;
- double *v;
- double *w;
-
- a = new double[3*n];
-
- u = r8vec_uniform ( n-1, 0.0, 1.0, seed );
- v = r8vec_uniform ( n, 0.0, 1.0, seed );
- w = r8vec_uniform ( n-1, 0.0, 1.0, seed );
-
- a[0+0*3] = 0.0;
- for ( i = 1; i < n; i++ )
- {
- a[0+i*3] = u[i-1];
- }
- for ( i = 0; i < n; i++ )
- {
- a[1+i*3] = v[i];
- }
- for ( i = 0; i < n-1; i++ )
- {
- a[2+i*3] = w[i];
- }
- a[2+(n-1)*3] = 0.0;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- return a;
-}
-//****************************************************************************80
-
-void data_to_dif ( int ntab, double xtab[], double ytab[], double diftab[] )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// DATA_TO_DIF sets up a divided difference table from raw data.
-//
-// Discussion:
-//
-// Space can be saved by using a single array for both the DIFTAB and
-// YTAB dummy parameters. In that case, the difference table will
-// overwrite the Y data without interfering with the computation.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 04 September 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int NTAB, the number of pairs of points
-// (XTAB[I],YTAB[I]) which are to be used as data.
-//
-// Input, double XTAB[NTAB], the X values at which data was taken.
-// These values must be distinct.
-//
-// Input, double YTAB[NTAB], the corresponding Y values.
-//
-// Output, double DIFTAB[NTAB], the divided difference coefficients
-// corresponding to the input (XTAB,YTAB).
-//
-{
- int i;
- int j;
-//
-// Copy the data values into DIFTAB.
-//
- for ( i = 0; i < ntab; i++ )
- {
- diftab[i] = ytab[i];
- }
-//
-// Make sure the abscissas are distinct.
-//
- for ( i = 0; i < ntab; i++ )
- {
- for ( j = i+1; j < ntab; j++ )
- {
- if ( xtab[i] - xtab[j] == 0.0 )
- {
- cout << "\n";
- cout << "DATA_TO_DIF - Fatal error!\n";
- cout << " Two entries of XTAB are equal!\n";
- cout << " XTAB[%d] = " << xtab[i] << "\n";
- cout << " XTAB[%d] = " << xtab[j] << "\n";
- exit ( 1 );
- }
- }
- }
-//
-// Compute the divided differences.
-//
- for ( i = 1; i <= ntab-1; i++ )
- {
- for ( j = ntab-1; i <= j; j-- )
- {
- diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] );
- }
- }
-
- return;
-}
-//****************************************************************************80
-
-double dif_val ( int ntab, double xtab[], double diftab[], double xval )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// DIF_VAL evaluates a divided difference polynomial at a point.
-//
-// Discussion:
-//
-// DATA_TO_DIF must be called first to set up the divided difference table.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 05 September 2004
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, integer NTAB, the number of divided difference
-// coefficients in DIFTAB, and the number of points XTAB.
-//
-// Input, double XTAB[NTAB], the X values upon which the
-// divided difference polynomial is based.
-//
-// Input, double DIFTAB[NTAB], the divided difference table.
-//
-// Input, double XVAL, a value of X at which the polynomial
-// is to be evaluated.
-//
-// Output, double DIF_VAL, the value of the polynomial at XVAL.
-//
-{
- int i;
- double value;
-
- value = diftab[ntab-1];
- for ( i = 2; i <= ntab; i++ )
- {
- value = diftab[ntab-i] + ( xval - xtab[ntab-i] ) * value;
- }
-
- return value;
-}
-//****************************************************************************80
-
-int i4_max ( int i1, int i2 )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// I4_MAX returns the maximum of two I4's.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 October 1998
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int I1, I2, are two integers to be compared.
-//
-// Output, int I4_MAX, the larger of I1 and I2.
-//
-//
-{
- if ( i2 < i1 )
- {
- return i1;
- }
- else
- {
- return i2;
- }
-
-}
-//****************************************************************************80
-
-int i4_min ( int i1, int i2 )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// I4_MIN returns the smaller of two I4's.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 13 October 1998
-//
-// Author:
-//
-// John Burkardt
-//
-// Parameters:
-//
-// Input, int I1, I2, two integers to be compared.
-//
-// Output, int I4_MIN, the smaller of I1 and I2.
-//
-//
-{
- if ( i1 < i2 )
- {
- return i1;
- }
- else
- {
- return i2;
- }
-
-}
-//****************************************************************************80
-
-void least_set ( int point_num, double x[], double f[], double w[],
- int nterms, double b[], double c[], double d[] )
-
-//****************************************************************************80
-//
-// Purpose:
-//
-// LEAST_SET defines a least squares polynomial for given data.
-//
-// Discussion:
-//
-// This routine is based on ORTPOL by Conte and deBoor.
-//
-// The polynomial may be evaluated at any point X by calling LEAST_VAL.
-//
-// Thanks to Andrew Telford for pointing out a mistake in the form of
-// the check that there are enough unique data points, 25 June 2008.
-//
-// Licensing:
-//
-// This code is distributed under the GNU LGPL license.
-//
-// Modified:
-//
-// 25 June 2008
-//
-// Author:
-//
-// Original FORTRAN77 version by Samuel Conte, Carl deBoor.
-// C++ version by John Burkardt
-//
-// Reference:
-//
-// Samuel Conte, Carl deBoor,
-// Elementary Numerical Analysis,
-// Second Edition,
-// McGraw Hill, 1972,
-// ISBN: 07-012446-4,
-// LC: QA297.C65.
-//
-// Parameters:
-//
-// Input, int POINT_NUM, the number of data values.
-//
-// Input, double X[POINT_NUM], the abscissas of the data points.
-// At least NTERMS of the values in X must be distinct.
-//
-// Input, double F[POINT_NUM], the data values at the points X(*).
-//
-// Input, double W[POINT_NUM], the weights associated with
-// the data points. Each entry of W should be positive.
-//
-// Input, int NTERMS, the number of terms to use in the
-// approximating polynomial. NTERMS must be at least 1.
-// The degree of the polynomial is NTERMS-1.
-//
-// Output, double B[NTERMS], C[NTERMS], D[NTERMS], are quantities
-// defining the least squares polynomial for the input data,
-// which will be needed to evaluate the polynomial.
-//
-{
- int i;
- int j;
- int k;
- double p;
- double *pj;
- double *pjm1;
- double *s;
- double tol = 0.0;
- int unique_num;
-//
-// Make sure at least NTERMS X values are unique.
-//
- unique_num = r8vec_unique_count ( point_num, x, tol );
-
- if ( unique_num < nterms )
- {
- cout << "\n";
- cout << "LEAST_SET - Fatal error!\n";
- cout << " The number of distinct X values must be\n";
- cout << " at least NTERMS = " << nterms << "\n";
- cout << " but the input data has only " << unique_num << "\n";
- cout << " distinct entries.\n";
- return;
- }
-//
-// Make sure all W entries are positive.
-//
- for ( i = 0; i < point_num; i++ )
- {
- if ( w[i] <= 0.0 )
- {
- cout << "\n";
- cout << "LEAST_SET - Fatal error!\n";
- cout << " All weights W must be positive,\n";
- cout << " but weight " << i << "\n";
- cout << " is " << w[i] << "\n";
- return;
- }
- }
-
- s = new double[nterms];
-//
-/...
[truncated message content] |