@@ -3835,92 +3835,91 @@ def __mutual_inductance(
38353835 """
38363836
38373837 jj = np .r_ [ii + 1 : self .total_elements_current_carriers ]
3838- # ll = lmod[jj]
3839- # mm = lmod[ii]
3840- # len_jj = self.total_elements_current_carriers - (ii + 1)
3841- # rr = dict(
3842- # end_end=np.zeros(len_jj),
3843- # end_start=np.zeros(len_jj),
3844- # start_start=np.zeros(len_jj),
3845- # start_end=np.zeros(len_jj),
3846- # )
3838+ ll = lmod [jj ]
3839+ mm = lmod [ii ]
3840+ len_jj = self .total_elements_current_carriers - (ii + 1 )
3841+ rr = dict (
3842+ end_end = np .zeros (len_jj ),
3843+ end_start = np .zeros (len_jj ),
3844+ start_start = np .zeros (len_jj ),
3845+ start_end = np .zeros (len_jj ),
3846+ )
38473847
3848- # for key in rr.keys():
3849- # rr[key] = self.__vertex_to_vertex_distance(key, ii, jj)
3850- # # End for key
3848+ for key in rr .keys ():
3849+ rr [key ] = self .__vertex_to_vertex_distance (key , ii , jj )
3850+ # End for key
38513851
3852- # # Additional parameters
3853- # alpha2 = (
3854- # rr["start_end"] ** 2
3855- # - rr["start_start"] ** 2
3856- # + rr["end_start"] ** 2
3857- # - rr["end_end"] ** 2
3858- # )
3852+ # Additional parameters
3853+ alpha2 = (
3854+ rr ["start_end" ] ** 2
3855+ - rr ["start_start" ] ** 2
3856+ + rr ["end_start" ] ** 2
3857+ - rr ["end_end" ] ** 2
3858+ )
38593859
3860- # cos_eps = np.minimum(np.maximum(alpha2 / (2 * ll * mm), -1.0), 1.0)
3861- # sin_eps = np.sin(np.arccos(cos_eps))
3862-
3863- # dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2
3864- # mu = (
3865- # ll
3866- # * (
3867- # 2 * mm ** 2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2)
3868- # + alpha2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2)
3869- # )
3870- # / dd
3871- # )
3872- # nu = (
3873- # mm
3874- # * (
3875- # 2 * ll ** 2 * (rr["start_end"] ** 2 - rr["start_start"] ** 2 - mm ** 2)
3876- # + alpha2 * (rr["end_start"] ** 2 - rr["start_start"] ** 2 - ll ** 2)
3877- # )
3878- # / dd
3879- # )
3880- # d2 = rr["start_start"] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps
3881-
3882- # # avoid rounding for segments in a plane
3883- # d2[d2 < abstol ** 2] = 0
3884- # d0 = np.sqrt(d2)
3885-
3886- # # solid angles
3887- # omega = (
3888- # np.arctan(
3889- # (d2 * cos_eps + (mu + ll) * (nu + mm) * sin_eps ** 2)
3890- # / (d0 * rr["end_end"] * sin_eps)
3891- # )
3892- # - np.arctan(
3893- # (d2 * cos_eps + (mu + ll) * nu * sin_eps ** 2)
3894- # / (d0 * rr["end_start"] * sin_eps)
3895- # )
3896- # + np.arctan(
3897- # (d2 * cos_eps + mu * nu * sin_eps ** 2)
3898- # / (d0 * rr["start_start"] * sin_eps)
3899- # )
3900- # - np.arctan(
3901- # (d2 * cos_eps + mu * (nu + mm) * sin_eps ** 2)
3902- # / (d0 * rr["start_end"] * sin_eps)
3903- # )
3904- # )
3905- # omega[d0 == 0.0] = 0.0
3860+ cos_eps = np .minimum (np .maximum (alpha2 / (2 * ll * mm ), - 1.0 ), 1.0 )
3861+ sin_eps = np .sin (np .arccos (cos_eps ))
39063862
3907- # # contribution
3908- # pp = np.zeros((len_jj, 5), dtype=float)
3909- # pp[:, 0] = (ll + mu) * np.arctanh(mm / (rr["end_end"] + rr["end_start"]))
3910- # pp[:, 1] = -nu * np.arctanh(ll / (rr["end_start"] + rr["start_start"]))
3911- # pp[:, 2] = (mm + nu) * np.arctanh(ll / (rr["end_end"] + rr["start_end"]))
3912- # pp[:, 3] = -mu * np.arctanh(mm / (rr["start_start"] + rr["start_end"]))
3913- # pp[:, 4] = d0 * omega / sin_eps
3863+ dd = 4 * ll ** 2 * mm ** 2 - alpha2 ** 2
3864+ mu = (
3865+ ll
3866+ * (
3867+ 2 * mm ** 2 * (rr ["end_start" ] ** 2 - rr ["start_start" ] ** 2 - ll ** 2 )
3868+ + alpha2 * (rr ["start_end" ] ** 2 - rr ["start_start" ] ** 2 - mm ** 2 )
3869+ )
3870+ / dd
3871+ )
3872+ nu = (
3873+ mm
3874+ * (
3875+ 2 * ll ** 2 * (rr ["start_end" ] ** 2 - rr ["start_start" ] ** 2 - mm ** 2 )
3876+ + alpha2 * (rr ["end_start" ] ** 2 - rr ["start_start" ] ** 2 - ll ** 2 )
3877+ )
3878+ / dd
3879+ )
3880+ d2 = rr ["start_start" ] ** 2 - mu ** 2 - nu ** 2 + 2 * mu * nu * cos_eps
39143881
3915- # # filter odd cases (e.g. consecutive segments)
3916- # pp[np.isnan(pp)] = 0.0
3917- # pp[np.isinf(pp)] = 0.0
3882+ # avoid rounding for segments in a plane
3883+ d2 [d2 < abstol ** 2 ] = 0
3884+ d0 = np .sqrt (d2 )
3885+
3886+ # solid angles
3887+ omega = (
3888+ np .arctan (
3889+ (d2 * cos_eps + (mu + ll ) * (nu + mm ) * sin_eps ** 2 )
3890+ / (d0 * rr ["end_end" ] * sin_eps )
3891+ )
3892+ - np .arctan (
3893+ (d2 * cos_eps + (mu + ll ) * nu * sin_eps ** 2 )
3894+ / (d0 * rr ["end_start" ] * sin_eps )
3895+ )
3896+ + np .arctan (
3897+ (d2 * cos_eps + mu * nu * sin_eps ** 2 )
3898+ / (d0 * rr ["start_start" ] * sin_eps )
3899+ )
3900+ - np .arctan (
3901+ (d2 * cos_eps + mu * (nu + mm ) * sin_eps ** 2 )
3902+ / (d0 * rr ["start_end" ] * sin_eps )
3903+ )
3904+ )
3905+ omega [d0 == 0.0 ] = 0.0
3906+
3907+ # contribution
3908+ pp = np .zeros ((len_jj , 5 ), dtype = float )
3909+ pp [:, 0 ] = (ll + mu ) * np .arctanh (mm / (rr ["end_end" ] + rr ["end_start" ]))
3910+ pp [:, 1 ] = - nu * np .arctanh (ll / (rr ["end_start" ] + rr ["start_start" ]))
3911+ pp [:, 2 ] = (mm + nu ) * np .arctanh (ll / (rr ["end_end" ] + rr ["start_end" ]))
3912+ pp [:, 3 ] = - mu * np .arctanh (mm / (rr ["start_start" ] + rr ["start_end" ]))
3913+ pp [:, 4 ] = d0 * omega / sin_eps
3914+
3915+ # filter odd cases (e.g. consecutive segments)
3916+ pp [np .isnan (pp )] = 0.0
3917+ pp [np .isinf (pp )] = 0.0
39183918
39193919 # Mutual inductances
39203920 matrix [ii , jj ] = (
3921- 5.0e-8 # Mutual inductance imposed by Zappatore input file
3922- # 2 * cos_eps * (pp[:, 0] + pp[:, 1] + pp[:, 2] + pp[:, 3])
3923- # - cos_eps * pp[:, 4]
3921+ 2 * cos_eps * (pp [:, 0 ] + pp [:, 1 ] + pp [:, 2 ] + pp [:, 3 ])
3922+ - cos_eps * pp [:, 4 ]
39243923 )
39253924 return matrix
39263925
@@ -3995,24 +3994,23 @@ def __self_inductance_mode1(self, lmod: np.ndarray) -> np.ndarray:
39953994
39963995 for ii , obj in enumerate (self .inventory ["StrandComponent" ].collection ):
39973996 self_inductance [ii :: self .inventory ["StrandComponent" ].number ] = (
3998- 1.0e-7
3999- # 2
4000- # * lmod[ii :: self.inventory["StrandComponent"].number]
4001- # * (
4002- # np.arcsinh(
4003- # lmod[ii :: self.inventory["StrandComponent"].number]
4004- # / obj.radius
4005- # )
4006- # - np.sqrt(
4007- # 1.0
4008- # + (
4009- # obj.radius
4010- # / lmod[ii :: self.inventory["StrandComponent"].number]
4011- # )
4012- # ** 2
4013- # )
4014- # + obj.radius / lmod[ii :: self.inventory["StrandComponent"].number]
4015- # )
3997+ 2
3998+ * lmod [ii :: self .inventory ["StrandComponent" ].number ]
3999+ * (
4000+ np .arcsinh (
4001+ lmod [ii :: self .inventory ["StrandComponent" ].number ]
4002+ / obj .radius
4003+ )
4004+ - np .sqrt (
4005+ 1.0
4006+ + (
4007+ obj .radius
4008+ / lmod [ii :: self .inventory ["StrandComponent" ].number ]
4009+ )
4010+ ** 2
4011+ )
4012+ + obj .radius / lmod [ii :: self .inventory ["StrandComponent" ].number ]
4013+ )
40164014 )
40174015 return self_inductance
40184016
@@ -4029,26 +4027,25 @@ def __self_inductance_mode2(self, lmod: np.ndarray) -> np.ndarray:
40294027 for ii , obj in enumerate (self .inventory ["StrandComponent" ].collection ):
40304028
40314029 self_inductance [ii :: self .inventory ["StrandComponent" ].number ] = (
4032- 1.0e-7
4033- # 2 * (
4034- # lmod[ii :: self.inventory["StrandComponent"].number]
4035- # * np.log(
4036- # (
4037- # lmod[ii :: self.inventory["StrandComponent"].number]
4038- # + np.sqrt(
4039- # lmod[ii :: self.inventory["StrandComponent"].number] ** 2
4040- # + obj.radius ** 2
4041- # )
4042- # )
4043- # / obj.radius
4044- # )
4045- # - np.sqrt(
4046- # lmod[ii :: self.inventory["StrandComponent"].number] ** 2
4047- # + obj.radius ** 2
4048- # )
4049- # + lmod[ii :: self.inventory["StrandComponent"].number] / 4
4050- # + obj.radius
4051- )
4030+ 2 * (
4031+ lmod [ii :: self .inventory ["StrandComponent" ].number ]
4032+ * np .log (
4033+ (
4034+ lmod [ii :: self .inventory ["StrandComponent" ].number ]
4035+ + np .sqrt (
4036+ lmod [ii :: self .inventory ["StrandComponent" ].number ] ** 2
4037+ + obj .radius ** 2
4038+ )
4039+ )
4040+ / obj .radius
4041+ )
4042+ - np .sqrt (
4043+ lmod [ii :: self .inventory ["StrandComponent" ].number ] ** 2
4044+ + obj .radius ** 2
4045+ )
4046+ + lmod [ii :: self .inventory ["StrandComponent" ].number ] / 4
4047+ + obj .radius
4048+ ))
40524049
40534050 return self_inductance
40544051
0 commit comments