;;;PROGRAM NAME: hAVERSINE.lsp
;;;Date Started: Oct 20,2004
;;;     Finished Nov 02,2004
;;;
;;;Note:Final shading has to be done manually. don't know why.
;;;Author:	Takaya Iwamoto
;;;Reference:	Bob Chamberlain 	www.census.gov/cgi-bin/geo/gisfaq?Q5.1
;;;					"What is the best way to calculate the distance between 2 points ?"
;;;		Latitude & Longitude of Tokyo and San Francisco
;;;				Tokyo	35:41:00N(35.683)	139:44:00E(139.733)
;;;				SFO	37:47:36N(37.793)	122:33:17W(-122.555)
;;;				New York	40:40:11N(40.670)	73:56:38W(-73.944)
;;;				Seattle	47:37:18N(47.622)	122:33:17W(-122.350)
;;;				LA	34:06:44N(35.683)	118:24:40W(-118.411)
;;;				Honolulu	21:24N(21.4)	157:54W(-157.9)
;;;				Anchorage	61:10:42N(61.178)	149:11:11W(-149.186)
;;;				Misawa	41:41N(41.683)	141:24E(141.24)
;;;				Wenatchee	47:25N(47.417)	120:19W(-120.317)
;;;				Paris	48:52N(48.867)	2:20E(2.333)
;;;				London	51:30N(51.5)	0:10W(-0.167)
;;;				??Greenwich	37:47:36N(37.793)	122:33:17W(-122.555)
;;;				Peking	39:55N(39.917)	116:25E(116.417)
;;;				SFO	37:47:36N(37.793)	122:33:17W(-122.555)
;;;				ALexandria,Egypt	31.2N	29.9E
;;;				Syene(Aswan),Egypt	24.088N	32.899E
;;;					http://geography.about.com/cs/latitudelongitude/
;;;
;;;
(defun C:Haversine()
;Intuitive solution	
;	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 (6367 km  or  3956 mile)
;The distance as a crow flies is computed as follows:
;	Compute the linear distance between two places using the cartesian coordinates values,x,y,& z.
;	Usng the Radius and thi distance, compute the arc angle in radian.
;	The surface distance is the product of this arc angle and the radius.
;
	(setup_sysvar)
	(setq Erad 6367.)	;;Radius of the Earth  (Min 6336 km(3937 mile)  Max 6399 km(3976 mile))
	(setq dir_ref (* 1.05 Erad))	;reference length to define 4 directions (NSEW)
	(setq pnt_n (list 0 0 dir_ref) pnt_s (list 0 0 (- dir_ref))
		pnt_e (list 0 dir_ref 0) pnt_w (list 0 (- dir_ref) 0 )
		pnt_o '(0 0 0) pnt_g (list dir_ref 0 0)
	)
	(command "cecolor" "green")
	(command "torus" pnt_o Erad 50)
	(command "rotate3d" (entlast) "" "X" "" 90.)
	;(setq long_line (entlast))
	(command "cecolor" "BYLAYER")
	(plot_cities)
	(setvar "shadedge" 3)
	(Setvar "cecolor" "blue")
	(command "sphere"  pnt_o Erad)
	(setq whole_earth (entlast))
	(setvar "CECOLOR" "red")
	;(make_circle_1 "0" 1 pnt_o dir_ref)
	(command "_.torus" pnt_o Erad 50)
	(setq equator (entlast))
	;(command "region" (entlast))
	(command "cecolor" "BYLAYER")
	;(command "_.vpoint" "6367<-157.9<21.4")	;;Honolulu,Hawaii
	
	(command "_.zoom" "E")
	(command "cecolor" 2)	;;yellow
	(alert "Input 2 points on this sphere.")
		
	(setq phi_1 (getreal "\nLatitude angle of the point A (35.683)"))
	(if (= phi_1 nil) (setq phi_1 35.683))
	(setq theta_1 (getreal "\nLongitude angle of the point A (139.733)"))
	(if (= theta_1 nil) (setq theta_1 139.733))
	(setq phi_a (dtr phi_1) theta_a (dtr theta_1)
	     	x_a (* Erad (cos phi_a) (cos theta_a))
		y_a (* Erad (cos phi_a) (sin theta_a))
		z_a (* Erad (sin phi_a))
		pnt_a (list x_a y_a z_a)
	)
	(command "_.sphere" pnt_a 150.)
	(setq sphere_a (entlast))
	(setq phi_2 (getreal "\nLatitude angle of the point B (37.793)"))
	(if (= phi_2 nil) (setq phi_2 37.793))
	(setq theta_2 (getreal "\nLongitude angle of the point B (-122.555)"))
	(if (= theta_2 nil) (setq theta_2 -122.555))
	(setq	phi_b (dtr phi_2) theta_b (dtr theta_2)
		x_b (* Erad (cos phi_b) (cos theta_b))
		y_b (* Erad (cos phi_b) (sin theta_b))
		z_b (* Erad (sin phi_b))
		pnt_b (list x_b y_b z_b)
		dist_ab (distance pnt_a pnt_b)
	)
	(command "_.sphere" pnt_b 150.)
	(setq sphere_b (entlast))
	(command "_.cecolor" "BYLAYER")
	(alert "\nInput View point ")
	(setq phi_3 (getreal "\nLatitude angle of the point C (21.4)"))
	(if (= phi_3 nil) (setq phi_3 21.4))
	(setq theta_3 (getreal "\nLongitude angle of the point C (-157.9)"))
	(if (= theta_3 nil) (setq theta_3 -157.9))
	(setq phi_c (dtr phi_3) theta_c (dtr theta_3)
	     	x_c (* Erad (cos phi_c) (cos theta_c))
		y_c (* Erad (cos phi_c) (sin theta_c))
		z_c (* Erad (sin phi_c))
		pnt_c (list x_c y_c z_c)
	)
	(command "_.sphere" pnt_a 150.)
	(alert "slice the Earth by a great circle \npassing through 2 points and the center.")
;	
	(setq sin_alpha (/ (/ dist_ab 2.0) Erad)
		alpha (asin sin_alpha)
		arc_ab (* 2. alpha Erad)
	)
;distance between A and C
	(setq   dist_ac (distance pnt_a pnt_c)
		sin_beta (/ (/ dist_ac 2.0) Erad)
		beta (asin sin_beta)
		arc_ac (* 2. beta Erad)
	)
;distance between C and B
	(setq   dist_cb (distance pnt_c pnt_b)
		sin_gamma (/ (/ dist_cb 2.0) Erad)
		gamma (asin sin_gamma)
		arc_cb (* 2. gamma Erad)
	)

	;(princ arc_len)(terpri)
	
	
;
;Haversine Formula
;Ref. R.W.Sinnot,"Virtues of the Haversine",Sky and Telescope,vol.68,no.2,1984,p.159
;_
;	dlon = lon2 - lon1
;	dlat = lat2 - lat1
;	a = (sin(dlat/2))^2 + cos(lat1)*cos(lat2)*(sin(dlon/2))^2
;	c = 2 * arcsin(min(1,sqrt(a)))
;	d = R * c
;or using arctan instead of arcsin
;	c = 2 * arctan(sqrt(a)/sqrt(1-a))	;make sure (1-a) is not zero.
;
	(setq dlon (- theta_b theta_a) 
		dlat  (- phi_b phi_a) a1 (expt (sin (* 0.5 dlat)) 2)
		a2 (expt (sin (* 0.5 dlon)) 2)
		a (+ a1 (* (cos phi_a) (cos phi_b) a2))
		c (* 2. (atan (/ (sqrt a) (sqrt (- 1. a)) ) ))
		d (* Erad c)
	)	
	
	(command "_.shade")
	(setq shade_earth (entlast))
	(command "cecolor" "red")
	(command "_.line" pnt_a pnt_b pnt_o "c")
	(command "_.line" pnt_a pnt_c pnt_o "c")
	(command "_.line" pnt_c pnt_b pnt_o "c")
	(command "_.cecolor" "BYLAYER")
	(command "_.slice" shade_earth whole_earth sphere_a sphere_b "" "" pnt_o pnt_a pnt_b "B")	
	
	(princ "Distance between A and B = ")(princ arc_ab) (princ "km")(terpri)
	(princ "Distance between A and C = ")(princ arc_ac) (princ "km")(terpri)
	(princ "Distance between C and B = ")(princ arc_cb) (princ "km")(terpri)
	
	;(command "cecolor" "red")
	;(command "region" equator "")
	;(command "cecolor" "green")
	;(command "region" long_line "")
	(alert "Set the viewpoint")
	(command "_.vpoint" pnt_c)	;;Honolulu,Hawaii
	(alert "Use AutoCAD SHADE command to shade the final status")
	;(command "shade")
	
	;(reset_sysvar)


);;;HAVERSINE
;;;
;;;Plot_cities
;;;Plot othe cities for reference
;;;
(defun plot_cities()
	
	(command "sphere" "6367.<-157.9<21.4" 150.)	;Honolulu
	(command "sphere" "6367.<-73.944<40.670" 150.)	;New York
	(command "sphere" "6367.<-122.350<47.622" 150.)	;Seattle
	(command "sphere" "6367.<-118.411<34.112" 150.)	;Los Angeles
	(command "sphere" "6367.<-149.186<61.178" 150.)	;Anchorage
	(command "sphere" "6367.<-0.167<51.5" 150.)	;London
	(command "sphere" "6367.<116.417<39.917" 150.)	;Peking
	(command "sphere" "6367.<-95.367<29.767" 150.)	;Houston
	(command "sphere" "6367.<-90.083<29.95" 150.)	;New Orleans
	(command "sphere" "6367.<-80.2<25.767" 150.)	;Miami,Fl
	(command "sphere" "6367.<-71.067<42.367" 150.)	;Boston
	(command "sphere" "6367.<-87.39<41.85" 150.)	;Chicago
	;(command "sphere" "6367.<-0.167<51.5" 150.)	;London
	;(command "sphere" "6367.<116.417<39.917" 150.)	;Peking


);;;PLot_cities
;;;
;;;Pythagorean theorem in spherical coordinate
;;;
(defun C:sphere_pythagorean()
	(setup_sysvar)
	(setq Erad 6367.
		chr_size 1000.
		E_center '(0 0 0)
		North (list 0 0 (* 1.2 Erad))
		South (list 0 0 (- Erad))
		East (list 0 (* 1.2 Erad) 0)
		West (list 0 (- Erad) 0)
	
	)
	(setvar "ISOLINES" 16)
;;point A	
	(setq phi_1 0. theta_1 110.)
	(setq phi_a (dtr phi_1) theta_a (dtr theta_1)
	     	x_a (* Erad (cos phi_a) (cos theta_a))
		y_a (* Erad (cos phi_a) (sin theta_a))
		z_a (* Erad (sin phi_a))
		pnt_a (list x_a y_a z_a)
	)
;;point B
	(setq phi_1 30. theta_1 180.)
	(setq phi_a (dtr phi_1) theta_a (dtr theta_1)
	     	x_a (* Erad (cos phi_a) (cos theta_a))
		y_a (* Erad (cos phi_a) (sin theta_a))
		z_a (* Erad (sin phi_a))
		pnt_b (list x_a y_a z_a)
	)
;;point C
	(setq phi_1 0. theta_1 180.)
	(setq phi_a (dtr phi_1) theta_a (dtr theta_1)
	     	x_a (* Erad (cos phi_a) (cos theta_a))
		y_a (* Erad (cos phi_a) (sin theta_a))
		z_a (* Erad (sin phi_a))
		pnt_c (list x_a y_a z_a)
	)
	(make_line_1 "0" 6 E_center East)
	(make_line_1 "0" 6 E_center North)
	(textdisplay "E" East chr_size -90.)
	(textdisplay "N" North chr_size -90.)

	(command "_.sphere" pnt_c 150.)	;Point C
	(command "_.sphere" pnt_b 150.)	;Point B
	(command "_.sphere" pnt_a 150.)   ;Point A
	(command "_.sphere" E_center 150.)   ;E_center

	(textdisplay "A" pnt_a chr_size -90.)
	(textdisplay "B" pnt_b chr_size -90.)
	(textdisplay "C" pnt_c chr_size -90.)
	(textdisplay "O" E_center chr_size -90.)

	(make_circle_1 "0" 1 E_center Erad)
	(make_circle_1 "0" 2 E_center Erad)
		(setq long_zero (entlast))
	(make_circle_1 "0" 3 E_center Erad)
		(setq long_90 (entlast))
	(command "_.rotate3d" long_zero ""  "X" E_center 90.)
	(command "_.rotate3d" long_90 ""  "Y" E_center 90.)
	(alert "\nDraw a sphere")
;;make a sphere with radius Erad
	(command "_.sphere" E_center Erad)
		(setq E_ball (entlast))
	(command "_.vpoint" '(-31315. 24955. 20252.))
	(command "_.zoom" "_E")
	(alert "\nSlice the whole sphere\nby OAB plane")
;;slice this by plane OAB
	(command "_.slice" E_ball "" E_center pnt_a pnt_b South)
		(setq half_ball (entlast))
	(alert "\nSlice the remaining half-sphere\nby OBC plane")
;;slice again by plane OBC
	(command "_.slice" half_ball "" E_center pnt_c pnt_b East)
		(setq part_sphere (entlast))
	(alert "\nSlice the remaining part_sphere\nby OBC plane")
;;Then finally by plane OAC
	(command "_.slice" part_sphere "" E_center pnt_c pnt_a North)

	(reset_sysvar)
);;;
;;;
;;;CHECK_THEOREM	check the spherical Pythagorean theorem given R, a,b,&c
;;;
(defun C:check_theorem()
	(setq len_a 3333.7534
		len_b 7778.7579
		len_c 8086.6416
		cos_c (cos (/ len_c Erad))
		cos_a (cos (/ len_a Erad))
		cos_b (cos (/ len_b Erad))
		cos_ab (* cos_a cos_b)
		ratio_1 (/ cos_c cos_ab)
		sqr_a (* len_a len_a)
		sqr_b (* len_b len_b)
		sqr_c (* len_c len_c)
		sqr_ab (+ sqr_a sqr_b)
		ratio_2 (/ sqr_c sqr_ab)
		
	)
);;;
;;check_theorem