2016-11-10 7 views
1

50層(d1、d2)の厚さに応じて、多層システムの反射係数を計算します。任意の2つの数字(d1、d2)を取ると、正しく動作します。しかし、私はワイヤフレームプロットを取得する必要があります.d1、d2はある範囲で意味を持ち、エラーが発生します。「ValueError:入力は13行目の四角形の配列でなければなりません。ワイヤフレームプロット

from math import pi 
import numpy as np 
import matplotlib.pyplot as plt 

def R(n1, n2, d1, d2, lamda): 
    phy1 = (-2*pi*n1*d1/lamda) 
    phy2 = (-2*pi*n2*d2/lamda) 
    DPD1 = 0.5*np.array([[2*np.cos(phy1), 2j*np.sin(phy1)/n1], [n1*2j*np.sin(phy1), 2*np.cos(phy1) ]]) 
    DPD2 = 0.5*np.array([[2*np.cos(phy2), 2j*np.sin(phy2)/n2], [n2*2j*np.sin(phy2), 2*np. cos(phy2) ]]) 
    D0 = 0.5 * np.array([[1, 1], [1, -1]]) 
    DS = np.array([[1, 1], [n1, -n1]]) 
    DPD = np.dot(DPD1, DPD2) 
    DPD = np.linalg.matrix_power(DPD, 50) 
    M = np.dot(D0, DPD) 
    M = np.dot(M, DS) 
    return(abs(M[1,0]/M[0,0])**2) 

x = np.arange(0, 10, 1) 
y = np.arange(0, 10, 1) 
X, Y = np.meshgrid(x, y) 
Z = R(0.99910053+0.00183184j, 0.92373900+0.00644652j, X, Y, 13.5) 
fig = plt.figure() 
ax = fig.gca(projection='3d') 
surf = ax.plot_wireframe(X, Y, Z, antialiased=True) 

答えて

0

私はforループなしで行う方法がわかりません。ここに私の解決策があります:

from math import pi 
import numpy as np 
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 

def R(n1, n2, d1, d2, lamda): 
    phy1 = (-2*pi*n1*d1/lamda) 
    phy2 = (-2*pi*n2*d2/lamda) 
    DPD1 = 0.5*np.array([[2*np.cos(phy1), 2j*np.sin(phy1)/n1], [n1*2j*np.sin(phy1), 2*np.cos(phy1) ]]) 
    DPD2 = 0.5*np.array([[2*np.cos(phy2), 2j*np.sin(phy2)/n2], [n2*2j*np.sin(phy2), 2*np. cos(phy2) ]]) 
    DPD = np.dot(DPD1, DPD2) 
    DPD = np.linalg.matrix_power(DPD, 50) 
    D0 = 0.5 * np.array([[1, 1], [1, -1]]) 
    DS = np.array([[1, 1], [n1, -n1]]) 
    M = np.dot(D0, DPD) 
    M = np.dot(M, DS) 
    return(abs(M[1,0]/M[0,0])**2) 

x = np.arange(0, 10, 1) 
y = np.arange(0, 10, 1) 
X, Y = np.meshgrid(x, y) 
Z = np.zeros(X.shape).ravel() 
for i, (x, y) in enumerate(zip(X.ravel(), Y.ravel())): 
    Z[i] = R(0.99910053+0.00183184j, 0.92373900+0.00644652j, x, y, 13.5) 
fig = plt.figure() 
ax = fig.gca(projection='3d') 
surf = ax.plot_wireframe(X, Y, Z.reshape(X.shape), antialiased=True) 
plt.show()