text:
function M = stiefelgeneralizedfactory(n, p, B)
% Returns a manifold structure of "scaled" orthonormal matrices.
%
% function M = stiefelgeneralizedfactory(n, p)
% function M = stiefelgeneralizedfactory(n, p, B)
%
% The generalized Stiefel manifold is the set of "scaled" orthonormal
% nxp matrices X such that X'*B*X is identity. B must be positive definite.
% If B is identity, then this is the standard Stiefel manifold.
%
% The generalized Stiefel manifold is endowed with a scaled metric
% by making it a Riemannian submanifold of the Euclidean space,
% again endowed with the scaled inner product.
%
% Some notions (not all) are from Section 4.5 of the paper
% "The geometry of algorithms with orthogonality constraints",
% A. Edelman, T. A. Arias, S. T. Smith, SIMAX, 1998.
%
% Paper link: http://arxiv.org/abs/physics/9806030.
%
% Note: egrad2rgrad and ehess2rhess involve solving linear systems in B. If
% this is a bottleneck for a specific application, then a way forward is to
% create a modified version of this file which preprocesses B to speed this
% up (typically, by computing a Cholesky factorization of it, then calling
% an appropriate solver).
%
% See also: stiefelfactory  grassmannfactory  grassmanngeneralizedfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, June 30, 2015.
% Contributors:
%
% Change log:
%

    if ~exist('B', 'var') || isempty(B)
        B = speye(n); % Standard Stiefel manifold.
    end

    M.name = @() sprintf('Generalized Stiefel manifold St(%d, %d)', n, p);

    M.dim = @() (n*p - .5*p*(p+1));

    M.inner = @(X, eta, zeta) trace(eta'*(B*zeta)); % Scaled metric.

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(X, Y) error('stiefelgeneralizedfactory.dist not implemented yet.');

    M.typicaldist = @() sqrt(p);

    % Orthogonal projection of an ambient vector U to the tangent space
    % at X.
    M.proj = @projection;
    function Up = projection(X, U)
        BX = B*X;

        % Projection onto the tangent space
        Up = U - X*symm(BX'*U);
    end

    M.tangent = M.proj;

    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)

        % First, scale egrad according the to the scaled metric in the
        % Euclidean space.
        egrad_scaled = B\egrad;

        % Second, project onto the tangent space.
        % rgrad = egrad_scaled - X*symm((B*X)'*egrad_scaled);
        %
        % Verify that symm(BX'*egrad_scaled) = symm(X'*egrad).

        rgrad = egrad_scaled - X*symm(X'*egrad);
    end

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)
        egraddot = ehess;
        Xdot = H;

        % Directional derivative of the Riemannian gradient.
        egrad_scaleddot = B\egraddot;
        rgraddot = egrad_scaleddot - Xdot*symm(X'*egrad)...
            - X*symm(Xdot'*egrad)...
            - X*symm(X'*egraddot);

        % Project onto the tangent space.
        rhess = M.proj(X, rgraddot);
    end

    M.retr = @retraction;
    function Y = retraction(X, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = guf(X + t*U); % Ensure that Y'*B*Y is identity.
    end

    M.exp = @exponential;
    function Y = exponential(X, Z, t)
        if nargin < 3
            t = 1.0;
        end
        Y = retraction(X, Z, t);
        warning('manopt:stiefelgeneralizedfactory:exp', ...
               ['Exponential for generalized Stiefel manifold ' ...
                'manifold not implemented yet. Used retraction instead.']);
    end

    M.hash = @(X) ['z' hashmd5(X(:))];

    M.rand = @random;
    function X = random()
        X = guf(randn(n, p)); % Ensure that X'*B*X is identity;
    end

    M.randvec = @randomvec;
    function U = randomvec(X)
        U = projection(X, randn(n, p));
        U = U / norm(U(:));
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(X) zeros(n, p);

    % This transport is compatible with the generalized polar retraction.
    M.transp = @(X1, X2, d) projection(X2, d);

    M.vec = @(X, u_mat) u_mat(:);
    M.mat = @(X, u_vec) reshape(u_vec, [n, p]);
    M.vecmatareisometries = @() false;

    % Some auxiliary functions
    symm = @(D) (D + D')/2;

    function X = guf(Y)
        % Generalized polar decomposition of an n-by-p matrix Y.
        % X'*B*X is identity.

        % Method 1
        [u, ~, v] = svd(Y, 0);

        % Instead of the following three steps, an equivalent, but an
        % expensive way is to do X = u*(sqrtm(u'*(B*u))\(v')).
        [q, ssquare] = eig(u'*(B*u));
        qsinv = q/sparse(diag(sqrt(diag(ssquare))));
        X = u*((qsinv*q')*v'); % X'*B*X is identity.

        % Another computation using restricted_svd
        % [u, ~, v] = restricted_svd(Y);
        % X = u*v'; % X'*B*X is identity.

    end

    function [u, s, v] = restricted_svd(Y)
        % We compute a thin svd-like decomposition of an n-by-p matrix Y
        % into matrices u, s, and v such that u is an n-by-p matrix
        % with u'*B*u being identity, s is a p-by-p diagonal matrix
        % with positive entries, and v is a p-by-p orthogonal matrix.
        % Y = u*s*v'.
        [v, ssquare] = eig(symm(Y'*(B*Y))); % Y*B*Y is positive definite
        ssquarevec = diag(ssquare);

        s = sparse(diag(abs(sqrt(ssquarevec))));
        u = Y*(v/s); % u'*B*u is identity.
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [z,uTx,uTarRx,uTarTx]=Cart2RuvStdRefrac(zC,useHalfRange,zTx,zRx,M,Ns,includeW,ce,rE,spherCent)
%%CART2RUVSTDREFRAC Convert points in Cartesian coordinates (either the
%         global system or a local system at the receiver) into local
%         bistatic r-u-v coordinates of the receiver, accounting for how a
%         standard exponential atmospheric model warps the measurements.
%         r-u-v coordinates consist of a bistatic range and direction
%         cosines at the receiver. The "direction cosines" u and v are just
%         the x and y coordinates of a unit vector from the receiver to
%         the target in the coordinate system at the receiver (Here, to the
%         refraction-corrupted apparent direction). This basically assumes
%         that the boresight direction of the receiver is the z axis.
%         Assuming the target is in front of the receiver, the third unit
%         vector coordinate is not needed. However, with the includeW
%         option, it can be provided, resulting in r-u-v-w coordinates. This
%         function is not suitable for computing refraction between
%         satellites, grazing the Earth's atmosphere. The algorithm might
%         have an error if the ray path goes too far underground.
%
%INPUT: zC  A 3XN matrix of Cartesian points (target locations) in global
%           [x;y;z] Cartesian coordinates.
% useHalfRange A boolean value specifying whether the bistatic range value
%           should be divided by two. This normally comes up when operating
%           in monostatic mode, so that the range reported is a one-way
%           range. The default if this parameter is not provided (or an
%           empty matrix is provided) is false.
%       zTx The 3X1 [x;y;z] location vector of the transmitter in global
%           Cartesian coordinates.  If this parameter is omitted or an
%           empty matrix is passed, then the receiver is placed at the origin.
%       zRx The 3X1 [x;y;z] location vector of the receiver in global
%           Cartesian coordinates. If this parameter is omitted or an empty
%           matrix is passed, then the receiver is placed at the origin.
%         M A 3X3 rotation matrix to go from the alignment of the global
%           coordinate system to the local alignment of the receiver. The z
%           vector of the local coordinate system of the receiver is the
%           pointing direction of the receiver. If this matrix is omitted,
%           then the identity matrix is used.
%        Ns The atmospheric refractivity reduced to the reference sphere.
%           Note that the refractivity is (n-1)*1e6, where n is the index
%           of refraction. The function reduceStdRefrac2Spher can be used
%           to reduce a refractivity to the surface of a reference
%           ellipsoid. This function does not allow different
%           refractivities to be used as the transmitter and receiver. If
%           this parameter is omitted or an empty matrix is passed, a
%           default value of 313 is used.
%  includeW An optional boolean value indicating whether a third direction
%           cosine component should be included. The u and v direction
%           cosines are two parts of a 3D unit vector. Generally, one might
%           assume that the target is in front of the sensor, so the third
%           component would be positive and is not needed. However, the
%           third component can be included if ambiguity exists. The
%           default if this parameter is omitted or an empty matrix is
%           passed is false.
%        ce The optional decay constant of the exponential model. The
%           refractivity N at height h is N=Ns*exp(-ce*(h-h0)) where h0 is
%           the reference height (in this function, the height of the
%           reference ellipsoid surface is used). ce is related to the
%           change in refractivity at an elevation of 1km based on the
%           refractivity at sea level as
%           ce=log(Ns/(Ns+DeltaN))/1000;%Units of inverse meters.
%           where the change in refractivity for a change in elevation of
%           1km is DeltaN=-multConst*exp(expConst*Ns); In [1], standard
%           values for the two constants are expConst=0.005577; and
%           multConst=7.32; If ce is omitted or an empty matrix is passed,
%           the value based on the standard model is used.
% rE,spherCent The radius of the Earth to use for the spherical Earth
%           approximation used in the model and also the offset between the
%           global model and the local spherical model. It is assumed that
%           zC,zTx,and zRx are all given in the global model and will need
%           to be transformed to the local model to the used. If rE is
%           omitted or an empty matrix is passed, then the default of
%           [rE,spherCent]=osculatingSpher4LatLon(Cart2Ellipse(zRx)) is
%           used. The defaults here mean that a WGS-84 reference ellipsoid
%           is approximated by the local osculating sphere.
%
%OUTPUTS: z The 3XN (or 4XN if includeW is true) matrix of location vectors
%           of the points in bistatic [r;u;v] coordinates. If
%           useHalfRange=true, then the r component is half the bistatic
%           range (half the round-trip range for a monostatic scenario).
%       uTx A 3XN set of unit vectors pointing from the transmitter to the
%           refraction-corrupted position (of the target as seen by the
%           transmitter. This is in the global coordinate system.
%    uTarRx A 3XN set of unit vectors pointing from the target to the
%           refraction-corrupted position of the receiver as seen by the
%           target. This is in the global coordinate system.
%    uTarTx A 3XN set of unit vectors pointing from the target to the
%           refraction-corrupted position of the transmitter as seen by the
%           target. This is in the global coordinate system.
%
%This function implements the refraction algorithm for the basic
%exponential atmosphere as described in [1] for the bistatic case. If the
%target is collocated with the transmitter or the receiver, then NaNs will
%be returned for some values. The basic exponential refraction model is in
%[2].
%
%The model is parameterized in terms of a height above a sphere. The Earth
%is more of an ellipsoid than a sphere. Thus, we use local spherical
%approximations about the transmitter and the receiver. That is, for
%computing the refraction from the transmitter to the target, we use the
%distance from the center of the Earth to the surface of the reference
%ellipsoid at the transmitter as the radius of an approximately spherical
%Earth. Similarly, the distance from the center of the Earth to the
%receiver is used in the approximation for the path from the target to the
%receiver.
%
%The algorithm in [1] performs ray tracing by solving a boundary value
%problem. here, the bvp5c function in Matlab is used to solve the problem.
%
%For paths that are nearly vertical, it is approximated that there is no
%bending in angle and an explicit solution to the integral over the index
%of refraction in the vertical direction is used to obtain the range.
%
%EXAMPLE:
%Here, we have two radars and one target near Hawaii.
% latLonRx=deg2rad([20.269202;-155.852051]);
% AltRx=0;
% latLonTx=deg2rad([20.724568;-155.978394]);
% AltTx=0;
% latLonTar=deg2rad([20.835390;-155.313721]);
% AltTar=8e3;%8km target altitude.
% %Convert locations to Cartesian.
% zRx=ellips2Cart([latLonRx;AltRx]);
% zTx=ellips2Cart([latLonTx;AltTx]);
% zTar=ellips2Cart([latLonTar;AltTar]);
%
% %The receiver faces 45 degrees East of North and 15 degrees up from the
% %local ellipsoidal level.
% M=findRFTransParam([latLonRx;AltRx],deg2rad(45),deg2rad(15));
% Ns=350;%Assumed refractivity at the sea surface.
% useHalfRange=false;
% includeW=true;%Include third dimension of unit vector.
% [z,uTx,uTarRx,uTarTx]=Cart2RuvStdRefrac(zTar,useHalfRange,zTx,zRx,M,Ns,includeW);
% zNoRefrac=Cart2Ruv(zTar,useHalfRange,zTx,zRx,M,includeW);
% z(1)-zNoRefrac(1)%Bistatic range difference of 31.0813 meters
% %Direction difference of 0.0927 degrees
% rad2deg(angBetweenVecs(z(2:end),zNoRefrac(2:end)))
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using 3D monostatic and bistatic
%    measurements in refractive environments," IEEE Aerospace and
%    Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 54-75, Aug.
%    2014.
%[2] B. R. Bean and G. D. Thayer, CRPL Exponential Reference Atmosphere.
%    Washington, D.C.: U. S. Department of Commerce, National Bureau of
%    Standards, Oct. 1959. [Online]. Available:
%    http://digicoll.manoa.hawaii.edu/techreports/PDF/NBS4.pdf
%
%June 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numMeas=size(zC,2);

%The assumes refractivity at sea level to use if none is provided.
if(nargin<7||isempty(includeW))
    includeW=false;
end

if(nargin<6||isempty(Ns))
    Ns=313;
end

if(nargin<5||isempty(M))
    M=eye(3);
end

if(nargin<4||isempty(zRx))
    zRx=zeros(3,1);
end

if(nargin<3||isempty(zTx))
    zTx=zeros(3,1);
end

if(nargin<2||isempty(useHalfRange))
    useHalfRange=false;
end

if(nargin<8||isempty(ce))
    expConst=0.005577;
    multConst=7.32;

    %The change in refractivity at an elevation of 1km based on the
    %refractivity on the surface of the Earth.
    DeltaN=-multConst*exp(expConst*Ns);
    ce=log(Ns/(Ns+DeltaN))/1000;%Units of inverse meters.
end

if(nargin<9||isempty(rE))
    %Use the radius of the Earth that is the radius of the osculating
    %sphere at the location of the observed. This will be the radius used
    %in the local spherical Earth approximation for computing atmospheric
    %refraction. This uses the WGS-84 reference ellipsoid.
    [rE,spherCent]=osculatingSpher4LatLon(Cart2Ellipse(zRx));
end

%Adjust all the Cartesian values based on the osculating sphere.
zC=zC-spherCent;
zTx=zTx-spherCent;
zRx=zRx-spherCent;

%Allocate space
if(includeW)
    z=zeros(4,numMeas);
else
    z=zeros(3,numMeas);
end
uTx=zeros(3,numMeas);
uTarTx=zeros(3,numMeas);
uTarRx=zeros(3,numMeas);

if(any(zRx~=zTx))%If the scenario is bistatic
    for curMeas=1:numMeas
        [r2,uArrive,uTarRx(:,curMeas)]=atmosRefracMeas(zRx,zC(:,curMeas),Ns,ce,rE);
        [r1,uTx(:,curMeas),uTarTx(:,curMeas)]=atmosRefracMeas(zTx,zC(:,curMeas),Ns,ce,rE);

        r=(r1+r2);
        u=M*uArrive;

        if(useHalfRange)
            r=r/2;
        end

        zCur=[r;u];

        if(includeW)
            z(:,curMeas)=zCur;
        else
            z(:,curMeas)=zCur(1:3);
        end
    end
else%The scenario is monostatic.
    for curMeas=1:numMeas
        [range,uArrive,uTarRx(:,curMeas)]=atmosRefracMeas(zTx,zC(:,curMeas),Ns,ce,rE);
        uTx=uArrive;
        r=2*range;%Round-trip range.
        u=M*uArrive;

        if(useHalfRange)
            r=r/2;
        end

        zCur=[r;u];

        if(includeW)
            z(:,curMeas)=zCur;
        else
            z(:,curMeas)=zCur(1:3);
        end

        uTx(:,curMeas)=uArrive;
        uTarTx(:,curMeas)=uTarRx(:,curMeas);
    end
end
end

function [range,uArrive,uDepart]=atmosRefracMeas(xObs,xObj,Ns,ce,rE)
%%ATMOSREFRACMEAS  Given the location of an observer and an object in the
%                  atmosphere of the Earth, find the delay and angle of
%                  arrival of a signal from the object to the observer,
%                  accounting for basic,standard  refraction. A
%                  low-fidelity exponential atmospheric model is used. This
%                  function is not suitable for computing refraction
%                  between satellites, grazing the Earth's atmosphere. The
%                  algorithm might have an error if the raypath
%                  goes too far underground.
%
%INPUTS: xObs The Cartesian location of the observer in ECEF coordinates in
%             meters as [x;y;z].
%        xObj The Cartesian location of the object being observed in
%             ECEF coordinates in meters as [x;y;z].
%          Ns The refractivity as reduced to at sea level.
% expConst,multConst The parameters of the refractivity model such that
%             increasing the height by 1km, the model is that the
%             refractivity changes by deltaN=-multConst*exp(expConst*N).
%             deltaN cannot be negative.
%          rE The radius of the Earth to use in the spherical Earth
%             approximation. An osculating sphere near a sensor is
%             suggested.
%
%OUTPUTS: range  The apparent one-way range (in meters) of the signal from
%                the transmitter to the target and back to the receiver.
%        uArrive A unit vector in ECEF coordinates pointing in the apparent
%                direction of the signal the observer received (as seen by
%                the observer).
%        uDepart The direction of the signal departing the object that
%                arrives at the observer. Put another way, if the observer
%                were to transmit a signal to the object, this is the
%                apparent direction of the observer as seen by the object.
%
%This function implements the refraction algorithm for the basic
%exponential atmosphere as described in [1] for the monostatic case. If the
%points are collocated, then NaNs are returned for uArrive and uDepart.
%
%The function will fail for paths that go too deep into the Earth.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using 3D monostatic and bistatic
%    measurements in refractive environments," IEEE Aerospace and
%    Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 54-75, Aug.
%    2014.
%
%May 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%We need the conversion from the 3D coordinate system of the observer and
%object into the 2D coordinate system used for raytracing. The 2D
%coordinate system has the center of the Earth as its origin and the x-y
%axes are in the plane of the vector from the observer to the target. One
%vector common to both coordinate systems in the local up vector, which
%will be the local y axis. The second vector common to both will be the
%local x vector, which will be the projection of xObj-xObs onto the local
%tangent plane. Here, the vertical is the spherical model vertical. Since
%the precision of the model is low enough that the difference between the
%spherical and gravitational verticals shouldn't matter.
uENU=getENUAxes(Cart2Ellipse(xObs,[],rE,0));
uVertGlobal=uENU(:,3);
uVertLocal=[0;1;0];

vec2TarGlobal=xObj-xObs;

%The projection of the xObj-xObs vector into the local tangent plane can be
%obtained by subtracting the component of the vector that is orthogonal to
%the plane.
uHorizGlobalOrig=vec2TarGlobal-dot(vec2TarGlobal,uENU(:,3))*uENU(:,3);

uHorizGlobal=uHorizGlobalOrig/norm(uHorizGlobalOrig);
uHorizLocal=[1;0;0];

%Find the rotation matrix from the global coordinate system into the local
%coordinate system.
ECEF2LocalRot=findTransParam([uVertLocal,uHorizLocal],[uVertGlobal,uHorizGlobal]);

%The third (z) coordinate in the local system should be zero after this
%transformation.
vec2TarLocal=ECEF2LocalRot*vec2TarGlobal;

%The location of the receiver in the local 2D coordinate system.
x0Init=0;
y0Init=norm(xObs);

%The location of the target in the local 2D coordinate system.
x1Init=vec2TarLocal(1);
y1Init=vec2TarLocal(2)+y0Init;

%If the two points are nearly vertical, then the ray tracing algorithm will
%fail. For nearly vertical points, the bending due to refraction in the
%model should be negligible, so we can perform an integral in the y
%direction to solve for the excess range instead of having to solve the
%more complicated general bistatic problem. To deal with x1Init not being
%exactly zero, we actually go to a full range of norm([x1Init;y1Init])
if(norm(uHorizGlobalOrig)<1e-3)
    %Integrating the index of refraction from y0Init to
    %norm([x1Init;y1Init]) at a constant x yields the following measured
    %range. It is not just the geometric range, because the index of
    %refraction is not a constant 1.
    yMax=norm([x1Init;y1Init]);

    range=((exp(ce*(rE-y0Init))-exp(ce*(rE-yMax)))*Ns)/(1e6*ce)+yMax-y0Init;

    uArrive=vec2TarGlobal/norm(vec2TarGlobal);
    uDepart=-uArrive;
    return;
end

%Now, set up the boundary-value problem to determine the path taken by
%light between the target and the receiver.

%The initial guess is just the linear solution. The solver requires a fixed
%number of steps. 20 is probably sufficient for things near the Earth. that
%is, up to distances of, say 400km. We can scale the number of steps as 20
%for every 400 kilometers with a minimum of, say 10.
%Things outside of the atmosphere should use the astronomical refraction
%routines.
numSteps=max(20,ceil(20*norm(vec2TarLocal)/400e3));
x=linspace(x0Init,x1Init,numSteps);
slope=(y1Init-y0Init)/(x1Init-x0Init);
b=y1Init-slope*x1Init;%The y-intercept.
%The initial estimate of the solution.
solInit=bvpinit(x,@(x)[x*slope+b;slope]);

%Now, solve the differential equation.
oldOpts=bvpset();
newOpts=bvpset(oldOpts,'RelTol',1e-8,'AbsTol',1e-8,'FJacobian',@(x,y)odefunJacob(x,y,Ns,rE,ce),'BCJacobian',@bcfunJacob);%Increase the accuracy.
sol=bvp5c(@(x,y)expDiffEq(x,y,Ns,rE,ce),@(y0,y1)bcfun(y0,y1,y0Init,y1Init),solInit,newOpts);

%Get the refraction-corrupted range measurement for a signal traveling from
%the object to the observer.
range=integral(@(x)pathFun2D(x,sol,Ns,rE,ce),x0Init,x1Init,'AbsTol',eps(1),'RelTol',1e-15);

if(nargout>1)
    %Get the angle of arrival for a signal traveling from the object to the
    %observer. The angle is determined by the slope at the initial point.
    thetaOrig=atan(sol.y(2,1));
    uLocal=[cos(thetaOrig);sin(thetaOrig);0];

    %The inverse rotation is given the transpose of the rotation
    %matrix. This is the apparent direction of the object as seen by the
    %observer.
    uArrive=ECEF2LocalRot'*uLocal;

    thetaEnd=atan(sol.y(2,end));
    uLocal=[cos(thetaEnd);sin(thetaEnd);0];
    %This is the apparent direction of the observer as seen by the object.
    uDepart=-ECEF2LocalRot'*uLocal;
end
end

function val=pathFun2D(x,sol,Ns,rE,ce)
    %This function is used to integrate the time taken
    y=deval(x,sol);
    val=(1+NRefracExp(x,y(1,:),Ns,rE,ce)).*sqrt(1+y(2,:).^2);
end

function res=bcfun(y0,y1,y0Init,y1Init)
    %The residue to define the boundary condition for the numeric
    %differential equation solver as applied to the 2D exponential
    %atmospheric model.

    res=[y0(1)-y0Init;
         y1(1)-y1Init];
end

function [dbcy0,dbcy1]=bcfunJacob(~,~)
    %The Jacobians of the boundary conditions for raytracing the 2D
    %exponential atmospheric refraction model.
    dbcy0=[1 0
           0 0];
    dbcy1=[0 0
           1 0];
end

function J=odefunJacob(x,y,Ns,rE,ce)
    %The Jacobian of the differential equation for raytracing the 2D
    %exponential atmospheric model.
    expVal=NRefracExp(x,y(1),Ns,rE,ce);

    J=zeros(2,2);
    J(1,2)=1;
    J(2,1)=ce*(1+y(2)^2)*(-expVal)*(ce*y(1)*(x*y(2)-y(1))*sqrt(x^2+y(1)^2)+x*(x+y(1)*y(2))*(expVal+1))/((x^2+y(1)^2)^(3/2)*(expVal+1)^2);
    J(2,2)=ce*(x-2*y(1)*y(2)+3*x*y(2)^2)*expVal/((expVal+1)*sqrt(x^2+y(1)^2));
end

function dxdy=expDiffEq(x,y,Ns,rE,ce)
    %Find the refractivity at location (x,y).
    expVal=NRefracExp(x,y(1),Ns,rE,ce);

    dxdy=[y(2)
          ce*(1+y(2)^2)*(x*y(2)-y(1))*expVal/((expVal+1)*sqrt(x^2+y(1)^2))];
end

function nRefrac=NRefracExp(x,y,Ns,rE,ce)
    %The refractivity. This is 10^6*(index of refraction-1)
    nRefrac=1e-6*Ns*exp(-ce*(sqrt(x.^2+y.^2)-rE));
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = offset(varargin)
% file_array's offset property
% For getting the value
% dat = offset(obj)
%
% For setting the value
% obj = offset(obj,dat)
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: offset.m 7147 2017-08-03 14:07:01Z spm $

if nargin==2
    varargout{1} = asgn(varargin{:});
elseif nargin==1
    varargout{1} = ref(varargin{:});
else
    error('Wrong number of arguments.');
end

%==========================================================================
% function dat = ref(obj)
%==========================================================================
function dat = ref(obj)
dat = obj.offset;

%==========================================================================
% function obj = asgn(obj,dat)
%==========================================================================
function obj = asgn(obj,dat)
if isnumeric(dat) && numel(dat)==1 && dat>=0 && rem(dat,1)==0
    obj.offset = double(dat);
else
    error('"offset" must be a positive integer.');
end

author: spm


text:
%Primitives Primitives enumeration for the brick
%
% Notes::
% - Primitives can be found in the EV3 documentation and source code
% (bytecodes.h 1584-1606)

classdef Primitives < uint8
    enumeration
        pshort (0) % PRIMPAR_SHORT = 0x00
        plong (128) % PRIMPAR_LONG = 0x80

        pconst(0) % PRIMPAR_CONST = 0x00
        pvariabel(64) % PRIMPAR_VARIABEL = 0x40
        plocal(0) % PRIMPAR_LOCAL = 0x00
        pglobal(32) % PRIMPAR_GLOBAL = 0x20
        phandle(16) % PRIMPAR_HANDLE = 0x10
        paddr(8) % PRIMPAR_ADDR = 0x08

        pindex(31) % PRIMPAR_INDEX = 0x1F
        pconst_sign(32) % PRIMPAR_CONST_SIGN = 0x20
        pvalue (63) % PRIMPAR_VALUE = 0x3F

        pbytes(7) % PRIMPAR_BYTES = 0x07

        pstring_old (0) % PRIMPAR_STIRNG_OLD = 0
        p1_byte (1) % PRIMPAR_1_BYTE = 1
        p2_byte (2) % PRIMPAR_2_BYTE = 2
        p4_byte (3) % PRIMPAR_4_BYTE = 3
        pstring (4) % PRIMPAR_STRING = 4

        plabel (32) % PRIMPAR_LABEL = 0x20
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
RadarName = 'staddon_clutter';
RadarCoords.lat = 50.346069;
RadarCoords.lon = -4.113670;

RadarCoords_lat = 50.339026;
RadarCoords_lon = -4.126307;

% Load dataset
load(strcat(RadarName,'.mat'));
load('gnn_llr_clutter_a_fakephd.mat');

lon_lim = [-4.176846509809199,-4.103030078951917];
lat_lim = [50.297532151330770,50.337022151330770];

lon_lim_obs = [-4.16737200000000,-4.14897200000000];
lat_lim_obs = [50.3537860000000,50.3617750000000];
% radar_meas_convert;

% Map plot
figure('units','normalized','outerposition',[0 0 .5 0.8])
ax(1) = gca;
plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
plots = [];
hold on;

N = size(MeasurementScans,2); % Simulation length

Tracks = [TrackList, DeletedTracks];

%% PLot detections
for k=2:N

%     MeasurementList = MeasurementScans(k);
%     if(MeasurementList.NumMeasurements>0)
%         % Convert measurements to LLA and plot them
%     %     [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
%     %                                MeasurementList.Vectors(1,:),...
%     %                                0,...
%     %                                RadarCoords.lat,...
%     %                                RadarCoords.lon,...
%     %                                0,...
%     %                                referenceEllipsoid('wgs84'));
%         latlon = [MeasurementList.Metadata.LatLon];
%         lat = latlon(1,:);
%         lon = latlon(2,:);
%         plots(end+1) = plot(ax(1), lon,lat,'y.','MarkerSize', 4);
%     end
end

%% Plot all existing tracks

for j=1:numel(Tracks)
    track = Tracks{j};
%     if find(t==track.Tag.ID,1)
%         continue
%     end
    % Convert track trajectory to LLA and plot it
    means = [track.Trajectory.Mean];
    [lat,lon,~] = ned2geodetic(means(3,:),...
                               means(1,:),...
                               0,...
                               RadarCoords.lat,...
                               RadarCoords.lon,...
                               0,...
                               referenceEllipsoid('wgs84'));
    traj_length = size(lon,2);

    plots(end+1) = plot(ax(1), lon(:,1),lat(:,1),'go','MarkerSize',15,'LineWidth',2);
    t = [51806, 6423, 8579, 82699, 37841];
    td = [28571, 57088, 46505, 50816];
    if (find(t == track.Tag.ID,1))
        plots(end+1) = plot(ax(1), lon, lat,'-','LineWidth',4,'Color', [0.9856, 0.7372, 0.2537]);
        plots(end+1) = plot(ax(1), lon, lat,'-.w','LineWidth',1);
    elseif(find(td==track.Tag.ID,1))
        % Delayed initiation
        plots(end+1) = plot(ax(1), lon, lat,'-','LineWidth',4,'Color', [0.0265, 0.6137, 0.8135]);
        plots(end+1) = plot(ax(1), lon, lat,'-.w','LineWidth',1);
    else
        plots(end+1) = plot(ax(1), lon, lat,'-.w','LineWidth',2);
    end
    plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'rx','MarkerSize',15,'LineWidth',2);
%     text(ax(1), lon(:,1), lat(:,1), num2str(track.Tag.ID),'FontSize',18,'Color','g');

    % Convert track velocity to LLA and plot it
%     [lat_vel,lon_vel,~] = ned2geodetic(track.State.Mean(4,end),...
%                                        track.State.Mean(2,end),...
%                                        0,...
%                                        lat(:,end),...
%                                        lon(:,end),...
%                                        0,...
%                                        referenceEllipsoid('wgs84'));
%     lat_vel = lat_vel-lat(:,end);
%     lon_vel = lon_vel-lon(:,end);
%     plots(end+1) = quiver(ax(1), lon(:,end),lat(:,end),20*lon_vel,20*lat_vel,'r','LineWidth',1.5);
%     if ShowTrackInfo
%     end

end

%% Plot region of high clutter
x = [-469.710602736631,-2572.92295458513,-2516.31234127506,-1320.74163107899,-835.307363807233, -469.710602736631];
y = [-4915.41679013513,-4065.50122858976,-2858.60918927224,-1380.19032924249,-1261.08679763585, -4915.41679013513];
%         x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
%         y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
[y,x,~] = ned2geodetic(y,...
                       x,...
                       0,...
                       RadarCoords.lat,...
                       RadarCoords.lon,...
                       0,...
                       referenceEllipsoid('wgs84'));
plots(end+1) = plot(ax(1), x, y, 'r-', 'LineWidth', 2);

plot(ax(1), RadarCoords_lon,RadarCoords_lat,...
                               '-s','MarkerSize',10,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);
%% Draw fragmented tracks
pos = [-4.129656720082580,50.304822261284470,0.004444050429163,0.002011064814816];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
pos = [-4.167167518620260,50.307747446469660,0.003464852876975,0.002513831018518];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
pos = [-4.169427205279155,50.312455166377070,0.003163561322454,0.002056770833327];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
pos = [-4.129656720082580,50.302994020543736,0.002485655324787,0.001508298611107];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
% x = [-4.12945334365527,-4.12917944224207,-4.12739908305628,-4.12753603376287,-4.12945334365527];
% y = [50.3032914000914,50.3022941778692,50.3023772797210,50.3034576037951,50.3032914000914];
% plot(ax(1), x,y, 'c-','LineWidth', 2);
%% PLot false tracks
pos = [-4.130334626080249,50.326669738136324,0.004293404651903,0.004753425925927];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'm');

%% Legend
h1 = plot(NaN, NaN, 'o','MarkerSize', 5, 'MarkerFaceColor', 'y', 'MarkerEdgeColor','k');
h2 = plot(NaN, NaN, 's','MarkerSize',10,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);
h3 = plot(NaN, NaN, 'ms', 'MarkerSize',15,'LineWidth', 2);
h4 = plot(NaN, NaN, 'r-', 'LineWidth', 2);
h5 = plot(NaN, NaN, 'go','MarkerSize',15,'LineWidth',2);
h6 = plot(NaN, NaN, 'rx','MarkerSize',15,'LineWidth',2);
h7 = plot(NaN, NaN, 'cs','MarkerSize',15,'LineWidth',2);
h8 = plot(NaN, NaN, 'b-','MarkerSize',15,'LineWidth',4);
h9 = plot(NaN, NaN, '-','LineWidth',4,'Color', [0.9856, 0.7372, 0.2537]);
h10 = plot(NaN, NaN, '-', 'LineWidth',4, 'Color', [0.0265, 0.6137, 0.8135]);
legend([h4, h5, h6, h3, h7, h9, h10], 'High-clutter Region', 'Track Birth', 'Track Death','False Tracks', 'Fragmentation Examples', 'Bonus Tracks', 'Delayed Initiation', 'Location','northeast');
% legend([h4, h5, h6], 'High-clutter Region', 'Track Birth', 'Track Death','Location','northeast');
% legend([h1, h3, h5, h6, h9], 'Detections', 'Obscured Region', 'Track Birth', 'Track Death', 'Delayed Initiation', 'Location','southeast');
% legend([h1, h2, h3, h4], 'Detections', 'Radar Position', 'Obscured Region', 'High-clutter region');
% legend([h1, h3, h5, h6, h7, h8], 'Detections', 'Obscured Region', 'Track Birth', 'Track Death', 'Fragmentation Examples', 'Jump-track','Location','southeast');
% legend([h1, h2, h3, h4], 'Detections', 'Radar Position', 'Obscured Region', 'High-clutter region');
% xlabel('Longitude');
% ylabel('Latitude');
set(ax(1),'YTickLabel',[]);
set(ax(1),'XTickLabel',[]);

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function varargout=lpcdb2pf(varargin)
%
% For calling details please see v_lpcdb2pf.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcdb2pf.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcdb2pf.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcdb2pf(varargin{:});
else
    v_lpcdb2pf(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


function [Phip,Sigmap] = reorderVAR(Phi,Sigma,reordering)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018

% Permute the autoregressive matrix and the variance coveriance of the
% shocks with the new ordering of the variables.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N           = size(Sigma,1);

if length(reordering) ~= N,
    error('Not enough permutations: the vector with new ordering of variables must have size N.')
end

Phip    = Phi(:,reordering);
Sigmap  = nan(N);

for jj  = 1 : length(reordering)
    for hh = 1 : length(reordering)

        tmp = Sigma(reordering(jj),reordering(hh));
        Sigmap(jj,hh) = tmp;
    end

end

end


text:
function f = randarb(x,y)
% RANDARB generates a random observation from any arbitrary PDF defined by x,y
% function f = hfarbrand(x,y) x and y are vectors of length N that describe a PDF
% to some precision implicitly dictated by the size of N.  Fhe returned scalar f
% is an observation from the set of x with a probability of (x/y(x))/sum(y)

%% Sanity checks
if nargin ~= 2
    error('Two input vectors are required');
end
if length(x) ~= length(y)
    error('x and y must be of the same length');
end
if length(x) < 30
    warning('The size of x and y is very low');
end

%% Compute a uniform random number between 0 and sum(y)
randx = sum(y)*rand();

%% Find where the number lies on a conceptual line comprised of "segments" of
% length y
i = 1;
while sum(y(1:i)) < randx
    i = i + 1;
end

%% Return the x value corresponding to the y value we "landed on."
f = x(i);

author: Sable

repo: mcbench-benchmarks


text:
function varargout = PID_Controller(varargin)

global sys_tf sys_ss sys_fd_P sys_fd_PD sys_fd_PI sys_fd_PID
global a_dou b_dou c_dou d_dou
global sys_final_P sys_final_PI sys_final_PD sys_final_PID

% PID_CONTROLLER M-file for PID_Controller.fig
%      PID_CONTROLLER, by itself, creates a new PID_CONTROLLER or raises the existing
%      singleton*.
%
%      H = PID_CONTROLLER returns the handle to a new PID_CONTROLLER or the handle to
%      the existing singleton*.
%
%      PID_CONTROLLER('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in PID_CONTROLLER.M with the given input arguments.
%
%      PID_CONTROLLER('Property','Value',...) creates a new PID_CONTROLLER or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before PID_Controller_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to PID_Controller_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help PID_Controller

% Last Modified by GUIDE v2.5 04-Sep-2008 09:14:05

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @PID_Controller_OpeningFcn, ...
                   'gui_OutputFcn',  @PID_Controller_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before PID_Controller is made visible.
function PID_Controller_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to PID_Controller (see VARARGIN)

% Choose default command line output for PID_Controller
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes PID_Controller wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = PID_Controller_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in ss_pushb.
function ss_pushb_Callback(hObject, eventdata, handles)
% hObject    handle to ss_pushb (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
ss_pushbutton = get (handles.ss_radiob,'Value');
tf_pushbutton = get (handles.tf_radiob,'Value');

%================= Implementing State Space Model ====================
if ss_pushbutton == get(handles.ss_radiob,'Max')

    a = get(handles.A_edit,'String');
    global a_dou
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    global b_dou
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    global c_dou
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    global d_dou
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ================

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D in correct format');
    end

    %======= Error If the State Space Model is not correct ================
    [r1 c1] = size(a_dou);
    [r2 c2] = size(c_dou);
    [r3 c3] = size(b_dou);
    [r4 c4] = size(d_dou);
        if c1 ~= c2
            errordlg('Please Enter Valid State Space Model');
        end
        if r1 ~= r3
            errordlg('Please Enter Valid State Space Model');
        end
        if c3 ~= c4
            errordlg('Please Enter Valid State Space Model');
        end

    global sys_ss
    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);            % State Space Model

     [s t]= step(sys_ss);
     axes(handles.axes1)
     plot(t,s),grid on

    %================= Error If State Space Button is not selected ========
else
    errordlg('Select "State Space"');
end

    %======================= Time Response ================================

%     sys_s = tf(a_dou,b_dou,c_dou,d_dou)
%     sys_s = tf(n,d)
    sys_s= tf(sys_ss);
    [num,den] = tfdata(sys_s);
    [r c] = size(den{1,1});

    if c ~= 3
        set(handles.MOS_text,'String','Approximation Not Valid');
        set(handles.PTS_text,'String','Approximation Not Valid');
        set(handles.STS_text,'String','Approximation Not Valid');
    end

    roots_sys = roots(den{1,1});
    roots_sys_str = num2str(roots_sys);

    set(handles.ROS_text,'String',roots_sys_str);

    den_sp = den{1,1}(1,3);
    wn = sqrt(den_sp);
    den_sp2 = den{1,1}(1,2);
    b = den_sp2/(2*wn);
    [r c] = size(den{1,1});

    if b<=1 && num{1,1}(1,end-1) == 0 && c == 3
        MOS = exp(-(b*pi)/(sqrt(1-b^2)))*100;
        MOS_str = num2str(MOS);
        MOS_real = strcat(MOS_str,' %');
        set(handles.MOS_text,'String',MOS_real);

        PTS = pi/(wn*(sqrt(1-b^2)));
        PTS_str = num2str(PTS);
        PTS_real = strcat(PTS_str,' Sec');
        set(handles.PTS_text,'String',PTS_real);

        STS = 4/(b*wn);
        STS_str = num2str(STS);
        STS_real = strcat(STS_str,' Sec');
        set(handles.STS_text,'String',STS_real);
    else
        set(handles.MOS_text,'String','Approximation Not Valid');
        set(handles.PTS_text,'String','Approximation Not Valid');
        set(handles.STS_text,'String','Approximation Not Valid');
    end

%==========================================================================

function A_edit_Callback(hObject, eventdata, handles)
% hObject    handle to A_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of A_edit as text
%        str2double(get(hObject,'String')) returns contents of A_edit as a double

% --- Executes during object creation, after setting all properties.
function A_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to A_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function B_edit_Callback(hObject, eventdata, handles)
% hObject    handle to B_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of B_edit as text
%        str2double(get(hObject,'String')) returns contents of B_edit as a double

% --- Executes during object creation, after setting all properties.
function B_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to B_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function C_edit_Callback(hObject, eventdata, handles)
% hObject    handle to C_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of C_edit as text
%        str2double(get(hObject,'String')) returns contents of C_edit as a double

% --- Executes during object creation, after setting all properties.
function C_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to C_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function D_edit_Callback(hObject, eventdata, handles)
% hObject    handle to D_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of D_edit as text
%        str2double(get(hObject,'String')) returns contents of D_edit as a double

% --- Executes during object creation, after setting all properties.
function D_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to D_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in tf_pushb.
function tf_pushb_Callback(hObject, eventdata, handles)
% hObject    handle to tf_pushb (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
ss_pushbutton = get (handles.ss_radiob,'Value');
tf_pushbutton = get (handles.tf_radiob,'Value');

%================= Implementing Transfer Function ====================

if tf_pushbutton == get(handles.tf_radiob,'Max')

    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty  ====================
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please enter in right format, For Example Num = [1 2] or Den = [1 2 3]');
    end
    % =====================================================================

    global sys_tf
    sys_tf = tf(num_dou,den_dou);            % Transfer Function

    [s t]= step(sys_tf);
    axes(handles.axes1)
    plot(t,s),grid on

    %================= Error If Transfer Function Button is not selected=======

else
    errordlg('Select "Transfer Function');
end

    %======================= Time Response ================================

    [num,den] = tfdata(sys_tf);
    [r c] = size(den{1,1});

    if c ~= 3
        set(handles.MOS_text,'String','Approximation Not Valid');
        set(handles.PTS_text,'String','Approximation Not Valid');
        set(handles.STS_text,'String','Approximation Not Valid');
    end

    roots_sys = roots(den{1,1});
    roots_sys_str = num2str(roots_sys);

    set(handles.ROS_text,'String',roots_sys_str);

    den_sp = den{1,1}(1,3);
    wn = sqrt(den_sp);
    den_sp2 = den{1,1}(1,2);
    b = den_sp2/(2*wn);
    [r c] = size(den{1,1});

    if b<=1 && num{1,1}(1,end-1) == 0 && c == 3
        MOS = exp(-(b*pi)/(sqrt(1-b^2)))*100;
        MOS_str = num2str(MOS);
        MOS_real = strcat(MOS_str,' %');
        set(handles.MOS_text,'String',MOS_real);

        PTS = pi/(wn*(sqrt(1-b^2)));
        PTS_str = num2str(PTS);
        PTS_real = strcat(PTS_str,' Sec');
        set(handles.PTS_text,'String',PTS_real);

        STS = 4/(b*wn);
        STS_str = num2str(STS);
        STS_real = strcat(STS_str,' Sec');
        set(handles.STS_text,'String',STS_real);
    else
        set(handles.MOS_text,'String','Approximation Not Valid');
        set(handles.PTS_text,'String','Approximation Not Valid');
        set(handles.STS_text,'String','Approximation Not Valid');
    end
%==========================================================================

function num_edit_Callback(hObject, eventdata, handles)
% hObject    handle to num_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of num_edit as text
%        str2double(get(hObject,'String')) returns contents of num_edit as a double

% --- Executes during object creation, after setting all properties.
function num_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to num_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function den_edit_Callback(hObject, eventdata, handles)
% hObject    handle to den_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of den_edit as text
%        str2double(get(hObject,'String')) returns contents of den_edit as a double

% --- Executes during object creation, after setting all properties.
function den_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to den_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in PID_popupmenu.
function PID_popupmenu_Callback(hObject, eventdata, handles)
% hObject    handle to PID_popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns PID_popupmenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from PID_popupmenu

% --- Executes during object creation, after setting all properties.
function PID_popupmenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to PID_popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Kp_edit_Callback(hObject, eventdata, handles)
% hObject    handle to Kp_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Kp_edit as text
%        str2double(get(hObject,'String')) returns contents of Kp_edit as a double

% --- Executes during object creation, after setting all properties.
function Kp_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Kp_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Ki_edit_Callback(hObject, eventdata, handles)
% hObject    handle to Ki_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Ki_edit as text
%        str2double(get(hObject,'String')) returns contents of Ki_edit as a double

% --- Executes during object creation, after setting all properties.
function Ki_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Ki_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Kd_edit_Callback(hObject, eventdata, handles)
% hObject    handle to Kd_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Kd_edit as text
%        str2double(get(hObject,'String')) returns contents of Kd_edit as a double

% --- Executes during object creation, after setting all properties.
function Kd_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Kd_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in PID_pushb.
function PID_pushb_Callback(hObject, eventdata, handles)
% hObject    handle to PID_pushb (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%===== Values from Controller popup menu and Responses menu ===============

value = get(handles.PID_popupmenu,'Value');

%====================== Values from Radio Button ======================

ss_pushbutton = get (handles.ss_radiob,'Value');
tf_pushbutton = get (handles.tf_radiob,'Value');

%==================== Selecting Only Proportional Controller %===============

if value == 1
    value_str = get(handles.Kp_edit,'String');
    value_dou = str2num (value_str);

    if tf_pushbutton == get(handles.tf_radiob,'Max');

        global sys_tf
        sys_final_P = sys_tf * value_dou;
        global sys_fd_P
        sys_fd_P = feedback(sys_final_P,1);

    end

    if ss_pushbutton == get(handles.ss_radiob,'Max')

        %================= Error If Either A,B,C or D is empty ================
    global a_dou b_dou c_dou d_dou
    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D in correct format');
    end

    %======= Error If the State Space Model is not correct ================
    [r1 c1] = size(a_dou);
    [r2 c2] = size(c_dou);
    [r3 c3] = size(b_dou);
    [r4 c4] = size(d_dou);
        if c1 ~= c2
            errordlg('Please Enter Valid State Space Model');
        end
        if r1 ~= r3
            errordlg('Please Enter Valid State Space Model');
        end
        if c3 ~= c4
            errordlg('Please Enter Valid State Space Model');
        end
        %==================================================================
        global a_dou b_dou c_dou d_dou
        [num den] = ss2tf(a_dou, b_dou ,c_dou ,d_dou);
        sys_tf_ss = tf(num,den);
        sys_final_P = sys_tf_ss * value_dou;
        global sys_fd_P
        sys_fd_P = feedback(sys_final_P,1);

    end

        [s t]= step(sys_final_P);
        axes(handles.axes5)
        plot(t,s),grid on
        set(handles.OLSR_text,'String','P Controller');

        global sys_fd_P
        [s t]= step(sys_fd_P);
        axes(handles.axes6)
        plot(t,s),grid on
        set(handles.CLSR_text,'String','P Controller');

end

%==================== Selecting Proportional Integral Controller ===============

if value == 2
    value_str = get(handles.Kp_edit,'String');
    value_dou = str2num (value_str);

    value_strPi = get(handles.Ki_edit,'String');
    value_douPi = str2num (value_strPi);
    value_tfpi = tf(value_douPi,[1 0]);

    if tf_pushbutton == get(handles.tf_radiob,'Max')

        global sys_tf
        PI = value_dou + value_tfpi;
        sys_final_PI = sys_tf * PI;
        global sys_fd_PI
        sys_fd_PI = feedback(sys_final_PI,1);

    end

    if ss_pushbutton == get(handles.ss_radiob,'Max')

         %================= Error If Either A,B,C or D is empty ================
    global a_dou b_dou c_dou d_dou
    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D in correct format');
    end

    %======= Error If the State Space Model is not correct ================
    [r1 c1] = size(a_dou);
    [r2 c2] = size(c_dou);
    [r3 c3] = size(b_dou);
    [r4 c4] = size(d_dou);
        if c1 ~= c2
            errordlg('Please Enter Valid State Space Model');
        end
        if r1 ~= r3
            errordlg('Please Enter Valid State Space Model');
        end
        if c3 ~= c4
            errordlg('Please Enter Valid State Space Model');
        end
        %==================================================================
        global a_dou b_dou c_dou d_dou
        [num den] = ss2tf(a_dou, b_dou ,c_dou ,d_dou);
        sys_tf_ss = tf(num,den);
        PI = value_dou + value_tfpi;
        sys_final_PI = sys_tf_ss * PI;
        global sys_fd_PI
        sys_fd_PI = feedback(sys_final_PI,1);
    end

     [s t]= step(sys_fd_PI);
     axes(handles.axes6)
     plot(t,s),grid on
     set(handles.CLSR_text,'String','PI Controller');

     [s t]= step(sys_final_PI);
     axes(handles.axes5)
     plot(t,s),grid on
     set(handles.OLSR_text,'String','PI Controller');

end

%==================== Selecting Proportional Derivative Controller ===============

if value == 3
    value_str = get(handles.Kp_edit,'String');
    value_dou = str2num (value_str);

    value_strPd = get(handles.Kd_edit,'String');
    value_douPd = str2num (value_strPd);
    value_tfpd = tf([value_douPd 0],1);

    if tf_pushbutton == get(handles.tf_radiob,'Max')

        global sys_tf
        PD = value_dou + value_tfpd;
        sys_final_PD = sys_tf * PD;
        global sys_fd_PD
        sys_fd_PD = feedback(sys_final_PD,1);
    end

    if ss_pushbutton == get(handles.ss_radiob,'Max')

         %================= Error If Either A,B,C or D is empty ================
    global a_dou b_dou c_dou d_dou
    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D in correct format');
    end

    %======= Error If the State Space Model is not correct ================
    [r1 c1] = size(a_dou);
    [r2 c2] = size(c_dou);
    [r3 c3] = size(b_dou);
    [r4 c4] = size(d_dou);
        if c1 ~= c2
            errordlg('Please Enter Valid State Space Model');
        end
        if r1 ~= r3
            errordlg('Please Enter Valid State Space Model');
        end
        if c3 ~= c4
            errordlg('Please Enter Valid State Space Model');
        end
        %==================================================================
        global a_dou b_dou c_dou d_dou
        [num den] = ss2tf(a_dou, b_dou ,c_dou ,d_dou);
        sys_tf_ss = tf(num,den);
        PD = value_dou + value_tfpd;
        sys_final_PD = sys_tf_ss * PD;
        global sys_fd_PD
        sys_fd_PD = feedback(sys_final_PD,1);
    end

    [s t]= step(sys_fd_PD);
    axes(handles.axes6)
    plot(t,s),grid on
    set(handles.CLSR_text,'String','PD Controller');

    [s t]= step(sys_final_PD);
    axes(handles.axes5)
    plot(t,s),grid on
    set(handles.OLSR_text,'String','PD Controller');
end

%==================== Selecting PID Controller ===============

if value == 4
    value_str = get(handles.Kp_edit,'String');
    value_dou = str2num (value_str);

    value_strPd = get(handles.Kd_edit,'String');
    value_douPd = str2num (value_strPd);
    value_tfpd = tf([value_douPd 0],1);

    value_strPi = get(handles.Ki_edit,'String');
    value_douPi = str2num (value_strPi);
    value_tfpi = tf(value_douPi,[1 0]);

     if tf_pushbutton == get(handles.tf_radiob,'Max')

        global sys_tf
        PID = value_dou + value_tfpi + value_tfpd;
        sys_final_PID = sys_tf * PID;
        global sys_fd_PID
        sys_fd_PID = feedback(sys_final_PID,1);
    end

    if ss_pushbutton == get(handles.ss_radiob,'Max')

         %================= Error If Either A,B,C or D is empty ================
    global a_dou b_dou c_dou d_dou
    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D in correct format');
    end

    %======= Error If the State Space Model is not correct ================
    [r1 c1] = size(a_dou);
    [r2 c2] = size(c_dou);
    [r3 c3] = size(b_dou);
    [r4 c4] = size(d_dou);
        if c1 ~= c2
            errordlg('Please Enter Valid State Space Model');
        end
        if r1 ~= r3
            errordlg('Please Enter Valid State Space Model');
        end
        if c3 ~= c4
            errordlg('Please Enter Valid State Space Model');
        end
        %==================================================================
        global a_dou b_dou c_dou d_dou
        [num den] = ss2tf(a_dou, b_dou ,c_dou ,d_dou);
        sys_tf_ss = tf(num,den);
        PID = value_dou + value_tfpi + value_tfpd;
        sys_final_PID = sys_tf_ss * PID;
        global sys_fd_PID
        sys_fd_PID = feedback(sys_final_PID,1);
    end

    [s t]= step(sys_fd_PID);
    axes(handles.axes6)
    plot(t,s),grid on
    set(handles.CLSR_text,'String','PID Controller');

    [s t]= step(sys_final_PID);
    axes(handles.axes5)
    plot(t,s),grid on
    set(handles.OLSR_text,'String','PID Controller');

end
% --------------------------------------------------------------------
function SRWC_Callback(hObject, eventdata, handles)
% hObject    handle to SRWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Error If both Options Selected ====================
if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Max')
    errordlg('Select One Option from "State Space" or "Transfer Function"');
end

%================= Error If both were not Selected ===================
if ss_pushbutton == get(handles.ss_radiob,'Min') && tf_pushbutton == get(handles.tf_radiob,'Min')
    errordlg('First Enter either "State Space" or "Transfer Function"');
end

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ==============
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    step(sys_tf),title('Step Response Without Controller'),grid on
end

%================= Implementing State Space Model ====================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ==============

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D');
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    step(sys_ss),title('Step Response Without Controller'),grid on
end

% --------------------------------------------------------------------
function IRWC_Callback(hObject, eventdata, handles)
% hObject    handle to IRWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Error If both were not Selected ===================
if ss_pushbutton == get(handles.ss_radiob,'Min') && tf_pushbutton == get(handles.tf_radiob,'Min')
    errordlg('First Enter either "State Space" or "Transfer Function"');
end

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ==============
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    Impulse(sys_tf),title('Impulse Response Without Controller'),grid on
end

%================= Implementing State Space Model ====================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ==============

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D')
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    impulse(sys_ss),title('Impulse Response Without Controller'),grid on
end
% --------------------------------------------------------------------
function PZMWC_Callback(hObject, eventdata, handles)
% hObject    handle to PZMWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ==============
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    pzmap(sys_tf),title('Pole Zero Map Without Controller'),grid on
end

%================= Implementing State Space Model ====================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ==============

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D')
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    pzmap(sys_ss),title('Pole Zero Map Without Controller'),grid on
end
%=====================================================================

% --------------------------------------------------------------------
function SRC_Callback(hObject, eventdata, handles)
% hObject    handle to SRC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,step(sys_fd_P),title('Step Response With Proportional Controller'),grid on
end

if value == 2
    global sys_fd_PI
    figure,step(sys_fd_PI),title('Step Response With Proportional Integral Controller'),grid on
end

if value == 3
    global sys_fd_PD
    figure,step(sys_fd_PD),title('Step Response With Proportional Derivative Controller'),grid on
end

if value == 4
    global sys_fd_PID
    figure,step(sys_fd_PID),title('Step Response With PID Controller'),grid on
end

% --------------------------------------------------------------------
function IRC_Callback(hObject, eventdata, handles)
% hObject    handle to IRC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,impulse(sys_fd_P),title('Impulse Response With Proportional Controller'),grid on
end

if value == 2
    global sys_fd_PI
    figure,impulse(sys_fd_PI),title('Impulse Response With Proportional Integral Controller'),grid on
end

if value == 3
    global sys_fd_PD
    figure,impulse(sys_fd_PD),title('Impulse Response With Proportional Derivative Controller'),grid on
end

if value == 4
    global sys_fd_PID
    figure,impulse(sys_fd_PID),title('Impulse Response With PID Controller'),grid on
end

% --------------------------------------------------------------------
function PZMC_Callback(hObject, eventdata, handles)
% hObject    handle to PZMC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,pzmap(sys_fd_P),title('Pole Zero Map With Proportional Controller'),grid on
end

if value == 2
    global sys_fd_PI
    figure,pzmap(sys_fd_PI),title('Pole Zero Map With Proportional Integral Controller'),grid on
end

if value == 3
    global sys_fd_PD
    figure,pzmap(sys_fd_PD),title('Pole Zero Map With Proportional Derivative Controller'),grid on
end

if value == 4
    global sys_fd_PID
    figure,pzmap(sys_fd_PID),title('Pole Zero Map With PID Controller'),grid on
end

% --------------------------------------------------------------------
function BPWC_Callback(hObject, eventdata, handles)
% hObject    handle to BPWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Error If both were not Selected ===================
if ss_pushbutton == get(handles.ss_radiob,'Min') && tf_pushbutton == get(handles.tf_radiob,'Min')
    errordlg('First Enter either "State Space" or "Transfer Function"');
end

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ================
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    bode(sys_tf),title('Bode Diagram Without Controller'),grid on
end

%================= Implementing State Space Model =========================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ================

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D');
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    bode(sys_ss),title('Bode Diagram Without Controller'),grid on
end
%=====================================================================
% --------------------------------------------------------------------
function NPWC_Callback(hObject, eventdata, handles)
% hObject    handle to NPWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Error If both were not Selected ===================
if ss_pushbutton == get(handles.ss_radiob,'Min') && tf_pushbutton == get(handles.tf_radiob,'Min')
    errordlg('First Enter either "State Space" or "Transfer Function"');
end

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ==============
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    nyquist(sys_tf),title('Nyquist Plot Without Controller'),grid on
end

%================= Implementing State Space Model ====================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ==============

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D');
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    nyquist(sys_ss),title('Nyquist Plot Without Controller'),grid on
end
%=====================================================================

% --------------------------------------------------------------------
function NCWC_Callback(hObject, eventdata, handles)
% hObject    handle to NCWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Error If both were not Selected ===================
if ss_pushbutton == get(handles.ss_radiob,'Min') && tf_pushbutton == get(handles.tf_radiob,'Min')
    errordlg('First Enter either "State Space" or "Transfer Function"');
end

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ==============
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    nichols(sys_tf),title('Nichols Chat Without Controller'),grid on
end

%================= Implementing State Space Model ====================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ==============

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D');
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    nichols(sys_ss),title('Nichols Chat Without Controller'),grid on
end
%=====================================================================

% --------------------------------------------------------------------
function BPC_Callback(hObject, eventdata, handles)
% hObject    handle to BPC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,bode(sys_fd_P),title('Bode Plot With Proportional Controller'),grid on
end

if value == 2
    global sys_fd_PI
    figure,bode(sys_fd_PI),title('Bode Plot With Proportional Integral Controller'),grid on
end

if value == 3
    global sys_fd_PD
    figure,bode(sys_fd_PD),title('Bode Plot With Proportional Derivative Controller'),grid on
end

if value == 4
    global sys_fd_PID
    figure,bode(sys_fd_PID),title('Bode Plot With PID Controller'),grid on
end

% --------------------------------------------------------------------
function NPC_Callback(hObject, eventdata, handles)
% hObject    handle to NPC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,nyquist(sys_fd_P),title('Nyquist Plot With Proportional Controller'),grid on
end

if value == 2
    global sys_fd_PI
    figure,nyquist(sys_fd_PI),title('Nyquist Plot With Proportional Integral Controller'),grid on
end

if value == 3
    global sys_fd_PD
    figure,nyquist(sys_fd_PD),title('Nyquist Plot With Proportional Derivative Controller'),grid on
end

if value == 4
    global sys_fd_PID
    figure,nyquist(sys_fd_PID),title('Nyquist Plot With PID Controller'),grid on
end

% --------------------------------------------------------------------
function NCC_Callback(hObject, eventdata, handles)
% hObject    handle to NCC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,nichols(sys_fd_P),title('Nichols ChartWith Proportional Controller'),grid on
end

if value == 2
    global sys_fd_PI
    figure,nichols(sys_fd_PI),title('Nichols Chart With Proportional Integral Controller'),grid on
end

if value == 3
    global sys_fd_PD
    figure,nichols(sys_fd_PD),title('Nichols Chart With Proportional Derivative Controller'),grid on
end

if value == 4
    global sys_fd_PID
    figure,nichols(sys_fd_PID),title('Nichols Chart With PID Controller'),grid on
end

% --------------------------------------------------------------------
function RLWC_Callback(hObject, eventdata, handles)
% hObject    handle to RLWC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%====================== Values from Radio Button ======================
tf_pushbutton = get (handles.tf_radiob,'Value');
ss_pushbutton = get(handles.ss_radiob,'Value');

%================= Error If both Options Selected ====================
if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Max')
    errordlg('Select One Option from "State Space" or "Transfer Function"');
end

%================= Error If both were not Selected ===================
if ss_pushbutton == get(handles.ss_radiob,'Min') && tf_pushbutton == get(handles.tf_radiob,'Min')
    errordlg('First Enter either "State Space" or "Transfer Function"');
end

%================= Implementing Transfer Function ====================
if tf_pushbutton == get(handles.tf_radiob,'Max') && ss_pushbutton == get(handles.ss_radiob,'Min')
    num= get (handles.num_edit,'string');
    num_dou = str2num(num);

    den= get (handles.den_edit,'string');
    den_dou = str2num(den);

    %================= Error If Either Num or Den is empty ==============
    if isempty(num_dou) || isempty(den_dou)
        errordlg('Please Enter "Transfer Function"');
    end

    sys_tf = tf(num_dou,den_dou);
    figure
    rlocus(sys_tf),title('Root Locus Without Controller')
end

%================= Implementing State Space Model ====================

if ss_pushbutton == get(handles.ss_radiob,'Max') && tf_pushbutton == get(handles.tf_radiob,'Min')

    a= get(handles.A_edit,'String');
    a_dou = str2num(a);

    b= get(handles.B_edit,'String');
    b_dou = str2num(b);

    c= get(handles.C_edit,'String');
    c_dou = str2num(c);

    d= get(handles.D_edit,'String');
    d_dou = str2num(d);

    %================= Error If Either A,B,C or D is empty ==============

    if isempty(a_dou) || isempty(b_dou) || isempty(c_dou) || isempty(d_dou)
        errordlg('Please enter all A,B,C and D');
    end

    sys_ss = ss(a_dou,b_dou,c_dou,d_dou);
    figure
    rlocus(sys_ss),title('Root Locus Without Controller'),grid on
end

% --------------------------------------------------------------------
function RLC_Callback(hObject, eventdata, handles)
% hObject    handle to RLC (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

value = get(handles.PID_popupmenu,'Value');

if value == 1
    global sys_fd_P
    figure,rlocus(sys_fd_P),title('Root Locus With Proportional Controller')
end

if value == 2
    global sys_fd_PI
    figure,rlocus(sys_fd_PI),title('Root Locus With Proportional Integral Controller')
end

if value == 3
    global sys_fd_PD
    figure,rlocus(sys_fd_PD),title('Root Locus With Proportional Derivative Controller')
end

if value == 4
    global sys_fd_PID
    figure,rlocus(sys_fd_PID),title('Root Locus With PID Controller')
end

% --------------------------------------------------------------------

function Untitled_2_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_3_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_4_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------

% --- Executes on selection change in TAFR_popupmenu.
function TAFR_popupmenu_Callback(hObject, eventdata, handles)
% hObject    handle to TAFR_popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns TAFR_popupmenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from TAFR_popupmenu

% --- Executes during object creation, after setting all properties.
function TAFR_popupmenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to TAFR_popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in TAFR_pushb.
function TAFR_pushb_Callback(hObject, eventdata, handles)
% hObject    handle to TAFR_pushb (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% =============== Values from Pop Up menus (Controller and Responses)======

value_response = get(handles.TAFR_popupmenu,'Value');
value = get(handles.PID_popupmenu,'Value');

% =========================================================================

% ========================Proportional Controller =========================
if value ==1

    % ==================== Time Response ==================================

    if value_response == 1
        set(handles.MOlabel_text,'String','Maximum Overshoot >');
        set(handles.PTlabel_text,'String','Peak Time >');
        set(handles.STlabel_text,'String','Settling Time >');

        global sys_fd_P
        order_sys = order(sys_fd_P);
        order_sys_str = num2str(order_sys);
        set(handles.Order_text,'String',order_sys_str);

        [num,den] = tfdata(sys_fd_P);
        [r c] = size(den{1,1});

        if c ~= 3
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end

        roots_sys = roots(den{1,1});
        roots_sys_str = num2str(roots_sys);

        set(handles.roots_text,'String',roots_sys_str);
        den_sp = den{1,1}(1,3);
        wn = sqrt(den_sp);
        den_sp2 = den{1,1}(1,2);
        b = den_sp2/(2*wn);

        if b<=1 && num{1,1}(1,end-1) == 0 && c==3
            MO = exp(-(b*pi)/(sqrt(1-b^2)))*100;
            MO_str = num2str(MO);
            MO_real = strcat(MO_str,' %');
            set(handles.MO_text,'String',MO_real);

            PT = pi/(wn*(sqrt(1-b^2)));
            PT_str = num2str(PT);
            PT_real = strcat(PT_str,' Sec');
            set(handles.PT_text,'String',PT_real);

            ST = 4/(b*wn);
            ST_str = num2str(ST);
            ST_real = strcat(ST_str,' Sec');
            set(handles.ST_text,'String',ST_real);

        else
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end
    end
    % ====================== Frequncy Response ============================

    if value_response ==2
        set(handles.MOlabel_text,'String','Gain Margin >');
        set(handles.PTlabel_text,'String','Phase Margin >');
        set(handles.STlabel_text,'String','Bandwidth >');

        global sys_fd_P
        order_sys = order(sys_fd_P);
        order_sys_str = num2str(order_sys);
        set(handles.Order_text,'String',order_sys_str);

        band = bandwidth(sys_fd_P);
        band_str = num2str(band);
        set(handles.ST_text,'String',band_str);

        [gm pm wcg wcp]=margin (sys_fd_P);
        gm_str = num2str(gm);
        set(handles.MO_text,'String',gm_str);

        pm_str = num2str(pm);
        set(handles.PT_text,'String',pm_str);
    end
end

% ================Proportional Integral Controller ========================

if value ==2

    % ==================== Time Response ==================================

    if value_response == 1
        set(handles.MOlabel_text,'String','Maximum Overshoot >');
        set(handles.PTlabel_text,'String','Peak Time >');
        set(handles.STlabel_text,'String','Settling Time >');

        global sys_fd_PI
        order_sys = order(sys_fd_PI);
        order_sys_str = num2str(order_sys);
        set(handles.Order_text,'String',order_sys_str);

        [num,den] = tfdata(sys_fd_PI);
        [r c] = size(den{1,1});

        if c ~= 3
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end

        roots_sys = roots(den{1,1});
        roots_sys_str = num2str(roots_sys);

        set(handles.roots_text,'String',roots_sys_str);
        den_sp = den{1,1}(1,3);
        wn = sqrt(den_sp);
        den_sp2 = den{1,1}(1,2);
        b = den_sp2/(2*wn);

        if b<=1 && num{1,1}(1,end-1) == 0 && c==3
            MO = exp(-(b*pi)/(sqrt(1-b^2)))*100;
            MO_str = num2str(MO);
            MO_real = strcat(MO_str,' %');
            set(handles.MO_text,'String',MO_real);

            PT = pi/(wn*(sqrt(1-b^2)));
            PT_str = num2str(PT);
            PT_real = strcat(PT_str,' Sec');
            set(handles.PT_text,'String',PT_real);

            ST = 4/(b*wn);
            ST_str = num2str(ST);
            ST_real = strcat(ST_str,' Sec');
            set(handles.ST_text,'String',ST_real);

        else
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end
    end

    % ==================== Frequency Response =============================

    if value_response ==2
        set(handles.MOlabel_text,'String','Gain Margin >');
        set(handles.PTlabel_text,'String','Phase Margin >');
        set(handles.STlabel_text,'String','Bandwidth >');

        global sys_fd_PI
        order_sys = order(sys_fd_PI);
        order_sys_str = num2str(order_sys);
        set(handles.Order_text,'String',order_sys_str);

        band = bandwidth(sys_fd_PI);
        band_str = num2str(band);
        set(handles.ST_text,'String',band_str);

        [gm pm wcg wcp]=margin (sys_fd_PI);
        gm_str = num2str(gm);
        set(handles.MO_text,'String',gm_str);

        pm_str = num2str(pm);
        set(handles.PT_text,'String',pm_str);
    end
end

% ================Proportional Integral Controller =========================

if value ==3

     % ==================== Time Response ==================================

     if value_response == 1
        set(handles.MOlabel_text,'String','Maximum Overshoot >');
        set(handles.PTlabel_text,'String','Peak Time >');
        set(handles.STlabel_text,'String','Settling Time >');

        global sys_fd_PD
        order_sys = order(sys_fd_PD);
        order_sys_str = num2str(order_sys);
        set(handles.Order_text,'String',order_sys_str);

        [num,den] = tfdata(sys_fd_PD);
        [r c] = size(den{1,1});

        if c ~= 3
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end

        roots_sys = roots(den{1,1});
        roots_sys_str = num2str(roots_sys);

        set(handles.roots_text,'String',roots_sys_str);
        den_sp = den{1,1}(1,3);
        wn = sqrt(den_sp);
        den_sp2 = den{1,1}(1,2);
        b = den_sp2/(2*wn);

        if b<=1 && num{1,1}(1,end-1) == 0 && c==3
            MO = exp(-(b*pi)/(sqrt(1-b^2)))*100;
            MO_str = num2str(MO);
            MO_real = strcat(MO_str,' %');
            set(handles.MO_text,'String',MO_real);

            PT = pi/(wn*(sqrt(1-b^2)));
            PT_str = num2str(PT);
            PT_real = strcat(PT_str,' Sec');
            set(handles.PT_text,'String',PT_real);

            ST = 4/(b*wn);
            ST_str = num2str(ST);
            ST_real = strcat(ST_str,' Sec');
            set(handles.ST_text,'String',ST_real);

        else
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end
     end

     % ==================== Frequency Response ============================

     if value_response ==2
         set(handles.MOlabel_text,'String','Gain Margin >');
         set(handles.PTlabel_text,'String','Phase Margin >');
         set(handles.STlabel_text,'String','Bandwidth >');

         global sys_fd_PD
         order_sys = order(sys_fd_PD);
         order_sys_str = num2str(order_sys);
         set(handles.Order_text,'String',order_sys_str);

         band = bandwidth(sys_fd_PD);
         band_str = num2str(band);
         set(handles.ST_text,'String',band_str);

         [gm pm wcg wcp]=margin (sys_fd_PD);
         gm_str = num2str(gm);
         set(handles.MO_text,'String',gm_str);

         pm_str = num2str(pm);
         set(handles.PT_text,'String',pm_str);

     end
end

% ================Proportional Integral Derivative Controller =============

if value ==4

     % ==================== Time Response =================================

     if value_response == 1
        set(handles.MOlabel_text,'String','Maximum Overshoot >');
        set(handles.PTlabel_text,'String','Peak Time >');
        set(handles.STlabel_text,'String','Settling Time >');

        global sys_fd_PID
        order_sys = order(sys_fd_PID);
        order_sys_str = num2str(order_sys);
        set(handles.Order_text,'String',order_sys_str);

        [num,den] = tfdata(sys_fd_PID);
        [r c] = size(den{1,1});

        if c ~= 3
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end

        roots_sys = roots(den{1,1});
        roots_sys_str = num2str(roots_sys);

        set(handles.roots_text,'String',roots_sys_str);
        den_sp = den{1,1}(1,3);
        wn = sqrt(den_sp);
        den_sp2 = den{1,1}(1,2);
        b = den_sp2/(2*wn);

        if b<=1 && num{1,1}(1,end-1) == 0 && c==3
            MO = exp(-(b*pi)/(sqrt(1-b^2)))*100;
            MO_str = num2str(MO);
            MO_real = strcat(MO_str,' %');
            set(handles.MO_text,'String',MO_real);

            PT = pi/(wn*(sqrt(1-b^2)));
            PT_str = num2str(PT);
            PT_real = strcat(PT_str,' Sec');
            set(handles.PT_text,'String',PT_real);

            ST = 4/(b*wn);
            ST_str = num2str(ST);
            ST_real = strcat(ST_str,' Sec');
            set(handles.ST_text,'String',ST_real);

        else
            set(handles.MO_text,'String','Approximation Not Valid');
            set(handles.PT_text,'String','Approximation Not Valid');
            set(handles.ST_text,'String','Approximation Not Valid');
        end
     end

     % ==================== Frequency Response ============================

     if value_response == 2
         set(handles.MOlabel_text,'String','Gain Margin >');
         set(handles.PTlabel_text,'String','Phase Margin >');
         set(handles.STlabel_text,'String','Bandwidth >');

         global sys_fd_PID
         order_sys = order(sys_fd_PID);
         order_sys_str = num2str(order_sys);
         set(handles.Order_text,'String',order_sys_str);

         band = bandwidth(sys_fd_PID);
         band_str = num2str(band);
         set(handles.ST_text,'String',band_str);

         [gm pm wcg wcp]=margin (sys_fd_PID);
         gm_str = num2str(gm);
         set(handles.MO_text,'String',gm_str);

         pm_str = num2str(pm);
         set(handles.PT_text,'String',pm_str);
     end
end