Code covered by the BSD License  

Highlights from
Circuit Analysis Toolbox

from Circuit Analysis Toolbox by John McDermid
The circuit analysis toolbox allows you to perform an AC analysis of a circuit.

tableau
classdef tableau
    % Creates an object for the simplified tableau format.
    % The tableau matrix is system of equations contrived from two matrix
    % relations. The first is the relationship between elemental values
    % described by the simple equation ki*Ib-kv*Vb=s. Ib is the current
    % through component and Vb is the voltage across the component. The
    % second is A*IB=0 where A is the reduced incidence matrix and IB is a
    % vector of branch currents. A*IB=0 is a matrix statement of Kirkoff's
    % current law for all nodes.
    %
    % Values of ki, kv, and s are picked to match each component type.
    %       Resistor
    %           ki=R,  kv=1,  s=0,   making R=Vb/Ib
    %       Inductor
    %           ki=j*2*pi*f*L,  kv=1,  s=0,  making Vb/Ib=j*2*pi*f*L=j*XL
    %       Capacitor
    %           ki=1,  kv=j*2*pi*f*C,  s=0  making Vb/Ib=1/(j*2*pi*f*C)=j*XC
    %       Independent Current Source
    %           ki=1,  kv=0,  s=I   making Ib=I
    %       Independent Voltage Source
    %           ki=0,  kv=1,  s=V   making -Vb=V (note that the current in
    %                                             the voltage source flows in the
    %                                             opposite direction to the
    %                                             voltage accross it and
    %                                             hence the minus sign)
    %       Current Dependent Current Source
    %           ki=[1 -Beta],  kv=[0 0],  s=0, making [1 -Beta)*[Ib
    %           Ib_cntrl]'=0
    %                                       or Ib=Beta*Ib_cntrl
    %       Voltage Dependent Current Source
    %           ki=1,  kv=gm,  s=0  making Ib-gm*Vb_cntrl
    %                                     or Ib=gm*Vb_cntrl
    %       Voltage Dependent Voltage Source
    %           ki=[0 0],  kv=[1 Av],  s=0,  -[1 -Av]*[Vb VB_cntrl]'=0
    %                                            or -Vb=Av*Vb_cntrl
    %                                             (note that the current in
    %                                             the voltage source flows in the
    %                                             opposite direction to the
    %                                             voltage accross it and
    %                                             hence the minus sign)
    %
    % Kv and Ki (constructed from kv and ki) are essentially diagonal
    % matricies except for the dependent sources.
    %
    % The tableau matrix is constructed in the following way:
    %       tableau=[Kv -Ki*A';A Z] where Z is a matrix of zeros and of proper
    %                                dimension to make tableau a square matrix
    % The S column vector is constructed from the individual s and 0
    %
    % PROPERTIES:
    %         ckt                Structure of circuit properties
    %         typesAvailable     A list of the types available in this circuit
    %         typeIndx           An index to the list of types indexed by type name
    %         desigIndx          An index to the designators in A, Ki, and Kv by designator name
    %         nodeIndx           An index to nodes in Aa by designator name
    %         VnNodeIndx         An index to nodes in A by designator name
    %         TmIndx             An index to all elements in the tableau matrix
    %                              including designators and nodes by name
    %         desigList          The list of designators in matrix order
    %         nodeList           The list of nodes in matrix order
    %         branchIndx         Index to all branch currents in tableau matrix
    %         nodeVoltIndx       Index to all node voltages in the tableau matrix
    %         Aa                 The complete incidence matrix
    %         accInAa            A logical vector indicating which nodes in Aa have access
    %         accInA             A logical vector indicating which nodes in A have access
    %         A                  The reduced incidence matrix with the reference node removed (Sparse)
    %         Ki                 A sparse matrix holding impedance values
    %         Kv                 A sparse matrix holding conductance values
    %         Z                  A sparse array of zeros
    %         S                  A sparse vector of source voltages and currents
    %         deltaKi            A sparse matrix containing a "1" at each value location in Ki
    %         deltaKv            A sparse matrix containing a "1" at each value location in Kv
    %         deltaS             A sparse vector containing a "1" at each value location in S
    %
    %
    % EXAMPLE:
    %        tableauObject = tableau(circuitObject)
    
    % AUTHOR:
    %       John McDermid
    % CREATED:
    %       October 19, 2009
    
    properties (SetAccess = protected)
    % tableau Properties:  
        ckt               % Structure of circuit properties
        typesAvailable    % A list of the types availble in this circuit
        typeIndx          % An index to the list of types indexed by type name
        desigIndx         % An index to the designators in A, Ki, and Kv by designator name
        nodeIndx          % An index to nodes in Aa by designator name
        VnNodeIndx        % An index to nodes in A by designator name
        TmIndx            % An index to all elements in the tableau matrix 
                          %   including designators and nodes by name
        desigList         % The list of designators in matrix order
        nodeList          % The list of nodes in matrix order
        branchIndx        % Index to all branch currents in tableau matrix
        nodeVoltIndx      % Index to all node voltages in the tableau matrix
        accInAa           % A logical vector indicating which nodes in Aa have access
        accInA            % A logical vector indicating which nodes in A have access
        Aa                % The complete incidence matrix (Sparse)
        A                 % The reduced incidence matrix with the reference node removed (Sparse)
        Ki                % A sparse matrix holding impedance values
        Kv                % A sparse matrix holding conductance values
        Z                 % A sparse array of zeros
        S                 % A sparse vector of source voltages and currents
        deltaKi           % A sparse matrix containing a "1" at each value location in Ki
        deltaKv           % A sparse matrix containing a "1" at each value location in Kv
        deltaS            % A sparse vector containing a "1" at each value location in S
    end
    
    methods 
        function obj=tableau(cirObj)
            % Constructor for a tableau object.
            % Requires a circuit object as input.
            
            % Find the index of types in the cirDesc structure so that
            % components can be grouped according to type and place them in
            % a container so that the indexof types can be retrieved by
            % type name
            uniqueTypes=unique({cirObj.cirDesc.type});
            M=numel(uniqueTypes);            % Number of unique component types in the circuit
            typeIndexCirObj{M}=[];           % Pre-allocate for speed
            N=numel({cirObj.cirDesc.type});  % Number of components in the circuit
            for J=1:M                        % Assemble the component index for each unique type
                row=[];
                for I=1:N
                    if cirObj.cirDesc(I).type==uniqueTypes{J}
                        row=[row I];   %#ok<AGROW>
                    end
                end
                typeIndexCirObj{J}=row;
            end
            types=containers.Map(uniqueTypes,typeIndexCirObj); 
            % Build the ckt property (a structure) from the unique types
            for I=1:M
                typeName=char(uniqueTypes(I));
                obj.ckt.(typeName).designator={cirObj.cirDesc(types(typeName)).designator};
                obj.ckt.(typeName).fromNode={cirObj.cirDesc(types(typeName)).fromNode};
                obj.ckt.(typeName).toNode={cirObj.cirDesc(types(typeName)).toNode};
                obj.ckt.(typeName).value=[cirObj.cirDesc(types(typeName)).value];
                obj.ckt.(typeName).plusTol=[cirObj.cirDesc(types(typeName)).plusTol];
                obj.ckt.(typeName).negTol=[cirObj.cirDesc(types(typeName)).negTol];
                obj.ckt.(typeName).cntrlBranch={cirObj.cirDesc(types(typeName)).cntrlBranch};
            end
            % Find the order of unique types that matches the order of the type names
            %  available
            typeNamesAvailable={'V' 'I' 'R' 'L' 'C' 'IDI' 'VDI' 'VDV'};
            currentOrder(M)=0;
            % Created an order for the unique types frm the order in type
            % names available
            for I=1:M
                currentOrder(I)=find(strcmp(uniqueTypes(I),typeNamesAvailable));
            end
            [~,typeOrder]=sort(currentOrder);
            % Make a container for designators that indexes the branch
            % dimension in the A matrix. Construct it in the order defined
            % by "obj.typesAvailable".
             obj.typesAvailable=uniqueTypes(typeOrder);
            % Makes a container for designators and a container for types
            % and store them in the object
            typeIndex=containers.Map;
            desigIndex=containers.Map;
            count=0;
            for I=obj.typesAvailable
                indexMap=[];
                for J=1:numel(obj.ckt.(I{:}).designator)
                    count=count+1;
                    indexMap=[indexMap count]; %#ok<AGROW>
                    keyName=char(obj.ckt.(I{:}).designator(J));
                    desig_List(count)={keyName}; %#ok<AGROW>
                    desigIndex(keyName)=count;
                end
                typeIndex(char(I{:}))=indexMap;
            end
            obj.typeIndx=typeIndex;
            obj.desigIndx=desigIndex;
            obj.desigList=desig_List;
            % Make a container for the index of nodes in the Aa matrix
            nNodes=numel(cirObj.nodeList);
            nodeIndex=containers.Map;
            for I=1:nNodes
                node=cirObj.nodeList(I).list;
                nodeIndex(node)=I;
            end
            obj.nodeIndx=nodeIndex;
            % Make a logical vector for accessible nodes in Aa
            obj.accInAa=[cirObj.nodeList.accessibleFlag];
            % Make a logical vector for for the accessible nodes in A
            obj.accInA=obj.accInAa(~[cirObj.nodeList.refNodeIndex]);
            % Create the incidence matrix. The unreduced incidence
            % matrix is created first and then the reference (ground) node is
            % removed to generate the reduced incidence matrix
            A_unReduced(obj.nodeIndx.Count,obj.desigIndx.Count)=0;
            for I=obj.typesAvailable
                typeName=char(I);
                N=numel(obj.ckt.(typeName).fromNode);
                %Define the "FROM" node value
                for J=1:N
                    indexName=char(obj.ckt.(typeName).fromNode(J));
                    rowIndexA=nodeIndex(indexName);
                    indexDesig=char(obj.ckt.(typeName).designator(J));
                    colIndexA=obj.desigIndx(indexDesig);
                    A_unReduced(rowIndexA,colIndexA)=1;
                end
            end
            for I=obj.typesAvailable
                typeName=char(I);
                N=numel(obj.ckt.(typeName).toNode);
                % Define the "TO" node value
                for J=1:N
                    indexName=char(obj.ckt.(typeName).toNode(J));
                    rowIndexA=nodeIndex(indexName);
                    indexDesig=char(obj.ckt.(typeName).designator(J));
                    colIndexA=obj.desigIndx(indexDesig);
                    A_unReduced(rowIndexA,colIndexA)=-1;
                end
            end
            obj.Aa=sparse(A_unReduced);
            % Create the reduced A matrix by removing the row defined by
            % the reference node
            A_unReduced(nodeIndex(cirObj.refNode),:)=[];
            obj.A=sparse(A_unReduced);
            % Create the index for the reduced incidence matrix
            VnNodeIndex=containers.Map; 
            nCount=0;
            for I=1:nNodes
                node=cirObj.nodeList(I).list;
                if logical(cirObj.nodeList(I).refNodeIndex)==0
                    nCount=nCount+1;
                    VnNodeIndex(node)=nCount;
                    node_List(nCount)={node};  %#ok<AGROW>
                end
            end
            obj.VnNodeIndx=VnNodeIndex;
            obj.nodeList=node_List;
            % Construct Ki, Kv, S, and Z matrices
            Ki_(obj.desigIndx.Count,obj.desigIndx.Count)=0;
            Kv_(obj.desigIndx.Count,obj.desigIndx.Count)=0;
            deltaKi_=Ki_;
            deltaKv_=Kv_;
            S_(double(obj.desigIndx.Count)+double(obj.VnNodeIndx.Count))=0;
            S_=S_';         % Make "S" a column Vector
            deltaS_=S_;
            for I=obj.typesAvailable
                typeName=char(I);
                Indx=typeIndex(typeName);
                switch typeName
                    case 'V'
                        S_(Indx)=obj.ckt.V.value;
                        deltaS_(Indx)=ones(numel(Indx),1);
                        Kv_(Indx,Indx)=eye(numel(Indx)); 
                    case 'I'
                        S_(Indx)=obj.ckt.I.value; 
                        deltaS_(Indx)=ones(numel(Indx),1);
                        Ki_(Indx)=eye(numel(Indx)); 
                    case 'R'
                        Ki_(Indx,Indx)=diag(obj.ckt.R.value);
                        deltaKi_(Indx,Indx)=eye(numel(Indx));
                        Kv_(Indx,Indx)=eye(numel(Indx)); 
                    case 'L'
                        Ki_(Indx,Indx)=diag(1j*obj.ckt.L.value);
                        deltaKi_(Indx,Indx)=eye(numel(Indx));
                        Kv_(Indx,Indx)=eye(numel(Indx)); 
                    case 'C'
                        Ki_(Indx,Indx)=eye(numel(Indx)); 
                        Kv_(Indx,Indx)=diag(1j*obj.ckt.C.value);
                        deltaKv_(Indx,Indx)=eye(numel(Indx));
                    case 'IDI'
                        Ki_(Indx,Indx)=eye(numel(Indx));
                        for J=1:numel(Indx)
                            dependentBranch=obj.desigIndx(char(obj.ckt.IDI.cntrlBranch(J)));
                            independentBranch=Indx(J);
                            Ki_(independentBranch,dependentBranch)=-obj.ckt.IDI.value(J);
                            deltaKi_(independentBranch,dependentBranch)=1;
                        end
                    case 'VDI'
                        Ki_(Indx,Indx)=eye(numel(Indx));
                        for J=1:numel(Indx)
                            dependentBranch=obj.desigIndx(char(obj.ckt.VDI.cntrlBranch(J)));
                            independentBranch=Indx(J);
                            Kv_(independentBranch,dependentBranch)=obj.ckt.VDI.value(J);
                            deltaKv_(independentBranch,dependentBranch)=1;
                        end
                    case 'VDV'
                        Kv_(Indx,Indx)=eye(numel(Indx));
                        for J=1:numel(Indx)
                            dependentBranch=obj.desigIndx(char(obj.ckt.VDV.cntrlBranch(J)));
                            independentBranch=Indx(J);
                            Kv_(independentBranch,dependentBranch)=obj.ckt.VDV.value(J);
                            deltaKv_(independentBranch,dependentBranch)=1;
                        end
                end            
            end
            obj.S=sparse(S_);  
            obj.Ki=sparse(Ki_);
            obj.Kv=sparse(Kv_);
            obj.deltaS=sparse(deltaS_);
            obj.deltaKi=sparse(deltaKi_);
            obj.deltaKv=sparse(deltaKv_);
            obj.Z=sparse(zeros(obj.VnNodeIndx.Count,obj.VnNodeIndx.Count)); 
            % Create an index for the completed tableau_matrix
            Tindex=containers.Map;
            for I=keys(obj.desigIndx)
                charI=char(I);
                Tindex(charI)=obj.desigIndx(charI);
                obj.branchIndx=1:double(obj.desigIndx.Count);
            end
            count=0;
            for I=keys(obj.VnNodeIndx)   
                count=count+1;
                charI=char(I);
                val=obj.VnNodeIndx(charI)+double(obj.desigIndx.Count);
                obj.nodeVoltIndx(count)=val;
                Tindex(charI)=val;
            end
            obj.TmIndx=Tindex;
        end
        
        function matrixT=tableau_matrix(T,f,varyKi,varyKv)
            % Creates a tableau matrix from a tableau object at the
            % specified frequency. Matrix T is a tableau object. Parameter
            % f is the frequency used to create the tableau matrix.
            % Variable "varyKi" and "varyKv" are parameters supplied by the
            % function "monteCarlo" and provide a random permutation of
            % parameters within the defined tolerance level for each
            % component in Kv and Ki. Parameters T and f must be supplied.
            % If varyKi and varyKv are not supplied, the tableau matrix is
            % created from the nominal values of the components. If they
            % are supplied, the tableau matrix will be created from a
            % random distribution of component values within the stated
            % component tolerance.
            % Example:
            %      matrixT=tableau_matrix(T,f)
            %
            %       or optionally
            %
            %      matrixT=tableau_matrix(T,f,varyKi,varyKv)
            %
            
            % Find all capacitors and set the appropriate parts of Kv to
            % the desired frequency.
            if any(strcmp('C',T.typesAvailable))
                Indx=T.typeIndx('C');
                T.Kv(Indx,Indx)=diag(1j*2*pi*f*T.ckt.C.value);
            end
            % Find all the inductors and set the appropriate parts of Ki to
            % the desired frequency.
            if any(strcmp('L',T.typesAvailable))
                Indx=T.typeIndx('L');
                T.Ki(Indx,Indx)=diag(1j*2*pi*f*T.ckt.L.value);
            end
            % If varyKi and varyKv are not supplied, use the nominal values
            % of the components.
            if nargin==2
                matrixT=[T.Ki -T.Kv*T.A';T.A T.Z];
            else
                matrixT=[T.Ki.*varyKi -T.Kv.*varyKv*T.A';T.A T.Z];
            end
        end
        
        function data=bodeplot(tableauObj,name,startFreq,stopFreq,nPoints)

            % BODEPLOT creates a bode plot for the node or branch specified
            % by name. The parameter "name" is a string specifying a branch
            % or node name. The frequency range is specified by the start
            % and stop frequencies. The number of points in the frequency
            % range may optionally be specified. (Default is 50) The
            % resulting plot is not docked in the desk top.
            %
            % Example:
            % bodeplot(tableauObj,name,startFreq,stopFreq)
            %
            %  or optionally
            %
            % bodeplot(tableauObj,name,startFreq,stopFreq,nPoints)
            %
            
            % Check for optional arguments, define the frequency range, and
            % predefine the output array IVmat for speed.
            if nargin==4
                nPoints=50;
            elseif nargin~=5
                error('bodeplot requires 4 or 5 parameters!')
            end
            freq=logspace(log10(startFreq),log10(stopFreq),nPoints);
            IVmat(numel(tableauObj.S),nPoints)=0;
            % Invert the matrix for each of the specified frequencies
            for n=1:nPoints
                T=tableau_matrix(tableauObj,freq(n));
                IVmat(:,n)=T\tableauObj.S;
            end
            % Create the bode plot as two subplots and label each plot
            indx=tableauObj.TmIndx(name);
            magnitude=abs(IVmat(indx,:));
            phase=180/pi*(angle(IVmat(indx,:)));
            % Make numeric data available in the workspace before plotting
            % if requested
            if nargout==1
                data=[freq;magnitude;phase]';
            end
            h=figure;
            set(h,'WindowStyle','normal')
            subplot(2,1,1);loglog(freq,magnitude);
                xlabel('Frequency (Hz)');
                if ismember(tableauObj.TmIndx(name),tableauObj.branchIndx)
                    ylabel(['Magnitude of Current ' name ' (Amps)']);
                    title(['Magnitude of the Current Through ' name]);
                else
                    ylabel(['Voltage Magnitude At ' name ' (Volts)']);
                    title(['Magnitude of the Voltage At Node ' name]);
                end
            subplot(2,1,2);semilogx(freq,phase);
                xlabel('Frequency (Hz)');
                if ismember(tableauObj.TmIndx(name),tableauObj.branchIndx)
                    ylabel(['Phase of Current ' name ' (Degrees)']);
                    title(['Phase of the Current Through ' name]);
                else
                    ylabel(['Voltage Phase At ' name ' (Degrees)']);
                    title(['Phase of the Voltage At Node ' name]);
                end 
            set(h,...
                'Position',[560 162 560 786],...
                'Name',['Bode Plot At ' name],...
                'NumberTitle','off');
        end
        
        function mCarloBodeplot(tableauObj,name,startFreq,stopFreq,nCkts,nPoints)
            % MCARLOBODEPLOT creates a bode plot for the node or branch
            % specified by string "name", one line on the graph for each
            % circuit simulated. The frequency range is specified by the
            % start and stop frequencies. The number of circuits simulated
            % is specified by "nCkts". The optional variable "nPoints"
            % (Default is 50) specifies the number of points in the
            % frequency sweep. The resulting plot is not docked in the desk
            % top.
            % 
            % Example:
            %      mCarloBodeplot(tableauObj,name,startFreq,stopFreq,nCkts)
            %
            %        or optionally
            %
            %      mCarloBodeplot(tableauObj,name,startFreq,stopFreq,nCkts,nPoints)
            %
            
            % Check for optional arguments, define the frequency range, and
            % initialize IVmat, magnitude, and phase for speed.
            if nargin==5
                nPoints=50;
            elseif nargin~=6
                error('mCarloBodeplot requires 5 or 6 parameters!')
            end
            freq=logspace(log10(startFreq),log10(stopFreq),nPoints);
            IVmat(numel(tableauObj.S),nPoints)=0;
            magnitude(nCkts,nPoints)=0;
            phase(nCkts,nPoints)=0;
            % Create all the circuit simulations
            indx=tableauObj.TmIndx(name);
            for m=1:nCkts
                % Create the matricies to vary component parameters
                [varyKi,varyKv,varyS]=monteCarlo(tableauObj);
                % Invert the matrix for each of the specified frequencies
                for n=1:nPoints
                    T=tableau_matrix(tableauObj,freq(n),varyKi,varyKv);
                    IVmat(:,n)=T\(tableauObj.S.*varyS);
                end
                % Create the bode plot as two subplots and label each plot
                magnitude(m,:)=abs(IVmat(indx,:));
                phase(m,:)=180/pi*(angle(IVmat(indx,:)));
            end
            % Create the gain and phase plots for all circuits
            h=figure;
            set(h,'WindowStyle','normal')
            subplot(2,1,1);loglog(freq,magnitude','color','blue');
                xlabel('Frequency (Hz)');
                if ismember(tableauObj.TmIndx(name),tableauObj.branchIndx)
                    ylabel(['Magnitude of Current ' name ' (Amps)']);
                    title(['Magnitude of the Current Through ' name]);
                else
                    ylabel(['Voltage Magnitude At ' name ' (Volts)']);
                    title(['Magnitude of the Voltage At Node ' name]);
                end
            subplot(2,1,2);semilogx(freq,phase','color','blue');
                xlabel('Frequency (Hz)');
                if ismember(tableauObj.TmIndx(name),tableauObj.branchIndx)
                    ylabel(['Phase of Current ' name ' (Degrees)']);
                    title(['Phase of the Current Through ' name]);
                else
                    ylabel(['Voltage Phase At ' name ' (Degrees)']);
                    title(['Phase of the Voltage At Node ' name]);
                end 
            set(h,...
                'Position',[560 162 560 786],...
                'Name',['Bode Plot At ' name],...
                'NumberTitle','off');
        end
        
        function [varyKi,varyKv,varyS]=monteCarlo(obj,Kinit,Sinit) 
            % Constructs a matrix that will randomly vary S, Ki, and Kv
            % within the defined tolerance range of pTol and nTol. The
            % input is a tableau object. The random variation is a uniform
            % distribution over the stated range. To vary S (which is a
            % sparse matrix) it is multiplied by 'varyS' which is also a
            % sparse matrix. The multiplication is S.*varyS which is an
            % element-by-element multiplication. Since both arrays are the
            % same size and sparse, only the non-zero elements are actually
            % multiplied. Each time the monteCarlo function is called a new
            % set of multiplying values is created.
            %
            % Kinit and Sinit are optional parameters and should be
            % initialized as sparse matricies (value of zero) to prevent
            % the overhead of initializing the matricies every time.
            %
            % Example:
            %      [varyKi,varyKv,varyS]=monteCarlo(obj)
            %
            %        or optionally
            %
            %      [varyKi,varyKv,varyS]=monteCarlo(obj,Kinit,Sinit)
            %
            
            % Initialize S_, Ki_, and Kv_ for performance
            switch nargin
                case 1
                    Ki_(obj.desigIndx.Count,obj.desigIndx.Count)=0;
                    Kv_(obj.desigIndx.Count,obj.desigIndx.Count)=0;
                    S_(double(obj.desigIndx.Count)+double(obj.VnNodeIndx.Count),1)=0;
                case 3
                    Ki_=Kinit;
                    Kv_=Kinit;
                    S_=Sinit;
                otherwise
                    error('monteCarlo requires 1 or 3 arguments')
            end
            for I=obj.typesAvailable
                typeName=char(I);
                Indx=obj.typeIndx(typeName);
                switch typeName
                    case 'V'
                        multVal=range(obj.ckt.V.negTol,obj.ckt.V.plusTol);
                        S_(Indx)=multVal;
                        Kv_(Indx,Indx)=eye(numel(Indx)); 
                    case 'I'
                        multVal=range(obj.ckt.I.negTol,obj.ckt.I.plusTol);
                        S_(Indx)=multVal;
                        Ki_(Indx)=eye(numel(Indx)); 
                    case 'R'
                        multVal=range(obj.ckt.R.negTol,obj.ckt.R.plusTol);
                        Ki_(Indx,Indx)=diag(multVal);
                        Kv_(Indx,Indx)=eye(numel(Indx)); 
                    case 'L'
                        multVal=range(obj.ckt.L.negTol,obj.ckt.L.plusTol);
                        Ki_(Indx,Indx)=diag(multVal);
                        Kv_(Indx,Indx)=eye(numel(Indx)); 
                    case 'C'
                        Ki_(Indx,Indx)=eye(numel(Indx)); 
                        multVal=range(obj.ckt.C.negTol,obj.ckt.C.plusTol);
                        Kv_(Indx,Indx)=diag(multVal);
                    case 'IDI'
                        Ki_(Indx,Indx)=eye(numel(Indx));
                        for J=1:numel(Indx)
                            dependentBranch=obj.desigIndx(char(obj.ckt.IDI.cntrlBranch(J)));
                            independentBranch=Indx(J);
                            multVal=range(obj.ckt.IDI.negTol,obj.ckt.IDI.plusTol);
                            Ki_(independentBranch,dependentBranch)=multVal(J);
                        end
                    case 'VDI'
                        Ki_(Indx,Indx)=eye(numel(Indx));
                        for J=1:numel(Indx)
                            dependentBranch=obj.desigIndx(char(obj.ckt.VDI.cntrlBranch(J)));
                            independentBranch=Indx(J);
                            multVal=range(obj.ckt.VDI.negTol,obj.ckt.VDI.plusTol);
                            Kv_(independentBranch,dependentBranch)=multVal(J);
                        end
                    case 'VDV'
                        Kv_(Indx,Indx)=eye(numel(Indx));
                        for J=1:numel(Indx)
                            dependentBranch=obj.desigIndx(char(obj.ckt.VDV.cntrlBranch(J)));
                            independentBranch=Indx(J);
                            multVal=range(obj.ckt.VDV.negTol,obj.ckt.VDV.plusTol);
                            Kv_(independentBranch,dependentBranch)=multVal(J);
                        end
                end            
            end
            % Make the matricies sparse
            varyS=sparse(S_);  
            varyKi=sparse(Ki_);
            varyKv=sparse(Kv_);
            % Function to create a random multiplier around unity that
            % covers the tolerance range for the component.
            function multVal=range(nTol,pTol)
                n=numel(nTol);
                multVal=1+(-nTol+(pTol+nTol).*rand(1,n))/100;
            end
        end
        
        function deltaT_IV=deltaTableau_matrix(T,f,IVn)
            % A circuit with its components at their nominal values can be
            % represented in tableau format as
            % 
            %      Tn*IVn = S
            % 
            % where Tn is the tableau matrix, IVn is the vector of branch
            % currents and node voltages, and S is the vector of
            % independent sources. The variation of branch currents and
            % node voltages can be explored by letting the components vary
            % from their nominal value. (Here we will assume that the
            % independent sources remain constant since source variation
            % would simply be a linear response to the stimulus.) In this
            % situation, the preturbed circuit represented in tableau
            % format is
            %
            %      (Tn+deltaT)*(IVn+deltaIV)=S
            %
            % When the first equation is subtracted from the second
            % equation, the result is
            %
            %      (Tn+deltaT)*(IVn+deltaIV)-Tn*IVn=0
            %
            % This equation can be simplified to
            %
            %      Tn*deltaIV=-deltaT*(IVn+deltaIV)
            %
            % Since the tableau matrix has an inverse, then
            %
            %      deltaIV=-inv(Tn)*deltaT*(IVn+deltaIV)
            %
            % The tableau system of equations has a property not shared by
            % all methods of formulating circuit equations; a component
            % value appears in one and only one row of the tableau format
            % at a single location within the row. Each location in a row
            % of deltaT has a number representing the change-in-value for
            % that component. For each row, the change in component value
            % can be extracted from deltaT and replaced by the appropriate
            % value of (IVn+delatIV) because there is only one matrix
            % element in a row. The equation becomes
            %
            %      deltaIV=-inv(Tn)*deltaT_IV*deltaComp
            %
            % where deltaT_IV is the matrix where the change in component
            % value has been replaced by (IVn+delatIV) and deltaComp is the
            % vector of changes in component values.
            %
            % In this form, the matrix 
            %
            %       -inv(Tn)*deltaT_IV
            %
            % maps the change-in-component-value space to the
            % change-in-current-and-voltage space.
            %
            % Example:
            % deltaT_IV=deltaTableau_matrix(T,f,IVn)
            %
            
            % Extract deltaKv and deltaKi to a local variable, preserving
            % the property value
            delta_Kv=T.deltaKv;
            delta_Ki=T.deltaKi;
            % Find all capacitors and set the appropriate parts of
            % T.deltaKv to the desired frequency.
            if any(strcmp('C',T.typesAvailable))
                Indx=T.typeIndx('C');
                delta_Kv(Indx,Indx)=1j*2*pi*f*T.deltaKv(Indx,Indx);
            end
            % Find all the inductors and set the appropriate parts of
            % T.deltaKi to the desired frequency.
            if any(strcmp('L',T.typesAvailable))
                Indx=T.typeIndx('L');
                delta_Ki(Indx,Indx)=1j*2*pi*f*T.deltaKi(Indx,Indx);
            end
            % Create the beginning of deltaT_IV. Frequency and imaginary
            %  components are present, (IV+deltaIV) are not yet added.
            [m,n]=size(T.A);
            Zero=sparse(m,n);
            deltaT_IV=[delta_Ki -delta_Kv*T.A';Zero T.Z];
            % Multiply each non-zero location by the appropriate value in
            % IVn 
            L=length(deltaT_IV);
            for I=1:L
                col=find(deltaT_IV(I,:));
                deltaT_IV(I,col)=IVn(I)*deltaT_IV(I,col); 
            end
        end
        
        function deltaComp=deltaComp_vector(T)
            % This function builds a structure that contains three vectors.
            % These vectors are deltaComp.nTol, deltaComp.pTol, and 
            % deltaComp.compVal. The vectors deltaComp.nTol and
            % deltaComp.pTol are the fractional limits
            % (tolerance_percent/100) of the components. The
            % deltaComp.compVal vector is the actual component value. The
            % vectors do not contain values for independent sources.
            % Example:
            %      deltaComp=deltaComp_vector(T)
            
            % Create various elements of deltaComp. Start with initializing
            % the values.
            L=T.TmIndx.Count;
            nTol(L,1)=0;
            pTol(L,1)=0;
            compVal(L,1)=0;
            desig{L,1}='';
            % Build the various elements from the types available
            for I=T.typesAvailable
                typeName=char(I);
                Indx=T.typeIndx(typeName);
                switch typeName
                    case 'R'
                        nTol(Indx)=T.ckt.R.negTol;
                        pTol(Indx)=T.ckt.R.plusTol;
                        compVal(Indx)=T.ckt.R.value;
                        desig(Indx)=T.ckt.R.designator;
                    case 'L'
                        nTol(Indx)=T.ckt.L.negTol;
                        pTol(Indx)=T.ckt.L.plusTol;
                        compVal(Indx)=T.ckt.L.value;
                        desig(Indx)=T.ckt.L.designator;
                    case 'C'
                        nTol(Indx)=T.ckt.C.negTol;
                        pTol(Indx)=T.ckt.C.plusTol;
                        compVal(Indx)=T.ckt.C.value;
                        desig(Indx)=T.ckt.C.designator;
                    case 'IDI'
                        for J=1:numel(Indx)
                            nTol(Indx(J))=T.ckt.IDI.negTol(J);
                            pTol(Indx(J))=T.ckt.IDI.plusTol(J);
                            compVal(Indx(J))=T.ckt.IDI.value(J);
                            desig(Indx(J))=T.ckt.IDI.designator(J);
                        end
                    case 'VDI'
                        for J=1:numel(Indx)
                            nTol(Indx(J))=T.ckt.VDI.negTol(J);
                            pTol(Indx(J))=T.ckt.VDI.plusTol(J);
                            compVal(Indx(J))=T.ckt.VDI.value(J);
                            desig(Indx(J))=T.ckt.VDI.designator(J);
                        end
                    case 'VDV'
                        for J=1:numel(Indx)
                            nTol(Indx(J))=T.ckt.VDV.negTol(J);
                            pTol(Indx(J))=T.ckt.VDV.plusTol(J);
                            compVal(Indx(J))=T.ckt.VDV.value(J);
                            desig(Indx(J))=T.ckt.VDV.designator(J);
                        end
                end
            end
            % Assign the elements of deltaComp to a structure
            deltaComp.nTol=nTol/100;
            deltaComp.pTol=pTol/100;
            deltaComp.compVal=compVal;
            deltaComp.desig=desig;
        end  
    end
end

Contact us