Source code for dyna.functions.Functions

# This file is part of Dyna.
#
# Dyna is a free open source software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
# Copyright (C) 2008-present CNRS, stephane.grenier@neel.cnrs.fr
"""
Functions.

@author: Stephane Grenier, Institut Néel, CNRS, France
TODO : keep M matrices in memory so as to not have to recalculate them for Transmittance or Kerr.
"""
import numpy as np
from scipy.constants import constants as cst

from dyna.functions.XrayFunc import Circ2Lin
from dyna.functions.Roughness import calc_roughness_array
import dyna.constants as dcst

[docs] def F0(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF_SIM) : r""" F0 : Parratt formalism. .. math:: \left({\frac{h}{\pi}}\right)^2 _Factor = (cst.h/np.pi)**2/(4*cst.epsilon_0)/cst.m_e * 1e30 / np.square(_EnergyRange) 2*2*np.pi/(cst.h*cst.c)*cst.e *1e-10 = 1.01355 e-3, in A-1 Parameters ---------- _Sequence : TYPE contains the sequence of the multilayer _Sample : TYPE DESCRIPTION. _XRange : TYPE _EnergyRange : TYPE DESCRIPTION. _AngleRangeRad : TYPE DESCRIPTION. _SF: TYPE Scattering Factors for the current simulation Shape is [len(Sample), len(energy), 2] Returns ------- Reflectivity ndarray shape (len(XRange), 2) """ # SAMPLE LOOP : PHASES, INDICES OF REFRACTION phase = np.zeros((len(_XRange), len(_Sample)), dtype=np.cfloat) Q = np.zeros((len(_XRange), 2, len(_Sample)), dtype=np.cfloat) _Factor = dcst.factor / np.square(_EnergyRange) for i in range(len(_Sample)): RefIndex = np.sqrt(1 - _Factor * _Sample[i].density[1] * _SF_SIM[i,:,0]) _twok_sin = 1.01355e-3 * _EnergyRange * np.sin(np.arccos(np.cos(_AngleRangeRad)/RefIndex)) Q[:, 0, i] = RefIndex * _twok_sin # SIGMA : 2kn i+1_th for i_th layer SIGMA Q[:, 1, i] = np.reciprocal(RefIndex) * _twok_sin # PI : 2k/n phase[:, i] = np.exp(1j * Q[:, 0, i] * _Sample[i].thick[1]) # SEQUENCE LOOP R = np.zeros((len(_XRange), 2), dtype=np.cfloat) for k in range(len(_Sequence)-1): # len(_Sequence)-1 : 0 is substrate _r = (Q[:, :, _Sequence[k+1]] - Q[:, :, _Sequence[k]]) / \ (Q[:, :, _Sequence[k+1]] + Q[:, :, _Sequence[k]]) * \ np.exp(-0.5 * Q[:, :, _Sequence[k+1]] * Q[:, :, _Sequence[k]] * np.square(_Sample[_Sequence[k]].rough[1])) R[:,0] = (_r[:,0] + R[:,0] * phase[:, _Sequence[k]]) / (1 + _r[:,0] * R[:,0] * phase[:, _Sequence[k]]) R[:,1] = (_r[:,1] + R[:,1] * phase[:, _Sequence[k]]) / (1 + _r[:,1] * R[:,1] * phase[:, _Sequence[k]]) return R, [] # empty array for RefCalcPhi
[docs] def F0F1SampleLoop(_Sample, dim_XRange, _EnergyRange, _AngleRangeRad, _SF): """ F0F1SampleLoop. calculates the propagation and boundary matrices, P and A, for each layer _Sample : Parameters ---------- _Sample : TYPE DESCRIPTION. dim_XRange : TYPE DESCRIPTION. _EnergyRange : TYPE DESCRIPTION. _AngleRangeRad : TYPE DESCRIPTION. _SF : TYPE = fp - i fpp with fpp positive Returns ------- A : TYPE Boundary Matrix. APhi : TYPE Boundary Matrix, magnetization inverted. P : TYPE Propagation Matrix. PPhi : TYPE Propagation Matrix, magnetization inverted. kz : TYPE DESCRIPTION. """ eps = np.zeros((dim_XRange, 3, 3), dtype=np.cfloat) u = np.zeros(3, dtype=float) eps_mag = np.zeros(dim_XRange, dtype=float) alphay = np.cos(_AngleRangeRad) alphaz = np.zeros(len(_AngleRangeRad), dtype=np.cfloat) kz = np.zeros((len(_Sample), dim_XRange), dtype=np.cfloat) # one array of _sample arrays with xrange elements A = np.zeros((len(_Sample), dim_XRange, 4, 4), dtype=np.cfloat) APhi = np.zeros((len(_Sample), dim_XRange, 4, 4), dtype=np.cfloat) P = np.zeros((len(_Sample), dim_XRange, 4, 4), dtype=np.cfloat) PPhi = np.zeros((len(_Sample), dim_XRange, 4, 4), dtype=np.cfloat) wave_k = 5.0677e-4 * _EnergyRange # 2*np.pi/(cst.h*cst.c)*cst.e *1e-10 = 1.01355/2 e-3, in A-1 _Factor = dcst.factor / np.square(_EnergyRange) # 4pi/k^2*r_0, E in eV for i in range(len(_Sample)): # i=0 is Vacuum, ... i=len(_sample) is substrate # CALCULATING THE DIELECRIC FUNCTION eps_mag = _Factor * _Sample[i].density[1] * _Sample[i].MMS[1] * _SF[i][:,1] u = [ np.sin(np.deg2rad(_Sample[i].phi[1])) * np.cos(np.deg2rad(_Sample[i].gamma[1])), np.sin(np.deg2rad(_Sample[i].phi[1])) * np.sin(np.deg2rad(_Sample[i].gamma[1])), np.cos(np.deg2rad(_Sample[i].phi[1]))] eps0 = 1 - _Factor * _Sample[i].density[1] * _SF[i][:,0] eps[:,0,0] = eps0 eps[:,0,1] = -1j * u[2] * eps_mag eps[:,0,2] = 1j * u[1] * eps_mag eps[:,1,0] = -eps[:,0,1] eps[:,1,1] = eps0 eps[:,1,2] = -1j * u[0] * eps_mag eps[:,2,0] = -eps[:,0,2] eps[:,2,1] = -eps[:,1,2] eps[:,2,2] = eps0 # CALCULATING THE REFRACTION ANGLES alphay = np.cos(_AngleRangeRad) / np.sqrt(eps[:,0,0]) alphaz = np.sqrt(1 - alphay**2) kz[i,:] = wave_k * np.sqrt(eps[:,0,0]) * alphaz # CALCULATING THE INDICES OF REFRACTION eps_mag_y = u[1] * eps_mag * alphay eps_mag_z = u[2] * eps_mag * alphaz n_right_down = np.sqrt( eps[:,0,0] + eps_mag_y - eps_mag_z ) n_left_down = np.sqrt( eps[:,0,0] - eps_mag_y + eps_mag_z ) n_right_up = np.sqrt( eps[:,0,0] + eps_mag_y + eps_mag_z ) n_left_up = np.sqrt( eps[:,0,0] - eps_mag_y - eps_mag_z ) # CALCULATES k_j PROJECTION ON Y (Snell's law) AND Z FROM Y alphay_right_down = np.cos(_AngleRangeRad)/n_right_down alphay_left_down = np.cos(_AngleRangeRad)/n_left_down alphay_right_up = np.cos(_AngleRangeRad)/n_right_up alphay_left_up = np.cos(_AngleRangeRad)/n_left_up alphaz_left_down = np.sqrt(1-alphay_left_down**2) alphaz_right_down = np.sqrt(1-alphay_right_down**2) alphaz_right_up = np.sqrt(1-alphay_right_up**2) alphaz_left_up = np.sqrt(1-alphay_left_up**2) # CALCULATING THE TRANSFER MATRIX A[i,:,0,0] = -1 - 1j * eps[:,0,1] * alphaz_right_down - 1j * eps[:,0,2] * alphay_right_down A[i,:,0,1] = 1 - 1j * eps[:,0,1] * alphaz_left_down - 1j * eps[:,0,2] * alphay_left_down A[i,:,0,2] = -1 + 1j * eps[:,0,1] * alphaz_right_up - 1j * eps[:,0,2] * alphay_right_up A[i,:,0,3] = 1 + 1j * eps[:,0,1] * alphaz_left_up - 1j * eps[:,0,2] * alphay_left_up A[i,:,1,0] = 1j * alphaz_right_down - eps[:,0,1] - 1j * eps[:,1,2] * alphay_right_down A[i,:,1,1] = 1j * alphaz_left_down + eps[:,0,1] - 1j * eps[:,1,2] * alphay_left_down A[i,:,1,2] = -1j * alphaz_right_up - eps[:,0,1] - 1j * eps[:,1,2] * alphay_right_up A[i,:,1,3] = -1j * alphaz_left_up + eps[:,0,1] - 1j * eps[:,1,2] * alphay_left_up A[i,:,2,0] = -1j * n_right_down * A[i,:,0,0] A[i,:,2,1] = 1j * n_left_down * A[i,:,0,1] A[i,:,2,2] = -1j * n_right_up * A[i,:,0,2] A[i,:,2,3] = 1j * n_left_up * A[i,:,0,3] A[i,:,3,0] = - alphaz_right_down * n_right_down * A[i,:,0,0] A[i,:,3,1] = - alphaz_left_down * n_left_down * A[i,:,0,1] A[i,:,3,2] = alphaz_right_up * n_right_up * A[i,:,0,2] A[i,:,3,3] = alphaz_left_up * n_left_up * A[i,:,0,3] APhi[i,:,0,0] = -1 + 1j * eps[:,0,1] * alphaz_left_down + 1j * eps[:,0,2] * alphay_left_down APhi[i,:,0,1] = 1 + 1j * eps[:,0,1] * alphaz_right_down + 1j * eps[:,0,2] * alphay_right_down APhi[i,:,0,2] = -1 - 1j * eps[:,0,1] * alphaz_left_up + 1j * eps[:,0,2] * alphay_left_up APhi[i,:,0,3] = 1 - 1j * eps[:,0,1] * alphaz_right_up + 1j * eps[:,0,2] * alphay_right_up APhi[i,:,1,0] = 1j * alphaz_left_down + eps[:,0,1] + 1j * eps[:,1,2] * alphay_left_down APhi[i,:,1,1] = 1j * alphaz_right_down - eps[:,0,1] + 1j * eps[:,1,2] * alphay_right_down APhi[i,:,1,2] = -1j * alphaz_left_up + eps[:,0,1] + 1j * eps[:,1,2] * alphay_left_up APhi[i,:,1,3] = -1j * alphaz_right_up - eps[:,0,1] + 1j * eps[:,1,2] * alphay_right_up APhi[i,:,2,0] = 1j * n_left_down * APhi[i,:,0,0] APhi[i,:,2,1] = -1j * n_right_down * APhi[i,:,0,1] APhi[i,:,2,2] = 1j * n_left_up * APhi[i,:,0,2] APhi[i,:,2,3] = -1j * n_right_up * APhi[i,:,0,3] APhi[i,:,3,0] = - alphaz_left_down * n_left_down * APhi[i,:,0,0] APhi[i,:,3,1] = - alphaz_right_down * n_right_down * APhi[i,:,0,1] APhi[i,:,3,2] = alphaz_left_up * n_left_up * APhi[i,:,0,2] APhi[i,:,3,3] = alphaz_right_up * n_right_up * APhi[i,:,0,3] A[i,:] = A[i,:] / (np.sqrt(2) * eps[:,0,0].reshape((dim_XRange,1,1))) APhi[i,:] = APhi[i,:] / (np.sqrt(2) * eps[:,0,0].reshape((dim_XRange,1,1))) # THE PHASE phase = wave_k * _Sample[i].thick[1] # CALCULATING THE PROPAGATION MATRIX P[i,:,0,0] = np.exp( 1j * phase * n_right_down * alphaz_right_down) P[i,:,1,1] = np.exp( 1j * phase * n_left_down * alphaz_left_down) P[i,:,2,2] = np.exp(-1j * phase * n_right_up * alphaz_right_up) P[i,:,3,3] = np.exp(-1j * phase * n_left_up * alphaz_left_up) PPhi[i,:,0,0] = P[i,:,1,1] PPhi[i,:,1,1] = P[i,:,0,0] PPhi[i,:,2,2] = P[i,:,3,3] PPhi[i,:,3,3] = P[i,:,2,2] return A, APhi, P, PPhi, kz
[docs] def F0F1XRangeLoop(_Sequence, _Sample, _XRange, A, APhi, P, PPhi, kz) : """ F0F1XRangeLoop. Runs through the sequence of layers and actually calculates the reflectivity final M matrix Parameters ---------- _Sequence : TYPE DESCRIPTION. _Sample : TYPE DESCRIPTION. _XRange : TYPE DESCRIPTION. A : TYPE DESCRIPTION. APhi : TYPE DESCRIPTION. P : TYPE DESCRIPTION. PPhi : TYPE DESCRIPTION. kz : TYPE DESCRIPTION. Returns ------- M : TYPE DESCRIPTION. MPhi : TYPE DESCRIPTION. """ rev_Sequence=list(reversed(_Sequence)) # in _rev_Seq: first is vacuum, last is substrate S = np.tile(np.identity(4, dtype=np.cfloat)[:,:], (len(_XRange),1,1)) SPhi = np.tile(np.identity(4, dtype=np.cfloat)[:,:], (len(_XRange),1,1)) for k in range(1, len(rev_Sequence)): # for all layers kk = rev_Sequence[k] # the mth layer; kk = _rev_Sequence[1] = 1 jj = rev_Sequence[k-1] # the previous layer _rev_Sequence[0] = 0 is vacuum F = np.matmul( np.linalg.inv(A[kk]), A[jj]) FPhi = np.matmul( np.linalg.inv(APhi[kk]), APhi[jj]) W = calc_roughness_array(_Sample[kk].rough[1], kz[kk], kz[jj]) F_rough = F * W FPhi_rough = FPhi * W if kk < len(_Sample)-1 : # to deal with limit conditions don't calculate for infinite layer S = np.matmul(P[kk,:], np.matmul(F_rough, S)) SPhi = np.matmul(PPhi[kk,:], np.matmul(FPhi_rough, SPhi)) else : M = np.matmul(F_rough, S) MPhi = np.matmul(FPhi_rough, SPhi) return M, MPhi
[docs] def F0F1Ref(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF) : """ F0F1Ref. Parameters ---------- _Sequence : TYPE DESCRIPTION. _Sample : TYPE DESCRIPTION. _XRange : TYPE DESCRIPTION. _EnergyRange : TYPE DESCRIPTION. _AngleRangeRad : TYPE DESCRIPTION. _SF: TYPE DESCRIPTION. Returns ------- TYPE DESCRIPTION. """ R = np.zeros((len(_XRange), 2, 2), dtype=np.cfloat) RPhi = np.zeros((len(_XRange), 2, 2), dtype=np.cfloat) # For each layer, A and P are calculated A, APhi, P, PPhi, kz = F0F1SampleLoop(_Sample, len(_XRange), _EnergyRange, _AngleRangeRad, _SF) # Calculations of the sequence of layers M, MPhi = F0F1XRangeLoop(_Sequence, _Sample, _XRange, A, APhi, P, PPhi, kz) tmp = np.divide(1, M[:,3,3] * M[:,2,2] - M[:,3,2] * M[:,2,3]) R[:, 0, 0] = (- M[:,3,3] * M[:,2,0] + M[:,2,3] * M[:,3,0]) * tmp R[:, 0, 1] = (- M[:,3,3] * M[:,2,1] + M[:,2,3] * M[:,3,1]) * tmp R[:, 1, 0] = ( M[:,3,2] * M[:,2,0] - M[:,2,2] * M[:,3,0]) * tmp R[:, 1, 1] = ( M[:,3,2] * M[:,2,1] - M[:,2,2] * M[:,3,1]) * tmp tmpPHI = np.divide(1, MPhi[:,3,3] * MPhi[:,2,2] - MPhi[:,3,2] * MPhi[:,2,3]) RPhi[:, 0, 0] = (- MPhi[:,3,3] * MPhi[:,2,0] + MPhi[:,2,3] * MPhi[:,3,0]) * tmpPHI RPhi[:, 0, 1] = (- MPhi[:,3,3] * MPhi[:,2,1] + MPhi[:,2,3] * MPhi[:,3,1]) * tmpPHI RPhi[:, 1, 0] = ( MPhi[:,3,2] * MPhi[:,2,0] - MPhi[:,2,2] * MPhi[:,3,0]) * tmpPHI RPhi[:, 1, 1] = ( MPhi[:,3,2] * MPhi[:,2,1] - MPhi[:,2,2] * MPhi[:,3,1]) * tmpPHI return Circ2Lin(R), Circ2Lin(RPhi)
[docs] def F0F1Trans(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF): """ F0F1Trans calculates the transmitivity matrix. Parameters ---------- _Sequence : TYPE DESCRIPTION. _Sample : TYPE DESCRIPTION. _XRange : TYPE DESCRIPTION. _EnergyRange : TYPE DESCRIPTION. _AngleRangeRad : TYPE DESCRIPTION. _SF: TYPE DESCRIPTION. Returns ------- TYPE DESCRIPTION. """ T = np.zeros((len(_XRange), 2, 2), dtype=np.cfloat) TPhi = np.zeros((len(_XRange), 2, 2), dtype=np.cfloat) A, APhi, P, PPhi, kz = F0F1SampleLoop(_Sample, len(_XRange), _EnergyRange, _AngleRangeRad, _SF) M, MPhi = F0F1XRangeLoop(_Sequence, _Sample, _XRange, A, APhi, P, PPhi, kz) R, RPhi = F0F1Ref(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF) T[:, 0, 0] = M[:,0,0] + M[:,0,2] * R[:, 0, 0] + M[:,0,3] * R[:, 1, 0] T[:, 0, 1] = M[:,0,1] + M[:,0,2] * R[:, 0, 1] + M[:,0,3] * R[:, 1, 1] T[:, 1, 0] = M[:,1,0] + M[:,1,2] * R[:, 0, 0] + M[:,1,3] * R[:, 1, 0] T[:, 1, 1] = M[:,1,1] + M[:,1,2] * R[:, 0, 1] + M[:,1,3] * R[:, 1, 1] TPhi[:, 0, 0] = MPhi[:,0,0] + MPhi[:,0,2] * RPhi[:, 0, 0] + MPhi[:,0,3] * RPhi[:, 1, 0] TPhi[:, 0, 1] = MPhi[:,0,1] + MPhi[:,0,2] * RPhi[:, 0, 1] + MPhi[:,0,3] * RPhi[:, 1, 1] TPhi[:, 1, 0] = MPhi[:,1,0] + MPhi[:,1,2] * RPhi[:, 0, 0] + MPhi[:,1,3] * RPhi[:, 1, 0] TPhi[:, 1, 1] = MPhi[:,1,1] + MPhi[:,1,2] * RPhi[:, 0, 1] + MPhi[:,1,3] * RPhi[:, 1, 1] return Circ2Lin(T), Circ2Lin(TPhi)
[docs] def F0F1Kerr(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF, _pi_incident) : """ F0F1Kerr calculates the Kerr effect. this is supposing that the rotation and the ellipticity are small enough to justifiy the expression. One needs to check this accordingly. Parameters ---------- sequence, Sample structure, xrange, energyrange, structure factors, pi_incident condition Returns ------- rotation angle and ellipticity """ R, RPhi = F0F1Ref(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF) if _pi_incident : _rotation = np.real(R[:, 0, 1]/R[:, 1, 1]) # Re (r_sp / r_pp) _ellipticity = np.imag(R[:, 0, 1]/R[:, 1, 1]) # Im (r_sp / r_pp else : _rotation = np.real(R[:, 1, 0]/R[:, 0, 0]) # Re (r_ps / r_ss) _ellipticity = np.imag(R[:, 1, 0]/R[:, 0, 0]) # Im (r_ps / r_ss) return 180/np.pi * _rotation, 180/np.pi * _ellipticity
[docs] def F0F1Kerr_exact(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF, _pi_incident) : """ F0F1Kerr_exact calculates the Kerr effect. calculations are exact with respect to the elliptic parameters Parameters ---------- sequence, Sample structure, xrange, energyrange, structure factors, pi_incident condition Returns ------- rotation angle and ellipticity """ R, RPhi = F0F1Ref(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF) if _pi_incident : _phi = np.angle(R[:, 0, 1]) - np.angle(R[:, 1, 1]) _A_pi = np.abs(R[:, 1, 1]) _A_sigma = np.abs(R[:, 0, 1]) _ratio = _A_sigma / _A_pi # _ratio = np.min(np.array([_A_sigma,_A_pi])) / np.max(np.array([_A_pi,_A_sigma])) _rotation = 0.5 * np.arctan( 2 * np.real(R[:, 0, 1]/R[:, 1, 1]) / (1 - np.square(_ratio))) # Re (r_ps / r_ss) _A_sigma_cos = _A_sigma * np.cos(_rotation) _A_pi_sin = _A_pi * np.sin(_rotation) _A_sigma_sin = _A_sigma * np.sin(_rotation) _A_pi_cos = _A_pi * np.cos(_rotation) _a = np.sqrt(np.square(_A_sigma_cos) + np.square(_A_pi_sin) + 2 * _A_sigma_cos * _A_pi_sin * np.cos(_phi) ) _b = np.sqrt(np.square(_A_sigma_sin) + np.square(_A_pi_cos) - 2 * _A_sigma_cos * _A_pi_sin * np.cos(_phi) ) _ellipticity = 0.5 * np.arcsin(2 * np.imag(R[:, 0, 1]/R[:, 1, 1]) / (1 + np.square(_ratio))) else : _phi = np.angle(R[:, 1, 0]) - np.angle(R[:, 0, 0]) _A_sigma = np.abs(R[:, 0, 0]) _A_pi = np.abs(R[:, 1, 0]) _ratio = _A_pi / _A_sigma _rotation = 0.5 * np.arctan( 2 * np.real(R[:, 1, 0]/R[:, 0, 0]) / (1 - np.square(_ratio))) # Re (r_ps / r_ss) _A_sigma_cos = _A_sigma * np.cos(_rotation) _A_pi_sin = _A_pi * np.sin(_rotation) _A_sigma_sin = _A_sigma * np.sin(_rotation) _A_pi_cos = _A_pi * np.cos(_rotation) _a = np.sqrt(np.square(_A_sigma_cos) + np.square(_A_pi_sin) + 2 * _A_sigma_cos * _A_pi_sin * np.cos(_phi) ) _b = np.sqrt(np.square(_A_sigma_sin) + np.square(_A_pi_cos) - 2 * _A_sigma_cos * _A_pi_sin * np.cos(_phi) ) _ellipticity = 0.5 * np.arcsin(2 * np.imag(R[:, 1, 0]/R[:, 0, 0]) / (1 + np.square(_ratio))) return 180/np.pi * _rotation, 180/np.pi * _ellipticity
[docs] def F0F1Farad(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF, _pi_incident) : """ F0F1Farad calculates the Faraday effect. This is supposing that the rotation and the ellipticity are small enough to justifiy the expression. One needs to check this accordingly. Parameters ---------- _Sequence : TYPE DESCRIPTION. _Sample : TYPE DESCRIPTION. _XRange : TYPE DESCRIPTION. _EnergyRange : TYPE DESCRIPTION. _AngleRangeRad : TYPE DESCRIPTION. _SF: TYPE DESCRIPTION. Returns ------- TYPE DESCRIPTION. """ T, TPhi = F0F1Trans(_Sequence, _Sample, _XRange, _EnergyRange, _AngleRangeRad, _SF) if _pi_incident : phi_p = np.real(T[:, 0, 1]/T[:, 1, 1]) phi_pp = np.imag(T[:, 0, 1]/T[:, 1, 1]) else : phi_p = np.real(T[:, 1, 0]/T[:, 0, 0]) phi_pp = np.imag(T[:, 1, 0]/T[:, 0, 0]) return 180/np.pi * phi_p, 180/np.pi * phi_pp