import matplotlib.pyplot as plt
import numpy as np
from src import coil_class as BC

import magpylib as magpy
import matplotlib.ticker as ticker


def circular_HH(HH, radius, distance, layers, windings,
                I = 1, #A
                wire_height = .5 *1e-3,
                wire_width = .5 *1e-3,
                insulation_thickness = (0.568-0.5)/2 *1e-3,
                print_info = True):
    """returns magpylib system for a circular HH-coil (parameters in m) and prints info on it if requested"""

    #### give out coil parameter s####
    HH_Coil = BC.BCoil(HH = HH, distance = distance*1e3, radius = radius*1e3, layers = layers, windings = windings,
                    wire_height = wire_height*1e3, wire_width = wire_width*1e3, insulation_thickness = insulation_thickness*1e3,
                    is_round = True, winding_scheme= 2)

    if print_info:
        print("Single coil parameters:")

        HH_Coil.print_info()
        print(f"N = {HH_Coil.get_N()}")
        print(f"wire length = {HH_Coil.get_wire_length()} m")

        HH_Coil.plot_raster()

        #### total inductance from Lenny's code ####
        def Self_Inductance_Abbot(alpha, beta, gamma, N):
            """
            Formula for self inductance of circular coil by J.J.Abbott
            """
            L = (7.9e-6 * alpha**2 * N**2) / (3*alpha + 9*beta + 10*gamma)
            return L

        alpha = 2*radius
        beta = HH_Coil.get_coil_width()
        gamma = HH_Coil.get_coil_height()
        N = HH_Coil.get_N()
        dist = distance

        L_pair = Self_Inductance_Abbot(alpha, dist + beta, gamma, N) +Self_Inductance_Abbot(alpha, dist - beta, gamma, N) - 2 * Self_Inductance_Abbot(alpha, dist, gamma, N) + 2 * Self_Inductance_Abbot(alpha, beta, gamma, N)
        R_pair = 2*HH_Coil.get_wire_length() * 0.08708

        print("Pair in series:")
        print(f"wire length = {2*HH_Coil.get_wire_length()} m")
        print(f"R = {2*HH_Coil.resistance(25)} Ohm")
        print(f"L = {L_pair*1e3} mH")
        print(f"tau = {L_pair/R_pair*1e3} ms")

        print(f"\nat {I} A:")
        print(f"V = {R_pair*I} V")
        print(f"P = {R_pair*I**2} W")

    #### calculate field properties ####

    #define coils
    winding_offsets = HH_Coil.full_raster(raster_value=1)[:,0]

    coil1 = magpy.Collection()
    for coord in winding_offsets:
        winding = magpy.current.Circle(
            current=1,
            diameter=2*coord[1],
            position=(0,0,coord[0]),
        )
        coil1.add(winding)

    coil2 = magpy.Collection()
    for coord in winding_offsets:
        winding = magpy.current.Circle(
            current=HH*1,
            diameter=2*coord[1],
            position=(0,0,-coord[0]),
        )
        coil2.add(winding)

    #move to test if slight changes in positioning have a large effect
    #coil1.move([0.5e-3,0,0])

    # Combine them into a collection (like multiple windings)
    system = magpy.Collection(coil1, coil2)

    return system, 2*HH_Coil.get_wire_length()


def rect_current(current, a,b,position, corner_radius):
    """returns rectangular mapy current with given dimensions (in m)"""
    #first quadrant
    angles = np.linspace(0,np.pi/2,50)
    radii = np.zeros((len(angles),3))
    radii[:,0] = np.cos(angles)
    radii[:,1] = np.sin(angles)
    radii *= corner_radius
    vertices = np.array([a/2-corner_radius, b/2-corner_radius, 0]) + radii

    #second quadrant
    angles = np.linspace(np.pi/2,np.pi,50)
    radii = np.zeros((len(angles),3))
    radii[:,0] = np.cos(angles)
    radii[:,1] = np.sin(angles)
    radii *= corner_radius
    vertices = np.append(vertices,np.array([-(a/2-corner_radius), b/2-corner_radius, 0]) + radii, axis=0)

    #third quadrant
    angles = np.linspace(np.pi,1.5*np.pi,50)
    radii = np.zeros((len(angles),3))
    radii[:,0] = np.cos(angles)
    radii[:,1] = np.sin(angles)
    radii *= corner_radius
    vertices = np.append(vertices,np.array([-(a/2-corner_radius), -(b/2-corner_radius), 0]) + radii, axis=0)

    #fourth quadrant
    angles = np.linspace(1.5*np.pi, 2*np.pi,50)
    radii = np.zeros((len(angles),3))
    radii[:,0] = np.cos(angles)
    radii[:,1] = np.sin(angles)
    radii *= corner_radius
    vertices = np.append(vertices,np.array([a/2-corner_radius, -(b/2-corner_radius), 0]) + radii, axis=0)

    #close loop
    vertices = np.append(vertices, [vertices[0]], axis=0)

    # Create a rectangular current coil
    coil1 = magpy.current.Polyline(position=position, vertices=vertices, current=current)
    return coil1


def rectangular_HH(HH, a, b, distance, corner_radius, layers, windings,
                    I = 1, #A
                    wire_height = .5 *1e-3,
                    wire_width = .5 *1e-3,
                    insulation_thickness = (0.568-0.5)/2 *1e-3,
                    print_info = True):
    """returns magpylib system for a rectangular HH-coil (parameters in m) and prints info on it if requested"""

    #### give out coil parameter s####
    HH_Coil = BC.BCoil(HH = HH, distance = 0, radius = 0, layers = layers, windings = windings,
                    wire_height = wire_height*1e3, wire_width = wire_width*1e3, insulation_thickness = insulation_thickness*1e3,
                    is_round = True, winding_scheme= 2)

    if print_info:
        N = HH_Coil.get_N()
        print("Single coil parameters: \n")
        #HH_Coil.print_info()
        print(f"coil width = {HH_Coil.get_coil_width()*1e3} mm")
        print(f"coil height = {HH_Coil.get_coil_height()*1e3} mm")
        print(f"N = {HH_Coil.get_N()}")

        wire_length = N*(2*a - 4*corner_radius + 2*b - 4*corner_radius + 2*np.pi*corner_radius)
        print(f"wire length = {wire_length} m")

        HH_Coil.plot_raster()


        #### total inductance from Lenny's code ####
        #convert to centimeters
        a *= 100
        b *= 100
        distance *= 100
        p1 = (a + np.sqrt(a ** 2 + distance ** 2)) * np.sqrt(b ** 2 + distance ** 2) / ((a + np.sqrt(a ** 2 + b ** 2 + distance ** 2)) * distance)
        p2 = (b + np.sqrt(b ** 2 + distance ** 2)) * np.sqrt(a ** 2 + distance ** 2) / ((b + np.sqrt(a ** 2 + b ** 2 + distance ** 2)) * distance)
        p3 = np.sqrt(a ** 2 + b ** 2 + distance ** 2) - np.sqrt(a ** 2 + distance ** 2) - np.sqrt(b ** 2 + distance ** 2) + distance
        M = 4 * (a * np.log(p1) + b * np.log(p2)) + 8 * p3
        M *= N**2 * 1e-9 #Conversion to SI and accounting for N wires

        r = (wire_width/2 + insulation_thickness) * 100
        diag = np.sqrt(a**2 + b**2)
        L_self = 4 * ( (a+b) * np.log(2*a*b / r) - a * np.log(a+diag) - b * np.log(b+diag) - 7/4 * (a+b) + 2 * (diag + r) )
        L_self *= N**2 * 1e-9 #Conversion to SI and accounting for N wires

        L_pair = 2 * (L_self+M)

        R_pair = 2*wire_length * 0.08708

        print("Pair in series:")
        print(f"wire length = {2*wire_length} m")
        print(f"R = {R_pair} Ohm")
        print(f"L = {L_pair*1e3} mH")
        print(f"tau = {L_pair/R_pair *1e3} ms")

        print(f"\nat {I} A:")
        print(f"V = {R_pair*I} V")
        print(f"P = {R_pair*I**2} W")

    #### calculate field properties ####
        #convert back to meters
        a /= 100
        b /= 100
        distance /= 100

    #define coils
    winding_offsets = HH_Coil.full_raster(raster_value=1)[:,0]

    coil1 = magpy.Collection()
    for coord in winding_offsets:
        winding = rect_current(
            current=1,
            a=a+coord[1], b=b+coord[1],
            position=[0,0,distance/2+coord[0]],
            corner_radius=corner_radius+coord[1]/2
        )
        coil1.add(winding)

    coil2 = magpy.Collection()
    for coord in winding_offsets:
        winding = rect_current(
            current=HH*1,
            a=a+coord[1], b=b+coord[1],
            position=[0,0,-(distance/2+coord[0])],
            corner_radius=corner_radius+coord[1]/2
        )
        coil2.add(winding)

    #move to test if slight changes in positioning have a large effect
    #coil1.move([0.5e-3,0,0])

    # Combine them into a collection (like multiple windings)
    system = magpy.Collection(coil1, coil2)

    return system, 2*wire_length