;;;Nikolas von Cusa 's angle trisection approximation
;;;This is based on the following question.
;let f(t) = (1+n)sin(t) / (n + cos(t))
;for what integer value of  n , f(t) will be closest to t ?
;
;;The answer is n = 2
;;
;;;voncusa.lsp
;;;
;;;May 15, 2006		Takaya Iwamoto
(prompt "#1")
;;;
(defun c:cusa_base( )

	(setup_cusa_base)	
	(setq n_max 8
		theta_start (dtr 0.)
		theta_end (dtr 90.)
		theta_step (dtr 5.)
		num_theta (fix (/ theta_end theta_step))
		n_cur 1 
	)
	
	(repeat n_max
		(setq n (float n_cur) nt_cur 0 n_color (itoa n_cur))
		(repeat (1+ num_theta)
			(setq t_cur (+ theta_start (* nt_cur theta_step))
				num (* (+ 1 n) (sin t_cur) )
				den (+ n (cos t_cur))
				ft (/ num den)
				dif (- ft t_cur)
				new_pt (list t_cur dif)
			)
			;(make_pt "0" n_cur new_pt)
			(if (= nt_cur 0)  (make_pt "0" n_cur new_pt)
							
					(make_line_1 "0" n_cur old_pt new_pt)
			)		
					
			(setq old_pt new_pt)
			(setq nt_cur (1+ nt_cur))
		);;;end of inner repeat loop
		(setq n_cur (1+ n_cur))
	);;;end of outer repeat loop
	(command "_.zoom" "_Extent")
	(command "_.regen")
	(reset_sysvar)
);;;CUSA_BASE
;;;
(prompt "#2")
;;;
(defun setup_cusa_base()
	(setup_sysvar)
	(setq half_pi (* 0.5 pi) x_right (list half_pi 0)  y_up '(0 0.5) y_down '(0 -0.5)
		 pnt_org '(0 0) pnt_ref (list half_pi half_pi) )

	(make_line_1 "0" 8 y_up y_down)
	(make_line_1 "0" 8 x_right pnt_org)
	(command "_.zoom" "_Extent")
	(command "_.regen")
);;;setup
;;;
(prompt "#2a")
;;;
(defun setup_cusa_base_2()
	(setup_sysvar)
	(setq half_pi (* 0.5 pi) x_right (list half_pi 0)  y_up (list 0 half_pi)
		 pnt_org '(0 0) pnt_ref (list half_pi half_pi) del_x (/ half_pi 9.) del_y del_x )
	(setq 	theta_start (dtr 0.)
		theta_end (dtr 90.)
		theta_step (dtr 10.)
		num_theta (fix (/ theta_end theta_step))
		n_cur 2 n_t 0 chr_size 0.05 
	)

	(make_line_1 "0" 8 y_up pnt_org)
		(setq vertical_line (entlast))
	(make_line_1 "0" 8 x_right pnt_org)
		(setq horizontal_line (entlast))
	;;;draw grid lines with del_x=del_y = half_pi/9 (10 degrees)
	(command "_.array" horizontal_line "" "R" 10 1 del_y)
	(command "_.array" vertical_line "" "R" 1 10 del_x)
	;;;add x & y axis values
	(repeat 10
		(setq 	x_t (* n_t 10.) 		y_t (* n_t 10.)
			x_axis (list (dtr x_t) -0.10) 	y_axis (list -0.2 (dtr y_t))
			x_text (itoa (fix x_t)) 	y_text (rtos (dtr y_t) 2 2)
		)
		(textdisplay x_text x_axis chr_size 0.)
		(textdisplay y_text y_axis chr_size 0.)
		(setq n_t (1+ n_t))
	);;end of repeat loop
		
	(make_line_1 "0" 8 pnt_ref pnt_org)
	(command "_.zoom" "_Extent")
	(command "_.regen")
);;;setup
(prompt "#2b")
;;;
(defun setup_cusa_base_3()
	(setup_sysvar)
	(setq half_pi (* 0.5 pi) x_right (list half_pi 0)  y_up (list 0 1.6)
		 pnt_org '(0 0) pnt_ref (list half_pi half_pi) del_x (/ half_pi 9.) del_y 0.4)
	(setq 	theta_start (dtr 0.)
		theta_end (dtr 90.)
		theta_step (dtr 10.)
		num_theta (fix (/ theta_end theta_step))
		n_cur 2 n_t 0 chr_size 0.07 
	)

	(make_line_1 "0" 8 y_up pnt_org)
		(setq vertical_line (entlast))
	(make_line_1 "0" 8 x_right pnt_org)
		(setq horizontal_line (entlast))
	;;;draw grid lines with del_x=del_y = half_pi/9 (10 degrees)
	(command "_.array" horizontal_line "" "R" 5 1 del_y)
	(command "_.array" vertical_line "" "R" 1 10 del_x)
	;;;add x & y axis values
	(repeat 10
		(setq 	x_t (* n_t 10.) 		
			x_axis (list (dtr x_t) -0.20) 	
			x_text (itoa (fix x_t)) 	
		)
		(textdisplay x_text x_axis chr_size 0.)
		;(textdisplay y_text y_axis chr_size 0.)
		(setq n_t (1+ n_t))
	);;end of repeat loop
	(setq n_t 0)
	(repeat 5
		(setq 		y_t  n_t
			 	y_axis (list -0.07 (* 0.4 y_t))
			 	y_text (itoa (fix y_t))
		)
		(textdisplay y_text y_axis chr_size 0.)
		(setq n_t (1+ n_t))
	);;end of repeat loop
		
	;(make_line_1 "0" 8 pnt_ref pnt_org)
	(command "_.zoom" "_Extent")
	(command "_.regen")
);;;setup
(prompt "#3")
;;;
(defun c:cusa_base_2( )

	(setup_cusa_base_2)	
	
	
	
		(setq n (float n_cur) nt_cur 0 n_color (itoa n_cur))
		(repeat (1+ num_theta)
			(setq t_cur (+ theta_start (* nt_cur theta_step))
				num (* (+ 1 n) (sin t_cur) )
				den (+ n (cos t_cur))
				ft (/ num den)
				;dif (- ft t_cur)
				new_pt (list t_cur ft)
			)
			;(make_pt "0" n_cur new_pt)
			(if (= nt_cur 0)  (make_pt "0" n_cur new_pt)
							
					(make_line_1 "0" n_cur old_pt new_pt)
			)		
					
			(setq old_pt new_pt)
			(setq nt_cur (1+ nt_cur))
		);;;end of inner repeat loop
		
	(command "_.zoom" "_Extent")
	(command "_.regen")
	(reset_sysvar)
);;;CUSA_BASE_2
;;;
(prompt "#4")
;;;
(defun c:cusa_error( )

	(setup_cusa_base_3)	
	
	
	
		(setq n (float n_cur) nt_cur 0 n_color (itoa n_cur))
		(repeat (1+ num_theta)
			(setq t_cur (+ theta_start (* nt_cur theta_step))
				num (sin t_cur) 
				den (+ n (cos t_cur))
				ft (/ num den)
				trisect (atan ft)
				 dif (* 0.4 (rtd (- (/ t_cur 3.) trisect)))
				new_pt (list t_cur dif)
			)
			(if (= nt_cur 0)  (make_pt "0" n_cur new_pt)
							
					(make_line_1 "0" n_cur old_pt new_pt)
			)		
			(setq old_pt new_pt)
			(setq nt_cur (1+ nt_cur))
		);;;end of inner repeat loop
		
	(command "_.zoom" "_Extent")
	(command "_.regen")
	(reset_sysvar)
);;;CUSA_ERROR
;;;
;;;Von Cusa - Snellius main
;;;
(defun C:Voncusa()
	(setup_voncusa)
	;connect points P and A, and extend it to meet BG at D.
	(setq pnt_d (inters pnt_p pnt_a pnt_b pnt_g nil))  (make_pt "0" 0 pnt_d)
	(make_line_1 "0" 0 pnt_p pnt_d)
	(mark_id pnt_d "D" 1 0.1)
	;find M on BD such that BM = BD/3
	(setq pnt_m (plt pnt_b pnt_d (/ 1. 3.))) (make_pt "0" 0 pnt_m)
	(mark_id pnt_m "M" 1 0.1)
	;
	;compute intersection of line PM and the unit circle
	(make_line_1 "0" 2 pnt_p pnt_m)
	(setq pnt_t (compute_intersection))
	(make_pt "0" 2 pnt_t) (mark_id pnt_t "T" 3 0.1)
	;;connect pnt_t to the origin This line trisects angle AOB
	(make_line_1 "0" 1  pnt_org  pnt_t)

	(command "_.zoom" "_Extent")
	(command "_.regen")
	(reset_sysvar)
);;;VONCUSA
;;;
;;;VONCUSA_MOD	modified Von Cusa by Takaya Iwamoto	May 22, 2006
;;;Instead of drawing PM, draw a line parallel to OB through M. Intersection of this line and the arc
;;;is the point T. For given angles below 90, this method gives much better accuracy.
;;;
(defun C:Voncusa_mod()
	(setup_voncusa)
	;connect points P and A, and extend it to meet BG at D.
	(setq pnt_d (inters pnt_p pnt_a pnt_b pnt_g nil))  (make_pt "0" 0 pnt_d)
	(make_line_1 "0" 0 pnt_p pnt_d)
	(mark_id pnt_d "D" 1 0.1)
	;find M on BD such that BM = BD/3
	(setq pnt_m (plt pnt_b pnt_d (/ 1. 3.))) (make_pt "0" 0 pnt_m)
	(mark_id pnt_m "M" 1 0.1)
	;
	;compute intersection of line PM and the unit circle
	;(make_line_1 "0" 2 pnt_p pnt_m)
	;(setq pnt_t (compute_intersection))
	(setq bm (cadr pnt_m) bm_sqr (* bm bm) pnt_t (list (sqrt (- 1 bm_sqr)) bm))
	(make_pt "0" 2 pnt_t) (mark_id pnt_t "T" 3 0.1)
	;;connect pnt_t to the origin This line trisects angle AOB
	(make_line_1 "0" 1  pnt_org  pnt_t)
	;;;compute angle TOB
	(setq rad_tob (angle pnt_org pnt_t) deg_tob (rtd rad_tob)
		dif_deg (- ( / deg_aob 3.) deg_tob) dif_sec (* 3600 dif_deg)
		dif_min (* 60. dif_deg)
	)
	
	(command "_.zoom" "_Extent")
	(command "_.regen")
	(reset_sysvar)
);;;VONCUSA
;;;
;;;
(defun setup_voncusa()
	(setup_sysvar)
	(setq pnt_org '(0 0) pnt_b '(1 0) pnt_e '(-1 0) pnt_p '(-2 0)
		pnt_g '( 1 1.75)
	)
	(setvar "PDMODE" 32)
	(setvar "PDSIZE" -5)
		;;;create points $ point id's
	(make_pt "0" 0 pnt_b) (make_pt "0" 0 pnt_org) (make_pt "0" 0 pnt_e)
	(make_pt "0" 0 pnt_p) (make_pt "0" 0 pnt_g)
	(mark_id pnt_b "B" 4 0.1) (mark_id pnt_org "O" 4 0.1) (mark_id pnt_e "E" 4 0.1)
	(mark_id pnt_p "P" 4 0.1) (mark_id pnt_g "G" 1 0.1)
		;;;draw ref lines
	(make_line_1 "0" 0 pnt_p pnt_b)

	(make_line_1 "0" 0 pnt_b pnt_g)
	;(command "_.arc" "C" pnt_org pnt_b pnt_e)
	(make_arc_cbe "0" 4 pnt_org pnt_b pnt_e)

	(setvar  "OSMODE" 512)	;; nearest mode
	(setq pnt_a (getpoint pnt_org "\nDefine an angle on the circle"))
	(make_line_1 "0" 3 pnt_a pnt_org)
	(make_pt "0" 1 pnt_a)  (mark_id pnt_a "A" 2 0.1)
	
	;;;compue angle AOB;;this will be used later in error estimate.
	(setq rad_aob (angle pnt_org pnt_a) deg_aob (rtd rad_aob))
	(setvar "OSMODE" 0)	;;reset Snap mode to zero


);;;SETUP_VONCUSA
;;;
;;;
(defun compute_intersection()
	(setq m (distance pnt_b pnt_m) m_sqr (* m m)
		denom (+ 9 m_sqr)
		num (- (* 3. (sqrt (- 9 (* 3 m_sqr)))) (* 2 m_sqr))
		x_pos (/ num denom) y_pos (* (/ m 3.) (+ x_pos 2.) )
		pnt_t (list x_pos y_Pos)
	)
);;;compute_intersection
;;;
;;;voncusa_error
;;;
(defun c:voncusa_error()
	(setup_error)


	(reset_sysvar)
);;;
;;;
;;;
(defun setup_error()
	(setup_sysvar)
	(setq start_angle 0. end_angle 90. step_angle 10.
		pnt_org '(0 0) pnt_a '(1 0) pnt_b '(1 1) pnt_c '(0 1)
		del_x 0.1 del_y 0.1
	)
	;;draw grid lines
	(make_line_1 "0" 8 pnt_org pnt_c)
		(setq v_line (entlast))
	(make_line_1 "0" 8 pnt_org pnt_a)
		(setq h_line (entlast))
	(command "_.array" h_line "" "R" 10 1 del_y)
	(command "_.array" v_line "" "R" 1 10 del_x) 

	(command "_.zoom" "_Extent")
	(command "_.regen")

);;setup_error
;;;
(prompt "error_table loading")
;;;
(defun C:error_table()
	(setup_sysvar)
	(setq start_val 0. end_val 90. step_val  10.)
	(setq cur_val 0. nstep 1)
	(setq er_list nil)
	(while (<= (setq cur_val (+ start_val (* (- nstep 1) step_val))) end_val)
		(setq th (dtr cur_val)
			num (sin th) denom (+ 2 (cos th))
			h (/ num denom)
			k (/ (* h h) 9.)
			num2 (- (sqrt (- 1. (* 3 k)))  (* 2 k) )
			denom2 (+ 1 k)
			x_0 (/ num2 denom2)
			x_by_y (* (/ h 3.) (/ (+ x_0 2) x_0) )
			tri_th (rtd (atan x_by_y))
			third_th (/ cur_val 3.)
			error_deg (- third_th tri_th)
			error_rad (dtr error_deg)
			error_in_minsec (angtos error_rad 1 4)
			output (list nstep cur_val third_th tri_th error_deg error_in_minsec)
			er_list (append   er_list (list (* 3600. error_deg)))
			nstep (1+ nstep)
		)
		(princ "output= ")(princ output)(terpri)
	);;;while loop end



	(reset_sysvar)
);;;error_table