;;;
;;;Rene Descartes's approach to find PAI geometrically
;;;	by Takaya Iwamoto
;;;Find a diameter of a circle the circumference of which is a given value (4.0)
;;;Start with a square with unit side length.
;;;FInd a rectangle with area 1/4 of the previous rectangle (or square for the first)
;;;
;;;Ref. "History of Pi" by Petr Beckmann
;;;Ref. "Geometrie"     by Rene Descartes
;;;Descart.lsp
;;; PI = 3.14159 26535 89793
;;;
;;;after execution 
;;;function  GET_PI will give an approximation value for pi.
;;;
;;;test1:   print the result in text windows up to 15 decimals.
;;;         in order to get precise matching
;;;         up to 5     9 iterations
;;;               10   18
;;;               15   25
;;;         This is based on the following equation.
;;;         x(k) ( x(k) - x(k-1) ) = x(0)*x(0) / 4 ^ k
;;;         
;;;test2:   Graphical presentations with different colors
;;;
;;;
;;;test3:   Graphical solution by Descartes
;;;         quadratic equation of the form 
;;;         x*x  - a*x  = b*b can be solved by a line and circle.
;;;         
;;;
(defun c:test1()
	(setup_descarte)
	(setq nstep 1)
	(setq step (getint "\nHow many steps?"))
	(while (<= nstep step)
		(prin1 x_km1)(terpri)
		(setq x_k (next_x)
			pnt_next (list x_k x_k)
			x_km1 x_k
			k (1+ k)
			pi_app (/ 4. x_k)
		)
		(prin1 (rtos pi_app 2 15))
	   (make_pt "0" 1 pnt_next)
		(setq nstep (1+ nstep))
		(prin1 nstep)(terpri)
	);while loop
	
);C:TEST1
;;;
;(defun c:descart_pi()
;	(setup_descarte)
;	(setq nstep 1)
;	(command "_.layer" "_set" "layer1" "")
;	(setq step (getint "\nHow many steps?"))
;	(while (<= nstep step)
;		(setq x_k (next_x)
;			pnt_h (list x_k x_k)
;			pnt_g (list x_km1 x_k)
;			pnt_i (list x_k 0)
;			pnt_f (list x_km1 x_km1)
;			x_km1 x_k
;			k (1+ k)
;		)
;		(command "_.pline" pnt_f pnt_g pnt_h pnt_i "")
;		(setq nstep (1+ nstep)
;			layer_name (strcat "layer" (itoa (rem nstep 8))))
;		(if (= layer_name "layer0") (setq layer_name "0"))
;		(command "_.layer" "_set" layer_name "")
;	);while loop
;	(c:get_pi)
;);C:DESCART_PI
;;;
;;;
(defun c:descartes_pi()
	(setup_descarte)
	(setq nstep 0)
	;(setvar "CECOLOR" "red")
	(alert "Pi approximation by Rene Descartes\n\n
		MAke a rectangle with an area equal to 1/4 of\n
		the immediate left")
	(while (= (Yes_or_No "Go to next step ?") "_Y")
		(setq x_k (next_x)
			pnt_h (list x_k x_k)
			pnt_g (list x_km1 x_k)
			pnt_i (list x_k 0)
			pnt_f (list x_km1 x_km1)
			x_km1 x_k
			k (1+ k)
		)
		(setq 	nstep (1+ nstep)
			col_id (itoa (rem nstep 7))
		)
		;(if (= col_id 7) (setq col_id "white"))
		(setvar "CECOLOR"  col_id)
		(command "_.pline" pnt_f pnt_g pnt_h pnt_i "")
		
		(C:get_pi)
	);while loop
	(c:get_pi)
	(result_display)
	(reset_sysvar)
);C:DESCARTES_PI
;;;
;;;
(defun c:test_2()
	(setup_descarte)
	(setq nstep 1)
	(setvar "CECOLOR" "red") 
	(setq step (getint "\nHow many steps?"))
	(while (<= nstep step)
		(setq x_k (next_x)
			pnt_h (list x_k x_k)
			pnt_g (list x_km1 x_k)
			pnt_i (list x_k 0)
			pnt_f (list x_km1 x_km1)
			x_km1 x_k
			k (1+ k)
		)
		(command "_.pline" pnt_f pnt_g pnt_h pnt_i "")
		(setq nstep (1+ nstep)
			col_id (itoa (rem nstep 7))
		)
		;(if (= col_id 7) (setq col_id "white"))
		(setvar "CECOLOR" col_id)
	);while loop
	(reset_sysvar)
	(c:get_pi)
);C:TEST_2
;;;
;;;
;;SETUP_DESCARTE
(defun setup_descarte()
	(setup_sysvar)
	(set_txstyle "arial")
	(setvar "PDMODE" 34)
	(setvar "PDSIZE" -3)
	(setq y_up '(0 2)
		x_right '(2 0)
		pnt_org '(0 0)
		pnt_b '(1 0)
		pnt_c '(1 1)
		pnt_d '(0 1)
		pnt_e '(2 2)
		chr_size 0.1
		x_km1 1.0
		k 	1
		x_0   1.0
		x_0_sqr 1.0
		nstep 1 x_loc 1.5 y_loc 1.1 ystep -0.07
	)
	(make_line_1 "0" 8 pnt_org y_up)
	(make_line_1 "0" 8 pnt_org x_right)
	(make_line_1 "0" 8 pnt_org pnt_e)
	
	(command "_.pline" pnt_org pnt_b pnt_c pnt_d "_C")

	(textdisplay "1.0" '(-0.02 0.475) 0.07 90.)
	(textdisplay "1.0" '(0.45 1.025)    0.07 0.)

	(textdisplay "Rene Descartes:"    '(1.7122  1.5732) 0.0889 0.)
	(textdisplay "Find a diameter of a circle whose" '(1.7122  1.4000) 0.0889 0.)
	;(textdisplay "a circle whose"     '(1.7122  1.0573) 0.0889 0.)
	(textdisplay "circumference is a given length (4)"   '(1.7122  1.2573) 0.0889 0.)
	;(textdisplay "a given length (4)" '(1.7122  0.7718) 0.0889 0.)

	(textdisplay "Step        Pi -- Approximation(4/Dia)              Exact Value" 
			'(1.45 1.1452) 0.0423 0.)


	(command "_.zoom" "_EXTENT")
	
);setup_descarte
;;;
;;;
;;;
(defun next_x()
	(setq denom (expt 4. (- k 1))
		x_km1_sqr (* x_km1 x_km1)
	)
	(* 0.5 (+ x_km1 (sqrt (+ x_km1_sqr (/ x_0_sqr denom)) )))
	
);NEXT_X
;;;
;;;
(defun c:get_pi()

	(setq pi_value (rtos (/ 4. x_k) 2 15))
	;(setq mes1 (strcat "Pi    ---    approximation\n" "\n (itoa nstep )  "\n\n" pi_value
	;"\nExact value is" "\n3.141592653589793" ) )
	(if (< nstep 10)
	(setq mes1  (strcat  (itoa nstep ) "\n\n\n" pi_value "\n     3.141592653589793" )) 
	(setq mes1  (strcat  (itoa nstep ) "\n\n"   pi_value "\n     3.141592653589793" )) 
	)
	(setq text_loc (list x_loc (+ y_loc (* nstep ystep))))
	;(alert mes1)
	(textdisplay mes1 text_loc 0.05 0.)
	;(command "_.polygon" "4" "_Edge" pnt_org pnt_i )
	(princ mes1)(terpri)

);C:GET_PI
;;;
;;;
(defun result_display()
	(setq pnt_j (list 0 x_k) pnt_k (mid_point pnt_j pnt_h)
		rad (distance pnt_j pnt_k)
	)
	(make_line_1 "0" 2 pnt_j pnt_h)
	(make_circle_1 "0" 2 pnt_k rad)
	(textdisplay "Pi x Diameter = 4" '(0.2  1.3264) 0.07 0.)
	(command "_.zoom" "_EXTENT")
);RESULT_DISPLAY
;;;
;;;
;;;