Sun position
Sun position algorithm
vectors3d.c
1 /*==============================================================================
2  * vectors3d.c - Three dimensional geometry, vectors and matrices
3  *
4  * Author: David Hoadley
5  *
6  * Description: (see vectors3d.h)
7  *
8  * Copyright (c) 2020, David Hoadley <vcrumble@westnet.com.au>
9  * All rights reserved.
10  *
11  * Redistribution and use in source and binary forms, with or without
12  * modification, are permitted provided that the following conditions are met:
13  *
14  * * Redistributions of source code must retain the above copyright notice, this
15  * list of conditions and the following disclaimer.
16  * * Redistributions in binary form must reproduce the above copyright notice,
17  * this list of conditions and the following disclaimer in the documentation
18  * and/or other materials provided with the distribution.
19  *
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
24  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  *
32  *==============================================================================
33  */
34 /*------------------------------------------------------------------------------
35  * Notes:
36  * Character set: UTF-8. (Non-ASCII characters appear in this file)
37  *----------------------------------------------------------------------------*/
38 
39 /* ANSI includes etc. */
40 #include "instead-of-math.h" /* for sincos() */
41 #include <math.h>
42 #include <stdlib.h> /* for NULL */
43 
44 /* Local and project includes */
45 #include "general.h"
46 #include "vectors3d.h"
47 
48 /*
49  * Local #defines and typedefs
50  */
51 DEFINE_THIS_FILE; /* For use by REQUIRE() - assertions */
52 
53 /*
54  *==============================================================================
55  *
56  * Implementation
57  *
58  *==============================================================================
59  *
60  * Global functions callable by other modules
61  *
62  *------------------------------------------------------------------------------
63  */
64 GLOBAL V3D_Vector *v3d_addV(V3D_Vector *destV,
65  const V3D_Vector *srcV1,
66  const V3D_Vector *srcV2)
67 /*! Add two vectors:
68  [destV] = [srcV1] + [srcV2]
69  \returns pointer to \a destV
70  \param[out] destV vector which will contain the result
71  \param[in] srcV1 first vector
72  \param[in] srcV2 second vector
73 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74 {
75  REQUIRE_NOT_NULL(srcV1);
76  REQUIRE_NOT_NULL(srcV2);
77  REQUIRE_NOT_NULL(destV);
78 
79  destV->a[0] = srcV1->a[0] + srcV2->a[0];
80  destV->a[1] = srcV1->a[1] + srcV2->a[1];
81  destV->a[2] = srcV1->a[2] + srcV2->a[2];
82  return destV;
83 }
84 
85 
86 
87 GLOBAL V3D_Vector *v3d_addToV(V3D_Vector *modV1, const V3D_Vector *srcV2)
88 /*! Modify first vector by adding a second one:
89  [modV1] += [srcV2]
90  \returns pointer to \a modV1
91  \param[in,out] modV1 first vector, to which the second will be added
92  \param[in] srcV2 second vector
93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94 {
95  REQUIRE_NOT_NULL(modV1);
96  REQUIRE_NOT_NULL(srcV2);
97 
98  modV1->a[0] += srcV2->a[0];
99  modV1->a[1] += srcV2->a[1];
100  modV1->a[2] += srcV2->a[2];
101  return modV1;
102 }
103 
104 
105 
106 GLOBAL V3D_Vector *v3d_subtractV(V3D_Vector *destV,
107  const V3D_Vector *srcV1,
108  const V3D_Vector *srcV2)
109 /*! Vector subtraction:
110  [destV] = [srcV1] - [srcV2]
111  \returns pointer to \a destV
112  \param[out] destV vector which will contain the result
113  \param[in] srcV1 first vector
114  \param[in] srcV2 second vector
115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116 {
117  REQUIRE_NOT_NULL(srcV1);
118  REQUIRE_NOT_NULL(srcV2);
119  REQUIRE_NOT_NULL(destV);
120 
121  destV->a[0] = srcV1->a[0] - srcV2->a[0];
122  destV->a[1] = srcV1->a[1] - srcV2->a[1];
123  destV->a[2] = srcV1->a[2] - srcV2->a[2];
124  return destV;
125 }
126 
127 
128 
129 GLOBAL V3D_Vector *v3d_subFromV(V3D_Vector *modV1, const V3D_Vector *srcV2)
130 /*! Modify first vector by subtracting a second one:
131  [modV1] -= [srcV2]
132  \returns pointer to \a modV1
133  \param[in,out] modV1 first vector, from which the second will be subtracted
134  \param[in] srcV2 second vector
135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 {
137  REQUIRE_NOT_NULL(modV1);
138  REQUIRE_NOT_NULL(srcV2);
139 
140  modV1->a[0] -= srcV2->a[0];
141  modV1->a[1] -= srcV2->a[1];
142  modV1->a[2] -= srcV2->a[2];
143  return modV1;
144 }
145 
146 
147 
148 GLOBAL double v3d_dotProductV(const V3D_Vector *srcV1, const V3D_Vector *srcV2)
149 /*! Return the dot product of the two vectors:
150  res = [srcV1] · [srcV2]
151  or res = ||srcV1|| * ||srcV2|| * cos(θ)
152  \returns scalar result of the dot product
153  \param[in] srcV1 first vector
154  \param[in] srcV2 second vector
155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156 {
157  REQUIRE_NOT_NULL(srcV1);
158  REQUIRE_NOT_NULL(srcV2);
159 
160  return srcV1->a[0] * srcV2->a[0]
161  + srcV1->a[1] * srcV2->a[1]
162  + srcV1->a[2] * srcV2->a[2];
163 }
164 
165 
166 
167 GLOBAL V3D_Vector *v3d_crossProductV(V3D_Vector *destV,
168  const V3D_Vector *srcV1,
169  const V3D_Vector *srcV2)
170 /*! Return the cross product of the two vectors:
171  [destV] = [srcV1] x [srcV2]
172  \returns pointer to vector result of the cross product (\a destV)
173  \param[out] destV vector which will contain the result
174  \param[in] srcV1 first vector
175  \param[in] srcV2 second vector
176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177 {
178  V3D_Vector temp; /* Allow caller to pass same vector for destV
179  * and srcV1 or srcV2 */
180  REQUIRE_NOT_NULL(srcV1);
181  REQUIRE_NOT_NULL(srcV2);
182  REQUIRE_NOT_NULL(destV);
183 
184  temp.a[0] = (srcV1->a[1] * srcV2->a[2]) - (srcV1->a[2] * srcV2->a[1]);
185  temp.a[1] = (srcV1->a[2] * srcV2->a[0]) - (srcV1->a[0] * srcV2->a[2]);
186  temp.a[2] = (srcV1->a[0] * srcV2->a[1]) - (srcV1->a[1] * srcV2->a[0]);
187  *destV = temp;
188  return destV;
189 }
190 
191 
192 
193 GLOBAL V3D_Vector *v3d_addToUV(V3D_Vector *modV1, const V3D_Vector *srcV2)
194 /*! Modify a unit length rectangular position vector \a modV1 by adding a
195  correction vector \a srcV2 to it and re-normalizing to unit length.
196  \returns Pointer to \a modV1
197  \param[in,out] modV1 First vector, to which \a srcV2 will be added, and
198  then re-normalized to unit length.
199  \param[in] srcV2 vector to be added to \a modV1
200 
201  \par When to call this function
202  If your vector \a modV1 is a unit vector, and you want it to be a unit
203  vector after the addition of \a srcV2, then you can use this function. But
204  if \a srcV2 is very small, you will be better off calling v3d_addToUVfast()
205  instead.
206 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207 {
208  double mag;
209 
210  REQUIRE_NOT_NULL(modV1);
211  REQUIRE_NOT_NULL(srcV2);
212 
213  v3d_addToV(modV1, srcV2);
214  mag = v3d_magV(modV1);
215  modV1->a[0] /= mag;
216  modV1->a[1] /= mag;
217  modV1->a[2] /= mag;
218  return modV1;
219 }
220 
221 
222 
223 GLOBAL V3D_Vector *v3d_addToUVfast(V3D_Vector *modV1, const V3D_Vector *srcV2)
224 /*! Modify a unit length rectangular position vector \a modV1 by adding a small
225  correction vector \a srcV2 to it and re-normalizing to unit length. The
226  scaling factor for the length is approximated by the expression
227  1 - [modV1] · [srcV2], where the dot indicates the Dot Product.
228  The correction vector \a srcV2 must be very small.
229  \returns Pointer to \a modV1
230  \param[in,out] modV1 First vector, to which \a srcV2 will be added, and
231  then re-normalized to unit length.
232  Valid range: |modV1| == 1
233  \param[in] srcV2 vector to be added to \a modV1.
234  Valid range: |srcV2| << 1
235 
236  \par When to call this function
237  This routine performs the same function as v3d_addToUV(), but without
238  calling \c sqrt(), or performing divisions. This can make it quicker than
239  v3d_addToUV(), depending on what kind of floating point processor you have,
240  if any. It is only valid if \a srcV2 is very small.
241 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
242 {
243  double scale;
244 
245  REQUIRE_NOT_NULL(modV1);
246  REQUIRE_NOT_NULL(srcV2);
247 
248  scale = 1.0 - v3d_dotProductV(modV1, srcV2);
249  v3d_addToV(modV1, srcV2);
250  modV1->a[0] *= scale;
251  modV1->a[1] *= scale;
252  modV1->a[2] *= scale;
253  return modV1;
254 }
255 
256 
257 
258 GLOBAL double v3d_magVSq(const V3D_Vector *srcV)
259 /*! Return the square of the magnitude of the specified vector
260 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
261 {
262  REQUIRE_NOT_NULL(srcV);
263 
264  return srcV->a[0] * srcV->a[0]
265  + srcV->a[1] * srcV->a[1]
266  + srcV->a[2] * srcV->a[2];
267 }
268 
269 
270 
271 GLOBAL double v3d_magV(const V3D_Vector *srcV)
272 /*! Return the magnitude of the specified vector
273 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274 {
275  return sqrt(v3d_magVSq(srcV));
276 }
277 
278 
279 
280 GLOBAL V3D_Vector *v3d_polarToRect(V3D_Vector *destV,
281  double alpha_rad,
282  double delta_rad)
283 /*! Converts polar (curvilinear) coordinates to equivalent rectangular
284  (Cartesian) coordinates. The vector returned is of unit length and in a
285  right or left-handed system according to the convention for measuring the
286  angle alpha_rad: for example left-handed for Hour Angle/Declination or
287  Azimuth/Elevation, and right-handed for Right Ascension/Declination.
288  \returns Pointer to \a destV, the resultant vector
289  \param[out] destV Resultant 3D unit vector in rectangular coordinates
290  \param[in] alpha_rad Secondary angle coordinate, e.g. RA or Azimuth (radians)
291  \param[in] delta_rad Primary angle coordinate, e.g. declination or elevation
292  (radians)
293  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
294 {
295  double sinA, cosA;
296  double sinD, cosD;
297 
298  REQUIRE_NOT_NULL(destV);
299 
300  sincos(alpha_rad, &sinA, &cosA);
301  sincos(delta_rad, &sinD, &cosD);
302 
303  destV->a[0] = cosD * cosA;
304  destV->a[1] = cosD * sinA;
305  destV->a[2] = sinD;
306  return destV;
307 }
308 
309 
310 
311 GLOBAL void v3d_rectToPolar(double *alpha_rad,
312  double *delta_rad,
313  const V3D_Vector *srcV)
314 /*! Converts rectangular (Cartesian) coordinates to the equivalent polar
315  (curvilinear) ones.
316  \param[out] alpha_rad Secondary angle coordinate (radians), range [-Pi, +Pi]
317  \param[out] delta_rad Primary angle coordinate (radians), range [-Pi/2, +Pi/2]
318  \param[in] srcV 3D unit vector in rectangular coordinates (direction
319  cosines)
320 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
321 {
322  REQUIRE_NOT_NULL(alpha_rad);
323  REQUIRE_NOT_NULL(delta_rad);
324  REQUIRE_NOT_NULL(srcV);
325 
326  *alpha_rad = atan2(srcV->a[1], srcV->a[0]);
327  *delta_rad = atan2(srcV->a[2],
328  sqrt(srcV->a[0] * srcV->a[0] + srcV->a[1] * srcV->a[1]));
329 }
330 
331 
332 
333 #if 0
334 void v3d_polar2Rect(double alpha_rad, double delta_rad, V3D_Vector *R)
335 /* Converts polar (curvilinear) coordinates to equivalent rectangular
336  (cartesian) coordinates. The vector returned is of unit length and in a
337  right or left-handed system according to the convention for measuring the
338  angle alpha_rad: for example left-handed for HA/Dec or Az/El and
339  right-handed for RA/Dec.
340  Inputs
341  alpha_rad - Secondary angle coordinate in radians
342  delta_rad - Primary angle coordinate in radians
343  Outputs
344  R - Resultant 3D unit vector in rectangular coordinates
345 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
346 {
347  R->a[0] = cos(delta_rad) * cos(alpha_rad);
348  R->a[1] = cos(delta_rad) * sin(alpha_rad);
349  R->a[2] = sin(delta_rad);
350 }
351 
352 
353 
354 void v3d_rect2Polar(const V3D_Vector *R,
355  bool alphaNonNegative,
356  double *alpha_rad,
357  double *delta_rad)
358 /* Converts rectangular (cartesian) coordinates to the equivalent polar
359  (curvilinear) ones.
360  Inputs
361  R - 3D unit vector in rectangular coordinates (direction cosines)
362  alphaNonNegative - if true, return alpha_rad in range [0, 2Pi)
363  Outputs
364  alpha_rad - Secondary angle result in radians, range (-Pi, +Pi]
365  (or in range [0, 2Pi) if alphaNonNegative == TRUE)
366  delta_rad - Primary angle result in radians, range [-Pi/2, +Pi/2]
367 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
368 {
369  *alpha_rad = atan2(R->a[1], R->a[0]);
370  if (alphaNonNegative && (*alpha_rad < 0.0)) { *alpha_rad += TWOPI; }
371  *delta_rad = atan2(R->a[2], sqrt(R->a[0] * R->a[0] + R->a[1] * R->a[1]));
372 }
373 #endif
374 
375 
376 
377 GLOBAL V3D_Matrix *v3d_createRotationMatrix(V3D_Matrix *destM,
378  V3D_AxisNames axis,
379  double angle_rad)
380 /*! Creates a matrix to rotate a coordinate system about an axis. This matrix
381  can then be used to convert a rectangular position vector from one
382  coordinate system to another
383  \returns pointer to \a destM, the resultant matrix
384  \param[out] destM 3x3 rotation matrix
385  \param[in] axis axis of rotation.
386  Xaxis (=0), Yaxis (=1), Zaxis (=2) No other values
387  allowed.
388  \param[in] angle_rad angle to rotate around the axis of rotation (radian).
389  When looking down the axis of rotation toward
390  the origin, positive angles of rotation are:
391  - clockwise for left-handed systems (like
392  Azimuth/Elevation)
393  - anticlockwise for right-handed systems (like
394  Right Ascension/Declination)
395 
396  \verbatim
397  ┌ 1 0 0 ┐
398  createRotationMatrix(M, Xaxis, θ) is R1(θ) = │ 0 cosθ sinθ│
399  └ 0 -sinθ cosθ┘
400 
401  ┌ cosθ 0 -sinθ┐
402  createRotationMatrix(M, Yaxis, θ) is R2(θ) = │ 0 1 0 │
403  └ sinθ 0 cosθ┘
404 
405  ┌ cosθ sinθ 0 ┐
406  createRotationMatrix(M, Zaxis, θ) is R3(θ) = │-sinθ cosθ 0 │
407  └ 0 0 1 ┘
408  \endverbatim
409  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
410 {
411  double cosA;
412  double sinA;
413  unsigned int i, j, k;
414 
415  REQUIRE((axis == Xaxis) || (axis == Yaxis) || (axis == Zaxis));
416  REQUIRE_NOT_NULL(destM);
417 
418  i = axis;
419  j = ((axis + 1) % 3);
420  k = ((axis + 2) % 3);
421 
422  destM->a[i][i] = 1.0;
423  destM->a[i][j] = 0.0;
424  destM->a[i][k] = 0.0;
425  destM->a[j][i] = 0.0;
426  destM->a[k][i] = 0.0;
427 
428  sincos(angle_rad, &sinA, &cosA);
429 
430  destM->a[j][j] = cosA;
431  destM->a[j][k] = sinA;
432  destM->a[k][j] = -sinA;
433  destM->a[k][k] = cosA;
434  return destM;
435 }
436 
437 
438 
439 GLOBAL V3D_Vector *v3d_multMxV(V3D_Vector *destV,
440  const V3D_Matrix *srcM,
441  const V3D_Vector *srcV)
442 /*! Multiply 3x3 matrix by 3x1 vector to give a new 3x1 vector, as per equation
443  [destV] = [srcM] * [srcV]
444  \returns pointer to \a destV
445  \param[out] destV vector which will contain the result
446  \param[in] srcM the matrix
447  \param[in] srcV vector which will be multiplied by the matrix
448 
449  \note
450  The same vector may be passed to \a destV and \a srcV, in which case the
451  contents of \a srcV will be replaced.
452 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
453 {
454  int i;
455  /* Allow the caller to pass same vector for destV and srcV. Also allows the
456  * compiler to optimise without us having to put "restrict" on the
457  * arguments. */
458  V3D_Vector temp;
459 
460  REQUIRE_NOT_NULL(srcV);
461  REQUIRE_NOT_NULL(srcM);
462  REQUIRE_NOT_NULL(destV);
463 
464  for (i = 0; i < 3; i++) {
465  temp.a[i] = srcM->a[i][0] * srcV->a[0]
466  + srcM->a[i][1] * srcV->a[1]
467  + srcM->a[i][2] * srcV->a[2];
468  }
469  *destV = temp;
470  return destV;
471 }
472 
473 
474 
475 GLOBAL V3D_Vector *v3d_multMtransxV(V3D_Vector *destV,
476  const V3D_Matrix *srcM,
477  const V3D_Vector *srcV)
478 /*! Multiply 3x1 vector by the transpose of the 3x3 matrix to give a new 3x1
479  vector, as per equation
480  [destV] = TRANSPOSE([srcM]) * [srcV]
481  \returns pointer to \a destV
482  \param[out] destV vector which will contain the result
483  \param[in] srcM the matrix
484  \param[in] srcV vector which will be multiplied by the matrix
485 
486  Note that for an orthogonal matrix which represents a rotation about
487  axis j by angle θ, i.e. Matrix = Rj(θ), the following identity holds:
488  TRANSPOSE(Matrix) = INVERSE(Matrix) = Rj(-θ).
489  This is why this is a useful routine
490 
491  \note
492  The same vector may be passed to \a destV and \a srcV, in which case the
493  contents of \a srcV will be replaced.
494 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
495 {
496  int i;
497  /* Allow the caller to pass same vector for destV and srcV. Also allows the
498  * compiler to optimise without us having to put "restrict" on the
499  * arguments. */
500  V3D_Vector temp;
501 
502  REQUIRE_NOT_NULL(srcV);
503  REQUIRE_NOT_NULL(srcM);
504  REQUIRE_NOT_NULL(destV);
505 
506  for (i = 0; i < 3; i++) {
507  temp.a[i] = srcM->a[0][i] * srcV->a[0]
508  + srcM->a[1][i] * srcV->a[1]
509  + srcM->a[2][i] * srcV->a[2];
510  }
511  *destV = temp;
512  return destV;
513 }
514 
515 
516 
517 GLOBAL V3D_Matrix *v3d_multMxM(V3D_Matrix *destM,
518  const V3D_Matrix *srcM1,
519  const V3D_Matrix *srcM2)
520 /*! Multiply 3x3 matrix by 3x3 matrix to give a new 3x3 matrix, as per
521  [destM] = [srcM1]*[srcM2]
522  \returns pointer to \a destM
523  \param[out] destM matrix which will contain the result
524  \param[in] srcM1 first matrix
525  \param[in] srcM2 second matrix
526 
527  \note
528  The same matrix may be passed to \a destM and either of the source matrices,
529  in which case the contents of that source matrix will be replaced.
530 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
531 {
532  int i, j;
533  /* Allow the caller to pass same matrix for destM and srcM1 or srcM2. Also
534  * allows the compiler to optimise without us having to put "restrict" on
535  * the arguments. */
536  V3D_Matrix temp;
537 
538  REQUIRE_NOT_NULL(srcM1);
539  REQUIRE_NOT_NULL(srcM2);
540  REQUIRE_NOT_NULL(destM);
541 
542  for (i = 0; i < 3; i++) {
543  for (j = 0; j < 3; j++) {
544  temp.a[i][j] = srcM1->a[i][0] * srcM2->a[0][j]
545  + srcM1->a[i][1] * srcM2->a[1][j]
546  + srcM1->a[i][2] * srcM2->a[2][j];
547  }
548  }
549  *destM = temp;
550  return destM;
551 }
552 
553 /*! \page page-why-struct-array Why define vector and matrix arrays as structs?
554  This has been done to assist compiler error checking, and reduce the chances
555  of bugs that can potentially be rather difficult to diagnose.
556 
557  Take for example the following code. It has a vector of rectangular
558  coordinates \c vec[], which contains the values one gets from converting
559  polar angles of 45° and 45° to rectangular form. This program simply
560  converts this rectangular vector back to polar form and writes out the two
561  angles, after converting from radian back to degrees.
562 
563  #include <math.h>
564  #include <stdio.h>
565 
566  void rectToPolar(double *alpha, double *delta, double rect[])
567  {
568  *alpha = atan2(rect[1], rect[0]);
569  *delta = atan2(rect[2], sqrt(rect[0] * rect[0] + rect[1] * rect[1]));
570  }
571 
572  int main(void)
573  {
574  double azimuth = 0.5; // Dummy initial value, will be replaced
575  double elevation = 0.5; // Ditto
576  double vec[3] = { 0.5, 0.5, 0.70710678119 };
577 
578  rectToPolar(vec, &azimuth, &elevation);
579  printf("Azimuth = %f, Elevation = %f\n",
580  azimuth * 180.0 / 3.141592653589793,
581  elevation * 180.0 / 3.141592653589793);
582  }
583 
584  There is a bug in the code, which in this short example is rather easy to
585  see, but in a large program might not be. The expected output from this
586  code is
587 
588  Azimuth = 45.000000, Elevation = 45.000000
589 
590  But what do you get if you compile and run it? On my running this program,
591  I got
592 
593  Azimuth = 48.002776, Elevation = 28.647890
594 
595  You may get something else. But in any case, the results may look
596  sufficiently plausible that you might not immediately notice that you have
597  got a garbage answer. Depending on your program, this might end up being
598  quite time-consuming to debug.
599 
600  Now instead, compile the following code, which uses our V3D_Vector struct
601  instead of using a simple array:
602 
603  #include <math.h>
604  #include <stdio.h>
605 
606  #include "vectors3d.h"
607 
608  int main(void)
609  {
610  double azimuth = 0.5; // Dummy initial value, will be replaced
611  double elevation = 0.5; // Ditto
612  V3D_Vector vec = {{ 0.5, 0.5, 0.70710678119 }};
613 
614  v3d_rectToPolar(&vec, &azimuth, &elevation);
615  printf("Azimuth = %f, Elevation = %f\n",
616  azimuth * 180.0 / 3.141592653589793,
617  elevation * 180.0 / 3.141592653589793);
618  }
619 
620  This time, you get two warning messages from the compiler. With the clang
621  compiler, the messages are as follows:
622 
623  vec2.c:12:21: warning: incompatible pointer types passing 'V3D_Vector *'
624  to parameter of type 'double *' [-Wincompatible-pointer-types]
625  v3d_rectToPolar(&vec, &azimuth, &elevation);
626  ^~~~
627  ./vectors3d.h:74:30: note: passing argument to parameter 'alpha_rad' here
628  void v3d_rectToPolar(double *alpha_rad,
629  ^
630  vec2.c:12:37: warning: incompatible pointer types passing 'double *' to
631  parameter of type 'const V3D_Vector *' [-Wincompatible-pointer-types]
632  v3d_rectToPolar(&vec, &azimuth, &elevation);
633  ^~~~~~~~~~
634  ./vectors3d.h:76:40: note: passing argument to parameter 'srcV' here
635  const V3D_Vector *srcV);
636 
637  (Actually, there is a case here for elevating this warning to an error that
638  stops compilation, using -Werror=incompatible-pointer-types)
639 
640  I believe that this justifies using this slightly more awkward way
641  of implementing the vectors and arrays needed for 3-D geometry.
642 
643  It also allows one to use a simple assignment statement to copy a vector or
644  matrix, as opposed to needing to write a function to do it. (Or in the case
645  of the IAU SOFA routines for matrix copying - two functions!)
646  */
v3d_magVSq
double v3d_magVSq(const V3D_Vector *srcV)
Return the square of the magnitude of the specified vector.
Definition: vectors3d.c:258
v3d_multMxM
V3D_Matrix * v3d_multMxM(V3D_Matrix *destM, const V3D_Matrix *srcM1, const V3D_Matrix *srcM2)
Multiply 3x3 matrix by 3x3 matrix to give a new 3x3 matrix, as per [destM] = [srcM1]*[srcM2].
Definition: vectors3d.c:517
general.h
general.h - definitions of general use to (standard) C programs
v3d_polarToRect
V3D_Vector * v3d_polarToRect(V3D_Vector *destV, double alpha_rad, double delta_rad)
Converts polar (curvilinear) coordinates to equivalent rectangular (Cartesian) coordinates.
Definition: vectors3d.c:280
REQUIRE
#define REQUIRE(test_)
Check preconditions.
Definition: general.h:273
v3d_magV
double v3d_magV(const V3D_Vector *srcV)
Return the magnitude of the specified vector.
Definition: vectors3d.c:271
V3D_AxisNames
V3D_AxisNames
Enumeration to be used by the v3d_createRotationMatrix() function.
Definition: vectors3d.h:62
v3d_rectToPolar
void v3d_rectToPolar(double *alpha_rad, double *delta_rad, const V3D_Vector *srcV)
Converts rectangular (Cartesian) coordinates to the equivalent polar (curvilinear) ones.
Definition: vectors3d.c:311
GLOBAL
#define GLOBAL
See above.
Definition: general.h:134
REQUIRE_NOT_NULL
#define REQUIRE_NOT_NULL(pointer_)
,
Definition: general.h:308
v3d_addToUV
V3D_Vector * v3d_addToUV(V3D_Vector *modV1, const V3D_Vector *srcV2)
Modify a unit length rectangular position vector modV1 by adding a correction vector srcV2 to it and ...
Definition: vectors3d.c:193
v3d_multMtransxV
V3D_Vector * v3d_multMtransxV(V3D_Vector *destV, const V3D_Matrix *srcM, const V3D_Vector *srcV)
Multiply 3x1 vector by the transpose of the 3x3 matrix to give a new 3x1 vector, as per equation [des...
Definition: vectors3d.c:475
TWOPI
#define TWOPI
,
Definition: general.h:160
v3d_crossProductV
V3D_Vector * v3d_crossProductV(V3D_Vector *destV, const V3D_Vector *srcV1, const V3D_Vector *srcV2)
Return the cross product of the two vectors: [destV] = [srcV1] x [srcV2].
Definition: vectors3d.c:167
v3d_subFromV
V3D_Vector * v3d_subFromV(V3D_Vector *modV1, const V3D_Vector *srcV2)
Modify first vector by subtracting a second one: [modV1] -= [srcV2].
Definition: vectors3d.c:129
V3D_Matrix
3x3 matrix.
Definition: vectors3d.h:51
v3d_addToV
V3D_Vector * v3d_addToV(V3D_Vector *modV1, const V3D_Vector *srcV2)
Modify first vector by adding a second one: [modV1] += [srcV2].
Definition: vectors3d.c:87
instead-of-math.h
instead-of-math.h - header to be included instead of math.h
v3d_subtractV
V3D_Vector * v3d_subtractV(V3D_Vector *destV, const V3D_Vector *srcV1, const V3D_Vector *srcV2)
Vector subtraction: [destV] = [srcV1] - [srcV2].
Definition: vectors3d.c:106
V3D_Vector
3x1 vector.
Definition: vectors3d.h:57
v3d_dotProductV
double v3d_dotProductV(const V3D_Vector *srcV1, const V3D_Vector *srcV2)
Return the dot product of the two vectors: res = [srcV1] · [srcV2] or res = ||srcV1|| * ||srcV2|| * c...
Definition: vectors3d.c:148
DEFINE_THIS_FILE
#define DEFINE_THIS_FILE
Leave DEFINE_THIS_FILE undefined.
Definition: general.h:252
v3d_addV
V3D_Vector * v3d_addV(V3D_Vector *destV, const V3D_Vector *srcV1, const V3D_Vector *srcV2)
Add two vectors: [destV] = [srcV1] + [srcV2].
Definition: vectors3d.c:64
v3d_createRotationMatrix
V3D_Matrix * v3d_createRotationMatrix(V3D_Matrix *destM, V3D_AxisNames axis, double angle_rad)
Creates a matrix to rotate a coordinate system about an axis.
Definition: vectors3d.c:377
vectors3d.h
vectors3d.h - Three dimensional geometry, vectors and matrices
v3d_addToUVfast
V3D_Vector * v3d_addToUVfast(V3D_Vector *modV1, const V3D_Vector *srcV2)
Modify a unit length rectangular position vector modV1 by adding a small correction vector srcV2 to i...
Definition: vectors3d.c:223
sincos
static void sincos(double angle_rad, double *sinA, double *cosA)
Calculates sine and cosine of an angle.
Definition: instead-of-math.h:114
v3d_multMxV
V3D_Vector * v3d_multMxV(V3D_Vector *destV, const V3D_Matrix *srcM, const V3D_Vector *srcV)
Multiply 3x3 matrix by 3x1 vector to give a new 3x1 vector, as per equation [destV] = [srcM] * [srcV]...
Definition: vectors3d.c:439