From: igb@fulcrum.co.uk (Ian G Batten)
Newsgroups: uk.misc
Subject: Re: Lat/Long Co-Ordinates Conversion to NGR source code wanted
Date: 31 Dec 92 09:56:07 GMT
Organization: Fulcrum Communications

In article <725770662snx@skewsby.demon.co.uk> sjh@skewsby.demon.co.uk writes:
> I now find myself in the position of needing C source code or even COBOL, to
> do this conversion.

Man of your talents can convert.  Get the leaflet.  This was written in
a hurry and works on CMU Common Lisp 16d.

;;; --- cut here ---

;;; Convert from lat/long pair to national grid.

;;; Reference: ``Transverse Mercator Projection Constants, Formulae
;;; and Methods'' Ordnance Survey March 1983.

;;; the scale factor

(defconstant f0 0.9996012717)

;;; The spheroid

(defconstant a (* f0 6377563.396))	; major semi-axis in metres 
(defconstant b (* f0 6356256.910))	; minor ditto

;;; handy constants

;; should be 0.006670539762

(defconstant e-squared (/ (- (* a a) (* b b)) (* a a)))

(defconstant n (/ (- a b) (+ a b)))

(defun nu (phi) 
  "radius of curvature at lattitude phi perpendicular to a meridian"
  (/ a 
     (expt (- 1 (* e-squared (expt (sin phi) 2))) 0.5)))

(defun rho (phi)
  "radious of curvature of a meridian at littitude phi"
  (/ (* (nu phi) (- 1 e-squared))
     (- 1 (* e-squared (expt (sin phi) 2)))))

(defun eta-squared (phi)
  (1- (/ (nu phi) (rho phi))))

(defun deg2rad (degrees &optional (minutes 0) (seconds 0))
  (/ (+ degrees (/ minutes 60) (/ seconds 3600))
     (/ 180 pi)))

;;; origin of national grid

(defconstant phi-zero (deg2rad 49))	; degrees north
(defconstant lambda-zero (deg2rad -2))	; degree east (neg is west)

;;; grid co-ordinates of origin

(defconstant E0 400000)			; metres east
(defconstant N0 -100000)			; metres north (neg is south)

(defun M (phi)
  "Developed arc of a meridian from phi to phi-zero"
  (let* ((n-squared (expt n 2))
	 (n-cubed (expt n 3))
	 (term-1 (* (+ 1 n (* (/ 5 4) n-squared) (* (/ 5 4) n-cubed)) 
		    (- phi phi-zero)))
	 (term-2 (* (+ (* 3 n) (* 3 n-squared) (* (/ 21 8) n-cubed))
		    (sin (- phi phi-zero))
		    (cos (+ phi phi-zero))))
	 (term-3 (* (+ (* (/ 15 8) n-squared) (* (/ 15 8) n-cubed))
		    (sin (* 2 (- phi phi-zero)))
		    (cos (* 2 (+ phi phi-zero)))))
	 (term-4 (* (/ 35 24) n-cubed
		    (sin (* 3 (- phi phi-zero)))
		    (cos (* 3 (+ phi phi-zero))))))
    (* b (+ term-1 (- term-2) term-3 (- term-4)))))

(defun grid (phi lambda)
  "(E N) from phi and lambda"
  (let* ((nu-phi (nu phi))
	 (sin-phi (sin phi))
	 (cos-phi (cos phi))
	 (tan-phi (tan phi))
	 (eta-squared-phi (eta-squared phi))
	 (I (+ (M phi) N0))
	 (II (* (/ nu-phi 2) sin-phi cos-phi))
	 (III (* (/ nu-phi 24) sin-phi (expt cos-phi 3)
		 (- 5 (expt tan-phi 2) (* 9 eta-squared-phi))))
	 (IIIA (* (/ nu-phi 720) sin-phi (expt cos-phi 5)
		  (+ 61 (- (expt tan-phi 2)) (expt tan-phi 4))))
	 (P (- lambda lambda-zero))
	 ;; at last
	 (north (+ I (* (expt P 2) II) (* (expt P 4) III) 
	       (* (expt P 6) IIIA)))
	 (IV (* nu-phi cos-phi))
	 (V (* (/ nu-phi 6) (expt cos-phi 3)
	       (- (/ nu-phi (rho phi)) (expt tan-phi 2))))
	 (VI (* (/ nu-phi 120) (expt cos-phi 5)
		(+ 5 (* -18 (expt tan-phi 2)) (expt tan-phi 4)
		   (* 14 eta-squared-phi) 
		   (* -58 (expt tan-phi 2) eta-squared-phi))))
	 ;; even more at last
	 (east (+ E0 (* P IV) (* (expt P 3) V) (* (expt P 5) VI))))
    (list east north)))

(defun pair-to-ref (pair)
  "Convert (east north) to a standard six-figure grid reference"
  (let ((east (car pair))
	(north (cadr pair)))
    (format nil "[~a~3,'0d~3,'0d]" 
	    (square east north) 
	    (round (/ (mod east small-square-size) 100))
	    (round (/ (mod north small-square-size) 100)))))

;; ABCDE
;; FGHJK
;; LMNOP
;; QRSTU
;; VWXYZ

(defconstant square-size 500000)	;500km in metres
(defconstant small-square-size 100000)	;100km in metres

(defun square-1 (x y)
  (char "VWXYZQRSTULMNOPFGHJKABCDE" (+ x (* 5 y))))

(defun square (east north)
  (let* (
	 ;; get the large square corners
	 (east-square (+ (floor (/ east square-size)) 2))
	 (north-square (+ (floor (/ north square-size)) 1))
	 ;; now the small ones
	 (small-east-square
	  (floor (/ (mod east square-size) small-square-size)))
	 (small-north-square
	  (floor (/ (mod north square-size) small-square-size))))
    (format t "# ~d ~d ~d ~d~%" east-square north-square
	    small-east-square small-north-square)
    (format nil "~c~c"
	    (square-1 east-square north-square)
	    (square-1 small-east-square small-north-square))))

(defun convert (lat long)
  (pair-to-ref (grid lat long)))