% 編寫建立自激過程微分方程的M—函數
% 將該M—函數定義為sh_ge_se_ex_ode(shunt_generator_self_excited_ode)
function dydt=tongbufadianji(t,y)
% 首先指定全局變量
%global   a1 a2 a3 a4 Rf  
% 下面輸入電機基本數據：3.0892E-04  3.2164E-03  9.7153E-04 
r=.26;Rfd=0.13;Rkd=0.0224;Rkq=0.02;ufd=17.8876;w=314;Ll=1.14e-3;Lmd=13.7e-3;Lmq=11.0e-3;
Llfd=2.1e-3;Llkd=1.4e-3;Llkq=1e-3;
Lfd=Llfd+Lmd;Lkd=Llkd+Lmd;Lkq=Llkq+Lmq;Ld=Lmd+Ll;Lq=Lmq+Ll;L0=Ll;
%4510116748510111/6/74951011"6"74:51011"6$;<5"017;=510/17;85"01(7;9510$67;:510'17;+<510%67;+:510)6
 %--------------------------------------------------------------------------------------------------------------------------------------   

%L=[  Ld,               0,                 Lmd,	        Lmd,         0,            0+Ld*sin(w*t)+0; 
%     0,                Lq,                0,	        0,           Lmq,          0+0*sin(w*t)+Lq*cos(w*t);
 %    Lmd,              0,                 Lfd,          Lmd,         0,            0+Lmd*sin(w*t)+0;
 %    Lmd,              0,                 Lmd,          Lkd,         0,            0+Lmd*sin(w*t)+0;
 %    0,                Lmq,               0,            0,           Lkq,          0+0*sin(w*t)+Lmq*cos(w*t);
 %   (Ld-L0)*sin(w*t),  (Lq-L0)*cos(w*t),  Lmd*sin(w*t), Lmd*sin(w*t),Lmq*cos(w*t), L0+(Ld-L0)*sin(w*t)*sin(w*t)+(Lq-L0)*cos(w*t)*cos(w*t) ] ;
 L=[  
        Lfd,                    -2*Lmd*cos(w*t);
        Lmd*cos(w*t),           -(Ld+Lq+L0)+(Ld-Lq)/2*cos(2*w*t)] ;

 
 
 
 %--------------------------------------------------------------------------------------------------------------------------------------   
    
G=[ 
         0,        2*Lmd*sin(w*t);
     -Lmd*sin(w*t),  2*(Ld-Lq)*sin(2*w*t)];

%G=[  0,     -Lq,    0,       0,	       Lmq,0+0-Lq*cos(w*t)+Ld*cos(w*t)+0*sin(w*t); 
%     Ld,     0,     Mafd0,   Makd0,	   0,  0+Ld*sin(w*t)+0+0*cos(w*t)+Lq*sin(w*t);
%     0,      0,     0,       0,        0,  0+0+0+Lmd*cos(w*t)+0*sin(w*t);
%     0,      0,     0,       0,        0,  0+0+0+Lmd*cos(w*t)+0*sin(w*t);
%     0,      0,     0,       0,        0,  0+0+0+0*cos(w*t)+Lmq*sin(w*t);
%    (Ld-L0)*cos(w*t),  -(Lq-L0)*sin(w*t),  Lmd*cos(w*t)+Ld*cos(w*t),
%    Lmd*cos(w*t),-Lmq*sin(w*t), 0+(Ld-L0)*cos(w*t)*sin(w*t)-(Lq-L0)*sin(w*t)*cos(w*t)+(Ld-L0)*sin(w*t)*cos(w*t)+(Lq-L0)*cos(w*t)*sin(w*t)];


 %--------------------------------------------------------------------------------------------------------------------------------------   
    
% 下面進行曲線擬合：

%R=[ r,   0,   0,    0,   0,  0+r*sin(w*t)+0;
 %   0,   r,   0,    0,   0,  0+0+r*cos(w*t);
%    0,   0,   Rfd,  0,   0,  0+0+0;
 %   0,   0,   0,    Rkd, 0,  0+0+0;
 %   0,   0,   0,    0,   Rkq,0+0+0;
 %   0,   0,   0,    0,   0,  r+0+0];
 
 R=[ 
       Rfd,    0;
      0,      -3*r];

 
 %--------------------------------------------------------------------------------------------------------------------------------------   
 
Uabc=[ufd,0]';

dydt=L\(Uabc-w*G*y-R*y);





% 下面列寫標準形式微分方程并將感應電勢用勵磁電流的擬合多項式函數表示：
% 這里的a(1)、a(2)、a(3)、a(4)是指上一步的多項式擬合系數a(i)
%dydt=[1/Lf*(a(1)*iff^3+a(2)*iff^2+a(3)*iff^1+a(4))-Rf/Lf*iff];
% 給全局變量賦值
%a1=a(1);a2=a(2);a3=a(3);a4=a(4);