text:
% Minimize transition bandwidth of a linear phase lowpass FIR filter
% "Filter design" lecture notes (EE364) by S. Boyd
% (figures are generated)
%
% Designs a linear phase FIR lowpass filter such that it:
% - minimizes the transition band width (i.e. minimize w_stop)
% - has a constraint on the maximum passband ripple
% - has a constraint on the maximum stopband attenuation
%
% This is a quasiconvex problem and is solved using a bisection.
%
%   minimize   w_stop
%       s.t.   1/delta <= H(w) <= delta     for w in the passband
%              |H(w)| <= atten_level        for w in the stopband
%
% where H is the frequency response function and variable is
% the filter impulse response h (and its order/length).
% Data is delta (max passband ripple) and atten_level (max stopband
% attenuation level).
%
% Written for CVX by Almir Mutapcic 02/02/06

%********************************************************************
% user's filter specifications
%********************************************************************
% starting point for the stopband (needs to be feasible)
wstop = 0.24*pi;        % stopband start freq (in radians)
TOL = 1e-3;             % precision to which we should run bisection

n = 10;                 % filter order (2n+1 is the full order)
wpass = 0.12*pi;        % passband cutoff freq (in radians)
delta = 1;              % max (+/-) passband ripple in dB
atten_level = -30;      % stopband attenuation level in dB

%********************************************************************
% create optimization parameters
%********************************************************************
m = 30*n; % freq samples (rule-of-thumb)
w = linspace(0,pi,m);

%*********************************************************************
% use bisection algorithm to solve the problem
%*********************************************************************

wstop_bot = wpass;
wstop_top  = wstop;

while( wstop_top - wstop_bot > TOL)
  % try to find a feasible design for given specs
  wstop_cur = (wstop_top + wstop_bot)/2;

  % create optimization matrices (matrix of cosines)
  A = [ones(m,1) 2*cos(kron(w',[1:n]))];

  % passband 0 <= w <= w_pass
  ind = find((0 <= w) & (w <= wpass));    % passband
  Ap  = A(ind,:);

  % transition band is not constrained (w_pass <= w <= w_stop)

  % stopband (w_stop <= w) (this is the changing constraint)
  ind = find((wstop_cur <= w) & (w <= pi));   % stopband
  As  = A(ind,:);

  % formulate and solve the feasibility linear-phase lp filter design
  cvx_begin quiet
    variable h_cur(n+1,1);
    % feasibility problem
    % passband bounds
    Ap*h_cur <= 10^(delta/20);
    Ap*h_cur >= 10^(-delta/20);
    % stopband bounds
    abs( As*h_cur ) <= 10^(atten_level/20);
  cvx_end

  % bisection
  if strfind(cvx_status,'Solved') % feasible
    fprintf(1,'Problem is feasible for stopband freq = %3.4f rads\n',wstop_cur);
    wstop_top = wstop_cur;
    % construct the full impulse response
    h = [flipud(h_cur(2:end)); h_cur];
  else % not feasible
    fprintf(1,'Problem is not feasible for stopband freq = %3.4f rads\n',wstop_cur);
    wstop_bot = wstop_cur;
  end
end

wstop = wstop_top;
fprintf(1,['\nOptimum stopband frequency for given specs is %3.4f*pi rads\n' ...
           'and the minimum transition width is %3.4f*pi radians.\n'],...
            wstop/pi, (wstop-wpass)/pi);

%********************************************************************
% plots
%********************************************************************
figure(1)
% FIR impulse response
plot([0:2*n],h','o',[0:2*n],h','b:')
xlabel('t'), ylabel('h(t)')

figure(2)
% frequency response
H = exp(-j*kron(w',[0:2*n]))*h;
% magnitude
subplot(2,1,1)
plot(w,20*log10(abs(H)),...
     [wstop pi],[atten_level atten_level],'r--',...
     [0 wpass],[delta delta],'r--',...
     [0 wpass],[-delta -delta],'r--');
axis([0,pi,-40,10])
xlabel('w'), ylabel('mag H(w) in dB')
% phase
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w'), ylabel('phase H(w)')

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ as, bs, cs ] = sphere_triangle_vertices_to_sides ( r, v1, v2, v3 )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_VERTICES_TO_SIDES computes spherical triangle sides.
%
%  Discussion:
%
%    We can use the ACOS system call here, but the ARC_COSINE routine
%    will automatically take care of cases where the input argument is
%    (usually slightly) out of bounds.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real V1(3), V2(3), V3(3), the vertices of the spherical
%    triangle.
%
%    Output, real AS, BS, CS, the (geodesic) length of the sides
%    of the triangle.
%
  dim_num = 3;

  as = r * r8_acos ( v2(1:dim_num)' * v3(1:dim_num) / r.^2 );
  bs = r * r8_acos ( v3(1:dim_num)' * v1(1:dim_num) / r.^2 );
  cs = r * r8_acos ( v1(1:dim_num)' * v2(1:dim_num) / r.^2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plotIPDF(ori,varargin)
% plot orientations into inverse pole figures
%
% Syntax
%   plotIPDF(ori,[r1,r2,r3])
%   plotIPDF(ori,[r1,r2,r3],'points',100)
%   plotIPDF(ori,[r1,r2,r3],'points','all')
%   plotIPDF(ori,[r1,r2,r3],'contourf')
%   plotIPDF(ori,[r1,r2,r3],'antipodal')
%   plotIPDF(ori,data,[r1,r2,r3])
%
% Input
%  ebsd - @EBSD
%  r   - @vector3d specimen directions
%
% Options
%  RESOLUTION - resolution of the plots
%  property   - user defined colorcoding
%  MarkerSize -
%  MarkerFaceColor -
%  MarkerEdgeColor -
%
% Flags
%  antipodal - include <VectorsAxes.html antipodal symmetry>
%  complete  - ignore fundamental region
%  upper     - restrict to upper hemisphere
%  lower     - restrict to lower hemisphere
%  filled    - fill the marker with current color
%
% See also
% S2Grid/plot savefigure Plotting Annotations_demo ColorCoding_demo PlotTypes_demo
% SphericalProjection_demo

[mtexFig,isNew] = newMtexFigure('datacursormode',@tooltip,varargin{:});

% maybe we should call this function with the option add2all
if ~isNew && ~check_option(varargin,'parent') && ...
    ((((ishold(mtexFig.gca) && nargin > 1 && isa(varargin{1},'vector3d') && length(varargin{1})>1))) || check_option(varargin,'add2all'))
  plot(ori,varargin{:},'add2all');
  return
end

% extract data
if check_option(varargin,'property')
  data = get_option(varargin,'property');
  data = reshape(data,[1,length(ori) numel(data)/length(ori)]);
elseif nargin > 1 && ~isa(varargin{1},'vector3d')
  [data,varargin] = extract_data(length(ori),varargin);
  data = reshape(data,[1,length(ori) numel(data)/length(ori)]);
else
  data = [];
end

% find inverse pole figure direction
r = [];
try r = getappdata(mtexFig.currentAxes,'inversePoleFigureDirection'); end
if isempty(r), r = varargin{1}; end
argin_check(r,'vector3d');

%  subsample if needed
if (length(ori)*numSym(ori.CS)*numSym(ori.SS) > 100000 || check_option(varargin,'points')) ...
    && ~check_option(varargin,{'all','contourf','smooth','contour','pcolor'})

  points = fix(get_option(varargin,'points',100000/numSym(ori.CS)/numSym(ori.SS)));
  disp(['  I''m plotting ', int2str(points) ,' random orientations out of ', int2str(length(ori)),' given orientations']);

  samples = discretesample(length(ori),points);
  ori = ori.subSet(samples);
  if ~isempty(data), data = data(:,samples,:); end

end

for ir = 1:length(r)

  if ir>1, mtexFig.nextAxis; end

  % the crystal directions
  rSym = symmetrise(r(ir),ori.SS);
  h = ori(:) \ rSym;

  %  plot
  [~,cax] = h.plot(repmat(data,1,length(rSym)),'symmetrised',...
    'fundamentalRegion','doNotDraw',varargin{:});
  if isNew, mtexTitle(cax(1),char(r(ir),'LaTeX')); end

  % plot annotations
  setappdata(cax,'inversePoleFigureDirection',r(ir));
  set(cax,'tag','ipdf');
  setappdata(cax,'CS',ori.CS);
  setappdata(cax,'SS',ori.SS);

  % TODO: unifyMarkerSize

end

if isNew || check_option(varargin,'figSize')
  mtexFig.drawNow('figSize',getMTEXpref('figSize'),varargin{:});
end

% --------------- Tooltip function ------------------
function txt = tooltip(varargin)

  [r_local,id,value] = getDataCursorPos(mtexFig,length(ori));

  %id = (id-1)/

  h_local = round(Miller(r_local,ori.CS));

  txt{1} = ['id = ' xnum2str(id)];
  txt{2} = ['(h,k,l) = ' char(h_local,'tolerance',3*degree,'commasep')];
  txt{3} = ['Euler = ' char(ori.subSet(id))];
  if ~isempty(value)
    txt{4} = ['value = ' xnum2str(value)];
  end

end

end

author: mtex-toolbox

repo: mtex

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


text:
function Population = STM(Population,W,z,znad)
% Selection based on STM model

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    N  = length(Population);
    NW = size(W,1);

    %% The modified Tchebycheff value of each solution on each subproblem
    g = zeros(N,NW);
    for i = 1 : N
        g(i,:) = max(repmat(abs(Population(i).obj-z),NW,1)./W,[],2)';
    end

    %% The perpendicular distance of each solution on each subproblem
    PopObj   = (Population.objs-repmat(z,N,1))./repmat(znad-z,N,1);
    Cosine   = 1 - pdist2(PopObj,W,'cosine');
    Distance = repmat(sqrt(sum(PopObj.^2,2)),1,NW).*sqrt(1-Cosine.^2);

    %% STM selection
    Fp  = zeros(1,NW);
    FX  = zeros(1,N);
    Phi = false(NW,N);
    while any(Fp==0)
        RemainW  = find(Fp==0);
        i        = RemainW(randi(length(RemainW)));
        RemainX  = find(~Phi(i,:));
        [~,best] = min(g(RemainX,i));
        j        = RemainX(best);
        Phi(i,j) = true;
        if FX(j) == 0
            Fp(i) = j;
            FX(j) = i;
        elseif Distance(j,i) < Distance(j,FX(j))
            Fp(i)     = j;
            Fp(FX(j)) = 0;
            FX(j)     = i;
        end
    end
    Population = Population(Fp);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function m = medianw( x, w, dim )
% Fast weighted median.
%
% Computes the weighted median of a set of samples.
%  http://en.wikipedia.org/wiki/Weighted_median
% "A weighted median of a sample is the 50% weighted percentile."
% For matrices computes median along each column (or dimension dim).
% If all weights are equal to 1 gives identical results to median.
%
% USAGE
%  m = medianw( x, w, [dim] )
%
% INPUTS
%  x      - vector or array of samples
%  w      - vector or array of weights
%  dim    - dimension along which to compute median
%
% OUTPUTS
%  m      - weighted median value of x
%
% EXAMPLE - simple toy example
%  x=[1 2 3]; w=[1 1 5]; medianw(x,w)
%
% EXAMPLE - comparison to median
%  n=randi(100); m=randi(100);
%  x=rand(n,m); w=ones(n,m);
%  m1=median(x); m2=medianw(x,w);
%  assert(isequal(m1,m2))
%
% See also median
%
% Piotr's Image&Video Toolbox      Version 3.24
% Copyright 2013 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if(nargin<3), dim=find(size(x)~=1,1); end
d=dim; nd=ndims(x); n=numel(x);

if( n==1 || size(x,d)==1 )
  m=x;
elseif( length(x)==n )
  [x,o]=sort(x); w=w(o); w=cumsum(w);
  w=w/w(end); [~,j]=min(w<=.5);
  if(j==1 || w(j-1)~=.5), m=x(j);
  else m=(x(j-1)+x(j))/2; end
else
  if(d>1), p=[d 1:d-1 d+1:nd]; x=permute(x,p); w=permute(w,p); end
  [x,o]=sort(x); w=w(o); w=cumsum(w); is={':'}; is=is(ones(1,nd-1));
  w=bsxfun(@rdivide,w,w(end,is{:})); [~,j]=min(w<=.5);
  s=size(x); s=reshape(((1:n/s(1))-1)*s(1),size(j));
  j0=max(1,j-1); j0=j0+s; j=j+s;
  same=w(j0)~=.5; j0(same)=j(same); m=(x(j0)+x(j))/2;
  if(d>1), p=[2:d 1 d+1:nd]; m=permute(m,p); end
end

end

author: Cloud-CV

repo: object-proposals


text:
function [out]=ImageRecover(y,D,CoefMatrix)

% method:  Sliding
% CoefMatrix  256*146405
% D 64*256

[N M]=size(y);
n=sqrt(size(D,1));
out=zeros(N,M);
weight=zeros(N,M);

i=1; j=1;
for k=1:(N-n+1)*(M-n+1)
    patch=reshape(D*CoefMatrix(:,k),[n,n]);
    out(i:i+n-1,j:j+n-1)=out(i:i+n-1,j:j+n-1)+patch;
    weight(i:i+n-1,j:j+n-1)=weight(i:i+n-1,j:j+n-1)+1;
    if i<N-n+1
        i=i+1;
    else
        i=1; j=j+1;
    end;
end;
out=out./weight;
return;

author: chongyangtao

repo: Color-Image-Inpainting


text:
function res = rcnn_test_bbox_regressor(imdb, rcnn_model, bbox_reg, suffix)
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

conf = rcnn_config('sub_dir', imdb.name);
image_ids = imdb.image_ids;

% assume they are all the same
feat_opts = bbox_reg.training_opts;
num_classes = length(rcnn_model.classes);

if ~exist('suffix', 'var') || isempty(suffix)
  suffix = '_bbox_reg';
else
  if suffix(1) ~= '_'
    suffix = ['_' suffix];
  end
end

try
  aboxes = cell(num_classes, 1);
  for i = 1:num_classes
    load([conf.cache_dir rcnn_model.classes{i} '_boxes_' imdb.name suffix]);
    aboxes{i} = boxes;
  end
catch
  aboxes = cell(num_classes, 1);
  box_inds = cell(num_classes, 1);
  for i = 1:num_classes
    load([conf.cache_dir rcnn_model.classes{i} '_boxes_' imdb.name]);
    aboxes{i} = boxes;
    box_inds{i} = inds;
    clear boxes inds;
  end

  for i = 1:length(image_ids)
    fprintf('%s: bbox reg test (%s) %d/%d\n', procid(), imdb.name, i, length(image_ids));
    d = rcnn_load_cached_pool5_features(feat_opts.cache_name, ...
        imdb.name, image_ids{i});
    if isempty(d.feat)
      continue;
    end

    d.feat = rcnn_pool5_to_fcX(d.feat, feat_opts.layer, rcnn_model);
    d.feat = rcnn_scale_features(d.feat, feat_opts.feat_norm_mean);

    if feat_opts.binarize
      d.feat = single(d.feat > 0);
    end

    for j = 1:num_classes
      I = box_inds{j}{i};
      boxes = aboxes{j}{i};
      if ~isempty(boxes)
        scores = boxes(:,end);
        boxes = boxes(:,1:4);
        assert(sum(sum(abs(d.boxes(I,:) - boxes))) == 0);
        boxes = rcnn_predict_bbox_regressor(bbox_reg.models{j}, d.feat(I,:), boxes);
        boxes(:,1) = max(boxes(:,1), 1);
        boxes(:,2) = max(boxes(:,2), 1);
        boxes(:,3) = min(boxes(:,3), imdb.sizes(i,2));
        boxes(:,4) = min(boxes(:,4), imdb.sizes(i,1));
        aboxes{j}{i} = cat(2, single(boxes), single(scores));

        if 0
          % debugging visualizations
          im = imread(imdb.image_at(i));
          keep = nms(aboxes{j}{i}, 0.3);
          for k = 1:min(10, length(keep))
            if aboxes{j}{i}(keep(k),end) > -0.9
              showboxes(im, aboxes{j}{i}(keep(k),1:4));
              title(sprintf('%s %d score: %.3f\n', rcnn_model.classes{j}, ...
                  k, aboxes{j}{i}(keep(k),end)));
              pause;
            end
          end
        end
      end
    end
  end

  for i = 1:num_classes
    save_file = [conf.cache_dir rcnn_model.classes{i} '_boxes_' imdb.name suffix];
    boxes = aboxes{i};
    inds = box_inds{i};
    save(save_file, 'boxes', 'inds');
    clear boxes inds;
  end
end

% ------------------------------------------------------------------------
% Peform AP evaluation
% ------------------------------------------------------------------------
for model_ind = 1:num_classes
  cls = rcnn_model.classes{model_ind};
  try
    ld = load([conf.cache_dir cls '_pr_' imdb.name suffix]);
    fprintf('!!! %s : %.4f %.4f\n', cls, ld.res.ap, ld.res.ap_auc);
    res(model_ind) = ld.res;
  catch
    res(model_ind) = imdb.eval_func(cls, aboxes{model_ind}, imdb, suffix);
  end
end

fprintf('\n~~~~~~~~~~~~~~~~~~~~\n');
fprintf('Results (bbox reg):\n');
aps = [res(:).ap]';
disp(aps);
disp(mean(aps));
fprintf('~~~~~~~~~~~~~~~~~~~~\n');

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function MaxDD = maxdrawdown(Data)
%MAXDRAWDOWN Calculate maximum drawdown for one or more total return price series

[T, N] = size(Data);

MaxDD = zeros(1,N);

for i = (T - 1):-1:1
	for j = 1:N
		DD(j) = (min(Data(i:T,j)) - Data(i,j))/Data(i,j);
	end
	p = find(DD < MaxDD);
	MaxDD(p) = DD(p);
end

author: Sable

repo: mcbench-benchmarks


text:
function kds1 = kdelsum1QPOTTS(IESM)
kd=zeros(1,8);
if IESM(2,2)==IESM(1,1);kd(1,1)=1;end;
if IESM(2,2)==IESM(1,2);kd(1,2)=1;end;
if IESM(2,2)==IESM(1,3);kd(1,3)=1;end;
if IESM(2,2)==IESM(2,3);kd(1,4)=1;end;
if IESM(2,2)==IESM(3,3);kd(1,5)=1;end;
if IESM(2,2)==IESM(3,2);kd(1,6)=1;end;
if IESM(2,2)==IESM(3,1);kd(1,7)=1;end;
if IESM(2,2)==IESM(2,1);kd(1,8)=1;end;
kds1=sum(sum(kd));

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34784-monte-carlo-simulation-of-two-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 2D square-lattice - microstructure/kdelsum1QPOTTS.m


text:
%compute length of inflhead vector

function [data,units]=compute_inflheadmag(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
inflheadmag=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    inflheadmag{1,i}=bsxfun(@hypot,trx(larva).xhead_mm-trx(larva).xinflection_mm,trx(larva).yhead_mm-trx(larva).yinflection_mm);
end

units=parseunits('mm');
data=inflheadmag;

author: kristinbranson

repo: JAABA


text:
function main(varargin)
%GUI for viewing the Mandelbrot and the monic quadratic Julia sets.
%
%   MAIN ('PropertyName',PropertyValue,...)
%
%   List of properties
%   FractalType
%      {'mandelbrot','julia'}  default 'mandelbrot'
%      Choose between the Mandelbrot or the monic quadratic Julia set
%   JuliaConstant
%      any complex value       default 0.285+0.01i
%      This is the constant 'c' used for generating the monic quadratic
%      Julia set with
%      z_n+1 = z_n + c
%   MaxIter
%      positive value          default 256
%      The maximum number of iterations for determining if a certain point
%      has diverged or not
%   EscapeRadius
%      positive value          default 2.0
%      This is used to determine when a point has diverged
%   Width
%      positive value          default 400
%      Determines the number of datapoints to use width-wise
%   Height
%      positive value          default 400
%      Determines the number of datapoints to use height-wise
%   Subfunction
%      {'c','matlab'}          default 'c'
%      Type of fractal generator subfunction to use. Either the C or the
%      matlab versions.
%   XLim
%      2-element vector        default [-2,2]
%      The initial region boundaries
%   YLim
%      2-element vector        default [-2,2]
%      The initial region boundaries
%   ImageCacheSize
%      positive value          default 8
%      This program stores the previous ImageCacheSize-2 images. Must be at
%      least 2.
%   NonDivergingSetColor
%      1x3 matrix              default [0,0,0]
%      Color of the non-diverging points in RGB components. Must be between
%      0.0 and 1.0.
%   Colormap
%      function handle         default @hsv
%      Function for determining the colormap. For example, @hsv. The
%      function must have the form out=function(in) where in is the number
%      of colors in the colormap and out is an [in x 3] matrix giving the
%      RGB components of the colormap. The RGB values must be scaled
%      between 0.0 and 1.0.
%
%   NumColor
%      positive value          default 255
%      Number of colors to use when coloring the diverging points.
%
%  This version is from October 28, 2010
%  October 28, 2010, fixed non-complex constant bug for makejulia_c.c
%  October 27, 2010, added inputparser and inhibited rotate3d
%  October 25, 2010, improved image cache and added more comments
%  October 23, 2010, first version
%  Written by Christopher Leung
%  christopher.leung@mail.mcgill.ca
%
%  In memory to Benoit Mandelbrot

% Parse input
p = inputParser;
p.addParamValue('FractalType', 'mandelbrot',...
    @(x)any(strcmpi(x,{'mandelbrot','m','julia','j'})));
p.addParamValue('JuliaConstant', 0.285+0.01i, @isnumeric);
p.addParamValue('MaxIter', 256, @(x)isnumeric(x)&&isreal(x)&&x>0);
p.addParamValue('EscapeRadius', 2, @(x)isnumeric(x)&&isreal(x)&&x>0);
p.addParamValue('Width', 400, @(x)isnumeric(x)&&isreal(x)&&x>0);
p.addParamValue('Height', 400, @(x)isnumeric(x)&&isreal(x)&&x>0);
p.addParamValue('Subfunction', 'c',...
    @(x)any(strcmpi(x,{'c','matlab','m'})));
p.addParamValue('XLim', [-2,2],...
    @(x)isnumeric(x)&&isreal(x)&&numel(x)>=2&&(x(2)-x(1))>0);
p.addParamValue('YLim', [-2,2],...
    @(x)isnumeric(x)&&isreal(x)&&numel(x)>=2&&(x(2)-x(1))>0);
p.addParamValue('ImageCacheSize', 8, @(x)isnumeric(x)&&isreal(x)&&x>2);
p.addParamValue('NonDivergingSetColor', [0,0,0],...
    @(x)isnumeric(x)&&isreal(x)&&all(x>=0)&&all(x<=1)&&size(x,1)==1);
p.addParamValue('Colormap', @hsv, @(x)isa(x,'function_handle'));
p.addParamValue('NumColor', 255, @(x)isnumeric(x)&&isreal(x)&&x>0);
p.parse(varargin{:});

% Colormap, windows has size constraint of 256
% Keep 1 slot for mandelbrot set points hence the 255
cm = [p.Results.NonDivergingSetColor(1:3); colormap(...
      p.Results.Colormap(min(p.Results.MaxIter,p.Results.NumColor)))];

% Determine if using C version or Matlab version for generating the fractal
if any(strcmpi(p.Results.Subfunction,{'c'}))
    if any(strcmpi(p.Results.FractalType,{'j','julia'}))
        if exist('makejulia_c','file') == 3
            hmakefractal = ...
                @(varargin)makejulia_c(varargin{:},p.Results.JuliaConstant);
        else
            warning('main:unknownfunction',...
                ['Must compile "makejulia_c.c" before using it. ',...
                'Reverted to matlab subfunction.']);
            hmakefractal = ...
                @(varargin)makejulia(varargin{:},p.Results.JuliaConstant);
        end
    else
        if exist('makemandel_c','file') == 3
            hmakefractal = @makemandel_c;
        else
            warning('main:unknownfunction',...
                ['Must compile "makemandel_c.c" before using it. ',...
                'Reverted to matlab subfunction.']);
            hmakefractal = @makemandel;
        end
    end
else
    if any(strcmpi(p.Results.FractalType,{'j','julia'}))
        hmakefractal = ...
            @(varargin)makejulia(varargin{:},p.Results.JuliaConstant);
    else
        hmakefractal = @makemandel;
    end
end

% Initialize figure
% Two axes, one for mandelbrot image, the other with initial limits to
% allow zooming and panning outside the image boundary
hfig = figure(1);
hax = zeros(2,1);
for k = 1:2
    hax(k) = axes(...
        'Parent', hfig,...
        'XLim', p.Results.XLim,...
        'YLim', p.Results.YLim,...
        'DataAspectRatio', [1 1 1],...
        'DataAspectRatioMode', 'manual',...
        'Box','on',...
        'XTick',[],...
        'YTick',[]);
end
set(hax(1),'Color',cm(2,:)); % bg color to 1 iteration for escape
set(hax(2),'Color','none');  % no color to these axes
xlabel(hax(1),[num2str(mean(p.Results.XLim),'%0.12f'),' \pm ',...
    num2str((p.Results.XLim(2)-p.Results.XLim(1))/2,'%0.12f')]);
ylabel(hax(1),[num2str(mean(p.Results.YLim),'%0.12f'),' \pm ',...
    num2str((p.Results.YLim(2)-p.Results.YLim(1))/2,'%0.12f')]);

% Display the image in hax(1) and let hax(2) be on top
% Have 'cache' versions of the image for better panning with himage(1)
% corresponding to the topmost image and himage(cache) is the background,
% i.e., the initial image. The other himage are previous generated image.
set(hfig,'CurrentAxes',hax(1));
[x,y] = meshgrid(...
    linspace(p.Results.XLim(1),p.Results.XLim(2),p.Results.Width),...
    linspace(p.Results.YLim(1),p.Results.YLim(2),p.Results.Height));
mu = hmakefractal(x,y,p.Results.EscapeRadius,p.Results.MaxIter);
CData = ones(size(mu));
CData(mu~=p.Results.MaxIter) = ...
    mod(mu(mu~=p.Results.MaxIter)-1,p.Results.NumColor)+2;
himage = zeros(p.Results.ImageCacheSize,1);
for k = p.Results.ImageCacheSize:-1:1
    if (k~=p.Results.ImageCacheSize)
        set(hax(1),'NextPlot','add');
    end
    himage(k) = image(...
        'XData',p.Results.XLim,...
        'YData',p.Results.YLim,...
        'CData',CData);
end
colormap(cm);

% Link the axes for automatic XLim and YLim resizing
linkaxes(hax,'xy');

% Anonymous function for generating and redrawing the fractal with
% the new boundaries
redraw = @redrawmandel;
    function redrawmandel
        % Get new plot limits from the axes
        XLim = get(hax(1),'XLim');
        YLim = get(hax(1),'YLim');
        % Update labels
        xlabel(hax(1),[num2str(mean(XLim),'%0.12f'),' \pm ',...
            num2str((XLim(2)-XLim(1))/2,'%0.12f')]);
        ylabel(hax(1),[num2str(mean(YLim),'%0.12f'),' \pm ',...
            num2str((YLim(2)-YLim(1))/2,'%0.12f')]);
        % Generate new image
        [x,y] = meshgrid(...
            linspace(XLim(1),XLim(2),p.Results.Height),...
            linspace(YLim(1),YLim(2),p.Results.Width));
        mu = hmakefractal(x,y,p.Results.EscapeRadius,p.Results.MaxIter);
        CData = ones(size(mu));
        CData(mu~=p.Results.MaxIter) = ...
            mod(mu(mu~=p.Results.MaxIter)-1,p.Results.NumColor)+2;
        % Update image cache
        for m = p.Results.ImageCacheSize-1:-1:2
            set(himage(m),...
                'XData',get(himage(m-1),'XData'),...
                'YData',get(himage(m-1),'YData'),...
                'CData',get(himage(m-1),'CData'));
        end
        % Update topmost image
        set(himage(1),...
            'XData',XLim,...
            'YData',YLim,...
            'CData',CData);
    end

% Set the zoom and pan callback functions
set(zoom,'ActionPostCallback',{@postcallback,redraw});
set(pan,'ActionPostCallback',{@postcallback,redraw});

% Eliminate rotate3d
set(rotate3d,'ButtonDownFilter',@rotatebuttondown);
end

function postcallback(obj,evd,redraw)
% Callback function for zoom and pan
redraw();
end

function res = rotatebuttondown(obj,evd)
% Return true to inhibit rotate3d
res = true(1);
end

function mu = makemandel(x,y,escape,maxiter)
% Generate mandelbrot set for xy coordinates given a certain escape radius
% and a maximum number of iterations
c = x(:)+1i*y(:);     % Put the coordinates in complex form
z = zeros(size(c));   % Initialize z_n
mu = ones(size(x));   % Initialize mu, the number of iterations before
                      % divergence
z_index = 1:numel(z); % Index of non-divergent z_n

% Loop until max iter reached
% Assuming maxiter will always be reached so no need to check for empty
% z_index (this assumes that there is nothing to see otherwise)
% Content of loop is taken from 'Mandelbrot set vectorized' by Lucio Cetto
% from Matlab Central
for iter = 1:maxiter
    z(z_index) = z(z_index).^2 + c(z_index);
    z_index = z_index(abs(z(z_index))<escape);
    mu(z_index) = iter;
end
end

function mu = makejulia(x,y,escape,maxiter,c)
% Generate monic quadratic julia set for xy coordinates given a certain
% escape radius and a maximum number of iterations
z = x(:)+1i*y(:);     % Put the coordinates in complex form
mu = ones(size(x));   % Initialize mu, the number of iterations before
                      % divergence
z_index = 1:numel(z); % Index of non-divergent z_n

% Loop until max iter reached
% Assuming maxiter will always be reached so no need to check for empty
% z_index (this assumes that there is nothing to see otherwise)
% Content of loop is derived from 'Mandelbrot set vectorized' by Lucio
% Cetto from Matlab Central
for iter = 1:maxiter
    z(z_index) = z(z_index).^2 + c;
    z_index = z_index(abs(z(z_index))<escape);
    mu(z_index) = iter;
end
end


text:
% Figure 5.40 Feedback Control of Dynamic Systems, 5e
%                   Franklin, Powell, Emami
%
% attitude hold auto pilot for Dakota
% root locus with integral control

clf
Wsp=6.3;  % short period
Zsp=.4;
Wph=.24;  % phugoid mode
Zph=.06;
Kg=160;    % gain so that 1" of control input produces 5 deg

NUMsp=[1 Zsp*Wsp];
NUMph=[1 3*Wph];
DENsp=[1 2*Zsp*Wsp Wsp^2];
DENph=[1 2*Zph*Wph Wph^2];

numG=Kg*conv(NUMsp,NUMph);
denG=conv(DENsp,DENph);
zeroG=roots(numG);

% lead compensation

numD=[1 3];
denD=[1 20];

num=conv(numG,numD);
den=conv(denG,denD);
zeroD=roots(numD);

K=1.5;
[r2]=rlocus(num,den,K);

% now add integral control

polesI=[r2';0];  % add integral control pole
denI=poly(polesI);

axis('square')
subplot(121)
axis([-30 0 -15 15])
plot(polesI,'x'),grid
hold on
plot(zeroD,0,'o')
plot(zeroG(1),0,'o',zeroG(2),0,'o')

ki0=.01:.05:.51;
ki1=1:.5:10;
ki2=12:5:102;
ki3=120:20:1000;
ki=[ki0 ki1 ki2 ki3 5000];
[R]=rlocus(K*num,denI,ki);
plot(R,'-')
axis([-25 2 -15 15]);
title ('Fig. 5.40 Root Locus vs. K_I for autopilot')
xlabel('Re(s)')
ylabel('Im(s)')

KI=.15;
[R1]=rlocus(K*num,denI,KI);
damp(R1')
plot(R1,'*')
hold off

% now zoom in to origin

subplot(122)
plot(polesI,'x'),grid
axis([-4 0 -2 2])
hold on
xlabel('Re(s)')
ylabel('Im(s)')
plot(zeroD,0,'o')
plot(zeroG(1),0,'o',zeroG(2),0,'o')
plot(R,'-')
plot(R1,'*')

hold off
axis('normal')

subplot(111)


text:
function pass = test_roots( )
% Test the chebfun2/roots command.

tol = 1e-12;  % rank 1 curves
tol2 = 1e-8;  % rank >1 curves, disjoint
tol3 = 1e-4;  % rank >1 curves, noncircular
tol4 = 1e-2;  % intersecting curves, not properly disentangled

% First we check separable functions, i.e., rank 1:

f = chebfun2( @(x,y) 1/2-y);
r = roots( f );
exact = chebfun( @(x) x + 1i/2 );
pass(1) = ( norm( r - exact ) < tol );

f = chebfun2( @(x,y) 1/2-x);
r = roots( f );
exact = chebfun( @(x) 1i*x + 1/2 );
pass(2) = ( norm( r - exact ) < tol );

f = chebfun2( @(x,y) x.*y);
r = roots( f );
exact = chebfun( @(x) [x, 1i*x]  );
pass(3) = ( norm( r - exact ) < tol );

f = chebfun2( @(x,y) cos(5*pi*x));
s = roots( squeeze( f ) );
r = roots( f );
exact = chebfun;
for j = 1:numel(s)
    exact = [ exact chebfun(@(x) s(j) + 1i*x )];
end
pass(4) = ( norm( r - exact ) < tol );

f = chebfun2( @(x,y) x.*(y-1/2), [-2 2 -3 5]);
r = roots( f );
exact = [chebfun( @(x) 2*x + 1i/2 ) chebfun(@(x) 1i*(4*(x+1)-3) )];
pass(5) = ( norm( r - exact ) < tol );

% Next we move to functions of rank > 1

f = chebfun2( @(x,y) x.^2 + y.^2 - 1/4 );
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - pi;
pass(6) = abs(err) < tol2;
%disp([int2str(pass(6)) '  circle  ' int2str(length(c))])

area = abs(sum(real(c).*diff(imag(c))));
err = area - pi/4;
pass(7) = abs(err) < tol2;
%disp([int2str(pass(7)) '  area of circle  ' int2str(length(c))])

f = chebfun2( @(x,y) x.^2 + (10*y).^2 - 1/4 );
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - 2.031987090050447;
pass(8) = abs(err) < tol3;
%disp([int2str(pass(8)) '  ellipse  ' int2str(length(c))])

area = abs(sum(real(c).*diff(imag(c))));
err = area - pi/40;
pass(9) = abs(err) < tol3;
%disp([int2str(pass(9)) '  area of ellipse  ' int2str(length(c))])

f = chebfun2( @(x,y) (x-1).^2 + y.^2 - 1/4 );
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - pi/2;
pass(10) = abs(err) < tol2;
%disp([int2str(pass(10)) '  semicircle at boundary  ' int2str(length(c))])

f = chebfun2( @(x,y) (x-1).^2 + (y+1).^2 - 1/4 );
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - pi/4;
pass(11) = abs(err) < tol2;
%disp([int2str(pass(11)) '  quarter-circle at corner  ' int2str(length(c))])

r = 1/pi;  % this radius gives circumference 2
x = chebfun2(@(x,y) x);
y = chebfun2(@(x,y) y);

%% One circle
f = x.^2 + y.^2 - r^2;
c = roots(f); arclength = sum(abs(diff(c)));
err = arclength - 2;
pass(12) = abs(err) < tol2;
%disp([int2str(pass(12)) '  another circle  ' int2str(length(c))])

%% Two circles
f = x.^2 + y.^2 - r^2;
g = (x-.6).^2 + (y-.5).^2 - r^2;
c = roots(f.*g); arclength = sum(abs(diff(c)));
err = arclength - [2 2];
pass(13) = norm(err,inf) < tol2;
%disp([int2str(pass(13)) '  two circles  ' int2str(length(c))])

%% Intersecting circles
f = x.^2 + y.^2 - r^2;
g = (x-.2).^2 + (y-.5).^2 - r^2;
c = roots(f.*g); arclength = sum(abs(diff(c)));
err = sum(arclength) - 4;
pass(14) = abs(err) < tol4;
%disp([int2str(pass(14)) '  intersecting circles  ' int2str(length(c))])

%% A circle and a line
f = x.^2 + y.^2 - r^2;
w = exp(x)-1.2;
c = roots(f.*w); arclength = sum(abs(diff(c)));
err = sum(arclength) - 4;
pass(15) = abs(err) < tol4;
%disp([int2str(pass(15)) '  intersecting circle and line  ' int2str(length(c))])

%% Concentric circles (well separated)
f = chebfun2(@(x,y) (x.^2 + y.^2 - .25).*(x.^2+y.^2-0.300));
c = roots(f); arclength = sum(abs(diff(c)));
err = arclength - pi*[sqrt(.300/.250) 1];
pass(16) = norm(err,inf) < tol2;
%disp([int2str(pass(16)) '  concentric circles  ' int2str(length(c))])

%% Concentric circles (less well separated)
f = chebfun2(@(x,y) (x.^2 + y.^2 - .25).*(x.^2+y.^2-0.260));
c = roots(f); arclength = sum(abs(diff(c)));
err = arclength - pi*[sqrt(.260/.250) 1];
pass(17) = norm(err,inf) < tol2;
%disp([int2str(pass(17)) '  close concentric circles  ' int2str(length(c))])

f = chebfun2( @(x,y) (x.^2+y.^2-1/4).*((x-1.1).^2+(y-3).^2-1/4),[-1.3 1.1 -.9 3]);
c = roots(f);
arclength = sum(abs(diff(c)));
err = sum(arclength) - 1.25*pi;
pass(18) = abs(err) < tol2;
%disp([int2str(pass(18)) '  different domain  ' int2str(length(c))])

f = chebfun2( @(x,y) 1e100*(x.^2 + y.^2 - 1/4) );
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - pi;
pass(19) = abs(err) < tol2;
%disp([int2str(pass(19)) '  vertical scaling 1  ' int2str(length(c))])

f = chebfun2( @(x,y) 1e100*(x.^2 + y.^2 - 1/4) );
area = abs(sum(real(c).*diff(imag(c))));
err = area - pi/4;
pass(20) = abs(err) < tol2;
%disp([int2str(pass(20)) '  vertical scaling 2  ' int2str(length(c))])

f = chebfun2( @(x,y) y-sin(x),[-pi pi,-1.1 1.3]);
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - 7.640395578055424035;
pass(21) = abs(err) < tol2;

f = chebfun2( @(x,y) x-tanh(y));
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - 2.53182746833076;
pass(22) = abs(err) < tol2;

f = chebfun2( @(x,y) y-x.^2);
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - 2.957885715089195;
pass(23) = abs(err) < tol2;

f = chebfun2( @(x,y) y.^2-x.^2-.1,[-1 1 -1.5 1.5]);  % hyperbola
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - [1 1]*2.51829399795912;
pass(24) = norm(err,inf) < tol2;

f = chebfun2( @(x,y) y.^2-x.^2);  % cross
c = roots(f);
arclength = sum(abs(diff(c)));
err = arclength - [1 1]*sqrt(8);
pass(25) = norm(err,inf) < tol4;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% connect to Arduino
try
    a = arduino;
    disp(a)
catch
    warning('Unable to connect, user input required')
    disp('For Windows:')
    disp('  Open device manager, select "Ports (COM & LPT)"')
    disp('  Look for COM port of Arduino such as COM4')
    disp('For MacOS:')
    disp('  Open terminal and type: ls /dev/*.')
    disp('  Search for /dev/tty.usbmodem* or /dev/tty.usbserial*. The port number is *.')
    disp('For Linux')
    disp('  Open terminal and type: ls /dev/tty*')
    disp('  Search for /dev/ttyUSB* or /dev/ttyACM*. The port number is *.')
    disp('')
    com_port = input('Specify port (e.g. COM4 for Windows or /dev/ttyUSB0 for Linux): ','s');
    a = arduino(com_port,'Uno');
    disp(a)
end

% voltage read functions
v1 = @() readVoltage(a, 'A0');
v2 = @() readVoltage(a, 'A2');

% temperature calculations as a function of voltage for TMP36
TC = @(V) (V - 0.5)*100.0;          % Celsius
TK = @(V) TC(V) + 273.15;           % Kelvin
TF = @(V) TK(V) * 9.0/5.0 - 459.67; % Fahrenhiet

% temperature read functions
T1C = @() TC(v1());
T2C = @() TC(v2());

% LED function (0 <= level <= 1)
led = @(level) writePWMDutyCycle(a,'D9',max(0,min(1,level)));  % ON

% heater output (0 <= heater <= 100)
% limit to 0-0.9 (0-100%)
h1 = @(level) writePWMDutyCycle(a,'D3',max(0,min(100,level))*0.9/100);
% limit to 0-0.5 (0-100%)
h2 = @(level) writePWMDutyCycle(a,'D5',max(0,min(100,level))*0.5/100);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function M=Find_LCM(x)

%%%%%%%%%% This fucntion can find least common multiple between several
%%%%%%%%%% numbers.

temp=x(1);

for i=2:length(x)
     temp=lcm(temp,x(i));
end

M=temp;


text:
% block_fatrix_test.m
% Test the block_fatrix object

if 1 || ~isvar('A5'), printm 'setup'
	rng(0)
%	A1 = Gsparse(sparse(rand(10,20)));
	A1 = rand(10,20);
	A2 = magic(20);
	A3 = rand(10,30);
	A4 = Gdft('mask', true(5,4), 'samp', rand(5,4) > 0.5); % for gram
	A5 = rand(size(A4));

	Ac = block_fatrix({A1, A2}, 'type', 'col');
	Ad = block_fatrix({A1, A4}); % diag
	Ak = block_fatrix({A1}, 'type', 'kron', 'Mkron', 2);
	Ar = block_fatrix({A1, A3}, 'type', 'row');
	As = block_fatrix({A4, A5}, 'type', 'sum');
end

if 1, printm 'basic Fatrix tests'
	tester = @(A, mask, name) ...
		Fatrix_test_basic(A, mask, 'complex', 1, 'name', name)
	tester(Ac, true(20,1), 'A:col')
	tester(Ad, true(40,1), 'A:diag')
	tester(Ak, true(40,1), 'A:kron')
	tester(Ar, true(50,1), 'A:row')
	tester(As, true(20,1), 'A:sum')
end

if 1, printm 'adjoint tests'
	tester = @(A) test_adjoint(A, 'complex', 1);
	tester(Ac);
	tester(Ad);
	tester(Ak);
	tester(Ar);
	tester(As);
end

if 1 % test with a "tomo" object
	ig = image_geom('nx', 8, 'ny', 7, 'dx', 1);
	ig.mask = ig.circ > 0;
	sg0 = sino_geom('par', 'nb', 9, 'na', 10, 'dr', 1);
	sg1 = sg0; sg1.na = 6; sg1.orbit = sg0.orbit / sg0.na * sg1.na;
	sg2 = sg1; sg2.na = sg0.na - sg1.na; sg2.orbit_start = sg1.orbit;
		sg2.orbit = sg0.orbit - sg1.orbit;
	tmp = {'square/strip', 'Ltab', [1000], 'strip_width', sg0.dr};
	At0 = Gtomo2_table(sg0, ig, tmp);
	At1 = Gtomo2_table(sg1, ig, tmp);
	At2 = Gtomo2_table(sg2, ig, tmp);
	Atb = block_fatrix({At1, At2}, 'type', 'col', 'tomo', true);

	if 1
		xx = ig.circ;
		yt = cat(2, At1 * xx, At2 * xx);
		yb = Atb * xx;
		jf_equal(yb, yt)
		jf_equal(yb, At0 * xx)
	end

	if 1
		yy = Atb * ig.circ;
		yc = mat2cell(yy, [sg1.nb], [sg1.na sg2.na]);
		xx = At1' * yc{1} + At2' * yc{2};
		x0 = At0' * yy;
		equivs(x0, xx)
		xb = Atb' * yy;
		jf_equal(xb, xx)
	end

%	equivs(full(At0), full(Atb))
	tester_tomo2(Atb, ig.mask, 'A2', At0, 'halt', 0, 'nblock', 0)
end

if 1 % test 'col'
	x = [1:ncol(A1)]';
	y1 = A1 * x;
	y2 = A2 * x;
	yy = Ac * x;
%	printm('col forw error %g', max_percent_diff([y1; y2], yy))
	jf_equal([y1; y2], yy)

	x1 = A1' * y1;
	x2 = A2' * y2;
	xx = Ac' * [y1; y2];
%	printm('col back error %g', max_percent_diff(x1+x2, xx))
	jf_equal(x1+x2, xx)
end

if 1 % test 'diag'
	x1 = [1:ncol(A1)]';
	x2 = [1:ncol(A4)]';
	x = [x1; x2];
	y1 = A1 * x1;
	y2 = A4 * x2;
	yy = Ad * x;
%	printm('diag forw error %g', max_percent_diff([y1; y2], yy))
	jf_equal([y1; y2], yy)

	x1 = A1' * y1;
	x2 = A4' * y2;
	xx = Ad' * yy;
%	printm('diag back error %g', max_percent_diff([x1; x2], xx))
	jf_equal([x1; x2], xx)

	Td = build_gram(Ad, [], 0);
	y1 = Td * xx;
	y2 = [A1' * A1 * x1; A4' * (A4 * x2)];
%	printm('diag gram error %g%%', max_percent_diff(y1, y2))
	jf_equal(y1, y2)
end

if 1 % test 'kron'
	xx = rand(ncol(A1), Ak.arg.Mkron);
	clear xt yt
	for ic=1:Ak.arg.Mkron
		yt(:,ic) = A1 * xx(:,ic);
	end
	yy = Ak * xx(:);
%	printm('kron forw error %g', max_percent_diff(yt(:), yy))
	jf_equal(yt(:), yy)

	yy = reshapee(yy, [], Ak.arg.Mkron);
	for ic=1:Ak.arg.Mkron
		xt(:,ic) = A1' * yy(:,ic);
	end
	xx = Ak' * yy(:);
%	printm('kron back error %g', max_percent_diff(xt(:), xx))
	jf_equal(xt(:), xx)
end

if 1 % test 'row'
	x1 = [1:ncol(A1)]';
	x2 = [1:ncol(A3)]';
	y1 = A1 * x1;
	y2 = A3 * x2;
	yy = Ar * [x1; x2];
%	printm('row forw error %g', max_percent_diff(y1+y2, yy))
	equivs(y1+y2, yy)

	x1 = A1' * yy;
	x2 = A3' * yy;
	xx = Ar' * yy;
%	printm('row back error %g', max_percent_diff([x1; x2], xx))
	equivs([x1; x2], xx)
end

if 1 % test 'sum'
	x = [1:ncol(A4)]';
	y1 = A4 * x;
	y2 = A5 * x;
	yy = As * x;
%	printm('sum forw error %g', max_percent_diff(y1+y2, yy))
	jf_equal(y1+y2, yy)
	x1 = A4' * yy;
	x2 = A5' * yy;
	xx = As' * yy;
%	printm('sum back error %g', max_percent_diff(x1+x2, xx))
	jf_equal(x1+x2, xx)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [km] = mi2km(mi)
% Convert length from miles to kilometers.
% Chad A. Greene 2012
km = mi*1.609344;


text:
function [sens] = ft_read_sens(filename, varargin)

% FT_READ_SENS read sensor positions from various manufacturer specific files. See
% further down for the list of file types that are supported.
%
% Use as
%   elec = ft_read_sens(filename, 'senstype', 'eeg', ...)  % for EEG electrodes
%   grad = ft_read_sens(filename, 'senstype', 'meg', ...)  % for MEG gradiometers
%   opto = ft_read_sens(filename, 'senstype', 'nirs', ...) % for NIRS optodes
%
% Additional options should be specified in key-value pairs and can be
%   'fileformat'     = string, see the list of supported file formats (the default is determined automatically)
%   'senstype'       = string, can be 'eeg', 'meg' or 'nirs', specifies which type of sensors to read from the file (default = 'eeg')
%   'coordsys'       = string, 'head' or 'dewar' (default = 'head')
%   'coilaccuracy'   = scalar, can be empty or a number (0, 1 or 2) to specify the accuracy (default = [])
%   'readbids'       = string, 'yes', no', or 'ifmakessense', whether to read information from the BIDS sidecar files (default = 'ifmakessense')
%
% The electrode, gradiometer and optode structures are defined in more detail
% in FT_DATATYPE_SENS.
%
% Files from the following acquisition systems and analysis platforms file formats
% are supported.
%   asa_elc besa_elp besa_pos besa_sfp yokogawa_ave yokogawa_con yokogawa_raw 4d
%   4d_pdf 4d_m4d 4d_xyz ctf_ds ctf_res4 itab_raw itab_mhd netmeg neuromag_fif
%   neuromag_mne neuromag_mne_elec neuromag_mne_grad polhemus_fil polhemus_pos
%   zebris_sfp spmeeg_mat eeglab_set localite_pos artinis_oxy3 artinis_oxy4
%   artinis_oxy5 artinis_oxyproj yorkinstruments_hdf5 matlab
%
% See also FT_READ_HEADER, FT_DATATYPE_SENS, FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD,

% Copyright (C) 2005-2022 Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% optionally get the data from the URL and make a temporary local copy
filename = fetch_url(filename);

% get the options
fileformat     = ft_getopt(varargin, 'fileformat', ft_filetype(filename));
senstype       = ft_getopt(varargin, 'senstype');         % can be eeg/meg/nirs, default is automatic and eeg when both meg+eeg are present
coordsys       = ft_getopt(varargin, 'coordsys', 'head'); % this is used for ctf and neuromag_mne, it can be head or dewar
coilaccuracy   = ft_getopt(varargin, 'coilaccuracy');     % empty, or a number between 0 to 2
readbids       = ft_getopt(varargin, 'readbids', 'ifmakessense');

realtime = any(strcmp(fileformat, {'fcdc_buffer', 'ctf_shm', 'fcdc_mysql'}));

if realtime
  % skip the rest of the initial checks to increase the speed for realtime operation
else
  % test whether the file exists
  if ~exist(filename, 'file')
    ft_error('file ''%s'' does not exist', filename);
  end
end

% start with an empty electrode, gradiometer or optode definition
sens = [];

% deal with data that is organized according to BIDS
if strcmp(readbids, 'yes') || strcmp(readbids, 'ifmakessense')
  [p, f, x] = fileparts(filename);
  % check whether it a BIDS dataset
  isbids = startsWith(f, 'sub-');
  if isbids
    tsvfile = bids_sidecar(filename, 'electrodes');
    if ~isempty(tsvfile) && (isempty(senstype) || strcmp(senstype, 'eeg'))
      % read the electrodes.tsv file
      electrodes_tsv = ft_read_tsv(tsvfile);
      sens         = [];
      sens.label   = electrodes_tsv.name;
      sens.elecpos = [electrodes_tsv.x electrodes_tsv.y electrodes_tsv.z];

      % also read the electrodes.json file
      [p, f] = fileparts(tsvfile);
      jsonfile = fullfile(p, [f '.json']);
      if exist(jsonfile, 'file')
        electrodes_json = ft_read_json(jsonfile);
        ft_warning('the content of the electrodes.json is not used')
        % FIXME do something with the content
      end

      % also read the coordsystem.json file
      coordsysfile = bids_sidecar(filename, 'coordsystem');
      if exist(coordsysfile, 'file')
        coordsys_json = ft_read_json(coordsysfile);
        ft_warning('the content of the coordsystem.json is not used')
        % FIXME do something with the content
      end
    end
  end
end

switch fileformat
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % gradiometer information is always stored in the header of the MEG dataset
  % hence we use the standard fieldtrip/fileio ft_read_header function
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case {'ctf_ds', 'ctf_res4', 'ctf_old', 'neuromag_fif', 'neuromag_mne', '4d', '4d_pdf', '4d_m4d', '4d_xyz', 'yokogawa_ave', 'yokogawa_con', 'yokogawa_raw', 'ricoh_ave', 'ricoh_con', 'itab_raw' 'itab_mhd', 'netmeg'}
    hdr = ft_read_header(filename, 'headerformat', fileformat, 'coordsys', coordsys, 'coilaccuracy', coilaccuracy, 'readbids', readbids);
    % sometimes there can also be electrode position information in the header
    if isfield(hdr, 'elec') && isfield(hdr, 'grad')
      if isempty(senstype)
        % set the default
        ft_warning('both electrode and gradiometer information is present, returning the electrode information by default');
        senstype = 'eeg';
      end
      switch lower(senstype)
        case 'eeg'
          sens = hdr.elec;
        case 'meg'
          sens = hdr.grad;
        otherwise
          ft_error('incorrect specification of senstype');
      end
    elseif isfield(hdr, 'grad')
      sens = hdr.grad;
    elseif isfield(hdr, 'elec')
      sens = hdr.elec;
    else
      ft_error('there is no electrode nor gradiometer information present in the header');
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % optode information is mostly stored in the header of the NIRS dataset
    % hence we use the standard fieldtrip/fileio ft_read_header function
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case {'homer_nirs', 'snirf', 'artinis_oxy3', 'artinis_oxy4', 'artinis_oxy5', 'artinis_oxyproj', 'nirx_wl1', 'nirx_wl2', 'nirx_tpl'}
    hdr = ft_read_header(filename, 'headerformat', fileformat, 'coordsys', coordsys, 'readbids', readbids);
    if isfield(hdr, 'opto')
      sens = hdr.opto;
    else
      ft_error('there is no optode information present in the header');
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % read the content from various files that contain EEG electrode positions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case 'asa_elc'
    sens = read_asa_elc(filename);

  case 'polhemus_pos'
    sens = read_polhemus_pos(filename);

  case 'besa_elp'
    fid = fopen_or_error(filename);
    % these files seem to come in different formats with 3, 4 or 5 columns
    % see http://wiki.besa.de/index.php?title=Channel_Definition_File_Formats
    % read the first two lines to determine the number of columns
    columns1 = length(strsplit(strtrim(fgetl(fid))));
    columns2 = length(strsplit(strtrim(fgetl(fid))));
    fseek(fid, 0, 'bof');
    if columns1==1 && columns2>1
      % EEGLAB includes ELP files that start with a line with the number of electrodes
      % skip the first line
      columns1 = columns2;
      fgetl(fid);
    end
    switch columns1
      case 3
        % 3-column: label, azimuth, elevation
        tmp = textscan(fid, '%s%f%f');
        type      = repmat({'EEG'}, size(tmp{1}));
        label     = tmp{1};
        theta     = tmp{2};
        phi       = tmp{3};
        radius    = repmat(85, size(tmp{1}));
      ft_warning('assuming a head radius of 85 mm');
      case 4
        % 4-column: type, label, azimuth, elevation
        tmp = textscan(fid, '%s%s%f%f');
        type      = tmp{1};
        label     = tmp{2};
        theta     = tmp{3};
        phi       = tmp{4};
        radius    = repmat(85, size(label));
      ft_warning('assuming a head radius of 85 mm');
      case 5
        % 5-column: type, label, azimuth, elevation, radius
        tmp = textscan(fid, '%s%s%f%f%f');
        type      = tmp{1};
        label     = tmp{2};
        theta     = tmp{3};
        phi       = tmp{4};
        radius    = tmp{5};
      otherwise
        ft_error('unsupported file format for .elp');
    end
    fclose(fid);

    radians = @(degree) degree*pi/180;
    x = radius .* cos(radians(phi))   .* sin(radians(theta));
    y = radius .* sin(radians(theta)) .* sin(radians(phi));
    z = radius .* cos(radians(theta));
    sel = strcmpi(type, 'EEG') | strcmpi(type, 'SCP') | strcmpi(type, 'REF');
    sens.elecpos = [x(sel) y(sel) z(sel)];
    sens.chanpos = [x(sel) y(sel) z(sel)];
    sens.label   = label(sel);
    sens.unit    = 'mm';
    sel = strcmpi(type, 'FID');
    if any(sel)
      sens.fid.pos    = [x(sel) y(sel) z(sel)];
      sens.fid.label  = label(sel);
    end

  case 'besa_pos'
    tmp = importdata(filename);
    if ~isnumeric(tmp)
      ft_error('unexpected file format for fileformat=besa_pos')
    end
    [nchan,nrow] = size(tmp);
    if nrow==3
      sens.pnt = tmp;
    elseif nrow==9
      pnt1 = tmp(:,1:3);  % bottom coil
      pnt2 = tmp(:,4:6);  % top coil
      ori  = tmp(:,7:9);  % orientation of bottom coil
      sens.pnt = [pnt1; pnt2];
      sens.ori = [ori; ori];
      sens.tra = [eye(nchan) -eye(nchan)];
    else
      ft_error('unexpected file format for fileformat=besa_pos')
    end
    [p, f, x] = fileparts(filename);
    elpfile = fullfile(p, [f '.elp']);
    elafile = fullfile(p, [f '.ela']);
    if exist(elpfile, 'file')
      ft_warning('reading channel labels from %s', elpfile);
      % read the channel names from the accompanying ELP file
      lbl = importdata(elpfile);
      sens.label = strrep(lbl.textdata(:,2) ,'''', '');
    elseif exist(elafile, 'file')
      ft_warning('reading channel labels from %s', elafile);
      % read the channel names from the accompanying ELA file
      lbl = importdata(elafile);
      lbl = strrep(lbl, 'MEG ', ''); % remove the channel type
      lbl = strrep(lbl, 'EEG ', ''); % remove the channel type
      sens.label = lbl;
    else
      % the file does not have channel labels in it
      ft_warning('creating fake channel names for besa_pos');
      for i=1:nchan
        sens.label{i} = sprintf('%03d', i);
      end
    end

  case 'besa_sfh'
    sfh = readBESAsfh(filename);
    sens.label   = sfh.SurfacePointsLabels(:);
    sens.elecpos = sfh.SurfacePointsCoordinates(:,1:3);
    sel = true(sfh.NrOfPoints, 1);
    for i=1:sfh.NrOfPoints
      tok = tokenize(sens.label{i}, '_');
      sens.label{i} = tok{2};
      sel(i) = ~strcmp(tok{1}, 'Fid');
    end
    sens.label   = sens.label(sel);
    sens.elecpos = sens.elecpos(sel,:);

  case 'besa_sfp'
    [lab, pos] = read_besa_sfp(filename);
    sens.label   = lab;
    sens.elecpos = pos;

  case 'bioimage_mgrid'
    sens = read_bioimage_mgrid(filename);

  case {'curry_dat', 'curry_cdt'}
    hdr = ft_read_header(filename);
    if ~isempty(hdr.orig.sensorpos)
      sens.elecpos = hdr.orig.sensorpos';
      sens.label   = hdr.label(1:size(sens.elecpos, 1));
    end

  case 'fcdc_buffer'
    % the online header should have a binary blob with the sensor information
    hdr = ft_read_header(filename, 'headerformat', fileformat);
    if isfield(hdr, 'orig') && isfield(hdr.orig, 'neuromag_header')
      hdr = decode_fif(hdr.orig.neuromag_header);
    end
    if isfield(hdr, 'orig') && isfield(hdr.orig, 'ctf_res4')
      hdr = decode_res4(hdr.orig.ctf_res4);
    end
    if isempty(senstype)
      % set the default
      ft_warning('both electrode and gradiometer information is present, returning the electrode information by default');
      senstype = 'eeg';
    end
    switch lower(senstype)
      case 'eeg'
        sens = hdr.elec;
      case 'meg'
        sens = hdr.grad;
      otherwise
        ft_error('incorrect specification of senstype');
    end

  case 'neuromag_mne_grad'
    % the file can contain both, force reading the gradiometer info
    % note that this functionality overlaps with senstype=eeg/meg
    hdr = ft_read_header(filename, 'headerformat', 'neuromag_mne', 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    sens = hdr.grad;

  case 'neuromag_mne_elec'
    % the file can contain both, force reading the electrode info
    % note that this functionality overlaps with senstype=eeg/meg
    hdr = ft_read_header(filename, 'headerformat', 'neuromag_mne', 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    sens = hdr.elec;

  case {'spmeeg_mat', 'eeglab_set'}
    % this is for EEG formats where electrode positions can be stored with the data
    hdr = ft_read_header(filename, 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    if isfield(hdr, 'grad')
      sens = hdr.grad;
    elseif isfield(hdr, 'elec')
      sens = hdr.elec;
    else
      ft_error('no electrodes or gradiometers found in the file')
    end

  case 'polhemus_fil'
    % these are created at the FIL in London with a polhemus tracker
    [sens.fid.pnt, sens.pnt, sens.fid.label] = read_polhemus_fil(filename, 0);
    % the file does not have channel labels in it
    ft_warning('no channel names in polhemus file, using numbers instead');
    for i=1:size(sens.pnt, 1)
      sens.label{i} = sprintf('%03d', i);
    end

  case 'matlab'
    % MATLAB files can contain all sensor arrays
    matfile = filename;   % this solves a problem with the MATLAB compiler v3
    var = whos('-file', filename);
    sel = intersect({'elec', 'grad', 'opto', 'sens', 'elc'}, {var(:).name});
    if numel(sel)==1
      % read the specific variable
      sens = loadvar(matfile, sel{1});
    else
      % read whatever variable is in the file, this will error if the file contains multiple variables
      sens = loadvar(matfile);
    end
    if ft_datatype(sens, 'layout')
      ft_error('"%s" contains a 2D layout for plotting, not 3D electrode positions', filename);
    end

  case 'zebris_sfp'
    % these are created by a Zebris tracker, at CRC in Liege at least.
    [sens.fid.pnt, sens.chanpos, sens.fid.label, sens.label] = read_zebris(filename, 0);
    % convert to columns
    sens.label = sens.label(:);
    sens.fid.label = sens.fid.label(:);

  case '4d_el_ascii'
    fid = fopen_or_error(filename, 'rt');
    c = textscan(fid, '%s%s%f%f%f');
    l = c{:,1}; % label
    s = c{:,2}; % status, it can be 'Collected' or empty
    x = c{:,3};
    y = c{:,4};
    z = c{:,5};
    % shift the columns with one where the status is not specified
    sel = isnan(z);
    z(sel) = y(sel);
    y(sel) = x(sel);
    x(sel) = str2double(s(sel));
    s(sel) = {''};
    fclose(fid);
    if false
      % return all positions, including the ones that do not correspond to
      % electrodes per see, such as the fiducials and localizer coils
      sens          = [];
      sens.label    = l;
      sens.elecpos  = [x y z];
    else
      % split the electrodes and fiducials
      % this is consistent with zebris_sfp and with the output of ft_read_headshape
      sens            = [];
      sens.label      = l(~sel);
      sens.elecpos    = [x(~sel) y(~sel) z(~sel)];
      sens.fid.label  = l(sel);
      sens.fid.pnt    = [x(sel) y(sel) z(sel)];
    end

  case {'localite_pos', 'localite_ins'}
    if ~usejava('jvm') % Using xml2struct requires java
      fid = fopen_or_error(filename);

      % Read marker-file and store contents in cells of strings
      tmp = textscan(fid,'%s');

      fclose(fid);

      % Search for cells that contain coordinates
      selx = strncmp('data0',tmp{1},5);
      sely = strncmp('data1',tmp{1},5);
      selz = strncmp('data2',tmp{1},5);
      sellab = strncmp('description',tmp{1},5);

      % Extract cells that contain coordinates
      xtemp  = tmp{1}(selx);
      ytemp  = tmp{1}(sely);
      ztemp  = tmp{1}(selz);
      labtemp = tmp{1}(sellab);

      % Determine which channels are set. In localite channels that are not set
      % automatically receive coordinates [0, 0, 0] and should therefore
      % be discarded.
      settemp = tmp{1}(strncmp('set',tmp{1},3));
      selset = strncmp('set="f',settemp,6);

      % Remove channels that are not set
      xtemp(selset) = [];
      ytemp(selset) = [];
      ztemp(selset) = [];
      labtemp(selset) = [];

      % Convert cells that contain coordinates from string to double
      x = [];
      y = [];
      z = [];
      lbl = [];

      for i=1:numel(xtemp)
        x(i,1) = str2double(xtemp{i}(8:end-1));
        y(i,1) = str2double(ytemp{i}(8:end-1));
        z(i,1) = str2double(ztemp{i}(8:end-3));
        lbl{i,1} = labtemp{i}(14:end-1);
      end

      % Create and fill sens structure
      sens = [];
      sens.elecpos = [x y z];
      sens.chanpos = sens.elecpos;
      sens.label = lbl;
    else
      tmp = xml2struct(filename);

      sens = [];

      % Loop through structure obtained from xml-file and store
      % coordinate information into sens structure of channels that are
      % set.
      for i=1:numel(tmp)
        if strcmp(tmp(i).Marker.set, 'true')
          sens.elecpos(i,1) = str2double(tmp(i).Marker.ColVec3D.data0);
          sens.elecpos(i,2) = str2double(tmp(i).Marker.ColVec3D.data1);
          sens.elecpos(i,3) = str2double(tmp(i).Marker.ColVec3D.data2);
          sens.label{i} = tmp(i).Marker.description;
        end
      end

      sens.chanpos = sens.elecpos;
    end

  case 'easycap_txt'
    % Read the file and store all contents in cells of strings
    fid = fopen_or_error(filename);
    tmp = textscan(fid,'%s%s%s%s');
    fclose(fid);

    sens = [];
    if all(cellfun(@isempty, tmp{4}))
      % it contains theta and phi
      sens.label = tmp{1}(2:end);
      theta = cellfun(@str2double, tmp{2}(2:end));
      phi   = cellfun(@str2double, tmp{3}(2:end));
      radians = @(degree) degree*pi/180;
      ft_warning('assuming a head radius of 85 mm');
      x = 85 * cos(radians(phi))   .* sin(radians(theta));
      y = 85 * sin(radians(theta)) .* sin(radians(phi));
      z = 85 * cos(radians(theta));
      sens.unit = 'mm';
      sens.elecpos = [x y z];
      sens.chanpos = [x y z];
    else
      % it contains X, Y, Z
      sens.label = tmp{1}(2:end);
      x = cellfun(@str2double, tmp{2}(2:end));
      y = cellfun(@str2double, tmp{3}(2:end));
      z = cellfun(@str2double, tmp{4}(2:end));
      sens.elecpos = [x y z];
      sens.chanpos = [x y z];
    end

  case 'neuromag_iso'
    ft_hastoolbox('mne', 1);
    FIFF = fiff_define_constants();
    [fid, tree, dir] = fiff_open(filename);
    isotrak = fiff_dir_tree_find(tree, FIFF.FIFFB_ISOTRAK);
    sel = find([isotrak.dir.kind]==FIFF.FIFF_DIG_POINT);
    sens = [];
    sens.elecpos = nan(numel(sel),3);
    sens.chanpos = nan(numel(sel),3);
    coordsys     = nan(numel(sel),1);
    for i=sel
      tag = fiff_read_tag(fid,isotrak.dir(i).pos);
      sens.elecpos(i,:) = tag.data.r;
      sens.chanpos(i,:) = tag.data.r;
      coordsys(i)       = tag.data.coord_frame;
    end
    fclose(fid);

    if all(coordsys==FIFF.FIFFV_COORD_DEVICE)
      sens.coordsys = 'device';
    elseif all(coordsys==FIFF.FIFFV_COORD_ISOTRAK)
      sens.coordsys = 'isotrak';
    elseif all(coordsys==FIFF.FIFFV_COORD_HPI)
      sens.coordsys = 'hpi';
    elseif all(coordsys==FIFF.FIFFV_COORD_HEAD)
      sens.coordsys = 'head';
    else
      sens.coordsys = 'unknown';
    end

    ft_warning('creating fake channel names for neuromag_iso');
    for i=1:size(sens.chanpos,1)
      sens.label{i} = sprintf('%d', i);
    end

  case 'neuromag_cal'
    dat = cell(1,14);
    [dat{:}] = textread(filename, '%s%f%f%f%f%f%f%f%f%f%f%f%f%f');
    label = dat{1};
    x = dat{2};
    y = dat{3};
    z = dat{4};
    rot = cat(2, dat{5:13});
    cal = dat{14}; % ??
    % construct the sensor structucture
    % it seems that the channel positions are expressed in dewar cordinates
    % it would be possible to use coil_def.dat to construct the coil positions
    sens.label = label;
    sens.chanpos = [x y z];

  case '3dslicer_fscv'
    csvData = readtable(filename, 'FileType', 'text');
    sens.label = csvData.label;
    sens.elecpos = [csvData.x,csvData.y,csvData.z];

  case 'brainsight_txt'
    ws = warning('off', 'MATLAB:table:ModifiedAndSavedVarnames');
    txtData = readtable(filename, 'FileType', 'text', 'Delimiter', '\t', 'ReadVariableNames', true);
    warning(ws); % revert to the previous warning state
    sens.label   = txtData{:,1};
    sens.elecpos = [txtData.Loc_X txtData.Loc_Y txtData.Loc_Z];

  case 'yorkinstruments_hdf5'
    acquisition='default';
    if isempty(senstype)
      % set the default
      ft_warning('both electrode and gradiometer information is present, returning the electrode information by default');
      senstype = 'eeg';
    end
    hdr=ft_read_header(filename);
    %i will be the channel index, sens_i is the sensor index
    sens_i=0;
    for i=1:hdr.nChans
      if string(hdr.chantype{i})==upper(senstype)
        sens_i = sens_i+1;
        sens.chantype{sens_i,1} = hdr.chantype{i};
        try
          sens.chanpos(sens_i,1:3) =  h5read(filename,['/config/channels/' hdr.label{i} '/position']);
          sens.chanori(sens_i,1:3) =  h5read(filename,['/config/channels/' hdr.label{i} '/orientation']);
          sens.chanunit{sens_i,1}  =  hdr.chanunit{i};
          sens.coilori(sens_i,1:3) =  sens.chanori(sens_i,1:3);
          sens.coilpos(sens_i,1:3) =  sens.chanpos(sens_i,1:3);
          sens.label{sens_i,1}     =  hdr.label{i};
        catch
          error('Error reading channel %i sensor details.',i);
        end
      else
        continue
      end
    end
    if sens_i<1
      error('No data corresponding to the chosen sensor type (%s) found.',senstype);
    end
    sens.tra  = eye(sens_i);
    sens.type= 'yorkinstruments248';
    if isempty(coordsys)
      coordsys = 'dewar';
    end
    if strcmp(coordsys, 'head')
      try
        tCCStoMegscanScs = h5read(filename, [strcat('/acquisitions/', char(string(acquisition))) '/ccs_to_scs_transform']);
        T = maketform('affine', tCCStoMegscanScs);
        sens.coilpos = tformfwd(T, sens.chanpos(:,1), sens.chanpos(:,2), sens.coilpos(:,3));
        R = tCCStoMegscanScs(1:3,1:3); % (mm)
        sens.coilori =  sens.coilori * R;
        sens.chanpos=sens.coilpos;
        sens.chanori=sens.coilori;
      catch
        error('No dewar to head transform available in hdf5 file');
      end
    end

  otherwise
    if ~isempty(sens)
      % the electrode or optode information has been read from the BIDS sidecar file
    else
      ft_error('unknown fileformat for electrodes or gradiometers');
    end

end % switch fileformat

% ensure that the sensor description is up-to-date
% this will also add chantype and units to the sensor array if missing
sens = ft_datatype_sens(sens);

% ensure that the output is consistent with the type requested by the user
if strcmpi(senstype, 'meg')
  assert(isfield(sens,'coilpos'), 'cannot read gradiometer information from %s', filename);
elseif strcmpi(senstype, 'eeg')
  assert(isfield(sens,'elecpos'), 'cannot read electrode information from %s', filename);
elseif strcmpi(senstype, 'nirs')
  assert(isfield(sens,'optopos'), 'cannot read optode information from %s', filename);
else
  % it is empty if not specified by the user, in that case either one is fine
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function fem1d ( )

%*****************************************************************************80
%
%% MAIN is the main program for FEM1D.
%
%  Discussion:
%
%    FEM1D solves a one dimensional ODE using the finite element method.
%
%    The differential equation has the form:
%
%      -d/dx ( p(x) du/dx ) + q(x) * u = f(x)
%
%    The finite-element method uses piecewise linear basis functions.
%
%    Here U is an unknown scalar function of X defined on the
%    interval [XL,XR], and P, Q and F are given functions of X.
%
%    The values of U or U' at XL and XR are also specified.
%
%    The interval [XL,XR] is "meshed" with NSUB+1 points,
%
%    XN(0) = XL, XN(1)=XL+H, XN(2)=XL+2*H, ..., XN(NSUB)=XR.
%
%    This creates NSUB subintervals, with interval number 1
%    having endpoints XN(0) and XN(1), and so on up to interval
%    NSUB, which has endpoints XN(NSUB-1) and XN(NSUB).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2008
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    real F(NU).
%    ASSEMBLE stores into F the right hand side of the linear
%    equations.
%    SOLVE replaces those values of F by the solution of the
%    linear equations.
%
%    real H(N), the length of the subintervals.
%
%    integer IBC. declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    integer INDX(1:N+1).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    integer NL, the number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    integer NQUAD.
%    The number of quadrature points used in a subinterval.
%    This code uses NQUAD = 1.
%
%    integer NSUB, the number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    integer NU, the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    real XL, the left endpoint of the interval over which the
%    differential equation is being solved.
%
%    real XN(1:N+1).
%    XN(I) is the location of the I-th node.  XN(1) is XL,
%    and XN(N+1) is XR.
%
%    real XQUAD(N)
%    XQUAD(I) is the location of the single quadrature point
%    in interval I.
%
%    real XR, the right endpoint of the interval over which the
%    differential equation is being solved.
%
  n = 100;
  nl = 2;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the two-point boundary value problem:\n' );
  fprintf ( 1, '    -d/dx (p(x) du/dx) + q(x)*u  =  f(x)\n' );
  fprintf ( 1, '  on an interval [xl,xr], with the values of\n' );
  fprintf ( 1, '  u or u'' specified at xl and xr.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The interval is broken into %d subintervals.\n', n );
  fprintf ( 1, '  The number of basis functions per element is %d\n', nl );
%
%  Initialize variables that define the problem.
%
  [ ibc, nquad, ul, ur, xl, xr ] = init ( );
%
%  Compute the quantities that describe the geometry of the problem.
%
  [ h, indx, node, nu, xn, xquad ] = geometry ( ibc, nl, n, xl, xr );
%
%  Assemble the matrix.
%
  [ adiag, aleft, arite, f ] = assemble ( h, indx, nl, node, ...
    nu, nquad, n, ul, ur, xn, xquad );
%
%  Print out the linear system.
%
  system_print ( adiag, aleft, arite, f, nu );
%
%  Solve the linear system.
%
  u = solve ( adiag, aleft, arite, f, nu );
%
%  Print the current solution.
%
  output ( u, ibc, indx, n, nu, ul, ur, xn );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ adiag, aleft, arite, f ] = assemble ( h, indx, nl, ...
  node, nu, nquad, n, ul, ur, xn, xquad )

%*****************************************************************************80
%
%% ASSEMBLE assembles the matrix and right hand side of the linear system.
%
%  Discussion:
%
%    Note that a 1 point quadrature rule, which is sometimes used to
%    assemble the matrix and right hand side, is just barely accurate
%    enough for simple problems.  If you want better results, you
%    should use a quadrature rule that is more accurate.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real H(N), the length of the subintervals.
%
%    Input, integer INDX(1:N+1).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Input, integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    Input, integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    Input, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Input, integer NQUAD, the number of quadrature points in a subinterval.
%    This code uses NQUAD = 1.
%
%    Input, integer N, the number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real XN(1:N+1), the location of the I-th node.  XN(1) is XL,
%    and XN(N+1) is XR.
%
%    Input, real XQUAD(N), the location of the single quadrature point
%    in interval I.
%
%    Output, real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    Output, real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    Output, real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    Output, real F(NU), the right hand side of the linear
%    equations.
%
  f(1:nu) = 0.0;
  adiag(1:nu) = 0.0;
  aleft(1:nu) = 0.0;
  arite(1:nu) = 0.0;
%
%  For element IE...
%
  for ie = 1 : n

    he = h(ie);
    xleft = xn(node(1,ie)+1);
    xrite = xn(node(2,ie)+1);
%
%  For quadrature point IQ...
%
    for iq = 1 : nquad

      xqe = xquad(ie);
%
%  For basis function IL...
%
      for il = 1 : nl

        ig = node(il,ie);
        iu = indx(ig+1);

        if ( 0 < iu )

          [ phii, phiix ] = phi ( il, xqe, xleft, xrite );

          f(iu) = f(iu) + he * ff ( xqe ) * phii;
%
%  Handle boundary conditions.
%
          if ( ig == 0 )

            x = xn(1);
            f(iu) = f(iu) - pp ( x ) * ul;

          elseif ( ig == n )

            x = xn(n+1);
            f(iu) = f(iu) + pp ( x ) * ur;

          end
%
%  For basis function JL...
%
          for jl = 1 : nl

            jg = node(jl,ie);
            ju = indx(jg+1);

            [ phij, phijx ] = phi ( jl, xqe, xleft, xrite );

            aij = he * ( pp ( xqe ) * phiix * phijx ...
                  + qq ( xqe ) * phii * phij );

            if ( ju <= 0 )

              if ( jg == 0 )
                f(iu) = f(iu) - aij * ul;
              elseif ( jg == n )
                f(iu) = f(iu) - aij * ur;
              end

            elseif ( iu == ju )
              adiag(iu) = adiag(iu) + aij;
            elseif ( ju < iu )
              aleft(iu) = aleft(iu) + aij;
            else
              arite(iu) = arite(iu) + aij;
            end

          end

        end

      end

    end

  end

  return
end
function value = ff ( x )

%*****************************************************************************80
%
%% FF returns the right hand side of the differential equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = 0.0;

  return
end
function [ h, indx, node, nu, xn, xquad ] = geometry ( ibc, nl, nsub, ...
  xl, xr )

%*****************************************************************************80
%
%% GEOMETRY sets up the geometry for the interval [XL,XR].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Input, integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    Input, integer NSUB.
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, real XL.
%    XL is the left endpoint of the interval over which the
%    differential equation is being solved.
%
%    Input, real XR.
%    XR is the right endpoint of the interval over which the
%    differential equation is being solved.
%
%    Output, real H(NSUB), the length of the subintervals.
%
%    Output, integer INDX(1:NSUB+1).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Output, integer NODE(NL,NSUB).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    Output, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be NSUB-1,
%    NSUB, or NSUB+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Output, real XN(1:NSUB+1).
%    XN(I) is the location of the I-th node.  XN(1) is XL,
%    and XN(N+1) is XR.
%
%    Output, real XQUAD(NSUB)
%    XQUAD(I) is the location of the single quadrature point
%    in interval I.
%

%
%  Set the value of XN, the locations of the nodes.
%
  xn = zeros ( nsub + 1, 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node      Location\n' );
  fprintf ( 1, '\n' );
  for i = 0 : nsub
    xn(i+1) =  ( ( nsub - i ) * xl   ...
               + (        i ) * xr ) ...
               / ( nsub     );
  end
  r8vec_print_some ( nsub + 1, xn, 1, 10, '  First 10 nodes:' );
%
%  Set the lengths of each subinterval.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Subint    Length\n' );
  fprintf ( 1, '\n' );
  h = zeros ( nsub, 1 );
  for i = 1 : nsub
    h(i) = xn(i+1) - xn(i);
  end
  r8vec_print_some ( nsub, h, 1, 10, '  First 10 interval widths:' );
%
%  Set the quadrature points, each of which is the midpoint of its subinterval.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Subint    Quadrature point\n' );
  fprintf ( 1, '\n' );
  xquad = zeros ( nsub );
  for i = 1 : nsub
    xquad(i) = 0.5 * ( xn(i) + xn(i+1) );
  end
  r8vec_print_some ( nsub, xquad, 1, 10, '  First 10 quadrature points:' );
%
%  Set the value of NODE, which records, for each interval,
%  the node numbers at the left and right.
%
  node = zeros ( 2, nsub );
  for i = 1 : nsub
    node(1,i) = i - 1;
    node(2,i) = i;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  First 10 pairs of nodes defining intervals:\n' );
  fprintf ( 1, 'Subint  Left Node  Right Node\n' );
  fprintf ( 1, '\n' );
  for i = 1 : min ( nsub, 10 )
    fprintf ( 1, '  %6d  %6d  %6d\n', i, node(1,i), node(2,i) );
  end
%
%  Starting with node 0, see if an unknown is associated with
%  the node.  If so, give it an index.
%
  nu = 0;
%
%  Handle first node.
%
  i = 0;
  if ( ibc == 1 || ibc == 3 )
    indx(i+1) = -1;
  else
    nu = nu + 1;
    indx(i+1) = nu;
  end
%
%  Handle nodes 1 through nsub-1
%
  for i = 1 : nsub-1
    nu = nu + 1;
    indx(i+1) = nu;
  end
%
%  Handle the last node.
%
  i = nsub;
  if ( ibc == 2 || ibc == 3 )
    indx(i+1) = -1;
  else
    nu = nu + 1;
    indx(i+1) = nu;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  First 10 unknown indices\n' );
  fprintf ( 1, '    Node  Unknown\n' );
  fprintf ( 1, '\n' );
  for i = 0 : min ( nsub, 9 )
    fprintf ( 1, '  %6d  %6d\n', i, indx(i+1) );
  end

  return
end
function [ ibc, nquad, ul, ur, xl, xr ] = init  ( )

%*****************************************************************************80
%
%% INIT initializes variables that define the problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Output, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Output, integer NQUAD.
%    The number of quadrature points used in a subinterval.
%    This code uses NQUAD = 1.
%
%    Output, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Output, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Output, real XL.
%    XL is the left endpoint of the interval over which the
%    differential equation is being solved.
%
%    Output, real XR.
%    XR is the right endpoint of the interval over which the
%    differential equation is being solved.
%
  ibc = 1;
  nquad = 1;
  ul = 0.0;
  ur = 1.0;
  xl = 0.0;
  xr = 1.0;
%
%  Print out the values that have been set.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'The equation is to be solved for\n' );
  fprintf ( 1, 'X greater than XL = %f\n', xl );
  fprintf ( 1, ' and less than XR = %f\n', xr );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'The boundary conditions are:\n' );
  fprintf ( 1, '\n' );

  if ( ibc == 1 || ibc == 3 )
    fprintf ( 1, '  At X = XL, U = %f\n', ul );
  else
    fprintf ( 1, '  At X = XL, U'' = %f\n', ul );
  end

  if ( ibc == 2 || ibc == 3 )
    fprintf ( 1, '  At X = XR, U = %f\n', ur );
  else
    fprintf ( 1, '  At X = XR, U'' = %f\n', ur );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Number of quadrature points per element is %d\n', nquad );

  return
end
function output ( f, ibc, indx, nsub, nu, ul, ur, xn )

%*****************************************************************************80
%
%% OUTPUT prints out the computed solution at the nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real F(NU), the solution of the linear equations.
%
%    Input, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Input, integer INDX(1:N+1).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    integer NSUB.
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real XN(1:N+1).
%    XN(I) is the location of the I-th node.  XN(1) is XL,
%    and XN(N+1) is XR.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'First 10 entries of computed solution:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node        X(I)          U(I)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : min ( nsub, 9 )

    if ( i == 0 )
      if ( ibc == 1 || ibc == 3 )
        u = ul;
      else
        u = f(indx(i+1));
      end
    elseif ( i == nsub )
      if ( ibc == 2 || ibc == 3 )
        u = ur;
      else
        u = f(indx(i+1));
      end
    else
      u = f(indx(i+1));
    end

    fprintf ( 1, '  %6d  %12f  %12f\n', i, xn(i+1), u );

  end

  return
end
function [ phii, phiix ] = phi ( il, x, xleft, xrite )

%*****************************************************************************80
%
%% PHI evaluates a linear basis function and its derivative.
%
%  Discussion:
%
%    In any interval, there are just two basis functions.  The first
%    basis function is a line which is 1 at the left endpoint
%    and 0 at the right.  The second basis function is 0 at
%    the left endpoint and 1 at the right.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer IL, the index of the basis function.
%    1, the function which is 1 at XLEFT and 0 at XRITE.
%    2, the function which is 0 at XLEFT and 1 at XRITE.
%
%    Input, real X, the evaluation point.
%
%    Input, real XLEFT, XRITE, the left and right
%    endpoints of the interval.
%
%    Output, real PHII, PHIIX, the value of the
%    basis function and its derivative at X.
%
  if ( xleft <= x && x <= xrite )

    if ( il == 1 )
      phii = ( xrite - x ) / ( xrite - xleft );
      phiix = -1.0 / ( xrite - xleft );
    else
      phii = ( x - xleft ) / ( xrite - xleft );
      phiix = 1.0 / ( xrite - xleft );
    end
%
%  If X is outside of the interval, then the basis function
%  is always zero.
%
  else

    phii = 0.0;
    phiix = 0.0;

  end

  return
end
function value = pp ( x )

%*****************************************************************************80
%
%% PP returns the value of the coefficient function P(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of P(X).
%
  value = 1.0;

  return
end
function value = qq ( x )

%*****************************************************************************80
%
%% QQ returns the value of the coefficient function Q(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of Q(X).
%
  value = 0.0;

  return
end
function r8vec_print_some ( n, a, i_lo, i_hi, title )

%*****************************************************************************80
%
%% R8VEC_PRINT_SOME prints "some" of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines to print.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for i = max ( 1, i_lo ) : min ( n, i_hi )
    fprintf ( 1, '  %8d: %12f\n', i, a(i) );
  end

  return
end
function u = solve ( adiag, aleft, arite, f, nu )

%*****************************************************************************80
%
%% SOLVE solves a tridiagonal matrix system of the form A*x = b.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real ADIAG(NU), ALEFT(NU), ARITE(NU).
%    the diagonal, left and right entries of the equations.
%    Note that for the first equation, there is no ALEFT
%    coefficient, and for the last, there is no ARITE.
%    So there is no need to store a value in ALEFT(1), nor
%    in ARITE(NU).
%
%    Input, real F(NU),  the right hand side of the linear
%    system to be solved.
%
%    Input, integer NU.
%    NU is the number of equations to be solved.
%
%    Output, real U(NU), the solution of the linear system.
%
  arite(1) = arite(1) / adiag(1);
  for i = 2 : nu-1
    adiag(i) = adiag(i) - aleft(i) * arite(i-1);
    arite(i) = arite(i) / adiag(i);
  end
  adiag(nu) = adiag(nu) - aleft(nu) * arite(nu-1);

  u = zeros ( nu, 1 );

  u(1) = f(1) / adiag(1);
  for i = 2 : nu
    u(i) = ( f(i) - aleft(i) * u(i-1) ) / adiag(i);
  end

  for i = nu-1 : -1 : 1
    u(i) = u(i) - arite(i) * u(i+1);
  end

  return
end
function system_print ( adiag, aleft, arite, f, nu )

%*****************************************************************************80
%
%% SYSTEM_PRINT prints out the tridiagonal linear system to be solved.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ADIAG(NU), ALEFT(NU), ARITE(NU),
%    the diagonal, left and right entries of the equations.
%
%    Input, real F(NU), the right hand side of the linear system.
%
%    Input, integer NU, the number of equations to be solved.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'First 10 rows of tridiagonal linear system:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Equation   ALEFT         ADIAG         ARITE         RHS\n' );
  fprintf ( 1, '\n' );

  for i = 1 : min ( nu, 10 )

    fprintf ( 1, '%3d', i );

    if ( i == 1 )
      fprintf ( 1, '              ' );
    else
      fprintf ( 1, '  %12f', aleft(i) );
    end

    fprintf ( 1, '  %12f', adiag(i) );

    if ( i < nu )
      fprintf ( 1, '  %12f', arite(i) );
    else
      fprintf ( 1, '              ' );
    end

    fprintf ( 1, '  %12f\n', f(i) );

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function [ a, b ] = p46_lim ( )

%*****************************************************************************80
%
%% P46_LIM returns the integration limits for problem 46.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 2.0 * pi;

  return
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% KERNEL/REGULARIZERS
%
% Contents of FAIR's Reguarization Toolbox
%
% general purpose tools:
%   contents                - this file
%   regularizer             - the specific regularizer used in FAIR
%     initialize: regularizer('reset','regularizer','mbElastic',...
%                             'alpha',1e3,'mu',1,'lambda',0);
%     use:        [Sc,dS,d2S] = regularizer(yc-yRef,omega,m);
%
%  The regularizers for NPIR are:
%  ------------------------------
%  curvature                curvature regularizer (cell-centered grid, requires parameter alpha)
%                           see, e.g., E10_2Ddisc2C_curvature
%  curvaturemexC            matrix free C version of curvature regularizer
%  curvatureDiagMex         matrix free access to the diagonal of the operator
%  curvatureHessianMex      matrix free access to the second derivative of the operator

%  elastic                  (linear) elastic regularizer (staggered grid)
%                             initialize: 'alpha', 'mu', 'lambda'
%                             see, e.g., E10_2Ddisc2C_elastic
%  hyperelastic             hyperelastic regularizer    (nodal grid)
%                             initialize: 'alphaLength', 'alphaArea', 'alphaVolume'
%                             see, e.g., E10_2Ddisc2C_hyperElastic
%
%  (linear) Differential operators are built in:
%
%  getCurvatureMatrix       generates curvature regularizer matrix (cell-centered grid)
%  getElasticMatrixNodal    generates elastic regularizer matrix (nodal grid)
%  getElasticMatrixStg      generates elastic regularizer matrix (staggeres grid)
%  getGradientNodal         generates gradient operator matrix (nodal grid)
%
%  For (nonlinear) hyperelastic regularization area and volume of tetrahedral partition
%  are computed in
%  geometry                contains areas and volumes functions
%                          used in hyperElastic
%  geometrymexC.h / .cpp   MEX version of area and volume computation
%
%
% =============================================================================

function debit = contents
if nargout == 0, help(mfilename); return; end;

debit = {
          'contents.m'

          'regularizer.m'

          'elastic.m'
          'getElasticMatrixNodal.m'
          'getElasticMatrixStg.m'

          'curvature.m'
          'curvatureMexC.cpp'
          'getCurvatureMatrix.m'
          'curvatureDiagMex.cpp'
          'curvatureHessianMex.cpp'

          'geometry.m'
          'geometryMexC.cpp'
          'geometryMexC.h'
          'geometryC.cpp'
          'geometryC.h'
          'hyperElastic.m'

          'getGradientNodal.m'
          'testRegularizers.m'
  };
%==============================================================================


text:
report_this_filefun(mfilename('fullpath'));

% drap a colormap onto topography
l = isnan(tmap);
tmap(l) = 1;

%l = tmap< 0.1;
%tmap(l) = 0;

[lat,lon] = meshgrat(tmap,tmapleg);
[smap,smapleg] = country2mtx('switzerland',100);
[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

sw = interp2(lon0,lat0,smap,lon,lat);

l =  isnan(sw) == 1 | sw == 2 ;

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
%l =  isnan(ren);
ren(l) = mi-0.01;

figure

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[45 48.5],'MapLonLimit',[5.3 11.5])

ll = tmap < 0 & ren < 0;
ren(ll) = ren(ll)*0 + 20;
meshm(ren,tmapleg,size(tmap),tmap);
meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',6);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

%   load worldlo
%   h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
%   h2 = displaym(PPpoint);
%   h = displaym(PPtext); trimcart(h);

pl = plotm(ma(:,2),ma(:,1),'hw');
set(pl,'LineWidth',1.5,'MarkerSize',12,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')
%load coast.mat
%c = coast;
% plotm(c(:,1), c(:,2),'k','Linewidth',1);
zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface

j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.9 0.9 0.9 ] ; j; [ 0.5 0.5 0.5] ];
caxis([ min(min(re4)) max(max(re4)) ]);

colormap(j); brighten(0.1);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',3);

setm(gca,'mlabellocation',2)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',2)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',12,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.8 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',12);
set(gcf,'Inverthardcopy','off');

author: CelsoReyes