用Python对不规则导体几何均距、自感进行近似计算。

from dataclasses import dataclass
from typing import List
import math


@dataclass
class vec:
    """
    几何 点或向量
    """
    x: float
    y: float

    def __add__(self, other):
        return vec(self.x+other.x, self.y+other.y)

    def __sub__(self, other):
        return vec(self.x-other.x, self.y-other.y)

    def __mul__(self, n):
        x, y = self.x, self.y
        return vec(x*n, y*n)

    def __truediv__(self, n):
        x, y = self.x, self.y
        return vec(x/n, y/n)

    def unit(self):
        return self/self.length

    @property
    def length(self):
        x, y = self.x, self.y
        return (x*x+y*y)**0.5


def lower_bound(ls, key):
    '''
    二分查找
    '''
    le = 0
    ri = len(ls)-1
    ans = None
    while(le <= ri):
        mid = (le+ri) >> 1
        if(ls[mid] >= key):
            ri = mid-1
            ans = mid
        else:
            le = mid+1
    return ans


class polygon:
    '''
    多边形
    '''
    @property
    def n(self):
        return len(self.points)

    def __init__(self, points: List[vec]):
        self.points = points
        self.vecs = []
        self.start_length = []
        self.length = 0

        for i in range(self.n):
            vec = self.points[(i+1) % self.n]-self.points[i]
            self.vecs.append(vec)
            self.start_length.append(self.length)
            self.length += vec.length

    def get_point_at(self, length):
        '''
        沿周长走length长度后点的位置
        '''
        idx = lower_bound(self.start_length, length)    # 查找在哪一线段上
        if(idx is None):
            idx = len(self.start_length)-1
        else:
            idx -= 1
        le = length - self.start_length[idx]
        vec = self.vecs[idx]
        return self.points[idx]+vec.unit()*le   # 线段起点再移动剩余长度


def f(poly: polygon, step=None, stepN=1e6):
    '''
    ln G = 1/l² ∫∫ ln(dist(u,v))
    '''
    if(step is None):
        step = poly.length/(stepN**0.5)
        print(step)
    sum = 0
    tot = 0
    i = 0
    while(i < poly.length):
        u = poly.get_point_at(i)
        j = step/2
        while(j < i):

            v = poly.get_point_at(j)
            length = (u-v).length
            try:
                l = math.log(length)
            except Exception as e:

                j += step
                continue
            sum += l
            tot += 1
            j += step
        i += step
    return sum/tot


if(True and __name__ == "__main__"):
    # H型钢筋
    H = 600/1000
    B = 600/1000
    T1 = 20/1000
    T2 = 32/1000
    points = []
    points.append(vec(0, 0))
    points.append(vec(T2, 0))
    points.append(vec(T2, (B-T1)/2))
    points.append(vec(H-T2, (B-T1)/2))
    points.append(vec(H-T2, 0))
    points.append(vec(H, 0))
    points.append(vec(H, B))
    points.append(vec(H-T2, B))
    points.append(vec(H-T2, (B+T1)/2))
    points.append(vec(T2, (B+T1)/2))
    points.append(vec(T2, B))
    points.append(vec(0, B))
    poly = polygon(points)
    print(poly.start_length)
    import time
    t = time.time()
    print(f(poly, stepN=1e6))
    print(time.time()-t)