angle-component.lisp 15.7 KB
Newer Older
1
2
;; Angles by component parts.
;; Liam Healy Thu May 30 2002 - 11:46
3
;; Time-stamp: <2014-10-06 22:29:48EDT angle-component.lisp>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.

(in-package :antik)

(export '(angle asine acosine atangent arise arun sine consine
	  rise run normalize-angle))
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

;;; Angles may be defined by "components", that is, trigonometric functions of
;;; the angle, halfplanes, etc.  Components may be defined without the whole
;;; angle being known, with the rest of the information to be supplied later. 
;;; At the point enough information is supplied to define the angle,
;;; the value may be extracted with pq-magnitude, which will happen
;;; automatically if it is printed.

;;; Examples:
;;; 1) A normal definition with a magnitude
;;; (setf normal-def (make-instance 'angle :magnitude 1.234))

;;; 2) Define "rise" and "run" (fixed multiples of sine and cosine);
;;;   note that a symbol label is required to confirm the same multiplier:
;;; (setf (rise true-anomaly 'eccentricity) 0.02)
;;; (setf (run true-anomaly 'eccentricity) 0.04)
;;; true-anomaly is now #_0.46365_RADIANS.

;;; 3) Define an angle with the cosine and the vertical halfplane:
;;; (setf (cosine ang) 0.34)
;;; (setf (vertical-halfplane ang) -1)
;;; Is #_-1.2239_RADIANS.

;;;;****************************************************************************
;;;; Angle object
;;;;****************************************************************************
healy's avatar
healy committed
51

Liam M. Healy's avatar
Liam M. Healy committed
52
53
#|
;;; EXPERIMENTAL
54
(defclass angle (physical-quantity)
55
56
  ((magnitude :initform nil)
   (dimension :initform (dimension 'angle))
57
   (sine :initform nil :initarg :sine :reader sine
58
	 :writer sinecomp :type (or nil number))
59
   (cosine :initform nil :initarg :cosine :reader cosine
60
	   :writer cosinecomp :type (or nil number))
61
62
63
64
   (tangent :initform nil :initarg :tangent
	    :accessor tangent :type (or nil list number))
   (vertical-halfplane :initform nil :initarg :vertical-halfplane
		       :accessor vertical-halfplane
65
		       :documentation "Sign of sine"
66
67
		       :type (or nil (integer -1 1)))
   (horizontal-halfplane :initform nil
68
			 :initarg :horizontal-halfplane
69
			 :accessor horizontal-halfplane
70
			 :documentation "Sign of cosine"
71
			 :type (or nil (integer -1 1)))
72
   (normalized :initform nil :initarg :normalized :accessor normalized
73
74
	       :documentation "Fractional part of whole revolution."
	       ;; may be -pi to pi, 0 to 2pi, etc.
75
	       :type t)
76
77
   (whole-revolutions :initform 0 :initarg :whole-revolutions
		      :accessor whole-revolutions
78
79
80
		      :documentation
		      "Number of complete revolutions, which when added to normalized,
                       gives the angle."
81
		      :type integer))
healy's avatar
healy committed
82
83
  (:documentation "Angles with component parts."))

84
85
86
87
;;;;****************************************************************************
;;;; Angle consistency, definition, component extraction, printing
;;;;****************************************************************************

88
89
90
91
92
93
94
95
(defmacro with-angle (angle &body body)
  `(with-slots
       (magnitude dimension sine cosine tangent
	vertical-halfplane horizontal-halfplane
	normalized whole-revolutions)
       ,angle
     ,@body))

96
97
98
99
100
101
(defun check-angle-consistency (angle)
  "Check that the components of the angle are consistent if defined,
   and return the angle."
  (check-type angle angle)
  (with-angle angle
    (when (and sine cosine)
102
      (assert (antik:zerop (-nz (cl:+ (cl:expt sine 2) (cl:expt cosine 2)) 1.0))
103
104
	  (sine cosine)
	"Sine and cosine do not both belong to the same angle."))
105
    (when (and sine cosine tangent (physical-quantity-p tangent))
106
      (assert (antik:zerop (-nz (cl:/ sine cosine) tangent))
107
108
109
110
111
112
113
114
115
116
117
118
	  (sine cosine tangent)
	"Tangent is not the ratio of sine and cosine."))
    (when (and sine vertical-halfplane)
      (assert (eql (signum sine) vertical-halfplane)
	  (sine vertical-halfplane)
	"Sine and vertical halfplane do not agree."))
    (when (and cosine horizontal-halfplane)
      (assert (eql (signum cosine) horizontal-halfplane)
	  (cosine horizontal-halfplane)
	"Cosine and horizontal halfplane do not agree.")))
  angle)

119
120
121
(defun make-angle (&rest args)
  (check-angle-consistency (apply #'make-instance 'angle args)))

122
(defmethod print-object ((angle angle) stream)
123
  (if (magnitude-or-nil angle)
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
      (call-next-method)
    (print-unreadable-object (angle stream :identity t)
      (with-angle angle
	(format
	 stream
	 ;;"Partially known angle:~#[ none~; ~S~; ~S and ~S~:;~@{~#[~; and~] ~S~^,~}~]"
	 "Partially known angle: ~:{~(~s~)=~d~:^, ~}"
	 (remove nil
		 `((sine ,sine) (cosine ,cosine) (tangent ,tangent)
				(vertical-halfplane ,vertical-halfplane)
				(horizontal-halfplane ,horizontal-halfplane)
				(normalized ,normalized)
				(whole-revolutions ,whole-revolutions))
		 :key #'second))))))

;;;;****************************************************************************
;;;; Trigonometric functions
;;;;****************************************************************************

143
144
#+(or) ; uses pythagorean-complement
(defmethod cos ((angle angle))
145
146
147
148
  (with-angle angle
    (or cosine
	(setf cosine
	  (cond
149
	   ((magnitude-or-nil angle) (antik:cos (pq-magnitude angle)))
150
	   ((and sine horizontal-halfplane)
151
	    (antik:* horizontal-halfplane (pythagorean-complement sine)))
152
	   ((and tangent vertical-halfplane)
153
154
	    (antik:* vertical-halfplane
		(/ (sqrt (1+ (expt (antik:tan angle) 2))))))))
155
	(error "Insufficient information to compute cosine."))))
156

157
158
#+(or) ; uses pythagorean-complement
(defmethod sin ((angle angle))
159
160
161
162
  (with-angle angle
    (or sine
	(setf sine
	  (cond
163
	   ((magnitude-or-nil angle) (antik:sin (pq-magnitude angle)))
164
	   ((and cosine vertical-halfplane)
165
	    (antik:* vertical-halfplane (pythagorean-complement cosine)))
166
	   ((and tangent vertical-halfplane)
167
168
	    (antik:* horizontal-halfplane
		(/ (antik:tan angle) (sqrt (1+ (expt (antik:tan angle) 2))))))))
169
	(error "Insufficient information to compute sine."))))
170

171
(defmethod tan ((angle angle))
172
  (with-angle angle
173
    (if (and tangent (physical-quantity-p tangent))
174
175
176
	tangent
      (setf tangent
	(cond
177
	 ((magnitude-or-nil angle) (antik:tan (pq-magnitude angle)))
178
	 (tangent
179
	  (cl:/ (first tangent) (second tangent)))
180
	 ((or sine cosine)
181
	  (cl:/ (antik:sin angle) (antik:cos angle)))
182
183
	 (t (error "Insufficient information to compute tangent.")))))))

184
;;;;****************************************************************************
185
;;;; Defining an angle by setting components
186
187
188
;;;;****************************************************************************

(defun acosine (cosine &optional angle)
189
190
  "Start defining an angle by specifying its cosine.
   Returns the angle object."
191
192
  (check-angle-consistency
   (if angle
193
       (progn (cosinecomp angle cosine) angle)
194
     (make-angle :cosine cosine))))
195

196
(defun asine (sine &optional angle)
197
198
199
200
  "Start defining an angle by specifying its sine.
   Returns the angle object."
  (check-angle-consistency
   (if angle
201
       (progn (sinecomp angle sine) angle)
202
     (make-angle :sine sine))))
203
204

(defun atangent (tangent &optional angle)
205
  "Define an angle by specifying the tangent.  Returns the angle object."
206
  (if angle
207
208
      (progn (setf (tangent angle) tangent)
	     angle)
209
210
    (make-angle :tangent tangent)))

211
(defun arise (rise &optional label angle)
212
  "Make or modify an angle by components, label is a symbol
213
   naming the multiplier of sine.  Returns the angle object."
214
215
216
217
218
219
  (if (and angle (tangent angle)
	   ;; If the run is already defined, 
	   (listp (tangent angle))
	   ;; check label
	   (eql (first (tangent angle)) label))
      ;; and replace
220
221
222
223
      (progn
	(setf (first (tangent angle)) rise)
	angle)
    ;; otherwise, make a new angle object
224
225
    (atangent (list rise label) nil)))

226
(defun arun (run &optional label angle)
227
  "Make or modify an angle by components, label is a symbol
228
   naming the multiplier of cosine.  Returns the angle object." 
229
230
231
232
233
234
  (if (and angle (tangent angle)
	   ;; If the run is already defined, 
	   (listp (tangent angle))
	   ;; check label
	   (eql (second (tangent angle)) label))
      ;; and replace
235
236
237
238
      (progn
	(setf (second (tangent angle)) run)
	angle)
    ;; otherwise, make a new angle object
239
240
    (atangent (list label run) nil)))

241
242
243
244
245
;;;;****************************************************************************
;;;; Setf expansions to set an angle's components
;;;;****************************************************************************

(defmacro place-unbound-p (place)
healy's avatar
healy committed
246
  "Is the place unbound?"
247
248
249
250
251
252
253
254
255
256
  `(typep (nth-value 1 (ignore-errors ,place)) 'cell-error))

(defun angle-component-writer-form (form angle sv newval setter getter)
  "The writer-form in the setf-expansion of the angle."
  `(progn
     ;; If ,getter is not bound, set to nil so that it may be
     ;; reset to the angle object.
     ;; boundp does not work on non-symbols (general places)
     ;; boundp will not detect lexical bindings.
     (when (place-unbound-p ,getter)
healy's avatar
healy committed
257
       (eval (list 'let (list (list (first ',newval) nil)) ',setter)))
258
259
260
     (setf ,angle ,form)
     ,sv))

261
(defmacro set-angle-component (angle setform env)
262
263
264
265
  "Define an expansion to setf an angle component.  The setform should
   have an argument ,sv which is the value to set."
  `(multiple-value-bind (dummies vals newval setter getter)
       (get-setf-expansion angle ,env)
266
     (alexandria:with-gensyms (sv)
267
       (values
268
269
270
271
272
273
274
	 dummies
	 vals
	 (list sv)
	 (angle-component-writer-form
	  ,setform
	  ,angle sv newval setter getter)
	 nil))))
275

276
(define-setf-expander rise (angle label &environment env)
277
278
279
280
  "Set the rise of an angle, e.g.
    (setf (rise true-anomaly 'eccentricity) 0.02).
   The label is a symbol naming the multiplier of sine.
   It is used to determine that the same multiplier was used for the run."
281
  (set-angle-component angle `(arise ,sv ,label ,angle) env))
282

283
(define-setf-expander run (angle label &environment env)
284
  "Set the run of an angle, e.g.
285
    (setf (run true-anomaly 'eccentricity) 0.02).
286
287
   The label is a symbol naming the multiplier of cosine.
   It is used to determine that the same multiplier was used for the rise."
288
  (set-angle-component angle `(arun ,sv ,label ,angle) env))
289

290
291
292
293
;;; These setf expanders are better than the :accessor because it is
;;; not necessary to define an angle object first.  2004-06-17 These
;;; are now the only ones; changed slot definition to be :reader
;;; because sbcl objected.
294
(define-setf-expander cosine (angle &environment env)
295
  "Set the cosine of an angle."
296
  (set-angle-component angle `(acosine ,sv ,angle) env))
297

298
(define-setf-expander sine (angle &environment env)
299
  "Set the sine of an angle."
300
  (set-angle-component angle `(asine ,sv ,angle) env))
healy's avatar
healy committed
301

302
;;;;****************************************************************************
303
;;;; Angle value
304
305
306
307
308
309
310
311
312
313
314
;;;;****************************************************************************

;;; Conditional relations for finding the angle
;;; angle = whole-revolutions + normalized
;;; normalized = (atan tangent)
;;; tangent = (list sine cosine)
;;; normalized = (* horizontal-halfplane (acos cosine))
;;; normalized =
;;;  (let ((asin (asin sine)))
;;;     (if (minusp vertical-halfplane) (- pi asin) (if (minusp asin) (+ 2pi asin) asin)))

315
(defparameter *half-rev* (make-pq 1/2 'revolution))
healy's avatar
healy committed
316
317

(defmethod pq-magnitude :before ((angle angle))
318
  ;; I'm undecided on whether to use antik:atan etc. here vs. atan;
319
320
  ;; I'd like to work strictly in radians, but it might be necessary
  ;; to process generic numbers like bigfloats.
healy's avatar
healy committed
321
  (with-angle angle
322
    ;; Try to solve for the angle when component information is present
healy's avatar
healy committed
323
    (when (null magnitude)
324
325
      (unless normalized
	(setf normalized
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
	      (cond ((or tangent (and sine cosine))
		     (unless tangent (setf tangent (list sine cosine)))
		     (if (and (listp tangent) (not (some #'symbolp tangent)))
			 ;; two numerical values for rise/run in tangent slot; compute
			 (apply #'antik:atan tangent)
			 (cond (vertical-halfplane
				(antik:+ (if (minusp vertical-halfplane)
					     *half-rev* 0)
					 (antik:atan tangent)))
			       (horizontal-halfplane
				(let ((atan (antik:atan tangent)))
				  (antik:+
				   (if (cl:minusp (cl:* atan horizontal-halfplane))
				       *half-rev* 0)
				   atan)))
			       (t nil))))
		    ((and cosine vertical-halfplane)
		     (antik:* vertical-halfplane (antik:acos cosine)))
		    ((and sine horizontal-halfplane)
		     (let ((asin (antik:asin sine)))
		       (if (minusp horizontal-halfplane)
			   (antik:- *half-rev* asin) ; result between pi/2, 3pi/2
			   asin)))	; result between -pi/2, pi/2
		    (t (error "Insufficient information to find angle.")))))
350
351
352
      (setf
       magnitude
       (pqval (antik:+ normalized (make-pq whole-revolutions 'revolution)))))))
353
354
355

(defun magnitude-or-nil (angle)
  (ignore-errors (pq-magnitude angle)))
Liam M. Healy's avatar
Liam M. Healy committed
356
|#
357

358
359
360
361
;;;;****************************************************************************
;;;; Angle normalization
;;;;****************************************************************************

362
363
364
(defgeneric normalize-angle (angle &optional positive)
  (:documentation "Normalizes the angle to (in radians) -pi to pi
   or 0 to 2pi (positive = t).
365
366
   Two values are returned: the normalized angle, and the number
   of full revolutions that should be added to the normalized angle
367
368
369
   to recreate the original angle.")
  (:method ((angle number) &optional positive)
	   (multiple-value-bind (revs norm)
370
	       (funcall (if positive #'antik:floor #'antik:round) 
371
372
373
			angle 2pi)
	     (values norm revs))))

Liam M. Healy's avatar
Liam M. Healy committed
374
#|
375
(defmethod normalize-angle ((angle physical-quantity) &optional positive)
376
377
378
379
  (if (check-dimension angle 'angle nil)
      (let ((*zero-is-dimensionless* nil))
	(with-pq ((angle angle))
	  (unless (typep angle 'angle)
380
	    ;; Object may be a pdq with angle dimension, but not yet
381
382
383
384
	    ;; an angle object, so change it.
	    (change-class angle 'angle))
	  (with-angle angle
	    (multiple-value-bind (norm revs)
385
		(normalize-angle (pq-magnitude angle) positive)
386
387
	      (setf whole-revolutions revs
		    normalized (make-pq norm 'radian))
Liam M. Healy's avatar
Liam M. Healy committed
388
	      (values norm revs)))))
389
    angle))
Liam M. Healy's avatar
Liam M. Healy committed
390
391
392
393
394
395
396
|#

(defmethod normalize-angle ((angle physical-quantity) &optional positive)
  (if (check-dimension angle 'angle nil)
      (let ((*zero-is-dimensionless* nil))
	(with-pq ((angle angle))
	  (unless (typep angle 'angle)
397
	    ;; Object may be a pdq with angle dimension, but not yet
Liam M. Healy's avatar
Liam M. Healy committed
398
399
400
	    ;; an angle object, so change it.
	    (change-class angle 'angle))
	  (multiple-value-bind (norm revs)
401
	      (normalize-angle (pq-magnitude angle) positive)
402
	    (values (make-pq norm 'radian) revs))))
Liam M. Healy's avatar
Liam M. Healy committed
403
      angle))
404

405
406
407
408
409
410
411
412
(defmethod normalize-angle ((seq sequence) &optional positive)
  (map (type-of seq)
    (lambda (angle)
      (if (check-dimension angle 'angle nil)
	  (normalize-angle angle positive)
	angle))
    seq))

Liam M. Healy's avatar
Liam M. Healy committed
413
#|
414
415
416
417
418
419
420
;;; Formerely norm-denorm-angle
(defun apply-function-to-normalized-angle (angle function &optional positive)
  "Normalize the angle, apply the function, and denormalize it back to the
   original cycle count.  This is useful for e.g. Kepler's equation,
   where the normalized angle is needed to for the algorithm to
   converge properly, but the full cycles are needed for f & g functions."
  (normalize-angle angle positive)
421
422
  (antik:+ (funcall function (normalized angle))
	   (make-pq (whole-revolutions angle) 'revolution)))
Liam M. Healy's avatar
Liam M. Healy committed
423
|#