Vh=150;

%olc=[ P1  P2   Ih      n    ]
olc=[346 44 3.76 -528;
     336 49 3.61 -279;
     312 47 3.33 0;
     304 57 3.20 205;
     282 59 2.95 437;
     250 59 2.60 685;
     209 53 2.18 890;
     182 48 1.90 1015;
     126 33 1.31 1205;
     97 22 1.01 1299;
     56 4 0.65 1412;
     42 -2 0.55 1450;
     26 -11 0.48 1500;
     10 -21 0.55 1549;
     -7 -34 0.73 1601;
     -38 -63 1.26 1705;
     -60 -95 1.90 1803;
     -64 -126 2.57 1910;
     -46 -152 3.20 2034;
     -8 -170 3.80 2166
    ];

% HESAPLAMALAR
P1=olc(:,1); P2=olc(:,2); P=P1+P2; Q=sqrt(3)*(P1-P2);
I1=olc(:,3);    % Y bağlı
V1=Vh/sqrt(3);    % Y bağlı
n=olc(:,4);
[I1min,yuksuz]=min(I1); % Senkron hızın indisi = yuksuz
PFe=P(yuksuz) ; % Demir kaybi, voltaj sabit olduğu için hep aynı
                % Bakır kaybı ise yüksüz (senkron hızda) çalışmada
                % ihmal edilebilecek kadar küçük.

kilitli=3 % 3. değer takımı kilitli rotor durumu olduğu için indisi = kilitli
rT=(P(kilitli)-PFe)/3/I1(kilitli)^2; % Kilitliyken de demir kaybı var ama
                                     % Pm yok. PFe çıkınca PCu kalıyor.
                                     % Buradan rT=r1+r2' hesaplanıyor.

PCu=3*(rT)*I1.^2;     % Tüm ölçüm takımları için bakır kaybı akımın karesiyle orantılı
Pm=(P-PCu-PFe);   % Elektromekanik güç
w=2*pi*n/60;
Tm=Pm./w;   % Böyle hesaplanan elektromekanik tork, w=0 için hata veya uyarı verir
            % w=0 için tork formülünden s=1 için hesaplanacak

ws=2*pi*50/2; % Senkron hız rad/s olarak (4 kutuplu olduğu için 2'ye bölündü)
r2ussu=rT-8; % r1=8 ohm ölçülmüştü. s=1 için tork formülünde r2ussu gerekiyor.
Tm(kilitli)=3*I1(kilitli)^2*r2ussu / ws % s=1 için tork formülü 
                                        % pay'daki 1-s paydadaki (1-s)ws'inkiyle sadeleşir

figure(1),set(1,'Name','Tm  hıza(w) karşı'),plot(w,Tm,w,0*Tm)
figure(2),set(2,'Name','Q hıza karşı'),plot(n,Q)
fi1=atan2(Q,P)*180/pi; % atan2 fonksiyonu 4 çeyrek bölge için de açıyı doğru hesaplar;
                      % ama önce sanal sonra gerçel kısım yazılır ve sonuç
                      % radyan çıkar. Biz bunu dereceye çevirdik.
figure(3),set(3,'Name','fi1 hıza karşı (I1 ''in V1 ''den GERİ olma açısı)'),plot(n,fi1)