% 編寫建立自激過程微分方程的M—函數
% 將該M—函數定義為sh_ge_se_ex_ode(shunt_generator_self_excited_ode)
function dydt=tongbufadianji(t,y)
% 下面輸入電機基本數據 
r=.26;Rfd=0.13;Rkd=0.0224;Rkq=0.02;ufd=17.8876;w=377;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;ud=100;L0=Ll;
 %--------------------------------------------------------------------------------------------------------------------------------------   
 L=[  -Ld*(cos(w*t)+sin(w*t))-.5*L0*cos(w*t),  -Lq*(cos(w*t)-sin(w*t))+.5*L0*sin(w*t),  Lmd*(cos(w*t)+sin(w*t)),  Lmd*(cos(w*t)+sin(w*t)), Lmq*(cos(w*t)-sin(w*t));   
        -Ld*sin(w*t),                          -Lq*cos(w*t),                            Lmd*sin(w*t),             Lmd*sin(w*t),            Lmq*cos(w*t);
        -Lmd,                                  0,                                       Lfd,                      Lmd,                     0;
        -Lmd,                                  0,                                       Lmd,                      Lkd,                     0;
         0,                                   -Lmq,                                     0,                        0,                       Lkq] ;
G=[  -Ld*(cos(w*t)-sin(w*t))+.5*L0*sin(w*t),  Lq*(cos(w*t)+sin(w*t))+.5*L0*cos(w*t),   Lmd*(cos(w*t)-sin(w*t)),  Lmd*(cos(w*t)-sin(w*t)), -Lmq*(cos(w*t)+sin(w*t));   
     -Ld*cos(w*t),                            Lq*sin(w*t),                             Lmd*cos(w*t),             Lmd*cos(w*t),            -Lmq*sin(w*t);
      0,                                      0,                                       0,                        0,                        0;
      0,                                      0,                                       0,                        0,                        0;
      0,                                      0,                                       0,                        0,                        0] ;
 %--------------------------------------------------------------------------------------------------------------------------------------   
 R=[ -r*(3/2*cos(w*t)+sin(w*t)),         -r*(cos(w*t)-3/2*sin(w*t)),      0 ,           0,          0;
     -r*sin(w*t),                         -r*cos(w*t),                      0 ,           0,          0; 
      0,                                   0,                               Rfd,          0,          0 ;
      0,                                   0,                               0,            Rkd,        0;
      0,                                   0,                               0,            0,          Rkq];
 %--------------------------------------------------------------------------------------------------------------------------------------   
Uabc=[ud,0,ufd,0,0]';
dydt=L\(Uabc-w*G*y-R*y);