text:
function varargout = ihashsum(varargin)
% VL_IHASHSUM  Accumulate integer labels into a hash table
%   [H,ID,NEXT] = VL_IHASHSUM(H,ID,NEXT,K,X) counts the number of
%   occurences of the columns of X, accumulating these to the hash
%   table represented by the tripled H,ID,NEXT.
%
%   X is a D x N array of class UINT8 each row of which defines an D
%   dimensional label. Labels cannot be all zeros.
%
%   H and NEXT are 1 x C arrays of class UINT32 and ID is a D x C
%   array of class UINT8. H is a vector of counts, ID stores, for each
%   element of H, the corresponding label, and NEXT is a vector of
%   indexes.
%
%   Once constructed, the hash table can be searched by means of the
%   VL_IHASHFIND() function.
%
%   The hash table uses double hashing [1] with an initial size equal
%   to K (so that C >= K). Given a label X, this is first hashed by
%   using the FNV algorithm [2] to one of K bucket. If this bucket is
%   free, it is assigned to label X and the count is incremented.  If
%   the bucket is already assigned to the same label X, the count is
%   incremented. If the bucket is already assigned to a different
%   label, a second hash is used to scan (probe) the table for a free
%   bucket.
%
%   If no free/matching bucket is found (because the hash table is
%   full) an overflow area containing extra buckets is used. This is
%   visited by reading off indexe from the NEXT vector, until a
%   matching bucket is found or the overflow area is enlarged.
%
%   Example::
%     The following example counts integer bi-dimensional label
%     occurences:
%
%       K = 5 ;
%       h = zeros(1,K,'uint32') ;
%       id = zeros(2,K,'uint8');
%       next = zeros(1,K,'uint32') ;
%       X = uint8([1 1 ; 1 2 ; 2 1 ; 1 1]') ;
%       [h,id,next] = vl_ihashsum(h,id,next,K,X) ;
%
%     resulting in
%
%       h = [1 0 1 2 0]
%       id = [1    0    2    1    0
%             2    0    1    1    0]
%       next = [0 0 0 0 0]
%
%     For example, [1;2] has a count of 1 and [1;1] has a count of
%     2. NEXT is zero because there have been no collisions.
%
%   REFERENCES::
%   [1] http://en.wikipedia.org/wiki/Double_hashing
%   [2] http://www.isthe.com/chongo/tech/comp/fnv
%
%   See also: VL_IHASHFIND().
[varargout{1:nargout}] = vl_ihashsum(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function value = r8_huge ( )

%*****************************************************************************80
%
%% R8_HUGE returns a "huge" real number.
%
%  Discussion:
%
%    The value returned by this function is NOT required to be the
%    maximum representable R8.  This value varies from machine to machine,
%    from compiler to compiler, and may cause problems when being printed.
%    We simply want a "very large" but non-infinite number.
%
%    MATLAB provides a built-in symbolic constant "inf" that can be used
%    if a huge number is really what you want!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a huge number.
%
  value = 1.0E+30;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function rand = getNoise(app, mu, noise, num_points)
    rand = zeros(1, num_points);
    for i = 1:num_points
        rand(i) = normrnd(mu, noise);
    end
end

author: UMich-BipedLab


text:
function dist = line_exp_point_dist_3d ( p1, p2, p )

%*****************************************************************************80
%
%% LINE_EXP_POINT_DIST_3D: distance ( explicit line, point ) in 3D.
%
%  Discussion:
%
%    The explicit form of a line in 3D is:
%
%      ( P1, P2 ) = ( (X1,Y1,Z1), (X2,Y2,Z2) ).
%
%    Thanks to Francois Struempfer for pointing out a parenthesis mistake in the
%    computation of the Euclidean norm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3,1), P2(3,1), two points on the line.
%
%    Input, real P(3,1), the point whose distance from the line is
%    to be measured.
%
%    Output, real DIST, the distance from the point to the line.
%
  dim_num = 3;

  bot = sum ( ( p2(1:dim_num,1) - p1(1:dim_num,1) ).^2 );

  if ( bot == 0.0 )

    pn(1:dim_num,1) = p1(1:dim_num,1);
%
%  (P-P1) dot (P2-P1) = Norm(P-P1) * Norm(P2-P1) * Cos(Theta).
%
%  (P-P1) dot (P2-P1) / Norm(P-P1)**2 = normalized coordinate T
%  of the projection of (P-P1) onto (P2-P1).
%
  else

    t = sum ( ( p(1:dim_num,1) - p1(1:dim_num,1) )' ...
            * ( p2(1:dim_num,1) - p1(1:dim_num,1) ) ) / bot;

    pn(1:dim_num,1) = p1(1:dim_num,1) + t * ( p2(1:dim_num,1) - p1(1:dim_num,1) );

  end
%
%  Now compute the distance between the projection point and P.
%
  dist = sqrt ( sum ( ( p(1:dim_num,1) - pn(1:dim_num,1) ).^2 ) );

  return
end


text:
function out_m = updownsample( in_m,out_x_sz,out_y_sz,is_fourier_flag,is_real_flag )
%
% updownsample - up-sample or down-sample an input series using fourier domain
%                input series needs to be continuous of a high degree
%
% format:   out_m = updownsample( in_m,out_x_sz,out_y_sz,is_fourier_flag,is_real_flag )
%
% input:    in_m                - input matrix for up/down sampling. can be in
%                                 space domain OR in fourier domain, in such
%                                 case, needs to be in matlab format !!!
%                                 (matlab-format = the save as given from fft/fft2)
%           out_x_sz,out_y_sz   - desired number of pixels in the output image
%           is_fourier_flag     - 1: the input is given in the fourier domain
%                                 0: the input is given in the space domain
%                                    (we need to use fft2 to convert to fourier domain)
%           is_real_flag        - 0: the input is a complex matrix -> don't use
%                                    abs() at the output, perform complex
%                                    up/down sampling
%                                 1: the input is real BUT has negative values ->
%                                    use real() at the output
%                                 2: the input is real and positive -> using
%                                    abs() at the output
%
% output:   out_m               - up/down sampled image
%
% NOTE: it is important to specify if the image is REAL or COMPLEX, since
%       if the image is REAL -> we have to use ABS() on the inverse fourier
%       transform (because of roundoff errors of the transform).
%
% NOTE: since a desired amount of pixels is needed at the output, there is
%       no attempt to use matrices which are in size of power of 2. this
%       optimization can not be used in this case
%
% NOTE: input series needs to be CONTINUOUS OF A HIGH DEGREE, since the
%       upsampling is done in the frequency domain, which samples the output
%       grid with SINE-like (harmonic) functions
%
%
% Example:  type "updownsample" for an example. it will load the matlab's
%           child image to a temporary variable, as upsample it
%
%               out = updownsample( A,300,300,0,1 );
%               figure;
%               colormap gray;
%               subplot( 1,2,1 );
%               imagesc( A );
%               subplot( 1,2,2 );
%               imagesc( out );
%

%
% Theory:   the upsampling is done by zero-padding in the input domain BETWEEN the samples,
%           then at the fourier domain, taking the single spectrum (out of the repetition of spectrums)
%           i.e. low pass with Fcutoff=PI/upsample_factor, zeroing the rest of the spectrum
%           and then doing ifft to the distribution.
%           since we have a zero padding operation in time, we need to multiply by the fourier gain.
%
%              +-----------+     +-------+     +---------+     +--------+     +--------+
%   y[n,m] --> | up-sample | --> |  FFT  | --> |   LPF   | --> | * Gain | --> |  IFFT  | --> interpolated
%              | factor M  |     +-------+     | Fc=PI/M |     +--------+     +--------+
%              +-----------+                   +---------+
%
%           this operation is the same as the following one (which has less operations):
%
%              +-------+     +--------+     +--------------+     +--------+
%   y[n,m] --> |  FFT  | --> | * Gain | --> | Zero Padding | --> |  IFFT  | --> interpolated
%              +-------+     +--------+     +--------------+     +--------+
%
%           NOTE THAT, the zero-padding must be such that the D.C. ferquency remains the D.C. frequency
%           and that the zero padding is applied to both positive and negative frequencies.
%           The zero padding actually condences the frequency -> which yields a longer series in the
%           image domain, but without any additional information, thus the operation must be an interpolation.

% Copyright (c) 2003, Ohad Gal
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% ==============================================
%   check input parameters
% ==============================================
if (nargin==0)
    hFig    = figure( 'visible','off' );
    hImage  = image;
    Child   = get( hImage,'cdata' );
    close( hFig );
	out     = updownsample( Child,300,300,0,1 );
	figure;
	colormap gray;
	subplot( 1,2,1 );
	imagesc( Child );
	subplot( 1,2,2 );
	imagesc( out );
    return
elseif (nargin<5)
    error( 'UpDownSample - insufficient input parameters' );
end

% ==============================================
% get input image size, and calculate the gain
% ==============================================
[in_y_sz,in_x_sz] = size( in_m );
gain_x = out_x_sz/in_x_sz;
gain_y = out_y_sz/in_y_sz;

% ==============================================
% check if up/down sampling is needed at all
% ==============================================
if (gain_x == 1) & (gain_y==1)

    % same gain -> do not change sampling rate
    % ==========================================
    if is_fourier_flag
        switch is_real_flag
            case 0, out_m = ifft2( in_m );
            case 1, out_m = real( ifft2( in_m ) );
            case 2, out_m = abs( ifft2( in_m ) );
        end
    else
        out_m = in_m;
    end

else

    % upsample or downsample as needed
    % ==================================

    % convert to fourier domain, if input is given in the space domain
    if ~is_fourier_flag
        in_m = fft2(in_m);
    end

    % build grid vectors for the up/down sampling
    % ============================================
    % if the input is even & output is odd-> use floor for all
    % if the output is even & input is odd -> use ceil for all
    % other cases - don't care
    % for downsampling -> the opposite
    if (~mod( in_x_sz,2 ) & (out_x_sz>in_x_sz)) | (mod( in_x_sz,2 ) & (out_x_sz<in_x_sz))
        x_output_space  = max(floor((out_x_sz-in_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
        x_input_space   = max(floor((in_x_sz-out_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
    else
        x_output_space  = max(ceil((out_x_sz-in_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
        x_input_space   = max(ceil((in_x_sz-out_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
    end
    if (~mod( in_y_sz,2 ) & (out_y_sz>in_y_sz)) | (mod( in_y_sz,2 ) & (out_y_sz<in_y_sz))
       y_output_space  = max(floor((out_y_sz-in_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
       y_input_space   = max(floor((in_y_sz-out_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
   else
       y_output_space  = max(ceil((out_y_sz-in_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
       y_input_space   = max(ceil((in_y_sz-out_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
   end

    % perform the up/down sampling
    padded_out_m    = zeros( out_y_sz,out_x_sz );
    in_m            = fftshift(in_m);
    padded_out_m( y_output_space,x_output_space ) = in_m(y_input_space,x_input_space);
    out_m           = (gain_x*gain_y)*ifft2(ifftshift(padded_out_m));

    % check the output format real or complex
    switch is_real_flag
        case 0, % do nothing
        case 1, out_m   = real( out_m );
        case 2, out_m   = abs( out_m );
    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ point, seed ] = p08_sample ( m, n, seed )

%*****************************************************************************80
%
%% P08_SAMPLE samples points from the region in problem 08.
%
%  Discussion:
%
%    A rejection method is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real POINT(M,N), the coordinates
%    of the points.
%
  batch = 1000;

  [ lo, hi ] = p08_box ( m );

  have = 0;
%
%  We are going to generate batches of sample points.
%
  sample_num = min ( batch, n );

  reject = 0;

  while ( 1 )
%
%  Generate a batch of points in the bounding box.
%
    sample(1:m,1:sample_num) = rand ( m, sample_num );
%
%  Remap the points to the box.
%
    sample(1,1:sample_num) = lo(1) + sample(1,1:sample_num) * ( hi(1) - lo(1) );
    sample(2,1:sample_num) = lo(2) + sample(2,1:sample_num) * ( hi(2) - lo(2) );

    inside(1:sample_num) = p08_inside ( m, sample_num, sample );
%
%  Accept those points which are inside the region.
%
    for j = 1 : sample_num

      if ( inside(j) )

        have = have + 1;
        point(1:m,have) = sample(1:m,j);

        if ( have == n )
          return
        end

      else

        reject = reject + 1;

      end

    end

    if ( 10 * n < reject )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P08_SAMPLE - Fatal error!\n' );
      fprintf ( 1, '  Too many points rejected!\n' );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Total number of accepted points = %d\n', have );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Number generated on this sweep =  %d\n', sample_num );
      fprintf ( 1, '  Number accepted =                 %d\n', ...
        sample_num - reject );
      fprintf ( 1, '  Number rejected =                 %d\n', reject );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Something appears to be wrong!\n' );
      error ( 'P08_SAMPLE - Fatal error!' );
    end

  end

  return
end


text:
urls = { ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/GoodbyeOxford/camera/Picture004.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/GoodbyeOxford/camera/Picture200.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/MMOxford/camera/Picture400.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/MMOxford/camera/Picture509.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/StoLed/camera/Picture060.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/OxfordLateSpring/camera/Picture073.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/ParisAgain/camera/168_6814.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Barcelona/camera/167_6738.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Barcelona/camera/167_6730.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Lyon/camera/165_6532.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Lyon/camera/165_6535.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/OxfordSpring/camera/164_6422.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/MoreOxfordWinter/camera/162_6244.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/MoreOxfordWinter/camera/164_6409.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/MoreOxfordWinter/camera/162_6245.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/OxfordWinter/camera/IMG_5622.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/OxfordWinter/camera/IMG_5630.JPG', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Home/camera/img_5400.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Home/camera/img_5358.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Paris/camera/IMG_5134.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Paris/camera/IMG_5072.jpg', ...
'http://www.cs.berkeley.edu/~efros/photos/Albums/Places/camera/138_3807.JPG'};
alyosha_dir = 'D:\Projects\Databases\context_images\alyosha_images';
for i = 1:length(urls)
    im = wwwimread(urls{i});
    bn = urls{i};
    while(bn(end-3)~='.')
        [rem, bn] = strtok(bn);
    end
    disp(bn)
    imwrite(im, [alyosha_dir '/' bn]);
end

author: zouchuhang

repo: LayoutNet


text:
function [nx, ny, sJ] = Normals2D()

% function [nx, ny, sJ] = Normals2D()
% Purpose : Compute outward pointing normals at elements faces and surface Jacobians

Globals2D;
xr = Dr*x; yr = Dr*y; xs = Ds*x; ys = Ds*y; J = xr.*ys-xs.*yr;

% interpolate geometric factors to face nodes
fxr = xr(Fmask, :); fxs = xs(Fmask, :); fyr = yr(Fmask, :); fys = ys(Fmask, :);

% build normals
nx = zeros(3*Nfp, K); ny = zeros(3*Nfp, K);
fid1 = (1:Nfp)'; fid2 = (Nfp+1:2*Nfp)'; fid3 = (2*Nfp+1:3*Nfp)';

% face 1
nx(fid1, :) =  fyr(fid1, :);
ny(fid1, :) = -fxr(fid1, :);

% face 2
nx(fid2, :) =  fys(fid2, :)-fyr(fid2, :);
ny(fid2, :) = -fxs(fid2, :)+fxr(fid2, :);

% face 3
nx(fid3, :) = -fys(fid3, :);
ny(fid3, :) =  fxs(fid3, :);

% normalise
sJ = sqrt(nx.*nx+ny.*ny); nx = nx./sJ; ny = ny./sJ;
return;

repo: Aero-matlab


text:
function configPolytopeSimplifier()
path = fileparts(mfilename('fullpath'));
addpath(fullfile(path, '..'));
config;

loadModular('CMatrix');

author: opencobra

repo: cobratoolbox


text:
% TODO delete this. looks like someone's play file. has error, undocumented
report_this_filefun(mfilename('fullpath'));

ttcat = newt2;
mati = maepi(1,3);
%[p,sdp] = mypval(3,mati);
% [p,sdp,c,sdc,dk,sdk,aa,bb]=mypval2(3, mati);

tmin1 = 0.05;

p0 = 1.44;
b0 = 0.96;
c0 = 1.1;
A0 = -1.62;

save_aspar2;
do = [ ' ! '  hodi '/aspar/myaspar' ];
eval(do)

load aspar3.out
re = aspar3;

bv = re(2,3);
p =  re(1,3)
c =  re(4,3);
A =  re(3,3);

%[bv magco stan av me mer me2,  pr] =  bvalca3(newt2,1,1);
%[me, bv, si, av] = bmemag(newt2) ;

%A = log10(bv/av)
la = 0;
m0 = maepi(1,6);

m = min(newt2.Magnitude);
dt = 0.5;

%t0 = ( max(newt2.Date) - mati)*365;
t0 = tlen;

la = [];ti = [];
for t = c:dt:t0
    la = [la (10^(A + bv*(m0-m)) * (t + c)^(-p))*dt ];
    ti = [ti mati+t/365];
end

la2 = [];ti2 = [];
for t = c0:dt:t0
    la2 = [la2 (10^(A + b0*(m0-m)) * (t + c0)^(-p0))*dt ];
    ti2 = [ti2 mati+t/365];
end
pla = 0; pla2 = 0;
M = max(a.Magnitude) - 5.;
for t = t0:dt:t0+365
    pla  = pla + (10^(A + bv*(M)) * (t + c)^(-p))  *dt;
    pla2 = pla2 + (10^(A + b0*(M)) * (t + c0)^(-p0))  *dt;
end

P = 1 - exp(-pla)
P0  = 1 - exp(-pla2);
figure_w_normalized_uicontrolunits(cum)

l = newt2.Date >= mati + c/365;
tmpn = newt2(l,:);
nu = (1:length(tmpn(:,3))+1); nu(length(tmpn(:,3))+1) = length(tmpn(:,3));
try delete(plc); catch ME, error_handler(ME,@do_nothing);end
try delete(plc2); catch ME, error_handler(ME,@do_nothing);end
try delete(plc3); catch ME, error_handler(ME,@do_nothing);end

hold on
plc = plot([tmpn(:,3) ; teb],nu,'k');

plc2 = plot(ti,cumsum(la),'r');
%plc3 = plot(ti2,cumsum(la2),'g');

author: CelsoReyes


text:
% Add to matlab Path
addpath(genpath('NodalDG'),'SBBGK','Quad','TVD','WENO')


text:
function [x,y,z,s,w,flag] = squadsolve(Q,c,A,b,C),
%SQUADSOLVE
%
% USAGE:   [x,y,z,s,w,flag] = squadsolve(Q,c,A,b,C)
%
% PARAMETERS:  Q -> (n,n) symetric matrix (definite positive)
%              c -> (n,1) vector
%              A -> (m,n) matrix
%              b -> (m,1) vector
%              C -> scalar
%
%            x -> primal variables
%            y -> lagrangian coeff of equality constraints
%            z -> dual variables of x
%            s -> primal auxiliary variable (only if C < Inf)
%            w -> dual variable of s
%            flag -> set to 0 => no problem, set to 1 => problem
%
% DESCRIPTION: Primal-dual method for quadratic programming
%
%            minimize c'*x + 0.5*x'*Q*x
%
%            subject to  A*x=b
%                        0<= x <= C
%            The method used here is a primal dual method with a predictor-corrector
%            approach and a logarithmic barrier. I used the heuristic from two
%            existing methods LOQO and HOPDM. The method is an iterative method. The
%            maximal number of iteration is stored in the variable 'max_iter'.
%
% ERRORS AND BUGS: there is no test about the conditionning of the matrix Q. If the iteration
%                  50 has been reached, then the optimization may not be finished and the output
%                  may be wrong.
%
% NOTES: 50 iterations have always been sufficient to solve all problems.
%        This code should be read with the tech. report:
%                "Regularized Symmetric Indefinite Systems in Interior Point
%                 Methods for Linear and Quadratic Optimization",
%                 A. Altman and J. Gondzio, Logilab Tech. Report 1998.6
%
% Andre Elisseeff, Dec. 1999
% aelissee@eric.univ-lyon2.fr

 disp('toto');
% init
  verbose = 1;
  n = size(Q,1);
  m = size(A,1);
  H = zeros(n+m,n+m);
  flag=1;

% Values of the original HOPDM of Gondzio and Altmann
  dinf = 10^(-14);
  smallz = 10^(-14);
  smallt = 2.3*10^(-16);
  opttol = 10^(-5)/C;
  itref = 1;
  mu = 1;
  maxiter = 50;

% init values of the primal and dual variables
  x=ones(n,1);
  z=ones(n,1);
  y=ones(m,1);
  if C < Inf,
    s=ones(n,1);
    w=ones(n,1);
  else
    s=[];w=[];
  end;
% Description of variables:
%
%    x,s     -> primal variables
%    z,w     -> dual variables
%    n       -> number of variables in the initial pb (size of x)
%    m       -> number of constraints in A
%
%    dinf     -> smallest value for all variables
%    smallz   -> smallest value of z
%    smallt   -> smallest value for t in the computation of the matrix theta
%    opttol   -> acceptable tolerance for optimality conditions
%    itref    -> iteration counter
%    maxiter  -> maximum number of iteration
%    DEBUG    -> 0 => no debug, 1 => debug

% Analyze the constraints...
  if verbose,
   disp(sprintf('Analyzing the equality constraints...\n'));
  end;
  [QQ,RR]=qr(A',0);
  [mm,nb] = size(QQ); %% number of eq constraints
  ind = 1:1:nb;
  for i=1:nb,
    if abs(RR(i,i)) < 10^(-10),
      disp(sprintf('Constraints %d removed because of dependence\n',i));
      ind(i)=0;
    end;
  end;
  indice = find(ind >0);

  if (isempty(indice))
    disp(sprintf('No equality constraints... \n'));
    A=[];
    m=0;
  else
    A = A(indice,:); %% new independent eq constraints
    b = b(indice);
    y = y(indice);
    m=length(indice);
  end;
  clear QQ;clear RR;
%%%% AE : scale the problem, has been tested and seems to work
%%%% better : more stable...
  if C < Inf,
    u = ones(n,1);
    Q = Q*C;
    b = b/C;
    c = c;
  end;
% init values before looping
  cont = 1;
  objQ = 0.5*x'*Q*x;
% init values of primal and dual objective functions
  pobjo = abs(c'*x+objQ) + 1;
  if C < inf,
    dobjo = abs(b'*y - u'*w - objQ);
  else
    dobjo = abs(b'*y - objQ);
  end;
%%%%%%%%%%%%%%%
%% MAIN LOOP
%%%%%%%%%%%%%%%
  while (cont)&(itref<=maxiter)
  % Compute the primal objective function
    objQ = 0.5*x'*Q*x;
    pobj = c'*x+objQ;
    if (C<Inf)
        if ~isempty(A),
            dobj = b'*y - u'*w - objQ;
        else
            dobj = - u'*w - objQ;
        end;
    else
        if ~isempty(A),
            dobj = b'*y - objQ;
        else
            dobj = - objQ;
        end;
    end;
    dlgap = pobj - dobj;
    dp = abs(pobj)/(abs(pobjo)+1);
    dd = abs(dobj)/(abs(dobjo)+1);
    dobjo = dobj;
    pobjo = pobj;
    if verbose,
        disp(sprintf('%d - pobj : %f - dobj : %f\n',itref,pobj,dobj));
    end;
  % Check if the solution are bounded
    if (dp > 10^6)
      disp(sprintf('Solution not bounded in the primal. Exit.\n'));
      return;
    end;
    if (dd > 10^6)
      disp(sprintf('Solution not bounded in the dual. Exit.\n'));
      return;
    end;

  % test if optimality
    oldgap = dlgap;
    dp = abs(dobj) + 1;
    if ((abs(dlgap)/dp) <= opttol)
      if verbose,
       disp(sprintf('Optimal solution found. Exit.\n'));
      end;
      cont = 0;
      break;
    end;

    dp = dp + abs(pobj);
    T = abs(dlgap)/dp;

  % put the variables away from zero (from HOPDM)
    if (itref <= 3)
      ax = 2*10^(-3);
      az = 10^(-3);
    elseif (T >= 0.8)
      ax = 2*10^(-4);
      az = 10^(-4);
    elseif (T >= 0.1)
      ax = 2*10^(-5);
      az = 10^(-5);
    elseif (T >=0.01)
      ax = 2*10^(-6);
      az = 10^(-6);
    elseif (T>=0.001)
      ax = 2*10^(-7);
      az = 10^(-7);
    elseif (T>=0.0001)
      ax = 2*10^(-8);
      az = 10^(-7);
    elseif (T>=0.00001)
      ax = 2*10^(-9);
      az = 10^(-9);
    else
      ax = T*10^(-5);
      az = ax;
    end;

  % consider only variables that can be changed
    x = x + ax;
    z = z + az;
    if C < Inf
      s = s + ax;
      w = w + az;
    end;
  % Compute the values of xi_b, xi_c and xi_u
    if ~isempty(A),
     xi_b = -A*x + b;
     xi_c = c - A'*y - z + Q*x;
    else
      xi_b = [];
      xi_c = c - z + Q*x;
     end;

    xi_z =   - x.*z;
    if C < Inf,
      xi_c=xi_c + w;
      xi_u = u - x - s;
      xi_w =  - s.*w;
    end;
  % Compute theta = (z/x + w/s)

  % for bounded variables
    if C<Inf,
      dp = x;
      if (max(abs(dp))<= smallz)
        disp(sprintf('Conditioning problem to invert theta. Abort.\n'));
        return;
      end;
      dpp= s;
      if (max(abs(dpp))<= smallz)
        disp(sprintf('Conditioning problem to invert theta. Abort.\n'));
        return;
      end;
      theta=z./dp + w./dpp;
    end;
  % for unbounded variables
    if C==Inf,
      dp = x;
      if (max(abs(dp))<= smallz)
        disp(sprintf('Conditioning problem to invert theta. Abort.\n'));
        return;
      end;
      theta = z./dp;
    end;

  % neglect small elements of theta array

    neglect = find(theta < smallt);
    if ~isempty(neglect)
      theta(neglect)=zeros(size(neglect));
    end;

  % and control large elements of theta

    neglect = find(theta >= 10^8);
    if ~isempty(neglect)
      theta(neglect)=(10^4)*sqrt(theta(neglect));
    end;
  % factorize H = [-Q-theta^(-1)   A^T]
  %               [ A               0 ]

    H = zeros(n+m,n+m);
    H(1:n,1:n) = -Q-diag(theta);
    H(n+1:n+m,1:n) = A;
    H(1:n,n+1:n+m) = A';

  % Compute the predictor step
    if C < Inf,
      f = xi_c-xi_z./x+(xi_w - xi_u.*w)./s ;
      h = xi_b;
    else
      f = xi_c - xi_z./x;
      h = xi_b;
    end;
    delta=H\[f;h];
    dx = delta(1:n);
    dy = delta(n+1:n+m);
    dz = (xi_z-z.*dx)./x;
    if C<Inf,
      ds = xi_u - dx;
      dw = (xi_w-w.*ds)./s;
    end;

 % determine the maximum step size alpha_p (primal) and
 % alpha_d (dual) to stay in feasible region
 % (x,s,z,w must be positive and greater than dinf)
      indz = find(dz<0);
      indx = find(dx<0);
      inds=[];mins=1;
      indw=[];minw=1;
      if C < Inf,
        inds = find(ds<0);
        indw = find(dw<0);
        if ~isempty(inds)
          mins = min(-(s(inds)-dinf)./ds(inds));
        else
         mins = 1;
        end;
        if ~isempty(indw)
          minw = min(-(w(indw)-dinf)./dw(indw));
        else
         minw = 1;
        end;
      end;
      if ~isempty(indx)
        minx = min(-(x(indx)-dinf)./dx(indx));
      else
        minx = 1;
      end;
      apk = min([minx,mins,1]);
      if ~isempty(indz),
        minz = min(-(z(indz)-dinf)./dz(indz));
      else
        minz = 1;
      end;
      adk = min([minw,minz,1]);

      ax = sum(x.*z);
      as = sum((x+apk*dx).*(z+adk*dz));
      az = sum(dx.^2+dz.^2);
      if C < Inf,
        ax = ax + sum(s.*w);
        as = as + sum((s+apk*ds).*(w+adk*dw));
        az = az + sum(ds.^2+dw.^2);
      end;
    % check if complementary gap is less than opttol
      if (as <= opttol)
        if verbose,
            disp(sprintf('Complementary gap is less than %f\n',opttol));
        end;
        cont = 0;
        x = x + apk*dx;
        y = y + adk*dy;
        z = z + adk*dz;
        if C < Inf,
         s = s + apk*ds;
         w = w + adk*dw;
        end;
        break;
      end;

    % Set the barrier parameter : LOQO's heuristic
      ap = min(apk,adk);
      mu = (ax/(2*n))*(0.95*(1/ap) -1)^2/(0.95*(1/ap)+10)^2;

    % Compute the new direction (algo. of Mehrotra) of order 1 (corrector step)
      xi_z =  - x.*z + mu*ones(size(x)) - dx.*dz;
      f = xi_c - xi_z./x;
      if C < Inf
        xi_w = - s.*w + mu*ones(size(s)) - ds.*dw;
        f=f+ xi_w./s - (w.*xi_u)./s;
      end;
      h = xi_b;

      delta=H\[f;h];
      dx = delta(1:n);
      dy = delta(n+1:n+m);
      dz = (xi_z-z.*dx)./x;
      if C<Inf,
        ds = xi_u - dx;
        dw = (xi_w-w.*ds)./s;
      end;

    % determine the maximum step size alpha_p (primal) and
    % alpha_d (dual) to stay in feasible region
    % (x,s,z,w must be positiv)
      indz = find(dz<0);
      indx = find(dx<0);
      if C < Inf,
        inds = find(ds<0);
        indw = find(dw<0);
        if ~isempty(inds)
          mins = min(-s(inds)./ds(inds));
        else
      mins=1;
        end;
        if ~isempty(indw)
          minw = min(-w(indw)./dw(indw));
        else
      minw=1;
        end;
      end;
      if ~isempty(indx)
        minx = min(-x(indx)./dx(indx));
      else
        minx=1;
      end;
      alpha_p = min([minx,mins,1]);
      if ~isempty(indz)
        minz = min(-z(indz)./dz(indz));
      else
        minz=1;
      end;
      alpha_d = min([minw,minz,1]);

    % Compute step factors
      fp =0.9*min(alpha_p,alpha_d);
      fd =0.9*min(alpha_d,alpha_p);
      x = x +fp*dx;
      y = y +fd*dy;
      z= z + fd*dz;
      if C < Inf,
        w = w + fd*dw;
        s = s +fp*ds;
      end;
      itref=itref+1;
  end;
%%%%%%%%%%%%%%%%%%%%%
%%% End of main loop
%%%%%%%%%%%%%%%%%%%%%

    if (cont==0)
        if verbose,
          disp(sprintf('Optimal Solution found after %d iteration.\n',itref));
          disp(sprintf('Value of the objective : %f.\n',pobj));
        end;
      flag=0;
    end;
  % rescale x and y
    if C<Inf,
      x= C*x;
      %y = C*y;
    end;

end;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% tHis subroutine assigns creates a grid with
% spacing dx,dy (in degreees). The size will
% be selected interactiVELY. The bvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/95

report_this_filefun(mfilename('fullpath'));

global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dd = 1.00;
    dx = 1.00 ;
    ni = 500;

    def = {num2str(ZG.maepi.Date(1))};
    ni2 = inputdlg('Input Time of Mainshock ?','Input',1,def);
    l = ni2{:};
    mati = str2double(l);

    def = {'30'};
    ni2 = inputdlg('Time length considered in days?','Input',1,def);
    l = ni2{:};
    tlen = str2double(l);

    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Grid Input Parameter',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ ZG.wex+200 ZG.wey-200 550 300]);
    axis off

    %
    freq_field=uicontrol('Style','edit',...
        'Position',[.60 .50 .22 .10],...
        'Units','normalized','String',num2str(ni),...
        'Callback','ni=str2double(freq_field.String); freq_field.String=num2str(ni);');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.60 .40 .22 .10],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(freq_field2.String); freq_field2.String=num2str(dx);');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.60 .30 .22 .10],...
        'Units','normalized','String',num2str(dd),...
        'Callback','dd=str2double(freq_field3.String); freq_field3.String=num2str(dd);');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .15 .12 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'Callback','close,sel =''ca''; pcross',...
        'String','Go');

    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.30 0.65 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String',' Grid Parameter');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.42 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing along projection [km]');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.32 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in depth in km:');

    txt1 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.53 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m,...
        'FontWeight','bold',...
        'String','Number of Events (Ni):');

    set(gcf,'visible','on');
    watchoff

end   % if sel == in

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'

    figure_w_normalized_uicontrolunits(xsec_fig)
    hold on

    messtext=...
        ['To select a polygon for a grid.       '
        'Please use the LEFT mouse button of   '
        'or the cursor to the select the poly- '
        'gon. Use the RIGTH mouse button for   '
        'the final point.                      '
        'Mac Users: Use the keyboard "p" more  '
        'point to select, "l" last point.      '
        '                                      '];

    zmap_message_center.set_message('Select Polygon for a grid',messtext);

    ax = findobj('Tag','main_map_ax');
    [x,y, mouse_points_overlay] = select_polygon(ax);
    zmap_message_center.set_info('Message',' Thank you .... ')

    plos2 = plot(x,y,'b-','era','xor');        % plot outline
    sum3 = 0.;
    pause(0.3)

    %create a rectangular grid
    xvect=[min(x):dx:max(x)];
    yvect=[min(y):dd:max(y)];
    gx = xvect;gy = yvect;
    tmpgri=zeros((length(xvect)*length(yvect)),2);
    n=0;
    for i=1:length(xvect)
        for j=1:length(yvect)
            n=n+1;
            tmpgri(n,:)=[xvect(i) yvect(j)];
        end
    end
    %extract all gridpoints in chosen polygon
    XI=tmpgri(:,1);
    YI=tmpgri(:,2);

    ll = polygon_filter(x,y, XI, YI, 'inside');
    %grid points in polygon
    newgri=tmpgri(ll,:);

    % Plot all grid points
    plot(newgri(:,1),newgri(:,2),'+k')

    if length(xvect) < 2  ||  length(yvect) < 2
        errordlg('Selection too small! (not a matrix)');
        return
    end

    itotal = length(newgri(:,1));

    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(newa.Date)  ;
    n = newa.Count;
    teb = newa.Date(n) ;
    tdiff = round(days(teb-t0b)/ZG.bin_days);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    bvg = [];
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','p-value grid - percent done');;
    drawnow
    %
    % loop
    dm = 0.1;
    dt = 0.01;
    M = max(ZG.newt2.Magnitude) - 5.;

    %
    for i= 1:length(newgri(:,1))
        x = newgri(i,1);y = newgri(i,2);
        allcount = allcount + 1.;
        i2 = i2+1;

        % calculate distance from center point and sort wrt distance
        l = sqrt(((xsecx' - x)).^2 + ((xsecy + y)).^2) ;
        [s,is] = sort(l);
        b = newa(is(:,1),:) ;       % re-orders matrix to agree row-wise

        % take first ni points
        b = b(1:ni,:);      % new data per grid point (b) is sorted in distance
        l2 = sort(l); di = l2(ni);

        [st,ist] = sort(b);   % re-sort wrt time for cumulative count
        b = b(ist(:,3),:);

        % call the p-value function
        % first find out what magco is
        l = b.Date > mati + days(3);
        %[bv magco stan av me mer me2,  pr] =  bvalca3(b(l,:),1,1);
        %l = b.Magnitude > magco+0.1;
        ZG.newt2 = b;
        tmin1 =0.05;

        calcp
        load aspar3.out
        re = aspar3;

        bvg = [bvg ; re(1,2) re(1,4) x y re(2,2) re(2,4) bv 0 di];
        waitbar(allcount/itotal)
    end  % for  newgri

    % save data
    %
    %  set(txt1,'String', 'Saving data...')
    drawnow
    gx = xvect;gy = yvect;

    %   catSave3 =...
    % [ 'zmap_message_center.set_info(''Save Grid'',''  '');think;',...
    %  '[file1,path1] = uiputfile(fullfile(ZmapGlobal.Data.data_dir, ''*.mat''), ''Grid Datafile Name?'') ;',...
    %  ' sapa2 = [''save '' path1 file1 '' ll tmpgri bvg xvect yvect gx gy dx dd ZG.bin_days ni newa maex maey maix maiy ''];',...
    %  ' if length(file1) > 1, eval(sapa2),end , done']; eval(catSave3)

    close(wai)
    watchoff

    % reshape a few matrices
    %
    normlap2=ones(length(tmpgri(:,1)),1)*nan;
    normlap2(ll)= bvg(:,1);
    re3=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,2);
    dp =reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,6);
    db =reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,9);
    r =reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,5);
    bv=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,7);
    bv2=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,8);
    pmap =reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,1);
    Pv=reshape(normlap2,length(yvect),length(xvect));

    old = re3;

    % View the b-value map
    view_pv2

end   %  if sel = ca

% Load exist b-grid
if sel == 'lo'
    load_existing_bgrid_version_A
end


text:
function [ Pos ] = getEEFCartesianOrientation( t )
%% This function is used to get the endeffector cartizian orientation of the KUKA iiwa 7 R 800.

%% Syntax:
% [ Ori ] = getEEFCartisianOrientation( t )

%% About:
% This function is used to get the endeffector cartizian orientation of the KUKA iiwa 7 R 800.
% The orientation is of the media flange of the robot relative to the base
% frame of the robot.

%% Arreguments:
% t: is the TCP/IP connection
% Ori: is 1x3 cell array,  alpha,beta and gama angles of the endeffector,
% the angles are measured in radians.

% Copyright, Mohammad SAFEEA, 3rd of May 2017

theCommand='Eef_pos';
fprintf(t, theCommand);
message=fgets(t);
%disp(message)
[P1,N]=getDoubleFromString(message);
ttt=0;
for i=4:6
    ttt=ttt+1;
    Pos{ttt}=P1{i};
end

end

function [jPos,j]=getDoubleFromString(message)
n=max(max(size(message)));
j=0;
numString=[];
for i=1:n
    if message(i)=='_'
        j=j+1;
        jPos{j}=str2num(numString);
        numString=[];
    else
        numString=[numString,message(i)];
    end
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function edge_length = triangle_edge_length_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_EDGE_LENGTH_2D returns edge lengths of a triangle in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real EDGE_LENGTH(3,1), the length of the edges.
%
  for j1 = 1 : 3
    j2 = i4_wrap ( j1 + 1, 1, 3 );
    edge_length(j1,1) = sqrt ( ( t(1,j2) - t(1,j1) )^2 ...
                             + ( t(2,j2) - t(2,j1) )^2 );
  end

  return
end


text:
function [THD_U,THD_I,K_U,K_I,MaxHarm_I,MaxIHarm_I,MaxHFHarm_I,MaxHarm_U,MaxIHarm_U,MaxHFHarm_U,Order_Harm,Order_IHarm,Order_HFHarm] = IEC_Harmonics_final(f,Vn,Inom,U1,U2,U3,I1,I2,I3,U1_time,U2_time,U3_time,I1_time,I2_time,I3_time,figures)
% This function will calculate the Harmonics spectrum values according to
% the IEC 61000-4-7 Norm which spesify the need for 10 min Voltage and
% current measured signals of the system under study and calulate the
% Harmonics mean values every 10 Periods of the measured signal (In this case 50 Hz signal)
% This Fucntion inputs are:
%           1- The nominal values of the measured signal: Frequency,
%           Voltage, and the current. For example f = 50%Hz
%                                                 Vn = 400%V (Phase to Phase Voltage of the grid)
%                                                 Inom = 475%A (The nominal
%                                                 current of the system under measurement)
%           2- The measured Voltages and current signals with there time
%           measurement signals (U1,U2,U3,I1,I2,I3,U1_time,U2_time,U3_time,I1_time,I2_time,I3_time)
%           3- figures is an input of 1 or 0 in order to activate showing the figures of the different  harmonics spectrums figures calculated in this function.
%PS: This function assume that you have the same sampling rate of the
%different measured signals and that all measured signals has the same
%length.
Fn = f;
Un = Vn;
In = Inom;
Delta_Time_Orginal = I1_time(2)-I1_time(1);
F_Sample = abs(1/Delta_Time_Orginal);
%------------------------------------------------------------------------------------------------------
%---------------------- The Harmonics Spectrum of the complet input signal ----------------------------
%------------------------------------------------------------------------------------------------------
Tp = length(I3_time);
T_S = abs(Tp*Delta_Time_Orginal);
[nyquist_number_I1,Bezugsfrequenz_I1,Signal_amp_I1] = Fast_Harmonics (I1,F_Sample,T_S);
[nyquist_number_I2,Bezugsfrequenz_I2,Signal_amp_I2] = Fast_Harmonics (I2,F_Sample,T_S);
[nyquist_number_I3,Bezugsfrequenz_I3,Signal_amp_I3] = Fast_Harmonics (I3,F_Sample,T_S);
[nyquist_number_U1,Bezugsfrequenz_U1,Signal_amp_U1] = Fast_Harmonics (U1,F_Sample,T_S);
[nyquist_number_U2,Bezugsfrequenz_U2,Signal_amp_U2] = Fast_Harmonics (U2,F_Sample,T_S);
[nyquist_number_U3,Bezugsfrequenz_U3,Signal_amp_U3] = Fast_Harmonics (U3,F_Sample,T_S);
if figures
    figure;
    subplot(2,1,1);
    plot([0:nyquist_number_I1]*Bezugsfrequenz_I1,Signal_amp_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
    plot([0:nyquist_number_I2]*Bezugsfrequenz_I2,Signal_amp_I2,'-.r');hold on
    plot([0:nyquist_number_I3]*Bezugsfrequenz_I3,Signal_amp_I3,'--k');hold off
    subplot(2,1,2);
    plot([0:nyquist_number_U1]*Bezugsfrequenz_U1,Signal_amp_U1);xlabel('frequency [Hz]'),ylabel('amplitude Voltages [V]');hold on
    plot([0:nyquist_number_U2]*Bezugsfrequenz_U2,Signal_amp_U2,'-.r');hold on
    plot([0:nyquist_number_U3]*Bezugsfrequenz_U3,Signal_amp_U3,'.k','MarkerSize',0.1);hold off
end
%**********************************************************************
%---------------------- Creating a grouping matrix --------------------
%**********************************************************************
Sampling_Time = Delta_Time_Orginal;
T_10_period = 10/Fn;
F_Sample_Orginal = 1/Sampling_Time;
if abs(F_Sample_Orginal) < 20000
    F_Sample_Complete_Data = 20000;
    [p, q] = rat(abs(Delta_Time_Orginal/(1/F_Sample_Complete_Data)));
    I1_Resampled = resample(I1,p,q);
    I2_Resampled = resample(I2,p,q);
    I3_Resampled = resample(I3,p,q);
    U1_Resampled = resample(U1,p,q);
    U2_Resampled = resample(U2,p,q);
    U3_Resampled = resample(U3,p,q);
    t_P = I1_time(1):((q/p)/F_Sample_Orginal):(I1_time(end)+((q/p)/F_Sample_Orginal));
    tt = length(t_P);
    if length(t_P) == length(I1_Resampled)
        I1_Resampled_time = t_P;
        I2_Resampled_time = t_P;
        I3_Resampled_time = t_P;
        U1_Resampled_time  = t_P;
        U2_Resampled_time  = t_P;
        U3_Resampled_time  = t_P;
    else if length(t_P) < length(I1_Resampled)
            for kk = 1:abs(length(I1_Resampled)-length(t_P))
                t_P(tt+kk) = ((q/p)/F_Sample_Orginal)+t_P(tt+kk-1);
            end
            I1_Resampled_time = t_P;
            I2_Resampled_time = t_P;
            I3_Resampled_time = t_P;
            U1_Resampled_time  = t_P;
            U2_Resampled_time  = t_P;
            U3_Resampled_time  = t_P;
        elseif length(t_P) > length(I1_Resampled)
            t_P = t_P(1:length(I1_Resample));
            I1_Resampled_time = t_P;
            I2_Resampled_time = t_P;
            I3_Resampled_time = t_P;
            U1_Resampled_time  = t_P;
            U2_Resampled_time  = t_P;
            U3_Resampled_time  = t_P;
        end
    end
    Delta_Time_Complete_Data = 1/F_Sample_Complete_Data;
    Sampling_Time_New = Delta_Time_Complete_Data;
    Sampling_frq_New = (1/Sampling_Time_New)/1000;
    Sample_New = T_10_period*Sampling_frq_New*1000;
    ii= 1;
    Points = 2^ii;
    while Points < Sample_New
        Points = 2^ii;
        ii = ii+1;
    end
else if abs(F_Sample_Orginal) > 20000
        F_Sample_Complete_Data = 20000;
        [p, q] = rat(abs(Delta_Time_Orginal/(1/F_Sample_Complete_Data)));
        I1_Resampled = resample(I1,p,q);
        I2_Resampled = resample(I2,p,q);
        I3_Resampled = resample(I3,p,q);
        U1_Resampled = resample(U1,p,q);
        U2_Resampled = resample(U2,p,q);
        U3_Resampled = resample(U3,p,q);
        t_Q = I1_time(1):((q/p)/F_Sample):I1_time(end);
        I1_Resampled_time = t_Q;
        I2_Resampled_time = t_Q;
        I3_Resampled_time = t_Q;
        U1_Resampled_time  = t_Q;
        U2_Resampled_time  = t_Q;
        U3_Resampled_time  = t_Q;
        Delta_Time_Complete_Data = 1/F_Sample_Complete_Data;
        %Sampling_Time_New = Delta_Time_Complete_Data;
        Sampling_frq_New = (1/Sampling_Time_Complete_Data)/1000;
        Sample_New = T_10_period*Sampling_frq_New*1000;
        ii= 1;
        Points = 2^ii;
        while Points < Sample_New
            Points = 2^ii;
            ii = ii+1;
        end
    else
        I1_Resampled = I1;
        I2_Resampled = I2;
        I3_Resampled = I3;
        U1_Resampled = U1;
        U2_Resampled = U2;
        U3_Resampled = U3;
        I1_Resampled_time = I1_time;
        I2_Resampled_time = I2_time;
        I3_Resampled_time = I3_time;
        U1_Resampled_time  = U1_time;
        U2_Resampled_time  = U2_time;
        U3_Resampled_time  = U3_time;
        F_Sample_Complete_Data = 20000;
        Delta_Time_Complete_Data = I1_Resampled_time(2)-I1_Resampled_time(1);
        Sampling_frq = 1/Delta_Time_Complete_Data;
        T_10_period = 10/Fn;
        Sample = T_10_period*Sampling_frq;
        ii= 1;
        Points = 2^ii;
        while Points < Sample
            Points = 2^ii;
            ii = ii+1;
        end
    end
end
%----------------------------------------------------------------------
%--------------- Zeros Results Array Definition -----------------------
%----------------------------------------------------------------------
%--------------- Try --------------------
LOWPASS_CUTOFF = 100;% Filter cutoff frequency is 100 Hz
LOWPASS_ORDER = 4;% Filter order is 4
[b_bw, a_bw] = butter(LOWPASS_ORDER, LOWPASS_CUTOFF / (F_Sample_Complete_Data / 2), 'low');
Filter = filter(b_bw, a_bw, U1_Resampled);
%--------------- End --------------------
null_L1 = find(diff(Filter < 0));
null_L1_x = Zeros_finding(Filter,I1_Resampled_time');
Tp = length(U1_Resampled_time);
len = length(null_L1_x);
%len_1_period = (1/Fn)/Delta_Time_Complete_Data;
%len_10_period = (10*(1/Fn))/Delta_Time_Complete_Data;
Results_U_H1 = zeros(50,length(0:20:len-30));
Results_U_H2 = zeros(50,length(0:20:len-30));
Results_U_H3 = zeros(50,length(0:20:len-30));
Results_U_IH1 = zeros(39,length(0:20:len-30));
Results_U_IH2 = zeros(39,length(0:20:len-30));
Results_U_IH3 = zeros(39,length(0:20:len-30));
Results_U_HF1 = zeros(35,length(0:20:len-30));
Results_U_HF2 = zeros(35,length(0:20:len-30));
Results_U_HF3 = zeros(35,length(0:20:len-30));
Results_P = zeros(length(0:20:len-30));
Results_Q = zeros(length(0:20:len-30));
%----------------------------------------------------------------------
%--------------------------- End --------------------------------------
%----------------------------------------------------------------------
%**********************************************************************
%----------------------------------END---------------------------------
%**********************************************************************
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%**********************************************************************
%--------------------------- Gruoping ---------------------------------
%**********************************************************************
j = 0;
jj = 1;
while j < len-30 % t < t2+0.1
    data_I1 = I1_Resampled(null_L1(j+3):null_L1(j+23)+1);
    data_I2 = I2_Resampled(null_L1(j+3):null_L1(j+23)+1);
    data_I3 = I3_Resampled(null_L1(j+3):null_L1(j+23)+1);
    data_U1 = U1_Resampled(null_L1(j+3):null_L1(j+23)+1);
    data_U2 = U2_Resampled(null_L1(j+3):null_L1(j+23)+1);
    data_U3 = U3_Resampled(null_L1(j+3):null_L1(j+23)+1);
    data_I1_time = I1_Resampled_time(null_L1(j+3):null_L1(j+23)+1);
    %------------------- Cut End --------------------------------------
    %------------------------------------------------------------------
    %--------- Data Resampling with the new Delta_10_Periods ----------
    %------------------------------------------------------------------
    delta = (null_L1_x(j+23)-null_L1_x(j+3))/Points;
    %Rampe = null_L1_x(j+3):delta:(((Points-1)*delta)+null_L1_x(j+3));
    iii=2;
    Rampe = zeros();
    Rampe(1) = null_L1_x(j+3);
    while iii <= Points
        Rampe(iii) = delta+Rampe(iii-1);
        iii = iii+1;
    end
    F_Sample_10_Periods = 1/delta;
    if abs(F_Sample_10_Periods) < F_Sample_Complete_Data
        F_Sample_New = F_Sample_10_Periods;
        [p, q] = rat(abs(Delta_Time/(1/F_Sample_New)));
        I1_Resampled_10_Period = resample(data_I1,p,q);
        I2_Resampled_10_Period = resample(data_I2,p,q);
        I3_Resampled_10_Period = resample(data_I3,p,q);
        U1_Resampled_10_Period = resample(data_U1,p,q);
        U2_Resampled_10_Period = resample(data_U2,p,q);
        U3_Resampled_10_Period = resample(data_U3,p,q);
        t_P = data_I1_time(1):(1/F_Sample_New):(data_I1_time(end)+(1/F_Sample_New));
        I1_Resampled_time_10_Period = t_P;
        I2_Resampled_time_10_Period = t_P;
        I3_Resampled_time_10_Period = t_P;
        U1_Resampled_time_10_Period  = t_P;
        U2_Resampled_time_10_Period  = t_P;
        U3_Resampled_time_10_Period  = t_P;
    else if abs(F_Sample_10_Periods) > F_Sample_Complete_Data
            F_Sample_New = F_Sample_10_Periods;
            [p, q] = rat(abs(Delta_Time_Complete_Data/(1/F_Sample_New)));
            I1_Resampled_10_Period = resample(data_I1,p,q);
            I2_Resampled_10_Period = resample(data_I2,p,q);
            I3_Resampled_10_Period = resample(data_I3,p,q);
            U1_Resampled_10_Period = resample(data_U1,p,q);
            U2_Resampled_10_Period = resample(data_U2,p,q);
            U3_Resampled_10_Period = resample(data_U3,p,q);
            for kk=1:length(I1_Resampled_10_Period)-1
                I1_Resample_10_Period_Final(kk) =  ((I1_Resampled_10_Period(kk)+I1_Resampled_10_Period(kk+1))/2);
                I2_Resample_10_Period_Final(kk) =  ((I2_Resampled_10_Period(kk)+I2_Resampled_10_Period(kk+1))/2);
                I3_Resample_10_Period_Final(kk) =  ((I3_Resampled_10_Period(kk)+I3_Resampled_10_Period(kk+1))/2);
                U1_Resample_10_Period_Final(kk) =  ((U1_Resampled_10_Period(kk)+U1_Resampled_10_Period(kk+1))/2);
                U2_Resample_10_Period_Final(kk) =  ((U2_Resampled_10_Period(kk)+U2_Resampled_10_Period(kk+1))/2);
                U3_Resample_10_Period_Final(kk) =  ((U3_Resampled_10_Period(kk)+U3_Resampled_10_Period(kk+1))/2);
            end
            if length(I1_Resample_10_Period_Final) > length(Rampe)
                I1_Resample_10_Period_Final = I1_Resample_10_Period_Final(1:length(Rampe));
                I2_Resample_10_Period_Final = I2_Resample_10_Period_Final(1:length(Rampe));
                I3_Resample_10_Period_Final = I3_Resample_10_Period_Final(1:length(Rampe));
            else if length(I1_Resample_10_Period_Final) < length(Rampe)
                    error('Please check the Resampled data something wrong there :(')
                end
            end
            t_Q = Rampe(1):((q/p)/F_Sample_Complete_Data):Rampe(end);
            if length(t_Q) == length(I1_Resampled_10_Period)
                I1_Resampled_time_10_Period = t_Q;
                I2_Resampled_time_10_Period = t_Q;
                I3_Resampled_time_10_Period = t_Q;
                U1_Resampled_time_10_Period  = t_Q;
                U2_Resampled_time_10_Period  = t_Q;
                U3_Resampled_time_10_Period  = t_Q;
            else if length(t_Q) > length(I1_Resampled_10_Period)
                    I1_Resampled_time_10_Period = t_Q(1:length(I1_Resampled_10_Period));
                    I2_Resampled_time_10_Period = t_Q(1:length(I2_Resampled_10_Period));
                    I3_Resampled_time_10_Period = t_Q(1:length(I3_Resampled_10_Period));
                    U1_Resampled_time_10_Period = t_Q(1:length(U1_Resampled_10_Period));
                    U2_Resampled_time_10_Period = t_Q(1:length(U2_Resampled_10_Period));
                    U3_Resampled_time_10_Period = t_Q(1:length(U3_Resampled_10_Period));
                else if length(t_Q) < length(I1_Resampled_10_Period)
                        I1_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
                        I2_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
                        I3_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
                        U1_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
                        U2_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
                        U3_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
                    end
                end
            end
        else
            error('Please check the Delta Value of the 10 Period signal')
        end
    end
    %------------------------------------------------------------------
    %---------------- Data Re-sampling End ----------------------------
    %------------------------------------------------------------------
    %----------------Spectrum Calculation------------------------------
    %------------------------------------------------------------------
    Tp_10_Period_Resampled = length(I1_Resampled_time_10_Period);
    Delta_Time_10_Period_Resampled = I1_Resampled_time_10_Period(2)-I1_Resampled_time_10_Period(3);
    F_Sample_10_Period_Resampled = abs(1/Delta_Time_10_Period_Resampled);
    T_S_10_Period_Resampled = abs(Tp*Delta_Time_10_Period_Resampled);
    %----------------------------------------------------------------------
    [nyquist_number_I1_Resampled_10_Period,Bezugsfrequenz_I1_Resampled_10_Period,Signal_amp_I1_Resampled_10_Period] = Fast_Harmonics_1 (I1_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
    [nyquist_number_I2_Resampled_10_Period,Bezugsfrequenz_I2_Resampled_10_Period,Signal_amp_I2_Resampled_10_Period] = Fast_Harmonics_1 (I2_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
    [nyquist_number_I3_Resampled_10_Period,Bezugsfrequenz_I3_Resampled_10_Period,Signal_amp_I3_Resampled_10_Period] = Fast_Harmonics_1 (I3_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
    [nyquist_number_U1_Resampled_10_Period,Bezugsfrequenz_U1_Resampled_10_Period,Signal_amp_U1_Resampled_10_Period] = Fast_Harmonics_1 (U1_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
    [nyquist_number_U2_Resampled_10_Period,Bezugsfrequenz_U2_Resampled_10_Period,Signal_amp_U2_Resampled_10_Period] = Fast_Harmonics_1 (U2_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
    [nyquist_number_U3_Resampled_10_Period,Bezugsfrequenz_U3_Resampled_10_Period,Signal_amp_U3_Resampled_10_Period] = Fast_Harmonics_1 (U3_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
    Figure_Show = 0;
    if Figure_Show == 1
        figure;
        subplot(2,1,1);
        plot([0:nyquist_number_data_I1]*Bezugsfrequenz_data_I1,Signal_amp_data_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
        plot([0:nyquist_number_data_I2]*Bezugsfrequenz_data_I2,Signal_amp_data_I2,'-.r');hold on
        plot([0:nyquist_number_data_I3]*Bezugsfrequenz_data_I3,Signal_amp_data_I3,'--k');hold off
        subplot(2,1,2);
        plot([0:nyquist_number_data_U1]*Bezugsfrequenz_data_U1,Signal_amp_data_U1);xlabel('frequency [Hz]'),ylabel('amplitude Voltages [A]');hold on
        plot([0:nyquist_number_data_U2]*Bezugsfrequenz_data_U2,Signal_amp_data_U2,'-.r');hold on
        plot([0:nyquist_number_data_U3]*Bezugsfrequenz_data_U3,Signal_amp_data_U3,'.k','MarkerSize',0.1);hold off
    end
    % Grouping U1
    data_specM_x_H = [0:nyquist_number_U1_Resampled_10_Period]*Bezugsfrequenz_U1_Resampled_10_Period;
    deta_specM_x_V = data_specM_x_H';
    data_specM_Y_V = Signal_amp_U1_Resampled_10_Period';
    data_specM = data_specM_Y_V;
    i = 11;
    ii = 1;
    while i < 502
        harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
        iTxt = (i-1)/10;
        Results_U_H1(ii,jj) = harm_calc;
        if i < 392
            iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
            Results_U_IH1(ii,jj) = iharm_calc;
        end
        if i <= 360
            j1 = (i-1)*4+360+1;
            j_end = j1+40;
            hfharm_calc = 0;
            while j1 <= j_end
                hfharm_calc = hfharm_calc+power(data_specM(j1),2);
                j1 =j1+1;
            end
            hfharm_calc = sqrt(hfharm_calc);
            j1 = 200*(i-1)/10+1900;
            Results_U_HF1(ii,jj) = hfharm_calc;
        end
        i = i+10;
        ii = ii+1;
    end
    % Grouping U2
    data_specM_x_H = [0:nyquist_number_U2_Resampled_10_Period]*Bezugsfrequenz_U2_Resampled_10_Period;
    deta_specM_x_V = data_specM_x_H';
    data_specM_Y_V = Signal_amp_U2_Resampled_10_Period';
    data_specM = data_specM_Y_V;
    i = 11;
    ii = 1;
    while i < 502
        harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
        iTxt = (i-1)/10;
        Results_U_H2(ii,jj) = harm_calc;
        if i < 392
            iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
            Results_U_IH2(ii,jj) = iharm_calc;
        end
        if i <= 360
            j1 = (i-1)*4+360+1;
            j_end = j1+40;
            hfharm_calc = 0;
            while j1 <= j_end
                hfharm_calc = hfharm_calc+power(data_specM(j1),2);
                j1 =j1+1;
            end
            hfharm_calc = sqrt(hfharm_calc);
            j1 = 200*(i-1)/10+1900;
            Results_U_HF2(ii,jj) = hfharm_calc;
        end
        i = i+10;
        ii = ii+1;
    end
    % Grouping U3
    data_specM_x_H = [0:nyquist_number_U3_Resampled_10_Period]*Bezugsfrequenz_U3_Resampled_10_Period;
    deta_specM_x_V = data_specM_x_H';
    data_specM_Y_V = Signal_amp_U3_Resampled_10_Period';
    data_specM = data_specM_Y_V;
    i = 11;
    ii = 1;
    while i < 502
        harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
        iTxt = (i-1)/10;
        Results_U_H3(ii,jj) = harm_calc;
        if i < 392
            iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
            Results_U_IH3(ii,jj) = iharm_calc;
        end
        if i <= 360
            j1 = (i-1)*4+360+1;
            j_end = j1+40;
            hfharm_calc = 0;
            while j1 <= j_end
                hfharm_calc = hfharm_calc+power(data_specM(j1),2);
                j1 =j1+1;
            end
            hfharm_calc = sqrt(hfharm_calc);
            j1 = 200*(i-1)/10+1900;
            Results_U_HF3(ii,jj) = hfharm_calc;
        end
        i = i+10;
        ii = ii+1;
    end
    % Grouping I1
    data_specM_x_H = [0:nyquist_number_I1_Resampled_10_Period]*Bezugsfrequenz_I1_Resampled_10_Period;
    %deta_specM_x_V = data_specM_x_H';
    data_specM_Y_V = Signal_amp_I1_Resampled_10_Period';
    data_specM = data_specM_Y_V;
    i = 11;
    ii = 1;
    while i < 502
        harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
        iTxt = (i-1)/10;
        Results_I_H1(ii,jj) = harm_calc;
        if i < 392
            iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
            Results_I_IH1(ii,jj) = iharm_calc;
        end
        if i <= 360
            j1 = (i-1)*4+360+1;
            j_end = j1+40;
            hfharm_calc = 0;
            while j1 <= j_end
                hfharm_calc = hfharm_calc+power(data_specM(j1),2);
                j1 =j1+1;
            end
            hfharm_calc = sqrt(hfharm_calc);
            j1 = 200*(i-1)/10+1900;
            Results_I_HF1(ii,jj) = hfharm_calc;
        end
        i = i+10;
        ii = ii+1;
    end
    % Grouping I2
    data_specM_x_H = [0:nyquist_number_I2_Resampled_10_Period]*Bezugsfrequenz_I2_Resampled_10_Period;
    %deta_specM_x_V = data_specM_x_H';
    data_specM_Y_V = Signal_amp_I2_Resampled_10_Period';
    data_specM = data_specM_Y_V;
    i = 11;
    ii = 1;
    while i < 502
        harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
        iTxt = (i-1)/10;
        Results_I_H2(ii,jj) = harm_calc;
        if i < 392
            iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
            Results_I_IH2(ii,jj) = iharm_calc;
        end
        if i <= 360
            j1 = (i-1)*4+360+1;
            j_end = j1+40;
            hfharm_calc = 0;
            while j1 <= j_end
                hfharm_calc = hfharm_calc+power(data_specM(j1),2);
                j1 =j1+1;
            end
            hfharm_calc = sqrt(hfharm_calc);
            j1 = 200*(i-1)/10+1900;
            Results_I_HF2(ii,jj) = hfharm_calc;
        end
        i = i+10;
        ii = ii+1;
    end
    % Grouping I3
    data_specM_x_H = [0:nyquist_number_I3_Resampled_10_Period]*Bezugsfrequenz_I3_Resampled_10_Period;
    %deta_specM_x_V = data_specM_x_H';
    data_specM_Y_V = Signal_amp_I3_Resampled_10_Period';
    data_specM = data_specM_Y_V;
    i = 11;
    ii = 1;
    while i < 502
        harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
        iTxt = (i-1)/10;
        Results_I_H3(ii,jj) = harm_calc;
        if i < 392
            iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
            Results_I_IH3(ii,jj) = iharm_calc;
        end
        if i <= 360
            j1 = (i-1)*4+360+1;
            j_end = j1+40;
            hfharm_calc = 0;
            while j1 <= j_end
                hfharm_calc = hfharm_calc+power(data_specM(j1),2);
                j1 =j1+1;
            end
            hfharm_calc = sqrt(hfharm_calc);
            j1 = 200*(i-1)/10+1900;
            Results_I_HF3(ii,jj) = hfharm_calc;
        end
        i = i+10;
        ii = ii+1;
    end
    if 0 == 1
        v_alpha = sqrt(2)/3*(data_U1-0.5*data_U2-0.5*data_U3);
        v_beta = sqrt(2)/3*(sqrt(3)/2*data_U2-sqrt(3)/2*data_U3);
        I_alpha = sqrt(2)/3*(data_I1-0.5*data_I2-0.5*data_I3);
        I_beta = sqrt(2)/3*(sqrt(3)/2*data_I2-sqrt(3)/2*data_I3);
        P = 3*(v_alpha*I_alpha+v_beta*I_beta);
        Q = 3*(-v_alpha*I_beta+v_beta*I_alpha);
        Results_P(jj) = P;
        Results_Q(jj) = Q;
    end
    j = j+20; % Calculation every 10 Periods
    jj = jj+1;
end
%**********************************************************************
%----------------------------------END---------------------------------
%**********************************************************************
%-------------------------- Normalization -----------------------------
%**********************************************************************
Results_U_H1_pu = Results_U_H1/sqrt(2)/Un;
Results_U_H2_pu = Results_U_H2/sqrt(2)/Un;
Results_U_H3_pu = Results_U_H3/sqrt(2)/Un;
Results_I_H1_pu = Results_I_H1/sqrt(2)/In;
Results_I_H2_pu = Results_I_H2/sqrt(2)/In;
Results_I_H3_pu = Results_I_H3/sqrt(2)/In;
%----------------------------------------------------------------------
Results_U_IH1_pu = Results_U_IH1/sqrt(2)/Un;
Results_U_IH2_pu = Results_U_IH2/sqrt(2)/Un;
Results_U_IH3_pu = Results_U_IH3/sqrt(2)/Un;
Results_I_IH1_pu = Results_I_IH1/sqrt(2)/In;
Results_I_IH2_pu = Results_I_IH2/sqrt(2)/In;
Results_I_IH3_pu = Results_I_IH3/sqrt(2)/In;
%----------------------------------------------------------------------
Results_U_HF1_pu = Results_U_HF1/sqrt(2)/Un;
Results_U_HF2_pu = Results_U_HF2/sqrt(2)/Un;
Results_U_HF3_pu = Results_U_HF3/sqrt(2)/Un;
Results_I_HF1_pu = Results_I_HF1/sqrt(2)/In;
Results_I_HF2_pu = Results_I_HF2/sqrt(2)/In;
Results_I_HF3_pu = Results_I_HF3/sqrt(2)/In;
%**********************************************************************
%------------------------ End -----------------------------------------
%**********************************************************************
% ------------------- Mean and Max Value Calculation ------------------
% *********************************************************************
[~, colum] = size(Results_I_H1_pu);
if colum > 1
MeanHarm_U1 = mean(Results_U_H1_pu');
MeanHarm_U2 = mean(Results_U_H2_pu');
MeanHarm_U3 = mean(Results_U_H3_pu');
MeanHarm_I1 = mean(Results_I_H1_pu');
MeanHarm_I2 = mean(Results_I_H2_pu');
MeanHarm_I3 = mean(Results_I_H3_pu');
%----------------------------------------------------------------------
MeanIHarm_U1 = mean(Results_U_IH1_pu');
MeanIHarm_U2 = mean(Results_U_IH2_pu');
MeanIHarm_U3 = mean(Results_U_IH3_pu');
MeanIHarm_I1 = mean(Results_I_IH1_pu');
MeanIHarm_I2 = mean(Results_I_IH2_pu');
MeanIHarm_I3 = mean(Results_I_IH3_pu');
%----------------------------------------------------------------------
MeanHFHarm_U1 = mean(Results_U_HF1_pu');
MeanHFHarm_U2 = mean(Results_U_HF2_pu');
MeanHFHarm_U3 = mean(Results_U_HF3_pu');
MeanHFHarm_I1 = mean(Results_I_HF1_pu');
MeanHFHarm_I2 = mean(Results_I_HF2_pu');
MeanHFHarm_I3 = mean(Results_I_HF3_pu');
%----------------------------------------------------------------------
MaxHarm_U = max(max(MeanHarm_U1,MeanHarm_U2),MeanHarm_U3);
MaxHarm_I = max(max(MeanHarm_I1,MeanHarm_I2),MeanHarm_I3);
MaxIHarm_U = max(max(MeanIHarm_U1,MeanIHarm_U2),MeanIHarm_U3);
MaxIHarm_I = max(max(MeanIHarm_I1,MeanIHarm_I2),MeanIHarm_I3);
MaxHFHarm_U = max(max(MeanHFHarm_U1,MeanHFHarm_U2),MeanHFHarm_U3);
MaxHFHarm_I = max(max(MeanHFHarm_I1,MeanHFHarm_I2),MeanHFHarm_I3);
else
MeanHarm_I1 = Results_I_H1_pu;
MeanHarm_I2 = Results_I_H2_pu;
MeanHarm_I3 = Results_I_H3_pu;
MeanHarm_U1 = Results_U_H1_pu;
MeanHarm_U2 = Results_U_H2_pu;
MeanHarm_U3 = Results_U_H3_pu;
%----------------------------------------------------------------------
MeanIHarm_U1 = Results_U_IH1_pu;
MeanIHarm_U2 = Results_U_IH2_pu;
MeanIHarm_U3 = Results_U_IH3_pu;
MeanIHarm_I1 = Results_I_IH1_pu;
MeanIHarm_I2 = Results_I_IH2_pu;
MeanIHarm_I3 = Results_I_IH3_pu;
%----------------------------------------------------------------------
MeanHFHarm_U1 = Results_U_HF1_pu;
MeanHFHarm_U2 = Results_U_HF2_pu;
MeanHFHarm_U3 = Results_U_HF3_pu;
MeanHFHarm_I1 = Results_I_HF1_pu;
MeanHFHarm_I2 = Results_I_HF2_pu;
MeanHFHarm_I3 = Results_I_HF3_pu;
%----------------------------------------------------------------------
MaxHarm_U = max(max(MeanHarm_U1,MeanHarm_U2),MeanHarm_U3);
MaxHarm_I = max(max(MeanHarm_I1,MeanHarm_I2),MeanHarm_I3);
MaxIHarm_U = max(max(MeanIHarm_U1,MeanIHarm_U2),MeanIHarm_U3);
MaxIHarm_I = max(max(MeanIHarm_I1,MeanIHarm_I2),MeanIHarm_I3);
MaxHFHarm_U = max(max(MeanHFHarm_U1,MeanHFHarm_U2),MeanHFHarm_U3);
MaxHFHarm_I = max(max(MeanHFHarm_I1,MeanHFHarm_I2),MeanHFHarm_I3);
end
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
%--------------- Order (x-axie) Calculation and Ploting ---------------
%**********************************************************************
Order_Harm = 1:1:50;
Order_IHarm = 1:1:39;
%Order_HFHarm = 1:1:35;
%Order_HFHarm = 2100:1:2134;
Order_HFHarm = 2100:200:8900;
if figures
    figure;
    subplot(2,3,1)
    bar(Order_Harm,MaxHarm_I);xlabel('Order Number');ylabel('Current Hamonics Levels I/In in pu');title('Current Harmonics spectrom')
    subplot(2,3,2)
    bar(Order_IHarm,MaxIHarm_I);xlabel('Order Number');ylabel('Current Hamonics Levels I/In in pu');title('Current interharmonics spectrom')
    subplot(2,3,3)
    bar(Order_HFHarm,MaxHFHarm_I);xlabel('Frequency in Hz');ylabel('Current Hamonics Levels I/In in pu');title('Current Higher Frequencies spectrom')
    subplot(2,3,4)
    bar(Order_Harm,MaxHarm_U);xlabel('Order Number');ylabel('Voltage Hamonics Levels U/Un in pu');title('Voltage Harmonics spectrom')
    subplot(2,3,5)
    bar(Order_IHarm,MaxIHarm_U);xlabel('Order Number');ylabel('Voltage Hamonics Levels U/Un in pu');title('Voltage interharmonics spectrom')
    subplot(2,3,6)
    bar(Order_HFHarm,MaxHFHarm_U);xlabel('Frequency in Hz');ylabel('Voltage Hamonics Levels U/Un in pu');title('Voltage Higher Frequencies spectrom')
end
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
%-------------------------- K, THD Factor  ----------------------------
%**********************************************************************
THD_U = sqrt(sum((MaxHarm_U(2:end).^2)'))/(MaxHarm_U(1));
THD_I = sqrt(sum((MaxHarm_I(2:end).^2)'))/(MaxHarm_I(1));
K_U = sqrt(sum((MaxHarm_U(2:end).^2)')/(MaxHarm_U(1)+sum((MaxHarm_U(2:end).^2)')));
K_I = sqrt(sum((MaxHarm_I(2:end).^2)')/(MaxHarm_I(1)+sum((MaxHarm_I(2:end).^2)')));
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
%------------------------ Saving Results ------------------------------
%**********************************************************************
save('Harmonics','THD_U','THD_I','K_U','K_I','MaxHarm_I','MaxHarm_U','MaxIHarm_I','MaxIHarm_U','MaxHFHarm_I','MaxHFHarm_U','MaxHarm_I','Order_IHarm','Order_Harm','Order_HFHarm')
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = surf(m,cdata,varargin)
%
% Syntax
%   surf(m,cdata)
%
% Input
%
% Output
%
% Options
%
% See also
%

varargin = [m.CS.plotOptions,varargin];
[varargout{1:nargout}] = surf@vector3d(m,cdata,varargin{:},m.CS);

author: mtex-toolbox

repo: mtex

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


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function [X_clean, clean_up, remaining_indices, edges] = cleanLiDARTargetWithOneDataSetWithIndices(X_original, target_len, opt)
    N = 4; % clean up using N std for x axis
    M = 3; % clean up using M std for y, z axis
    opt = optimizeCost(opt, X_original, target_len, 0.001);
    if size(X_original, 2) ~= 4
        X = [X_original(1:3, :);
             ones(1 ,size(X_original, 2))];
    end

    X_ref = opt.H_opt * X;
%     distance = sum(X_ref(1:3,:), 1); % L1
    L_infinity = max(abs(X_ref(1:3,:)));  % L infinity
    K = find(L_infinity < (target_len/2)*1.025);
    X_ref_clean_yz = X_ref(:, K); % clean up y and z axis
    L_infinity = max(X_ref_clean_yz(1:3,:));
    K_center = find(L_infinity < target_len/4);
    X_std = std(X_ref_clean_yz(:, K_center), 1, 2);
    Q = find(abs(X_ref(1,:)) < N*(X_std(1))); % clean up x with 2 std
    remaining_indices = intersect(K, Q);
    X_ref_clean = X_ref(:, remaining_indices);
    X_clean = inv(opt.H_opt) * X_ref_clean;
    clean_up.std = N*(X_std(1));
    clean_up.L_infinity = L_infinity;
    clean_up.L_1 = sum(X_ref(1:3,:), 1);
    X_ref_with_all_info = X_original;
    X_ref_with_all_info(1:3, :) = X_ref(1:3,:);
    edges = findEdgePointsInIdealFrame(opt.H_opt, X_ref_with_all_info, target_len);
%     figure(200);
%     clf('reset')
%     scatter3(X_ref(1,:), X_ref(2,:), X_ref(3,:))
%     xlabel('x')
%     ylabel('y')
%     zlabel('z')
%     axis equal
%     hold on;
% %     scatter3(X_ref_clean_yz(1,:), X_ref_clean_yz(2,:), X_ref_clean_yz(3,:))
%     scatter3(X_ref_clean(1,:), X_ref_clean(2,:), X_ref_clean(3,:))
%     axis equal
%     view(90,0)
%     hold off;
end


text:
function [ height, ival, val ] = p09_gx ( option )

%*****************************************************************************80
%
%% P09_GX is used by problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real HEIGHT, the height of the structure.
%
%    Output, integer IVAL, the index of the variable being fixed.
%
%    Output, real VAL, the value of the fixed variable.
%
  height = 2.0;
  ival = 4;

  if ( option == 1 | option == 5 | option == 9 )
    val = 0.00;
  elseif ( option == 2 | option == 6 | option == 10 )
    val = 0.25;
  elseif ( option == 3 | option == 7 | option == 11 )
    val = 0.50;
  elseif ( option == 4 | option == 8 | option == 12 )
    val = 1.00;
  elseif ( option == 13 )
    val = 0.00;
  end

  return
end


text:
function [ft] = um2ft(um)
% Convert length from micrometers (or microns) to feet.
% Chad A. Greene 2012
ft = um*0.000003280839895013;


text:
function [is,os,str] = strengths_dir(CIJ)
%STRENGTHS_DIR      Instrength and outstrength
%
%   [is,os,str] = strengths_dir(CIJ);
%
%   Node strength is the sum of weights of links connected to the node. The
%   instrength is the sum of inward link weights and the outstrength is the
%   sum of outward link weights.
%
%   Input:      CIJ,    directed weighted connection matrix
%
%   Output:     is,     node instrength
%               os,     node outstrength
%               str,    node strength (instrength + outstrength)
%
%   Notes:  Inputs are assumed to be on the columns of the CIJ matrix.
%
%
%   Olaf Sporns, Indiana University, 2002/2006/2008

% compute strengths
is = sum(CIJ,1);    % instrength = column sum of CIJ
os = sum(CIJ,2)';   % outstrength = row sum of CIJ
str = is+os;        % strength = instrength+outstrength

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function [decentr_state, opt_handle, dgs_stats, opt_increment] = ...
    manageAsyncGaussSeidel(...
    decentr_state, opt_handle, distributed_mapper_location, ...
    wait_finished, max_iters, renew)

n_robots = numel(decentr_state);
dgs_stats = {};
opt_increment = zeros(n_robots);

% If an optimization has previously been launched...
if (opt_handle.running)
    % If optimization still running, don't do anything.
    for i = 1:numel(opt_handle.futures)
        if (~isempty(opt_handle.futures{i}))
            if (wait_finished)
                wait(opt_handle.futures{i});
            else
                if (~strcmp(opt_handle.futures{i}.State, 'finished'))
                    return;
                end
            end
        end
    end
    % Catch any failures.
    for i = 1:numel(opt_handle.futures)
        if (~isempty(opt_handle.futures{i}))
            if (~isempty(opt_handle.futures{i}.Error))
                opt_handle.futures{i}.Error
                opt_handle.futures{i}.Diary
                opt_handle.futures{i}.Error.stack(1)
                opt_handle.futures{i}.Error.stack(end)
                assert(false);
            end
        end
    end

    % Optimization finished successfully: Get stats, update state.
    update_state = cell(n_robots, 1);
    for i = 1:numel(opt_handle.futures)
        if (~isempty(opt_handle.futures{i}))
            update_state(opt_handle.group_assignment == i) = ...
                opt_handle.futures{i}.OutputArguments{1};
            dgs_stats_i = opt_handle.futures{i}.OutputArguments{2};
            dgs_stats = [dgs_stats; dgs_stats_i];

            % Reconstruct communication between pairs of robots.
            dgs_stats_i.comm_link_count
            opt_handle.between_robot_links{i}
            links = opt_handle.between_robot_links{i};
            links = links + links';
            linksum = sum(links(:));
            assert(dgs_stats_i.comm_link_count == linksum);

            group = find(opt_handle.group_assignment == i);
            group_data_exchange = ...
                links * dgs_stats_i.exchange_gs / linksum;
            opt_increment(group, group) = ...
                opt_increment(group, group) + group_data_exchange;
        end
    end

    decentr_state = consistentlyApplyState(update_state, decentr_state);
    decentr_state = setConvergedWhereApplicable(...
        update_state, decentr_state);

    % Indicate which groups are co-optimized, for accuracy eval:
    opt_handle.optimized_groups = [];
    groups = unique(opt_handle.group_assignment);
    for group_i = 1:numel(groups)
        opt_handle.optimized_groups(group_i).members = find(...
            opt_handle.group_assignment == groups(group_i));
    end
    disp (['Have ' num2str(numel(opt_handle.optimized_groups)) ...
        ' optimization groups.']);
end

if renew
    % Launch new optimization
    %opt_future = parfeval(gcp(), @runSyncGaussSeidel, 1, decentr_state, ...
    %    distributed_mapper_location);

    opt_handle.group_assignment = 1:n_robots;
    for i = 1:n_robots
        opt_handle.group_assignment(i) = ...
            min([decentr_state{i}.grouped_with i]);
    end

    opt_groups = unique(opt_handle.group_assignment);

    opt_handle.futures = {};
    for group_i = opt_groups
        group_mask = (opt_handle.group_assignment == group_i);
        group_size = nnz(group_mask);
        if group_size > 1
            members_converged = cellfun(@(x) x.converged, ...
                decentr_state(group_mask));
            if ~all(members_converged)
                group_reindexing = zeros(n_robots, 1);
                group_reindexing(group_mask) = 1:group_size;
                opt_handle.futures{group_i} = parfeval(gcp(), ...
                    @runSyncGaussSeidel, 2, decentr_state(group_mask), ...
                    distributed_mapper_location, group_i, group_reindexing, ...
                    max_iters);
                opt_handle.between_robot_links{group_i} = ...
                    countBetweenRobotLinks(decentr_state(group_mask), ...
                    group_reindexing);
            end
        end
    end
    opt_handle.running = true;
end

end

author: uzh-rpg


text:
function [nVertices, nFaces, BstMat] = out_tess_tri( BstFile, OutputFile, isOpenMEEG )
% OUT_TESS_TRI: Exports a surface to a BrainVISA ASCII .tri file.
%
% USAGE:  [nVertices, nFaces] = out_tess_tri( BstFile, OutputFile )
%         [nVertices, nFaces] = out_tess_tri( BstFile/TessMat, OutputFile, isOpenMEEG=0 )
%
% INPUT:
%    - BstFile    : full path to Brainstorm file to export
%    - OutputFile : full path to output file (with '.tri' extension)
%    - isOpenMEEG : if flag set to 1, write the positions in meters instead of millimeters

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2011-2019

if (nargin < 3) || isempty(isOpenMEEG)
    isOpenMEEG = 0;
end

% ===== LOAD BRAINSTORM SURFACE =====
if ischar(BstFile)
    BstMat = in_tess_bst(BstFile);
else
    BstMat = BstFile;
end

% ===== PREPARE VALUES ======
% Vertices (=> in millimeters)
if ~isOpenMEEG
    BstMat.Vertices = BstMat.Vertices * 1000;
end
% Faces : remove 1 (convert to 0-based indices)
Faces = BstMat.Faces - 1;
% Normals
VertNormals = BstMat.VertNormals;
% Return surface sizes
nVertices = length(BstMat.Vertices);
nFaces = length(BstMat.Faces);

% ===== SAVE FILE =====
% Open file
[fid, message] = fopen(OutputFile, 'w');
if (fid < 0)
    error(['Could not create file : ' message]);
end
% Write vertices and normals
fprintf(fid, '- %g\n', size(BstMat.Vertices,1));
fprintf(fid, '%g %g %g %g %g %g\n', [BstMat.Vertices, VertNormals]');
% Write faces
nfaces = size(Faces,1);
fprintf(fid, '- %g %g %g\n', [nfaces nfaces nfaces]);
fprintf(fid, '%g %g %g\n', Faces');
% Close file
fclose(fid);

author: brainstorm-tools


text:
function [LE,J]=CylNeumanLePerp_PGSE(d, R, G, delta, smalldel, roots)
% Substrate: Parallel, impermeable cylinders with one radius in an empty
%            background.
% Pulse sequence: Pulsed gradient spin echo
% Signal approximation: Gaussian phase distribution.
%
% [LE,J] = CylNeumanLePerp_PGSE(d, R, G, delta, smalldel, roots)
%
% returns the log signal attenuation in perpendicular direction (LePerp) for
% EACH RADIUS specified in R according to the Neuman model and the Jacobian J
% of LePerp with respect to the parameters.
%
% The Jacobian DOES NOT include derivates with respect to the fibre direction.
%
% d is the diffusivity of the material inside the cylinders.
%
% R is the list of the radii of the cylinders. It has size [1 m] where m is the
% number of radii.
%
% G, delta and smalldel are the gradient strength, pulse separation and
% pulse length of each measurement in the protocol.  Each has
% size [N 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%         Gary Hui Zhang     (gary.zhang@ucl.ac.uk)
%

% When R=0, no need to do any calculation
if (R == 0.00)
    LE = zeros(size(G,1), size(R,2));
    J = zeros([size(LE), 2]);
	 J(:,:,:) = 0;
    return;
end

% Check the roots array is correct
if(abs(roots(1) - 1.8412)>0.0001)
    error('Looks like the roots array is wrong.  First value should be 1.8412, but is %f', roots(1));
end

% Radial wavenumbers
GAMMA = 2.675987E8;

% number of gradient directions, i.e. number of measurements
l_q=size(G,1);
l_a=numel(R);
k_max=numel(roots);

R_mat=repmat(R,[l_q 1]);
R_mat=R_mat(:);
R_mat=repmat(R_mat,[1 1 k_max]);
R_matSq=R_mat.^2;

root_m=reshape(roots,[1 1 k_max]);
alpha_mat=repmat(root_m,[l_q*l_a 1 1])./R_mat;
amSq=alpha_mat.^2;
amP6=amSq.^3;

deltamx=repmat(delta,[1,l_a]);
deltamx_rep = deltamx(:);
deltamx_rep = repmat(deltamx_rep,[1 1 k_max]);

smalldelmx=repmat(smalldel,[1,l_a]);
smalldelmx_rep = smalldelmx(:);
smalldelmx_rep = repmat(smalldelmx_rep,[1 1 k_max]);

Gmx=repmat(G,[1,l_a]);
GmxSq = Gmx.^2;

% Perpendicular component (Neuman model)
sda2 = smalldelmx_rep.*amSq;
bda2 = deltamx_rep.*amSq;
emdsda2 = exp(-d*sda2);
emdbda2 = exp(-d*bda2);
emdbdmsda2 = exp(-d*(bda2 - sda2));
emdbdpsda2 = exp(-d*(bda2 + sda2));

sumnum1 = 2*d*sda2;
% the rest can be reused in dE/dR
sumnum2 = - 2 + 2*emdsda2 + 2*emdbda2;
sumnum2 = sumnum2 - emdbdmsda2 - emdbdpsda2;
sumnum = sumnum1 + sumnum2;

sumdenom = d^2*amP6.*(R_matSq.*amSq - 1);

% Check for zeros on top and bottom
%sumdenom(find(sumnum) == 0) = 1;
sumterms = sumnum./sumdenom;

testinds = find(sumterms(:,:,end)>0);
test = sumterms(testinds,1)./sumterms(testinds,end);
if(min(test)<1E4)
    warning('Ratio of largest to smallest terms in Neuman model sum is <1E4.  May need more terms.');
end

s = sum(sumterms,3);
s = reshape(s,[l_q,l_a]);
if(min(s)<0)
    warning('Negative sums found in Neuman sum.  Setting to zero.');
    s(find(s<0))=0;
end

LE = -2*GAMMA^2*GmxSq.*s;

% Compute the Jacobian matrix
if(nargout>1)

    % dLE/dd
    sumnumD = 2*sda2;
    sumnumD = sumnumD - 2*sda2.*emdsda2;
    sumnumD = sumnumD - 2*bda2.*emdbda2;
    sumnumD = sumnumD + (bda2 - sda2).*emdbdmsda2;
    sumnumD = sumnumD + (bda2 + sda2).*emdbdpsda2;
    sumtermsD = sumnumD./sumdenom;

    sD = sum(sumtermsD,3);
    sD = reshape(sD,[l_q,l_a]);

    dLEdd = -2*GAMMA^2*GmxSq.*(sD - 2*s/d);

    % dLE/dR
    sumtermsR = (6*sumterms - 2*d*sumtermsD)./R_mat;

    sR = sum(sumtermsR,3);
    sR = reshape(sR,[l_q,l_a]);

    dLEdr = -2*GAMMA^2*GmxSq.*sR;

    % Construct the jacobian matrix.
    J = zeros([size(LE), 2]);
    J(:,:,1) = dLEdd;
    J(:,:,2) = dLEdr;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function medHi=medianHi(x)
%%MEDIANHI Find the high (upper) median of a set of real values. That is if
%          the n values were sorted, then return the value in index
%          fix(n/2)+1.
%
%INPUTS: x An nX1 or 1Xn vector. The elements do not need to be in any
%          particular order.
%
%OUTPUTS: medHi The value in x that holds the high median.
%
%The function kthOrderStat is used to obtain the value without sorting.
%
%EXAMPLE:
%We look at what one obtains with an odd and an even number of values.
% medHi10=medianHi(1:10)
% medHi11=medianHi(1:11)
%One gets medHi10=6 and medHi11=6.
%
%August 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=length(x);

k=fix(n/2)+1;
medHi=kthOrderStat(x,k);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% demo_wavedet_3D_detector.m
% OVERVIEW This script demonstrates the usage of the wavedet_3D_detector function.
% The wavedet_3D_detector function detects QRS onset, Q, R, S, J fiducial points for ECG required to perform MV analysis for ECG.
%
% The inputs and outputs to the main wavedet_3D_detector function are listed below.
% INPUTS        MANDATORY           DESCRIPTION
%               sig_all (uV)            Single channel ECG in mV.
%
%               ext_anot                The reference R peaks generated using jsqi. Provided to assist in beat detection.
%
%               heasig                  Structure variable which contains
%                                       the following fields.
%               .nsamp                  Length of ecg data in sig_all in
%                                       samples.
%               .freq                   Sampling frequency for data in
%                                       sig_all.
%               .nsig                   Number of ecg channels in sig_all. Set to 1.
%
%               messages                            Strucuture variable with following field
%               .setup.wavedet.QRS_detection_only   Peforms fiducial point detection for the above listed
%                                                   points of interest when set to 1.
%
%
% OUTPUTS (of interest)
%               position                Structure variable which contains
%                                       the following fields of interest.
%               .QRSon                  QRS onset location for each beat in
%                                       sig_all variable.
%               .Q                      Q point locations for each beat in
%                                       the sig_all variable.
%               .R                      R point locations for each beat in
%                                       the sig_all variable.
%               .S                      S point locations for each beat in
%                                       the sig_all variable.
%               .QRSoff                 QRS offset locations for each beat
%                                       in the sig_all variable.
%               .Toff                   T offset locations for each beat in
%                                       the sig_all variable.
%
% The detected fiducial points may be plotted along with the ecg to verify
% the demo executed correctly.
%
%   REF:
%       The wavedet_3D_detector algorithm has been adapted from the
%       opensource ecgkit toolbox referenced below,
%       Demski AJ, Llamedo Soria M. ecg-kit a Matlab Toolbox for Cardiovascular Signal Processing.
%       Journal of Open Research Software. 2016;4(1):e8. DOI: http://doi.org/10.5334/jors.86
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.
%

% add toolbox dependencies
addpath(genpath('../../../../../PhysioNet-Cardiovascular-Signal-Toolbox-master/'));
addpath('../Tools/Annotation_generator/');
% load ecg data
sig = load('./testdata/mvm/sample_ecg_datam');
siginfo = readheader('./testdata/mvm/sample_ecg_datam.hea');
Fs = siginfo.freq;
ecg = (sig.val - siginfo.adczero)./siginfo.gain; % Adjust signal according to gain and dc offset
clear siginfo sig;

% generate annotations
HRVparams = InitializeHRVparams('MVanalysis');
[qrs_pos,sign,en_thres] = jqrs(ecg,HRVparams);
ECG_header.nsig = 1; ECG_header.freq = Fs; ECG_header.nsamp = length(ecg);
wavedet_config.setup.wavedet.QRS_detection_only = 0;
[ann,~,~] = wavedet_3D_detector(ecg, qrs_pos', ECG_header, wavedet_config );
% improve fiducial point detection
[detection_flg,ann] = improvfiducials(ann, Fs, ecg);

if(~detection_flg)
    disp('Warning: Unable to improve fiducial point detection.')
end

% Plotting annotations for verification
figure(2); plot(ecg); hold on;
stem(ann.QRSon, ecg(ann.QRSon)); stem(ann.Q, ecg(ann.Q));
stem(ann.R, ecg(ann.R)); stem(ann.S, ecg(ann.S));
stem(ann.QRSoff, ecg(ann.QRSoff)); stem(ann.Toff(~isnan(ann.Toff)), ecg(ann.Toff(~isnan(ann.Toff)))); hold off
legend('ECG','QRSon','Q','R','S','QRSoff','Toff')

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox