text:
%kGenRecMat 'generates a matrix, counting peaks from left to right '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros GenRecMat.pane file
%
% Parameters:
% Integer: sx 'Pattern steps X', default: 5: 'number of equally devided steps the pattern was shifted along X'
% Integer: sy 'Pattern steps Y', default: 1: 'number of equally devided steps the pattern was shifted along Y'
% Integer: nx 'Recon Orders X', default: 1: 'number of orders to be included (in + and - direction)'
% Integer: ny 'Recon Orders Y', default: 0: 'number of orders to be included (in + and - direction)'
% Toggle: diagscan 'Diagonal Scan', default: 0: 'For scans without the zero-order, this has to be active'
% Toggle: nozeroorder 'Zero Order Blocked', default: 0: 'If selected, the zero order is assumed to not be present. One full pattern movement corresponds now to two unit cells.'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kGenRecMat( {'sx',5;'sy',1;'nx',1;'ny',0;'diagscan',0;'nozeroorder',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% GenRecMat - generates a matrix, counting peaks from left to right
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kGenRecMat(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kGenRecMat(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'sx', 5;'sy', 1;'nx', 1;'ny', 0;'diagscan', 0;'nozeroorder', 0;'o', '__output'};
maxval={2,2,1,1,0,0,0};
minval={2,2,1,1,0,0,0};
istoggle=[0,0,0,0,1,1,0];
was_set=istoggle * 0;
paramtype={'Integer','Integer','Integer','Integer','Toggle','Toggle','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'genrecmat"  -k'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef SLERP < handle

    properties (Access = public)
        tau
    end

    properties (Access = private)
        phi
        theta
        scalar_product
    end

    methods (Access = public)

        function obj = SLERP(cParams)
            obj.init(cParams);
        end

        function x = update(obj,g,~)
            obj.computeTheta(g);
            x = obj.computeNewLevelSet(g);
        end

        function computeFirstStepLength(obj,~,~,~)
            obj.tau = 1;
        end

        function is = isTooSmall(obj)
            is = obj.tau < 1e-10;
        end

        function increaseStepLength(obj,f)
            obj.tau = min(f*obj.tau,1);
        end

        function decreaseStepLength(obj)
            obj.tau = obj.tau/1.5;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.phi            = cParams.designVar;
            obj.scalar_product = cParams.uncOptimizerSettings.scalarProductSettings;
        end

        function computeTheta(obj,g)
            pN   = obj.normalizeFunction(obj.phi.value);
            gN   = obj.normalizeFunction(g);
            phiG = obj.scalar_product.computeSP(pN,gN);
            obj.theta = max(acos(phiG),1e-14);
        end

        function p = computeNewLevelSet(obj,g)
            k  = obj.tau;
            t  = obj.theta;
            pN = obj.normalizeFunction(obj.phi.value);
            gN = obj.normalizeFunction(g);
            a  = sin((1-k)*t)*pN;
            b  = sin(k*t)*gN;
            p  = (a + b)/sin(t);
            p  = obj.normalizeFunction(p);
        end

        function x = normalizeFunction(obj,x)
            norm2 = obj.scalar_product.computeSP(x,x);
            xNorm = sqrt(norm2);
            x = x/xNorm;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function X = sigm(P)
    X = 1./(1+exp(-P));
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function [GAmp,GTime]=GzCartesian(p)

global VCtl;
global VObj;
global VVar;

t1Start=p.t1Start; %ms
t1End=p.t1End; %ms
t2Start=p.t2Start;
t2End=p.t2End;
tRamp=p.tRamp;
Gz1Sign=p.Gz1Sign;
Gz2Sign=p.Gz2Sign;
% GzOrder=p.GzOrder;
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

% if strcmp(GzOrder,'centric')
%     if VVar.SliceCount==1 & VVar.PhaseCount==1
%        SliceCount = (-1)^rem(VVar.SliceCount-1,2)*(VVar.SliceCount-1);
%        VVar.SliceCountTmp = 0;
%     elseif VVar.PhaseCount~=1
%        SliceCount = VVar.SliceCountTmp;
%     else
%        SliceCount = (-1)^rem(VVar.SliceCount-1,2)*(VVar.SliceCount-1) + VVar.SliceCountTmp;
%        VVar.SliceCountTmp = SliceCount;
%     end
% elseif strcmp(GzOrder,'sequential')
%     SliceCount = VVar.SliceCount-floor(VCtl.SliceNum/2)-1;
% end

% phase encoding
GzdG=1/((VObj.Gyro/(2*pi))*(t1End-t1Start)*VCtl.FOVSlice);
[GAmp1,GTime1]=StdTrap(t1Start-tRamp, ...
                       t1End+tRamp,   ...
                       t1Start,               ...
                       t1End,                 ...
                       (VVar.SliceCount-floor(VCtl.SliceNum/2)-1)*GzdG*Gz1Sign,2,2,2); % use floor to make Kz = 0 when VCtl.SliceNum=1
% rephasing
GzdG=1/((VObj.Gyro/(2*pi))*(t2End-t2Start)*VCtl.FOVSlice);
[GAmp2,GTime2]=StdTrap(t2Start-tRamp, ...
                       t2End+tRamp,   ...
                       t2Start,               ...
                       t2End,                 ...
                       (VVar.SliceCount-floor(VCtl.SliceNum/2)-1)*GzdG*Gz2Sign,2,2,2);

GAmp=GAmp1;
GTime=GTime1;

if Gz2Sign~=0
    GAmp=[GAmp, GAmp2];
    GTime=[GTime, GTime2];
end

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function result = square_unit_sum ( func, norder, xtab, ytab, weight )

%*****************************************************************************80
%
%% SQUARE_UNIT_SUM carries out a quadrature rule over the unit square.
%
%  Integration region:
%
%    Points (X,Y) such that
%
%      -1 <= X <= 1, and
%      -1 <= Y <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, external FUNC, the name of the function to be
%    integrated.  The user must declare the name an EXTERNAL
%    parameter in the calling program, pass the name of the
%    function in FUNC, and write a function of the form
%      function value = func ( x, y )
%    which evaluates the function at the point (X,Y).
%
%    Input, integer NORDER, the order of the rule.
%
%    Input, real XTAB(NORDER), YTAB(NORDER), the abscissas of
%    the rule.
%
%    Input, real WEIGHT(NORDER), the weights of the rule.
%
%    Output, real RESULT, the approximate integral of the function.
%
  quad = 0.0E+00;
  for i = 1 : norder
    quad = quad + weight(i) * feval ( func, xtab(i), ytab(i) ) / 4.0E+00;
  end

  volume = 1.0E+00;
  result = quad * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% GRIFFIN_LIM The Griffin & Lim phase retrieval algorithm
%
% Usage
%    x = griffin_lim(x_init, x_phi, x_psi_mod, filters, options)
%
% Input
%    x_init (numeric): The inital guess for x.
%    x_phi (numeric): The lowpass coefficients
%    x_psi_mod (cell): The wavelet modulus coefficients.
%    filters (struct): The filter bank used to calculate the wavelet transform
%        in WAVELET_1D.
%    options (struct, optional): Different parameters to the algorithm and to
%        WAVELET_1D and INVERSE_WAVELET_1D that are called. Parameters for
%        GRIFFIN_LIM are:
%            gl_iter (numeric): The number of iterations (default 32).
%            verbose (boolean): If true, shows computation information (default
%               false).
%            x_phi_resolution (numeric): The resolution of the x_phi input,
%               with respect to the maximum resolution allowed by filters
%               (default 0).
%
% Output
%     x (numeric): The result of the algorithm.
%
% Description
%    Given the lowpass coefficients and the modulus of the wavelet
%    coefficients, the Griffin & Lim algorithm attempts to reconstruct the
%    original signal by recovering the phase of the wavelet coefficients. To
%    do this, it performs an alternating projection on the wavelet reproducing
%    kernel (the set of coefficients that are valid wavelet transforms) and
%    the set of coefficients of the desired modulus. Since the latter is a
%    non-convex set, the algorithm will not converge. In some applications,
%    however, the approximation that it provides is sufficient. For more
%    details, see [1].
%
% References
%    [1] D. W. Grifﬁn and J. S. Lim, “Signal estimation from modiﬁed short-
%        time fourier transform,” IEEE Trans. Acoust., Speech, Signal
%        Process., vol. 32, no. 2, pp. 236–243, 1984.
%
% See also
%    WAVELET_1D, INVERSE_WAVELET_1D

function x = griffin_lim(x_init, x_phi, x_psi_mod, filters, options)
	if nargin < 5
		options = struct();
	end

	options = fill_struct(options, 'gl_iter', 32);
	options = fill_struct(options, 'verbose', 0);
	options = fill_struct(options, 'x_phi_resolution', 0);

	dual_filters = dual_filter_bank(filters);

	x0 = x_init;
	[x0_phi, x0_psi, meta_phi, meta_psi] = wavelet_1d(x0, filters, options);

	for k = 1:length(filters.psi.filter)
		if isempty(x_psi_mod{k})
			continue;
		end
		ds = length(x_psi_mod{k})/length(x0_psi{k});
		if ds > 1
			x_psi_mod{k} = x_psi_mod{k}(1:ds:end)*sqrt(ds);
		else
			x_psi_mod{k} = upsample(x_psi_mod{k}, length(x0_psi{k}));
		end
	end

	for iter = 1:options.gl_iter
		if options.verbose, fprintf('%d...',iter); end
		for k = 1:length(filters.psi.filter)
			if isempty(x_psi_mod{k})
				x1_psi{k} = x0_psi{k};
			else
				x1_psi{k} = x0_psi{k}.*x_psi_mod{k}./abs(x0_psi{k});
			end
		end

		meta_phi.resolution = options.x_phi_resolution;

		x1 = inverse_wavelet_1d(length(x_init), x_phi, x1_psi, meta_phi, ...
			meta_psi, dual_filters, options);
		[x0_phi, x0_psi, meta_phi, meta_psi] = wavelet_1d(x1, filters, ...
			options);
	end
	if options.verbose, fprintf('\n'); end

	x = x1;
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% This demo compares dynamic Bayesian GPLVM with and without image pyramids
% to see if we get a better reconstruction of the dog tail motion using
% pyramids.
%
% Created by Andreas and Teo on 31 July 2013, based on demosDynamics.m

%------------ GENERATION -------
% Train

%% Setting paths
addpath(genpath('../../../GPmat/matlab/')); % GPmat
% addpath(genpath('../../../vargplvm/vargplvm/matlab/')); % Bayesian GPLVM
addpath('../../../GPmatlab/netlab3_3/'); % NetLab, which is not included in gitHub

%% Preparing pyramid dataset:
% This only needs to be done once.
load dog
% To view it:
% for f=1:600, image(reshape(Y(mod(f-1,61)+1,:),[360 640])); colormap gray(256); pause(0.033); end
levels = 5;
Y = im2pyramid(Y, height, width, levels);
save('dogPyramid', 'Y', 'height', 'width', 'levels');
clear

%% Selecting dataset -- To run experiments, change the dataset below and run the test again:
% function doAll(dataSetName)
dataSetName = 'dog';
% Running the rest for each dataset.
fprintf(1,'\n\n#-----  DOG DEMO: Generation ----#\n');
experimentNo=61;
indPoints=-1;
latentDim=35;
fixedBetaIters= 0;
initVardistIters = 800;
mappingKern = 'rbfardjit';
reconstrIters = 1; % no reconstruction needed here
itNo=repmat(1000, 1, 16); %16000
periodicPeriod = 4.3983; % Recalculated for dataToKeep=60
dynamicKern={'rbfperiodic','whitefixed','bias','rbf'};
vardistCovarsMult=0.8;
whiteVar = 1e-6;
dataToKeep = 60;
dataSetSplit = 'custom';
indTr = [1:60];
indTs = 60; % We don't really reconstruct in this experiment
learnSecondVariance = 0;
DgtN = 1; % Uncomment for the faster version
predWithMs = 0;
%%
demHighDimVargplvm3

%%
%-- Then generate:
% experimentNo=61;
% dataToKeep = 60;
% dataSetSplit = 'custom';
% indTr = [1:60];
% indTs = 60;
futurePred = 40;
doSampling = 0;
demHighDimVargplvmTrained
%clear Yts; clear YtsOriginal; clear Testmeans2; clear Testcovars2;
%playMov(height, width, [], [Ytr(end-5:end,:); Varmu2]);

%%

% Produce plots
bar(prunedModelInit.kern.comp{1}.inputScales)
print -depsc ../diagrams/dog_scalesInit.eps; system('epstopdf ../diagrams/dog_scalesInit.eps');
bar(model.kern.comp{1}.inputScales)
print -depsc ../diagrams/dog_scalesOpt.eps; system('epstopdf ../diagrams/dog_scalesOpt.eps');

nPixels = height*width; % This will be used in case Y was built using pyramids.
fr=reshape(Ytr(end,1:nPixels),height,width); imagesc(fr); colormap('gray'); % Last training image
print -depsc ../diagrams/dogGeneration_lastOfTraining.eps; system('epstopdf ../diagrams/dogGeneration_lastOfTraining.eps');
fr=reshape(Varmu2(1,1:nPixels),height,width); imagesc(fr); colormap('gray');  % First predicted
print -depsc ../diagrams/dogGeneration_firstOfTest.eps; system('epstopdf ../diagrams/dogGeneration_firstOfTest.eps');
fr=reshape(Varmu2(13,1:nPixels),height,width); imagesc(fr); colormap('gray'); % A subsequent frame
print -depsc ../diagrams/dogGeneration_frame14.eps; system('epstopdf ../diagrams/dogGeneration_frame14.eps');

for i=1:size(Varmu2,1)
    img = reshape(Varmu2(i,1:nPixels), height, width);
    % imagesc(img); colormap('gray');
    imwrite(uint8(img), sprintf('../diagrams/%s_%02d.png',dataSetName, i))
    % pause(0.1);
end

%% Compare the dogs
diffs = zeros(40,2);
for f=1:40,
    imgFlat = imread(sprintf('../diagrams/dog_%02d.png',f));
    imgPyramid = imread(sprintf('../diagrams/dogPyramid_%02d.png',f));
    diffImg = abs(imgFlat - imgPyramid);
    imagesc(diffImg);
    colorbar
    diffs(f,1) = mean(diffImg(:));
    diffs(f,2) = max(diffImg(:));
    fprintf('f=%d, mean_difference=%f, max_difference=%f\n', f, diffs(f,1), diffs(f,2));
    pause(0.1);
end

author: SheffieldML

repo: vargplvm


text:
function [result, why] = structeq(struct1, struct2, funh2string, ignorenan)
% STRUCTEQ performs an equality comparison between two structures by
% recursively comparing the elements of the struct array, their fields and
% subfields. This function requires companion function CELLEQ to compare
% two cell arrays.
%
% USAGE:
%
% structeq(struct1, struct2)
%       Performs a comparison and returns true if all the subfields and
%       properties of the structures are identical. It will fail if
%       subfields include function handles or other objects which don't
%       have a defined eq method.
%
% [iseq, info] = structeq(struct1, struct2)
%       This syntax returns a logical iseq and a second output info which
%       is a structure that contains a field "Reason" which gives you a
%       text stack of why the difference occurred as well as a field
%       "Where" which contains the indices and subfields of the structure
%       where the comparison failed. If iseq is true, info contains empty
%       strings in its fields.
%
% [...] = structeq(struct1, struct2, funh2string, ignorenan)
%       Illustrates an alternate syntax for the function with additional
%       input arguments. See the help for CELLEQ for more information on the
%       meaning of the arguments
%
% METHOD:
% 1. Compare sizes of struct arrays
% 2. Compare numbers of fields
% 3. Compare field names of the arrays
% 4. For every element of the struct arrays, convert the field values into
% a cell array and do a cell array comparison recursively (this can result
% in multiple recursive calls to CELLEQ and STRUCTEQ)
%
% EXAMPLE:
% % Compare two handle graphics hierarchies
% figure;
% g = surf(peaks(50));
% rotate3d
% hg1 = handle2struct(gcf);
% set(g,'XDataMode', 'manual');
% hg2 = handle2struct(gcf);
%
% structeq(hg1, hg2)
% [iseq, info] = structeq(hg1, hg2)
% [iseq, info] = structeq(hg1, hg2, true)

if nargin < 3
    funh2string = false;
end
if nargin < 4
    ignorenan = false;
end

why = struct('Reason','','Where','');

if any(size(struct1) ~= size(struct2))
    result = false;
    why = struct('Reason','Sizes are different',Where,'');
    return
end

fields1 = fieldnames(struct1);
fields2 = fieldnames(struct2);

% Check field lengths
if length(fields1) ~= length(fields2)
    result = false;
    why = struct('Reason','Number of fields are different','Where','');
    return
end

% Check field names
result = celleq(fields1,fields2);
result = all(result);
if ~result
    why = struct('Reason','Field names are different','Where','');
    return
end

for i = 1:numel(struct1)
    props1 = struct2cell(struct1(i));
    props2 = struct2cell(struct2(i));
    [result, subwhy] = celleq(props1,props2,funh2string,ignorenan);
    result = all(result);
    if ~result

        [fieldidx, subwhy.Where] = strtok(subwhy.Where, '}');
        fieldidx = str2double(fieldidx(2:end));
        %str2double(regexp(subwhy.Where,'{([0-9]+)}','tokens','once'));
        where = sprintf('(%d).%s%s',i,fields1{fieldidx},subwhy.Where(2:end));
        why = struct('Reason',sprintf('Properties are different <- %s',subwhy.Reason),'Where',where);
        return
    end
end

author: opencobra

repo: cobratoolbox


text:
% Program to form Admittance And Impedance Bus Formation....

function ybus = ybus();  % Returns ybus

linedata = linedata6(); % Calling "linedata6.m" for Line Data...
fb = linedata(:,1);     % From bus number...
tb = linedata(:,2);     % To bus number...
r = linedata(:,3);      % Resistance, R...
x = linedata(:,4);      % Reactance, X...
b = linedata(:,5);      % Ground Admittance, B/2...
z = r + i*x;            % Z matrix...
y = 1./z;               % To get inverse of each element...
b = i*b;                % Make B imaginary...

nbus = max(max(fb),max(tb));    % no. of buses...
nbranch = length(fb);           % no. of branches...
ybus = zeros(nbus,nbus);        % Initialise YBus...

 % Formation of the Off Diagonal Elements...
 for k=1:nbranch
     ybus(fb(k),tb(k)) = -y(k);
     ybus(tb(k),fb(k)) = ybus(fb(k),tb(k));
 end

 % Formation of Diagonal Elements....
 for m=1:nbus
     for n=1:nbranch
         if fb(n) == m | tb(n) == m
             ybus(m,m) = ybus(m,m) + y(n) + b(n);
         end
     end
 end
 ybus;                  % Bus Admittance Matrix
 zbus = inv(ybus);      % Bus Impedance Matrix

author: Sable

repo: mcbench-benchmarks


text:
function writeFiberSelection (sel, filename)

% function writeFiberSelection (sel, filename)
% Author: DA

fid = fopen (filename, 'wb');

fwrite (fid, length(sel.selected), 'uint');
fwrite (fid, sel.selected, 'uint');

fclose(fid);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [uv, out] = surface_param(x, param)
% surface_param  2D parametrization of a scattered set of 3D points or of
% the vertices of a triangular mesh, for open and closed surfaces.
%
%   surface_param takes a scattered data set X in 3D that is supposed to
%   represent a sampling of a surface, and finds a 2D parameterisation U, V
%   for them. The surface can be open or closed, in which case the
%   parametrization domain will be the R^2 plane or the unit sphere,
%   respectively.
%
%   surface_param provides several parametrization methods, some
%   implemented by us, some by other authors.
%
% [UV, OUT] = surface_param(X, PARAM)
%
%   X is a 3-column matrix. Each row has the coordinates of a point that
%   belongs to the surface we want to interpolate.
%
%   UV is a 2-column matrix with the parameterisation of X, (U, V)->X. In
%   planar parameterisations, UV has units of meters. In spherical
%   parameterisations, UV=[LAT, LON], in units of radians.
%
%   OUT is a struct with the extra output from some of the parametrization
%   methods.
%
%   PARAM is a string or struct that selects the parametrization method
%   and possibly arguments for the method.
%
%     Parametrization methods summary (PARAM.type):
%
%     XY plane:
%                 'xy':       Simple projection on XY.
%                 'pca':      Projection on main PCA plane.
%                 'isomap':   Isomap (with Dijkstra's shortest path).
%                 'cmdsmap':  MDS mapping, global neighbourhood, classical
%                             MDS.
%                 'lmdscale': MDS mapping, local neighbourhood, local MDS.
%
%     Unit sphere:
%                 'sphproj':  Simple projection on sphere around centroid.
%                 'cald':     CALD method by Shen & Makedon 2006.
%                 'smdscale': Spherical extension of lmdscale.
%
%=========================== OPEN SURFACES ================================
%
%   * 'xy': Projection on the xy-plane, i.e. xy-coordinates of the points
%   in X.
%
%--------------------------------------------------------------------------
%
%   * 'pca': Projection on the plane defined by the two eigenvectors with
%     the largest eigenvalues.
%
%--------------------------------------------------------------------------
%
%   * 'isomap': Isometric mapping (Tenenbaum et al. [1][2]). This is a
%     non-linear method that can be applied to complex open surfaces that
%     bend on themselves or are rolled-up. This option uses the
%     third-party function IsomapII. Because isomap uses Dijkstra's
%     shortest-path algorithm internally, it suffers strongly from
%     metrication errors in regular meshes. On the other hand, it can be
%     applied to any adjacency graph, not only meshes.
%
%     Small neighbourhoods capture higher curvature in the surface, but
%     too small neighbourhoods create several connected components and
%     will produce an error.
%
%     PARAM.d: [Opt] Distance matrix between points. The parameter can have
%       the following formats:
%       (1) a full N x N matrix with distances between all pairs of points.
%       In this case, you can define the local neighbourhoods using
%       PARAM.neigh and PARAM.size, so that some of the connections in the
%       full matrix will be deleted and replaced by Disjkstra's
%       shortest-path lengths.
%       (2) a sparse N x N matrix (missing entries are treated as Inf). In
%       this case, you usually want to make PARAM.size=Inf, because the
%       local neighbourhood is already given by the sparse matrix. The
%       missing connections will be replaced internally by Disjkstra's
%       shortest-path lengths.
%       (3) the name of a function (e.g. 'd_fun') that takes one argument,
%       i, and returns a row vector containing the distances from all N
%       points to point i.
%
%       If not provided, PARAM.d is the full Euclidean matrix of distances
%       between every pair of points.
%
%       As a particular case, note that in case of a triangular mesh, the
%       sparse local distance that corresponds to the mesh can be computed
%       using dmatrix_mesh().
%
%     PARAM.neigh: [Opt] Type of neighbourhood ('epsilon' or 'k').
%     PARAM.size:  [Opt] Size of neighbourhood.
%
%       Each point in X is directly connected only to neighbours within a
%       distance PARAM.size (for 'epsilon'), or the PARAM.size closest
%       neighbours (for 'k'), according to the distances in PARAM.d. This
%       allows to further restrict the neighbourhoods in PARAM.d.
%
%     PARAM.options: [Opt] Extra arguments passed to IsomapII. Most users
%       won't need to change this. For details, see help to IsomapII().
%
%--------------------------------------------------------------------------
%
%   * 'cmdsmap': Open surface classical MDS mapping improved with a Fast
%     Marching method (Zigelman et al. [3]). The fast marching algorithm
%     doesn't suffer from metrication errors on regular meshes as badly as
%     Dijkstra's. On the other hand, the current implementation is limited
%     to triangular meshes, while Dijkstra's will work with any adjacency
%     graph. This option uses classical MDS on a full distance matrix, so
%     it's not valid for local neighbourhoods or closed surfaces. For local
%     neighbourhoods, see 'lmdscale' below.
%
%     PARAM.d or PARAM.tri must be provided. If PARAM.d is provided,
%     PARAM.tri and PARAM.options are ignored. Otherwise, PARAM.d is
%     computed from PARAM.tri and PARAM.options.
%
%     PARAM.d: [Req/Opt] This is expected to be the full distance matrix
%       between mesh vertices that can be computed with
%         [~, param.d] = dmatrix_mesh(param.tri, x, 'fastmarching');
%       Note that for comparison purposes, Dijkstra's algorithm can be used
%       instead of Fast Marching
%         [~, param.d] = dmatrix_mesh(param.tri, x, 'dijkstra');
%       The latter is equivalent to the 'isomap' method above.
%
%     If PARAM.d is not provided:
%
%     PARAM.tri: [Req/Opt] 3-column matrix. Each row contains the 3 nodes
%       that form one triangular facet in the mesh.
%
%     PARAM.dmethod: [Opt] String with the method to compute distances
%       between non-neighbours: 'fastmarching' (default) or 'dijkstra'. The
%       latter is faster, but suffers much more from metrication errors.
%
%     PARAM.options: [Opt] Extra arguments passed to the fast marching
%       method. From the help to perform_fast_marching_mesh():
%       - options.W: non-uniform speed.
%       - options.heuristic: heuristic that tries to guess the distance
%         that remains from a given node to a given target. This is an
%         array of same size as W.
%       - parameters that create a local neighbourhood are ignored with a
%         warning.
%
%--------------------------------------------------------------------------
%
%   * 'lmdscale': Open surface local neighbourhood MDS mapping (Schwartz,
%     Wolfson and Shaw [4, 5]), improved with a Fast Marching method
%     (Zigelman et al. [3]). This is like 'cmdsmap' above, but with local
%     neighbourhoods. Local neighbourhoods were already proposed in the
%     first publications [4, 5] on what we are calling MDS mapping. We use
%     our own implementation of this idea, lmdscale().
%
%     PARAM.d or PARAM.tri must be provided. If PARAM.d is provided,
%     PARAM.tri and PARAM.options are ignored. Otherwise, PARAM.d is
%     computed from PARAM.tri and PARAM.options.
%
%     PARAM.d: [Req/Opt] Can be a sparse or full distance matrix. Distances
%       greater than 0 reflect that two vertices are within a common local
%       neighbourhood.
%         [~, param.d] = dmatrix_mesh(param.tri, x, 'fastmarching');
%
%     If PARAM.d is not provided:
%
%     PARAM.tri: [Req/Opt] 3-column matrix. Each row contains the 3 nodes
%       that form one triangular facet in the mesh.
%
%     PARAM.dmethod: [Opt] String with the method to compute distances
%       between non-neighbours: 'fastmarching' (default) or 'dijkstra'. The
%       latter is faster, but suffers much more from metrication errors,
%       and produces only full matrices in the current implementation, so
%       it's not valid for local neighbourhoods.
%
%     PARAM.options: [Opt] Extra arguments passed to the fast marching
%       method. From the help to perform_fast_marching_mesh():
%       - options.constraint_map: Exploration radius to reduce the set of
%         explored points. Only points with current distance smaller than L
%         will be expanded. Set some entries of L to -Inf to avoid any
%         exploration of these points.
%       - options.W: non-uniform speed.
%       - options.end_points: stop when these points are reached.
%       - options.nb_iter_max: stop when a given number of iterations is
%         reached.
%       - options.heuristic: heuristic that tries to guess the distance
%         that remains from a given node to a given target. This is an
%         array of same size as W.
%
%     PARAM.options2: [Opt] Extra arguments passed to the local
%       neighbourhood MDS algorithm. See help to lmdscale().
%       - opt.MaxIter: (default = 20) maximum number of iterations we allow
%         the optimisation algorithm. Each iteration consists of an entire
%         sweep of all points on the sphere.
%       - opt.MaxInc:  inc is the Euclidean distance that each output point
%         is moved in an iteration. The algorithm will stop if no point has
%         been moved more than MaxInc.
%
%     OUT.stopCondition: cell-array with the condition/s that made the
%       algorithm stop, in string form.
%
%     OUT.err:  struct with several error measures at each algorithm
%     iteration. See help to lmdscale() for details.
%
%
%========================== CLOSED SURFACES ===============================
%
%   * 'sphproj':   Direct projection on unit sphere centered on point set
%     centroid.
%
%     OUT.medrad:  Median value of the radii of all points in the point
%       set.
%
%--------------------------------------------------------------------------
%
%   * 'smdscale': Local neighbourhood Multidimensional scaling on a sphere.
%     Our extension to lmdscale() based on method by Agarwal et al. 2010
%     [6].
%
%     PARAM.d or PARAM.tri must be provided. If PARAM.d is provided,
%     PARAM.tri and PARAM.options are ignored. Otherwise, PARAM.d is
%     computed from PARAM.tri and PARAM.options.
%
%     PARAM.d: [Req/Opt] Can be a sparse or full distance matrix. Distances
%       greater than 0 reflect that two vertices are within a common local
%       neighbourhood.
%         [~, param.d] = dmatrix_mesh(param.tri, x, 'fastmarching');
%
%     If PARAM.d is not provided:
%
%     PARAM.tri: [Req/Opt] 3-column matrix. Each row contains the 3 nodes
%       that form one triangular facet in the mesh.
%
%     PARAM.dmethod: [Opt] String with the method to compute distances
%       between non-neighbours: 'fastmarching' (default) or 'dijkstra'. The
%       latter is faster, but suffers much more from metrication errors,
%       and produces only full matrices in the current implementation, so
%       it's not valid for local neighbourhoods.
%
%     PARAM.options: [Opt] Extra arguments passed to the fast marching
%       method. From the help to perform_fast_marching_mesh():
%       - options.constraint_map: Exploration radius to reduce the set of
%         explored points. Only points with current distance smaller than L
%         will be expanded. Set some entries of L to -Inf to avoid any
%         exploration of these points.
%       - options.W: non-uniform speed.
%       - options.end_points: stop when these points are reached.
%       - options.nb_iter_max: stop when a given number of iterations is
%         reached.
%       - options.heuristic: heuristic that tries to guess the distance
%         that remains from a given node to a given target. This is an
%         array of same size as W.
%
%     PARAM.init: [Opt] Parameterization initialisation:
%       - 'sphproj' (default): Direct projection on unit sphere centered on
%         point set centroid, as the option above.
%       - 'random': Random distribution of points on the sphere. This is
%         useful if 'sphproj' produces a fold-over of points that
%         corresponds to a local minimum.
%       - [lat lon]: A 2-column matrix with the initial latitude and
%         longitude coordinates of the parameterization.
%
%     PARAM.sphrad: [Opt] Parameterization sphere's radius. If not
%       provided, it is computed as the median of the distance of the
%       points in X to their centroid.
%
%     PARAM.options2: [Opt] Extra arguments passed to the spherical local
%       neighbourhood MDS algorithm. See help to smdscale().
%       - opt.MaxIter: (default = 20) maximum number of iterations we allow
%         the optimisation algorithm. Each iteration consists of an entire
%         sweep of all points on the sphere.
%       - opt.MaxAlpha:  inc is the angular distance in radians that each
%         output point is moved in an iteration. The algorithm will stop if
%         no point has been moved more than MaxAlpha.
%
%     OUT.stopCondition: cell-array with the condition/s that made the
%       algorithm stop, in string form.
%
%     OUT.err: Error measures vs. iterations:
%       - err.rawstress: Raw stress (sigma_r), as defined above.
%       - err.stress1:   Stress-1 (sigma_1).
%         sigma_1 = sqrt(sigma_r / sum_{i,j} (d_ij)^2)
%       - err.maxalpha:  Maximum angular displacement of any point in each
%         iteration.
%       - err.medalpha:  Median angular displacement of all points in each
%         iteration.
%
%--------------------------------------------------------------------------
%
%   * 'cald': Li Shen's Control Area and Length Distortions (CALD)
%     spherical parametrization [7, 8].
%
%     PARAM.tri: [Req] 3-column matrix. Each row contains the 3 nodes that
%       form one triangular facet in the mesh.
%
%     PARAM.options: [Opt] Extra arguments passed to CALD:
%       - options.MeshGridSize: The interpolation mesh used to smooth the
%         spherical parametrization has length 2*PARAM.MeshGridSize+1. By
%         default, MeshGridSize = 50.
%       - options.MaxSPHARMDegree: Degree of the spherical harmonics used
%         for interpolation to smooth the spherical parametrization. By
%         default, MaxSPHARMDegree = 6.
%       - options.Tolerance: Scalar. In the interpolation smoothing, grid
%         points with a height larger than PARAM.Tolerance*gmin, where gmin
%         is the lowest point in the grid, will be truncated. By default,
%         Tolerance = 2.
%       - options.Smoothing: Scalar. Before parametrization smoothing,
%         parametrization triangle areas relative to triangular mesh areas
%         are smoothed by elevating to (1/PARAM.Smoothing). By default,
%         Smoothing = 2.
%       - options.Iteration: Maximum number of smoothing iterations.
%         Smoothing will stop before reaching Iteration if the algorithm
%         converges to a solution. By default, Iteration = 100.
%       - options.LocalIteration: Number of smoothing iteration in the
%         local smoothing algorithm. By default, LocalIteration = 10.
%
%--------------------------------------------------------------------------
%
% [1] J.B. Tenenbaum, V. de Silva and J.C. Langford, "A Global Geometric
% Framework for Nonlinear Dimensionality Reduction", Science 290(5500):
% 2319-2323, 2000.
%
% [2] Isomap Homepage, http://isomap.stanford.edu/
%
% [3] G. Zigelman, R. Kimmel, and N. Kiryati, "Texture mapping using
% surface flattening via multidimensional scaling,” Visualization and
% Computer Graphics, IEEE Transactions on, vol. 8, no. 2, pp. 198–207,
% 2002.
%
% [4] E. Wolfson and E.L. Schwartz, "Computing minimal distances on
% polyhedral surfaces", IEEE Pattern Analysis and Machine Intelligence,
% 11(9):1001-1005, 1989.
%
% [5] E.L. Schwartz, A. Shaw, E. Wolfson, "A numerical solution to the
% generalized mapmaker’s problem: flattening nonconvex polyhedral", IEEE
% Pattern Analysis and Machine Intelligence, 11(9):1005-1008, 1989.
%
% [6] A. Agarwal et al., "Universal Multi-Dimensional Scaling", Procs. of
% the 16th ACM SIGKDD international conference on Knowledge discovery and
% data mining (KDD '10), p. 1149-1158.
%
% [7] Shen and Makedon, "Spherical mapping for processing of 3D closed
% surfaces", Image and Vision Computing, 24(7):743–761, 2006.
% http://www.sciencedirect.com/science/article/pii/S0262885606000485
%
% [8] CALD and SPHARM-MAT homepage,
% http://www.iupui.edu/~shenlab/software.html
%
% See also: surface_tridomain, surface_interp, scimat_surface2seg.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013-2014 University of Oxford
% Version: 0.6.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 2);

if (size(x, 2) ~= 3)
    error('X must be a 3-column matrix')
end

if ischar(param)
    aux = param; param = [];
    param.type = aux;
end

out = [];

% defaults
if (isempty(param) || ~isfield(param, 'type') ...
        || isempty(param.type))
    error('No PARAM.type provided');
end

% obtain a 2D parameterisation for the input 3D points
switch param.type

%--------------------------------------------------------------------------
% Open surface
%--------------------------------------------------------------------------

    case 'xy'

        % check arguments
        nargoutchk(0, 1);

        % (u,v) is simply (x,y)
        uv = x(:, 1:2);

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

    case 'pca'

        % this function doesn't provide the "out" output
        nargoutchk(0, 1);

        % rotate valve points to make the valve surface as horizontal as
        % possible
        m = mean(x, 1);
        uv = x - m(ones(1, size(x, 1)), :);
        eigv = pts_pca(uv');
        uv = uv * eigv;
        uv = uv(:, 1:2);

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

    case 'isomap'

        % check arguments
        nargoutchk(0, 1);

        % IsomapII takes the following distance matrices:
        %  (1) a full N x N matrix (as in isomap.m)
        %  (2) a sparse N x N matrix (missing entries are treated as INF)
        %  (3) the name of a function (e.g. 'd_fun') that takes
        %        one argument, i, and returns a row vector containing the
        %        distances from all N points to point i.

        % default parameters
        if (~isfield(param, 'neigh') || isempty(param.neigh))
            param.neigh = 'epsilon';
        end
        if (~isfield(param, 'size') || isempty(param.size))
            param.size = Inf;
        end

        % if distance matrix is not provided by the user, it is computed
        % internally as the full Euclidean distance matrix
        if (~isfield(param, 'd') || isempty(param.d))
            param.d = dmatrix(x', x', 'euclidean');
        end

        % compute 2-d projection of the 3-d data
        if (~isfield(param, 'options') || isempty(param.options))
            param.options.dims = 2;
            param.options.display = 0;
            param.options.overlay = 0;
            param.options.verbose = 0;
            param.options.dijkstra = 1;
        end
        uv = IsomapII(param.d, param.neigh, param.size, param.options);
        uv = uv.coords{1}';
        if (size(uv, 1) ~= size(x, 1))
            error('The neighbourhood size (param.size) is too small to connect all points')
        end

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

    case 'cmdsmap'

        % check arguments
        nargoutchk(0, 1);

        % if distance matrix is not provided we need to compute it from the
        % mesh description
        if (~isfield(param, 'd') || isempty(param.d))

            % but if the user hasn't provided the mesh triangulation
            % either, we cannot, so we give an error
            if (~isfield(param, 'tri') || isempty(param.tri))
                error('Either PARAM.d or PARAM.tri must be provided')
            end

            % default for fast marching method options
            if (~isfield(param, 'dmethod') || isempty(param.dmethod))
                param.dmethod = 'fastmarching';
            end
            if (~isfield(param, 'options') || isempty(param.options))
                param.options = [];
            end

            % remove fast marching options that don't apply to this method
            if (isfield(param.options, 'constraint_map'))
                warning('Gerardus:InvalidInputArg', 'PARAM.options.constraint_map ignored, this method only works with a full distance matrix')
                param.options = rmfield(param.options, 'constraint_map');
            end
            if (isfield(param.options, 'end_points'))
                warning('Gerardus:InvalidInputArg', 'PARAM.options.end_points ignored, this method only works with a full distance matrix')
                param.options = rmfield(param.options, 'end_points');
            end
            if (isfield(param.options, 'nb_iter_max'))
                warning('Gerardus:InvalidInputArg', 'PARAM.options.nb_iter_max ignored, this method only works with a full distance matrix')
                param.options = rmfield(param.options, 'nb_iter_max');
            end

            % compute full distance matrix from the mesh
            [~, param.d] ...
                = dmatrix_mesh(param.tri, x, param.dmethod, param.options);

            if (issparse(param.d))
                error('Assertion fail: param.d should be a ');
            end

        end

        % if the user provided a sparse distance matrix, we assume that
        % it's a mistake
        if (issparse(param.d))
            error('PARAM.d must be a full distance matrix, not a sparse mesh adjacency-distance matrix');
        end

        % use classic MDS for the optimization
        uv = cmdscale(param.d, 2);

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

    case 'lmdscale'

        % check arguments
        nargoutchk(0, 2);

        % if distance matrix is not provided we need to compute it from the
        % mesh description
        if (~isfield(param, 'd') || isempty(param.d))

            % but if the user hasn't provided the mesh triangulation
            % either, we cannot, so we give an error
            if (~isfield(param, 'tri') || isempty(param.tri))
                error('Either PARAM.d or PARAM.tri must be provided')
            end

            % default for distance method
            if (~isfield(param, 'options') || isempty(param.options))
                param.options = [];
            end

            % default for fast marching method options
            if (~isfield(param, 'dmethod') || isempty(param.dmethod))
                param.dmethod = 'fastmarching';
            end
            if (~isfield(param, 'options') || isempty(param.options))
                param.options = [];
            end

            % compute neighbourhood distance matrix from the mesh
            [~, param.d] ...
                = dmatrix_mesh(param.tri, x, param.dmethod, param.options);

        end

        % at this point, the distance matrix can be full or sparse, both
        % are valid, unlike in the cases above, where it had to be full

        % compute parameterization using local neighbourhood method
        [u, v, out.stopCondition, out.err] ...
            = lmdscale(param.d, [], [], param.options2);
        uv = [u, v];

%--------------------------------------------------------------------------
% Closed surface
%--------------------------------------------------------------------------

    case 'sphproj'

        % check arguments
        nargoutchk(0, 2);

        % init output
        uv = zeros(size(x, 1), 1);

        % project 3D Cartesian coordinates onto the surface of a sphere
        % lon, lat given in radians, centered around 0
        [uv(:, 2), uv(:, 1), r] = cart2sph(...
            x(:, 1) - mean(x(:, 1)), ...
            x(:, 2) - mean(x(:, 2)), ...
            x(:, 3) - mean(x(:, 3)));

        % we use the median radius of the points in the configuration as the radius
        % of the sphere that best can contains their projections
        out.medrad = median(r);

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

    case 'smdscale'

        % check arguments
        nargoutchk(0, 2);

        % number of points to embed
        N = size(x, 1);

        % if distance matrix is not provided we need to compute it from the
        % mesh description
        if (~isfield(param, 'd') || isempty(param.d))

            % but if the user hasn't provided the mesh triangulation
            % either, we cannot, so we give an error
            if (~isfield(param, 'tri') || isempty(param.tri))
                error('Either PARAM.d or PARAM.tri must be provided')
            end

            % if no options are provided for fastmarching, we need to have
            % at least an empty field
            if (~isfield(param, 'options') || isempty(param.options))
                param.options = [];
            end

            % default for fast marching method options
            if (~isfield(param, 'dmethod') || isempty(param.dmethod))
                param.dmethod = 'fastmarching';
            end
            if (~isfield(param, 'options') || isempty(param.options))
                param.options = [];
            end

            % compute neighbourhood distance matrix from the mesh
            [~, param.d] ...
                = dmatrix_mesh(param.tri, x, param.dmethod, param.options);

        end

        % defaults
        if (~isfield(param, 'init') || isempty(param.init))
            param.init = 'sphproj';
        end

        % sMDS initialisation
        if ischar(param.init) % initial parametrization as a string

            switch param.init
                case 'sphproj'

                    % initial guess for the sphere embedding by simply
                    % projecting each point along the radius to the
                    % centroid
                    [lat, lon, sphrad] = proj_on_sphere(x);

                    % if the user didn't provide an estimate for the
                    % parameterization sphere's radius, we can use sphrad
                    if (~isfield(param, 'sphrad') || isempty(param.sphrad))
                        param.sphrad = sphrad;
                    end

                case 'random'

                    % counter-intuitively, starting from a completely
                    % random distribution usually gives good results
                    lat = (rand(N, 1) - .5) * pi;
                    lon = 2 * (rand(N, 1) - .5) * pi;

                otherwise
                    error('Not valid initialisation for spherical Isomap parameterisation')
            end

        else % initial parametrization as a point set

            if (size(param.init, 1) ~= N || size(param.init, 1) ~= 2)
                error('Initial parametrization provided, but it is not a 2-column matrix with the same number of points as X')
            end

            lat = param.init(:, 1);
            lon = param.init(:, 2);

        end

        % if the user has not provided an estimate for the parameterization
        % sphere's radius, we have to compute it internally
        if (~isfield(param, 'sphrad') || isempty(param.sphrad))
            [~, ~, param.sphrad] = proj_on_sphere(x);
        end

        % if no options are provided for the MDS algorithm, we need to have
        % at least an empty field
        if (~isfield(param, 'options2') || isempty(param.options2))
            param.options = [];
        end

        % compute the parametrization that optimises isometry
        [lat, lon, out.stopCondition, out.err] = ...
            smdscale(param.d, param.sphrad, lat, lon, param.options2);

        % create output with the parameterisation values for each point
        % lat = em(1, :)
        % lon = em(2, :)
        uv = [lat lon];

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

    case 'cald'

        % check arguments
        nargoutchk(0, 1);

        % required parameters
        if (~isfield(param, 'tri') || isempty(param.tri))
            error('PARAM.tri must be provided')
        end

        % set up parameters for the CALD algorithm. We have to create a
        % structure with the right format. We can then replace default
        % values by values provided by the user
        param.options.vars = {'MeshGridSize', 'MaxSPHARMDegree', 'Tolerance', ...
            'Smoothing', 'Iteration', 'LocalIteration', 't_major', ...
            'SelectDiagonal', 'OutDirectory'};
        param.options.args = [1 1 1 1 1 1 10 10 200];
        param.options.inFilter = {
            '*_bim.mat;*_fix.mat;*_obj.mat' 'Binary and Mesh Objects (*_bim.mat;*_fix.mat;*_obj.mat)'
            '*_obj.mat'                     'Mesh Objects (*_obj.mat)'
            '*_bim.mat;*_fix.mat'           'Binary Objects (*_bim.mat, *_fix.mat)'
            '*.mat'                         'All Objects (*.mat)'
            };
        param.options.default = {'50'  '6'  '2'  '2'  '100'  '10'  ''  ''  './'};
        if (~isfield(param.options, 'MeshGridSize') || isempty(param.options.MeshGridSize))
            param.options.MeshGridSize = 50;
        end
        if (~isfield(param.options, 'MaxSPHARMDegree') || isempty(param.options.MaxSPHARMDegree))
            param.options.MaxSPHARMDegree = 6;
        end
        if (~isfield(param.options, 'Tolerance') || isempty(param.options.Tolerance))
            param.options.Tolerance = 2;
        end
        if (~isfield(param.options, 'Smoothing') || isempty(param.options.Smoothing))
            param.options.Smoothing = 2;
        end
        if (~isfield(param.options, 'Iteration') || isempty(param.options.Iteration))
            param.options.Iteration = 100;
        end
        if (~isfield(param.options, 'LocalIteration') || isempty(param.options.LocalIteration))
            param.options.LocalIteration = 10;
        end
        param.options.t_major = 'x';
        param.options.SelectDiagonal = 'ShortDiag';
        param.options.OutDirectory = '.';
        param.options = orderfields(param.options, ...
            {'vars', 'args', 'inFilter', 'default', ...
            'MeshGridSize', 'MaxSPHARMDegree', 'Tolerance', ...
            'Smoothing', 'Iteration', 'LocalIteration', 't_major', ...
            'SelectDiagonal', 'OutDirectory'});

        % initial parametrization
        [sph_verts, name3] ...
            = initParamCALD(x, param.tri, '', param.options);

        if (any(isnan(sph_verts)))
            error('Initial parameterization produced NaN values')
        end

        % smooth the parameterization
        [~, ~, sph_verts, new_name] ...
            = smootheCALD(x, param.tri, sph_verts, name3, param.options);

        % convert xyz coordinates to spherical coordinates
        [lon, lat] ...
            = cart2sph(sph_verts(:, 1), sph_verts(:, 2), sph_verts(:, 3));

        % delete files created by the SPHARM functions
        if exist(new_name, 'file')
            delete(new_name) % CALD_smo.mat
        end
        if exist('initParamCALD', 'file')
            rmdir('initParamCALD', 's')
        end

        % create output
        uv = [lat lon];

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

    otherwise

        error('Parametrization method not implemented')
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
clear
close all
clc;

%time
dt=0.01;

%coeffs
Cg=100.;
Cd=.1;
minSpeed=.5;

name='greenspikes.png.tiff'
sanitizedname=[strrep(name,'.','') '_']
mkdir([sanitizedname 'outputdir'])
%green
%greenheight=double(imread('testgreen.tiff'))/256*1.7;
%greenheight=double(imread('greencone.tiff'))/256*.3;
%greenheight=double(imread('greenconeinv.tiff'))/256*.3;
greenheight=fliplr(rot90(rot90(double(imread(name))/256*.8)));
greenheightPLOT=fliplr(rot90(rot90(double(imread(name))/256*.8)));
[gradX,gradY]=gradient(greenheight(:,:,1));
forceFun=@(R,RL) (RL-R)/dt*Cd + [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))]*Cg;
%forceFun=@(R,RL) [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))];
l=size(greenheight,1);

%pin
holeLoc=size(greenheight)'/2;
holeRadius=3;

%graphics
%sph=sphere(30)
alt=30.85;
az=-39+90;
cmapgreen=flipud([181,228,138; 160,220,104; 124,216,87; 81,212,76; 52,188,67; 37,167,60; 32,154,61; 1,114,56; 0,86,19]/255);
stickX=[holeLoc(1),holeLoc(1)];
stickY=[holeLoc(2),holeLoc(2)];
stickZ=[greenheight(floor(holeLoc(1)),floor(holeLoc(1))),greenheight(floor(holeLoc(1)),floor(holeLoc(1)))+100];
flagX=[stickX;stickX+40];
flagY=[stickY;stickY];
flagZ=[stickZ(2),stickZ(2)-30;stickZ(2),stickZ(2)-30];
flagC=zeros(2,2,3);
flagC(:,:,1)=.9;

% balls
startLoc=[64;64];

%angles=(45)*pi/180;
%speeds=400;
%angles=linspace((45-15)*pi/180,(45+15)*pi/180,10);%-pi/4;
angles=linspace((45-90)*pi/180,(45+90)*pi/180,100);%-pi/4;
speeds=linspace(minSpeed,150,100);

numParticles=length(speeds)*length(angles);

speeds=reshape(speeds,1,1,length(speeds));
angleVectors=[cos(angles);sin(angles)];
startConditions=bsxfun(@times,speeds,angleVectors);
startConditions=reshape(startConditions,2,length(speeds)*length(angles));

r=repmat(startLoc,1,numParticles);
rl=r-startConditions*dt;

rstart=r;
rlstart=rl;

%sequentialhits staging
%each ball laucnhes after 1/4 second (300 iters/4) = 80 ticks
launchTicks=80

i=0;
while sum(sum(r~=rl))
    %energy verlet + drag
    i=i+1;
    rn=2*r-rl+(forceFun(r,rl))*dt^2;
    rl=r;
    r=rn;

    %%%LAUNCH CONDITIONS
    if i<50
        r=rstart;
        rl=rlstart;
    end

    %%%STOP CONDITIONS
    %static friction
    dr=rl-r;
    s=sqrt(dr(1,:).^2+dr(2,:).^2)/dt;
    %walls
    haltBallsEdge=(r>l)|(r<1);
    %in the hole
    distToHole=bsxfun(@minus,holeLoc,rl);%use rl not r to check if the ball is in the hole so it can't escape
    distToHole=sqrt(distToHole(1,:).^2+distToHole(2,:).^2);
    %stop those in stop conditions
    holeBalls=(distToHole<holeRadius);
    rl(:,holeBalls)=repmat(holeLoc,1,sum(holeBalls));
    haltBalls=(haltBallsEdge(1,:)|haltBallsEdge(2,:))|(s<minSpeed)|(distToHole<holeRadius);
    r(:,haltBalls)=rl(:,haltBalls);

    if mod(i,10)==0
        figure(1)
        surf(greenheightPLOT,'edgecolor','none')
        colormap(cmapgreen)
        set(gca,'DataAspectRatio',[1 1 1])
        hold on
        scatter3(r(1,:),r(2,:),3+greenheightPLOT(sub2ind(size(greenheightPLOT), floor(r(2,:)),floor(r(1,:)))),'wo','filled')
        plot3(stickX,stickY,stickZ,'k-')
        surf(flagX,flagY,flagZ,flagC,'edgecolor','none')
        hold off
        az=az+.1;
        view([az alt])
%         imshow(greenheight,[min(min(greenheight)),max(max(greenheight))],'colormap',colormap('parula'))
        hold on
        %rectangle('Position',[holeLoc(1)-holeRadius holeLoc(2)-holeRadius holeRadius*2 holeRadius*2],'Curvature',[1,1],'FaceColor','black');
        %scatter(r(1,:),r(2,:),'w.')
        hold off
        xlim([0,512])
        ylim([0,512])
        %title(num2str(s(1)))
        set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
        saveas(gcf,[[sanitizedname 'outputdir'] '\' num2str(i/10,'%05i') '.png'])
    end
end

%phase space

figure(2)
maxDist=256;
maxColor=reshape([1,0,0],1,1,3);
minColor=reshape([1,1,1],1,1,3);
dists=fliplr(flipud(reshape(distToHole,length(angles),length(speeds))'));
imshow(dists,[holeRadius,maxDist],'colormap',colormap('jet'))
hold on
h=imshow(repmat(ones(size(dists)),1,1,3));
hold off
set(h,'AlphaData',dists<holeRadius);
ylabel('Putt speed')
xlabel('putt angle')

author: BrianHaidet

repo: AlphaPhoenix


text:
classdef LaunchVehicleViewProfile < matlab.mixin.SetGet
    %LaunchVehicleViewProfile Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %profile properties
        name(1,:) char = 'Untitled View Profile';
        desc(1,1) string = "";

        %axes properties
        backgroundColor(1,1) ColorSpecEnum = ColorSpecEnum.White;
        gridType(1,1) ViewGridTypeEnum = ViewGridTypeEnum.Major;
        majorGridColor(1,1) ColorSpecEnum = ColorSpecEnum.DarkGrey;
        minorGridColor(1,1) ColorSpecEnum = ColorSpecEnum.DarkGrey;
        gridTransparency(1,1) double = 0.15;
        meshEdgeAlpha(1,1) double = 0.1;

        %render mode
        renderer(1,1) FigureRendererEnum = FigureRendererEnum.OpenGL;

        %axis properties
        dispXAxis(1,1) logical = false;
        dispYAxis(1,1) logical = false;
        dispZAxis(1,1) logical = false;

        %trajectory view options
        trajEvtsViewType(1,1) ViewEventsTypeEnum = ViewEventsTypeEnum.All %either chunked by event/SoI or all
        frame AbstractReferenceFrame
        projType ViewProjectionTypeEnum = ViewProjectionTypeEnum.Perspective;

        %body fixed options
        showLongLatAnnotations(1,1) logical = true;

        %thrust vectors
        showThrustVectors(1,1) logical = false;
        thrustVectColor(1,1) ColorSpecEnum  = ColorSpecEnum.Red;
        thrustVectLineType(1,1) LineSpecEnum = LineSpecEnum.SolidLine;
        thrustVectScale(1,1) double = 1;
        thrustVectEntryIncr(1,1) double = 1;

        %SoI and other body options
        showSoIRadius(1,1) logical = false;
        bodiesToPlot(1,:) KSPTOT_BodyInfo = KSPTOT_BodyInfo.empty(1,0);
        bodyPlotStyle(1,1) ViewProfileBodyPlottingStyle = ViewProfileBodyPlottingStyle.Dot;

        %show sc body axes
        showScBodyAxes(1,1) logical = false;
        scBodyAxesScale(1,1) double = 100; %km

        %lighting
        showLighting(1,1) logical = false;
        showSunVect(1,1) logical = false;

        %show atmosphere
        showAtmosphere(1,1) logical = false;

        %ground objects
        groundObjsToPlot(1,:) LaunchVehicleGroundObject
        showGndTracks(1,1) logical = true;
        showGrdObjLoS(1,1) logical = true;

        %geometric points
        pointsToPlot(1,:) AbstractGeometricPoint

        %geometric vectors
        vectorsToPlot(1,:) AbstractGeometricVector

        %geometric reference frames
        refFramesToPlot(1,:) AbstractGeometricRefFrame

        %geometric angles
        anglesToPlot(1,:) AbstractGeometricAngle

        %geometric planes
        planesToPlot(1,:) AbstractGeometricPlane

        %sensors
        sensorsToPlot(1,:) AbstractSensor

        %sensor targets
        sensorTgtsToPlot(1,:) AbstractSensorTarget

        %central body transform
        hCBodySurfXForm = []; %matlab.graphics.GraphicsPlaceholder()

        %view properties (set by user indirectly through UI controls)
        updateViewAxesLimits(1,1) logical = false;
        orbitNumToPlot(1,1) double = 1;
        viewAzEl(1,2) = [-37.5, 30]; %view(3)
        viewZoomAxLims(3,2) double = NaN(3,2);
        viewCameraPosition(1,3) double = NaN(1,3);
        viewCameraTarget(1,3) double = NaN(1,3);
        viewCameraUpVector(1,3) double = NaN(1,3);
        viewCameraViewAngle(1,1) double = NaN(1,1);
    end

    properties(Transient)
        markerTrajData(1,:) LaunchVehicleViewProfileTrajectoryData = LaunchVehicleViewProfileTrajectoryData.empty(1,0);
        markerBodyData(1,:) LaunchVehicleViewProfileBodyData = LaunchVehicleViewProfileBodyData.empty(1,0);
        markerTrajAxesData(1,:) LaunchVehicleViewProfileBodyAxesData = LaunchVehicleViewProfileBodyAxesData.empty(1,0);
        markerGrdObjData(1,:) LaunchVehicleViewProfileGroundObjData = LaunchVehicleViewProfileGroundObjData.empty(1,0);
        sunLighting(1,:) LaunchVehicleViewProfileSunLighting = LaunchVehicleViewProfileSunLighting.empty(1,0);
        centralBodyData(1,:) LaunchVehicleViewProfileCentralBodyData = LaunchVehicleViewProfileCentralBodyData.empty(1,0);

        pointData(1,:) LaunchVehicleViewProfilePointData = LaunchVehicleViewProfilePointData.empty(1,0);
        vectorData(1,:) LaunchVehicleViewProfileVectorData = LaunchVehicleViewProfileVectorData.empty(1,0);
        refFrameData(1,:) LaunchVehicleViewProfileRefFrameData = LaunchVehicleViewProfileRefFrameData.empty(1,0);
        angleData(1,:) LaunchVehicleViewProfileAngleData = LaunchVehicleViewProfileAngleData.empty(1,0);
        planeData(1,:) LaunchVehicleViewProfilePlaneData = LaunchVehicleViewProfilePlaneData.empty(1,0);

        sensorData(1,:) LaunchVehicleViewProfileSensorData = LaunchVehicleViewProfileSensorData.empty(1,0);
        sensorTgtData(1,:) LaunchVehicleViewProfileSensorTargetData = LaunchVehicleViewProfileSensorTargetData.empty(1,0);

        userDefinedRefFrames(1,:)
    end

    properties(Access=private)
        generic3DTrajView(1,1) Generic3DTrajectoryViewType = Generic3DTrajectoryViewType();
    end

    methods
        function obj = LaunchVehicleViewProfile()

        end

        function removeGrdObjFromList(obj, grdObj)
            obj.groundObjsToPlot([obj.groundObjsToPlot] == grdObj) = [];
        end

        function removeGeoPointFromList(obj, point)
            obj.pointsToPlot([obj.pointsToPlot] == point) = [];
        end

        function removeGeoVectorFromList(obj, vector)
            obj.vectorsToPlot([obj.vectorsToPlot] == vector) = [];
        end

        function removeGeoRefFrameFromList(obj, refFrame)
            obj.refFramesToPlot([obj.refFramesToPlot] == refFrame) = [];
        end

        function removeGeoAngleFromList(obj, angle)
            obj.anglesToPlot([obj.anglesToPlot] == angle) = [];
        end

        function removeGeoPlaneFromList(obj, plane)
            obj.planesToPlot([obj.planesToPlot] == plane) = [];
        end

        function removeSensorFromList(obj, sensor)
            obj.sensorsToPlot([obj.sensorsToPlot] == sensor) = [];
        end

        function removeSensorTargetFromList(obj, target)
            obj.sensorTgtsToPlot([obj.sensorTgtsToPlot] == target) = [];
        end

        function plotTrajectory(obj, lvdData, handles, app)
            obj.generic3DTrajView.plotStateLog(obj.orbitNumToPlot, lvdData, obj, handles, app);
        end

        function createTrajectoryMarkerData(obj, subStateLogs, evts)
            obj.clearAllTrajData();
            trajMarkerData = obj.createTrajData();

            for(j=1:length(subStateLogs)) %#ok<*NO4LP>
                if(size(subStateLogs{j},1) > 0)
                    [times, rVects, evtColor] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, j, evts);

                    if(length(unique(times)) > 0)
                        trajMarkerData.addData(times, rVects, evtColor);
                    end
                end
            end
        end

        function createBodyMarkerData(obj, dAxes, subStateLogs, viewInFrame, showSoI, meshEdgeAlpha, evts)
            obj.clearAllBodyData();

            timesArr = {};
            rVectArr = {};
            bodyColorsArr = {};
            bodyMarkerDataObjs = {};

            for(i=1:length(obj.bodiesToPlot))
                bodyToPlot = obj.bodiesToPlot(i);

                if(bodyToPlot == viewInFrame.getOriginBody())
                    continue;
                end

                bColorRGB = bodyToPlot.getBodyRGB();

                bodyOrbitPeriod = [];
                if(bodyToPlot.sma > 0)
                    gmu = bodyToPlot.getParentGmuFromCache();

                    if(not(isempty(gmu)) && not(isnan(gmu)) && isfinite(gmu))
                        bodyOrbitPeriod = computePeriod(bodyToPlot.sma, gmu);
                    end
                else
                    bodyOrbitPeriod = Inf;
                end

                bodyMarkerData = obj.createBodyData(bodyToPlot, viewInFrame, obj.bodyPlotStyle, showSoI, meshEdgeAlpha);

                times = [];
                for(j=1:length(subStateLogs)) %#ok<NO4LP>
                    if(size(subStateLogs{j},1) > 0)
                        times = [times; subStateLogs{j}(:,1)]; %#ok<AGROW>
                    end
                end

                if(isfinite(bodyOrbitPeriod))
                    numPeriods = (max(times) - min(times))/bodyOrbitPeriod;
                    times = linspace(min(times), max(times), max(10*numPeriods,length(times)));
                else
                    times = linspace(min(times), max(times), length(times));
                end

                states = bodyToPlot.getElementSetsForTimes(times);
                if(numel(states) >= 1)
                    doConversion = states(1).frame ~= viewInFrame;
                else
                    doConversion = true;
                end

                if(doConversion)
                    states = convertToFrame(states, viewInFrame);
                end

                rVects = [states.rVect];

                [times,ia,~] = unique(times,'stable');
                rVects = rVects(:,ia);

                [times,I] = sort(times);
                rVects = rVects(:,I);

                timesInner = times;
                rVectsInner = rVects;

                timesArr{i} = timesInner; %#ok<AGROW>
                rVectArr{i} = rVectsInner; %#ok<AGROW>

                bodyColorsArr{i} = bColorRGB; %#ok<AGROW>
                bodyMarkerDataObjs{i} = bodyMarkerData; %#ok<AGROW>
            end

            for(i=1:length(timesArr))
                timesInner = timesArr{i};
                rVectsInner = rVectArr{i};
                bColorRGB = bodyColorsArr{i};
                bodyMarkerData = bodyMarkerDataObjs{i};

                if(not(isempty(timesInner)))
                    obj.markerBodyData(end+1) = bodyMarkerData;

                    times = timesInner;
                    rVects = rVectsInner;

                    if(length(unique(times)) > 1)
                        plot3(dAxes, rVects(1,:), rVects(2,:), rVects(3,:), '-', 'Color',bColorRGB, 'LineWidth',1.5);
                        bodyMarkerData.addData(times, rVects);
                    end
                end
            end
        end

        function createSunLightSrc(obj, dAxes, viewInFrame)
            obj.sunLighting = LaunchVehicleViewProfileSunLighting(dAxes, viewInFrame, obj.showLighting, obj.showSunVect);
        end

        function updateLightPosition(obj, time)
            obj.sunLighting.updateSunLightingPosition(time);
        end

        function createBodyAxesData(obj, vehPosVelData, vehAttData) %lvdStateLogEntries, evts, viewInFrame
            obj.clearAllBodyAxesData();
            obj.markerTrajAxesData = LaunchVehicleViewProfileBodyAxesData(vehPosVelData, vehAttData, obj.scBodyAxesScale, obj.showScBodyAxes);

%             if(obj.showScBodyAxes)
%                 for(i=1:length(evts))
%                     evt = evts(i);
%                     evtStateLogEntries = lvdStateLogEntries([lvdStateLogEntries.event] == evt);
%                     evtStateLogEntries = evtStateLogEntries(:)';
%
%                     cartElem = convertToFrame(getCartesianElementSetRepresentation(evtStateLogEntries), viewInFrame);
%
%                     times = [evtStateLogEntries.time];
%
%                     rVects = [cartElem.rVect];
%                     rotMatsBodyToView = NaN(3, 3, length(evtStateLogEntries));
%                     for(j=1:length(evtStateLogEntries))
%                         %get body position in view frame
%                         entry = evtStateLogEntries(j);
%
%                         %get body axes in view frame
%                         rotMatBodyToInertial = entry.steeringModel.getBody2InertialDcmAtTime(entry.time, entry.position, entry.velocity, entry.centralBody);
%
%                         [~, ~, ~, rotMatToInertial12] = viewInFrame.getOffsetsWrtInertialOrigin(entry.time, cartElem(j));
%                         [~, ~, ~, rotMatToInertial32] = entry.centralBody.getBodyCenteredInertialFrame().getOffsetsWrtInertialOrigin(entry.time, cartElem(j));
%
%                         rotMatsBodyToView(:,:,j) = rotMatToInertial12' * rotMatToInertial32 * rotMatBodyToInertial; %body to inertial -> inertial to inertial -> inertial to view frame
%                     end
%
%                     [times,ia,~] = unique(times,'stable');
%                     rVects = rVects(:,ia);
%                     rotMatsBodyToView = rotMatsBodyToView(:,:,ia);
%
%                     [times,I] = sort(times);
%                     rVects = rVects(:,I);
%                     rotMatsBodyToView = rotMatsBodyToView(:,:,I);
%
%                     switch(evt.plotMethod)
%                         case EventPlottingMethodEnum.PlotContinuous
%                             %nothing
%
%                         case EventPlottingMethodEnum.SkipFirstState
%                             times = times(2:end);
%                             rVects = rVects(:,2:end);
%                             rotMatsBodyToView = rotMatsBodyToView(:,:,2:end);
%
%                         case EventPlottingMethodEnum.DoNotPlot
%                             times = [];
%                             rVects = [];
%                             rotMatsBodyToView = [];
%
%                         otherwise
%                             error('Unknown event plotting method: %s', EventPlottingMethodEnum.DoNotPlot.name);
%                     end
%
%                     if(length(times) >= 2 && all(diff(times)>0))
%                         obj.markerTrajAxesData.addData(times, rVects, rotMatsBodyToView);
%                     end
%                 end
%             end
        end

        function createGroundObjMarkerData(obj, dAxes, lvdStateLogEntries, vehPosVelData, evts, viewInFrame, celBodyData)
            obj.clearAllGrdObjData();

            for(i=1:length(obj.groundObjsToPlot))
                grdObj = obj.groundObjsToPlot(i);

                if(ismember(grdObj.centralBodyInfo, obj.bodiesToPlot) || grdObj.centralBodyInfo == viewInFrame.getOriginBody())
                    grdObjData = LaunchVehicleViewProfileGroundObjData(grdObj, vehPosVelData, celBodyData);
                    obj.markerGrdObjData(end+1) = grdObjData;

                    for(j=1:length(evts))
                        evt = evts(j);
                        evtStateLogEntries = lvdStateLogEntries([lvdStateLogEntries.event] == evt);
                        numEntries = length(evtStateLogEntries);

%                         scCartElem = convertToFrame(getCartesianElementSetRepresentation(evtStateLogEntries), viewInFrame);
                        allTimes = [evtStateLogEntries.time];

                        elemSet = GeographicElementSet.empty(1,0);

                        for(k=1:numEntries)
                            ge = grdObj.getStateAtTime(allTimes(k));
                            if(not(isempty(ge)))
                                elemSet(end+1) = ge; %#ok<AGROW>
                            end
                        end
                        elemSet = convertToFrame(convertToCartesianElementSet(elemSet), viewInFrame);

                        times = [elemSet.time];
                        rVectsGrdObj = [elemSet.rVect];
%                         rVectsSc= [scCartElem.rVect];

                        [times,ia,~] = unique(times);
                        rVectsGrdObj = rVectsGrdObj(:,ia);
%                         rVectsSc = rVectsSc(:,ia);

                        switch(evt.plotMethod)
                            case EventPlottingMethodEnum.PlotContinuous
                                %nothing

                            case EventPlottingMethodEnum.SkipFirstState
                                times = times(2:end);
                                rVectsGrdObj = rVectsGrdObj(2:end,:);
%                                 rVectsSc = rVectsSc(2:end,:);

                            case EventPlottingMethodEnum.DoNotPlot
                                times = [];
                                rVectsGrdObj = [];
%                                 rVectsSc = [];

                            otherwise
                                error('Unknown event plotting method: %s', EventPlottingMethodEnum.DoNotPlot.name);
                        end

                        if(length(times) >= 2 && all(diff(times)>0))
                            grdObjData.addData(times, rVectsGrdObj, viewInFrame, obj.showGrdObjLoS);

                            if(obj.showGndTracks)
                                hold(dAxes,'on');
                                plot3(dAxes, rVectsGrdObj(1,:), rVectsGrdObj(2,:), rVectsGrdObj(3,:), 'Color',grdObj.grdTrkLineColor.color, 'LineStyle',grdObj.grdTrkLineSpec.linespec);
                                hold(dAxes,'off');
                            end
                        end
                    end
                end
            end
        end

        function configureTimeSlider(obj, minTime, maxTime, subStateLogs, handles, app)
            timeSlider = app.DispAxesTimeSlider;
            curSliderTime = timeSlider.Value;
            if(not(isfinite(minTime) && isfinite(maxTime)) || ...
               minTime == maxTime)
                onlyTime = subStateLogs{1}(1,1);

                minTime = onlyTime;
                maxTime = onlyTime + 1;
            end
            timeSlider.Limits = [minTime maxTime];
            timeSlider.MajorTicks = linspace(minTime, maxTime, 10);
            timeSlider.MinorTicks = linspace(minTime, maxTime, 100);
            timeSlider.MajorTickLabels = "";

            if(curSliderTime > maxTime)
                timeSlider.Value = maxTime;
            elseif(curSliderTime < minTime)
                timeSlider.Value = minTime;
            end

            lvdData = getappdata(handles.ma_LvdMainGUI,'lvdData');
            timeSliderCb = @(src,evt) timeSliderStateChanged(src,evt, lvdData, handles, app);
            timeSlider.ValueChangingFcn = timeSliderCb;

            timeSlider.ValueChangingFcn(timeSlider, true);
        end

        function trajData = createTrajData(obj)
            trajData = LaunchVehicleViewProfileTrajectoryData();
            obj.markerTrajData = trajData;
        end

        function createCentralBodyData(obj, bodyInfo, hCBodySurfXForm, viewFrame)
            obj.centralBodyData = LaunchVehicleViewProfileCentralBodyData(bodyInfo, hCBodySurfXForm, viewFrame);
        end

        function bodyData = createBodyData(obj, bodyInfo, viewInFrame, bodyPlotStyle, showSoI,meshEdgeAlpha)
            bodyData = LaunchVehicleViewProfileBodyData(bodyInfo, viewInFrame, bodyPlotStyle, showSoI, meshEdgeAlpha);
            obj.markerBodyData(end+1) = bodyData;
       end

        function createPointData(obj, viewFrame, subStateLogs, evts)
            obj.clearPointData();
            points = obj.pointsToPlot;

            for(i=1:length(points))
                obj.pointData(end+1) = LaunchVehicleViewProfilePointData(points(i), viewFrame);
            end

            for(i=1:length(subStateLogs))
                if(size(subStateLogs{i},1) > 0)
                    [times, rVects, ~, vVects] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, i, evts);

                    if(length(unique(times)) > 1)
                        for(j=1:length(obj.pointData))
                            obj.pointData(j).addData(times, rVects, vVects);
                        end
                    end
                end
            end
        end

        function createVectorData(obj, viewFrame, subStateLogs, evts)
            obj.clearVectorData();
            vectors = obj.vectorsToPlot;

            for(i=1:length(vectors))
                obj.vectorData(end+1) = LaunchVehicleViewProfileVectorData(vectors(i), viewFrame);
            end

            for(i=1:length(subStateLogs))
                if(size(subStateLogs{i},1) > 0)
                    [times, rVects, ~, vVects] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, i, evts);

                    if(length(unique(times)) > 1)
                        for(j=1:length(obj.vectorData))
                            obj.vectorData(j).addData(times, rVects, vVects);
                        end
                    end
                end
            end
        end

        function createRefFrameData(obj, viewFrame, subStateLogs, evts)
            obj.clearRefFrameData();
            refFrames = obj.refFramesToPlot;

            for(i=1:length(refFrames))
                obj.refFrameData(end+1) = LaunchVehicleViewProfileRefFrameData(refFrames(i), viewFrame);
            end

            for(i=1:length(subStateLogs))
                if(size(subStateLogs{i},1) > 0)
                    [times, rVects, ~, vVects] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, i, evts);

                    if(length(unique(times)) > 1)
                        for(j=1:length(obj.refFrameData))
                            obj.refFrameData(j).addData(times, rVects, vVects);
                        end
                    end
                end
            end
        end

        function createAngleData(obj, viewFrame, subStateLogs, evts)
            obj.clearAngleData();
            angles = obj.anglesToPlot;

            for(i=1:length(angles))
                obj.angleData(end+1) = LaunchVehicleViewProfileAngleData(angles(i), viewFrame);
            end

            for(i=1:length(subStateLogs))
                if(size(subStateLogs{i},1) > 0)
                    [times, rVects, ~, vVects] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, i, evts);

                    if(length(unique(times)) > 1)
                        for(j=1:length(obj.angleData))
                            obj.angleData(j).addData(times, rVects, vVects);
                        end
                    end
                end
            end
        end

        function createPlaneData(obj, viewFrame, subStateLogs, evts)
            obj.clearPlaneData();
            planes = obj.planesToPlot;

            for(i=1:length(planes))
                obj.planeData(end+1) = LaunchVehicleViewProfilePlaneData(planes(i), viewFrame);
            end

            for(i=1:length(subStateLogs))
                if(size(subStateLogs{i},1) > 0)
                    [times, rVects, ~, vVects] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, i, evts);

                    if(length(unique(times)) > 1)
                        for(j=1:length(obj.planeData))
                            obj.planeData(j).addData(times, rVects, vVects);
                        end
                    end
                end
            end
        end

        function createSensorData(obj, lvdStateLogEntries, vehPosVelData, vehAttData, viewFrame)
            obj.clearSensorData();

            for(i=1:length(obj.sensorsToPlot))
                sensor = obj.sensorsToPlot(i);
                obj.sensorData(i) = LaunchVehicleViewProfileSensorData(sensor, obj.sensorTgtsToPlot, lvdStateLogEntries, vehPosVelData, vehAttData, viewFrame);
            end
        end

        function createSensorTargetData(obj, viewFrame)
            obj.clearSensorTgtData();

            for(i=1:length(obj.sensorTgtsToPlot))
                target = obj.sensorTgtsToPlot(i);
                obj.sensorTgtData(i) = LaunchVehicleViewProfileSensorTargetData(target, viewFrame);
            end
        end

        function clearAllTrajData(obj)
            obj.markerTrajData = LaunchVehicleViewProfileTrajectoryData.empty(1,0);
        end

        function clearAllBodyData(obj)
            obj.markerBodyData = LaunchVehicleViewProfileBodyData.empty(1,0);
        end

        function clearAllBodyAxesData(obj)
            obj.markerTrajAxesData = LaunchVehicleViewProfileBodyAxesData.empty(1,0);
        end

        function clearAllGrdObjData(obj)
            obj.markerGrdObjData = LaunchVehicleViewProfileGroundObjData.empty(1,0);
        end

        function clearPointData(obj)
            obj.pointData = LaunchVehicleViewProfilePointData.empty(1,0);
        end

        function clearVectorData(obj)
            obj.vectorData = LaunchVehicleViewProfileVectorData.empty(1,0);
        end

        function clearRefFrameData(obj)
            obj.refFrameData = LaunchVehicleViewProfileRefFrameData.empty(1,0);
        end

        function clearAngleData(obj)
            obj.angleData = LaunchVehicleViewProfileAngleData.empty(1,0);
        end

        function clearPlaneData(obj)
            obj.planeData = LaunchVehicleViewProfilePlaneData.empty(1,0);
        end

        function clearSensorData(obj)
            obj.sensorData = LaunchVehicleViewProfileSensorData.empty(1,0);
        end

        function clearSensorTgtData(obj)
            obj.sensorTgtData = LaunchVehicleViewProfileSensorTargetData.empty(1,0);
        end

    end

    methods(Static)
        function vehPosVelData = createVehPosVelData(subStateLogs, evts, viewFrame)
            vehPosVelData = LaunchVehicleViewPosVelInterp(viewFrame);

            for(i=1:length(subStateLogs))
                if(size(subStateLogs{i},1) > 0)
                    [times, rVects, ~, vVects] = LaunchVehicleViewProfile.parseTrajDataFromSubStateLogs(subStateLogs, i, evts);

                    if(length(unique(times)) > 1)
                        vehPosVelData.addData(times, rVects, vVects);
                    end
                end
            end
        end

        function vehAttData = createVehAttitudeData(vehPosVelData, lvdStateLogEntries, evts, viewInFrame)
            vehAttData = LaunchVehicleViewProfileAttitudeData();

            for(i=1:length(evts))
                evt = evts(i);
                bool = [lvdStateLogEntries.event] == evt;
                evtStateLogEntries = lvdStateLogEntries(bool);
                evtStateLogEntries = evtStateLogEntries(:)';

                if(isempty(evtStateLogEntries))
                    continue;
                end

                cartElem = convertToFrame(getCartesianElementSetRepresentation(evtStateLogEntries), viewInFrame);

                times = [evtStateLogEntries.time];

                rotMatsBodyToView = NaN(3, 3, length(evtStateLogEntries));
                for(j=1:length(evtStateLogEntries))
                    %get body position in view frame
                    entry = evtStateLogEntries(j);

                    %get body axes in view frame
                    rotMatBodyToInertial = entry.steeringModel.getBody2InertialDcmAtTime(entry.time, entry.position, entry.velocity, entry.centralBody);

                    [~, ~, ~, rotMatToInertial12] = viewInFrame.getOffsetsWrtInertialOrigin(entry.time, cartElem(j));
                    [~, ~, ~, rotMatToInertial32] = entry.centralBody.getBodyCenteredInertialFrame().getOffsetsWrtInertialOrigin(entry.time, cartElem(j));

                    rotMatsBodyToView(:,:,j) = rotMatToInertial12' * rotMatToInertial32 * rotMatBodyToInertial; %body to inertial -> inertial to inertial -> inertial to view frame
                end

                switch(evt.plotMethod)
                    case EventPlottingMethodEnum.PlotContinuous
                        %nothing

                    case EventPlottingMethodEnum.SkipFirstState
                        times = times(2:end);
                        rotMatsBodyToView = rotMatsBodyToView(:,:,2:end);

                    case EventPlottingMethodEnum.DoNotPlot
                        times = [];
                        rotMatsBodyToView = [];

                    otherwise
                        error('Unknown event plotting method: %s', evt.plotMethod);
                end

                [times,ia,~] = unique(times,'stable');
                rotMatsBodyToView = rotMatsBodyToView(:,:,ia);

                [times,I] = sort(times);
                rotMatsBodyToView = rotMatsBodyToView(:,:,I);

                if(length(times) >= 2 && all(diff(times)>0))
                    vehAttData.addData(times, rotMatsBodyToView);
                end
            end
        end
    end

    methods(Static, Access=private)
        function [times, rVects, evtColor, vVects] = parseTrajDataFromSubStateLogs(subStateLogs, j, evts)
            times = subStateLogs{j}(:,1);
            rVects = subStateLogs{j}(:,2:4);
            vVects = subStateLogs{j}(:,5:7);

            evtNum = subStateLogs{j}(1,13);
            if(isempty(evtNum) || isnan(evtNum))
                times = [];
                rVects = [];
                evtColor = [];
                vVects = [];

                return;
            end
            evt = evts(evtNum);
            evtColor = evt.colorLineSpec.color;

            switch(evt.plotMethod)
                case EventPlottingMethodEnum.PlotContinuous
                    %nothing

                case EventPlottingMethodEnum.SkipFirstState
                    times = times(2:end);
                    rVects = rVects(2:end,:);
                    vVects = vVects(2:end,:);

                case EventPlottingMethodEnum.DoNotPlot
                    times = [];
                    rVects = [];
                    vVects = [];

                otherwise
                    error('Unknown event plotting method: %s', event.plotMethod);
            end

            if(not(isempty(times)))
                [times,ia,~] = unique(times,'stable','rows');
                rVects = rVects(ia,:);
                vVects = vVects(ia,:);

                [times,I] = sort(times);
                rVects = rVects(I,:);
                vVects = vVects(I,:);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function circle_segment_test13 ( )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_TEST13 demonstrates GAUSS for quadrature rules.
%
%  Discussion:
%
%    Some recursion coefficients ALPHA and BETA are listed in Kautsky
%    and Elhay.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference
%
%    Jaroslav Kautsky, Sylvan Elhay,
%    Calculation of the Weights of Interpolatory Quadratures,
%    Numerische Mathematik,
%    Volume 40, Number 3, October 1982, pages 407-422.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_SEGMENT_TEST13\n' );
  fprintf ( 1, '  GAUSS computes the points and weights for a\n' );
  fprintf ( 1, '  Gauss quadrature rule, given the ALPHA and BETA\n' );
  fprintf ( 1, '  recursion coefficients.\n' );
%
%  Legendre rule.
%
  n = 10;

  alpha = zeros ( n, 1 );
  beta = zeros ( n, 1 );

  for i = 1 : n
    alpha(i) = 0.0;
    if ( i == 1 )
      beta(i) = 2.0;
    else
      beta(i) = 1.0 / ( 4.0 - 1.0 / ( i - 1.0 )^2 );
    end
  end

  ab(1:n,1) = alpha(1:n);
  ab(1:n,2) = beta(1:n);

  xw = gauss ( n, ab );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEGENDRE RULE\n' );
  fprintf ( 1, '  Point   Weight\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', xw(i,1), xw(i,2) );
  end
%
%  Hermite rule.
%
  n = 10;

  alpha = zeros ( n, 1 );
  beta = zeros ( n, 1 );

  for i = 1 : n
    alpha(i) = 0.0;
    if ( i == 1 )
      beta(i) = sqrt ( pi );
    else
      beta(i) = ( i - 1 ) / 2;
    end
  end

  ab(1:n,1) = alpha(1:n);
  ab(1:n,2) = beta(1:n);

  xw = gauss ( n, ab );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  HERMITE RULE\n' );
  fprintf ( 1, '  Point   Weight\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', xw(i,1), xw(i,2) );
  end
%
%  Laguerre rule.
%
  n = 10;

  alpha = zeros ( n, 1 );
  beta = zeros ( n, 1 );

  for i = 1 : n
    alpha(i) = 2.0 * i - 1.0;
    if ( i == 1 )
      beta(i) = 1;
    else
      beta(i) = ( i - 1 )^2;
    end
  end

  ab(1:n,1) = alpha(1:n);
  ab(1:n,2) = beta(1:n);

  xw = gauss ( n, ab );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LAGUERRE RULE\n' );
  fprintf ( 1, '  Point   Weight\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12g  %12g\n', xw(i,1), xw(i,2) );
  end

  return
end


text:
function [ pp, normal, seed ] = plane_normal_uniform_nd ( dim_num, seed )

%*****************************************************************************80
%
%% PLANE_NORMAL_UNIFORM_ND generates a random normal plane in ND.
%
%  Discussion:
%
%    The normal form of a plane is:
%
%      PP is a point on the plane,
%      N is a normal vector to the plane.
%
%    The point PP will be chosen at random inside the unit sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer SEED, a seed for the random number generator.
%
%    Output, real PP(DIM_NUM), a point on the plane.
%
%    Output, real NORMAL(DIM_NUM), the unit normal vector.
%

%
%  Pick PP as a random point inside the unit sphere in ND.
%
  [ pp, seed ] = ball_unit_sample_nd ( dim_num, seed );
%
%  Get values from a standard normal distribution.
%
  [ normal, seed ] = r8vec_normal_01 ( dim_num, seed );
%
%  Compute the length of the vector.
%
  norm = sqrt ( sum ( normal(1:dim_num).^2 ) );
%
%  Normalize the vector.
%
  normal(1:dim_num) = normal(1:dim_num) / norm;

  return
end