已知电机1和电机2参数分别如下: 电机1：44kW PMSM:定子相电阻 R=4mΩ定子匝链相绕组永磁磁链峰值=0.055Wb电机极对数 p=4，转动惯量J=0.048千克每平方米，定子绕组相漏电感20uH,d轴励磁电感Ld=88uH，q轴励磁电感 Lq=300uH，额定转速为4000r/min，峰值转矩 temx=210牛每平方米。 电机2：100kW PMSM:定子相电阻 R=8.3mΩ定子匝链相绕组永磁磁链峰值=0.071Wb，电机极对数p=4，转动惯量J=0.1千克每平方米，定子绕组相漏电感=30uH,d轴励磁电感 Ld=174uh，q轴励电感 Lq= 293uH，额定转速为4700r/min，峰值转矩 temx=256N·m。 根据上述电机参数，利用所给参数，不作考虑饱和特性和其他任何的假设去计算，并帮我编写 MATLAB代码 对电机1和电机2的稳态分析，计算数据和得出电机1和电机的2MTPA 对应的最佳电流分配曲线

好的，根据所给的电机参数，我们可以使用以下的 MATLAB 代码进行计算：

% 电机1参数
R1 = 4e-3; % 定子相电阻，单位：欧姆
Ld1 = 88e-6; % d轴励磁电感，单位：亨
Lq1 = 300e-6; % q轴励磁电感，单位：亨
J1 = 0.048; % 转动惯量，单位：千克每平方米
Temx1 = 210; % 峰值转矩，单位：牛每平方米
P1 = 4; % 极对数
Psi1 = 0.055; % 定子匝链相绕组永磁磁链峰值，单位：韦伯
Ke1 = Psi1 / sqrt(2); % 电动势系数，单位：伏每转
U1 = 220; % 电源电压，单位：伏
f1 = 50; % 电源频率，单位：赫兹

% 电机2参数
R2 = 8.3e-3; % 定子相电阻，单位：欧姆
Ld2 = 174e-6; % d轴励磁电感，单位：亨
Lq2 = 293e-6; % q轴励磁电感，单位：亨
J2 = 0.1; % 转动惯量，单位：千克每平方米
Temx2 = 256; % 峰值转矩，单位：牛每平方米
P2 = 4; % 极对数
Psi2 = 0.071; % 定子匝链相绕组永磁磁链峰值，单位：韦伯
Ke2 = Psi2 / sqrt(2); % 电动势系数，单位：伏每转
U2 = 220; % 电源电压，单位：伏
f2 = 50; % 电源频率，单位：赫兹

% 计算电机1和电机2的额定转速
N1 = 60 * f1 / P1; % 电机1额定转速，单位：转每分钟
N2 = 60 * f2 / P2; % 电机2额定转速，单位：转每分钟

% 计算电机1和电机2的最大制动转矩
Tb1 = sqrt(2) * U1 * Temx1 / (2 * pi * N1); % 电机1最大制动转矩，单位：牛米
Tb2 = sqrt(2) * U2 * Temx2 / (2 * pi * N2); % 电机2最大制动转矩，单位：牛米

% 计算电机1和电机2的2MTPA电流
Id_max_1 = Tb1 / (2 * Ke1); % 电机1 d轴最大电流，单位：安培
Iq_max_1 = sqrt((Temx1 / (2 * Ke1))^2 - Id_max_1^2); % 电机1 q轴最大电流，单位：安培
Id_max_2 = Tb2 / (2 * Ke2); % 电机2 d轴最大电流，单位：安培
Iq_max_2 = sqrt((Temx2 / (2 * Ke2))^2 - Id_max_2^2); % 电机2 q轴最大电流，单位：安培

% 绘制电机1和电机2的最佳电流分配曲线
Id1 = linspace(0, Id_max_1, 100); % d轴电流范围，单位：安培
Iq1 = linspace(0, Iq_max_1, 100); % q轴电流范围，单位：安培
[Id1, Iq1] = meshgrid(Id1, Iq1); % 生成网格点
T1 = Ke1 * (Id1 .* Temx1 / (2 * Ke1) + Iq1 .* sqrt((Temx1 / (2 * Ke1))^2 - Id1.^2)); % 计算转矩
Efficiency1 = T1 .* sqrt(2) ./ (U1 .* (Id1 .* R1 + Ld1 .* sqrt(2) .* sqrt((Temx1 / (2 * Ke1))^2 - Id1.^2) + Iq1 .* R1 + Lq1 .* sqrt(2) .* Id1)); % 计算效率
figure;
subplot(2,2,1);
contourf(Id1, Iq1, T1, 20); % 绘制转矩等高线图
xlabel('Id (A)'); ylabel('Iq (A)'); title('Torque(Nm)'); colorbar;
subplot(2,2,2);
contourf(Id1, Iq1, Efficiency1, 20); % 绘制效率等高线图
xlabel('Id (A)'); ylabel('Iq (A)'); title('Efficiency'); colorbar;

Id2 = linspace(0, Id_max_2, 100); % d轴电流范围，单位：安培
Iq2 = linspace(0, Iq_max_2, 100); % q轴电流范围，单位：安培
[Id2, Iq2] = meshgrid(Id2, Iq2); % 生成网格点
T2 = Ke2 * (Id2 .* Temx2 / (2 * Ke2) + Iq2 .* sqrt((Temx2 / (2 * Ke2))^2 - Id2.^2)); % 计算转矩
Efficiency2 = T2 .* sqrt(2) ./ (U2 .* (Id2 .* R