# include # include # include # include # include # include "sphere_triangle_quad.h" /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_power ( int i, int j ) /******************************************************************************/ /* Purpose: I4_POWER returns the value of I^J. Licensing: This code is distributed under the MIT license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input, int I, J, the base and the power. J should be nonnegative. Output, int I4_POWER, the value of I^J. */ { int k; int value; if ( j < 0 ) { if ( i == 1 ) { value = 1; } else if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" ); exit ( 1 ); } else { value = 0; } } else if ( j == 0 ) { if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" ); exit ( 1 ); } else { value = 1; } } else if ( j == 1 ) { value = i; } else { value = 1; for ( k = 1; k <= j; k++ ) { value = value * i; } } return value; } /******************************************************************************/ 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 r8_pi = 3.141592653589793; if ( c <= -1.0 ) { angle = r8_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 r8_pi = 3.141592653589793; double theta; double theta_0; /* Special cases: */ if ( x == 0.0 ) { if ( 0.0 < y ) { theta = r8_pi / 2.0; } else if ( y < 0.0 ) { theta = 3.0 * r8_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 = r8_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 = r8_pi - theta_0; } else if ( x < 0.0 && y < 0.0 ) { theta = r8_pi + theta_0; } else if ( 0.0 < x && y < 0.0 ) { theta = 2.0 * r8_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. */ { const 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; } /******************************************************************************/ double r8vec_dot_product ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's. Licensing: This code is distributed under the MIT license. Modified: 26 July 2007 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], A2[N], the two vectors to be considered. Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a1[i] * a2[i]; } return value; } /******************************************************************************/ double r8vec_norm ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_NORM returns the L2 norm of an R8VEC. Discussion: The vector L2 norm is defined as: R8VEC_NORM = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ). Licensing: This code is distributed under the MIT license. Modified: 01 March 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in A. Input, double A[N], the vector whose L2 norm is desired. Output, double R8VEC_NORM, the L2 norm of A. */ { int i; double v; v = 0.0; for ( i = 0; i < n; i++ ) { v = v + a[i] * a[i]; } v = sqrt ( v ); return v; } /******************************************************************************/ void r8vec_polarize ( int n, double a[], double p[], double a_normal[], double a_parallel[] ) /******************************************************************************/ /* Purpose: R8VEC_POLARIZE decomposes an R8VEC into normal and parallel components. Discussion: An R8VEC is a vector of R8's. The (nonzero) vector P defines a direction. The vector A can be written as the sum A = A_normal + A_parallel where A_parallel is a linear multiple of P, and A_normal is perpendicular to P. Licensing: This code is distributed under the MIT license. Modified: 22 August 2010 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double A[N], the vector to be polarized. Input, double P[N], the polarizing direction. Output, double A_NORMAL[N], A_PARALLEL[N], the normal and parallel components of A. */ { double a_dot_p; int i; double p_norm; p_norm = 0.0; for ( i = 0; i < n; i++ ) { p_norm = p_norm + pow ( p[i], 2 ); } p_norm = sqrt ( p_norm ); if ( p_norm == 0.0 ) { for ( i = 0; i < n; i++ ) { a_normal[i] = a[i]; } for ( i = 0; i < n; i++ ) { a_parallel[i] = 0.0; } return; } a_dot_p = 0.0; for ( i = 0; i < n; i++ ) { a_dot_p = a_dot_p + a[i] * p[i]; } a_dot_p = a_dot_p / p_norm; for ( i = 0; i < n; i++ ) { a_parallel[i] = a_dot_p * p[i] / p_norm; } for ( i = 0; i < n; i++ ) { a_normal[i] = a[i] - a_parallel[i]; } return; } /******************************************************************************/ void r8vec_transpose_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_TRANSPOSE_PRINT prints an R8VEC "transposed". Discussion: An R8VEC is a vector of R8's. Example: A = (/ 1.0, 2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8, 10.9, 11.0 /) TITLE = 'My vector: ' My vector: 1.0 2.1 3.2 4.3 5.4 6.5 7.6 8.7 9.8 10.9 11.0 Licensing: This code is distributed under the MIT license. Modified: 12 November 2010 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; int ihi; int ilo; printf ( "\n" ); printf ( "%s\n", title ); if ( n <= 0 ) { printf ( " (Empty)\n" ); return; } for ( ilo = 0; ilo < n; ilo = ilo + 5 ) { ihi = i4_min ( ilo + 5, n ); for ( i = ilo; i < ihi; i++ ) { printf ( " %12g", a[i] ); } printf ( "\n" ); } return; } /******************************************************************************/ double sphere01_distance_xyz ( double xyz1[3], double xyz2[3] ) /******************************************************************************/ /* Purpose: SPHERE01_DISTANCE_XYZ computes great circle distances on a unit sphere. Discussion: XYZ coordinates are used. We assume the points XYZ1 and XYZ2 lie on the unit sphere. This computation is a special form of the Vincenty formula. It should be less sensitive to errors associated with very small or very large angular separations. Licensing: This code is distributed under the MIT license. Modified: 23 September 2010 Author: John Burkardt Reference: "Great-circle distance", Wikipedia. Parameters: Input, double XYZ1[3], the coordinates of the first point. Input, double XYZ2[3], the coordinates of the second point. Output, double DIST, the great circle distance between the points. */ { double bot; double dist; double lat1; double lat2; double lon1; double lon2; double top; lat1 = r8_asin ( xyz1[2] ); lon1 = r8_atan ( xyz1[1], xyz1[0] ); lat2 = r8_asin ( xyz2[2] ); lon2 = r8_atan ( xyz2[1], xyz2[0] ); top = pow ( cos ( lat2 ) * sin ( lon1 - lon2 ), 2 ) + pow ( cos ( lat1 ) * sin ( lat2 ) - sin ( lat1 ) * cos ( lat2 ) * cos ( lon1 - lon2 ), 2 ); top = sqrt ( top ); bot = sin ( lat1 ) * sin ( lat2 ) + cos ( lat1 ) * cos ( lat2 ) * cos ( lon1 - lon2 ); dist = atan2 ( top, bot ); return dist; } /******************************************************************************/ double *sphere01_sample ( int n, int *seed ) /******************************************************************************/ /* Purpose: SPHERE01_SAMPLE picks random points on a unit sphere. Licensing: This code is distributed under the MIT license. Modified: 23 September 2010 Author: John Burkardt Parameters: Input, int N, the number of samples. Input/output, int *SEED, a seed for the random number generator. Output, double SPHERE01_SAMPLE[3*N], the sample points. */ { int j; double phi; double r8_pi = 3.141592653589793; double theta; double vdot; double *x; x = ( double * ) malloc ( 3 * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { /* Pick a uniformly random VDOT, which must be between -1 and 1. This represents the dot product of the random vector with the Z unit vector. Note: this works because the surface area of the sphere between Z and Z + dZ is independent of Z. So choosing Z uniformly chooses a patch of area uniformly. */ vdot = r8_uniform_01 ( seed ); vdot = 2.0 * vdot - 1.0; phi = r8_acos ( vdot ); /* Pick a uniformly random rotation between 0 and 2 Pi around the axis of the Z vector. */ theta = r8_uniform_01 ( seed ); theta = 2.0 * r8_pi * theta; x[0+j*3] = cos ( theta ) * sin ( phi ); x[1+j*3] = sin ( theta ) * sin ( phi ); x[2+j*3] = cos ( phi ); } return x; } /******************************************************************************/ double sphere01_triangle_angles_to_area ( double a, double b, double c ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_ANGLES_TO_AREA: area of spherical triangle on unit sphere. Discussion: A unit sphere centered at 0 in 3D satisfies the equation: X^2 + Y^2 + Z^2 = 1 A spherical triangle is specified by three points on the surface of the sphere. The area formula is known as Girard's formula. The area of a spherical triangle on the unit sphere is: AREA = A + B + C - PI where A, B and C are the (surface) angles of the triangle. Licensing: This code is distributed under the MIT license. Modified: 25 September 2010 Author: John Burkardt Parameters: Input, double A, B, C, the angles of the triangle. Output, double SPHERE01_TRIANGLE_ANGLES_TO_AREA, the area of the spherical triangle. */ { double area; double r8_pi = 3.141592653589793; area = a + b + c - r8_pi; return area; } /******************************************************************************/ double *sphere01_triangle_project ( double a_xyz[3], double b_xyz[3], double c_xyz[3], int f1, int f2, int f3 ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_PROJECT projects from plane to spherical triangle. Discussion: We assume that points A, B and C lie on the unit sphere, and they thus define a spherical triangle. They also, of course, define a planar triangle. Let (F1,F2,F3) be the barycentric coordinates of a point in this planar triangle. This function determines the coordinates of the point in the planar triangle identified by the barycentric coordinates, and returns the coordinates of the projection of that point onto the unit sphere. Licensing: This code is distributed under the MIT license. Modified: 23 September 2010 Author: John Burkardt Parameters: Input, double A_XYZ[3], B_XYZ[3], C_XYZ[3], the coordinates of the points A, B, and C. Input, int F1, F2, F3, the barycentric coordinates of a point in the triangle ABC. Normally, these coordinates would be real numbers, and would sum to 1. For convenience, we allow these to be integers which must be divided by F1+F2+F3. Output, double NODE_XYZ[3], the coordinates of the point on the unit sphere which is the projection of the point on the plane whose barycentric coordinates with respect to A, B, and C is (F1,F2,F3)/(F1+F2+F3). */ { int i; double *node_xyz; double norm; node_xyz = ( double * ) malloc ( 3 * sizeof ( double ) ); for ( i = 0; i < 3; i++ ) { node_xyz[i] = ( ( double ) ( f1 ) * a_xyz[i] + ( double ) ( f2 ) * b_xyz[i] + ( double ) ( f3 ) * c_xyz[i] ) / ( double ) ( f1 + f2 + f3 ); } norm = r8vec_norm ( 3, node_xyz ); for ( i = 0; i < 3; i++ ) { node_xyz[i] = node_xyz[i] / norm; } return node_xyz; } /******************************************************************************/ double *sphere01_triangle_project2 ( double a_xyz[3], double b_xyz[3], double c_xyz[3], int f1, int f2, int f3 ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_PROJECT2 projects from plane to spherical triangle. Discussion: We assume that points A, B and C lie on the unit sphere, and they thus define a spherical triangle. They also, of course, define a planar triangle. Let (F1,F2,F3) be the barycentric coordinates of a point in this planar triangle. This function determines the coordinates of the point in the planar triangle identified by the barycentric coordinates, and returns the coordinates of the projection of that point onto the unit sphere. Licensing: This code is distributed under the MIT license. Modified: 23 September 2010 Author: John Burkardt Parameters: Input, double A_XYZ(3), B_XYZ(3), C_XYZ(3), the coordinates of the points A, B, and C. Input, int F1, F2, F3, the barycentric coordinates of a point in the triangle ABC. Normally, these coordinates would be real numbers, and would sum to 1. For convenience, we allow these to be integers which must be divided by F1+F2+F3. Output, double SPHERE01_TRIANGLE_PROJECT2[3], the coordinates of the point on the unit sphere which is the projection of the point on the plane whose barycentric coordinates with respect to A, B, and C is (F1,F2,F3)/(F1+F2+F3). */ { double ab[3]; double ac[3]; double acn[3]; double acp[3]; double angle; double bn[3]; double bp[3]; double cn[3]; double cp[3]; int i; double *node_xyz; double norm; double theta_ab; double theta_ac; double theta_bc; node_xyz = ( double * ) malloc ( 3 * sizeof ( double ) ); /* This check avoids 0/0 calculations later. */ if ( f2 == 0 && f3 == 0 ) { for ( i = 0; i < 3; i++ ) { node_xyz[i] = a_xyz[i]; } return node_xyz; } else if ( f1 == 0 && f3 == 0 ) { for ( i = 0; i < 3; i++ ) { node_xyz[i] = b_xyz[i]; } return node_xyz; } else if ( f1 == 0 && f2 == 0 ) { for ( i = 0; i < 3; i++ ) { node_xyz[i] = c_xyz[i]; } return node_xyz; } /* Determine the angular distances (A,B) and (A,C). */ theta_ab = sphere01_distance_xyz ( a_xyz, b_xyz ); theta_ac = sphere01_distance_xyz ( a_xyz, c_xyz ); /* Polarize B = BP + BN Normalize BN, Same for C. */ r8vec_polarize ( 3, b_xyz, a_xyz, bn, bp ); norm = r8vec_norm ( 3, bn ); for ( i = 0; i < 3; i++ ) { bn[i] = bn[i] / norm; } r8vec_polarize ( 3, c_xyz, a_xyz, cn, cp ); norm = r8vec_norm ( 3, cn ); for ( i = 0; i < 3; i++ ) { cn[i] = cn[i] / norm; } /* Determine AB and AC that use cos ( ( F2 + F3 ) / ( F1 + F2 + F3 ) ) of A and cos ( F1 / ( F1 + F2 + F3 ) ) of B or C. */ angle = ( ( double ) ( f2 + f3 ) * theta_ab ) / ( double ) ( f1 + f2 + f3 ); for ( i = 0; i < 3; i++ ) { ab[i] = cos ( angle ) * a_xyz[i] + sin ( angle ) * bn[i]; } angle = ( ( double ) ( f2 + f3 ) * theta_ac ) / ( double ) ( f1 + f2 + f3 ); for ( i = 0; i < 3; i++ ) { ac[i] = cos ( angle ) * a_xyz[i] + sin ( angle ) * cn[i]; } /* Determine the angular distance between AB and AC. */ theta_bc = sphere01_distance_xyz ( ab, ac ); /* Polarize AC = ACP + ACN, normalize ACN. */ r8vec_polarize ( 3, ac, ab, acn, acp ); norm = r8vec_norm ( 3, acn ); for ( i = 0; i < 3; i++ ) { acn[i] = acn[i] / norm; } /* The interval between AB and AC is marked by F2+F3+1 vertices 0 through F2+F3. */ angle = ( ( double ) ( f3 ) * theta_bc ) / ( double ) ( f2 + f3 ); for ( i = 0; i < 3; i++ ) { node_xyz[i] = cos ( angle ) * ab[i] + sin ( angle ) * acn[i]; } return node_xyz; } /******************************************************************************/ double sphere01_triangle_quad_00 ( int n, double v1[3], double v2[3], double v3[3], double f ( double x[] ), int *seed ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_00: quadrature over a triangle on the unit sphere. Discussion: This is a Monte Carlo approach. The integral is approximated by averaging the values at N random points, multiplied by the area. Licensing: This code is distributed under the MIT license. Modified: 27 September 2010 Author: John Burkardt Parameters: Input, integer N, the number of sample points. Input, real V1[3], V2[3], V3[3], the XYZ coordinates of the vertices of the triangle. Input, double F ( double x[] ), evaluates the integrand. Input/output, integer *SEED, a seed for the random number generator. Output, double SPHERE01_TRIANGLE_QUAD_00, the approximate integral. */ { double area; int j; double quad; double result; double *vc; area = sphere01_triangle_vertices_to_area ( v1, v2, v3 ); vc = sphere01_triangle_sample ( n, v1, v2, v3, seed ); quad = 0.0; for ( j = 0; j < n; j++ ) { quad = quad + f ( vc+3*j ); } result = quad * area / ( double ) ( n ); free ( vc ); return result; } /******************************************************************************/ double sphere01_triangle_quad_01 ( double v1[3], double v2[3], double v3[3], double f ( double x[] ) ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_01: quadrature over a triangle on the unit sphere. Discussion: The integral is approximated by the value at the centroid, multiplied by the area. Licensing: This code is distributed under the MIT license. Modified: 27 September 2010 Author: John Burkardt Parameters: Input, real V1[3], V2[3], V3[3], the XYZ coordinates of the vertices of the triangle. Input, double F ( double x[] ), evaluates the integrand. Output, double SPHERE01_TRIANGLE_QUAD_01, the approximate integral. */ { double area; double quad; double result; double *vc; area = sphere01_triangle_vertices_to_area ( v1, v2, v3 ); vc = sphere01_triangle_vertices_to_centroid ( v1, v2, v3 ); quad = f ( vc ); result = quad * area; free ( vc ); return result; } /******************************************************************************/ double sphere01_triangle_quad_02 ( double v1[3], double v2[3], double v3[3], double f ( double x[] ) ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_02: quadrature over a triangle on the unit sphere. Discussion: The integral is approximated by the average of the vertex values, multiplied by the area. Licensing: This code is distributed under the MIT license. Modified: 27 September 2010 Author: John Burkardt Parameters: Input, real V1[3], V2[3], V3[3], the XYZ coordinates of the vertices of the triangle. Input, double F ( double x[] ), evaluates the integrand. Output, double SPHERE01_TRIANGLE_QUAD_02, the approximate integral. */ { double area; double quad; double result; area = sphere01_triangle_vertices_to_area ( v1, v2, v3 ); quad = ( f ( v1 ) + f ( v2 ) + f ( v3 ) ) / 3.0; result = quad * area; return result; } /******************************************************************************/ double sphere01_triangle_quad_03 ( double v1[3], double v2[3], double v3[3], double f ( double x[] ) ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_03: quadrature over a triangle on the unit sphere. Discussion: The integral is approximated by the average of the midside values, multiplied by the area. Licensing: This code is distributed under the MIT license. Modified: 27 September 2010 Author: John Burkardt Parameters: Input, real V1[3], V2[3], V3[3], the XYZ coordinates of the vertices of the triangle. Input, double F ( double x[] ), evaluates the integrand. Output, double SPHERE01_TRIANGLE_QUAD_03, the approximate integral. */ { double area; double quad; double result; double v4[3]; double v5[3]; double v6[3]; area = sphere01_triangle_vertices_to_area ( v1, v2, v3 ); sphere01_triangle_vertices_to_midpoints ( v1, v2, v3, v4, v5, v6 ); quad = ( f ( v4 ) + f ( v5 ) + f ( v6 ) ) / 3.0; result = quad * area; return result; } /******************************************************************************/ double sphere01_triangle_quad_icos1c ( double a_xyz[3], double b_xyz[3], double c_xyz[], int factor, double fun ( double x[] ), int *node_num ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_ICOS1C: centroid rule, subdivide then project. Discussion: This function estimates an integral over a spherical triangle on the unit sphere. This function subdivides each edge of the triangle into FACTOR subedges. These edges define a grid within the triangle. The centroids of these triangles can be determined. All of these calculations are done, essentially, on the FLAT faces of the planar triangle. Only then are the triangle vertices and centroids projected to the sphere. Licensing: This code is distributed under the MIT license. Modified: 28 September 2010 Author: John Burkardt Parameters: Input, double A_XYZ[3], B_XYZ[3], C_XYZ[3], the vertices of the spherical triangle. Input, int FACTOR, the subdivision factor, which must be at least 1. Input, double FUN ( double x[] ), evaluates the integrand. Output, int *NODE_NUM, the number of evaluation points. Output, double SPHERE01_TRIANGLE_QUAD_ICOS1C, the estimated integral. */ { double *a2_xyz; double area; double area_total; double *b2_xyz; double *c2_xyz; int f1; int f2; int f3; double *node_xyz; double result; double v; /* Initialize the integral data. */ result = 0.0; area_total = 0.0; *node_num = 0; /* Some subtriangles will have the same direction as the face. Generate each in turn, by determining the barycentric coordinates of the centroid (F1,F2,F3), from which we can also work out the barycentric coordinates of the vertices of the subtriangle. */ for ( f3 = 1; f3 <= 3 * factor - 2; f3 = f3 + 3 ) { for ( f2 = 1; f2 <= 3 * factor - f3 - 1; f2 = f2 + 3 ) { f1 = 3 * factor - f3 - f2; node_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 + 2, f2 - 1, f3 - 1 ); b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 - 1, f2 + 2, f3 - 1 ); c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 - 1, f2 - 1, f3 + 2 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); v = fun ( node_xyz ); *node_num = *node_num + 1; result = result + area * v; area_total = area_total + area; free ( a2_xyz ); free ( b2_xyz ); free ( c2_xyz ); free ( node_xyz ); } } /* The other subtriangles have the opposite direction from the face. Generate each in turn, by determining the barycentric coordinates of the centroid (F1,F2,F3), from which we can also work out the barycentric coordinates of the vertices of the subtriangle. */ for ( f3 = 2; f3 <= 3 * factor - 4; f3 = f3 + 3 ) { for ( f2 = 2; f2 <= 3 * factor - f3 - 2; f2 = f2 + 3 ) { f1 = 3 * factor - f3 - f2; node_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 - 2, f2 + 1, f3 + 1 ); b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 + 1, f2 - 2, f3 + 1 ); c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 + 1, f2 + 1, f3 - 2 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); v = fun ( node_xyz ); *node_num = *node_num + 1; result = result + area * v; area_total = area_total + area; free ( a2_xyz ); free ( b2_xyz ); free ( c2_xyz ); free ( node_xyz ); } } return result; } /******************************************************************************/ double sphere01_triangle_quad_icos1m ( double a_xyz[3], double b_xyz[3], double c_xyz[], int factor, double fun ( double x[] ), int *node_num ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_ICOS1M: midside rule, subdivide then project. Discussion: This function estimates an integral over a spherical triangle on the unit sphere. This function subdivides each edge of the triangle into FACTOR subedges. These edges define a grid within the triangle. The midsides of the edges of these triangles can be determined. All of these calculations are done, essentially, on the FLAT faces of the planar triangle. Only then are the triangle vertices and midsides projected to the sphere. Licensing: This code is distributed under the MIT license. Modified: 29 September 2010 Author: John Burkardt Parameters: Input, double A_XYZ[3], B_XYZ[3], C_XYZ[3], the vertices of the spherical triangle. Input, int FACTOR, the subdivision factor, which must be at least 1. Input, double FUN ( double x[] ), evaluates the integrand. Output, int *NODE_NUM, the number of evaluation points. Output, double SPHERE01_TRIANGLE_QUAD_ICOS1M, the estimated integral. */ { double *a2_xyz; double *a3_xyz; double area; double area_total; double *b2_xyz; double *b3_xyz; double *c2_xyz; double *c3_xyz; int f1; int f2; int f3; double result; double va; double vb; double vc; /* Initialize the integral data. */ result = 0.0; area_total = 0.0; *node_num = 0; /* Some subtriangles will have the same direction as the face. */ for ( f1 = 0; f1 <= factor - 1; f1++ ) { for ( f2 = 0; f2 <= factor - f1 - 1; f2++ ) { f3 = factor - f1 - f2; a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 + 1, f2, f3 - 1 ); b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2 + 1, f3 - 1 ); c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); a3_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, 2 * f1 + 1, 2 * f2 + 1, 2 * f3 - 2 ); b3_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, 2 * f1, 2 * f2 + 1, 2 * f3 - 1 ); c3_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, 2 * f1 + 1, 2 * f2, 2 * f3 - 1 ); *node_num = *node_num + 3; va = fun ( a3_xyz ); vb = fun ( b3_xyz ); vc = fun ( c3_xyz ); result = result + area * ( va + vb + vc ) / 3.0; area_total = area_total + area; free ( a2_xyz ); free ( a3_xyz ); free ( b2_xyz ); free ( b3_xyz ); free ( c2_xyz ); free ( c3_xyz ); } } /* The other subtriangles have the opposite direction from the face. */ for ( f3 = 0; f3 <= factor - 2; f3++ ) { for ( f2 = 1; f2 <= factor - f3 - 1; f2++ ) { f1 = factor - f2 - f3; a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 - 1, f2, f3 + 1 ); b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2 - 1, f3 + 1 ); c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); a3_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, 2 * f1 - 1, 2 * f2 - 1, 2 * f3 + 2 ); b3_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, 2 * f1, 2 * f2 - 1, 2 * f3 + 1 ); c3_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, 2 * f1 - 1, 2 * f2, 2 * f3 + 1 ); *node_num = *node_num + 3; va = fun ( a3_xyz ); vb = fun ( b3_xyz ); vc = fun ( c3_xyz ); result = result + area * ( va + vb + vc ) / 3.0; area_total = area_total + area; free ( a2_xyz ); free ( a3_xyz ); free ( b2_xyz ); free ( b3_xyz ); free ( c2_xyz ); free ( c3_xyz ); } } return result; } /******************************************************************************/ double sphere01_triangle_quad_icos1v ( double a_xyz[3], double b_xyz[3], double c_xyz[], int factor, double fun ( double x[] ), int *node_num ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_ICOS1V: vertex rule, subdivide then project. Discussion: This function estimates an integral over a spherical triangle on the unit sphere. This function subdivides each edge of the triangle into FACTOR subedges. These edges define a grid within the triangle. The vertices of these triangles can be determined. All of these calculations are done, essentially, on the FLAT faces of the planar triangle. Only then are the triangle vertices projected to the sphere. Licensing: This code is distributed under the MIT license. Modified: 29 September 2010 Author: John Burkardt Parameters: Input, double A_XYZ[3], B_XYZ[3], C_XYZ[3], the vertices of the spherical triangle. Input, int FACTOR, the subdivision factor, which must be at least 1. Input, double FUN ( double x[] ), evaluates the integrand. Output, int *NODE_NUM, the number of evaluation points. Output, double SPHERE01_TRIANGLE_QUAD_ICOS1V, the estimated integral. */ { double *a2_xyz; double area; double area_total; double *b2_xyz; double *c2_xyz; int f1; int f2; int f3; double result; double va; double vb; double vc; /* Initialize the integral data. */ result = 0.0; area_total = 0.0; *node_num = 0; /* Some subtriangles will have the same direction as the face. */ for ( f1 = 0; f1 <= factor - 1; f1++ ) { for ( f2 = 0; f2 <= factor - f1 - 1; f2++ ) { f3 = factor - f1 - f2; a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 + 1, f2, f3 - 1 ); b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2 + 1, f3 - 1 ); c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); *node_num = *node_num + 1; va = fun ( a2_xyz ); vb = fun ( b2_xyz ); vc = fun ( c2_xyz ); result = result + area * ( va + vb + vc ) / 3.0; area_total = area_total + area; free ( a2_xyz ); free ( b2_xyz ); free ( c2_xyz ); } } /* The other subtriangles have the opposite direction from the face. */ for ( f3 = 0; f3 <= factor - 2; f3++ ) { for ( f2 = 1; f2 <= factor - f3 - 1; f2++ ) { f1 = factor - f2 - f3; a2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1 - 1, f2, f3 + 1 ); b2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2 - 1, f3 + 1 ); c2_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); *node_num = *node_num + 1; va = fun ( a2_xyz ); vb = fun ( b2_xyz ); vc = fun ( c2_xyz ); result = result + area * ( va + vb + vc ) / 3.0; area_total = area_total + area; free ( a2_xyz ); free ( b2_xyz ); free ( c2_xyz ); } } return result; } /******************************************************************************/ double sphere01_triangle_quad_icos2v ( double a_xyz[3], double b_xyz[3], double c_xyz[], int factor, double fun ( double x[] ), int *node_num ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_QUAD_ICOS2V: vertex rule, subdivide then project. Discussion: This function estimates an integral over a spherical triangle on the unit sphere. This function subdivides each edge of the triangle into FACTOR subedges. These edges define a grid within the triangle. The vertices of these triangles can be determined. All of these calculations are done, essentially, on the FLAT faces of the planar triangle. Only then are the triangle vertices projected to the sphere. This function uses a more sophisticated projection scheme than that used by SPHERE01_TRIANGLE_QUAD_ICOS1V. Licensing: This code is distributed under the MIT license. Modified: 29 September 2010 Author: John Burkardt Parameters: Input, double A_XYZ[3], B_XYZ[3], C_XYZ[3], the vertices of the spherical triangle. Input, int FACTOR, the subdivision factor, which must be at least 1. Input, double FUN ( double x[] ), evaluates the integrand. Output, int *NODE_NUM, the number of evaluation points. Output, double SPHERE01_TRIANGLE_QUAD_ICOS2V, the estimated integral. */ { double *a2_xyz; double area; double area_total; double *b2_xyz; double *c2_xyz; int f1; int f2; int f3; double result; double va; double vb; double vc; /* Initialize the integral data. */ result = 0.0; area_total = 0.0; *node_num = 0; /* Some subtriangles will have the same direction as the face. */ for ( f1 = 0; f1 <= factor - 1; f1++ ) { for ( f2 = 0; f2 <= factor - f1 - 1; f2++ ) { f3 = factor - f1 - f2; a2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1 + 1, f2, f3 - 1 ); b2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1, f2 + 1, f3 - 1 ); c2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); *node_num = *node_num + 1; va = fun ( a2_xyz ); vb = fun ( b2_xyz ); vc = fun ( c2_xyz ); result = result + area * ( va + vb + vc ) / 3.0; area_total = area_total + area; free ( a2_xyz ); free ( b2_xyz ); free ( c2_xyz ); } } /* The other subtriangles have the opposite direction from the face. */ for ( f3 = 0; f3 <= factor - 2; f3++ ) { for ( f2 = 1; f2 <= factor - f3 - 1; f2++ ) { f1 = factor - f3 - f2; a2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1 - 1, f2, f3 + 1 ); b2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1, f2 - 1, f3 + 1 ); c2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1, f2, f3 ); area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz ); *node_num = *node_num + 1; va = fun ( a2_xyz ); vb = fun ( b2_xyz ); vc = fun ( c2_xyz ); result = result + area * ( va + vb + vc ) / 3.0; area_total = area_total + area; free ( a2_xyz ); free ( b2_xyz ); free ( c2_xyz ); } } return result; } /******************************************************************************/ double *sphere01_triangle_sample ( int n, double v1[3], double v2[3], double v3[3], int *seed ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_SAMPLE: sample points from triangle on unit sphere. Discussion: The sphere has center 0 and radius 1. A spherical triangle on the surface of the unit sphere contains those points with radius R = 1, bounded by the vertices V1, V2, V3. Licensing: This code is distributed under the MIT license. Modified: 24 September 2010 Author: John Burkardt Reference: James Arvo, Stratified sampling of spherical triangles, Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH '95, pages 437-438, 1995. Parameters: Input, int N, the number of points. Input, double V1[3], V2[3], V3[3], the XYZ coordinates of the vertices of the spherical triangle. Input/output, int *SEED, a seed for the random number generator. Output, double SPHERE01_TRIANGLE_SAMPLE[3*N], the XYZ coordinates of the sample points. */ { double a; double alpha; double area; double area_hat; double b; double beta; double c; double gamma; int i; int j; double norm; double q; double s; double t; double u; double v; double v31[3]; double v4[3]; double v42[3]; double w; double *x; double xsi1; double xsi2; double z; sphere01_triangle_vertices_to_sides ( v1, v2, v3, &a, &b, &c ); sphere01_triangle_sides_to_angles ( a, b, c, &alpha, &beta, &gamma ); area = sphere01_triangle_angles_to_area ( alpha, beta, gamma ); x = ( double * ) malloc ( 3 * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { /* Select the new area. */ xsi1 = r8_uniform_01 ( seed ); area_hat = xsi1 * area; /* Compute the sine and cosine of the angle phi. */ s = sin ( area_hat - alpha ); t = cos ( area_hat - alpha ); /* Compute the pair that determines beta_hat. */ u = t - cos ( alpha ); v = s + sin ( alpha ) * cos ( c ); /* Q is the cosine of the new edge length b_hat. */ q = ( ( v * t - u * s ) * cos ( alpha ) - v ) / ( ( v * s + u * t ) * sin ( alpha ) ); /* We very occasionally get a Q value out of bounds. */ q = fmax ( q, - 1.0 ); q = fmin ( q, + 1.0 ); /* V31 = normalized ( V3 - ( V3 dot V1 ) * V1 ) */ w = r8vec_dot_product ( 3, v3, v1 ); for ( i = 0; i < 3; i++ ) { v31[i] = v3[i] - w * v1[i]; } norm = r8vec_norm ( 3, v31 ); for ( i = 0; i < 3; i++ ) { v31[i] = v31[i] / norm; } /* V4 is the third vertex of the subtriangle V1, V2, V4. */ for ( i = 0; i < 3; i++ ) { v4[i] = q * v1[i] + sqrt ( 1.0 - q * q ) * v31[i]; } /* Select cos theta, which will sample along the edge from V2 to V4. */ xsi2 = r8_uniform_01 ( seed ); z = 1.0 - xsi2 * ( 1.0 - r8vec_dot_product ( 3, v4, v2 ) ); /* V42 = normalized ( V4 - ( V4 dot V2 ) * V2 ) */ w = r8vec_dot_product ( 3, v4, v2 ); for ( i = 0; i < 3; i++ ) { v42[i] = v4[i] - w * v2[i]; } norm = r8vec_norm ( 3, v42 ); for ( i = 0; i < 3; i++ ) { v42[i] = v42[i] / norm; } /* Construct the point. */ for ( i = 0; i < 3; i++ ) { x[i+j*3] = z * v2[i] + sqrt ( 1.0 - z * z ) * v42[i]; } } return x; } /******************************************************************************/ void sphere01_triangle_sides_to_angles ( double as, double bs, double cs, double *a, double *b, double *c ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_SIDES_TO_ANGLES: spherical triangle angles on unit sphere. Licensing: This code is distributed under the MIT license. Modified: 25 September 2010 Author: John Burkardt Parameters: Input, double AS, BS, CS, the (geodesic) length of the sides of the triangle. Output, double *A, *B, *C, the spherical angles of the triangle. Angle A is opposite the side of length AS, and so on. */ { double asu; double bsu; double csu; double ssu; double tan_a2; double tan_b2; double tan_c2; asu = as; bsu = bs; csu = cs; ssu = ( asu + bsu + csu ) / 2.0; tan_a2 = sqrt ( ( sin ( ssu - bsu ) * sin ( ssu - csu ) ) / ( sin ( ssu ) * sin ( ssu - asu ) ) ); *a = 2.0 * atan ( tan_a2 ); tan_b2 = sqrt ( ( sin ( ssu - asu ) * sin ( ssu - csu ) ) / ( sin ( ssu ) * sin ( ssu - bsu ) ) ); *b = 2.0 * atan ( tan_b2 ); tan_c2 = sqrt ( ( sin ( ssu - asu ) * sin ( ssu - bsu ) ) / ( sin ( ssu ) * sin ( ssu - csu ) ) ); *c = 2.0 * atan ( tan_c2 ); return; } /******************************************************************************/ double sphere01_triangle_vertices_to_area ( double v1[3], double v2[3], double v3[3] ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_VERTICES_TO_AREA: area of spherical triangle on unit sphere. Discussion: A sphere centered at 0 in 3D satisfies the equation: X*X + Y*Y + Z*Z = 1 A spherical triangle is specified by three points on the surface of the sphere. The area formula is known as Girard's formula. The area of a spherical triangle on the unit sphere is: AREA = A + B + C - PI where A, B and C are the (surface) angles of the triangle. Licensing: This code is distributed under the MIT license. Modified: 25 September 2010 Author: John Burkardt Parameters: Input, double V1[3], V2[3], V3[3], the vertices of the triangle. Output, double SPHERE_TRIANGLE_VERTICES_TO_AREA_3D, the area of the spherical triangle. */ { double area; double a; double as; double b; double bs; double c; double cs; /* Compute the lengths of the sides of the spherical triangle. */ sphere01_triangle_vertices_to_sides ( v1, v2, v3, &as, &bs, &cs ); /* Get the spherical angles. */ sphere01_triangle_sides_to_angles ( as, bs, cs, &a, &b, &c ); /* Get the area */ area = sphere01_triangle_angles_to_area ( a, b, c ); return area; } /******************************************************************************/ double *sphere01_triangle_vertices_to_centroid ( double v1[3], double v2[3], double v3[3] ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_VERTICES_TO_CENTROID: spherical triangle centroid on unit sphere. Discussion: A unit sphere centered at 0 in 3D satisfies the equation: X*X + Y*Y + Z*Z = 1 A spherical triangle is specified by three points on the sphere. The (true) centroid of a spherical triangle is the point VT = (XT,YT,ZT) = Integral ( X, Y, Z ) dArea / Integral 1 dArea Note that the true centroid does NOT, in general, lie on the sphere. The "flat" centroid VF is the centroid of the planar triangle defined by the vertices of the spherical triangle. The "spherical" centroid VS of a spherical triangle is computed by the intersection of the geodesic bisectors of the triangle angles. The spherical centroid lies on the sphere. VF, VT and VS lie on a line through the center of the sphere. We can easily calculate VF by averaging the vertices, and from this determine VS by normalizing. (Of course, we still will not have actually computed VT, which lies somewhere between VF and VS!) Licensing: This code is distributed under the MIT license. Modified: 25 September 2010 Author: John Burkardt Parameters: Input, double V1[3], V2[3], V3[3], the vertices of the triangle. Output, double SPHERE01_TRIANGLE_VERTICES_TO_CENTROID[3], the coordinates of the "spherical centroid" of the spherical triangle. */ { # define DIM_NUM 3 int i; double norm; double *vs; vs = ( double * ) malloc ( 3 * sizeof ( double ) ); for ( i = 0; i < DIM_NUM; i++ ) { vs[i] = ( v1[i] + v2[i] + v3[i] ) / 3.0; } norm = r8vec_norm ( DIM_NUM, vs ); for ( i = 0; i < DIM_NUM; i++ ) { vs[i] = vs[i] / norm; } return vs; # undef DIM_NUM } /******************************************************************************/ void sphere01_triangle_vertices_to_midpoints ( double v1[3], double v2[3], double v3[3], double m1[3], double m2[3], double m3[3] ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_VERTICES_TO_MIDPOINTS gets the midsides of a spherical triangle. Discussion: The points are assumed to lie on the unit sphere. Licensing: This code is distributed under the MIT license. Modified: 23 September 2010 Author: John Burkardt Parameters: Input, double V1[3], V2[3], V3[3], the vertices of the triangle. Output, double M1[3], M2[3], M3[3], the coordinates of the midpoints of the sides of the spherical triangle. */ { int i; double norm; for ( i = 0; i < 3; i++ ) { m1[i] = ( v1[i] + v2[i] ) / 2.0; } norm = r8vec_norm ( 3, m1 ); for ( i = 0; i < 3; i++ ) { m1[i] = m1[i] / norm; } for ( i = 0; i < 3; i++ ) { m2[i] = ( v2[i] + v3[i] ) / 2.0; } norm = r8vec_norm ( 3, m2 ); for ( i = 0; i < 3; i++ ) { m2[i] = m2[i] / norm; } for ( i = 0; i < 3; i++ ) { m3[i] = ( v3[i] + v1[i] ) / 2.0; } norm = r8vec_norm ( 3, m3 ); for ( i = 0; i < 3; i++ ) { m3[i] = m3[i] / norm; } return; } /******************************************************************************/ void sphere01_triangle_vertices_to_sides ( double v1[3], double v2[3], double v3[3], double *as, double *bs, double *cs ) /******************************************************************************/ /* Purpose: SPHERE01_TRIANGLE_VERTICES_TO_SIDES_3D: spherical triangle sides on the unit sphere. Discussion: We can use the ACOS system call here, but the ARC_COSINE routine will automatically take care of cases where the input argument is (usually slightly) out of bounds. Licensing: This code is distributed under the MIT license. Modified: 25 September 2010 Author: John Burkardt Parameters: Input, double V1[3], V2[3], V3[3], the vertices of the spherical triangle. Output, double *AS, *BS, *CS, the (geodesic) length of the sides of the triangle. */ { *as = r8_acos ( r8vec_dot_product ( 3, v2, v3 ) ); *bs = r8_acos ( r8vec_dot_product ( 3, v3, v1 ) ); *cs = r8_acos ( r8vec_dot_product ( 3, v1, v2 ) ); return; } /******************************************************************************/ 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 }