diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..7c65259 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 Umberto Zanovello + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..df09040 --- /dev/null +++ b/README.md @@ -0,0 +1,140 @@ +# CoSimPy + +CoSimPy is an open source Pyhton library aiming to combine results from electromagnetic (EM) simulation with circuits analysis through a co-simulation environment. + +## Summary + + - [Getting Started](#getting-started) + - [Deployment](#deployment) + - [License](#license) + - [Acknowledgments](#acknowledgments) + +## Getting Started + +The library has been developed with Python 3.7. and tested on previous versions down to Python 3.5. + +### Prerequisites + +The library uses the follwong additional packages: + +- [numpy](https://numpy.org/) (>=1.15.2) +- [matplotlib](https://matplotlib.org/) (>=3.0.0) +- [h5py](https://www.h5py.org/) (>=2.8.0) +- [scipy](https://www.scipy.org/) (>=1.1.0) + +The package versions reported in brackets represent the oldest releases with which the library has been succesfully tested. + +### Installing + +With [pip](https://pypi.org/project/pip/) (coming soon...): +``` +pip install cosimpy +``` + +With [anaconda](https://www.anaconda.com/products/individual): +``` +conda install --channel umbertopy cosimpy +``` + +## Deployment + +After installation, the library can be imported as: + +```python +import cosimpy +``` + +### An Example + +In the following example, a 1-port RF coil is modeled as a 5 ohm resistance in series with a 300 nH inductance. The RF coil is supposed to generate a 0.1 μT magnetic flux density oriented along the y-direction when it is supplied with 1 W incident power at 128 MHz. The coil is connected to a tuning/matching network through a 5 cm long lossless transmission line. The network is designed to transform the impedance at its output to 50 ohm at 128 MHz. + +```python +import numpy as np +import cosimpy + +L_coil = 300e-9 #Coil inductance +R_coil = 5 #Coil resistance + +#Frequency values at which the S parameters are evaluated +frequencies = np.linspace(50e6,250e6,1001) + +#Number of points along x-, y-, z-direction where the magnetic flux density is evaluated +nPoints = [20,20,20] + +#b_field is evaluated at one frequency (128 MHz) at one port +b_field = np.zeros((1,1,3,np.prod(nPoints))) +#Only the y-component is different from zero +b_field[:,:,1,:] = 0.1e-6 + +#S_Matrix instance to be associated with the RF coil instance +s_coil = cosimpy.S_Matrix.sMatrixRLseries(R_coil,L_coil,frequencies) +#EM_Field instance defined at 128 MHz to be associated with the RF coil instance +em_coil = cosimpy.EM_Field([128e6], nPoints, b_field) + +#RF_Coil instance +rf_coil = cosimpy.RF_Coil(s_coil,em_coil) + +#The average value of the y-component of the magnetic flux density +np.average(np.abs(rf_coil.em_field.b_field[0,0,1,:])).round(10) + +''' +Out: + 1e-07 +''' + +#5 cm lossless transmission line +tr_line = cosimpy.S_Matrix.sMatrixTrLine(5e-2,frequencies) + +#Connection between the RF coil and the transmission line +rf_coil_line = rf_coil.singlePortConnRFcoil([tr_line],True) + +#To design the tuning/matching network, I need to know the impedance value at 128 MHz +rf_coil_line.s_matrix[128e6].getZMatrix() + +''' +Out: + array([[[41.66705459+708.46385311j]]]) +''' + +#The impedance can be transormed to 50 ohm at 128 MHz deploying a T-network made of two capacitors and one inductor with the following values: + +Ca = 1.87e-12 #farad +Cb = 27.24e-12 #farad +L = 56.75e-9 #henry + +#I create the S_Matrix instances associated with Ca, Cb and L +S_Ca = cosimpy.S_Matrix.sMatrixRCseries(0,Ca,frequencies) +S_Cb = cosimpy.S_Matrix.sMatrixRCseries(0,Cb,frequencies) +S_L = cosimpy.S_Matrix.sMatrixRLseries(0,L,frequencies) + +#I create the S_Matrix instance of the tuning/matching network. +tun_match_network = cosimpy.S_Matrix.sMatrixTnetwork(S_Ca,S_L,S_Cb) + +#The RF coil is connected to the matching network. The capacitor Ca will be in series with the transmission line +rf_coil_line_matched = rf_coil_line.singlePortConnRFcoil([tun_match_network], True) + +#The average value of the y-component of the magnetic flux density +np.average(np.abs(rf_coil_line_matched.em_field.b_field[0,0,1,:])).round(10) + +''' +Out: + 7.825e-07 +''' + +rf_coil_line_matched.s_matrix.plotS(["S1-1"]) +``` +![](./docs/images/example_S.png) + +## License + +This project is licensed under the MIT +License - see the [LICENSE](LICENSE) file for +details. + +## Acknowledgments + +The library has been developed in the framework of the Researcher Mobility Grant (RMG) associated with the european project 17IND01 MIMAS. This RMG: 17IND01-RMG1 MIMAS has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme. + +[![](./docs/images/EMPIR_logo.jpg)](https://www.euramet.org/research-innovation/research-empir/) +[![](./docs/images/MIMAS_logo.png)](https://www.ptb.de/mimas/home/) + diff --git a/cosimpy/EM_Field.py b/cosimpy/EM_Field.py new file mode 100644 index 0000000..ad47a9f --- /dev/null +++ b/cosimpy/EM_Field.py @@ -0,0 +1,534 @@ +import numpy as np +from matplotlib import pyplot as plt +import h5py +from scipy.io import loadmat + +class EM_Field(): + + def __init__(self, freqs, nPoints, b_field=None, e_field=None, **kwargs): + + if not isinstance(freqs, np.ndarray) and not isinstance(freqs, list): + raise TypeError("Frequencies can only be an Nf list") + elif not isinstance(nPoints, np.ndarray) and not isinstance(nPoints, list): + raise TypeError("nPoints can only be a list with length equal to 3") + elif len(nPoints) != 3: + raise TypeError("nPoints can only be a list with length equal to 3") + + if e_field is None and b_field is None: + raise ValueError("At least one among e_field and b_field arguments has to be different from None") + + if e_field is not None: + if not isinstance(e_field, np.ndarray): + raise TypeError("e_field can only be np.ndarray") + elif len(e_field.shape) != 4: + raise ValueError("e_field can only be an Nf x Np x 3 x Nn matrix") + elif e_field.shape[0] != len(freqs): + raise ValueError("The frequencies list is not compatible with the e_field matrix first dimension") + elif e_field.shape[2] != 3: + raise ValueError("The third dimension of e_field is expected to be 3 (the number of field components)") + elif e_field.shape[3] != np.prod(nPoints): + raise ValueError("The fourth dimension of e_field is expected to be equal to nPoints[0]*nPoints[1]*nPoints[2]") + if b_field is not None: + if not isinstance(b_field, np.ndarray): + raise TypeError("b_field can only be np.ndarray") + elif len(b_field.shape) != 4: + raise ValueError("b_field can only be an Nf x Np x 3 x Nn matrix") + elif b_field.shape[0] != len(freqs): + raise ValueError("The frequencies list is not compatible with the b_field matrix first dimension") + elif b_field.shape[2] != 3: + raise ValueError("The third dimension of b_field is expected to be 3 (the number of field components)") + elif b_field.shape[3] != np.prod(nPoints): + raise ValueError("The fourth dimension of b_field is expected to be equal to nPoints[0]*nPoints[1]*nPoints[2]") + + for arg in kwargs.values(): + if not isinstance(arg, list) and not isinstance(arg, np.ndarray): + raise ValueError("All the additional properties have to be list or np.ndarray of floats") + elif np.array(arg).dtype != np.dtype(np.float): + raise ValueError("All the additional properties have to be list or np.ndarray of floats") + elif arg.size != np.prod(nPoints): + raise ValueError("At least one of the kwargs arguments has a length different from nPoints[0]*nPoints[1]*nPoints[2]") + + if e_field is not None: + self.__e_field = np.array(e_field,dtype = "complex") + self.__nPorts = e_field.shape[1] + else: + self.__e_field = None + if b_field is not None: + self.__b_field = np.array(b_field,dtype = "complex") + self.__nPorts = b_field.shape[1] + else: + self.__b_field = None + + self.__freqs = np.array(freqs) + self.__nPoints = nPoints + self.__n_f = len(freqs) + self.__prop = kwargs + + + @property + def e_field(self): + return self.__e_field + + + @property + def b_field(self): + return self.__b_field + + + @property + def nPoints(self): + return self.__nPoints + + + @property + def n_f(self): + return self.__n_f + + + @property + def nPorts(self): + return self.__nPorts + + + @property + def frequencies(self): + return self.__freqs + + + @property + def properties(self): + return self.__prop + + + def __repr__(self): + string = '"""""""""""""""\n EM FIELD\n"""""""""""""""\n\n' + string += "Number of frequency values = %d\nNumber of ports = %d\nNumber of point (nx, ny, nz) = %d, %d, %d\n\n"%(self.__n_f,self.__nPorts, self.__nPoints[0], self.__nPoints[1],self.__nPoints[2]) + if self.e_field is None: + string += "E field not defined\n\n" + elif self.b_field is None: + string += "B field not defined\n\n" + if self.__prop != {}: + for key in self.__prop: + string += "'%s' additional property defined\n" %key + return string + + + def compSensitivities(self): + + if self.__b_field is None: + raise ValueError("No b field property is specifed for the EM_Field instance") + + sens = np.copy(self.__b_field) + sens = np.delete(sens, 2, axis = 2) + + b1p = 0.5*(sens[:,:,0,:] + 1j*sens[:,:,1,:]) + b1m = 0.5*np.conj(sens[:,:,0,:] - 1j*sens[:,:,1,:]) + + sens[:,:,0,:] = b1p + sens[:,:,1,:] = b1m + + return sens + + + def plotB(self, comp, freq, port, plane, sliceIdx, vmin=None, vmax=None): + + if self.__b_field is None: + raise ValueError("No b field property is specifed for the EM_Field instance") + + f_idx = np.where(self.__freqs==freq)[0][0] + if f_idx is None: + raise ValueError("No B field for the specified frequency") + + if port not in np.arange(self.__nPorts) + 1: + raise ValueError("No B field for the specified port") + + if comp not in ['b1+', 'b1-']: + b = self.__b_field[f_idx, port-1,:,:] + + if comp.lower() == "x": + b = np.abs(b[0,:]) + elif comp.lower() == 'y': + b = np.abs(b[1,:]) + elif comp.lower() == 'z': + b = np.abs(b[2,:]) + elif comp.lower() == 'mag': + b = np.linalg.norm(b,axis=0) + else: + raise ValueError("comp must take one of the following values: 'mag', 'x', 'y', 'z', 'b1+', 'b1-'") + else: + b = self.compSensitivities() + b = b[f_idx, port-1,:,:] + + if comp == "b1+": + b = np.abs(b[0,:]) + else: + b = np.abs(b[1,:]) + + b = np.reshape(b, self.__nPoints, order='F') + + plt.figure("B field, Frequency %.2f MHz, Port %d, Plane %s, Index %d" %(freq*1e-6, port, plane, sliceIdx)) + + if plane == 'xy': + plt.imshow(1e6*b[:,:,sliceIdx].T,vmin=vmin,vmax=vmax) + plt.xlabel("x") + plt.ylabel("y") + elif plane == 'xz': + plt.imshow(1e6*b[:,sliceIdx,:].T,vmin=vmin,vmax=vmax) + plt.xlabel("x") + plt.ylabel("z") + elif plane == 'yz': + plt.imshow(1e6*b[sliceIdx,:,:].T,vmin=vmin,vmax=vmax) + plt.xlabel("y") + plt.ylabel("z") + + plt.title("Port %d: %s component" %(port, comp)) + cbar = plt.colorbar() + cbar.ax.set_ylabel("B field ($\mu$T)") + + + def plotE(self, comp, freq, port, plane, sliceIdx, vmin=None, vmax=None): + + if self.__e_field is None: + raise ValueError("No e field property is specifed for the EM_Field instance") + + f_idx = np.where(self.__freqs==freq)[0][0] + if f_idx is None: + raise ValueError("No E field for the specified frequency") + + if port not in np.arange(self.__nPorts) + 1: + raise ValueError("No E field for the specified port") + + if comp != "mag": + e = self.__e_field[f_idx, port-1,:,:] + + if comp == "x": + e = np.abs(e[0,:]) + elif comp == 'y': + e = np.abs(e[1,:]) + elif comp == 'z': + e = np.abs(e[2,:]) + else: + raise ValueError("comp must take one of the following values: 'mag', 'x', 'y', 'z', 'mag'") + else: + e = np.sqrt(np.abs(self.__e_field[f_idx, port-1,0,:])**2 + np.abs(self.__e_field[f_idx, port-1,1,:])**2 + np.abs(self.__e_field[f_idx, port-1,2,:])**2) + + e = np.reshape(e, self.__nPoints, order='F') + + plt.figure("E field, Frequency %.2f MHz, Port %d, Plane %s, Index %d" %(freq*1e-6, port, plane, sliceIdx)) + + if plane == 'xy': + plt.imshow(e[:,:,sliceIdx].T,vmin=vmin,vmax=vmax) + plt.xlabel("x") + plt.ylabel("y") + elif plane == 'xz': + plt.imshow(e[:,sliceIdx,:].T,vmin=vmin,vmax=vmax) + plt.xlabel("x") + plt.ylabel("z") + elif plane == 'yz': + plt.imshow(e[sliceIdx,:,:].T,vmin=vmin,vmax=vmax) + plt.xlabel("y") + plt.ylabel("z") + + plt.title("Port %d: %s component" %(port, comp)) + cbar = plt.colorbar() + cbar.ax.set_ylabel("E field (V/m)") + + + + def exportXMF(self, filename): + + dx=dy=dz=1 + nx = self.__nPoints[0] + ny = self.__nPoints[1] + nz = self.__nPoints[2] + + n_elem = (nx*ny*nz) + n_nodes = (nx+1) * (ny+1) * (nz+1) + + nodes = np.arange(n_nodes).reshape((nx+1,ny+1,nz+1),order="F") + + connections = np.empty((n_elem,8)) + + print("Arranging connection data...\n") + + #Fortan order + r = 0 + for k in range(nz): + for j in range(ny): + for i in range(nx): + print("\r%.2f %%" %(r/n_elem*100), end='') + connections[r] = nodes[i:i+2,j:j+2,k:k+2].ravel('F') + r += 1 + + #Swap elements + connections[:,[2,3]] = connections[:,[3,2]] + connections[:,[6,7]] = connections[:,[7,6]] + + x = dx*np.tile(np.tile(np.arange(nx+1),ny+1),nz+1) + y = dy*np.tile(np.arange(ny+1).repeat(nx+1),nz+1) + z = dz*np.arange(nz+1).repeat((nx+1)*(ny+1)) + + print("\n\n\nCompiling .xmf file...\n\n") + + with open(filename+".xmf", "w") as f: + f.write('\n') + f.write('\n') + f.write('\n') + f.write('\n') + f.write('\n') + f.write('\n'%n_elem) + f.write('\n'%n_elem) + f.write('%s.h5:/Mesh/Connections\n'%(filename)) + f.write('\n') + f.write('\n') + f.write('\n') + f.write('\n'%n_nodes) + f.write('%s.h5:/Mesh/Nodes_X\n'%(filename)) + f.write('\n') + f.write('\n'%n_nodes) + f.write('%s.h5:/Mesh/Nodes_Y\n'%(filename)) + f.write('\n') + f.write('\n'%n_nodes) + f.write('%s.h5:/Mesh/Nodes_Z\n'%(filename)) + f.write('\n') + f.write('\n') + + for freq in self.__freqs: + for port in range(self.__nPorts): + if self.__b_field is not None: + f.write('\n'%(freq/1e6,port+1)) + f.write('\n'%n_elem) + f.write('%s.h5:/%d_MHz/Port_%d/Breal\n'%(filename, freq/1e6, port+1)) + f.write('\n') + f.write('\n') + f.write('\n'%(freq/1e6,port+1)) + f.write('\n'%n_elem) + f.write('%s.h5:/%d_MHz/Port_%d/Bimag\n'%(filename, freq/1e6, port+1)) + f.write('\n') + f.write('\n') + if self.__e_field is not None: + f.write('\n'%(freq/1e6,port+1)) + f.write('\n'%n_elem) + f.write('%s.h5:/%d_MHz/Port_%d/Ereal\n'%(filename, freq/1e6, port+1)) + f.write('\n') + f.write('\n') + f.write('\n'%(freq/1e6,port+1)) + f.write('\n'%n_elem) + f.write('%s.h5:/%d_MHz/Port_%d/Eimag\n'%(filename, freq/1e6, port+1)) + f.write('\n') + f.write('\n') + + for p in self.__prop: + f.write('\n'%p) + f.write('\n'%n_elem) + f.write('%s.h5:/Properties/%s\n'%(filename, p)) + f.write('\n') + f.write('\n') + + f.write('\n') + f.write('\n') + f.write('\n') + + print("Compiling .h5py file...\n\n") + + with h5py.File(filename+".h5", "w") as f: + f["Mesh/Connections"] = connections + f["Mesh/Nodes_X"] = x + f["Mesh/Nodes_Y"] = y + f["Mesh/Nodes_Z"] = z + + for i,freq in enumerate(self.__freqs): + for port in range(self.__nPorts): + if self.__b_field is not None: + #Substitute nan values with zero + b_field = np.copy(self.__b_field) + b_field[np.isnan(b_field)] = 0+0j + f["%d_MHz/Port_%d/Breal"%(freq/1e6, port+1)] = np.real(b_field[i,port]).T + f["%d_MHz/Port_%d/Bimag"%(freq/1e6, port+1)] = np.imag(b_field[i,port]).T + if self.__e_field is not None: + #Substitute nan values with zero + e_field = np.copy(self.__e_field) + e_field[np.isnan(e_field)] = 0+0j + f["%d_MHz/Port_%d/Ereal"%(freq/1e6, port+1)] = np.real(e_field[i,port]).T + f["%d_MHz/Port_%d/Eimag"%(freq/1e6, port+1)] = np.imag(e_field[i,port]).T + for p in self.__prop: + f["Properties/%s"%p] = self.__prop[p] + + + def _newFieldComp(self, p_incM, phaseM): + + if p_incM.shape != phaseM.shape: + raise ValueError("p_incM and phaseM arrays are not coherent") + elif p_incM.shape[2] != self.__nPorts or phaseM.shape[2] != self.__nPorts: + raise ValueError("The number of ports has to be equal to p_incM and phase third dimension") + elif p_incM.shape[0] != self.__n_f or phaseM.shape[0] != self.__n_f: + raise ValueError("The number of frequencies of self has to be equal to p_incM and phaseM first dimension") + + norm = np.sqrt(p_incM)*np.exp(phaseM*1j) + #Move norm axis to obtain norm.shape = [self.n_f, self.nPorts, output.nPorts] + norm = np.moveaxis(norm,1,-1) + norm = np.expand_dims(norm,1)#norm.shape = [self.n_f,1, self.nPorts, output.nPorts] + norm = np.repeat(norm,3,1)#norm.shape = [self.n_f,3, self.nPorts, output.nPorts] + + efield_new = None + bfield_new = None + + if self.e_field is not None: + + #Temporary axis change so field.shape = [self.n_f, 3, self.nPoints, self.nPorts] + efield_new = np.moveaxis(self.e_field,1,-1) + efield_new = efield_new @ norm + efield_new = np.moveaxis(efield_new,-1,1) + + if self.b_field is not None: + + #Temporary axis change so field.shape = [self.n_f, 3, self.nPoints, self.nPorts] + bfield_new = np.moveaxis(self.b_field,1,-1) + bfield_new = bfield_new @ norm + bfield_new = np.moveaxis(bfield_new,-1,1) + + return EM_Field(self.__freqs, self.__nPoints, bfield_new, efield_new, **self.__prop) + + + @staticmethod + def importFields_cst(directory, freqs, nPorts, nPoints=None, Pinc_ref=1, b_multCoeff=1, pkORrms='pk', imp_efield=True, imp_bfield=True, **kwargs): + + if not imp_efield and not imp_bfield: + raise ValueError("At least one among imp_efield and imp_bfield has to be True") + elif nPoints is not None and len(nPoints) != 3: + raise TypeError("nPoints can only be None or a list with length equal to 3") + + if nPoints is None: #I try to evaluate nPoints + + if imp_efield: + x,y,z = np.loadtxt(directory+"/efield_%s_port1.txt"%(freqs[0]), skiprows=2, unpack=True, usecols=(0,1,2)) + else: + x,y,z = np.loadtxt(directory+"/bfield_%s_port1.txt"%(freqs[0]), skiprows=2, unpack=True, usecols=(0,1,2)) + + orig_len = len(x) #Total number of points + + x = np.unique(x) + y = np.unique(y) + z = np.unique(z) + + nPoints = [len(x), len(y), len(z)] + + assert np.prod(nPoints) == orig_len, "nPoints evaluation failed. Please specify its value in the method argument" + + n = np.prod(nPoints) + + if imp_efield: + e_field = np.empty((0,nPorts,3,n), dtype="complex") + else: + e_field = None + if imp_bfield: + b_field = np.empty((0,nPorts,3,n), dtype="complex") + else: + b_field = None + + if pkORrms.lower() == "pk": + rmsCoeff = 1/np.sqrt(2) + else: + rmsCoeff = 1 + + for f in freqs: + print("Importing %s MHz fields\n"%f) + if imp_efield: + e_f = np.empty((1,0,3,n), dtype="complex") + if imp_bfield: + b_f = np.empty((1,0,3,n), dtype="complex") + + for port in range(nPorts): + print("\tImporting port%d fields\n\n"%(port+1)) + if imp_efield: + e_real = np.loadtxt(directory+"/efield_%s_port%d.txt"%(f,port+1), skiprows=2, usecols=(3,4,5)) + e_imag = np.loadtxt(directory+"/efield_%s_port%d.txt"%(f,port+1), skiprows=2, usecols=(6,7,8)) + assert e_real.shape[0] == n, "At least one of e_field files is not compatible with the evaluated or passed nPoints" + e_p = np.empty((1,1,3,n), dtype="complex") + e_p[0,0,:,:] = (e_real+1j*e_imag).T + e_f = np.append(e_f,e_p,axis=1) + if imp_bfield: + b_real = np.loadtxt(directory+"/bfield_%s_port%d.txt"%(f,port+1), skiprows=2, usecols=(3,4,5)) + b_imag = np.loadtxt(directory+"/bfield_%s_port%d.txt"%(f,port+1), skiprows=2, usecols=(6,7,8)) + assert b_real.shape[0] == n, "At least one of b_field files is not compatible with the evaluated or passed nPoints" + b_p = np.empty((1,1,3,n), dtype="complex") + b_p[0,0,:,:] = (b_real+1j*b_imag).T + b_f = np.append(b_f,b_p,axis=1) + + if imp_efield: + e_field = np.append(e_field,e_f,axis=0) + if imp_bfield: + b_field = np.append(b_field,b_f,axis=0) + + if imp_efield: + e_field = np.sqrt(1/Pinc_ref) * rmsCoeff * e_field #cst exported field values are peak values + if imp_bfield: + b_field = b_multCoeff * np.sqrt(1/Pinc_ref) * rmsCoeff * b_field #cst exported field values are peak values + + freqs = 1e6*np.array(freqs).astype(np.float) + + return EM_Field(freqs, nPoints, b_field, e_field, **kwargs) + + + @staticmethod + def importFields_s4l(directory, freqs, nPorts, Pinc_ref=1, b_multCoeff=1, pkORrms='rms', imp_efield=True, imp_bfield=True, **kwargs): + + if not imp_efield and not imp_bfield: + raise ValueError("At least one among imp_efield and imp_bfield has to be True") + + if pkORrms.lower() == "pk": + rmsCoeff = 1/np.sqrt(2) + else: + rmsCoeff = 1 + + if not imp_efield: + e_field = None + + if not imp_bfield: + b_field = None + + for port in range(nPorts): + + print("\tImporting port%d fields\n\n"%(port+1)) + + if imp_efield: + + data = loadmat(directory+"/efield_port%d.mat"%(port+1)) + + if port == 0: + if data["Axis0"].shape[-1] * data["Axis1"].shape[-1] * data["Axis2"].shape[-1] == data["Snapshot0"].shape[0]: + n = data["Axis0"].shape[-1] * data["Axis1"].shape[-1] * data["Axis2"].shape[-1] + nPoints = [data["Axis0"].shape[-1], data["Axis1"].shape[-1], data["Axis2"].shape[-1]] + else: + n = (data["Axis0"].shape[-1]-1) * (data["Axis1"].shape[-1]-1) * (data["Axis2"].shape[-1]-1) + nPoints = [data["Axis0"].shape[-1]-1, data["Axis1"].shape[-1]-1, data["Axis2"].shape[-1]-1] + + e_field = np.empty((len(freqs),nPorts,3,n), dtype="complex") + + for f in range(freqs): + e_field[f,port,:,:] = np.moveaxis(data["Snapshot%d"%port],-1,0) + + if imp_bfield: + + data = loadmat(directory+"/bfield_port%d.mat"%(port+1)) + + if port == 0: + if data["Axis0"].shape[-1] * data["Axis1"].shape[-1] * data["Axis2"].shape[-1] == data["Snapshot0"].shape[0]: + n = data["Axis0"].shape[-1] * data["Axis1"].shape[-1] * data["Axis2"].shape[-1] + nPoints = [data["Axis0"].shape[-1], data["Axis1"].shape[-1], data["Axis2"].shape[-1]] + else: + n = (data["Axis0"].shape[-1]-1) * (data["Axis1"].shape[-1]-1) * (data["Axis2"].shape[-1]-1) + nPoints = [data["Axis0"].shape[-1]-1, data["Axis1"].shape[-1]-1, data["Axis2"].shape[-1]-1] + + b_field = np.empty((len(freqs),nPorts,3,n), dtype="complex") + + for f in range(len(freqs)): + b_field[f,port,:,:] = np.moveaxis(data["Snapshot%d"%f],-1,0) + + if imp_efield: + e_field = np.sqrt(1/Pinc_ref) * rmsCoeff * e_field + if imp_bfield: + b_field = b_multCoeff * np.sqrt(1/Pinc_ref) * rmsCoeff * b_field + + + return EM_Field(freqs, nPoints, b_field, e_field, **kwargs) diff --git a/cosimpy/RF_Coil.py b/cosimpy/RF_Coil.py new file mode 100644 index 0000000..e0d3ccf --- /dev/null +++ b/cosimpy/RF_Coil.py @@ -0,0 +1,85 @@ +from .S_Matrix import S_Matrix +from .EM_Field import EM_Field +import numpy as np + +class RF_Coil(): + + def __init__(self, s_matrix, em_field=None): + + if not isinstance(s_matrix, S_Matrix): + raise TypeError("s_matrix must be an instance of S_Matrix class") + + if em_field is not None: + if not isinstance(em_field, EM_Field): + raise TypeError("em_field must either be None or an instance of EM_Field class") + if em_field.nPorts != s_matrix.nPorts: + raise ValueError("The S matrix and the em_field are not compatible. There is not a e_field and b_field distribution for each port of the S matrix") + if not all(elem in s_matrix.frequencies for elem in em_field.frequencies): + raise ValueError("One or more frequencies at which the em_field is computed, differ from the frequencies at which the S matrix is defined") + + self.__s_matrix = s_matrix + self.__em_field = em_field + + + @property + def em_field(self): + return self.__em_field + + + @property + def s_matrix(self): + return self.__s_matrix + + + def __repr__(self): + string = ' """""""""""""""\n RF COIL\n """""""""""""""\n\n\n\n' + string+= '"""""""""""""""\n S MATRIX\n"""""""""""""""\n\n' + string += "|V-| = |S||V+|\n|%d x 1| = |%d x %d||%d x 1|\n\nNumber of frequency values = %d\n\n"%(self.s_matrix.nPorts,self.s_matrix.nPorts,self.s_matrix.nPorts,self.s_matrix.nPorts, self.s_matrix.n_f) + string+= '"""""""""""""""\n EM FIELD\n"""""""""""""""\n\n' + if self.__em_field is not None: + string += "Number of frequency values = %d\n\nNumber of points = (%d, %d, %d)\n\n"%(self.__em_field.n_f,self.__em_field.nPoints[0], self.__em_field.nPoints[1], self.__em_field.nPoints[2]) + if self.__em_field.e_field is None: + string += "E field not defined\n" + elif self.__em_field.b_field is None: + string += "B field not defined\n" + else: + string += 'Not defined\n' + return string + + + def singlePortConnRFcoil(self, networks, comp_Pinc=False): + + if not comp_Pinc or self.em_field is None: + + S_l = self.s_matrix._singlePortConnSMatrix(networks, False) + + return RF_Coil(S_l, None) + + else: + + S_l, p_inc, phase = self.s_matrix._singlePortConnSMatrix(networks, True) + idxs=np.where(self.s_matrix.frequencies[:,None]==self.em_field.frequencies)[0] + p_inc = p_inc[idxs] + phase = phase[idxs] + em_field_new = self.em_field._newFieldComp(p_inc, phase) + + return RF_Coil(S_l, em_field_new) + + + def fullPortsConnRFcoil(self, s_matrix_other, comp_Pinc=False): + + if not comp_Pinc or self.em_field is None: + + S_l = self.s_matrix._fullPortsConnSMatrix(s_matrix_other, False) + + return RF_Coil(S_l, None) + + else: + + S_l, p_inc, phase = self.s_matrix._fullPortsConnSMatrix(s_matrix_other, True) + idxs=np.where(self.s_matrix.frequencies[:,None]==self.em_field.frequencies)[0] + p_inc = p_inc[idxs] + phase = phase[idxs] + em_field_new = self.em_field._newFieldComp(p_inc, phase) + + return RF_Coil(S_l, em_field_new) diff --git a/cosimpy/S_Matrix.py b/cosimpy/S_Matrix.py new file mode 100644 index 0000000..5d273d0 --- /dev/null +++ b/cosimpy/S_Matrix.py @@ -0,0 +1,870 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from matplotlib import pyplot as plt +from scipy.interpolate import make_interp_spline +from copy import copy + +eps = 1e-10 + +class S_Matrix(): + + def __init__(self, S, freqs, z0=50): + + if not isinstance(S, np.ndarray): + raise TypeError("S matrix can only be np.ndarray") + elif not isinstance(freqs, np.ndarray) and not isinstance(freqs, list): + raise TypeError("Frequencies can only be an Nf list") + elif (isinstance(z0, np.ndarray) or isinstance(z0, list)) and len(z0) != S.shape[1]: + raise TypeError("The port impedances list is not compatible with the S matrix ports number") + elif len(S.shape) != 3: + raise TypeError("S matrix can only be an Nf x Np x Np matrix (also 1 x 1 x 1 is accepted in case of single load at single frequency") + elif S.shape[1] != S.shape[2]: + raise TypeError("S matrix can only be a Nf square matrices") + elif len(freqs) != S.shape[0]: + raise TypeError("The frequencies list is not compatible with the S matrix first dimension") + elif (np.round(np.abs(S),6) > 1).any(): + raise ValueError("At least one S parameter is higher than one") + + self.__S = np.array(S,dtype = "complex") + self.__f = np.array(freqs) #Hz + self.__n_f = len(freqs) + self.__nPorts = S.shape[-1] + + if not isinstance(z0, np.ndarray) and not isinstance(z0, list): + self.__z0 = z0 * np.ones(self.__nPorts) + else: + self.__z0 = np.array(z0) + + + @property + def S(self): + return self.__S + + + @property + def frequencies(self): + return self.__f + + + @property + def n_f(self): + return self.__n_f + + + @property + def nPorts(self): + return self.__nPorts + + @property + def z0(self): + return self.__z0 + + + def __repr__(self): + string = '"""""""""""""""\n S MATRIX\n"""""""""""""""\n\n' + string += "|V-| = |S||V+|\n|%d x 1| = |%d x %d||%d x 1|\n\nNumber of frequency values = %d\n"%(self.__nPorts,self.__nPorts,self.__nPorts,self.__nPorts, self.__n_f) + return string + + + def __getitem__(self, key): + + if isinstance(key,int) or isinstance(key,float): + idx = self.__findFreqIndex(key) + ret_S = np.expand_dims(self.__S[idx], axis = 0) + ret_frequencies = [self.__f[idx]] + + elif isinstance(key,tuple) or isinstance(key,list) or isinstance(key,np.ndarray): + if isinstance(key,np.ndarray) and len(key.shape) > 1: + raise SyntaxError + idxs = list(map(self.__findFreqIndex, key)) + ret_S = self.__S[idxs] + ret_frequencies = self.__f[idxs] + + elif isinstance(key,slice): + if key.start is None: + idx0 = 0 + else: + idx0 = self.__findFreqIndex(key.start) + if key.stop is None: + idx1 = self.__n_f + else: + idx1 = self.__findFreqIndex(key.stop) + + ret_S = self.__S[idx0:idx1] + ret_frequencies = self.__f[idx0:idx1] + + return S_Matrix(ret_S, ret_frequencies) + + + def __add__(self, other): + + if not isinstance(other, S_Matrix): + raise TypeError("Series operation can be performed only between two 1-port S_Matrix instances") + elif self.__nPorts != 1 or other.__nPorts != 1: + raise TypeError("Series operation can be performed only between two 1-port S_Matrix instances") + elif not np.array_equal(self.__f, other.__f): + raise ValueError("The S_Matrix instances must be defined over the same frequency values") + elif not np.array_equal(self.__z0, other.__z0): + raise ValueError("The port impedances of the S_Matrix instances must be the same") + + z0 = self.__z0[0] + + z_s = z0 * (1 + self.__S) / (1 - self.__S) + z_o = z0 * (1 + other.__S) / (1 - other.__S) + + z_r = z_s + z_o + s_r = (z_r - z0) / (z_r + z0) + + return S_Matrix(s_r, self.__f, z0) + + + def __mul__(self, other): + + if not isinstance(other, S_Matrix): + raise TypeError("Series operation can be performed only between two 1-port S_Matrix instances") + elif self.__nPorts != 1 or other.__nPorts != 1: + raise TypeError("Series operation can be performed only between two 1-port S_Matrix instances") + elif not np.array_equal(self.__f, other.__f): + raise ValueError("The S_Matrix instances must be defined over the same frequency values") + elif not np.array_equal(self.__z0, other.__z0): + raise ValueError("The port impedances of the S_Matrix instances must be the same") + + z0 = self.__z0[0] + + z_s = z0 * (1 + self.__S) / (1 - self.__S) + z_o = z0 * (1 + other.__S) / (1 - other.__S) + + z_r = (z_s * z_o) / (z_s + z_o) + s_r = (z_r - z0) / (z_r + z0) + + return S_Matrix(s_r, self.__f, z0) + + + def plotS(self, parameters, dB=True, smooth=True): + + if not isinstance(parameters,list) and not isinstance(parameters, np.ndarray): + raise TypeError("Parameters should be list or numpy array") + + eligibleParams = [] + smooth_limit = 1000 #If self.__n_f > smooth_limit it does not perform smooth + + for r in range(self.__S.shape[1]): + for c in range(self.__S.shape[2]): + eligibleParams.append('S%d-%d'%(r+1,c+1)) + for param in parameters: + if not param in eligibleParams: + raise TypeError("One of the parameters is not compliant with the S matrix") + + if self.n_f <= 2: + smooth = False + elif self.n_f == 3: + interpOrder = 2 + else: + interpOrder = 3 + + if dB: + plt.figure() + plt.title("S parameters") + plt.xlabel("Frequency (MHz)") + plt.ylabel("Magnitude (dB)") + + for param in parameters: + i = int(param[1:].split("-")[0]) - 1 + j = int(param[1:].split("-")[1]) - 1 + + if not smooth or self.__n_f > smooth_limit: + x = self.__f + y = self.S[:,i,j] + else: + x = np.linspace(self.__f[0], self.__f[-1], 1000) + spl = make_interp_spline(self.__f, self.__S[:,i,j], k=interpOrder) + y = spl(x) + + #if y = 0 I set -50 dB + dB_values = -50 * np.ones_like(y, dtype=np.float) + dB_values[y!=0] = 20*np.log10(np.abs(y[y!=0])) + + if self.__n_f == 1: + plt.scatter(x*1e-6, dB_values, label=param) + else: + plt.plot(x*1e-6, dB_values, label=param) + + plt.legend() + + else: + fig, axs = plt.subplots(2,1) + fig.suptitle("S parameters") + axs[0].set_ylabel("Magnitude") + axs[0].set_xlabel("Frequency (MHz)") + axs[1].set_ylabel("Phase (rad)") + axs[1].set_xlabel("Frequency (MHz)") + + + for param in parameters: + i = int(param[1:].split("-")[0]) - 1 + j = int(param[1:].split("-")[1]) - 1 + + if not smooth or self.__n_f > smooth_limit: + x = self.__f + y = self.__S[:,i,j] + else: + x = np.linspace(self.__f[0], self.__f[-1], 1000) + spl = make_interp_spline(self.__f, self.__S[:,i,j], k=interpOrder) + y = spl(x) + + if self.__n_f == 1: + axs[0].scatter(x*1e-6, np.abs(y), label=param) + axs[1].scatter(x*1e-6, np.angle(y), label=param) + else: + axs[0].plot(x*1e-6, np.abs(y), label=param) + axs[1].plot(x*1e-6, np.angle(y), label=param) + + axs[0].legend() + axs[1].legend() + + + def getZMatrix(self): + + z0 = np.diag(self.__z0) + + S = np.copy(self.__S) + + z0sqrt = np.sqrt(z0) + + Z = z0sqrt @ (np.eye(self.__nPorts) + S) @ np.linalg.inv(np.eye(self.__nPorts) - S) @ z0sqrt + + return Z + + + def getYMatrix(self): + + y0 = np.diag(1/self.__z0) + + S = np.copy(self.__S) + + y0sqrt = np.sqrt(y0) + + Y = y0sqrt @ np.linalg.inv(np.eye(self.__nPorts) + S) @ (np.eye(self.__nPorts) - S) @ y0sqrt + + return Y + + + def _singlePortConnSMatrix(self, networks, comp_Pinc=False): + + if len(networks) != self.__nPorts: + raise ValueError("Invalid networks length value") + + input_ports = self.__nPorts #Ports to be connected to self.S + output_ports = 0 #Ports of the returned S Matrix + output_ports_impedances = np.empty(0) + for p, network in enumerate(networks): + if network is None: + output_ports += 1 + output_ports_impedances = np.append(output_ports_impedances, self.__z0[p]) + elif not isinstance(network, S_Matrix): + raise TypeError("All the elements of networks have either to be None or an instance of S_Matrix") + elif not np.array_equal(self.__f, network.__f): + raise TypeError("All the S_Matrix instances in networks have to be defined over the same frequency values of self.S") + else: + if network.__z0[0] != self.__z0[p]: + raise TypeError("The input port impedances of the S Matrix elements of networks must be equal to the relevant port impedances of self.S") + output_ports += network.__nPorts - 1 + output_ports_impedances = np.append(output_ports_impedances, network.__z0[1:]) + + if output_ports == 0: + raise TypeError("The returned S_Matrix has to result at least into a 1 port S_Matrix") + + #The s matrix that will be used as other by the __resSMatrix method + s_forComp = np.zeros((self.__n_f, input_ports+output_ports, input_ports+output_ports), dtype='complex') + + out_idx = 0 + for in_idx, network in enumerate(networks): + if network is None: + s_forComp[:, in_idx, input_ports+out_idx] = 1 + s_forComp[:, input_ports+out_idx, in_idx] = 1 + out_idx += 1 + else: + n_c = network.__nPorts + s_forComp[:, in_idx, in_idx] = network.__S[:, 0, 0] + s_forComp[:, in_idx, input_ports+out_idx:input_ports+out_idx+n_c-1] = network.__S[:, 0, 1:] + s_forComp[:, input_ports+out_idx:input_ports+out_idx+n_c-1, in_idx] = network.__S[:, 1:, 0] + s_forComp[:, input_ports+out_idx:input_ports+out_idx+n_c-1, input_ports+out_idx:input_ports+out_idx+n_c-1] = network.__S[:, 1:, 1:] + out_idx += n_c - 1 + + s_forComp = S_Matrix(s_forComp, self.__f, np.concatenate((self.__z0, output_ports_impedances))) + + S_res = self.__resSMatrix(s_forComp) + + + if comp_Pinc: + #s matrix for closing s_forComp output ports on matched loads + s_loads_Pinc = np.zeros((self.__n_f, 2*input_ports+output_ports+1, 2*input_ports+output_ports+1)) + + for port in range(input_ports+1): + s_loads_Pinc[:, port, input_ports+output_ports+port] = 1 + s_loads_Pinc[:, input_ports+output_ports+port, port] = 1 + + p_incM = np.empty((self.__n_f, 0, self.__nPorts)) + phaseM = np.empty_like(p_incM) + + if output_ports == 1: #If there is only an output_port, no ports need to be closed on a matched load + p_inc, phase = self.__load_Pinc(s_forComp) + p_incM = np.append(p_incM, p_inc, axis=1) + phaseM = np.append(phaseM, phase, axis=1) + else: + for port in range(output_ports): + if port != 0: + s_forComp = S_Matrix.__movePort(s_forComp, input_ports, -1) + + s_loads_Pinc_M = S_Matrix(s_loads_Pinc, self.__f, np.concatenate((s_forComp.__z0,s_forComp.__z0[:input_ports+1]))) + + s_Pinc = s_forComp.__resSMatrix(s_loads_Pinc_M) + + p_inc, phase = self.__load_Pinc(s_Pinc) + p_incM = np.append(p_incM, p_inc, axis=1) + phaseM = np.append(phaseM, phase, axis=1) + + return S_res, p_incM, phaseM + + return S_res + + + def _fullPortsConnSMatrix(self, other, comp_Pinc=False): + + S_res = self.__resSMatrix(other) + + if comp_Pinc: + + #s matrix for closing the output ports of other on matched loads + s_loads_Pinc = np.zeros((self.__n_f, self.__nPorts+other.__nPorts+1, self.__nPorts+other.__nPorts+1)) + + for port in range(self.__nPorts+1): + s_loads_Pinc[:, port, other.__nPorts+port] = 1 + s_loads_Pinc[:, other.__nPorts+port, port] = 1 + + p_incM = np.empty((self.__n_f, 0, self.__nPorts)) + phaseM = np.empty_like(p_incM) + + output_ports = other.__nPorts - self.__nPorts + if output_ports == 1: #If there is only an output_port, no ports need to be closed on a matched load + p_inc, phase = self.__load_Pinc(other) + p_incM = np.append(p_incM, p_inc, axis=1) + phaseM = np.append(phaseM, phase, axis=1) + else: + s_forComp = copy(other) + + for port in range(output_ports): + if port != 0: + s_forComp = S_Matrix.__movePort(s_forComp, self.__nPorts, -1) + + s_loads_Pinc_M = S_Matrix(s_loads_Pinc, self.__f, np.concatenate((s_forComp.__z0,s_forComp.__z0[:self.__nPorts+1]))) + + s_Pinc = s_forComp.__resSMatrix(s_loads_Pinc_M) + + p_inc, phase = self.__load_Pinc(s_Pinc) + p_incM = np.append(p_incM, p_inc, axis=1) + phaseM = np.append(phaseM, phase, axis=1) + + return S_res, p_incM, phaseM + + return S_res + + + def __resSMatrix(self, other): + + if not isinstance(other, S_Matrix): + raise TypeError("The S matrix can be connected only to an instance of S_Matrix") + elif other.__nPorts < self.__nPorts + 1: + raise ValueError("The number of ports of other must be higher than (self.nPorts + 1)") + elif not np.array_equal(other.__f, self.__f): + raise ValueError("The other S matrix has to be defined on the same frequency values as the S matrix to be loaded") + elif (other.__z0[:self.__nPorts] != self.__z0).any(): + raise ValueError("The first self.nPorts of other must have the same impedances of the ports of self") + + S_C_11 = other.__S[:,:self.nPorts,:self.nPorts] + S_C_12 = other.__S[:,:self.nPorts,self.nPorts:] + S_C_21 = other.__S[:,self.nPorts:,:self.nPorts] + S_C_22 = other.__S[:,self.nPorts:,self.nPorts:] + + S_0 = np.copy(self.__S) + + + #Management of singular matrices: Look for singular matrices + + #Look for singular matrices + mask_sing = np.isclose(np.linalg.det(S_0), 0) + + if mask_sing.any(): + #Mask of manageble (A procedure) singular matrices: Array of shape equal to S_0.shape[0],S_0.shape[1]. rc_zero_mask[i,j] is True if S_0[i,:,j] == S_0[i,j,:] == 0 + rc_zero_mask = np.logical_and((S_0==0).all(axis=1),(S_0==0).all(axis=2)) + #Indices of S_0 first dimension containing singular matrices which can be managed + inv_S_0_ixds = np.where((rc_zero_mask).any(axis=1))[0] + + #I turn into np.nan S_0 at frequencies at which S_0 = is singualar + S_0[mask_sing] = np.nan * np.ones((self.nPorts, self.nPorts)) + + S_0_inv = np.linalg.inv(S_0) + + S_res = S_C_21 @ np.linalg.inv(S_0_inv - S_C_11) @ S_C_12 + S_C_22 + + + #Management of singular matrices: procedure A + + if mask_sing.any(): + S_0 = np.copy(self.__S) + for idx in inv_S_0_ixds: + #Indices of the potentially invertible sub-matrix + rc_idxs = np.where(np.logical_not(rc_zero_mask[idx]))[0] + sub_S_0 = S_0[idx][rc_idxs][:,rc_idxs] + + if np.isclose(np.linalg.det(sub_S_0), 0): #If the sub-matrix is not invertible + rc_zero_mask[idx] = False #S matrix is singular and not manageable with A procedure + else: + sub_S_0_inv = np.linalg.inv(sub_S_0) + sub_S_C_11 = S_C_11[idx][rc_idxs][:,rc_idxs] + sub_S_C_12 = S_C_12[idx][rc_idxs] + sub_S_C_21 = S_C_21[idx][:,rc_idxs] + + S_res[idx] = sub_S_C_21 @ np.linalg.inv(sub_S_0_inv - sub_S_C_11) @ sub_S_C_12 + S_C_22[idx] + + #Management of singular matrices: procedure B (I try using Z matrices) + + if mask_sing.any(): + #Mask which is True where S0 are singular and not manageable with A procedure + mask_notSolved = np.logical_and(mask_sing,np.logical_not((rc_zero_mask).any(axis=1))) + + if mask_notSolved.any(): + notSolved_idx = np.where(mask_notSolved)[0] + flagAcc = False #It becomes True if one or more original values should be changed by eps to avoid errors + flagEff = False #It becomes True if procedure B is not effective for at least one time + for idx in notSolved_idx: + S_0_extr = np.expand_dims(S_0[idx],axis=0) + S_C_extr = np.copy(other.__S)[idx] + + #I look for couple of rows i and j where row i is maed of zeros except for 1 at j column and + #row j is zeros except 1 at column i. This causes the S_Matrix.getZMatrix() method to return error + idxs_one = np.vstack(np.where(S_C_extr==1)).T + prev_idx = np.array([None, None]) + for idx_one in idxs_one: + if (idx_one[::-1] != prev_idx).all(): + r1 = np.zeros((other.__nPorts)) + r2 = np.zeros((other.__nPorts)) + r1[idx_one[1]] = 1 + r2[idx_one[0]] = 1 + if (S_C_extr[idx_one[0]] == r1).all() and (S_C_extr[idx_one[1]] == r2).all(): + S_C_extr[idx_one[0], idx_one[1]] = 1 - eps + S_C_extr[idx_one[1], idx_one[0]] = 1 - eps + if not flagAcc: + print("__resSMatrix WARNING: The matrices are singular at least at one frequency value. Results can be inaccurate at those frequencies") + prev_idx = idx_one + + S_0_extr = S_Matrix(S_0_extr, [None], self.__z0) #Fictitious S_Matrix instance to use getZMatrix method + S_C_extr = np.expand_dims(S_C_extr, axis=0) + S_C_extr = S_Matrix(S_C_extr, [None], other.__z0) #Fictitious S_Matrix instance to use getZMatrix method + + try: + Z_0 = S_0_extr.getZMatrix() + Z_C = S_C_extr.getZMatrix() + + + Z_C_11 = Z_C[0,:self.nPorts,:self.nPorts] + Z_C_12 = Z_C[0,:self.nPorts,self.nPorts:] + Z_C_21 = Z_C[0,self.nPorts:,:self.nPorts] + Z_C_22 = Z_C[0,self.nPorts:,self.nPorts:] + + Z_0_inv = np.linalg.inv(Z_0[0]) + + Z_res = Z_C_22 - Z_C_21 @ Z_0_inv @ np.linalg.inv(np.eye(self.nPorts) + Z_C_11 @ Z_0_inv) @ Z_C_12 + + S_res[idx] = S_Matrix.fromZtoS(np.expand_dims(Z_res,axis=0), [None], other.__z0[self.nPorts:]).__S[0] + + except np.linalg.LinAlgError: + + if not flagEff: + print("__resSMatrix WARNING: The matrices are singular at least at one frequency value. NaN is returned at those frequencies") + flagEff = True + + S_res = S_Matrix(S_res, self.__f, other.__z0[self.nPorts:]) + + return S_res + + + def __load_Pinc(self, other): + + if not isinstance(other, S_Matrix): + raise TypeError("The S matrix can be loaded only with an instance of S_Matrix") + elif other.__nPorts != self.__nPorts + 1: + raise ValueError("other.nPorts must be equal to self.nPorts+1") + elif not np.array_equal(other.__f, self.__f): + raise ValueError("The load S matrix has to be defined on the same frequency values as the S matrix to be loaded") + elif (self.__z0 != other.__z0[:self.__nPorts]).any(): + raise ValueError("The first self.nPorts of other must have the same impedances of the ports of self") + + z0 = self.__z0 + + S_O = self.__S + S_CL_11 = other.__S[:,:-1,:-1] + S_CL_12 = other.__S[:,:-1,-1:] + + V_O_p = (np.linalg.inv(np.eye(self.__nPorts) - S_CL_11 @ S_O) @ S_CL_12) * np.sqrt(np.real(other.__z0[-1])) + + #For testing + # V_O_m = S_O @ V_O_p + # V_O = V_O_p + V_O_m + # Z_O = self.getZMatrix() + # I_O = np.linalg.inv(Z_O) @ V_O + # Vg = V_O + z0m @ I_O + + p_inc = np.abs(V_O_p)**2 / (np.repeat(np.expand_dims(np.real(z0),axis=(0,2)), self.__n_f, axis=0)) + phase = np.angle(V_O_p) + + return np.swapaxes(p_inc, 1, 2), np.swapaxes(phase, 1, 2) + + + def __findFreqIndex(self, freq): + + if freq in self.__f: + idx = np.where(self.__f == freq)[0][0] + else: + idx = np.argmin(np.abs(self.__f - freq)) + print("WARNING: %e Hz is not contained in the frequencies list. %e Hz is returned instead" %(freq, self.__f[idx])) + + return idx + + + @staticmethod + def fromZtoS(Z, freqs, z0 = 50): + + if isinstance(z0, list) or isinstance(z0, np.ndarray): + if len(z0) != Z.shape[-1]: + raise ValueError("z0 has to be a list or numpy array with a length equal to the number of ports associated with the S matrix") + + y0 = np.diag(1./np.array(z0)) + + else: + + y0 = 1./z0 * np.eye(Z.shape[-1]) + + y0sqrt = np.sqrt(y0) + + S = (y0sqrt @ Z @ y0sqrt - np.eye(Z.shape[-1])) @ np.linalg.inv(y0sqrt @ Z @ y0sqrt + np.eye(Z.shape[-1])) + + return S_Matrix(S,freqs,z0) + + + @staticmethod + def fromYtoS(Y, freqs, y0 = 0.02): + + if isinstance(y0, list) or isinstance(y0, np.ndarray): + if len(y0) != Y.shape[-1]: + raise ValueError("y0 has to be a list or numpy array with a length equal to the number of ports associated with the S matrix") + + z0 = np.diag(1./np.array(y0)) + + else: + + z0 = 1./y0 * np.eye(Y.shape[-1]) + + z0sqrt = np.sqrt(z0) + + S = (np.eye(Y.shape[-1]) - z0sqrt @ Y @ z0sqrt) @ np.linalg.inv(np.eye(Y.shape[-1]) + z0sqrt @ Y @ z0sqrt) + + return S_Matrix(S,freqs,np.diag(z0)) + + + @staticmethod + def importTouchstone(filename, fmt='R_I', freqUnit='MHz', fix=False, n_f=None): + + if fmt not in ['R_I', 'Mag_Deg', 'dB_Deg', 'Mag_Rad', 'dB_Rad']: + raise ValueError("fmt is not valid") + elif freqUnit not in ['Hz', 'MHz']: + raise ValueError("freqUnit is not valid") + + if fix: + if n_f is None: + raise ValueError("Please specify the number of frequency values 'n_f' in order to fix the touchstone file") + with open(filename, 'r') as f: + data = f.readlines() + + comm_rows = 0 + fc = data[0][0] + while fc in ["!", "#"]: + comm_rows += 1 + fc = data[comm_rows][0] + + comments = "".join(data[:comm_rows]) + data = data[comm_rows:] + + if len(data) % n_f != 0: #Check that each frequency value cooresponds to the same number of rows + raise ValueError("Impossible to fix the touchstone file") + else: + n_r = int(len(data) / n_f) + + new_data = np.array(data[0::n_r]) + + for i in range(n_r-1): + new_data = np.core.defchararray.add(new_data, data[i+1::n_r]) + + new_data = np.chararray.split(new_data) + new_data = np.stack(new_data,axis=0) + new_data = new_data.astype(np.float) + + filename = filename.split(".") + filename = ".".join(filename[:-1]) + "_fixed." + filename[-1] + + np.savetxt(filename, new_data, header=comments, fmt="%.6f", delimiter="\t") + + data = np.loadtxt(filename, comments=["#","!"]) + if freqUnit == 'Hz': + freq = data[:,0] + else: + freq = 1e6*data[:,0] + + s_parameters_A = data[:,1:][:,::2] + s_parameters_B = data[:,2:][:,::2] + + if fmt == 'R_I': + s_parameters = s_parameters_A + 1j*s_parameters_B + elif fmt == 'Mag_Deg': + s_parameters = s_parameters_A * np.exp(1j*np.deg2rad(s_parameters_B)) + elif fmt == 'Mag_Rad': + s_parameters = s_parameters_A * np.exp(1j*s_parameters_B) + elif fmt == 'dB_Deg': + s_parameters = (10**(s_parameters_A/20)) * np.exp(1j*np.deg2rad(s_parameters_B)) + elif fmt == 'dB_Rad': + s_parameters = (10**(s_parameters_A/20)) * np.exp(1j*s_parameters_B) + + s_parameters = np.reshape(s_parameters, (freq.size, np.int(np.sqrt(s_parameters_A.shape[-1])),-1)) + + return S_Matrix(s_parameters, freq) + + + @staticmethod + def sMatrixRCseries(Rvalue, Cvalue, freqs, z0=50): + + if Rvalue < 0 or Cvalue <= 0: + raise ValueError("Rvalue and Cvalue must be equal or higher than zero") + + + f = np.array(freqs) + n_f = len(f) + + z = (Rvalue - 1.j/(2*np.pi*f*Cvalue)).reshape(n_f,1,1) + + s = (z - z0) / (z + z0) + + return S_Matrix(s, f, z0) + + + @staticmethod + def sMatrixRLseries(Rvalue, Lvalue, freqs, z0=50): + + if Rvalue < 0 or Lvalue < 0: + raise ValueError("Rvalue and Lvalue must be equal or higher than zero") + + f = np.array(freqs) + n_f = len(f) + + z = (Rvalue + 1.j*2*np.pi*f*Lvalue).reshape(n_f,1,1) + + s = (z - z0) / (z + z0) + + return S_Matrix(s, f, z0) + + + @staticmethod + def sMatrixRCparallel(Rvalue, Cvalue, freqs, z0=50): + + if Rvalue < 0 or Cvalue <= 0: + raise ValueError("Rvalue and Cvalue must be equal or higher than zero") + + f = np.array(freqs) + n_f = len(f) + + tau = Rvalue*Cvalue + omega = 2*np.pi*f + z = (Rvalue*(1. - 1j*omega*tau) / ((omega*tau)**2 + 1)).reshape(n_f,1,1) + + s = (z - z0) / (z + z0) + + return S_Matrix(s, f, z0) + + + @staticmethod + def sMatrixRLparallel(Rvalue, Lvalue, freqs, z0=50): + + if Rvalue < 0 or Lvalue < 0: + raise ValueError("Rvalue and Lvalue must be equal or higher than zero") + + + f = np.array(freqs) + n_f = len(f) + + tau = Lvalue / Rvalue + omega = 2*np.pi*f + + z = ((omega*Lvalue) * (omega*tau + 1.j) / ((omega*tau)**2 + 1)).reshape(n_f,1,1) + + s = (z - z0) / (z + z0) + + return S_Matrix(s, f, z0) + + + @staticmethod + def sMatrixTnetwork(SL_1, SL_2, ST, z0=50): + + if not isinstance(SL_1, S_Matrix) or not isinstance(SL_2, S_Matrix) or not isinstance(ST, S_Matrix): + raise TypeError("SL_1, SL_2 and ST must be S_Matrix instances") + elif not np.array_equal(SL_1.__f, SL_2.__f) or not np.array_equal(SL_1.__f, ST.__f): + raise ValueError("SL_1, SL_2 and ST must be defined over the same frequency values") + elif (SL_1.nPorts != 1) or (SL_2.nPorts != 1) or (ST.nPorts != 1): + raise ValueError("SL_1, SL_2 and ST must be one port S_Matrix instances") + + if isinstance(z0,list) or isinstance(z0,np.ndarray): + if len(z0) != 2: + raise ValueError("z0 must be a complex value or a list with length equal to 2") + z0 = np.array(z0) + else: + z0 = z0*np.ones(2) + + f = SL_1.__f + + #Try to avoid some possible Singular Matrix errors + if (SL_1.__S==1).any() or (SL_2.__S==1).any() or (ST.__S==1).any(): + print("sMatrixTnetwork WARNING: Singular matrix error detected at least at one frequency value. Inaccurate values can be obtained at those frequencies") + SL_1.__S[SL_1.__S==1] = 1 - eps + SL_2.__S[SL_2.__S==1] = 1 - eps + ST.__S[ST.__S==1] = 1 - eps + + ZL_1 = SL_1.getZMatrix() + ZL_2 = SL_2.getZMatrix() + ZT = ST.getZMatrix() + + Z11 = ZL_1 + ZT + Z12 = ZT + Z21 = ZT + Z22 = ZL_2 + ZT + + Z1 = np.concatenate((Z11, Z12), axis=2) + Z2 = np.concatenate((Z21, Z22), axis=2) + + z = np.concatenate((Z1, Z2), axis=1) + + return S_Matrix.fromZtoS(z, f, z0) + + + @staticmethod + def sMatrixPInetwork(ST_1, ST_2, SL, z0=50): + + if not isinstance(ST_1, S_Matrix) or not isinstance(ST_2, S_Matrix) or not isinstance(SL, S_Matrix): + raise TypeError("ST_1, ST_2 and SL must be S_Matrix instances") + elif not np.array_equal(ST_1.__f, ST_2.__f) or not np.array_equal(ST_1.__f, SL.__f): + raise ValueError("ST_1, ST_2 and SL must be defined over the same frequency values") + + if isinstance(z0,list) or isinstance(z0,np.ndarray): + if len(z0) != 2: + raise ValueError("z0 must be a complex value or a list with length equal to 2") + y0 = 1. / np.array(z0) + else: + y0 = (1./z0) * np.ones(2) + + f = ST_1.__f + + #Try to avoid some possible Singular Matrix errors + if (ST_1.__S==-1).any() or (ST_2.__S==-1).any() or (SL.__S==-1).any(): + print("sMatrixPInetwork WARNING: Singular matrix error detected at least at one frequency value. Inaccurate values can be obtained at those frequencies") + ST_1.__S[ST_1.__S==-1] = -1 + eps + ST_2.__S[ST_2.__S==-1] = -1 + eps + SL.__S[SL.__S==-1] = -1 + eps + + YT_1 = ST_1.getYMatrix() + YT_2 = ST_2.getYMatrix() + YL = SL.getYMatrix() + + Y11 = YT_1 + YL + Y12 = -YL + Y21 = -YL + Y22 = YT_2 + YL + + Y1 = np.concatenate((Y11, Y12), axis=2) + Y2 = np.concatenate((Y21, Y22), axis=2) + + y = np.concatenate((Y1, Y2), axis=1) + + s = S_Matrix.fromYtoS(y, f, y0) + + return s + + + @staticmethod + def sMatrixTrLine(l, freqs, z0=50, c_f = 1, alpha=0): + + f = np.array(freqs) + wl = 3e8 * c_f * 1/f + beta = 2 * np.pi / wl + + gamma = alpha + 1j*beta + + S11 = np.zeros((len(f),1,1)) + S12 = np.moveaxis([[np.exp(-gamma*l)]],-1,0) + S21 = S12 + S22 = S11 + + S1 = np.concatenate((S11, S12), axis=2) + S2 = np.concatenate((S21, S22), axis=2) + + s = np.concatenate((S1, S2), axis=1) + + return S_Matrix(s, f, z0) + + + @staticmethod + def __movePort(Smat, idx0, idx1): + + if not isinstance(Smat, S_Matrix): + raise TypeError("Smat must be an S_Matrix instance") + + if idx0 < 0: + idx0 = Smat.__nPorts+idx0 + if idx1 < 0: + idx1 = Smat.__nPorts+idx1 + + tempS = np.copy(Smat.__S) + tempz0 = np.copy(Smat.__z0) + + if idx0 < idx1: + + #move row + tempS[:,idx0:idx1,:] = Smat.__S[:,idx0+1:idx1+1,:] + tempS[:,idx1,:] = Smat.__S[:,idx0,:] + tempS[:,idx1+1:,:] = Smat.__S[:,idx1+1:,:] + + #move column + S_new = np.copy(tempS) + S_new[:,:,idx0:idx1] = tempS[:,:,idx0+1:idx1+1] + S_new[:,:,idx1] = tempS[:,:,idx0] + + #fix z0 + tempz0[idx0:idx1+1] = np.roll(tempz0[idx0:idx1+1],-1) + + elif idx0 > idx1: + + #move row + tempS[:,idx1,:] = Smat.__S[:,idx0,:] + tempS[:,idx1+1:idx0+1,:] = Smat.__S[:,idx1:idx0,:] + tempS[:,idx0+1:,:] = Smat.__S[:,idx0+1:,:] + + #move column + S_new = np.copy(tempS) + S_new[:,:,idx1] = tempS[:,:,idx0] + S_new[:,:,idx1+1:idx0+1] = tempS[:,:,idx1:idx0] + + #fix z0 + tempz0[idx0:idx1+1] = np.roll(tempz0[idx0:idx1+1],1) + + else: + + S_new = tempS + + return S_Matrix(S_new, Smat.__f, tempz0) diff --git a/cosimpy/__init__.py b/cosimpy/__init__.py new file mode 100644 index 0000000..e65fbc1 --- /dev/null +++ b/cosimpy/__init__.py @@ -0,0 +1,3 @@ +from .RF_Coil import RF_Coil +from .S_Matrix import S_Matrix +from .EM_Field import EM_Field diff --git a/docs/Documentation.md b/docs/Documentation.md new file mode 100644 index 0000000..f437973 --- /dev/null +++ b/docs/Documentation.md @@ -0,0 +1,1117 @@ +# CoSimPy +CoSimPy is an open source Pyhton library aiming to combine results from electromagnetic (EM) simulation with circuit analysis through a co-simulation environment. +The library is developed specifically to deal with Magnetic Resonance Imaging (MRI) radiofrequency (RF) coils. Nevertheless, it is sufficiently general to be adopted in any context involving EM multiport simulations. + +The library is composed of three different classes: +* [S_Matrix class](#s_matrix-class) +* [EM_Field class](#em_field-class) +* [RF_Coil class](#rf_coil-class) + +The **S_Matrix** class is defined to manage the scattering parameters associated with any specific circuit. These can be either the scattering parameters associated with the simulated RF coil or with any other circuit which has to be connected. Specific methods are defined to handle the connection between different instances of its class. + +The **EM_Field** class is defined to manage the EM field generated by any port of the simulated device. + +The **RF_Coil** class represents the highest level class of the library. An instance of the RF_Coil class virtually represents the simulated device being, its properties, an *S_Matrix* instance managing the scattering parameters associated with the RF coil and an *EM_Field* instance managing the EM field generated by any port of the RF coil. This class is designed to handle the interaction between the scattering parameters and the EM field allowing the computation of the latter when the RF coil is connected to any other passive device. +___ +## S_Matrix class +This class is defined to manage the scattering matrix (S matrix) associated with any specific circuit. These can be either the scattering parameters associated with the simulated RF coil or with any other circuit which has to be connected. + +### Properties +* `S_Matrix.S` : *numpy ndarray*
+Nf 🞩 NP 🞩 NP *numpy ndarray* representing an NP ports S matrix defined over Nf frequency values +* `S_Matrix.frequencies` : *numpy ndarray*
+Nf *numpy ndarray* representing the frequency values, in hertz, over which the scattering parameters are defined +* `S_Matrix.n_f` : *int*
+number of frequency values over which the scattering parameters are defined (Nf) +* `S_Matrix.nPorts` : *int*
+number of ports associated with the S Matrix (NP) + +### Methods + +#### `__init__(self, S, freqs, z0=50)` + +It creates an *S_Matrix* instance. + +Parameters + +* S : *numpy ndarray*
+Nf 🞩 NP 🞩 NP *numpy ndarray* representing an NP ports S Matrix defined over Nf frequency values +* freqs : *list* or *numpy ndarray*
+Nf *list* or *numpy ndarray* representing the frequency values, in hertz, over which the scattering parameters are defined +* z0 : *float*, *list* or *numpy ndarray*, *optional*
+port impedances in ohm. These can be given as a NP *list* or *numpy ndarray*. If all the ports share the same impedances, a *float* value can be passed as parameter. Default is 50 ohm + +Returns + +* S_Matrix : *S_Matrix* + +Example +```python +import numpy as np +from cosimpy import * + +n_p = 5 #Number of ports +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False) + +s_real = np.random.rand(n_f,n_p,n_p) +s_imag = np.random.rand(n_f,n_p,n_p) +s = (s_real + 1j*s_imag) / np.max(np.abs(s_real + 1j*s_imag)) + +S_Matrix(s, frequencies) + +''' +Out: + + """"""""""""""" + S MATRIX + """"""""""""""" + + |V-| = |S||V+| + |5 x 1| = |5 x 5||5 x 1| + + Number of frequency values = 10 + +''' +``` + +#### `__repr__(self)` + +Method for returning a printable representation of the *S_Matrix* instance. + +Parameters + +* self : *S_Matrix* + +Returns + +* string : *string*
+The string identifies the class of the instance. It gives a textual representation of the S matrix shape and reports the number of frequency values over which the scattering parameters are defined + +#### `__getitem__(self, key)` + +Indexing method for *S_Matrix* instances. Indices are interpreted as frequency values. + +Parameters + +* key : *tuple*, *list*, *numpy ndarray*, *int*, *float*, *slice*
+frequency values used to index the `self` *S_Matrix* + +Returns + +* S_Matrix : *S_Matrix* + +Example +```python +import numpy as np +from cosimpy import * + +n_p = 5 #Number of ports +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +S_Mat = S_Matrix(np.random.uniform(size=(n_f,n_p,n_p)), frequencies) + +S_Mat.frequencies + +''' +Out: + + array([5.0e+07, 6.0e+07, 7.0e+07, 8.0e+07, 9.0e+07, 1.0e+08, 1.1e+08, + 1.2e+08, 1.3e+08, 1.4e+08]) +''' + +S_Mat[90e6:].frequencies + +''' +Out: + + array([9.0e+07, 1.0e+08, 1.1e+08, 1.2e+08, 1.3e+08, 1.4e+08]) +''' + +S_Mat[90e6:130e6].frequencies + +''' +Out: + + array([9.0e+07, 1.0e+08, 1.1e+08, 1.2e+08]) +''' +S_Mat[95e6].frequencies + +''' +Out: + +'WARNING: 9.500000e+07 Hz is not contained in the frequencies list. 9.000000e+07 Hz is returned instead' + array([90000000.]) + +''' +``` + +#### `__add__(self, other)` + +\_\_add__ method for *S_Matrix* instances. It computes the *S_Matrix* resultant from the series among 1-port `self` and `other` *S_Matrix* instances. `self` and `other` must be defined over the same frequency values and must share the same port impedance. + +Parameters + +* self : *S_Matrix*
+1-port *S_Matrix* +* other : *S_Matrix*
+1-port *S_Matrix* + +Returns + +* S_Matrix : *S_Matrix*
+ 1-port *S_Matrix* + +Example +```python +import numpy as np +from cosimpy import * + +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +S_Mat1 = S_Matrix(np.random.uniform(size=(n_f,1,1)), frequencies) +S_Mat2 = S_Matrix(-1*np.ones((n_f,1,1)), frequencies) #S_Mat2: short circuit + +S_res = S_Mat1 + S_Mat2 + +(np.round(S_res.S,6) == np.round(S_Mat1.S,6)).all() + +''' +Out: + + True + +''' +``` + +#### `__mul__(self, other)` + +\_\_mul__ method for *S_Matrix* instances. It computes the *S_Matrix* resultant from the parallel among 1-port `self` and `other` *S_Matrix* instances. `self` and `other` must be defined over the same frequency values and must share the same port impedance. + +Parameters + +* self : *S_Matrix*
+1-port *S_Matrix* +* other : *S_Matrix*
+1-port *S_Matrix* + +Returns + +* S_Matrix : *S_Matrix*
+ 1-port *S_Matrix* + +Example +```python +import numpy as np +from cosimpy import * + +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +S_Mat1 = S_Matrix(np.random.uniform(size=(n_f,1,1)), frequencies) +S_Mat2 = S_Matrix(-1*np.ones((n_f,1,1)), frequencies) #S_Mat2: short circuit + +S_res = S_Mat1 * S_Mat2 + +(S_res.S == -1).all() + +''' +Out: + + True + +''' +``` + +#### `plotS(self, parameters, dB=True, smooth=True)` + +The method is used to plot the scattering parameters as a function of the frequency values over which they are defined. + +Parameters + +* parameters : *list* or *numpy ndarray*
+*list* or *numpy ndarray* identifying the scattering parapeters to plot. Each element of parameters must be a string formatted as "S\-\" where \ and \ are the relevant port numbers +* dB : *bool*, *optional*
+if `True`, the scattering parameters are plotted in Decibel. Otherwise, they are plotted as magnitude and phase. Default is `True` +* smooth : *bool*, *optional*
+if `True`, the scattering parameters are interpolated along their frequency range to produce smoother plots. Default is `True` + +Returns + +* None + +Example +```python +import numpy as np +from cosimpy import * + +n_p = 5 #Number of ports +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +S_Mat = S_Matrix(np.random.uniform(size=(n_f,n_p,n_p)), frequencies) + +S_Mat.plot(["S1-1, S1-2, S3-2"]) +``` + +#### `getZMatrix(self)` + +it returns the Z matrix associated with the *S_Matrix* `self` + +Parameters + +* self : *S_Matrix* + +Returns + +* Z : *numpy ndarray*
+ Nf 🞩 NP 🞩 NP *numpy ndarray* of impedance parameters. Nf is the number of frequency values over which `self` is defined and NP is the number of ports + +#### `getYMatrix(self)` + +it returns the Y matrix associated with the *S_Matrix* `self` + +Parameters + +* self : *S_Matrix* + +Returns + +* Y : *numpy ndarray*
+ Nf 🞩 NP 🞩 NP *numpy ndarray* of admittance parameters. Nf is the number of frequency values over which `self` is defined and NP is the number of ports + +#### `_singlePortConnSMatrix(self, networks, comp_Pinc=False)` + +The method performs the cascade connection between `self` and the *S_Matrix* instances contained in `networks` *list*. + +Parameters + +* self : *S_Matrix* +* networks : *list*
+*list*, with length equal to the number of ports of `self`, either containing `None` values or *S_Matrix* instances defined over the same frequency values of `self`. If the n-th element of `networks` is `None`, the n-th port of `self` is not connected to any other network in the returned *S_Matrix*. If the n-th element of `networks` is an N-ports *S_Matrix*, its first port must share the same impedance of the n-th port of `self` and the two ports are connected together. In the returned *S_Matrix*, the n-th port of `self` is therefore expanded into (N-ports - 1) ports. The example proposed in the image below should clarify the workflow + +![image](./images/singlePortConnSMatrix.png) + +* comp_Pinc : *bool*, *optional*
+ if `True` the method returns the magnitude and 'phase' of the incident powers across the ports of `self` when each port of the returned S matrix is supplied with 1 W of incident power and all its other ports are closed on characteristic loads. Default is `False` + +Returns + +* S_res : *S_Matrix*
+ The *S_Matrix* which results from the connection. The *S_Matrix* is defined over the same frequency values of `self` and has a number of ports equal to NP - NT + ∑i=1NT(Np,i - 1) where NP is the number of ports of `self`, NT is the number of *S_Matrix* instances in `networks` and Np,i is the number of ports of the i-th *S_Matrix* of `networks`. The n-th port impedance of S_res is equal to the impedance of the relevant port of `self` if it derives from a `None` element of `networks`. Otherwise, it is equal to the impedance of the relevant port of the *S_Matrix* of `networks` from which the port derives +* p_incM : *numpy ndarray*, *optional*
+ Nf 🞩 NPout 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined, NPout is the number of ports of the returned *S_Matrix* and NP is the number of ports of `self`. PincM[i,j,k] represents the magnitude of the incident power at the k-port of `self` when the j-port of the output *S_Matrix* is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` +* phaseM : *numpy ndarray*, *optional*
+ Nf 🞩 NPout 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined, NPout is the number of ports of the returned *S_Matrix* and NP is the number of ports of `self`. + phaseM[i,j,k] represents the 'phase' of the incident power at the k-port of `self` when the j-port of the output *S_Matrix* is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` + +Example +```python +import numpy as np +from cosimpy import * + +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +s = np.repeat(np.array([[[0, 1], [1, 0]]]), n_f, axis=0) #Straight connection between port 1 and port 2 +S_Mat = S_Matrix(s, frequencies) +S_ch = S_Matrix(np.zeros((n_f,1,1)), frequencies) #Characteristic load + +S_res, p_incM, phaseM = S_Mat._singlePortConnSMatrix([None, S_ch], True) + +np.array_equal(S_res.S,np.repeat([[[0]]],10,axis=0)) + +''' +Out: + + True +''' + +np.array_equal(np.round(p_incM,10),np.repeat([[[1., 0.]]],n_f,axis=0)) + +''' +Out: + + True +''' + +np.equal(phaseM, np.zeros_like(phaseM)).all() + +''' +Out: + + True + +''' +``` + +#### `_fullPortsConnSMatrix(self, other, comp_Pinc=False)` + +The method performs the cascade connection between `self` and `other` *S_Matrix* instances. + +Parameters + +* self : *S_Matrix* +* other : *S_Matrix*
+*S_Matrix* instance defined over the same frequency values of `self` and with a number of ports higher than those of `self`. In the returned *S_Matrix*, the first NP ports of `other` are connected with the NP ports of `self` and must share the same port impedances +* comp_Pinc : *bool*, *optional*
+if `True` the method returns the magnitude and 'phase' of the incident powers across the ports of `self` when each port of the returned S matrix is supplied with 1 W of incident power and all its other ports are closed on characteristic loads. Default is `False` + +Returns + +* S_res : *S_Matrix*
+ The *S_Matrix* which results from the connection. The *S_Matrix* is defined over the same frequency values of `self` and has a number of ports equal to NPo - NP being NP the number of ports of `self` and NPo the number of ports of `other`. The port impedances of S_res will be equal to those of the last NPo ports of other +* p_incM : *numpy ndarray*, *optional*
+ Nf 🞩 NPout 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined, NPout is the number of ports of the returned *S_Matrix* and NP is the number of ports of `self`. + PincM[i,j,k] represents the magnitude of the incident power at the k-port of `self` when the j-port of the output *S_Matrix* is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` +* phaseM : *numpy ndarray*, *optional*
+ Nf 🞩 NPout 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined, NPout is the number of ports of the returned *S_Matrix* and NP is the number of ports of `self`. + phaseM[i,j,k] represents the 'phase' of the incident power at the k-port of `self` when the j-port of the output *S_Matrix* is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` + + Example +```python +import numpy as np +from scipy.linalg import circulant +from cosimpy import * + +n_p = 2 #Number of ports +n_f = 10 #Number of frequency values + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +S_Mat = S_Matrix(np.random.uniform(size=(n_f,n_p,n_p)), frequencies) + +s_other = circulant([0,0,1,0]).T #Straight connection between ports 1 with port 3 and port 2 with port 4 + +s_other + +''' +Out: + + array([[0, 0, 1, 0], + [0, 0, 0, 1], + [1, 0, 0, 0], + [0, 1, 0, 0]]) +''' + +other = S_Matrix(np.repeat([s_other], n_f, axis=0), frequencies) + +S_res, p_incM, phaseM = S_Mat._fullPortsConnSMatrix(other, True) + +(np.round(S_res.S,10) == np.round(S_Mat.S,n_f)).all() + +''' +Out: + + True +''' + +np.equal(np.round(p_incM, 10),np.repeat([[[1,0],[0,1]]], n_f, axis=0)).all() + +''' +Out: + + True +''' + +np.equal(phaseM, np.zeros_like(phaseM)).all() + +''' +Out: + + True + +''' +``` + +#### `__resSMatrix(self, other)` + +Private method used to compute the *S_Matrix* instance resulting from the connection of the NP ports of `self` with the first NP ports of `other`. + +Parameters + +* self : *S_Matrix* +* other : *S_Matrix*
+S_Matrix instance* defined over the same frequency values of `self` and with a number of ports higher than those of `self`. In the returned *S_Matrix*, the first NP ports of `other` are connected with the NP ports of `self` and must share the same port impedances + + +Returns + +* S_res : *S_Matrix*
+ The *S_Matrix* which results from the connection of `self` with `other`. The *S_Matrix* is defined over the same frequency values of `self` and has a number of ports equal to NPo - NP being NP the number of ports of `self` and NPo the number of ports of `other`. The port impedances of S_res will be equal to those of the last NPo ports of other + +#### `__load_Pinc(self, other, comp_Pinc=False)` + +Private method used to compute the magnitude and 'phase' of the powers incident to the ports of `self` when the last port of `other` is supplied with 1 W incident power. + +Parameters + +* self : *S_Matrix* +* other : *S_Matrix*
+S_Matrix instance* defined over the same frequency values of `self` and with a number of ports equal to (NP + 1) where NP is equal to the number of ports of `self`. The first NP ports of `other` must share the same port impedances of the ports of `self` + +Returns + +* p_incM : *numpy ndarray*
+ Nf 🞩 1 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined and NP is the number of ports of `self`. + PincM[i,0,k] represents the magnitude of the incident power at the k-port of `self` when the last port of `other` is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` +* phaseM : *numpy ndarray*
+ Nf 🞩 1 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined and NP is the number of ports of `self`. + phaseM[i,0,k] represents the 'phase' of the incident power at the k-port of `self` when the last port of `other` is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` + +#### `__findFreqIndex(self, freq)` + +Private method which returns the index of `self.frequencies` corresponding to the frequency `freq`. + +Parameters + +* self : *S_Matrix* +* freq : *float*
+ +Returns + +* idx : *int*
+Index of of `self.frequencies` corresponding to the frequency `freq`. In case `freq` is not in `self.frequencies`, the lower index of the nearest element of `self.frequencies` to `freq` is returned + +#### `fromZtoS(Z, freqs, z0 = 50)` + +*static method* used to return an *S_Matrix* starting from an impedance matrix defined over the frequency values listed in `freqs`. + +Parameters + +* Z : *numpy ndarray*
+Nf 🞩 NP 🞩 NP *numpy ndarray* of impedance parameters. Nf is the number of frequency values over which the impedance parameters are defined. NP is the number of ports +* freqs : *list*, *numpy ndarray*
+Nf *list* or *numpy ndarray* containing the frequency values, in hertz, over which the impedance parameters are defined +* z0 : *int*, *float*, *list* or *numpy ndarray*, *optional*
+port impedances in ohm. These can be given as a NP *list* or *numpy ndarray*. If all the ports share the same impedances, an *int* or *float* value can be passed as parameter. Default is 50 ohm + +Returns + +* S_Matrix : *S_Matrix*
+*S_Matrix* obtained from the conversion of the impedance parameters and defined over the frequency values listed in `freqs`. The port impedances of the returned *S_Matrix* are given by the relevant method parameter + +#### `fromYtoS(Z, freqs, y0 = 0.02)` + +*static method* used to return an *S_Matrix* starting from an admittance matrix defined over the frequency values listed in `freqs`. + +Parameters + +* Z : *numpy ndarray*
+Nf 🞩 NP 🞩 NP *numpy ndarray* of admittance parameters. Nf is the number of frequency values over which the admittance parameters are defined. NP is the number of ports +* freqs : *list*, *numpy ndarray*
+Nf *list* or *numpy ndarray* containing the frequency values, in hertz, over which the admittance parameters are defined +* y0 : *int*, *float*, *list* or *numpy ndarray*, *optional*
+port admittances. These can be given as a NP *list* or *numpy ndarray*. If all the ports share the same admittances, an *int* or *float* value can be passed as parameter. Default is 0.02 S + +Returns + +* S_Matrix : *S_Matrix*
+*S_Matrix* obtained from the conversion of the admittance parameters and defined over the frequency values listed in `freqs`. The port impedances of the returned *S_Matrix* are obtained from the port admittances given by the relevant method parameter + +#### `importTouchstone(filename, fmt='R_I', freqUnit='MHz', fix=False, n_f=None)` + +*static method* used to return an *S_Matrix* importing the scattering parameters from a *Touchstone* file. + +Parameters + +* filename : *string*
+*string* reporting the *touchstone* file name +* fmt : *string*, *optional*
+*string* defining the format of the scattering parameters reported in the *touchstone* file. The string can be equal to: + * 'R_I' : Real and Imaginary parts + * 'Mag_Deg' : Magnitude and phase in degrees + * 'dB_Deg' : Magnitude in decibel and phase in degrees + * 'Mag_Rad' : Magnitude and phase in radians + * 'dB_Rad' : Magnitude in decibel and phase in radians +Default is 'R_I' +* freqUnit : *string*, *optional*
+*string* defining the measurement unit of the frequency values reported in the *touchstone* file. It can be equal to 'Hz' or 'MHz'. Default is 'MHz' +* fix : *bool*, *optional*
+if `True` the method tries to fix the *touchstone* file where more than one row of scattering parameters corresponds to each single frequency value. If `True` the method returns a new *touchstone* file in the same folder of the original one and named as 'filename_fixed.*'. Default is `False` +* n_f : *int*, *optional*
+number of frequency values over which the scattering parameters reported in the *touchstone* file are defined. This parameter is mandatory if `fix=True` otherwise it can be `None` (default value) + +Returns + +* S_Matrix : *S_Matrix*
+*S_Matrix* obtained from the *touchstone* file + +#### `sMatrixRCseries(Rvalue, Cvalue, freqs, z0=50)` + +*static method* used to return an *S_Matrix* of a 1-port Resistance-Capacitance series circuit. + +Parameters + +* Rvalue : *float*
+The resistance value in ohm +* Cvalue : *float*
+The capacitance value in farad +* freqs : *list*, *numpy ndarray*
+*list*, *numpy ndarray* reporting the frequency values, in hertz, over which the returned *S_Matrix* is defined +* z0 : *float*, *optional*
+The impedance, in ohm, of the port of the returned *S_Matrix*. Default is 50 ohm + +Returns + +* S_Matrix : *S_Matrix*
+1-port *S_Matrix*, defined over the frequency values passed as method parameter, of a Resistance-Capacitance series. The port impedance is specified as method parameter + +#### `sMatrixRLseries(Rvalue, Lvalue, freqs, z0=50)` + +*static method* used to return an *S_Matrix* of a 1-port Resistance-Inductance series circuit. + +Parameters + +* Rvalue : *float*
+The resistance value in ohm +* Lvalue : *float*
+The inductance value in henry +* freqs : *list*, *numpy ndarray*
+*list*, *numpy ndarray* reporting the frequency values, in hertz, over which the returned *S_Matrix* is defined +* z0 : *float*, *optional*
+The impedance, in ohm, of the port of the returned *S_Matrix*. Default is 50 ohm + +Returns + +* S_Matrix : *S_Matrix*
+1-port *S_Matrix*, defined over the frequency values passed as method parameter, of a Resistance-Inductance series. The port impedance is specified as method parameter + +#### `sMatrixRCparallel(Rvalue, Cvalue, freqs, z0=50)` + +*static method* used to return an *S_Matrix* of a 1-port Resistance-Capacitance parallel circuit. + +Parameters + +* Rvalue : *float*
+The resistance value in ohm +* Cvalue : *float*
+The capacitance value in farad +* freqs : *list*, *numpy ndarray*
+*list*, *numpy ndarray* reporting the frequency values, in hertz, over which the returned *S_Matrix* is defined +* z0 : *float*, *optional*
+The impedance, in ohm, of the port of the returned *S_Matrix*. Default is 50 ohm + +Returns + +* S_Matrix : *S_Matrix*
+1-port *S_Matrix*, defined over the frequency values passed as method parameter, of a Resistance-Capacitance parallel. The port impedance is specified as method parameter + +#### `sMatrixRLparallel(Rvalue, Lvalue, freqs, z0=50)` + +*static method* used to return an *S_Matrix* of a 1-port Resistance-Inductance parallel circuit. + +Parameters + +* Rvalue : *float*
+The resistance value in ohm +* Lvalue : *float*
+The inductance value in henry +* freqs : *list*, *numpy ndarray*
+*list*, *numpy ndarray* reporting the frequency values, in hertz, over which the returned *S_Matrix* is defined +* z0 : *float*, *optional*
+The impedance, in ohm of the port of the returned *S_Matrix*. Default is 50 ohm + +Returns + +* S_Matrix : *S_Matrix*
+1-port *S_Matrix*, defined over the frequency values passed as method parameter, of a Resistance-Inductance parallel. The port impedance is specified as method parameter + +#### `sMatrixTnetwork(SL_1, SL_2, ST, z0=50)` + +*static method* used to return an *S_Matrix* of a 2-port T-network. + +Parameters + +* SL_1 : *S_Matrix*
+The 1-port S_Matrix of the longitudinal parameter of the T-network towards the first port +* SL_2 : *S_Matrix*
+The 1-port S_Matrix of the longitudinal parameter of the T-network towards the second port +* ST : *S_Matrix*
+The 1-port S_Matrix of the transversal parameter of the T-network +* z0 : *int*, *float*, *list* or *numpy ndarray*, *optional*
+port impedances in ohm. These can be given as a NP *list* or *numpy ndarray*. If all the ports share the same impedances, an *int* or *float* value can be passed as parameter. Default is 50 ohm + +![image](./images/sMatrixTnetwork.png) + +Returns + +* S_Matrix : *S_Matrix*
+2-port *S_Matrix*, of the T-network defined over the same frequency values over which SL_1, SL_2 and ST are defined. The impedances of the two ports are specified by the relevant method parameter + +#### `sMatrixPInetwork(ST_1, ST_2, SL, z0=50)` + +*static method* used to return an *S_Matrix* of a 2-port PI-network. + +Parameters + +* ST_1 : *S_Matrix*
+The 1-port S_Matrix of the transversal parameter of the PI-network towards the first port +* ST_2 : *S_Matrix*
+The 1-port S_Matrix of the transversal parameter of the PI-network towards the second port +* SL : *S_Matrix*
+The 1-port S_Matrix of the longitudinal parameter of the PI-network +* z0 : *int*, *float*, *list* or *numpy ndarray*, *optional*
+port impedances. These can be given as a NP *list* or *numpy ndarray*. If all the ports share the same impedances, an *int* or *float* value can be passed as parameter. Default is 50 ohm + +![image](./images/sMatrixPInetwork.png) + +Returns + +* S_Matrix : *S_Matrix*
+2-port *S_Matrix*, of the PI-network defined over the same frequency values over which ST_1, ST_2 and SL are defined. The impedances of the two ports are specified by the relevant method parameter + +#### `sMatrixTrLine(l, freqs, z0=50, c_f = 1, alpha=0)` + +*static method* used to return an *S_Matrix* of a low-losses transmission line. + +Parameters + +* l : *float*
+length, in meters, of the transmission line +* freqs : *list*, *numpy ndarray*
+*list*, *numpy ndarray* reporting the frequency values, in hertz, over which the returned *S_Matrix* is defined +* z0 : *int*, *float*, *optional*
+The characteristic impedance, in ohm, of the transmission line. Default is 50 ohm +* c_f : *float*, *optional*
+attenuating coefficient of the transmission line. Default is 0 + +Returns + +* S_Matrix : *S_Matrix*
+2-port *S_Matrix*, of the transmission line defined over the frequency values given as method parameter. The impedances of the two ports are equal to the characteristic impedance of the transmission line given as method parameter + +#### `__movePort(Smat, idx0, idx1)` + +private *static method* used to return an *S_Matrix* where the `Smat` port at index `idx0` is exchanged with the port at index `idx1`. + +Parameters + +* Smat : *S_Matrix*
+* idx0 : *int*
+start point index. It can be included in the range 0 to (NP-1) where NP is the number of ports of `Smat`. The index can be also included in the range -NP to -1. In the last case `idx0 += N_P` +* idx1 : *int*
+end point index. It can be included in the range 0 to (NP-1) where NP is the number of ports of `Smat`. The index can be also included in the range -NP to -1. In the last case `idx1 += N_P` + +Returns + +* S_Matrix : *S_Matrix*
+*S_Matrix* where the `Smat` port at index `idx0` is exchanged with the port at index `idx1`. If `idx0 < idx1` it moves the port at `idx0` of `Smat` after the port at `idx1` of `Smat`. If `idx0 > idx1` it moves the port at `idx0` of `Smat` before the port at `idx1` of `Smat`
+If:
+`Smat_t = movePort(Smat, idx0, idx1)`
+it is always verified:
+`Smat == movePort(Smat_t, idx1, idx0)`
+___ +## EM_Field class +This class is defined to manage harmonic electric and magnetic flux density fields generated by any multiport device. +The EM field must be computed on a regular cartesian grid (Δx = Δy = Δz). + +### Properties +* `EM_Field.b_field` : *numpy ndarray* or *None*
+Nf 🞩 NP 🞩 3 🞩 NN *numpy ndarray* where Nf is the number of frequency values over which the magnetic flux density field, expressed in tesla, is defined. NP is the number of ports of the device and NN is the number of spatial points over which the magnetic flux density field is defined. `EM_Field.b_field[i,k,j,n]` represents the rms of the `j` component of the magnetic flux density field at the frequency identified by the index `i`, when the port `k` is supplied with 1 W incident power in the point identified by index `n`. The values along the last axis of the array follow a fortran-like index ordering (first index changing faster) considering a Cartesian coordinate system. It can be `None` if the magnetic flux density field is not defined. +* `EM_Field.e_field` : *numpy ndarray* or *None*
+Nf 🞩 NP 🞩 3 🞩 NN *numpy ndarray* where Nf is the number of frequency values over which the electric field, expressed in volts per meter, is defined. NP is the number of ports of the device and NN is the number of spatial points over which the electric field is defined. `EM_Field.e_field[i,k,j,n]` represents the rms of the `j` component of the electric field at the frequency identified by the index `i`, when the port `k` is supplied with 1 W incident power in the point identified by index `n`. The values along the last axis of the array follow a fortran-like index ordering (first index changing faster) considering a Cartesian coordinate system. It can be `None` if the electric field is not defined. +* `EM_Field.frequencies` : *numpy ndarray*
+Nf *numpy ndarray* representing the frequency values, in hertz, over which electric and magnetic flux density fields are defined +* `EM_Field.nPoints` : *numpy ndarray*
+*numpy ndarray* with a length equal to three reporting the number of spatial points over which the electric and magnetic flux density fields are defined: [Nx, Ny, Nz] +* `EM_Field.n_f` : *int*
+number of frequency values over which the electric and magnetic flux density fields are defined (Nf) +* `EM_Field.nPorts` : *int*
+number of ports associated with the electric and magnetic flux density fields (NP) +* `EM_Field.properties` : *dict*
+*dictionary* which include additional properties related to the electric and magnetic flux density fields. These can be, for example, sample electric permittivity, electrical conductivity. Each additional property has to be defined over the same spatial points over which the electric and magnetic flux density fields are defined and must be represented by a *real* value. + +### Methods + +#### `__init__(self, freqs, nPoints, b_field=None, e_field=None, **kwargs)` + +It creates an *EM_Field* instance. At least one among b_field and e_field must be different from `None` + +Parameters + +* freqs : *list* or *numpy ndarray*
+Nf *list* or *numpy ndarray* representing the frequency values, in hertz, over which the electric and magnetic flux density fields are defined +* nPoints : *list* or *numpy ndarray*
+*list* or *numpy ndarray* with a length equal to three reporting the number of spatial points over which the electric and magnetic flux density fields are defined: [Nx, Ny, Nz] +* b_field : *numpy ndarray*, *optional*
+Nf 🞩 NP 🞩 3 🞩 NN *numpy ndarray* where Nf is the number of frequency values over which the magnetic flux density field, expressed in tesla, is defined. NP is the number of ports of the device and NN is the number of spatial points over which the magnetic flux density field is defined. `b_field[i,k,j,n]` represents the rms of the `j` component of the magnetic flux density field at the frequency identified by the index `i`, when the port `k` is supplied with 1 W incident power in the point identified by index `n`. The values along the last axis of the array follow a fortran-like index ordering (first index changing faster) considering a Cartesian coordinate system. It can be `None` (default value) if the magnetic flux density field is not defined. +* e_field : *numpy ndarray*, *optional*
+Nf 🞩 NP 🞩 3 🞩 NN *numpy ndarray* where Nf is the number of frequency values over which the electric field, expressed in volts per meter, is defined. NP is the number of ports of the device and NN is the number of spatial points over which the electric field is defined. `EM_Field.e_field[i,k,j,n]` represents the rms of the `j` component of the electric field at the frequency identified by the index `i`, when the port `k` is supplied with 1 W incident power in the point identified by index `n`. The values along the last axis of the array follow a fortran-like index ordering (first index changing faster) considering a Cartesian coordinate system. It can be `None` (default value) if the electric field is not defined. +* **kwargs : *list*, *numpy ndarray*, *optional*
+additional properties defined over the same NN spatial points over which the electric and magnetic flux density fields are defined. They can be *list* or *numpy ndarray* with a length equal to NN + +Returns + +* EM_Field : *EM_Field* + +Example +```python +import numpy as np +from cosimpy import * + +n_p = 5 #Number of ports +n_f = 10 #Number of frequency values +nPoints = [3,4,5] #Number of points along the three Cartesian directions + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +bfield = np.random.rand(n_f,n_p,3,np.prod(nPoints)) + 1j*np.random.rand(n_f,n_p,3,np.prod(nPoints)) + +cond = np.random.uniform(size=np.prod(nPoints)) + +EM_Field(frequencies, nPoints, b_field=bfield, conductivity=cond) + +''' +Out: + + """"""""""""""" + EM FIELD + """"""""""""""" + + Number of frequency values = 10 + Number of ports = 5 + Number of point (nx, ny, nz) = 3, 4, 5 + + E field not defined + + 'conductivity' additional property defined + +''' +``` + +#### `__repr__(self)` + +Method for returning a printable representation of the *EM_Field* instance. + +Parameters + +* self : *EM_Field* + +Returns + +* string : *string*
+The string identifies the class of the instance. It reports the number of frequency values, the number of ports and the number of spatial points over which the EM field is defined + +#### `compSensitivities(self)` + +If the magnetic flux density field is defined, the method returns the complex B1+ and B1- values. + +Parameters + +* self : *EM_Field* + +Returns + +* sens : *numpy ndarray*
+Nf 🞩 NP 🞩 2 🞩 NN *numpy ndarray* where Nf is the number of frequency values over which the sensitivity, expressed in tesla, are defined. NP is the number of ports of the device and NN is the number of spatial points over which the sensitivities are defined. `sens[i,k,j,n]` represents the B1+ (if `j` is equal to 0) or B1- (if `j` is equal to 1) complex value at the frequency identified by the index `i`, generated by 1 W incident power at port `k` in the point identified by index `n`. The values along the last axis of the array follow a fortran-like index ordering (first index changing faster) considering a Cartesian coordinate system + +#### `plotB(self, comp, freq, port, plane, sliceIdx, vmin=None, vmax=None)` + +Method for plotting the magnetic flux density field over a specific slice + +![image](./images/plotB.png) + +Parameters + +* self : *EM_Field* +* comp : *string*
+component of the magnetic flux density field, in microtesla, to be plotted. It can take the following values: + * 'x' : x-component + * 'y' : y-component + * 'z' : z-component + * 'mag' : magnitude of the magnetic flux density field vector (rms value) + * 'b1+' : magnitude of B1+ + * 'b1-' : magnitude of B1- +* freq : *float*
+frequency, in hertz, of the magnetic flux density field to be plotted +* port : *int*
+port number which is supplied with 1 W incident power to generate the magnetic flux density field. All the other ports are closed on matched loads. It can take values from 1 to NP being NP the number of ports of the device +* plane : *string*
+slice to be plotted. It can take the following values: + * 'xy' : axial slice + * 'xz' : coronal slice + * 'yz' : sagittal slice +* sliceIdx : *int*
+index of the slice to be plotted. It can take negative values. In that case `sliceIdx = Ns + sliceIdx` where Ns is the total number of slices in the relevant direction +* vmin *float*, *optional*
+minimum value, in microtesla, of the chromatic bar. Default is `None` +* vmax *float*, *optional*
+maximum value, in microtesla, of the chromatic bar. Default is `None` + +Returns + +* None + +#### `plotE(self, comp, freq, port, plane, sliceIdx, vmin=None, vmax=None)` + +Method for plotting the melectric field over a specific slice + +![image](./images/plotE_n.png) + +Parameters + +* self : *EM_Field* +* comp : *string*
+component of the electric field, in volts per meter, to be plotted. It can take the following values: + * 'x' : x-component + * 'y' : y-component + * 'z' : z-component + * 'mag' : magnitude of the electric field vector (rms value) +* freq : *float*
+frequency, in hertz, of the electric field to be plotted +* port : *int*
+port number which is supplied with 1 W incident power to generate the electric field. All the other ports are closed on matched loads. It can take values from 1 to NP being NP the number of ports of the device +* plane : *string*
+slice to be plotted. It can take the following values: + * 'xy' : axial slice + * 'xz' : coronal slice + * 'yz' : sagittal slice +* sliceIdx : *int*
+index of the slice to be plotted. It can take negative values. In that case `sliceIdx = Ns + sliceIdx` where Ns is the total number of slices in the relevant direction +* vmin *float*, *optional*
+minimum value, in volts per meter, of the chromatic bar. Default is `None` +* vmax *float*, *optional*
+maximum value, in volts per meter, of the chromatic bar Default is `None` + +Returns + +* None + +#### `exportXMF(self, filename)` + +Method for generating a .xmf and a .h5 file to be imported in data analysis and visualization applications such as [ParaView](https://www.paraview.org/). In the .h5 files, the electric field and magnetic flux density field are stored as MHz-p-_ where is the frequency value in megahertz, is the port number, is 'E' (electric field) or 'B' (magnetic flux density field) and can be 'real' or 'imag' identifying the real or imaginary part of the field phasor. Additional properties, defined along the *EM_Field* instance, are stored with their full name in the .h5 file. + +![image](./images/ParaView.png) + +Parameters + +* self : *EM_Field* +* filename : *string*
+name of the .xmf and .h5 files + +Returns + +* None + +#### `_newFieldComp(self, p_inc, phase)` + +It returns an *EM_Field* instance where the EM fields are updated according to new incident powers whose magnitude and 'phase' are contained in `p_inc` and `phase` *numpy ndarray* + +Parameters + +* self : *EM_Field* +* p_incM : *numpy ndarray*
+Nf 🞩 NPout 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined, NPout is the number of ports of the returned *EM_Field* instance and NP is the number of ports of `self`. +PincM[i,j,k] represents the magnitude of the incident power at the k-port of `self` when the j-port of the output *EM_Field* is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` +* phaseM : *numpy ndarray*
+Nf 🞩 NPout 🞩 NP *numpy ndarray* where Nf is the number of frequency values over which `self` is defined, NPout is the number of ports of the returned *EM_Field* instance and NP is the number of ports of `self`. +PincM[i,j,k] represents the magnitude of the incident power at the k-port of `self` when the j-port of the output *EM_Field* is supplied with 1 W incident power at a frequency equal to `self.frequencies[i]` + +Returns + +* EM_Field : *EM_Field*
+*EM_Field* defined over the same frequency values of `self` *EM_Field*. The returned *EM_Field* is characterised by NPout different distributions of EM fields corresponding to the number of rows of `p_incM` and `phaseM` method parameters + +#### `importFields_cst(directory, freqs, nPorts, nPoints=None, Pinc_ref=1, b_multCoeff=1, pkORrms='pk', imp_efield=True, imp_bfield=True, **kwargs)` + +Static method which returns an *EM_Field* instance importing the data from [CST® STUDIO SUITE](https://www.3ds.com/products-services/simulia/products/cst-studio-suite/) standard ASCII export files whose columns are formatted as: + +|x, y, z, xRe, yRe, zRe, xIm, yIm, zIm| + +Results files must be collected in a dedicated directory and named as __port.txt where can be either 'efield', for the electric field result, or 'bfield' for magnetic flux density field results, is a string identifying the frequency value in megahertz (e.g. '128.4') and is the number of the port supplied, in the simulation environment, to generate the relevant EM fields. All the EM quantities must be exported on the same regular grid. + +Parameters + +* directory : *string*
+*string* reporting the path of the directory which contains the text files +* freqs : *list*, *numpy ndarray*
+*list* or *numpy ndarray* of *string* values reporting the frequency values as they are indicated in the results filenames () +* nPorts : *int*
+number of ports for which the results are available +* nPoints : *list*, *numpy ndarray*, *optional*
+*list* or *numpy ndarray* with a length equal to three reporting the number of spatial points over which the electric and magnetic flux density fields are defined: [Nx, Ny, Nz]. If `None` (default value) the method deduces them from the result files at the expenses of a slight longer import process +* Pinc_ref : *float*, *optional*
+magnitude of the power incident at the ports which generate the EM fields. Default is 1 W. +* b_multCoeff : *float*, *optional*
+multiplicative factor for the magnetic results. In case the exported results are magnetic field values (A/m), it has to be equal to 4π10-7 to obtain magnetic flux density values in tesla. Default is 1 +* pkORrms : *string*, *optional*
+if 'pk' the values contained in the result files are interpreted as peak values, if 'rms' they are interpreted as rms values. Default is 'pk' +* imp_efield : *bool*, *optional*
+if `True` the method imports the results relevant to the electric field. Default is `True` +* imp_bfield : *bool*, *optional*
+if `True` the method imports the results relevant to the magnetic flux density field. Default is `True` +* **kwargs : *list*, *numpy ndarray*, *optional*
+additional properties defined over the same NN spatial points over which the electric and magnetic flux density fields are defined. They can be *list* or *numpy ndarray* with a length equal to NN + +Returns + +* EM_Field : *EM_Field*
+*EM_Field* obtained from the EM results exported from CST® STUDIO SUITE. `EM_Field.properties` is a dictionary based on the `**kwargs`parameters passed to the method + +#### `importFields_s4l(directory, freqs, nPorts, Pinc_ref=1, b_multCoeff=1, pkORrms='rms', imp_efield=True, imp_bfield=True, **kwargs)` + +Static method which returns an *EM_Field* instance importing the data from [Sim4Life](https://zmt.swiss/sim4life/) standard export .mat files. Results files must be collected in a dedicated directory and named as _port.mat where can be either 'efield', for the electric field results, or 'bfield' for magnetic flux density field results and is the number of the ports supplied, in the simulation environment, to generate the relevant EM field. All the EM quantities must be exported on the same regular grid. + +Parameters + +* directory : *string*
+*string* reporting the path of the directory which contains the .mat files +* freqs : *list*, *numpy ndarray*
+*list* or *numpy ndarray* of *float* values reporting the sorted frequency values at which the results have been exported from Sim4Life +* nPorts : *int*
+number of ports for which the results are available +* nPoints : *list*, *numpy ndarray*, *optional*
+*list* or *numpy ndarray* with a length equal to three reporting the number of spatial points over which the electric and magnetic flux density fields are defined: [Nx, Ny, Nz]. If `None` (default value) the method deduces them from the result files at the expenses of a slight longer import process +* Pinc_ref : *float*, *optional*
+magnitude of the power incident at the ports which generate the EM fields. Default is 1 W. +* b_multCoeff : *float*, *optional*
+multiplicative factor for the magnetic results. In case the exported results are magnetic field values (A/m), it has to be equal to 4π10-7 to obtain magnetic flux density values in tesla. Default is 1 +* pkORrms : *string*, *optional*
+if 'pk' the values contained in the result files are interpreted as peak values, if 'rms' they are interpreted as rms values. Default is 'rms' +* imp_efield : *bool*, *optional*
+if `True` the method imports the results relevant to the electric field. Default is `True` +* imp_bfield : *bool*, *optional*
+if `True` the method imports the results relevant to the magnetic flux density field. Default is `True` +* **kwargs : *list*, *numpy ndarray*, *optional*
+additional properties defined over the same NN spatial points over which the electric and magnetic flux density fields are defined. They can be *list* or *numpy ndarray* with a length equal to NN + +Returns + +* EM_Field : *EM_Field*
+*EM_Field* obtained from the EM results exported from Sim4Life. `EM_Field.properties` is a dictionary based on the `**kwargs` parameters passed to the method + +___ +## RF_Coil class +This class represents the link between the *S_Matrix* and *EM_Field* classes. An instance of this class represents the simulated device being, its properties, an *S_Matrix* and an *EM_Field* instance. The *S_Matrix* must be defined over all the frequency values at which the *EM_Field* is defined. + +### Properties +* `RF_Coil.s_matrix` : *S_Matrix*
+*S_Matrix* instance defined over Nf,s frequency values and with NP ports +* `RF_Coil.em_field` : *EM_Field*
+*EM_Field* instance defined over Nf,em frequency values and with NP ports + +### Methods + +#### `__init__(self, s_matrix, em_field=None)` + +It creates an *RF_Coil* instance. + +Parameters + +* s_matrix : *S_Matrix*
+*S_Matrix* instance defined over Nf,s frequency values and with NP ports +* em_field : *EM_Field*, optional
+*EM_Field* instance defined over Nf,em frequency values and with NP ports. Default is `None` + +Returns + +* RF_Coil : *RF_Coil* + +Example +```python +import numpy as np +from cosimpy import * + +n_p = 5 #Number of ports +n_f = 10 #Number of frequency values +nPoints = [3,4,5] #Number of points for EM field along the three Cartesian directions + +frequencies = np.linspace(50e6, 150e6, n_f, endpoint=False)) + +s_matrix = S_Matrix(np.random.uniform(size=(n_f,n_p,n_p)), freqs=frequencies) + +bfield = np.random.rand(n_f,n_p,3,np.prod(nPoints)) + +em_field = EM_Field(frequencies, nPoints, b_field=bfield) + +''' +Out: + + """"""""""""""" + RF COIL + """"""""""""""" + + + +""""""""""""""" + S MATRIX +""""""""""""""" + +|V-| = |S||V+| +|5 x 1| = |5 x 5||5 x 1| + +Number of frequency values = 10 + +""""""""""""""" + EM FIELD +""""""""""""""" + +Number of frequency values = 10 + +Number of points = (3, 4, 5) + +E field not defined + +''' +``` + +#### `__repr__(self)` + +Method for returning a printable representation of the *RF_Coil* instance. + +Parameters + +* self : *RF_Coil* + +Returns + +* string : *string*
+The string identifies the class of the instance. It reports information about the *S_Matrix* and *EM_Field* properties + +#### `singlePortConnRFcoil(self, networks, comp_Pinc=False)` + +The method performs the cascade connection between `self` and the *S_Matrix* instances contained in `networks` *list*. + +Parameters + +* self : *RF_Coil* +* networks : *list*
+*list* with length equal to the number of port of `self.s_matrix` either containing `None` values or *S_Matrix* instances defined over the same frequency values of `self.s_matrix`. If the n-th element of `networks` is `None`, the n-th port of `self` is not connected to any other network in the returned *RF_Coil*. If the n-th element of `networks` is an N-ports *S_Matrix*, its first port must share the same impedance of the n-th port of `self` and the two ports are connected together. In the returned *RF_Coil*, the n-th port of `self` is therefore expanded into (N-ports - 1) ports. The example proposed in the image below should clarify the workflow + +![image](./images/singlePortConnSMatrix.png) + +* comp_Pinc : *bool*, *optional*
+if `True` the method returns an *RF_Coil* instance whose `em_field` property is updated according to the new connections. If `False` (default value), the `em_field` property of the returned *RF_Coil* is `None` + +Returns + +* RF_Coil : *RF_Coil*
+The new *RF_Coil* resulting form the connection of the original `self` *RF_Coil* with the *S_Matrix* instances contained in `networks`. It is worth noting that, if any of the network in the `networks` *list*, generates, under the given supplying conditions, an EM field which is not negligible with respect to that generated by the unconnected RF coil, the EM field associated with the returned *RF_Coil* can be inaccurate + +#### `fullPortsConnRFcoil(self, s_matrix_other, comp_Pinc=False)` + +The method performs the cascade connection between `self` and `other` *S_Matrix* instance. + +Parameters + +* self : *RF_Coil* +* other : *S_Matrix*
+*S_Matrix* instance defined over the same frequency values of `self.s_matrix` and with a number of ports higher than those of `self`. In the returned *RF_Coil*, the first NP ports of `other` are connected with the NP ports of `self` and must share the same port impedances +* comp_Pinc : *bool*, *optional*
+if `True` the method returns an *RF_Coil* instance whose `em_field` property is updated according to the new connections. If `False` (default value), the `em_field` property of the returned *RF_Coil* is `None` + +Returns + +* RF_Coil : *RF_Coil*
+The new *RF_Coil* resulting form the connection of the original `self` *RF_Coil* with the *S_Matrix* `other`. It is worth noting that, if the network, described by the *S_Matrix* `other`, generates, under the given supplying conditions, an EM field which is not negligible with respect to that generated by the unconnected RF coil, the EM field associated with the returned *RF_Coil* can be inaccurate \ No newline at end of file diff --git a/docs/images/EMPIR_logo.jpg b/docs/images/EMPIR_logo.jpg new file mode 100644 index 0000000..27cde6a Binary files /dev/null and b/docs/images/EMPIR_logo.jpg differ diff --git a/docs/images/MIMAS_logo.png b/docs/images/MIMAS_logo.png new file mode 100644 index 0000000..5d9ddad Binary files /dev/null and b/docs/images/MIMAS_logo.png differ diff --git a/docs/images/ParaView.png b/docs/images/ParaView.png new file mode 100644 index 0000000..79f9aa2 Binary files /dev/null and b/docs/images/ParaView.png differ diff --git a/docs/images/example_S.png b/docs/images/example_S.png new file mode 100644 index 0000000..ac9a9bd Binary files /dev/null and b/docs/images/example_S.png differ diff --git a/docs/images/plotB.png b/docs/images/plotB.png new file mode 100644 index 0000000..3910969 Binary files /dev/null and b/docs/images/plotB.png differ diff --git a/docs/images/plotE.png b/docs/images/plotE.png new file mode 100644 index 0000000..9de54c7 Binary files /dev/null and b/docs/images/plotE.png differ diff --git a/docs/images/plotE_n.png b/docs/images/plotE_n.png new file mode 100644 index 0000000..007266f Binary files /dev/null and b/docs/images/plotE_n.png differ diff --git a/docs/images/sMatrixPInetwork.png b/docs/images/sMatrixPInetwork.png new file mode 100644 index 0000000..accd41b Binary files /dev/null and b/docs/images/sMatrixPInetwork.png differ diff --git a/docs/images/sMatrixTnetwork.png b/docs/images/sMatrixTnetwork.png new file mode 100644 index 0000000..fe53c54 Binary files /dev/null and b/docs/images/sMatrixTnetwork.png differ diff --git a/docs/images/singlePortConnSMatrix.png b/docs/images/singlePortConnSMatrix.png new file mode 100644 index 0000000..ccbedc6 Binary files /dev/null and b/docs/images/singlePortConnSMatrix.png differ diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..30250b1 --- /dev/null +++ b/setup.py @@ -0,0 +1,6 @@ +import setuptools + +with open(".//docs//Documentation.md", "r", encoding="utf-8") as fh: + long_description = fh.read() + +setuptools.setup(name="cosimpy", version="1.0.0", licence="MIT", url='https://github.com/umbertozanovello/CoSimPy', packages=setuptools.find_packages(), author="Umberto Zanovello", description="Python electromagnetic co-simulation library", long_description=long_description, long_description_content_type="text/markdown", classifiers=["Programming Language :: Python :: 3", "License :: OSI Approved :: MIT License", "Operating System :: OS Independent",], python_requires='>=3.5',)