Skip to content

Commit

Permalink
Merge pull request #3 from umbertozanovello/Develop
Browse files Browse the repository at this point in the history 
Develop
  • Loading branch information
@umbertozanovello
umbertozanovello authored Sep 6, 2021
2 parents c461adf + 3ee748d commit 658684b
Show file tree
Hide file tree
Showing 7 changed files with 535 additions and 73 deletions.
73 changes: 61 additions & 12 deletions cosimpy/EM_Field.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,16 +10,16 @@ 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")
raise TypeError("nPoints can only be a list or numpy ndarray with length equal to 3")
elif len(nPoints) != 3:
raise TypeError("nPoints can only be a list with length equal to 3")
raise TypeError("nPoints can only be a list or numpy ndarray 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")
raise TypeError("e_field can only be numpy 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):
Expand All @@ -30,7 +30,7 @@ def __init__(self, freqs, nPoints, b_field=None, e_field=None, **kwargs):
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")
raise TypeError("b_field can only be numpy 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):
Expand All @@ -42,10 +42,10 @@ def __init__(self, freqs, nPoints, b_field=None, e_field=None, **kwargs):

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("All the additional properties have to be list or numpy ndarray of float or int")
elif np.array(arg).dtype != np.dtype(np.float) and np.array(arg).dtype != np.dtype(np.int):
raise ValueError("All the additional properties have to be list or numpy ndarray of float or int")
elif np.array(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:
Expand Down Expand Up @@ -116,7 +116,7 @@ def __repr__(self):
def compSensitivities(self):

if self.__b_field is None:
raise ValueError("No b field property is specifed for the EM_Field instance")
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)
Expand All @@ -130,10 +130,59 @@ def compSensitivities(self):
return sens


def compPowDens(self, elCond=None, p_inc=None):

if elCond is None and not "elCond" in self.__prop.keys():
raise ValueError("No 'elCond' key is found in self.properties. Please, provide the electrical conductivity as argument of the method")
elif elCond is not None:
if not isinstance(elCond, list) and not isinstance(elCond, np.ndarray):
raise ValueError("elCond has to be a list or a numpy ndarray of float or int")
elif np.array(elCond).dtype != np.dtype(np.float) and np.array(elCond).dtype != np.dtype(np.int):
raise ValueError("elCond has to be a list or an numpy ndarray of float or int")
elif np.array(elCond).size != np.prod(self.__nPoints):
raise ValueError("elCond has a length different from nPoints[0]*nPoints[1]*nPoints[2]")
else:
elCond = self.__prop["elCond"]

if self.__e_field is None:
raise ValueError("No e field property is specifed for the EM_Field instance. Power density cannot be computed")

if p_inc is not None: #Power density is computed for a defined supply configuration
if not isinstance(p_inc, np.ndarray) and not isinstance(p_inc, list):
raise TypeError("S matrix can only be numpy ndarray or a list")
else:
p_inc = np.array(p_inc)
if p_inc.size != self.__nPorts:
raise TypeError("p_inc has to be a self.nPorts length list or numpy ndarray")

norm = np.sqrt(np.abs(p_inc))*np.exp(1j*np.angle(p_inc))
efield_new = np.moveaxis(self.e_field,1,-1) #Temporary axis change so field.shape = [self.n_f, 3, self.nPoints, self.nPorts]
efield_new = efield_new @ norm
else:
efield_new = self.__e_field

powDens = np.linalg.norm(efield_new, axis=-2)**2 * elCond

return powDens


def compDepPow(self, voxVols, elCond=None, p_inc=None):


if not isinstance(voxVols, int) and not isinstance(voxVols, float):
raise ValueError("voxVols has to be int or float")

powDens = self.compPowDens(elCond, p_inc)

depPow = np.nansum(powDens*voxVols, axis=-1)

return depPow


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")
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:
Expand Down Expand Up @@ -189,7 +238,7 @@ def plotB(self, comp, freq, port, plane, sliceIdx, vmin=None, vmax=None):
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")
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:
Expand Down Expand Up @@ -531,4 +580,4 @@ def importFields_s4l(directory, freqs, nPorts, Pinc_ref=1, b_multCoeff=1, pkORrm
b_field = b_multCoeff * np.sqrt(1/Pinc_ref) * rmsCoeff * b_field


return EM_Field(freqs, nPoints, b_field, e_field, **kwargs)
return EM_Field(freqs, nPoints, b_field, e_field, **kwargs)
79 changes: 67 additions & 12 deletions cosimpy/RF_Coil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,37 +49,92 @@ def __repr__(self):

def singlePortConnRFcoil(self, networks, comp_Pinc=False):

rets = self.s_matrix._singlePortConnSMatrix(networks, comp_Pinc)

if not comp_Pinc or self.em_field is None:
#If comp_Pinc is True but the em_field property is None,
#the S_Matrix method returns the Pinc and phase matrices. These will be discared

S_l = self.s_matrix._singlePortConnSMatrix(networks, False)

return RF_Coil(S_l, None)
return RF_Coil(rets[0], 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]
idxs = np.where(self.s_matrix.frequencies[:,None]==self.em_field.frequencies)[0]
p_inc, phase = rets[1:]
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)
return RF_Coil(rets[0], em_field_new)


def fullPortsConnRFcoil(self, s_matrix_other, comp_Pinc=False):

rets = self.s_matrix._fullPortsConnSMatrix(s_matrix_other, comp_Pinc)

if not comp_Pinc or self.em_field is None:
#If comp_Pinc is True but the em_field property is None,
#the S_Matrix method returns the Pinc and phase matrices. These will be discared

S_l = self.s_matrix._fullPortsConnSMatrix(s_matrix_other, False)

return RF_Coil(S_l, None)
return RF_Coil(rets[0], 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]
idxs = np.where(self.s_matrix.frequencies[:,None]==self.em_field.frequencies)[0]
p_inc, phase = rets[1:]
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)
return RF_Coil(rets[0], em_field_new)


def powerBalance(self, p_inc, voxVols=None, elCond=None, printReport=False):

powBal = {}
tot_p_inc = np.sum(np.abs(p_inc))
powBal["P_inc_tot"] = tot_p_inc

pows = self.s_matrix.powerBalance(p_inc)

if self.em_field is None or self.em_field.e_field is None:

frequencies = self.s_matrix.frequencies

if len(pows) == 1: # No available data for external circuitries losses
powBal["P_refl"] = tot_p_inc - pows[0]
powBal["P_other"] = pows[0]
elif len(pows) == 2:
powBal["P_refl"] = tot_p_inc - pows[0]
powBal["P_circ_loss"] = pows[0] - pows[1]
powBal["P_other"] = tot_p_inc - powBal["P_refl"] - powBal["P_circ_loss"]

else:

idxs = np.where(self.s_matrix.frequencies[:,None]==self.em_field.frequencies)[0]

depPow = self.em_field.compDepPow(voxVols, elCond, p_inc)

if len(pows) == 1: # No available data for external circuitries losses
powBal["P_refl"] = tot_p_inc - pows[0][idxs]
powBal["P_dep"] = depPow
powBal["P_other"] = pows[0][idxs] - depPow
elif len(pows) == 2:
powBal["P_refl"] = tot_p_inc - pows[0][idxs]
powBal["P_dep"] = depPow
powBal["P_circ_loss"] = pows[0] - pows[1]
powBal["P_other"] = tot_p_inc - powBal["P_refl"] - powBal["P_circ_loss"] - powBal["P_dep"]

frequencies = self.s_matrix.frequencies[idxs]

if printReport:
print("\n\nTotal incident power: %.2e W" %(tot_p_inc))
for i, freq in enumerate(frequencies):
print("\n\n"+"-"*10+"\n")
print("Frequency: %.2e Hz" %(freq))
print("\n"+"-"*10+"\n\n")

for key in powBal:
if key != "P_inc_tot":
print("%s: %.2e W\t(%.2f %%)\n"%(key, powBal[key][i], powBal[key][i]/tot_p_inc*100))
return powBal
Loading

0 comments on commit 658684b

Please sign in to comment.