%Motor2.m  -A Ruina,     Feb 16, 2010
%This is NOT adequately commented
disp('************************************')
disp(['START Calculation at   '  datestr(now)])
clear all  %added on 2/21/10

syms omega V I  T k c R x D 


%The basic motor equations
eq1 = T- k*I + c*omega;
eq2 = V- k*omega - I*R;

%Find I,V from T and omega
S=solve(eq1,eq2, I, V);

%Write power in and out in terms of T and omega
Pin = S.I*S.V;  %S.I and S.V are in terms of T and omega
Pout = T*omega;

eff1 = simplify(Pout/Pin);


eff2 = subs(eff1, T, x*omega);  % that is, x is defined to be T/omega

eff3=simplify(eff2)


%Use first derivative to find the local max
bestx = solve(diff(eff3, x))  %eff replaced by eff3 on 2/21/10
bestx = bestx(1) %  just keep the positive root, added on 2/21/10

% Substitute back in to get best efficiency
besteff1 = simple(subs(eff3, x, bestx));

% Try to tidy up and simplify
besteff2=simple(subs(besteff1, k, D*sqrt(R*c))); % defines D = k/sqrt(R*c)
pretty(besteff2)


%Specifications for Faulhaver 2657 - 012 CR
    k= 0.0173; %in consistant SI units
    R=0.71; %ohms. actually, the spec is wrong and R = 1.3
    c=3e-6; %in consistant SI units
    D = k/sqrt(R*c)   %The key motor parameter

besteffnumerical = subs(besteff2)
bestxnumerical = subs(bestx) % Remember, this is the ratio T/omega

%Motor use specifications, given by design team, say
Pout = 30; % in watts  
omega_out = 5;  % in rad/s

% We know the bestx = T/omega. We know Pout = T*omega. So
omegamotor = sqrt(Pout/bestxnumerical)
G = omegamotor/omega_out
Tmotornum = bestxnumerical * omegamotor
Vbat = subs(S.V); Vbat = subs(Vbat,  [T omega],[Tmotornum omegamotor])

disp(['END Calculation at   '  datestr(now)])