text:
function Answer=inputdlg2(Prompt, Title, NumLines, DefAns, Resize)
%INPUTDLG Input dialog box.
%  ANSWER = INPUTDLG(PROMPT) creates a modal dialog box that returns user
%  input for multiple prompts in the cell array ANSWER. PROMPT is a cell
%  array containing the PROMPT strings.
%
%  INPUTDLG uses UIWAIT to suspend execution until the user responds.
%
%  ANSWER = INPUTDLG(PROMPT,NAME) specifies the title for the dialog.
%
%  ANSWER = INPUTDLG(PROMPT,NAME,NUMLINES) specifies the number of lines for
%  each answer in NUMLINES. NUMLINES may be a constant value or a column
%  vector having one element per PROMPT that specifies how many lines per
%  input field. NUMLINES may also be a matrix where the first column
%  specifies how many rows for the input field and the second column
%  specifies how many columns wide the input field should be.
%
%  ANSWER = INPUTDLG(PROMPT,NAME,NUMLINES,DEFAULTANSWER) specifies the
%  default answer to display for each PROMPT. DEFAULTANSWER must contain
%  the same number of elements as PROMPT and must be a cell array of
%  strings.
%
%  ANSWER = INPUTDLG(PROMPT,NAME,NUMLINES,DEFAULTANSWER,OPTIONS) specifies
%  additional options. If OPTIONS is the string 'on', the dialog is made
%  resizable. If OPTIONS is a structure, the fields Resize, WindowStyle, and
%  Interpreter are recognized. Resize can be either 'on' or
%  'off'. WindowStyle can be either 'normal' or 'modal'. Interpreter can be
%  either 'none' or 'tex'. If Interpreter is 'tex', the prompt strings are
%  rendered using LaTeX.
%
%  Examples:
%
%  prompt={'Enter the matrix size for x^2:','Enter the colormap name:'};
%  name='Input for Peaks function';
%  numlines=1;
%  defaultanswer={'20','hsv'};
%
%  answer=inputdlg(prompt,name,numlines,defaultanswer);
%
%  options.Resize='on';
%  options.WindowStyle='normal';
%  options.Interpreter='tex';
%
%  answer=inputdlg(prompt,name,numlines,defaultanswer,options);
%
%  See also DIALOG, ERRORDLG, HELPDLG, LISTDLG, MSGBOX,
%    QUESTDLG, TEXTWRAP, UIWAIT, WARNDLG .

%  Copyright 1994-2017 The MathWorks, Inc.

%%%%%%%%%%%%%%%%%%%%
%%% Nargin Check %%%
%%%%%%%%%%%%%%%%%%%%
narginchk(0,5);
nargoutchk(0,1);

if nargin > 0
    Prompt = convertStringsToChars(Prompt);
end

if nargin > 1
    Title = convertStringsToChars(Title);
end

if nargin > 3
    if isstring(DefAns)
        DefAns = cellstr(DefAns);
    end
end

if nargin > 4
    Resize = convertStringsToChars(Resize);
end

%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Handle Input Args %%%
%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<1
    Prompt=getString(message('MATLAB:uistring:popupdialogs:InputDlgInput'));
end
if ~iscell(Prompt)
    Prompt={Prompt};
end
NumQuest=numel(Prompt);

if nargin<2
    Title=' ';
end

if nargin<3
    NumLines=1;
end

if nargin<4
    DefAns=cell(NumQuest,1);
    for lp=1:NumQuest
        DefAns{lp}='';
    end
end

if nargin<5
    Resize = 'off';
end
WindowStyle='modal';
Interpreter='none';

Options = struct([]); %#ok
if nargin==5 && isstruct(Resize)
    Options = Resize;
    Resize  = 'off';
    if isfield(Options,'Resize'),      Resize=Options.Resize;           end
    if isfield(Options,'WindowStyle'), WindowStyle=Options.WindowStyle; end
    if isfield(Options,'Interpreter'), Interpreter=Options.Interpreter; end
end

[rw,cl]=size(NumLines);
OneVect = ones(NumQuest,1);
if (rw == 1 & cl == 2) %#ok Handle []
    NumLines=NumLines(OneVect,:);
elseif (rw == 1 & cl == 1) %#ok
    NumLines=NumLines(OneVect);
elseif (rw == 1 & cl == NumQuest) %#ok
    NumLines = NumLines';
elseif (rw ~= NumQuest | cl > 2) %#ok
    error(message('MATLAB:inputdlg:IncorrectSize'))
end

if ~iscell(DefAns)
    error(message('MATLAB:inputdlg:InvalidDefaultAnswer'));
end

%%%%%%%%%%%%%%%%%%%%%%%
%%% Create InputFig %%%
%%%%%%%%%%%%%%%%%%%%%%%
FigWidth=175;
FigHeight=100;
FigPos(3:4)=[FigWidth FigHeight];  %#ok
FigColor=get(0,'DefaultUicontrolBackgroundColor');

InputFig=dialog(                     ...
    'Visible'          ,'off'      , ...
    'KeyPressFcn'      ,@doFigureKeyPress, ...
    'Name'             ,Title      , ...
    'Pointer'          ,'arrow'    , ...
    'Units'            ,'pixels'   , ...
    'UserData'         ,'Cancel'   , ...
    'Tag'              ,Title      , ...
    'HandleVisibility' ,'callback' , ...
    'Color'            ,FigColor   , ...
    'NextPlot'         ,'add'      , ...
    'WindowStyle'      ,WindowStyle, ...
    'Resize'           ,Resize       ...
    );

%%%%%%%%%%%%%%%%%%%%%
%%% Set Positions %%%
%%%%%%%%%%%%%%%%%%%%%
DefOffset    = 5;
DefBtnWidth  = 53;
DefBtnHeight = 23;

TextInfo.Units              = 'pixels'   ;
TextInfo.FontSize           = get(0,'FactoryUicontrolFontSize');
TextInfo.FontWeight         = get(InputFig,'DefaultTextFontWeight');
TextInfo.HorizontalAlignment= 'left'     ;
TextInfo.HandleVisibility   = 'callback' ;

StInfo=TextInfo;
StInfo.Style              = 'text'  ;
StInfo.BackgroundColor    = FigColor;

EdInfo=StInfo;
EdInfo.FontWeight      = get(InputFig,'DefaultUicontrolFontWeight');
EdInfo.Style           = 'edit';
EdInfo.BackgroundColor = 'white';

BtnInfo=StInfo;
BtnInfo.FontWeight          = get(InputFig,'DefaultUicontrolFontWeight');
BtnInfo.Style               = 'pushbutton';
BtnInfo.HorizontalAlignment = 'center';

% Add VerticalAlignment here as it is not applicable to the above.
TextInfo.VerticalAlignment  = 'bottom';
TextInfo.Color              = get(0,'FactoryUicontrolForegroundColor');

% adjust button height and width
btnMargin=1.4;
ExtControl=uicontrol(InputFig   ,BtnInfo     , ...
    'String'   ,getString(message('MATLAB:uistring:popupdialogs:Cancel'))        , ...
    'Visible'  ,'off'         ...
    );

% BtnYOffset  = DefOffset;
BtnExtent = get(ExtControl,'Extent');
BtnWidth  = max(DefBtnWidth,BtnExtent(3)+8);
BtnHeight = max(DefBtnHeight,BtnExtent(4)*btnMargin);
delete(ExtControl);

% Determine # of lines for all Prompts
TxtWidth=FigWidth-2*DefOffset;
ExtControl=uicontrol(InputFig   ,StInfo     , ...
    'String'   ,''         , ...
    'Position' ,[ DefOffset DefOffset 0.96*TxtWidth BtnHeight ] , ...
    'Visible'  ,'off'        ...
    );

WrapQuest=cell(NumQuest,1);
QuestPos=zeros(NumQuest,4);

for ExtLp=1:NumQuest
    if size(NumLines,2)==2
        [WrapQuest{ExtLp},QuestPos(ExtLp,1:4)]= ...
            textwrap(ExtControl,Prompt(ExtLp),NumLines(ExtLp,2));
    else
        [WrapQuest{ExtLp},QuestPos(ExtLp,1:4)]= ...
            textwrap(ExtControl,Prompt(ExtLp),80);
    end
end % for ExtLp

delete(ExtControl);
QuestWidth =QuestPos(:,3);
QuestHeight=QuestPos(:,4);
if ismac % Change Edit box height to avoid clipping on mac.
    editBoxHeightScalingFactor = 1.4;
else
    editBoxHeightScalingFactor = 1;
end
TxtHeight=QuestHeight(1)/size(WrapQuest{1,1},1) * editBoxHeightScalingFactor;
EditHeight=TxtHeight*NumLines(:,1);
EditHeight(NumLines(:,1)==1)=EditHeight(NumLines(:,1)==1)+4;

FigHeight=(NumQuest+2)*DefOffset    + ...
    BtnHeight+sum(EditHeight) + ...
    sum(QuestHeight);

TxtXOffset=DefOffset;

QuestYOffset=zeros(NumQuest,1);
EditYOffset=zeros(NumQuest,1);
QuestYOffset(1)=FigHeight-DefOffset-QuestHeight(1);
EditYOffset(1)=QuestYOffset(1)-EditHeight(1);

for YOffLp=2:NumQuest
    QuestYOffset(YOffLp)=EditYOffset(YOffLp-1)-QuestHeight(YOffLp)-DefOffset;
    EditYOffset(YOffLp)=QuestYOffset(YOffLp)-EditHeight(YOffLp);
end % for YOffLp

QuestHandle=[];
EditHandle=[];

AxesHandle=axes('Parent',InputFig,'Position',[0 0 1 1],'Visible','off');

inputWidthSpecified = false;

for lp=1:NumQuest
    if ~ischar(DefAns{lp})
        delete(InputFig);
        error(message('MATLAB:inputdlg:InvalidInput'));
    end

    EditHandle(lp)=uicontrol(InputFig    , ...
        EdInfo      , ...
        'Max'        ,NumLines(lp,1)       , ...
        'Position'   ,[ TxtXOffset EditYOffset(lp) TxtWidth EditHeight(lp)], ...
        'String'     ,DefAns{lp}           , ...
        'Tag'        ,'Edit'                ...
        );
if strfind(WrapQuest{lp}{1},'save as')
    set(EditHandle(lp), 'Enable'     ,'inactive' );
    set(EditHandle(lp),'ButtonDownFcn',@(a,b) uiputfileeditbox(a,b))
elseif strfind(WrapQuest{lp}{1},'browse')
    set(EditHandle(lp), 'Enable'     ,'inactive' );
    set(EditHandle(lp),'ButtonDownFcn',@(a,b) uigetfileeditbox(a,b))
end

    QuestHandle(lp)=text('Parent'     ,AxesHandle, ...
        TextInfo     , ...
        'Position'   ,[ TxtXOffset QuestYOffset(lp)], ...
        'String'     ,WrapQuest{lp}                 , ...
        'Interpreter',Interpreter                   , ...
        'Tag'        ,'Quest'                         ...
        );

    MinWidth = max(QuestWidth(:));
    if (size(NumLines,2) == 2)
        % input field width has been specified.
        inputWidthSpecified = true;
        EditWidth = setcolumnwidth(EditHandle(lp), NumLines(lp,1), NumLines(lp,2));
        MinWidth = max(MinWidth, EditWidth);
    end
    % Get the extent of the text object. See g1008152
    questExtent = get(QuestHandle(lp), 'Extent');
    MinWidth = max(MinWidth, questExtent(3));
    FigWidth=max(FigWidth, MinWidth+2*DefOffset);

end % for lp

% fig width may have changed, update the edit fields if they dont have user specified widths.
if ~inputWidthSpecified
    TxtWidth=FigWidth-2*DefOffset;
    for lp=1:NumQuest
        set(EditHandle(lp), 'Position', [TxtXOffset EditYOffset(lp) TxtWidth EditHeight(lp)]);
    end
end

FigPos=get(InputFig,'Position');

FigWidth=max(FigWidth,2*(BtnWidth+DefOffset)+DefOffset);
FigPos(1)=0;
FigPos(2)=0;
FigPos(3)=FigWidth;
FigPos(4)=FigHeight;

set(InputFig,'Position',getnicedialoglocation(FigPos,get(InputFig,'Units')));

OKHandle=uicontrol(InputFig     ,              ...
    BtnInfo      , ...
    'Position'   ,[ FigWidth-2*BtnWidth-2*DefOffset DefOffset BtnWidth BtnHeight ] , ...
    'KeyPressFcn',@doControlKeyPress , ...
    'String'     ,getString(message('MATLAB:uistring:popupdialogs:OK'))        , ...
    'Callback'   ,@doCallback , ...
    'Tag'        ,'OK'        , ...
    'UserData'   ,'OK'          ...
    );

setdefaultbutton(InputFig, OKHandle);

CancelHandle=uicontrol(InputFig     ,              ...
    BtnInfo      , ...
    'Position'   ,[ FigWidth-BtnWidth-DefOffset DefOffset BtnWidth BtnHeight ]           , ...
    'KeyPressFcn',@doControlKeyPress            , ...
    'String'     ,getString(message('MATLAB:uistring:popupdialogs:Cancel'))    , ...
    'Callback'   ,@doCallback , ...
    'Tag'        ,'Cancel'    , ...
    'UserData'   ,'Cancel'       ...
    ); %#ok

handles = guihandles(InputFig);
handles.MinFigWidth = FigWidth;
handles.FigHeight   = FigHeight;
handles.TextMargin  = 2*DefOffset;
guidata(InputFig,handles);
set(InputFig,'ResizeFcn', {@doResize, inputWidthSpecified});

% make sure we are on screen
movegui(InputFig)

% if there is a figure out there and it's modal, we need to be modal too
if ~isempty(gcbf) && strcmp(get(gcbf,'WindowStyle'),'modal')
    set(InputFig,'WindowStyle','modal');
end

set(InputFig,'Visible','on');
drawnow;

if ~isempty(EditHandle)
    uicontrol(EditHandle(1));
end

if ishghandle(InputFig)
    % Go into uiwait if the figure handle is still valid.
    % This is mostly the case during regular use.
    c = matlab.ui.internal.dialog.DialogUtils.disableAllWindowsSafely();
    uiwait(InputFig);
    delete(c);
end

% Check handle validity again since we may be out of uiwait because the
% figure was deleted.
if ishghandle(InputFig)
    Answer={};
    if strcmp(get(InputFig,'UserData'),'OK')
        Answer=cell(NumQuest,1);
        for lp=1:NumQuest
            Answer(lp)=get(EditHandle(lp),{'String'});
        end
    end
    delete(InputFig);
else
    Answer={};
end
drawnow; % Update the view to remove the closed figure (g1031998)

function doFigureKeyPress(obj, evd) %#ok
switch(evd.Key)
    case {'return','space'}
        set(gcbf,'UserData','OK');
        uiresume(gcbf);
    case {'escape'}
        delete(gcbf);
end

function doControlKeyPress(obj, evd)
switch(evd.Key)
    case {'return'}
        if ~strcmp(get(obj,'UserData'),'Cancel')
            set(gcbf,'UserData','OK');
            uiresume(gcbf);
        else
            delete(gcbf)
        end
    case 'escape'
        delete(gcbf)
end

function doCallback(obj, evd) %#ok
if ~strcmp(get(obj,'UserData'),'Cancel')
    set(gcbf,'UserData','OK');
    uiresume(gcbf);
else
    delete(gcbf)
end

function doResize(FigHandle, evd, multicolumn) %#ok
% TBD: Check difference in behavior w/ R13. May need to implement
% additional resize behavior/clean up.

Data=guidata(FigHandle);

resetPos = false;

FigPos = get(FigHandle,'Position');
FigWidth = FigPos(3);
FigHeight = FigPos(4);

% the current and the target sizes are considered equal
% if the difference is less then 1 pixel
% (Position property values can be doubles due to other units,
%  and the non integer values show up due to conversions/round offs)
widthDiff = Data.MinFigWidth - FigWidth;
if widthDiff >= 1
    FigWidth  = Data.MinFigWidth;
    FigPos(3) = Data.MinFigWidth;
    resetPos = true;
end

% make sure edit fields use all available space if
% number of columns is not specified in dialog creation.
if ~multicolumn
    for lp = 1:length(Data.Edit)
        EditPos = get(Data.Edit(lp),'Position');
        EditPos(3) = FigWidth - Data.TextMargin;
        set(Data.Edit(lp),'Position',EditPos);
    end
end

% the current and the target sizes are considered equal
% if the difference is less then 1 pixel
% (Position property values can be doubles due to other units,
%  and the non integer values show up due to conversions/round offs)
heightDiff = abs(FigHeight - Data.FigHeight);
if heightDiff >= 1
    FigPos(4) = Data.FigHeight;
    resetPos = true;
end

if resetPos
    set(FigHandle,'Position',FigPos);
end

% set pixel width given the number of columns
function EditWidth = setcolumnwidth(object, rows, cols)
% Save current Units and String.
old_units = get(object, 'Units');
old_string = get(object, 'String');
old_position = get(object, 'Position');

set(object, 'Units', 'pixels')
set(object, 'String', char(ones(1,cols)*'x'));

new_extent = get(object,'Extent');
if (rows > 1)
    % For multiple rows, allow space for the scrollbar
    new_extent = new_extent + 19; % Width of the scrollbar
end
new_position = old_position;
new_position(3) = new_extent(3) + 1;
set(object, 'Position', new_position);

% reset string and units
set(object, 'String', old_string, 'Units', old_units);

EditWidth = new_extent(3);

function uiputfileeditbox(a,b)
[file, path] = uiputfile('*.nii;*.nii.gz','save as','Mask.nii.gz');
if file
    set(a,'String',fullfile(path,file));
end

function uigetfileeditbox(a,b)
[fileref,pathref] = uigetfile('*.nii;*.nii.gz','Select a reference','Reference.nii.gz');
if fileref
    set(a,'String',fullfile(pathref,fileref));
end

function figure_size = getnicedialoglocation(figure_size, figure_units)
% adjust the specified figure position to fig nicely over GCBF
% or into the upper 3rd of the screen

%  Copyright 1999-2010 The MathWorks, Inc.

parentHandle = gcbf;
convertData.destinationUnits = figure_units;
if ~isempty(parentHandle)
    % If there is a parent figure
    convertData.hFig = parentHandle;
    convertData.size = get(parentHandle,'Position');
    convertData.sourceUnits = get(parentHandle,'Units');
    c = [];
else
    % If there is no parent figure, use the root's data
    % and create a invisible figure as parent
    convertData.hFig = figure('visible','off');
    convertData.size = get(0,'ScreenSize');
    convertData.sourceUnits = get(0,'Units');
    c = onCleanup(@() close(convertData.hFig));
end

% Get the size of the dialog parent in the dialog units
container_size = hgconvertunits(convertData.hFig, convertData.size ,...
    convertData.sourceUnits, convertData.destinationUnits, get(convertData.hFig,'Parent'));

delete(c);

figure_size(1) = container_size(1)  + 1/2*(container_size(3) - figure_size(3));
figure_size(2) = container_size(2)  + 2/3*(container_size(4) - figure_size(4));

function setdefaultbutton(figHandle, btnHandle)
% WARNING: This feature is not supported in MATLAB and the API and
% functionality may change in a future release.

%SETDEFAULTBUTTON Set default button for a figure.
%  SETDEFAULTBUTTON(BTNHANDLE) sets the button passed in to be the default button
%  (the button and callback used when the user hits "enter" or "return"
%  when in a dialog box.
%
%  This function is used by inputdlg.m, msgbox.m, questdlg.m and
%  uigetpref.m.
%
%  Example:
%
%  f = figure;
%  b1 = uicontrol('style', 'pushbutton', 'string', 'first', ...
%       'position', [100 100 50 20]);
%  b2 = uicontrol('style', 'pushbutton', 'string', 'second', ...
%       'position', [200 100 50 20]);
%  b3 = uicontrol('style', 'pushbutton', 'string', 'third', ...
%       'position', [300 100 50 20]);
%  setdefaultbutton(b2);
%

%  Copyright 2005-2017 The MathWorks, Inc.

%--------------------------------------- NOTE ------------------------------------------
% This file was copied into matlab/toolbox/local/private.
% These two files should be kept in sync - when editing please make sure
% that *both* files are modified.

% Nargin Check
narginchk(1,2)

if (usejava('awt') == 1)
    % We are running with Java Figures
    useJavaDefaultButton(figHandle, btnHandle)
else
    % We are running with Native Figures
    useHGDefaultButton(figHandle, btnHandle);
end
function useJavaDefaultButton(figH, btnH)
% Get a UDD handle for the figure.
fh = handle(figH);
% Call the setDefaultButton method on the figure handle
fh.setDefaultButton(btnH);

function useHGDefaultButton(~, btnHandle)
% First get the position of the button.
btnPos = getpixelposition(btnHandle);

% Next calculate offsets.
leftOffset   = btnPos(1) - 1;
bottomOffset = btnPos(2) - 2;
widthOffset  = btnPos(3) + 3;
heightOffset = btnPos(4) + 3;

% Create the default button look with a uipanel.
% Use black border color even on Mac or Windows-XP (XP scheme) since
% this is in natve figures which uses the Win2K style buttons on Windows
% and Motif buttons on the Mac.
h1 = uipanel(get(btnHandle, 'Parent'), 'HighlightColor', 'black', ...
    'BorderType', 'etchedout', 'units', 'pixels', ...
    'Position', [leftOffset bottomOffset widthOffset heightOffset]);

% Make sure it is stacked on the bottom.
uistack(h1, 'bottom');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function cont=progress(status,rate)
%PROGRESS   Text progress bar
%   Similar to waitbar but without the figure display.
%
%   Start:
%      C = PROGRESS('init',TITLE), where the default TITLE is 'please wait'
%
%   On progress:
%      C = PROGRESS(C,RATE);
%
%   Examples:
%      n=20;
%      for i=1:n
%        if i==1, c=progress('init');
%        else     c=progress(c,i/n); end
%
%        % computing something ...
%        pause(.1)
%      end
%
%      % inside a script you may use:
%       c=progress('init','wait for ... whatever');
%       for i=1:n
%         c=progress(c,i/n);
%         ...
%       end
%

lmax=50;
if isequal(status,'init')
  if nargin > 1
    title = rate;
    fprintf(1,'\n   %s\n',title);
    str = repmat(' ',1,lmax-4);
  else
    str = '                please wait                   ';
    fprintf(1,'\n');
  end
  fprintf(1,'  |-%s-|\n',str);
  fprintf(1,'%s','  ');
  cont=0;
else
  cont = status;
  n=ceil(rate*lmax);
  dif = n - cont;
  if dif > 0 & n <=lmax
    N=dif;
  else
    N=0;
  end
  cont=cont+N;
  str = repmat('*',1,N);
  fprintf(1,'%s',str);
  if rate==1
   fprintf(1,'%s\n',' done');
  end
end

author: peterhcharlton

repo: RRest


text:
function [YY,XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V);

Ap(2,:) = 0;

for ii = 1:1:length(Ap)-1
    Ap(ii+1,ii) = 1;
end

inx = 1;
UUk = [Uk(inx); 0; 0; 0; 0];
PPk = (Xint_v*Xint_v');
VVk = [Vk(inx); 0; 0; 0; 0];
Qv = V*V';

for ii = 1:1:length(Xint_v)

XKk(ii,1) = Xint_v(ii)^2;                                             % FIRST STEP

end

PPk = Ap*PPk*Ap';                                                   % SECOND STEP

Kk = PPk*C'*inv( (C*PPk*C') + (V*Qv*V') );                          % THIRD STEP

for ii = 1:1:length(Xint_v)

XUPK(ii,1) = XKk(ii)^2 + UUk(ii);                                     % UPPER EQUATIONS.

Zk(ii,1) = cos(XUPK(ii)) +  VVk(ii);                                  % UPPER EQUATIONS.

end

for ii = 1:1:length(XKk)

XBARk(ii,1) = XKk(ii) + Kk(ii)*(Zk(ii) - (cos(XKk(ii)))) ;            % FOURTH STEP

end

II = eye(5,5);

Pk = ( II -  Kk*C)*PPk;                                             % FIFTH STEP

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

for ii = 1:1:n

UUk = [Uk(ii+1); 0; 0; 0; 0];
PPk = XBARk*XBARk';
VVk = [Vk(ii+1); 0; 0; 0; 0];

XKk = exp(-XBARk);                                                % FIRST STEP

PPkM = Ap*PPk*Ap';                                                 % SECOND STEP

Kk = PPkM*C'*inv( (C*PPkM*C') + (V*Qv*V') );                      % THIRD STEP

for nn = 1:1:length(XBARk)

XUPK(nn) = exp(-XKk(nn)) + UUk(nn);                              % UPPER EQUATIONS.

Zk(nn) = cos(XUPK(nn)) +  VVk(nn);                               % UPPER EQUATIONS.

end

for in = 1:1:length(XUPK)

XNEW(in) = XBARk(in) + Kk(in)*(Zk(in) - cos(XBARk(in)));           % FOURTH STEP

end

II = eye(5,5);

Pk = (II -  Kk*C)*PPkM;                                            % FIFTH STEP

XBARk = XNEW;

OUTX(ii) = XBARk(1,1);
OUTY(ii) = Zk(1,1);

end

YY = OUTY;
XX = OUTX;

author: Sable

repo: mcbench-benchmarks


text:
function test_bug2220(datainfo)

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY ft_preprocessing ft_preproc_padding preproc

%load('C:\Users\jorhor\Downloads\bugdatafile.mat')
cfg = [];
cfg.baselinewindow = [-.5 0];
cfg.demean = 'yes';

data = [];
data.label  = {'chan'};
data.fsample = 512;

for tr=1:10
  data.time{tr}    = linspace(-1, 2, 3*data.fsample+1);
  data.trial{tr}   = rand(1, length(data.time{tr})) .^ 2;
end

tmpdata = ft_preprocessing(cfg, data);

tbase = tmpdata.time{1} >= cfg.baselinewindow(1) & tmpdata.time{end} <= cfg.baselinewindow(2);

[ok, msg] = isalmostequal(sum(tmpdata.trial{1}(:, tbase), 2), zeros(length(data.label), 1),'abstol',eps*100);

if ~ok
  error('baseline not zero: %s', msg{1});
end

%% test with data on disk
if nargin<1
  datainfo = ref_datasets;
end

for k = 1:numel(datainfo)
  dataset = datainfo(k);

  cfg            = [];
  cfg.dataset    = fullfile(dataset.origdir,'original',dataset.type,dataset.datatype,'/',dataset.filename);

  % get header and event information
  if ~isempty(dataset.dataformat)
    hdr   = ft_read_header(cfg.dataset, 'headerformat', dataset.dataformat);
    event = ft_read_event(cfg.dataset, 'eventformat', dataset.dataformat);

    cfg.dataformat   = dataset.dataformat;
    cfg.headerformat = dataset.dataformat;
  else
    hdr   = ft_read_header(cfg.dataset);
    event = ft_read_event(cfg.dataset);
  end

  % create 10 1-second trials to be used as test-case
  begsample = ((1:10)-1)*round(hdr.Fs) + 1;
  endsample = ((1:10)  )*round(hdr.Fs);
  offset    = zeros(1,10);

  cfg.trl   = [begsample(:) endsample(:) offset(:)];
  sel       = cfg.trl(:,2)<=hdr.nSamples*hdr.nTrials;
  cfg.trl   = cfg.trl(sel,:);

  cfg.continuous     = 'yes';
  cfg.demean         = 'yes';
  cfg.baselinewindow = [0 0.5];
  if strcmp(dataset.type, 'meg') || strcmp(dataset.type, 'eeg')
    cfg.channel      = dataset.type;
  else
    cfg.channel      = 'all';
  end
  data               = ft_preprocessing(cfg);

  tbase = data.time{1} >= cfg.baselinewindow(1) & data.time{end} <= cfg.baselinewindow(2);
 [ok, msg] = isalmostequal(sum(data.trial{1}(1:end, tbase), 2), zeros(length(data.label), 1),'abstol',abs(.0001 * mean(data.trial{1}(:))));

  if ~ok
    error('baseline not zero: %s', msg{1});
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%--- help for dsge_solver_h ---
%
%  H1 line
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%
%    Other functions named dsge_solver_h
%
%       dsge/dsge_solver_h
%

author: jmaih


text:
function OutputFile = export_ssp(Projectors, ChannelNames, OutputFile)
% IMPORT_SSP: Saves SSP projectors to a file.
%
% USAGE:  SspFile = export_ssp(Projectors, ChannelNames, OutputFile=Ask)

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

% Ask for filename if not defined
if (nargin < 2) || isempty(OutputFile)
    % Get default directories and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    % Build a default file name
    OutputFile = bst_process('GetNewFilename', LastUsedDirs.ExportData, 'proj.mat');
    [fPath,fBase,fExt] = bst_fileparts(OutputFile);
    OutputFile = bst_fullfile(LastUsedDirs.ExportData, [fBase, fExt]);
    % Get filename where to store the filename
    [OutputFile, ProjFormat] = java_getfile('save', 'Save projectors', OutputFile, 'single', 'files', ...
                            {{'.fif'}, 'Elekta-Neuromag/MNE (*.fif)',        'FIF'; ...
                            {'_proj'}, 'Brainstorm SSP (*proj*.mat)', 'BST'}, 2);
    if isempty(OutputFile)
        return;
    end
    % Save new default export path
    LastUsedDirs.ExportData = bst_fileparts(OutputFile);
    bst_set('LastUsedDirs', LastUsedDirs);
end

% Get file extension
[fPath, fBase, fExt] = bst_fileparts(OutputFile);
% Save file
switch (fExt)
    case '.fif'
        out_projector_fif(OutputFile, ChannelNames, Projectors);
    case '.mat'
        NewMat.Projector = Projectors;
        NewMat.RowNames = ChannelNames;
        bst_save(OutputFile, NewMat, 'v7');
    otherwise
        error(['Unknown file extension "' fExt '"']);
end

author: brainstorm-tools


text:
%% Microstrip branch line-coupler analysis using 3D FDTD code with UPML
%% absorbing borders (ABC)
%
% Again we use FDTD 3D with UPML to calculate scattering coefficients S_{11},
% S_{21}, S_{31} and S_{41} for planar microstrip branch line-coupler
% following by the original paper by D. Sheen, S. Ali, M. Abouzahra, J. Kong
% "Application of the Three-Dimensional Finite-Difference Time-Domain Method
% to the Analysis of planar Microstrip Circuits", IEEE Trans. on Microwave
% Theory and Techniques (http://dx.doi.org/10.1109/22.55775).
%% Current code includes some improvements in comparison with the original
%% calculations:
% 1) imposing UPML instead of Mur ABCs;
% 2) using real metal (copper) as a patch conductor material instead of PEC;
% 3) applying matched load at the ends of filter's transmission microstrip
%    lines to prevent physical reflections;
% 4) no "magnetic wall" or "electric wall" conditions at the Ez source
%    plane;
% 5) using real properties for dielectric with losses (Duroid).
%
function FDTD_3D_Branch_Line_Coupler
close all; clear; clc;
%% Physical constants
   epsilon0 = 8.85418782e-12; mu0 = 1.25663706e-6;
   c = 1.0/sqrt(mu0*epsilon0);

%% Gaussian half-width
   t_half = 15.0e-12;

%% Microstrip transmission lines parameters
   lineW = 2.413e-3;
   lineH = 1.0e-3;
   % Roger's 5880 Duroid parameters
   lineEr = 2.2;     % eps_r
   lineTan = 0.0009; % loss tangent
   Z0 = 49.2526; % Matched load transmission line

%% End time
   t_end = 1.5e-9;

%% Total mesh dimensions and grid cells sizes (without PML)
   nx = 60; ny = 100; nz = 16;
   dx = 0.406e-3; dy = 0.406e-3; dz = 0.265e-3;

%% Number of PML layers
   PML = 5;

%% Matrix of material's constants
   number_of_materials = 4;
   % For material of number x = 1,2,3... :
   % Material(x,1) - relative permittivity, Material(x,2) - relative permeability,
   % Material(x,3) - specific conductivity
   % Vacuum
   Material(1,1) = 1.0;   Material(1,2) = 1.0;   Material(1,3) = 0.0;
   % Metal (Copper)
   Material(2,1) = 1.0;   Material(2,2) = 1.0;   Material(2,3) = 5.88e+7;
   % Substrate material (RT/Duroid 5880)
   Material(3,1) = lineEr;   Material(3,2) = 1.0;
   % Calculate conductivity of Duroid at 20 GHz from loss tangent and eps_r
   Material(3,3) = 2*pi*20e9*lineTan*lineEr*epsilon0;
   % Matched load material is calculated from transmission line parameters
   Material(4,1) = 1.0;   Material(4,2) = 1.0;   Material(4,3) = lineH/(Z0*lineW*dy);

% Add PML layers
   nx = nx + 2*PML; ny = ny + 2*PML; nz = nz + 2*PML;
% Calculate dt
   dt = (1.0/c/sqrt( 1.0/(dx^2) + 1.0/(dy^2) + 1.0/(dz^2)))*0.9999;
   number_of_iterations = ceil(t_end/dt);

%% 3D array for geometry
   Index = ones(nx, ny, nz);

%% Define of low-pass filter geometry
   % Ground plane
   Index((1+PML):(nx-PML-1), (1+PML):(ny-PML-1), PML+1) = 2;
   % Rectangular patch (thickness is equal to one cell)
   Index((nx/2-14):(nx/2+15), (ny/2-14):(ny/2-9), PML+5) = 2;
   Index((nx/2-14):(nx/2+15), (ny/2+9):(ny/2+14), PML+5) = 2;
   Index((nx/2-14):(nx/2-9), (1+PML):(ny-PML-1), PML+5) = 2;
   Index((nx/2+10):(nx/2+15), (1+PML):(ny-PML-1), PML+5) = 2;
   Index((nx/2-16):(nx/2-7), (ny/2-12):(ny/2+12), PML+5) = 2;
   Index((nx/2+8):(nx/2+17), (ny/2-12):(ny/2+12), PML+5) = 2;
   % Dielectric substrate between ground plane and filter patch
   Index((1+PML):(nx-PML-1), (1+PML):(ny-PML-1), (PML+2):(PML+4)) = 3;
   % Matched load before port 1 and 2
   Index((nx/2-14):(nx/2-9), PML+1, (PML+2):(PML+4)) = 4;
   Index((nx/2+10):(nx/2+15), PML+1, (PML+2):(PML+4)) = 4;
   % Matched load after port 2
   Index((nx/2-14):(nx/2-9), ny-PML-1, (PML+2):(PML+4)) = 4;
   Index((nx/2+10):(nx/2+15), ny-PML-1, (PML+2):(PML+4)) = 4;

%% 3D FDTD physical (fields) and additional arrays are defined as 'single'
%% to increase performance
   Ex = zeros(nx, ny+1, nz+1, 'single');
   Gx = zeros(nx, ny+1, nz+1, 'single');
   Fx = zeros(nx, ny+1, nz+1, 'single');
   Ey = zeros(nx+1, ny, nz+1, 'single');
   Gy = zeros(nx+1, ny, nz+1, 'single');
   Fy = zeros(nx+1, ny, nz+1, 'single');
   Ez = zeros(nx+1, ny+1, nz, 'single');
   Gz = zeros(nx+1, ny+1, nz, 'single');
   Fz = zeros(nx+1, ny+1, nz, 'single');
   Hx = zeros(nx+1, ny, nz, 'single');
   Bx = zeros(nx+1, ny, nz, 'single');
   Hy = zeros(nx, ny+1, nz, 'single');
   By = zeros(nx, ny+1, nz, 'single');
   Hz = zeros(nx, ny, nz+1, 'single');
   Bz = zeros(nx, ny, nz+1, 'single');

%% FDTD PML coefficients arrays. Here they are already filled with values
%% corresponding to free space
   m = 4; ka_max = 1.0; R_err = 1.0e-16;
   eta = sqrt(mu0/epsilon0*Material(1,1)/Material(1,2));
   k_Ex_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;
   k_Ex_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);
   k_Ey_a = ones(nx+1, ny, nz, 'single');
   k_Ey_b = ones(nx+1, ny, nz, 'single')/(2.0*epsilon0);
   k_Gz_a = ones(nx+1, ny, nz, 'single');
   k_Gz_b = ones(nx+1, ny, nz, 'single');
   k_Hy_a = ones(nx, ny, nz, 'single');
   k_Hy_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);
   k_Hx_c = ones(nx+1, ny, nz, 'single')*2.0*epsilon0/mu0;
   k_Hx_d = ones(nx+1, ny, nz, 'single')*(-2.0*epsilon0/mu0);
   k_Bz_a = ones(nx, ny, nz, 'single');
   k_Bz_b = ones(nx, ny, nz, 'single')*dt;
   k_Gx_a = ones(nx, ny+1, nz, 'single');
   k_Gx_b = ones(nx, ny+1, nz, 'single');
   k_Ey_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;
   k_Ey_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);
   k_Ez_a = ones(nx, ny+1, nz, 'single');
   k_Ez_b = ones(nx, ny+1, nz, 'single')/(2.0*epsilon0);
   k_Bx_a = ones(nx, ny, nz, 'single');
   k_Bx_b = ones(nx, ny, nz, 'single')*dt;
   k_Hy_c = ones(nx, ny+1, nz, 'single')*2.0*epsilon0/mu0;
   k_Hy_d = ones(nx, ny+1, nz, 'single')*(-2.0*epsilon0/mu0);
   k_Hz_a = ones(nx, ny, nz, 'single');
   k_Hz_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);
   k_Ex_a = ones(nx, ny, nz+1, 'single');
   k_Ex_b = ones(nx, ny, nz+1, 'single')/(2.0*epsilon0);
   k_Gy_a = ones(nx, ny, nz+1, 'single');
   k_Gy_b = ones(nx, ny, nz+1, 'single');
   k_Ez_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;
   k_Ez_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);
   k_Hx_a = ones(nx, ny, nz, 'single');
   k_Hx_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);
   k_By_a = ones(nx, ny, nz, 'single');
   k_By_b = ones(nx, ny, nz, 'single')*dt;
   k_Hz_c = ones(nx, ny, nz+1, 'single')*2.0*epsilon0/mu0;
   k_Hz_d = ones(nx, ny, nz+1, 'single')*(-2.0*epsilon0/mu0);

%% General FDTD coefficients
   I = 1:number_of_materials;
   K_a(I) = (2.0*epsilon0*Material(I,1) - Material(I,3)*dt)./...
            (2.0*epsilon0*Material(I,1) + Material(I,3)*dt);
   K_b(I) = 2.0*dt./(2.0*epsilon0*Material(I,1) + Material(I,3)*dt);
   K_c(I) = Material(I,2);
   Ka = single(K_a(Index)); Kb = single(K_b(Index)); Kc = single(K_c(Index));

%% PML coefficients along x-axis
   sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dx);
   for I=0:(PML-1)
        sigma_x = sigma_max*((PML - I)/PML)^m;
        ka_x = 1.0 + (ka_max - 1.0)*((PML - I)/PML)^m;
        k_Ey_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Ey_a(nx-I,:,:) = k_Ey_a(I+1,:,:);
        k_Ey_b(I+1,:,:) = 1.0/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Ey_b(nx-I,:,:) = k_Ey_b(I+1,:,:);
        k_Gz_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Gz_a(nx-I,:,:) = k_Gz_a(I+1,:,:);
        k_Gz_b(I+1,:,:) = 2.0*epsilon0/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Gz_b(nx-I,:,:) = k_Gz_b(I+1,:,:);
        k_Hx_c(I+1,:,:) = (2.0*epsilon0*ka_x + sigma_x*dt)/mu0;
        k_Hx_c(nx-I,:,:) = k_Hx_c(I+1,:,:);
        k_Hx_d(I+1,:,:) = -(2.0*epsilon0*ka_x - sigma_x*dt)/mu0;
        k_Hx_d(nx-I,:,:) = k_Hx_d(I+1,:,:);

        sigma_x = sigma_max*((PML - I - 0.5)/PML)^m;
        ka_x = 1.0 + (ka_max - 1.0)*((PML - I - 0.5)/PML)^m;
        k_Ex_c(I+1,:,:) = 2.0*epsilon0*ka_x + sigma_x*dt;
        k_Ex_c(nx-I-1,:,:) = k_Ex_c(I+1,:,:);
        k_Ex_d(I+1,:,:) = -(2.0*epsilon0*ka_x - sigma_x*dt);
        k_Ex_d(nx-I-1,:,:) = k_Ex_d(I+1,:,:);
        k_Hy_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Hy_a(nx-I-1,:,:) = k_Hy_a(I+1,:,:);
        k_Hy_b(I+1,:,:) = 1.0/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Hy_b(nx-I-1,:,:) = k_Hy_b(I+1,:,:);
        k_Bz_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Bz_a(nx-I-1,:,:) = k_Bz_a(I+1,:,:);
        k_Bz_b(I+1,:,:) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Bz_b(nx-I-1,:,:) = k_Bz_b(I+1,:,:);
   end

%% PML coefficients along y-axis
   sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dy);
   for J=0:(PML-1)
        sigma_y = sigma_max*((PML - J)/PML)^m;
        ka_y = 1.0 + (ka_max - 1.0)*((PML - J)/PML)^m;
        k_Gx_a(:,J+1,:) = (2.0*epsilon0*ka_y - sigma_y*dt)/...
                          (2.0*epsilon0*ka_y + sigma_y*dt);
        k_Gx_a(:,ny-J,:) = k_Gx_a(:,J+1,:);
        k_Gx_b(:,J+1,:) = 2.0*epsilon0/(2.0*epsilon0*ka_y + sigma_y*dt);
        k_Gx_b(:,ny-J,:) = k_Gx_b(:,J+1,:);
        k_Ez_a(:,J+1,:) = (2.0*epsilon0*ka_y - sigma_y*dt)/...
                          (2.0*epsilon0*ka_y + sigma_y*dt);
        k_Ez_a(:,ny-J,:) = k_Ez_a(:,J+1,:);
        k_Ez_b(:,J+1,:) = 1.0/(2.0*epsilon0*ka_y + sigma_y*dt);
        k_Ez_b(:,ny-J,:) = k_Ez_b(:,J+1,:);
        k_Hy_c(:,J+1,:) = (2.0*epsilon0*ka_y + sigma_y*dt)/mu0;
        k_Hy_c(:,ny-J,:) = k_Hy_c(:,J+1,:);
        k_Hy_d(:,J+1,:) = -(2.0*epsilon0*ka_y - sigma_y*dt)/mu0;
        k_Hy_d(:,ny-J,:) = k_Hy_d(:,J+1,:);

        sigma_y = sigma_max*((PML - J - 0.5)/PML)^m;
        ka_y = 1.0 + (ka_max - 1.0)*((PML - J - 0.5)/PML)^m;
        k_Ey_c(:,J+1,:) = 2.0*epsilon0*ka_y+sigma_y*dt;
        k_Ey_c(:,ny-J-1,:) = k_Ey_c(:,J+1,:);
        k_Ey_d(:,J+1,:) = -(2.0*epsilon0*ka_y-sigma_y*dt);
        k_Ey_d(:,ny-J-1,:) = k_Ey_d(:,J+1,:);
        k_Bx_a(:,J+1,:) = (2.0*epsilon0*ka_y-sigma_y*dt)/...
                          (2.0*epsilon0*ka_y+sigma_y*dt);
        k_Bx_a(:,ny-J-1,:) = k_Bx_a(:,J+1,:);
        k_Bx_b(:,J+1,:) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_y+sigma_y*dt);
        k_Bx_b(:,ny-J-1,:) = k_Bx_b(:,J+1,:);
        k_Hz_a(:,J+1,:) = (2.0*epsilon0*ka_y-sigma_y*dt)/...
                          (2.0*epsilon0*ka_y+sigma_y*dt);
        k_Hz_a(:,ny-J-1,:) = k_Hz_a(:,J+1,:);
        k_Hz_b(:,J+1,:) = 1.0/(2.0*epsilon0*ka_y+sigma_y*dt);
        k_Hz_b(:,ny-J-1,:) = k_Hz_b(:,J+1,:);
   end

%% PML coefficients along z-axis
   sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dz);
   for K=0:(PML-1)
        sigma_z = sigma_max*((PML - K)/PML)^m;
        ka_z = 1.0 + (ka_max - 1.0)*((PML-K)/PML)^m;
        k_Ex_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_Ex_a(:,:,nz-K) = k_Ex_a(:,:,K+1);
        k_Ex_b(:,:,K+1) = 1.0/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_Ex_b(:,:,nz-K) = k_Ex_b(:,:,K+1);
        k_Gy_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_Gy_a(:,:,nz-K) = k_Gy_a(:,:,K+1);
        k_Gy_b(:,:,K+1) = 2.0*epsilon0/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_Gy_b(:,:,nz-K) = k_Gy_b(:,:,K+1);
        k_Hz_c(:,:,K+1) = (2.0*epsilon0*ka_z + sigma_z*dt)/mu0;
        k_Hz_c(:,:,nz-K) = k_Hz_c(:,:,K+1);
        k_Hz_d(:,:,K+1) = -(2.0*epsilon0*ka_z - sigma_z*dt)/mu0;
        k_Hz_d(:,:,nz-K) = k_Hz_d(:,:,K+1);

        sigma_z = sigma_max*((PML - K - 0.5)/PML)^m;
        ka_z = 1.0 + (ka_max - 1.0)*((PML - K - 0.5)/PML)^m;
        k_Ez_c(:,:,K+1) = 2.0*epsilon0*ka_z + sigma_z*dt;
        k_Ez_c(:,:,nz-K-1) = k_Ez_c(:,:,K+1);
        k_Ez_d(:,:,K+1) = -(2.0*epsilon0*ka_z - sigma_z*dt);
        k_Ez_d(:,:,nz-K-1) = k_Ez_d(:,:,K+1);
        k_Hx_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_Hx_a(:,:,nz-K-1) = k_Hx_a(:,:,K+1);
        k_Hx_b(:,:,K+1) = 1.0/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_Hx_b(:,:,nz-K-1) = k_Hx_b(:,:,K+1);
        k_By_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_By_a(:,:,nz-K-1) = k_By_a(:,:,K+1);
        k_By_b(:,:,K+1) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_By_b(:,:,nz-K-1) = k_By_b(:,:,K+1);
   end

%% Main 3D FDTD+UPML routine (operates with 'singles' to increase speed)
   hhh = waitbar(0, 'Calculations in progress...');
   tic;
   for T=0:(number_of_iterations-1)
        %% Calculate Fx -> Gx -> Ex
        I = 1:nx; J = 2:ny; K = 2:nz;
        Fx_r = Fx(I,J,K);
        Fx(I,J,K) = Ka(I,J,K).*Fx(I,J,K) + Kb(I,J,K).*...
                    ((Hz(I,J,K) - Hz(I,J-1,K))/dy - (Hy(I,J,K) - Hy(I,J,K-1))/dz);
        Gx_r = Gx(I,J,K);
        Gx(I,J,K) = k_Gx_a(I,J,K).*Gx(I,J,K) + k_Gx_b(I,J,K).*(Fx(I,J,K) - Fx_r);
        Ex(I,J,K) = k_Ex_a(I,J,K).*Ex(I,J,K) + k_Ex_b(I,J,K).*...
                    (k_Ex_c(I,J,K).*Gx(I,J,K) + k_Ex_d(I,J,K).*Gx_r);

        %% Calculate Fy -> Gy -> Ey
        I = 2:nx; J = 1:ny; K = 2:nz;
        Fy_r = Fy(I,J,K);
        Fy(I,J,K) = Ka(I,J,K).*Fy(I,J,K) + Kb(I,J,K).*...
                    ((Hx(I,J,K) - Hx(I,J,K-1))/dz - (Hz(I,J,K) - Hz(I-1,J,K))/dx);
        Gy_r = Gy(I,J,K);
        Gy(I,J,K) = k_Gy_a(I,J,K).*Gy(I,J,K) + k_Gy_b(I,J,K).*(Fy(I,J,K) - Fy_r);
        Ey(I,J,K) = k_Ey_a(I,J,K).*Ey(I,J,K) + k_Ey_b(I,J,K).*...
                    (k_Ey_c(I,J,K).*Gy(I,J,K) + k_Ey_d(I,J,K).*Gy_r);

        %% Calculate Fz -> Gz -> Ez
        I = 2:nx; J = 2:ny; K = 1:nz;
        Fz_r = Fz(I,J,K);
        Fz(I,J,K) = Ka(I,J,K).*Fz(I,J,K) + Kb(I,J,K).*...
                    ((Hy(I,J,K) - Hy(I-1,J,K))/dx - (Hx(I,J,K) - Hx(I,J-1,K))/dy);
        Gz_r = Gz(I,J,K);
        Gz(I,J,K) = k_Gz_a(I,J,K).*Gz(I,J,K) + k_Gz_b(I,J,K).*(Fz(I,J,K) - Fz_r);
        Ez(I,J,K) = k_Ez_a(I,J,K).*Ez(I,J,K) + k_Ez_b(I,J,K).*...
                    (k_Ez_c(I,J,K).*Gz(I,J,K) + k_Ez_d(I,J,K).*Gz_r);

        %% Source of vertical electric field Ez applied to the whole face
        %% plane at the port 1
        if (T*dt<=10.0*t_half)
            Ez((nx/2-14):(nx/2-9), PML+2, (PML+2):(PML+4)) = Source(t_half, 1.0, T*dt);
        end

        %% Save reflected Ez at ports 1 and 2, as well as passed through ports 3 and 4
        Incident(T+1) = Source(t_half, 1.0, T*dt);
        % Remove incident signal from Port 1
        Port_1(T+1) = Ez(nx/2-11, PML+2, PML+3) - Incident(T+1);
        Port_2(T+1) = Ez(nx/2+13, PML+2, PML+3);
        Port_3(T+1) = Ez(nx/2+13, ny-PML-2, PML+3);
        Port_4(T+1) = Ez(nx/2-11, ny-PML-2, PML+3);
        Ez_saved(1:(ny-2*PML-1), 1:(nx-2*PML-1), T+1) = ...
                Ez((PML+1):(nx-PML-1), (PML+1):(ny-PML-1), PML+3)';

        %% Calculate Bx -> Hx
        I = 2:nx; J = 1:ny; K = 1:nz;
        Bx_r = Bx(I,J,K);
        Bx(I,J,K) = k_Bx_a(I,J,K).*Bx(I,J,K) + k_Bx_b(I,J,K).*...
                    ((Ey(I,J,K+1) - Ey(I,J,K))/dz - (Ez(I,J+1,K) - Ez(I,J,K))/dy);
        Hx(I,J,K) = k_Hx_a(I,J,K).*Hx(I,J,K) + k_Hx_b(I,J,K).*...
                    (k_Hx_c(I,J,K).*Bx(I,J,K) + k_Hx_d(I,J,K).*Bx_r)./Kc(I,J,K);

        %% Calculate By -> Hy
        I = 1:nx; J = 2:ny; K = 1:nz;
        By_r = By(I,J,K);
        By(I,J,K) = k_By_a(I,J,K).*By(I,J,K) + k_By_b(I,J,K).*...
                    ((Ez(I+1,J,K) - Ez(I,J,K))/dx - (Ex(I,J,K+1) - Ex(I,J,K))/dz);
        Hy(I,J,K) = k_Hy_a(I,J,K).*Hy(I,J,K) + k_Hy_b(I,J,K).*...
                    (k_Hy_c(I,J,K).*By(I,J,K) + k_Hy_d(I,J,K).*By_r)./Kc(I,J,K);

        %% Calculate Bz -> Hz
        I = 1:nx; J = 1:ny; K = 2:nz;
        Bz_r = Bz(I,J,K);
        Bz(I,J,K) = k_Bz_a(I,J,K).*Bz(I,J,K) + k_Bz_b(I,J,K).*...
                    ((Ex(I,J+1,K) - Ex(I,J,K))/dy - (Ey(I+1,J,K) - Ey(I,J,K))/dx);
        Hz(I,J,K) = k_Hz_a(I,J,K).*Hz(I,J,K) + k_Hz_b(I,J,K).*...
                    (k_Hz_c(I,J,K).*Bz(I,J,K) + k_Hz_d(I,J,K).*Bz_r)./Kc(I,J,K);

        %% Progress bar updates at every time percent
        if ( mod(T, ceil(number_of_iterations/100)) == 0 )
            waitbar((T+1)/number_of_iterations, hhh);
        end

    end
    toc
    close(hhh);

%% Plots: signals, scattering parameters and Ez field between plates
   figure('units','normalized','outerposition',[0 0 1 1]);
   set(gcf, 'doublebuffer', 'on');
   % Input signal, reflected (at port 1) and passed through (port 2)
     subplot(2,2,1);
     plot(1e9*dt*(0:T), Incident, 'm-', ...
          1e9*dt*(0:T), Port_1, 'k-', ...
          1e9*dt*(0:T), Port_2, 'r-', ...
          1e9*dt*(0:T), Port_3, 'b-', ...
          1e9*dt*(0:T), Port_4, 'g-', 'LineWidth', 2);
     title('Input and Port 1-4 signals', 'FontSize', 18);
     xlabel('Time, [ns]', 'FontSize', 16);
     ylabel('Strength, [V/cm]', 'FontSize', 16);
     legend('Input' ,'Port 1', 'Port 2', 'Port 3', 'Port 4');
     grid on;
   % |S_{1,2,3,4}|
     NFFT = 2^nextpow2(T);
     S_11 = fft(Port_1, NFFT)./fft(Incident, NFFT);
     S_21 = fft(Port_2, NFFT)./fft(Incident, NFFT);
     S_31 = fft(Port_3, NFFT)./fft(Incident, NFFT);
     S_41 = fft(Port_4, NFFT)./fft(Incident, NFFT);
     subplot(2,2,3);
     f = (0.5/dt)*linspace(0, 1, NFFT/2);
     plot(f*1e-9, 20*log10(abs(S_11(1:NFFT/2))), 'k', ...
          f*1e-9, 20*log10(abs(S_21(1:NFFT/2))), 'r', ...
          f*1e-9, 20*log10(abs(S_31(1:NFFT/2))), 'b', ...
          f*1e-9, 20*log10(abs(S_41(1:NFFT/2))), 'g', 'LineWidth', 2);
     title('Scattering parameters - |S_{1,2,3,4}|', 'FontSize', 18);
     xlabel('Frequency, [GHz]', 'FontSize', 16);
     ylabel('Magnitude, [dB]', 'FontSize', 16);
     legend('Port 1', 'Port 2', 'Port 3', 'Port 4', 'Location', 'SouthEast');
     xlim([0 10]);
     grid on;
   % Electric field strength Ez between electrodes
     for T=1:20:number_of_iterations
         subplot(2,2,[2 4]);
         contour(1e3*dx*(1:(nx-2*PML-1)), 1e3*dy*(1:(ny-2*PML-1)), ...
         Index((PML+1):(nx-PML-1), (PML+1):(ny-PML-1), PML+5)', 1);
         hold on;
         subplot(2,2,[2 4]);
         pcolor(1e3*dx*(1:(nx-2*PML-1)), 1e3*dy*(1:(ny-2*PML-1)), double(Ez_saved(:, :, T)));
         title(['Time t = ', num2str(round(T*dt*1e12)),' ps'], 'FontSize', 18);
         xlabel('[mm]', 'FontSize', 16);
         ylabel('[mm]', 'FontSize', 16);
         shading interp;
         axis image;
         caxis([-1 1]);
         colorbar;
         drawnow;
     end

%% Gauss function for voltage source
function [res] = Source(t_half, amplitude, t)
% Pulse delay
t0 = 3.0*t_half;
res = amplitude*exp(-0.5*((t-t0)/(t_half/2.35482))^2.0);


text:
function [s, vw] = openGrayWindow
%
%   [s, vw] = openGrayWindow
%
% Calls openRawVolumeWindow to set up the VOLUME data structure,
% opens and initializes the volume window.  Then calls switch2Gray.
%
% djh, 2/99
% Junjie 2004.02.14 let function return s

mrGlobals

s = openRawVolumeWindow;
VOLUME{s} = switch2Gray(VOLUME{s});

selectView(VOLUME{s});

% If user requested a view output, give it to them:
if nargout > 1,	vw = VOLUME{s}; end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% WriteSegy : writes data to disk using SEGY REV 1 standard.
%
% EX
% WriteSegy('datacube.segy',data,'dt',.004,'Inline3D',Inline,'Crossline3D',Crossline,'cdpX',X,'cdpY',Y);
%
% to use a specific SEG revision use :
% WriteSegy('test.segy',seisdata,'revision',0); % SEG-Y Revision 0
% WriteSegy('test.segy',seisdata,'revision',1); % SEG-Y Revision 1
%
% to use a specific Data Sampling Format use :
% WriteSegy('test.segy',seisdata,'dsf',1); % IBM FLAOTING POINT
%
% Forice Revision 1 and IEEE Floating point :
% WriteSegy('test.segy',seisdata,'dsf',5,'revision',1);
%
% See also : WriteSegyStructure, WriteSu, WriteSuStructure
%

%
% (C) 2001-2007, Thomas Mejer Hansen, tmh@gfy.ku.dk/thomas.mejer.hansen@gmail.com
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function WriteSegy(filename,data,varargin);

    [ns,ntraces]=size(data);

    for i=1:2:length(varargin)
        var=varargin{i};
        val=varargin{i+1};
        eval([var,'=[',num2str(val),'];']);
    end

    if exist('dt')~=1
        dt=0.004;
        SegymatVerbose([mfilename,' - No dt set. Using dt=',num2str(dt)],0)
    end

    %%%%%%%%%%
    % SET UP SegyHeader structure.
    if exist('TextualFileHeader'), SegyHeader.TextualFileHeader=TextualFileHeader; end
    if exist('Job')==1, SegyHeader.Job=Job; end;
    if exist('Line')==1, SegyHeader.Line=Line; end
    if exist('Reel')==1, SegyHeader.Reel=Reel; end
    if exist('DataTracePerEnsemble')==1, SegyHeader.DataTracePerEnsemble=DataTracePerEnsemble; end
    if exist('AuxiliaryTracePerEnsemble')==1, SegyHeader.AuxiliaryTracePerEnsemble=AuxiliaryTracePerEnsemble; end
    if exist('dt')==1, SegyHeader.dt=dt.*1e+6; end
    if exist('dtOrig')==1, SegyHeader.dtOrig=dtOrig; end
    if exist('ns')==1, SegyHeader.ns=ns; end
    if exist('nsOrig')==1, SegyHeader.nsOrig=nsOrig; end
    if exist('EnsembleFold')==1, SegyHeader.EnsembleFold=EnsembleFold;  end
    if exist('TraceSorting')==1, SegyHeader.TraceSorting=TraceSorting; end
    if exist('VerticalSumCode')==1, SegyHeader.VerticalSumCode=VerticalSumCode; end

    %%%%%%%%%%
    %% CHOOSE WHICH SEGY REVISION AND DATATYPE TO USE
    %% DEFAULT IS REVISION 1, 4-byte IEEE

    % IF A SPECFIC REVISION HAS BEEN CHOSEN, USE THAT
    if exist('revision')==1,
        if revision==1,
            SegyHeader.SegyFormatRevisionNumber=256;
        else
            SegyHeader.SegyFormatRevisionNumber=0;
        end
        SegymatVerbose([mfilename,' :  Using user specified SEG Y revision : ',num2str(revision)],1)
    end

    % IF A SPECFIC DATA SAMPLING FORMAT HAS BEEN SELECTED USE THAT
    if exist('dsf')==1,
        SegyHeader.DataSampleFormat=dsf;
        SegymatVerbose([mfilename,' :  Using user specified Data Sample Format : ',num2str(revision)],1)
    end

    % OPEN SEGY FILE HANDLE
    segyid = fopen(filename,'w','b');   % ALL DISK FILES ARE IN BIG
                                        % ENDIAN FORMAT, ACCORDING TO SEG
                                        % Y rev 1

    % Write SEGY HEADER
    SegyHeader=PutSegyHeader(segyid,SegyHeader);

    for i=1:ntraces;
        if (i/100)==round(i/100),
            SegymatVerbose(['writing trace ',num2str(i),' of ',num2str(ntraces)],0);
        end
        % Basic TraceHeader information

        % INITALIZE SEGY TRACE HEADER
        if exist('SegyTraceHeader')==0;
            SegyTraceHeader=InitSegyTraceHeader(ns,dt*1e+6);
        end

        if exist('TraceNumber')==0
            SegyTraceHeader.TraceNumber=i;
            SegyTraceHeader.TraceSequenceFile=i;
        end

        % Update TraceHeader information if available
        if exist('cdpX')==1,SegyTraceHeader.cdpX = cdpX(i);end
        if exist('offset')==1,SegyTraceHeader.offset = offset(i);end
        if exist('cdpY')==1,SegyTraceHeader.cdpY = cdpY(i);end
        if exist('Inline3D')==1,SegyTraceHeader.Inline3D = Inline3D(i);end
        if exist('Crossline3D')==1,SegyTraceHeader.Crossline3D=Crossline3D(i);end

        if exist('YearDataRecorded')==1,SegyTraceHeader.YearDataRecorded=YearDataRecorded(i);end
        if exist('DayOfYear')==1,SegyTraceHeader.DayOfYear=DayOfYear(i);end
        if exist('HourOfDay')==1,SegyTraceHeader.HourOfDay=HourOfDay(i);end
        if exist('MinuteOfOur')==1,SegyTraceHeader.MinuteOfOur=MinuteOfOur(i);end
        if exist('SecondOfMinute')==1,SegyTraceHeader.SecondOfMinute=SecondOfMinute(i);end
        if exist('TimeBaseCode')==1,SegyTraceHeader.TimeBaseCode=TimeBaseCode(i);end

        % Write the Trace
        PutSegyTrace(segyid,data(:,i),SegyTraceHeader,SegyHeader);
    end

    % CLOSE SEGY FILE HANDLE

    fclose(segyid);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function eciVector = rotateVectorFromRsw2Eci(rswVector, rVect, vVect)
    R = normVector(rVect);
    W = normVector(crossARH(rVect, vVect));
    S = normVector(crossARH(W,R));

    RSW2ECIRotMat = [R,S,W]';

    eciVector = RSW2ECIRotMat \ rswVector;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function qr = rightQuaternionProductMatrix(q)
% q 4x1 w,x,y,z
% p * q = rightQuaternionProductMatrix(q) * p
% eq 18 in Sola Quaternion kinematics for the error-state Kalman filter
w = q(1);
x = q(2);
y = q(3);
z = q(4);
qr = [w, -x, -y, -z;
    x, w, z, -y;
    y, -z, w, x;
    z, y, -x, w];
end

repo: ekfmonoslam


text:
function pickStatrStopNodes(src,evt,hAxis,axisType,~)

global airwayStateS

% set(src,'Marker','o','MarkerSize',12,'MarkerEdgeColor',[0,0,0])
xV = get(src,'XData');
yV = get(src,'YData');
zV = get(src,'ZData');
%x = mean(xV);
%y = mean(yV);
%z = mean(zV);
distV = (xV-evt.IntersectionPoint(1)).^2 + (yV-evt.IntersectionPoint(2)).^2 + ...
    (zV-evt.IntersectionPoint(3)).^2;
[minDist,indMin] = min(distV);
if minDist > 0.5
    return;
end
x = xV(indMin(1));
y = yV(indMin(1));
z = zV(indMin(1));

if strcmpi(axisType,'base')
    if airwayStateS.selectBaseStartNode % size(airwayStateS.baseSegmentStartNodes,1) > size(airwayStateS.baseSegmentEndNodes,1)
        %airwayStateS.baseSegmentEndNodes(end+1,:) = [x,y,z];
        airwayStateS.baseSegmentS(airwayStateS.currentSegment).startNode = indMin(1);
        airwayStateS.selectBaseStartNode = 0;
        set(airwayStateS.hBaseStart,'XData',x,'YData',y,'zData',z,'visible','on')
        %set(airwayStateS.hBaseEnd,'visible','off')
%         airwayStateS.basehEndNodes = [airwayStateS.basehEndNodes,...
%             plot3(x,y,z,'parent',hAxis,'marker','o','MarkerFaceColor','r',...
%             'markerSize',8,'hittest','off')];
    else
        %airwayStateS.baseSegmentStartNodes(end+1,:) = [x,y,z];
        %set(airwayStateS.hBaseEnd,'XData',x,'YData',y,'zData',z,'visible','on')
        airwayStateS.baseSegmentS(airwayStateS.currentSegment).endNode = indMin(1);
        airwayGraph = airwayStateS.baseTreeS.airwayGraph;
        startNode = airwayStateS.baseSegmentS(airwayStateS.currentSegment).startNode;
        endNode = airwayStateS.baseSegmentS(airwayStateS.currentSegment).endNode;
        allNodesV = shortestpath(airwayGraph,startNode,endNode);
        airwayStateS.baseSegmentS(airwayStateS.currentSegment).allNodes = allNodesV;
        xyzM = airwayStateS.baseTreeS.nodeXyzM(allNodesV,:);
        if ishandle(airwayStateS.baseSegmentS(airwayStateS.currentSegment).hPlot)
            delete(airwayStateS.baseSegmentS(airwayStateS.currentSegment).hPlot)
        end
        airwayStateS.baseSegmentS(airwayStateS.currentSegment).hPlot = ...
            plot3(airwayStateS.hAxisBase,xyzM(:,1),xyzM(:,2),xyzM(:,3),'linewidth',5,...
            'color',[0,0,1,1],'hittest','off');
        airwayStateS.selectBaseStartNode = 1;
        set(airwayStateS.hBaseStart,'visible','off')
        airwaySegmentCallback('SEGMENT_CLICKED')
%         airwayStateS.basehStartNodes = [airwayStateS.basehStartNodes,...
%             plot3(x,y,z,'parent',hAxis,'marker','o','MarkerFaceColor','g',...
%             'markerSize',8,'hittest','off')];
    end
else
    if airwayStateS.selectFollowupStartNode  %size(airwayStateS.followupSegmentStartNodes,1) > size(airwayStateS.followupSegmentEndNodes,1)
        %airwayStateS.followupSegmentEndNodes(end+1,:) = [x,y,z];
        airwayStateS.followupSegmentS(airwayStateS.currentSegment).startNode = indMin(1);
        airwayStateS.selectFollowupStartNode = 0;
        set(airwayStateS.hFollowupStart,'XData',x,'YData',y,'zData',z,'visible','on')
        %set(airwayStateS.hFollowupEnd,'visible','off')
%         airwayStateS.followuphEndNodes = [airwayStateS.followuphEndNodes,...
%             plot3(x,y,z,'parent',hAxis,'marker','o','MarkerFaceColor','r',...
%             'markerSize',8,'hittest','off')];
    else
        %         airwayStateS.followupSegmentStartNodes(end+1,:) = [x,y,z];
        airwayStateS.followupSegmentS(airwayStateS.currentSegment).endNode = indMin(1);
        airwayGraph = airwayStateS.followupTreeS.airwayGraph;
        startNode = airwayStateS.followupSegmentS(airwayStateS.currentSegment).startNode;
        endNode = airwayStateS.followupSegmentS(airwayStateS.currentSegment).endNode;
        allNodesV = shortestpath(airwayGraph,startNode,endNode);
        airwayStateS.followupSegmentS(airwayStateS.currentSegment).allNodes = allNodesV;
        %set(airwayStateS.hFollowupEnd,'XData',x,'YData',y,'zData',z,'visible','on')
        xyzM = airwayStateS.followupTreeS.nodeXyzM(allNodesV,:);
        if ishandle(airwayStateS.followupSegmentS(airwayStateS.currentSegment).hPlot)
            delete(airwayStateS.followupSegmentS(airwayStateS.currentSegment).hPlot)
        end
        airwayStateS.followupSegmentS(airwayStateS.currentSegment).hPlot = ...
            plot3(airwayStateS.hAxisFollowup,xyzM(:,1),xyzM(:,2),xyzM(:,3),'linewidth',5,...
            'color',[0,0,1,1],'hittest','off');
        airwayStateS.selectFollowupStartNode = 1;
        set(airwayStateS.hFollowupStart,'visible','off')
        airwaySegmentCallback('SEGMENT_CLICKED')

%         airwayStateS.followuphStartNodes = [airwayStateS.followuphStartNodes,...
%             plot3(x,y,z,'parent',hAxis,'marker','o','MarkerFaceColor','g',...
%             'markerSize',8,'hittest','off')];
    end
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function xy = bud_radial ( n, xy )

%*****************************************************************************80
%
%% BUD_RADIAL sets the coordinates for a "radial" bud.
%
%  Discussion:
%
%    A radial bud is created by a generator on the boundary.  The initial
%    position of the bud should be created by copying the generator's
%    coordinates, reducing the radial component slightly.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of buds.
%
%    Input, real P_XY(2,N), the coordinates of the bud generators.
%
%    Output, real P_XY(2,N), the coordinates of the buds.
%
  r(1,1:n) = 0.999 * sqrt ( sum ( xy.^2, 1 ) );
  t(1,1:n) = atan2 ( xy(2,:), xy(1,:) );

  xy(1,:) = r(1,1:n) .* cos ( t(1,1:n) );
  xy(2,:) = r(1,1:n) .* sin ( t(1,1:n) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2014-2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym symsum (@var{f}, @var{n}, @var{a}, @var{b})
%% @defmethodx @@sym symsum (@var{f}, @var{n}, [@var{a} @var{b}])
%% @defmethodx @@sym symsum (@var{f}, @var{a}, @var{b})
%% @defmethodx @@sym symsum (@var{f}, [@var{a} @var{b}])
%% @defmethodx @@sym symsum (@var{f}, @var{n})
%% @defmethodx @@sym symsum (@var{f})
%% Symbolic summation.
%%
%% The sum of the expression @var{f} for @var{n} from @var{a} to
%% @var{b}.  When @var{n} is omitted it is determined using
%% @code{symvar} and defaults to @code{x} if @var{f} is constant. The
%% limits @var{a} and @var{b} default to 0 and @var{n} - 1
%% respectively.
%%
%% @example
%% @group
%% syms n m
%% symsum(1/n^2, n, 1, m)
%%   @result{} (sym) harmonic(m, 2)
%%
%% symsum(exp(2*n)/sin(n), n, 2*m, 6*m)
%%   @result{} (sym)
%%           6⋅m
%%           ____
%%           ╲
%%            ╲      2⋅n
%%             ╲    ℯ
%%             ╱   ──────
%%            ╱    sin(n)
%%           ╱
%%           ‾‾‾‾
%%         n = 2⋅m
%% @end group
%% @end example
%%
%% @seealso{@@sym/symprod, @@sym/sum}
%% @end defmethod

function S = symsum(f, n, a, b)

  if (nargin > 4)
    print_usage ();
  end

  idx1.type = '()';
  idx1.subs = {1};
  idx2.type = '()';
  idx2.subs = {2};

  if (nargin == 1)
    n = symvar(f, 1);
    if (isempty(n))
      n = sym('x');
    end
    a = sym(0);
    b = n - 1;
  elseif (nargin == 2) && (length(n) == 2)
    f = sym(f);
    %a = n(1);  % issue #17
    %b = n(2);
    a = subsref(n, idx1);
    b = subsref(n, idx2);
    n = symvar(f, 1);
    if (isempty(n))
      n = sym('x');
    end
  elseif (nargin == 2)
    f = sym(f);
    n = sym(n);
    a = sym(0);
    b = n - 1;
  elseif (nargin == 3) && (length(a) == 2)
    f = sym(f);
    n = sym(n);
    %b = a(2);  % issue #17
    %a = a(1);
    b = subsref(a, idx2);
    a = subsref(a, idx1);
  elseif (nargin == 3)
    f = sym(f);
    b = a;
    a = n;
    n = symvar(f, 1);
    if (isempty(n))
      n = sym('x');
    end
  else
    f = sym(f);
    n = sym(n);
    a = sym(a);
    b = sym(b);
  end

  cmd = { '(f, n, a, b) = _ins'
           'S = sp.summation(f, (n, a, b))'
           'return S,' };

  S = pycall_sympy__ (cmd, sym(f), sym(n), sym(a), sym(b));

end

%!error symsum (sym(1), 2, 3, 4, 5)

%!test
%! % finite sums
%! syms n
%! assert (isequal (symsum(n,n,1,10), 55))
%! assert(isa(symsum(n,n,1,10), 'sym'))
%! assert (isequal (symsum(n,n,sym(1),sym(10)), 55))
%! assert (isequal (symsum(n,n,sym(1),sym(10)), 55))
%! assert (isequal (symsum(1/n,n,1,10), sym(7381)/2520))

%!test
%! % negative limits
%! syms n
%! assert (isequal (symsum(n,n,-3,3), sym(0)))
%! assert (isequal (symsum(n,n,-3,0), sym(-6)))
%! assert (isequal (symsum(n,n,-3,-1), sym(-6)))

%!test
%! % one input
%! syms n
%! f = symsum (n);
%! g = n^2/2 - n/2;
%! assert (isequal (f, g))
%! f = symsum (2*n);
%! g = n^2 - n;
%! assert (isequal (f, g))

%!test
%! % constant input
%! f = symsum (sym(2));
%! syms x
%! g = 2*x;
%! assert (isequal (f, g))

%!test
%! % two inputs
%! syms n
%! f = symsum (2*n, n);
%! g = n^2 - n;
%! assert (isequal (f, g))

%!test
%! % two inputs, second is range
%! syms n
%! f = symsum (n, [1 6]);
%! g = 21;
%! assert (isequal (f, g))
%! f = symsum (n, [sym(1) 6]);
%! g = 21;
%! assert (isequal (f, g))
%! f = symsum (2*n, [1 6]);
%! g = 2*21;
%! assert (isequal (f, g))

%!test
%! % three inputs, last is range
%! syms n
%! f = symsum (2*n, n, [1 4]);
%! g = sym(20);
%! assert (isequal (f, g))
%! f = symsum (2*n, n, [sym(1) 4]);
%! g = sym(20);
%! assert (isequal (f, g))
%! f = symsum (2, n, [sym(1) 4]);
%! g = sym(8);
%! assert (isequal (f, g))

%!test
%! % three inputs, no range
%! syms n
%! f = symsum (2*n, 1, 4);
%! g = sym(20);
%! assert (isequal (f, g))
%! f = symsum (5, sym(1), 3);
%! g = sym(15);
%! assert (isequal (f, g))

%!test
%! % ok to use double's for arguments in infinite series
%! syms n oo
%! assert(isequal(symsum(1/n^2,n,1,oo), sym(pi)^2/6))
%! assert(isequal(symsum(1/n^2,n,1,inf), sym(pi)^2/6))

%!test
%! % should be oo because 1 is real but seems to be
%! % zoo/oo depending on sympy version
%! syms n oo
%! zoo = sym('zoo');
%! assert (isequal (symsum(1/n,n,1,oo), oo) || ...
%!         isequal (symsum(1/n,n,1,oo), zoo))

repo: octsympy


text:
clear

FV = sphere_tri('ico',2);

facecolor = [.7 .7 .7];
figure
Hp = patch('faces',FV.faces,'vertices',FV.vertices,...
    'facecolor',facecolor,'facealpha',1,'edgecolor',[.8 .8 .8]);
camlight('headlight','infinite'); daspect([1 1 1]); axis vis3d; axis off
material dull; rotate3d
hold on

f = 25;
vertex_index1 = FV.faces(f,1);
vertex_index2 = FV.faces(f,2);
vertex_index3 = FV.faces(f,3);
vertex1 = FV.vertices(vertex_index1,:);
vertex2 = FV.vertices(vertex_index2,:);
vertex3 = FV.vertices(vertex_index3,:);

plot3(vertex1(1),vertex1(2),vertex1(3),'ro')
plot3(vertex2(1),vertex2(2),vertex2(3),'go')
plot3(vertex3(1),vertex3(2),vertex3(3),'bo')

vertNormals = get(Hp,'vertexnormals');
quiver3(vertex1(1),vertex1(2),vertex1(3),...
    vertNormals(vertex_index1,1),vertNormals(vertex_index1,2),vertNormals(vertex_index1,3),0);

[faceNormals,faceNormalsUnit,centroids] = mesh_face_normals(FV);
quiver3(centroids(f,1),centroids(f,2),centroids(f,3),...
    faceNormals(f,1),faceNormals(f,2),faceNormals(f,3),0);

[MYvertNormals,MYvertNormalsUnit] = mesh_vertex_normals(FV);
MYvertNormalsMag = vector_magnitude(MYvertNormals);
quiver3(vertex1(1),vertex1(2),vertex1(3),...
    MYvertNormals(vertex_index1,1),MYvertNormals(vertex_index1,2),MYvertNormals(vertex_index1,3),0);

view(3)
rotate3d

vertNormalsMag = vector_magnitude(vertNormals);
vertNormalsUnit = vertNormals ./ repmat(vertNormalsMag,1,3);

% these comparisons are accurate to within 10^-12 or less, not sure why
% there is even these very small differences.
for v = 1:size(FV.vertices,1),
    dotprod(v,1) = dot( vertNormalsUnit(v,:), MYvertNormalsUnit(v,:) );
end

for v = 1:size(FV.vertices,1),
    crossprod(v,:) = cross( vertNormalsUnit(v,:), MYvertNormalsUnit(v,:) );
end

author: PatternRecognition

repo: OpenBMI


text:
function [ result, seed ] = sphere01_quad_mc ( f, h, seed, n )

%*****************************************************************************80
%
%% SPHERE01_QUAD_MC uses the Monte Carlo rule for sphere quadrature.
%
%  Discussion:
%
%    A number of points N are chosen at random on the sphere, with N
%    being determined so that, if the points were laid out on a regular
%    grid, the average spacing would be no more than H.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, function v = F ( n, x ), evaluates the integrand.
%
%    Input, real H, the maximum length of a side of the spherical
%    quadrilaterals.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Input, integer N, the number of points to use.
%
%    Output, real RESULT, the approximate integral.
%
  sphere_area = 4.0 * pi;

  [ x, seed ] = sphere01_sample_3d ( n, seed );

  v = f ( n, x )

  result = sphere_area * sum ( v(1:n) ) / n;

  return
end


text:
function [] = SLC_export(slcstack,slclist, Path, extention, InSAR_processor, reference_index)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file is part of TomoSAR.
%
%   TomoSAR is distributed in the hope that it will be useful,
%   but without warranty of any kind; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%   See the Apache License for more details.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author : Dinh Ho Tong Minh (INRAE) and Yen Nhi Ngo, Jan. 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if not(exist('InSAR_processor', 'var'))
     InSAR_processor = 'snap';
end

if not(exist('reference_index', 'var'))
     reference_index = 1;
end

[nlines,nwidths,n_slc] = size(slcstack);

real_index = 1:2:nwidths*2-1;
imag_index = 2:2:nwidths*2;
line_cpx = zeros(2*nwidths, 1);

for i = 1:n_slc
     switch InSAR_processor
        case 'snap' %
            filename = [Path,'/',num2str(slclist(i)),extention];
            fid = fopen(filename, 'wb', 'ieee-be');
        case 'isce'
            if i == reference_index
                filename = [Path,'/reference/','reference.slc',extention];
            else
                filename = [Path,'/',num2str(slclist(i)),'/','secondary.slc',extention];
            end
            fid = fopen(filename, 'wb');
        otherwise
            disp('not yet support')
    end

    data = squeeze(slcstack(:,:,i));
    for k=1:nlines
        line_cpx(real_index) = real(data(k,:));
        line_cpx(imag_index) = imag(data(k,:));
        line_count = fwrite(fid, line_cpx, 'float32');
    end
    fclose(fid);
end

author: DinhHoTongMinh

repo: TomoSAR


text:
function cvx_optval = norm_largest( x, k )

%NORM_LARGEST Sum of the k largest magnitudes of a vector.
%   NORM_LARGEST( X, k ) computes the 'largest-k' norm; that is, it computes
%   the sum of the magnitudes of the k largest elements in X. X must
%   be a vector, and k must be a real scalar.
%
%   Disciplined convex programming information:
%       NORM_LARGEST is convex and nonmonotonic of X, though it is monotonic for
%       positive values of X. So when used in CVX expressions, X must be affine,
%       monomial, or posynomial. k must be a real scalar constant.

error( nargchk( 2, 2, nargin ) );
if ~any( size( x ) ~= 1 ),
    error( 'The first argument must be a vector.' );
elseif ~isnumeric( k ) || ~isreal( k ) || length( k ) ~= 1,
    error( 'Third argument must be a scalar.' );
else
    cvx_optval = sum_largest( abs( x ), k );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function days = month_length_bahai ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_BAHAI returns the number of days in a Bahai month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of
%    days in the month.
%

%
%  Copy the input.
%
  m2 = m;
  y2 = y;
%
%  Check the input.
%
  [ y2, m2, ierror ] = ym_check_bahai ( y2, m2 );

  if ( ierror ~= 0 )
    days = 0;
    return
  end

  if ( m2 <= 18 || m2 == 20 )
    days = 19;
  elseif ( year_is_leap_bahai ( y2 ) )
    days = 5;
  else
    days = 4;
  end

  return
end


text:
function out = cfg_run_file_fplist(job)

% function out = cfg_run_file_fplist(job)
%
% Select files non-interactively using cfg_getfile('FPList',...) or
% cfg_getfile('FPListRec',...).
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_file_fplist.m 6918 2016-11-02 14:33:11Z guillaume $

rev = '$Rev: 6918 $'; %#ok

[out.files, out.dirs] = cfg_getfile(job.rec, job.dir{1}, job.filter);
if numel(out.files) < 2, plural = ''; else plural = 's'; end
if strcmpi(job.rec,'FPListRec'), rec = ' (and subdirectories)'; else rec = ''; end
fprintf('Found %d file%s matching filter ''%s'' in %s%s.\n',...
    numel(out.files),plural,job.filter,job.dir{1},rec);

author: spm


text:
% Copyright 2013 The MathWorks, Inc
%% Beamforming Example
% This example shows how to model antenna and target motion, and monitor
% array sidelobes while beamforming.
%% Model
% A 4-panel, ground based antenna array is mounted on a vehicle moving on a
% circular path at constant speed.
fc    = 3e9;
sAnt  = lowProfileArray('FrequencyRange',[2/3*fc 4/3*fc],'ViewArray',false);
[sSAnt, sSAntPlat]= splitSubArrays(sAnt);
for i=1:4,
    sCol{i}  = phased.Collector('Sensor',sSAnt{i},'OperatingFrequency',fc);
    sBF{i}   = phased.PhaseShiftBeamformer('OperatingFrequency',fc,...
                        'SensorArray',sSAnt{i},'DirectionSource','Input port',...
                        'WeightsOutputPort',true); %#ok<*SAGROW>
    sResp{i} = phased.ArrayResponse('SensorArray',sSAnt{i},'WeightsInputPort',true);
end
% A target is flying over the array.
sTgtMotion = phased.Platform('InitialPosition',[10e4; 0; 5e4],...
                             'Velocity',[-400; -200; -40]);
s=sin(2*pi*5*(0:.01:1))'; % Target Echo

%% Simulate
N = 500;
scanAz = -180:180;
for ii=1:N
    [AzEl,face,antpos,tgtpos,range] = updatePos(sSAntPlat,sTgtMotion,1/N,ii);
    viewTrajectories                               % Visualize trajectories
    y = []; resp = []; chan = 0;
    for k=find(face),
        x = step(sCol{k},s,AzEl);                  % Collect signal
        noise = 0.01*(randn(size(x)));             % Add noise
        chan = chan + 1;
        [y(:,chan),w] = step(sBF{k},x+noise,AzEl); % Beamform
        ang = [scanAz;ones(size(scanAz))*AzEl(2)];
        resp(:,chan) = step(sResp{k},fc,ang,w);    % Array response
    end
    y = mean(y,2);                                 % Combine active array outputs
    resp = mean(resp,2);
    plotBeam                                       % Visualize array response
end


text:
%WEIGHTEDMEDIANFILTER  Applies weighted median filter to an image
%
%     dst = cv.weightedMedianFilter(src, joint)
%     dst = cv.weightedMedianFilter(src, joint, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Source 8-bit or floating-point, 1-channel or 3-channel image.
% * __joint__ Joint 8-bit, 1-channel or 3-channel image.
%
% ## Output
% * __dst__ Destination image of the same size and type as `src`.
%
% ## Options
% * __Radius__ Radius of filtering kernel, should be a positive integer.
%   Default 7
% * __Sigma__ Filter range standard deviation for the joint image.
%   Default 25.5
% * __WeightType__ The type of weight definition. Specifies weight types of
%   weighted median filter, default 'EXP'. One of:
%   * __EXP__ `exp(-|I1-I2|^2 / (2*sigma^2))`
%   * __IV1__ `(|I1-I2| + sigma)^-1`
%   * __IV2__ `(|I1-I2|^2 + sigma^2)^-1`
%   * __COS__ `dot(I1,I2) / (|I1|*|I2|)`
%   * __JAC__ `(min(r1,r2) + min(g1,g2) + min(b1,b2)) / (max(r1,r2) + max(g1,g2) + max(b1,b2))`
%   * __OFF__ unweighted
% * __Mask__ A 0-1 mask that has the same size with `I`. This mask is used to
%   ignore the effect of some pixels. If the pixel value on mask is 0, the
%   pixel will be ignored when maintaining the joint-histogram. This is useful
%   for applications like optical flow occlusion handling. Not set by default.
%
% For more details about this implementation, please see [zhang2014100+].
%
% ## References
% [zhang2014100+]:
% > Qi Zhang, Li Xu, and Jiaya Jia. "100+ Times Faster Weighted Median Filter
% > (WMF)". In Computer Vision and Pattern Recognition (CVPR), 2014 IEEE
% > Conference on, pages 2830-2837. IEEE, 2014.
%
% See also: cv.medianBlur, cv.jointBilateralFilter
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef L2N < dagnn.Filter
% L2N layer that is L2-normalizing input vectors
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

  methods
    function outputs = forward(self, inputs, params)
      outputs{1} = vl_nnnormalizelp(inputs{1}, 'p', 2, 'epsilon', 1e-6);
    end

    function [derInputs, derParams] = backward(self, inputs, params, derOutputs)
      derInputs{1} = vl_nnnormalizelp(inputs{1}, derOutputs{1}, 'p', 2, 'epsilon', 1e-6);
      derParams = {} ;
    end

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

author: filipradenovic

repo: cnnimageretrieval


text:
function test_bug1902

% WALLTIME 02:00:00
% MEM 2gb
% DEPENDENCY ft_volumesegment ft_prepare_sourcemodel volumesmooth

mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/mri/nifti/single_subj_T1.nii'));
mri.coordsys = 'spm';

cfg = [];
tpm = ft_volumesegment(cfg,mri);  % gray, white, csf tissue prob. map

% there are some spm functions which change their own input

cfg = [];
cfg.mri = tpm;
grid01 = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.mri = tpm;
grid02 = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.mri = tpm;
grid03 = ft_prepare_sourcemodel(cfg);

grid01 = rmfield(grid01, 'cfg');
grid02 = rmfield(grid02, 'cfg');
grid03 = rmfield(grid03, 'cfg');

% make sure that they are the same, and not super-smoohted
isequal(grid01,grid02);
isequal(grid01,grid03);


text:
function dStates = rhs(~,States,Phase,P)

%This function is a wrapper for the controller and dynamics functions that
%makes them easier to call from ode45.

Actuators = Controller_PD(States,P.Control,Phase);

switch Phase
    case 'F'
        dStates = dynamics_flight(States,Actuators,P.Dyn);
    case 'D'
        dStates = dynamics_doubleStance(States,Actuators,P.Dyn);
    case 'S1'
        dStates = dynamics_singleStanceOne(States,Actuators,P.Dyn);
    case 'S2'
        dStates = dynamics_singleStanceTwo(States,Actuators,P.Dyn);
    otherwise
        error('Invalid phase, must be element of {''F'',''D'',''S1'',''S2''}');
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function path_name = PTKSaveAs(image_data, patient_name, path_name, is_secondary_capture, dicom_metadata, reporting)
    % MimSaveAs. Prompts the user for a filename and file type, and saves the image
    %
    %     Syntax
    %     ------
    %
    %         PTKSaveAs(image_data, patient_name, path_name, is_secondary_capture, dicom_metadata, reporting)
    %
    %             image_data      is a PTKImage (or PTKDicomImage) class containing the image
    %                             to be saved
    %             patient_name    specifies the patient name to be stored in the image (only
    %                             used when there is no metadata available in the image)
    %             path_name       specify the location to save the DICOM data. One 2D file
    %                             will be created for each image slice in the z direction.
    %                             Each file is numbered, starting from 0.
    %                             So if filename is 'MyImage.DCM' then the files will be
    %                             'MyImage0.DCM', 'MyImage1.DCM', etc.
    %             is_secondary_capture   true if the image is derived, false if the pixel data
    %                             directly corresponds to the original image pixel data
    %             dicom_metadata  a structure containing additional manufacturer tags
    %                             used to construct Dicom images
    %             reporting       an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if ~isa(image_data, 'PTKImage')
        error('Requires a PTKImage as input');
    end

    if nargin < 4 || isempty(is_secondary_capture)
        is_secondary_capture = false;
    end

    if nargin < 5 || isempty(dicom_metadata)
        dicom_metadata = PTKDicomMetadata();
    end

    if nargin < 6
        reporting = CoreReportingDefault();
    end

    path_name = MimSaveAs(image_data, patient_name, path_name, is_secondary_capture, dicom_metadata, reporting);
end

author: tomdoel

repo: pulmonarytoolkit


text:
figure1=gcf;
% Create textarrow
annotation(figure1,'textarrow',[0.281132075471698 0.226415094339623],...
    [0.845386533665835 0.800498753117207],'TextEdgeColor','none','FontSize',12,...
    'String',{'1st iteration'});

% Create textarrow
annotation(figure1,'textarrow',[0.516981132075472 0.518490566037736],...
    [0.389027431421446 0.324039900249375],'TextEdgeColor','none','FontSize',12,...
    'String',{'5th (final) iteration'});

% Create textarrow
annotation(figure1,'textarrow',[0.443396226415094 0.44188679245283],...
    [0.244389027431421 0.309127182044887],'TextEdgeColor','none','FontSize',12,...
    'String',{'4th iteration'});

% Create textarrow
annotation(figure1,'textarrow',[0.369811320754717 0.372452830188679],...
    [0.448877805486284 0.386334164588527],'TextEdgeColor','none','FontSize',12,...
    'String',{'3rd iteration'});

% Create textarrow
annotation(figure1,'textarrow',[0.284905660377358 0.286037735849057],...
    [0.321695760598504 0.416408977556109],'TextEdgeColor','none','FontSize',12,...
    'String',{'2nd iteration'});

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [q,Qv] = v2q(v)

% V2Q Rotaiton vector to quaternion conversion.
%   [Q,Qv] = V2Q(V) returns the quaternion Q correscponding to the rotation
%   encoded in rotation vector V, and the associated Jacobian Qv = dQ/dV.

if nargout == 1

    [a,u] = v2au(v);
    q = au2q(a,u);

else
    a = sqrt(dot(v,v));

    if isnumeric(a) && a < 1e-6

        % Use small signal approximation:
        q = [1-norm(v)^2/8
            v(:)/2];

        Qv = [-1/4*v(:)'
            0.5*eye(3)];

    else

        [a,u,Av,Uv] = v2au(v);
        [q,Qa,Qu] = au2q(a,u);
        Qv = Qa*Av + Qu*Uv;

    end
end

return

%%
syms r s t real
v = [r;s;t];
[q,Qv] = v2q(v)

simplify(Qv - jacobian(q,v))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [E]=dtiComputeApproximateEmbeddingVectors(fgfile, NfibersInA, method, npoints, kernelsigma, nvec)

%[E]=dtiComputeApproximateEmbeddingVectors(fgfile, [NfibersInA=1500', [method='pairwise_dist'], [npoints=15], [kernelsigma=30   ], [nvec=15])

%Computes embedded coordinates on a large fiberset
%Parameters: full fiber set (fgfile)
% NfibersInA: how many fibers form the "A" marix (full affinity matrix).
% Recommended/max my laptom with 2GB RAM can handle is 1500
%npoints (nodes in a resampled fiber)
%Note that your sample for A matrix will be picked from the first
%NfibersInA fibers. I is therefore crucial that the original fiberset is
%reshuffled using  fg = dtiShuffleFibers(fg) before it is passed into this
%function.

%Kernel is a parameter used for gaussian transformation from distances to
%affinities  (sigma of affinity=-distance^2/sigmasquare)

%ER 03/2008
%ER 11/2008 added default parameters

%TODO: = Specify recommended params
%= Add reshuffling by default
%= Add an option of using NfibersInA=size(fg.fibers, 2)

%%%%%%%%%%%%%%%%%%%%%%%%%
if(~exist('npoints','var')||isempty(npoints))
    npoints = 15;
end

if(~exist('kernelsigma','var')||isempty(kernelsigma))
    kernelsigma = 30; %30 for short distances, like clustering 1/2 brain; 60 for full brain distances
end

if(~exist('nvec','var')||isempty(nvec))
    nvec = 15;
end

if(~exist('method','var')||isempty(method))
    method = 'pairwise_dist';
end

if(~exist('NfibersInA','var')||isempty(NfibersInA))
    NfibersInA = 1500;
end

load(fgfile);
NumFibersTotal=size(fg.fibers, 1);
clear fg;
if (NumFibersTotal<NfibersInA)
NfibersInA=NumFibersTotal;
end

%Compute distances for A. Compute distances for B.
range1start=1; range2start=1; range1end=NfibersInA; range2end=NumFibersTotal;

compute_interfiber_distances(fgfile, npoints, method, range1start, range1end, range2start,  range2end, 1);
%If you want, use compute_interfiber_distances(fgfile, npoints, method, range1start, range1end, range2start,  range2end, 1);
%This will save some intermediate results on disk.

outfile=[prefix(fgfile) 'dist' num2str(range1start) 'to' num2str(range1end) 'vs' num2str(range2start) 'to' num2str(range2end) method '.mat'];

%Transform distances into proximities. Use kernel sigma=30; That makes sigmasquare=900.
load(outfile);

if (strcmp(method,'frenet'))

    for i=1:size(distmsr(:));
    distmsr(i)=log(distmsr(i));
    end

    distmsr(distmsr<-3)=-3; %rescale so that zero is the smalles distance
    distmsr= distmsr+3;
    distmsr(isinf(distmsr))=0;
end

%figure; hist(distmsr(:));
%kernelsigma=mean(distmsr(:))    ;

distmsr(1:NfibersInA, 1:NfibersInA)=(distmsr(1:NfibersInA, 1:NfibersInA)'+distmsr(1:NfibersInA, 1:NfibersInA))./2;
%just to fix assymmetry which shldnt be there  btw.

affinities=dtiComputeAffinitiesFromDistances(distmsr, kernelsigma);	clear distmsr;
clear fibergroup1 fibergroup2

%figure; hist(affinities(:));

%Perform estimation such that the size of B matrix is 5xSizeOfA which makes
%A about 20%

NSamples=min((NumFibersTotal-NfibersInA), NfibersInA*5);
%display(['Nsamples ' num2str(NSamples) ' NumFibersTotal ' num2str(NumFibersTotal) ' NFibersInA ' num2str(NfibersInA)]);
A=affinities(1:NfibersInA, 1:NfibersInA);
B=affinities(1:NfibersInA, NfibersInA+(1:NSamples));
clear affinities;

%Compute embedding vectors;
[E, embbasis]=dtiApproximateEmbeddingVectors(A, B, nvec);
display('Embedded space basis computed');

clear B;

maxlast=NfibersInA+NSamples;

while maxlast<NumFibersTotal

display(['Embedding fibers ' num2str(maxlast+1) ' to' num2str(min(maxlast+NSamples, NumFibersTotal))]);

load(outfile);
affinities=dtiComputeAffinitiesFromDistances(distmsr, kernelsigma);	clear distmsr;
clear fibergroup1 fibergroup2
S=affinities(1:NfibersInA, (maxlast+1):min(maxlast+NSamples, NumFibersTotal));
clear affinities;

E=[E; dtiNewDataOntoEmbeddingVectors(A, S, embbasis, nvec)];
maxlast=maxlast+NSamples;
end

%The next step will be computing actual clustering which we will put into a
%different function in case we wanted to stop here and use current embedded
%vectors for atlas creation purposes.

save([prefix(fgfile) 'EV' method], 'A', 'E', 'embbasis', 'npoints', 'kernelsigma', 'method');


text:
function [fevd_record,fevd_estimates]=panel5fevd(N,n,struct_irf_record,gamma_record,It,Bu,IRFperiods,FEVDband)

% run the Gibbs sampler to obtain FEVD draws
[fevd_record]=bear.fevd(struct_irf_record,gamma_record,It,Bu,N*n,IRFperiods,FEVDband);
% obtain point estimates and credibility intervals
[fevd_estimates]=bear.fevdestimates(fevd_record,N*n,IRFperiods,FEVDband);

author: european-central-bank

repo: BEAR-toolbox


text:
function [x,dx] = grid1d(a,b,nx)
%% Description:

% INPUTS
% [a,b]: x range, nx: the desided number of nodes.

% OUTPUTS
% uniformly distributed mesh grid for computation.
% x: vector x's coordinates, dx: cell size.

%% Compute cell size
dx = (b-a)/(nx-1);

%% Compute vector arrays
x = a:dx:b;

repo: Aero-matlab


text:
clc
clear
close all

%% Execute the configuration file to read parameters for data paths
addpath('config');
addpath(genpath('functions'));
configFile;

%% Starting parallel pooling (requires Parallel Processing Toolbox)
% This section takes a while to load for the first time
% To shutdown, run: delete(gcp('nocreate'));
if (isempty(gcp) && data_params.use_multithreads)
    parpool();
end

%% Read directories containing images
img_files1 = dir(strcat(data_params.path1,'*.png'));
img_files2 = dir(strcat(data_params.path2,'*.png'));
num_of_images = length(img_files1);

%% Read camera parameters
[P1, P2] = createCamProjectionMatrices(cam_params);

%% Read ground truth file if flag is true
if data_params.show_gt_flag
  ground_truth = load(data_params.gt_file);
  gt_x_max = max(ground_truth(:, end - 8));
  gt_x_min = min(ground_truth(:, end - 8));
  gt_z_max = max(ground_truth(:, end));
  gt_z_min = min(ground_truth(:, end));
end

%% Initialize variables for odometry
pos = [0;0;0];
Rpos = eye(3);

%% Start Algorithm
start = 0;
for t = 1 : num_of_images
    %% Read images for time instant t
    I2_l = imread([img_files1(t+1).folder, '/', img_files1(t).name]);
    I2_r = imread([img_files2(t+1).folder, '/', img_files2(t).name]);
    fprintf('Frame: %i\n', t);

    %% Bootstraping for initialization
    if (start == 0)
        vo_previous.pts1_l = computeFeatures(I2_l, vo_params.feature);
        vo_previous.pts1_r = computeFeatures(I2_r, vo_params.feature);
        start = 1;
        I1_l = I2_l;
        I1_r = I2_r;
        fprintf('\n---------------------------------\n');
        continue;
    end

    %% Implement SOFT for time instant t+1
    [R, tr, vo_previous] = visualSOFT(t, I1_l, I2_l, I1_r, I2_r, P1, P2, vo_params, vo_previous);

    %% Estimated pose relative to global frame at t = 0
    pos = pos + Rpos * tr';
    Rpos = R * Rpos;

    %% Prepare frames for next iteration
    I1_l = I2_l;
    I1_r = I2_r;

    %% Plot the odometry transformed data
    subplot(2, 2, [2, 4]);

    % Read ground truth pose if flag is true
    if data_params.show_gt_flag
      axis([gt_x_min gt_x_max gt_z_min gt_z_max])
      T = reshape(ground_truth(t, :), 4, 3)';
      pos_gt = T(:, 4);
      scatter(pos_gt(1), pos_gt(3), 'r', 'filled');
      hold on;
    end
    scatter( - pos(1), pos(3), 'b', 'filled');
    title(sprintf('Odometry plot at frame %d', t))
    xlabel('x-axis (in meters)');
    ylabel('z-axis (in meters)');

    if data_params.show_gt_flag
        legend('Ground Truth Pose', 'Estimated Pose')
    else
        legend('Estimated Pose')
    end

    %% Pause to visualize the plot
    pause(0.0001);
    fprintf('\n---------------------------------\n');
end

repo: Stereo-Odometry-SOFT


text:
function demo_RecPF

close all
clear all

addpath('utilities');
addpath('solver');
addpath('images');

% load images
Ims = cell(6,1);
for k = 1:6
    filename = ['ortho' int2str(k) '.jpg'];
    I = single(imread(filename))/255;
    if ndims(I) == 3
        Ims{k} = rgb2gray(I);
    end
end

% call tester
tester(Ims);

function tester(Ims)

%
% min aTV*TV(u) + aL1*||Phi*u||_1 + 0.5*||F_p*u - f_p||_2^2
%
aTV = 1.e-6;
aL1 = 0.e-6;

h = figure(1);
set(h,'units','normalized','outerposition',[0 .4 .9 .6]);

for nimg = 1:length(Ims)

    % load an image
    I = single(Ims{nimg});
    [m n] = size(I);
    N = m*n;

    % generate mask: two methods
    mth = 1;
    if mth == 1
        acc = 3.6;
        picks = selectPF(n/acc,n);
        K = length(picks);
    else
        Ls = 45;
        picks = fftshift(MRImask(n,Ls));
        picks = union(find(picks~=0),1);
        K = length(picks);
    end

    % image information display
    fprintf('Image %i: size %3i by %3i, N/M = %6.2f\n',nimg,m,n,N/K);
    snr(2*I,I);
    subplot(2,6,nimg); imshow(I,[]);
    title(['Original ' int2str(m) ' x ' int2str(n)]);
    drawnow;

    % add noise
    sigma = 0.01;
    noise = sigma*(randn(K,1) + sqrt(-1)*randn(K,1));

    % DWT and IDWT
    if exist('midwt','file')
        wav = daubcqf(2);
        Psi = @(x) midwt(x,wav);
        PsiT = @(x) mdwt(x,wav);
    elseif exist('wavedec2','file')
        Psi = @(x) Wavedb1Phi(x,1);
        PsiT = @(x) Wavedb1Phi(x,0);
    else
        Psi = @(x) x; PsiT = Psi;
    end

    % observation B
    FI = fft2(I);
    B = FI(picks) + noise;

    % run RecPF
    opts = [];
    t = cputime;
    U = RecPF(m,n,aTV,aL1,picks,B,PsiT,Psi,opts);
    t = cputime - t;

    % plot reconstructed images
    subplot(2,6,6+nimg);
    imshow(U,[]);
    title(sprintf('%4.2fdB, %4.2fs',snr(U,I),t));
    drawnow;

end


text:
function mbr = mbrot()
    i_list = 1:1024;
    mbr = nan(1024,1024,3);
    for i = i_list(:)'
        for j = i_list(:)'

            mbr(i,j,1) = red(i,j);
            mbr(i,j,2) = green(i,j);
            mbr(i,j,3) = blue(i,j);
        end
    end

    figure; image(mbr./max(mbr(:)));

    function out = red(i, j)
        a = 0;
        b = 0;
        d = 0;
        n = 0;
        while (a*a+(b*b)<4 && n<1024)
            b=2*a*b+(j-1)/5e4+.06;
            a=a*a-d+(i-1)/5e4+.34;
            d = b*b;
            n = n+1;
        end
        out = n/4;
    end

    function out = green(i, j)
        out = 4 * red(i, j);
    end

    function out = blue(i, j)
        out = 6 * red(i, j);
    end
end

author: goGPS-Project


text:
fid = fopen('results_20130124.token');
fid2 = fopen('flickr30k-train&val.txt','w');
load('train_val_test_split.mat');
test_index = find(set==3);
% uid sid tweet
tline = fgetl(fid);
raw_txt = [];
count = 1;
while ischar(tline)
    %skip test
    if ~isempty(find(ceil(count/5)==test_index))
        disp(count);
        count = count+1;
        tline = fgetl(fid);
        continue;
    end
    split_tline = strsplit(tline);
    s1 = numel(split_tline{1})+1;
    tline = tline(s1+1:end);
    fprintf(fid2,'%s\n',tline);
    %disp(tline)
    raw_txt{count,1} = tline;
    tline = fgetl(fid);
    count = count+1;
end
fclose(fid);
fclose(fid2);

author: layumi

repo: Image-Text-Embedding


text:
function [anat, inplanes] = GetAnatomyFromAnalyze2(fileName)
% [anat, inplanes] = GetAnatomyFromAnalyze2(fileName)
%
% Similar to getAnatomy but you pass in a 3d Analyze / nifti file that you
% already made. Requires SPM5
% ARW 061407
% SOD 090520

anat = [];

if ~exist(fileName, 'file');
    disp(sprintf('[%s]:File (%s) does not exist!',mfilename,fileName));
    return;
end

V=spm_vol(fileName);

[anat,XYZ] = spm_read_vols(V);

% Do a rot90 on each slice of the anat. That's just the way it is...
for thisSlice=1:size(anat,3)
    anat2(:,:,thisSlice)=rot90(anat(:,:,thisSlice),1);
end
anat=anat2;

% Use imatrix to recover slice spacing etc.
[inTrans]=spm_imatrix(V(1).mat);

sz = size(anat);
mm=abs(inTrans(7:9));

inplanes.FOV = mm(1)*sz(1);
inplanes.fullSize = sz(1:2);
inplanes.voxelSize = mm;
inplanes.spacing = 0;
% We already checked hdr.image.slquant==nList, not again -- Junjie
inplanes.nSlices = sz(3);

inplanes.examNum = V(1).descrip;
inplanes.crop = [];
inplanes.cropSize = [];

author: vistalab

repo: vistasoft


text:
function tri = triangulatePolygon(poly)
%TRIANGULATEPOLYGON Compute a triangulation of the polygon
%
%   TRI = triangulatePolygon(POLY)
%   Computes a triangulation TRI of the polygon defined by POLY
%   POLY contains the polygon vertices, as a Nv-by-2 array of double.
%   TRI is a Nt-by-3 array containing indices of vertices forming the
%   triangles.
%
%   Example
%     % creates a simple polygon and computes its Delaunay triangulation
%     poly = [0 0 ; 10 0;5 10;15 15;5 20;-5 10];
%     figure;drawPolygon(poly); axis equal
%     tri = triangulatePolygon(poly);
%     figure;
%     % patch('Faces', tri, 'Vertices', poly, 'facecolor', 'c');
%     drawMesh(poly, tri, 'facecolor', 'c');
%     axis equal
%
%     % Another example for which constrains were necessary
%     poly2 = [10 10;80 10; 140 20;30 20; 80 30; 140 30; 120 40;10 40];
%     tri2 = triangulatePolygon(poly2);
%     figure; drawMesh(poly2, tri2);
%     hold on, drawPolygon(poly2, 'linewidth', 2);
%     axis equal
%     axis([0 150 0 50])
%
%   See also
%     delaunayTriangulation, drawMesh, patch
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-11-25,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% compute constraints
nv = size(poly, 1);
cons = [(1:nv)' [2:nv 1]'];

if verLessThan('matlab', '8.1')
    % Code for versions before R2013a

    % delaunay triangulation
    dt = DelaunayTri(poly(:,1), poly(:, 2), cons); %#ok<DDELTRI>

    % find which triangles are contained in polygon
    centers = incenters(dt);
    inds = isPointInPolygon(centers, poly);

    % keep selected triangles
    tri = dt.Triangulation(inds, :);

else
    % Code for versions R2013a and later

    % delaunay triangulation
    % dt = DelaunayTri(poly(:,1), poly(:, 2), cons);
    dt = delaunayTriangulation(poly(:,1), poly(:, 2), cons);

    % find which triangles are contained in polygon
    % centers = incenters(dt);
    centers = incenter(dt);
    inds = isPointInPolygon(centers, poly);

    % keep selected triangles
    % tri = dt.Triangulation(inds, :);
    tri = dt.ConnectivityList(inds, :);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function N = congruenceblocks(exponent_m,exponent_p,options,csclasses)
%CONGRUENCEBLOCKS Partitions monomials based on sign symmetry
%
% V = CONGRUENCEBLOCKS(P)
%
% Input
%  V : Vector with SDPVAR objects
%  P : Scalar SDPVAR object
%
% Output
%  V : Cell with SDPVAR objects
%
% Example:
%
% sdpvar x y z
% p = 1+x*y+x^4+y^4+z+z^6;
% v = newtonmonoms(p);
% v = congruenceblocks(v,p);
% sdisplay(v{1}) % Even w.r.t (x,y)
% sdisplay(v{2}) % Odd w.r.t (x,y)
%
% See also NEWTONREDUCE, NEWTONMONOMS, CONSISTENT

sdpvarout = 0;
if isa(exponent_m,'sdpvar')
    z = depends(exponent_p);
    z = recover(unique([depends(exponent_p) depends(exponent_m)]));
    [exponent_p,p_base] = getexponentbase(exponent_p,z);
    [m,m_base] = getexponentbase(exponent_m,z);
    exponent_m = cell(1);exponent_m{1} = m;
    sdpvarout = 1;
end

if nargin < 3
    options.verbose = 0;
    options.sos.congruence = 2;
end

if nargin < 4
    csclasses = 1;
end

if ~isempty(exponent_m{1}) & options.sos.congruence>0 & ((size(exponent_p,2)<=16)  | options.sos.congruence==1)

    % **********************************************
    % DEFINE CONGRUENCE CLASSES
    % **********************************************
    if options.verbose>0;fprintf('Finding symmetries..............');end;
    n = size(exponent_p,2);
    t = cputime;
    switch options.sos.congruence
        case 1
            Htemp = eye(n);                 % CHEAP VERSION; ONLY CHECK IF IT IS EVEN WRT x_i
        case 2
            Htemp = dec2decbin(1:2^n-1,n)'; % ALL POSSIBLE COMBINATIONS
        otherwise
            error('sos.congruence should be 0, 1 or 2')
    end
    %try
        H = Htemp(:,find(~any(rem(exponent_p*Htemp,2),1))); % Find "even" rows
%     catch
%         i = [];
%         % Loop instead
%         for j = 1:size(Htemp,2)
%             if ~any(rem(exponent_p*Htemp(:,j),2))
%                 i = [i j];
%             end
%         end
%     end
    if isempty(H)
        N = exponent_m;
        if options.verbose>0;disp(['Found no symmetries (' num2str(cputime-t) 'sec)']);end
        return
    end

    t = cputime-t;
    if size(H,2)>=1
        if options.verbose>0
            if size(H,2)>1
                disp(['Found ' num2str(size(H,2)) ' symmetries  (' num2str(t) 'sec)']);
            else
                disp(['Found ' num2str(size(H,2)) ' symmetry  (' num2str(t) 'sec)']);
            end
        end
    else
        if options.verbose>0;disp(['Found no symmetries  (' num2str(t) 'sec)']);end;
    end

    % **********************************************
    % CLASSIFY MONMS ACCORDING TO CONGRUENCE CLASSES
    % **********************************************
    if size(H,2)>=1
        the_text = 'Partitioning using symmetry.....';
    end
    N = cell(0,1);
    for cs = 1:length(csclasses)
        [ur,j,k]=uniquesafe(mod(exponent_m{cs}*H,2),'rows');
        Ntemp = cell(size(ur,1),1);
        temp = [];
        for i = 1:length(k)
            Ntemp{k(i),1} = [Ntemp{k(i)};exponent_m{cs}(i,:)];
            temp = [temp;size(exponent_m{cs}(i,:),1)];
        end

        for i = 1:length(Ntemp)
            N{end+1,1} = Ntemp{i};
        end
    end
    % **********************************************
    % PRINT SOME RESULTS
    % **********************************************
    if size(H,2)>=1
        [uu,ii,oo] = uniquesafe(cellfun('prodofsize',N)/size(N{1},2));
        for i = 1:length(uu)
            n_this = length(find(oo==i));

            the_text = [the_text num2str(uu(i)) 'x' num2str(uu(i)) '(' num2str(n_this) ')' ' '];
        end
    end
    if options.verbose>0;;disp(the_text);end;
else
    N = exponent_m;
end

if sdpvarout
    for i = 1:length(N)
        N{i} = recovermonoms(N{i},z);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP