simulation.py

from datetime import datetime
import math
import os

time = f'{datetime.today().strftime('%Y-%m-%d__%H.%M.%S')}'

script_dir = os.path.dirname(os.path.abspath(__file__))
filePath = os.path.join(script_dir, 'results')
fileName = f'{time}.txt'

targetFile = os.path.join(filePath, fileName)

file = open(targetFile, 'a')

file.write(f'New experiment: {time}')

def logger(message):
    file.write(f'\n{message}')

dictionary = { }

def listDictionary(where, what):
    try:
        dictionary[where] # see whether where exists
    except:
        dictionary[where] = what # creates a new entry for where if it does not exist
        print(f'new dict value for {where}: {dictionary[where]}')
        return

    currentvalue = dictionary[where] # stores the currently attributed value for where

    if (type(what) == list):
        for value in what:
            listDictionary(where, value)
    else:
        if (type(currentvalue) == list):
            currentvalue.append(what)
            dictionary[where] = currentvalue
        else:    
            newList = [currentvalue, what]
            dictionary[where] = newList
    return

def writeDataDict():
    logger('dataDict = {')
    for x in dictionary:
        logger(f'    {x}: {dictionary[x]},')
    logger('}')


materialList = { # resistivity => soortelijke weerstand
	"copper": {"resistivity": 17E-9}
}

μ0 = 1.25663706212E-6

π = math.pi

class Magnet:
    def __init__(self, mass, area, magnetic_flux, length):
        self.mass = mass
        self.area = area
        self.magnetic_flux = magnetic_flux
        self.length = length

    def getVectorDipoleMoment(self):
        # Alternatively called magnetic moment
        return self.length * self.magnetic_flux
        # Definition used is weird, I can find this being said here and there, but no real source
        # Wikipedia, and most other sources, use moment = external magnetic field / torque acting on dipole, which is not usefull in this case


class Coil:
    def __init__(self, material, windings, length, area, height, Δt):
        if (material in materialList):
            self.material = material
        else:
            self.material = "copper"
        self.windings = windings
        self.length = length
        self.area = area
        self.height = height
        self.Δt = Δt

    def inducedCurrent(self, Δφ):
        current = -1 *  self.windings * Δφ / self.Δt
        listDictionary('U', current)
        return current


class FallingObject:
    def __init__(self, mass, area, height, velocity, Δt):
        self.mass = mass
        self.area = area
        self.velocity = velocity
        self.height = height
        self.Δt = Δt
        self.time = 0
        self.steps = 0

    def step(self, steps, skip=False):
        g = -9.81
        rho = 1.293 # BiNaS 12 Lucht - dichtheid
        cw = 0.82 # BiNaS 28A Cilinder Lang - Luchtsweerstandcoëfficient

        A = self.area
        m  = self.mass
        dt = self.Δt

        v = self.velocity
        t = self.time
        h = self.height

        i = 0

        while (i < steps):
            if (v < 0.0):
                Fw = 0.5 * rho * A * cw * v**2
            else:
                Fw = 0
            Fz = m * g
            Fres = Fz + Fw
            a = Fres / m
            dv = a * dt
            v = v + dv
            dh = v * dt
            h = h + dh
            t = t + dt

            self.steps = self.steps + 1
            i += 1
            listDictionary('t', t)
            listDictionary('h', h)
            listDictionary('v', v)
            if (skip): continue
            logger(f'========[Step {self.steps}]========\nt: { t }; h: { h }; dh: { dh }; v: { v }; dv: { dv }; a: { a }; Fres: { Fres }; Fz: { Fz }; Fw: { Fw };')

        self.time = t
        self.height = h
        self.velocity = v


class MagneticFlux:
    def __init__(self, moment):
        self.moment = moment
        self.oldFlux = 0

    def getΔφ(self, radius, angle=0):
        oldFlux = self.oldFlux
        newFlux = self.__calculateφ(radius, angle)
        Δφ = newFlux-oldFlux
        logger(f'radius: {radius}; oldFlux: {oldFlux}; newFlux: {newFlux};')
        listDictionary('r', radius)
        self.oldFlux = newFlux
        listDictionary('dP', Δφ)
        return Δφ

    def __calculateφ(self, radius, angle, vector = []):
        if (radius == 0): return 0
        Flux = (μ0 / 4*π) * (self.moment / (radius**3) ) * math.sqrt(1+3*math.sin(angle)**2)
        # Wikipedia > Dipole > Field of a static magnetic dipole > Magnitude
        return Flux
    
    # both deprecated due to using magnitude instead of vector form
    def __getUnitVector(self, vector, absVal):
        i = 0
        while i < len(vector):
            vector[i] = vector[i] / absVal
            i +=1
        return vector
    
    def __getAbsoluteValue(self, vector):
        absoluteValue = 0
        for value in vector:
            absolute += value**2

        return math.sqrt(absoluteValue)


def main(Δt, h):
    print('initiated main function')
    mass = 0.009; magnetArea = 0.0004; flux = 0.1; magnetLength = 0.05
    material = "copper"; windings = 200; coilLength = 0.1; coilArea = 0; coilHeight = 0.36
    velocity = 0

    magnet = Magnet(mass, magnetArea, flux, magnetLength)
    coil = Coil(material, windings, coilLength, coilArea, coilHeight, Δt)

    moment = magnet.getVectorDipoleMoment()

    logger(f'starting height: {h}; Delta t: {Δt}')
    logger(f'Magnet stats:\n- mass: {mass}\n- area: {magnetArea}\n- magnetic field: {flux}\n- length: {magnetLength}\n- moment: {moment}')
    logger(f'Coil stats:\n- material: {material}\n- windings: {windings}\n- length: {coilLength}\n- area: {coilArea}\n- height: {coilHeight}')
    logger('='*50)

    fallingMagnet = FallingObject(mass, magnetArea, h, velocity, Δt)
    fluxField = MagneticFlux(moment)

    # fallingMagnet.step(40)

    while fallingMagnet.height > 0:
        fallingMagnet.step(1)
        Δφ = fluxField.getΔφ((fallingMagnet.height - coil.height), 0)
        U = coil.inducedCurrent(Δφ)
        logger(f'Delta Phi (flux): {Δφ}; U: {U};')

    # fallingMagnet.step(10)
    # print(fluxField.getΔφ(0, 0.1))
    # print((μ0 / 4*π) * ( ( 3 * ( moment * 1 ) * 1 - moment ) / (0.1**3) ))

    with open(targetFile) as f:
        print(f.read())

    writeDataDict()

main(0.001, 1)