classdef BrushlessDCWheelMotor < SoModel
    
    methods
        
        % declarative part
        
        function obj = BrushlessDCWheelMotor
            
            % create object and declare name and version
            obj = obj@SoModel('Brushless DC Wheel Motor','Sophemis 4 Model');
            
            % declare inputs of model
            obj.addInput('C',struct('name','C','comment','Couple electromagnétique','unit','N.m','type','continuous'));
            obj.addInput('omega',struct('name','omega','comment','vitesse de rotation nominale','unit','tr/min','type','continuous'));
            obj.addInput('omegamax',struct('name','omegamax','comment','vitesse de rotation maxinale','unit','tr/min','type','continuous'));
            obj.addInput('kr',struct('name','kr','comment','coefficient de remplissage','unit','-','type','continuous'));
            obj.addInput('kfoi',struct('name','kfoi','comment','coefficient de foisonnement','unit','-','type','continuous'));
            obj.addInput('kfui',struct('name','kfui','comment','coefficient de fuite entre aimants','unit','-','type','continuous'));
            obj.addInput('rocu',struct('name','rocu','comment','resistivite du conducteur à 0°C','unit','ohm.m','type','continuous'));
            obj.addInput('Br',struct('name','Br','comment','induction magnétique rémanente','unit','T','type','continuous'));
            obj.addInput('Bc',struct('name','Bc','comment','induction magnétique de démagnétisation','unit','T','type','continuous'));
            obj.addInput('dt',struct('name','dt','comment','densité du stator','unit','kg/m3','type','continuous'));
            obj.addInput('da',struct('name','da','comment','densité des aimants','unit','kg/m3','type','continuous'));
            obj.addInput('dcu',struct('name','dcu','comment','densité du conducteur','unit','kg/m3','type','continuous'));
            obj.addInput('dcr',struct('name','dcr','comment','densité de la culasse rotor','unit','kg/m3','type','continuous'));
            obj.addInput('qt',struct('name','qt','comment','pertes fer specifiques stator','unit','W/kg','type','continuous'));
            obj.addInput('ft',struct('name','ft','comment','frequence des pertes spécifiques','unit','Hz','type','continuous'));
            obj.addInput('Bt',struct('name','Bt','comment','induction magnétique spécifique','unit','T','type','continuous'));
            obj.addInput('mua',struct('name','mua','comment','perméabilité magnétique relative des aimants','unit','-','type','continuous'));
            obj.addInput('Pm',struct('name','Pm','comment','pertes mécaniques','unit','W','type','continuous'));
            obj.addInput('alphaa',struct('name','alphaa','comment','variation de resistivite des aimants','unit','°C','type','continuous'));
            obj.addInput('alphacu',struct('name','alphacu','comment','variation de resistivite du conducteur','unit','°C','type','continuous'));
            obj.addInput('Text',struct('name','Text','comment','température extérieure','unit','°C','type','continuous'));
            obj.addInput('h',struct('name','h','comment','coefficient de convection','unit','W/m2/kg','type','continuous'));
            obj.addInput('e',struct('name','e','comment','épaisseur d''entrefer','unit','m','type','continuous'));
            obj.addInput('Udc',struct('name','Udc','comment','tension bus continu','unit','V','type','continuous'));
            obj.addInput('Bcr',struct('name','Bcr','comment','induction magnétique culasse rotor','unit','T','type','continuous'));
            obj.addInput('rrs',struct('name','rrs','comment','rapport longueurs rotor / stator','unit','°C','type','continuous'));
            obj.addInput('p',struct('name','p','comment','nombre de paires de pôles','unit','-','type','continuous'));
            obj.addInput('Bd',struct('name','Bd','comment','induction magnétique dents','unit','T','type','continuous'));
            obj.addInput('Be',struct('name','Be','comment','induction magnétique entrefer','unit','T','type','continuous'));
            obj.addInput('Bcs',struct('name','Bcs','comment','induction magnétique culasse stator','unit','T','type','continuous'));
            obj.addInput('Lm',struct('name','Lm','comment','longueur magnétique stator','unit','m','type','continuous'));
            obj.addInput('Ds',struct('name','Ds','comment','diamètre alésage stator','unit','m','type','continuous'));
            obj.addInput('J',struct('name','J','comment','densité de courant','unit','A/m2','type','continuous'));
            
            % declare outputs of model
            obj.addOutput('Mtot',struct('name','Mtot','comment','masse totale','unit','kg'));
            obj.addOutput('Ta',struct('name','Temp','comment','temperature moyenne','unit','°C'));
            obj.addOutput('Dint',struct('name','Dint','comment','diamètre intérieur','unit','m'));
            obj.addOutput('Dext',struct('name','Dext','comment','diamètre extérieur','unit','m'));
            obj.addOutput('Ltot',struct('name','Ltot','comment','longueur totale','unit','m'));
            obj.addOutput('eta',struct('name','rend','comment','rendement','unit','-'));
            obj.addOutput('discr',struct('name','discr','comment','discriminant (doit être >0)','unit','-'));
            obj.addOutput('Pf',struct('name','Pfer','comment','pertes stator','unit','W'));
            obj.addOutput('Pj',struct('name','Pcond','comment','pertes conducteurs','unit','W'));
            obj.addOutput('Imax',struct('name','Idemag','comment','courant de démagnétisation','unit','A'));
            obj.addOutput('t1',struct('name','t1','comment','temps relatif de montée du courant (doit être <1)','unit','-'));
            obj.addOutput('hd',struct('name','hd','comment','hauteur dent','unit','m'));
            obj.addOutput('ha',struct('name','ha','comment','hauteur aimant','unit','m'));
            obj.addOutput('eb',struct('name','hp','comment','hauteur pièce polaire','unit','m'));
            
        end
        
        % procedural part
        
        function output = compute(obj,input)
            
            % convert input structure to variables
            fnames = fieldnames(input);
            for i=1:length(fnames)
                eval([fnames{i} '=input.' fnames{i} ';']);
            end
            
            omega = omega*pi/30;
            omegamax = omegamax*pi/30;
            
            % equations for model
            Ne = 3/2*p;
            muo = 4*pi*1e-7;
            alpha=pi/p;
            beta=pi/p;
            alphai=alpha/5;
            kv=omegamax/omega;
            E=Udc/2/kv;
            I=C*omega/2/E;
            n=4*E/Be/Ds/Lm/omega;
            li=Be/Bd*alphai*Ds/2;
            ld=Be/Bd*alpha*Ds/2;
            eb=(Be/Bd*(alpha/2-asin(ld/Ds))*Ds/2)+Ds/2*(1-cos(asin(ld/Ds)));
            hcs=1/2*Bd/Bcs*ld;
            Ba=Be*alpha*Ds/2/beta/(Ds/2+e)/rrs/kfui;
            hcr=1/2*Ba/Bcr*beta*(Ds/2+e);
            Senc=3/2*n*I/J/kr;
            a=pi;
            b=-(2*pi*(Ds/2-eb)-Ne*(li+ld));
            c=Senc;
            discr=b^2-4*a*c;
            hd=(-b-sqrt(abs(discr)))/2/a;
            Dint=Ds-2*(eb+hd+hcs);
            f=p*omega/2/pi;
            hc=(eb/cos(alpha/2)-Ds/2*(1/cos(alpha/2)-1));
            hi=(Ds/2*(1-cos(alphai/2))+hc*cos(alphai/2));
            Rtb=((ld-li)/4+(Ds/2-eb-hd/2)*pi/2/Ne);
            Lds=(Lm/kfoi+pi*Rtb);
            Ltot=(Lm/kfoi+2*((Ds/2-eb-hd/2)*pi/Ne-li/2));
            Mcs=dt*pi*hcs*(-hcs+2*(Ds/2-eb-hd))*Lm;
            Mds=dt*Ne*((ld+li)*hd+Ds/2*(alphai*(hi+hc)/2+alpha*(eb+hc)/2))*Lm;
            Pf=qt*((f/ft)^1.5)*(Mcs*((Bcs/Bt)^2)+Mds*((Bd/Bt)^2));
            Ta=100; Tcu=110;
            while abs(Ta-Tcu)/Ta>1e-6
                Ta=Tcu;
                ha=-mua*Be*e/(Ba-Br*(1+alphaa*Ta));
                Dext=Ds+2*(e+ha+hcr);
                Rph=rocu*(1+alphacu*Tcu)*n/2*Lds*J/I;
                Pj=2*Rph*I^2;
                Sext=pi/2*(Dext^2)+pi*Dext*Ltot;
                Tcu=Text+(Pj+Pf+Pm)/h/Sext;
            end
            Imax=-4*p/n*(Bc/muo*beta/alpha*(1+2*e/Ds)*rrs*kfui*e+(Bc-Br*(1+alphaa*Ta))*ha/muo/mua);
            Ma=da*p*beta*ha*(ha+2*(Ds/2+e))*Lm*rrs;
            Mcr=dcr*pi*hcr*(hcr+2*(Ds/2+e+ha))*Lm*rrs;
            Mcu=dcu*3/2*n*I/J*Lds;
            Mtot=Ma+Mcr+Mcs+Mds+Mcu;
            eta=(C*omega-Pm)/(C*omega+Pj+Pf);
            FencIrc=muo/16*(n^2)/Ne*hd*Lm/((Ds/2-eb-hd/2)*pi/Ne-(li+ld)/2);
            FentIrc=muo*3/16*(n^2)/Ne*Lm/(e+ha)*alpha*Ds/2;
            FbecIrc=muo*3/16*(n^2)/Ne*hc*Lm/(Ds/2-hc/2)/(pi/Ne-(alpha+alphai)/2);
            Lph=3/2*FentIrc+2*(FencIrc+FbecIrc);
            t1=-Lph/Rph*log(abs(1-3*Rph*I/Udc/(2-1/kv)));
            
            % convert output variables to structure
            fnames = fieldnames(obj.outputs);
            for i=1:length(fnames)
                output.(fnames{i}) = eval(fnames{i});
            end
            
        end
        
    end
end
