# include # include # include # include # include # include "edge.h" /******************************************************************************/ double fx1 ( double x ) /******************************************************************************/ /* Purpose: FX1 is 1D example #1, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX1, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx1_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fx2 ( double x ) /******************************************************************************/ /* Purpose: FX2 is 1D example #2, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX2, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx2_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fx3 ( double x ) /******************************************************************************/ /* Purpose: FX3 is 1D example #3, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX3, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx3_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fx4 ( double x ) /******************************************************************************/ /* Purpose: FX4 is 1D example #4, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX4, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx4_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fx5 ( double x ) /******************************************************************************/ /* Purpose: FX5 is 1D example #5, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX5, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx5_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fx6 ( double x ) /******************************************************************************/ /* Purpose: FX6 is 1D example #6, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 23 September 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX6, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx6_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fx7 ( double x ) /******************************************************************************/ /* Purpose: FX7 is 1D example #7, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 23 September 2014 Author: John Burkardt Parameters: Input, double X, the argument. Output, double FX7, the function value. */ { double value; double *value_vec; double x_vec[1]; x_vec[0] = x; value_vec = fx7_vec ( 1, x_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fxy1 ( double x, double y ) /******************************************************************************/ /* Purpose: FXY1 is the first 2D example, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt */ { double value; double *value_vec; double x_vec[1]; double y_vec[1]; x_vec[0] = x; y_vec[0] = y; value_vec = fxy1_vec ( 1, x_vec, y_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fxy2 ( double x, double y ) /******************************************************************************/ /* Purpose: FXY2 is the second 2D example, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt */ { double value; double *value_vec; double x_vec[1]; double y_vec[1]; x_vec[0] = x; y_vec[0] = y; value_vec = fxy2_vec ( 1, x_vec, y_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fxy3 ( double x, double y ) /******************************************************************************/ /* Purpose: FXY3 is the third 2D example, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 19 September 2014 Author: John Burkardt */ { double value; double *value_vec; double x_vec[1]; double y_vec[1]; x_vec[0] = x; y_vec[0] = y; value_vec = fxy3_vec ( 1, x_vec, y_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fxy4 ( double x, double y ) /******************************************************************************/ /* Purpose: FXY4 is the fourth 2D example, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 21 September 2014 Author: John Burkardt */ { double value; double *value_vec; double x_vec[1]; double y_vec[1]; x_vec[0] = x; y_vec[0] = y; value_vec = fxy4_vec ( 1, x_vec, y_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fxy5 ( double x, double y ) /******************************************************************************/ /* Purpose: FXY5 is the fifth 2D example, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 21 September 2014 Author: John Burkardt */ { double value; double *value_vec; double x_vec[1]; double y_vec[1]; x_vec[0] = x; y_vec[0] = y; value_vec = fxy5_vec ( 1, x_vec, y_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double fxyz1 ( double x, double y, double z ) /******************************************************************************/ /* Purpose: FXYZ1 is the first 3D example, scalar version. Discussion: This function allows the user a more convenient interface when only a single input argument is supplied. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt */ { double value; double *value_vec; double x_vec[1]; double y_vec[1]; double z_vec[1]; x_vec[0] = x; y_vec[0] = y; z_vec[0] = z; value_vec = fxyz1_vec ( 1, x_vec, y_vec, z_vec ); value = value_vec[0]; free ( value_vec ); return value; } /******************************************************************************/ double *fx1_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX1_VEC is 1D example #1, vector version. Discussion: This is example 3.1 in the reference. The function should be plotted over [-1.0,+1.0]. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX1_VEC[N], the function values. Local parameters: Local, real STEEP, controls the steepness of the slope. The default value is a moderate 5. For a sharp rise, use 25 instead. */ { double *f; int i; const double r8_pi = 3.141592653589793; const double steep = 5.0; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { if ( x[i] < 0.0 ) { f[i] = cos ( 3.0 * r8_pi * x[i] ); } else if ( 0.0 <= x[i] ) { f[i] = - 1.0 + 2.0 / ( 1.0 + 3.0 * exp ( - steep * ( 2.0 * x[i] - 1.0 ) ) ); } } return f; } /******************************************************************************/ double *fx2_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX2_VEC is 1D example #2, vector version. Discussion: The function should be plotted over [-1,+1]. The "internal" coordinate range will be [-2.0,6.0*pi]. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX2_VEC[N], the function values. */ { double *f; int i; const double r8_pi = 3.141592653589793; double x2; f = ( double * ) malloc ( n * sizeof ( double ) ); /* Map from the convenient range [-1,+1] to the physical range [-2,6pi]. */ for ( i = 0; i < n; i++ ) { x2 = ( ( 1.0 - x[i] ) * ( - 2.0 ) + ( 1.0 + x[i] ) * 6.0 * r8_pi ) / 2.0; if ( x2 < 0.0 ) { f[i] = exp ( x2 ); } else if ( 0.0 <= x2 && x2 < 3.0 * r8_pi / 2.0 ) { f[i] = - exp ( - x2 ); } else if ( 3.0 * r8_pi / 2.0 <= x2 ) { f[i] = -1.5 * sin ( x2 ); } } return f; } /******************************************************************************/ double *fx3_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX3_VEC is 1D example #3, vector version. Discussion: The function should be plotted over [-1.0,+1.0]. Internally, this range is mapped to [-3.0,+3.0]. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX3_VEC[N], the function values. */ { double *f; int i; double x2; f = ( double * ) malloc ( n * sizeof ( double ) ); /* Map from the convenient range [-1,+1] to the physical range [-3,+3]. */ for ( i = 0; i < n; i++ ) { x2 = ( ( 1.0 - x[i] ) * ( -3.0 ) + ( 1.0 + x[i] ) * ( +3.0 ) ) / 2.0; if ( -2.0 <= x2 && x2 <= -1.0 ) { f[i] = 1.0; } else if ( -0.5 <= x2 && x2 <= 0.5 ) { f[i] = 0.5 + 4.0 * pow ( x2 + 0.5, 2 ); } else if ( 1.25 <= x2 && 3.0 * x2 <= 7.0 ) { f[i] = 3.0 * ( 2.0 - x2 ); } else { f[i] = 0.0; } } return f; } /******************************************************************************/ double *fx4_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX4_VEC is 1D example #4, vector version. Discussion: The function should be plotted over [0.0,+1.0]. The function is continuous, but the derivative has a discontinuity at 0.5. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX4_VEC[N], the function values. */ { double *f; int i; const double r8_pi = 3.141592653589793; double x2; f = ( double * ) malloc ( n * sizeof ( double ) ); /* Convert from -1 <= x <= 1 to 0 <= x <= 1: */ for ( i = 0; i < n; i++ ) { x2 = ( x[i] + 1.0 ) / 2.0; if ( x2 <= 0.5 ) { f[i] = - ( x2 - 0.5 ) + sin ( 4.0 * r8_pi * x2 ) / 6.0; } else if ( 0.5 < x2 ) { f[i] = ( x2 - 0.5 ) + sin ( 4.0 * r8_pi * x2 ) / 6.0; } } return f; } /******************************************************************************/ double *fx5_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX5_VEC is 1D example #5, vector version. Discussion: The function should be plotted over [-1.0,+1.0]. The function actually has no discontinuities, but does have a steep rise. The local parameter S controls the steepness of the rise. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX5_VEC[N], the function values. */ { double *f; int i; const double steep = 20.0; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { f[i] = tanh ( steep * x[i] ); } return f; } /******************************************************************************/ double *fx6_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX6_VEC is 1D example #6, vector version. Discussion: This is example 2.1 in the reference. The function should be plotted over [0.0,+1.0]. The function has a discontinuous first derivative at 1/2. Licensing: This code is distributed under the MIT license. Modified: 23 September 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Determining the location of discontinuities in the derivatives of functions, Applied Numerical Mathematics, Volume 58, 2008, pages 577-592. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX6_VEC[N], the function values. */ { double *f; int i; const double r8_pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { f[i] = sin ( 2.0 * r8_pi * x[i] ) / 6.0; if ( x[i] < 0.5 ) { f[i] = f[i] - ( x[i] - 0.5 ); } else { f[i] = f[i] + ( x[i] - 0.5 ); } } return f; } /******************************************************************************/ double *fx7_vec ( int n, double x[] ) /******************************************************************************/ /* Purpose: FX7_VEC is 1D example #7, vector version. Discussion: This is example 2.1 in the reference. The function should be plotted over [0.0,+1.0]. The function has a discontinuous second derivative at 1/2. Licensing: This code is distributed under the MIT license. Modified: 23 September 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Determining the location of discontinuities in the derivatives of functions, Applied Numerical Mathematics, Volume 58, 2008, pages 577-592. Parameters: Input, int N, the number of points. Input, double X[N], the arguments. Output, double FX7_VEC[N], the function values. */ { double *f; int i; const double r8_pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { f[i] = sin ( 2.0 * r8_pi * x[i] ) / 6.0; if ( x[i] < 0.5 ) { f[i] = f[i] - 0.5 * pow ( x[i] - 0.5, 2 ); } else { f[i] = f[i] + 0.5 * pow ( x[i] - 0.5, 2 ); } } return f; } /******************************************************************************/ double *fxy1_vec ( int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: FXY1_VEC is the first 2D example, vector version. Discussion: This is example 4.1 in the reference. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Parameters: Input, int N, the number of points. Input, double X[N], Y[N], the arguments. Output, double FXY1_VEC[N], the function values. */ { double *f; int i; const double r8_pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { f[i] = x[i] * y[i] + cos ( 2.0 * r8_pi * x[i] * x[i] ) - sin ( 2.0 * r8_pi * x[i] * x[i] ); if ( 0.25 < x[i] * x[i] + y[i] * y[i] ) { f[i] = f[i] + 10.0 * x[i] - 5.0; } } return f; } /******************************************************************************/ double *fxy2_vec ( int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: FXY2_VEC is the second 2D example, vector version. Discussion: This is example 4.2 in the reference. It is known as the Shepp-Logan phantom. It should be plotted on [-1,+1] x [-1,+1]. Note that the Archibald reference incorrectly describes the divisor of x in the second ellipse as 0.06624, when it should be 0.6624. Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, Volume 43, Number 1, 2006, pages 259-279. Larry Shepp, Ben Logan, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science, Volume NS-21, June 1974, pages 21-43. Parameters: Input, int N, the number of points. Input, double X[N], Y[N], the arguments. Output, double FXY2_VEC[N], the function values. Local parameters: Local, integer CHOICE: 1, use Archibald's (and Shepp and Logan's) level values; 2, use Matlab's level values; 3, use Matlab's enhanced contrast level values. */ { double *c; double c1[4] = { 2.0, -0.98, -0.02, +0.01 }; double c2[4] = { 1.0, -0.98, -0.02, +0.01 }; double c3[4] = { 1.0, -0.8, -0.2, +0.1 }; int choice; double eta1; double eta2; double *f; int i; const double r8_pi = 3.141592653589793; double xi1; double xi2; f = ( double * ) malloc ( n * sizeof ( double ) ); choice = 3; if ( choice == 1 ) { c = r8vec_copy_new ( 4, c1 ); } else if ( choice == 2 ) { c = r8vec_copy_new ( 4, c2 ); } else { c = r8vec_copy_new ( 4, c3 ); } for ( i = 0; i < n; i++ ) { f[i] = 0.0; xi1 = ( x[i] - 0.22 ) * cos ( 0.4 * r8_pi ) + y[i] * sin ( 0.4 * r8_pi ); eta1 = - ( x[i] - 0.22 ) * sin ( 0.4 * r8_pi ) + y[i] * cos ( 0.4 * r8_pi ); xi2 = ( x[i] + 0.22 ) * cos ( 0.6 * r8_pi ) + y[i] * sin ( 0.6 * r8_pi ); eta2 = - ( x[i] + 0.22 ) * sin ( 0.6 * r8_pi ) + y[i] * cos ( 0.6 * r8_pi ); if ( pow ( x[i] / 0.69, 2 ) + pow ( y[i] / 0.92, 2 ) <= 1.0 ) { f[i] = f[i] + c[0]; } if ( pow ( x[i] / 0.6624, 2 ) + pow ( ( y[i] + 0.0184 ) / 0.874, 2 ) <= 1.0 ) { f[i] = f[i] + c[1]; } if ( pow ( xi1 / 0.31, 2 ) + pow ( eta1 / 0.11, 2 ) <= 1.0 || pow ( xi2 / 0.41, 2 ) + pow ( eta2 / 0.16, 2 ) <= 1.0 ) { f[i] = f[i] + c[2]; } if ( ( pow ( ( x[i] - 0.35 ) / 0.3, 2 ) + pow ( y[i] / 0.6, 2 ) <= 1.0 ) || ( pow ( x[i] / 0.21, 2 ) + pow ( ( y[i] - 0.35 ) / 0.25, 2 ) <= 1.0 ) || ( pow ( x[i] / 0.046, 2 ) + pow ( ( y[i] - 0.1 ) / 0.046, 2 ) <= 1.0 ) || ( pow ( x[i] / 0.046, 2 ) + pow ( ( y[i] + 0.1 ) / 0.046, 2 ) <= 1.0 ) || ( pow ( ( x[i] + 0.08 ) / 0.046, 2 ) + pow ( ( y[i] + 0.605 ) / 0.023, 2 ) <= 1.0 ) || ( pow ( x[i] / 0.023, 2 ) + pow ( ( y[i] + 0.605 ) / 0.023, 2 ) <= 1.0 ) || ( pow ( ( x[i] - 0.06 ) / 0.023, 2 ) + pow ( ( y[i] + 0.605 ) / 0.023, 2 ) <= 1.0 ) ) { f[i] = f[i] + c[3]; } } free ( c ); return f; } /******************************************************************************/ double *fxy3_vec ( int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: FXY3_VEC is the third 2D example, vector version. Discussion: This is example 3.2 in the reference. It is known as the modified two-dimensional Harten function. It should be plotted on [-1,+1] x [-1,+1]. Licensing: This code is distributed under the MIT license. Modified: 19 September 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Determining the locations and discontinuities in the derivatives of functions, Applied Numerical Mathematics, Volume 58, 2008, pages 577-592. Parameters: Input, int N, the number of points. Input, double X[N], Y[N], the arguments. Output, double FXY3_VEC[N], the function values. */ { double *f; int i; double r; const double r8_pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { r = ( 4.0 * x[i] * x[i] + 4.0 * y[i] * y[i] - 1.0 ) / 6.0; if ( 3.0 * r <= -1.0 ) { f[i] = - r * sin ( 0.5 * r8_pi * r * r ); } else if ( 3.0 * r < 1.0 ) { f[i] = fabs ( sin ( 2.0 * r8_pi * r ) ); } else { f[i] = 2.0 * r - 1.0 - sin ( 3.0 * r8_pi * r ) / 6.0; } } return f; } /******************************************************************************/ double *fxy4_vec ( int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: FXY4_VEC is the fourth 2D example, vector version. Discussion: This is example 3.1 in the reference. It is known as the discontinuous medium wave function. Here, we are computing the first component of the solution, P(X,Y). It should be plotted on (x,y) in [-1,0]x[0,0.1]. The second variable y actually represents time. Note that in the reference, the formula reads: f(i) = 2.0D+00 * rhor * cr / ( rhol * cl + rhor * cr ) * sin ( r8_pi * omega * ( y(i) - ( x(i) + 0.5D+00 ) / cr ) ) but I believe this should be: f(i) = 2.0D+00 * rhor * cr / ( rhol * cl + rhor * cr ) * sin ( r8_pi * omega * ( y(i) - ( x(i) + 0.5D+00 ) / cl ) ) Licensing: This code is distributed under the MIT license. Modified: 21 September 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Determining the location of discontinuities in the derivatives of functions, Applied Numerical Mathematics, Volume 58, 2008, pages 577-592. Parameters: Input, int N, the number of points. Input, double X[N], Y[N], the arguments. Output, double FXY4_VEC[N], the function values. */ { const double cl = 0.87879; const double cr = 1.0; double *f; int i; const double omega = 12.0; const double rhol = 0.55556; const double rhor = 1.0; const double r8_pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { if ( x[i] <= -0.5 ) { f[i] = sin ( r8_pi * omega * ( y[i] - ( x[i] + 0.5 ) / cl ) ) - ( rhol * cl - rhor * cr ) / ( rhol * cl + rhor * cr ) * sin ( r8_pi * omega * ( y[i] + ( x[i] + 0.5 ) / cl ) ); } else { f[i] = 2.0 * rhor * cr / ( rhol * cl + rhor * cr ) * sin ( r8_pi * omega * ( y[i] - ( x[i] + 0.5 ) / cl ) ); } } return f; } /******************************************************************************/ double *fxy5_vec ( int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: FXY5_VEC is the fifth 2D example, vector version. Discussion: This is example 3.1 in the reference. It is known as the discontinuous medium wave function. Here, we are computing the second component of the solution, U(X,Y). It should be plotted on (x,y) in [-1,0]x[0,0.1]. The second variable y actually represents time. Licensing: This code is distributed under the MIT license. Modified: 21 September 2014 Author: John Burkardt Reference: Rick Archibald, Anne Gelb, Jungho Yoon, Determining the location of discontinuities in the derivatives of functions, Applied Numerical Mathematics, Volume 58, 2008, pages 577-592. Parameters: Input, int N, the number of points. Input, double X[N], Y[N], the arguments. Output, double FXY4_VEC[N], the function values. */ { const double cl = 0.87879; const double cr = 1.0; double *f; int i; const double omega = 12.0; const double rhol = 0.55556; const double rhor = 1.0; const double r8_pi = 3.141592653589793; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { if ( x[i] <= -0.5 ) { f[i] = sin ( r8_pi * omega * ( y[i] - ( x[i] + 0.5 ) / cl ) ) + ( rhol * cl - rhor * cr ) / ( rhol * cl + rhor * cr ) / ( rhol * cl ) * sin ( r8_pi * omega * ( y[i] + ( x[i] + 0.5 ) / cl ) ); } else { f[i] = 2.0 / ( rhol * cl + rhor * cr ) * sin ( r8_pi * omega * ( y[i] - ( x[i] + 0.5 ) / cl ) ); } } return f; } /******************************************************************************/ double *fxyz1_vec ( int n, double x[], double y[], double z[] ) /******************************************************************************/ /* Purpose: FXYZ1_VEC is the first 3D example, vector version. Discussion: This example is known as the 3D Shepp-Logan phantom. It should be plotted on [-1,+1] x [-1,+1.5] x [-1.5,+1.5]. Seventeen objects are modeled by ellipses of various gray levels, including: 1: Outer skull 2: Inner skull 3: Left eye 4: Right eye 5: Nose 6: Mouth 7: Left ear 8: Right ear 9: Left small tumor 10: Center small tumor 11: Right small tumor 12: Old f 13: Old g 14: Old e 15: Right ventricle 16: Left ventricle 17: Blood clot Licensing: This code is distributed under the MIT license. Modified: 14 February 2014 Author: John Burkardt Reference: Larry Shepp, Computerized tomography and nuclear magnetic resonance, Journal of Computer Assisted Tomography, Volume 4, Number 1, February 1980, pages 94-107. Parameters: Input, int N, the number of points. Input, double X[N], Y[N], Z[N], the arguments. Output, double FXYZ1_VEC[N], the function values. */ { double a1[17] = { 0.7233, 0.7008, 0.1270, 0.1270, 0.1270, 0.4575, 0.0635, 0.0635, 0.0460, 0.0230, 0.0230, 0.0460, 0.2100, 0.1100, 0.1600, 0.1600, 0.0300 }; double a2[17] = { 0.9644, 0.9246, 0.1270, 0.1270, 0.3400, 0.6099, 0.3175, 0.3175, 0.0230, 0.0230, 0.0460, 0.0460, 0.2581, 0.2500, 0.3100, 0.4100, 0.2000 }; double a3[17] = { 1.2700, 1.2241, 0.1270, 0.1270, 0.1700, 0.5080, 0.3175, 0.3175, 0.0230, 0.0460, 0.0230, 0.0460, 0.2581, 0.2300, 0.2540, 0.3810, 0.2000 }; double c; double *f; double g[17] = { 2.0000, -0.9800, -1.0000, -1.0000, 1.5000, -1.0000, 1.0000, 1.0000, 0.0100, 0.0100, 0.0100, 0.0100, 0.0100, 0.0100, -0.0200, -0.0200, 0.0300 }; int e; int i; double v11[17] = { 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9903, -0.9903, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9511, -0.9511, 0.9192 }; double v12[17] = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.1085, -0.1085, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.3090, -0.3090, -0.3381 }; double v13[17] = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.0865, -0.0865, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.2020 }; double v21[17] = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.1089, -0.1089, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.3090, -0.3090, 0.3452 }; double v22[17] = { 1.0000, 1.0000, 1.0000, 1.0000, 0.5446, 1.0000, 0.9941, 0.9941, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9511, 0.9511, 0.9385 }; double v23[17] = { 0.0000, 0.0000, 0.0000, 0.0000, -0.8387, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 }; double v31[17] = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0860, -0.0860, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.1896 }; double v32[17] = { 0.0000, 0.0000, 0.0000, 0.0000, 0.8387, 0.0000, -0.0094, -0.0094, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, -0.0697 }; double v33[17] = { 1.0000, 1.0000, 1.0000, 1.0000, 0.5446, 1.0000, 0.9963, 0.9963, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, -0.9794 }; double x0[17] = { 0.0000, 0.0000, 0.2583, -0.2583, 0.0000, 0.0000, 0.7076, -0.7076, -0.0800, 0.0000, 0.0600, 0.0000, 0.0000, 0.0000, 0.2200, -0.2200, 0.5600 }; double y0[17] = { 0.0000, -0.0184, 0.7534, 0.7534, 1.1398, 0.0000, -0.1378, -0.1378, -0.6050, -0.6050, -0.6050, 0.1000, -0.1000, 0.3500, 0.0000, 0.0000, -0.4000 }; double z0[17] = { 0.0000, -0.0185, 0.0000, 0.0000, -0.1957, -0.7620, -0.1905, -0.1905, 0.3810, 0.3810, 0.3810, 0.3810, 0.1270, 0.3810, 0.3810, 0.3810, 0.3810 }; f = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { f[i] = 0.0; for ( e = 0; e < 17; e++ ) { c = pow ( ( ( x[i] - x0[e] ) * v11[e] + ( y[i] - y0[e] ) * v12[e] + ( z[i] - z0[e] ) * v13[e] ) / a1[e], 2 ) + pow ( ( ( x[i] - x0[e] ) * v21[e] + ( y[i] - y0[e] ) * v22[e] + ( z[i] - z0[e] ) * v23[e] ) / a2[e], 2 ) + pow ( ( ( x[i] - x0[e] ) * v31[e] + ( y[i] - y0[e] ) * v32[e] + ( z[i] - z0[e] ) * v33[e] ) / a3[e], 2 ); if ( c <= 1.0 ) { f[i] = f[i] + g[e]; } } } return f; } /******************************************************************************/ double *r8vec_copy_new ( int n, double a1[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY_NEW copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Output, double R8VEC_COPY_NEW[N], the copy of A1. */ { double *a2; int i; a2 = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return a2; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: R8VEC_LINSPACE_NEW creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A, B, the first and last entries. Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ double *r8vec_uniform_ab_new ( int n, double a, double b, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_AB_NEW returns a scaled pseudorandom R8VEC. Discussion: Each dimension ranges from A to B. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 30 January 2005 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int N, the number of entries in the vector. Input, double A, B, the lower and upper limits of the pseudorandom values. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_AB_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_AB_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = a + ( b - a ) * ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }