text:
function varargout = ma_OtherSpacecraftGUI(varargin)
% MA_OTHERSPACECRAFTGUI MATLAB code for ma_OtherSpacecraftGUI.fig
%      MA_OTHERSPACECRAFTGUI, by itself, creates a new MA_OTHERSPACECRAFTGUI or raises the existing
%      singleton*.
%
%      H = MA_OTHERSPACECRAFTGUI returns the handle to a new MA_OTHERSPACECRAFTGUI or the handle to
%      the existing singleton*.
%
%      MA_OTHERSPACECRAFTGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MA_OTHERSPACECRAFTGUI.M with the given input arguments.
%
%      MA_OTHERSPACECRAFTGUI('Property','Value',...) creates a new MA_OTHERSPACECRAFTGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ma_OtherSpacecraftGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ma_OtherSpacecraftGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ma_OtherSpacecraftGUI

% Last Modified by GUIDE v2.5 29-Mar-2019 17:03:34

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ma_OtherSpacecraftGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @ma_OtherSpacecraftGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ma_OtherSpacecraftGUI is made visible.
function ma_OtherSpacecraftGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ma_OtherSpacecraftGUI (see VARARGIN)

% Choose default command line output for ma_OtherSpacecraftGUI
handles.output = hObject;
handles.ma_MainGUI = varargin{1};
handles.ksptotMainGUI = varargin{2};

populateBodiesCombo(getappdata(handles.ma_MainGUI,'celBodyData'), handles.centralBodyCombo);
bodyComboValue = findValueFromComboBox('Kerbin', handles.centralBodyCombo);
if(isempty(bodyComboValue))
    bodyComboValue = 1;
end
set(handles.centralBodyCombo,'Value',bodyComboValue);

setupGUI(handles);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes ma_OtherSpacecraftGUI wait for user response (see UIRESUME)
uiwait(handles.ma_OtherSpacecraftGUI);

function setupGUI(handles)
    maData = getappdata(handles.ma_MainGUI,'ma_data');

    if(~isfield(maData.spacecraft,'otherSC'))
        maData.spacecraft.otherSC = cell(0,1);
        setappdata(handles.ma_MainGUI,'ma_data',maData);
    end

    if(isempty(maData.spacecraft.otherSC))
        set(handles.otherSCListbox,'String','');
        updateOSCListbox(handles);
    else
        oSCValue = get(handles.otherSCListbox,'Value');
        oSC = maData.spacecraft.otherSC{oSCValue};

        populateSpacecraft(handles,oSC);
        updateOSCListbox(handles);
        otherSCListbox_Callback(handles.otherSCListbox, [], handles);
    end

function populateSpacecraft(handles,oSC)
    set(handles.scNameText,'String',oSC.name);

    set(handles.epochText,'String',num2str(oSC.epoch));
    set(handles.smaText,'String',num2str(oSC.sma));
    set(handles.eccText,'String',num2str(oSC.ecc));
    set(handles.incText,'String',num2str(oSC.inc));
    set(handles.raanText,'String',num2str((oSC.raan)));
    set(handles.argText,'String',num2str(oSC.arg));
    set(handles.truText,'String',num2str(rad2deg(computeTrueAnomFromMean(deg2rad(oSC.mean), oSC.ecc))));
    set(handles.maxCommRangeText,'String',num2str(oSC.maxCommRange));

    bodyComboValue = findValueFromComboBox(oSC.parent, handles.centralBodyCombo);
    if(isempty(bodyComboValue))
        bodyComboValue = 1;
    end
    set(handles.centralBodyCombo,'Value',bodyComboValue);

    colorComboStr = getStringFromLineSpecColor(oSC.color);
	value = findValueFromComboBox(colorComboStr, handles.otherSCColor);
    set(handles.otherSCColor,'Value',value);

    if(isfield(oSC,'linestyle'))
        lineStyle = getLineStyleFromString(oSC.linestyle);

        value = findValueFromComboBox(lineStyle, handles.otherSCLineStyle);
        set(handles.otherSCLineStyle,'Value',value);
    else
        set(handles.otherSCLineStyle,'Value',2);
    end

    if(isfield(oSC,'lineWidth'))
        lineWidth = oSC.lineWidth;

        contents = handles.lineWidthCombo.String;
        contentsDouble = str2double(contents);
        ind = find(contentsDouble == lineWidth);
        set(handles.lineWidthCombo,'Value',ind);
    else
        set(handles.lineWidthCombo,'Value',2);
    end

% --- Outputs from this function are returned to the command line.
function varargout = ma_OtherSpacecraftGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    varargout{1} = [];

% --- Executes on selection change in otherSCListbox.
function otherSCListbox_Callback(hObject, eventdata, handles)
% hObject    handle to otherSCListbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns otherSCListbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from otherSCListbox
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    otherSC = maData.spacecraft.otherSC;
    oSCValue = get(handles.otherSCListbox,'Value');

    if(~isempty(otherSC) && ~isempty(oSCValue))
        oSC = otherSC{oSCValue};
        populateSpacecraft(handles,oSC);

        if(strcmpi(get(handles.ma_OtherSpacecraftGUI,'SelectionType'),'open'))
            celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
            bodyInfo = getBodyInfoByNumber(oSC.parentID, celBodyData);
            tru = computeTrueAnomFromMean(deg2rad(oSC.mean), oSC.ecc);
            [rVect,vVect]=getStatefromKepler(oSC.sma, oSC.ecc, deg2rad(oSC.inc), deg2rad(oSC.raan), deg2rad(oSC.arg), tru, bodyInfo.gm);

            UT = oSC.epoch;

            rX = rVect(1);
            rY = rVect(2);
            rZ = rVect(3);

            vX = vVect(1);
            vY = vVect(2);
            vZ = vVect(3);

            state = [UT, rX, rY, rZ, vX, vY, vZ, oSC.parentID, 0 0 0 0 1];
            viewSpacecraftStatePopupGUI(maData.spacecraft.propellant.names, state, 1, celBodyData);
        end
    end

% --- Executes during object creation, after setting all properties.
function otherSCListbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to otherSCListbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in updateSCButton.
function updateSCButton_Callback(hObject, eventdata, handles)
% hObject    handle to updateSCButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);
    if(~isempty(errMsg))
        msgbox(errMsg,'Errors were found while updating a spacecraft.','error');
        return;
    end

    maData = getappdata(handles.ma_MainGUI,'ma_data');
    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');

    if(isempty(get(handles.otherSCListbox,'String')))
        return;
    end
    oscLBValue = get(handles.otherSCListbox,'Value');

    oldOtherSC = maData.spacecraft.otherSC{oscLBValue};
    oldID = oldOtherSC.id;

    ecc = str2double(get(handles.eccText,'String'));
    tru = deg2rad(str2double(get(handles.truText,'String')));
    mean = rad2deg(computeMeanFromTrueAnom(tru, ecc));

    contents = cellstr(get(handles.otherSCColor,'String'));
    oSCColorStr = contents{get(handles.otherSCColor,'Value')};
    oscColor = getLineSpecColorFromString(oSCColorStr);

    contents = cellstr(get(handles.otherSCLineStyle,'String'));
	oscLineStyle = contents{get(handles.otherSCLineStyle,'Value')};
    lineStyle = getLineStyleStrFromText(oscLineStyle);

    contents = cellstr(get(handles.lineWidthCombo,'String'));
	oscLineWidth = str2double(contents{get(handles.lineWidthCombo,'Value')});

    otherSC = struct();

    otherSC.name = get(handles.scNameText,'String');

    otherSC.epoch = str2double(get(handles.epochText,'String'));
    otherSC.sma = str2double(get(handles.smaText,'String'));
    otherSC.ecc = ecc;
    otherSC.inc = str2double(get(handles.incText,'String'));
    otherSC.raan = str2double(get(handles.raanText,'String'));
    otherSC.arg = str2double(get(handles.argText,'String'));
    otherSC.mean = mean;
    otherSC.maxCommRange = str2double(get(handles.maxCommRangeText,'String'));

    contents = cellstr(get(handles.centralBodyCombo,'String'));
    cbName = strtrim(contents{get(handles.centralBodyCombo,'Value')});
    bodyInfo = celBodyData.(strtrim(lower(cbName)));
    otherSC.parent = cbName;
    otherSC.parentID = bodyInfo.id;

    otherSC.id = oldID;
    otherSC.color = oscColor;
    otherSC.linestyle = lineStyle;
    otherSC.lineWidth = oscLineWidth;

    maData.spacecraft.otherSC{oscLBValue} = otherSC;
    setappdata(handles.ma_MainGUI,'ma_data',maData);

    updateOSCListbox(handles);

% --- Executes on button press in addSCButton.
function addSCButton_Callback(hObject, eventdata, handles)
% hObject    handle to addSCButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');

    oscColor = 'r';
    lineStyle = '-';

    otherSC = struct();

    otherSC.name = 'New Spacecraft';

    otherSC.epoch = 0.0;
    otherSC.sma = 7000;
    otherSC.ecc = 0.0;
    otherSC.inc = 0.0;
    otherSC.raan = 0.0;
    otherSC.arg = 0.0;
    otherSC.mean = 0.0;
    otherSC.maxCommRange = Inf;

    contents = cellstr(get(handles.centralBodyCombo,'String'));
    cbName = contents{1};
    bodyInfo = celBodyData.(strtrim(lower(cbName)));
    otherSC.parent = cbName;
    otherSC.parentID = bodyInfo.id;

    otherSC.id = rand();
    otherSC.color = oscColor;
    otherSC.linestyle = lineStyle;
    otherSC.lineWidth = 0.5;

    maData.spacecraft.otherSC{end+1} = otherSC;
    setappdata(handles.ma_MainGUI,'ma_data',maData);

    updateOSCListbox(handles);
    set(handles.otherSCListbox,'Value',length(maData.spacecraft.otherSC));
    otherSCListbox_Callback(handles.otherSCListbox, [], handles);

function updateOSCListbox(handles)
	maData = getappdata(handles.ma_MainGUI,'ma_data');

    if(isempty(maData.spacecraft.otherSC))
        enableUI('off', maData, handles);
    else
        enableUI('on', maData, handles);
    end

    otherSC = maData.spacecraft.otherSC;
    otherSCStr = cell(0,1);
    for(i=1:length(otherSC)) %#ok<*NO4LP>
        otherSCStr{end+1} = otherSC{i}.name; %#ok<AGROW>
    end

    if(isempty(otherSC))
        enableUI('off', maData, handles)
    else
        enableUI('on', maData, handles)
    end

	if(isempty(get(handles.otherSCListbox,'Value')))
        set(handles.otherSCListbox,'Value',1);
    elseif(get(handles.otherSCListbox,'Value') > length(otherSC))
        set(handles.otherSCListbox,'Value',length(otherSC));
    elseif(get(handles.otherSCListbox,'Value') == 0 && ~isempty(otherSC))
        set(handles.otherSCListbox,'Value',1);
	end

    set(handles.otherSCListbox,'String',otherSCStr);

% --- Executes on button press in removeSCButton.
function removeSCButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeSCButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    oSCValue = get(handles.otherSCListbox,'Value');

    otherSC = maData.spacecraft.otherSC;

    script = maData.script;
    for(i=1:length(script))
        event = script{i};

        if(strcmpi(event.type,'Docking'))
            if(event.oScId == otherSC{oSCValue}.id)
                warnstr = sprintf('Cannot delete other spacecraft "%s".  This spacecraft is in use by event %u "%s".', otherSC{oSCValue}.name, i, event.name);
                warndlg(warnstr,'Cannot Delete Spacecraft', 'modal');

                return;
            end
        end
    end

    otherSC(oSCValue) = [];
    maData.spacecraft.otherSC = otherSC;
    setappdata(handles.ma_MainGUI,'ma_data',maData);

    updateOSCListbox(handles);
    otherSCListbox_Callback(handles.otherSCListbox, [], handles);

function enableUI(enable, maData, handles)
    set(handles.otherSCListbox,'Enable',enable)
    set(handles.scNameText,'Enable',enable)
    set(handles.otherSCColor,'Enable',enable)
    set(handles.otherSCLineStyle,'Enable',enable)
    set(handles.centralBodyCombo,'Enable',enable)
    set(handles.epochText,'Enable',enable)
    set(handles.smaText,'Enable',enable)
    set(handles.eccText,'Enable',enable)
    set(handles.incText,'Enable',enable)
    set(handles.raanText,'Enable',enable)
    set(handles.scNameText,'Enable',enable)
    set(handles.argText,'Enable',enable)
    set(handles.truText,'Enable',enable)
    set(handles.maxCommRangeText,'Enable',enable)

	if(isempty(maData.spacecraft.otherSC))
        set(handles.removeSCButton,'Enable','off')
    else
        set(handles.removeSCButton,'Enable','on')
	end

function scNameText_Callback(hObject, eventdata, handles)
% hObject    handle to scNameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of scNameText as text
%        str2double(get(hObject,'String')) returns contents of scNameText as a double
    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function scNameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to scNameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in centralBodyCombo.
function centralBodyCombo_Callback(hObject, eventdata, handles)
% hObject    handle to centralBodyCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns centralBodyCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from centralBodyCombo
    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function centralBodyCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to centralBodyCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function epochText_Callback(hObject, eventdata, handles)
% hObject    handle to epochText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of epochText as text
%        str2double(get(hObject,'String')) returns contents of epochText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function epochText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to epochText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function smaText_Callback(hObject, eventdata, handles)
% hObject    handle to smaText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of smaText as text
%        str2double(get(hObject,'String')) returns contents of smaText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function smaText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to smaText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function eccText_Callback(hObject, eventdata, handles)
% hObject    handle to eccText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of eccText as text
%        str2double(get(hObject,'String')) returns contents of eccText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function eccText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eccText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function incText_Callback(hObject, eventdata, handles)
% hObject    handle to incText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of incText as text
%        str2double(get(hObject,'String')) returns contents of incText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function incText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to incText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function raanText_Callback(hObject, eventdata, handles)
% hObject    handle to raanText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of raanText as text
%        str2double(get(hObject,'String')) returns contents of raanText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function raanText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to raanText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function argText_Callback(hObject, eventdata, handles)
% hObject    handle to argText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of argText as text
%        str2double(get(hObject,'String')) returns contents of argText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function argText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to argText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function truText_Callback(hObject, eventdata, handles)
% hObject    handle to truText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of truText as text
%        str2double(get(hObject,'String')) returns contents of truText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function truText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to truText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function errMsg = validateInputs(handles)
    errMsg = {};

    eEcc = str2double(get(handles.eccText,'String'));
    enteredStr = get(handles.eccText,'String');
    numberName = 'Orbit Eccentricity';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(eEcc, numberName, lb, ub, isInt, errMsg, enteredStr);

    if(not(isempty(errMsg)))
        msgbox(errMsg,'Errors were found while creating the spacecraft orbit.','error');
        return;
    end

    epoch = str2double(get(handles.epochText,'String'));
    enteredStr = get(handles.epochText,'String');
    numberName = 'Orbit Epoch';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(epoch, numberName, lb, ub, isInt, errMsg, enteredStr);

    eSMA = str2double(get(handles.smaText,'String'));
    enteredStr = get(handles.smaText,'String');
    if(eEcc < 1.0)
        lb = 1;
        ub = Inf;
        numberName = 'Eccentric Orbit Semi-Major Axis';
    else
        lb = -Inf;
        ub = -1;
        numberName = 'Hyperbolic Orbit Semi-Major Axis';
    end
    isInt = false;
    errMsg = validateNumber(eSMA, numberName, lb, ub, isInt, errMsg, enteredStr);

    eInc = str2double(get(handles.incText,'String'));
    enteredStr = get(handles.incText,'String');
    numberName = 'Orbit Inclination';
    lb = 0;
    ub = 180;
    isInt = false;
    errMsg = validateNumber(eInc, numberName, lb, ub, isInt, errMsg, enteredStr);

    eRAAN = str2double(get(handles.raanText,'String'));
    enteredStr = get(handles.raanText,'String');
    numberName = 'Orbit RAAN';
    lb = 0;
    ub = 360;
    isInt = false;
    errMsg = validateNumber(eRAAN, numberName, lb, ub, isInt, errMsg, enteredStr);

    eArg = str2double(get(handles.argText,'String'));
    enteredStr = get(handles.argText,'String');
    numberName = 'Orbit Argument of Periapse';
    lb = 0;
    ub = 360;
    isInt = false;
    errMsg = validateNumber(eArg, numberName, lb, ub, isInt, errMsg, enteredStr);

    eTru = str2double(get(handles.truText,'String'));
    enteredStr = get(handles.truText,'String');
    numberName = 'Orbit True Anomaly';
    lb = -180;
    ub = 360;
    isInt = false;
    errMsg = validateNumber(eTru, numberName, lb, ub, isInt, errMsg, enteredStr);

    maxCommRange = str2double(get(handles.maxCommRangeText,'String'));
    enteredStr = get(handles.maxCommRangeText,'String');
    numberName = 'Max Comm. Range';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(maxCommRange, numberName, lb, ub, isInt, errMsg, enteredStr);

% --------------------------------------------------------------------
function getOrbitFromSFSFileContextMenu_Callback(hObject, eventdata, handles)
% hObject    handle to getOrbitFromSFSFileContextMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    [refBodyID, scName] = orbitPanelGetOrbitFromSFSContextCallBack(handles.ksptotMainGUI, handles.smaText, handles.eccText, handles.incText, handles.raanText, handles.argText, handles.truText, handles.epochText);
    tru = computeTrueAnomFromMean(deg2rad(str2double(get(handles.truText,'String'))), str2double(get(handles.eccText,'String')));
    set(handles.truText,'String',num2str(rad2deg(tru)));

    if(~isempty(refBodyID) && isnumeric(refBodyID))
        celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
        bodyInfo = getBodyInfoByNumber(refBodyID, celBodyData);
        value = findValueFromComboBox(bodyInfo.name, handles.centralBodyCombo);
        set(handles.centralBodyCombo,'Value',value);
    end

    if(~isempty(scName))
        if(strcmpi(handles.scNameText.String,'New Spacecraft'))
            set(handles.scNameText,'String',scName);
        else
            qstring = sprintf('Do you wish overwrite the name of the spacecraft with the name of the imported spacecraft (%s)?', scName);
            title = 'Update S/C Name?';
            button = questdlg(qstring,title,'Yes','No','No');

            if(strcmpi(button,'Yes'))
                set(handles.scNameText,'String',scName);
            end
        end
    end

    updateSCButton_Callback([], [], handles);

% --------------------------------------------------------------------
function getOrbitFromKSPTOTConnectContextMenu_Callback(hObject, eventdata, handles)
% hObject    handle to getOrbitFromKSPTOTConnectContextMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    [refBodyID, scName] = orbitPanelGetOrbitFromKSPTOTConnectCallBack(handles.smaText, handles.eccText, handles.incText, handles.raanText, handles.argText, handles.truText, handles.epochText);
    tru = computeTrueAnomFromMean(deg2rad(str2double(get(handles.truText,'String'))), str2double(get(handles.eccText,'String')));
    set(handles.truText,'String',num2str(rad2deg(tru)));

    if(~isempty(refBodyID) && isnumeric(refBodyID))
        celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
        bodyInfo = getBodyInfoByNumber(refBodyID, celBodyData);
        value = findValueFromComboBox(bodyInfo.name, handles.centralBodyCombo);
        set(handles.centralBodyCombo,'Value',value);
    end

    if(~isempty(scName))
        if(strcmpi(handles.scNameText.String,'New Spacecraft'))
            set(handles.scNameText,'String',scName);
        else
            qstring = sprintf('Do you wish overwrite the name of the spacecraft with the name of the imported spacecraft (%s)?', scName);
            title = 'Update S/C Name?';
            button = questdlg(qstring,title,'Yes','No','No');

            if(strcmpi(button,'Yes'))
                set(handles.scNameText,'String',scName);
            end
        end
    end

    updateSCButton_Callback([], [], handles);

% --- Executes on selection change in otherSCColor.
function otherSCColor_Callback(hObject, eventdata, handles)
% hObject    handle to otherSCColor (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns otherSCColor contents as cell array
%        contents{get(hObject,'Value')} returns selected item from otherSCColor
    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function otherSCColor_CreateFcn(hObject, eventdata, handles)
% hObject    handle to otherSCColor (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function copyOrbitToClipboardMenu_Callback(hObject, eventdata, handles)
% hObject    handle to copyOrbitToClipboardMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    oscLBValue = get(handles.otherSCListbox,'Value');
    oldOtherSC = maData.spacecraft.otherSC{oscLBValue};

    copyOrbitToClipboardFromText(handles.epochText, handles.smaText, handles.eccText, ...
                                 handles.incText, handles.raanText, handles.argText, ...
                                 handles.truText, true, oldOtherSC.parentID);

% --------------------------------------------------------------------
function pasteOrbitFromClipboardMenu_Callback(hObject, eventdata, handles)
% hObject    handle to pasteOrbitFromClipboardMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
    pasteOrbitFromClipboard(handles.epochText, handles.smaText, handles.eccText, ...
                                 handles.incText, handles.raanText, handles.argText, ...
                                 handles.truText, true, handles.centralBodyCombo, celBodyData);

    updateSCButton_Callback([], [], handles);

% --- Executes on key press with focus on ma_OtherSpacecraftGUI or any of its controls.
function ma_OtherSpacecraftGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to ma_OtherSpacecraftGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'escape'
            close(handles.ma_OtherSpacecraftGUI);
    end

% --- Executes on selection change in otherSCLineStyle.
function otherSCLineStyle_Callback(hObject, eventdata, handles)
% hObject    handle to otherSCLineStyle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns otherSCLineStyle contents as cell array
%        contents{get(hObject,'Value')} returns selected item from otherSCLineStyle
    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function otherSCLineStyle_CreateFcn(hObject, eventdata, handles)
% hObject    handle to otherSCLineStyle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function getOrbitFromKSPActiveVesselMenu_Callback(hObject, eventdata, handles)
% hObject    handle to getOrbitFromKSPActiveVesselMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    [refBodyID, scName] = orbitPanelGetOrbitFromKSPTOTConnectActiveVesselCallBack(handles.smaText, handles.eccText, handles.incText, handles.raanText, handles.argText, handles.truText, handles.epochText);
    tru = computeTrueAnomFromMean(deg2rad(str2double(get(handles.truText,'String'))), str2double(get(handles.eccText,'String')));
    set(handles.truText,'String',fullAccNum2Str(rad2deg(tru)));

    if(~isempty(refBodyID) && isnumeric(refBodyID))
        celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
        bodyInfo = getBodyInfoByNumber(refBodyID, celBodyData);
        value = findValueFromComboBox(bodyInfo.name, handles.centralBodyCombo);
        set(handles.centralBodyCombo,'Value',value);
    end

    if(~isempty(scName))
        if(strcmpi(handles.scNameText.String,'New Spacecraft'))
            set(handles.scNameText,'String',scName);
        else
            qstring = sprintf('Do you wish overwrite the name of the spacecraft with the name of the imported spacecraft (%s)?', scName);
            title = 'Update S/C Name?';
            button = questdlg(qstring,title,'Yes','No','No');

            if(strcmpi(button,'Yes'))
                set(handles.scNameText,'String',scName);
            end
        end
    end

    updateSCButton_Callback([], [], handles);

function maxCommRangeText_Callback(hObject, eventdata, handles)
% hObject    handle to maxCommRangeText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of maxCommRangeText as text
%        str2double(get(hObject,'String')) returns contents of maxCommRangeText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function maxCommRangeText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to maxCommRangeText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in lineWidthCombo.
function lineWidthCombo_Callback(hObject, eventdata, handles)
% hObject    handle to lineWidthCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns lineWidthCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from lineWidthCombo
    updateSCButton_Callback([], [], handles);

% --- Executes during object creation, after setting all properties.
function lineWidthCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lineWidthCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
function stage_transition_type

before_contact = ocl.Stage([], @before_contact_vars, @before_contact_ode, 'N', 3, 'd', 2);
after_contact = ocl.Stage(1, @after_contact_vars, @after_contact_ode, ...
  @after_contact_cost, 'N', 5, 'd', 2);

before_contact.setInitialBounds('s', 1);
before_contact.setInitialBounds('v', 0);
before_contact.setEndBounds('s', 0);

before_contact.setBounds('energy_loss', 2);

after_contact.setBounds('s', 0, 2);

after_contact.setEndBounds('s', 1);

ocp = ocl.MultiStageProblem({before_contact, after_contact}, {@stage_transition},...
  'transition_type', 2);

[sol,times] = ocp.solve(ocp.getInitialGuess());

stage_1 = sol{1}.states(:,[1,end]).value;
ocl.utils.assertAlmostEqual(stage_1, [1 0;0 -4.47214], 'Feature test bouncing ball failed.');

end

function before_contact_vars(sh)
sh.addState('s');
sh.addState('v');

sh.addParameter('energy_loss');
end

function before_contact_ode(sh,x,~,~,~)
sh.setODE('s', x.v);
sh.setODE('v', -10);
end

function after_contact_vars(sh)
sh.addState('s');
sh.addState('v');
sh.addControl('F');
end

function after_contact_ode(dh,x,~,u,~)
dh.setODE('s', x.v);
dh.setODE('v', -10 + 10*u.F);

ocl.utils.assert(isempty(dh.userdata()));
end

function after_contact_cost(ch,~,~,u,~)
ch.add( u.F^2 );
end

function stage_transition(ch, x0_cur, xF_prev, p_cur, p_prev)
% x0 current stage
% xF previous stage
ch.add(x0_cur.s, '==', xF_prev.s);
ch.add(x0_cur.v, '==', -xF_prev.v / p_prev.energy_loss);
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function value = r8vec_i4vec_dot_product ( n, r8vec, i4vec )

%*****************************************************************************80
%
%% R8VEC_I4VEC_DOT_PRODUCT finds the dot product of an R8VEC and an I4VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    An I4VEC is a vector of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vectors.
%
%    Input, real R8VEC(N), the first vector.
%
%    Input, integer I4VEC(N), the second vector.
%
%    Output, real R8VEC_I4VEC_DOT_PRODUCT, the dot product.
%
  value = r8vec * i4vec';

  return
end

author: johannesgerer

repo: jburkardt-m


%==============================================================================
% This code is part of the Finite Element Method app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/FAIRFEM
%==============================================================================
%
% function [Sc,dS,d2S] = elasticFEM(uc,xc,tri,varargin)
%
% linear elastic regularizer for tetrahedral finite element discretization
% with linear basis functions
%
% S(u) = \int || B * u  ||^2 dx
%
% where B is the navier lame operator
%
% Input:
%   uc    - coefficients of current displacement
%   yRef  - coefficients of reference transformation, yc = uc + yRef
%   Mesh  - representation of mesh
%
% Output:
%    Sc    - hyperelastic energy
%    dS    - first derivative
%    d2S   - Hessian
%
% see also elastic.m (staggered grid version)
%==============================================================================

function [Sc,dS,d2S] = elasticFEM(uc,~,Mesh,varargin)

if nargin == 0, help(mfilename); runMinimalExample;return;end

persistent A alphaOld
if ~exist('A','var'),        A        = []; end;
if ~exist('alphaOld','var'), alphaOld = []; end;

% default parameter
matrixFree  = 0;
alpha       = 1;
mu          = 1;
lambda      = 0;

for k=1:2:length(varargin), % overwrites default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

if not(matrixFree), % matrix-based
    build = isempty(alphaOld) || any(alphaOld~=[alpha mu lambda]) ||...
        isempty(A) || size(A,2)~=numel(uc);
    if build,
        alphaOld = [alpha mu lambda];
        A = getElasticMatrixFEM(Mesh,mu,lambda);
        A = alpha * (A'* A);
    end;
    dS  = uc'*A;
    Sc  = 0.5*dS*uc;
    d2S = A;
else % matrix-free
    d2S.regularizer = regularizer;
    d2S.alpha  = alpha;
    d2S.By     = @(u,Mesh) elasticOperator(u,Mesh,mu,lambda,'By');
    d2S.BTy    = @(u,Mesh) elasticOperator(u,Mesh,mu,lambda,'BTy');
    d2S.B      = @(Mesh) getElasticMatrixFEM(Mesh,mu,lambda);
    d2S.diag   = @(Mesh) getDiag(Mesh,mu,lambda,alpha);
    d2S.solver = @FEMMultiGridsolveGN;
    d2S.d2S  = @(uc,Mesh) ...
        alpha * ...
        elasticOperator(elasticOperator(uc,Mesh,mu,lambda,'By'),Mesh,mu,lambda,'BTy');

    dS   = d2S.d2S(uc,Mesh)';
    Sc   = .5*dS*uc;
end

function By = elasticOperator(u,Mesh,mu,lambda,flag)
vol  = sqrt(Mesh.vol);
a    = sdiag(sqrt(mu)*vol);
b    = sdiag(sqrt(mu+lambda)*vol);
flag = [num2str(Mesh.dim) 'D-' flag];
switch flag
    case '2D-By'
        u  = reshape(u,[],2);
        dx1 = Mesh.mfdx1.D ; dx2 = Mesh.mfdx2.D;

        By = [...
            a*dx1(u(:,1));...
            a*dx2(u(:,1));...
            a*dx1(u(:,2));...
            a*dx2(u(:,2));...
            b*dx1(u(:,1))+b*dx2(u(:,2))...
            ];
    case '2D-BTy'
        u = reshape(u,[],5);
        dx1 = Mesh.mfdx1.Dadj ; dx2 = Mesh.mfdx2.Dadj;

        By = [dx1(a*u(:,1))+dx2(a*u(:,2))+dx1(b*u(:,5));...
            dx1(a*u(:,3))+dx2(a*u(:,4))+dx2(b*u(:,5))];
    case '3D-By'
        u  = reshape(u,[],3);
        dx1 = Mesh.mfdx1.D ; dx2 = Mesh.mfdx2.D; dx3 = Mesh.mfdx3.D;

        By = [...
                a*dx1(u(:,1));...
                a*dx2(u(:,1));...
                a*dx3(u(:,1));...
                a*dx1(u(:,2));...
                a*dx2(u(:,2));...
                a*dx3(u(:,2));...
                a*dx1(u(:,3));...
                a*dx2(u(:,3));...
                a*dx3(u(:,3));...
                b*dx1(u(:,1))+b*dx2(u(:,2))+b*dx3(u(:,3))...
            ];
    case '3D-BTy'
        u = reshape(u,[],10);
        dx1 = Mesh.mfdx1.Dadj ; dx2 = Mesh.mfdx2.Dadj; dx3 = Mesh.mfdx3.Dadj;

        By = [dx1(a*u(:,1))+dx2(a*u(:,2))+dx3(a*u(:,3))+dx1(b*u(:,10));...
            dx1(a*u(:,4))+dx2(a*u(:,5))+dx3(a*u(:,6))+dx2(b*u(:,10));...
            dx1(a*u(:,7))+dx2(a*u(:,8))+dx3(a*u(:,9))+dx3(b*u(:,10))...
            ];

    otherwise
        error('unknown flag: %s',flag);
end

% get diagonal of d2S (interesting in matrix free mode)
function D = getDiag(Mesh,mu,lambda,alpha)
dim = Mesh.dim;
vol = Mesh.vol;
if dim==2
    xn = Mesh.xn;
    x1 =  Mesh.mfPi(xn,1);
    x2 =  Mesh.mfPi(xn,2);
    x3 =  Mesh.mfPi(xn,3);
    e1 =  x1 - x3;
    e2 =  x2 - x3;

    % compute gradients of basis functions
    dphi1 =   [ e2(:,2) -e2(:,1)]./[2*vol,2*vol];
    dphi2 =   [-e1(:,2)  e1(:,1)]./[2*vol,2*vol];
    dphi3 = -dphi1 - dphi2;

    % get boundaries
    Dxi = @(i)   Mesh.mfPi(vol.*dphi1(:,i).^2,1) ...
        + Mesh.mfPi(vol.*dphi2(:,i).^2,2) ...
        + Mesh.mfPi(vol.*dphi3(:,i).^2,3); % diagonal of Dx1'*Dx1

    D1 = (2*mu+lambda)*Dxi(1)+            mu*Dxi(2);
    D2 = mu*Dxi(1)           + (2*mu+lambda)*Dxi(2);
    D = [D1;D2];
else
    xn = Mesh.xn;
    % compute edges
    v1   = Mesh.mfPi(xn,1);
    v2   = Mesh.mfPi(xn,2);
    v3   = Mesh.mfPi(xn,3);
    v4   = Mesh.mfPi(xn,4);
    e1   =  v1 - v4;
    e2   =  v2 - v4;
    e3   =  v3 - v4;
    % compute inverse transformation to reference element
    Ainv =  [
        e2(:,2).*e3(:,3)-e2(:,3).*e3(:,2), ...
        -(e1(:,2).*e3(:,3)-e1(:,3).*e3(:,2)), ...
        e1(:,2).*e2(:,3)-e1(:,3).*e2(:,2), ...
        -(e2(:,1).*e3(:,3)-e2(:,3).*e3(:,1)), ...
        e1(:,1).*e3(:,3)-e1(:,3).*e3(:,1), ...
        -(e1(:,1).*e2(:,3)-e1(:,3).*e2(:,1)), ...
        e2(:,1).*e3(:,2)-e2(:,2).*e3(:,1), ...
        -(e1(:,1).*e3(:,2)-e1(:,2).*e3(:,1)), ...
        e1(:,1).*e2(:,2)-e1(:,2).*e2(:,1), ...
        ];
    detA = e1(:,1).*Ainv(:,1) + e2(:,1).*Ainv(:,2)+ e3(:,1).*Ainv(:,3);

    % compute gradients of basis functions
    dphi1 =   Ainv(:,1:3)./repmat(detA,[1 3]);
    dphi2 =   Ainv(:,4:6)./repmat(detA,[1 3]);
    dphi3 =   Ainv(:,7:9)./repmat(detA,[1 3]);
    dphi4 =   1 - dphi1 - dphi2 - dphi3;

    Dxi = @(i) Mesh.mfPi(vol.*dphi1(:,i).^2,1) +  Mesh.mfPi(vol.*dphi2(:,i).^2,2)  ...
        + Mesh.mfPi(vol.*dphi3(:,i).^2,3) +  Mesh.mfPi(vol.*dphi4(:,i).^2,4);

    D1 = (2*mu+lambda)*Dxi(1)+            mu*Dxi(2) + mu *Dxi(3);
    D2 = (2*mu+lambda)*Dxi(2)+            mu*Dxi(1) + mu *Dxi(3);
    D3 = (2*mu+lambda)*Dxi(3)+            mu*Dxi(1) + mu *Dxi(2);
    D = [D1;D2;D3];
end
D = alpha*reshape(D,[],1);

function B = getElasticMatrixFEM(Mesh,mu,lambda)
vol = sqrt(Mesh.vol);
a   = sdiag(sqrt(mu)*vol);
b   = sdiag(sqrt(mu+lambda)*vol);
switch Mesh.dim
    case 2
        dx1 = Mesh.dx1; dx2 = Mesh.dx2;
        B = [a*dx1;a*dx2];
        B = [ blkdiag(B,B); b*dx1, b*dx2];
    case 3
        dx1 = Mesh.dx1; dx2 = Mesh.dx2; dx3 = Mesh.dx3;
        B = [a*dx1;a*dx2;a*dx3];
        B = [ blkdiag(B,B,B); b*dx1, b*dx2, b*dx3];
end

% shortcut for sparse diagonal matrix
function A = sdiag(v)
A = spdiags(v(:),0,numel(v),numel(v));

function runMinimalExample
% ====== 2D
omega = [0,10,0,8]; m = [17,16]; p = [5,6];
w = zeros([p,2]);  w(3,3,1) = 0.06; w(3,4,2) = -0.05;
xn = reshape(getNodalGrid(omega,m),[],2);
yn = splineTransformation2D(w(:),xn(:),'omega',omega,'m',m+1,'p',p,'Q',[]);

Mesh = TriMesh1(omega,m);

regOptn = {'alpha',1,'mu',1,'lambda',0};

fctn = @(y) feval(mfilename,y-xn(:),xn(:),Mesh,regOptn{:},'matrixFree',0);
[Sc, dS, d2S] = fctn(yn(:));
checkDerivative(fctn,yn(:));

[Scmf, dSmf, d2Smf] = feval(mfilename,yn(:)-xn(:),xn(:),Mesh,regOptn{:},'matrixFree',1);

err1 = abs(Sc-Scmf)/abs(Sc);
err2 = norm(dS-dSmf)/norm(dS);

y = randn(5*size(Mesh.tri,1),1); x = rand(numel(xn),1);
err3 = norm(d2S*x-d2Smf.d2S(x,Mesh))/norm(d2S*x);
err4 = abs(y'*d2Smf.By(x,Mesh)-x'*d2Smf.BTy(y,Mesh));
fprintf('%s, MFerror(S,dS,d2S): [%1.2e,%1.2e,%1.2e] OK? %d\n\t\tAdjointError: %1.2e, OK? %d\n',...
    mfilename,err1,err2,err3,max([err1;err2;err3])<1e-14,err4,err4<1e-13);
err5 = norm(full(diag(d2S))-d2Smf.diag(Mesh))/norm(full(diag(d2S)));
fprintf('\t\tDiagErrorMF: %1.2e OK? %d\n',err5,err5<1e-13);

% ====== 3D
omega = [0,10,0,8,0 3]; m = [3,5,6];
Mesh = TetraMesh1(omega,m);

xn = Mesh.xn;
yn = xn + randn(size(xn));

regOptn = {'alpha',1,'mu',1,'lambda',0};

fctn = @(y) feval(mfilename,y-xn(:),xn(:),Mesh,regOptn{:},'matrixFree',0);
[Sc, dS, d2S] = fctn(yn(:));
checkDerivative(fctn,yn(:));

[Scmf, dSmf, d2Smf] = feval(mfilename,yn(:)-xn(:),xn(:),Mesh,regOptn{:},'matrixFree',1);

err1 = abs(Sc-Scmf)/abs(Sc);
err2 = norm(dS-dSmf)/norm(dS);

y = randn(10*size(Mesh.tri,1),1); x = rand(numel(xn),1);
err3 = norm(d2S*x-d2Smf.d2S(x,Mesh))/norm(d2S*x);
err4 = abs(y'*d2Smf.By(x,Mesh)-x'*d2Smf.BTy(y,Mesh));
fprintf('%s, MFerror(S,dS,d2S): [%1.2e,%1.2e,%1.2e] OK? %d\n\t\tAdjointError: %1.2e, OK? %d\n',...
    mfilename,err1,err2,err3,max([err1;err2;err3])<1e-14,err4,err4<1e-13);
err5 = norm(full(diag(d2S))-d2Smf.diag(Mesh))/norm(full(diag(d2S)));
fprintf('\t\tDiagErrorMF: %1.2e OK? %d\n',err5,err5<1e-13);


text:
function [flow_out, time_shifted_points_out] = em1_flow_with_prior(...
    events, ...
    feature_pos, ...
    prev_shifted_points, ...
    prev_shifted_weights, ...
    flow_init, ...
    params, ...
    fig)
%EM1_FLOW_WITH_PRIOR Estimates optical flow for an event feature.
%
% EM1_FLOW_WITH_PRIOR estimates the optical flow of a set of events
% combined to form a 'feature'. This EM step uses the time shifted points
% from the previous iteration as a template, as in:
% Alex Zihao Zhu, Nikolay Atanasov and Kostas Daniilidis.
% "Event-based Visual Inertial Odometry",
% IEEE International Conference on Computer Vision and Pattern Recognition (CVPR), 2017.
%
% Syntax:  [flow_out, time_shifted_points_out] = EM1_FLOW_WITH_PRIOR(...
%            events, ...
%            feature_pos, ...
%            prev_shifted_points, ...
%            prev_shifted_weights, ...
%            flow_init, ...
%            params, ...
%            fig)
%
% Inputs:
%    events               - 4xN, each column is (x,y,t,p).
%    feature_pos          - 2x1, pixel position of the feature.
%    prev_shifted_points  - 2xM, time shifted points from the previous
%                           iteration, used as a template.
%    prev_shifted_weights - 1xM, weights for each shifted point.
%    flow_init            - 2x1, initialization for the flow.
%    params               - parameters, defined in get_params().
%    fig                  - figure handle for plotting.
%
% Outputs:
%    flow_out                - 2x1, estimated flow.
%    time_shifted_points_out - 2xN, input events in the feature window
%                              shifted by the flow [x; y] + dt * flow.
%
% See also GET_PARAMS
%
% Author: Alex Zihao Zhu, University of Pennsylvania
% Email: alexzhu(at)seas.upenn.edu
% Copyright 2018 University of Pennsylvania
% Alex Zihao Zhu, Nikolay Atanasov, Kostas Daniilidis
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS, CONTRIBUTORS, AND THE
% TRUSTEES OF THE UNIVERSITY OF PENNSYLVANIA "AS IS" AND ANY EXPRESS OR
% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT OWNER, CONTRIBUTORS OR THE TRUSTEES OF
% THE UNIVERSITY OF PENNSYLVANIA BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
% NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
% THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%% Initialization
flow = flow_init;
prev_flow = flow;
% Estimated flow
flow_out = [nan; nan];
% Estimated canonical events
time_shifted_points_out = [];
delta_flows = zeros(params.em1_params.max_iters, 1);
num_iter = 0;
scatter_plot_handle = [];
event_window = [];

target_time = events(3, 1);

centered_events = events;
centered_events(1:2, :) = bsxfun(@minus, events(1:2, :), feature_pos);

kdtree = KDTreeSearcher(...
    prev_shifted_points' / (sqrt(2) * params.em1_params.sigma), ...
    'Distance', 'euclidean');

if params.debug
    set(0, 'CurrentFigure', fig)
    subplot(2,1,2)
    if ~isempty(prev_shifted_points)
        scatter(prev_shifted_points(1, :), prev_shifted_points(2, :), prev_shifted_weights * 10, 'b.')
    end
end

%% Main EM Loop
while true
    if num_iter > params.em1_params.max_iters
        return
    end

    time_shifted_points = centered_events(1:2,:) + ...
        bsxfun(@times, flow,(target_time-centered_events(3,:)));

    if isempty(event_window)
        event_window = ...
            time_shifted_points(1, :) >= -params.window_size/2 &...
            time_shifted_points(2, :) >= -params.window_size/2 & ...
            time_shifted_points(1, :) <= params.window_size/2 & ...
            time_shifted_points(2, :) <= params.window_size/2;

        n_in_window = sum(event_window);
        if n_in_window < params.min_events_for_em
            return
        end

        centered_events = centered_events(:, event_window);
        time_shifted_points = time_shifted_points(:, event_window);

        weights = zeros(size(prev_shifted_points, 2), size(time_shifted_points, 2));
    end

    [neighbors_cell, distances_cell] = rangesearch(...
        kdtree, ...
        time_shifted_points' / (sqrt(2) * params.em1_params.sigma), ...
        params.em1_params.max_distance);

    distancesstacked = cell2mat(cellfun(...
        @transpose, distances_cell, 'UniformOutput', false))';

    if isempty(distancesstacked)
        return
    end

    % NOTE: 'length' is not the same as @length
    num_neighbors = cellfun('length', neighbors_cell);
    prev_correspondences = cell2mat(cellfun(...
        @transpose, neighbors_cell, 'UniformOutput', false))';
    curr_correspondences = repelem(1:size(time_shifted_points, 2), num_neighbors);

    % It's cheaper to multiply to 0 to reset the weights matrix than to
    % reinitialize it using zeros.
    weightsstacked = exp(-distancesstacked);
    weights = weights * 0;

    valid_inds = sub2indc(...
        curr_correspondences, ...
        prev_correspondences, ...
        size(weights));

    weights(valid_inds) = weightsstacked;
    weights = bsxfun(@times, prev_shifted_weights, weights);
    weights_sum = sum(weights, 1);
    valid_weights = weights_sum > 0;
    weights = bsxfun(@rdivide, weights, weights_sum + 1e-10);

    % Simultaneously minimize over flow and translation
    weighted_prior_points = prev_shifted_points * weights;
    weighted_prior_points = weighted_prior_points(:, valid_weights);

    valid_centered_events = centered_events(:, valid_weights);
    dx = weighted_prior_points(1:2, :) - valid_centered_events(1:2, :);
    dt = target_time - valid_centered_events(3, :);

    flow = (dx*dt') / (dt*dt');

    %% Calculate change in flow, plot debug information.
    if (norm(flow - prev_flow) < params.em1_params.min_err)
        break;
    end

    delta_flows(num_iter+1) = norm(flow - prev_flow);
    prev_flow = flow;

    if params.debug
        set(0, 'CurrentFigure', fig)
        subplot(2,1,1)
        plot(delta_flows(delta_flows > 0),'b')
        title('EM1 with prior change in flow (convergence criterion)')
        xlim([0 params.em1_params.max_iters])

        subplot(2,1,2)
        if (~isempty(scatter_plot_handle))
            delete(scatter_plot_handle)
        end

        if ~isempty(prev_shifted_points)
            hold on
        end

        scatter_plot_handle = scatter(time_shifted_points(1, :), time_shifted_points(2, :),'r.');
        hold off
        axis equal
        axis([-params.window_size/2-5 params.window_size/2+5 -params.window_size/2-5 params.window_size/2+5])
        axis ij

        title('EM1 with prior time shifted events')
        pause(0.01)
    end

    num_iter = num_iter + 1;
end

if params.debug
    pause(0.5)
end

flow_out = flow;

dt = events(3, end) - events(3, 1);
centered_events = bsxfun(@minus, events(1:2, :), feature_pos + flow * dt);
time_shifted_points = centered_events(1:2, :) + ...
    bsxfun(@times, flow, (events(3, end) - events(3,:)));

% Make the window a little bigger.
window_size = round(params.window_size * 1.5);

event_window = time_shifted_points(1, :) >= -window_size/2 & ...
    time_shifted_points(2, :) >= -window_size/2 & ...
    time_shifted_points(1, :) <= window_size/2 & ...
    time_shifted_points(2, :) <= window_size/2;

time_shifted_points_out = time_shifted_points(:, event_window);
end

author: daniilidis-group


text:
% Chapter 11 - Hamiltonian Systems, Lyapunov Functions, and Stability.
% Programs_11a - Surface and Contour Plots. Finding Vdot.
% Copyright Birkhauser 2013. Stephen Lynch.

% The double-well potential (Fig. 11.5(b)).
% Symbolic Math toolbox required.
ezsurfc('-x^2/2+x^4/4+y^2/2',[-1.5,1.5]);

% Calculating Vdot (Exercise 6).
clear
syms x y V
V=(1-4*x^2-y^2)^2;
xdot=-y*(1+x)/2+x*(1-4*x^2-y^2);
ydot=2*x*(1+x)+y*(1-4*x^2-y^2);
Vdot=factor(diff(V,x)*xdot+diff(V,y)*ydot)

% Phase portraits may be plotted using the M-files in Chapter 8.
% End of Programs_11a.

author: Sable

repo: mcbench-benchmarks


text:
function [ h, hmin, hmax ] = p01_stepsize ( option )

%*****************************************************************************80
%
%% P01_STEPSIZE returns step sizes for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =    0.30000;
  hmin = 0.03125;
  hmax = 4.00000;

  return
end


text:
% std_readersp() - load ERSP measures for data channels or  for all
%                  components of a specified cluster. This function is
%                  also being used to read ITC and ERPimage data.
% Usage:
%   >> [STUDY, erspdata, times, freqs, erspbase] = ...
%                   std_readersp(STUDY, ALLEEG, varargin);
% Inputs:
%       STUDY - studyset structure containing some or all files in ALLEEG
%      ALLEEG - vector of loaded EEG datasets
%
% Optional inputs:
%  'design'    - [integer] read files from a specific STUDY design. Default
%                is empty (use current design in STUDY.currentdesign).
%  'channels'  - [cell] list of channels to import {default: none}
%  'clusters'  - [integer] list of clusters to import {[]|default: all but
%                the parent cluster (1) and any 'NotClust' clusters}
%  'singletrials' - ['on'|'off'] load single trials spectral data (if
%                available). Default is 'off'.
%  'forceread' - ['on'|'off'] Force rereading data from disk.
%                Default is 'off'.
%  'subject'   - [string] select a specific subject {default:all}
%  'component' - [integer] select a specific component in a cluster
%                 {default:all}
%  'datatype'  - {'ersp'|'itc'|'erpim'} This function is used to read all
%                2-D STUDY matrices stored on disk (not only ERSP). It may
%                read ERSP ('ersp' option), ITC ('itc' option) or ERPimage
%                data ('erpim' option).
%
% ERSP specific options:
%  'timerange' - [min max] time range {default: whole measure range}
%  'freqrange' - [min max] frequency range {default: whole measure range}
%  'subbaseline' - ['on'|'off'] subtract the ERSP baseline for paired
%                conditions. The conditions for which baseline is removed
%                are indicated on the command line. See help
%                std_studydesign for more information about paired and
%                unpaired variables.
%
% ERPimage specific option:
%   This function is used to read all 2-D STUDY matrices stored on disk
%   (this includes ERPimages). It therefore takes as input specific
%   ERPimage options. Note that the 'singletrials' optional input is
%   irrelevant for ERPimages (which are always stored as single trials).
%   'concatenate' - ['on'|'off'] read concatenated ERPimage data ('on') or
%                   stacked ERPimage data. See help std_erpimage for more
%                   information.
%   'timerange'   - [min max] time range {default: whole measure range}
%   'trialrange'  - [min max] read only a specific range of the ERPimage
%                   output trials {default: whole measure range}
% Output:
%  STUDY    - updated studyset structure
%  erspdata - [cell array] ERSP data (the cell array size is
%             condition x groups). This may also be ITC data or ERPimage
%             data (see above).
%  times    - [float array] array of time points
%  freqs    - [float array] array of frequencies. For ERPimage this
%             contains trial indices.
%  erspbase - [cell array] baseline values
%  events   - [cell array] events (ERPimage only).
%
% Author: Arnaud Delorme, CERCO, 2006-

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% 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 [STUDY, erspdata, alltimes, allfreqs, erspbase, events, unitPower] = std_readersp(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_readersp;
    return;
end
events = {};
unitPower = 'dB';
if ~isstruct(ALLEEG) % old calling format
    dataset = ALLEEG;
    EEG = STUDY(dataset);
    comp = varargin{1};
    if length(varargin) > 1, timelim = varargin{2}; else timelim = []; end;
    if length(varargin) > 2, freqlim = varargin{3}; else freqlim = []; end;
    filebase = fullfile(EEG.filepath, EEG.filename(1:end-4));
    if comp < 0
        error('Old format function call, channel reading not supported');
        [ersp params alltimes allfreqs] = std_readfile(filebase, 'channels', -comp, 'timelimits', timelim, 'freqlimits', freqlim);
    else
        [ersp     params alltimes allfreqs] = std_readfile(filebase, 'components', comp, 'timelimits', timelim, 'freqlimits', freqlim, 'measure', 'ersp');
        [erspbase params alltimes allfreqs] = std_readfile(filebase, 'components', comp, 'timelimits', timelim, 'freqlimits', freqlim, 'measure', 'erspbase');
    end;
    STUDY = ersp;
    erspdata = allfreqs;
    allfreqs = params;
    return;
end;

STUDY = pop_erspparams( STUDY, 'default');
STUDY = pop_erpimparams(STUDY, 'default');
[opt moreopts] = finputcheck( varargin, { ...
    'design'        'integer' []             STUDY.currentdesign;
    'channels'      'cell'    []             {};
    'clusters'      'integer' []             [];
    'trialrange'    'real'    []             STUDY.etc.erpimparams.trialrange;
    'freqrange'     'real'    []             STUDY.etc.erspparams.freqrange;
    'timerange'     'real'    []             NaN;
    'singletrials'  'string'  { 'on','off' } 'off';
    'forceread'     'string'  { 'on','off' } 'off';
    'concatenate'   'string'  { 'on','off' } STUDY.etc.erpimparams.concatenate;
    'subbaseline'   'string'  { 'on','off' } STUDY.etc.erspparams.subbaseline;
    'component'     'integer' []             [];
    'infotype'      'string'  { 'ersp','itc','erpim' } 'ersp'; ... % deprecated
    'datatype'      'string'  { 'ersp','itc','erpim' } 'ersp'; ...
    'subject'       'string'  []             '' }, ...
    'std_readersp', 'ignore');
if isstr(opt), error(opt); end;
if ~strcmpi(opt.infotype, 'ersp'), opt.datatype = opt.infotype; end;
if strcmpi(opt.datatype, 'erpim'), if isnan(opt.timerange), opt.timerange = STUDY.etc.erpimparams.timerange; end;
                                   opt.freqrange = opt.trialrange;
                                   ordinate      = 'trials';
else                               if isnan(opt.timerange) opt.timerange = STUDY.etc.erspparams.timerange; end;
                                   ordinate      = 'freqs';
end;
nc = max(length(STUDY.design(opt.design).variable(1).value),1);
ng = max(length(STUDY.design(opt.design).variable(2).value),1);
paired1 = STUDY.design(opt.design).variable(1).pairing;
paired2 = STUDY.design(opt.design).variable(2).pairing;

dtype = opt.datatype;

% find channel indices
% --------------------
if ~isempty(opt.channels)
     allChangrp = lower({ STUDY.changrp.name });
     finalinds = std_chaninds(STUDY, opt.channels);
else finalinds = opt.clusters;
end;

for ind = 1:length(finalinds)
    erspbase = cell( nc, ng );

    % find indices
    % ------------
    if ~isempty(opt.channels)
        tmpstruct = STUDY.changrp(finalinds(ind));
        allinds       = tmpstruct.allinds;
        setinds       = tmpstruct.setinds;
    else
        tmpstruct = STUDY.cluster(finalinds(ind));
        allinds       = tmpstruct.allinds;
        setinds       = tmpstruct.setinds;
    end;

    % check if data is already here
    % -----------------------------
    if strcmpi(dtype, 'erpim') % ERP images
        dataread = 0;
        eqtf  = isequal( STUDY.etc.erpimparams.timerange , opt.timerange) && ...
                isequal( STUDY.etc.erpimparams.trialrange, opt.trialrange);
        if isfield(tmpstruct, 'erpimdata') && eqtf && ~isempty(tmpstruct.erpimdata)
            dataread = 1;
        end;
    else
        dataread = 0;
        eqtf  = isequal( STUDY.etc.erspparams.timerange, opt.timerange) && ...
                isequal( STUDY.etc.erspparams.freqrange, opt.freqrange);
        eqtfb = isequal( STUDY.etc.erspparams.subbaseline, opt.subbaseline) && eqtf;
        if strcmpi(opt.singletrials,'off')
            if isfield(tmpstruct, [ dtype 'data']) && eqtfb && ~isempty(getfield(tmpstruct, [ dtype 'data']))
                dataread = 1;
            end;
        else
            if isfield(tmpstruct, [ dtype 'datatrials']) && eqtf && ~isempty(getfield(tmpstruct, [ dtype 'datatrials']))
                if ~isempty(opt.channels) && strcmpi(getfield(tmpstruct, [ dtype 'trialinfo' ]), opt.subject)
                    dataread = 1;
                elseif isempty(opt.channels) && isequal(getfield(tmpstruct, [ dtype 'trialinfo' ]), opt.component)
                    dataread = 1;
                end;
            end;
        end;
    end;

    if dataread && strcmpi(opt.forceread, 'off')
        disp('Using pre-loaded data. To force rereading data from disk use the ''forceread'' flag');
    else

        if strcmpi(dtype, 'erpim') % ERP images

            if strcmpi(opt.singletrials, 'on')
                error( [ 'Single trial loading option not supported with STUDY ERP-image' 10 '(there is no such thing as a single-trial ERPimage)' ]);
            end;
            % read the data and select channels
            % ---------------------------------
            setinfo = STUDY.design(opt.design).cell;
            erpim    = cell( nc, ng );
            events   = cell( nc, ng );
            for c = 1:nc
                for g = 1:ng
                    if ~isempty(setinds{c,g})
                        if ~isempty(opt.channels), opts = { 'channels',    allChangrp(allinds{c,g}(:)), 'timelimits', opt.timerange, 'triallimits', opt.trialrange, 'concatenate', opt.concatenate };
                        else                       opts = { 'components',  allinds{c,g}(:), 'timelimits', opt.timerange, 'triallimits', opt.trialrange, 'concatenate', opt.concatenate };
                        end;
                        [erpim{c, g} tmpparams alltimes alltrials events{c,g}] = std_readfile( setinfo(setinds{c,g}(:)), 'measure', 'erpim', opts{:});
                        fprintf('.');
                    end;
                end;
            end;
            if strcmpi(opt.concatenate, 'on'), alltrials = []; end;
            fprintf('\n');
            ersp     = erpim;
            allfreqs = alltrials;

        else % ERSP/ITC

            % reserve arrays
            % --------------
            events   = {};
            ersp     = cell( nc, ng );
            erspinds = cell( nc, ng );

            % find total nb of trials
            % THIS CODE IS NOT NECESSARY ANY MORE (SEE BUG 1170)
            % -----------------------
            %setinfo = STUDY.design(opt.design).cell;
            %tottrials = cell( nc, ng );
            %if strcmpi(opt.singletrials, 'on')
            %    for indSet = 1:length(setinfo)
            %        condind = std_indvarmatch( setinfo(indSet).value{1}, STUDY.design(opt.design).variable(1).value );
            %        grpind  = std_indvarmatch( setinfo(indSet).value{2}, STUDY.design(opt.design).variable(2).value );
            %        if isempty(tottrials{condind, grpind}), tottrials{condind, grpind} = sum(cellfun(@length, setinfo(indSet).trials));
            %        else       tottrials{condind, grpind} = tottrials{condind, grpind} + sum(cellfun(@length, setinfo(indSet).trials));
            %        end;
            %    end;
            %end;

            % read the data and select channels
            % ---------------------------------
            fprintf('Reading all %s data:', upper(dtype));
            setinfo = STUDY.design(opt.design).cell;
            if strcmpi(opt.singletrials, 'on')
                for c = 1:nc
                    for g = 1:ng
                        if ~isempty(opt.channels)
                            allsubjects = { STUDY.design(opt.design).cell.case };
                            if ~isempty(opt.subject), inds = strmatch( opt.subject, allsubjects(setinds{c,g}));
                            else inds = 1:length(allinds{c,g}); end;
                        else
                            if ~isempty(opt.component) inds = find( allinds{c,g} == STUDY.cluster(finalinds(ind)).comps(opt.component));
                            else inds = 1:length(allinds{c,g}); end;
                        end;
                        if ~isempty(inds)
                            count{c, g} = 1;
                            for indtmp = 1:length(inds)
                                setindtmp = STUDY.design(opt.design).cell(setinds{c,g}(inds(indtmp)));
                                tmpopts = { 'measure', 'timef' 'timelimits', opt.timerange, 'freqlimits', opt.freqrange };
                                if ~isempty(opt.channels), [ tmpersp tmpparams alltimes allfreqs] = std_readfile(setindtmp, 'channels', allChangrp(allinds{c,g}(inds(indtmp))), tmpopts{:});
                                else                       [ tmpersp tmpparams alltimes allfreqs] = std_readfile(setindtmp, 'components',          allinds{c,g}(inds(indtmp)),  tmpopts{:});
                                end;
                                indices = [count{c, g}:count{c, g}+size(tmpersp,3)-1];
                                if indtmp == 1
                                    ersp{c, g} = permute(tmpersp, [2 1 3]);
                                else ersp{c, g}(:,:,indices) = permute(tmpersp, [2 1 3]);
                                end;
                                erspinds{c, g}(1:2,indtmp) = [ count{c, g} count{c, g}+size(tmpersp,3)-1 ];
                                count{c, g} = count{c, g}+size(tmpersp,3);
                                if size(tmpersp,3) ~= sum(cellfun(@length, setindtmp.trials))
                                    error( sprintf('Wrong number of trials in datafile for design index %d\n', setinds{c,g}(inds(indtmp))));
                                end;
                            end;
                        end;
                    end;
                end;
            else
                for c = 1:nc
                    for g = 1:ng
                        if ~isempty(setinds{c,g})
                            if ~isempty(opt.channels), opts = { 'channels',  allChangrp(allinds{c,g}(:)), 'timelimits', opt.timerange, 'freqlimits', opt.freqrange };
                            else                       opts = { 'components',           allinds{c,g}(:) , 'timelimits', opt.timerange, 'freqlimits', opt.freqrange };
                            end;
                            if strcmpi(dtype, 'ersp')
                                 erspbase{c, g}                             = std_readfile( setinfo(setinds{c,g}(:)), 'measure', 'erspbase', opts{:});
                                 [ ersp{c, g} tmpparams alltimes allfreqs ] = std_readfile( setinfo(setinds{c,g}(:)), 'measure', 'ersp'    , opts{:});
                            else [ ersp{c, g} tmpparams alltimes allfreqs ] = std_readfile( setinfo(setinds{c,g}(:)), 'measure', 'itc'     , opts{:});
                                 ersp{c, g} = ersp{c, g};
                            end;
                            fprintf('.');
                            %ersp{c, g}      = permute(ersp{c, g}           , [3 2 1]);
                            %erspbase{c, g}  = 10*log(permute(erspbase{c, g}, [3 2 1]));
                        end;
                    end;
                end;
            end;
            fprintf('\n');

            % compute ERSP or ITC if trial mode
            % (since only the timef have been loaded)
            % ---------------------------------------
            if strcmpi(opt.singletrials, 'on')
                tmpparams2      = fieldnames(tmpparams); tmpparams2 = tmpparams2';
                tmpparams2(2,:) = struct2cell(tmpparams);
                precomp.times = alltimes;
                precomp.freqs = allfreqs;
                precomp.recompute = dtype;
                for c = 1:nc
                    for g = 1:ng
                        if ~isempty(ersp{c,g})
                            precomp.tfdata = permute(ersp{c,g}, [2 1 3]);
                            if strcmpi(dtype, 'itc')
                                [tmp ersp{c,g}] = newtimef(zeros(ALLEEG(1).pnts,2), ALLEEG(1).pnts, [ALLEEG(1).xmin ALLEEG(1).xmax]*1000, ...
                                    ALLEEG(1).srate, [], tmpparams2{:}, 'precomputed', precomp, 'verbose', 'off');
                            elseif strcmpi(dtype, 'ersp')
                                [ersp{c,g} tmp] = newtimef(zeros(ALLEEG(1).pnts,2), ALLEEG(1).pnts, [ALLEEG(1).xmin ALLEEG(1).xmax]*1000, ...
                                    ALLEEG(1).srate, [], tmpparams2{:}, 'precomputed', precomp, 'verbose', 'off');
                            end;
                            ersp{c,g} = permute(ersp{c,g}, [2 1 3]);
                        end;
                    end;
                end;
            end;

            % compute average baseline across groups and conditions
            % -----------------------------------------------------
            if strcmpi(opt.subbaseline, 'on') && strcmpi(dtype, 'ersp')
                if strcmpi(opt.singletrials, 'on')
                    disp('WARNING: no ERSP baseline may not be subtracted when using single trials');
                else
                    disp('Recomputing baseline...');
                    if strcmpi(paired1, 'on') && strcmpi(paired2, 'on')
                        disp('Removing ERSP baseline for both indep. variables');
                        meanpowbase = computeerspbaseline(erspbase(:), opt.singletrials);
                        ersp        = removeerspbaseline(ersp, erspbase, meanpowbase);
                    elseif strcmpi(paired1, 'on')
                        disp('Removing ERSP baseline for first indep. variables (second indep. var. is unpaired)');
                        for g = 1:ng        % ng = number of groups
                            meanpowbase = computeerspbaseline(erspbase(:,g), opt.singletrials);
                            ersp(:,g)   = removeerspbaseline(ersp(:,g), erspbase(:,g), meanpowbase);
                        end;
                    elseif strcmpi(paired2, 'on')
                        disp('Removing ERSP baseline for second indep. variables (first indep. var. is unpaired)');
                        for c = 1:nc        % ng = number of groups
                            meanpowbase = computeerspbaseline(erspbase(c,:), opt.singletrials);
                            ersp(c,:)   = removeerspbaseline(ersp(c,:), erspbase(c,:), meanpowbase);
                        end;
                    else
                        disp('Not removing ERSP baseline (both indep. variables are unpaired');
                    end;
                end;
            end;
        end;

%         if strcmpi(opt.statmode, 'common')
%             % collapse the two last dimensions before computing significance
%             % i.e. 18 subject with 4 channels -> same as 4*18 subjects
%             % --------------------------------------------------------------
%             disp('Using all channels for statistics...');
%             for c = 1:nc
%                 for g = 1:ng
%                     ersp{c,g} = reshape( ersp{c,g}, size(ersp{c,g},1), size(ersp{c,g},2), size(ersp{c,g},3)*size(ersp{c,g},4));
%                 end;
%             end;
%         end;

        % copy data to structure
        % ----------------------
        if ~isempty(events)
            tmpstruct = setfield(tmpstruct, [ dtype 'events' ], events);
        end;
        tmpstruct = setfield(tmpstruct, [ dtype ordinate ], allfreqs);
        tmpstruct = setfield(tmpstruct, [ dtype 'times'  ], alltimes);
        if strcmpi(opt.singletrials, 'on')
            tmpstruct = setfield(tmpstruct, [ dtype 'datatrials' ], ersp);
            tmpstruct = setfield(tmpstruct, [ dtype 'subjinds' ], erspinds);
            tmpstruct = setfield(tmpstruct, [ dtype 'times' ], alltimes);
            if ~isempty(opt.channels)
                 tmpstruct = setfield(tmpstruct, [ dtype 'trialinfo' ], opt.subject);
            else tmpstruct = setfield(tmpstruct, [ dtype 'trialinfo' ], opt.component);
            end;
        else
            tmpstruct = setfield(tmpstruct, [ dtype 'data' ], ersp);
            if strcmpi(dtype, 'ersp')
                tmpstruct = setfield(tmpstruct, [ dtype 'base' ], erspbase);
            end;
        end;

        % copy results to structure
        % -------------------------
        fields = { [ dtype 'data'     ] [ dtype 'events' ] [ dtype ordinate ] [ dtype 'datatrials' ] ...
                   [ dtype 'subjinds' ] [ dtype 'base'   ] [ dtype 'times'  ] [ dtype 'trialinfo'  ]  'allinds' 'setinds' };
        for f = 1:length(fields)
            if isfield(tmpstruct, fields{f}),
                tmpdata = getfield(tmpstruct, fields{f});
                if ~isempty(opt.channels)
                     STUDY.changrp = setfield(STUDY.changrp, { finalinds(ind) }, fields{f}, tmpdata);
                else STUDY.cluster = setfield(STUDY.cluster, { finalinds(ind) }, fields{f}, tmpdata);
                end;
            end;
        end;
    end;
end;

% output units
% -----------
if exist('tmpparams') ~= 1
    tmpparams = []; ctmp = 1;
    while isempty(tmpparams) && ctmp <= length(STUDY.design(opt.design).cell)
        try
            if ~isempty(opt.channels), [tmpersp tmpparams] = std_readfile(STUDY.design(opt.design).cell(ctmp), 'channels', {STUDY.changrp(1).name}, 'measure', opt.datatype);
            else                       [tmpersp tmpparams] = std_readfile(STUDY.design(opt.design).cell(ctmp), 'components', STUDY.cluster(finalinds(end)).allinds{1,1}(1), 'measure', opt.datatype);
            end;
        catch
        end
        ctmp = ctmp + 1;
    end
end;
if ~isfield(tmpparams, 'baseline'), tmpparams.baseline = 0;     end;
if ~isfield(tmpparams, 'scale'   ), tmpparams.scale    = 'log'; end;
if ~isfield(tmpparams, 'basenorm'), tmpparams.basenorm = 'off'; end;
if strcmpi(tmpparams.scale, 'log')
    if strcmpi(tmpparams.basenorm, 'on')
        unitPower = '10*log(std.)'; % impossible
    elseif isnan(tmpparams.baseline)
        unitPower = '10*log10(\muV^{2}/Hz)';
    else
        unitPower = 'dB';
    end;
else
    if strcmpi(tmpparams.basenorm, 'on')
        unitPower = 'std.';
    elseif isnan(tmpparams.baseline)
        unitPower = '\muV^{2}/Hz';
    else
        unitPower = '% of baseline';
    end;
end;

% return structure
% ----------------
if nargout <2
    return
end
allinds   = finalinds;
if ~isempty(opt.channels)
    structdat = STUDY.changrp;
    erspdata  = cell(nc, ng);
    events    = {};
    for ind =  1:length(erspdata(:))
        if strcmpi(opt.singletrials, 'on')
             tmpdat = getfield(structdat(allinds(1)), [ dtype 'datatrials' ]);
        else tmpdat = getfield(structdat(allinds(1)), [ dtype 'data' ]);
        end;
        if ndims(tmpdat{ind}) == 2, erspdata{ind} = zeros([ size(tmpdat{ind}) 1 length(allinds)]);
        else                        erspdata{ind} = zeros([ size(tmpdat{ind}) length(allinds)]);
        end;
        for index = 1:length(allinds)
            if strcmpi(opt.singletrials, 'on')
                 tmpdat = getfield(structdat(allinds(index)), [ dtype 'datatrials' ]);
            else tmpdat = getfield(structdat(allinds(index)), [ dtype 'data' ]);
            end;
            erspdata{ind}(:,:,:,index) = tmpdat{ind};
            allfreqs = getfield(structdat(allinds(index)), [ dtype ordinate ]);
            alltimes = getfield(structdat(allinds(index)), [ dtype 'times'  ]);
            if isfield(structdat, [ dtype 'events' ])
                events   = getfield(structdat(allinds(index)), [ dtype 'events' ]);
            end;
            compinds = structdat(allinds(index)).allinds;
            setinds  = structdat(allinds(index)).setinds;
        end;
        erspdata{ind} = permute(erspdata{ind}, [1 2 4 3]); % time freqs elec subjects
    end;
    if ~isempty(opt.subject) && strcmpi(opt.singletrials,'off')
        erspdata = std_selsubject(erspdata, opt.subject, setinds, { STUDY.design(opt.design).cell.case }, 2);
    end;
else
    if strcmpi(opt.singletrials, 'on')
         erspdata = getfield(STUDY.cluster(allinds(1)), [ dtype 'datatrials' ]);
    else erspdata = getfield(STUDY.cluster(allinds(1)), [ dtype 'data' ]);
    end;
    allfreqs = getfield(STUDY.cluster(allinds(1)), [ dtype ordinate ]);
    alltimes = getfield(STUDY.cluster(allinds(1)), [ dtype 'times'  ]);
    if isfield(STUDY.cluster, [ dtype 'events' ])
        events   = getfield(STUDY.cluster(allinds(1)), [ dtype 'events' ]);
    end;
    compinds = STUDY.cluster(allinds(1)).allinds;
    setinds  = STUDY.cluster(allinds(1)).setinds;
    if ~isempty(opt.component) && length(allinds) == 1 && strcmpi(opt.singletrials,'off')
        erspdata = std_selcomp(STUDY, erspdata, allinds, setinds, compinds, opt.component);
    end;
end;

% compute ERSP baseline
% ---------------------
function meanpowbase = computeerspbaseline(erspbase, singletrials)

    try
        len = length(erspbase(:));
        count = 0;
        for index = 1:len
            if ~isempty(erspbase{index})
                if strcmpi(singletrials, 'on')
                    if count == 0, meanpowbase = mean(erspbase{index},3);
                    else           meanpowbase = meanpowbase + mean(erspbase{index},3);
                    end;
                else
                    if count == 0, meanpowbase = erspbase{index};
                    else           meanpowbase = meanpowbase + erspbase{index};
                    end;
                end;
                count = count+1;
            end;
        end;
        meanpowbase = reshape(meanpowbase  , [size(meanpowbase,1) 1 size(meanpowbase,2)])/count;
    catch,
        error([ 'Problem while subtracting common ERSP baseline.' 10 ...
                'Common baseline subtraction is performed based on' 10 ...
                'pairing settings in your design. Most likelly, one' 10 ...
                'independent variable should not have its data paired.' ]);
    end;

% remove ERSP baseline
% ---------------------
function ersp = removeerspbaseline(ersp, erspbase, meanpowbase)
    convert2log = 0;
    for g = 1:size(ersp,2)        % ng = number of groups
        for c = 1:size(ersp,1)
            if ~isempty(erspbase{c,g}) && ~isempty(ersp{c,g})
                erspbasetmp = reshape(erspbase{c,g}, [size(meanpowbase,1) 1 size(meanpowbase,3)]);
                if any(erspbasetmp(:) > 1000)
                    convert2log = 1;
                end;
                tmpmeanpowbase = repmat(meanpowbase, [1 size(ersp{c,g},2) 1]);
                if convert2log
                     ersp{c,g} = ersp{c,g} - repmat(10*log10(erspbasetmp), [1 size(ersp{c,g},2) 1 1]) + 10*log10(tmpmeanpowbase);
                else ersp{c,g} = ersp{c,g} - repmat(erspbasetmp, [1 size(ersp{c,g},2) 1 1]) + tmpmeanpowbase;
                end;
            end;
        end;
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [Rob,Sen,Lmk,Obs] = correctKnownLmks(Rob, Sen, Raw, Lmk, Obs, Opt)

% CORRECTKNOWNLMKS  Correct known landmarks.
%   [Rob,Sen,Lmk,Obs] = correctKnownLmks(Rob, Sen, Raw, Lmk, Obs, Opt)
%   performs EKF corrections to a selection of observations, updating Map,
%   Rob, Sen, Lmk and Obs structures.
%
%   Input/output structures:
%       Rob:  the robot
%       Sen:  the sensor
%       Raw:  the raw data issued from Sen
%       Lmk:  the set of landmarks
%       Obs:  the set of observations for the sensor Sen
%       Opt:  the algorithm options
%
%   The selection of the landmarks to observe is based on active-search
%   procedures and individual compatibility tests to reject outliers.
%
%   This function takes care of:
%       1. Selection of observations to correct
%       2. Feature matching
%       3. EKF correction
%       4. Correction of landmark parameters out of the EKF
%       5. Landmark re-parametrization
%       6. Landmark deletion in case of corruption
%
%   The algorithm is configurable through a number of options in structure
%   Opt.correct. Edit USERDATA to access and modify these options.
%
%   See also PROJECTLMK, PROJEUCPNTINTOPINHOLEONROB, SMARTDELETELMK.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% steps in this function
% 0- update Rob and Sen info from Map
% 1- project all landmarks
% 2- select landmarks to correct. For each one:
% 3- do feature matching. If feature found:
% 4- compute innovation.
% 5- perform consistency test. If it is OK:
% 6- do correction
% 7- eventually reparametrize landmark
% 8- eventually delete corrupted landmarks
% 9- force covariance symmetry

global Map

% 0. UPDATE ROB AND SEN INFO FROM MAP
Rob = map2rob(Rob);
Sen = map2sen(Sen);

% 1. PROJECT ALL LMKS - get all expectations
for lmk = find([Lmk.used])

    Obs(lmk) = projectLmk(Rob,Sen,Lmk(lmk),Obs(lmk),Opt);

end  % --- all landmarks are now projected.

vis = [Obs.vis];

if any(vis) % Consider only visible observations

    % 2. SELECT LMKS TO CORRECT
    [lmksToObs,lmksToSkip] = selectLmksToObserve(Obs(vis),Opt.correct.nUpdates);

    % lmks to skip, update Obs info
    [Obs(lmksToSkip).measured] = deal(false);
    [Obs(lmksToSkip).matched]  = deal(false);
    [Obs(lmksToSkip).updated]  = deal(false);

    for lmk = lmksToObs % for each landmark to correct

        % Update Lmk search counter
        Lmk(lmk).nSearch = Lmk(lmk).nSearch + 1;

        % 3. TRY TO MATCH FEATURE
        Obs(lmk) = matchFeature(Sen,Raw,Obs(lmk));

        if Obs(lmk).matched

            % Update Lmk matched counter
            Lmk(lmk).nMatch = Lmk(lmk).nMatch + 1;

            % 4. COMPUTE INNOVATIONS
            Rob = map2rob(Rob);
            Sen = map2sen(Sen);

            if Opt.correct.reprojectLmks
                % re-project landmark for improved Jacobians
                Obs(lmk) = projectLmk(Rob,Sen,Lmk(lmk),Obs(lmk),Opt);
            end

            Obs(lmk) = observationInnovation(Obs(lmk),Opt.correct.lines.innType);

            % 5. CHECK CONSISTENCE
            if Obs(lmk).inn.MD2 < Opt.correct.MD2th

                % Update Lmk inlier counter
                Lmk(lmk).nInlier = Lmk(lmk).nInlier + 1;

                % 6. LANDMARK CORRECTION AND 7. REPARAMETRIZATION
                % fully correct landmark - EKF, reparam. and off-filter
                [Rob,Sen,Lmk(lmk),Obs(lmk)] = correctLmk(...
                    Rob,      ...
                    Sen,      ...
                    Lmk(lmk), ...
                    Obs(lmk), ...
                    Opt          );

            else % obs is inconsistent - do not update

                Obs(lmk).updated = false;

            end % if consistent

        end % if matched

    end % for lmk = lmkList

    % 8. LANDMARK DELETION -- loop all visible
    for lmk = find([Obs.vis])

        [Lmk(lmk),Obs(lmk)] = smartDeleteLmk(Lmk(lmk),Obs(lmk));

    end

    % 9. FORCE COVARIANCE SYMMETRY
    Map.P(Map.used,Map.used) = (Map.P(Map.used,Map.used) + Map.P(Map.used,Map.used)')/2;

end % if any(vis)

end % function

% ========== 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:
classdef DesignCapturerOfMaterialDesignExperiment < handle

    properties (Access = public)

    end

    properties (Access = private)
        outputFigureName
        inputFile
        experimentPath
        pathInput
        iterations
    end

    properties (Access = private)

    end

    methods (Access = public)

        function obj = DesignCapturerOfMaterialDesignExperiment()
            obj.init()
            obj.captureAllImages();
        end

    end

    methods (Access = private)

        function init(obj)
            path = '/home/alex/Dropbox/MaterialDesign/CC/Poisson';
            obj.outputFigureName = fullfile(path,'Poisson');
            obj.inputFile = 'HorizontalMaterialDesign';
            %obj.experimentPath = '/media/alex/MyPassport/MaterialDesign/CStar/LargeSquareWithPerimeter10_1';
            %obj.iterations = round(linspace(1,133,15));

         %   obj.experimentPath = '/media/alex/MyPassport/MaterialDesign/CStar/LargeCircleWithPerimeter10_12c';
         %   obj.iterations = round(linspace(1,67,15));

           % obj.experimentPath = '/media/alex/MyPassport/MaterialDesign/CStar/LargeCircleWithPerimeter10_1c';
           % obj.iterations = round(linspace(1,154,15));

            %obj.experimentPath = '/media/alex/MyPassport/MaterialDesign/CStar/HorizontalFromCircleWithPerimeter10_8_LargeCircle';
           % obj.iterations = round(linspace(1,27,9));

            obj.experimentPath = '/media/alex/MyPassport/MaterialDesign/CStar/NegPoissonNoPerimeter5x6';
            obj.iterations = round(linspace(1,162,10));

            %obj.experimentPath = '/home/alex/git-repos/Swan/Output/HorizontalMaterialDesign';
            %obj.iterations = 0;

        end

        function captureAllImages(obj)
            nIter = length(obj.iterations);
            for iter = 1:nIter
                iterImage = obj.iterations(iter);
                obj.captureImage(iterImage);
            end
        end

        function captureImage(obj,iter)
            i = iter;
            f = obj.inputFile;
            outPutNameWithIter = [obj.outputFigureName,'Iter',num2str(i)];
            inputFileName = fullfile(obj.experimentPath,[f,num2str(i),'.flavia.res']);
            s.fileName = f;
            s.outPutImageName = outPutNameWithIter;
            s.inputFileName   = inputFileName;
            imageCapturer     = GiDImageCapturer(s);
            imageCapturer.capture();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ n_data, x, fx ] = arcsinh_values ( n_data )

%*****************************************************************************80
%
%% ARCSINH_VALUES returns some values of the hyperbolic arc sine function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      ArcSinh[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
    -2.3124383412727526203, ...
    -0.88137358701954302523, ...
     0.00000000000000000000, ...
     0.099834078899207563327, ...
     0.19869011034924140647, ...
     0.29567304756342243910, ...
     0.39003531977071527608, ...
     0.48121182505960344750, ...
     0.56882489873224753010, ...
     0.65266656608235578681, ...
     0.73266825604541086415, ...
     0.80886693565278246251, ...
     0.88137358701954302523, ...
     1.4436354751788103425, ...
     1.8184464592320668235, ...
     2.0947125472611012942, ...
     2.3124383412727526203, ...
     2.9982229502979697388, ...
     5.2983423656105887574, ...
     7.6009027095419886115 ];

  x_vec = [ ...
       -5.0, ...
       -1.0, ...
        0.0, ...
        0.1, ...
        0.2, ...
        0.3, ...
        0.4, ...
        0.5, ...
        0.6, ...
        0.7, ...
        0.8, ...
        0.9, ...
        1.0, ...
        2.0, ...
        3.0, ...
        4.0, ...
        5.0, ...
       10.0, ...
      100.0, ...
     1000.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function iList = MapGroupTree2List(iGroup, iSubj, iSess, iRun)
global maingui

% Function to convert from processing element 3-tuple index to
% a linear index used to select a listbox entry

viewSetting = maingui.listboxGroupTreeParams.viewSetting;
idxs = maingui.listboxGroupTreeParams.listMaps(viewSetting).idxs;

% Convert processing element tuple index to a scalar
scalar0 = index2scalar(iGroup, iSubj, iSess, iRun);

% Find closest match to processing element argument in the listMap
for ii = 1:size(idxs, 1)
    ig = idxs(ii,1);
    is = idxs(ii,2);
    ie = idxs(ii,3);
    ir = idxs(ii,4);
    scalar1 = index2scalar(ig, is, ie, ir);
    if scalar0<=scalar1
        break;
    end
end
iList = ii;

% -----------------------------------------------------
function scalar = index2scalar(ig, is, ie, ir)
scalar = ig*1000 + is*100 + ie*10 + ir;

author: BUNPC