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

cabinsummer

　　上一季大挑战中的圆周率吸引了不少人的关注，但是参与者很少。到目前为止，也只有高飞鸟的程序最高效。可能是由于太偏理论了，都是数值计算，而这并不是lisp的强项，所以很多人望而却步。
　　这一季的lisp大挑战是图形方面的，图形可是CAD的强项，应该能吸引不少人直接参与吧？但愿这次更火爆。
　　好了，唠叨这么久，风该把题目拿出来吧。
　　呵呵，本次挑战的题目是光路计算。图片中有两条光线通过透镜形成的折射光路，红色光和绿色光。
　　这次要求计算任意折射曲面（一个自定义曲面，图片中是三个曲面难度大了点）的折射光路。曲面是空气和某种玻璃的分界面，假定折射率是1.33吧（这个折射率应该能在程序中设定），求空气中穿过某点（任意点）的光线（任意光线）穿过空气-玻璃分界面后形成的光路，只计算方向。最终输出结果画出光路。冠军将是生成1000条光路速度最快的那位。

　　折射率的公式大家都知道吧？n=sinθ/sinφ。θ是光线的入射角，φ是光线的折射角。这个不再普及了，不清楚的就Google一下。
　　祝大家好运！

tigcat

baitang36 发表于 2020-12-4 08:16
已经够50楼了，gu版应该公布源码了。如果他不公布，你找我私聊，我发给你源码。

好的代码就像美酒一样，多少年过去，仍芬香四溢。等一个月看看顾版有空上来不，不行大侠发我邮箱，tigcat@126.com，先谢谢了。

fanpaulm

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

G版的东西，确实不一般！

baitang36

tigcat 发表于 2020-12-1 21:09
不知道gu版说的话还有效没？现在来帮忙顶一把。

已经够50楼了，gu版应该公布源码了。如果他不公布，你找我私聊，我发给你源码。

cabinsummer

首先声明一点，风不是学光学的，从事的专业和光也没有半点关系。

LLXXZZ

貌似,我不会.有点难,数学不好.

cabinsummer

提供思路
求入射光与曲面的交点，做曲面在交点处的法线。这样得出入射角，再根据折射定律计算折射角。有了交点和折射角，就能画出折射光线。

yjr111

我估计会无疾而终，因为。。。物理的书找不到了

cabinsummer

这是几何光学里的，根本不需要太多的物理知识，再说，谷歌百度维基都是很好的学习平台

prolee2008

写的比较繁琐，效率也不高，大家不要笑话。先画100条，还只是一个反射面的.

1. (defun c:zheshe(/ po_list n qiexian_list);;__LiXuDong
   
2.   (vl-load-com)
   
3.   (setvar "OSMODE" 0)
   
4.   (setq n 1.33)
   
5.   (setq curve_n1 (car (entsel)))
   
6.   (setq curve_n1 (vlax-ename->vla-object curve_n1))
   
7.   (setq startSpline (vlax-curve-getStartParam curve_n1))
   
8.   (setq endSpline (vlax-curve-getEndParam curve_n1))
   
9.   (setq po_list '())
   
10.   (setq qiexian_list '())
    
11.   (setq i 1)
    
12.   (repeat 100
    
13.     (setq po_curve (vlax-curve-getpointatparam curve_n1 (* i (/ (- endSpline startSpline) 100.))))
    
14. ;    (print po_curve)
    
15.     (setq po_list (cons  po_curve po_list))
    
16.     (setq i (1+ i))
    
17.    
    
18.     )
    
19.     (setq i 1)
    
20.     (repeat 100
    
21.     (setq po_qiexiang_curve (vlax-curve-getfirstderiv curve_n1 (* i (/ (- endSpline startSpline) 100.))))
    
22.     (setq po_qiexiang_curve (mapcar '(lambda(x) (* x 1000)) po_qiexiang_curve))
    
23.     (setq po_qiexiang_curve (mapcar '+ (nth (- 100 i) po_list) po_qiexiang_curve))
    
24.     (setq qiexian_list (cons  po_qiexiang_curve qiexian_list))
    
25.     (setq i (1+ i))
    
26.       )
    
27. ;(print qiexian_list)
    
28.     (setq i 0)
    
29.     (setq in_po (getpoint))
    
30.     (repeat 100
    
31.       (command "line" in_po (nth i po_list) "")
    
32.       (setq i (1+ i))
    
33.       )
    
34.     (setq i 0)
    
35.     (repeat 100
    
36.       (setq angle1 (jaj in_po (nth i po_list) (nth i qiexian_list)))
    
37.       (setq shurejiao (abs (- 90 angle1)))
    
38.       (setq chushejiao_sin (/ (sin (* PI (/ shurejiao 180.))) n))
    
39.       (setq chushejiao_cos (sqrt (- 1 (* chushejiao_sin chushejiao_sin))))
    
40.       (setq chushejiao (getc chushejiao_cos))
    
41.       (command "line" (nth i po_list) (nth i qiexian_list) "")
    
42.       (setq obj (entlast))
    
43.       (command "rotate" obj "" (nth i po_list) "0.5" "")
    
44.       (setq obj_duandian(assoc 11 (entget obj)))
    
45.       (setq obj_duandian (list (nth 1 obj_duandian) (nth 2 obj_duandian) (nth 3 obj_duandian)))
    
46.       (setq angle2 (jaj in_po (nth i po_list) obj_duandian))
    
47.       (command "undo" "1" "")
    
48.       (if(< angle2 angle1)
    
49.         (if(< angle1 90.)
    
50.           (command "rotate" obj "" (nth i po_list) (- 270 chushejiao)"")
    
51.           (command "rotate" obj "" (nth i po_list) (+ 270 chushejiao)"")
    
52.          
    
53.           )
    
54.         (if(< angle1 90.)
    
55.           (command "rotate" obj "" (nth i po_list) (+ 90 chushejiao)"")
    
56.           (command "rotate" obj "" (nth i po_list) (- 90 chushejiao)"")
    
57.           )
    
58.         
    
59.         
    
60.         )
    
61.       (setq i (1+ i))
    
62.       
    
63.       )
    
64.   )
    
65.   
    
66.   
    
67.   
    
68.   (defun jaj (pt1 pt2 pt3)
    
69. ;摘自modernyoung的博客
    
70.   (setq a (distance pt1 pt3)
    
71.          b (distance pt1 pt2)
    
72.          c (distance pt2 pt3)
    
73.   )
    
74.   (setq cos_a (/ (- (+ (* b b) (* c c)) (* a a))
    
75.    (* b c 2)
    
76.    ;;求余弦
    
77.        )
    
78.   )
    
79.   (if (= cos_a 0)
    
80.       (setq jad 90)
    
81. 
    
82. 
    
83.     (progn (setq jad1 (atan (/ (sqrt (- 1 (* cos_a cos_a))
    
84.          ;;转换成正切
    
85.           )
    
86.           cos_a
    
87.        )
    
88.         )
    
89.     )
    
90.     (if (<= jad1 0)
    
91.       ;;当值是负值时，转换角度
    
92. 
    
93.       (setq jad (* (/ (+ pi jad1) pi)
    
94.      180
    
95.          )
    
96.       )
    
97.       (setq jad (* (/ jad1 pi)
    
98.      180
    
99.          )
    
100.       )
     
101.     )
     
102.     )
     
103.   )
     
104. 
     
105.   jad
     
106. )
     
107. 
     
108. (defun GetC (cosc / C)
     
109. 
     
110.   (setq C (atan (sqrt (1- (/ 1 (* cosc cosc))))))
     
111. 
     
112.   (setq c (* 180.0 (/ C PI)))
     
113.   (if(> c 90)
     
114.     (setq c (- C 90)))
     
115. 
     
116.   C
     
117. 
     
118. )
     

(defun c:zheshe(/ po_list n qiexian_list);;__LiXuDong<br />   (vl-load-com)<br />   (setvar "OSMODE" 0)<br />   (setq n 1.33)<br />   (setq curve_n1 (car (entsel)))<br />   (setq curve_n1 (vlax-ename->vla-object curve_n1))<br />   (setq startSpline (vlax-curve-getStartParam curve_n1))<br />   (setq endSpline (vlax-curve-getEndParam curve_n1))<br />   (setq po_list '())<br />   (setq qiexian_list '())<br />   (setq i 1)<br />   (repeat 100<br />     (setq po_curve (vlax-curve-getpointatparam curve_n1 (* i (/ (- endSpline startSpline) 100.))))<br /> ;    (print po_curve)<br />     (setq po_list (cons  po_curve po_list))<br />     (setq i (1+ i))<br />     <br />     )<br />     (setq i 1)<br />     (repeat 100<br />     (setq po_qiexiang_curve (vlax-curve-getfirstderiv curve_n1 (* i (/ (- endSpline startSpline) 100.))))<br />     (setq po_qiexiang_curve (mapcar '(lambda(x) (* x 1000)) po_qiexiang_curve))<br />     (setq po_qiexiang_curve (mapcar '+ (nth (- 100 i) po_list) po_qiexiang_curve))<br />     (setq qiexian_list (cons  po_qiexiang_curve qiexian_list))<br />     (setq i (1+ i))<br />       )<br /> ;(print qiexian_list)<br />     (setq i 0)<br />     (setq in_po (getpoint))<br />     (repeat 100<br />       (command "line" in_po (nth i po_list) "")<br />       (setq i (1+ i))<br />       )<br />     (setq i 0)<br />     (repeat 100<br />       (setq angle1 (jaj in_po (nth i po_list) (nth i qiexian_list)))<br />       (setq shurejiao (abs (- 90 angle1)))<br />       (setq chushejiao_sin (/ (sin (* PI (/ shurejiao 180.))) n))<br />       (setq chushejiao_cos (sqrt (- 1 (* chushejiao_sin chushejiao_sin))))<br />       (setq chushejiao (getc chushejiao_cos))<br />       (command "line" (nth i po_list) (nth i qiexian_list) "")<br />       (setq obj (entlast))<br />       (command "rotate" obj "" (nth i po_list) "0.5" "")<br />       (setq obj_duandian(assoc 11 (entget obj)))<br />       (setq obj_duandian (list (nth 1 obj_duandian) (nth 2 obj_duandian) (nth 3 obj_duandian)))<br />       (setq angle2 (jaj in_po (nth i po_list) obj_duandian))<br />       (command "undo" "1" "")<br />       (if(< angle2 angle1)<br />         (if(< angle1 90.) <br />           (command "rotate" obj "" (nth i po_list) (- 270 chushejiao)"")<br />           (command "rotate" obj "" (nth i po_list) (+ 270 chushejiao)"")<br />           <br />           )<br />         (if(< angle1 90.)<br />           (command "rotate" obj "" (nth i po_list) (+ 90 chushejiao)"")<br />           (command "rotate" obj "" (nth i po_list) (- 90 chushejiao)"")<br />           )<br />         <br />         <br />         )<br />       (setq i (1+ i))<br />       <br />       )<br />   )<br />   <br />   <br />   <br />   (defun jaj (pt1 pt2 pt3)<br /> ;摘自modernyoung的博客 <br />   (setq a (distance pt1 pt3)<br />          b (distance pt1 pt2)<br />          c (distance pt2 pt3)<br />   )<br />   (setq cos_a (/ (- (+ (* b b) (* c c)) (* a a))<br />    (* b c 2)<br />    ;;求余弦<br />        )<br />   )<br />   (if (= cos_a 0)<br />       (setq jad 90)<br /> <br /> <br />     (progn (setq jad1 (atan (/ (sqrt (- 1 (* cos_a cos_a))<br />          ;;转换成正切<br />           )<br />           cos_a<br />        )<br />         )<br />     )<br />     (if (<= jad1 0)<br />       ;;当值是负值时，转换角度<br /> <br />       (setq jad (* (/ (+ pi jad1) pi)<br />      180<br />          )<br />       )<br />       (setq jad (* (/ jad1 pi)<br />      180<br />          )<br />       )<br />     )<br />     )<br />   )<br /> <br />   jad<br /> )<br /> <br /> (defun GetC (cosc / C)<br /> <br />   (setq C (atan (sqrt (1- (/ 1 (* cosc cosc))))))<br /> <br />   (setq c (* 180.0 (/ C PI)))<br />   (if(> c 90) <br />     (setq c (- C 90)))<br /> <br />   C<br /> <br /> )<br />

prolee2008

任意曲面。

cabinsummer

prolee2008 发表于 2011-11-7 21:52
写的比较繁琐，效率也不高，大家不要笑话。先画100条，还只是一个反射面的.

先对prolee2008的程序表示赞赏，不过程序有点小小的瑕疵。在上面的图片中，曲线为折射面，入射点位于左下方，但是程序运行结果是穿过曲线找第二交点折射，这明显与我们的常识不符。希望能加以改进。不过，prolee2008的程序目前暂列挑战榜首位。