function transF_2prodConv_Estim_xAsU(filename,savename,DOF1, DOF2)
t1 = 65.87; % deg C, boiling temperature of Methanol 
t2 = 119.4; % deg C, boiling temperature of Butanol 
% General Transfer Function
global serv Col HyCase D B nominal Specs J_opt xDs xBs y_nom

serv = actxserver('Hysys.Application.2006.5');
HyCase = serv.SimulationCases.Open([pwd,filename]);
HyCase.Activate;
% TWOPRODUCTCOLUMN_PREF.HSC
Col = HyCase.Flowsheet.Operations.Item('T-100');
y_nom = Col.ColumnFlowsheet.TemperaturesValue;
Specs = Col.ColumnFlowsheet.Specifications;

D     = Col.Flowsheet.MaterialStreams.Item('D');
B     = Col.Flowsheet.MaterialStreams.Item('B');

fid = fopen('Fconv_all_final_pref.txt');
rl = fgetl(fid);
i = 1;
Fmat_2est = [];
while 1
    rl = fgetl(fid);
    if ~ischar(rl), break, end
    Fmat_2est(i,:) = str2num(rl);
    i = i+1;
end
fclose(fid);
% this piece of code is to get the F matrix from the data in txt file
save Fmat_2est.mat Fmat_2est

nominal = Fmat_2est(1,:);

for i = 1:Specs.count
    Specs.Item(i-1).goalValue = Specs.Item(i-1).CurrentValue;
end

bests(1) = Specs.Item(DOF1).CurrentValue;
bests(2) = Specs.Item(DOF2).CurrentValue;
bests(3) = HyCase.Flowsheet.MaterialStreams.Item('F').MolarFlowValue;
load_bests_2est_xAsU
Col.ColumnFlowsheet.Run

xDs = D.ComponentMolarFractionValue;
xBs = B.ComponentMolarFractionValue;

CV1s = Specs.Item(DOF1).CurrentValue;
CV2s = Specs.Item(DOF2).CurrentValue;

J_opt = nominal(22);


%% Deriving F matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% stp_F stp_q stp_za stp_zb stp_zc stp_zd stp_pf stp_pt stp_pb 
obj = Fmat_2est(:,22);


Findex_m = [2 4 6 8 10 12 14 16];
Fmat_m = Fmat_2est(Findex_m,:);

Fi_m = (Fmat_m(1:end,23:end) - repmat(nominal(23:end), [size(Fmat_m,1),1])) ./ ...
    ((diag(Fmat_m) - nominal(1:8)') * ones(1,size(Fmat_m,2)-22));
F_m = Fi_m';

ally = Fmat_m(1:end,23:end);
Y0 = (nominal(23:end) - t1) ./ (t2 - nominal(23:end));
Y = (ally - t1) ./ (t2-ally);
Fiscal_m = (Y - repmat(Y0, [size(Fmat_m,1),1])) ./ ...
    ((diag(Fmat_m) - nominal(1:8)') * ones(1,size(Fmat_m,2)-22));
Fscal_m = Fiscal_m';

%plus direction for disturbances
Findex_p = [3 5 7 9 11 13 15 17];
Fmat_p = Fmat_2est(Findex_p,:);

Fi_p = (Fmat_p(1:end,23:end) - repmat(nominal(23:end), [size(Fmat_p,1),1])) ./ ...
    ((diag(Fmat_p) - nominal(1:8)') * ones(1,size(Fmat_p,2)-22));
F_p = Fi_p';

ally = Fmat_p(1:end,23:end);
Y0 = (nominal(23:end) - t1) ./ (t2 - nominal(23:end));
Y = (ally - t1) ./ (t2-ally);
Fiscal_p = (Y - repmat(Y0, [size(Fmat_p,1),1])) ./ ...
    ((diag(Fmat_p) - nominal(1:8)') * ones(1,size(Fmat_p,2)-22));
Fscal_p = Fiscal_p';
% This part is written to check if the step size for perturbation was ok
% and if we are in the linear region

%% Plus step in u%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%u(1)    =   yD

stp_x1 = .08;  
stp_x2 = .08;

disp('u1 - plus')
load_bests_2est_xAsU
Col.ColumnFlowsheet.Run
x1 = Specs.Item(DOF1).CurrentValue;
X1_scal = -log(x1/(1-x1));
Specs.Item(DOF1).GoalValue = x1 * (1 + stp_x1);
X1_scalo = -log((x1/(1-x1)) * (1 + stp_x1));
Col.ColumnFlowsheet.Run
G_x1_p = MGain_2est(x1 * stp_x1);
Gscal_x1_p = MGain_2est(X1_scalo-X1_scal);
Gscalxy_x1_p = MGain_2est_xyScal(X1_scalo-X1_scal,t1,t2);
G1_x1_p = Gain_G1_2est(x1 * stp_x1,DOF1,DOF2,CV1s,CV2s);
G1scal_x1_p = Gain_G1scal_2est(X1_scalo-X1_scal,3,2); %not correct


% u(2)        : xB
disp('u2 - plus')
load_bests_2est_xAsU
x2 = Specs.Item(DOF2).CurrentValue;
X2_scal = log(x2/(1-x2));
Specs.Item(DOF2).GoalValue = x2 * (1 + stp_x2);
X2_scalo = log((x2/(1-x2)) * (1 + stp_x2));
Col.ColumnFlowsheet.Run
G_x2_p = MGain_2est(x2 * stp_x2);
Gscal_x2_p = MGain_2est(X2_scalo-X2_scal);
Gscalxy_x2_p = MGain_2est_xyScal(X2_scalo-X2_scal,t1,t2);
G1_x2_p = Gain_G1_2est(x2 * stp_x2,DOF1,DOF2,CV1s,CV2s);
G1scal_x2_p = Gain_G1scal_2est(X2_scalo-X2_scal,3,2);

G_p = [G_x1_p; G_x2_p];
G1_p = [G1_x1_p; G1_x2_p];
G_p = G_p';

Gscal_p = [Gscal_x1_p; Gscal_x2_p];
Gscalxy_p = [Gscalxy_x1_p; Gscalxy_x2_p]';
G1scal_p = [G1scal_x1_p; G1scal_x2_p];
Gscal_p = Gscal_p';

%% Minus step in u%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('u1 - minus')
load_bests_2est_xAsU
Col.ColumnFlowsheet.Run
x1 = Specs.Item(DOF1).CurrentValue;
X1_scal = -log(x1/(1-x1));
Specs.Item(DOF1).GoalValue = x1 * (1 - stp_x1);
X1_scalo = -log((x1/(1-x1)) * (1 - stp_x1));
Col.ColumnFlowsheet.Run
G_x1_m = MGain_2est(x1 * stp_x1);
Gscal_x1_m = MGain_2est(X1_scalo-X1_scal);
Gscalxy_x1_m = MGain_2est_xyScal(X1_scalo-X1_scal,t1,t2);
G1_x1_m = Gain_G1_2est(-x1 * stp_x1,DOF1,DOF2,CV1s,CV2s);
G1scal_x1_m = Gain_G1scal_2est(X1_scalo-X1_scal,3,2);

disp('u2 - minus')
load_bests_2est_xAsU
x2 = Specs.Item(DOF2).CurrentValue;
X2_scal = log(x2/(1-x2));
Specs.Item(DOF2).GoalValue = x2 * (1 - stp_x2);
X2_scalo = log((x2/(1-x2)) * (1 - stp_x2));
Col.ColumnFlowsheet.Run
G_x2_m = MGain_2est(x2 * stp_x2);
Gscal_x2_m = MGain_2est(X2_scalo-X2_scal);
Gscalxy_x2_m = MGain_2est_xyScal(X2_scalo-X2_scal,t1,t2);
G1_x2_m = Gain_G1_2est(-x2 * stp_x2,DOF1,DOF2,CV1s,CV2s);
G1scal_x2_m = Gain_G1scal_2est(X2_scalo-X2_scal,3,2);

load_bests_2est_xAsU

G_m = [G_x1_m; G_x2_m];
G1_m = [G1_x1_m; G1_x2_m];

G_m = G_m';

Gscal_m = [Gscal_x1_m; Gscal_x2_m];
G1scal_m = [G1scal_x1_m; G1scal_x2_m];
Gscal_m = Gscal_m';
Gscalxy_m = [Gscalxy_x1_m; Gscalxy_x2_m]';

checkG = (abs(G_m) - abs(G_p))./abs(G_m);
checkG1 = (abs(G1_m) - abs(G1_p))/abs(G1_m);

Gy = G_p;
G1 = G1_p;

Gyscal = Gscal_p;
G1scal = G1scal_p;

%% Saving the data

Gy_orig = Gy;
Gyscal_orig = Gyscal;
Gscalxy_orig = Gscalxy_p;
Gy = Gy_orig(2:end-1,:);
Gyscal = Gyscal_orig(2:end-1,:);
F = F_p(2:end-1,:);
Fscal = Fscal_p(2:end-1,:);

Gyscalxy = Gscalxy_orig(2:end-1,:);
%d       = [stp_F;stp_za;stp_zb;stp_zc;stp_zd; stp_q;stp_pf; stp_pt];
d       = [10;.05;.05;.05;.05; .1;10; 10];

% Scaling matrices Wd and Wn
Wd      = diag(d);

% implementation error
e_n     = 0.1; % degrees Kelvin
Wn      = diag(ones(length(Gy),1)*e_n);

save(savename,'Gy','Gyscal','Gyscalxy','G1','G1scal','F','Fscal','Wd','Wn','y_nom','CV1s','CV2s')
