
function res = DistGain_2est(J_1u, J_1d, step1, step21, step22, step23, step24, step25, step26, step_pf, step_pt, step_pb, flag)
global serv HyCase J_opt Col xDs xBs
    % Disturbance in Feed flow rate
    Feed = HyCase.Flowsheet.MaterialStreams.Item('F').MolarFlowValue;
    HyCase.Flowsheet.MaterialStreams.Item('F').MolarFlowValue = Feed * (1 + step21);
    Col.ColumnFlowsheet.Run
    if isempty(J_1u) && isempty(step1) 
        Gd_sub(:,1) = MGain_2est(Feed * step21)';
        J_1d(1) = Goal_2est(xDs,xBs);
    else 
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(1) = (J_1u2d + J_opt - J_1u - J_1d(1)) / (step1 * step21);
    end
    
    HyCase.Flowsheet.MaterialStreams.Item('F').MolarFlowValue = Feed;
    
   

    % Disturbance in Feed Composition
    % +5% in C1OH
    Compfrac = HyCase.Flowsheet.MaterialStreams.Item('F').componentMolarFractionValue;
    compositionDisturbance(1, Compfrac, step22);
  
    
    HyCase.Application.PlayScript([pwd,'\compChange.SCP'])
    if isempty(J_1u) && isempty(step1)
        Gd_sub(:,2) = MGain_2est(step22)';
        J_1d(2) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(2) = (J_1u2d + J_opt - J_1u - J_1d(2)) / (step1 * step22);
    end    

    HyCase.Application.PlayScript([pwd,'\compnormal.SCP'])
    Col.ColumnFlowsheet.Run

    % +5% in C2OH
    
    compositionDisturbance(2, Compfrac, step23);

    HyCase.Application.PlayScript([pwd,'\compChange.SCP'])
    if isempty(J_1u) && isempty(step1)
        Gd_sub(:,3) = MGain_2est(step23)';
        J_1d(3) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(3) = (J_1u2d + J_opt - J_1u - J_1d(3)) / (step1 * step23);
    end    

    HyCase.Application.PlayScript([pwd,'\compnormal.SCP'])
    Col.ColumnFlowsheet.Run

    % -5% in C3OH
    compositionDisturbance(3, Compfrac, step24);

    HyCase.Application.PlayScript([pwd,'\compChange.SCP'])
    if isempty(J_1u) && isempty(step1)
        Gd_sub(:,4) = MGain_2est(step24)';
        J_1d(4) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(4) = (J_1u2d + J_opt - J_1u - J_1d(4)) / (step1 * step24);
    end    

    HyCase.Application.PlayScript([pwd,'\compnormal.SCP'])
    Col.ColumnFlowsheet.Run

    % -5% in C4OH
    compositionDisturbance(4, Compfrac, step25);

    HyCase.Application.PlayScript([pwd,'\compChange.SCP'])
    if isempty(J_1u) && isempty(step1)
        Gd_sub(:,5) = MGain_2est(step25)';
        J_1d(5) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(5) = (J_1u2d + J_opt - J_1u - J_1d(5)) / (step1 * step25);
    end    

    HyCase.Application.PlayScript([pwd,'\compnormal.SCP'])
    Col.ColumnFlowsheet.Run

    
     % Disturbance in Feed quality

    vFrac = HyCase.Flowsheet.MaterialStreams.Item('F').VapourFraction.value;
    HyCase.Flowsheet.MaterialStreams.Item('F').VapourFraction.Erase
    temptoset = fsolve(@reachq,100,[],HyCase,1+step26);
    HyCase.Flowsheet.MaterialStreams.Item('F').Temperature.Value = temptoset;
    Col.ColumnFlowsheet.Run
    
    if isempty(J_1u) && isempty(step1) 
        Gd_sub(:,6) = MGain_2est(step26)';
        J_1d(6) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(6) = (J_1u2d + J_opt - J_1u - J_1d(6)) / (step1 * step26);
    end
        
    HyCase.Flowsheet.MaterialStreams.Item('F').Temperature.Erase
    HyCase.Flowsheet.MaterialStreams.Item('F').VapourFraction.value = vFrac;
    Col.ColumnFlowsheet.Run
    
    
     % Disturbance in Feed pressure

    Pnom_Feed = HyCase.Flowsheet.MaterialStreams.Item('F').Pressure.value;
    HyCase.Flowsheet.MaterialStreams.Item('F').Pressure.value = Pnom_Feed*(1+step_pf);
    Col.ColumnFlowsheet.Run
    if isempty(J_1u) && isempty(step1) 
        Gd_sub(:,7) = MGain_2est(Pnom_Feed*step_pf)';
        J_1d(7) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(7) = (J_1u2d + J_opt - J_1u - J_1d(7)) / (step1 * step_pf);
    end
    HyCase.Flowsheet.MaterialStreams.Item('F').Pressure.value = Pnom_Feed;
    

    % Disturbance in top pressure

    Pnom_top = Col.ColumnFlowsheet.Pressures.value(1);
    PressChange('PconChange.scp',serv,Col,Pnom_top*(1+step_pt));
    Col.ColumnFlowsheet.Run
    
    if isempty(J_1u) && isempty(step1) 
        Gd_sub(:,8) = MGain_2est(step_pt*Pnom_top)';
        J_1d(8) = Goal_2est(xDs,xBs);
    else
        J_1u2d = Goal_2est(xDs,xBs);
        Jud_sub(8) = (J_1u2d + J_opt - J_1u - J_1d(8)) / (step1 * step_pt);
    end
    PressChange('PconChange.scp',serv,Col,Pnom_top);
    
   
     % Disturbance in top pressure

%     Pnom_reb = Col.ColumnFlowsheet.Pressures.value(end);
%     PressChange('PrebChange.scp',serv,Col,Pnom_reb*(1+step_pb));
%     Col.ColumnFlowsheet.Run
%     
%     if isempty(J_1u) && isempty(step1) 
%         Gd_sub(:,9) = MGain_2est(Pnom_reb*step_pb)';
%         J_1d(9) = Goal_2est(xDs,xBs);
%     else
%         J_1u2d = Goal_2est(xDs,xBs);
%         Jud_sub(9) = (J_1u2d + J_opt - J_1u - J_1d(9)) / (step1 * step_pb);
%     end
%     PressChange('PrebChange.scp',serv,Col,Pnom_reb);
%    
    
    
    
   if isempty(J_1u) && isempty(step1) 
       if flag == 1
           res = J_1d;
       else
          res = Gd_sub; 
       end
    else
        res = Jud_sub;
    end  

    
