C 代码 执行椭圆的几何计算，包括 面积、到点的距离、偏心率、周长、沿 周长、随机抽样.rar

# include <float.h> # include <math.h> # include <stdbool.h> # include <stdio.h> # include <stdlib.h> # include "ellipse.h" /******************************************************************************/ double ellipse_area1 ( double a[], double r ) /******************************************************************************/ /* Purpose: ellipse_area1() returns the area of an ellipse in quadratic form. Discussion: The points X in the ellipse are described by a 2 by 2 positive definite symmetric matrix A, and a "radius" R, such that X' * A * X <= R * R Licensing: This code is distributed under the GNU LGPL license. Modified: 13 August 2014 Author: John Burkardt Input: double A[2*2], the matrix that describes the ellipse. A must be symmetric positive definite. double R, the "radius" of the ellipse. Output: double ELLIPSE_AREA1, the area of the ellipse. */ { const double r8_pi = 3.141592653589793; double value; value = r * r * r8_pi / sqrt ( ( a[0+0*2] * a[1+1*2] - a[1+0*2] * a[0+1*2] ) ); return value; } /******************************************************************************/ double ellipse_area2 ( double a, double b, double c, double d ) /******************************************************************************/ /* Purpose: ellipse_area2() returns the area of an ellipse defined by an equation. Discussion: The ellipse is described by the formula a x^2 + b xy + c y^2 = d Licensing: This code is distributed under the GNU LGPL license. Modified: 08 November 2016 Author: John Burkardt Input: double A, B, C, coefficients on the left hand side. double D, the right hand side. Output: double ELLIPSE_AREA2, the area of the ellipse. */ { const double r8_pi = 3.141592653589793; double value; value = 2.0 * d * d * r8_pi / sqrt ( 4.0 * a * c - b * b ); return value; } /******************************************************************************/ double ellipse_area3 ( double r1, double r2 ) /******************************************************************************/ /* Purpose: ellipse_area3() returns the area of an ellipse in standard form. Discussion: An ellipse in standard form has a center at the origin, and axes aligned with the coordinate axes. Any point P on the ellipse satisfies ( X / A )^2 + ( Y / B )^2 == 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 17 May 2010 Author: John Burkardt Input: double A, B, the lengths of the major and minor semi-axes. Output: double ELLIPSE_AREA3, the area of the ellipse. */ { double area; const double r8_pi = 3.141592653589793; area = r8_pi * r1 * r2; return area; } /******************************************************************************/ double ellipse_eccentricity ( double a, double b ) /******************************************************************************/ /* Purpose: ellipse_eccentricity() computes the eccentricity of an ellipse. Discussion: An ellipse in standard form has a center at the origin, and axes aligned with the coordinate axes. Any point P on the ellipse satisfies ( X / A )^2 + ( Y / B )^2 == 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 24 March 2021 Author: John Burkardt Input: double A, B, the major and minor semi-axes. Output: double ELLIPSE_ECCENTRICITY, the eccentricity of the ellipse. */ { double ecc; double t; a = fabs ( a ); b = fabs ( b ); if ( a < b ) { t = a; a = b; b = t; } ecc = 1.0 - pow ( b / a, 2 ); return ecc; } /******************************************************************************/ double ellipse_point_dist_2d ( double r1, double r2, double p[2] ) /******************************************************************************/ /* Purpose: ellipse_point_dist_2d() finds the distance from a point to an ellipse in 2D. Discussion: An ellipse in standard form has a center at the origin, and axes aligned with the coordinate axes. Any point P on the ellipse satisfies ( X / A )^2 + ( Y / B )^2 == 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 17 May 2010 Author: John Burkardt Reference: Dianne O'Leary, Elastoplastic Torsion: Twist and Stress, Computing in Science and Engineering, July/August 2004, pages 74-76. September/October 2004, pages 63-65. Input: double R1, R2, the ellipse parameters. Normally, these are both positive quantities. Generally, they are also distinct. double P[2], the point. Output: double ELLIPSE_POINT_DIST_2D, the distance to the ellipse. */ { double dist; int i; double *pn; pn = ellipse_point_near_2d ( r1, r2, p ); dist = 0.0; for ( i = 0; i < 2; i++ ) { dist = dist + pow ( p[i] - pn[i], 2 ); } dist = sqrt ( dist ); free ( pn ); return dist; } /******************************************************************************/ double *ellipse_point_near_2d ( double r1, double r2, double p[2] ) /******************************************************************************/ /* Purpose: ellipse_point_near_2d() finds the nearest point on an ellipse in 2D. Discussion: An ellipse in standard form has a center at the origin, and axes aligned with the coordinate axes. Any point P on the ellipse satisfies ( X / A )^2 + ( Y / B )^2 == 1 The nearest point PN on the ellipse has the property that the line from PN to P is normal to the ellipse. Points on the ellipse can be parameterized by T, to have the form ( R1 * cos ( T ), R2 * sin ( T ) ). The tangent vector to the ellipse has the form ( -R1 * sin ( T ), R2 * cos ( T ) ) At PN, the dot product of this vector with ( P - PN ) must be zero: - R1 * sin ( T ) * ( X - R1 * cos ( T ) ) + R2 * cos ( T ) * ( Y - R2 * sin ( T ) ) = 0 This nonlinear equation for T can be solved by Newton's method. Licensing: This code is distributed under the GNU LGPL license. Modified: 17 May 2010 Author: John Burkardt Input: double R1, R2, the ellipse parameters. Normally, these are both positive quantities. Generally, they are also distinct. Input, double P[2], the point. Output: double ELLIPSE_POINT_NEAR_2D[2], the point on the ellipse which is closest to P. */ { double ct; double f; double fp; int iteration; int iteration_max = 100; const double r8_pi = 3.141592653589793; double *pn; double st; double t; double x; double y; x = fabs ( p[0] ); y = fabs ( p[1] ); if ( y == 0.0 && r1 * r1 - r2 * r2 <= r1 * x ) { t = 0.0; } else if ( x == 0.0 && r2 * r2 - r1 * r1 <= r2 * y ) { t = r8_pi / 2.0; } else { if ( y == 0.0 ) { y = sqrt ( DBL_EPSILON ) * fabs ( r2 ); } if ( x == 0.0 ) { x = sqrt ( DBL_EPSILON ) * fabs ( r1 ); } /* Initial parameter T: */ t = atan2 ( y, x ); iteration = 0; for ( ; ; ) { ct = cos ( t ); st = sin ( t ); f = ( x - fabs ( r1 ) * ct ) * fabs ( r1 ) * st - ( y - fabs ( r2 ) * st ) * fabs ( r2 ) * ct; if ( fabs ( f ) <= 100.0 * DBL_EPSILON ) { break; } if ( iteration_max <= iteration ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "ELLIPSE_POINT_NEAR_2D - Warning!\n" ); fprintf ( stderr, " Reached iteration limit.\n" ); fprintf ( stderr, " T = %f\n", t ); fprintf ( stderr, " F = %f\n", f ); break; } iteration = iteration + 1; fp = r1 * r1 * st * st + r2 * r2 * ct * ct + ( x - fabs ( r1 ) * ct ) * fabs ( r1