text:
 function mat_write(file, data, varargin)
%function mat_write(file, data,
%		['-check', '-nocheck', '-type', datatype, '-name', varname])
%
% Save matlab "data" to "file", first checking if "file" exists!
% If file suffix is .fld, then save as AVS .fld file in xdr_float.
%
% datatype can be 'xdr_float' 'xdr_int' ...
% Jeff Fessler, The University of Michigan

if nargin == 1 && streq(file, 'test'), mat_write_test, return, end
if nargin < 2, ir_usage, end

% use .fld file if filename ends in .fld, or if filename is t0, t1, ...
isfld = 0;
if length(file) > 4
	isfld = streq(file((end-3):end), '.fld');
elseif length(file) == 2 && file(1) == 't'
	isfld = true;
end

check1st = true;
if isfld
	varname = '';
	datatype = 'xdr_float';
else
	varname = 'data';
	datatype = '';
end

while length(varargin)
	if streq(varargin{1}, '-nocheck')
		check1st = 0;
		varargin = {varargin{2:end}};
	elseif streq(varargin{1}, '-check')
		check1st = 1;
		varargin = {varargin{2:end}};
	elseif streq(varargin{1}, '-type')
		if length(varargin) < 2
			error '-type needs arg'
		else
			datatype = varargin{2};
			varargin = {varargin{3:end}};
		end
	elseif streq(varargin{1}, '-name')
		if length(varargin) < 2
			error '-name needs arg'
		else
			varname = varargin{2};
			varargin = {varargin{3:end}};
		end
	else
		error(sprintf('bad arg "%s"', varargin{1}))
	end
end

if isfld && ~isempty(varname)
	warning('variable name for .fld file ignored')
end
if ~isfld && ~isempty(datatype)
	warning('datatype for .mat file ignored')
end

if 2==exist(file) && check1st
	if has_aspire
		if isfld
			os_run(['op range ' file])
		else
			os_run(['op -chat 0 matls ' file])
		end
	end
	t = sprintf('file "%s" exists.  overwrite? [y|n]: ', file);
	t = input(t, 's');
	if ~streq(t, 'y')
		return
	end
end

% avs file
if isfld
	if check1st
		fld_write(file, data, '-check', datatype)
	else
		fld_write(file, data, '-nocheck', datatype)
	end

% matlab file
else
	error 'writing to matlab file for reading by aspire no longer supported because mathworks changes the file format too often.  use fld_write'
	if ~streq(varname, 'data')
		eval([varname '= data'])	% caution: dangerous!
	end

	t = version;
	if t(1) == '4'
		save(file, varname)
	elseif t(1) == '5' || t(1) == '6'
		save(file, varname, '-v4')
	else
		error version
	end
end

printm('file "%s" written ok', file)

% mat_write_test()
function mat_write_test
dat = [5:8 1:4];
file = 'tmp.fld';
mat_write(file, dat, '-check', '-name', 'tester')
%d = load(file);
d = fld_read(file);
if any(d.tester ~= dat), error 'bug', else, printf('%s test ok', mfilename), end
delete(file)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function value = i4_gcdb ( i, j, k )

%*****************************************************************************80
%
%% I4_GCDB finds the greatest common divisor of the form K**N of two numbers.
%
%  Discussion:
%
%    Note that if J is negative, I4_GCDB will also be negative.
%    This is because it is likely that the caller is forming
%    the fraction I/J, and so any minus sign should be
%    factored out of J.
%
%    If I and J are both zero, I4_GCDB is returned as 1.
%
%    If I is zero and J is not, I4_GCDB is returned as J,
%    and vice versa.
%
%    If I and J are nonzero, and have no common divisor of the
%    form K**N, I4_GCDB is returned as 1.
%
%    Otherwise, I4_GCDB is returned as the largest common divisor
%    of the form K**N shared by I and J.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, two numbers whose greatest common divisor K**N
%    is desired.
%
%    Input, integer K, the possible divisor of I and J.
%
%    Output, integer VALUE, the greatest common divisor of
%    the form K^N shared by I and J.
%
  value = 1;
%
%  If both I and J are zero, I4_GCDB is 1.
%
  if ( i == 0 && j == 0 )
    value = 1;
    return
  end
%
%  If just one of I and J is zero, I4_GCDB is the other one.
%
  if ( i == 0 )
    value = j;
    return
  elseif ( j == 0 )
    value = i;
    return
  end
%
%  Divide out K as long as you can.
%
  if ( 0 < j )
    value = 1;
  else
    value = -1;
  end

  while

    if ( mod ( i, k ) ~= 0 || mod ( j, k ) ~= 0 )
      break
    end

    value = value * k;
    i = i / k;
    j = j / k;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_noiseless_s_4.mat';
N = 256;
S = 4;
pta = spx.pursuit.PhaseTransitionAnalysis(N);

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function w = set(w, varargin)
   %SET Set properties for waveform object(s)
   %   w = Set(w,prop_name, val, [prop_name, val])
   %   SET is one of the two gateway functions of an object, such as waveform.
   %   Properties that are changed through SET are typechecked and otherwise
   %   scrutinized before being stored within the waveform object.  This
   %   ensures that the other waveform methods are all retrieving valid data,
   %   thereby increasing the reliability of the code.
   %
   %   Another strong advantage to using SET and GET to change  and retrieve
   %   properties, rather than just assigning them to the waveform directly,
   %   is that the underlying data structure can change and grow without
   %   harming the code that is written based on the waveform object.
   %
   %   Valid property names:
   %
   %       STATION, CHANNEL, FREQ, DATA, UNITS: self explanatory
   %       START :         Set the start-time using matlab format or string
   %       START_EPOCH:    Set the start-time using epoch (antelope) time
   %       SAMPLE_LENGTH : adjusts the length of the data to length VAL.
   %                       If this number is shorter than existing data, SET
   %                       truncates after sample VAL.  If it is longer, then
   %                       SET will pad the end with zeroes to get length VAL.
   %
   %       If user-defined fields were added to the waveform (ie, through
   %       addField), these fieldnames are also available through set.
   %
   %       for example
   %           % create a waveform object, and add a field called DUMMY with
   %           % the value 12
   %           w = addfield(waveform,'DUMMY',12);
   %
   %           % change the value of the DUMMY field
   %           w = set(w,'DUMMY',32);
   %
   %           % change the station and channel at once
   %           w = set(w,'Station','YCK', 'Channel', 'BHZ');
   %
   %   Batch changes can be made if input w is a matrix (use with care!)
   %
   %   NOTE: There is no set support for "version".  This is handled entirely
   %         within the waveform constructor and the loadobj function.
   %
   % NOTE: To avoid logging the history, append the argument 'nohist'
   %
   %  See also WAVEFORM/GET

   % VERSION: 1.1 of waveform objects
   % AUTHOR: Celso Reyes (celso@gi.alaska.edu)
   % LASTUPDATE: March 15, 2009

   %4/22/08 added synonym for sample_length of data_length

   %global dep2mep

   % check to see if history is bypassed

   idx = numel(varargin);
   while idx > 1
      prop_name = upper(varargin{idx-1});
      val = varargin{idx};
      idx = idx - 2;
      %{
Vidx = 1 : numel(varargin);

while numel(Vidx) >= 2
  prop_name = upper(varargin{Vidx(1)});
  val = varargin{Vidx(2)};
      %}

      %for n = 1 : numel(w);
      %out = w(n);

      switch prop_name
         case 'SCNLOBJECT'
            switch class(val)
               case 'scnlobject'
                  [w.cha_tag] = deal(get(val,'channeltag'));
               case 'ChannelTag'
                  [w.cha_tag] = deal(val);
               case 'char'
                  [w.cha_tag] = deal(ChannelTag.parse(val));
               otherwise
                  % try it!
                  [w.cha_tag] = deal(ChannelTag(val));
                  % error('Waveform:set:propertyTypeMismatch','Expected a SCNLOBJECT or CHANNELTAG');
            end
            % whine while we refractor
            warning('Use set(ChannelTag) instead') % TODO: Decide if better generic name exists

         case {'CHANNELINFO', 'CHANNELTAG'}
            switch class(val)
               case 'scnlobject'
                  [w.cha_tag] = deal(ChannelTag(char(val)));
               case 'ChannelTag'
                  [w.cha_tag] = deal(val);
               case 'char'
                  [w.cha_tag] = deal(ChannelTag(val));
               otherwise
                  % try it!
                  [w.cha_tag] = deal(ChannelTag(val));
                  % error('Waveform:set:propertyTypeMismatch','Expected a SCNLOBJECT or CHANNELTAG');
            end
         case {'STATION','NETWORK', 'CHANNEL','LOCATION'}
            if ~ischar(val)
               error('Waveform:set:propertyTypeMismatch',...
                  '%s should be a string not a %s', prop_name,class(val));
            end
            all_chantags = set([w.cha_tag],prop_name,val);
            for n=1:numel(w)
               w(n).cha_tag = all_chantags(n);
            end

         case {'FS', 'FREQ'}
            if ~isnumeric(val)
               error('Waveform:set:propertyTypeMismatch',...
                  'Frequency should be numeric, not %s', class(val));
            end
            if numel(w)==1
               w.Fs = val(1);
            else
               [w.Fs] = deal(val(1));
            end

         case {'START', 'START_MATLAB'}
            %AUTOMATICALLY figures out whether date is antelope or matlab
            %format;
            if numel(val) == 1
               %expecting a number
               if numel(w) == 1
                  w.start = val;
               else
                  [w.start] = deal(val);
               end
            else
               %might be char, might be vector. give it a shot.
               if numel(w) == 1
                  w.start = datenum(val);
               else
                  [w.start] = deal(datenum(val));
               end
               % else
               %   error('Waveform:set:propertyTypeMismatch',...
               %    'Start time not assigned... Unknown value type: %s', class(val));
            end

         case {'START_ANTELOPE', 'START_EPOCH'}
            if ~isnumeric(val),
               error('Waveform:set:propertyTypeMismatch',...
                  'Epoch time should be numeric, not %s', class(val));
            end
            [w.start] = deal(datenum(dep2mep(val)));

         case 'DATA'
            if ~isnumeric(val)
               error('Waveform:set:propertyTypeMismatch',...
                  'DATA should be numeric, not %s', class(val));
            end
            if numel(w) == 1
               w.data = val(:);
            else
               [w.data] = deal(val(:)); %always a column
            end

         case 'UNITS'
            if ~ischar(val)
               error('Waveform:set:propertyTypeMismatch',...
                  'UNITS should be a string, not %s',class(val));
            end
            [w.units] = deal(val);

         case {'SAMPLE_LENGTH','DATA_LENGTH'} % chop data or zero pad to specific length
            if ~isnumeric(val),
               error('Waveform:set:propertyTypeMismatch',...
                  'SAMPLE_LENGTH should be numeric, not %s',class(val));
            end
            for n=1:numel(w);
               if numel(w(n).data) >= val
                  w(n).data = w(n).data(1:val,1);  % shrink the data
               else
                  w(n).data(numel(w(n).data)+1 :val , 1) = 0; %zero pad
               end
            end

         case 'HISTORY'
            w = addhistory(w, val);

         otherwise
            for n=1:numel(w)
               switch prop_name
                  case w(n).misc_fields
                     %mask = ismember(w(n).misc_fields, prop_name);
                     mask = strcmp(prop_name,w(n).misc_fields);
                     w(n).misc_values(mask) = {val};
                  otherwise
                     error('Waveform:set:unknownProperty',...
                        'can''t understand property name : %s', prop_name);
               end %switch
            end %n
      end %switch

      % Vidx(1:2) = []; %done with those parameters, move to the next ones...

   end
end

author: geoscience-community-codes

repo: GISMO


text:
function [F, p, resid, df_model, df_error] = F_test_full_vs_red(y, X, Xred, px, pxred)
% :Usage:
% ::
%
%     [F, p, resid, df_model, df_error] = F_test_full_vs_red(y, X, Xred, px, pxred)
%
% :Examples:
% ::
%
%    X = randn(100, 3); Xred = X(:,1); y = X(:,2) + randn(100, 1);
%    px = pinv(X); pxred = pinv(Xred);
%    [F, p, resid] = F_test_full_vs_red(y, X, Xred, px, pxred);
%
%
%    % Test full-model F-value against regress.m
%    Xred = X(:,end); % intercept only
%    px = pinv(X);
%    pxred = pinv(Xred);
%    [F, p, resid, dfm, dfe] = F_test_full_vs_red(y, X, Xred, px, pxred); % full model F-test
%    [b, bint, r, rint, stats] = regress(y, X);
%
% ..
%    Tested OK on 11/27/07, by tor
% ..

    T = length(y);      % Length of time course

    k = size(px, 1);        % predictors: full model
    kred = size(pxred, 1);  % predictors: reduced model

    % Degrees of freedom: model: Full - reduced
    df_model = k - kred;                        % degrees of freedom for model (param - 1)
    df_error = T - k;                           % error DF, full model

    % Step 1: Find the OLS solution for the FULL model
    % ---------------------------------------------------
    beta = px * y;                             % Beta values
    resid = y - X * beta;                       % Residuals

    % Sums of squares
    %SST = y' * y;
    SSE = resid' * resid;
    %SSfull = SST - SSE;
    var_est = SSE / df_error;                   % Estimate of Sigma^2

    % Step 2: Find the OLS solution for the REDUCED model
    % F-test for full vs. reduced
    % ---------------------------------------------------
    betared = pxred * y;                          % Beta values
    residred = y - Xred * betared;
    SSEred = residred' * residred;
    %SSred = betared' * Xred' * y;                      % Full model sum of squares

    % F stat
    % (SSred - SSfull) ./ (var * df_model)
    F = (SSEred - SSE) / (var_est * df_model);         % F-statistic - compare with F-distribution with (param-1, df) degrees of freedom

    p = 1 - fcdf(F, df_model, df_error);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function fiff_write_complex_matrix(fid,kind,mat)
%
% fiff_write_complex_matrix(fid,kind,mat)
%
% Writes a single-precision complex matrix tag
%
%     fid           An open fif file descriptor
%     kind          The tag kind
%     mat           The data matrix
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.1  2006/09/23 14:43:56  msh
%   Added routines for writing complex and double complex matrices.
%   Added routine for writing double-precision real matrix.
%
%

me='MNE:fiff_write_complex_matrix';

if nargin ~= 3
        error(me,'Incorrect number of arguments');
end

FIFFT_COMPLEX_FLOAT = 20;
FIFFT_MATRIX  = bitshift(1,30);
FIFFT_MATRIX_COMPLEX = bitor(FIFFT_COMPLEX_FLOAT,FIFFT_MATRIX);
FIFFV_NEXT_SEQ=0;

datasize = 2*4*numel(mat) + 4*3;

count = fwrite(fid,int32(kind),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFT_MATRIX_COMPLEX),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(datasize),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFV_NEXT_SEQ),'int32');
if count ~= 1
   error(me,'write failed');
end
nrow = size(mat,1);
ncol = size(mat,2);
for j = 1:nrow
   for k = 1:ncol
      count = fwrite(fid,real(mat(j,k)),'single');
      if count ~= 1
          error(me,'write failed');
      end
      count = fwrite(fid,imag(mat(j,k)),'single');
      if count ~= 1
          error(me,'write failed');
      end
   end
end
dims(1) = size(mat,2);
dims(2) = size(mat,1);
dims(3) = 2;
count = fwrite(fid,int32(dims),'int32');
if count ~= 3
    error(me,'write failed');
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% SnirfSave( filenm, snirf )
%
% This function will save a snirf object to file 'filenm'.
% Recall that you create a snirf object with
%   snirf = SnirfClass();
% and then populate it with the required fields indicated by the SNIRF
% specification.

function SnirfSave( filenm, snirf )

snirf.Save( filenm )

author: BUNPC


text:
clear;
data = importdata('road1_fiber_ms50_r10.csv');
data = data.data;

figure('DefaultAxesFontSize',16);
ax = gca;

hold on;plot(1:144,data(1:144,4),'k-.','linewidth',0.5);
hold on;plot(144+1:2*144,data(144+1:2*144,4),'color',[0,0.4470,0.7410],'linewidth',2);
hold on;plot(2*144+1:3*144,data(2*144+1:3*144,4),'k-.','linewidth',0.5);
hold on;plot(3*144+1:5*144,data(3*144+1:5*144,4),'color',[0,0.4470,0.7410],'linewidth',2);
hold on;plot(5*144+1:6*144,data(5*144+1:6*144,4),'k-.','linewidth',0.5);
hold on;plot(6*144+1:7*144,data(6*144+1:7*144,4),'color',[0,0.4470,0.7410],'linewidth',2);

ax.XTick = [144,2*144,3*144,4*144,5*144,6*144,7*144];
ax.XTickLabel = [' ';' ';' ';' ';' ';' ';' '];
ax.YTick = [0,60];
ax.YTickLabel = [' ';' '];
x = [144,2*144,2*144,144];
y = [10,10,50,50];

v1 = [144,0; 2*144,0; 2*144,60; 144,60];
v2 = v1; v2(:,1)=v1(:,1)+2*144;
v3 = v2; v3(:,1)=v2(:,1)+144;
v4 = v3; v4(:,1)=v3(:,1)+2*144;
f = [1,2,3,4];
patch('faces',f,'vertices',v1,'edgecolor',[0.5 0.5 0.5],'facecolor','green','facealpha',0.1);
patch('faces',f,'vertices',v2,'edgecolor',[0.5 0.5 0.5],'facecolor','green','facealpha',0.1);
patch('faces',f,'vertices',v3,'edgecolor',[0.5 0.5 0.5],'facecolor','green','facealpha',0.1);
patch('faces',f,'vertices',v4,'edgecolor',[0.5 0.5 0.5],'facecolor','green','facealpha',0.1);
xlim([10,50]);
xlim([0,7*144]);

grid on;
box off;

set(gcf, 'PaperSize', [6 2]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 6 2]);
saveas(gcf,'curve1','pdf');

author: xinychen

repo: academic-drawing


text:
classdef InstanceSegMMLoss < dagnn.Loss
    properties
        marginAlpha_
    end
    properties (Transient)
        curSimMat
        weightMat
        mass_
        size_
    end

    methods
        function outputs = forward(obj, inputs, params)
            [h, w, ch, bs] = size(inputs{1});
            sz = [h, w, ch, bs];
            mass = sz(1) * sz(2) + 1;

            obj.size_ = sz;
            obj.mass_ = mass;

%             gpuMode = isa(inputs{1}, 'gpuArray');
%             if gpuMode
%                 grndLabel = gpuArray(zeros(sz, 'single'));
%             else
%                 grndLabel = zeros(sz, 'single');
%             end
%
%             for j = 1:sz(4)
%                 C = inputs{2}(:,:,:,j);
%                 C = reshape(C, [sz(1), 1]);
%                 C = repmat(C, 1, sz(2));
%                 C = (C==C');
%                 grndLabel(:,:,:,j) = C;
%             end
%             obj.curSimMat = grndLabel;

            layerName = 'obj_instSeg';
            if isnan(obj.net.getLayerIndex(layerName))
                layerName = 'obj_instSeg_reg';
            end
            obj.curSimMat = obj.net.layers(obj.net.getLayerIndex(layerName)).block.curSimMat;
            obj.weightMat = obj.net.layers(obj.net.getLayerIndex(layerName)).block.weightMat;

            outputs{1} = vl_nnloss(inputs{1}, obj.curSimMat, 'marginAlpha_', obj.marginAlpha_, ...
                'loss', obj.loss, ...
                'instanceWeights', obj.weightMat ) ; % 1./mass
            n = obj.numAveraged ;
            m = n + size(inputs{1},4) ;
            obj.average = (n * obj.average + double(gather(outputs{1}))) / m ;
            obj.numAveraged = m ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
%             sz = obj.size_;
%             mass = obj.mass_;
%             grndLabel = obj.curSimMat;

            derInputs{1} = vl_nnloss(inputs{1}, obj.curSimMat, derOutputs{1}, 'marginAlpha_', obj.marginAlpha_, ...
                'loss', obj.loss, ...
                'instanceWeights', obj.weightMat ) ; % 1./mass
            derInputs{2} = [] ;
            derParams = {} ;
        end

        function obj = InstanceSegMMLoss(varargin)
            obj.load(varargin) ;
        end
    end
end

author: aimerykong


text:
%FINDHOMOGRAPHY  Finds a perspective transformation between two planes
%
%     H = cv.findHomography(srcPoints, dstPoints)
%     [H, mask] = cv.findHomography(...)
%     [...] = cv.findHomography(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __srcPoints__ Coordinates of the points in the original plane, a numeric
%   array of size Nx2/1xNx2/Nx1x2 or cell array of 2-elements vectors
%   `{[x,y], ...}` (single floating-point precision).
% * __dstPoints__ Coordinates of the points in the target plane, of same size
%   and type as `srcPoints`.
%
% ## Output
% * __H__ 3x3 Homography matrix.
% * __mask__ Nx1 mask array of same length as input points, indicates inliers
%   (which points were actually used in the best computation of `H`).
%
% ## Options
% * __Method__ Method used to compute a homography matrix. The following
%   methods are possible:
%   * __0__ a regular method using all the points, i.e. the least squares
%     method (default)
%   * __Ransac__ RANSAC-based robust method.
%   * __LMedS__ Least-Median of squares robust method.
%   * __Rho__ PROSAC-based robust method, introduced in [Bazargani15].
%     (weighted RANSAC modification, faster in the case of many outliers).
% * __RansacReprojThreshold__ Maximum allowed reprojection error to treat a
%   point pair as an inlier (used in the RANSAC and RHO methods only). That
%   is, if
%   `|| dstPoints_i - convertPointsToHomogeneous(H*srcPoints_i) ||_2 > RansacReprojThreshold`
%   then the point `i` is considered as an outlier. If `srcPoints` and
%   `dstPoints` are measured in pixels, it usually makes sense to set this
%   parameter somewhere in the range of 1 to 10. default 3.0.
% * __MaxIters__ The maximum number of RANSAC iterations. default 2000
% * __Confidence__ Confidence level, between 0 and 1. default 0.995
%
% The function finds and returns the perspective transformation `H` between
% the source and the destination planes:
%
%     s_i * [x_i'; y_i'; 1] ~ H * [x_i; y_i; 1]
%
% so that the back-projection error:
%
%     sum_{i} (x_i' - (h11*x_i + h12*y_i + h13)/(h31*x_i + h32*y_i + h33))^2 +
%             (y_i' - (h21*x_i + h22*y_i + h23)/(h31*x_i + h32*y_i + h33))^2
%
% is minimized. If the parameter method is set to the default value 0, the
% function uses all the point pairs to compute an initial homography estimate
% with a simple least-squares scheme.
%
% However, if not all of the point pairs `(srcPoints_i, dstPoints_i)` fit the
% rigid perspective transformation (that is, there are some outliers), this
% initial estimate will be poor. In this case, you can use one of the three
% robust methods. The methods RANSAC, LMedS and RHO try many different
% random subsets of the corresponding point pairs (of four pairs each,
% collinear pairs are discarded), estimate the homography matrix using this
% subset and a simple least-squares algorithm, and then compute the
% quality/goodness of the computed homography (which is the number of inliers
% for RANSAC or the least median re-projection error for LMedS). The best
% subset is then used to produce the initial estimate of the homography matrix
% and the mask of inliers/outliers.
%
% Regardless of the method, robust or not, the computed homography matrix
% is refined further (using inliers only in case of a robust method) with
% the Levenberg-Marquardt method to reduce the re-projection error even
% more.
%
% The methods RANSAC and RHO handle practically any ratio of outliers but
% need a threshold to distinguish inliers from outliers. The method LMedS
% does not need any threshold but it works correctly only when there are
% more than 50% of inliers. Finally, if there are no outliers and the noise
% is rather small, use the default method (`Method=0`).
%
% The function is used to find initial intrinsic and extrinsic matrices.
% Homography matrix is determined up to a scale. Thus, it is normalized so
% that `h33 = 1`. Note that whenever an `H` matrix cannot be estimated, an
% empty one will be returned.
%
% ## References
% [Fischler81]:
% > Fischler, M. A., and R. C. Bolles. "Random sample consensus: A paradigm
% > for model fitting with applications to image analysis and automated
% > cartography", Communications of the Association for Computing Machinery
% > 24 (1981): 381-395.
%
% [Rousseeuw84]:
% > Rousseeuw, P. J. "Least median of squares regression", Journal of the
% > American Statistical Association, 79 (1984): 871-880.
%
% [Inui03]:
% > Inui, K., S. Kaneko, and S. Igarashi. "Robust Line Fitting using LMedS
% > Clustering", Systems and Computers in Japan 34 (2003): 92-100.
%
% [Bazargani15]:
% > Hamid Bazargani, Olexa Bilaniuk, and Robert Laganiere. "A fast and robust
% > homography scheme for real-time planar target detection". Journal of
% > Real-Time Image Processing (2015): 1-20.
%
% See also: cv.getAffineTransform, cv.estimateAffine2D,
% cv.estimateAffinePartial2D, cv.getPerspectiveTransform,
% cv.estimateRigidTransform, cv.warpPerspective, cv.perspectiveTransform
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function vw = refreshFlatLevelView(vw,recomputeFlag)
%
% vw = refreshFlatLevelView(vw,recomputeFlag);
%
% This is a version of refreshView that works on the
% Flat Across-Gray-Levels View. There are 2 main reasons
% I'm breaking it off into a separate function:
%
%   1) It allows for working off a mosaic of several layers
%      at once, which may be useful for things like inplane
%      views but would also make the refreshView code confusing;
%
%   2) It's still being tested out, so it's structure is likely
%      to change.
%
% But basically, it refreshes the flat multi view interface
% based off the settings in the UI controls.
%
% ras, 08/04, off refreshView

if ~exist('recomputeFlag','var'), recomputeFlag=1; end

% Make this one the selectedInplane, selectedVolume, or selectedFlat
if iscell(vw), selectView(vw); end

% Set window title
% Check for sessionCode - this should have been set in openxxxxWindow but
% some other programs (like mrAlign) don't beshave properly.
if (~isfield(vw,'sessionCode'))
    vw.sessionCode=pwd;
end

ui = viewGet(vw,'ui');

% Get colormap, numGrays, numColors and clipMode
str = [vw.ui.displayMode,'Mode'];

% Make sure the phMode, or ampMode structure exists.  If not, then set them
% to the default values.
if ~checkfields(vw,'ui',str), vw=resetDisplayModes(vw); end

% Here is the mode, get the relevant parameters.
modeStr=['vw.ui.',str];
mode = eval(modeStr);
cmap = mode.cmap;
numGrays = mode.numGrays;
numColors = mode.numColors;
clipMode = mode.clipMode;

% get/set flat rotations
[rotations,flipLR] = getFlatRotations(vw);

if (~(recomputeFlag==2)) % Only set when the hemisphere changes
    % Retreive the imageRotation from the UI
    thisRotation = getImageRotate(vw);

    % Set the corresponding value in degrees in FLAT.imageRotateDegrees
    rotations(viewGet(vw, 'Current Slice')) = thisRotation*(180/pi);
end

vw = setFlatRotations(vw,rotations,flipLR);

% recompute the image if selected
if (recomputeFlag | isempty(vw.ui.image))
    vw = recomputeFlatLevelImage(vw);
end

% Select the window
set(0,'CurrentFigure',vw.ui.windowHandle);

% Update annotation string
set(vw.ui.annotationHandle,'string',annotation(vw,getCurScan(vw)));

% Display final image
dispMode = sprintf('%sMode',viewGet(vw,'displayMode'));
cmap = vw.ui.(dispMode).cmap;
imshow(vw.ui.image);
colormap(cmap);

% Draw colorbar
if isempty(vw.ui.cbarRange)
    setColorBar(vw,'off',[numGrays+1,numGrays+numColors]);
else
    setColorBar(vw,'on',[numGrays+1,numGrays+numColors]);
end

% Draw ROIs
if vw.ui.showROIs > 0,     drawROIsMontage(vw);
elseif vw.ui.showROIs < 0, drawROIsPerimMontage(vw);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%  Figure 10.59      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_59.m is a script for disk read/write head control design
% the design is based on 1/s^2 with robustness
% with respect to a resonance of damping zeta = z at wm.
% the design calls for a single lead selected to maximize wc
% the compensator has a second order roll-off filter
% the resonance is 'gain stabilized' with GM
% the phase margin is set by alpha = a
% the time is scaled to milliseconds
wm= 5*pi
% input the nominal plant with a gain of 1.
numGo=1;
z=.05;
GM = 4;
a=0.1;
denGo= [1 0 0];
sysGo=tf(numGo,denGo);
% input the plant with the resonance
numG= [1/(50*pi) 1];
denG=[1/(25*pi^2) 1/(50*pi) 1  0 0];
sysG=tf(numG,denG);
% Define the design constant
A = ((a*2*z)/(GM))^.333;
% PM >50; we have selected alpha = a = 0.1
Tinv= A*wm*sqrt(a);
T=1/Tinv;
wc =  A*wm
% K =  wc^2;
K=  wc^2;
% Input the lead compensation
numD =K*sqrt(a)*[T 1];
denD = [a*T 1];
sysD=tf(numD,denD);
% Now design the roll-off filter
w1=wm*sqrt(A/sqrt(a));
% Design the roll-off damping
z1=.3 ;
numF=1;
denF=[1/w1^2  2*z1/w1 1];
sysF=tf(numF,denF);
% now compute the entire compensation
sysDc=sysD*sysF;
sysOLo=sysD*sysGo;
sysOL=sysDc*sysG;
clf;
% bode(sysOLo);
% hold on
w=logspace(-1,2);
[mag,ph]=bode(sysOL,w);
magdb=20*log10(mag(:,:));
magdb1=[magdb; zeros(size(magdb))];
subplot(211);
semilogx(w,magdb1);
grid on;
title('Fig. 10.59 Disk drive control with gain-stabilized resonance and roll-off filter');
xlabel('\omega (rad/sec)');
ylabel('Magnitude (db)');

subplot(212);
semilogx(w,[ph(:,:); -180*ones(size(ph(:,:)))]);
grid on;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
[Gm,Pm,Wcg,Wcp] = margin(mag,ph,w)

author: Sable

repo: mcbench-benchmarks


text:
function [success ips] = genInterestPoints(moviename,varargin)

[ipnames,doaskleftright] = myparse(varargin,'ipnames',{'food','mouth','perch'},...
  'doaskleftright',true);

annObj = HandleObj; % handle obj to store annotation info
hFig = playfmf('moviename',moviename,'annDataObj',annObj,'ipnames',ipnames,'doaskleftright',doaskleftright);
waitfor(hFig);
if isempty(annObj.data)
  success = false;
  ips = [];
  return;
end

fprintf(1,'Interest points generated:\n');
disp(annObj.data);

ips = struct();
for i = 1:numel(ipnames),
  ips.(ipnames{i}) = annObj.data.(ipnames{i});
end
if doaskleftright
  ips.face = annObj.data.lr;
end
success = true;

author: kristinbranson

repo: JAABA


text:
function [cc] = in32cc(in3)
% Convert volume from cubic inches to cubic centimeters*.
% * Not to be confused with cubic cubits.
% Chad Greene 2012
cc = in3*16.387064;


text:
function [curv] = write_curv(fname, curv, fnum)
% [curv] = write_curv(fname, curv, fnum)
%
% writes a curvature vector into a binary file
%				fname - name of file to write to
%				curv  - vector of curvatures
%				fnum  - # of faces in surface.
%

%
% write_curv.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

% open it as a big-endian file
fid = fopen(fname, 'wb', 'b') ;
vnum = length(curv) ;
NEW_VERSION_MAGIC_NUMBER = 16777215;
fwrite3(fid, NEW_VERSION_MAGIC_NUMBER ) ;
fwrite(fid, vnum,'int32') ;
fwrite(fid, fnum,'int32') ;
fwrite(fid, 1, 'int32');
fwrite(fid, curv, 'float') ;
fclose(fid) ;


text:
function Y = spm_data_read(V,varargin)
% Read data from disk [Y = V(I)]
% FORMAT Y = spm_data_read(V)
% V        - a structure array (see spm_data_hdr_read)
% Y        - an array of data values; the last dimension indexes numel(V)
%
% FORMAT Y = spm_data_read(V,'slice',S)
% V        - a structure array of image volumes (see spm_data_hdr_read)
% S        - an array of slice indices
% Y        - an array of data values with dimensions (x,y,s,v)
%
% FORMAT Y = spm_data_read(V,'xyz',XYZ)
% V        - a structure array (see spm_data_hdr_read)
% XYZ      - a [n x m] array of m coordinates {voxel (n=3 or 4)/vertex (n=1)}
% Y        - an array of data values with dimensions (v,m)
%
% FORMAT Y = spm_data_read(V,I1,I2,...)
% V        - a structure array (see spm_data_hdr_read)
% I1,I2,...- subscript arrays
% Y        - an array of data values with dimensions (v,m)
%__________________________________________________________________________
% Copyright (C) 2012-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_data_read.m 7577 2019-04-24 08:59:56Z guillaume $

if ~isstruct(V)
    V = spm_data_hdr_read(V);
end

cl = class(V(1).private);
if isfield(V(1),'dat') && ~isequal(cl,'gifti'), cl = 'nifti'; end

switch cl
    case 'nifti'
        if isempty(varargin)
            % Y = full(V.private.dat); % if numel(V)==1, is faster
            Y = spm_read_vols(V);
        elseif ischar(varargin{1}) && ~isequal(varargin{1},':')
            switch lower(varargin{1})
                case 'slice'
                    for i=1:numel(V), for p=1:numel(varargin{2})
                        Y(:,:,p,i) = spm_slice_vol(V(i),spm_matrix([0 0 varargin{2}(p)]),V(i).dim(1:2),0);
                    end, end
                    if numel(V)==1, Y=Y(:,:,:,1); end
                case 'xyz'
                    Y = spm_get_data(V,varargin{2});
                otherwise
                    error('Unknown input option.');
            end
        else
            indices = varargin;
            n = get_ndata(V(1).dim,indices{:});
            Y = zeros(numel(V),prod(n));
            for i=1:numel(V)
                if numel(indices) == 1
                    ind = {indices{1} + (V(i).n(1)-1)*prod(V(i).dim)};
                else
                    ind = indices;
                end
                Y(i,:) = reshape(V(i).private.dat(ind{:}),1,[]);
            end
        end

    case 'gifti'
        indices = varargin;
        if isempty(indices)
            indices = repmat({':'},1,ndims(V));
        elseif strcmpi(indices{1},'xyz')
            indices = {indices{2}(1,:)};
        end
        n = get_ndata(V(1).dim,indices{:});
        Y = zeros(numel(V),prod(n));
        for i=1:numel(V)
            Y(i,:) = reshape(V(i).private.cdata(indices{:}),1,[]);
        end
        if isempty(varargin), Y = Y'; end % to be coherent with spm_read_vols

    otherwise
        error('Unknown data type.');
end

%==========================================================================
function n = get_ndata(dim,varargin)
n = zeros(1,numel(varargin));
for i=1:numel(varargin)
    if isequal(varargin{i},':')
        if i==numel(varargin)
            n(i) = dim(i); %prod(dim(i:end));
        else
            n(i) = dim(i);
        end
    else
        n(i) = numel(varargin{i});
    end
end

author: spm


text:
function varargout = lvd_EditActionSetPwrSinkStateGUI(varargin)
% LVD_EDITACTIONSETPWRSINKSTATEGUI MATLAB code for lvd_EditActionSetPwrSinkStateGUI.fig
%      LVD_EDITACTIONSETPWRSINKSTATEGUI, by itself, creates a new LVD_EDITACTIONSETPWRSINKSTATEGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITACTIONSETPWRSINKSTATEGUI returns the handle to a new LVD_EDITACTIONSETPWRSINKSTATEGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITACTIONSETPWRSINKSTATEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITACTIONSETPWRSINKSTATEGUI.M with the given input arguments.
%
%      LVD_EDITACTIONSETPWRSINKSTATEGUI('Property','Value',...) creates a new LVD_EDITACTIONSETPWRSINKSTATEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditActionSetPwrSinkStateGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditActionSetPwrSinkStateGUI_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 lvd_EditActionSetPwrSinkStateGUI

% Last Modified by GUIDE v2.5 06-Aug-2020 15:23:10

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditActionSetPwrSinkStateGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditActionSetPwrSinkStateGUI_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 lvd_EditActionSetPwrSinkStateGUI is made visible.
function lvd_EditActionSetPwrSinkStateGUI_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 lvd_EditActionSetPwrSinkStateGUI (see VARARGIN)

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

    centerUIFigure(hObject);

    action = varargin{1};
    setappdata(hObject,'action',action);

    lv = varargin{2};
    setappdata(hObject,'lv',lv);

    populateGUI(handles, action, lv);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, action, lv)
    [powerSinkListStr, powerSinks] = lv.getPowerSinksListBoxStr();
    set(handles.refPwrCompCombo,'String',powerSinkListStr);

    if(not(isempty(action.sink)))
        ind = find(powerSinks == action.sink,1,'first');
    else
        ind = [];
    end

    if(not(isempty(ind)))
        set(handles.refPwrCompCombo,'Value',ind);
    end

    if(action.activeStateToSet)
        set(handles.stateCombo,'Value',1);
    else
        set(handles.stateCombo,'Value',2);
    end

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditActionSetPwrSinkStateGUI_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
    if(isempty(handles))
        varargout{1} = false;
    else
        action = getappdata(hObject,'action');
        lv = getappdata(hObject,'lv');

        [~, powerSinks] = lv.getPowerSinksListBoxStr();
        ind = get(handles.refPwrCompCombo,'Value');
        action.sink = powerSinks(ind);

        contents = cellstr(get(handles.stateCombo,'String'));
        stateStr = contents{get(handles.stateCombo,'Value')};

        switch stateStr
            case 'Active'
                state = true;
            case 'Inactive'
                state = false;
            otherwise
                error('Invalid state string: %s', stateStr);
        end

        action.activeStateToSet = state;

        varargout{1} = true;
        close(handles.lvd_EditActionSetPwrSinkStateGUI);
    end

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

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

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

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

% --- Executes during object creation, after setting all properties.
function refPwrCompCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to refPwrCompCombo (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 selection change in stateCombo.
function stateCombo_Callback(hObject, eventdata, handles)
% hObject    handle to stateCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

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

% --- Executes during object creation, after setting all properties.
function stateCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to stateCombo (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 key press with focus on lvd_EditActionSetPwrSinkStateGUI or any of its controls.
function lvd_EditActionSetPwrSinkStateGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditActionSetPwrSinkStateGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditActionSetPwrSinkStateGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function LnX = coefficient(x,y,L)
% Lagrange form for the polynomial interpolation
% this function produces elements of the form y_k*L_k(x)
n = length(x);
Num = '';
Den = '';
syms X

for i = 1:n
    if i ~= L
    TempNum = strcat('(X-x(',num2str(i),')',')','*');
    Num = strcat(Num,TempNum);
    TempDen = strcat('(x(',num2str(L),')','-','x(',num2str(i),')',')','*');
    Den = strcat(Den,TempDen);
    end
end

Num(end) = []; Den(end) = [];
LnX = (eval(strcat('(',Num,')','/','(',Den,')')))*y(L);


text:
function dispfrac(mode)
% DISPFRAC Display terms in laplacian format.
%          DISPFRAC displays the current terms in fraction format.

% Author: Craig Borghesani
% Date: 8/23/94
% Revised: 10/18/94
% Copyright (c) 1999, Prentice-Hall

% obtain handle information
f=gcf;
ui_data     = get(f,'userdata');
ui_han      = ui_data{1};
ui_obj      = ui_data{2};
%popup       = ui_han(1);
popup       = ui_obj.SysUI(5);
term_axis   = ui_han(10);
plant       = ui_han(30);
for_cont    = ui_han(31);
bac_cont    = ui_han(32);
set_to_one1 = ui_han(33);
set_to_one2 = ui_han(34);
cont_type   = ui_han(35);
num_den     = ui_han(36);
pid_1       = ui_han(37);
pid_2       = ui_han(38);
component   = ui_han(44);
open_loop   = ui_han(45);
closed_loop = ui_han(46);
factored    = ui_han(47);
coefficient = ui_han(48);
inverse     = ui_han(49);
neg_loc     = ui_han(68);
cur_sys     = get(ui_han(30),'userdata');
cur_page    = get(ui_han(31),'userdata');
sys_state   = get(ui_han(33),'userdata');
page_setup  = get(ui_han(34),'userdata');
stat_bar    = get(ui_han(43),'userdata');
lab_txt     = get(ui_han(44),'userdata');
ers = 'norm';
fontname = 'times';

% define various shades of grey
grey = get(0,'defaultuicontrolbackground');
ltgrey = [0.5,0.5,0.5]*1.5;
dkgrey = [0.5,0.5,0.5]*0.5;

if nargin == 0, mode = 0; end

% update checked properties of menu items
if mode,
 if any(abs(mode) == [1,2,3]), % switching systems

  if cur_sys == mode, return; end

  eraser1 = ui_obj.eraser1;
  tab_uis = ui_obj.SysTabs;

   %sel_dims = [-3,0,+6,+2];
   %unsel_dims = [3,0,-6,-2];

  sel_dims = [0,0,0,+2];
  unsel_dims = [0,0,0,-2];

  sel_tab = findobj(f,'tag','selected_tab1');
  cur_obj = tab_uis(abs(mode));
  if mode > 0 & sel_tab == cur_obj, return; end

  if length(sel_tab),
     tab_pos = get(sel_tab,'pos');
     set(sel_tab,'pos',tab_pos+unsel_dims,'tag','',...
                 'back',ltgrey);
  end

  cur_pos = get(cur_obj,'pos');
  set(cur_obj,'pos',cur_pos+sel_dims,'tag','selected_tab1',...
              'back',grey);
  eraser_pos = cur_pos+sel_dims;
  set(eraser1,'pos',[eraser_pos(1:3)+[2,0,-4],2]);

  set(ui_han(25+cur_sys),'checked','off');
  set(ui_han(25+abs(mode)),'checked','on');
  cur_sys = abs(mode);
  set(ui_han(30),'userdata',cur_sys);
  set(lab_txt(3),'callback',['dispfrac(',int2str(rem(cur_sys,3)+1),')'],...
                 'string',['System ',int2str(cur_sys),' Info']);

% restore settings for system
% Plant option
  if sys_state(1,cur_sys),
   popup_str = ['Gain|Integrator|Differentiator|Real Pole|Real Zero|',...
      'Complex Pole|Complex Zero|Lead/Lag|Num/Den|Delay Time|Design Point'];
   set(popup,'string',popup_str,'callback','termmang(0,1,1)');
   set(plant,'checked','on','enable','on');
   set(cont_type,'enable','off');
   set([for_cont,bac_cont],'checked','off');
   set(ui_obj.SysUI(1),'value',1);
   set(ui_obj.SysUI(2:3),'value',0);
   set(lab_txt(1),'callback','dispfrac(5)');
   loc_type = 1;

% set to one option
   set(set_to_one1,'label',['Set G(s) = 1']);
   set(set_to_one2,'label',['Set H(s) = 1']);
   if sys_state(4,cur_sys),
    set(set_to_one1,'checked','on');
    set(for_cont,'enable','off');
   else
    set(set_to_one1,'checked','off');
    set(for_cont,'enable','on');
   end
   if sys_state(5,cur_sys),
    set(set_to_one2,'checked','on');
    set(bac_cont,'enable','off');
   else
    set(set_to_one2,'checked','off');
    set(bac_cont,'enable','on');
   end
  end

% Forward Controller option
  if sys_state(2,cur_sys),
   set(for_cont,'checked','on','enable','on');
   set([plant,bac_cont],'checked','off');
   set(ui_obj.SysUI(2),'value',1);
   set(ui_obj.SysUI([1,3]),'value',0);
   set(cont_type,'enable','on');
   set(lab_txt(1),'callback','dispfrac(6)');

% set to one option
   set(set_to_one1,'label',['Set P(s) = 1']);
   set(set_to_one2,'label',['Set H(s) = 1']);
   if sys_state(4,cur_sys),
    set(set_to_one1,'checked','on');
    set(plant,'enable','off');
   else
    set(set_to_one1,'checked','off');
    set(plant,'enable','on');
   end
   if sys_state(5,cur_sys),
    set(set_to_one2,'checked','on');
    set(bac_cont,'enable','off');
   else
    set(set_to_one2,'checked','off');
    set(bac_cont,'enable','on');
   end

% what type of forward controller
   loc_type = find(sys_state(6:8,cur_sys)==1);
   set(ui_han(36:38),'checked','off');
   set(ui_han(35+loc_type),'checked','on');

   if loc_type == 1,
    popup_str = ['Gain|Integrator|Differentiator|Real Pole|Real Zero|',...
      'Complex Pole|Complex Zero|Lead/Lag|Num/Den|Delay Time|Design Point'];
    set(popup,'string',popup_str,'callback','termmang(0,1,1)');
   else
    popup_str = ['Proportional Gain|Integral Gain|Derivative Gain|Design Poin'];
    set(popup,'string',popup_str,'callback','pidmang(0,1)');
   end

  end

% Feedback Controller option
  if sys_state(3,cur_sys),
   set(bac_cont,'checked','on','enable','on');
   set([plant,for_cont],'checked','off');
   set(cont_type,'enable','on');
   set(ui_obj.SysUI(3),'value',1);
   set(ui_obj.SysUI(1:2),'value',0);
   set(lab_txt(1),'callback','dispfrac(6)');

% set to one option
   set(set_to_one1,'label',['Set P(s) = 1']);
   set(set_to_one2,'label',['Set G(s) = 1']);
   if sys_state(4,cur_sys),
    set(set_to_one1,'checked','on');
    set(plant,'enable','off');
   else
    set(set_to_one1,'checked','off');
    set(plant,'enable','on');
   end
   if sys_state(5,cur_sys),
    set(set_to_one2,'checked','on');
    set(for_cont,'enable','off');
   else
    set(set_to_one2,'checked','off');
    set(for_cont,'enable','on');
   end

% what type of feedback controller
   loc_type = find(sys_state(9:11,cur_sys)==1);
   set(ui_han(36:38),'checked','off');
   set(ui_han(35+loc_type),'checked','on');

   if loc_type == 1,
    popup_str = ['Gain|Integrator|Differentiator|Real Pole|Real Zero|',...
      'Complex Pole|Complex Zero|Lead/Lag|Num/Den|Delay Time|Design Point'];
    set(popup,'string',popup_str,'callback','termmang(0,1,1)');
   else
    popup_str = ['Proportional Gain|Integral Gain|Derivative Gain|Design Poin'];
    set(popup,'string',popup_str,'callback','pidmang(0,1)');
   end

  end

  if loc_type == 1,
   termmang(1,2,1);
  else
   pidmang(1,2)
  end

% clear root locus storage
  set(ui_han([37,69:75]),'userdata',[]);

  pageplot(0,13);

 elseif mode == 4, % switching to plant
  if strcmp(get(plant,'enable'),'off'),
   dispfrac(5);
   return
  end
  set(stat_bar,'string',['Switching to Plant ',int2str(cur_sys),' transfer function']);
  set(lab_txt(1),'callback','dispfrac(5)');
  set(plant,'checked','on');
  set([for_cont,bac_cont],'checked','off');
  set(ui_obj.SysUI(1),'value',1);
  set(ui_obj.SysUI(2:3),'value',0);

  set(set_to_one1,'label',['Set G(s) = 1']);
  set(set_to_one2,'label',['Set H(s) = 1']);

  if strcmp(get(for_cont,'enable'),'off'),
   set(set_to_one1,'checked','on');
   temp1 = get(ui_han(12+cur_sys*2-1),'userdata');
   temp2 = get(ui_han(12+cur_sys*2),'userdata');
   set(ui_han(12+cur_sys*2-1),'userdata',temp2);
   set(ui_han(12+cur_sys*2),'userdata',temp1);
  end

  if strcmp(get(bac_cont,'enable'),'off'),
   set(set_to_one2,'checked','on');
  end

  set(cont_type,'enable','off');
  popup_str = ['Gain|Integrator|Differentiator|Real Pole|Real Zero|',...
            'Complex Pole|Complex Zero|Lead/Lag|Num/Den|Delay Time|Design Point'];
  set(popup,'string',popup_str,'callback','termmang(0,1,1)');
  termmang(1,2,1)

% update sys_state matrix
  sys_state(2:3,cur_sys) = [0;0];
  sys_state(1,cur_sys) = 1;

% need to update gain plots because of transfer function change
  pageplot(4,[]);

 elseif mode == 5, % switching to forward controller
  if strcmp(get(for_cont,'enable'),'off'),
   dispfrac(6);
   return
  end
  set(stat_bar,'string',['Switching to Forward Controller ',int2str(cur_sys),' transfer function']);
  set(lab_txt(1),'callback','dispfrac(6)');
  set(for_cont,'checked','on');
  set([plant,bac_cont],'checked','off');
  set(ui_obj.SysUI(2),'value',1);
  set(ui_obj.SysUI([1,3]),'value',0);

  set(set_to_one1,'label',['Set P(s) = 1']);
  set(set_to_one2,'label',['Set H(s) = 1']);

  if strcmp(get(plant,'enable'),'off'),
   set(set_to_one1,'checked','on');
  end
  if strcmp(get(bac_cont,'enable'),'off'),
   set(set_to_one2,'checked','on');
  end

  set(cont_type,'enable','on');

% what type of forward controller
  loc_type = find(sys_state(6:8,cur_sys)==1);
  set(ui_han(36:38),'checked','off');
  set(ui_han(35+loc_type),'checked','on');

  if strcmp(get(num_den,'checked'),'on'),
   popup_str = ['Gain|Integrator|Differentiator|Real Pole|Real Zero|',...
     'Complex Pole|Complex Zero|Lead/Lag|Num/Den|Delay Time|Design Point'];
   set(popup,'string',popup_str,'callback','termmang(0,1,1)');
   termmang(1,2,1)
  else
   popup_str = ['Proportional Gain|Integral Gain|Derivative Gain|Design Point'];
   set(popup,'string',popup_str,'callback','pidmang(0,1)');
   pidmang(1,2)
  end

% update sys_state matrix
  sys_state([1,3],cur_sys) = [0;0];
  sys_state(2,cur_sys) = 1;

% need to update gain plots because of transfer function change
  pageplot(4,[]);

 elseif mode == 6, % switching to feedback controller
  if strcmp(get(bac_cont,'enable'),'off'),
   dispfrac(4);
   return
  end
  set(stat_bar,'string',['Switching to Feedback Controller ',int2str(cur_sys),' transfer function']);
  set(lab_txt(1),'callback','dispfrac(4)');
  set(bac_cont,'checked','on');
  set([plant,for_cont],'checked','off');
  set(ui_obj.SysUI(3),'value',1);
  set(ui_obj.SysUI(1:2),'value',0);

  set(set_to_one1,'label',['Set P(s) = 1']);
  set(set_to_one2,'label',['Set G(s) = 1']);

  if strcmp(get(plant,'enable'),'off'),
   set(set_to_one1,'checked','on');
  end
  if strcmp(get(for_cont,'enable'),'off'),
   set(set_to_one2,'checked','on');
   temp1 = get(ui_han(12+cur_sys*2-1),'userdata');
   temp2 = get(ui_han(12+cur_sys*2),'userdata');
   set(ui_han(12+cur_sys*2-1),'userdata',temp2);
   set(ui_han(12+cur_sys*2),'userdata',temp1);
  end

  set(cont_type,'enable','on');

% what type of feedback controller
  loc_type = find(sys_state(9:11,cur_sys)==1);
  set(ui_han(36:38),'checked','off');
  set(ui_han(35+loc_type),'checked','on');

  if strcmp(get(num_den,'checked'),'on'),
   popup_str = ['Gain|Integrator|Differentiator|Real Pole|Real Zero|',...
     'Complex Pole|Complex Zero|Lead/Lag|Num/Den|Delay Time|Design Point'];
   set(popup,'string',popup_str,'callback','termmang(0,1,1)');
   termmang(1,2,1)
  else
   popup_str = ['Proportional Gain|Integral Gain|Derivative Gain|Design Point'];
   set(popup,'string',popup_str,'callback','pidmang(0,1)');
   pidmang(1,2)
  end

% update sys_state matrix
  sys_state(1:2,cur_sys) = [0;0];
  sys_state(3,cur_sys) = 1;

% need to update gain plots because of transfer function change
  pageplot(4,[]);

 elseif mode == 7, % component
  set(stat_bar,'string','Displaying individual transfer function');
  set(ui_han([45:46]),'checked','off')
  set(component,'checked','on');

 elseif mode == 8, % open-loop
  set(stat_bar,'string','Displaying open-loop transfer function');
  set(ui_han([44,46]),'checked','off')
  set(open_loop,'checked','on');

 elseif mode == 9, % closed-loop
  set(stat_bar,'string','Displaying closed-loop transfer function');
  set(ui_han([44:45]),'checked','off');
  set(closed_loop,'checked','on');

 elseif mode == 10, % factored
  set(stat_bar,'string','Displaying pole-zero factored form');
  set(ui_han([48:49]),'checked','off');
  set(factored,'checked','on');

 elseif mode == 11, % coefficient
  set(stat_bar,'string','Displaying polynomial form');
  set(ui_han([47,49]),'checked','off');
  set(coefficient,'checked','on');

 elseif mode == 12, % inverse laplace transform
  set(stat_bar,'string','Displaying inverse Laplace form');
  set(ui_han([47:48]),'checked','off');
  set(inverse,'checked','on');

 end
 set(ui_han(33),'userdata',sys_state);

end

% determine states of display

% which component is displayed
if strcmp(get(plant,'checked'),'on'), mode1 = 1;
elseif strcmp(get(for_cont,'checked'),'on'), mode1 = 2;
else mode1 = 3; end

% what type of controller is displayed
if strcmp(get(num_den,'checked'),'on'), mode2 = 1;
elseif strcmp(get(pid_1,'checked'),'on'), mode2 = 2;
else mode2 = 3; end

% what configuration is desired
if strcmp(get(component,'checked'),'on'), mode3 = 1;
elseif strcmp(get(open_loop,'checked'),'on'), mode3 = 2;
else mode3 = 3; end

% what form is desired
if strcmp(get(factored,'checked'),'on'), mode4 = 1;
elseif strcmp(get(coefficient,'checked'),'on'), mode4 = 2;
else mode4 = 3; end

if mode1 == 1,
 set(lab_txt(1),'string',['Plant ',int2str(cur_sys),' Entry']);
elseif mode1 == 2,
 set(lab_txt(1),'string',['Forward ',int2str(cur_sys),' Entry']);
else
 set(lab_txt(1),'string',['Feedback ',int2str(cur_sys),' Entry']);
end

plant_mat = get(ui_han(3+cur_sys),'userdata');
for_mat = get(ui_han(6+cur_sys),'userdata');
bac_mat = get(ui_han(9+cur_sys),'userdata');

if mode3 == 1,

 if mode1 == 1,
  tf_str = ['P',int2str(cur_sys)];
  term_mat = plant_mat;
 elseif mode1 == 2,
  tf_str = ['G',int2str(cur_sys)];
  term_mat = for_mat;
 else
  tf_str = ['H',int2str(cur_sys)];
  term_mat = bac_mat;
 end

elseif mode3 == 2,

 tf_str = ['GPH',int2str(cur_sys)];
 term_mat = termjoin(plant_mat,for_mat,bac_mat);

else

 if any(mode == [4:6]), return; end
 tf_str = ['CL',int2str(cur_sys)];
 [nump,denp] = termextr(plant_mat);
 [numg,deng] = termextr(for_mat);
 [numh,denh] = termextr(bac_mat);
 neg_locus = strcmp(get(neg_loc,'checked'),'on');
 numcl = conv(conv(nump,numg),denh);
 num = conv(conv(nump,numg),numh);
 den = conv(conv(denp,deng),denh);
 ln = length(num); ld = length(den);
 num = [zeros(1,ld-ln),num];
 den = [zeros(1,ln-ld),den];
 dencl = den + num;
 term_mat = termpars(numcl,dencl);

end

if mode4 < 3 , tf_str = [tf_str,'(s)='];
else tf_str = [tf_str,'(t)=']; end

if mode1 > 1 & mode2 > 1 & mode3 == 1, % pid controller
 if mode4 > 1,
  if term_mat(4) == 40, term_mat = pid2pid(term_mat,1); end
  term_mat = termpars(term_mat([3,1,2]),[1,0]);
 end
end

% destroy present transfer function
axes(term_axis);
cla;

line_yloc = 0.5;
T = [];
txt = text(0,0.5,tf_str,'erase',ers,'horizontalalignment','left',...
        'verticalalignment','middle','fontname',fontname,'color','k');
tf_ext = get(txt,'extent');

% initial locations of text for the numerator and denominator
num_btm = line_yloc+0.01;
den_btm = line_yloc-tf_ext(4)-0.01;
txt_lft = sum(tf_ext([1,3]))+0.01;
term_val = [1:8];

% mode1 - transfer function
% mode2 - controller type
% mode3 - component, open-loop, closed-loop
% mode4 - factored, coefficient, inverse laplace

if mode4 == 1,
 if mode1 > 1 & mode2 == 2 & mode3 == 1, % pid form #1
  p_gain = num2str(term_mat(1));
  i_gain = num2str(term_mat(2));
  d_gain = num2str(term_mat(3));
  txt = text(txt_lft,0.5,[p_gain,' + '],'erase',ers,...
            'horizontalalignment','left','verticalalignment','middle',...
            'userdata',[1,0],'fontname',fontname,...
            'tag','0','color','k');
%            'buttondownfcn','pagemous(0)');
  txt_ext = get(txt,'extent');
  txt_lft = txt_lft+txt_ext(3);
  txt = text(txt_lft,num_btm,i_gain,'erase',ers,...
            'horizontalalignment','left','verticalalignment','bottom',...
            'userdata',[2,0],'fontname',fontname,...
            'tag','0','color','k');
%            'buttondownfcn','pagemous(0)');
  txt_ext = get(txt,'extent');
  new_lft = txt_lft+txt_ext(3);
  lin_lft = txt_lft; lin_rht = new_lft;
  lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');
  txt_mid = txt_lft + (lin_rht-lin_lft)/2;
  txt = text(txt_mid,den_btm,'s','erase',ers,...
           'horizontalalignment','center','verticalalignment','bottom',...
           'userdata',[2,0],'fontname',fontname,...
           'tag','0','color','k');
%           'buttondownfcn','pagemous(0)');
  text(new_lft,0.5,[' + ',d_gain,'s'],'erase',ers,...
            'horizontalalignment','left','verticalalignment','middle',...
            'userdata',[3,0],'fontname',fontname,...
            'tag','0','color','k');
%            'buttondownfcn','pagemous(0)');

 elseif mode1 > 1 & mode2 == 3 & mode3 == 1, % pid form #2

  p_gain = num2str(term_mat(1));
  i_gain = num2str(abs(term_mat(2)));
  d_gain = num2str(abs(term_mat(3)));

  if term_mat(2) >= 0, p_gain = [p_gain,' + '];
  else p_gain = [p_gain,' - ']; end

  if term_mat(3) >= 0, d_gain = [' + ',d_gain];
  else d_gain = [' - ',d_gain]; end

  txt = text(txt_lft,0.5,['(1',d_gain,'s)'],'erase',ers,...
            'horizontalalignment','left','verticalalignment','middle',...
            'userdata',[3,0],'fontname',fontname,...
            'tag','0','color','k');
%            'buttondownfcn','pagemous(0)');
  txt_ext = get(txt,'extent');
  txt_lft = txt_lft+txt_ext(3);

  txt = text(txt_lft,0.5,['(',p_gain],'erase',ers,...
            'horizontalalignment','left','verticalalignment','middle',...
            'userdata',[1,0],'fontname',fontname,...
            'tag','0','color','k');
%            'buttondownfcn','pagemous(0)');
  txt_ext = get(txt,'extent');
  txt_lft = txt_lft+txt_ext(3);

  txt = text(txt_lft,num_btm,i_gain,'erase',ers,...
            'horizontalalignment','left','verticalalignment','bottom',...
            'userdata',[2,0],'fontname',fontname,...
            'tag','0','color','k');
%            'buttondownfcn','pagemous(0)');
  txt_ext = get(txt,'extent');
  new_lft = txt_lft+txt_ext(3);
  lin_lft = txt_lft; lin_rht = new_lft;
  lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');
  txt_mid = txt_lft + (lin_rht-lin_lft)/2;
  txt = text(txt_mid,den_btm,'s','erase',ers,...
           'horizontalalignment','center','verticalalignment','bottom',...
           'userdata',[2,0],'fontname',fontname,...
           'tag','0','color','k');
%           'buttondownfcn','pagemous(0)');
  text(new_lft,0.5,')','erase',ers,...
            'horizontalalignment','left','verticalalignment','middle',...
            'userdata',[2,0],'fontname',fontname,...
            'tag','0','color','k');

 else
  for k = term_val,
   term_loc=find(k==term_mat(:,4));
   if length(term_loc),
    for h=1:length(term_loc),

     if k==1, % Gain

      gain = num2str(term_mat(1,1));
      txt = text(txt_lft,0.5,gain,'erase',ers,...
                 'horizontalalignment','left','verticalalignment','middle',...
                 'fontname',fontname,'color','k');
      txt_ext = get(txt,'extent');

      if mode3 == 1 & mode4 == 1,
       set(txt,'userdata',term_loc(h),'tag','0');
%       set(txt,'userdata',term_loc(h),'buttondownfcn','pagemous(0)');
      else
       set(txt,'tag','9');
      end

      txt_lft = txt_lft + txt_ext(3)+0.01;

     elseif k==2,

      integ = term_mat(term_loc(h),1);
      if integ~=0,
       if integ > 0,
        txt1_btm = den_btm; txt2_btm = num_btm;
       else
        txt1_btm = num_btm; txt2_btm = den_btm;
       end

       coeffs = [1,zeros(1,abs(integ))];
       [txt_han,new_lft] = coefdisp(txt_lft,txt1_btm,coeffs,0,term_loc(h));

       lin_lft = txt_lft; lin_rht = new_lft;
       lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');
       txt_mid = txt_lft + (lin_rht-lin_lft)/2;
       txt = text(txt_mid,txt2_btm,'1','erase',ers,...
            'horizontalalignment','center','verticalalignment','bottom',...
            'fontname',fontname,'color','k');
       set(txt,'userdata',term_loc(h));
       set(lin,'userdata',[txt_han,txt]);

       if mode3 == 1 & mode4 == 1,
        set([txt,txt_han],'tag','0');
%        set([txt,txt_han],'buttondownfcn','pagemous(0)');
       else
        set([txt,txt_han],'tag','9');
       end

       txt_lft = new_lft + 0.01;
      end

     elseif any(k==[4,5]), % real pole/zero

      if k==4,
       txt1_btm = den_btm; txt2_btm = num_btm;
      else
       txt1_btm = num_btm; txt2_btm = den_btm;
      end

      coeffs = [1,term_mat(term_loc(h),1)];
      [txt_han,new_lft] = coefdisp(txt_lft,txt1_btm,coeffs,1,term_loc(h));

      lin_lft = txt_lft; lin_rht = new_lft;
      lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');
      txt_mid = txt_lft + (lin_rht-lin_lft)/2;
      txt = text(txt_mid,txt2_btm,'1','erase',ers,...
              'horizontalalignment','center','verticalalignment','bottom',...
              'fontname',fontname,'color','k');
      set(txt,'userdata',term_loc(h));

      if mode3 == 1 & mode4 == 1,
       set([txt,txt_han],'tag','0');
%       set([txt,txt_han],'buttondownfcn','pagemous(0)');
      else
       set([txt,txt_han],'tag','9');
      end

      set(lin,'userdata',[txt_han,txt]);
      txt_lft = new_lft + 0.01;

     elseif any(k==[6,7]), % complex pole/zero

      if k==6,
       txt1_btm = den_btm; txt2_btm = num_btm;
      else
       txt1_btm = num_btm; txt2_btm = den_btm;
      end

      zeta = term_mat(term_loc(h),1);
      wn = term_mat(term_loc(h),2);
      coeffs = [1,2*zeta*wn,wn^2];
      [txt_han,new_lft] = coefdisp(txt_lft,txt1_btm,coeffs,1,term_loc(h));

      lin_lft = txt_lft; lin_rht = new_lft;
      lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');
      txt_mid = txt_lft + (lin_rht-lin_lft)/2;
      txt = text(txt_mid,txt2_btm,'1','erase',ers,...
             'horizontalalignment','center','verticalalignment','bottom',...
             'fontname',fontname,'color','k');
      set(txt,'userdata',term_loc(h));

      if mode3 == 1 & mode4 == 1,
       set([txt,txt_han],'tag','0');
%       set([txt,txt_han],'buttondownfcn','pagemous(0)');
      else
       set([txt,txt_han],'tag','9');
      end

      set(lin,'userdata',[txt_han,txt]);
      txt_lft = new_lft + 0.01;

     elseif k==8, % lead/lag

      phase = term_mat(term_loc(h),1);
      freq = term_mat(term_loc(h),2);
      [jk,zero,pole]=leadlag(phase,freq,[],T);

      zero_str = num2str(zero);
      pole_str = num2str(pole);

      if zero > 0,
       lead = ['(s+',zero_str,')'];
      else
       lead = ['(s',zero_str,')'];
      end

      if pole > 0,
       lag = ['(s+',pole_str,')'];
      else
       lag = ['(s',pole_str,')'];
      end

      txt1 = text(txt_lft,num_btm,lead,'erase',ers,...
            'horizontalalignment','left','verticalalignment','bottom',...
            'fontname',fontname,'color','k');
      set(txt1,'userdata',term_loc(h));

      txt_ext = get(txt1,'extent');
      lin_lft = txt_lft; lin_rht = txt_lft+txt_ext(3);
      lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');

      txt2 = text(txt_lft,den_btm,lag,'erase',ers,...
            'horizontalalignment','left','verticalalignment','bottom',...
            'fontname',fontname,'color','k');
      set(txt2,'userdata',term_loc(h));
      set(lin,'userdata',[txt1,txt2]);

      if mode3 == 1 & mode4 == 1,
       set([txt1,txt2],'tag','0');
%       set([txt1,txt2],'buttondownfcn','pagemous(0)');
      else
       set([txt1,txt2],'tag','9');
      end

      txt_lft = txt_lft + txt_ext(3) + 0.01;

     end % end else-ifs
    end % end for h = loop
   end % end of if length(term_loc)
  end % end of for k = loop

  if term_mat(1,2)~=0,
   e_str = 'e';
   txt1=text(txt_lft,line_yloc,e_str,'erase',ers,...
             'horizontalalignment','left','verticalalignment','middle',...
             'fontname',fontname,'color','k');
   ext_txt = get(txt1,'extent');
   txt_lft = txt_lft + ext_txt(3);
   if term_mat(1,2) > 0,
    delay_str = ['-',num2str(term_mat(1,2)),'s'];
   else
    delay_str = [num2str(-term_mat(1,2)),'s'];
   end
   txt2=text(txt_lft,line_yloc,delay_str,'erase',ers,...
             'horizontalalignment','left','verticalalignment','bottom',...
             'fontname',fontname,'color','k');
   if mode3 == 1 & mode4 == 1,
    set([txt1,txt2],'userdata',[0,1],'tag','0');
%       set([txt1,txt2],'buttondownfcn','pagemous(0)');
   else
    set([txt1,txt2],'tag','9');
   end
  end
 end

elseif mode4 == 2, % coefficient

 [num,den] = termextr(term_mat);
 num = num/term_mat(1,1);

 gain = num2str(term_mat(1,1));
 txt = text(txt_lft,0.5,gain,'erase',ers,...
            'horizontalalignment','left','verticalalignment','middle',...
            'fontname',fontname,'color','k');
 txt_ext = get(txt,'extent');
 set(txt,'userdata',1);
 txt_lft = txt_lft + txt_ext(3)+0.01;

 [txt_han1,new_lft1] = coefdisp(txt_lft,num_btm,num,0,[0,0]);
 [txt_han2,new_lft2] = coefdisp(txt_lft,den_btm,den,0,[0,0]);
 min_lft = min(new_lft1,new_lft2); max_lft = max(new_lft1,new_lft2);
 shft = (max_lft - min_lft)/2;
 if new_lft1 < new_lft2,
  for k = length(txt_han1):-1:1,
   txt_pos = get(txt_han1(k),'pos');
   set(txt_han1(k),'pos',txt_pos+[shft,0,0]);
  end
 elseif new_lft2 < new_lft1,
  for k = length(txt_han2):-1:1,
   txt_pos = get(txt_han2(k),'pos');
   set(txt_han2(k),'pos',txt_pos+[shft,0,0]);
  end
 end
 lin_lft = txt_lft; lin_rht = max_lft;
 lin = line([lin_lft,lin_rht],[line_yloc,line_yloc],'erase',ers,'color','k');
 set(lin,'userdata',[txt_han1,txt_han2]);
 set([txt_han1,txt_han2],'tag','9');

elseif mode4 == 3, % Inverse Laplace Transform

 residisp(term_mat,tf_ext);

end


text:
function [] = showSuperquadrics(x, varargin)

R = eul2rotm(x(6 : 8));
t = x(9 : 11);

% with tapering or not
taper = false;
color = 'r';
ViewAxis = [0 0];
CamRoll = 0;
ShowAxis = 0;
arclength = 0.02;
FaceAlpha = 1;
FaceLighting = 'flat';
lighting = false;

for k = 1 : size(varargin, 2)
    if strcmp(varargin{k}, 'Taper')
        taper = varargin{k + 1};
    end
    if strcmp(varargin{k}, 'Color')
        color = varargin{k + 1};
    end
    if strcmp(varargin{k}, 'ViewAxis')
        ViewAxis = varargin{k + 1};
    end
    if strcmp(varargin{k}, 'CamRoll')
        CamRoll = varargin{k + 1};
    end
    if strcmp(varargin{k}, 'ShowAxis')
        ShowAxis= varargin{k + 1};
    end
    if strcmp(varargin{k}, 'Arclength')
        arclength= varargin{k + 1};
    end
    if strcmp(varargin{k}, 'FaceAlpha')
        FaceAlpha= varargin{k + 1};
    end
    if strcmp(varargin{k}, 'FaceLighting')
        FaceLighting= varargin{k + 1};
    end
    if strcmp(varargin{k}, 'Light')
        lighting= varargin{k + 1};
    end
end

% validate dimensionality
if taper == true
    if size(x, 2) ~= 13
        error('Input parameters should have dimension (:, 13) for taperred SQ.')
    end
else
    if size(x, 2) ~= 11
        error('Input parameters should have dimension (:, 11) for taperred SQ.')
    end
end

% avoiding numerical instability of points sampling on superquadrics
if x(1) < 0.01 %0.007
    x(1) = 0.01;
end
if x(2) < 0.01
    x(2) = 0.01;
end

[point_eta] = uniformSampledSuperellipse(x(1), [1, x(5)], arclength);
[point_omega] = uniformSampledSuperellipse(x(2), [x(3), x(4)], arclength);

x_mesh = ones(size(point_omega, 2), size(point_eta, 2));
y_mesh = ones(size(point_omega, 2), size(point_eta, 2));
z_mesh = ones(size(point_omega, 2), size(point_eta, 2));

for m = 1 : size(point_omega, 2)
    for n = 1 : size(point_eta, 2)
        point_temp = [point_omega(:, m) * point_eta(1, n); point_eta(2, n)];

        if taper == true
            fx = x(12) * point_temp(3) / x(5) + 1;
            fy = x(13) * point_temp(3) / x(5) + 1;
            fz = 1;

            point_temp(1) = point_temp(1) * fx;
            point_temp(2) = point_temp(2) * fy;
            point_temp(3) = point_temp(3) * fz;
        end

        point_temp = R * point_temp + t';

        x_mesh(m, n) = point_temp(1);
        y_mesh(m, n) = point_temp(2);
        z_mesh(m, n) = point_temp(3);
    end
end

mesh(x_mesh, y_mesh, z_mesh, 'FaceAlpha', FaceAlpha, 'facecolor', color, ...
    'LineStyle', 'none', 'FaceLighting', FaceLighting)
if lighting == 1
    light
    material dull
end

axis equal
view(ViewAxis)
camroll(CamRoll)

if ShowAxis == 0
    axis off
end

hold off
% ---------------------------------utility functions ----------------------
    function [point, theta] = uniformSampledSuperellipse(epsilon, scale, arclength)
        threshold = 1e-2;
        num_limit = 10000;
        theta = zeros(1, num_limit);
        theta(1) = 0;

        for i = 2 : num_limit
            dt = dtheta(theta(i - 1), arclength, threshold, scale, epsilon);
            theta_temp = theta(i - 1) + dt;

            if theta_temp > pi/4
                break
            else
                if i < num_limit
                    theta(i) = theta_temp;
                else
                    error(['The number of the sampled points exceeds the limit of ', ...
                        num2str(num_limit * 4),...
                        '. Please increase the arclength or raise the limit'])
                end
            end
        end
        critical = i;

        for j = critical + 1 : num_limit
            dt = dtheta(theta(j - 1), arclength, threshold, flip(scale), epsilon);
            theta_temp = theta(j - 1) + dt;

            if theta_temp > pi/4
                break
            else
                if j < num_limit
                    theta(j) = theta_temp;
                else
                    error(['The number of the sampled points exceeds the limit of ', ...
                        num2str(num_limit * 4),...
                        '. Please increase the arclength or raise the limit'])
                end
            end
        end

        num_pt = j - 1;
        theta = theta(1 : num_pt);

        points_fw = angle2points(theta(1 : critical - 1), scale, epsilon);
        points_bw = flip(angle2points(theta(critical : end), flip(scale), epsilon), 2);
        point = [points_fw, [points_bw(2, :); points_bw(1, :)]];

        point = [point, flip([-point(1, 1 : num_pt - 1); point(2, 1 : num_pt - 1)], 2), ...
            [-point(1, 2 : end); -point(2, 2 : end)], flip([point(1, 1 : num_pt - 1); ...
            -point(2, 1 : num_pt - 1)], 2)];

    end

    function [dt] = dtheta(theta, arclength, threshold, scale, sigma)
        if theta < threshold
            dt = abs((arclength / scale(2) + (theta)^(sigma))^(1 / sigma) ...
                - (theta));
        else
            dt = arclength / sigma * ((cos(theta) ^ 2 * sin(theta) ^ 2) / ...
                (scale(1) ^ 2 * cos(theta) ^ (2 * sigma) * sin(theta) ^ 4 + ...
                scale(2) ^ 2 * sin(theta) ^ (2 * sigma) * cos(theta) ^ 4))^(1 / 2);
        end
    end

    function [point] = angle2points(theta, scale, sigma)
        point = zeros(2, size(theta, 2));
        point(1, :) = scale(1) .* sign(cos(theta)) .* abs(cos(theta)).^sigma;
        point(2, :) = scale(2) .* sign(sin(theta)) .* abs(sin(theta)).^sigma;
    end

end

author: ChirikjianLab

repo: Marching-Primitives


text:
%input is defined in "b", ref is defined in "n" (e.g inp=[g*sin(45) g*cos(45) 0], ref=[0 0 1])
%Cbn_est=(I-S(e))Cbn, e=E[inp_err]

function [Cnb E]=align_wander(inp, ref)

%Compute the dcm matrix
Cnb=align_vector(inp,ref);

if (nargout==2) %Compute the error matrix
    E=zeros(3);
    ep=0.001;

    for i=1:3
        %ith perturbation
        inp_err=inp;
        inp_err(i)=inp_err(i)+ep;
        Cnb_err=align_vector(inp_err,ref);
        err=Cnb_err*Cnb'-eye(3);
        E(:,i)=[-err(2,3);err(1,3);-err(1,2)]/ep;
    end

end

function Cnb=align_vector(inp, ref)
%normalize
inp=inp/norm(inp);
ref=ref/norm(ref);

%angles
dotp=dot(ref,inp);
crosv=cross(ref,inp);
crosp=norm(crosv);

if (crosp==0) %inp and ref is already aligned
    if (dotp>=0) %no rotation is necessary assert(dotp==1)
        Cnb=eye(3);
    else %assert(dotp==-1)
        %Find a vector that is ortogonal to ref
        if (ref(1)^2>ref(2)^2)
            vr_a=[-ref(3);0;ref(1)];
        else
            vr_a=[0;-ref(3);ref(2)];
        end
        vr_a=vr_a/norm(vr_a)*pi;

        %rotate around computed orthogonal vector
        Cnb=rot2dcm_v000(vr_a);
    end
else
    Cnb=eye(3)+skew(crosv)+((1-dotp)/crosp^2)*skew(crosv)*skew(crosv);
end

return;

repo: ekfmonoslam


text:
function result = spsym (A, quick)					    %#ok
%SPSYM determine if a sparse matrix is symmetric, Hermitian, or skew-symmetric.
%   If so, also determine if its diagonal has all positive real entries.
%   A must be sparse.
%
%   Example:
%   result = spsym (A) ;
%   result = spsym (A,quick) ;
%
%   If quick = 0, or is not present, then this routine returns:
%
%       1: if A is rectangular
%       2: if A is unsymmetric
%       3: if A is symmetric, but with one or more A(j,j) <= 0
%       4: if A is Hermitian, but with one or more A(j,j) <= 0 or with
%           nonzero imaginary part
%       5: if A is skew symmetric (and thus the diagonal is all zero as well)
%       6: if A is symmetric with real positive diagonal
%       7: if A is Hermitian with real positive diagonal
%
%   If quick is nonzero, then the function can return more quickly, as soon as
%   it finds a diagonal entry that is <= 0 or with a nonzero imaginary part.
%   In this case, it returns 2 for a square matrix, even if the matrix might
%   otherwise be symmetric or Hermitian.
%
%   Regardless of the value of "quick", this function returns 6 or 7 if A is
%   a candidate for sparse Cholesky.
%
%   For an MATLAB M-file function that computes the same thing as this
%   mexFunction (but much slower), see the get_symmetry function by typing
%   "type spsym".
%
%   This spsym function does not compute the transpose of A, nor does it need
%   to examine the entire matrix if it is unsymmetric.  It uses very little
%   memory as well (just size-n workspace, where n = size (A,1)).
%
%   Examples:
%       load west0479
%       A = west0479 ;
%       spsym (A)
%       spsym (A+A')
%       spsym (A-A')
%       spsym (A+A'+3*speye(size(A,1)))
%
%   See also mldivide.

%       function result = get_symmetry (A,quick)
%       %GET_SYMMETRY: does the same thing as the spsym mexFunction.
%       % It's just a lot slower and uses much more memory.  This function
%       % is meant for testing and documentation only.
%       [m n] = size (A) ;
%       if (m ~= n)
%           result = 1 ;            % rectangular
%           return
%       end
%       if (nargin < 2)
%           quick = 0 ;
%       end
%       d = diag (A) ;
%       posdiag = all (real (d) > 0) & all (imag (d) == 0) ;
%       if (quick & ~posdiag)
%           result = 2 ;            % Not a candidate for sparse Cholesky.
%       elseif (~isreal (A) & nnz (A-A') == 0)
%           if (posdiag)
%               result = 7 ;        % complex Hermitian, with positive diagonal
%           else
%               result = 4 ;        % complex Hermitian, nonpositive diagonal
%           end
%       elseif (nnz (A-A.') == 0)
%           if (posdiag)
%               result = 6 ;        % symmetric with positive diagonal
%           else
%               result = 3 ;        % symmetric, nonpositive diagonal
%           end
%       elseif (nnz (A+A.') == 0)
%           result = 5 ;            % skew symmetric
%       else
%           result = 2 ;            % unsymmetric
%       end

% With additional outputs, spsym computes the following for square matrices:
% (in this case "quick" is ignored, and set to zero):
%
% [result xmatched pmatched nzoffdiag nnzdiag] = spsym(A)
%
%   xmatched is the number of nonzero entries for which A(i,j) = conj(A(j,i)).
%   pmatched is the number of entries (i,j) for which A(i,j) and A(j,i) are
%   both in the pattern of A (the value doesn't matter).  nzoffdiag is the
%   total number of off-diagonal entries in the pattern.  nzdiag is the number
%   of diagonal entries in the pattern.  If the matrix is rectangular,
%   xmatched, pmatched, nzoffdiag, and nzdiag are not computed (all of them are
%   returned as zero).  Note that a matched pair, A(i,j) and A(j,i) for i != j,
%   is counted twice (once per entry).

%   Copyright 2006-2007, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

error ('spsym mexFunction not found') ;


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function r8vec_indicator1_test ( )

%*****************************************************************************80
%
%% R8VEC_INDICATOR1_TEST tests R8VEC_INDICATOR1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_INDICATOR1_TEST\n' );
  fprintf ( 1, '  R8VEC_INDICATOR1 returns a 1-based indicator vector.\n' );

  n = 10;
  a = r8vec_indicator1 ( n );
  r8vec_print ( n, a, '  The indicator1 vector:' );

  return
end