;;;Pi approxiomation by Archimedes ,then Viete
;;;		  by  Takaya Iwamoto   Jan 26,2007
;;;Ref. "A history of Pi"  by  Petr Beckmann page 62.
;;;
;;;********************************************************************
;;;		output text file name is pi_number.txt
;;;********************************************************************
;;;

;;;
;;;
(defun c:List_pi_number()
	(setup_pi_approx)
	(graph_note)
	(make_jpg)
	(setq fil (open "pi_number.txt" "w"))
	(setq k_times (getint "\nNumber of cycles ?:(def = 17)"))
	(if (= k_times nil) (setq k_times 17))

	(princ (setq title1 "\t\tPi Number List\n\n") fil)
	;(princ (setq title2 (strcat  	"\t " "k"  "\t" "lower_pi" "\t"
	;				"upper_pi" "\n")) fil)
	(setq n 6
		theta (/ pi 6.)
		nk 0
	)
	(repeat k_times
		
		(setq two_k (* 2 nk)
			exp_2k (expt 2 nk)
			theta_2k (/ theta exp_2k)
			lower_pi (* exp_2k n (sin theta_2k))
			upper_pi (* exp_2k n (tan theta_2k))
			delta_low (- pi lower_pi)
			delta_up (- upper_pi pi)
			y_low (- (+ 12 (log10 delta_low)))
			y_up  (+ 12 (log10 delta_up))
			pnt_low (list nk y_low)
			pnt_up (list nk y_up)
			app_low (rtos lower_pi 2 16)
			app_up (rtos upper_pi 2 16)
		)
		(princ "*** ")(princ nk)(princ " : ")(princ delta_low)(princ " : ")(princ delta_up)(terpri)
		;(princ y_low)(princ " : ")(princ y_up)(terpri)
		(princ app_low)(princ " : ")(princ app_up)(terpri)
		(princ " ")(terpri)
;;Archimedes result
		(if (= nk 4)
			(progn
				(textdisplay (rtos lower_pi 2 6) pnt_low 0.5 0.)
				(textdisplay (rtos upper_pi 2 6) pnt_up 0.5 0.)
			)
		)
;;Vieta result
		(if (and (> nk 7) (= (rem nk 4) 0))
			(progn
				(textdisplay (rtos lower_pi 2 11) pnt_low 0.5 0.)
				(textdisplay (rtos upper_pi 2 11) pnt_up 0.5 0.)
			)
		)
			
		(if (= nk 0) 
			(progn
				(setq pnt_low_old pnt_low 
					pnt_up_old pnt_up)
				(make_pt "0" 1 pnt_low)
				(make_pt "0" 2 pnt_up)
			)

			(progn
				(setq pnt_low_new pnt_low
					pnt_up_new pnt_up)
				(make_pt "0" 1 pnt_low)
				(make_pt "0" 2 pnt_up)
				(make_line_1 "0" 1 pnt_low_old pnt_low_new)
				(make_line_1 "0" 2 pnt_up_old pnt_up_new)
				(setq pnt_low_old pnt_low_new
					pnt_up_old pnt_up_new
				)
			);progn
		);if
			
		

		;(princ (setq final_str (strcat  "\t"
		;				nk  "  ,  " lower_pi  "  ,  "
		;				upper_pi  "\n"))	fil)
		(setq nk (1+ nk))
	(make_jpg)
	(command "_.delay" 1000)
		
	);;;repeat loop
	(make_jpg)
	;(final_result)
	(close fil)
);c:List_pi_number
;;;
;;;log with base 10--AutoLISP's log is base "e".
;;;
(defun log10( x / inv_log10)
	(setq inv_log10 (/ 1. (log 10))	)
	(princ (* (log x) inv_log10))
);;
;;;
;;;
(defun setup_pi_approx()
	(setup_sysvar)
	(setvar "PDMODE" 32)
	(setvar "PDSIZE" -2)
	(setq pnt_org '(0 0)
		x_right '(18 0)
		yaxis_up '(0 12)
		yaxis_down '(0 -12)
		upper_right '(20.5 12.5)
		lower_left '(-3.0 -12.5)
	)
	(make_line_1 "0" 8 pnt_org x_right)
		(setq horiz_line (entlast))
	(make_line_1 "0" 8 yaxis_down yaxis_up)
		(setq vert_line (entlast))
	(command "_.array" vert_line "" "_R" 1 19 1.0)
	(command "_.array" horiz_line "" "_R" 13 1 1.0)
	(command "_.array" horiz_line "" "_R" 13 1 -1.0)
	(make_pt "0" 0 pnt_org)
	(make_line_1 "0" 0 pnt_org x_right)
		
	(make_line_1 "0" 0 yaxis_down yaxis_up)
	
	(graph_note)	
	(command "_.zoom" "_W" lower_left upper_right)
	(command "_.regen")


)
;;;
;;;
(defun graph_note()
	(setq chr_size 0.5 icnt 1)
	
	(repeat 12
		(setq num (itoa (+ icnt -12))
			y_plus (list -1 icnt)
			y_minus (list -1 (- icnt))
		)
		(textdisplay num y_plus chr_size 0.)
		(textdisplay num y_minus chr_size 0.)
		(setq icnt (1+ icnt))
	)
	(textdisplay "-12" '(-1 0) chr_size 0.)
	(textdisplay "Deviation from exact ( 10E-n )" '(-2.5 -9) 1.0 90.)

	(textdisplay "Lower" '(	5.0 -4) 0.75 0.)
	(textdisplay "Bound" '(	5.0 -5) 0.75 0.)
	(textdisplay "Upper" '(	5.0 5.0	) 0.75 0.)
	(textdisplay "Bound" '(	5.0 4.0	) 0.75 0.)

	(setq icnt 2)
	(repeat 8
		(setq num (itoa icnt)
			x_loc (list icnt -0.6)
		)
		(textdisplay num x_loc chr_size 0.)
		(setq icnt (+ 2 icnt))
	)
;;color id
	(make_Line_1 "0" 3 '(4 -9) '(4 9))
	(make_line_1 "0" 4 '(16 -4) '(16 4))
	(textdisplay "Archimedes" '(5.5 10.5) 1.0 0.)
	(textdisplay "Viete" '(14.5 5.0) 1.0 0.)
	

)
;;;
;;;
(defun C:Archimedes_pi()
	(setup_pi_approx)
	(setq fil (open "pi_number.txt" "w"))
	(setq circle_center '(-9 0)
		circle_radius 5.0
		graph_center '(9 0)
		nk 0
	)
	(textdisplay "N = " '(-10.3 0.5) 0.7 0.)
	(textdisplay "< pi <" '(-10.1 -1.0) 0.5 0.)
	(make_circle_1 "0" 8 circle_center circle_radius)
	(command "_.vports" 2 "_V")
	(setvar "CVPORT" 2)
	
	(command "_.zoom" "_W" '(-15 -5.5) '(-3 5.5))
	(command "_.regen")
	 
	(setvar "CVPORT" 3)
	
	(command "_.zoom" "_W" lower_left upper_right)
	(setq chr_size 0.5 icnt 1)
	
	

	(graph_note)
	(make_jpg)
;;;*********************************************************************
    (while (= (Yes_or_No "Go to next step ?") "_Y")
    (setq npoly (getint "\nNumber of sides in regular polygon(def =6)"))
	(if (= npoly nil) (setq npoly 6))
	(if (/= nk 0) (progn
			(entdel in_scribed)
			(entdel circum_scribed)
			(entdel npoly_text)
			(entdel lower_pi_text)
			(entdel upper_pi_text)
			)
	);;endif
	(textdisplay (itoa npoly) '(-8.366 0.486) 0.75 0.)
		(setq npoly_text (entlast))
	(setvar "CECOLOR" "1")
	(command "_.polygon" npoly circle_center "_I" circle_radius)
		(setq in_scribed (entlast))
	(setvar "CECOLOR" "2")
	(command "_.polygon" npoly circle_center "_C" circle_radius)
		(setq circum_scribed (entlast))
	(command "_.list" in_scribed "")
		(setq lower_pi (* 0.1 (getvar "PERIMETER")))
	(command "_.list" circum_scribed "")
		(setq upper_pi (* 0.1 (getvar "PERIMETER")))
	(setvar "CECOLOR" "0")
	(textdisplay (rtos lower_pi 2 5) '(-12.77 -1.00) 0.5 0.)
		(setq lower_pi_text (entlast))
	(textdisplay (rtos upper_pi 2 5) '(-8.05 -1.00) 0.5 0.)
		(setq upper_pi_text (entlast))
	(setq result (strcat (itoa npoly) "," (rtos lower_pi 2 8) "," (rtos upper_pi 2 8) ))
	(princ result)(terpri)
	(setq delta_low (- pi lower_pi)
		delta_up (- upper_pi pi)
		y_low (- (+ 12 (log10 delta_low)))
		y_up  (+ 12 (log10 delta_up))
		pnt_low (list nk y_low)
		pnt_up (list nk y_up)
		app_low (rtos lower_pi 2 16)
		app_up (rtos upper_pi 2 16)
	)
	(textdisplay (rtos lower_pi 2 6) pnt_low 0.5 0.)
	(textdisplay (rtos upper_pi 2 6) pnt_up 0.5 0.)
	(if (= nk 0) 
			(progn
				(setq pnt_low_old pnt_low 
					pnt_up_old pnt_up)
				(make_pt "0" 1 pnt_low)
				(make_pt "0" 2 pnt_up)
			)

			(progn
				(setq pnt_low_new pnt_low
					pnt_up_new pnt_up)
				(make_pt "0" 1 pnt_low)
				(make_pt "0" 2 pnt_up)
				(make_line_1 "0" 1 pnt_low_old pnt_low_new)
				(make_line_1 "0" 2 pnt_up_old pnt_up_new)
				(setq pnt_low_old pnt_low_new
					pnt_up_old pnt_up_new
				)
			);progn
	);if
	(setq nk (1+ nk))
    	(make_jpg)
  );;while loop	

	(reset_sysvar)
)
;;;
;;;
;;;Huygens
;;;
(defun C:Huygens()
	(setq n_repeat 10
		pi_by_6 (/ pi 6.)
		k_cnt 0
	)
	(repeat n_repeat
		(setq c1 (expt 2 k_cnt)
			c2 (expt 2 (1+ k_cnt))
			f1 (* 6 c1)
			f2 (* 6 c2)
			theta_1 (/ pi_by_6 c1)
			theta_2 (/ pi_by_6 c2)
			sp1 (* f1 (sin theta_1))
			sp2 (* f2 (sin theta_2))
			bp1 (* f1 (tan theta_1))
			bp2 (* f2 (tan theta_2))
			x_low (/ (* 3. f1 (sin theta_1)) (+ 2. (cos theta_1)))
		)
		(setq lower_pi (/ (- (* 4. sp2) sp1) 3.)
			upper_pi (/ (+ (* 2 sp1) bp1) 3.)
			low_pi (rtos lower_pi 2 13)
			up_pi (rtos upper_pi 2 13)
			low_x (rtos x_low 2 13)
		)
		(setq k_cnt (1+ k_cnt))
		(princ k_cnt)(terpri)
		(princ low_pi)(princ " :: ")(princ up_pi)(princ " :: ")(princ low_x)(terpri)
		
	);;end of loop

)
;;;
;;;
;;;Snell_Huygens_LB
;;;
(defun c:Snell_Huygens_LB( )
	(setup_snell_huyg_lb)
	(reset_sysvar)
);;
;;;
;;;
(defun setup_snell_huyg_lb()
	(setup_sysvar)
	(setvar "PDMODE" 32)
	(setvar "PDSIZE" -4)
;;lower bound
	(setq pnt_e '(-2 0)
		pnt_b '(1 0)
		pnt_a '(-1 0)
		pnt_org '(0 0)
		pnt_f (polar pnt_org (/ pi 3.) 1.0)
		pnt_tan '(1 0.75)
		pnt_g (inters pnt_e pnt_f pnt_b pnt_tan nil)
		dist_bg (distance pnt_b pnt_g)
		lower_left '(-2.1 -0.5)
		upper_right '(1.25 1.25)
	)
;;
	(make_pt "0" 0 pnt_org)
	(make_pt "0" 0 pnt_g)
	(make_pt "0" 0 pnt_e)
	(make_pt "0" 0 pnt_a)
	(make_pt "0" 0 pnt_b)
	(make_pt "0" 0 pnt_f)

	(mark_id pnt_org "O" 4 0.1)
	(mark_id pnt_b "B" 4 0.1)
	(mark_id pnt_a "A" 4 0.1)
	(mark_id pnt_e "E" 4 0.1)
	(mark_id pnt_g "G1" 1 0.1)
	(mark_id pnt_f "F" 2 0.1)
;;
	(command "_.arc" "_C" pnt_org pnt_b pnt_e)
	(make_line_1 "0" 8 pnt_e pnt_b)
	(make_line_1 "0" 8 pnt_e pnt_g)
	(make_line_1 "0" 2 pnt_b pnt_g)
	(make_line_1 "0" 8 pnt_org pnt_f)
;;
	(textdisplay "r" '(-1.5 0) 0.1 0.)
	(textdisplay "r" '(-0.5 0) 0.1 0.)
	(textdisplay "r" '( 0.5 0) 0.1 0.)
	(textdisplay "arc BF > BG1" '(-0.75 -0.48) 0.2 0.)
	(textdisplay "Lower" '(-1.68 -0.33) 0.15 0.)
	(textdisplay "Bound" '(-1.68 -0.57) 0.15 0.)
;;
	(mark_angle pnt_g pnt_e pnt_b 0.75 "b" 0.15 1)
	(mark_angle pnt_f pnt_org pnt_b 0.35 "q" 0.15 1)
;;
	(command "_.zoom" "_W" lower_left upper_right)
	(command "_.regen")
	
)
;;;setup_snell_huyg_lb
;;;
;;;
;;;Snell_Huygens_UB
;;;
(defun c:Snell_Huygens_UB( )
	(setup_snell_huyg_ub)
	(reset_sysvar)
);;
;;;
;;;
(defun setup_snell_huyg_ub()
	(setup_sysvar)
	(setvar "PDMODE" 32)
	(setvar "PDSIZE" -4)
;;lower bound
	(setq 	pnt_b '(1 0)
		pnt_a '(-1 0)
		pnt_org '(0 0)
		pnt_f (polar pnt_org (/ pi 3.) 1.0)
		pnt_d (polar pnt_org (/ (* 8. pi) 9.) 1.0)
		pnt_tan '(1 0.75)
		pnt_e (inters pnt_d pnt_f pnt_a pnt_b nil)
		pnt_g (inters pnt_d pnt_f pnt_b pnt_tan nil) 
		dist_bg (distance pnt_b pnt_g)
		lower_left '(-2.1 -0.5)
		upper_right '(1.25 1.25)
	)
;;
	(make_pt "0" 0 pnt_org)
	(make_pt "0" 0 pnt_g)
	(make_pt "0" 0 pnt_e)
	(make_pt "0" 0 pnt_a)
	(make_pt "0" 0 pnt_b)
	(make_pt "0" 0 pnt_f)
	(make_pt "0" 0 pnt_d)

	(mark_id pnt_org "O" 4 0.1)
	(mark_id pnt_b "B" 4 0.1)
	(mark_id pnt_a "A" 4 0.1)
	(mark_id pnt_e "E" 4 0.1)
	(mark_id pnt_g "G2" 1 0.1)
	(mark_id pnt_f "F" 2 0.1)
	(mark_id pnt_d "D" 3 0.1)
;;
	(command "_.arc" "_C" pnt_org pnt_b pnt_e)
	(make_line_1 "0" 8 pnt_e pnt_b)
	(make_line_1 "0" 8 pnt_e pnt_g)
	(make_line_1 "0" 2 pnt_b pnt_g)
	(make_line_1 "0" 8 pnt_org pnt_f)
	(make_line_1 "0" 8 pnt_org pnt_d)
;;
	(textdisplay "r" '(-1.5 0.21) 0.1 0.)
	(textdisplay "r" '(-0.47 0.23) 0.1 0.)
	(textdisplay "r" '( 0.5 0) 0.1 0.)
	(textdisplay "arc BF < BG2" '(-0.75 -0.48) 0.2 0.)
	(textdisplay "Upper" '(-1.68 -0.33) 0.15 0.)
	(textdisplay "Bound" '(-1.68 -0.57) 0.15 0.)
;;
	(mark_angle pnt_g pnt_e pnt_b 0.50 "q/3" 0.15 1)
	(mark_angle pnt_f pnt_org pnt_b 0.35 "q" 0.15 1)
;;
	(command "_.zoom" "_W" lower_left upper_right)
	(command "_.regen")
	
)
;;;setup_snell_huyg_ub
;;;
;;;SNELL_PI
;;;
(defun c:Snell_pi()
	(setq k_times (getint "\nNumber of cycles ?:(def = 10)"))
	(if (= k_times nil) (setq k_times 10))
	(setq 	hexa (/ pi 6.)
		nk 0
	)
	(repeat k_times
		
		(setq two_k (* 2 nk)
			exp_2k (expt 2 nk)
			theta (/ hexa exp_2k)
			alpha (/ theta 3.)
			lbf (* 18. exp_2k)
			ubf (* 6. exp_2k)
			lb_num (sin theta)
			lb_denom (+ 2 (cos theta))
			ub_1 (+ (* 2. (cos alpha)) 1)
			ub_2 (tan alpha)
			lower_pi (/ (* lbf lb_num) lb_denom)
			upper_pi (* ubf ub_1 ub_2)

			
			delta_low (- pi lower_pi)
			delta_up (- upper_pi pi)
			
			app_low (rtos lower_pi 2 16)
			app_up (rtos upper_pi 2 16)
		)
		(princ "** ")(princ nk)(princ " **")(terpri)
		(princ app_low)(princ " : ")(princ app_up)(terpri)
		(princ " : ")(princ delta_low)(princ " : ")(princ delta_up)(terpri)
		(princ " ")(terpri)
		(setq nk (1+ nk))
	)


)