Autolisp高斯反算，CGCS2000大地坐标系平面坐标转换为经纬度坐标

CAD制图中，经常要用到根据图上某点坐标标注经纬度的功能，之前在网上找了一个西安80坐标系下标注的，现在用2000的比较多，网上没有找到合适的，所以自己照着三调的培训资料写了一下，三调的资料内有高斯投影的计算公式以及图幅分幅和随球面积等计算内容。网上有很多.net下写的2000坐标系下高斯反算的代码，但是我没看懂，没理解需要迭代循环的一部分，如果有懂的，可以给科普一下。谢谢。
五个常数自行翻阅三调培训资料，我不确定是不是涉密。

(defun xytobl2000 (y x cm / a bb f e2 e12 p k0 k1 k2 k3	k4 E Bf	tt n2 C
		   V N B L)
					;椭球参数及常数
  (setq	a   6378137.0
	bb  6356752.31414036
	f   0.003352810681182
	e2  0.006694380022901
	e12 0.006739496775479
	p   206264.806247097
	k0  0.00000015***
	k1  0.005052****
	k2  0.00002***
	k3  0.00000024*****
	k4  0.0000*****
  )
  (setq	E  (* k0 x)
	Bf (+ E
	      (* (cos E)
		 (- (+ (- (* k1 (sin E)) (* k2 (expt (sin E) 3)))
		       (* k3 (expt (sin E) 5))
		    )
		    (* k4 (expt (sin E) 7))
		 )
	      )
	   )
	tt (/ (sin Bf) (cos Bf))
	n2 (* e12 (expt (cos Bf) 2))
	C  (/ (expt a 2) bb)
	V  (sqrt (+ 1 n2))
	N  (/ C V)
  )
  (setq	B (-
	    (+
	      (- Bf
		 (* 0.5
		    (* (expt V 2) tt (expt (/ y N) 2))
		 )
	      )
	      (* (/ 1 24.0)
		 (- (+ 5 (* 3 (expt tt 2)) n2) (* 9 n2 (expt tt 2)))
		 (* (expt V 2) tt)
		 (expt (/ y N) 4)
	      )
	    )
	    (*
	      (/ 1 720.0)
	      (+
		61
		(* 90 (expt tt 2))
		(* 45 (expt tt 4))
	      )
	      (* (expt V 2) tt)
	      (expt (/ y N) 6)
	    )
	  )
  )					;endB
  (setq	L (+
	    (-
	      (*
		(/ 1 (cos Bf))
		(/ y N)
	      )
	      (*
		(/ 1 6.0)
		(+
		  1
		  (* 2 (expt tt 2))
		  n2
		)
		(/ 1 (cos Bf))
		(expt (/ y N) 3)
	      )
	    )
	    (*
	      (/ 1 120.0)
	      (+
		5
		(* 28 (expt tt 2))
		(* 24 (expt tt 4))
		(* 6 n2)
		(* 6 n2 (expt tt 2))
	      )
	      (/ 1 (cos Bf))
	      (expt (/ y N) 5)
	    )
	    cm
	  )

  )					;endL

  (list (dfm L) (dfm B))
)					;end defun

实际上没有什么技术含量，就是照书上的公式写出来就行，有了公式和简化后的常量，EXCEL也可以代入计算。