(prompt "\nEratosthenes Model in ACAD - Eratosthenes_MODEL2.LSP")
(prompt "\nCommands:  Measure_earth ")
(prompt "\nFunctions: (SETUP_Eratosthenes2) (Add_pole_well) ")
(prompt "\nFunctions: (three_ports) (shadow_ref) (add_text)")
(prompt "\n           ")
;----------------------------------------------------------
;;;**************************************************************
;;;	Copyright	July 2006	Takaya Iwamoto
;;;**************************************************************
;;;Greek Astronomy
;;;By Takaya Iwamoto   April 30, 2005
;;;Eratosthenes_MODEL2.LSP
;;;Dimensions of sundial & well is increased by 1.e+3 times.
;;;1.e+1 yields a reasonable rendering result. We must use a high precision rendering
;;;scheme to make the real scale model possible.--->Maybe  3D studio or Viz is a solution.
;;;The longitude and latitude of Syene(present day Aswan)  &  Alexandria are
;;;		ALexandria	31.2 N    	29.9 E
;;;		Syene		24.088 N	32.899 E
;;;But in this model data for Syene is changed to
;;;				23.5 N		29.9E
;;;
;;;coordinate value in Cartesian (X,Y,Z)
;	x = R * cos(phi) * cos(theta)
;	y = R * cos(phi) * sin(theta)
;	z = R * sin(phi)
;  where
;	phi is latitude angle (north is plus, south is minus.
;	theta is longitude angle (east is plus, west is minus)
;	R is the radius of the Earth (6378 km  or  3956 mile)
;;;
;;;MEASURE_EARTH
;;;main trial routine
(defun c:measure_earth()
	(setup_eratosthenes2)		;setup parameters
	(set_UCS)
	(alert "\nAdd Pole and well")
	(add_pole_well)
	(alert "\nSplit into 3 View ports")
	(three_ports)
	(alert "\nMeasure the angle\nThis is equal to\nthe difference in \nlatitude")
	(shadow_ref)
	;(add_text)	;;add descriptions
);measure_earth
;;;

;;;
(defun set_UCS()
	;;;define UCS named UCS_A
	(setq pnt_x1 (list (/ earth_rad (cos phi_a)) 0 0)
		pnt_y1 (list (/ earth_rad (cos phi_a)) earth_rad 0)
	)
	(alert "\nDefine UCS named UCS_A\n
		at ALexandria,Egypt")
	(command "_.ucs" "_New" "3point" pnt_alex pnt_x1 pnt_y1)
	(command "_.ucsicon" "_ON")
	(command "_.ucs" "_Save" "UCS_A")
	(alert "\nDefine another UCS named UCS_S\n at Syene,Egypt")
	(command "_.ucs" "_World")
	;;;define UCS named UCS_S
	(setq pnt_x1 (list (/ earth_rad (cos phi_s)) 0 0)
		pnt_y1 (list (/ earth_rad (cos phi_s)) earth_rad 0)
	)
	(command "_.ucs" "_New" "3point" pnt_syen pnt_x1 pnt_y1)
	
	(command "_.ucs" "_Save" "UCS_S")
	(alert "\nGo back to world coordinate")
	(command "_.ucs" "_World")

);;;set_UCS
;;;
;;;Add post at Alexandria and drill a well at Syene
(prompt "\nLoading add_pole_well")
;;;
(defun Add_pole_well()
	(setvar "SURFTAB1" 36)
	(setvar "SURFTAB2" 36)
	(setvar "FACETRES" 10)
	(alert "\nMake a sundial at Alexandria")
	(command "_.ucs" "_Restore" "UCS_A")
	;;;Create a sundial at Alexandria
	(command "_.pline" '(0 0) '(3 0) '(3 1) '(0.1 1) '(0.1 3) '(0 3) "_cl")
	(setq sund_line (entlast))
	(command "_.revolve" sund_line "" "_Object" line_alex 360)
	(setq sundial (entlast))

	(alert "\nMake a well at Syene")
	(command "_.ucs" "_Restore" "UCS_S")
	;;;Create a well at Syene
	;;;Create a well ID=2.0  OD=3.0 Depth 5 m,rim height=1 m
	(command "_.pline" '(0 -5) '(1 -5) '(1 1) '(1.5 1) '(1.5 -6) '(0 -6) "_cl")
	(setq well_line (entlast))
	(command "_.revolve" well_line "" "_Object" line_syen 360)
	(setq well (entlast))

	;;;water in the well displayed in CYAN color
	(setvar "CECOLOR" "CYAN")
	(command "_.pline" '(0 -5) '(1 -5) '(1 -4) '(0 -4) "_CL")
	(command "_.revolve" (entlast) "" "_Object" line_syen 360)
	;;;rim around the well is covered by mud
	(setvar "CECOLOR" "42")
	(command "_.pline" '(1.5 -1) '(4 -1) '(4 0) '(1.5 0) "_CL")
	(command "_.revolve" (entlast) "" "_Object" line_syen 360)
	(setvar "CECOLOR" "BYLAYER")
	;(command "_.sphere" '(0 50 0) 0.5 )

	(command "_.ucs" "_World")


	;;;define UCS named UCS_A
	;(setq pnt_x1 (list (/ earth_rad (cos phi_a)) 0 0)
	;	pnt_y1 (plt_3d earth_ctr pnt_alex 1.0)
	;)
	;(command "_.ucs" "3point" pnt_alex pnt_x1 pnt_y1)
	;(command "_.ucs" 
	;(command "_.cylinder" pnt_alex pole_base_rad "C" pole_base_top)
	;	(setq pole_base (entlast))
	;(command "_.cylinder" pole_base_top pole_rad "C" pole_top)
	;	(setq pole (entlast))
	
	;(command "_.cylinder" in_shell_bot in_shell_rad  "C" in_shell_top)
	;	(setq in_shell (entlast))
	;(command "_.cylinder" out_shell_bot out_shell_rad  "C" out_shell_top)
	;	(setq out_shell (entlast))
	;(command "_.subtract" out_shell "" in_shell "")
	;	(setq well (entlast))

);;;ADD_POLE_WELL
;;;
;;;
(defun shadow_ref()
	;;Alexandria is vport #3
	;;;VIEWDIR = 5.5913,5.0841,16.1236
	(setvar "CVPORT" 3)
	(command "_.ucs" "_Restore" "UCS_A")
	
	(command "_.trace" 0.05   '(-0.1 3 0) '(-0.37 1 0) "")	;;;shadow ref line
	
	(command "_.regen")

)
;;;
;;;
(prompt "\nLoading three_ports")
;;;Split views into 3 
;;;
(defun three_ports()
	(setq alex_low_left '(5455.5091 3303.97450 0)
	      alex_up_right '(5455.5167 3303.9788 0)
	      syen_up_right '(5826.5963 2594.1720 0)
	      syen_low_left '(5826.5671 2594.1529 0)
	      all_low_left  '(-8329.4751 -6909.1081 0)
	      all_up_right  '( 7456.7067 10493.3440 0)
	);;;each view port's limit

	(command "_.vports" 3 "_R")
	;;Alexandria is vport #3
	;;;VIEWDIR = 5.5913,5.0841,16.1236
	(setvar "CVPORT" 3)
	(command "_.viewres" "_N" 10000)
	(command "_.zoom" "_W" alex_low_left alex_up_right)
	(command "_.trace" 0.05 '(-0.1 3 0) '(-0.37 1 0) "")	;;;shadow ref line
	;(command "_.zoom" "_Object" sundial "")
	;(command "_.vpoint" "5.5913,5.0841,16.1236")
	(command "_.regen")

	;;Syene is vport #4
	(setvar "CVPORT" 4)
	(command "_.viewres" "_N" 10000)
	(command "_.zoom" "_W" syen_low_left syen_up_right)
	;(command "_.zoom" "_Object" well "")
	;(command "_.vpoint" "17.2469,3.4575,3.1594")
	(command "_.regen")

	;;Whole picture is vport #2
	(setvar "CVPORT" 2)
	(command "_.zoom" "_W" all_low_left all_up_right)
	(command "_.regen")

	(alert "\ndisplay shaded object in View port 3 & 4")
	;;;display shaded object in View port 3 & 4
	;;Alexandria is vport #3
	(setvar "CVPORT" 3)
	(command "_.zoom" "_Object" sundial "")
	(command "_.regen")

	;;Syene is vport #4
	(setvar "CVPORT" 4)
	(command "_.zoom" "_Object" well "")
	(command "_.regen")
	
)
;;;
;;;SETUP_Eratosthenes2
(prompt "\nLoading setup_eratosthenes2")
;;;setup for Eratosthenes Model for estimating Earth's Equatorial perimeter
(defun setup_Eratosthenes2()
	(setup_sysvar)
	(setvar "SURFTAB1" 36)
	(setvar "SURFTAB2" 36)
	;(setvar "PDMODE" 32)
	(setq sun_rad 6.95e5 moon_rad 1.734e3 Earth_rad 6.378e3
		sun_dist 1.496e8 moon_dist 3.844e5 deg_23_5 (dtr 23.5)
		earth_ctr '(0 0 0) deg_20 (dtr 20.) deg_30 (dtr 30.)
		lat_20_rad (* earth_rad (cos deg_20))		;;Latitude  20 deg
		lat_20_ctr (list 0 0 (* earth_rad (sin deg_20)) )
		lat_30_rad (* earth_rad (cos deg_30))		;;Latitude  30 deg
		lat_30_ctr (list 0 0 (* earth_rad (sin deg_30)) )
		tropic_cancer_rad (* earth_rad (cos deg_23_5))	;;Tropic of Cancer
		tropic_cancer_ctr (list 0 0 (* earth_rad (sin deg_23_5)) )
		torus_rad 15 point_rad 50
		phi_a (dtr 31.2) phi_s (dtr 23.5) theta_a (dtr 29.9) theta_s (dtr 29.9)
		x_a (* earth_rad  (cos phi_a))
		y_a (* earth_rad  (sin phi_a))
		x_s (* earth_rad  (cos phi_s))
		y_s (* earth_rad  (sin phi_s))
		pnt_alex (list x_a y_a 0 )
		pnt_syen (list x_s y_s 0)
		in_shell_top (pld_3d pnt_syen earth_ctr -1.e-3)
		in_shell_bot (pld_3d pnt_syen earth_ctr 1.9e-2)
		out_shell_top in_shell_top
		out_shell_bot (pld_3d pnt_syen earth_ctr 2.4e-2)
		in_shell_rad	 1.0e-3  
		out_shell_rad	 1.5e-3
		
		pole_base_bot pnt_alex
		pole_base_top (pld_3d pnt_alex earth_ctr -1.0e-3)
		pole_bot pole_base_top
		pole_top (pld_3d pnt_alex earth_ctr -3.0e-3)
		pole_rad	 2.5e-5
		pole_base_rad	 2.0e-3
		
		
	)
	
	(make_pt "0" 5 earth_ctr)
	(make_line_1 "0" 8 earth_ctr pnt_alex)
		(setq line_alex (entlast))
	(make_line_1 "0" 8 earth_ctr pnt_syen)
		(setq line_syen (entlast))	
	(make_line_1 "0" 0 earth_ctr (list earth_rad 0 0))
	(make_line_1 "0" 0 earth_ctr (list 0 earth_rad 0))

	
	
	;(setvar "shadedge" 3)
	
	(make_circle_1 "0" 5 earth_ctr earth_rad)
	(setq earth (entlast))
	
	(command "_.zoom" "_EXTENT" )
	;(command "_.shade")
	(command "_.regen")
	
);SETUP_ERATOSTHENES2
;;;
(prompt "\nLoading add_text")
;;;
(defun add_text()
	(command "_.ucs" "_World")
	;;Earth is vport #2
	(setvar "CVPORT" 2)
	(textdisplay "Alexandria"  	'( 1630.7278  3416.8348  0)  629.0043 0.)
	(textdisplay "Syene"  		'( 5058.4761  1390.0431  0)  629.0043 0.)
	(textdisplay "Eratosthenes Model"  		'( -5959.2865  8658.5372  0)  1118.7768 0.)
	(textdisplay "N"  		'( -432.9167  6841.4137  0)  663.9490 0.)
	(textdisplay "Equator Plane"  	'( 3764.3263  -846.4166  0)  663.9490 0.)
	(textdisplay "Earth"  		'( -2115.2814  -3676.7250  0)  1327.8980 0.)
	(command "_.zoom" "_Extent")

	;;Alexandria is vport #3
	(setvar "CVPORT" 3)
	(textdisplay "Sundial"  	'( 5450.3212  3305.8278  0)  0.6196 0.)
	(textdisplay "at"  		'( 5450.3212  3304.8325  0)  0.6196 0.)
	(textdisplay "Alexandria"  	'( 5450.3212  3303.8371  0)  0.6196 0.)

	;;Syene is vport #4
	(setvar "CVPORT" 4)
	(textdisplay "Well"  		'( 5840.9478  2544.8444  0)  0.5657 0.)
	(textdisplay "at"  		'( 5840.9478  2543.9356  0)  0.5657 0.)
	(textdisplay "Syene"  		'( 5840.9478  2543.0268  0)  0.5657 0.)

	

);;add_text

(prompt "\nLoading pld_3d")
;;;
;;;************************************************************************
;;; Function: PLT_3D
;;;
;;; get a point off the line segment defined by two end points and angle.
;;;	
; utility routine 
;; get a point along the line segment defined by two end points 
;;
;; pnt_a pnt_b      end points
;; loc_from_a       ratio of distance from a to length AB
;; output           point data (x,y)
(defun PLT_3D
	(pnt_a 
	pnt_b 
	loc_from_a  
	/ 
	xa xb ya yb za zb dist_ab dx dy dz xm ym zm
	)

	(setq xa (car pnt_a) ya (cadr pnt_a) za (caddr pnt_a)
		xb (car pnt_b) yb (cadr pnt_b) zb (caddr pnt_b)
		dx (- xb xa) dy (- yb ya) dz (- zb za)
		xm (+ xa (* loc_from_a dx))
		ym (+ ya (* loc_from_a dy))
		zm (+ za (* loc_from_a dz))
	)
	(list xm ym zm)
);PLT_3D
;;;************************************************************************
;;;************************************************************************
;;; Function: PLD_3d
;;;
;;; get a point off the line segment defined by two end points and angle.
;;;	
; utility routine 
;; get a point along the line segment defined by two end points 
;;
;; pnt_a pnt_b      end points
;; dist_from_a      distance from a to b
;;
;; output           point data (x,y)
;;
(defun PLD_3D
	(pnt_a 
	pnt_b 
	dist_from_a  
	/ 
	xa xb ya yb za zb dx dy dz xm ym zm loc_from_a
	)

	(setq xa (car pnt_a) ya (cadr pnt_a) za (caddr pnt_a)
		xb (car pnt_b) yb (cadr pnt_b) zb (caddr pnt_b)
		dx (- xb xa) dy (- yb ya) dz (- zb za)
		loc_from_a (/ dist_from_a (distance pnt_a pnt_b))
		xm (+ xa (* loc_from_a dx))
		ym (+ ya (* loc_from_a dy))
		zm (+ za (* loc_from_a dz))
	)
	;(princ zm)
	(list xm ym zm)
	
);PLD_3D
(princ)