DriCy_Selo.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Created on Thu Mar  4 08:53:28 2021

@author:    Michael Schlüter \n
            Technische Universität Berlin \n
            https://www.pe.tu-berlin.de/ \n
            m.schlueter@tu-berlin.de
"""


from unittest.loader import VALID_MODULE_NAME
import toml #used for conifg data (vehicles, semiconductor, ...)
import pickle #used to save output data to plot later...
import click
import numpy as np
import numpy.polynomial.polynomial as poly
import matplotlib.pyplot as plt
import matplotlib.font_manager
import os
from csv import reader
import sympy as sy
from scipy import signal
from pathlib import Path
import pyfiglet
import datetime

def f(question):
    """Answer to the Ultimate Question of Life, the Universe, and Everything.
    tbd: Find the inverse function..."""  
    return 42


def read_toml(file):
    """This function loads the data of the vehicle."""
    vehicle = toml.load(file)
    Felge_d = vehicle['basic']['rim_dia']*0.0254
    H = vehicle['basic']['H_B'] * vehicle['basic']['B'] * 0.01
    vehicle['basic']['wheel_radius'] = (Felge_d + (2*H))/2
    vehicle['basic']['i_g'] = vehicle['basic']['n_max'] / 60 *(vehicle['basic']['wheel_radius'] * 2 * np.pi) / (vehicle['basic']['v_max']/3.6)

    return(vehicle)

def read_csv(file):
    """Read csv export from automeris.io
    Column Separator: comma [,]
    decimal separator: dot[.]"""
    with open(file, 'r') as read_obj: # read csv file as a list of lists
        csv_reader = reader(read_obj) # pass the file object to reader() to get the reader object
        list_of_rows = list(csv_reader) # Pass reader object to list() to get a list of lists
        Data = np.asarray(list_of_rows).astype(float)
        return Data

def rolling_fiction(vehicle, v, alpha):
    """This function calculates the rolling friction acording to Mitschke, Wallentowitz 2014 - Dynamik der Kraftfahzeuge p. 14.:
    
    .. math::
        
        F_{R} = (M_{w} + M_{z}) \\cdot g \\cdot cos(\\alpha) \\cdot \\left( f_{R0} + f_{R1} \\cdot \\frac{v}{100 km/h} + f_{R4} \\cdot  {\\left(\\frac{v}{100 km/h}\\right)}^4\\right)
    
    
    """
    if v == 0:
        return(0)
    else:
        return((vehicle['basic']['m_w'] + vehicle['basic']['m_z']) * vehicle['basic']['g_earth'] * np.cos(alpha /180 * np.pi) * (vehicle['basic']['f_B'] + vehicle['basic']['f_R0'] + (vehicle['basic']['f_R1'] * (v*3.6)/100) + (vehicle['basic']['f_R4'] * ((v*3.6)/100)**4)))

def aerodynamic_resistance(vehicle, v):
    """This function calculates the aerodynamic resistance acording to:
    
    .. math::
        
        F_{L} = c_{w} \\cdot A_{front} \\cdot \\frac{p_{air}}{2} \\cdot v^2
    
    
    """
    return(0.5 * vehicle['basic']['cw'] * vehicle['basic']['A_front'] * vehicle['basic']['p_air'] * v**2)

def acceleration_resistance(a, vehicle):
    """This function calculates the acceleration resistance acording to:
    
   .. math::
        
        F_{A} = a \\cdot (M_{w} \\cdot \\lambda + M_{z})
     
    
    """
    return(a * (vehicle['basic']['m_w'] * vehicle['basic']['lambda'] + vehicle['basic']['m_z']))

def slope_resistance(vehicle, alpha):
    """This function calculates the slope resistance acording to:
    
    .. math::
    
        F_{S} = sin (\\alpha) \\cdot (M_{w} + M_{z}) \\cdot g
    
    """
    return(np.sin (alpha * np.pi/180) * (vehicle['basic']['m_w'] + vehicle['basic']['m_z']) * vehicle['basic']['g_earth'])

def total_force(vehicle, a, v, alpha):
    """This function calculates the total force for a vehicle in one operation point:
    
    .. math::
    
        F_{Total} = F_{S} + F_{A} + F_{L} + F_{R}
        
    """
    return(slope_resistance(vehicle, alpha) + acceleration_resistance(a, vehicle) + aerodynamic_resistance(vehicle, v) + rolling_fiction(vehicle, v, alpha))

def dq0_to_abc(d, q, wt, z=0, delta=0):
    """Inverse Park transform"""
    a = d*np.cos(wt+delta) - q*np.sin(wt+delta) + z
    b = d*np.cos(wt-(2*np.pi/3)+delta) - q*np.sin(wt-(2*np.pi/3)+delta) + z
    c = d*np.cos(wt+(2*np.pi/3)+delta) - q*np.sin(wt+(2*np.pi/3)+delta) + z
    return a, b, c

def calc_PMSM(vehicle, T, n, verbose = False):
    """Calculation of the inverter set point for a given vehicle, tourque and speed of a permanent magnet synchronos motor (PMSM).
    ToDo: speed up -> add complete numerical calculation
    """
    erg = {}
    if T == 0:
        erg["status"] = 'Status: Basic Speed Range'
        erg["i_d"] = 0
        erg["i_q"] = 0
        erg["v_d"] = 0
        erg["v_q"] = 0
        return (erg)
     
 
    if n > vehicle['basic']['n_max']/60:
        erg["status"] = 'Error: Rotational limit reached!'
        return(erg)

    i_d = sy.symbols('i_d')
    i_dlist = (sy.solveset(sy.Eq(sy.diff(( i_d**2 + ( (2*T) / ( 3* vehicle['electric']['zp'] * (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d) ) )**2) ,i_d),0), i_d, sy.Interval(-vehicle['electric']['I_max'],0)))
    if len(list(i_dlist)) == 0:
        erg["status"] = 'Error: Current limit reached!'
        return(erg)
    i_ds = list(i_dlist)[0]
    erg["i_d"] = i_ds
    i_q = (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])* i_ds ))
    erg["i_q"] = i_q
    if verbose:
        erg["fuc_T"] = sy.lambdify(i_d, (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d )) )
        erg["MTPC"] = (i_ds, i_q)
    if i_ds**2 + i_q**2  > vehicle['electric']['I_max']**2:
        erg["status"] = 'Error: Current limit reached!'
        return(erg)
    omega = 2*sy.pi* n*vehicle['electric']['zp']
    v_max = vehicle['electric']['V_DC'] / np.sqrt(3) * vehicle['electric']['V_DC_res']
    v_d = (i_ds * vehicle['electric']['Rs'] - omega*vehicle['electric']['Lq'] * i_q).evalf()
    erg["v_d"] = v_d
    v_q = (i_q * vehicle['electric']['Rs'] + (vehicle['electric']['Ld'] * i_ds + vehicle['electric']['PSI'])* omega).evalf()
    erg["v_q"] = v_q
    
    if verbose:
        if T > 0: 
            erg['fuc_Imax'] = sy.lambdify(i_d, sy.sqrt(vehicle['electric']['I_max']**2 - i_d**2))
        elif T < 0:
            erg['fuc_Imax'] = sy.lambdify(i_d, -sy.sqrt(vehicle['electric']['I_max']**2 - i_d**2))

    if sy.sqrt( v_d**2 + v_q**2) > v_max:
        if T > 0:                
            func = sy.lambdify(i_d, (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d ) ) - sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ))
            with np.errstate(invalid='ignore'):
                xx = np.linspace (-vehicle['electric']['I_max'],0,10000)
                if np.nanmin(func(xx)) > 0:
                    # erg["fuc_V"] = func
                    erg["fuc_V"] = sy.lambdify(i_d, sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ) )
                    erg["status"] = 'Error: Voltage limit reached!'
                    return(erg)                    

            with np.errstate(invalid='ignore'):
                A = np.diff(np.sign(func(xx)))
            A[np.isnan(A)] = 0
            test = xx[np.where(A)[0]]
            
            i_ds = test[0]
            erg["i_d"] = i_ds
            i_q = ( (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*list(test)[0] )) ) #.evalf() )
            erg["i_q"] = i_q
            erg["status"] = 'Status: Field-Weakening!'
            v_d = (i_ds * vehicle['electric']['Rs'] - omega*vehicle['electric']['Lq'] * i_q).evalf()
            erg["v_d"] = v_d
            v_q = (i_q * vehicle['electric']['Rs'] + (vehicle['electric']['Ld'] * i_ds + vehicle['electric']['PSI'])* omega).evalf()
            erg["v_q"] = v_q
            if verbose:
                erg["fuc_V"] = sy.lambdify(i_d, sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ) )
            if i_ds**2 + i_q**2  > vehicle['electric']['I_max']**2:
                erg["status"] = 'Error: Current limit under voltage limit reached!'
                return(erg)
            return(erg)
        elif T < 0:                
            func = sy.lambdify(i_d, (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d ) ) + sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ))
            with np.errstate(invalid='ignore'):
                xx = np.linspace (-vehicle['electric']['I_max'],0,10000)
                if np.nanmax(func(xx)) < 0:
                    # erg["fuc_V"] = func
                    erg["fuc_V"] = sy.lambdify(i_d, sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ) )
                    erg["status"] = 'Error: Voltage limit reached!'
                    return(erg)            
            with np.errstate(invalid='ignore'):
                A = np.diff(np.sign(func(xx)))
            A[np.isnan(A)] = 0
            test = xx[np.where(A)[0]]
            
            i_ds = test[0]
            erg["i_d"] = i_ds
            i_q = ( (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*list(test)[0] )) )
            erg["i_q"] = i_q
            erg["status"] = 'Status: Field-Weakening!'
            v_d = (i_ds * vehicle['electric']['Rs'] - omega*vehicle['electric']['Lq'] * i_q).evalf()
            erg["v_d"] = v_d
            v_q = (i_q * vehicle['electric']['Rs'] + (vehicle['electric']['Ld'] * i_ds + vehicle['electric']['PSI'])* omega).evalf()
            erg["v_q"] = v_q
            if verbose:
                erg["fuc_V"] = sy.lambdify(i_d, -sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ) )
            if i_ds**2 + i_q**2  > vehicle['electric']['I_max']**2:
                erg["status"] = 'Error: Current limit under voltage limit reached!'
                return(erg)
            return(erg)

    erg["status"] = 'Status: Basic Speed Range'
    return(erg)

def plot_PMSM(vehicle, T,n):
    """Plot the PMSM set point in the dq-current plane."""
    erg = calc_PMSM(vehicle, T, n, True)
    xx = np.linspace (-700,-1,1000)
    fig, ax = plt.subplots(figsize=(6,5))
    with np.errstate(invalid='ignore'):
        ax.plot(xx,erg["fuc_T"](xx),'g', alpha = 1, label='Tourque Curve')
        ax.plot(erg['MTPC'][0], erg['MTPC'][1],'ro', label='MTPC')
        if erg['status'] == 'Status: Field-Weakening!' or erg['status'] == 'Error: Voltage limit reached!' or erg['status'] == 'Error: Current limit under voltage limit reached!':
                ax.plot(erg['i_d'], erg['i_q'],'bo', label = 'MTPV')
                ax.plot(xx, erg['fuc_V'](xx), 'b', label = 'Voltage Boundary')
        
        ax.plot(xx, erg['fuc_Imax'](xx) ,'r', label='I_max Boundary' )
        ax.set_xlim(-700, 0)
        #ax.set_ylim(0,700)
    ax.grid(True, axis='both', ls="-", color='0.8')
    ax.set_ylabel('i_q')
    ax.set_xlabel('i_d')
    ax.legend()
    title_text = vehicle['name'] + ' @ ' + str(int(T)) + ' Nm ' + str(int(n*60)) + ' rmp ' + '\n' + erg['status']
    ax.set_title( title_text )
    # plt.show()
    
def calc_SVPWM(vd, vq, n, vehicle):
    """Calculation of the Space Vetor Pulse Width Modulation (SVPWM) for the inverter."""
    f = n * vehicle['electric']['zp']
    t_dead = vehicle['electric']['t_dead']
    Vdc = vehicle['electric']['V_DC']
    fs = vehicle['electric']['fs']
    betrag = np.sqrt(vd**2 + vq**2)
    M = betrag / Vdc * 2

    if betrag > Vdc/np.sqrt(3):
        Status = 'Warning: reached maximum Voltage -> overmodulation or worse'
    else:
        Status = 'Status: Normal'

    if t_dead == 0:
        t = np.arange(0,1/f, 1/fs/100)
    else:
        t = np.arange(0,1/f, t_dead)

    w = 2 * np.pi * f
    wt = w*t
    va,vb,vc = dq0_to_abc(vd,vq,wt)

    vaz = va - ((np.maximum.reduce([va,vb,vc]) + np.minimum.reduce([va,vb,vc]))/2) #Holmes (6.33)
    vbz = vb - ((np.maximum.reduce([va,vb,vc]) + np.minimum.reduce([va,vb,vc]))/2) #Holmes (6.33)
    vcz = vc - ((np.maximum.reduce([va,vb,vc]) + np.minimum.reduce([va,vb,vc]))/2) #Holmes (6.33)

    triangle = Vdc/2 * signal.sawtooth(2 * np.pi * fs * t, 0.5)

    B0wf = np.greater(vaz, triangle).astype(int)
    B0wf[np.where(B0wf == 0)] = -1
    if t_dead != 0:
        B0wf[1:-1:2][np.equal(B0wf[:-2:2], -B0wf[2::2])] = 0

    B1wf = np.greater(vbz, triangle).astype(int)
    B1wf[np.where(B1wf == 0)] = -1
    if t_dead != 0: 
        B1wf[1:-1:2][np.equal(B1wf[:-2:2], -B1wf[2::2])] = 0

    B2wf = np.greater(vcz, triangle).astype(int)
    B2wf[np.where(B2wf == 0)] = -1
    if t_dead != 0: 
        B2wf[1:-1:2][np.equal(B2wf[:-2:2], -B2wf[2::2])] = 0
        
    return(t, M, B0wf, B1wf, B2wf, Status)

def MOSFET_losses(vehicle, f, time, B0wf, B1wf, B2wf, i_d, i_q):

    """Calculate the MOSFET losses: 
    Conduction losses including parallel conduction of MOSFET channel and body diode in reverese conduction; 
    dead/blanking time conduction losses of diodes; 
    switching losses bases on Eon and Eoff including Gate ON and OFF resistors;
    Reverse recovery losses of body diode"""
   
    Temp = vehicle['electric']['Tj']
    Rg_on = vehicle['electric']['Rg_on']
    Rg_off = vehicle['electric']['Rg_off']
    order = 2
    Vdcref = vehicle['electric']['Vdcref']
    data_folder = Path('datasheets/' + vehicle['electric']['foulder_MOSFET'])
    t_dead = vehicle['electric']['t_dead']
    
    fs = vehicle['electric']['fs']
    Vdc = vehicle['electric']['V_DC']
    w = 2 * np.pi * f
    wt = w*time
    ia,ib,ic = dq0_to_abc(i_d, i_q, wt)
    
    file = data_folder /'Rds_on_400A_T.csv'
    Rds_on_400A_T = read_csv(file)
    file = data_folder / 'Rds_on_25C_Id.csv'
    Rds_on_25C_Id = read_csv(file)
    file = data_folder /'Isd_T25_Vsd.csv'
    Isd_T25_Vsd = read_csv(file)
    Isd_T25_Vsdp = poly.Polynomial(poly.polyfit(Isd_T25_Vsd[:,1],Isd_T25_Vsd[:,0],order))
    # maybe an exponential fit like: https://swharden.com/blog/2020-09-24-python-exponential-fit/#exponential-fit-with-python would be better for extrapolation, interpolation should be fine
    Rds_on_400A_T = poly.Polynomial(poly.polyfit(Rds_on_400A_T[:,0],Rds_on_400A_T[:,1],order))
    Rds_on_25C_Idp = poly.Polynomial(poly.polyfit(Rds_on_25C_Id[:,0],Rds_on_25C_Id[:,1],order-1))
    
    file = data_folder / 'Eoff_T125_Rg.csv'
    Eoff_Rg = read_csv(file)
    file = data_folder / 'Eon_T125_Rg.csv'
    Eon_Rg = read_csv(file)

    Eoff_Rg_p = poly.Polynomial(poly.polyfit(Eoff_Rg[:,0],Eoff_Rg[:,1],order))
    Eon_Rg_p = poly.Polynomial(poly.polyfit(Eon_Rg[:,0],Eon_Rg[:,1],order))
    
    # Erec
    file = data_folder / 'Erec_T125.csv'
    Erec_T125 = read_csv(file)
    file = data_folder / 'Erec_Rg.csv'
    Erec_Rg = read_csv(file)

    Erec_Rg_p = poly.Polynomial(poly.polyfit(Erec_Rg[:,0],Erec_Rg[:,1],14))
    Erec_T125_p = poly.Polynomial(poly.polyfit(Erec_T125[:,0],Erec_T125[:,1],order))
    
    A = Rds_on_25C_Idp.coef[1]
    B = Rds_on_25C_Idp.coef[0]
    C = Isd_T25_Vsdp.coef[2]
    D = Isd_T25_Vsdp.coef[1]
    E = Isd_T25_Vsdp.coef[0]

    def calculation(B0wf):
        if t_dead != 0:
            Wcon_MOSh = np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )])**2 * Rds_on_25C_Idp(np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * t_dead
            Wcon_MOSl = np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )])**2 * Rds_on_25C_Idp(np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25)* t_dead
            
            Wcon_MOSh = np.append(Wcon_MOSh, np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == 1) & (ia < 0) )])**2 
                                  * Rds_on_25C_Idp(np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == 1) & (ia < 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * t_dead, axis=0)
            I_MOShp = np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) > Isd_T25_Vsdp(0) ) & (B0wf == 1) & (ia < 0) )]) 
            I_MOSh = -((2*C*I_MOShp + D+ B) / (A-C))/2 + np.sqrt((((2*C*I_MOShp + D+ B) / (A-C))/2)**2 - (-(C * I_MOShp**2 + E + D*I_MOShp)  / (A-C)))
            I_MOS_dh = I_MOShp - I_MOSh
            Wcon_MOSh = np.append(Wcon_MOSh, I_MOSh**2 *  Rds_on_25C_Idp(I_MOSh) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * t_dead)
            Wcon_MOS_dh = I_MOS_dh * Isd_T25_Vsdp(I_MOS_dh)
            Wcon_MOSl = np.append(Wcon_MOSl, np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == -1) & (ia > 0) )])**2 \
                                  * Rds_on_25C_Idp(np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == -1) & (ia > 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * t_dead, axis=0)
            I_MOSlp = np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) > Isd_T25_Vsdp(0) ) & (B0wf == -1) & (ia > 0) )])
            I_MOSl = -((2*C*I_MOSlp + D+ B) / (A-C))/2 + np.sqrt((((2*C*I_MOSlp + D+ B) / (A-C))/2)**2 - (-(C * I_MOSlp**2 + E + D*I_MOSlp)  / (A-C)))
            I_MOS_dl = I_MOSlp - I_MOSl
            Wcon_MOSl = np.append(Wcon_MOSl, I_MOSl**2 *  Rds_on_25C_Idp(I_MOSl) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * t_dead)
            Wcon_MOS_dl = I_MOS_dl * Isd_T25_Vsdp(I_MOS_dl)
        else:
            Wcon_MOSh = np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )])**2 * Rds_on_25C_Idp(np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * 1/fs/100
            Wcon_MOSl = np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )])**2 * Rds_on_25C_Idp(np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25)* 1/fs/100
            
            Wcon_MOSh = np.append(Wcon_MOSh, np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == 1) & (ia < 0) )])**2 
                                  * Rds_on_25C_Idp(np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == 1) & (ia < 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * 1/fs/100, axis=0)
            I_MOShp = np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) > Isd_T25_Vsdp(0) ) & (B0wf == 1) & (ia < 0) )]) 
            I_MOSh = -((2*C*I_MOShp + D+ B) / (A-C))/2 + np.sqrt((((2*C*I_MOShp + D+ B) / (A-C))/2)**2 - (-(C * I_MOShp**2 + E + D*I_MOShp)  / (A-C)))
            I_MOS_dh = I_MOShp - I_MOSh
            Wcon_MOSh = np.append(Wcon_MOSh, I_MOSh**2 *  Rds_on_25C_Idp(I_MOSh) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * 1/fs/100)
            Wcon_MOS_dh = I_MOS_dh * Isd_T25_Vsdp(I_MOS_dh)
            Wcon_MOSl = np.append(Wcon_MOSl, np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == -1) & (ia > 0) )])**2 \
                                  * Rds_on_25C_Idp(np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) <= Isd_T25_Vsdp(0) ) & (B0wf == -1) & (ia > 0) )])) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * 1/fs/100, axis=0)
            I_MOSlp = np.abs(ia[np.where( ( np.abs(ia) * Rds_on_25C_Idp(ia) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) > Isd_T25_Vsdp(0) ) & (B0wf == -1) & (ia > 0) )])
            I_MOSl = -((2*C*I_MOSlp + D+ B) / (A-C))/2 + np.sqrt((((2*C*I_MOSlp + D+ B) / (A-C))/2)**2 - (-(C * I_MOSlp**2 + E + D*I_MOSlp)  / (A-C)))
            I_MOS_dl = I_MOSlp - I_MOSl
            Wcon_MOSl = np.append(Wcon_MOSl, I_MOSl**2 *  Rds_on_25C_Idp(I_MOSl) * Rds_on_400A_T(Temp)/Rds_on_400A_T(25) * 1/fs/100)
            Wcon_MOS_dl = I_MOS_dl * Isd_T25_Vsdp(I_MOS_dl)
          
        if t_dead != 0:
            Wcon_d = np.abs(ia[np.where(B0wf == 0)]) * Isd_T25_Vsdp(np.abs(ia[np.where(B0wf == 0)])) * t_dead 
        else: 
            Wcon_d = 0
        
        P_cond = ( np.sum(Wcon_MOSh) + np.sum(Wcon_MOSl) + np.sum(Wcon_MOS_dl) + np.sum(Wcon_MOS_dh) + np.sum(Wcon_d))*f
        
        file = data_folder / 'Eoff_ Tvj _ 125_C.csv'
        Eoff = read_csv(file)
        file = data_folder / 'Eon_ Tvj _ 125_C.csv'
        Eon = read_csv(file)
    
        if Eon[0,0] != 0:
            Eon = np.insert(Eon, [0], [0], axis=0)
        if Eoff[0,0] != 0:
            Eoff = np.insert(Eoff, [0], [0], axis=0)
    
   
        wEoff = np.ones(Eoff[:,1].shape)
        wEoff[0] = 10 # force Eoffp(0) closer to zero
        wEon = np.ones(Eon[:,1].shape)
        wEon[0] = 10 # force Eonp(0) closer to zero
    
    
        Eoffp = poly.Polynomial(poly.polyfit(Eoff[:,0],Eoff[:,1],order,w = np.sqrt(wEoff)))
        Eonp = poly.Polynomial(poly.polyfit(Eon[:,0],Eon[:,1],order,w = np.sqrt(wEon)))
        
        if t_dead != 0:
            T1hon = Eonp(ia[np.where((B0wf[1::] ==0) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Eon_Rg_p(Rg_on) / Eon_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T1h off
            T1hoff = Eoffp(ia[np.where((B0wf[1::] ==1) & (B0wf[:-1:] == 0)& (ia[1::] > 0))]) * Eoff_Rg_p(Rg_off) / Eoff_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T1h on
            T2loff = Eoffp(np.abs(ia[np.where((B0wf[1::] ==0) & (B0wf[:-1:] == -1)& (ia[1::] < 0))])) * Eoff_Rg_p(Rg_off) / Eoff_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T2l off
            T2lon = Eonp(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 0) & (ia[1::] < 0))])) * Eon_Rg_p(Rg_on) / Eon_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T2l on
            #Erec
            D1hoff = Erec_T125_p(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 0) & (ia[1::] < 0))])) * Erec_Rg_p(Rg_on)/Erec_Rg_p(5.1)
            D2loff = Erec_T125_p(ia[np.where((B0wf[1::] ==0) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Erec_Rg_p(Rg_on)/Erec_Rg_p(5.1)
        else:
            T1hon = Eonp(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Eon_Rg_p(Rg_on) / Eon_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T1h off
            T1hoff = Eoffp(ia[np.where((B0wf[1::] ==1) & (B0wf[:-1:] == -1)& (ia[1::] > 0))]) * Eoff_Rg_p(Rg_off) / Eoff_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T1h on
            T2loff = Eoffp(np.abs(ia[np.where((B0wf[1::] ==1) & (B0wf[:-1:] == -1)& (ia[1::] < 0))])) * Eoff_Rg_p(Rg_off) / Eoff_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T2l off
            T2lon = Eonp(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] < 0))])) * Eon_Rg_p(Rg_on) / Eon_Rg_p(5.1) * (Vdc/Vdcref)**1.4 # T2l on
            #Erec
            D1hoff = Erec_T125_p(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] < 0))])) * Erec_Rg_p(Rg_on)/Erec_Rg_p(5.1)
            D2loff = Erec_T125_p(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Erec_Rg_p(Rg_on)/Erec_Rg_p(5.1)
        
        P_sw = ( np.sum(T1hon) + np.sum(T1hoff) + np.sum(T2loff) + np.sum(T2lon) + np.sum(D1hoff) + np.sum(D2loff) ) *f
    
      
        P_tot = P_cond + P_sw
        
        return(P_tot, P_cond, P_sw)
    
    (P_tot, P_cond, P_sw) = [sum(x) for x in zip( calculation(B0wf), calculation(B1wf), calculation(B2wf))]
    return(P_tot, P_cond, P_sw)

def IGBT_losses(vehicle, f, time, B0wf, B1wf, B2wf, i_d, i_q):
    """Calculation of IGBT losses: Conduction losses of IGBT and diode;
    Switching losses based on Eon and Eoff depenting on Ron and Roff;
    reverse recovery losses of diode"""

    #Temp = vehicle['electric']['Tj'] # tbt include temperature dependency to IGBT losses
    Rg_on = vehicle['electric']['Rg_on']
    Rg_off = vehicle['electric']['Rg_off']
    order = 2
    Vdcref = vehicle['electric']['Vdcref']
    data_folder = Path('datasheets/' + vehicle['electric']['foulder_IGBT'])
    t_dead = vehicle['electric']['t_dead']
    

    fs = vehicle['electric']['fs']
    Vdc = vehicle['electric']['V_DC']
    w = 2 * np.pi * f
    wt = w*time
    ia,ib,ic = dq0_to_abc(i_d, i_q, wt)
    
    file = data_folder / 'Vce_Ic_T25.csv'
    # Vce_Ic_T150.csv
    # Vce_Ic_T125.csv
    Vce_Ic = read_csv(file)
    Vce_Ic_p = poly.Polynomial(poly.polyfit(Vce_Ic[:,1],Vce_Ic[:,0],order))
    file = data_folder / 'If_Vf_T25.csv'
    If_Vf = read_csv(file)
    Vf_If_p = poly.Polynomial(poly.polyfit(If_Vf[:,1],If_Vf[:,0],order))
    
    file = data_folder / 'Eon_Tvj125.csv'
    Eon_Tvj125 = read_csv(file)
    if Eon_Tvj125[0,0] != 0:
        Eon_Tvj125 = np.insert(Eon_Tvj125, [0], [0], axis=0)
    wEon = np.ones(Eon_Tvj125[:,1].shape)
    wEon[0] = 10 # force Eonp(0) closer to zero
    Eon_Tvj125_p = poly.Polynomial(poly.polyfit(Eon_Tvj125[:,0],Eon_Tvj125[:,1],order, w = np.sqrt(wEon)))
    file = data_folder / 'Eon_Rg_Tvj125.csv'
    Eon_Rg_Tvj125 = read_csv(file)
    Eon_Rg_Tvj125_p = poly.Polynomial(poly.polyfit(Eon_Rg_Tvj125[:,0],Eon_Rg_Tvj125[:,1],order))
    file = data_folder / 'Eoff_Tvj125.csv'
    Eoff_Tvj125 = read_csv(file)
    if Eoff_Tvj125[0,0] != 0:
        Eoff_Tvj125 = np.insert(Eoff_Tvj125, [0], [0], axis=0)
    wEon = np.ones(Eoff_Tvj125[:,1].shape)
    wEon[0] = 10 # force Eonp(0) closer to zero
    Eoff_Tvj125_p = poly.Polynomial(poly.polyfit(Eoff_Tvj125[:,0],Eoff_Tvj125[:,1],order, w = np.sqrt(wEon)))
    file = data_folder / 'Eoff_Rg_Tvj125.csv'
    Eoff_Rg_Tvj125 = read_csv(file)
    Eoff_Rg_Tvj125_p = poly.Polynomial(poly.polyfit(Eoff_Rg_Tvj125[:,0],Eoff_Rg_Tvj125[:,1],order))
    
    # Erec
    file = data_folder / 'Erec_Tvj125.csv'
    Erec_Tvj125 = read_csv(file)
    if Erec_Tvj125[0,0] != 0:
        Erec_Tvj125 = np.insert(Erec_Tvj125, [0], [0], axis=0)
    wEon = np.ones(Erec_Tvj125[:,1].shape)
    wEon[0] = 10 # force Eonp(0) closer to zero
    Erec_Tvj125_p = poly.Polynomial(poly.polyfit(Erec_Tvj125[:,0],Erec_Tvj125[:,1],order, w = np.sqrt(wEon)))
    file = data_folder / 'Erec_Rg_Tvj125.csv'
    Erec_Rg_Tvj125 = read_csv(file)
    Erec_Rg_Tvj125_p = poly.Polynomial(poly.polyfit(Erec_Rg_Tvj125[:,0],Erec_Rg_Tvj125[:,1],order))

    def calculation(B0wf):

        if t_dead != 0:
            Wcon_igbth = np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )]) * Vce_Ic_p(np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )])) * t_dead 
            Wcon_dl = np.abs(ia[np.where( (B0wf == -1) & (ia > 0) )]) * Vf_If_p(np.abs(ia[np.where((B0wf == -1) & (ia > 0) )])) * t_dead 
            Wcon_igbtl = np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )]) * Vce_Ic_p(np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )])) * t_dead 
            Wcon_dh = np.abs(ia[np.where( (B0wf == 1) & (ia < 0) )]) * Vf_If_p(np.abs(ia[np.where((B0wf == 1) & (ia < 0) )])) * t_dead 
        else: 
            Wcon_igbth = np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )]) * Vce_Ic_p(np.abs(ia[np.where( (B0wf == 1) & (ia > 0) )])) * 1/fs/100 
            Wcon_dl = np.abs(ia[np.where( (B0wf == -1) & (ia > 0) )]) * Vf_If_p(np.abs(ia[np.where((B0wf == -1) & (ia > 0) )])) * 1/fs/100 
            Wcon_igbtl = np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )]) * Vce_Ic_p(np.abs(ia[np.where( (B0wf == -1) & (ia < 0) )])) * 1/fs/100
            Wcon_dh = np.abs(ia[np.where( (B0wf == 1) & (ia < 0) )]) * Vf_If_p(np.abs(ia[np.where((B0wf == 1) & (ia < 0) )])) * 1/fs/100 
         
        # conduction losses of diodes during dead time    
        if t_dead != 0:
            Wconb2_igbt_d = np.abs(ia[np.where(B0wf == 0)]) * Vf_If_p(np.abs(ia[np.where(B0wf == 0)])) * t_dead 
        else: 
            Wconb2_igbt_d = 0
    
        Pcon_igbt = ( np.sum(Wcon_igbth) + np.sum(Wcon_igbtl) )*f
        Pcon_igbt_d = ( np.sum(Wcon_dl) + np.sum(Wcon_dh) )*f
        Pcon_igbt_dd = np.sum(Wconb2_igbt_d)*f
        Pcon_igbt_ges = Pcon_igbt + Pcon_igbt_d + Pcon_igbt_dd
        
        if t_dead !=0:
            T1hon = Eon_Tvj125_p(ia[np.where((B0wf[1::] ==0) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Eon_Rg_Tvj125_p(Rg_on) / Eon_Rg_Tvj125_p(2.2) * (Vdc/Vdcref)**1.4 # T1h off
            T1hoff = Eoff_Tvj125_p(ia[np.where((B0wf[1::] ==1) & (B0wf[:-1:] == 0)& (ia[1::] > 0))]) * Eoff_Rg_Tvj125_p(Rg_off) / Eoff_Rg_Tvj125_p(2.2) * (Vdc/Vdcref)**1.4 # T1h on
            T2loff = Eoff_Tvj125_p(np.abs(ia[np.where((B0wf[1::] ==0) & (B0wf[:-1:] == -1)& (ia[1::] < 0))])) * Eoff_Rg_Tvj125_p(Rg_off) / Eoff_Rg_Tvj125_p(2.2) * (Vdc/Vdcref)**1.4 # T2l off
            T2lon = Eon_Tvj125_p(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 0) & (ia[1::] < 0))])) * Eon_Rg_Tvj125_p(Rg_on) / Eon_Rg_Tvj125_p(5.1) * (Vdc/Vdcref)**1.4 # T2l on
            D1hoff = Erec_Tvj125_p(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 0) & (ia[1::] < 0))])) * Erec_Rg_Tvj125_p(Rg_on)/Erec_Rg_Tvj125_p(2.2)
            D2loff = Erec_Tvj125_p(ia[np.where((B0wf[1::] ==0) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Erec_Rg_Tvj125_p(Rg_on)/Erec_Rg_Tvj125_p(2.2)
        else:
            T1hon = Eon_Tvj125_p(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Eon_Rg_Tvj125_p(Rg_on) / Eon_Rg_Tvj125_p(2.2) * (Vdc/Vdcref)**1.4 # T1h off
            T1hoff = Eoff_Tvj125_p(ia[np.where((B0wf[1::] ==1) & (B0wf[:-1:] == -1)& (ia[1::] > 0))]) * Eoff_Rg_Tvj125_p(Rg_off) / Eoff_Rg_Tvj125_p(2.2) * (Vdc/Vdcref)**1.4 # T1h on
            T2loff = Eoff_Tvj125_p(np.abs(ia[np.where((B0wf[1::] ==1) & (B0wf[:-1:] == -1)& (ia[1::] < 0))])) * Eoff_Rg_Tvj125_p(Rg_off) / Eoff_Rg_Tvj125_p(2.2) * (Vdc/Vdcref)**1.4 # T2l off
            T2lon = Eon_Tvj125_p(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] < 0))])) * Eon_Rg_Tvj125_p(Rg_on) / Eon_Rg_Tvj125_p(5.1) * (Vdc/Vdcref)**1.4 # T2l on
            D1hoff = Erec_Tvj125_p(np.abs(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] < 0))])) * Erec_Rg_Tvj125_p(Rg_on)/Erec_Rg_Tvj125_p(2.2)
            D2loff = Erec_Tvj125_p(ia[np.where((B0wf[1::] ==-1) & (B0wf[:-1:] == 1) & (ia[1::] > 0))]) * Erec_Rg_Tvj125_p(Rg_on)/Erec_Rg_Tvj125_p(2.2)
       
        Psw = ( np.sum(T1hon) + np.sum(T1hoff) + np.sum(T2loff) + np.sum(T2lon) + np.sum(D1hoff) + np.sum(D2loff) ) *f
    
        P_tot = Pcon_igbt_ges + Psw
        
        return(P_tot, Pcon_igbt_ges, Psw)
    (P_tot, P_cond, P_sw) = [sum(x) for x in zip( calculation(B0wf), calculation(B1wf), calculation(B2wf))]
    return(P_tot, P_cond, P_sw)

def read_drive_cycle(file, vehicle):
    """This fuction reads the drive cycle of the cycle folder and calculates acceloration, total force, power, rouns per second, torque...:
    
    .. math::
    
        \\omega_{1/s} = \\frac{v_{\\mathrm{m/s}}}{wheel_{radius} \\cdot 2 \\pi} \\cdot  i_{g}
        
        T = \\frac{Force \\cdot wheel_{radius}}{i_{g}} \cdot \\eta_{i_{g}} \\;   \\mathrm{for \\; M => 0}
        
        T = \\frac{Force \\cdot wheel_{radius}}{i_{g}} \cdot \\frac{1}{\\eta_{i_{g}}} \\;   \\mathrm{for \\; M < 0}
    
    """
    with open(file, 'r') as read_obj: # read csv file as a list of lists
        csv_reader = reader(read_obj) # pass the file object to reader() to get the reader object
        list_of_rows = list(csv_reader) # Pass reader object to list() to get a list of lists
        cycle = np.delete(np.asarray(list_of_rows), (0), axis = 0).astype(float)
        a = np.zeros((cycle[:,0].shape[0],24))
        SVPWM = {}
        items = list(range(0,cycle[:,0].shape[0]-1))
        #for j in progressBar(items, prefix = 'Drive Cycle calculation:', suffix = 'Complete', length = 50):
        with click.progressbar(items, label='calculate drive cycle') as bar:
            for j in bar:
                if j == 0:
                    a[j,0] = (cycle[j+1,1] - cycle[j,1])/3.6/(cycle[j+1,0] - cycle[j,0]) # acceleration [m/s**2]
                else:
                    a[j,0] = (cycle[j+1,1] - cycle[j-1,1])/3.6/(cycle[j+1,0] - cycle[j-1,0]) # acceleration [m/s**2]
                a[j,1] = cycle[j,1]/3.6 # speed [m/s]
                a[j,2] = total_force(vehicle, a[j,0], a[j,1], cycle[j,2]) # Force [Nm]
                a[j,3] = a[j,1] * a[j,2] # Power [W]
                a[j,4] = a[j,1] /( vehicle['basic']['wheel_radius'] * 2 * np.pi) * vehicle['basic']['i_g'] # rotational speed / rounds per second [1/s]
                if a[j,0] >= 0:
                    a[j,5] = a[j,2] * vehicle['basic']['wheel_radius'] / vehicle['basic']['i_g'] * vehicle['basic']['my_i_g'] # torque [Nm]
                else:
                    a[j,5] = a[j,2] * vehicle['basic']['wheel_radius'] / vehicle['basic']['i_g'] / vehicle['basic']['my_i_g'] # torque [Nm]
                a[j,6] = 2* np.pi * a[j,4] * a[j,5] # Power M1 [W]
                PMSM = calc_PMSM(vehicle,a[j,5], a[j,4])
                # print(j)
                # print('M [Nm]: ', str(a[j,5]))
                # print('n [1/s]: ', str(a[j,4]))
                if PMSM['status'][0:6] != 'Status':
                    a[j,7] = np.nan
                    a[j,8] = np.nan
                    a[j,9] = np.nan
                    a[j,10] = np.nan
                else:
                    a[j,7] = PMSM['i_d'] # i_d current setpoint of PMSM
                    a[j,8] = PMSM['i_q'] # i_q current setpoint of PMSM
                    a[j,9] = PMSM['v_d'] # v_d current setpoint of PMSM
                    a[j,10] = PMSM['v_q'] # v_q current setpoint of PMSM
                    #a[j,11] = PMSM['status'] # status of PMSM -> I sould store this somewhere...
                if a[j,4] == 0 or a[j,5] == 0:
                    (time, M, B0wf, B1wf, B2wf,Status) = (0,0,0,0,0,'Status: M or n equals null')
                else:
                    (time, M, B0wf, B1wf, B2wf, Status) = calc_SVPWM(float(PMSM['v_d']), float(PMSM['v_q']), a[j,4], vehicle)
                # SVPWM[j] = {}
                # SVPWM[j]['time'] = time
                # SVPWM[j]['M'] = M
                # SVPWM[j]['SV0'] = B0wf
                # SVPWM[j]['SV1'] = B1wf
                # SVPWM[j]['SV2'] = B2wf
                # SVPWM[j]['StatusSVPWM'] = Status
                # SVPWM[j]['StatusPSM'] = PMSM['status']
                if (a[j,4] == 0 or a[j,5] == 0):
                    MOSFET = (0,0,0)
                    IGBT = (0,0,0)
                else:
                    MOSFET = MOSFET_losses(vehicle, a[j,4] * vehicle['electric']['zp'], time, B0wf, B1wf, B2wf, a[j,7], a[j,8])
                    IGBT = IGBT_losses(vehicle, a[j,4] * vehicle['electric']['zp'], time, B0wf, B1wf, B2wf, a[j,7], a[j,8])
                a[j,11] = (np.abs(a[j,7] + a[j,8]*1j) / np.sqrt(2)) **2 * vehicle['electric']['Rs'] # PMSM Loss Stator [W]
                a[j,12] = MOSFET[0] # MOSFET total losses
                a[j,13] = MOSFET[1] # MOSFET conduction losses
                a[j,14] = MOSFET[2] # MOSFET switching losses
                if np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) == 0:
                    a[j,15] = np.nan
                elif ((np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) -a[j,12] )  / np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) >= 1 or (np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) -a[j,12] )  / np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) < 0 ):
                    a[j,15] = np.nan
                else:
                    a[j,15] =  (np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) -a[j,12] )  / np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j))  # Energy conversion efficiency
                a[j,17] = IGBT[0] #  IGBT total losses
                a[j,18] = IGBT[1] # IGBT conduction losses
                a[j,19] = IGBT[2] # IGBT switching losses
                if np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) == 0:
                    a[j,20] = np.nan
                elif ((np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) -a[j,17] )  / np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) >= 1 or (np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) -a[j,17] )  / np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) < 0 ):
                    a[j,20] = np.nan
                else:
                    a[j,20] =  (np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j)) -a[j,17] )  / np.abs( np.real( a[j,9] + a[j,10]*1j) * np.conjugate(a[j,7] + a[j,8]*1j))  # Energy conversion efficiency
                
            a[:,16] = np.nancumsum(a[:,12])/3600000 #E_total_MOSFET [kWh]
            a[:,21] = np.nancumsum(a[:,17])/3600000 # E_total_IGBT [kWh]
            a[:,22] = np.nancumsum(np.abs(a[:,3]))/3600000 #abs(E_total) [kWh]
            a[:,23] = np.nancumsum(a[:,11])/3600000 #E_total_loss PMSM Stator [kWh]
        cycle = np.concatenate((cycle, a), axis=1 )
        variable_names = ['time','speed [km/h]','slope [°]','acceleration [m/s**2]','speed [m/s]','Force [Nm]','Total Power [W]', 'rounds per second [1/s]',
        'torque [Nm]','Power M1 [W]','i_d [A]','i_q [A]','v_d [V]','v_q [V]', 'PMSM_Loss Stator [W]', 'P_total_MOSFET [W]','P_cond_MOSFET [W]','P_sw_MOSFET [W]', 'ETA_MOSFET', 'E_total_MOSFET [kWh]',
        'P_total_IGBT [W]', 'P_cond_IGBT [W]', 'P_sw_IGBT [W]', 'ETA_IGBT', 'E_total_IGBT [kWh]', 'abs(E_total) [kWh]','E_total_loss PMSM Stator [kWh]']
    return(cycle, variable_names)

def calc_con(vehicle, qns, qTs):
    """Calculation of conture plots in the torque-speed plane."""
    ns = np.linspace(1,float(vehicle['basic']['n_max'])/60,qns) 
    T_max = 3/2 * vehicle['electric']['zp'] * ( (vehicle['electric']['PSI'] * (vehicle['electric']['I_max']/np.sqrt(2)) ) + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*-(vehicle['electric']['I_max']/np.sqrt(2))**2 )
    T_max = np.ceil( T_max/100)*100
    Ts = np.linspace(-T_max, T_max,qTs)

    con_OP = np.empty((len(ns),len(Ts),),dtype=object)
    con_OP.fill(np.nan)
    con_i_d = np.empty((len(ns),len(Ts),))
    con_i_d.fill(np.nan)
    con_i_q = np.empty((len(ns),len(Ts),))
    con_i_q.fill(np.nan)
    con_v_d = np.empty((len(ns),len(Ts),))
    con_v_d.fill(np.nan)
    con_v_q = np.empty((len(ns),len(Ts),))
    con_v_q.fill(np.nan)
    con_P = np.empty((len(ns),len(Ts),))
    con_P.fill(np.nan)

    con_MOSFET_T = np.empty((len(ns),len(Ts),))
    con_MOSFET_T.fill(np.nan)
    con_MOSFET_C = np.empty((len(ns),len(Ts),))
    con_MOSFET_C.fill(np.nan)
    con_MOSFET_S = np.empty((len(ns),len(Ts),))
    con_MOSFET_S.fill(np.nan)
    con_IGBT_T = np.empty((len(ns),len(Ts),))
    con_IGBT_T.fill(np.nan)
    con_IGBT_C = np.empty((len(ns),len(Ts),))
    con_IGBT_C.fill(np.nan)
    con_IGBT_S = np.empty((len(ns),len(Ts),))
    con_IGBT_S.fill(np.nan)


    with click.progressbar(range(len(ns)),label='Tn plot calculating: ') as bar:
        for n in bar:
            # print(ns[n]*60)
            for T in range(len(Ts)):
                # print(ns[n]*60)
                # print(Ts[T])
                erg = calc_PMSM(vehicle, Ts[T], ns[n], False)
                # print(erg['status'])
                if not erg['status'][0:5] == 'Error':
                    con_OP[n,T] = (ns[n],Ts[T])
                    con_i_d[n,T] = float(erg['i_d'])
                    con_i_q[n,T] = float(erg['i_q'])
                    con_v_d[n,T] = float(erg['v_d'])
                    con_v_q[n,T] = float(erg['v_q'])
                    con_P[n,T] = 2*np.pi*Ts[T]*ns[n]
                    (time, M, B0wf, B1wf, B2wf, Status) = calc_SVPWM(float(erg['v_d']),float(erg['v_q']),ns[n],vehicle)
                    # print(Status)
                    if not Status[0:5] == 'Error':
                        (M_T,M_C,M_S) = MOSFET_losses(vehicle, ns[n] * vehicle['electric']['zp'], time, B0wf, B1wf, B2wf, float(erg['i_d']), float(erg['i_q']))
                        con_MOSFET_T[n,T] = M_T
                        con_MOSFET_C[n,T] = M_C
                        con_MOSFET_S[n,T] = M_S
                        (I_T,I_C,I_S) = IGBT_losses(vehicle, ns[n] * vehicle['electric']['zp'], time, B0wf, B1wf, B2wf, float(erg['i_d']), float(erg['i_q']))
                        con_IGBT_T[n,T] = I_T
                        con_IGBT_C[n,T] = I_C
                        con_IGBT_S[n,T] = I_S
                        # if Status == 'Warning: reached maximum Voltage -> overmodulation or worse':
                        #     print(ns[n]*60)
                        #     print(Ts[T])
                        #     print(Status)
    con_values = (ns, Ts, con_OP, con_i_d, con_i_q, con_v_d, con_v_q, con_P, con_MOSFET_T, con_MOSFET_C, con_MOSFET_S, con_IGBT_T, con_IGBT_C, con_IGBT_S)
    con_names = ['ns', 'Ts', 'con_OP', 'con_i_d', 'con_i_q', 'con_v_d', 'con_v_q', 'con_P', 'con_MOSFET_T', 'con_MOSFET_C', 'con_MOSFET_S', 'con_IGBT_T', 'con_IGBT_C', 'con_IGBT_S']
    return con_values, con_names


class Config(object):
    
    def __init__(self):
        self.vebose=False


pass_config= click.make_pass_decorator(Config, ensure=True)

@click.group()
@click.option('--verbose',is_flag=True)
@click.option('--home-directory', type=click.Path())
@pass_config
def cli(config, verbose, home_directory):
    click.clear()
    """This is the main function."""
    config.verbose = verbose
    if home_directory is None:
        home_directory = os.getcwd()
    config.home_directory = home_directory
    click.clear()
    click.echo(click.style(r"""        __-------__
      / _---------_ \
     / /           \ \
     | |           | |
  _  |_|___________|_|  _
 /-\|                 |/-\
| _ |\       4       /| _ |
|(_)| \      2      / |(_)|
|___|__\_____!_____/__|___|
[________|SiCWell|________] 
 ||||    ~~~~~~~~     ||||
 `--'                 `--'""", fg = 'green'))
    



@cli.command()
@click.option('-v', '--vehicle', type=click.File('r'), default = None,
            help ='Here you can specify the vehicle of the vehicles foulder. If none is specified you will be asked to choose.')
@click.option('-c','--cycle', type=click.File('r'), default = None,
            help = 'Here you can specify the drive cycle of the cycle foulder. If none is specified you will be asked to choose.')
@pass_config
def drive_cycle(config, vehicle, cycle):
    """Calculate a drive cycle."""
    if config.verbose:
        click.echo('We are in verbose mode.')
    click.echo('The Working directory is %s' % config.home_directory)
    
    if vehicle == None:
        vehicles = list(os.listdir('vehicles'))
        n = 0
        click.echo('Here is the list of available vehicles.')
        for vehicle in vehicles:
            click.echo(str(n) + ': ' + str(vehicle))
            n = n+1
        nr = int(click.prompt('Please select vehicle [nr]', type=click.Choice([str(i) for i in range(n)])))
        vehicle = read_toml(os.path.join('vehicles', vehicles[nr]))

    #click.echo('You choose: ' + str(vehicle['name']))

    if cycle == None:
        cycles = list(os.listdir('cycles'))
        n = 0
        click.echo('Here is the list of available cycles.')
        for cycle in cycles:
            click.echo(str(n) + ': ' + str(cycle))
            n = n+1
        nur = int(click.prompt('Please select cycle [nr]', type=click.Choice([str(i) for i in range(n)])))
        click.echo('You choose: ' + str(vehicle['name']) + ' and ' + str(cycles[nur] ) )
        
        current_time = datetime.datetime.now() 
        savename = '_' + str(current_time.year) + '_' + str(current_time.month) + '_' + str(current_time.day) + '_' + str(current_time.hour) + '-' + str(current_time.minute)

        filename = 'save/' + str(vehicle['name'])+ '_' + str(cycles[nur])[0:-4] + savename + '.p'
        save_files = [save_file for save_file in list(os.listdir('save')) if str(vehicle['name']) + '_' + str(cycles[nur] )[0:-4] in save_file]

        if save_files:
            click.echo('There are saved files for your configuration, choose one or start a new calculation.')
            n = 0
            for save_file in save_files:
                click.echo(str(n) + ': ' + str(save_file))
                n = n+1
            click.echo(str(n) + ': ' + 'new calculation')
            nr = int(click.prompt('Please select number [nr]', type=click.Choice([str(i) for i in range(n+1)])))
            if int(nr) == n:
                (cycle, variable_names) = read_drive_cycle(os.path.join('cycles', cycles[nur]), vehicle)
                pickle.dump([cycle, variable_names], open(filename,'wb'))
            else:
                [cycle, variable_names] = pickle.load( open(os.path.join('save', save_files[nr]),'rb'))
        else:
            (cycle, variable_names) = read_drive_cycle(os.path.join('cycles', cycles[nur]), vehicle)
            pickle.dump([cycle, variable_names], open(filename,'wb'))

        #plt.ion()
        fonts = pyfiglet.FigletFont.getFonts()
        ascii_banner = pyfiglet.figlet_format("Plotting",font=fonts[100])
        click.clear()
        click.echo(ascii_banner)
        variable_names.append('Info')
        n = 0
        click.echo(str(n) + ': ' + 'Exit and close all plots')
        n = n+1
        for v_name in variable_names[1:-1]:
            click.echo(str(n) + ': ' + v_name + ' / time(s)' )
            n = n+1
        click.echo(str(n) + ': ' + variable_names[-1])
        nr = int(click.prompt('Please select number of plot', type=click.Choice([str(i) for i in range(n+1)])))
        plt.ion()
        while nr != 0:

        #for nr in range(len(variable_names)-1):
            if variable_names[nr] != 'Info':
                fig, ax = plt.subplots(figsize=(6,5))
                p1 = ax.plot(cycle[:,0], cycle[:,nr], label = variable_names[nr])
                # p2 = ax.bar(labels, V_inv_sim, label = 'Ploss inverter sim', alpha = .6)
                # p3 = ax.bar(labels, V_s_sim, label='Ploss static sim', alpha = .6)
                # p4 = ax.bar(labels, V_d_sim, bottom = V_s_sim, label='Ploss dynamic sim', alpha = .6)
                ax.legend(loc=2)
                # if (nr == 15 or nr == 20):
                #     ax.set_ylim( (0, max( np.nanmax(cycle[:,15]), np.nanmax(cycle[:,20]) ) ) )
                if (nr == 16 or nr == 21):
                    ax.set_ylim( (0, max( np.nanmax(cycle[:,16]), np.nanmax(cycle[:,21]) ) ) )
                if (nr == 17 or nr == 22):
                    ax.set_ylim( (0, max( np.nanmax(cycle[:,17]), np.nanmax(cycle[:,22]) ) ) )
                if (nr == 19 or nr == 24 or nr==26):
                    ax.set_ylim( (0, max( np.nanmax(cycle[:,19]), np.nanmax(cycle[:,24]), np.nanmax(cycle[:,26]) ) ) )
                #ax.set_xlim(43.5,72.5)
                #ax.bar_label(p1,fmt='%0.4f', label_type='edge')
                # ax.bar_label(p2,fmt='%0.0f', label_type='edge')
                # ax.bar_label(p3,fmt='%0.0f', label_type='center')
                # ax.bar_label(p4,fmt='%0.0f', label_type='center')
                ax.set_ylabel(variable_names[nr])
                ax.set_xlabel('Time [s]')
                ax.grid(True, axis='both', ls="-", color='0.8')
                fig.tight_layout()
                #plt.show()
                click.clear()
                click.echo(ascii_banner)
            else:
                Info_banner = pyfiglet.figlet_format("Info",font=fonts[100])
                click.clear()
                click.echo(Info_banner)
                click.echo('You choose: ' + str(vehicle['name']) + ' and ' + str(cycles[nur] ))
                click.echo('The cycle takes: ' + str(cycle[:,0][-1]) + ' s')
                click.echo('The distance is: ' + str(np.cumsum(cycle[:,4])[-1]/1000) + ' km')
                click.echo('The average speed is: ' + str(np.average(abs(cycle[:,1]))) + ' km/h')
                click.echo('The max/min acceleration is: ' + str(np.max((cycle[:,3]))) + ' / ' + str(np.min((cycle[:,3]))) + ' m/s**2')
                click.echo('The max slope is: ' + str(np.max((cycle[:,2]))) + ' °')
                click.echo('The average power is: ' + str(np.average(abs(cycle[:,6]))/1000) + ' kW')
                click.echo('The max speed is: ' + str(np.max(cycle[:,1])) + ' km/h')
                click.echo('The max force is: ' + str(np.max(cycle[:,5])) + ' Nm')
                click.echo('The max speed of the PMSM is: ' +  str(np.max(cycle[:,7])*60) + ' rmp')
                click.echo('The max tourque of the PMSM is: ' + str(np.max(cycle[:,8])) + ' Nm')
                click.echo('The drive cycle efficiency of the SiC-MOSFET inverter is: ' + str(cycle[:,25][-1] / (cycle[:,25][-1] + cycle[:,19][-1])))
                click.echo('The total losses of the SiC-MOSFET: ' + str(cycle[:,19][-1]) + ' kWh')
                click.echo('The drive cycle efficiency of the IGBT inverter is: ' + str(cycle[:,25][-1] / (cycle[:,25][-1] + cycle[:,24][-1])))
                click.echo('The total losses of the Si-IGBT: ' + str(cycle[:,24][-1]) + ' kWh')
                click.echo('Delta Eta SiC-MOSFET vs Si-IGBT: ' + str(cycle[:,25][-1] / (cycle[:,25][-1] + cycle[:,19][-1]) - cycle[:,25][-1] / (cycle[:,25][-1] + cycle[:,24][-1])))
                click.pause('Press any key to continue.')
            n = 0
            click.echo(str(n) + ': ' + 'Exit and close all plots')
            n = n+1
            for v_name in variable_names[1:-1]:
                click.echo(str(n) + ': ' + v_name + ' / time(s)')
                n = n+1
            click.echo(str(n) + ': ' + variable_names[-1])
            nr = int(click.prompt('Please select number of plot', type=click.Choice([str(i) for i in range(n+1)])))
        #click.echo(variable_names)
        click.echo('simulation done')
        plt.close('all')    

@cli.command()
@click.option('-v', '--vehicle', type=click.File('r'), default = None,
            help ='Here you can specify the vehicle of the vehicles foulder. If none is specified you will be asked to choose.')
@click.option('-T', '--torque', type=float, default = None,
            help = 'Torque value for the operation point. If none is specified you will be asked to choose.')
@click.option('-n', '--speed', type=float, default = None,
            help = 'Speed in rounds per minute for the operation point. If none is specified you will be asked to choose.')
#@pass_config
def plot_OP(vehicle, torque, speed):
    """Here you can plot one operation point for the permanet-magnet synchronous motor (PMSM)."""
    if vehicle == None:
        vehicles = list(os.listdir('vehicles'))
        n = 0
        click.echo('Here is the list of available vehicles.')
        for vehicle in vehicles:
            click.echo(str(n) + ': ' + str(vehicle))
            n = n+1
        nr = int(click.prompt('Please select vehicle [nr]', type=click.Choice([str(i) for i in range(n)])))
        vehicle = read_toml(os.path.join('vehicles', vehicles[nr]))
    click.echo('You choose: ' + str(vehicle['name']))
    
    if torque == None:
        torque = click.prompt('Enter torque [Nm] for operation point', type = float)
    if speed == None:
        speed = click.prompt('Enter speed [rpm] for operation point', type = float)
    click.echo('Close plot windows to resume.')
    plot_PMSM(vehicle, torque, speed/60)
    plt.show()

@cli.command()
@click.option('-v', '--vehicle', type=click.File('r'), default = None,
            help ='Here you can specify the vehicle of the vehicles foulder. If none is specified you will be asked to choose.')
@click.option('-qts', '--quantity_t', type=int, default = 90,
            help ='Here you can specify the quantity of torque steps. The default value is 50.')
@click.option('-qns', '--quantity_n', type=int, default = 100,
            help ='Here you can specify the quatity of speed steps. The default value is 100.')
#@pass_config
def plot_Tn(vehicle, quantity_n, quantity_t):
    """Here you can plot the semiconductor losses (total, conduction and switching) on the speed-torque plane of the permanet-magnet synchronous motor (PMSM) (contourplot)."""
    fonts = pyfiglet.FigletFont.getFonts()
    ascii_banner = pyfiglet.figlet_format("s-T Plots",font=fonts[100])
    click.clear()
    click.echo(ascii_banner)
    if vehicle == None:
        vehicles = list(os.listdir('vehicles'))
        n = 0
        click.echo('Here is the list of available vehicles.')
        for vehicle in vehicles:
            click.echo(str(n) + ': ' + str(vehicle))
            n = n+1
        nr = int(click.prompt('Please select vehicle [nr]', type=click.Choice([str(i) for i in range(n)])))
        vehicle = read_toml(os.path.join('vehicles', vehicles[nr]))
    click.echo('You choose: ' + str(vehicle['name']))

    current_time = datetime.datetime.now() 
    savename = '_' + str(current_time.year) + '_' + str(current_time.month) + '_' + str(current_time.day) + '_' + str(current_time.hour) + '-' + str(current_time.minute)

    filename = 'save/' + str(vehicle['name'])+ '_T_' + str(quantity_t) + '_n_' + str(quantity_n) + savename + '.p'
    save_files = [save_file for save_file in list(os.listdir('save')) if str(vehicle['name']) + '_T_' + str(quantity_t) + '_n_' + str(quantity_n) in save_file]

    if save_files:
        click.echo('There are saved files for your configuration, choose one or start a new calculation.')
        n = 0
        for save_file in save_files:
            click.echo(str(n) + ': ' + str(save_file))
            n = n+1
        click.echo(str(n) + ': ' + 'new calculation')
        nr = int(click.prompt('Please select number [nr]', type=click.Choice([str(i) for i in range(n+1)])))
        if int(nr) == n:
            (con_values, con_names) = calc_con(vehicle, quantity_n, quantity_t)
            pickle.dump((con_values, con_names), open(filename,'wb'))
        else:
            (con_values, con_names) = pickle.load( open(os.path.join('save', save_files[nr]),'rb'))
    else:
        (con_values, con_names) = calc_con(vehicle, quantity_n, quantity_t)
        pickle.dump((con_values, con_names), open(filename,'wb'))
    # calc_Tn(vehicle, qTs, qns)    
    # (ns, Ts, con_OP, con_i_d, con_i_q, con_v_d, con_v_q, con_P, con_MOSFET_T, con_Mosfet_C, con_Mosfet_S, con_IGBT_T, con_IGBT_C, con_IGBT_S) = con_values
    n = 0
    click.echo(str(n) + ': ' + 'Exit and close all plots')
    n = n+1
    for v_name in con_names[3:]:
        click.echo(str(n) + ': ' + v_name )
        n = n+1
    click.echo(str(n) + ': Efficiency MOSFET')
    n = n+1
    click.echo(str(n) + ': Efficiency IGBT')
    n = n+1
    click.echo(str(n) + ': Delta Efficiency MOSFET/IGBT')
    n = n+1
    # click.echo(str(n) + ': ' + variable_names[-1])
    nr = int(click.prompt('Please select number of plot', type=click.Choice([str(i) for i in range(n)])))
    plt.ion()
    while nr != 0:

        if nr == 12:
            fig,ax = plt.subplots()
            levels1 = np.linspace(0.9,1,101)
            c1 = ax.contourf(con_values[0]*60, con_values[1], np.abs(con_values[7].T) / (np.abs(con_values[7].T) + np.abs(con_values[8].T)), levels=levels1, cmap='gnuplot', extend = 'min'  )
            c1.cmap.set_under('black')
            bar = fig.colorbar(c1, label='Eta')
            bar.set_ticks(list(np.linspace(0.9,1,10)))
            ax.set_xlabel('n [rpm]')
            ax.set_ylabel('T [Nm]')
            ax.set_title('Efficiency MOSFET')
        elif nr == 13:
            fig,ax = plt.subplots()
            levels2 = np.linspace(0.9,1,101)