text:
function Keff = electrolyteConductivity(ce,T,param,batterySection)
% electrolyteConductivity  Evaluates the conductivity coefficients for the
% electrolyte phase. The measurement unit is [S/m]
%
%   Keff = electrolyteConductivity(ce,T,param) evaluates
%   the conductivity coefficients for the anode, separator and cathode electrolyte phase of the
%   battery. You can modify the script to meet your particular needs.
%
%   The conductivity coefficients can be evaluated in isothermal case
%   (param.TemperatureEnabled=0) or adiabatic case
%   (param.TemperatureEnabled=1 or 2).
%
%   You can modify the way that the conductivity coefficients are computed, as
%   function of electrolyte concentration and temperature. The main script
%   will pass also the param array.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

if(param.TemperatureEnabled>=1)
    switch(batterySection)
        case'p'
            Keff = param.eps_p^param.brugg_p *(1e-4*ce.*((-10.5+0.668*1e-3*ce+0.494*1e-6*ce.^2) +...
                (0.074  -1.78*1e-5*ce -8.86*1e-10*ce.^2).*T + (-6.96*1e-5+2.8*1e-8*ce).*T.^2).^2);
        case 's'
            Keff = param.eps_s^param.brugg_s *(1e-4*ce.*((-10.5+0.668*1e-3*ce+0.494*1e-6*ce.^2) +...
                (0.074  -1.78*1e-5*ce -8.86*1e-10*ce.^2).*T + (-6.96*1e-5+2.8*1e-8*ce).*T.^2).^2);
        case 'n'
            Keff = param.eps_n^param.brugg_n *(1e-4*ce.*((-10.5+0.668*1e-3*ce+0.494*1e-6*ce.^2) +...
                (0.074  -1.78*1e-5*ce -8.86*1e-10*ce.^2).*T + (-6.96*1e-5+2.8*1e-8*ce).*T.^2).^2);
    end
else
        switch(batterySection)
        case'p'
            Keff = param.eps_p^param.brugg_p *(4.1253*1e-2 + 5.007*1e-4*ce - 4.7212*1e-7*ce.^2 +1.5094*1e-10*ce.^3 -1.6018*1e-14*ce.^4);
        case 's'
            Keff = param.eps_s^param.brugg_s *(4.1253*1e-2 + 5.007*1e-4*ce - 4.7212*1e-7*ce.^2 +1.5094*1e-10*ce.^3 -1.6018*1e-14*ce.^4);
        case 'n'
            Keff = param.eps_n^param.brugg_n *(4.1253*1e-2 + 5.007*1e-4*ce - 4.7212*1e-7*ce.^2 +1.5094*1e-10*ce.^3 -1.6018*1e-14*ce.^4);
        end
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [m,q]=qrabs(q1)
%QRABS absoloute value and normalization of a real quaternions [m,q]=[q1]
%
% Inputs:
%
%     q1(4,1)  A real quaternion in the form [r, i, j, k]' where i^2=j^2=k^2=ijk=-1
%
% Outputs:
%
%     m   Magnitude of the quaternion, q1.
%     q   Normalized version of q1 with unit magnitude.

%      Copyright (C) Mike Brookes 2000-2008
%      Version: $Id: qrabs.m,v 1.1 2008/12/03 10:07:53 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=sqrt(q1'*q1);
if m>0
    q=q1/m;
else
    q=[1 0 0 0]';
end

author: decouples


text:
classdef trigspec < coeffsDiscretization
%TRIGSPEC    Fourier spectral method in coefficient space.
%   TRIGSPEC is an implementation of OPDISCRETIZATION that implements a
%   Fourier spectral method in coefficient space.
%
% See also TRIGCOLLOC.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        coeffs        % Coefficients of the operator.
        outputSpace   % The range of the operator.
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function disc = trigspec(varargin)
            disc = disc@coeffsDiscretization(varargin{:});
            % No dimension adjustment are required for TRIGSPEC.
            disc.dimAdjust = 0;
            disc.projOrder = 0;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        function tech = returnTech()
            %RETURNTECH    Return the appropriate tech to use for TRIGSPEC.
            tech = @trigtech;
        end

        % Differentiation matrices for TRIGSPEC.
        D = diffmat(N, m, flag)

        % Multiplication matrices for TRIGSPEC.
        D = multmat(N, f)

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ n_unique, a_unique ] = i4vec_sorted_unique ( n, a )

%*****************************************************************************80
%
%% I4VEC_SORTED_UNIQUE finds the unique elements in a sorted I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements in A.
%
%    Input, integer A(N), the sorted integer array.
%
%    Output, integer N_UNIQUE, the number of unique elements in A.
%
%    Output, integer A_UNIQUE[N_UNIQUE], the unique elements.
%
  n_unique = 0;
  a_unique = [];

  if ( n <= 0 )
    return;
  end

  n_unique = 1;
  a_unique = [ a_unique a(1) ];

  for i = 2 : n

    if ( a(i) ~= a_unique(n_unique) )
      n_unique = n_unique + 1;
      a_unique = [ a_unique a(i) ];
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,options] = lsqnonneg_fast(C,d,options)
%LSQNONNEG Linear least squares with nonnegativity constraints.
%
%function [x,resnorm,resid,exitflag,output,lambda] =
%lsqnonneg(C,d,x0,options)
%
%   X = LSQNONNEG(C,d) returns the vector X that minimizes NORM(C*X - d)
%   subject to X >= 0. C and d must be real.
%
%   X = LSQNONNEG(C,d,X0) uses X0 as the starting point if all(X0 > 0);
%   otherwise the default is used. The default start point is the
%   origin (the default is used when X0==[] or when only two input
%   arguments are provided).
%
%   X = LSQNONNEG(C,d,X0,OPTIONS) minimizes with the default optimization
%   parameters replaced by values in the structure OPTIONS, an argument
%   created with the OPTIMSET function.  See OPTIMSET for details.  Used
%   options are Display and TolX. (A default tolerance TolX of
%   10*MAX(SIZE(C))*NORM(C,1)*EPS is used).
%
%   [X,RESNORM] = LSQNONNEG(...) also returns the value of the squared 2-norm of
%   the residual: norm(C*X-d)^2.
%
%   [X,RESNORM,RESIDUAL] = LSQNONNEG(...) also returns the value of the
%   residual: C*X-d.
%
%   [X,RESNORM,RESIDUAL,EXITFLAG] = LSQNONNEG(...) returns an EXITFLAG that
%   describes the exit condition of LSQNONNEG.
%   If EXITFLAG is:
%     1 then LSQNONNEG converged with a solution X.
%     0 then the iteration count was exceeded. Increasing the tolerance
%       (OPTIONS.TolX) may lead to a solution.
%
%   [X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT] = LSQNONNEG(...) returns a structure
%   OUTPUT with the number of steps taken in OUTPUT.iterations and the type
%   of algorithm used in OUTPUT.algorithm.
%
%   [X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA] = LSQNONNEG(...) returns
%   the dual vector LAMBDA  where LAMBDA(i) <= 0 when X(i) is (approximately) 0
%   and LAMBDA(i) is (approximately) 0 when X(i) > 0.
%
%   See also LSCOV, SLASH.

%   L. Shure 5-8-87
%   Revised, 12-15-88,8-31-89 LS, 5-26-98 MAB.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.6 $  $Date: 1998/08/25 15:18:02 $

% Reference:
%  Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974.

if nargin<3
    %disp('Calculating optimset');

    defaultopt = optimset('display','final','TolX','10*eps*norm(C,1)*length(C)');
    if ~isreal(C) || ~isreal(d), error('C and d must be real.'); end
    options = [];
    options = optimset(defaultopt,options);
    printtype = optimget(options,'display');
    tol = optimget(options,'tolx');

    % In case the defaults were gathered from calling: optimset('fminsearch'):
    c=C;
    if ischar(tol)
        tol = eval(tol);
    end

    switch printtype
        case {'none','off'}
            verbosity = 0;
        case 'iter'
            warning('''iter'' value not valid for ''Display'' parameter for lsqnonneg');
            verbosity = 2;
        case 'final'
            verbosity = 1;
        otherwise
            error('Bad value for options parameter: ''Display''');
    end
end
tol = optimget(options,'tolx');
c=C;
if ischar(tol)
    tol = eval(tol);
end

[m,n] = size(C);
P = zeros(1,n);
Z = 1:n;
x = P';

ZZ=Z;
resid = d-C*x;
w = C'*(resid);

% set up iteration criterion
outeriter = 0;
iter = 0;
itmax = 3*n;
exitflag = 1;

% outer loop to put variables into set to hold positive coefficients
while any(Z) && any(w(ZZ) > tol)
   outeriter = outeriter + 1;
   [wt,t] = max(w(ZZ));
   t = ZZ(t);
   P(1,t) = t;
   Z(t) = 0;
   PP = find(P);
   ZZ = find(Z);
   nzz = size(ZZ);
   CP(1:m,PP) = C(:,PP);
   CP(:,ZZ) = zeros(m,nzz(2));
   z = pinv(CP)*d;
   z(ZZ) = zeros(nzz(2),nzz(1));
   % inner loop to remove elements from the positive set which no longer belong
   while any((z(PP) <= tol))
      iter = iter + 1;
      if iter > itmax
         if verbosity
            warnstr = sprintf('Exiting: Iteration count is exceeded, exiting LSQNONNEG.', ...
               '\n','Try raising the tolerance (OPTIONS.TolX).');
            disp(warnstr);
         end
         exitflag = 0;
         output.iterations = outeriter;
         resnorm = sum(resid.*resid);
         x = z;
         lambda = w;
         return
      end
      QQ = find((z <= tol) & P');
      alpha = min(x(QQ)./(x(QQ) - z(QQ)));
      x = x + alpha*(z - x);
      ij = find(abs(x) < tol & P' ~= 0);
      Z(ij)=ij';
      P(ij)=zeros(1,length(ij));
      PP = find(P);
      ZZ = find(Z);
      nzz = size(ZZ);
      CP(1:m,PP) = C(:,PP);
      CP(:,ZZ) = zeros(m,nzz(2));
      z = pinv(CP)*d;
      z(ZZ) = zeros(nzz(2),nzz(1));
   end
   x = z;
   resid = d-C*x;
   w = C'*(resid);
end

lambda = w;
resnorm = sum(resid.*resid);
output.iterations = outeriter;
output.algorithm = 'active-set using svd';

verbosity =0;

if verbosity > 0
   disp('Optimization terminated successfully.');
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function plotOutput(output)

FigNum = 1;

sol = output.result.solution;

idx = output.result.setup.auxdata.index;

Time = sol.phase.time;
th = sol.phase.state(:,idx.THETA);
phi = sol.phase.state(:,idx.PHI);
Dth = sol.phase.state(:,idx.DTHETA);
Dphi = sol.phase.state(:,idx.DPHI);

uHip = sol.phase.control(:,idx.HIP);
uAnk = sol.phase.control(:,idx.ANK);

%Print things:
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
disp(['Push-off impulse: ' num2str(sol.parameter) ' Nm']);
disp(['Objective value: ' num2str(output.result.objective)...
    ' (' output.result.setup.auxdata.cost.model ')']);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
figure(FigNum); clf;

TitleFontSize = 16;
LabelFontSize = 14;
AxisFontSize = 10;

    subplot(3,2,1)
    plot(Time,th,'LineWidth',3)
    title('Stance Angle','FontSize',TitleFontSize)
    xlabel('Time (s)','FontSize',LabelFontSize)
    ylabel('Angle (rad)','FontSize',LabelFontSize)
    set(gca,'fontsize',AxisFontSize);

    subplot(3,2,2)
    plot(Time,phi,'LineWidth',3)
    title('Swing Angle','FontSize',TitleFontSize)
    xlabel('Time (s)','FontSize',LabelFontSize)
    ylabel('Angle (rad)','FontSize',LabelFontSize)
    set(gca,'fontsize',AxisFontSize);

    subplot(3,2,3)
    plot(Time,Dth,'LineWidth',3)
    title('Stance Rate','FontSize',TitleFontSize)
    xlabel('Time (s)','FontSize',LabelFontSize)
    ylabel('Rate (rad/s)','FontSize',LabelFontSize)
    set(gca,'fontsize',AxisFontSize);

    subplot(3,2,4)
    plot(Time,Dphi,'LineWidth',3)
    title('Swing Rate','FontSize',TitleFontSize)
    xlabel('Time (s)','FontSize',LabelFontSize)
    ylabel('Rate (rad/s)','FontSize',LabelFontSize)
    set(gca,'fontsize',AxisFontSize);

    subplot(3,2,5)
    plot(Time,uHip,'LineWidth',3)
    title('Hip Torque','FontSize',TitleFontSize)
    xlabel('Time (s)','FontSize',LabelFontSize)
    ylabel('Torque (Nm)','FontSize',LabelFontSize)
    set(gca,'fontsize',AxisFontSize);

    subplot(3,2,6)
    plot(Time,uAnk,'LineWidth',3)
    title('Ankle Torque','FontSize',TitleFontSize)
    xlabel('Time (s)','FontSize',LabelFontSize)
    ylabel('Torque (Nm)','FontSize',LabelFontSize)
    set(gca,'fontsize',AxisFontSize);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function varargout = tsvread( varargin )
%[data, header, raw] = tsvread( file ) reads in text file with tab-seperated variables. default value for data is nan.
%alternative input/output option is suppluying header strings
%[col1, col2, col3, ..., header, raw] = tsvread( file, header1, header2, header3, ... )
%header is the first row (assumed to have header names) and raw is the imported text
%if a vector is supplied, this specifies the number rows to be imported.
%examples:
%[col1, col2,col3,header,raw] = tsvread( 'example.tsv', 'header1', 'header2', 'header3', 1:5 )
%will import data from example.tsv, and cols corresponding to header1,
%header2, header3, cols 1 t0 5.
%if no outputs are requested, then a portion of the rquested table is
%displayed -- good idea to see how the import and header requests are
%working!
%If there is a tsv file in local directory, then just running "tsvread" at
%the command line will read the newest tsv file and display first ten rows to screen.
%sak 9/2/11
if nargin == 0
    disp( 'importing most recent tsv file in local directory' );
    d = dir( '*.tsv' );
    [~,i] = sort( [d.datenum], 'descend' );
    fprintf( 'tsvread( ''%s'', 1:10 )\n', d(i(1)).name );
    eval( sprintf( 'tsvread( ''%s'', 1:10 )', d(i(1)).name ) );
    return;
end;

fid = fopen( varargin{1}, 'r' );
if nargout == 0
    fprintf( 'loading %s, and displaying sideways\n', varargin{1} );
end
varargin(1) = [];
stuff = textscan( fid, '%s', 'delimiter', '\n');
stuff = stuff{1};
fclose(fid);

numrows = 1:size(stuff,1);
ind = cellfun( 'isclass', varargin, 'double' );
if any( ind )
    numrows = varargin{ind};
    varargin(ind) = [];
    if numel(numrows) == 1
        numrows = 1:min(numrows, size( stuff, 1) );
    end
end
numrows = intersect( numrows, 1:size( stuff, 1 ) );
header = regexprep( regexp( stuff{1}, '[^\t]*\t', 'match' ), '\t', '' );
raw = repmat( {}, numel(numrows), numel(header) );
for i=numrows
    stuff{i}(end+1) = 9;
    tmp = regexprep( regexp( stuff{i}, '[^\t]*\t', 'match' ), '\t', '');
    raw(i,1:numel(tmp)) = tmp;
end
header = raw(1,:);
data = nan(size(raw));
for i=numrows
    for j=1:size( raw, 2 )
        if ~isempty( raw{i,j} )
            [a, count, errmsg] = sscanf( raw{i,j}, '%f' );
            if ~isempty( a )
                data(i,j) = a;
            end
        end
    end
end
%%
if numel( varargin ) == 0
    j=1:size(data, 2);
else
    j = [];
    for i=1:numel(varargin)
        j = [j, find( strncmp( header, varargin{i}, numel(varargin{i})) )];
    end
end
if nargout == 0
    disp( [ strvcat( header(j)), num2str( data(numrows,j)' ) ] );
    return;
end
if numel(varargin)==0
    varargout = {data, header, raw};
    varargout = varargout(1:nargout);
    return;
end
varargout = {};
for i=j
    varargout{end+1} = data( :, i );
end
varargout{end+1} = header(:,j);
varargout{end+1} = raw(:,j);

author: Sable

repo: mcbench-benchmarks


text:
%CDATS Support routine for checking datasets
%
% [B,C,LABLIST,P] = CDATS(A,REDUCE)
%
% INPUT
%   A        Dataset or double
%   REDUCE   0/1, Reduce A to labeled samples (1) or not (0, default), optional.
%
% OUTPUT
%   B        Dataset
%   C        Number of classes
%   LABLIST  Label list of A
%   P        Priors
%
% DESCRIPTION
% This routine supports dataset checking and conversion for mappings
% and density estimators. If A is double it is converted to a one-class
% dataset with all labels set to 1. The same holds if A is an entirely
% unlabeled dataset. The label list of A is returned in LABLIST, but
% for multi-class datasets LABLIST is set to 1. C is the number of classes
% in A. The priors are returned in P (length C), the dataset itself in B.
% If REDUCE is 1, all unlabeled objects are removed.
%
% If A has soft labels, B = A.
% If A does not have labels but targets: C = 1, LABLIST = [], P = 1.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS

function [a,c,lablist,p] = cdats(a,red)

	if (nargin < 2), red = 0; end

	if (isdataset(a) & islabtype(a,'targets'))
		c = 1; lablist = []; p = 1;
	elseif (isdataset(a) & islabtype(a,'soft'))
		c = getsize(a,3);
		lablist = getlablist(a);
		if nargout > 3, p = getprior(a); end
	else
		if ~isvaldfile(a)
			a = prdataset(a);
      % we decided to allow doubles instead of datasets where possible.
      % However, if the user calls cdats we assume he needs a dataset
      % definitely.
		end
		c = getsize(a,3);
		if nargout > 3 % avoid unnecessary warnings
			p = getprior(a);
		end
		if (c == 0)
			a = setlabels(a,1);
			p = 1;
			c = 1;
			lablist = 1;
			prwarning(4,'Dataset unlabeled: all objects will be used')
		elseif (c == 1)
			p = 1;
			if (red == 1)
				a = seldat(a); % remove unlabeled objects
			end
			lablist = getlablist(a);
		else
			prwarning(4,'Multiclass dataset will be combined using priors')
			if (red == 1)
				a = seldat(a); % remove unlabeled objects
			end
			lablist = 1;
		end
		isvaldfile(a,1,1); % at least one object per class
	end;
return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function indices = matRad_convRtssContours2Indices(structure,ct)
% matRad function to convert a polygon segmentation from an rt structure
% set into a binary segmentation as required within matRad's cst struct
%
% call
%   indices = matRad_convRtssContours2Indices(contPoints,ct)
%
% input
%   structure:      information about a single structure
%   ct:             matRad ct struct where the binary segmentations will
%                   be aligned to
%
% output
%   indicies:       indices of voxels of the ct cube that are inside the
%                   contour
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

voiCube = zeros(ct.cubeDim);

% loop over all closed contour items
for i = 1:size(structure.item,2)

    if ~isempty(structure.item(i).points)

        dicomCtSlicePos = unique(structure.item(i).points(:,3));

        if numel(dicomCtSlicePos) > 1
            error('Contour defined over multiple planes\n');
        end

        round2 = @(a,b) round(a*10^b)/10^b;
        dicomCtSliceThickness = ct.dicomInfo.SliceThickness(round2(ct.dicomInfo.SlicePositions,1)==round2(dicomCtSlicePos,1));

        coords1 = interp1(ct.x,1:ct.cubeDim(2),structure.item(i).points(:,1),'linear','extrap');
        coords2 = interp1(ct.y,1:ct.cubeDim(1),structure.item(i).points(:,2),'linear','extrap');

        binIn = poly2mask(coords1,coords2,ct.cubeDim(1),ct.cubeDim(2));

        slicesInMatradCt = find(dicomCtSlicePos+dicomCtSliceThickness/2 > ct.z & dicomCtSlicePos-dicomCtSliceThickness/2 <= ct.z);

        % loop over all slices in matRad ct
        for j = 1:numel(slicesInMatradCt)
            voiCube(:,:,slicesInMatradCt(j)) = voiCube(:,:,slicesInMatradCt(j)) | binIn;
        end

    end

end

indices = find(voiCube(:));

repo: matRad


text:
%GETCOST Get classification cost matrix
%
%  [COST,LABLIST] = GETCOST(A)
%
% Returns the classification cost matrix as defined for the dataset A.
% An empty cost matrix is interpreted as equal costs for misclassification.
% In that case COST = ONES(C+1) - EYE(C+1) is returned, if C is the number
% of classes. Row C+1 and column C+1 in COST refer to unlabeld objects.
% In LABLIST the class labels are returned.
%
% If A has target labels an error is returned as in that case no classes
% and thereby no costs are defined.
%
% See DATASETS


text:
function [ lambda, x ] = jacobi_iterate ( n, a )

%*****************************************************************************80
%
%% JACOBI_ITERATE applies the Jacobi eigenvalue iteration to a symmetric matrix.
%
%  Discussion:
%
%    I had to modify the code, in order to avoid cases where the
%    off-diagonal element was not exactly zero, but very very close.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real A(N,N), a symmetric matrix.
%    On output, the matrix has been overwritten by an approximately
%    diagonal matrix, with the eigenvalues on the diagonal.
%
%    Output, real LAMBDA(N), the computed eigenvalues.
%
%    Output, real X(N,N), the computed eigenvector matrix.
%
  eps = 0.00001;
  maxit = 100;

  error_frobenius = r8mat_is_symmetric ( n, n, a );

  if ( eps < error_frobenius )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_ITERATE - Fatal error!\n' );
    fprintf ( 1, '  The input matrix is not symmetric.\n' );
    error ( 'JACOBI_ITERATE - Fatal error!' );
  end

  b(1:n,1:n) = a(1:n,1:n);

  norm_fro = r8mat_norm_fro ( n, n, b );

  x = identity ( n, n );

  for it = 1 : maxit

    for i = 1 : n
      for j = 1 : i - 1

        if ( eps * norm_fro < abs ( b(i,j) ) + abs ( b(j,i) ) )

          u = ( b(j,j) - b(i,i) ) / ( b(i,j) + b(j,i) );

          t = r8_sign ( u ) / ( abs ( u ) + sqrt ( u * u + 1.0 ) );
          c = 1.0 / sqrt ( t * t + 1.0 );
          s = t * c;
%
%  A -> A * Q.
%
          for k = 1 : n
            t1 = b(i,k);
            t2 = b(j,k);
            b(i,k) = t1 * c - t2 * s;
            b(j,k) = t1 * s + t2 * c;
          end
%
%  A -> Q' * A
%
          for k = 1 : n
            t1 = b(k,i);
            t2 = b(k,j);
            b(k,i) = c * t1 - s * t2;
            b(k,j) = s * t1 + c * t2;
          end
%
%  X -> Q' * X
%
          for k = 1 : n
            t1 = x(k,i);
            t2 = x(k,j);
            x(k,i) = c * t1 - s * t2;
            x(k,j) = s * t1 + c * t2;
          end

        end

      end
    end
%
%  Test the size of the off-diagonal elements.
%
    sum2 = 0.0;
    for i = 1 : n
      sum2 = sum2 + sum ( abs ( b(i,1:i-1) ) );
    end

    if ( sum2 <= eps * ( norm_fro + 1.0 ) )
      break
    end

  end

  lambda = diag ( b );

  return
end


text:
% setup_paths() sets up paths for the matlabserver code to include subdirectories, etc.
% This file is sourced by the serverInit() and slaveInit() functions.  You can add your
% own paths here if you like to ensure that they're added to the server and slaves.
%
% The returned variable 'rootpath' is the path of the setup_paths.m file
% (i.e., the matlabserver/ directory) which can be used to reference other
% internal paths by the caller.
function rootpath = setup_paths()
p=which('setup_paths');
i=strfind(p, 'setup_paths.m');

if (isempty(i))
  error('Couldn''t determine path of setup_paths.m');
end

p=p(1:i-1);
addpath(p); % add root path in case cwd changes.

rootpath=p;

author: huashiyiqike

repo: LSTM-MATLAB


text:
function makenoiseband(mrSESSION, fdrift, fresp, fcard)
% makenoiseband(mrSESSION, fdrift, fresp, fcard)
% PURPOSE: Provides input for mrVista function AddNoiseBand
% Works in Hz rather than cycles per scan : converts between the two
% using the correct TR taken from mrSESSION.functionals(1).framePeriod
% and nFrames
%
% Thomas Ferree @ UCSF
% Created 1/23/2005

% extract necessary parameters from mrSESSION
TR = mrSESSION.functionals(1).framePeriod;
nFrames = mrSESSION.functionals(1).nFrames;
df = 1/(TR*nFrames);

% compute discrete frequency indices from input

if length(fdrift) > 0
    drift1 = round(fdrift(1)/df) + 1;
    drift2 = round(fdrift(2)/df) + 1;
    driftband = [drift1:drift2];
    if drift1 < 1 | drift2 < 1 | drift1 > (fix(nFrames/2)+1) | drift2 > (fix(nFrames/2)+1)
        error('Drift band is out of range.');
    end
else
    driftband = [];
end

if length(fresp) > 0
    resp1 = round(fresp(1)/df) + 1;
    resp2 = round(fresp(2)/df) + 1;
    respband = [resp1:resp2];
    if resp1 < 1 | resp2 < 1 | resp1 > (fix(nFrames/2)+1) | resp2 > (fix(nFrames/2)+1)
        error('Respiration band is out of range.');
    end
else
    respband = [];
end

if length(fcard) > 0
    card1 = round(fcard(1)/df) + 1;
    card2 = round(fcard(2)/df) + 1;
    cardband = [card1:card2];
    if card1 < 1 | card2 < 1 | card1 > (fix(nFrames/2)+1) | card2 > (fix(nFrames/2)+1)
        error('Cardiac band is out of range.');
    end
else
    cardband = [];
end

% out of sequence errors

if length(fdrift) > 0 & length(fresp) > 0 & drift2 > resp1
    error('Definition of drift and respiration bands is mixed up.');
end

if length(fresp) > 0 & length(fcard) > 0 & resp2 > card1
    error('Definition of respiration and cardiac bands is mixed up.');
end

% compute noise band (frequencies to be kept)
allband = 1:nFrames/2+1;
dropband = union(union(driftband,respband),cardband);
noiseband = setdiff(allband,dropband);

% call mrVista function AddNoiseBand
AddNoiseBand(noiseband);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function d=sqdist(a,b)
% SQDIST - computes squared Euclidean distance matrix
%          computes a rectangular matrix of pairwise distances
% between points in A (given in columns) and points in B

% NB: very fast implementation taken from Roland Bunschoten

aa = sum(a.*a,1); bb = sum(b.*b,1); ab = a'*b;
d = abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);

repo: REST


text:
function [ xdp, ydp ] = dif_deriv_table ( nd, xd, yd )

%*****************************************************************************80
%
%% DIF_DERIV_TABLE computes the divided difference table for a derivative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer ND, the size of the input table.
%
%    Input, real XD(ND), the abscissas for the divided
%    difference table.
%
%    Input, real YD(ND), the divided difference table.
%
%    Output, real XDP(ND-1), the abscissas for the divided
%    difference table for the derivative.
%
%    Output, real YDP(ND-1), the divided difference
%    table for the derivative.
%

%  Using a temporary copy of the difference table, shift the
%  abscissas to zero.
%
  xd_temp(1:nd) = xd(1:nd);
  yd_temp(1:nd) = yd(1:nd);

  [ xd_temp, yd_temp ] = dif_shift_zero ( nd, xd_temp, yd_temp );
%
%  Construct the derivative.
%
  xdp(1:nd-1) = 0.0;

  for i = 1 : nd - 1
    ydp(i) = i * yd_temp(i+1);
  end

  return
end


text:
%GOBJ_ELLIPSE Create ellipse geometry object.
%
%   [ GOBJ ] = GOBJ_ELLIPSE( P, RX, RY, TAG ) Creates an ellipse
%   geometry object. Accepts the following input parameters.
%
%       Parameter   Value/{Default}           Description
%       -----------------------------------------------------------------------------------
%       p           array  {[0 0]}            Coordinates of center point
%       rx          scalar {1}                Radius along x-axis
%       ry          scalar {0.5}              Radius along y-axis
%       tag         string {E1}               Geometry object tag/name

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function dohyst(indir,outdir,hmult,nthresh)
% function dohyst(indir,outdir,hmult,nthresh)
%
% Read pb files from indir, apply hysteresis thresholding, and
% write the resulting pb files into outdir.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% April 2003

if nargin<3, hmult=1/3; end
if nargin<4, nthresh=100; end

iids = imgList('test');
unused = mkdir(outdir);
for i = 1:numel(iids),
  iid = iids(i);
  pbold = double(imread(sprintf('%s/%d.bmp',indir,iid)))/255;
  thresh = linspace(1/nthresh,1-1/nthresh,nthresh);
  fwrite(2,sprintf('%3d/%d ',i,numel(iids)));
  progbar(0,nthresh);
  pbnew = zeros(size(pbold));
  for i = 1:nthresh,
    progbar(i,nthresh);
    [r,c] = find(pbold>=thresh(i));
    if numel(r)==0, continue; end
    b = bwselect(pbold>hmult*thresh(i),c,r,8);
    pbnew = max(pbnew,b*thresh(i));
  end
  imwrite(pbnew,sprintf('%s/%d.bmp',outdir,iid),'bmp');
end

author: Cloud-CV

repo: object-proposals


text:
function resaveHTK2KaldiArk(flist, uttID, outdir, nUttPerArk, archive_name_prefix)
my_mkdir(outdir);
nUtt = length(flist);
nArk = ceil(nUtt/nUttPerArk);

for ai = 1:nArk
    PrintProgress(ai, nArk, 1, 'Resave features in Kaldi archive format');

    startIdx = (ai-1)*nUttPerArk + 1;
    stopIdx = min(nUtt, ai*nUttPerArk);
    ark_name = [outdir '/' archive_name_prefix '.' num2str(ai)];
    saveHTK2KaldiArk(flist(startIdx:stopIdx), uttID(startIdx:stopIdx), ark_name);
end
end

function saveHTK2KaldiArk(flist, uttID, ark_name)
for i=1:length(flist)
    spec{i} = readHTK(flist{i});
end
writeKaldiArchiveText(ark_name, uttID, spec);
gzip(ark_name);        delete(ark_name);
end

author: singaxiong

repo: SignalGraph


text:
%                                                 y = triumtriu(r,u,K)
% TRIUMTRIU  Computes y = r * u
%   Both r and u should be upper triangular.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

function y = triumtriu(r,u,K) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 2 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.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')


text:
function TestEnhanceNet()
addpath('local');

dnn_files{1} = 'nnet/EnhanceRegression.noCMN.DeltaByEqn.MbSize40.U28539.771-LSTM-2048-771.L2_3E-4.LR_3E-3/nnet.itr37.LR1E-5.CV2453.828.mat';
dnn_files{2} = 'nnet/EnhanceRegression.noCMN.DeltaByEqn.MbSize20.U28539.771-LSTM-2048-771.L2_3E-4.LR_1E-4/nnet.itr62.LR4.97E-8.CV2439.245.mat';
dnn_files{3} = 'nnet/EnhanceGaussian.noCMN.init1.expDecay.Decay0.999.DeltaByEqn.MbSize40.U28539.771-LSTM-2048-771.L2_0.LR_3E-3/nnet.itr26.LR3.62E-5.CV756.273.mat';
IsMaskNet = [0 0 0];

testset = 2;
T60 = [0.01 0.3 0.6];
noise = {'buccaneer2', 'pink', 'm109', 'hfchannel', 'factory2', 'factory1', 'destroyerops', 'babble', 'volvo', 'leopard', 'buccaneer1', 'destroyerengine', 'f16', 'machinegun', 'white'};
% noise = {'f16', 'machinegun', 'white'};
SNR = [10 0 -10];
DEBUG = 0;
measures = {'PESQ', 'FWSEQ', 'CD', 'IS', 'LLR', 'WSS'};
useGPU = 0;

for t60_i = 2%:length(T60)
    for n_i = 12%:length(noise)
        for s_i = 3%:length(SNR)
            TestEnhanceNetByCategory(dnn_files, testset, T60(t60_i), noise{n_i}, SNR(s_i), IsMaskNet, measures, useGPU, DEBUG)
        end
    end
end

end


text:
% plot the ROC curve with associated level set samples and true positive
% and false positive rate calculations

% Zoya Bylinskii, April 2016
% linked to: "What do different evaluation metrics tell us about saliency models?"

function visualize_AUC(salMap,fixations)
% salMap is the saliency map
% fixations is a binary map of fixation locations

npoints = 10; % number of points to sample on ROC curve

% prepare the color map for the correctly detected and missed fixations
colmap = fliplr(colormap(jet(npoints)));
G = linspace(0.5,1,20)';
tpmap = horzcat(zeros(size(G)),G,zeros(size(G))); % green color space
fpmap = horzcat(G,zeros(size(G)),zeros(size(G))); % red color space

% compute AUC-Judd
heatmap = im2double(salMap);
[score,tp,fp,allthreshes] = AUC_Judd(heatmap, fixations);

N = ceil(length(allthreshes)/npoints);
allthreshes_samp = allthreshes(1:N:end);

heatmap_norm = (heatmap-min(heatmap(:)))/(max(heatmap(:))-min(heatmap(:)));

salMap_col = makeLevelSets(heatmap, allthreshes_samp, colmap);
figure; subplot(1,2,1); imshow(salMap_col);

% plot the ROC curve
tp1 = tp(1:N:end); fp1 = fp(1:N:end);
subplot(1,2,2); plot(fp,tp,'b'); hold on;
for ii = 1:npoints
     plot(fp1(ii),tp1(ii),'.','color',colmap(ii,:),'markersize',20); hold on; axis square;
end
title(sprintf('AUC: %2.2f',score),'fontsize',14);
xlabel('FP rate','fontsize',14); ylabel('TP rate','fontsize',14);

% plot the level sets, one per subplot
figure('name','saliency map level sets');
nplot = floor(npoints/2); % plot every other level set
for ii = 1:nplot
   %temp = heatmap_norm>=allthreshes_samp(ii);   % plot every level set
   temp = heatmap_norm>=allthreshes_samp(2*ii);  % plot every other level set
   temp2 = zeros(size(temp,1),size(temp,2),3);
   temp2(:,:,1) = temp*colmap(2*ii,1); temp2(:,:,2) = temp*colmap(2*ii,2); temp2(:,:,3) = temp*colmap(2*ii,3);
   subplottight(1,nplot,ii,0.05); imshow(temp2);
end

% plot the fixations falling within and outside of each level set
figure('name','true positives and false negatives');
for ii = 1:nplot
   % temp = heatmap_norm>=allthreshes_samp(ii); % plot every level set
   temp = heatmap_norm>=allthreshes_samp(2*ii); % plot every other level set
   res = fixations.*temp;
   res_neg = fixations.*(1-temp);
   subplottight(1,nplot,ii,0.05); imshow(temp); hold on;
   [J,I] = ind2sub(size(res),find(res==1));
   scatterPlotHeatmap(I,J,tpmap);
   [J,I] = ind2sub(size(res_neg),find(res_neg==1));
   scatterPlotHeatmap(I,J,fpmap);
end

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function identify_redcells_sourcery(db, ops0)

ops = build_ops3(db, ops0);

redcells = [];

for i = 1:length(ops.planesToProcess)

    iplane  = ops.planesToProcess(i);
    %%
    try
        fname = sprintf('%s/F_%s_%s_plane%d.mat', ops.ResultsSavePath, ops.mouse_name, ops.date, iplane);
        dat = load(fname);
        while isfield(dat, 'dat')
            dat = dat.dat;
        end
    catch
        fname = sprintf('%s/F_%s_%s_plane%d_proc.mat', ops.ResultsSavePath, ops.mouse_name, ops.date, iplane, ops.Nk);
        dat = load(fname);
        while isfield(dat, 'dat')
            dat = dat.dat;
        end
    end

    % subtract bleedthrough of green into red channel
    % if there's a specific green channel (from a red+green short session)
    % using non-rigid linear regression
    if isfield(dat.ops,'mimgGREEN')
        mimgG = dat.ops.mimgGREEN;
    else
        mimgG = dat.ops.mimg1;
    end
    mimgR = dat.ops.mimgRED;
    mimgR = mimgR(dat.ops.yrange, dat.ops.xrange);
    mimgG = mimgG(dat.ops.yrange, dat.ops.xrange);
    [Ny, Nx] = size(mimgR);

    % regression
    nblks               = 3;
    yB                  = round(linspace(1,Ny,nblks+1));
    xB                  = round(linspace(1,Nx,nblks+1));
    [yBL,xBL] = MakeQuadrants(yB,xB);
    msk                 = zeros(Ny,Nx,length(yB));
    for j = 1:length(yBL)
        xg              = mimgG(yBL{j},xBL{j});
        xr              = mimgR(yBL{j},xBL{j});
        A               = polyfit(xg,xr,1);
        msk(:,:,j)      = QuadrantMask(yBL,xBL,Ny,Nx,j);
        gw(j)           = A(1);
    end
    msk                 = bsxfun(@times,msk,1./sum(msk,3));
    pixweight           = zeros(Ny,Nx);
    for j = 1:length(yBL)
        pixweight = pixweight + msk(:,:,j)*gw(j);
    end
    mimgR0              = mimgR - pixweight.*mimgG;

    % save to dat for GUI
    dat.ops.mimgREDcorrected = zeros(dat.ops.Ly, dat.ops.Lx);
    dat.ops.mimgREDcorrected(dat.ops.yrange, dat.ops.xrange) = mimgR0;
    %

    %%
    %%%% compute red pixels in cell and in area around cell
    % (exclude pixels from other cells (cellPix))
    % (use surround neuropil masks to compute red around cell)
    ops.ratioNeuropil     = 4;
    ops.minNeuropilPixels = 80;
    [dat.stat, cellPix, ~]       = createCellMasks(dat.stat, Ny, Nx);
    [~, neuropMasks] = createNeuropilMasks(ops, dat.stat, cellPix);

    %
    redSum = zeros(length(dat.stat),2);
    for k = 1:numel(dat.stat)
        ipix                 = dat.stat(k).ipix;
        rpix                 = mimgR0(ipix) .* dat.stat(k).lam/sum(dat.stat(k).lam);
        % rescale by lam? .*(dat.stat(j).lambda'/sum(dat.stat(j).lambda));
        extpix               = squeeze(neuropMasks(k,:,:));
        extpix               = find(extpix~=0);
        ext_rpix             = mimgR0(extpix);
        redSum(k,:)          = [sum(rpix) mean(ext_rpix)];
    end
    redSum = redSum - min(redSum(:));

    % set threshold for redpix
    ops.redthres = getOr(ops, 'redthres', 1.5);
    ops.redmax   = getOr(ops, 'redmax', 1.25);

    rrat = redSum(:,1)./(redSum(:,2)+redSum(:,1));
    redcell  = rrat > nanmean(rrat) + ops.redthres*nanstd(rrat);
    notred   = rrat <= nanmean(rrat) + ops.redmax*nanstd(rrat);

    fprintf('plane %d  reds %d\n',iplane,sum(redcell(:)));

    %%
    for j = 1:length(dat.stat)
        dat.stat(j).redcell = redcell(j);
        dat.stat(j).redprob = rrat(j);
        dat.stat(j).notred  = notred(j);
    end

    dat.ops.redthres = ops.redthres;
    dat.ops.redmax   = ops.redmax;
    save(fname, '-struct', 'dat')
end

%%
if 0
    [~,ix] = sort(rrat, 'descend');
    for k = ix(:)'
        clf;
        cent = round(dat.stat(k).med);
        yl = cent(1) + [-20:20];
        yl(yl<1 | yl > Ny) = [];
        xl = cent(2) + [-20:20];
        xl(xl<1 | xl > Nx) = [];
        imagesc(mimgR0)
        hold all;
        plot(cent(2),cent(1),'k*');
        axis([xl(1) xl(end) yl(1) yl(end)]);
        title([rrat(k) redcell(k)]);
        drawnow;
        pause;
    end
end
    %

author: cortex-lab


text:
function [dist1, dist2] = distance2(i, bgevent, ac, newcat)
    % calculates the distance in [km] between two eqs
    % precise version based on Raesenbergs Program
    % the calculation is done simultaniously for the biggest event in the
    % cluster and for the current event
    %
    % A. Allmann

    global err derr

    if isempty(err)
        err=1.5; % from InitVariables --epicenter error
    end

    if isempty(derr)
        derr=2; % from InitVariables --depth error
    end

    % columns
    LAT = 2;
    LON = 1;
    DEP = 3;

    % if any of these fail, then we are dealing with multiple events simultaneously, and must take that into account
    assert(numel(i)==1)
    assert(numel(bgevent) == 1)

    evs = newcat.XYZ([i,bgevent,reshape(ac,1,[])],:);

    evs(:,[LAT,LON]) = deg2rad(evs(:,[LAT,LON])); % convert lats & lons to radians

    i_ev = evs(1,:);
    big_ev = evs(2,:);
    ac_ev = evs(3,:);

    pi2 = pi/2;
    flat= 0.993231;

    tana = flat .* tan([i_ev(LAT); big_ev(LAT)]);
    acol = pi2 - atan(tana);

    tanb = flat * tan(ac_ev(LAT));
    bcol = pi2 - atan(tanb);

    diflon = (ac_ev(LON) - [i_ev(LON),  big_ev(LON)]);

    cosdel = [ ...
        (sin(acol(1))*sin(bcol)) .* cos(diflon(1)) + (cos(acol(1))*cos(bcol)) , ...
        (sin(acol(2))*sin(bcol)) .* cos(diflon(2)) + (cos(acol(2))*cos(bcol)) ];

    colat = pi2 - (ac_ev(LAT) + [i_ev(LAT), big_ev(LAT)] ./2);

    radius = 6371.227*(1+(3.37853e-3)*(1/3-((cos(colat)).^2)));
    r = acos(cosdel) .* radius;            % epicenter distance

    r = r - 1.5 * err;               %influence of epicenter error

    r(r<0) = 0;

    % depth distance
    z = abs(ac_ev(DEP) - [i_ev(DEP), big_ev(DEP)] );
    z = z - derr;
    z(z<0)=0;

    r = sqrt(z .^ 2 + r.^2);                   %hypocenter distance
    dist1 = r(:,1);           %distance between eqs
    dist2 = r(:,2);

author: CelsoReyes


text:
function mean = arcsin_mean ( a )

%*****************************************************************************80
%
%% ARCSIN_MEAN returns the mean of the Arcsin PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the CDF.
%    A must be positive.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = 0.0;

  return
end


text:
function coef=ref_dgtns_1(f,gamma,V)
%REF_DGT_1 Reference DGTNS using P.Prinz algorithm
%   Usage:  c=ref_dgtns_1(f,gamma,V);
%

a=V(1,1);
b=V(2,2);
r=-V(2,1);
L=size(gamma,1);
M=L/b;
N=L/a;
W=size(f,2);

c=gcd(a,M);
p=a/c;
q=M/c;
d=N/q;

if r==0
  % The grid is rectangular. Call another reference algorithm.
  coef=zeros(M*N,W);
  coef(:)=ref_dgt(f,gamma,a,M);

  return;
end;

% We can now assume that the grid is truly nonseperable,
% and so d>1.

% Conjugate.
gammac=conj(gamma);

% Level 2: Block diagonalize, and use that some blocks are
% the same up to permutations.

p1=stridep(M,L);
p2=stridep(N,M*N);

% Get shift offsets for stage 2 of the algorithm.
[mll,all]=shiftoffsets(a,M);

% Step 1: Permute
s1 = f(p1,:);

% Step 2: Multiply by DG'
s2=zeros(M*N,W);

% Do interpreter-language-optimized indexing.
[n_big,m_big]=meshgrid(0:N-1,0:b-1);
base=m_big*M-n_big*a+L;
base=base.';

% Work arrays.
work=zeros(b,M/c*W);
wk=zeros(N,b);
wkrect=zeros(N,b);

% Create fixed modulation matrix (Does not change with k)
fixedmod=zeros(N,b);
for n=0:N-1
  fixedmod(n+1,:)=exp(2*pi*i*r*n/L*(0:M:L-1));
end;

% This loop iterates over the number of truly different wk's.
for ko=0:c-1

  % Create the wk of the rectangular-grid case.
  wkrect(:)=gammac(mod(base+ko,L)+1);

  % Create wk of skewed case.
  wk=(fixedmod.*wkrect);

  % Setup work array.
  for l=0:M/c-1
    k=ko+l*c;
    rowmod=exp(2*pi*i*r*(0:b-1)/b*all(l+1)).';

    work(:,l*W+1:(l+1)*W)=circshift(rowmod.*s1(1+(ko+l*c)*b:(ko+l*c+1)*b,:),-mll(l+1));
  end;

  % Do the actual multiplication,
  work2=wk*work;

  % Place the result correctly.
  for l=0:M/c-1
    k=ko+l*c;
    kmod=exp(2*pi*i*r*(0:N-1)*k/L).';
    colmod=exp(2*pi*i*r*(0:N-1)/b*mll(l+1)).';
    doublefac=exp(-2*pi*i*r/b*all(l+1)*mll(l+1));

    s2(1+(ko+l*c)*N:(ko+l*c+1)*N,:)=doublefac*colmod.*kmod.*circshift(work2(:,l*W+1:(l+1)*W),all(l+1));
  end;

end;

% Step 3: Permute again.
coef = s2(p2,:);

% Apply fft.
for n=1:N
  coef((n-1)*M+1:n*M,:)=fft(coef((n-1)*M+1:n*M,:));
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function hRD = plugin_rd_tx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.fmcomms5.common.plugin_rd('ZC702', 'Tx');
AnalogDevices.fmcomms5.zc702.add_io(hRD,'Tx');

author: analogdevicesinc


text:
function fd2d_heat_steady_test ( )

%*****************************************************************************80
%
%% FD2D_HEAT_STEADY_TEST tests the FD2D_HEAT_STEADY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 August 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD2D_HEAT_STEADY_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FD2D_HEAT_STEADY library.'

  fd2d_heat_steady_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD2D_HEAT_STEADY_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function p = jaccard_similarity(A,B)
%% jaccard similarity
% This function calculates jaccard similarity index of inputs arrays A and
% B. The formula to find the Index is (number of entries in both sets) / (number of entries in either set) * 100
% The higher the percentage, the more similar the two arrays.
% For this, each of input arrays is modified by removing its same entries
% (except on them), then number of common entries between two new arrays is
% calculated by comparing them.

modified_A = unique(A);
modified_B = unique(B);

length_mA = length(modified_A);
length_mB = length(modified_B);
common_number = 0;                    %initialize the number of common entries

if length_mA <= length_mB
    X = modified_A;
    Y = modified_B;
else
    X = modified_B;
    Y = modified_A;
end

for i = 1:length(X)
    for j = 1:length(Y)
        if X(i) == Y(j)
           common_number = common_number + 1;
        end
    end
end

total_number = length_mA + length_mB - common_number;
p = (common_number/total_number)*100;
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
%% EBSD Analysis
%
%%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function str = cell2str(theCell)
% CELL2STR - Convert 1-D cells to a string
%
%   STR = CELL2STR(THECELL)
%
% Converts a 1-D cell to a string.
%
% Example:
%   A = {'test','test2','test3'};
%   str = cell2str(A)
%
%       produced str =
%    '{ 'test','test2','test3' }'
%
 %1-dim cells only, only chars and matricies

if isempty(theCell), str = '{}'; return; end;

str = '{ ';
for i=1:length(theCell),
        if ischar(theCell{i})
                str = [str '''' theCell{i} ''', '];
        elseif isnumeric(theCell{i}),
                str = [str mat2str(theCell{i}) ', '];
        end;
end;
str = [str(1:end-2) ' }'];

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function v = ksvdver(history)
%KSVDVER K-SVD toolbox version information.
%  KSVDVER displays the current KSVD toolbox version information.
%
%  KSVDVER('history') also displays history information about the previous
%  versions of the KSVD toolbox and their change logs.
%
%  V = KSVDVER returns the version number of the current KSVD toolbox, and
%  does not display any information.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  August 2009

ver = 12;

if (nargout>0)
  if (nargin>0)
    error('Invalid number of parameters.');
  end
  v = ver;

else

  if (nargin==0 || (nargin==1 && strcmpi(history,'history')))

    disp(' ');
    disp('---------------------------------');
    printf('    KSVD Toolbox version %d       ',ver);
    disp('---------------------------------');
    disp(' ');

  else
    error('Unknown parameters.');
  end

  if (nargin>0)

    %% changes only displayed since first public release %%

%     disp(' ');
%     disp('KSVD Toolbox version 1, 24.10.07');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Initial release.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('KSVD Toolbox version 2, 13.12.07');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. KSVD: support for sparsity and error based coding.');
%     disp('  2. New OMPDENOISE function for generic OMP-based denoising.');
%     disp('  3. ''memusage'' parameter introduced for controlling memory usage.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. KSVD: when initializing the dictionary with a set of signals, ');
%     disp('     only selects non-zero ones.');
%     disp('  2. KSVD: only negates atoms with first element negative, instead of');
%     disp('     multiplying each atom by its sign (which nullified atoms whose');
%     disp('     first element was zero).');
%     disp('  3. KSVD: tracks training signals used to replace unused atoms, so');
%     disp('     the same one is not used twice.');
%     disp(' ');
%     disp('Internal changes:');
%     disp('  1. Eliminated the use of function cleardict in KSVD, following bugfix (3).');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('KSVD Toolbox version 3, 16.01.08');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. Option for specifying a list of initial signals for the dictionary using');
%     disp('     the parameter initdict.')
%     disp('  2. Changed error computation from Frobenius norm to MSE.');
%     disp('  3. Added the VERBOSE parameter for controlling printed output.');
%     disp('  4. Added the muthresh parameter to KSVD.');
%     disp('  5. Added parameter maxval to OMPDENOISE and KSVDDENOISE.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. Added the sparsecode method for cleaner code.');
%     disp('  2. Re-introduced the use of function cleardict.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. Function cleardict now checks correlation between atoms in absolute value.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('KSVD Toolbox version 4, 27.01.08');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. Added functions OMPDENOISE2 and OMPDENOISE3 for faster denoising of');
%     disp('     2-D and 3-D signals.');
%     disp('  2. Added gain and lambda parameters to OMPDENOISE and KSVDDENOISE.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('KSVD Toolbox version 5, 03.02.08');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. KSVD: atom update is now in random order rather than sequential.');
%     disp('  2. KSVD: improved error computation: error is now computed after the');
%     disp('     dictionary update step, eliminating the need for an additional sparse');
%     disp('     coding step at the end of the algorithm for error computation.');
%     disp('     Computed error is now RMSE for sparsity-based minimization, and average');
%     disp('     atom count for error-based minimization.');
%     disp('  3. KSVD: cleardict now replaces every atom not used by at least 4 training');
%     disp('     signals. Atoms already replaced in the dictionary update step are');
%     disp('     excluded from this criterion.');
%     disp('  4. Added the ''exact'' parameter to KSVD/KSVDDENOISE, to allow using the exact');
%     disp('     SVD solution rather than the default quicker approximation method.');
%     disp('  5. Changed behavior of VERBOSE=3 mode in KSVD, and added VERBOSE=4 mode.');
%     disp('  6. VERBOSE in OMPDENOISE/2/3 now prints a message every VERBOSE seconds,');
%     disp('     not every VERBOSE blocks.');
%     disp('  7. OMPDENOISE: faster computation of the matrix cnt (a counter matrix');
%     disp('     containing the number of blocks covering each sample in the signal.');
%     disp('  8. KSVDDENOISE: new default values for gain and lambda.');
%     disp('  9. KSVDDENOISE: reverts to function OMPDENOISE when memusage is ''low''');
%     disp(' 10. Set a limit on the number of atoms used to code each block in OMPDENOISE');
%     disp('     and KSVDDENOISE. Added the maxatoms parameter to control this limit.');
%     disp(' 11. Improved documentation.');
%     disp(' 12. Added the KSVDVER function.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. Support for the new call syntax of OMP Toolbox v6.');
%     disp('  2. OMPDENOISE3 now computes the matrix cnt using much less memory. Instead');
%     disp('     of allocating three matrices the size of the input signal, it only');
%     disp('     allocates one such matrix.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. OMPDENOISE2/3 now compute the matrix cnt correctly.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('KSVD Toolbox version 6, 14.02.08');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. New stepsize parameter in OMPDENOISE/2/3 and KSVDDENOISE.');
%     disp('  2. ''blocksize'' parameter in OMPDENOISE/2/3 and KSVDDENOISE can now be a scalar.');
%     disp('  3. Added the Contents.m file.');
%     disp('  4. Revised VERBOSE settings in KSVDDENOISE, and slightly improved display format');
%     disp('     in KSVD.');
%     disp(' ');
%     disp('Bug fixes:');
%     disp(' ');
%     disp('  1. OMPDENOISE/2/3: fixed a bug in the computation of the matrix cnt when');
%     disp('     the signal is not square.');
%     disp(' ');
%     disp(' ');
%     disp(' ');
%
%     disp('KSVD Toolbox version 7, 18.04.08');
%     disp('---------------------------------');
%     disp(' ');
%     disp('Features:');
%     disp(' ');
%     disp('  1. OMPDENOISE now prints the estimated remaining time.');
%     disp('  2. New ''msgdelta'' parameter in KSVD and KSVDDENOISE for printing');
%     disp('     periodic status updates.');
%     disp(' ');
%     disp('Internal changes:');
%     disp(' ');
%     disp('  1. Compatible with OMP Toolbox v7.');
%     disp(' ');
%     disp(' ');
%     disp(' ');

    disp(' ');
    disp(' ');
    disp('KSVD Toolbox version 7, 18.04.08');
    disp('---------------------------------');
    disp(' ');
    disp('First public release.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('KSVD Toolbox version 8, 29.05.08');
    disp('--------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. Added ''sigma'' and ''noisemode'' parameters to OMPDENOISE and KSVDDENOISE.');
    disp('  2. Added the overcomplete DCT dictionary as an optional initial dictionary');
    disp('     in KSVDDENOISE, by setting the initdict parameter to the string ''odct''.');
    disp('  3. Improved default gain and lambda parameters in OMPDENOISE.');
    disp(' ');
    disp('Bug fixes:');
    disp(' ');
    disp('  1. KSVDDENOISE: corrected a bug which caused incorrect gain and lambda');
    disp('     parameters to be sent to OMPDENOISE in some cases.');
    disp(' ');
    disp('Internal changes:');
    disp(' ');
    disp('  1. When replacing an unused atom, the search for an example with maximal');
    disp('     error is limited to a fixed number of signals.');
    disp('  2. Improved implementation of KSVDVER.');
    disp('  3. Uses the TIMER package for time estimation.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('KSVD Toolbox version 9, 30.11.08');
    disp('--------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. OMPDENOISE and KSVDDENOISE can now return the average number of');
    disp('     non-zeros in the denoised block representations.');
    disp(' ');
    disp('Bug fixes:');
    disp(' ');
    disp('  1. Corrected an error when parameter MSGDELTA was not specified.');
    disp(' ');
    disp('Internal changes:');
    disp(' ');
    disp('  1. The COUNTCOVER function counts the number of times a pixel in a');
    disp('     multi-dimensional signal is covered by overlapping blocks.');
    disp('  2. Function CLEARDICT now computes the error norm in blocks to reduce');
    disp('     memory consumption.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('KSVD Toolbox version 10, 21.05.09');
    disp('---------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. Source code for supporting C functions now available for compilation');
    disp('     on any Matlab platform.');
    disp('  2. Added KSVDDEMO and KSVDDENOISEDEMO.');
    disp('  3. Improved performance of OMPDENOISE2, OMPDENOISE3 and KSVD.');
    disp('  4. New SHOWDICT function for displaying a dictionary of image patches.');
    disp('  5. KSVD: improved VERBOSE/MSGDELTA parameter syntax.');
    disp('  6. KSVD: messages are now displayed by default.');
    disp('  7. KSVD: only one of the parameters dictsize/initdict needs to be specified.');
    disp('  8. KSVD: eliminated the data saving option.');
    disp('  9. OMPDENOISE: messages are now displayed by default.');
    disp(' 10. KSVDDENOISE: simplified MSGDELTA parameter.');
    disp(' 11. KSVDDENOISE: default initial dictionary is now overcomplete DCT.');
    disp(' 12. KSVDDENOISE: parameter initdict can be the string ''data'' to indicate');
    disp('     initialization from training examples.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('KSVD Toolbox version 11, 03.08.09');
    disp('---------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. Added OMPDENOISE1 for fast OMP denoising of 1-D signals.');
    disp('  2. Added support for 1-D signals in KSVDDENOISE.');
    disp(' ');
    disp('Bug fixes:');
    disp('  1. Corrected a bug in OMPDENOISE/2/3 when memusage is low.');
    disp(' ');
    disp('Internal changes:');
    disp('  1. KSVD: only compute sum(data.*data) when memusage is high.');
    disp(' ');
    disp(' ');
    disp(' ');

    disp('KSVD Toolbox version 12, 24.08.09');
    disp('---------------------------------');
    disp(' ');
    disp('Features:');
    disp(' ');
    disp('  1. KSVD: significantly faster atom update computation via optimized');
    disp('     sparse-matrix code.');
    disp('  2. KSVD: reduced memory consumption by computing one column of G');
    disp('     at a time when clearing the dictionary, and using block processing');
    disp('     to compute representation error.');
    disp(' ');
    disp('Internal changes:');
    disp('  1. KSVD: when replacing an atom in cleardict, the search for an example');
    disp('     with maximal error is limited to a fixed number of signals.');
    disp('  2. KSVDDENOISEDEMO: memusage is set to ''high'' for faster performance.');
    disp(' ');
    disp(' ');

   end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
classdef ThrottlePolyModel < AbstractThrottleModel
    %AbstractSteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        throttleModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);
    end

    methods
        function throttle = getThrottleAtTime(obj, ut, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~)
            throttle = obj.throttleModel.getValueAtTime(ut);

            if(throttle < 0)
                throttle = 0.0;
            elseif(throttle > 1)
                throttle = 1.0;
            end
        end

        function enum = getThrottleModelTypeEnum(~)
            enum = ThrottleModelEnum.PolyModel;
        end

        function initThrottleModel(obj, initialStateLogEntry)
            t0 = initialStateLogEntry.time;
            obj.setT0(t0)

            if(obj.throttleContinuity)
                throttle = initialStateLogEntry.throttle;
                obj.throttleModel.constTerm = throttle;
            end
        end

        function setInitialThrottleFromState(obj, stateLogEntry, tOffsetDelta)
            t0 = stateLogEntry.time;
            obj.setT0(t0);

            obj.throttleModel.tOffset = obj.throttleModel.tOffset + tOffsetDelta;
        end

        function t0 = getT0(obj)
            t0 = obj.throttleModel.t0;
        end

        function setT0(obj, newT0)
            obj.throttleModel.t0 = newT0;
        end

        function setPolyTerms(obj, const, linear, accel)
            obj.throttleModel.constTerm = const;
            obj.throttleModel.linearTerm = linear;
            obj.throttleModel.accelTerm = accel;
        end

        function setTimeOffsets(obj, timeOffset)
            obj.throttleModel.tOffset = timeOffset;
        end

        function optVar = getNewOptVar(obj)
            optVar = SetPolyThrottleModelActionOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function [addActionTf, throttleModel] = openEditThrottleModelUI(obj, lv, useContinuity)
            output = AppDesignerGUIOutput({false});
            lvd_EditActionSetThrottleModelGUI_App(obj, lv, useContinuity, output);
            addActionTf = output.output{1};
            throttleModel = output.output{2};
        end
    end

    methods(Access=private)
        function obj = ThrottlePolyModel(throttleModel)
            obj.throttleModel = throttleModel;
        end
    end

    methods(Static)
        function model = getDefaultThrottleModel()
            throttlePoly = PolynominalModel(0,0,0,0);
            model = ThrottlePolyModel(throttlePoly);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function formatPlot(fig,titleString,subtitleStrings,xString,yString,X,Y,legendStrings,notes)
% Creates a set of 1x3 subplots for the given input data. Adds formatting
% based on the given x and y axes labels, legend, title, subtitles and
% notes.
%
% Parameters
% ----------
%     fig : [1 1] handle
%         Figure handle that the plots are created on
%
%     titleString : [1 1] string
%         Overall title shown above all subplots
%
%     subTitleStrings : [1 3] cell of strings
%         Subplot titles shown above each subplot
%
%     xString : [1 3] cell of strings
%         X axis label shown on each subplot
%
%     yString : [1 3] cell of strings
%         Y axis label shown on each subplot
%
%     X : [1 1] struct
%         Struct of X data. Fields numbered sequentially. The number of
%         fields is the number of datasets being plotted. Each field must
%         correspond to one field of Y (e.g. X.X1 will be the x-axis for
%         all values in Y.Y1)
%
%     Y : [1 1] struct
%         Struct of Y data. Fields numbered sequentially. The number of
%         fields is the number of datasets being plotted. Each field must
%         correspond to one field of Y and have a size of 3 in dimension 1
%
%     legendStrings : [: :] cell of strings
%         Strings shown in the legend. The size of dimension 1 must be the
%         size of Y.Y1 dimension 2. The size of dimension 2 must be the
%         number of fields in Y
%
%     notes : [1 :] cell of strings
%         Annotations shown below all subplots. Each column corresponds to
%         a new line of text

% Overall title for all subplots
sgtitle(titleString,'Interpreter','latex','FontWeight','bold','FontSize',12);

legendStrings1D = legendStrings(:);
legendStrings1D = legendStrings1D(~cellfun(@isempty, legendStrings1D));

% Figures
nHydro = length(fieldnames(X));
nBodiesTotal = 0;
for iHydro = 1:nHydro
    tmp1 = strcat('X',num2str(iHydro));
    tmp2 = strcat('Y',num2str(iHydro));
    [~,nBodies,~] = size(Y.(tmp2));
    nBodiesTotal = nBodiesTotal + nBodies;

    % Surge
    subplot('Position',[0.0731 0.3645 0.2521 0.4720])
    hold('on');
    box('on');
    for iBody = 1:nBodies
        plot(X.(tmp1),squeeze(Y.(tmp2)(1,iBody,:)),'LineWidth',1)
    end
    if iHydro == nHydro
        legend(legendStrings1D,'location','best','Box','off','Interpreter','none')
        title(subtitleStrings(1));
        xlabel(xString(1),'Interpreter','latex');
        ylabel(yString(1),'Interpreter','latex');
    end

    % Heave
    subplot('Position',[0.3983 0.3645 0.2521 0.4720]);
    hold('on');
    box('on');
    for iBody = 1:nBodies
        plot(X.(tmp1),squeeze(Y.(tmp2)(2,iBody,:)),'LineWidth',1);
    end
    if iHydro == nHydro
        legend(legendStrings1D,'location','best','Box','off','Interpreter','none')
        title(subtitleStrings(2));
        xlabel(xString(2),'Interpreter','latex');
        ylabel(yString(2),'Interpreter','latex');
    end

    % pitch
    subplot('Position',[0.7235 0.3645 0.2521 0.4720]);
    hold('on');
    box('on');
    for iBody = 1:nBodies
        plot(X.(tmp1),squeeze(Y.(tmp2)(3,iBody,:)),'LineWidth',1);
    end
    if iHydro == nHydro
        legend(legendStrings1D,'location','best','Box','off','Interpreter','none')
        title(subtitleStrings(3));
        xlabel(xString(3),'Interpreter','latex');
        ylabel(yString(3),'Interpreter','latex');
    end
end

% Footer
annotation(fig,'textbox',[0.0 0.0 1.0 0.2628],...
    'String',notes,...
    'Interpreter','latex',...
    'FitBoxToText','off',...
    'EdgeColor','none');

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the following exercises:
%   A) Parabolic. with friction.
%
% Help: Function prototype
% [y, t] = runge_kutta(f, y0, [t0 tfinal], timestep)
% where f is the function being integrated as dy/dt = f(t, y).
% y0 are the initial conditions
% t0: initial time
% tfinal: final time.
% h: time step for the calculations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function parabolic_with_friction()
global g k m
close all;
g = 9.81; %m/s^ 2. gravity at equator
% datos de ajustados a 1000 m de https://en.wikipedia.org/wiki/Ballistic_table#/media/File:Ballistic_table_for_7.62x51_mm_NATO_(mil_and_moa).png
% 7.62x51 mm NATO ammunition
theta0 = 0; %initial angle
k = 0.0000005; % N*s/m friction coefficient
m=9.7/1000; %mass of the bullet kg
v=860; % m/s initial speed

theta0 = 0; %initial angle
k = 0.0000005; % N*s/m friction coefficient
m=0.47; %mass of the bullet kg
v=609; % m/s initial speed
 % 914

t0 = 0;
tfinal = 2.2;
x0 = 0;
y0 = 0;
vx0 = v*cos(theta0);
vy0 = v*sin(theta0);

[t, x] = runge_kutta(@parabolic_friction_square, [x0 y0 vx0 vy0]', [t0 tfinal], 0.01);
x = x(:, 1:length(t));
plot(t, x(1,:), 'r'), hold
plot(t, x(2,:), 'g')
plot(t, x(3,:), 'b')
plot(t, x(4,:), 'c')
legend('Position X (m) RK4', 'Position Y (m) RK4', 'Speed X (m/s) RK4', 'Speed Y (m/s) RK4')

figure,
plot(x(1,:), x(2,:), 'r')
xlabel('X Position (m)')
ylabel('Y Position (m)')

x(2,end)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to solve a constant acceleration on x axis on a ramp.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = parabolic_friction_square(t, x)
global k g m
% We must return the solution of
% [dx1/dt; dx2/dt; dx3/dt dx4/dt]
vx = x(3);
vy = x(4);

%forces in both axes
sq = sqrt(vx^2+vy^2);
xd(1) = vx;
xd(2) = vy;
xd(3) = -(k/m)*vx*sq;
xd(4) = -(k/m)*vy*sq - g;
xd = xd(:);

repo: ARTE


text:
clear all; close all; clear classes; clc;

%% Set flags.
inspect_only = true;

%% Solve the system.
gray = [0.5 0.5 0.5];  % [r g b]
solveropts.method = 'aws';
[E, H, obj_array, src_array, J] = maxwell_run(...
	'OSC', 1e-9, 1550, ...
	'DOM', {'Palik/SiO2', 'none'}, [0, 700; 0, 600; -200, 1700], 20, [BC.e, BC.m, BC.p], [0 200; 0 200; 200 200], ...
	'OBJ', ...
		{'Palik/SiO2', 'none'}, Box([0, 50; 0, 50; -200, 1700], [2, 2, 20]), ...
		{'CRC/Ag', gray}, Box([25, 700; 0, 25; -200, 1700], 20), ...
	'SRC', ModalSrc(Axis.z, 200, 2.0), ...
	solveropts, inspect_only);

%% Visualize the solution.
if ~inspect_only
	figure;
	clear opts;
	opts.withabs = true;
	opts.withgrid = true;
	visall(E{Axis.x}, obj_array, src_array, opts);
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd