import numpy,scipy,scipy.integrate
from numpy import pi,sqrt
from scipy.constants import mu_0 as u0
from scipy.special import jv as besselj
from scipy.special import ellipe,ellipk
import warnings
import math  
from scipy.interpolate import interp1d  
factor = 1 #*0.96,0.96为层间交错导致的层间距修正系数

def coil_calculate_O(outer_Radius_mm, wire_size, Paint, inner_Radius_mm, coil_length, Coil_material):
              #输入 外径,铜线直径,涂漆厚度,支架直径,线圈长度,材料选择（铜/铝）
    Coil_material_dict={"Cu":[8920,5.8*1e+7],"Al":[2700,3.7*1e+7]}      #线圈材料(密度，导电率)
    
    Enameled_wire_size = wire_size+2*Paint*1e-3 #漆包线外径
    
    q=int(coil_length/Enameled_wire_size)  #q每层缠绕数
    #qend=N%q      #未满的最后一层的缠绕数(p+1层)
    #p=int(N/q)    #p满层层数
    Length=0  
    D_layer=inner_Radius_mm - Enameled_wire_size
    p=0
    for layer in range(100):   #循环计算每层直径，每层线长
        D_layer=D_layer+2*Enameled_wire_size*factor #*0.96,0.96为层间交错导致的层间距修正系数
        #Length=Length+( sqrt( (D_layer* pi)**2 + Enameled_wire_size**2 ))*q    #累加计算线长
        Length=Length+( pi*D_layer)*q
        
        if(D_layer+2*Enameled_wire_size*factor+Enameled_wire_size > outer_Radius_mm):
            p=layer+1
            break
    #print("p:",p)        
           
           
    #外径计算(单位 mm)
    #outer_Radius_mm = D_layer + (qend/q)*2*Enameled_wire_size*factor +Enameled_wire_size  #计算外径 并对最后一层修正
    qend=(outer_Radius_mm-D_layer-Enameled_wire_size)/(2*Enameled_wire_size*factor/q)
    #print("qend:",qend)   
    N=p*q+qend  
    #print("N:",N)  
    Length=Length+( pi*D_layer)*qend  #计算导线长度 并对最后一层修正
    
    #线圈电阻计算(单位 mohm)
    Resistivity = (1 / (Coil_material_dict[Coil_material][1]))  #Ω.m
    Res_Metre = (Resistivity)/((1e-3*wire_size/2)**2* pi) #每米电阻 (欧)
    Res_mohm=(Length/1000)*Res_Metre*1e+3       #总电阻=每米电阻×线长 (单位毫欧)
    
    #线圈重量计算(单位 g)
    Coil_Weigh_g=(Length/10)*((wire_size/20)**2* pi)*Coil_material_dict[Coil_material][0]*1e-3    
    
    #根据Wheeler公式 计算电感量(单位 uH)
    M=(inner_Radius_mm+outer_Radius_mm)/2  #线圈平均直径
    B=coil_length                          #线圈长度
    C=(outer_Radius_mm-inner_Radius_mm)/2  #线圈径向长度
    L_uh=1e-3*(7.87*N**2*M**2)/(3*M+9*B+10*C)
    
    #线圈时间常数计算(单位 us)
    coil_time_us=1000*L_uh/Res_mohm
    #print(N,Res_mohm,wire_size,outer_Radius_mm/inner_Radius_mm,coil_time_us)
    #输出 每层缠绕数，满层的层数，未满层的缠绕数，外径，线圈阻抗，电感，线圈重量，时间常数
    return q,p,qend,N,outer_Radius_mm*1e-3,Res_mohm*1e-3,L_uh*1e-6,Coil_Weigh_g,coil_time_us

def coil_calculate_N(N,wire_size,Paint,inner_Radius_mm,coil_length,Coil_material):  
              #输入 缠绕圈数,铜线直径,涂漆厚度,支架直径,线圈长度,材料选择（铜/铝）
    Coil_material_dict={"Cu":[8920,5.8*1e+7],"Al":[2700,3.7*1e+7]}      #线圈材料(密度，导电率)
    
    Enameled_wire_size = wire_size+2*Paint*1e-3 #漆包线外径
    
    q=int(coil_length/Enameled_wire_size)  #q每层缠绕数
    qend=N%q      #未满的最后一层的缠绕数(p+1层)
    p=int(N/q)    #p满层层数
    Length=0  
    D_layer=inner_Radius_mm - Enameled_wire_size
    for layer in range(p):   #循环计算每层直径，每层线长
        D_layer=D_layer+2*Enameled_wire_size*1 #*0.96,0.96为层间交错导致的层间距修正系数
        #Length=Length+( sqrt( (D_layer* pi)**2 + Enameled_wire_size**2 ))*q    #累加计算线长
        Length=Length+( pi*D_layer)*q
    Length=Length+( pi*D_layer)*qend  #计算导线长度 并对最后一层修正

    #外径计算(单位 mm)
    outer_Radius_mm = D_layer + (qend/q)*2*Enameled_wire_size*1 +Enameled_wire_size  #计算外径 并对最后一层修正
    
    #线圈电阻计算(单位 mohm)
    Resistivity = (1 / (Coil_material_dict[Coil_material][1]))  #Ω.m
    Res_Metre = (Resistivity)/((1e-3*wire_size/2)**2* pi) #每米电阻 (欧)
    Res_mohm=(Length/1000)*Res_Metre*1e+3       #总电阻=每米电阻×线长 (单位毫欧)
    
    '''
    print(f'Length:{Length}')
    wire_size_yinxian = 1
    Length_yinxian = 2*390
    Res_Metre_yinxian = (Resistivity)/((1e-3*wire_size_yinxian/2)**2* pi)
    Res_mohm_yinxian= (Length_yinxian/1000)*Res_Metre_yinxian*1e+3
    Coil_Weigh_g_yinxian=(Length_yinxian/10)*((wire_size_yinxian/20)**2* pi)*Coil_material_dict[Coil_material][0]*1e-3    
    print(f'Coil_Weigh_g_yinxian:{Coil_Weigh_g_yinxian}')

    Res_mohm +=Res_mohm_yinxian'''
    
    #线圈重量计算(单位 g)
    Coil_Weigh_g=(Length/10)*((wire_size/20)**2* pi)*Coil_material_dict[Coil_material][0]*1e-3    
    
    #根据Wheeler公式 计算电感量(单位 uH)
    M=(inner_Radius_mm+outer_Radius_mm)/2  #线圈平均直径
    B=coil_length                          #线圈长度
    C=(outer_Radius_mm-inner_Radius_mm)/2  #线圈径向长度
    L_uh=1e-3*(7.87*N**2*M**2)/(3*M+9*B+10*C)
    
    #线圈时间常数计算(单位 us)
    coil_time_us=1000*L_uh/Res_mohm
    
    ##输出 序号，每层缠绕数：，满层的层数，未满层的缠绕数，匝数，自感，电阻
    #print(q,p,qend,N,outer_Radius_mm,Res_mohm,L_uh,Coil_Weigh_g,coil_time_us)
    return q,p,qend,N,outer_Radius_mm*1e-3,Res_mohm*1e-3,L_uh*1e-6,Coil_Weigh_g,coil_time_us


def coil_calculate_R(Res_mohm,wire_size,Paint,inner_Radius_mm,coil_length,Coil_material):  
              #输入 电阻,铜线直径,涂漆厚度,支架直径,线圈长度,材料选择（铜/铝）
    Coil_material_dict={"Cu":[8920,5.8*1e+7],"Al":[2700,3.7*1e+7]}      #线圈材料(密度，导电率)
    
    Enameled_wire_size = wire_size+2*Paint*1e-3 #漆包线外径
    
    #线圈电阻计算(单位 mohm)
    Resistivity = (1 / (Coil_material_dict[Coil_material][1]))  #Ω.m
    Res_Metre = (Resistivity)/((1e-3*wire_size/2)**2* pi) #每米电阻 (欧)
    
    Length = 1000*Res_mohm/(Res_Metre*1e+3) #线长
    
    q=int(coil_length/Enameled_wire_size)  #q每层缠绕数
    p=0
    
    D_layer=inner_Radius_mm - Enameled_wire_size
    lengthend=0
    Length_lins=Length
    for layer in range(100):   #循环计算每层直径，每层线长
        D_layer_lins=D_layer+2*Enameled_wire_size*1 #*0.96,0.96为层间交错导致的层间距修正系数
        Length2=Length_lins-( pi*D_layer_lins)*q
        if Length2 >0:
            D_layer =  D_layer_lins  
            Length_lins = Length2
        else:
            lengthend =Length_lins
            p= layer+1 -1
            break
            
    qend = lengthend/( pi*D_layer)
    
    N=p*q+qend

    #外径计算(单位 mm)
    outer_Radius_mm = D_layer + (qend/q)*2*Enameled_wire_size*1 +Enameled_wire_size  #计算外径 并对最后一层修正
    
    #线圈重量计算(单位 g)
    Coil_Weigh_g=(Length/10)*((wire_size/20)**2* pi)*Coil_material_dict[Coil_material][0]*1e-3    
    
    #根据Wheeler公式 计算电感量(单位 uH)
    M=(inner_Radius_mm+outer_Radius_mm)/2  #线圈平均直径
    B=coil_length                          #线圈长度
    C=(outer_Radius_mm-inner_Radius_mm)/2  #线圈径向长度
    L_uh=1e-3*(7.87*N**2*M**2)/(3*M+9*B+10*C)
    
    #线圈时间常数计算(单位 us)
    coil_time_us=1000*L_uh/Res_mohm
    
    ##输出 序号，每层缠绕数：，满层的层数，未满层的缠绕数，匝数，自感，电阻
    #print(q,p,qend,N,outer_Radius_mm,Res_mohm,L_uh,Coil_Weigh_g,coil_time_us)
    return q,p,qend,N,outer_Radius_mm*1e-3,Res_mohm*1e-3,L_uh*1e-6,Coil_Weigh_g,coil_time_us




 
def calculate_M(inner_Radius,A_outer_Radius,A_Length,A_N,B_outer_Radius,B_Length,B_N ,spacer): # 计算互感
    warnings.filterwarnings("ignore")#忽略警告
    
    A_Ir=(inner_Radius/2)    #线圈A内半径
    A_Or=(A_outer_Radius/2)  #线圈A外半径
    #A_Length=               #线圈A长度
    
    B_Ir= A_Ir               #线圈B内半径
    B_Or=(B_outer_Radius/2)  #线圈B外半径
    #B_Length=               #线圈B长度
    
    #spacer=        #隔板距离(俩线圈距离)

    r11=A_Ir
    r12=A_Or
    r21=B_Ir
    r22=B_Or
    z11=0
    z12=A_Length
    z21=z12+spacer
    z22=z21+B_Length
    N1=A_N
    N2=B_N

    INTEG_STEPS = 5
    def Ml(a,z1,b,z2):
        h=z1-z2
        k=sqrt(4*a*b/((a+b)**2+h**2))
        return u0*numpy.sqrt(a*b)*((2/k-k)*ellipk(k**2)-2*ellipe(k**2)/k)
    S1 = (r12-r11) * (z12-z11)
    S2 = (r22-r21) * (z22-z21)
    # M = (N1*N2/S1/S2) * scipy.integrate.nquad(Ml,[[r11,r12],[z11,z12],[r21,r22],[z21,z22]],opts={"limit":3})[0]
    M = 0
    temp4 = []
    for r1 in numpy.linspace(r11,r12,INTEG_STEPS):
        temp3 = []
        for z1 in numpy.linspace(z11,z12,INTEG_STEPS):
            temp2 = []
            for r2 in numpy.linspace(r21,r22,INTEG_STEPS):
                temp1 = []
                for z2 in numpy.linspace(z21,z22,INTEG_STEPS):
                    temp1.append(Ml(r1,z1,r2,z2))
                M = scipy.integrate.simpson(temp1, numpy.linspace(z21,z22,INTEG_STEPS))
                temp2.append(M)
            M = scipy.integrate.simpson(temp2, numpy.linspace(r21,r22,INTEG_STEPS))
            temp3.append(M)
        M = scipy.integrate.simpson(temp3, numpy.linspace(z11,z12,INTEG_STEPS))
        temp4.append(M)
    M = scipy.integrate.simpson(temp4, numpy.linspace(r11,r12,INTEG_STEPS))

    M = (N1*N2/S1/S2) * M 
    return M
    
    
def y_function(x, attenuation_coefficient=0.8):   #
    # 我们可以使用x的对数来模拟这种衰减，但结果会是一个平滑的衰减，而不是离散的。  
    return 1 * (attenuation_coefficient ** (math.log2(x))  )  
  

#[["Cu",0.6,3,2.8,10],["Al",0.75,3,2.8,30]]
def calculate_all(scanning_select,Mcaulate_select,Inner_diameteri,L_coil1,coil_layout):
    wires=[]
    Ns=[]
    Ls=[]
    Ms=[]
    Rs=[]

    Outer_diameter=[]
    Coil_length=[]
    pw=[]
    pwk=[]
    
    Coil_weight_list=[]
    
    Coil_numi=0 #线圈编号
    
    Mk2=0 #互感系数
    
    for coil in coil_layout:
        #print("+++70:::",coil)

        
        
        mytype = coil['匝数或外径模式选择(0or1)']
        Coil_material_type=["Cu","Al"]
        Coil_material=Coil_material_type[ int(coil['漆包线材料(CuorAl)']) ]
        

        
        wire_size =coil['漆包线线径(mm)']
        thinkness_um = coil['漆层厚度(um)']
        L_coil0   =coil['线圈长度(mm)']
        Ostart    =coil['匝径start']
        attenuation_coefficient      =coil['匝数幂衰减系数']
        Coil_num      =coil['线圈数量']
        pwii      =coil['脉宽(pw)']
        pw_coefficient = coil['脉宽幂递增系数']
        pwki =coil['关断策略偏移(mm)']

        Ostartnum=  Coil_numi
        

        
        if(Ostart==0): #使用匝数延续功能
            Ostart= (Ns[-1]/y_function(Coil_numi,attenuation_coefficient) )*y_function(Coil_numi+1,attenuation_coefficient) 
        
        else:
            if mytype == 1: #选择输入O,外径 均转化为 输入N,匝数
                q,p,qend,Nstart,Outer_diameteri,Rsi,Lsi,Coil_weighti,coil_time_usi = coil_calculate_O(Ostart,wire_size,thinkness_um,Inner_diameteri*1000,L_coil0,Coil_material)
                Ostart = Nstart
            elif mytype == 2: #选择输入R,电阻 均转化为 输入N,匝数  
                Rmes = (480/Ostart )*1000        
                q,p,qend,Nstart,Outer_diameteri,Rsi,Lsi,Coil_weighti,coil_time_usi = coil_calculate_R(Rmes,wire_size,thinkness_um,Inner_diameteri*1000,L_coil0,Coil_material)
                Ostart = Nstart
                  
        pwendnum=Coil_numi+Coil_num-1
        for i in range(int(Coil_num)):

            pwi = float( (float(pwii)/y_function(pwendnum+1,pw_coefficient) )*y_function(Coil_numi+1,pw_coefficient) )  
  
            N_O = int( (Ostart/y_function(Ostartnum +1,attenuation_coefficient) )*y_function(Coil_numi+1,attenuation_coefficient) )
            
            q,p,qend,N,Outer_diameteri,Rsi,Lsi,Coil_weighti,coil_time_usi = coil_calculate_N(N_O,wire_size,thinkness_um,Inner_diameteri*1000,L_coil0,Coil_material)

            if scanning_select== '非参数化扫描':
                print(f"线圈编号：{Coil_numi} |匝数：{int(N)} ={q}*{p}+{int(qend)} |自感：{int(Lsi*1e+6)}uh |阻值：{int(Rsi*1000)}mΩ |外径：{int(Outer_diameteri*10000)/10}mm |重量：{Coil_weighti:.2f}g ")
            
            
            if i==0:
                wires.append(round(wire_size,2))
            else:
                wires.append("_")
            Ns.append(N)
            Ls.append(Lsi)
            
            Rs.append(Rsi)
            Coil_length.append(L_coil0*1e-3) 
            Outer_diameter.append(Outer_diameteri)
            
            Coil_weight_list.append(Coil_weighti)
            
            pw.append(pwi*1e-3)
            pwk.append(pwki*1e-3)
            
            Coil_numi+=1 
            
           
    #计算线圈位置
    coil_number=sum([coil_l['线圈数量'] for coil_l in coil_layout]    ) #线圈总数
    coil_location=[0]
    for i in range(1,coil_number) :
        coil_location.append(0.5*Coil_length[0]+ sum(Coil_length[1:i]) + i*L_coil1   +0.5*Coil_length[i]  )
    

    
    for i in range(coil_number):
        #计算线圈互感
        M_lin={'-2':0,'-1':0,'1':0,'2':0}
        Mk_lin={'-2':0,'-1':0,'1':0,'2':0}
        if(Mcaulate_select =='精细计算互感'):
        
            if i > 2:
                M_lin['-2']= calculate_M(Inner_diameteri,Outer_diameter[i],Coil_length[i] ,Ns[i],Outer_diameter[i-2],Coil_length[i-2] ,Ns[i-2]  ,2*L_coil1+Coil_length[i-1]     )
            if i > 1:
                M_lin['-1']= calculate_M(Inner_diameteri,Outer_diameter[i],Coil_length[i] ,Ns[i],Outer_diameter[i-1],Coil_length[i-1] ,Ns[i-1]  ,L_coil1    )
     
            if i < coil_number-1:
                M_lin['1']= calculate_M(Inner_diameteri,Outer_diameter[i],Coil_length[i] ,Ns[i],Outer_diameter[i+1],Coil_length[i+1] ,Ns[i+1]  ,L_coil1    )
     
            if i < coil_number-2:
                M_lin['2']= calculate_M(Inner_diameteri,Outer_diameter[i],Coil_length[i] ,Ns[i],Outer_diameter[i+2],Coil_length[i+2] ,Ns[i+2]  ,2*L_coil1+Coil_length[i+1]     )
        
        elif(Mcaulate_select =='粗略计算互感') :   #默认粗略计算互感
            M_lin['-1'] = calculate_M(Inner_diameteri,Outer_diameter[i],Coil_length[i] ,Ns[i],Outer_diameter[i],Coil_length[i] ,Ns[i]  ,L_coil1    )
            M_lin['1']  = M_lin['-1']
            
            M_lin['-2'] = M_lin['-1']*0.2  #默认系数0.2 ,默认互感分布类似于：[0.08,0.4,1,0.4,0.08]
            M_lin['2']  = M_lin['-2']
            
        elif(Mcaulate_select =='不计算互感') :   #默认粗略计算互感
            M_lin['-1'] = Ls[i]*0.2
            M_lin['1']  = M_lin['-1']
            
            M_lin['-2'] = M_lin['-1']*0.2  #默认系数0.2 ,默认互感分布类似于：[0.08,0.4,1,0.4,0.08]
            M_lin['2']  = M_lin['-2']
        
        
        Ms.append(M_lin)
    #print(Ms)    

    #Coil_weight_all = sum(Coil_weight_list)
    return wires,Ns,Ls,Ms,Rs,Coil_length,Outer_diameter,coil_location,pw,pwk ,Coil_weight_list
    
    
    
'''
wire_size,Paint,inner_Radius_mm,coil_length,Coil_material =(0.38,25,8.5,10,'Cu')

q2,p2,qend2,N2,Outer_diameteri2,Rsi2,Lsi2,Coil_weighti2,coil_time_usi2 = coil_calculate_N(200,wire_size,Paint,inner_Radius_mm,coil_length,Coil_material)

print(f"线圈编号：{0} |匝数：{int(N2)} ={q2}*{p2}+{int(qend2)} |自感：{int(Lsi2*1e+6)}uh |阻值：{int(Rsi2*1000)}mΩ |外径：{int(Outer_diameteri2*10000)/10}mm |重量：{int(Coil_weighti2)}g ")

q2,p2,qend2,N2,Outer_diameteri2,Rsi2,Lsi2,Coil_weighti2,coil_time_usi2 = coil_calculate_O(15.9,wire_size,Paint,inner_Radius_mm,coil_length,Coil_material)

print(f"线圈编号：{1} |匝数：{int(N2)} ={q2}*{p2}+{int(qend2)} |自感：{int(Lsi2*1e+6)}uh |阻值：{int(Rsi2*1000)}mΩ |外径：{int(Outer_diameteri2*10000)/10}mm |重量：{int(Coil_weighti2)}g ")


q2,p2,qend2,N2,Outer_diameteri2,Rsi2,Lsi2,Coil_weighti2,coil_time_usi2 = coil_calculate_R(1163,wire_size,Paint,inner_Radius_mm,coil_length,Coil_material)

print(f"线圈编号：{2} |匝数：{int(N2)} ={q2}*{p2}+{int(qend2)} |自感：{int(Lsi2*1e+6)}uh |阻值：{int(Rsi2*1000)}mΩ |外径：{int(Outer_diameteri2*10000)/10}mm |重量：{int(Coil_weighti2)}g ")
'''