
%%%%%%%%%%%%%% This is called from the main program%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%% LOAD FLOW DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


load lfl.dat;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%%%%%%%%%%%%%%%%%%%%%%%INITIAL CONDITIONS FOR GENERATOR %%%%%%%%%%%%%%%%%%%                  
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

xd=gen(:,2)';
xdd=gen(:,3)';
xddd=gen(:,4)';
Td0d=gen(:,5)';
Td0dd=gen(:,6)';
xq=gen(:,7)';
xqd=gen(:,8)';
xqdd=gen(:,9)';
Tq0d=gen(:,10)';
Tq0dd=gen(:,11)';
H=gen(:,12)';
Dm=gen(:,13)';
Vg0=lfl(gen(:,1),2)';
thetag0=pi/180*lfl(gen(:,1),3)';

Vg0bar=(Vg0.*cos(thetag0))+j*Vg0.*sin(thetag0);
Ig0bar=conj((lfl(gen(:,1),4)'+j*lfl(gen(:,1),5)')./Vg0bar);

Eqbar=(Vg0bar+(j*xq).*Ig0bar);
delta0=angle(Eqbar);
Eq=abs(Eqbar);

VQg0=real(Vg0bar);
VDg0=imag(Vg0bar);

Ig0bardq=Ig0bar.*(cos(delta0)-j*sin(delta0));
Vg0bardq=Vg0bar.*(cos(delta0)-j*sin(delta0));

iq0=real(Ig0bardq);
id0=imag(Ig0bardq);

iDg0=imag(Ig0bar);
iQg0=real(Ig0bar);

vqg0=real(Vg0bardq);
vdg0=imag(Vg0bardq);

siq0=-vdg0;
sid0=vqg0;
Efd0=Eq-(xd-xq).*id0;

% d-axis states
sih0=sid0;
sif0=sid0 + (xdd./(xd-xdd)).*Efd0;

% q-axis states
sig0=siq0;
sik0=siq0;

% dummy coil
Edummydd0=-(xqdd-xddd).*iq0;
Tdummy=0.01;

%~~~~~~~~~~~~~~~~~~~~~~field current~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
xfl=(xd.*xdd)./(xd-xdd);
ifd0= (sif0-sid0 )./xfl;
IFD0=ifd0.*xd;

%~~~~~~~~~~~~~~~~ To determine Tdd and Tddd from Td0d and Td0dd~~~~~~~~~%
   a=(1- xd./xdd + xd./xddd);
   b=-(Td0d + Td0dd);
   c=(xddd./xdd).*Td0d.*Td0dd;
   Tddd1= (-b + sqrt(b.*b - 4*a.*c))./(2*a);
   Tddd2= (-b - sqrt(b.*b - 4*a.*c))./(2*a);
   Tddd= min(Tddd1,Tddd2);
   Tdd = Td0d.*Td0dd.*(xddd./xd)./Tddd;

   
%~~~~~~~~~~~~~~~~~~~~~To determine Tqd and Tqdd from Tq0d and Tq0dd~~~~~~~~~%
   a=(1- xq./xqd + xq./xqdd);
   b=-(Tq0d + Tq0dd);
   c=(xqdd./xqd).*Tq0d.*Tq0dd;
   Tqdd1= (-b + sqrt(b.*b - 4*a.*c))./(2*a);
   Tqdd2= (-b - sqrt(b.*b - 4*a.*c))./(2*a);
   Tqdd= min(Tqdd1,Tqdd2);
   Tqd = Tq0d.*Tq0dd.*(xqdd./xq)./Tqdd;

%---------------------------------------------------------------------------%
Eqdd0=((xdd-xddd)./xdd).*sih0 + ((xd-xdd)./xd).*(xddd./xdd).*sif0;
Eddd0=-(((xqd-xqdd)./xqd).*sik0 + ((xq-xqd)./xq).*(xqdd./xqd).*sig0);

Tm0=Eqdd0.*iq0+Eddd0.*id0+id0.*iq0.*(xddd-xqdd);
%----------------------------------------------------------------------------%
clear iq0,iQg0;
clear id0,iDg0;
clear Ig0bardq;
clear Vg0bardq;


[mm order_gen]=sort([setxor(1:nb,gen(:,1)),gen(:,1)']);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Exciter Initial Condition calculations
EFD01=[zeros(1, nb-ngen) Efd0];
EFD0=EFD01(order_gen);
IFD01=[zeros(1, nb-ngen) IFD0];
IFD0_a=IFD01(order_gen);
Vg0bar1=[zeros(1, nb-ngen) Vg0bar];
Vg0bar_a=Vg0bar1(order_gen);
Ig0bar1=[zeros(1, nb-ngen) Ig0bar];
Ig0bar_a=Ig0bar1(order_gen);

%----------------------Single time constant Static type exciters-----~~~--

static_exciter

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


%----------------------- IEEE DC1A -Type%Exciters------------------------

DC1A_exciter

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%----------------IEEE AC4A -Type exciter~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

AC4A_exciter

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


%-----------Main Selector for excitation systems.

%AVR = 1's (DISABLED); =0's (ENABLED)
%AVR=ones(1,nb);
AVR=zeros(1,nb);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SL_static(1:nb)=1;

SL_DC1A(1:nb)=1;

SL_AC4A(1:nb)=1;

%~~~~~~~~~~~~~Indicate the generator number on which a specfic type of ~~%
%~~~~~~~~~~~~~~~~~~~~exciter is to be enabled, otherwise null~~~~~~~~~~~~%

%load ng_static.dat
ng_static=[1,2,3,4];

%load ng_DC1A.dat
ng_DC1A=[];

%load ng_AC4A.dat
ng_AC4A=[];
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SL_static(ng_static)=0;

SL_DC1A(ng_DC1A)=0;

SL_AC4A(ng_AC4A)=0;


%-------------------------Load Modelling----------------------------------

%~~~~~~~~~~~~~~Load varables~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
VL0=lfl(ld(:,1),2)';
theL0=lfl(ld(:,1),3)'*pi/180;
PL0=ld(:,2)';
QL0=ld(:,3)';

%%%%%%%%%%%%%%%%% set p1, p2 and p3, r1, r2 and r3  in load_zip_model.m

load_zip_model

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%------Turbine and Speed Governers- Initial condition calculations.------
H1=[zeros(1, nb-ngen) H];
H_a=H1(order_gen);

Tm01=[zeros(1, nb-ngen) Tm0];
TmA0=Tm01(order_gen);

%speed Governer system with Hydro Turbine.

hydro_turbine

%Speed governer with reheat type turbine

Reheat_turbine

% ---------------------Main Selector for Turbine models---------------%
%TURB = 1's (DISABLED); =0's (ENABLED)
TURB=ones(1,nb);
%TURB=zeros(1,nb);
%TURB(4)=1;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SL_HYDRO(1:nb)=1; 

SL_RHST(1:nb)=1;

%-----Indicate the generator number on which a specfic type of --------%
%------- speed-governor-turbine is to be enabled, otherwise null-------%

%load ng_hydro.dat
ng_hydro=[1,2,3,4];

%load ng_rht.dat
ng_rht=[];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SL_HYDRO(ng_hydro)=0;

SL_RHST(ng_rht)=0;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%---Main Selector for PSS----
%PSS = 1's (DISABLED); =0's (ENABLED)
PSS=ones(1,nb);
%PSS=zeros(1,nb);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SL_slip_pss(1:nb)=1;

SL_delPw_pss(1:nb)=1;

SL_power_pss(1:nb)=1;
%-----Indicate the generator number on which a specfic type of 
%------- PSS is to be enabled, otherwise nullmatrix---

%load ng_slip_pss.dat
ng_slip_pss=[1];

%load ng_delPw_pss.dat
ng_delPw_pss=[];

%load ng_power_pss.dat
ng_power_pss=[];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SL_slip_pss(ng_slip_pss)=0;

SL_delPw_pss(ng_delPw_pss)=0;

SL_power_pss(ng_power_pss)=0;

%---------------PSS models----------------------
%------Slip signal based PSS----m


pss_slip_signal


%----Delta Power-Omega signal based PSS----
pss_delPw_signal


%----Power signal based PSS----
pss_power_signal
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%---------------------------------------------------------------------
busang=lfl(:,3)*pi/180;
Vpre=lfl(:,2).*(cos(busang)+j*sin(busang));
%---------------------------------------------------------------------


