# include # include # include # include # include # include # include "circle_segment.h" /******************************************************************************/ double circle_segment_angle_from_chord ( double r, double c[2], double p1[2], double p2[2] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_ANGLE_FROM_CHORD computes the angle of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the ends of the chord. Output, double CIRCLE_SEGMENT_ANGLE_FROM_CHORD, the value of THETA, the angle of the circle segment. 0 <= THETA < 2 * PI. */ { const double pi = 3.141592653589793; double theta; double v1[2]; double v2[2]; /* Compute the radial vectors V1 and V2. */ v1[0] = p1[0] - c[0]; v1[1] = p1[1] - c[1]; v2[0] = p2[0] - c[0]; v2[1] = p2[1] - c[1]; /* The arc cosine will only give us an answer between 0 and PI. */ theta = r8_atan ( v2[1], v2[0] ) - r8_atan ( v1[1], v1[0] ); /* Force 0 <= THETA < 2 * PI. */ while ( theta < 0.0 ) { theta = theta + 2.0 * pi; } while ( 2.0 * pi <= theta ) { theta = theta - 2.0 * pi; } return theta; } /******************************************************************************/ double circle_segment_angle_from_chord_angles ( double omega1, double omega2 ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_ANGLE_FROM_CHORD_ANGLES computes angle of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double OMEGA1, OMEGA2, the angles of the points P1 and P2. OMEGA1 <= OMEGA2. Output, double CIRCLE_SEGMENT_ANGLE_FROM_CHORD_ANGLES, the angle THETA of the circle segment. Essentially, THETA = OMEGA2 - OMEGA1. */ { const double pi = 3.141592653589793; double theta; while ( omega2 < omega1 ) { omega2 = omega2 + 2.0 * pi; } theta = omega2 - omega1; return theta; } /******************************************************************************/ double circle_segment_angle_from_height ( double r, double h ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_ANGLE_FROM_HEIGHT computes the angle of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double H, the "height" of the circle segment. 0 <= H <= 2 * R. Output, double CIRCLE_SEGMENT_ANGLE_FROM_HEIGHT, the angle THETA of the circle segment. */ { const double pi = 3.141592653589793; double theta; if ( h <= 0.0 ) { theta = 0.0; } else if ( h <= r ) { theta = 2.0 * r8_acos ( ( r - h ) / r ); theta = 2.0 * r8_asin ( sqrt ( r * r - ( r - h ) * ( r - h ) ) / r ); } else if ( h <= 2.0 * r ) { theta = 2.0 * r8_acos ( ( r - h ) / r ); theta = 2.0 * r8_asin ( sqrt ( r * r - ( r - h ) * ( r - h ) ) / r ); theta = 2.0 * pi - theta; } else { theta = 2.0 * pi; } return theta; } /******************************************************************************/ double circle_segment_area_from_angle ( double r, double theta ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_AREA_FROM_ANGLE computes the area of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double THETA, the angle of the circle segment. Output, double CIRCLE_SEGMENT_AREA_FROM_ANGLE, the area of the circle segment. */ { double area; area = r * r * ( theta - sin ( theta ) ) / 2.0; return area; } /******************************************************************************/ double circle_segment_area_from_chord ( double r, double c[2], double p1[2], double p2[2] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_AREA_FROM_CHORD computes the area of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the ends of the chord. Output, double CIRCLE_SEGMENT_AREA_FROM_CHORD, the area of the circle segment. */ { double area; double theta; theta = circle_segment_angle_from_chord ( r, c, p1, p2 ); area = r * r * ( theta - sin ( theta ) ) / 2.0; return area; } /******************************************************************************/ double circle_segment_area_from_height ( double r, double h ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_AREA_FROM_HEIGHT computes the area of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double H, the height of the circle segment. 0 <= H <= 2 * R. Output, double CIRCLE_SEGMENT_AREA_FROM_HEIGHT, the area of the circle segment. */ { double area; const double pi = 3.141592653589793; double theta; if ( h <= 0.0 ) { area = 0.0; } else if ( h <= r ) { theta = 2.0 * r8_asin ( sqrt ( r * r - ( r - h ) * ( r - h ) ) / r ); area = r * r * ( theta - sin ( theta ) ) / 2.0; } else if ( h <= 2.0 * r ) { theta = 2.0 * r8_asin ( sqrt ( r * r - ( r - h ) * ( r - h ) ) / r ); theta = 2.0 * pi - theta; area = r * r * ( theta - sin ( theta ) ) / 2.0; } else { area = pi * r * r; } return area; } /******************************************************************************/ double circle_segment_area_from_sample ( double r, double c[2], double p1[2], double p2[2], int n, int *seed ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_AREA_FROM_SAMPLE computes the area of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the ends of the chord. Input, int N, the number of sample points. Input/output, int *SEED, a seed for the random number generator. Output, double CIRCLE_SEGMENT_AREA_FROM_SAMPLE, the area of the circle segment. */ { double *angle; double area; int i; int m; double omega1; double omega2; double p[2]; const double pi = 3.141592653589793; double *r2; double rmh; double *vdotp; double *x; double *y; /* Determine the angles of the chord endpoints. */ omega1 = r8_atan ( p1[1] - c[1], p1[0] - c[0] ); while ( omega1 < 0.0 ) { omega1 = omega1 + 2.0 * pi; } omega2 = r8_atan ( p2[1] - c[1], p2[0] - c[0] ); while ( omega2 < omega1 ) { omega2 = omega2 + 2.0 * pi; } /* Get N random points in the circle. To simplify angle measurement, take OMEGA1 as your smallest angle. That way, the check OMEGA1 <= ANGLE <= OMEGA2 will be legitimate. */ angle = r8vec_uniform_01_new ( n, seed ); for ( i = 0; i < n; i++ ) { angle[i] = omega1 + 2.0 * pi * angle[i]; } r2 = r8vec_uniform_01_new ( n, seed ); for ( i = 0; i < n; i++ ) { r2[i] = sqrt ( r2[i] ); } x = ( double * ) malloc ( n * sizeof ( double ) ); y = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { x[i] = c[0] + r2[i] * cos ( angle[i] ); y[i] = c[0] + r2[i] * sin ( angle[i] ); } /* Determine the vector that touches the circle segment base. */ p[0] = 0.5 * ( p1[0] + p2[0] ) - c[0]; p[1] = 0.5 * ( p1[1] + p2[1] ) - c[1]; rmh = sqrt ( p[0] * p[0] + p[1] * p[1] ); p[0] = p[0] / rmh; p[1] = p[1] / rmh; if ( pi < omega2 - omega1 ) { p[0] = - p[0]; p[1] = - p[1]; rmh = - rmh; } /* Compute the projection of each point onto P. */ vdotp = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { vdotp[i] = ( x[i] - c[0] ) * p[0] + ( y[i] - c[1] ) * p[1]; } /* Points in the segment lie in the sector, and project at least RMH onto P. */ m = 0; for ( i = 0; i < n; i++ ) { if ( omega1 < angle[i] && angle[i] < omega2 && rmh < vdotp[i] ) { m = m + 1; } } /* The area of the segment is its relative share of the circle area. */ area = pi * r * r * ( double ) ( m ) / ( double ) ( n ); free ( angle ); free ( r2 ); free ( vdotp ); free ( x ); free ( y ); return area; } /******************************************************************************/ double circle_segment_cdf ( double r, double h, double h2 ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_CDF computes a CDF related to a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Now, suppose we want to assign a cumulative density function or CDF based on a variable H2 which measures the height of the circle segment formed by an arbitrary point in the circle segment. CDF(H2) will measure the probability that a point drawn uniformly at random from the circle segment defines a (smaller) circle segment of height H2. If we can define this CDF, then we will be able to sample uniformly from the circle segment, since our "Y" value can be determined from H2, and our X value is chosen uniformly over the horizontal chord through (0,Y). Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double H, the "height" of the circle segment. 0 <= H <= 2 * R. Input, double H2, the "height" of the new circle segment defined by a given point in the circle segment. 0 <= H2 <= H. Output, double CDF, the cumulative density function for H2, the probability that a point chosen at random in the circle segment would define a smaller circle segment of height H2 or less. */ { double a; double a2; double cdf; if ( h2 <= 0.0 ) { cdf = 0.0; } else if ( h <= h2 ) { cdf = 1.0; } else { a = circle_segment_area_from_height ( r, h ); a2 = circle_segment_area_from_height ( r, h2 ); cdf = a2 / a; } return cdf; } /******************************************************************************/ double *circle_segment_centroid_from_chord ( double r, double c[2], double p1[2], double p2[2] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_CENTROID_FROM_CHORD computes the centroid of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. For this function, we assume that the center of the circle is at (0,0), that the chord is horizontal, and that the circle segment is at the top. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the coordinates of the endpoints of the chord. Output, double CIRCLE_SEGMENT_CENTROID_FROM_CHORD[2], the coordinates of the centroid. */ { double *d; double s; double theta; double thetah; double v1[2]; /* Get the angle subtended by P1:P2. */ theta = circle_segment_angle_from_chord ( r, c, p1, p2 ); /* Construct V1, the vector from C to P1. */ v1[0] = p1[0] - c[0]; v1[1] = p1[1] - c[1]; /* Rotate V1 through THETA / 2. */ thetah = theta / 2.0; d = ( double * ) malloc ( 2 * sizeof ( double ) ); d[0] = cos ( thetah ) * v1[0] - sin ( thetah ) * v1[1]; d[1] = sin ( thetah ) * v1[0] + cos ( thetah ) * v1[1]; /* Scale this vector so it represents the distance to the centroid relative to R. */ s = 4.0 * pow ( sin ( theta / 2.0 ), 3 ) / 3.0 / ( theta - sin ( theta ) ); d[0] = s * d[0]; d[1] = s * d[1]; /* Add the center. */ d[0] = d[0] + c[0]; d[1] = d[1] + c[1]; return d; } /******************************************************************************/ double *circle_segment_centroid_from_height ( double r, double h ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_CENTROID_FROM_HEIGHT computes centroid of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. For this function, we assume that the center of the circle is at (0,0). Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double H, the "height" of the circle segment. 0 <= H <= 2 * R. Output, double CIRCLE_SEGMENT_CENTROID_FROM_HEIGHT[2], the coordinates of the centroid. */ { double *d; double theta; theta = circle_segment_angle_from_height ( r, h ); d = ( double * ) malloc ( 2 * sizeof ( double ) ); d[0] = 0.0; d[1] = 4.0 * r * pow ( sin ( theta / 2.0 ), 3 ) / 3.0 / ( theta - sin ( theta ) ); return d; } /******************************************************************************/ double *circle_segment_centroid_from_sample ( double r, double c[2], double p1[2], double p2[2], int n, int *seed ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_CENTROID_FROM_SAMPLE estimates a circle segment centroid. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the ends of the chord. Input, int N, the number of sample points. Input/output, int *SEED, a seed for the random number generator. Output, double CIRCLE_SEGMENT_CENTROID_FROM_SAMPLE[2], the estimated centroid of the circle segment. */ { double *d; double *x; double *y; x = ( double * ) malloc ( n * sizeof ( double ) ); y = ( double * ) malloc ( n * sizeof ( double ) ); circle_segment_sample_from_chord ( r, c, p1, p2, n, seed, x, y ); d = ( double * ) malloc ( 2 * sizeof ( double ) ); d[0] = r8vec_sum ( n, x ) / ( double ) ( n ); d[1] = r8vec_sum ( n, y ) / ( double ) ( n ); free ( x ); free ( y ); return d; } /******************************************************************************/ int circle_segment_contains_point ( double r, double c[2], double omega1, double omega2, double xy[2] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_CONTAINS_POINT reports whether a point is in a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. In this function, we allow the circle to have an arbitrary center C, arbitrary radius R, and we describe the points P1 and P2 by specifying their angles OMEGA1 and OMEGA2. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double OMEGA1, OMEGA2, the angles of the two points on the circumference of the circle that define the circle segment. OMEGA1 < OMEGA2 <= OMEGA1 + 2 * PI Input, double XY[2], a point. Output, int CIRCLE_SEGMENT_CONTAINS_POINT, is TRUE if the point is inside the circle segment. */ { double h; double omegah; const double pi = 3.141592653589793; double theta; double v[2]; double v_omega; double v_project; double v_r; int value; if ( r <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_CONTAINS_POINT - Fatal error!\n" ); fprintf ( stderr, " R <= 0.0.\n" ); exit ( 1 ); } while ( omega2 < omega1 ) { omega2 = omega2 + 2.0 * pi; } /* Compute the vector V = XY - C: */ v[0] = xy[0] - c[0]; v[1] = xy[1] - c[1]; /* a: Point must be inside the circle, so ||V|| <= R. */ v_r = sqrt ( v[0] * v[0] + v[1] * v[1] ); if ( r < v_r ) { value = 0; return value; } /* b: Angle made by the vector V must be between OMEGA1 and OMEGA2. */ v_omega = atan2 ( v[1], v[0] ); while ( omega1 <= v_omega + 2.0 * pi ) { v_omega = v_omega - 2.0 * pi; } while ( v_omega + 2.0 * pi <= omega1 ) { v_omega = v_omega + 2.0 * pi; } if ( omega2 < v_omega ) { value = 0; return value; } /* c: Projection of V onto unit centerline must be at least R-H. */ omegah = 0.5 * ( omega1 + omega2 ); v_project = v[0] * cos ( omegah ) + v[1] * sin ( omegah ); theta = omega2 - omega1; h = circle_segment_height_from_angle ( r, theta ); if ( v_project < r - h ) { value = 0; return value; } value = 1; return value; } /******************************************************************************/ double circle_segment_height_from_angle ( double r, double angle ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_HEIGHT_FROM_ANGLE: height of a circle segment from angle. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. This function is given the radius R and angle of the segment, and determines the corresponding height. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double ANGLE, the angle of the circle segment. 0 <= ANGLE <= 2.0 * PI. Output, double CIRCLE_SEGMENT_HEIGHT_FROM_ANGLE, the height of the circle segment. */ { double h; const double pi = 3.141592653589793; if ( angle < 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_HEIGHT_FROM_ANGLE - Fatal error!\n" ); fprintf ( stderr, " ANGLE < 0.0.\n" ); exit ( 1 ); } if ( angle == 0.0 ) { h = 0.0; return h; } if ( angle == 2.0 * pi ) { h = 2.0 * r; return h; } if ( 2.0 * pi < angle ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_HEIGHT_FROM_ANGLE - Fatal error!\n" ); fprintf ( stderr, " 2.0 * pi < ANGLE.\n" ); exit ( 1 ); } if ( r <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_HEIGHT_FROM_ANGLE - Fatal error!\n" ); fprintf ( stderr, " R <= 0.0.\n" ); exit ( 1 ); } if ( angle <= pi ) { h = r * ( 1.0 - cos ( angle / 2.0 ) ); } else { h = r * ( 1.0 + cos ( ( 2.0 * pi - angle ) / 2.0 ) ); } return h; } /******************************************************************************/ double circle_segment_height_from_area ( double r, double area ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_HEIGHT_FROM_AREA: height of a circle segment from area. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. This function is given the radius R and area of the segment, and determines the corresponding height. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double AREA, the area of the circle segment. 0 <= AREA <= 2.0 * PI * R^2. Output, double CIRCLE_SEGMENT_HEIGHT_FROM_AREA, the height of the circle segment. */ { double a; double area_circle; double eps; double h; double h1; double h2; int it; const double pi = 3.141592653589793; if ( area < 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!\n" ); fprintf ( stderr, " AREA < 0.0.\n" ); exit ( 1 ); } area_circle = 2.0 * pi * r * r; if ( area == 0.0 ) { h = 0.0; return h; } if ( area == area_circle ) { h = 2.0 * r; return h; } if ( area_circle < area ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error!\n" ); fprintf ( stderr, " 2.0 * pi * r^2 < AREA.\n" ); exit ( 1 ); } if ( r <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CIRCLE_SEGMENT_HEIGHT_FROM_AREA - Fatal error\n" ); fprintf ( stderr, " R <= 0.0.\n" ); exit ( 1 ); } h1 = 0.0; h2 = 2.0 * r; it = 0; eps = DBL_EPSILON; while ( it < 30 ) { h = 0.5 * ( h1 + h2 ); a = circle_segment_area_from_height ( r, h ); it = it + 1; if ( fabs ( a - area ) < sqrt ( eps ) * area ) { break; } if ( a < area ) { h1 = h; } else { h2 = h; } } return h; } /******************************************************************************/ double circle_segment_height_from_chord ( double r, double c[2], double p1[2], double p2[2] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_HEIGHT_FROM_CHORD: height of a circle segment from chord. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the coordinates of the circle center. Input, double P1[2], P2[2], the coordinates of the chord endpoints. Output, double CIRCLE_SEGMENT_HEIGHT_FROM_CHORD, the height of the circle segment. */ { double h; double theta; theta = circle_segment_angle_from_chord ( r, c, p1, p2 ); h = circle_segment_height_from_angle ( r, theta ); return h; } /******************************************************************************/ double circle_segment_rotation_from_chord ( double r, double c[], double p1[], double p2[] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_ROTATION_FROM_CHORD computes the rotation of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 16 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the ends of the chord. Warning! If P1 = P2, we can't tell whether the segment is the whole circle or none of it! Output, double CIRCLE_SEGMENT_ROTATION_FROM_CHORD, the rotation of the circle segment. 0 <= ALPHA < 2 * PI. */ { double alpha; double pi = 3.141592653589793; double rho1; double rho2; double theta; double v1[2]; double v2[2]; /* Compute the radial vectors V1 and V2. */ v1[0] = p1[0] - c[0]; v1[1] = p1[1] - c[1]; v2[0] = p2[0] - c[0]; v2[1] = p2[1] - c[1]; /* Use R8_ATAN to guarantee that 0 <= RHO1, RHO2 <= 2 * PI. */ rho1 = r8_atan ( v1[1], v1[0] ); rho2 = r8_atan ( v2[1], v2[0] ); /* Force RHO2 to be bigger than RHO1. */ while ( rho2 <= rho1 ) { rho2 = rho2 + 2.0 * pi; } /* Compute THETA. */ theta = rho2 - rho1; /* ALPHA is RHO1, plus half of the angular distance between P1 and P2. */ alpha = rho1 + 0.5 * theta; while ( 2.0 * pi <= alpha ) { alpha = alpha - 2.0 * pi; } return alpha; } /******************************************************************************/ void circle_segment_sample_from_chord ( double r, double c[2], double p1[2], double p2[2], int n, int *seed, double x[], double y[] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_SAMPLE_FROM_CHORD samples points from a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double C[2], the center of the circle. Input, double P1[2], P2[2], the endpoints of the chord. Input, int N, the number of sample points. Input/output, int *SEED, a seed for the random number generator. Output, double X[N], Y[N], the sample points. */ { double c2[2]; double *eta; double h; int i; double t; double vc[2]; double vr[2]; double *xi; /* Determine unit vectors VR and VC. VR points to the center of the chord from the radius. VC points along the chord, from P1 to P2. */ vr[0] = 0.5 * ( p1[0] + p2[0] ) - c[0]; vr[1] = 0.5 * ( p1[1] + p2[1] ) - c[1]; t = sqrt ( vr[0] * vr[0] + vr[1] * vr[1] ); vr[0] = vr[0] / t; vr[1] = vr[1] / t; vc[0] = p2[0] - p1[0]; vc[1] = p2[1] - p1[1]; t = sqrt ( vc[0] * vc[0] + vc[1] * vc[1] ); vc[0] = vc[0] / t; vc[1] = vc[1] / t; /* Get the height of the circle segment. */ c2[0] = 0.0; c2[1] = 0.0; h = circle_segment_height_from_chord ( r, c2, p1, p2 ); /* Sample (xi,eta) in the reference coordinates, where the chord is horizontal. */ xi = ( double * ) malloc ( n * sizeof ( double ) ); eta = ( double * ) malloc ( n * sizeof ( double ) ); circle_segment_sample_from_height ( r, h, n, seed, xi, eta ); /* XI is the left/right coordinate along VC. ETA is the distance along VR. */ for ( i = 0; i < n; i++ ) { x[i] = c[0] + eta[i] * vr[0] + xi[i] * vc[0]; y[i] = c[1] + eta[i] * vr[1] + xi[i] * vc[1]; } free ( eta ); free ( xi ); return; } /******************************************************************************/ void circle_segment_sample_from_height ( double r, double h, int n, int *seed, double x[], double y[] ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_SAMPLE_FROM_HEIGHT samples points from a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double H, the height of the circle segment. 0 <= H <= 2 * R. Input, int N, the number of sample points. Input/output, int *SEED, a seed for the random number generator. Output, double X[N], Y[N], the sample points. */ { double area; double *area2; double *h2; int i; double *u; double *wh; area = circle_segment_area_from_height ( r, h ); /* Pick CDF's randomly. */ u = r8vec_uniform_01_new ( n, seed ); /* Choose points randomly by choosing ordered areas randomly. */ area2 = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { area2[i] = u[i] * area; } /* Each area corresponds to a height H2. Find it. */ h2 = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { h2[i] = circle_segment_height_from_area ( r, area2[i] ); } /* Determine the half-width WH of the segment for each H2. */ wh = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { wh[i] = sqrt ( h2[i] * ( 2.0 * r - h2[i] ) ); } /* Choose an X position randomly in [-WH,+WH]. */ free ( u ); u = r8vec_uniform_01_new ( n, seed ); for ( i = 0; i < n; i++ ) { x[i] = ( 2.0 * u[i] - 1.0 ) * wh[i]; } /* Our circle center is at (0,0). Our height of H2 is subtracted from the height R at the peak of the circle. Determine the Y coordinate using this fact. */ for ( i = 0; i < n; i++ ) { y[i] = r - h2[i]; } free ( area2 ); free ( h2 ); free ( u ); free ( wh ); return; } /******************************************************************************/ double circle_segment_width_from_height ( double r, double h ) /******************************************************************************/ /* Purpose: CIRCLE_SEGMENT_WIDTH_FROM_HEIGHT computes the width of a circle segment. Discussion: Begin with a circle of radius R. Choose two points P1 and P2 on the circle, and draw the chord P1:P2. This chord divides the circle into two pieces, each of which is called a circle segment. Consider one of the pieces. The "angle" of this segment is the angle P1:C:P2, where C is the center of the circle. Let Q be the point on the chord P1:P2 which is closest to C. The "height" of the segment is the distance from Q to the perimeter of the circle. The "width" of the circle segment is the length of P1:P2. This function is given the radius R and height H of the segment, and determines the corresponding width W. Licensing: This code is distributed under the MIT license. Modified: 06 July 2013 Author: John Burkardt Parameters: Input, double R, the radius of the circle. 0 < R. Input, double H, the height of the circle segment. 0 <= H <= 2 * R. Output, double CIRCLE_SEGMENT_WIDTH_FROM_HEIGHT, the width of the circle segment. */ { double w; w = 2.0 * sqrt ( h * ( 2.0 * r - h ) ); return w; } /******************************************************************************/ void filename_inc ( char *filename ) /******************************************************************************/ /* Purpose: FILENAME_INC increments a partially numeric file name. Discussion: It is assumed that the digits in the name, whether scattered or connected, represent a number that is to be increased by 1 on each call. If this number is all 9's on input, the output number is all 0's. Non-numeric letters of the name are unaffected. If the name is empty, then the routine stops. If the name contains no digits, the empty string is returned. Example: Input Output ----- ------ "a7to11.txt" "a7to12.txt" (typical case. Last digit incremented) "a7to99.txt" "a8to00.txt" (last digit incremented, with carry.) "a9to99.txt" "a0to00.txt" (wrap around) "cat.txt" " " (no digits to increment) " " STOP (error) Licensing: This code is distributed under the MIT license. Modified: 22 November 2011 Author: John Burkardt Parameters: Input/output, char *FILENAME, the filename to be incremented. */ { int change; int n; char *t; n = s_len_trim ( filename ); if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILENAME_INC - Fatal error!\n" ); fprintf ( stderr, " The input string is empty.\n" ); exit ( 1 ); } change = 0; t = filename + n - 1; while ( 0 < n ) { if ( '0' <= *t && *t <= '9' ) { change = change + 1; if ( *t == '9' ) { *t = '0'; } else { *t = *t + 1; return; } } t--; n--; } /* No digits were found. Return blank. */ if ( change == 0 ) { n = s_len_trim ( filename ); t = filename + n - 1; while ( 0 < n ) { *t = ' '; t--; n--; } } return; } /******************************************************************************/ void gauss ( int n, double alpha[], double beta[], double x[], double w[] ) /******************************************************************************/ /* Purpose: GAUSS computes a Gauss quadrature rule. Discussion: Given a weight function W encoded by the first N recurrence coefficients ALPHA and BETA for the associated orthogonal polynomials, the call call gauss ( n, alpha, beta, x, w ) generates the nodes and weights of the N-point Gauss quadrature rule for the weight function W. Licensing: This code is distributed under the MIT license. Modified: 16 July 2013 Author: Original MATLAB version by Walter Gautschi. C version by John Burkardt. Reference: Walter Gautschi, Orthogonal Polynomials: Computation and Approximation, Oxford, 2004, ISBN: 0-19-850672-4, LC: QA404.5 G3555. Parameters: Input, int N, the order of the desired quadrature rule. Input, double ALPHA[N], BETA[N], the alpha and beta recurrence coefficients for the othogonal polynomials associated with the weight function. Output, double X[N], W[N], the nodes and weights of the desired quadrature rule. The nodes are listed in increasing order. */ { double *a; int i; int it_max; int it_num; int j; int rot_num; double *v; /* Define the tridiagonal Jacobi matrix. */ a = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == j ) { a[i+j*n] = alpha[i]; } else if ( i == j - 1 ) { a[i+j*n] = sqrt ( beta[j] ); } else if ( i - 1 == j ) { a[i+j*n] = sqrt ( beta[i] ); } else { a[i+j*n] = 0.0; } } } /* Get the eigenvectors and eigenvalues. */ it_max = 100; v = ( double * ) malloc ( n * n * sizeof ( double ) ); jacobi_eigenvalue ( n, a, it_max, v, x, &it_num, &rot_num ); for ( j = 0; j < n; j++ ) { w[j] = beta[0] * v[0+j*n] * v[0+j*n]; } free ( a ); free ( v ); return; } /******************************************************************************/ int i4vec_sum ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SUM sums the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: A = ( 1, 2, 3, 4 ) Output: I4VEC_SUM = 10 Licensing: This code is distributed under the MIT license. Modified: 29 May 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int A[N], the vector to be summed. Output, int I4VEC_SUM, the sum of the entries of A. */ { int i; int sum; sum = 0; for ( i = 0; i < n; i++ ) { sum = sum + a[i]; } return sum; } /******************************************************************************/ void jacobi_eigenvalue ( int n, double a[], int it_max, double v[], double d[], int *it_num, int *rot_num ) /******************************************************************************/ /* Purpose: JACOBI_EIGENVALUE carries out the Jacobi eigenvalue iteration. Discussion: This function computes the eigenvalues and eigenvectors of a real symmetric matrix, using Rutishauser's modfications of the classical Jacobi rotation method with threshold pivoting. Licensing: This code is distributed under the MIT license. Modified: 15 July 2013 Author: C version by John Burkardt Parameters: Input, int N, the order of the matrix. Input, double A[N*N], the matrix, which must be square, real, and symmetric. Input, int IT_MAX, the maximum number of iterations. Output, double V[N*N], the matrix of eigenvectors. Output, double D[N], the eigenvalues, in descending order. Output, int *IT_NUM, the total number of iterations. Output, int *ROT_NUM, the total number of rotations. */ { double *bw; double c; double g; double gapq; double h; int i; int j; int k; int l; int m; int p; int q; double s; double t; double tau; double term; double termp; double termq; double theta; double thresh; double w; double *zw; for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == j ) { v[i+j*n] = 1.0; } else { v[i+j*n] = 0.0; } } } for ( i = 0; i < n; i++ ) { d[i] = a[i+i*n]; } bw = ( double * ) malloc ( n * sizeof ( double ) ); zw = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { bw[i] = d[i]; zw[i] = 0.0; } *it_num = 0; *rot_num = 0; while ( *it_num < it_max ) { *it_num = *it_num + 1; /* The convergence threshold is based on the size of the elements in the strict upper triangle of the matrix. */ thresh = 0.0; for ( j = 0; j < n; j++ ) { for ( i = 0; i < j; i++ ) { thresh = thresh + a[i+j*n] * a[i+j*n]; } } thresh = sqrt ( thresh ) / ( double ) ( 4 * n ); if ( thresh == 0.0 ) { break; } for ( p = 0; p < n; p++ ) { for ( q = p + 1; q < n; q++ ) { gapq = 10.0 * fabs ( a[p+q*n] ); termp = gapq + fabs ( d[p] ); termq = gapq + fabs ( d[q] ); /* Annihilate tiny offdiagonal elements. */ if ( 4 < *it_num && termp == fabs ( d[p] ) && termq == fabs ( d[q] ) ) { a[p+q*n] = 0.0; } /* Otherwise, apply a rotation. */ else if ( thresh <= fabs ( a[p+q*n] ) ) { h = d[q] - d[p]; term = fabs ( h ) + gapq; if ( term == fabs ( h ) ) { t = a[p+q*n] / h; } else { theta = 0.5 * h / a[p+q*n]; t = 1.0 / ( fabs ( theta ) + sqrt ( 1.0 + theta * theta ) ); if ( theta < 0.0 ) { t = - t; } } c = 1.0 / sqrt ( 1.0 + t * t ); s = t * c; tau = s / ( 1.0 + c ); h = t * a[p+q*n]; /* Accumulate corrections to diagonal elements. */ zw[p] = zw[p] - h; zw[q] = zw[q] + h; d[p] = d[p] - h; d[q] = d[q] + h; a[p+q*n] = 0.0; /* Rotate, using information from the upper triangle of A only. */ for ( j = 0; j < p; j++ ) { g = a[j+p*n]; h = a[j+q*n]; a[j+p*n] = g - s * ( h + g * tau ); a[j+q*n] = h + s * ( g - h * tau ); } for ( j = p + 1; j < q; j++ ) { g = a[p+j*n]; h = a[j+q*n]; a[p+j*n] = g - s * ( h + g * tau ); a[j+q*n] = h + s * ( g - h * tau ); } for ( j = q + 1; j < n; j++ ) { g = a[p+j*n]; h = a[q+j*n]; a[p+j*n] = g - s * ( h + g * tau ); a[q+j*n] = h + s * ( g - h * tau ); } /* Accumulate information in the eigenvector matrix. */ for ( j = 0; j < n; j++ ) { g = v[j+p*n]; h = v[j+q*n]; v[j+p*n] = g - s * ( h + g * tau ); v[j+q*n] = h + s * ( g - h * tau ); } *rot_num = *rot_num + 1; } } } for ( i = 0; i < n; i++ ) { bw[i] = bw[i] + zw[i]; d[i] = bw[i]; zw[i] = 0.0; } } /* Restore upper triangle of input matrix. */ for ( j = 0; j < n; j++ ) { for ( i = 0; i < j; i++ ) { a[i+j*n] = a[j+i*n]; } } /* Ascending sort the eigenvalues and eigenvectors. */ for ( k = 0; k < n - 1; k++ ) { m = k; for ( l = k + 1; l < n; l++ ) { if ( d[l] < d[m] ) { m = l; } } if ( m != k ) { t = d[m]; d[m] = d[k]; d[k] = t; for ( i = 0; i < n; i++ ) { w = v[i+m*n]; v[i+m*n] = v[i+k*n]; v[i+k*n] = w; } } } free ( bw ); free ( zw ); return; } /******************************************************************************/ void r_jacobi ( int n, double a, double b, double alpha[], double beta[] ) /******************************************************************************/ /* Purpose: R_JACOBI computes recurrence coefficients for monic Jacobi polynomials. Discussion: This function generates the first N recurrence coefficients for monic Jacobi polynomials with parameters A and B. These polynomials are orthogonal on [-1,1] relative to the weight w(x) = (1.0-x)^A * (1.0+x)^B. Licensing: This code is distributed under the MIT license. Modified: 17 July 2013 Author: Original MATLAB version by Dirk Laurie, Walter Gautschi. C version by John Burkardt. Reference: Walter Gautschi, Orthogonal Polynomials: Computation and Approximation, Oxford, 2004, ISBN: 0-19-850672-4, LC: QA404.5 G3555. Parameters: Input, int N, the number of coefficients desired. Input, double A, B, the parameters for the Jacobi polynomial. -1.0 < A, -1.0 < B. Output, double ALPHA[N], BETA[N], the first N recurrence coefficients. */ { int i; double i_r8; double mu; double nab; double nu; if ( a <= -1.0 ) { printf ( "\n" ); printf ( "R_JACOBI - Fatal error!\n" ); printf ( " Illegal value of A.\n" ); exit ( 1 ); } if ( b <= -1.0 ) { printf ( "\n" ); printf ( "R_JACOBI - Fatal error!\n" ); printf ( " Illegal value of B.\n" ); exit ( 1 ); } nu = ( b - a ) / ( a + b + 2.0 ); mu = pow ( 2.0, a + b + 1.0 ) * tgamma ( a + 1.0 ) * tgamma ( b + 1.0 ) / tgamma ( a + b + 2.0 ); alpha[0] = nu; beta[0] = mu; if ( n == 1 ) { return; } for ( i = 1; i < n; i++ ) { i_r8 = ( double ) ( i + 1 ); alpha[i] = ( b - a ) * ( b + a ) / ( 2.0 * ( i_r8 - 1.0 ) + a + b ) / ( 2.0 * i_r8 + a + b ); } beta[1] = 4.0 * ( a + 1.0 ) * ( b + 1.0 ) / ( a + b + 2.0 ) / ( a + b + 2.0 ) / ( a + b + 3.0 ); for ( i = 2; i < n; i++ ) { i_r8 = ( double ) ( i + 1 ); nab = 2.0 * ( i_r8 - 1.0 ) + a + b; beta[i] = 4.0 * ( i_r8 - 1.0 + a ) * ( i_r8 - 1.0 + b ) * ( i_r8 - 1.0 ) * ( i_r8 - 1.0 + a + b ) / nab / nab / ( nab + 1.0 ) / ( nab - 1.0 ); } return; } /******************************************************************************/ double r8_acos ( double c ) /******************************************************************************/ /* Purpose: R8_ACOS computes the arc cosine function, with argument truncation. Discussion: If you call your system ACOS routine with an input argument that is outside the range [-1.0, 1.0 ], you may get an unpleasant surprise. This routine truncates arguments outside the range. Licensing: This code is distributed under the MIT license. Modified: 13 June 2002 Author: John Burkardt Parameters: Input, double C, the argument, the cosine of an angle. Output, double R8_ACOS, an angle whose cosine is C. */ { double angle; const double pi = 3.141592653589793; if ( c <= -1.0 ) { angle = pi; } else if ( 1.0 <= c ) { angle = 0.0; } else { angle = acos ( c ); } return angle; } /******************************************************************************/ double r8_asin ( double s ) /******************************************************************************/ /* Purpose: R8_ASIN computes the arc sine function, with argument truncation. Discussion: If you call your system ASIN routine with an input argument that is outside the range [-1.0, 1.0 ], you may get an unpleasant surprise. In particular, you may get the value NaN returned. This routine truncates arguments outside the range, avoiding the problem. Licensing: This code is distributed under the MIT license. Modified: 04 March 2008 Author: John Burkardt Parameters: Input, double S, the argument. Output, double R8_ASIN, an angle whose sine is S. */ { double value; s = fmax ( s, -1.0 ); s = fmin ( s, +1.0 ); value = asin ( s ); return value; } /******************************************************************************/ double r8_atan ( double y, double x ) /******************************************************************************/ /* Purpose: R8_ATAN computes the inverse tangent of the ratio Y / X. Discussion: R8_ATAN returns an angle whose tangent is ( Y / X ), a job which the built in functions ATAN and ATAN2 already do. However: * R8_ATAN always returns a positive angle, between 0 and 2 PI, while ATAN and ATAN2 return angles in the interval [-PI/2,+PI/2] and [-PI,+PI] respectively; * R8_ATAN accounts for the signs of X and Y, (as does ATAN2). The ATAN function by contrast always returns an angle in the first or fourth quadrants. Licensing: This code is distributed under the MIT license. Modified: 09 May 2003 Author: John Burkardt Parameters: Input, double Y, X, two quantities which represent the tangent of an angle. If Y is not zero, then the tangent is (Y/X). Output, double R8_ATAN, an angle between 0 and 2 * PI, whose tangent is (Y/X), and which lies in the appropriate quadrant so that the signs of its cosine and sine match those of X and Y. */ { double abs_x; double abs_y; const double pi = 3.141592653589793; double theta; double theta_0; /* Special cases: */ if ( x == 0.0 ) { if ( 0.0 < y ) { theta = pi / 2.0; } else if ( y < 0.0 ) { theta = 3.0 * pi / 2.0; } else if ( y == 0.0 ) { theta = 0.0; } } else if ( y == 0.0 ) { if ( 0.0 < x ) { theta = 0.0; } else if ( x < 0.0 ) { theta = pi; } } /* We assume that ATAN2 is correct when both arguments are positive. */ else { abs_y = fabs ( y ); abs_x = fabs ( x ); theta_0 = atan2 ( abs_y, abs_x ); if ( 0.0 < x && 0.0 < y ) { theta = theta_0; } else if ( x < 0.0 && 0.0 < y ) { theta = pi - theta_0; } else if ( x < 0.0 && y < 0.0 ) { theta = pi + theta_0; } else if ( 0.0 < x && y < 0.0 ) { theta = 2.0 * pi - theta_0; } } return theta; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 26 June 2013 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14f", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8mat_uniform_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_UNIFORM_01_NEW fills an R8MAT with pseudorandom values scaled to [0,1]. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. 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 June 2009 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns. Input/output, int *SEED, the "seed" value. Normally, this value should not be 0, otherwise the output value of SEED will still be 0, and R8_UNIFORM will be 0. On output, SEED has been updated. Output, double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; int k; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10; } } return r; } /******************************************************************************/ 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_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of 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 vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: 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: 19 August 2004 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/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_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] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Licensing: This code is distributed under the MIT license. Modified: 26 April 2003 Author: John Burkardt Parameters: Input, char *S, a pointer to a string. Output, int S_LEN_TRIM, the length of the string to the last nonblank. If S_LEN_TRIM is 0, then the string is entirely blank. */ { int n; char *t; n = strlen ( s ); t = s + strlen ( s ) - 1; while ( 0 < n ) { if ( *t != ' ' ) { return n; } t--; n--; } return n; } /******************************************************************************/ void timestamp ( void ) /******************************************************************************/ /* 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 } /******************************************************************************/ double *tridisolve ( int n, double a[], double b[], double c[], double d[] ) /******************************************************************************/ /* Purpose: TRIDISOLVE solves a tridiagonal system of linear equations. Discussion: We can describe an NxN tridiagonal matrix by vectors A, B, and C, where A and C are of length N-1. In that case, a linear system can be represented as b(1) * x(1) + c(1) * x(2) = d(1), a(j-1) * x(j-1) + b(j) * x(j) + c(j) * x(j+1) = d(j), j = 2:n-1, a(n-1) * x(n-1) + b(n) * x(n) = d(n) This function produces the solution vector X. This function is derived from Cleve Moler's Matlab suite. Licensing: This code is distributed under the MIT license. Modified: 19 May 2013 Author: John Burkardt. Parameters: Input, int N, the order of the linear system. Input, double A(N-1), B(N), C(N-1), the matrix entries. Input, double D(N), the right hand side. Output, double TRIDISOLVE[N], the solution. */ { double *bi; int j; double mu; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { x[j] = d[j]; } bi = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { bi[j] = 1.0 / b[j]; } for ( j = 0; j < n - 1; j++ ) { mu = a[j] * bi[j]; b[j+1] = b[j+1] - mu * c[j]; x[j+1] = x[j+1] - mu * x[j]; } x[n-1] = x[n-1] * bi[n-1]; for ( j = n - 2; 0 <= j; j-- ) { x[j] = ( x[j] - c[j] * x[j+1] ) * bi[j]; } free ( bi ); return x; }