[风之影][Lisp大挑战第二季]光路计算

Gu_xl · 发表于 2011-11-15 09:04:29

本帖最后由 Gu_xl 于 2011-11-24 16:28 编辑

参与的人好像不是太多啊！我先发一个编译好的程序，生成1万条光路不到5秒种！希望有兴趣的积极参与！人气足够多时再放源码！请楼主测试下程序是否符合题意！
2011.11.24更新加了反射线效果，并优化了代码，速度更快，生成1万条反射光线和1万条折射光线不到2秒！

cabinsummer · 发表于 2011-11-15 18:42:19

Gu_xl 发表于 2011-11-15 09:04
参与的人好像不是太多啊！我先发一个编译好的程序，生成1万条光路不到5秒种！希望有兴趣的积极参与！人气足 ...

感谢G版的大力支持。程序的效果和速度令人叹为观止。直线圆弧样条线作为分界面都没有问题，即使光线与曲线的交点不止一处时，折射点也会自动找到最近点。不过有点小小的瑕疵。见图片：

用多义线做折射面（图中紫色圆弧组成的一个整体透镜）时，折射光线的方向反了。一开始我还以为G版连反射都做出了，后来一量角度，发现不是反射线。不过，除了对多义线有点不足，其它部分十分完美。期待公布源码。

如果G版有兴趣，可以挑战一个稍微复杂一点得光线追踪，就是给定光线的长度，无论经过多少次的折射，总长度不变，折射面也不一定就是一个了，每个折射面的折射率也不一定是一样了。

再次感谢G版支持。

trjrtj · 发表于 2011-11-20 21:26:54

为了源码，顶起

cabinsummer · 发表于 2011-11-20 21:48:09

trjrtj 发表于 2011-11-20 21:26
为了源码，顶起

再多找一些人，一起顶

trjrtj · 发表于 2011-11-20 22:44:03

再顶一次算否？

cabinsummer · 发表于 2011-11-21 07:07:03

trjrtj 发表于 2011-11-20 22:44
再顶一次算否？

当然算。盖满50楼，gu版共享源码。所以要号召大家一起顶

x_s_s_1 · 发表于 2011-11-21 08:02:27

cabinsummer 发表于 2011-11-21 07:07
当然算。盖满50楼，gu版共享源码。所以要号召大家一起顶

反对没有实质内容的顶帖

cabinsummer · 发表于 2011-11-21 19:47:23

x_s_s_1 发表于 2011-11-21 08:02
反对没有实质内容的顶帖

大侠可以尝试一下挑战，看大侠的程序确实很精彩

qjchen · 发表于 2011-11-23 10:39:49

本帖最后由 qjchen 于 2011-11-24 17:20 编辑

用矢量计算编了一下，本来希望是空间可用的，不过其中有一个矢量绕轴旋转会麻烦些，就只是采用了绕z轴旋转的方法。
没有仔细测试和gu版程序的速度差别，可能用了矢量计算会慢些。没有仔细测试vla-addLine和entmake的差别，10000根都似乎都不太慢。顺手就编了一下反射的，这个比较简单。
由于vla-intersecthwith有精度问题，在处理复杂的spline的时候会有些小问题（刚刚测试了gu版的，也似乎有这个问题），主要是应该有2个交点的被判断为只有一个交点。

;;;;vec plus
(defun q:vec:+(v1 v2)
(mapcar '+ v1 v2)
)
;;;;vec substract
(defun q:vec:-(v1 v2)
(mapcar '- v1 v2)
)
;;;;vec plus constant
(defun q:vec:*c(v a)
(mapcar '(lambda(x) (* x a)) v)
)
;;;;vec dot product
(defun q:vec:dot*(v1 v2)
(apply '+ (mapcar '* v1 v2))
)
;;;;Normalize a vec
(defun q:vec:Norm(v / l) (if (not (zerop (setq l (distance '(0 0 0) v)))) (mapcar '(lambda(x) (/ x l)) v)))
;;;;rotate a vec along z axis
(defun q:vec:rot90z(v) (list (- (cadr v))(car v)(caddr v)))
;;http://en.wikipedia.org/wiki/Snell%27s_law，注意，其中第二条公式网页是错误的
;;n1 是第一个介质; n2是第二个介质
(defun q:refreact(l n n1 n2 / costheta1 costheta2 n12 vrefract)
(setq costheta1 (q:vec:dot* n (q:vec:*c l -1)) n12 (/ n1 n2))
(setq costheta2 (sqrt (- 1 (* (expt n12 2) (- 1 (expt costheta1 2))))))
(if (> costheta1 0)
(setq vrefract (q:vec:+ (q:vec:*c l n12) (q:vec:*c n (- (* n12 costheta1) costheta2))))
(setq vrefract (q:vec:+ (q:vec:*c l n12) (q:vec:*c n (- costheta2 (* -1 n12 costheta1) ))))
)
(setq vreflect (q:vec:+ l (q:vec:*c n (* 2 costheta1))))
(list vrefract vreflect)
)
;; 曲线上取点
(defun q:curvedivide(curve divn / dpar endpar i res startpar)
(setq startPar (vlax-curve-getStartParam curve))
(setq endPar (vlax-curve-getEndParam curve))
(setq i startPar dpar (/ (- endPar startPar) divn))
(repeat divn (setq res (cons (vlax-curve-getPointAtParam curve i) res))(setq i (+ i dpar)))
res
)
(defun q:entmake:line (pt1 pt2 color)
(entmake (list (cons 0 "LINE") ;***
(cons 6 "BYLAYER") (cons 10 pt1) ;***
(cons 11 pt2) ;***
(cons 39 0.0) (cons 62 color) (cons 210 (list 0.0 0.0 1.0))
)
)
)
;;计算两点与curve的交点，非常遗憾的是，IntersectWith存在着计算误差，在高精度时会有问题
(defun q:ent:2pintersect(p1 p2 curve / tempLine pub)
(setq tempLine (vla-addline mSpace (vlax-3d-point p1) (vlax-3d-point p2)))
(setq pub (vlax-invoke tempLine 'IntersectWith curve acExtendOtherEntity))
(vla-delete tempLine)
pub
)
;;;;光线折射和反射的写法
(defun c:test( / a b c curve d f tan temp vec1 vec2 x x1)
;;Vla的写法，速度也比较快，没有仔细测试和entmake的速度
(setq mSpace (vla-get-ModelSpace (vla-get-ActiveDocument (vlax-get-acad-object))))
(if @_d (if (not (setq d @_d @_d (getreal (strcat "\n折射率:<" (rtos @_d 2 2) ">:"))))(setq @_d d)) (setq @_d (getreal "\n折射率:")))
(if @_f (if (not (setq f @_f @_f (getint (strcat "\n根数:<" (rtos @_f 2 0) ">:"))))(setq @_f f)) (setq @_f (getint "\n根数:")))
(setq a (car (entsel "\n折射曲线:")))
(setq curve (vlax-ename->vla-object a))
(setq b (q:curvedivide a f))
(setq c (getpoint "\n 选择基准点:"))
(setq dis (* 0.1 (distance c (car b))))
(foreach x b
(if (= (length (setq temp (q:ent:2pintersect x c curve))) 3)
;;;这句话是迫不得已写的，因为Intersectwith有误差问题，不过即使这样写，在spline有时候也还会出问题
(if (< (distance temp x) 1e-5)
(progn
;(vla-addline mSpace (vlax-3D-point x) (vlax-3D-point c))
(q:entmake:line x c 1)
(setq tan (q:vec:Norm (q:vec:rot90z (vlax-curve-getFirstDeriv a (vlax-curve-getParamAtPoint a x))) ))
(setq vec1 (q:vec:Norm (q:vec:- x c)))
(setq vec2 (q:vec:*c (car (q:refreact vec1 tan 1 d)) dis))
(setq vec3 (q:vec:*c (cadr (q:refreact vec1 tan 1 d)) dis))
(setq x2 (q:vec:+ x vec2) x3 (q:vec:+ x vec3))
;(vla-addline mSpace (vlax-3D-point x) (vlax-3D-point x1))
(q:entmake:line x x2 2)
(q:entmake:line x x3 5)
)
)
)
)
)
(princ "\n By qjchen@gmail.com, 光路计算，命令test")
(princ)

;;;;vec plus 
(defun q:vec:+(v1 v2) 
  (mapcar '+ v1 v2) 
) 
;;;;vec substract 
(defun q:vec:-(v1 v2) 
  (mapcar '- v1 v2) 
) 
;;;;vec plus constant 
(defun q:vec:*c(v a) 
  (mapcar '(lambda(x) (* x a)) v) 
) 
;;;;vec dot product 
(defun q:vec:dot*(v1 v2) 
  (apply '+ (mapcar '* v1 v2)) 
) 
 
;;;;Normalize a vec 
(defun q:vec:Norm(v / l) (if (not (zerop (setq l (distance '(0 0 0) v)))) (mapcar '(lambda(x) (/ x l)) v))) 
 
;;;;rotate a vec along z axis 
(defun q:vec:rot90z(v) (list (- (cadr v))(car v)(caddr v))) 
 
;;http://en.wikipedia.org/wiki/Snell%27s_law，注意，其中第二条公式网页是错误的 
;;n1 是第一个介质; n2是第二个介质 
(defun q:refreact(l n n1 n2 / costheta1 costheta2 n12 vrefract) 
  (setq costheta1 (q:vec:dot* n (q:vec:*c l -1)) n12 (/ n1 n2)) 
  (setq costheta2 (sqrt (- 1 (* (expt n12 2) (- 1 (expt costheta1 2)))))) 
  (if (> costheta1 0) 
  (setq vrefract (q:vec:+ (q:vec:*c l n12)  (q:vec:*c n (- (* n12 costheta1) costheta2)))) 
  (setq vrefract (q:vec:+ (q:vec:*c l n12)  (q:vec:*c n (- costheta2 (* -1 n12 costheta1) )))) 
  ) 
  (setq vreflect (q:vec:+ l (q:vec:*c n (* 2 costheta1)))) 
   
  (list vrefract vreflect) 
) 
 
;; 曲线上取点 
(defun q:curvedivide(curve divn / dpar endpar i res startpar) 
 (setq startPar (vlax-curve-getStartParam curve)) 
 (setq endPar (vlax-curve-getEndParam curve)) 
 (setq i startPar dpar (/ (- endPar startPar) divn)) 
 (repeat divn (setq res (cons (vlax-curve-getPointAtParam curve i) res))(setq i (+ i dpar))) 
 res 
) 
 
(defun q:entmake:line (pt1 pt2 color) 
  (entmake (list (cons 0 "LINE") ;*** 
                 (cons 6 "BYLAYER") (cons 10 pt1) ;*** 
                 (cons 11 pt2) ;*** 
                 (cons 39 0.0) (cons 62 color) (cons 210 (list 0.0 0.0 1.0)) 
           ) 
  ) 
) 
 
;;计算两点与curve的交点，非常遗憾的是，IntersectWith存在着计算误差，在高精度时会有问题 
(defun q:ent:2pintersect(p1 p2 curve / tempLine pub) 
 (setq tempLine (vla-addline mSpace (vlax-3d-point p1) (vlax-3d-point p2))) 
 (setq pub (vlax-invoke tempLine 'IntersectWith curve acExtendOtherEntity)) 
 (vla-delete tempLine) 
 pub 
) 
 
;;;;光线折射和反射的写法 
(defun c:test( / a b c curve d f tan temp vec1 vec2 x x1) 
  ;;Vla的写法，速度也比较快，没有仔细测试和entmake的速度 
  (setq mSpace (vla-get-ModelSpace (vla-get-ActiveDocument (vlax-get-acad-object)))) 
  (if @_d (if (not (setq d @_d @_d (getreal (strcat "\n折射率:<" (rtos @_d 2 2) ">:"))))(setq @_d d)) (setq @_d (getreal "\n折射率:"))) 
  (if @_f (if (not (setq f @_f @_f (getint (strcat "\n根数:<" (rtos @_f 2 0) ">:"))))(setq @_f f)) (setq @_f (getint "\n根数:"))) 
  (setq a (car (entsel "\n折射曲线:"))) 
  (setq curve (vlax-ename->vla-object a)) 
  (setq b (q:curvedivide a f)) 
  (setq c (getpoint "\n 选择基准点:")) 
  (setq dis (* 0.1 (distance c (car b)))) 
  (foreach x b 
    (if (= (length (setq temp (q:ent:2pintersect x c curve))) 3) 
      ;;;这句话是迫不得已写的，因为Intersectwith有误差问题，不过即使这样写，在spline有时候也还会出问题 
      (if (< (distance temp x) 1e-5) 
        (progn 
         ;(vla-addline mSpace (vlax-3D-point x) (vlax-3D-point c)) 
         (q:entmake:line x c 1) 
         (setq tan (q:vec:Norm (q:vec:rot90z (vlax-curve-getFirstDeriv a  (vlax-curve-getParamAtPoint a x))) )) 
         (setq vec1 (q:vec:Norm (q:vec:- x c))) 
         (setq vec2 (q:vec:*c (car (q:refreact vec1 tan 1 d)) dis)) 
         (setq vec3 (q:vec:*c (cadr (q:refreact vec1 tan 1 d)) dis)) 
         (setq x2 (q:vec:+ x vec2) x3 (q:vec:+ x vec3)) 
         ;(vla-addline mSpace (vlax-3D-point x) (vlax-3D-point x1)) 
         (q:entmake:line x x2 2) 
         (q:entmake:line x x3 5) 
         ) 
      ) 
    ) 
  ) 
) 
 
(princ "\n By qjchen@gmail.com, 光路计算，命令test") 
(princ)

zyhandw · 发表于 2011-11-23 11:18:06

为了G版的源码，顶啊！！！

账号		自动登录	找回密码
密码			注册

[【风之影】] [风之影][Lisp大挑战第二季]光路计算

本帖子中包含更多资源

评分

本帖子中包含更多资源

本帖子中包含更多资源

评分