diff --git a/wptherml/lightlib.py b/wptherml/lightlib.py index 9238af8..2a4e877 100644 --- a/wptherml/lightlib.py +++ b/wptherml/lightlib.py @@ -37,19 +37,6 @@ def normalized_power(lam, TE, BB): return numerator/denominator -def Lum_efficiency_prime(lam, TE, TE_prime): - upper = np.amax(lam) - ph = datalib.PhLum(lam) - num = ph*TE - num_prime = ph*TE_prime - numerator = numlib.Integrate(num, lam, 0, upper ) - numerator_prime = numlib.Integrate(num_prime, lam, 0, upper) - - denominator = numlib.Integrate(TE, lam, 0, upper) - denominator_prime = numlib.Integrate(TE_prime, lam, 0, upper) - - return (denominator*numerator_prime - numerator*denominator_prime)/(denominator**2) - def Lum_efficacy(lam, TE): le = Lum_efficiency(lam, TE) efficacy = 683*le @@ -64,29 +51,3 @@ def IdealSource(lam, T): TE_ideal = ph * rho return TE_ideal -''' - -def luminous_efficiency(lam, thermal_emission_array): - upper = np.amax(lam) - ph = datalib.PhLum(lam) - num = datalib.PhLum(lam)*thermal_emission_array - numerator = numlib.Integrate(num, lam, 0, upper ) - den = thermal_emission_array - denominator = numlib.Integrate(den, lam, 0, upper) - luminous_efficiency_value = (numerator/denominator) - return luminous_efficiency_value - -def luminous_efficacy(lam, thermal_emission_array): - le = luminous_efficiency_value(lam, thermal_emission_array) - luminous_efficacy_value = 683*le - return luminous_efficacy_value - -def thermal_emission_ideal_source(lam, T): - ### compute blackbody spectrum at current T - rho = datalib.BB(lam, T) - ### get photopic luminosity function - ph = datalib.PhLum(lam) - ### ideal thermal emission is product of the two - thermal_emission_ideal_source = ph * rho - return thermal_emission_ideal_source_array -''' \ No newline at end of file diff --git a/wptherml/tmm.py b/wptherml/tmm.py index 0fce1f9..1b28850 100644 --- a/wptherml/tmm.py +++ b/wptherml/tmm.py @@ -25,16 +25,6 @@ def BuildP(phil): ### this function will compute the derivative of the P matrix for a given ### layer l with respect to the thickness of layer l -def Build_dP_ds(kzl, dl): - P = np.zeros((2,2),dtype=complex) - ci = 0+1j - a = -1*ci*kzl*dl - b = ci*kzl*dl - P[0][1] = 0+0j - P[1][0] = 0+0j - P[0][0] = -ci*kzl*np.exp(a) - P[1][1] = ci*kzl*np.exp(b) - return P def BuildD(nl, ctheta,pol): D = np.zeros((2,2),dtype=complex) @@ -172,151 +162,7 @@ def tmm(k0, theta0, pol, nA, tA): return M -### analytically differentiates the transfer matrix -### with respect to the thickness of a layer li to be specified by -### the user! -def d_tmm_dsi(li, k0, theta0, pol, nA, tA): - t1 = np.zeros((2,2),dtype=complex) - t2 = np.zeros((2,2),dtype=complex) - Dl = np.zeros((2,2),dtype=complex) - Dli = np.zeros((2,2),dtype=complex) - Pl = np.zeros((2,2),dtype=complex) - M = np.zeros((2,2),dtype=complex) - - D1 = BuildD(nA[0], np.cos(theta0), pol) - ### Note it is actually faster to invert the 2x2 matrix - ### "By Hand" than it is to use linalg.inv - ### and this inv step seems to be the bottleneck for the TMM function - tmp = D1[0,0]*D1[1,1]-D1[0,1]*D1[1,0] - det = 1/tmp - M[0,0] = det*D1[1,1] - M[0,1] = -det*D1[0,1] - M[1,0] = -det*D1[1,0] - M[1,1] = det*D1[0,0] - #D1i = inv(D1) - #print("D1i is ") - #print(D1i) - - - ### This is the number of layers in the structure - L = len(nA) - kz = np.zeros(L,dtype=complex) - phil = np.zeros(L,dtype=complex) - ctheta = np.zeros(L,dtype=complex) - theta = np.zeros(L,dtype=complex) - - ### since kx is conserved through all layers, just compute it - ### in the upper layer (layer 1), for which you already known - ### the angle of incidence - kx = nA[0]*k0*np.sin(theta0) - kz[0] = np.sqrt((nA[0]*k0)**2 - kx**2) - kz[L-1] = np.sqrt((nA[L-1]*k0)**2 - kx**2) - - ### keeping consistent with K-R excitation - if np.real(kz[0])<0: - kz[0] = -1*kz[0] - if np.imag(kz[L-1])<0: - kz[L-1] = -1*kz[L-1] - - ''' For dM/dsi we have three cases: - we compute contributions to M as normal for layers 2-li-1 - we compute dPi/dsi for layer li - we compute contributions to M as normal for layers li+1 - L-1 - ''' - - ''' layers 2 - li-1 ''' - for i in range(1,li): - kz[i] = np.sqrt((nA[i]*k0)**2 - kx**2) - if np.imag(kz[i])<0: - kz[i] = -1*kz[i] - - ctheta[i] = kz[i]/(nA[i]*k0) - theta[i] = np.arccos(ctheta[i]) - - phil[i] = kz[i]*tA[i] - - Dl = BuildD(nA[i],ctheta[i], pol) - ## Invert Dl - tmp = Dl[0,0]*Dl[1,1]-Dl[0,1]*Dl[1,0] - det = 1/tmp - Dli[0,0] = det*Dl[1,1] - Dli[0,1] = -det*Dl[0,1] - Dli[1,0] = -det*Dl[1,0] - Dli[1,1] = det*Dl[0,0] - #Dli = inv(Dl) - ## form Pl - Pl = BuildP(phil[i]) - - t1 = np.matmul(M,Dl) - t2 = np.matmul(t1,Pl) - M = np.matmul(t2,Dli) - - ''' layer li ''' - kz[li] = np.sqrt((nA[li]*k0)**2 - kx**2) - if np.imag(kz[li])<0: - kz[li] = -1*kz[li] - - ctheta[li] = kz[li]/(nA[li]*k0) - theta[li] = np.arccos(ctheta[li]) - - Dl = BuildD(nA[li], ctheta[li], pol) - tmp = Dl[0,0]*Dl[1,1]-Dl[0,1]*Dl[1,0] - det = 1/tmp - Dli[0,0] = det*Dl[1,1] - Dli[0,1] = -det*Dl[0,1] - Dli[1,0] = -det*Dl[1,0] - Dli[1,1] = det*Dl[0,0] - - Pl = Build_dP_ds(kz[li], tA[li]) - - t1 = np.matmul(M, Dl) - t2 = np.matmul(t1, Pl) - M = np.matmul(t2, Dli) - - - ''' layers li+1 - L-1 ''' - for i in range(li+1,(L-1)): - kz[i] = np.sqrt((nA[i]*k0)**2 - kx**2) - if np.imag(kz[i])<0: - kz[i] = -1*kz[i] - - ctheta[i] = kz[i]/(nA[i]*k0) - theta[i] = np.arccos(ctheta[i]) - - phil[i] = kz[i]*tA[i] - - Dl = BuildD(nA[i],ctheta[i], pol) - ## Invert Dl - tmp = Dl[0,0]*Dl[1,1]-Dl[0,1]*Dl[1,0] - det = 1/tmp - Dli[0,0] = det*Dl[1,1] - Dli[0,1] = -det*Dl[0,1] - Dli[1,0] = -det*Dl[1,0] - Dli[1,1] = det*Dl[0,0] - #Dli = inv(Dl) - ## form Pl - Pl = BuildP(phil[i]) - - t1 = np.matmul(M,Dl) - t2 = np.matmul(t1,Pl) - M = np.matmul(t2,Dli) - - ### M is now the product of D_1^-1 .... D_l-1^-1... just need to - ### compute D_L and multiply M*D_L - kz[L-1] = np.sqrt((nA[L-1]*k0)**2 - kx**2) - ctheta[L-1]= kz[L-1]/(nA[L-1]*k0) - DL = BuildD(nA[L-1], ctheta[L-1], pol) - t1 = np.matmul(M,DL) - ### going to create a dictionary called M which will - ### contain the matrix elements of M as well as - ### other important quantities like incoming and outgoing angles - dM_ds = {"dM11_ds": t1[0,0], - "dM12_ds": t1[0,1], - "dM21_ds": t1[1,0], - "dM22_ds": t1[1,1] - } - return dM_ds ### This function will take a thickness array ### and a scalar thickness and will return ### the index of the layer in which you are in diff --git a/wptherml/wpml.py b/wptherml/wpml.py index d2cfb30..aff74f5 100644 --- a/wptherml/wpml.py +++ b/wptherml/wpml.py @@ -112,10 +112,6 @@ def __init__(self, args): self.valid_emiss_vs_theta = [] self.validation_option = 1 - ### derivative quantities - self.reflectivity_prime_array = np.zeros(len(self.lambda_array)) - self.emissivity_prime_array = np.zeros(len(self.lambda_array)) - ### In some cases the user may wish to compute ### R, T, or eps vs angle at a specific wavelength @@ -297,39 +293,7 @@ def fresnel(self): return 1 ### currently will return derivative of reflectivity and emissivity ### wrt to thickness of layer i - def fresnel_prime(self, layer_i): - nc = np.zeros(len(self.d),dtype=complex) - for i in range(0,len(self.lambda_array)): - for j in range(0,len(self.d)): - nc[j] = self.n[j][i] - - k0 = np.pi*2/self.lambda_array[i] - ### get transfer matrix for this k0, th, pol, nc, and d - M = tmm.tmm(k0, self.theta, self.pol, nc, self.d) - dM_ds = tmm.d_tmm_dsi(layer_i, k0, self.theta, self.pol, nc, self.d) - - ### store all relevant matrix elements in variable names - M21 = M["M21"] - M11 = M["M11"] - M21p = dM_ds["dM21_ds"] - M11p = dM_ds["dM11_ds"] - - ### get reflection amplitude - r = M["M21"]/M["M11"] - r_star = np.conj(r) - - ### get derivative of reflection amplitudes - r_prime = (M11*M21p - M21*M11p)/(M11*M11) - r_prime_star = np.conj(r_prime) - - ### get reflectivity - R_prime = r_prime * r_star + r * r_prime_star - ### get Reflectivity - self.reflectivity_prime_array[i] = np.real(R_prime) - ### get Transmissivity - self.emissivity_prime_array[i] = 1 - self.reflectivity_prime_array[i] - return 1 def angular_fresnel(self, lambda_0): ### create an array for RI of each layer at the @@ -463,11 +427,6 @@ def thermal_emission_ea(self): return 1 - ''' METHOD FOR J-AGG ENHANCEMENT ''' - def jagg_sd(self): - self.jagg_sd_val = numlib.Integrate(self.emissivity_array, self.lambda_array, 500e-9, 700e-9)/200e-9 - return 1 - ''' METHODS FOR STPVLIB!!! ''' ### Normal versions first - no explicit dependence on angle