# include # include # include # include # include "sphere_grid.h" /******************************************************************************/ double arc_cosine ( double c ) /******************************************************************************/ /* Purpose: ARC_COSINE 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 GNU LGPL license. Modified: 13 June 2002 Author: John Burkardt Parameters: Input, double C, the argument, the cosine of an angle. Output, double ARC_COSINE, an angle whose cosine is C. */ { double angle; 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 arc_sine ( double s ) /******************************************************************************/ /* Purpose: ARC_SINE 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. This routine truncates arguments outside the range. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 May 2010 Author: John Burkardt Parameters: Input, double S, the argument, the sine of an angle. Output, double ARC_SINE, an angle whose sine is S. */ { double angle; double pi = 3.141592653589793; if ( s <= -1.0 ) { angle = - pi / 2.0; } else if ( 1.0 <= s ) { angle = pi / 2.0; } else { angle = asin ( s ); } return angle; } /******************************************************************************/ double atan4 ( double y, double x ) /******************************************************************************/ /* Purpose: ATAN4 computes the inverse tangent of the ratio Y / X. Discussion: ATAN4 returns an angle whose tangent is ( Y / X ), a job which the built in functions ATAN and ATAN2 already do. However: * ATAN4 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; * ATAN4 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 GNU LGPL license. Modified: 14 May 2010 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 ATAN4, 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 pi = 3.141592653589793; /* Special cases: */ if ( x == 0.0 ) { if ( 0.0 < y ) { return ( pi / 2.0 ); } else if ( y < 0.0 ) { return ( 3.0 * pi / 2.0 ); } else if ( y == 0.0 ) { return ( 0.0 ); } } else if ( y == 0.0 ) { if ( 0.0 < x ) { return 0.0; } else if ( x < 0.0 ) { return pi; } } /* We assume that ATAN2 is reliable when both arguments are positive. */ if ( 0.0 < x && 0.0 < y ) { return atan2 ( y, x ); } else if ( x < 0.0 && 0.0 < y ) { return ( pi - atan2 ( y, -x ) ); } else if ( x < 0.0 && y < 0.0 ) { return ( pi + atan2 ( -y, -x ) ); } else if ( 0.0 < x && y < 0.0 ) { return ( 2.0 * pi - atan2 ( -y, x ) ); } return 0.0; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the GNU LGPL 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 GNU LGPL 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; } /******************************************************************************/ void i4mat_transpose_print ( int m, int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the GNU LGPL license. Modified: 31 January 2005 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, int A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void i4mat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 August 2010 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, int 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 10 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); /* Print the columns of the matrix, in strips of INCX. */ for ( i2lo = ilo; i2lo <= ihi; i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); fprintf ( stdout, "\n" ); /* For each row I in the current range... Write the header. */ fprintf ( stdout, " Row: " ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, "%6d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { /* Print out (up to INCX) entries in column J, that lie in the current strip. */ fprintf ( stdout, "%5d: ", j - 1 ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, "%6d ", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void i4vec_copy ( int n, int a1[], int a2[] ) /******************************************************************************/ /* Purpose: I4VEC_COPY copies an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 25 April 2007 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, int A1[N], the vector to be copied. Input, int A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ void icos_shape ( int point_num, int edge_num, int face_num, int face_order_max, double point_coord[], int edge_point[], int face_order[], int face_point[] ) /******************************************************************************/ /* Purpose: ICOS_SHAPE describes a icosahedron. Discussion: The input data required for this routine can be retrieved from ICOS_NUM. The vertices lie on the unit sphere. The dual of an icosahedron is the dodecahedron. The data has been rearranged from a previous assignment. The STRIPACK program refuses to triangulate data if the first three nodes are "collinear" on the sphere. Licensing: This code is distributed under the GNU LGPL license. Modified: 30 August 2010 Author: John Burkardt Parameters: Input, int POINT_NUM, the number of points (12). Input, int EDGE_NUM, the number of edges (30). Input, int FACE_NUM, the number of faces (20). Input, int FACE_ORDER_MAX, the maximum number of vertices per face (3). Output, double POINT_COORD[3*POINT_NUM], the point coordinates. Output, int EDGE_POINT[2*EDGE_NUM], the points that make up each edge, listed in ascending order of their indexes. Output, int FACE_ORDER[FACE_NUM], the number of vertices per face. Output, int FACE_POINT[FACE_ORDER_MAX*FACE_NUM]; FACE_POINT(I,J) contains the index of the I-th point in the J-th face. The points are listed in the counter clockwise direction defined by the outward normal at the face. The nodes of each face are ordered so that the lowest index occurs first. The faces are then sorted by nodes. */ { # define DIM_NUM 3 # define EDGE_NUM 30 # define EDGE_ORDER 2 # define FACE_NUM 20 # define POINT_NUM 12 double phi = 0.5 * ( sqrt ( 5.0 ) + 1.0 ); static int edge_point_save[EDGE_ORDER*EDGE_NUM] = { 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 1, 2, 1, 3, 1, 6, 1, 7, 2, 4, 2, 6, 2, 8, 3, 5, 3, 7, 3, 9, 4, 5, 4, 8, 4, 10, 5, 9, 5, 10, 6, 7, 6, 8, 6, 11, 7, 9, 7, 11, 8, 10, 8, 11, 9, 10, 9, 11, 10, 11 }; static int face_order_save[FACE_NUM] = { 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 }; static int face_point_save[3*FACE_NUM] = { 0, 1, 3, 0, 2, 1, 0, 3, 5, 0, 4, 2, 0, 5, 4, 1, 2, 6, 1, 6, 7, 1, 7, 3, 2, 4, 8, 2, 8, 6, 3, 7, 9, 3, 9, 5, 4, 5, 10, 4, 10, 8, 5, 9, 10, 6, 8, 11, 6, 11, 7, 7, 11, 9, 8, 10, 11, 9, 11, 10 }; double point_coord_save[DIM_NUM*POINT_NUM]; point_coord_save[0+0*3] = + phi / sqrt ( 1.0 + phi * phi ); point_coord_save[1+0*3] = + 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[2+0*3] = 0.0; point_coord_save[0+1*3] = + phi / sqrt ( 1.0 + phi * phi ); point_coord_save[1+1*3] = - 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[2+1*3] = 0.0; point_coord_save[0+2*3] = + 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[1+2*3] = 0.0; point_coord_save[2+2*3] = + phi / sqrt ( 1.0 + phi * phi ); point_coord_save[0+3*3] = + 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[1+3*3] = 0.0; point_coord_save[2+3*3] = - phi / sqrt ( 1.0 + phi * phi ); point_coord_save[0+4*3] = 0.0; point_coord_save[1+4*3] = + phi / sqrt ( 1.0 + phi * phi ); point_coord_save[2+4*3] = + 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[0+5*3] = 0.0; point_coord_save[1+5*3] = + phi / sqrt ( 1.0 + phi * phi ); point_coord_save[2+5*3] = - 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[0+6*3] = 0.0; point_coord_save[1+6*3] = - phi / sqrt ( 1.0 + phi * phi ); point_coord_save[2+6*3] = + 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[0+7*3] = 0.0; point_coord_save[1+7*3] = - phi / sqrt ( 1.0 + phi * phi ); point_coord_save[2+7*3] = - 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[0+8*3] = - 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[1+8*3] = 0.0; point_coord_save[2+8*3] = + phi / sqrt ( 1.0 + phi * phi ); point_coord_save[0+9*3] = - 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[1+9*3] = 0.0; point_coord_save[2+9*3] = - phi / sqrt ( 1.0 + phi * phi ); point_coord_save[0+10*3] = - phi / sqrt ( 1.0 + phi * phi ); point_coord_save[1+10*3] = + 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[2+10*3] = 0.0; point_coord_save[0+11*3] = - phi / sqrt ( 1.0 + phi * phi ); point_coord_save[1+11*3] = - 1.0 / sqrt ( 1.0 + phi * phi ); point_coord_save[2+11*3] = 0.0; r8vec_copy ( DIM_NUM * point_num, point_coord_save, point_coord ); i4vec_copy ( EDGE_ORDER * edge_num, edge_point_save, edge_point ); i4vec_copy ( face_num, face_order_save, face_order ); i4vec_copy ( face_order_max * face_num, face_point_save, face_point ); return; # undef DIM_NUM # undef EDGE_NUM # undef EDGE_ORDER # undef FACE_NUM # undef POINT_NUM } /******************************************************************************/ void icos_num ( int *point_num, int *edge_num, int *face_num, int *face_order_max ) /******************************************************************************/ /* Purpose: ICOS_NUM gives "sizes" for an icosahedron. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 July 2007 Author: John Burkardt Parameters: Output, int *POINT_NUM, the number of points. Output, int *EDGE_NUM, the number of edges. Output, int *FACE_NUM, the number of faces. Output, int *FACE_ORDER_MAX, the maximum order of any face. */ { *point_num = 12; *edge_num = 30; *face_num = 20; *face_order_max = 3; return; } /******************************************************************************/ double r8_abs ( double x ) /******************************************************************************/ /* Purpose: R8_ABS returns the absolute value of an R8. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, the quantity whose absolute value is desired. Output, double R8_ABS, the absolute value of X. */ { double value; if ( 0.0 <= x ) { value = + x; } else { value = - x; } return value; } /******************************************************************************/ double r8_modp ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MODP returns the nonnegative remainder of R8 division. Formula: If REM = R8_MODP ( X, Y ) RMULT = ( X - REM ) / Y then X = Y * RMULT + REM where REM is always nonnegative. Discussion: The MOD function computes a result with the same sign as the quantity being divided. Thus, suppose you had an angle A, and you wanted to ensure that it was between 0 and 360. Then mod(A,360.0) would do, if A was positive, but if A was negative, your result would be between -360 and 0. On the other hand, R8_MODP(A,360.0) is between 0 and 360, always. Example: I J MOD R8_MODP R8_MODP Factorization 107 50 7 7 107 = 2 * 50 + 7 107 -50 7 7 107 = -2 * -50 + 7 -107 50 -7 43 -107 = -3 * 50 + 43 -107 -50 -7 43 -107 = 3 * -50 + 43 Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, double X, the number to be divided. Input, double Y, the number that divides X. Output, double R8_MODP, the nonnegative remainder when X is divided by Y. */ { double value; if ( y == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_MODP - Fatal error!\n" ); fprintf ( stderr, " R8_MODP ( X, Y ) called with Y = %f\n", y ); exit ( 1 ); } value = x - ( ( double ) ( ( int ) ( x / y ) ) ) * y; if ( value < 0.0 ) { value = value + fabs ( y ); } return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a unit pseudorandom R8. 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 GNU LGPL 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 k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } /* Although SEED can be represented exactly as a 32 bit integer, it generally cannot be represented exactly as a 32 bit real number! */ r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8mat_transpose_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, char *TITLE, a title. */ { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, int ILO, JLO, the first row and column to print. Input, int IHI, JHI, the last row and column to print. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2; int i2hi; int i2lo; int inc; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); inc = i2hi + 1 - i2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, " Row:" ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, " %7d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j - 1 ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14f", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void r8mat_write ( char *output_filename, int m, int n, double table[] ) /******************************************************************************/ /* Purpose: R8MAT_WRITE writes an R8MAT file. Discussion: An R8MAT is an array of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 June 2009 Author: John Burkardt Parameters: Input, char *OUTPUT_FILENAME, the output filename. Input, int M, the spatial dimension. Input, int N, the number of points. Input, double TABLE[M*N], the data. */ { int i; int j; FILE *output; /* Open the file. */ output = fopen ( output_filename, "wt" ); if ( !output ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_WRITE - Fatal error!\n" ); fprintf ( stderr, " Could not open the output file.\n" ); exit ( 1 ); } /* Write the data. */ for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { fprintf ( output, " %24.16g", table[i+j*m] ); } fprintf ( output, "\n" ); } /* Close the file. */ fclose ( output ); return; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Input, double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ double r8vec_diff_norm ( int n, double a[], double b[] ) /******************************************************************************/ /* Purpose: R8VEC_DIFF_NORM returns the L2 norm of the difference of R8VEC's. Discussion: An R8VEC is a vector of R8's. The vector L2 norm is defined as: R8VEC_NORM_L2 = sqrt ( sum ( 1 <= I <= N ) A(I)**2 ). Licensing: This code is distributed under the GNU LGPL license. Modified: 01 March 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in A. Input, double A[N], B[N], the vectors. Output, double R8VEC_DIFF_NORM, the L2 norm of A - B. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + ( a[i] + b[i] ) * ( a[i] - b[i] ); } value = sqrt ( value ); return value; } /******************************************************************************/ 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 GNU LGPL 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; } /******************************************************************************/ int r8vec_eq ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_EQ is true if every pair of entries in two R8VEC's is equal. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 August 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], A2[N], two vectors to compare. Output, int R8VEC_EQ, is TRUE if every pair of elements A1(I) and A2(I) are equal, and FALSE otherwise. */ { int i; for ( i = 0; i < n; i++ ) { if ( a1[i] != a2[i] ) { return 0; } } return 1; } /******************************************************************************/ 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 GNU LGPL 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 GNU LGPL 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_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 April 2009 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; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } return; } /******************************************************************************/ void sphere_cubed_ijk_to_xyz ( int n, int i, int j, int k, double xyz[3] ) /******************************************************************************/ /* Purpose: SPHERE_CUBED_IJK_TO_XYZ: cubed sphere IJK to XYZ coordinates. Licensing: This code is distributed under the GNU LGPL license. Modified: 11 October 2012 Author: John Burkardt Parameters: Input, int N, the number of sections into which each face of the cube is to be divided. Input, int I, J, K, indices between 0 and N. Normally, at least one of the indices should have the value 0 or N. Output, double XYZ[3], coordinates of the point. */ { double pi = 3.141592653589793; double xc; double xyzn; double yc; double zc; if ( i == 0 ) { xc = -1.0; } else if ( i == n ) { xc = +1.0; } else { xc = tan ( ( double ) ( 2 * i - n ) * 0.25 * pi / ( double ) ( n ) ); } if ( j == 0 ) { yc = -1.0; } else if ( j == n ) { yc = +1.0; } else { yc = tan ( ( double ) ( 2 * j - n ) * 0.25 * pi / ( double ) ( n ) ); } if ( k == 0 ) { zc = -1.0; } else if ( k == n ) { zc = +1.0; } else { zc = tan ( ( double ) ( 2 * k - n ) * 0.25 * pi / ( double ) ( n ) ); } xyzn = sqrt ( xc * xc + yc * yc + zc * zc ); xyz[0] = xc / xyzn; xyz[1] = yc / xyzn; xyz[2] = zc / xyzn; return; } /******************************************************************************/ int sphere_cubed_line_num ( int n ) /******************************************************************************/ /* Purpose: SPHERE_CUBED_LINE_NUM counts lines on a cubed sphere grid. Licensing: This code is distributed under the GNU LGPL license. Modified: 11 October 2012 Author: John Burkardt Parameters: Input, int N, the number of sections into which each face of the cube is to be divided. Output, int LINE_NUM, the number of lines. */ { int line_num; line_num = 0; /* If N = 1, the corners form 12 lines. */ if ( n == 1 ) { line_num = 12; return line_num; } /* If 1 < N, each of 8 corners connects to three neighboring edges. */ else { line_num = line_num + 8 * 3; } /* If 2 < N, then each of the 12 edges includes lines. */ if ( 2 < n ) { line_num = line_num + 12 * ( n - 2 ); } /* Lines that belong to one of the six faces. */ if ( 1 < n ) { line_num = line_num + 6 * 2 * n * ( n - 1 ); } return line_num; } /******************************************************************************/ double *sphere_cubed_lines ( int n, int line_num ) /******************************************************************************/ /* Purpose: SPHERE_CUBED_LINES computes the lines on a cubed sphere grid. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 October 2012 Author: John Burkardt Parameters: Input, int N, the number of sections into which each face of the cube is to be divided. Input, int LINE_NUM, the number of lines. Output, double SPHERE_CUBED_LINES[3*2*LINE_NUM], distinct points on the unit sphere generated by a cubed sphere grid. */ { int i; int j; int l; double *line_data; line_data = ( double * ) malloc ( 3 * 2 * line_num ); l = 0; /* If N = 1, the corners form 12 lines. */ if ( n == 1 ) { sphere_cubed_ijk_to_xyz ( n, 0, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 0, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 0, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 0, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 0, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 0, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 0, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n, n, line_data+0+1*3+l*6 ); l = l + 1; return line_data; } /* If 1 < N, each of 8 corners connects to three neighboring edges. */ else { sphere_cubed_ijk_to_xyz ( n, 0, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 1, 0, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 1, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 0, 1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n-1, 0, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 1, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 0, 1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n-1, n, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n-1, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n, 1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 1, n, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n-1, 0, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n, 1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 1, 0, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 1, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 0, n-1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n-1, 0, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 1, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 0, n-1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n-1, n, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n-1, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, n, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n, n-1, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 1, n, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n-1, n, line_data+0+1*3+l*6 ); l = l + 1; sphere_cubed_ijk_to_xyz ( n, 0, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n, n-1, line_data+0+1*3+l*6 ); l = l + 1; } /* If 2 < N, then each of the 12 edges includes lines. */ if ( 2 < n ) { for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, i, 0, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i+1, 0, 0, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, n, i, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, i+1, 0, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, n-i, n, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n-i-1, n, 0, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, 0, n-i, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n-i-1, 0, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, i, 0, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i+1, 0, n, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, n, i, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, i+1, n, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, n-i, n, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n-i-1, n, n, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, 0, n-i, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n-i-1, n, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, 0, 0, i, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, 0, i+1, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, n, 0, i, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, 0, i+1, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, n, n, i, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, n, i+1, line_data+0+1*3+l*6 ); l = l + 1; } for ( i = 1; i <= n - 2; i++ ) { sphere_cubed_ijk_to_xyz ( n, 0, n, i, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, n, i+1, line_data+0+1*3+l*6 ); l = l + 1; } } /* Lines that belong to one of the six faces. */ if ( 1 < n ) { /* 000, nn0 */ for ( i = 1; i <= n - 1; i++ ) { for ( j = 0; j <= n - 1; j++ ) { sphere_cubed_ijk_to_xyz ( n, i, j, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i, j+1, 0, line_data+0+1*3+l*6 ); l = l + 1; } } for ( j = 1; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { sphere_cubed_ijk_to_xyz ( n, i, j, 0, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i+1, j, 0, line_data+0+1*3+l*6 ); l = l + 1; } } /* 00n, nnn */ for ( i = 1; i <= n - 1; i++ ) { for ( j = 0; j <= n - 1; j++ ) { sphere_cubed_ijk_to_xyz ( n, i, j, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i, j+1, n, line_data+0+1*3+l*6 ); l = l + 1; } } for ( j = 1; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { sphere_cubed_ijk_to_xyz ( n, i, j, n, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i+1, j, n, line_data+0+1*3+l*6 ); l = l + 1; } } /* 000:n0n */ for ( i = 1; i <= n - 1; i++ ) { for ( j = 0; j <= n - 1; j++ ) { sphere_cubed_ijk_to_xyz ( n, i, 0, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i, 0, j+1, line_data+0+1*3+l*6 ); l = l + 1; } } for ( j = 1; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { sphere_cubed_ijk_to_xyz ( n, i, 0, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i+1, 0, j, line_data+0+1*3+l*6 ); l = l + 1; } } /* 0n0:nnn */ for ( i = 1; i <= n - 1; i++ ) { for ( j = 0; j <= n - 1; j++ ) { sphere_cubed_ijk_to_xyz ( n, i, n, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i, n, j+1, line_data+0+1*3+l*6 ); l = l + 1; } } for ( j = 1; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { sphere_cubed_ijk_to_xyz ( n, i, n, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, i+1, n, j, line_data+0+1*3+l*6 ); l = l + 1; } } /* 000:0nn */ for ( i = 1; i <= n - 1; i++ ) { for ( j = 0; j <= n - 1; j++ ) { sphere_cubed_ijk_to_xyz ( n, 0, i, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, i, j+1, line_data+0+1*3+l*6 ); l = l + 1; } } for ( j = 1; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { sphere_cubed_ijk_to_xyz ( n, 0, i, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, 0, i+1, j, line_data+0+1*3+l*6 ); l = l + 1; } } /* n00:nnn */ for ( i = 1; i <= n - 1; i++ ) { for ( j = 0; j <= n - 1; j++ ) { sphere_cubed_ijk_to_xyz ( n, n, i, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, i, j+1, line_data+0+1*3+l*6 ); l = l + 1; } } for ( j = 1; j <= n - 1; j++ ) { for ( i = 0; i <= n - 1; i++ ) { sphere_cubed_ijk_to_xyz ( n, n, i, j, line_data+0+0*3+l*6 ); sphere_cubed_ijk_to_xyz ( n, n, i+1, j, line_data+0+1*3+l*6 ); l = l + 1; } } } if ( l != line_num ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SPHERE_CUBED_LINES - Fatal error!\n" ); fprintf ( stderr, " LINE_NUM = %d\n", line_num ); fprintf ( stderr, " L = %d\n", l ); exit ( 1 ); } return line_data; } /******************************************************************************/ double *sphere_cubed_points ( int n, int ns ) /******************************************************************************/ /* Purpose: SPHERE_CUBED_POINTS computes the points on a cubed sphere grid. Discussion: For a value of N = 3, for instance, each of the 6 cube faces will be divided into 3 sections, so that a single cube face will have (3+1)x(3+1) points: X---X---X---X | 1 | 4 | 7 | X---X---X---X | 2 | 5 | 8 | X---X---X---X | 3 | 6 | 9 | X---X---X---X Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2012 Author: John Burkardt Parameters: Input, int N, the number of sections into which each face of the cube is to be divided. Input, int NS, the number of points. Output, double SPHERE_CUBED_POINTS[3*NS], distinct points on the unit sphere generated by a cubed sphere grid. */ { int ns2; double *xyz; xyz = ( double * ) malloc ( 3 * ns * sizeof ( double ) ); ns2 = 0; /* Bottom full. */ sphere_cubed_points_face ( n, 0, 0, 0, n, n, 0, &ns2, xyz ); /* To avoid repetition, draw the middles as grids of n-2 x n-1 points. */ sphere_cubed_points_face ( n, 0, 0, 1, 0, n-1, n-1, &ns2, xyz ); sphere_cubed_points_face ( n, 0, n, 1, n-1, n, n-1, &ns2, xyz ); sphere_cubed_points_face ( n, n, 1, 1, n, n, n-1, &ns2, xyz ); sphere_cubed_points_face ( n, 1, 0, 1, n, 0, n-1, &ns2, xyz ); /* Top full. */ sphere_cubed_points_face ( n, 0, 0, n, n, n, n, &ns2, xyz ); if ( ns2 != ns ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SPHERE_CUBED_POINTS - Fatal error\n" ); fprintf ( stderr, " Expected to generated NS = %d points.\n", ns ); fprintf ( stderr, " Generated %d points.\n", ns2 ); exit ( 1 ); } return xyz; } /******************************************************************************/ void sphere_cubed_points_face ( int n, int i1, int j1, int k1, int i2, int j2, int k2, int *ns, double xyz[] ) /******************************************************************************/ /* Purpose: SPHERE_CUBED_POINTS_FACE: points on one face of a cubed sphere grid. Discussion: This routine starts with NS = 0, and is called repeatedly to add points for another face. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2012 Author: John Burkardt Parameters: Input, int N, the number of sections into which each face of the cube is to be divided. Input, int I1, J1, K1, I2, J2, K2, the logical indices, between 0 and N, of two corners of the face grid. It is guaranteed that I1 <= I2, J1 <= J2, and K1 <= K2. Input/output, int *NS, the number of points. Input/output, double XYZ[3*NS], distinct points on the unit sphere generated by a cubed sphere grid. */ { int i; int j; int k; double pi = 3.141592653589793; double xyzn; double xc; double yc; double zc; for ( i = i1; i <= i2; i++ ) { if ( i1 < i2 ) { xc = tan ( ( double ) ( 2 * i - n ) * 0.25 * pi / ( double ) ( n ) ); } else if ( i1 == 0 ) { xc = -1.0; } else if ( i1 == n ) { xc = +1.0; } else { xc = 0.0; } for ( j = j1; j <= j2; j++ ) { if ( j1 < j2 ) { yc = tan ( ( double ) ( 2 * j - n ) * 0.25 * pi / ( double ) ( n ) ); } else if ( j1 == 0 ) { yc = -1.0; } else if ( j1 == n ) { yc = +1.0; } else { yc = 0.0; } for ( k = k1; k <= k2; k++ ) { if ( k1 < k2 ) { zc = tan ( ( double ) ( 2 * k - n ) * 0.25 * pi / ( double ) ( n ) ); } else if ( k1 == 0 ) { zc = -1.0; } else if ( k1 == n ) { zc = +1.0; } else { zc = 0.0; } xyzn = sqrt ( xc * xc + yc * yc + zc * zc ); xyz[0+ *ns *3] = xc / xyzn; xyz[1+ *ns *3] = yc / xyzn; xyz[2+ *ns *3] = zc / xyzn; *ns = *ns + 1; } } } return; } /******************************************************************************/ int sphere_cubed_point_num ( int n ) /******************************************************************************/ /* Purpose: SPHERE_CUBED_POINT_NUM counts the points on a cubed sphere grid. Discussion: For a value of N = 3, for instance, each of the 6 cube faces will be divided into 3 sections, so that a single cube face will have (3+1)x(3+1) points: X---X---X---X | 1 | 4 | 7 | X---X---X---X | 2 | 5 | 8 | X---X---X---X | 3 | 6 | 9 | X---X---X---X The number of points is simply (N+1)^3 - (N-1)^3. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2012 Author: John Burkardt Parameters: Input, int N, the number of sections into which each face of the cube is to be divided. Output, int SPHERE_CUBED_POINT_NUM, the number of points. */ { int ns; ns = i4_power ( n + 1, 3 ) - i4_power ( n - 1, 3 ); return ns; } /******************************************************************************/ double sphere_distance_xyz ( double xyz1[3], double xyz2[3] ) /******************************************************************************/ /* Purpose: SPHERE_DISTANCE_XYZ computes great circle distances on a sphere. Discussion: XYZ coordinates are used. We assume the points XYZ1 and XYZ2 lie on the same 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 GNU LGPL license. Modified: 29 August 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 r; double top; r = r8vec_norm ( 3, xyz1 ); lat1 = arc_sine ( xyz1[2] ); lon1 = atan4 ( xyz1[1], xyz1[0] ); lat2 = arc_sine ( xyz2[2] ); lon2 = atan4 ( 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 = r * atan2 ( top, bot ); return dist; } /******************************************************************************/ int *sphere_grid_q4 ( int lat_num, int long_num ) /******************************************************************************/ /* Purpose: SPHERE_GRID_Q4: rectangular grid on a sphere. Discussion: The point numbering system is the same used in SPHERE_GRIDPOINTS, and that routine may be used to compute the coordinates of the points. A sphere in 3D satisfies the equation: sum ( ( P(1:DIM_NUM) - PC(1:DIM_NUM) )^2 ) = R * R Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2010 Author: John Burkardt Parameters: Input, int LAT_NUM, the number of "rows" of rectangles to be created. LAT_NUM must be at least 2. Input, int LONG_NUM, the number of "columns" of rectangles to be created. Output, int SPHERE_GRID_Q4[4*LAT_NUM*LONG_NUM], the indices of the nodes that make up each rectangle. */ { int i; int j; int n; int n_max; int n_min; int ne; int nw; int s; int s_max; int s_min; int se; int sw; int *rectangle_node; int rectangle_num; rectangle_node = ( int * ) malloc ( ( 4 * lat_num * long_num ) * sizeof ( int ) ); rectangle_num = 0; /* The first row. */ n = 0; sw = 1; se = sw + 1; s_min = 1; s_max = long_num; for ( j = 1; j <= long_num; j++ ) { rectangle_node[0+rectangle_num*4] = sw; rectangle_node[1+rectangle_num*4] = se; rectangle_node[2+rectangle_num*4] = n; rectangle_node[3+rectangle_num*4] = n; rectangle_num = rectangle_num + 1; sw = se; if ( se == s_max ) { se = s_min; } else { se = se + 1; } } /* The intermediate rows. */ for ( i = 2; i < lat_num; i++ ) { n_max = s_max; n_min = s_min; s_max = s_max + long_num; s_min = s_min + long_num; nw = n_min; ne = nw + 1; sw = s_min; se = sw + 1; for ( j = 1; j <= long_num; j++ ) { rectangle_node[0+rectangle_num*4] = sw; rectangle_node[1+rectangle_num*4] = se; rectangle_node[2+rectangle_num*4] = ne; rectangle_node[3+rectangle_num*4] = nw; rectangle_num = rectangle_num + 1; sw = se; nw = ne; if ( se == s_max ) { se = s_min; } else { se = se + 1; } if ( ne == n_max ) { ne = n_min; } else { ne = ne + 1; } } } /* The last row. */ n_max = s_max; n_min = s_min; s = n_max + 1; nw = n_min; ne = nw + 1; for ( j = 1; j <= long_num; j++ ) { rectangle_node[0+rectangle_num*4] = ne; rectangle_node[1+rectangle_num*4] = nw; rectangle_node[2+rectangle_num*4] = s; rectangle_node[3+rectangle_num*4] = s; rectangle_num = rectangle_num + 1; nw = ne; if ( ne == n_max ) { ne = n_min; } else { ne = ne + 1; } } return rectangle_node; } /******************************************************************************/ int *sphere_grid_t3 ( int lat_num, int long_num ) /******************************************************************************/ /* Purpose: SPHERE_GRID_T3 produces a triangle grid on a sphere. Discussion: The point numbering system is the same used in SPHERE_GRIDPOINTS, and that routine may be used to compute the coordinates of the points. A sphere in 3D satisfies the equation: sum ( ( P(1:DIM_NUM) - pc(1:DIM_NUM) )^2 ) = R * R Licensing: This code is distributed under the GNU LGPL license. Modified: 28 August 2010 Author: John Burkardt Parameters: Input, int LAT_NUM, LONG_NUM, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so LAT_NUM = 5, for instance, will result in points along 7 lines of latitude. Output, int SPHERE_GRID_T3[3*2*(LAT_NUM+1)*LONG_NUM], the triangle vertices. */ { int i; int j; int n; int n_max; int n_min; int ne; int nw; int s; int s_max; int s_min; int se; int sw; int *triangle_node; int triangle_num; triangle_node = ( int * ) malloc ( 3 * 2 * ( lat_num + 1 ) * long_num * sizeof ( int ) ); triangle_num = 0; /* The first row. */ /* Not working!: n = 0; sw = 1; se = sw + 1; s_min = 1; s_max = long_num; */ n = 1; sw = 2; se = sw + 1; s_min = 2; s_max = long_num + 1; for ( j = 0; j < long_num; j++ ) { triangle_node[0+triangle_num*3] = sw - 1; triangle_node[1+triangle_num*3] = se - 1; triangle_node[2+triangle_num*3] = n - 1; triangle_num = triangle_num + 1; sw = se; if ( se == s_max ) { se = s_min; } else { se = se + 1; } } /* The intermediate rows. */ for ( i = 1; i <= lat_num; i++ ) { n_max = s_max; n_min = s_min; s_max = s_max + long_num; s_min = s_min + long_num; nw = n_min; ne = nw + 1; sw = s_min; se = sw + 1; for ( j = 0; j < long_num; j++ ) { triangle_node[0+triangle_num*3] = sw - 1; triangle_node[1+triangle_num*3] = se - 1; triangle_node[2+triangle_num*3] = nw - 1; triangle_num = triangle_num + 1; triangle_node[0+triangle_num*3] = ne - 1; triangle_node[1+triangle_num*3] = nw - 1; triangle_node[2+triangle_num*3] = se - 1; triangle_num = triangle_num + 1; sw = se; nw = ne; if ( se == s_max ) { se = s_min; } else { se = se + 1; } if ( ne == n_max ) { ne = n_min; } else { ne = ne + 1; } } } /* The last row. */ n_max = s_max; n_min = s_min; s = n_max + 1; nw = n_min; ne = nw + 1; for ( j = 0; j < long_num; j++ ) { triangle_node[0+triangle_num*3] = ne - 1; triangle_node[1+triangle_num*3] = nw - 1; triangle_node[2+triangle_num*3] = s - 1; triangle_num = triangle_num + 1; nw = ne; if ( ne == n_max ) { ne = n_min; } else { ne = ne + 1; } } return triangle_node; } /******************************************************************************/ int sphere_icos_edge_num ( int factor ) /******************************************************************************/ /* Purpose: SPHERE_ICOS_EDGE_NUM sizes an icosahedral grid on a sphere. Discussion: With FACTOR = 1, the grid has 20 triangular faces, 30 edges, and 12 nodes. With FACTOR = 2, each triangle of the icosahedron is subdivided into 2x2 subtriangles, resulting in 80 faces, 120 edges, and 42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes. With FACTOR = 3, each triangle of the icosahedron is subdivided into 3x3 subtriangles, resulting in 180 faces, 270 edges and 92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes. In general, each triangle is subdivided into FACTOR*FACTOR subtriangles, resulting in 20 * FACTOR * FACTOR faces, 30 * FACTOR * FACTOR edges, and 12 + 20 * 3 * (FACTOR-1) / 2 + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2010 Author: John Burkardt Parameters: Input, int FACTOR, the subdivision factor, which must be at least 1. Output, int SPHERE_ICOS_EDGE_NUM, the number of edges. */ { int edge_num; edge_num = 30 * factor * factor; return edge_num; } /******************************************************************************/ int sphere_icos_face_num ( int factor ) /******************************************************************************/ /* Purpose: SPHERE_ICOS_FACE_NUM sizes an icosahedral grid on a sphere. Discussion: With FACTOR = 1, the grid has 20 triangular faces, 30 edges, and 12 nodes. With FACTOR = 2, each triangle of the icosahedron is subdivided into 2x2 subtriangles, resulting in 80 faces, 120 edges, and 42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes. With FACTOR = 3, each triangle of the icosahedron is subdivided into 3x3 subtriangles, resulting in 180 faces, 270 edges and 92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes. In general, each triangle is subdivided into FACTOR*FACTOR subtriangles, resulting in 20 * FACTOR * FACTOR faces, 30 * FACTOR * FACTOR edges, and 12 + 20 * 3 * (FACTOR-1) / 2 + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2010 Author: John Burkardt Parameters: Input, int FACTOR, the subdivision factor, which must be at least 1. Output, int SPHERE_ICOS_FACE_NUM, the number of triangles. */ { int face_num; face_num = 20 * factor * factor; return face_num; } /******************************************************************************/ int sphere_icos_point_num ( int factor ) /******************************************************************************/ /* Purpose: SPHERE_ICOS_POINT_NUM sizes an icosahedral grid on a sphere. Discussion: With FACTOR = 1, the grid has 20 triangular faces, 30 edges, and 12 nodes. With FACTOR = 2, each triangle of the icosahedron is subdivided into 2x2 subtriangles, resulting in 80 faces, 120 edges, and 42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes. With FACTOR = 3, each triangle of the icosahedron is subdivided into 3x3 subtriangles, resulting in 180 faces, 270 edges and 92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes. In general, each triangle is subdivided into FACTOR*FACTOR subtriangles, resulting in 20 * FACTOR * FACTOR faces, 30 * FACTOR * FACTOR edges, and 12 + 20 * 3 * (FACTOR-1) / 2 + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2010 Author: John Burkardt Parameters: Input, int FACTOR, the subdivision factor, which must be at least 1. Output, int SPHERE_ICOS_NODE_NUM, the number of nodes. */ { int node_num; node_num = 12 + 10 * 3 * ( factor - 1 ) + 10 * ( factor - 2 ) * ( factor - 1 ); return node_num; } /******************************************************************************/ double *sphere_icos1_points ( int factor, int node_num ) /******************************************************************************/ /* Purpose: SPHERE_ICOS1_POINTS returns icosahedral grid points on a sphere. Discussion: With FACTOR = 1, the grid has 20 triangular faces and 12 nodes. With FACTOR = 2, each triangle of the icosahedron is subdivided into 2x2 subtriangles, resulting in 80 faces and 42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes. With FACTOR = 3, each triangle of the icosahedron is subdivided into 3x3 subtriangles, resulting in 180 faces and 92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes. In general, each triangle is subdivided into FACTOR*FACTOR subtriangles, resulting in 20 * FACTOR * FACTOR faces and 12 + 20 * 3 * (FACTOR-1) / 2 + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 July 2007 Author: John Burkardt Parameters: Input, int FACTOR, the subdivision factor, which must be at least 1. Input, int NODE_NUM, the number of nodes, as reported by SPHERE_GRID_ICOS_NUM. Output, double SPHERE_GRIDPOINTS_ICOS1[3*NODE_NUM], the node coordinates. Local Parameters: POINT_NUM, EDGE_NUM, FACE_NUM and FACE_ORDER_MAX are counters associated with the icosahedron, and POINT_COORD, EDGE_POINT, FACE_ORDER and FACE_POINT are data associated with the icosahedron. We need to refer to this data to generate the grid. NODE counts the number of nodes we have generated so far. At the end of the routine, it should be equal to NODE_NUM. */ { int a; int b; int c; int dim; int edge; int edge_num; int *edge_point; int f; int f1; int f2; int face; int face_num; int *face_order; int *face_point; int face_order_max; int node; double node_norm; double *node_xyz; int point; double *point_coord; int point_num; node_xyz = ( double * ) malloc ( 3 * node_num * sizeof ( double ) ); /* Size the icosahedron. */ icos_num ( &point_num, &edge_num, &face_num, &face_order_max ); /* Set the icosahedron. */ point_coord = ( double * ) malloc ( 3 * point_num * sizeof ( double ) ); edge_point = ( int * ) malloc ( 2 * edge_num * sizeof ( int ) ); face_order = ( int * ) malloc ( face_num * sizeof ( int ) ); face_point = ( int * ) malloc ( face_order_max * face_num * sizeof ( int ) ); icos_shape ( point_num, edge_num, face_num, face_order_max, point_coord, edge_point, face_order, face_point ); /* Generate the point coordinates. A. Points that are the icosahedral vertices. */ node = 0; for ( point = 0; point < point_num; point++ ) { for ( dim = 0; dim < 3; dim++ ) { node_xyz[dim+node*3] = point_coord[dim+point*3]; } node = node + 1; } /* B. Points in the icosahedral edges, at 1/FACTOR, 2/FACTOR, ..., (FACTOR-1)/FACTOR. */ for ( edge = 0; edge < edge_num; edge++ ) { a = edge_point[0+edge*2]; b = edge_point[1+edge*2]; for ( f = 1; f < factor; f++ ) { for ( dim = 0; dim < 3; dim++ ) { node_xyz[dim+node*3] = ( ( double ) ( factor - f ) * point_coord[dim+a*3] + ( double ) ( f ) * point_coord[dim+b*3] ) / ( double ) ( factor ); } node_norm = r8vec_norm ( 3, node_xyz+node*3 ); for ( dim = 0; dim < 3; dim++ ) { node_xyz[dim+node*3] = node_xyz[dim+node*3] / node_norm; } node = node + 1; } } /* C. Points in the icosahedral faces. */ for ( face = 0; face < face_num; face++ ) { a = face_point[0+face*3]; b = face_point[1+face*3]; c = face_point[2+face*3]; for ( f1 = 1; f1 < factor; f1++ ) { for ( f2 = 1; f2 < factor - f1; f2++ ) { for ( dim = 0; dim < 3; dim++ ) { node_xyz[dim+node*3] = ( ( double ) ( factor - f1 - f2 ) * point_coord[dim+a*3] + ( double ) ( f1 ) * point_coord[dim+b*3] + ( double ) ( f2 ) * point_coord[dim+c*3] ) / ( double ) ( factor ); } node_norm = r8vec_norm ( 3, node_xyz+node*3 ); for ( dim = 0; dim < 3; dim++ ) { node_xyz[dim+node*3] = node_xyz[dim+node*3] / node_norm; } node = node + 1; } } } /* Discard allocated memory. */ free ( edge_point ); free ( face_order ); free ( face_point ); free ( point_coord ); return node_xyz; } /******************************************************************************/ double *sphere_icos2_points ( int factor, int node_num ) /******************************************************************************/ /* Purpose: SPHERE_ICOS2_POINTS returns icosahedral grid points on a sphere. Discussion: With FACTOR = 1, the grid has 20 triangular faces and 12 nodes. With FACTOR = 2, each triangle of the icosahedron is subdivided into 2x2 subtriangles, resulting in 80 faces and 42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes. With FACTOR = 3, each triangle of the icosahedron is subdivided into 3x3 subtriangles, resulting in 180 faces and 92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes. In general, each triangle is subdivided into FACTOR*FACTOR subtriangles, resulting in 20 * FACTOR * FACTOR faces and 12 + 20 * 3 * (FACTOR-1) / 2 + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2010 Author: John Burkardt Parameters: Input, int FACTOR, the subdivision factor, which must be at least 1. Input, int NODE_NUM, the number of nodes, as reported by SPHERE_IMP_GRID_ICOS_NUM. Output, double SPHERE_GRIDPOINTS_ICOS2[3*NODE_NUM], the node coordinates. Local Parameters: POINT_NUM, EDGE_NUM, FACE_NUM and FACE_ORDER_MAX are counters associated with the icosahedron, and POINT_COORD, EDGE_POINT, FACE_ORDER and FACE_POINT are data associated with the icosahedron. We need to refer to this data to generate the grid. NODE counts the number of nodes we have generated so far. At the end of the routine, it should be equal to NODE_NUM. */ { int a; double angle; double ab[3]; double ac[3]; double acn[3]; double acn_norm; double acp[3]; int b; double bn[3]; double bn_norm; double bp[3]; int c; double cn[3]; double cn_norm; double cp[3]; int edge; int edge_num; int *edge_point; int f; int fa; int fbc; int face; int face_num; int *face_order; int *face_point; int face_order_max; int i; int j; int node; double *node_xyz; double *point_coord; int point_num; double theta; double theta_ab; double theta_ac; double theta_bc; node_xyz = ( double * ) malloc ( 3 * node_num * sizeof ( double ) ); /* Size the icosahedron. */ icos_num ( &point_num, &edge_num, &face_num, &face_order_max ); /* Set the icosahedron. */ point_coord = ( double * ) malloc ( 3 * point_num * sizeof ( double ) ); edge_point = ( int * ) malloc ( 2 * edge_num * sizeof ( int ) ); face_order = ( int * ) malloc ( face_num * sizeof ( int ) ); face_point = ( int * ) malloc ( face_order_max * face_num * sizeof ( int ) ); icos_shape ( point_num, edge_num, face_num, face_order_max, point_coord, edge_point, face_order, face_point ); /* Generate the point coordinates. A. Points that are the icosahedral vertices. */ node = 0; for ( j = 0; j < point_num; j++ ) { for ( i = 0; i < 3; i++ ) { node_xyz[i+j*3] = point_coord[i+j*3]; } node = node + 1; } /* B. Points in the icosahedral edges, at 1/FACTOR, 2/FACTOR, ..., (FACTOR-1)/FACTOR. */ for ( edge = 0; edge < edge_num; edge++ ) { a = edge_point[0+edge*2]; b = edge_point[1+edge*2]; /* Determine the "distance" = angle between points A and B. */ theta = sphere_distance_xyz ( point_coord+a*3, point_coord+b*3 ); /* Polarize B into BP + BN and normalize BN. */ r8vec_polarize ( 3, point_coord+b*3, point_coord+a*3, bn, bp ); bn_norm = r8vec_norm ( 3, bn ); for ( i = 0; i < 3; i++ ) { bn[i] = bn[i] / bn_norm; } /* March from A to B, by taking equally spaced angles from 0 to THETA. F = 0 => ANGLE = 0 => A F = FACTOR => ANGLE = THETA => B */ for ( f = 1; f < factor; f++ ) { angle = ( ( double ) ( f ) * theta ) / ( double ) ( factor ); for ( i = 0; i < 3; i++ ) { node_xyz[i+node*3] = cos ( angle ) * point_coord[i+a*3] + sin ( angle ) * bn[i]; } node = node + 1; } } /* C. Points in the icosahedral faces. */ for ( face = 0; face < face_num; face++ ) { a = face_point[0+face*3]; b = face_point[1+face*3]; c = face_point[2+face*3]; /* Determine the "distance" = angle between points A and B, A and C. */ theta_ab = sphere_distance_xyz ( point_coord+a*3, point_coord+b*3 ); theta_ac = sphere_distance_xyz ( point_coord+a*3, point_coord+c*3 ); /* Polarize B = BP + BN and normalize BN, C = CP + CN, and normalize CN. */ r8vec_polarize ( 3, point_coord+b*3, point_coord+a*3, bn, bp ); bn_norm = r8vec_norm ( 3, bn ); for ( i = 0; i < 3; i++ ) { bn[i] = bn[i] / bn_norm; } r8vec_polarize ( 3, point_coord+c*3, point_coord+a*3, cn, cp ); cn_norm = r8vec_norm ( 3, cn ); for ( i = 0; i < 3; i++ ) { cn[i] = cn[i] / cn_norm; } /* March AB from A to B: FA = 0 => ANGLE = 0 => AB = A FA = FACTOR => ANGLE = THETA_AB => AB = B March AC from A to C: FA = 0 => ANGLE = 0 => AC = A FA = FACTOR => ANGLE = THETA_AC => AC = C */ for ( fa = 2; fa < factor; fa++ ) { /* Determine points AB and AC that use cos ( FA / FACTOR ) of A and cos ( ( FACTOR - FA ) / FACTOR ) of B or C. */ angle = ( (double ) ( fa ) * theta_ab ) / ( double ) ( factor ); for ( i = 0; i < 3; i++ ) { ab[i] = cos ( angle ) * point_coord[i+a*3] + sin ( angle ) * bn[i]; } angle = ( ( double ) ( fa ) * theta_ac ) / ( double ) ( factor ); for ( i = 0; i < 3; i++ ) { ac[i] = cos ( angle ) * point_coord[i+a*3] + sin ( angle ) * cn[i]; } /* Determine the "distance" = angle between points AB and AC. */ theta_bc = sphere_distance_xyz ( ab, ac ); /* Polarize AC into ACP + ACN and normalize ACN. */ r8vec_polarize ( 3, ac, ab, acn, acp ); acn_norm = r8vec_norm ( 3, acn ); for ( i = 0; i < 3; i++ ) { acn[i] = acn[i] / acn_norm; } /* The interval between AB and AC is broken into FA intervals. Go from 1 to FA - 1. */ for ( fbc = 1; fbc < fa; fbc++ ) { angle = ( ( double ) ( fbc ) * theta_bc ) / ( double ) ( fa ); for ( i = 0; i < 3; i++ ) { node_xyz[i+node*3] = cos ( angle ) * ab[i] + sin ( angle ) * acn[i]; } node = node + 1; } } } /* Discard allocated memory. */ free ( edge_point ); free ( face_order ); free ( face_point ); free ( point_coord ); return node_xyz; } /******************************************************************************/ int sphere_line_project ( double r, double pc[3], int n, double p[], int maxpnt2, double pp[], double thetamin, double thetamax ) /******************************************************************************/ /* Purpose: SPHERE_LINE_PROJECT projects a line onto an implicit sphere. Discussion: The implicit form of a sphere in 3D is: pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) + pow ( P[2] - PC[2], 2 ) = pow ( R, 2 ) The line to be projected is specified as a sequence of points. If two successive points subtend a small angle, then the second point is essentially dropped. If two successive points subtend a large angle, then intermediate points are inserted, so that the projected line stays closer to the sphere. Note that if any P coincides with the center of the sphere, then its projection is mathematically undefined. P will be returned as PC. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 July 2005 Author: John Burkardt Parameters: Input, double R, the radius of the sphere. If R is zero, PP will be returned as PC, and if R is negative, points will end up diametrically opposite from where you would expect them for a positive R. Input, double PC[3], the coordinates of the center of the sphere. Input, int N, the number of points on the line that is to be projected. Input, double P[3*N], the coordinates of the points on the line that is to be projected. Input, int MAXPNT2, the maximum number of points on the projected line. Even if the routine thinks that more points are needed, no more than MAXPNT2 will be generated. Output, double PP[3*N2], the coordinates of the points representing the projected line. The value N2 is returned as the function value of this routine. These points lie on the sphere. Successive points are separated by at least THETAMIN radians, and by no more than THETAMAX radians. Input, double THETAMIN, THETAMAX, the minimum and maximum angular projections allowed between successive projected points. If two successive points on the original line have projections separated by more than THETAMAX radians, then intermediate points will be inserted, in an attempt to keep the line closer to the sphere. If two successive points are separated by less than THETAMIN radians, then the second point is dropped, and the line from the first point to the next point is considered. Output, int SPHERE_LINE_PROJECT, the number of points on the projected line. This value can be zero, if the line has an angular projection of less than THETAMIN radians. */ { # define DIM_NUM 3 double alpha; double ang3d; double dot; int i; int j; int nfill; int n2; double tnorm; double p1[DIM_NUM]; double p2[DIM_NUM]; double pi[DIM_NUM]; /* Check the input. */ if ( r == 0.0 ) { n2 = 0; return n2; } r8vec_copy ( DIM_NUM, pc, p1 ); r8vec_copy ( DIM_NUM, pc, p2 ); n2 = 0; for ( i = 0; i < n; i++ ) { if ( r8vec_eq ( DIM_NUM, p, pc ) ) { } else { r8vec_copy ( DIM_NUM, p2, p1 ); alpha = sqrt ( pow ( p[0+i*3] - pc[0], 2 ) + pow ( p[1+i*3] - pc[1], 2 ) + pow ( p[2+i*3] - pc[2], 2 ) ); p2[0] = pc[0] + r * ( p[0+i*3] - pc[0] ) / alpha; p2[1] = pc[1] + r * ( p[1+i*3] - pc[1] ) / alpha; p2[2] = pc[2] + r * ( p[2+i*3] - pc[2] ) / alpha; /* If we haven't gotten any points yet, take this point as our start. /*/ if ( n2 == 0 ) { pp[0+n2*3] = p2[0]; pp[1+n2*3] = p2[1]; pp[2+n2*3] = p2[2]; n2 = n2 + 1; } /* Compute the angular projection of P1 to P2. */ else if ( 1 <= n2 ) { dot = ( p1[0] - pc[0] ) * ( p2[0] - pc[0] ) + ( p1[1] - pc[1] ) * ( p2[1] - pc[1] ) + ( p1[2] - pc[2] ) * ( p2[2] - pc[2] ); ang3d = arc_cosine ( dot / ( r * r ) ); /* If the angle is at least THETAMIN, (or it's the last point), then we will draw a line segment. */ if ( thetamin < r8_abs ( ang3d ) || i == n ) { /* Now we check to see if the line segment is too long. */ if ( thetamax < r8_abs ( ang3d ) ) { nfill = ( int ) ( r8_abs ( ang3d ) / thetamax ); for ( j = 1; j < nfill; j++ ) { pi[0] = ( ( double ) ( nfill - j ) * ( p1[0] - pc[0] ) + ( double ) ( j ) * ( p2[0] - pc[0] ) ); pi[1] = ( ( double ) ( nfill - j ) * ( p1[1] - pc[1] ) + ( double ) ( j ) * ( p2[1] - pc[1] ) ); pi[2] = ( ( double ) ( nfill - j ) * ( p1[2] - pc[2] ) + ( double ) ( j ) * ( p2[2] - pc[2] ) ); tnorm = r8vec_norm ( DIM_NUM, pi ); if ( tnorm != 0.0 ) { pi[0] = pc[0] + r * pi[0] / tnorm; pi[1] = pc[1] + r * pi[1] / tnorm; pi[2] = pc[2] + r * pi[2] / tnorm; pp[0+n2*3] = pi[0]; pp[1+n2*3] = pi[1]; pp[2+n2*3] = pi[2]; n2 = n2 + 1; } } } /* Now tack on the projection of point 2. */ pp[0+n2*3] = p2[0]; pp[1+n2*3] = p2[1]; pp[2+n2*3] = p2[2]; n2 = n2 + 1; } } } } return n2; # undef DIM_NUM } /******************************************************************************/ int *sphere_ll_lines ( int nlat, int nlong, int line_num ) /******************************************************************************/ /* Purpose: SPHERE_LL_LINES produces lines on a latitude/longitude grid. Discussion: The point numbering system is the same used in SPHERE_LL_POINTS, and that routine may be used to compute the coordinates of the points. The implicit form of a sphere in 3D is: pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) + pow ( P[2] - PC[2], 2 ) = pow ( R, 2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 08 October 2012 Author: John Burkardt Parameters: Input, int NLAT, NLONG, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so NLAT = 5, for instance, will result in points along 7 lines of latitude. Input, int LINE_NUM, the number of grid lines. Output, int SPHERE_LL_LINES[2*LINE_NUM], contains pairs of point indices for line segments that make up the grid. */ { int i; int j; int l; int *line; int next; int newcol; int old; line = ( int * ) malloc ( 2 * line_num * sizeof ( int ) ); l = 0; /* "Vertical" lines. */ for ( j = 0; j <= nlong - 1; j++ ) { old = 1; next = j + 2; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; for ( i = 1; i <= nlat-1; i++ ) { old = next; next = old + nlong; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } old = next; line[0+l*2] = old; line[1+l*2] = 1 + nlat * nlong + 1; l = l + 1; } /* "Horizontal" lines. */ for ( i = 1; i <= nlat; i++ ) { next = 1 + ( i - 1 ) * nlong + 1; for ( j = 0; j <= nlong-2; j++ ) { old = next; next = old + 1; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } old = next; next = 1 + ( i - 1 ) * nlong + 1; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } /* "Diagonal" lines. */ for ( j = 0; j <= nlong - 1; j++ ) { old = 1; next = j + 2; newcol = j; for ( i = 1; i <= nlat - 1; i++ ) { old = next; next = old + nlong + 1; newcol = newcol + 1; if ( nlong - 1 < newcol ) { newcol = 0; next = next - nlong; } line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } } return line; } /******************************************************************************/ int sphere_ll_line_num ( int lat_num, int long_num ) /******************************************************************************/ /* Purpose: SPHERE_LL_LINE_NUM counts lines for a latitude/longitude grid. Discussion: The number returned is the number of pairs of points to be connected. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 October 2012 Author: John Burkardt Parameters: Input, int LAT_NUM, LONG_NUM, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so LAT_NUM = 5, for instance, will result in points along 7 lines of latitude. Output, int SPHERE_LL_LINE_NUM, the number of grid lines. */ { int line_num; line_num = long_num * ( lat_num + 1 ) + lat_num * long_num + long_num * ( lat_num - 1 ); return line_num; } /******************************************************************************/ double *sphere_ll_points ( double r, double pc[3], int lat_num, int lon_num, int point_num ) /******************************************************************************/ /* Purpose: SPHERE_LL_POINTS produces points on a latitude/longitude grid. Discussion: The implicit form of a sphere in 3D is: pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) + pow ( P[2] - PC[2], 2 ) = pow ( R, 2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 08 October 2012 Author: John Burkardt Parameters: Input, double R, the radius of the sphere. Input, double PC[3], the coordinates of the center of the sphere. Input, int LAT_NUM, LON_NUM, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so LAT_NUM = 5, for instance, will result in points along 7 lines of latitude. Input, int POINT_NUM, the number of points. Output, double SPHERE_LL_POINTS[3*POINT_NUM], the coordinates of the grid points. */ { int lat; int lon; int n; double *p; double phi; double pi = 3.141592653589793; double theta; p = ( double * ) malloc ( 3 * point_num * sizeof ( double ) ); n = 0; /* The north pole. */ theta = 0.0; phi = 0.0; p[0+n*3] = pc[0] + r * sin ( phi ) * cos ( theta ); p[1+n*3] = pc[1] + r * sin ( phi ) * sin ( theta ); p[2+n*3] = pc[2] + r * cos ( phi ); n = n + 1; /* Do each intermediate ring of latitude. */ for ( lat = 1; lat <= lat_num; lat++ ) { phi = ( double ) ( lat ) * pi / ( double ) ( lat_num + 1 ); /* Along that ring of latitude, compute points at various longitudes. */ for ( lon = 0; lon < lon_num; lon++ ) { theta = ( double ) ( lon ) * 2.0 * pi / ( double ) ( lon_num ); p[0+n*3] = pc[0] + r * sin ( phi ) * cos ( theta ); p[1+n*3] = pc[1] + r * sin ( phi ) * sin ( theta ); p[2+n*3] = pc[2] + r * cos ( phi ); n = n + 1; } } /* The south pole. */ theta = 0.0; phi = pi; p[0+n*3] = pc[0] + r * sin ( phi ) * cos ( theta ); p[1+n*3] = pc[1] + r * sin ( phi ) * sin ( theta ); p[2+n*3] = pc[2] + r * cos ( phi ); n = n + 1; return p; } /******************************************************************************/ int sphere_ll_point_num ( int lat_num, int long_num ) /******************************************************************************/ /* Purpose: SPHERE_LL_POINT_NUM counts points for a latitude/longitude grid. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 October 2012 Author: John Burkardt Parameters: Input, int LAT_NUM, LONG_NUM, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so LAT_NUM = 5, for instance, will result in points along 7 lines of latitude. Output, int SPHERE_LL_POINT_NUM, the number of grid points. */ { int point_num; point_num = 2 + lat_num * long_num; return point_num; } /******************************************************************************/ int *sphere_llq_lines ( int nlat, int nlong, int line_num ) /******************************************************************************/ /* Purpose: SPHERE_LLQ_LINES: latitude/longitude quadrilateral grid lines. Discussion: The point numbering system is the same used in SPHERE_LL_POINTS, and that routine may be used to compute the coordinates of the points. The implicit form of a sphere in 3D is: pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) + pow ( P[2] - PC[2], 2 ) = pow ( R, 2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 08 October 2012 Author: John Burkardt Parameters: Input, int NLAT, NLONG, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so NLAT = 5, for instance, will result in points along 7 lines of latitude. Input, int LINE_NUM, the number of grid lines. Output, int SPHERE_LLQ_LINES[2*LINE_NUM], contains pairs of point indices for line segments that make up the grid. */ { int i; int j; int l; int *line; int next; int old; line = ( int * ) malloc ( 2 * line_num * sizeof ( int ) ); l = 0; /* "Vertical" lines. */ for ( j = 0; j <= nlong - 1; j++ ) { old = 1; next = j + 2; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; for ( i = 1; i <= nlat-1; i++ ) { old = next; next = old + nlong; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } old = next; line[0+l*2] = old; line[1+l*2] = 1 + nlat * nlong + 1; l = l + 1; } /* "Horizontal" lines. */ for ( i = 1; i <= nlat; i++ ) { next = 1 + ( i - 1 ) * nlong + 1; for ( j = 0; j <= nlong-2; j++ ) { old = next; next = old + 1; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } old = next; next = 1 + ( i - 1 ) * nlong + 1; line[0+l*2] = old; line[1+l*2] = next; l = l + 1; } return line; } /******************************************************************************/ int sphere_llq_line_num ( int lat_num, int long_num ) /******************************************************************************/ /* Purpose: SPHERE_LLQ_LINE_NUM counts lines for a latitude/longitude quadrilateral grid. Discussion: The number returned is the number of pairs of points to be connected. Licensing: This code is distributed under the GNU LGPL license. Modified: 09 October 2012 Author: John Burkardt Parameters: Input, int LAT_NUM, LONG_NUM, the number of latitude and longitude lines to draw. The latitudes do not include the North and South poles, which will be included automatically, so LAT_NUM = 5, for instance, will result in points along 7 lines of latitude. Output, int SPHERE_LLQ_LINE_NUM, the number of grid lines. */ { int line_num; line_num = long_num * ( lat_num + 1 ) + lat_num * long_num; return line_num; } /******************************************************************************/ double *sphere_spiralpoints ( double r, double pc[3], int n ) /******************************************************************************/ /* Purpose: SPHERE_SPIRALPOINTS produces spiral points on an implicit sphere. Discussion: The points should be arranged on the sphere in a pleasing design. The implicit form of a sphere in 3D is: pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) + pow ( P[2] - PC[2], 2 ) = pow ( R, 2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 03 July 2005 Author: John Burkardt Reference: Edward Saff, Arno Kuijlaars, Distributing Many Points on a Sphere, The Mathematical Intelligencer, Volume 19, Number 1, 1997, pages 5-11. Parameters: Input, double R, the radius of the sphere. Input, double PC[3], the coordinates of the center of the sphere. Input, int N, the number of points to create. Output, double SPHERE_SPIRALPOINTS[3*N], the coordinates of the grid points. */ { double cosphi = 0.0; int i; double *p; double pi = 3.141592653589793; double sinphi = 0.0; double theta = 0.0; p = ( double * ) malloc ( 3 * n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { cosphi = ( ( double ) ( n - i - 1 ) * ( -1.0 ) + ( double ) ( i ) * ( 1.0 ) ) / ( double ) ( n - 1 ); sinphi = sqrt ( 1.0 - cosphi * cosphi ); if ( i == 0 || i == n - 1 ) { theta = 0.0; } else { theta = theta + 3.6 / ( sinphi * sqrt ( ( double ) n ) ); theta = r8_modp ( theta, 2.0 * pi ); } p[0+i*3] = pc[0] + r * sinphi * cos ( theta ); p[1+i*3] = pc[1] + r * sinphi * sin ( theta ); p[2+i*3] = pc[2] + r * cosphi; } return p; } /******************************************************************************/ double *sphere_unit_sample ( int n, int *seed ) /******************************************************************************/ /* Purpose: SPHERE_UNIT_SAMPLE picks a random point on the unit sphere. Discussion: The unit sphere in 3D satisfies the equation: X * X + Y * Y + Z * Z = 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2010 Author: John Burkardt Parameters: Input, int N, the number of points to generate. Input/output, int *SEED, a seed for the random number generator. Output, double SPHERE_UNIT_SAMPLE[3*N], the sample point. */ { int j; double phi; double 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. 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 = 2.0 * r8_uniform_01 ( seed ) - 1.0; phi = arc_cosine ( vdot ); /* Pick a uniformly random rotation between 0 and 2 Pi around the axis of the Z vector. */ theta = 2.0 * pi * r8_uniform_01 ( seed ); x[0+j*3] = cos ( theta ) * sin ( phi ); x[1+j*3] = sin ( theta ) * sin ( phi ); x[2+j*3] = cos ( phi ); } return x; } /******************************************************************************/ 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 GNU LGPL license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_NUM 40 static char time_buffer[TIME_NUM]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_NUM, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_NUM }