% ***************************************************************************
% COPYRIGHT BY DAVE JEWELL AND T.E. MARLIN 1994
% FILNAME : INITIAL.M
% Version 2.0
% LAST UPDATED : DECEMBER 1994
% DISCRIPTION OF FILE :
% THIS M-FILE IS ACCESSED FROM THE MAIN SIMULATION AND CONTAINS THE GLOBAL
% VARIABLES FOR PLANT CALCULATIONS, MATRIX STORAGE AND PLANT CONSTANTS
% ***************************************************************************

% ***************************************************************************
% *************************** NOISE SETTINGS ********************************
% ***************************************************************************
% NOISE IS ADDED TO THE SENSORS (PLANT)
if NOISEYN == 1
		    ; % sets the random numbers to a normal distribution
     randn('seed',0); % resets the seed to 0 for identical repeat
end

% *********************  SIMULATION PARAMETERS  *****************************
SStime=max([max(max(DT)) max(DDT) max(max(DT_m))]); % INTIAL STEADY-STATE TIME FOR SIMULATION
nstart=round(SStime/delt)+3;                    % START OF SIMULATION
nmax=round(tend/delt)+nstart;                   % END OF SIMULATION
if nmax > maxstep
    nmax=maxstep;                               % SETS NUMBER OF ITERATIONS
    delt=(SStime+tend)/(maxstep-3);             % CALCULATES STEP SIZE
    nstart=round(SStime/delt)+3;                % START OF SIMULATION
end                                             % WHILE
gamma=round(DT./delt);
Dgamma=round(DDT./delt);
%gamma_m=round(DT_m./delt);
tSTEP=0.1*tend;
nSTEP=round(tSTEP/delt)+nstart;
% ****************  MATRIX CREATION FOR STORAGE  ****************************
cvp=zeros(2,2);
cvp1=zeros(2,2);
cvp2=zeros(2,2);
noise1=zeros(nmax,1);
noise2=zeros(nmax,1);
cv1=zeros(nmax,1);
cv2=zeros(nmax,1);
cvd=zeros(2,1);
cvd1=zeros(2,1);
cvd2=zeros(2,1);
dist=zeros(nmax,1);
dist(nSTEP:nmax,1)=ones(nmax-nSTEP+1,1).*MAGSTEP;

% ***************************************************************************
% STEADY-STATE INITIALIZATION FOR PV, MV and SP (data comes from base file)
% ***************************************************************************
% ADDS NOISE TO THE INITIAL VARIABLES IF NOISE IS ON THE COMPLETE SIMULATION
if NOISEYN==1 
cv1(1:nstart,1)=ones(nstart,1)*initcv(1)+Kn(1)*(rand(nstart,1));
cv2(1:nstart,1)=ones(nstart,1)*initcv(2)+Kn(2)*(rand(nstart,1));

else
cv1(1:nstart,1)=ones(nstart,1)*initcv(1);
cv2(1:nstart,1)=ones(nstart,1)*initcv(2);
end

% ***************************************************************************
% ********  EVALUATION OF CONSTANTS FOR THE DIFFERENCE EQUATIONS  ***********
% *****************  SECOND-ORDER LEAD-LAG PLUS DEADTIME  *******************
% *****************************  PROCESS  ***********************************
% ***************************************************************************
	tau12=tau1-tau2;
	tauA1=tauA-tau1;
	tau2A=tau2-tauA;
	exp1=zeros(2,2);
	exp2=zeros(2,2);
	a1=zeros(2,2);
	a2=zeros(2,2);
	b1=zeros(2,2);
	b2=zeros(2,2);

for x=1:2
	for y=1:2
	if Kp(x,y) ~= 0         
	if (tau12(x,y)==0) & (tau1(x,y)==0)
	    tau1(x,y)=1e-3;     
	elseif (tau12(x,y)==0) & (tau1(x,y)~=0)
	    tau1(x,y)=1.01*tau1(x,y);
	end

	tau12(x,y)=tau1(x,y)-tau2(x,y);
	tauA1(x,y)=tauA(x,y)-tau1(x,y);
	tau2A(x,y)=tau2(x,y)-tauA(x,y);

	if tau1(x,y)>0
		exp1(x,y)=exp(-delt/tau1(x,y));
	end                                     % IF
	
	if tau2(x,y)>0
	exp2(x,y)=exp(-delt/tau2(x,y));
	end                                     % IF
	    a1(x,y)=-exp1(x,y)-exp2(x,y);
	    a2(x,y)=exp1(x,y)*exp2(x,y);
	    b1(x,y)=Kp(x,y)*(1+tauA1(x,y)/tau12(x,y)*exp1(x,y)+...
	    tau2A(x,y)/tau12(x,y)*exp2(x,y));
	    b2(x,y)=Kp(x,y)*(a2(x,y)+tauA1(x,y)/tau12(x,y)*...
	    exp2(x,y)+tau2A(x,y)/tau12(x,y)*exp1(x,y));
	end     % Kp CHECK
  end           % y LOOP
end             % x LOOP

% ***************************************************************************
% ********  EVALUATION OF CONSTANTS FOR THE DIFFERENCE EQUATIONS  ***********
% *****************  SECOND-ORDER LEAD-LAG PLUS DEADTIME  *******************
% *****************************  DISTURBANCE  *******************************
% ***************************************************************************
Dtau12=Dtau1-Dtau2;
DtauA1=DtauA-Dtau1;
Dtau2A=Dtau2-DtauA;
Dexp1=zeros(2,1);
Dexp2=zeros(2,1);
Db1=zeros(2,1);
Db2=zeros(2,1);
Da1=zeros(2,1);
Da2=zeros(2,1);

for l=1:2
	if Kd(l) ~= 0
	if (Dtau12(l)==0) & (Dtau1(l)==0)
	    Dtau1(l)=1e-3;      
	elseif (Dtau12(l)==0) & (Dtau1(l)~=0)
	    Dtau1(l)=1.01*Dtau1(l);
	end

	Dtau12(l)=Dtau1(l)-Dtau2(l);
	DtauA1(l)=DtauA(l)-Dtau1(l);
	Dtau2A(l)=Dtau2(l)-DtauA(l);

	if Dtau1(l)>0
		Dexp1(l)=exp(-delt/Dtau1(l));
	end             % if
	
	if Dtau2(l)>0
		Dexp2(l)=exp(-delt/Dtau2(l));
	end             % if
	    Da1(l)=-Dexp1(l)-Dexp2(l);
	    Da2(l)=Dexp1(l)*Dexp2(l);
	    Db1(l)=Kd(l)*(1+DtauA1(l)/Dtau12(l)*Dexp1(l)+...
	    Dtau2A(l)/Dtau12(l)*Dexp2(l));
	    Db2(l)=Kd(l)*(Da2(l)+DtauA1(l)/Dtau12(l)*...
	    Dexp2(l)+Dtau2A(l)/Dtau12(l)*Dexp1(l));
	end     % GAIN CHECK
end     % l

% ***************************************************************************
% ****************  MATRIX CREATION FOR STORAGE  ****************************
% ***************************************************************************
contout1=zeros(nmax,1);
contout2=zeros(nmax,1);
deltacontout1=zeros(nmax,1);
deltacontout2=zeros(nmax,1);
initsp=zeros(2,1);
sp1=zeros(nmax,1);
sp2=zeros(nmax,1);
mv1=zeros(nmax,1);
mv2=zeros(nmax,1);
error1=zeros(nmax,1);
error2=zeros(nmax,1);

% ***************************************************************************
% ****************** INITIALIZATION OF INSTRUMENTS **************************
% ***************************************************************************
initsp=initcv;
mv1(1:nmax,1)=ones(nmax,1)*initmv(1);
mv2(1:nmax,1)=ones(nmax,1)*initmv(2);
sp1(1:nstart-1,1)=ones(nstart-1,1)*initsp(1);
sp2(1:nstart-1,1)=ones(nstart-1,1)*initsp(2);

dmc_init;