;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;									     ;;
;;	File:     gps.el						     ;;
;;	Author:   Wolfgang S. Rupprecht <wolfgang@charlotte.dontspam.wsrcc.com>		     ;;
;;	Created:  Sun Jun  9 07:27:22 PDT 1996				     ;;
;;	Contents: various gps calculation programs			     ;;
;;									     ;;
;;	Copyright (c) 1996 Wolfgang S. Rupprecht.				     ;;
;;	All rights reserved.						     ;;
;;									     ;;
;;	$Id: gps-el.txt,v 1.4 2003/10/20 08:27:03 wolfgang Exp $							     ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; WGS semi-major axis: 	6378137.0
;; in meters. 
;; WGS 1/f: 			    298.257223563

;; ave. radius 3,959 miles
;; km/mile 1.609344e+00
;; (* 3959.0 1.609344)
;; 6371.3928960000003

(defconst wgs84:a0 6378137.0 "Semimajor axis of the Earth in meters")
(defconst wgs84:1f0 298.257223563 "Inverse flattening")

(setq wgs84:f0 (/ 1 wgs84:1f0))
;; 0.0033528106647474805

(setq wgs84:b0 (* wgs84:a0 (- 1 wgs84:f0)))
;; 6356752.3142451793

;; (- wgs84:a0 wgs84:b0)
;;   21384.685754820704

;; dist = 	Arccos[Cos[lat1] Cos[lon1] Cos[lat2] Cos[lon2] + 
;; 	   	Cos[lat1] Sin[lon1] Cos[lat2] Sin[lon2] +
;; 		Sin[lat1] Sin[lat2]]/360 * 2Pi * r

(defun lat-lon-dist (lat1 lon1 lat2 lon2)
  "Distance between 2 pts on the earth (in meters)."
  (interactive "nLattitude 1\nnLongitude 1\nnLattitude 2\nnLongitude 2")
  (setq lat1 (degrees-to-radians lat1)
	lon1 (degrees-to-radians lon1)
	lat2 (degrees-to-radians lat2)
	lon2 (degrees-to-radians lon2))
  (* wgs84:a0
     (acos
      (+
       (* (cos lat1)
	  (cos lon1)
	  (cos lat2)
	  (cos lon2))
       (* (cos lat1)
	  (sin lon1)
	  (cos lat2)
	  (sin lon2))
       (* (sin lat1)
	  (sin lat2))))))

;; distance for 1 degree
;; (lat-lon-dist 37.0 121.0 38.0 121.0)
;; 111,319.49 m.

;; distance for 1 minute
;; (lat-lon-dist 37.0 121.0 (+ 37.0 (/ 1.0 60)) 121.0)
;; 1855.32 m.

;; distance for 1 second
;; (lat-lon-dist 37.0 121.0 (+ 37.0 (/ 1.0 3600)) 121.0)
;; 30.92 m.

;; beacon at: 
;; 37 11.22469 -122 23.39607

(setq latbeacon (hms2h 37.0 11.22469 0.0)
      lonbeacon (hms2h -122.0 23.39607 0.0))

(defun lat-lon-sa (lat lon latsa lonsa)
  "Given LATITUDE LONGITUDE and LATITUDANAL-NOISE LOGITUDINAL-NOISE
caluculate the ns and ew noise and the distance from the dgps beacon.
All figures are in meters."
  (list
   "ns:"
   (lat-lon-dist lat lon (+ lat latsa) lon)
   "ew:"
   (lat-lon-dist lat lon lat (+ lon lonsa))
   "dist:"
   (lat-lon-dist lat lon latbeacon lonbeacon)
   "meters"))


(defun decimalize-minutes ()
  "convert the previous number from minutes to fractional hours"
  (interactive)
  (save-excursion
    (if (not (looking-at "\\s_\\|\\sw")) ; if not in a symbol
	(re-search-backward "\\s_\\|\\sw" nil t)) ; go into prev. one
    (setq minutes
	  (string-to-int
	   (buffer-substring
	    (progn (backward-sexp 1) (point))
	    (setq end (progn (forward-sexp 1) (point)))
	    )))

    (setq hours
	  (string-to-int
	   (buffer-substring
	    (setq start (progn (backward-sexp 2) (point)))
	    (progn (forward-sexp 1) (point)))))

    ;; (error (cons hours minutes))
    
    (delete-region start end)
    
    (let ((float-output-format "%.3f"))
      (insert (format "%s" (+ hours (/ minutes 60.0)))))
    ))

(defun decimalize-hms ()
  "convert the previous number from minutes to fractional hours"
  (interactive)
  (save-excursion
    (if (not (looking-at "\\s_\\|\\sw")) ; if not in a symbol
	(re-search-backward "\\s_\\|\\sw" nil t)) ; go into prev. one
    (setq secs
	  (string-to-int
	   (buffer-substring
	    (progn (backward-sexp 1) (point))
	    (setq end (progn (forward-sexp 1) (point)))
	    )))

    (setq minutes
	  (string-to-int
	   (buffer-substring
	    (progn (backward-sexp 2) (point))
	    (setq end (progn (forward-sexp 1) (point)))
	    )))
    (setq hours
	  (string-to-int
	   (buffer-substring
	    (setq start (progn (backward-sexp 3) (point)))
	    (progn (forward-sexp 1) (point)))))

    ;; (error (cons hours minutes))
    
    (delete-region start end)
    
    (let ((float-output-format nil))
      (insert (format "%s" (+ hours (/ minutes 60.0) (/ secs 3600.0)))))))

(defun hms2h (h m s)
  (if (< h 0)
      (- h
	 (/ m 60.0)
	 (/ s 3600.0)
	 )
    (+ h
       (/ m 60.0)
       (/ s 3600.0)
       )))

(defun int-to-hms (x)
  (let (hour min sec (mult 1))
    (if (< x 0)
	(setq x (- x)
	      mult -1))
    (setq hours (truncate x))
    (setq x (* 60 (- x hours)))
    (setq min (truncate x))
    (setq sec (* 60 (- x min)))
    (list (* mult hours) min sec)
    ))

;;
;; end
;;