Autolisp - UTM para LatLon e viceversa

Bom, hoje vou voltar um pouco às origens do blog!!! Um pouco de autolisp pra relembrar os velhos tempos de programação em POG!!! A lisp abaixo na verdade são algumas subrotinas para conversão de coordenadas geográficas em UTM e viceversa. É bem fácil de usar se você souber o que é UTM e coordenada geográfica e está familiarizado com Georeferenciamento. Muitas pessoas me perguntam se eu tenho uma rotina pra converter e.... Bem, tenho!! Está aí!!

---

;| Conversão de UTM para GEOGRAFICA
baseado em http://recursos.gabrielortiz.com/index.asp?Info=058a
metodo: Coticchia-Surace
ay:     coorenada do semi eixo maior (m)
bx:     coorenada do semi eixo menor (m)
pt:     (cord_X coord_Y coord_Z)
fuso:   fuso, inteiro
hmsf:   hemisferio, "N" para norte e "S" para sul
se for conhecido f, temos: bx=(1-f)*ay
sad69 ->          ay = 6.378.160,000m e f = 1/298,25
corrego alegre -> ay = 6.378.388,000m e f = 1/297,00
(utm2geo (GETPOINT) 6378160.0  298.25 22 "S") => 25º25'51"15216  49º17'02"51881
|;
(defun utm2geo (pt ay f fuso hmsf / tmp e el el² c alpha lo phil nu A1 A2 J2 J4 re
                J6 beta gama Bo zeta xi eta sinhxi dl tau lat lon x y x1 y1 fe bx)
  (setq ay     (float ay)
        bx     (* ay (- 1 (/ 1.0 f)))
        x      (car pt)
        y      (cadr pt)
        re     6366197.724 ;raio da terra
        fe     0.9996      ;fator de escala
        tmp    (sqrt (- (expt ay 2) (expt bx 2)))
        e      (/ tmp ay) ;excentricidade
        el     (/ tmp bx) ;2ª excentricidade
        el²    (expt el 2)
        c      (/ (expt ay 2) bx);raio polar de curvatura
        x1     (- x 500000.0)
        y1     (if (or (= hmsf 'N) (= hmsf "N")) y (- y 10000000.0))
        phil   (/ y1 (* re fe))
        lo     (- (* 6 fuso) 183)
        nu     (/ (* c fe) (sqrt (1+ (* el² (expt (cos phil) 2)))))
        a      (/ x1 nu)
        A1     (sin (* 2 phil))
        A2     (* A1 (expt (cos phil) 2))
        J2     (+ phil (/ A1 2.0))
        J4     (/ (+ (* J2 3.0) A2) 4.0)
        J6     (/ (+ (* 5.0 J4) (* A2 (expt (cos phil) 2))) 3.0)
        alpha  (/ (* 3.0 el²) 4.0)
        beta   (* (/ 5.0 3.0) (expt alpha 2))
        gama   (* (/ 35.0 27.0) (expt alpha 3))
        Bo     (* fe c (+ phil (* (- alpha) J2) (* beta J4) (* (- gama) J6)))
        b      (/ (- y1 Bo) nu)
        zeta   (* (/ (* el² (expt a 2)) 2.0) (expt (cos phil) 2))
        xi     (* a (- 1.0 (/ zeta 3.0)))
        eta    (+ (* b (- 1.0 zeta)) phil)
        sinhxi (/ (- (exp xi) (exp (- xi))) 2.0)
        dl     (atan (/ sinhxi (cos eta)))
        tau    (atan (* (cos dl) (tan eta)))
        lon    (+ (* (/ 180.0 pi) dl) lo)
        lat    (* (/ 180.0 pi)
                  (+ phil (* (+ 1.0
                                (* el² (expt (cos phil) 2.0))
                                (* (/ -3.0 2.0) el² (sin phil) (cos phil) (- tau phil)))
                             (- tau phil)))))
  (if (caddr pt)
    (list lon lat (caddr pt))
    (list lon lat 0.0)))


;pt -> long lat
;a  -> semi eixo maior
;f  -> achatamento
(defun geo2utm (pt a f / b e el el² c lamb fi fuso lo deltal Am eps n v
                S A1 A2 J2 J4 J6 alfa beta gama bo)
  (setq a      (float a)
        b      (- a (/ a f))
        el     (/ (sqrt (- (expt a 2) (expt b 2))) b)
        el²    (expt el 2)
        c      (/ (expt a 2) b)
        fuso   (fix (+ (/ (car pt) 6.0) 31))
        lamb   (/ (* (car pt) pi) 180.0)
        fi     (/ (* (cadr pt) pi) 180.0)
        lo     (- (* fuso 6) 183) ;meridiano central
        deltal (- lamb (/ (* lo pi) 180.0))
        Am     (* (cos fi) (sin deltal))
        eps    (* 0.5 (log (/ (+ 1 Am) (- 1 Am))))
        n      (- (atan (/ (tan fi) (cos deltal))) fi)
        v      (/ (* c 0.9996) (sqrt (+ 1 (* el² (expt (cos fi) 2)))))
        S      (/ (expt (* el eps (cos fi)) 2) 2.0)
        A1     (sin (* 2.0 fi))
        A2     (* A1 (expt (cos fi) 2.0))
        J2     (+ fi (/ A1 2.0))
        J4     (/ (+ (* 3.0 J2) A2) 4.0)
        J6     (/ (+ (* 5 J4) (* A2 (expt (cos fi) 2))) 3.0)
        alfa   (/ (* 3.0 el²) 4.0)
        beta   (* (/ 5.0 3.0) (expt alfa 2))
        gama   (* (/ 35.0 27.0) (expt alfa 3))
        bo     (* 0.9996 c (+ fi (* (- alfa) J2) (* beta J4) (* (- gama) J6))))
  (list  (+ 500000.0 (* eps v (1+ (/ S 3.0)))) ;x
         (+ bo (* n v (1+ S)) (if (< lat 0.0) 10000000.0 0.0));y
         (caddr pt)
         ))

(defun LLA_wgs84->sad69 (pt)
   (geo2geo pt 6378137.0 298.257223563 6378160.0 298.25 66.87 -4.37  38.52))

(defun LLA_sad69->wgs84 (pt)
   (geo2geo pt 6378160.0 298.25  6378137.0 298.257223563 -66.87  4.37 -38.52))


      
(defun geo2geo (pt ;long_from lat_from h_from ;coordenadas geodesicas de origem
                a_from f_from             ;parametros geodesicos de origem
                a_to   f_to               ;parametros geodesicos de destino
                dx dy dz                  ;translação origem->destino
                / lat1 long1 f1 f2 a1 a2 e²1 e²2 N1 N2 Xw Yw Zw b2 p teta fi lamb hb ep²)
  (setq lat1   (/ (* (cadr pt) pi) 180.0)
        long1  (/ (* (car pt) pi) 180.0)
        f1     (/ 1.0 f_from)        ;wgs
        a1     a_from                ;wgs
        e²1    (* f1 (- 2.0 f1))
        N1     (/ a1 (sqrt (- 1.0 (* e²1 (expt (sin lat1) 2.0)))))
        ;coord carteziana no sistema de origem:
        Xw     (* (+ N1 (caddr pt)) (cos lat1) (cos long1))
        Yw     (* (+ N1 (caddr pt)) (cos lat1) (sin long1))
        Zw     (* (+ (* N1 (- 1.0 e²1)) (caddr pt)) (sin lat1))
        ;coord carteziana do ponto no novo sistema:
        Xb     (+ Xw dx)
        Yb     (+ Yw dy)
        Zb     (+ Zw dz)
        ;converter carteziana para geodesica:
        f2     (/ 1.0 f_to)
        a2     a_to
        e²2    (* f2 (- 2.0 f2))
        b2     (* a2 (- 1.0 f2))
        p      (sqrt (+ (expt Xb 2) (expt Yb 2)))
        teta   (atan (/ (* Zb a2) (* p b2)))
        ep²    (/ (- (expt a2 2.0) (expt b2 2.0)) (expt b2 2.0))
        fi     (atan (/ (+ Zb (* ep² b2 (expt (sin teta) 3.0))) (- p (* e²2 a2 (expt (cos teta) 3.0)))))
        lamb   (atan (/ Yb Xb))
        N2     (/ a2 (sqrt (- 1.0 (* e²2 (expt (sin fi) 2.0)))))
        hb     (- (/ p (cos fi)) N2))
  (list (/ (* 180.0 lamb) pi);long
        (/ (* 180.0 fi) pi)  ;lat
        hb)                  ;altitude
  )


Link(s) da(s) subrotina(s) usada(s): tan

---

A dificuldade nem está nos cálculos em si, mas na sintaxe das fórmulas, não acham??