Python和科学计算认证群组  - 讨论区

标题:我的小作业,欢迎吐槽。。。

2011年02月07日 星期一 06:43

话说这个程序捣鼓了一天。本来想用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])

 

 

2011年02月07日 星期一 20:18

基本上看懂了,不过最优解就是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)

 

 

2011年02月09日 星期三 17:07

恩 。。

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

这句话不太懂。。。

什么叫 reshape(-1,1)?

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

威武。。。

2011年02月09日 星期三 18:19

因为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]])

2011年02月09日 星期三 18:25

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

2011年02月19日 星期六 22:56

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

如下红色区域有误,请重新填写。

    你的回复:

    请 登录 后回复。还没有在Zeuux哲思注册吗?现在 注册 !

    Zeuux © 2020

    京ICP备05028076号