1

我制作了一个程序，它从一系列矩阵中计算 R（反射率）的值......但我知道 n 个值（n3 是反射指数）中的一个只是被猜到了。现在我想使用定义的函数来遍历不同的 n3 值以找到最适合我的数据的 R(l)（它们存储在单独的文件中）。我尝试使用 optimization.curve_fit 但老实说我不知道​​在这种情况下如何正确使用它，因为对应于“y”的 R 是由 2 个循环生成的，而不是真正的函数......

from scipy import*
from scipy import linalg
import math
import cmath
from numpy import linalg
import numpy
import matplotlib.pyplot as plt
import scipy.optimize as optimization

infile=open('C:\Users\...\Data.txt', 'r')

n0=1. 
n1=2.31
n2=2. 
n3=1.8
d1=20*10**(-9)
d2=20*10**(-9)

lmin=700*10**(-9)
lmax=2200*10**(-9)
dl=1*10**(-9)
l=lmin

Theta=0

D0=matrix([[math.cos(math.radians(Theta)), math.cos(math.radians(Theta))],[n0, -n0]], dtype=complex)
D1=matrix([[math.cos(math.radians(Theta)), math.cos(math.radians(Theta))],[n1, -n1]], dtype=complex)
D2=matrix([[math.cos(math.radians(Theta)), math.cos(math.radians(Theta))],[n2, -n2]], dtype=complex)



Fi1= 0+0j
Fi2= 0+0j


P1=matrix([[1, 0], [0, 1]], dtype=complex)
P2=matrix([[1, 0], [0, 1]], dtype=complex)

M=matrix([[1,0],[0,1]], dtype=complex)

R=0+0j
r=0+0j

#File reading
x,y = numpy.loadtxt('C:\Users\...\Data.txt',dtype=float).transpose()


DL = []
Res = []

#Funcion
def R(l):
        l=lmin
        D3=matrix([[math.cos(math.radians(Theta)), math.cos(math.radians(Theta))],[n3, -n3]], dtype=complex)
        while l<lmax:
                Fi1= 2.* pi * n1 *d1 * math.cos(math.radians(Theta))/l
                Fi2= 2.* pi * n2 *d2 * math.cos(math.radians(Theta))/l
                f1 = math.e**(1j*Fi1)
                f2 = math.e**(1j*Fi2)
                P1 = matrix([[1/f1, 0], [0, f1]])
                P2 = matrix([[1/f2, 0], [0, f2]])
                T=matrix([[1,0],[0,1]], dtype=complex)
                for i in range(0,9): 
                        T=D1*P1*(D1)**(-1)*D2*P2*(D2)**(-1)*T
                M=D0**(-1)*T*D3
                r=M[1,0]/M[0,0]
                R=abs(M[1,0]/M[0,0])**2
                Res.append(R)
                DL.append(l)
                l=l+dl



        return R

#print optimization.curve_fit(R(l), x, y)

R(l)


plt.plot(x,y)
plt.plot(DL,Res)
plt.show()
4

0 回答 0