// To Determine the rotor angle and angular frequency  using runga kutta and euler`s modified method

clc
clear
Pm=3;
r1Pm=1.2;
r2Pm=2;
H=3;
f=60;
Dt=.02;
Pe=1.5;
Do=asind(1.5/3);
do=Do/57.33;
wo=0;
d=0;
K10=0;
l10=62.83*(1.5-1.2*sin(do))*.02;
K20=(377.5574-376.992)*.02;
l20=62.83*(1.5-1.2*sin(do))*.02;
K30=(377.5574-376.992)*.02;
l30=62.83*(1.5-1.2*sin(.5296547))*.02;
K40=l30*0.02;
l40=62.83*(1.5-1.2*sin(.5353094))*.02;
d1=.53528;
Dwo=(3*1.13094+2*1.123045+1.115699)/6;
w1=wo+Dwo;
d1=.53528;
mprintf("Runga-Kutta method-\n")
mprintf("w1=%.6f \nd1=%.5f\n",w1,d1);
d7=1.026;
w7=6.501;
wp=376.992+6.501;
K17=(wp-376.992)*0.02;
l17=62.83*(1.5-1.2*sin(1.026))*.02;
K27=(6.501+.297638)*0.02;
l27=62.83*(1.5-1.2*sin(1.09101))*.02;
K37=(6.501+.2736169)*0.02;
l37=62.83*(1.5-1.2*sin(1.0939863))*.02;
K47=(6.501+.545168)*0.02;
l47=62.83*(1.5-1.2*sin(1.16149))*.02;
Dd7=(K17+2*K27+2*K37+K47)/6;
d8=d7+Dd7;
Dw7=(l17+2*l27+2*l37+l47)/6;
w8=w7+Dw7;
mprintf("d8=%.5f rad.\nw8=%.4frad/sec\n\n",d8,w8)
mprintf("using Euler`s Modified Method-\n");
d0=0;
d10=.524;
w=62.83*(1.5-1.2*sin(.524));
d11=d10+0;
w11=w*.02;
d=1.13094;
dav=(0+d)/2;
wav=(56.547+56.547)/2;
d01=.524+.56547*.02;
w11=0+56.547*0.02;
mprintf("d01=%.4f\nw11=%.5f",d01,w11);