话说这个程序捣鼓了一天。本来想用class的，但是一直出错。。后来某位大哥提点了我下，可是我已经把这个写的差不多了。

一开始生成一束光，类似于高斯光束。x——y面是高斯面，原点是焦点。然后用每道光和平面的交点到这些点的中心的平方和来评价这个面是否为最好的面。 r是面到原点的距离，dx，dy，dz是面的法向。

感谢张若愚大哥的提点。。

嗨皮兔爷！

# -*- coding: utf-8 -*-

"""

Created on Wed Feb 02 13:35:50 2011

@author: Praktikant

"""

from scipy.optimize import fmin,fmin_powell

import numpy as np

n = 200

def normalize(x):

     norm = np.sqrt(np.dot(x[:,0],x[:,0]))

     return x/norm

def twistbeam(n,r):

     z1 = np.ones(n)

     z2 = -np.ones(n)

     phi1 = 2*np.pi/n*(np.arange(n)+1)

     phi2 = phi1+np.pi/4

     x1 = r*np.cos(phi1)

     y1 = r*np.sin(phi1)

     P1 = np.array([x1,y1,z1])

     x2 = r*np.cos(phi2)

     y2 = r*np.sin(phi2)

     P2 = np.array([x2,y2,z2])

     V = P1-P2

     return P1,V

Beam_p,Beam_v = twistbeam(n,5)

print 'Beam\n',Beam_p,'\n',Beam_v

def note(P1,V,P2,N): #line cross plane: P1,V of line, P2,N of plane

     t = np.dot((P2-P1),N)/np.dot(V,N)

     note = P1 + t*V

     return note

def rms(P): #P: a set of point

     Pm = np.array([np.mean(P[0]),np.mean(P[1]),np.mean(P[2])])

     P1 = P.T-Pm

     P1 = P1.T

     rms = np.sqrt(np.sum(P1**2)/n)

     return rms

def evaluation(Plane):

     points = np.zeros([3,n])

     N = np.array([Plane[0],Plane[1],Plane[2]])

     N = N/np.dot(N,N)

     P0 = Plane[3]*N

     for i in range(n):

          P1,V = Beam_p[:,i],Beam_v[:,i]

          points[:,i]=note(P1,V,P0,N)

     rr = rms(points)

     return rr

r=0

dx = 0

dy = 0

dz = 1

test1 = evaluation([r,dx,dy,dz])

print test1

bestfit = fmin_powell(evaluation,(0.,0.,2.,1.),ftol = 0.000000001)

bestfit1 = fmin(evaluation,(0.,0.,2.,1.),ftol = 0.000000001)

print '\n\n',bestfit,'\nresult:',evaluation(bestfit)

print '\n\n',bestfit1,'\nresult1:',evaluation(bestfit1)

print evaluation([0,0,1,0.001])

基本上看懂了，不过最优解就是X-Y平面吧。是不是还有别的光线分布情况使得最优解在别的位置？

另外有几个地方可以改进：

N = np.array([Plane[0],Plane[1],Plane[2]])

改为

N = Plane[:3]

Pm = np.array([np.mean(P[0]),np.mean(P[1]),np.mean(P[2])])

P1 = P.T-Pm

P1 = P1.T

可改为

P1 = P - np.mean(P,axis=1).reshape(-1,1)

恩 。。

P1 = P - np.mean(P,axis=1).reshape(-1,1)

这句话不太懂。。。

什么叫 reshape（-1，1）?

楼上是做什么的？很神的样子。。

威武。。。

因为P1的shape为(3,200), np.mean(P, axis=1)所得到的数组的shape为(3,)，为了让P1中的每行减去其中的一个数，需要将数组的shape改为(3,1)。reshape完成shape变换，-1表示根据数组原来的长度计算，实际上就是3。

In [2]: np.array([1,2,3]).reshape(-1,1)

Out[2]:

array([[1], [2], [3]])

非常不错。。。谢谢了。。

张大哥,真的很认真阿,能耐心看一个程序,耐心可嘉.