text:
function saveas(this,filename,format)
% Save GIfTI object in external file format
% FORMAT saveas(this,filename,format)
% this      - GIfTI object
% filename  - name of file to be created [Default: 'untitled.vtk']
% format    - optional argument to specify encoding format, among
%             VTK (.vtk,.vtp), Collada (.dae), IDTF (.idtf), Wavefront OBJ
%             (.obj), JavaScript (.js), JSON (.json), FreeSurfer
%             (.surf,.curv), MZ3 (.mz3) [Default: VTK]
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: saveas.m 7676 2019-10-22 10:19:29Z guillaume $

% Check filename and file format
%--------------------------------------------------------------------------
ext = '.vtk';
if nargin == 1
    filename = ['untitled' ext];
else
    if nargin == 3 && strcmpi(format,'js')
        ext = '.js';
    end
    if nargin == 3 && strcmpi(format,'json')
        ext = '.json';
    end
    if nargin == 3 && strcmpi(format,'collada')
        ext = '.dae';
    end
    if nargin == 3 && strcmpi(format,'idtf')
        ext = '.idtf';
    end
    if nargin == 3 && strncmpi(format,'vtk',3)
        format = lower(format(5:end));
        ext = '.vtk';
    end
    if nargin == 3 && strcmpi(format,'obj')
        ext = '.obj';
    end
    if nargin == 3 && strcmpi(format,'ply')
        ext = '.ply';
    end
    if nargin == 3 && strcmpi(format,'curv')
        ext = '.curv';
    end
    if nargin == 3 && strcmpi(format,'surf')
        ext = '.surf';
    end
    if nargin == 3 && strcmpi(format,'mz3')
        ext = '.mz3';
    end
    [p,f,e] = fileparts(filename);
    if strcmpi(e,'.gii')
        warning('Use save instead of saveas.');
        save(this,filename);
        return;
    end
    if ~ismember(lower(e),{ext})
        warning('Changing file extension from %s to %s.',e,ext);
        e = ext;
    end
    filename = fullfile(p,[f e]);
end

% Write file
%--------------------------------------------------------------------------
s = struct(this);

switch ext
    case {'.js','.json'}
        save_js(s,filename);
    case '.dae'
        save_dae(s,filename);
    case '.idtf'
        save_idtf(s,filename);
    case {'.vtk','.vtp'}
        if nargin < 3, format = 'legacy-ascii'; end
        mvtk_write(s,filename,format);
    case '.obj'
        save_obj(s,filename);
    case '.ply'
        save_ply(s,filename);
    case {'.curv','.surf'}
        write_freesurfer_file(s,filename,format);
    case '.mz3'
        mz3_write(s,filename,true);
    otherwise
        error('Unknown file format.');
end

%==========================================================================
% function save_js(s,filename)
%==========================================================================
function save_js(s,filename)

% Vertices & faces
%--------------------------------------------------------------------------
trace = struct(...
    'type','mesh3d',...
    'x',s.vertices(:,1),...
    'y',s.vertices(:,2),...
    'z',s.vertices(:,3),...
    'i',s.faces(:,1)-1,...
    'j',s.faces(:,2)-1,...
    'k',s.faces(:,3)-1);
if isfield(s,'cdata') && ~isempty(s.cdata)
    if size(s.cdata,2) == 1
        trace.intensity = s.cdata(:,1);
    elseif size(s.cdata,2) == 3
        trace.vertexcolor = cell(1,size(trace.x,1));
        for i=1:numel(trace.vertexcolor)
            trace.vertexcolor{i} = sprintf('rgb(%d,%d,%d)',floor(s.cdata(i,:)*255));
        end
    end
end
data = {trace};

%-Save as JS variable or JSON file
%--------------------------------------------------------------------------
if filename(end) == 's' % JS
    if exist('jsonencode','builtin')
        data = jsonencode(data);
    else
        data = spm_jsonwrite(data);
    end
    fid = fopen(filename,'wt');
    if fid == -1
        error('Unable to write file %s: permission denied.',filename);
    end
    fprintf(fid,'/*\n * JavaScript file saved by %s\n */\n\n','SPM'); % spm('Version')
    fprintf(fid,'/*\n  <!DOCTYPE html><meta charset="utf-8"/>\n');
    fprintf(fid,'  \n','https://cdn.plot.ly/plotly-latest.min.js');
    fprintf(fid,'  \n');
    fprintf(fid,'

\n','plotly',650,800);
    fprintf(fid,'  \n');
    fprintf(fid,'  \n','plotly');
    fprintf(fid,'  \n*/\n');
    fclose(fid);
else % JSON
    if exist('jsonencode','builtin')
        fid = fopen(filename,'wt');
        if fid == -1
            error('Unable to write file %s: permission denied.',filename);
        end
        fprintf(fid,'%s',jsonencode(struct('data',{data},'layout',struct())));
        fclose(fid);
    else
        spm_jsonwrite(filename,struct('data',{data},'layout',struct()),struct('indent',' '));
    end
end

%==========================================================================
% function save_dae(s,filename)
%==========================================================================
function save_dae(s,filename)

o = @(x) blanks(x*3);

% Split the mesh into connected components
%--------------------------------------------------------------------------
try
    C = spm_mesh_label(s.faces);
    d = [];
    for i=1:numel(unique(C))
        d(i).faces    = s.faces(C==i,:);
        u             = unique(d(i).faces);
        d(i).vertices = s.vertices(u,:);
        a             = 1:max(d(i).faces(:));
        a(u)          = 1:size(d(i).vertices,1);
        %a = sparse(1,double(u),1:1:size(d(i).vertices,1));
        d(i).faces    = a(d(i).faces);
    end
    s = d;
end

% Open file for writing
%--------------------------------------------------------------------------
fid = fopen(filename,'wt');
if fid == -1
    error('Unable to write file %s: permission denied.',filename);
end

% Prolog & root of the Collada XML file
%--------------------------------------------------------------------------
fprintf(fid,'<?xml version="1.0"?>\n');
fprintf(fid,'<COLLADA xmlns="http://www.collada.org/2008/03/COLLADASchema" version="1.5.0">\n');

% Assets
%--------------------------------------------------------------------------
fprintf(fid,'%s<asset>\n',o(1));
fprintf(fid,'%s<contributor>\n',o(2));
fprintf(fid,'%s<author_website>%s</author_website>\n',o(3),...
    'http://www.fil.ion.ucl.ac.uk/spm/');
fprintf(fid,'%s<authoring_tool>%s</authoring_tool>\n',o(3),'SPM');
fprintf(fid,'%s</contributor>\n',o(2));
fprintf(fid,'%s<created>%s</created>\n',o(2),datestr(now,'yyyy-mm-ddTHH:MM:SSZ'));
fprintf(fid,'%s<modified>%s</modified>\n',o(2),datestr(now,'yyyy-mm-ddTHH:MM:SSZ'));
fprintf(fid,'%s<unit name="millimeter" meter="0.001"/>\n',o(2));
fprintf(fid,'%s<up_axis>Z_UP</up_axis>\n',o(2));
fprintf(fid,'%s</asset>\n',o(1));

% Image, Materials, Effects
%--------------------------------------------------------------------------
%fprintf(fid,'%s<library_images/>\n',o(1));

fprintf(fid,'%s<library_materials>\n',o(1));
for i=1:numel(s)
    fprintf(fid,'%s<material id="material%d" name="material%d">\n',o(2),i,i);
    fprintf(fid,'%s<instance_effect url="#material%d-effect"/>\n',o(3),i);
    fprintf(fid,'%s</material>\n',o(2));
end
fprintf(fid,'%s</library_materials>\n',o(1));

fprintf(fid,'%s<library_effects>\n',o(1));
for i=1:numel(s)
    fprintf(fid,'%s<effect id="material%d-effect" name="material%d-effect">\n',o(2),i,i);
    fprintf(fid,'%s<profile_COMMON>\n',o(3));
    fprintf(fid,'%s<technique sid="COMMON">\n',o(4));
    fprintf(fid,'%s<lambert>\n',o(5));
    fprintf(fid,'%s<emission>\n',o(6));
    fprintf(fid,'%s<color>%f %f %f %d</color>\n',o(7),[0 0 0 1]);
    fprintf(fid,'%s</emission>\n',o(6));
    fprintf(fid,'%s<ambient>\n',o(6));
    fprintf(fid,'%s<color>%f %f %f %d</color>\n',o(7),[0 0 0 1]);
    fprintf(fid,'%s</ambient>\n',o(6));
    fprintf(fid,'%s<diffuse>\n',o(6));
    fprintf(fid,'%s<color>%f %f %f %d</color>\n',o(7),[0.5 0.5 0.5 1]);
    fprintf(fid,'%s</diffuse>\n',o(6));
    fprintf(fid,'%s<transparent>\n',o(6));
    fprintf(fid,'%s<color>%d %d %d %d</color>\n',o(7),[1 1 1 1]);
    fprintf(fid,'%s</transparent>\n',o(6));
    fprintf(fid,'%s<transparency>\n',o(6));
    fprintf(fid,'%s<float>%f</float>\n',o(7),0);
    fprintf(fid,'%s</transparency>\n',o(6));
    fprintf(fid,'%s</lambert>\n',o(5));
    fprintf(fid,'%s</technique>\n',o(4));
    fprintf(fid,'%s</profile_COMMON>\n',o(3));
    fprintf(fid,'%s</effect>\n',o(2));
end
fprintf(fid,'%s</library_effects>\n',o(1));

% Geometry
%--------------------------------------------------------------------------
fprintf(fid,'%s<library_geometries>\n',o(1));
for i=1:numel(s)
    fprintf(fid,'%s<geometry id="shape%d" name="shape%d">\n',o(2),i,i);
    fprintf(fid,'%s<mesh>\n',o(3));
    fprintf(fid,'%s<source id="shape%d-positions">\n',o(4),i);
    fprintf(fid,'%s<float_array id="shape%d-positions-array" count="%d">',o(5),i,numel(s(i).vertices));
    fprintf(fid,'%f ',reshape(s(i).vertices',1,[]));
    fprintf(fid,'</float_array>\n');
    fprintf(fid,'%s<technique_common>\n',o(5));
    fprintf(fid,'%s<accessor count="%d" offset="0" source="#shape%d-positions-array" stride="3">\n',o(6),size(s(i).vertices,1),i);
    fprintf(fid,'%s<param name="X" type="float" />\n',o(7));
    fprintf(fid,'%s<param name="Y" type="float" />\n',o(7));
    fprintf(fid,'%s<param name="Z" type="float" />\n',o(7));
    fprintf(fid,'%s</accessor>\n',o(6));
    fprintf(fid,'%s</technique_common>\n',o(5));
    fprintf(fid,'%s</source>\n',o(4));
    fprintf(fid,'%s<vertices id="shape%d-vertices">\n',o(4),i);
    fprintf(fid,'%s<input semantic="POSITION" source="#shape%d-positions"/>\n',o(5),i);
    fprintf(fid,'%s</vertices>\n',o(4));
    fprintf(fid,'%s<triangles material="material%d" count="%d">\n',o(4),i,size(s(i).faces,1));
    fprintf(fid,'%s<input semantic="VERTEX" source="#shape%d-vertices" offset="0"/>\n',o(5),i);
    fprintf(fid,'%s
',o(5));
    fprintf(fid,'%d ',reshape(s(i).faces',1,[])-1);
    fprintf(fid,'
\n');
    fprintf(fid,'%s</triangles>\n',o(4));
    fprintf(fid,'%s</mesh>\n',o(3));
    fprintf(fid,'%s</geometry>\n',o(2));
end
fprintf(fid,'%s</library_geometries>\n',o(1));

% Scene
%--------------------------------------------------------------------------
fprintf(fid,'%s<library_visual_scenes>\n',o(1));
fprintf(fid,'%s<visual_scene id="VisualSceneNode" name="SceneNode">\n',o(2));
for i=1:numel(s)
    fprintf(fid,'%s<node id="node%d">\n',o(3),i);
    fprintf(fid,'%s<instance_geometry url="#shape%d">\n',o(4),i);
    fprintf(fid,'%s<bind_material>\n',o(5));
    fprintf(fid,'%s<technique_common>\n',o(6));
    fprintf(fid,'%s<instance_material symbol="material%d" target="#material%d"/>\n',o(7),i,i);
    fprintf(fid,'%s</technique_common>\n',o(6));
    fprintf(fid,'%s</bind_material>\n',o(5));
    fprintf(fid,'%s</instance_geometry>\n',o(4));
    fprintf(fid,'%s</node>\n',o(3));
end
fprintf(fid,'%s</visual_scene>\n',o(2));
fprintf(fid,'%s</library_visual_scenes>\n',o(1));
fprintf(fid,'%s<scene>\n',o(1));
fprintf(fid,'%s<instance_visual_scene url="#VisualSceneNode" />\n',o(2));
fprintf(fid,'%s</scene>\n',o(1));

% End of XML
%--------------------------------------------------------------------------
fprintf(fid,'</COLLADA>\n');

% Close file
%--------------------------------------------------------------------------
fclose(fid);

%==========================================================================
% function save_idtf(s,filename)
%==========================================================================
function save_idtf(s,filename)

o = @(x) blanks(x*3);

% Compute normals
%--------------------------------------------------------------------------
if ~isfield(s,'normals')
    try
        s.normals = spm_mesh_normals(...
            struct('vertices',s.vertices,'faces',s.faces),true);
    catch
        s.normals = [];
    end
end

% Split the mesh into connected components
%--------------------------------------------------------------------------
try
    C = spm_mesh_label(s.faces);
    d = [];
    try
        if size(s.cdata,2) == 1 && (any(s.cdata>1) || any(s.cdata<0))
            mi = min(s.cdata); ma = max(s.cdata);
            s.cdata = (s.cdata-mi)/ (ma-mi);
        else
        end
    end
    for i=1:numel(unique(C))
        d(i).faces    = s.faces(C==i,:);
        u             = unique(d(i).faces);
        d(i).vertices = s.vertices(u,:);
        d(i).normals  = s.normals(u,:);
        a             = 1:max(d(i).faces(:));
        a(u)          = 1:size(d(i).vertices,1);
        %a = sparse(1,double(u),1:1:size(d(i).vertices,1));
        d(i).faces    = a(d(i).faces);
        d(i).mat      = s.mat;
        try
            d(i).cdata = s.cdata(u,:);
            if size(d(i).cdata,2) == 1
                d(i).cdata = repmat(d(i).cdata,1,3);
            end
        end
    end
    s = d;
end

% Open file for writing
%--------------------------------------------------------------------------
fid = fopen(filename,'wt');
if fid == -1
    error('Unable to write file %s: permission denied.',filename);
end

% FILE_HEADER
%--------------------------------------------------------------------------
fprintf(fid,'FILE_FORMAT "IDTF"\n');
fprintf(fid,'FORMAT_VERSION 100\n\n');

% NODES
%--------------------------------------------------------------------------
for i=1:numel(s)
    fprintf(fid,'NODE "MODEL" {\n');
    fprintf(fid,'%sNODE_NAME "%s"\n',o(1),sprintf('Mesh%04d',i));
    fprintf(fid,'%sPARENT_LIST {\n',o(1));
    fprintf(fid,'%sPARENT_COUNT %d\n',o(2),1);
    fprintf(fid,'%sPARENT %d {\n',o(2),0);
    fprintf(fid,'%sPARENT_NAME "%s"\n',o(3),'<NULL>');
    fprintf(fid,'%sPARENT_TM {\n',o(3));
    I = s(i).mat; % eye(4);
    for j=1:size(I,2)
        fprintf(fid,'%s',o(4)); fprintf(fid,'%f ',I(:,j)'); fprintf(fid,'\n');
    end
    fprintf(fid,'%s}\n',o(3));
    fprintf(fid,'%s}\n',o(2));
    fprintf(fid,'%s}\n',o(1));
    fprintf(fid,'%sRESOURCE_NAME "%s"\n',o(1),sprintf('Mesh%04d',i));
    %fprintf(fid,'%sMODEL_VISIBILITY "BOTH"\n',o(1));
    fprintf(fid,'}\n\n');
end

% NODE_RESOURCES
%--------------------------------------------------------------------------
for i=1:numel(s)
    fprintf(fid,'RESOURCE_LIST "MODEL" {\n');
    fprintf(fid,'%sRESOURCE_COUNT %d\n',o(1),1);
    fprintf(fid,'%sRESOURCE %d {\n',o(1),0);
    fprintf(fid,'%sRESOURCE_NAME "%s"\n',o(2),sprintf('Mesh%04d',i));
    fprintf(fid,'%sMODEL_TYPE "MESH"\n',o(2));
    fprintf(fid,'%sMESH {\n',o(2));
    fprintf(fid,'%sFACE_COUNT %d\n',o(3),size(s(i).faces,1));
    fprintf(fid,'%sMODEL_POSITION_COUNT %d\n',o(3),size(s(i).vertices,1));
    fprintf(fid,'%sMODEL_NORMAL_COUNT %d\n',o(3),size(s(i).normals,1));
    if ~isfield(s(i),'cdata') || isempty(s(i).cdata)
        c = 0;
    else
        c = size(s(i).cdata,1);
    end
    fprintf(fid,'%sMODEL_DIFFUSE_COLOR_COUNT %d\n',o(3),c);
    fprintf(fid,'%sMODEL_SPECULAR_COLOR_COUNT %d\n',o(3),0);
    fprintf(fid,'%sMODEL_TEXTURE_COORD_COUNT %d\n',o(3),0);
    fprintf(fid,'%sMODEL_BONE_COUNT %d\n',o(3),0);
    fprintf(fid,'%sMODEL_SHADING_COUNT %d\n',o(3),1);
    fprintf(fid,'%sMODEL_SHADING_DESCRIPTION_LIST {\n',o(3));
    fprintf(fid,'%sSHADING_DESCRIPTION %d {\n',o(4),0);
    fprintf(fid,'%sTEXTURE_LAYER_COUNT %d\n',o(5),0);
    fprintf(fid,'%sSHADER_ID %d\n',o(5),0);
    fprintf(fid,'%s}\n',o(4));
    fprintf(fid,'%s}\n',o(3));

    fprintf(fid,'%sMESH_FACE_POSITION_LIST {\n',o(3));
    fprintf(fid,'%d %d %d\n',s(i).faces'-1);
    fprintf(fid,'%s}\n',o(3));

    fprintf(fid,'%sMESH_FACE_NORMAL_LIST {\n',o(3));
    fprintf(fid,'%d %d %d\n',s(i).faces'-1);
    fprintf(fid,'%s}\n',o(3));

    fprintf(fid,'%sMESH_FACE_SHADING_LIST {\n',o(3));
    fprintf(fid,'%d\n',zeros(size(s(i).faces,1),1));
    fprintf(fid,'%s}\n',o(3));

    if c
        fprintf(fid,'%sMESH_FACE_DIFFUSE_COLOR_LIST {\n',o(3));
        fprintf(fid,'%d %d %d\n',s(i).faces'-1);
        fprintf(fid,'%s}\n',o(3));
    end

    fprintf(fid,'%sMODEL_POSITION_LIST {\n',o(3));
    fprintf(fid,'%f %f %f\n',s(i).vertices');
    fprintf(fid,'%s}\n',o(3));

    fprintf(fid,'%sMODEL_NORMAL_LIST {\n',o(3));
    fprintf(fid,'%f %f %f\n',s(i).normals');
    fprintf(fid,'%s}\n',o(3));

    if c
        fprintf(fid,'%sMODEL_DIFFUSE_COLOR_LIST {\n',o(3));
        fprintf(fid,'%f %f %f\n',s(i).cdata');
        fprintf(fid,'%s}\n',o(3));
    end

    fprintf(fid,'%s}\n',o(2));
    fprintf(fid,'%s}\n',o(1));
    fprintf(fid,'}\n');
end

% Close file
%--------------------------------------------------------------------------
fclose(fid);

%==========================================================================
% function save_obj(s,filename)
%==========================================================================
function save_obj(s,filename)

% Open file for writing
%--------------------------------------------------------------------------
fid = fopen(filename,'wt');
if fid == -1
    error('Unable to write file %s: permission denied.',filename);
end

% Vertices & faces
%--------------------------------------------------------------------------
fprintf(fid,'# Wavefront OBJ file saved by %s\n','SPM'); % spm('Version')
fprintf(fid,'v %f %f %f\n',s.vertices');
fprintf(fid,'f %d %d %d\n',s.faces');

% Close file
%--------------------------------------------------------------------------
fclose(fid);

%==========================================================================
% function save_ply(s,filename)
%==========================================================================
function save_ply(s,filename)

% Open file for writing
%--------------------------------------------------------------------------
fid = fopen(filename,'wt');
if fid == -1
    error('Unable to write file %s: permission denied.',filename);
end

% Header
%--------------------------------------------------------------------------
fprintf(fid,'ply\n');
fprintf(fid,'format ascii 1.0\n');
fprintf(fid,'element vertex %d\n',size(s.vertices,1));
fprintf(fid,'property float32 x\n');
fprintf(fid,'property float32 y\n');
fprintf(fid,'property float32 z\n');
fprintf(fid,'element face %d\n',size(s.faces,1));
fprintf(fid,'property list uint8 int32 vertex_indices\n');
fprintf(fid,'end_header\n');

% Vertices & faces
%--------------------------------------------------------------------------
fprintf(fid,'%f %f %f\n',s.vertices');
fprintf(fid,'3 %d %d %d\n',s.faces'-1);

% Close file
%--------------------------------------------------------------------------
fclose(fid);

%==========================================================================
% function write_freesurfer_file(s,filename,format)
%==========================================================================
function write_freesurfer_file(s,filename,format)

switch lower(format)
    case 'surf'
        fid = fopen(filename,'wb','ieee-be');
        if fid == -1, error('Cannot open "%s".',filename); end
        fwrite(fid,[255 255 254],'uchar');
        fprintf(fid,'created by @gifti on %s\n\n',...
            datestr(now,'ddd mmm dd HH:MM:SS yyyy'));
        fwrite(fid,size(s.vertices,1),'int32'); % nv
        fwrite(fid,size(s.faces,1),'int32'); % nf
        fwrite(fid,s.vertices','float32');
        fwrite(fid,s.faces'-1,'int32');
        fclose(fid);
    case 'curv'
        fid = fopen(filename,'wb','ieee-be');
        if fid == -1, error('Cannot open "%s".',filename); end
        fwrite(fid,[255 255 255],'uchar');
        fwrite(fid,size(s.cdata,1),'int32'); % nv
        if isfield(s,'faces')
            nf = size(s.faces,1);
        else
            nf = 0;
        end
        fwrite(fid,nf,'int32'); % nf
        fwrite(fid,size(s.cdata,2),'int32'); % n
        fwrite(fid,s.cdata','float32');
        fclose(fid);
    otherwise
        error('File format not supported.');
end


text:
function [m,n] = size(a,dim)
%SIZE         Implements  size(a)  for gradients
%
%   [m,n] = size(a,dim)
%
% functionality as Matlab function size, second argument dim optional
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargout==2
    [m n] = size(a.x);
  else
    if nargin==1
      m = size(a.x);
    else
      m = size(a.x,dim);
    end
  end

repo: OpenMAS


text:
function option_num = p20_option_num ( )

%*****************************************************************************80
%
%% P20_OPTION_NUM returns the number of options for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer OPTION_NUM, the number of options.
%
  option_num = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function grad = B_beamforming_power(future_layer, curr_layer, prev_layer)
% we don't implement this function now, as it will be used with
% beamforming. We will implement the joint gradient of Power spectrum and
% frequency domain beamforming.
% X is the multichannel complex spectrum inputs
future_grad = future_layer{1}.grad;
[X, weight] =  prepareBeamforming(prev_layer);
[C,T,N] = size(X);
% Y is the beamforming's output
Y = curr_layer.a;

if 0
    for i=1:C
        cross_term = future_grad .* Y .* squeeze(conj(X(i,:,:)))';
        grad(:,i) = 2*sum(cross_term,2);
    end
else
    future_grad_Y = future_grad .* Y;
    X_tmp = permute(X, [3 2 1]);
    cross_term = bsxfun(@times, X_tmp, future_grad_Y);
    grad = 2*squeeze(sum(cross_term,2));
end

end

author: singaxiong

repo: SignalGraph


text:
function [cubR,cubS,cubW, Ncub] = Cubature2D(Corder)

% function [cubR,cubS,cubW, Ncub] = Cubature2D(Corder)
% Purpose: provide multidimensional quadrature (i.e. cubature)
%          rules to integrate up to Corder polynomials

CubatureData2D;

if(Corder<=28)
  cubR = cub2D{Corder}(:,1);
  cubS = cub2D{Corder}(:,2);
  cubW = cub2D{Corder}(:,3);
else
  cubNA = ceil( (Corder+1)/2);
  [cubA,cubWA] = JacobiGQ(0,0, cubNA-1);
  cubNB = ceil( (Corder+1)/2);
  [cubB,cubWB] = JacobiGQ(1,0, cubNB-1);

  cubA = ones(cubNB,1)*cubA';
  cubB = cubB*ones(1,cubNA);

  cubR = 0.5*(1+cubA).*(1-cubB)-1;
  cubS = cubB;
  cubW = 0.5*cubWB*(cubWA');

  cubR = cubR(:);
  cubS = cubS(:);
  cubW = cubW(:);
end
Ncub = length(cubW(:));
return

repo: Aero-matlab


text:
function N=nmultichoosek(n,k,type)
%NMULTICHOOSEK   Unordered samples WITH or WITHOUT repetition.
%   NMULTICHOOSEK(N,K) finds the number of multisets of length k on n
%   symbols. NMULTICHOOSEK can take vector or scalar input.
%
%   NMULTICHOOSEK(N,K,'single') is the same as NCHOOSEK (unordered samples WITHOUT
%    repetition), except that it accepts vector inputs for both n and k.
%
%   NMULTICHOOSEK(N,K,'multi') is the same as NMULTICHOOSEK(N,K).
%
%   Examples:
%      N = nmultichoosek(5,1:5)
%
%   finds the number of multisets of length 1 to 5 from a 5 symbol set
%
%      N = nmultichoosek(5,1:5,'single')
%
%   is the same as:
%
%      for k=1:5
%         N(k) = nchoosek(5,k);
%      end
%
%   See also nchoosek, perms.

% Author: Peter (PB) Bodin
% Email: pbodin@kth.se
% Created: 14-Nov-2005 19:48:49

msg = nargchk(2,3,nargin);
error(msg);
n = n(:);           % Convert N-D arrays to 1-D
k = k(:);           % Convert N-D arrays to 1-D
if numel(n)>1 & numel(k)>1
    if numel(n)~=numel(k)
        error('vectors must have the same number of elements')
    end
end

if nargin > 2
    switch type
        case 'single'
            N = round(exp(gammaln(n+1) - gammaln(k+1) - gammaln(n-k+1)));
            return;
        case 'multi'
        otherwise
            error('valid types are ''multi'' and ''single''');
    end
end

N = round(exp(gammaln(n+k) - gammaln(k+1) - gammaln(n)));
N(n<k)=NaN;

author: Sable

repo: mcbench-benchmarks


text:
function dimacs = computedimacs(b,c,A,xin,y,s,K);
% COMPUTEDIMACS
%
% min <C,X> s.t     AX = b, X > 0
% max b'y   s.t S-C+A'y =0, S > 0

% If no primal exist, fake till later
if isempty(xin)
    x = c*0;
else
    x = xin;
end

if isempty(s)
    s = c-A'*y;
end

xres = inf;
sres = inf;

% Not officially defined in DIMACS
if K.f>0
    sres = -min(norm(s(1:K.f),inf));
end

% Errors in linear cone
if K.l>0
    xres = min(x(1+K.f:K.f+K.l));
    sres = min(s(1+K.f:K.f+K.l));
end

% Errors in quadratic cone
if any(K.q)
    top = K.f+K.l;
    for i = 1:length(K.q)
        X = x(1+top:top+K.q(i));
        S = s(1+top:top+K.q(i));
        xres = min(xres,X(1)-norm(X(2:end)));
        sres = min(sres,S(1)-norm(S(2:end)));
        top = top + K.q(i);
    end
end

% Errors in semidefinite cone
if any(K.s)
    top = K.f+K.l+K.q+K.r;
    for i = 1:length(K.s)
        X = reshape(x(1+top:top+K.s(i)^2),K.s(i),K.s(i));
        S = reshape(s(1+top:top+K.s(i)^2),K.s(i),K.s(i));
        xres = min(xres,min(eig(full(X))));
        sres = min(sres,min(eig(full(S))));
        top = top + K.s(i)^2;
    end
end

err1 = norm(b-A*x)/(1 + norm(b,inf));
err2 = max(0,-xres)/(1 + norm(b,inf));
err3 = conenorm(s-(c-A'*y),K)/(1+norm(c,inf));
%err3 = norm(s-(c-A'*y))/(1+norm(c,inf)); % Used by some solvers
err4 = max(0,-sres)/(1+max(abs(c)));
err5 = (c'*x-b'*y)/(1+abs(c'*x)+abs(b'*y));
err6 = x'*(c-A'*y)/(1+abs(c'*x)+abs(b'*y));

% No primal was computed
if isempty(xin)
    err1 = nan;
    err2 = nan;
    err5 = nan;
    err6 = nan;
end
dimacs = [err1 err2 err3 err4 err5 err6];

function t = conenorm(s,K)

% Implementation of the norm described on
% http://plato.asu.edu/dimacs/node3.html

t = 0;

if K.f + K.l>0
    t = t + norm(s(1:K.f+K.l));
end

top = 1+K.f+K.l;
if any(K.q)
    for i = 1:length(K.q)
        t = t + norm(s(top:top+K.q(i)-1));
        top  = top + K.q(i);
    end
end

if any(K.s)
    for i = 1:length(K.s)
        S = reshape(s(top:top+K.s(i)^2-1),K.s(i),K.s(i));
        t = t + norm(S,'fro');
        top  = top + K.s(i)^2;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function y=rdct(x,n)
%RDCT     Discrete cosine transform of real data Y=(X,N)
% Data is truncated/padded to length N.
%
% This routine is equivalent to multiplying by the matrix
%
%   rdct(eye(n)) = diag([sqrt(2) 2*ones(1,n-1)]) * cos((0:n-1)'*(0.5:n)*pi/n)
%
% The rows and columns of the matrix are orthogonal but not unit modulus.
% Various versions of the DCT are obtained by pre-multiplying the above
% matrix by diag([b/a ones(1,n-1)/a]) and post-multiplying the
% inverse transform matrix by its inverse. A common choice is a=n and/or b=sqrt(2).
% Choose a=sqrt(2n) and b=1 to make the matrix orthogonal.
% If b~=1 then the columns are no longer orthogonal.
%
% see IRDCT for the inverse transform

%      Copyright (C) Mike Brookes 1998
%
%      Last modified Tue Apr 13 15:56:48 1999
%
%   VOICEBOX is a MATLAB toolbox for speech processing. Home page is at
%   http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fl=size(x,1)==1;
if fl x=x(:); end
[m,k]=size(x);
if nargin<2 n=m;
elseif n>m x=[x; zeros(n-m,k)];
elseif n<m x(n+1:m,:)=[];
end

x=[x(1:2:n,:); x(2*fix(n/2):-2:2,:)];
z=[sqrt(2) 2*exp((-0.5i*pi/n)*(1:n-1))].';
y=real(fft(x).*z(:,ones(1,k)));

if fl y=y.'; end

author: bastamon


text:
function varargout = FAULT(varargin)
% FAULT M-file for FAULT.fig
%      FAULT, by itself, creates a new FAULT or raises the existing
%      singleton*.
%
%      H = FAULT returns the handle to a new FAULT or the handle to
%      the existing singleton*.
%
%      FAULT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in FAULT.M with the given input arguments.
%
%      FAULT('Property','Value',...) creates a new FAULT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before FAULT_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to FAULT_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 FAULT

% Last Modified by GUIDE v2.5 28-Nov-2012 12:06:59

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @FAULT_OpeningFcn, ...
                   'gui_OutputFcn',  @FAULT_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 FAULT is made visible.
function FAULT_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 FAULT (see VARARGIN)

% Choose default command line output for FAULT
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = FAULT_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} = handles.output;

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (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 edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (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 edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (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 edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (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 edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (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 edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (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 edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (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 edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double

% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (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 edit5_Callback(hObject, eventdata, handles)
% hObject    handle to edit5 (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 edit5 as text
%        str2double(get(hObject,'String')) returns contents of edit5 as a double

% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (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 edit6_Callback(hObject, eventdata, handles)
% hObject    handle to edit6 (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 edit6 as text
%        str2double(get(hObject,'String')) returns contents of edit6 as a double

% --- Executes during object creation, after setting all properties.
function edit6_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit6 (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 edit8_Callback(hObject, eventdata, handles)
% hObject    handle to edit8 (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 edit8 as text
%        str2double(get(hObject,'String')) returns contents of edit8 as a double

% --- Executes during object creation, after setting all properties.
function edit8_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit8 (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 edit9_Callback(hObject, eventdata, handles)
% hObject    handle to edit9 (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 edit9 as text
%        str2double(get(hObject,'String')) returns contents of edit9 as a double

% --- Executes during object creation, after setting all properties.
function edit9_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit9 (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 button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
a=(get(handles.popupmenu1,'value'));
real=str2num(get(handles.edit1,'string'));
img=str2num(get(handles.edit2,'string'));
Z0= real+img*i;
real1=str2num(get(handles.edit3,'string'));
img1=str2num(get(handles.edit4,'string'));
Z1=real1+img1*i;
real2=str2num(get(handles.edit5,'string'));
img2=str2num(get(handles.edit6,'string'));
Z2=real2+img2*i;
line_voltage=str2num(get(handles.edit9,'string'));
phase_voltage=line_voltage/(sqrt(3));
kva=str2num(get(handles.edit8,'string'));
% for popup menu commond
a=(get(handles.popupmenu1,'value'));switch a
    case 1
       %LINE TO GROUND FAULT
Io=(phase_voltage)/(Z0+Z1+Z2);
fault_current=(3*(phase_voltage)/(Z0+Z1+Z2));
set(handles.edit10,'string',imag(fault_current));
    case 2
   %LINE TO LINE FAULT
fault_current2=(-i*sqrt(3)*phase_voltage)/(Z1+Z1);
set(handles.edit10,'string',(fault_current2));
end

function edit10_Callback(hObject, eventdata, handles)
% hObject    handle to edit10 (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 edit10 as text
%        str2double(get(hObject,'String')) returns contents of edit10 as a double

% --- Executes during object creation, after setting all properties.
function edit10_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit10 (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 edit11_Callback(hObject, eventdata, handles)
% hObject    handle to edit11 (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 edit11 as text
%        str2double(get(hObject,'String')) returns contents of edit11 as a double

% --- Executes during object creation, after setting all properties.
function edit11_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit11 (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 popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (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 popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1

% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (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 drawbtf (A, p, q, r)
%DRAWBTF plot the BTF form of a matrix
%
% A(p,q) is in BTF form, r the block boundaries
%
% Example:
%   [p,q,r] = dmperm (A)
%   drawbtf (A, p, q, r)
%
% See also btf, maxtrans, strongcomp, dmperm.

% Copyright 2004-2007, Tim Davis, University of Florida

nblocks = length (r) - 1 ;

hold off
spy (A (p,abs(q)))
hold on

for k = 1:nblocks
    k1 = r (k) ;
    k2 = r (k+1) ;
    nk = k2 - k1 ;
    if (nk > 1)
        plot ([k1 k2 k2 k1 k1]-.5, [k1 k1 k2 k2 k1]-.5, 'r') ;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = cyl2cart(varargin)
%CYL2CART  Convert cylindrical to cartesian coordinates
%
%   CART = cyl2cart(CYL)
%   convert the 3D cylindrical coordinates of points CYL (given by
%   [THETA R Z] where THETA, R, and Z have the same size) into cartesian
%   coordinates CART, given by [X Y Z].
%   The transforms is the following :
%   X = R*cos(THETA);
%   Y = R*sin(THETA);
%   Z remains inchanged.
%
%   CART = cyl2cart(THETA, R, Z)
%   provides coordinates as 3 different parameters
%
%   Example
%   cyl2cart([-1 0 2])
%   gives : 4.7124    1.0000     2.0000
%
%   See also angles3d, cart2pol, cart2sph2, cart2cyl
%
%
% ------
% Author: David Legland
% e-mail: david.legland@jouy.inra.fr
% Created: 2006-03-23
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% process input parameters
if length(varargin)==1
    var = varargin{1};
    theta = var(:,1);
    r = var(:,2);
    z = var(:,3);
elseif length(varargin)==3
    theta = varargin{1};
    r = varargin{2};
    z = varargin{3};
end

% convert coordinates
dim = size(theta);
x = reshape(r(:).*cos(theta(:)), dim);
y = reshape(r(:).*sin(theta(:)), dim);

% process output parameters
if nargout==0 ||nargout==1
    if length(dim)>2 || dim(2)>1
        varargout{1} = {x y z};
    else
        varargout{1} = [x y z];
    end
elseif nargout==3
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
classdef NomadOptimizer < AbstractOptimizer
    %NomadOptimizer Summary of this class goes here
    %   Detailed explanation goes here

    properties
        options(1,1) NomadOptions = NomadOptions();
    end

    methods
        function obj = NomadOptimizer()
            obj.options = NomadOptions();
        end

        function [exitflag, message] = optimize(obj, lvdOpt, writeOutput, callOutputFcn, hLvdMainGUI)
            [x0All, actVars, varNameStrs] = lvdOpt.vars.getTotalScaledXVector();
            [lbAll, ubAll, lbUsAll, ubUsAll] = lvdOpt.vars.getTotalScaledBndsVector();
            typicalX = lvdOpt.vars.getTypicalScaledXVector();
            lvdData = lvdOpt.lvdData;

            if(isempty(x0All) && isempty(actVars))
                exitflag = 0;
                message = 'No variables enabled on script.  Aborting optimization.';

                return;
            end

            evtNumToStartScriptExecAt = obj.getEvtNumToStartScriptExecAt(lvdOpt, actVars);
            evtToStartScriptExecAt = lvdOpt.lvdData.script.getEventForInd(evtNumToStartScriptExecAt);

            objFuncWrapper = @(x) lvdOpt.objFcn.evalObjFcn(x, evtToStartScriptExecAt);
            nonlcon = @(x, tfRunScript, stateLog) lvdOpt.constraints.evalConstraints(x, tfRunScript, evtToStartScriptExecAt, true, stateLog);

            opts = obj.options.getOptionsForOptimizer();
            constrTypeStr = obj.options.getConstrTypeStr();
            useParallel = obj.options.usesParallel();

            [c, ceq] = lvdOpt.constraints.evalConstraints(x0All, true, evtToStartScriptExecAt, false, []);
            numConstr = length(c) + 2*length(ceq);
            bboutput = horzcat({'OBJ'}, repmat({constrTypeStr},1,numConstr));
            numElemPerOutput = length(bboutput);
            bboutput = strjoin(bboutput,' ');
            nomadNonlconWrapper = @(x) NomadOptimizer.nomadConstrWrapper(x, nonlcon);

            pp = gcp('nocreate');
            if(isempty(pp) || useParallel == false) %
                numWorkers = 0;
                useParallel = false;
            else
                numWorkers = pp.NumWorkers;
                useParallel = true;
            end

            numVars = length(x0All);
            f = @(x) NomadOptimizer.nomadObjConstrWrapper(x, objFuncWrapper, nomadNonlconWrapper, numVars, numElemPerOutput, numWorkers, lvdData);
            opts = nomadset(opts, 'bb_output_type',bboutput);

            problem.objective = f; %f
            problem.x0 = x0All;
            problem.lb = lbAll;
            problem.ub = ubAll;
            problem.options = opts;

            problem.solver = 'nomad';

            problem.lvdData = lvdOpt.lvdData; %need to get lvdData in somehow

            %%% Run optimizer
            celBodyData = lvdOpt.lvdData.celBodyData;
            recorder = ma_OptimRecorder();

            if(callOutputFcn)
                propNames = lvdOpt.lvdData.launchVehicle.tankTypes.getFirstThreeTypesCellArr();
%                 handlesObsOptimGui = ma_ObserveOptimGUI(celBodyData, problem, true, writeOutput, [], varNameStrs, lbUsAll, ubUsAll);

                out = AppDesignerGUIOutput();
                ma_ObserveOptimGUI_App(out);
                handlesObsOptimGui = out.output{1};

%                 outputFnc = @(x, optimValues, state) ma_OptimOutputFunc(x, optimValues, state, handlesObsOptimGui, problem.objective, problem.lb, problem.ub, celBodyData, recorder, propNames, writeOutput, varNameStrs, lbUsAll, ubUsAll);
                hOptimStatusLabel = handlesObsOptimGui.optimStatusLabel;
                hFinalStateOptimLabel = handlesObsOptimGui.finalStateOptimLabel;
                hDispAxes = handlesObsOptimGui.dispAxesPanel;
                hCancelButton = handlesObsOptimGui.cancelButton;
                optimStartTic = tic();

                outputFnc = @(x, optimValues, state) NomadOptimizer.getOutputFunction(x, optimValues, state, hOptimStatusLabel, hFinalStateOptimLabel, hDispAxes, hCancelButton, ...
                                                                                problem.objective, problem.lb, problem.ub, celBodyData, recorder, propNames, writeOutput, varNameStrs, lbUsAll, ubUsAll, optimStartTic);

                nomadOutput1 = @(iter, fval, x, state) NomadOptimizer.nomadIterFunWrapper(iter, fval, x, outputFnc, state, numVars);
                nomadOutput2 = @(iter, fval, x) nomadOutput1(iter, fval, x, 'iter');
                problem.options = nomadset(problem.options, 'iterfun',nomadOutput2);
                nomadOutput1(0,NaN,x0All,'init');
            end

            problem.UseParallel = useParallel;
            [exitflag, message] = lvd_executeOptimProblem(celBodyData, writeOutput, problem, recorder, callOutputFcn);

            if(callOutputFcn)
                close(handlesObsOptimGui.ma_ObserveOptimGUI);
            end
        end

        function options = getOptions(obj)
            options = obj.options;
        end

        function tf = usesParallel(obj)
            tf = obj.getOptions().usesParallel();
        end

        function numWorkers = getNumParaWorkers(obj)
            numWorkers = obj.options.getNumParaWorkers();
        end

        function openOptionsDialog(obj)
%             lvd_editNomadOptionsGUI(obj);

            output = AppDesignerGUIOutput({false});
            lvd_editNomadOptionsGUI_App(obj, output);
        end
    end

    methods(Static, Access=private)
        function [f, stateLog] = nomadObjFuncWrapper(x, objFuncWrapper)
            global nomadCachedX nomadCachedStateLog

            [f, stateLog] = objFuncWrapper(x);

            nomadCachedX = x;
            nomadCachedStateLog = stateLog;
        end

        function c = nomadConstrWrapper(x, nonlcon)
            global nomadCachedX nomadCachedStateLog

            if(numel(x) ~= numel(nomadCachedX))
                nomadCachedX = NaN(size(x));
            end

            if(all(x(:) == nomadCachedX(:)))
                [cM, ceqM] = nonlcon(x, false, nomadCachedStateLog);
            else
                [cM, ceqM] = nonlcon(x, true, []);
            end

            c = [cM(:); ceqM(:); -1*ceqM(:)]';
            c = c(:);
        end

        function [fc, stateLogs] = nomadObjConstrWrapper(x, objFun, nonlcon, numVars, numElementsInEachOutput, M, lvdData)
            try
                if(numel(x) == numVars)
                    numEvals = 1;
                else
                    numEvals = size(x,1);
                end

                x = reshape(x, numEvals, numel(x)/numEvals);

                fc = NaN(numEvals,numElementsInEachOutput);

                for(i=1:numEvals) %#ok<NO4LP>
                    stateLogs(i) = LaunchVehicleStateLog(lvdData); %#ok<AGROW>
                end

                if(M > 0 && numEvals > 1)
                    pp = gcp();
                    opts = parforOptions(pp,'RangePartitionMethod','fixed','SubrangeSize',ceil(pp.NumWorkers/numEvals));
                    parfor(i=1:numEvals,opts)
                        xI = x(i,:);
                        [fcRow, stateLog] = NomadOptimizer.loopInternal(xI, objFun, nonlcon);

                        fc(i,:) = fcRow;
                        stateLogs(i) = stateLog;
                    end
                else
                    for(i=1:numEvals)
                        xI = x(i,:);
                        [fcRow, stateLog] = NomadOptimizer.loopInternal(xI, objFun, nonlcon);

                        fc(i,:) = fcRow;
                        stateLogs(i) = stateLog;
                    end
                end
            catch ME
                disp(ME.message);
            end
        end

        function [fcRow, stateLog] = loopInternal(xI, objFun, nonlcon)
            [f, stateLog] = objFun(xI);
            [c] = nonlcon(xI);

            fcRow = [f; c(:);]';
        end

        function stop = nomadIterFunWrapper(iter, fval, x, outputFnc, state, numVars)
            if(numel(x) == numVars)
                numEvals = 1;
            else
                numEvals = size(x,1);
            end

            x = reshape(x, numEvals, numel(x)/numEvals);

            f = fval(:,1);
            c = fval(:,2:end);

            if(isempty(c))
                [fval,I] = min(f);
                cViol = 0;
                xx = x(I,:);
            else
                c(c <= 0) = 0;
                cViol = sqrt(sum(c.^2,2));

                [minViolation,I] = min(cViol);

                if(sum(cViol == minViolation) > 1)
                    boolC = cViol == minViolation;
                    fBool = f(boolC);
                    [minFBool,~] = min(fBool);

                    boolF = f == minFBool;
                    II = find(boolF & boolC,1,'first');

                    fval = f(II);
                    cViol = minViolation;
                    xx = x(II,:);
                else
                    fval = f(I);
                    cViol = minViolation;
                    xx = x(I,:);
                end
            end

            optimValues.constrviolation = max(cViol, 0);
            optimValues.funccount = iter;
            optimValues.fval = fval;
            optimValues.iteration = iter;
            optimValues.stepsize = 0;
            optimValues.firstorderopt = 0;

            stop = outputFnc(xx, optimValues, state);
            stop = logical(stop);
        end

        function stop = getOutputFunction(x, optimValues, state, hOptimStatusLabel, hFinalStateOptimLabel, hDispAxes, hCancelButton, ...
                                                               objFcn, lb, ub, celBodyData, recorder, propNames, writeOutput, varLabels, lbUsAll, ubUsAll, optimStartTic)
            switch state
                case 'iter'
                    stop = get(hCancelButton,'Value');

                    recorder.iterNums(end+1) = optimValues.iteration;
                    recorder.xVals(end+1) = {x};
                    recorder.fVals(end+1) = optimValues.fval;
                    recorder.maxCVal(end+1) = optimValues.constrviolation;
                case {'init','interrupt','done'}
                    stop = get(hCancelButton,'Value');
            end

            if(stop == true)
                return;
            end

            [~, stateLog] = objFcn(x);

%             finalStateLogEntry = stateLog.getFinalStateLogEntry();
%             finalStateLogEntryMA = finalStateLogEntry.getMAFormattedStateLogMatrix(true);

            stateLogMA = stateLog.getMAFormattedStateLogMatrix(true);

            if(strcmpi(state,'init') || strcmpi(state,'iter'))
                NomadOptimizer.writeOptimStatus(hOptimStatusLabel, optimValues, state, writeOutput, optimStartTic);
                ma_UpdateStateReadout(hFinalStateOptimLabel, 'final', propNames, stateLogMA, celBodyData);
                NomadOptimizer.generatePlots(x, optimValues, state, hDispAxes, lb, ub, varLabels, lbUsAll, ubUsAll);
                drawnow;
            end
        end

        function writeOptimStatus(hOptimStatusLabel, optimValues, state, writeOutput, timer)
            elapTime = toc(timer);

            outStr = {};
            outStr{end+1} = ['State                = ', state];
            outStr{end+1} = '                        ';
            outStr{end+1} = ['Iterations           = ', num2str(optimValues.iteration)];
            outStr{end+1} = ['Function Evals       = ', num2str(optimValues.funccount)];
            outStr{end+1} = ['Objective Value      = ', num2str(optimValues.fval)];
            outStr{end+1} = ['Constraint Violation = ', num2str(optimValues.constrviolation)];
            outStr{end+1} = '                       ';
            outStr{end+1} = ['Elapsed Time         = ', num2str(elapTime), ' sec'];

            set(hOptimStatusLabel, 'String', outStr);

            switch state
                case 'iter'
                    formatstr = ' %- 12.1i %- 12.0i %- 12.6g %- 12.3g';

                    iter = optimValues.iteration;
                    fcnt = optimValues.funccount;
                    val  = optimValues.fval;
                    feas = optimValues.constrviolation;

                    hRow = sprintf(formatstr,iter,fcnt,val,feas);
                    writeOutput(hRow,'append');
                case 'init'
                    hdrStr = sprintf('%- 13s%- 13s%- 13s%- 13s', 'Iteration','Fcn-Count','f(x)-Value', 'Feasibility');
                    writeOutput(hdrStr,'append');
            end
        end

        function generatePlots(x, optimValues, state, hDispAxes, lb, ub, varLabels, lbUsAll, ubUsAll)
            persistent fValPlotIsLog tLayout hPlot1 hPlot2 hPlot3

            if(isempty(fValPlotIsLog))
                fValPlotIsLog = true;
            end

            switch state
                case 'init'
                    if(isvalid(hDispAxes))
%                         set(hDispAxes,'Visible','on');
%                         subplot(hDispAxes);
%                         axes(hDispAxes);
                        tLayout = tiledlayout(hDispAxes, 3,1);
                    end
                    fValPlotIsLog = true;
            end

            hPlot1 = nexttile(tLayout, 1);
            if(strcmpi(state,'init'))

                hPlot1.XTickLabel= [];
                hPlot1.YTickLabel= [];
                hPlot1.ZTickLabel= [];
%                 axes(hPlot1);
            else
%                 axes(hPlot1);
            end
            optimplotxKsptot(x, optimValues, state, lb, ub, varLabels, lbUsAll, ubUsAll);

            hPlot2 = nexttile(tLayout, 2);
            if(strcmpi(state,'init'))

                hPlot2.XTickLabel= [];
                hPlot2.YTickLabel= [];
                hPlot2.ZTickLabel= [];
                h = hPlot2;
            else
                h = hPlot2;
%                 axes(hPlot2);
            end
            if(optimValues.fval<=0)
                fValPlotIsLog = false;
                set(h,'yscale','linear');
            end
            optimplotfvalKsptot(x, optimValues, state);
            if(fValPlotIsLog)
                set(h,'yscale','log');
            else
                set(h,'yscale','linear');
            end
            grid on;
            grid minor;

            hPlot3 = nexttile(tLayout, 3);
            if(strcmpi(state,'init'))
                hPlot3.XTickLabel= [];
                hPlot3.YTickLabel= [];
                hPlot3.ZTickLabel= [];
                h = hPlot3;
            else
                h = hPlot3;
%                 axes(hPlot3);
            end
            optimplotconstrviolationKsptot(x, optimValues, state);

            if(not(isempty(h.Children)))
                hLine = h.Children(1);
                if(isa(hLine,'matlab.graphics.chart.primitive.Line'))
                    yDataLine = hLine.YData;
                    if(abs(max(yDataLine) / min(yDataLine)) >= 10 && all(yDataLine > 0))
                        set(h,'yscale','log');
                    else
                        set(h,'yscale','linear');
                    end
                else
                    set(h,'yscale','linear');
                end
            end

            grid on;
            grid minor;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% ----------------------------------------------------------------------------
% function hfssAddMaterial(fid, Name, Er, bSigma, tanDelta)
%
% Description :
% -------------
% Creates VB Script necessary to add a new material to the HFSS Materials
% Manager.
%
% Parameters :
% ------------
% fid      - file identifier of the HFSS script file.
% Name     - name of the material to be added.
% Er       - dielectric constant of the material to be added.
% bSigma   - bulk conductivity of the material to be added (Siemens)
% tanDelta - loss tangent of the material.
%
% Note :
% ------
% If a material with the given name already exists in HFSS, this code will
% have no effect and the script will continue running. A warning message will
% appear in the HFSS message log.
%
% Example :
% ---------
% fid = fopen('Dipole.vbs', 'wt');
% ...
% hfssAddMaterial(fid, 'CoaxDielectric', 2.07, 0, 0);
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% This file is part of HFSS-MATLAB-API.
%
% HFSS-MATLAB-API 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.
%
% HFSS-MATLAB-API 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
% Foobar; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place, Suite 330, Boston, MA  02111-1307  USA
%
% Copyright 2004, Vijay Ramasami (rvc@ku.edu)
% ----------------------------------------------------------------------------
function hfssAddMaterial(fid, Name, Er, bSigma, tanDelta)

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'oProject.AddMaterial _\n');

% Name.
fprintf(fid, 'Array("NAME:%s", _\n', Name);

% Dielectric Properties.
fprintf(fid, '"permittivity:=", "%f", _\n', Er);
fprintf(fid, '"conductivity:=", "%f", _\n', bSigma);
fprintf(fid, '"dielectric_loss_tangent:=", "%f")\n', tanDelta);

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
clc;
clear all;
warning off;
chos=0;
possibility=11;

while chos~=possibility,
    chos=menu('Digital 3-level image watermarking','select the cover image','select the watermark image','show 3-level coverimage','show 3-level watermarkimage','show watermarked image','show extracted image','Calculate MSE for embedding','Calculate PSNR for embedding','Calculate MSE for extraction','Calculate PSNR for extraction','exit');
    if chos==1
        [fname pname]=uigetfile('*.jpg','select the Cover Image');
        %eval('imageinput=imread(fname)');
        imageinput=imread(fname);

A=rgb2gray(imageinput);
P1=im2double(A);
P=imresize(P1,[2048 2048]);
%imshow(P);
%figure(1);
%title('original image');

[F1,F2]= wfilters('haar', 'd');
[LL,LH,HL,HH] = dwt2(P,'haar','d');
[LL1,LH1,HL1,HH1] = dwt2(LL,'haar','d');
[LL2,LH2,HL2,HH2] = dwt2(LL1,'haar','d');
%figure(2)
%imshow(LL2,'DisplayRange',[]), title('3 level dwt of cover image');
    end
    if chos==2
    [fname pname]=uigetfile('*.jpg','select the Watermark');
    %eval('imageinput=imread(fname)');
    %imageinput=imread(fname);
    imw2=imread(fname);
    imw=rgb2gray(imw2);
watermark=im2double(imw);
watermark=imresize(watermark,[2048 2048]);
%figure(3)
%imshow(uint8(watermark));title('watermark image')
[WF1,WF2]= wfilters('haar', 'd');
[L_L,L_H,H_L,H_H] = dwt2(watermark,'haar','d');
[L_L1,L_H1,H_L1,H_H1] = dwt2(L_L,'haar','d');
[L_L2,L_H2,H_L2,H_H2] = dwt2(L_L1,'haar','d');
%figure(4)
%imshow(L_L2,'DisplayRange',[]), title('3 level dwt of watermark image');
    end
    if chos==3
        imshow(LL2,'DisplayRange',[]), title('3 level dwt of cover image')
    end
    if chos==4
        imshow(L_L2,'DisplayRange',[]), title('3 level dwt of watermark image')
    end
    if chos==5
        Watermarkedimage=LL2+0.0001*L_L2;

%computing level-1 idwt2
Watermarkedimage_level1= idwt2(Watermarkedimage,LH2,HL2,HH2,'haar');
%figure(5)
%imshow(Watermarkedimage_level1,'DisplayRange',[]), title('Watermarkedimage level1');

%computing level-2 idwt2
Watermarkedimage_level2=idwt2(Watermarkedimage_level1,LH1,HL1,HH1,'haar');
%figure(6)
%imshow(Watermarkedimage_level2,'DisplayRange',[]), title('Watermarkedimage level2');

%computing level-3 idwt2
Watermarkedimage_final=idwt2(Watermarkedimage_level2,LH,HL,HH,'haar');
%figure(7)
imshow(Watermarkedimage_final,'DisplayRange',[]), title('Watermarkedimage final')
    end
    if chos==6
        [F11,F22]= wfilters('haar', 'd');
[a b c d]=dwt2(Watermarkedimage_final,'haar','d');
[aa bb cc dd]=dwt2(a,'haar','d');
[aaa bbb ccc ddd]=dwt2(aa,'haar','d');

recovered_image=aaa-LL2;
%figure(8)
imshow(recovered_image,[]);
%title('extracted watermark')
    end
   if chos==7

        pic1= P;
    pic2= Watermarkedimage_final;
       mse=MSE(pic1,pic2)
   end
       if chos==8
            pic1= P;
    pic2= Watermarkedimage_final;
psnr=PSNR(pic1,pic2)
       end
       if chos==9
           clear pic1;
           clear pic2;
           pic1=L_L2;
           pic2=recovered_image;
           mse_extraction=MSE(pic1,pic2)
       end
           if chos==10
           psnr_extraction=PSNR(pic1,pic2)
           end
end


text:
function a = toeplitz_5diag ( n, d1, d2, d3, d4, d5 )

%*****************************************************************************80
%
%% TOEPLITZ_5DIAG returns a pentadiagonal Toeplitz matrix.
%
%  Formula:
%
%    if ( I - J == 2 )
%      A(I,J) = D1
%    elseif ( I - J == 1 )
%      A(I,J) = D2
%    elseif ( I - J == 0 )
%      A(I,J) = D3
%    elseif ( I - J == -1 )
%      A(I,J) = D4
%    elseif ( I - J == -2 )
%      A(I,J) = D5
%    else
%      A(I,J) = 0.0
%
%  Example:
%
%    N = 5, D1 = 1, D2 = -10, D3 = 0, D4 = 10, D5 = 1
%
%      0  10   1   .   .
%    -10   0  10   1   .
%      1 -10   0  10   1
%      .   1 -10   0  10
%      .   .   1 -10   0
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is banded, with bandwidth 5.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    The special data D1 = 1, D2 = -10, D3 = 0, D4 = 10, D5 = 1 corresponds
%    to a matrix of Rutishauser.
%
%    The matrix has eigenvalues lying approximately on the complex line
%    segment 2 * cos ( 2 * t ) + 20 * I * sin ( t ).
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    RM Beam, RF Warming,
%    The asymptotic spectra of banded Toeplitz and quasi-Toeplitz matrices,
%    SIAM Journal on Scientific and Statistical Computing,
%    Volume 14, Number 4, 1993, pages 971-1006.
%
%    Heinz Rutishauser,
%    On test matrices,
%    Programmation en Mathematiques Numeriques,
%    Centre National de la Recherche Scientifique,
%    1966, pages 349-365.
%
%  Parameters:
%
%    Input, integer N, the order of A.  N should be at least 3.
%
%    Input, D1, D2, D3, D4, D5, values that define the nonzero diagonals
%    of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

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

      if ( i - j == 2 )
        a(i,j) = d1;
      elseif ( i - j == 1 )
        a(i,j) = d2;
      elseif ( i - j == 0 )
        a(i,j) = d3;
      elseif ( i - j == -1 )
        a(i,j) = d4;
      elseif ( i - j == -2 )
        a(i,j) = d5;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
classdef KSPTOT_BodyInfo < matlab.mixin.SetGet
    %KSPTOT_BodyInfo Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %Orbit
        epoch double
        sma double
        ecc double
        inc double
        raan double
        arg double
        mean double

        %Physical
        gm double
        radius double

        %Atmo
        atmohgt double
        atmopresscurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        atmotempcurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        atmotempsunmultcurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        lattempbiascurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        lattempsunmultcurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        axialtempsunbiascurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        axialtempsunmultcurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        ecctempbiascurve = griddedInterpolant([0 1]',[0 0]','spline', 'nearest');
        atmomolarmass double

        %Rotation
        rotperiod double
        rotini double

        %Display
        bodycolor char

        %Others
        canbecentral double
        canbearrivedepart double
        parent
        parentid double
        name char
        id double

        %body surface textures
        surftexturefile (1,:) char = '';
        surftexturezrotoffset(1,1) double = 0; %degrees

        %orientation of inertial frame (spin axis handling)
        bodyZAxis(3,1) double = [0;0;1];
        bodyXAxis(3,1) double = [1;0;0];
        bodyRotMatFromGlobalInertialToBodyInertial = [];

        %Propagation type
        proptype(1,:) char = 'analytic_two_body';
        propTypeEnum(1,1) BodyPropagationTypeEnum = BodyPropagationTypeEnum.TwoBody
        numIntStateCache CelestialBodySunRelStateDataCache
    end

    properties
        celBodyData

        parentBodyInfo KSPTOT_BodyInfo
        parentBodyInfoNeedsUpdate logical = true;

        childrenBodyInfo KSPTOT_BodyInfo
        childrenBodyNames cell
        childrenBodyInfoNeedsUpdate logical = true;
    end

    properties(Transient)
        doNotUseAtmoTempSunMultCurve(1,1) logical = false;
        doNotUseLatTempSunMultCurve(1,1) logical = false;
        doNotUseAxialTempSunMultCurve(1,1) logical = false;
        doNotUseAxialTempSunBiasCurveGI(1,1) logical = false;
        doNotUseEccTempBiasCurveGI(1,1) logical = false;
        doNotUseAtmoPressCurveGI(1,1) logical = false;

        orbitElemsChainCache cell = {};
        fixedFrameFromInertialFrameCache cell = {};

        surfTextureCache uint8 = [];

        lastComputedTime = NaN;
        lastComputedRVect = NaN(3,1);
        lastComputedVVect = NaN(3,1);

        lastComputedRVectVVect = [];

        propTypeIsTwoBody(1,1) logical = true;
        propTypeIsNumerical(1,1) logical = false;
    end

    properties(Access=private)
        bodyInertialFrameCache BodyCenteredInertialFrame = BodyCenteredInertialFrame.empty(1,0);
        bodyFixedFrameCache BodyFixedFrame = BodyFixedFrame.empty(1,0);

        atmoTempCache AtmoTempDataCache = AtmoTempDataCache.empty(1,0);
        sunDotNormalCache SunDotNormalDataCache = SunDotNormalDataCache.empty(1,0);
        parentGmuCache(1,1) double = NaN;

        soiRadiusCache(1,1) double = NaN;
    end

    methods
        function obj = KSPTOT_BodyInfo()
            obj.sunDotNormalCache = SunDotNormalDataCache(obj);
            obj.atmoTempCache = AtmoTempDataCache(obj);
            obj.numIntStateCache = CelestialBodySunRelStateDataCache(obj);

            obj.bodyInertialFrameCache = BodyCenteredInertialFrame(obj, obj.celBodyData);
            obj.bodyFixedFrameCache = BodyFixedFrame(obj, obj.celBodyData);

            if(obj.propTypeEnum == BodyPropagationTypeEnum.TwoBody)
                obj.propTypeIsTwoBody = true;
                obj.propTypeIsNumerical = false;
            elseif(obj.propTypeEnum == BodyPropagationTypeEnum.Numerical)
                obj.propTypeIsTwoBody = false;
                obj.propTypeIsNumerical = true;
            end
        end

        function set.atmotempsunmultcurve(obj,newValue)
            obj.atmotempsunmultcurve = newValue;
            obj.doNotUseAtmoTempSunMultCurve = isa(obj.atmotempsunmultcurve,'griddedInterpolant') && all(obj.atmotempsunmultcurve.Values == 0); %#ok<MCSUP>
        end

        function set.lattempsunmultcurve(obj,newValue)
            obj.lattempsunmultcurve = newValue;
            obj.doNotUseLatTempSunMultCurve = isa(obj.lattempsunmultcurve,'griddedInterpolant') && all(obj.lattempsunmultcurve.Values == 0); %#ok<MCSUP>
        end

        function set.axialtempsunmultcurve(obj,newValue)
            obj.axialtempsunmultcurve = newValue;
            obj.doNotUseAxialTempSunMultCurve = isa(obj.axialtempsunmultcurve,'griddedInterpolant') && all(obj.axialtempsunmultcurve.Values == 0); %#ok<MCSUP>
        end

        function set.axialtempsunbiascurve(obj,newValue)
            obj.axialtempsunbiascurve = newValue;
            obj.doNotUseAxialTempSunBiasCurveGI = isa(obj.axialtempsunbiascurve,'griddedInterpolant') && all(obj.axialtempsunbiascurve.Values == 0); %#ok<MCSUP>
        end

        function set.ecctempbiascurve(obj,newValue)
            obj.ecctempbiascurve = newValue;
            obj.doNotUseEccTempBiasCurveGI = isa(obj.ecctempbiascurve,'griddedInterpolant') && all(obj.ecctempbiascurve.Values == 0); %#ok<MCSUP>
        end

        function set.atmopresscurve(obj,newValue)
            obj.atmopresscurve = newValue;

            if(not(isempty(obj.atmohgt))) %#ok<MCSUP>
                obj.doNotUseAtmoPressCurveGI = logical(obj.atmohgt > 0 && isempty(obj.atmopresscurve)); %#ok<MCSUP>
            end
        end

        function set.propTypeEnum(obj, newValue)
            obj.propTypeEnum = newValue;

            if(obj.propTypeEnum == BodyPropagationTypeEnum.TwoBody)
                obj.propTypeIsTwoBody = true; %#ok<MCSUP>
                obj.propTypeIsNumerical = false; %#ok<MCSUP>
            elseif(obj.propTypeEnum == BodyPropagationTypeEnum.Numerical)
                obj.propTypeIsTwoBody = false; %#ok<MCSUP>
                obj.propTypeIsNumerical = true; %#ok<MCSUP>
            end
        end

        function parentBodyInfo = getParBodyInfo(obj, ~)
            if(obj.parentBodyInfoNeedsUpdate || (obj.parentid > 0 && isempty(obj.parentBodyInfo)))
                pBodyInfo = getParentBodyInfo(obj, obj.celBodyData);

                if(not(isempty(pBodyInfo)))
                    obj.parentBodyInfo = pBodyInfo;
                end

                obj.soiRadiusCache = NaN;
                obj.parentBodyInfoNeedsUpdate = false;
            end

            parentBodyInfo = obj.parentBodyInfo;
        end

        function [childrenBodyInfo, cbNames] = getChildrenBodyInfo(obj, celBodyData)
            if(obj.childrenBodyInfoNeedsUpdate)
                [cBodyInfo,cbNames] = getChildrenOfParentInfo(celBodyData, obj.name);

                if(not(isempty(cBodyInfo)))
                    for(i=1:length(cBodyInfo)) %#ok<*NO4LP>
                        obj.childrenBodyInfo(i) = cBodyInfo{i};
                    end
                    obj.childrenBodyNames = cbNames;
                end

                obj.childrenBodyInfoNeedsUpdate = false;
            end

            childrenBodyInfo = obj.childrenBodyInfo;
            cbNames = obj.childrenBodyNames;
        end

        function soiRadius = getCachedSoIRadius(obj)
            if(isnan(obj.soiRadiusCache))
                obj.soiRadiusCache = getSOIRadius(obj, obj.getParBodyInfo);
            end

            soiRadius = obj.soiRadiusCache;
        end

        function rotMat = get.bodyRotMatFromGlobalInertialToBodyInertial(obj)
            if(isempty(obj.bodyRotMatFromGlobalInertialToBodyInertial))
                obj.bodyZAxis = normVector(obj.bodyZAxis);
                obj.bodyXAxis = normVector(obj.bodyXAxis);

                rotY = normVector(crossARH(obj.bodyZAxis, obj.bodyXAxis));
                rotX = normVector(crossARH(rotY, obj.bodyZAxis));
                rotMat = [rotX, rotY, obj.bodyZAxis];

                obj.bodyRotMatFromGlobalInertialToBodyInertial = rotMat;
            else
                rotMat = obj.bodyRotMatFromGlobalInertialToBodyInertial;
            end
        end

        function frame = getBodyCenteredInertialFrame(obj)
%             if(isempty(obj.bodyInertialFrameCache))
%                 obj.bodyInertialFrameCache = BodyCenteredInertialFrame(obj, obj.celBodyData);
%             end

            frame = obj.bodyInertialFrameCache;
        end

        function frame = getBodyFixedFrame(obj)
%             if(isempty(obj.bodyFixedFrameCache))
%                 obj.bodyFixedFrameCache = BodyFixedFrame(obj, obj.celBodyData);
%             end

            frame = obj.bodyFixedFrameCache;
        end

        function states = getElementSetsForTimes(obj, times)
            parBodyInfo = obj.getParBodyInfo();
            if(not(isempty(parBodyInfo)))
                frame = parBodyInfo.getBodyCenteredInertialFrame();
                gmu = obj.getParentGmuFromCache();

                [rVects, vVects] = getStateAtTime(obj, times, gmu);

                states = CartesianElementSet(times, rVects, vVects, frame);
            else
                frame = obj.getBodyCenteredInertialFrame();

                for(i=1:length(times))
                    states(i) = CartesianElementSet(times(i), [0;0;0], [0;0;0], frame); %#ok<AGROW>
                end
            end
        end

        function [rVect, vVect] = getCachedPositionVelWrtSun(obj, time)
            if(numel(time) == 1)
                if(time == obj.lastComputedTime)
                    rVect = obj.lastComputedRVect;
                    vVect = obj.lastComputedVVect;
                else
                    [rVect, vVect] = getPositOfBodyWRTSun(time, obj, obj.celBodyData);

                    obj.lastComputedRVect = rVect;
                    obj.lastComputedVVect = vVect;
                end
            else
                [rVect, vVect] = getPositOfBodyWRTSun(time, obj, obj.celBodyData);
            end
        end

        function bColorRGB = getBodyRGB(obj)
            bColor = obj.bodycolor;

%             cmap = colormap(bColor);
            try
                cmap = feval(bColor);
            catch
                cmap = feval('gray');
                warning('Could not find a colormap for body color "%s".  Defaulting to grey.', bColor);
            end

            midRow = round(size(cmap,1)/2);
            bColorRGB = cmap(midRow,:);
        end

        function temperature = getBodyAtmoTemperature(obj, time, lat, long, alt)
            temperature = getTemperatureAtAltitude(obj, alt, lat, time, long);
        end

        function pressure = getBodyAtmoPressure(obj, altitude)
            pressure = getPressureAtAltitude(obj, altitude);
        end

        function sunDotNormal = getCachedSunDotNormal(obj, time, long)
            sunDotNormal = obj.sunDotNormalCache.getCachedBodyStateAtTimeAndLong(time, long);
        end

        function parentGM = getParentGmuFromCache(obj)
            if(isnan(obj.parentGmuCache) || obj.parentBodyInfoNeedsUpdate == true)
                obj.parentGmuCache = getParentGM(obj, obj.celBodyData);
                obj.parentBodyInfoNeedsUpdate = false;
            end

            parentGM = obj.parentGmuCache;
        end

        function chain = getOrbitElemsChain(obj)
            if(isempty(obj.orbitElemsChainCache))
                smas = [];
                eccs = [];
                incs = [];
                raans = [];
                args = [];
                means = [];
                epochs = [];
                parentGMs = [];

                loop = true;
                bodyInfo = obj;
                while(loop)
                    smas(end+1) = bodyInfo.sma; %#ok<AGROW>
                    eccs(end+1) = bodyInfo.ecc; %#ok<AGROW>
                    incs(end+1) = bodyInfo.inc; %#ok<AGROW>
                    raans(end+1) = bodyInfo.raan; %#ok<AGROW>
                    args(end+1) = bodyInfo.arg; %#ok<AGROW>
                    means(end+1) = bodyInfo.mean; %#ok<AGROW>
                    epochs(end+1) = bodyInfo.epoch; %#ok<AGROW>

                    try
                        thisParentBodyInfo = bodyInfo.getParBodyInfo(obj.celBodyData);
                    catch
                        thisParentBodyInfo = getParentBodyInfo(bodyInfo, obj.celBodyData);
                    end

                    if(isempty(thisParentBodyInfo))
                        parentGMs(end+1) = 0;

                        break;
                    else
                        parentGMs(end+1) = bodyInfo.getParentGmuFromCache(); %#ok<AGROW>

                        bodyInfo = thisParentBodyInfo;
                    end
                end

                obj.orbitElemsChainCache = {smas, eccs, incs, raans, args, means, epochs, parentGMs};
            end

            chain = obj.orbitElemsChainCache;
        end

        function inputs = getFixedFrameFromInertialFrameInputsCache(obj)
            if(isempty(obj.fixedFrameFromInertialFrameCache))
                obj.fixedFrameFromInertialFrameCache = {obj.rotperiod, obj.rotini};
            end

            inputs = obj.fixedFrameFromInertialFrameCache;
        end

        function [surfTexture] = getSurfaceTexture(obj)
            if(isempty(obj.surfTextureCache))
                if(not(isempty(obj.surftexturefile)))
                    try
                        [I,~] = imread(obj.surftexturefile);
                        I = flip(I, 1);
                    catch ME
                        I = NaN;
                        warning('Could not load surface texture for "%s".  Message: \n\n%s', obj.name, ME.message);
                    end

                    obj.surfTextureCache = I;
                else
                    obj.surfTextureCache = NaN;
                end
            end

            surfTexture = obj.surfTextureCache;
        end

        function tf = eq(A,B)
            tf = [A.id] == [B.id];

            if(isempty(tf))
                tf = false;
            end
        end

        function setPropTypeEnum(obj)
            switch obj.proptype
                case 'analytic_two_body'
                    obj.propTypeEnum = BodyPropagationTypeEnum.TwoBody;
                case 'numerical_integration'
                    obj.propTypeEnum = BodyPropagationTypeEnum.Numerical;
                otherwise
                    error('Unknown Body Prop Type: %s', obj.propType);
            end
        end

        function setStateCacheData(obj, times, stateVects, frame)
            obj.numIntStateCache.setStateData(times, stateVects, frame);
        end
    end

    methods
        function setBaseBodySurfaceTexture(obj)
            if(isempty(obj.surftexturefile) && strcmpi(obj.name,'Sun') && obj.id == 0)
                obj.surftexturefile = 'images/body_textures/surface/sunSurface.jpg';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Kerbin') && obj.id == 1)
                obj.surftexturefile = 'images/body_textures/surface/kerbinSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Mun') && obj.id == 2)
                obj.surftexturefile = 'images/body_textures/surface/munSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Minmus') && obj.id == 3)
                obj.surftexturefile = 'images/body_textures/surface/minmusSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Moho') && obj.id == 4)
                obj.surftexturefile = 'images/body_textures/surface/mohoSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Eve') && obj.id == 5)
                obj.surftexturefile = 'images/body_textures/surface/eveSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Gilly') && obj.id == 13)
                obj.surftexturefile = 'images/body_textures/surface/gillySurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Duna') && obj.id == 6)
                obj.surftexturefile = 'images/body_textures/surface/dunaSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Ike') && obj.id == 7)
                obj.surftexturefile = 'images/body_textures/surface/ikeSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Dres') && obj.id == 15)
                obj.surftexturefile = 'images/body_textures/surface/dresSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Jool') && obj.id == 8)
                obj.surftexturefile = 'images/body_textures/surface/joolSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Laythe') && obj.id == 9)
                obj.surftexturefile = 'images/body_textures/surface/laytheSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Vall') && obj.id == 10)
                obj.surftexturefile = 'images/body_textures/surface/vallSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Tylo') && obj.id == 12)
                obj.surftexturefile = 'images/body_textures/surface/tyloSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Bop') && obj.id == 11)
                obj.surftexturefile = 'images/body_textures/surface/bopSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Pol') && obj.id == 14)
                obj.surftexturefile = 'images/body_textures/surface/polSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Eeloo') && obj.id == 16)
                obj.surftexturefile = 'images/body_textures/surface/eelooSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Mercury') && obj.id == 199)
                obj.surftexturefile = 'images/body_textures/surface/mercurySurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Venus') && obj.id == 299)
                obj.surftexturefile = 'images/body_textures/surface/venusSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Earth') && obj.id == 399)
                obj.surftexturefile = 'images/body_textures/surface/earthSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Moon') && obj.id == 301)
                obj.surftexturefile = 'images/body_textures/surface/moonSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Mars') && obj.id == 499)
                obj.surftexturefile = 'images/body_textures/surface/marsSurface.png';

            elseif(isempty(obj.surftexturefile) && strcmpi(obj.name,'Pluto') && obj.id == 999)
                obj.surftexturefile = 'images/body_textures/surface/plutoSurface.png';
            end
        end
    end

	methods(Static)
        function obj = loadobj(obj)
            if(isempty(obj.epoch))
                return; %this should only happen if something is bugged out
            end

            obj.doNotUseAtmoTempSunMultCurve = isa(obj.atmotempsunmultcurve,'griddedInterpolant') && all(obj.atmotempsunmultcurve.Values == 0);
            obj.doNotUseLatTempSunMultCurve = isa(obj.lattempsunmultcurve,'griddedInterpolant') && all(obj.lattempsunmultcurve.Values == 0);
            obj.doNotUseAxialTempSunMultCurve = isa(obj.axialtempsunmultcurve,'griddedInterpolant') && all(obj.axialtempsunmultcurve.Values == 0);
            obj.doNotUseAxialTempSunBiasCurveGI = isa(obj.axialtempsunbiascurve,'griddedInterpolant') && all(obj.axialtempsunbiascurve.Values == 0);
            obj.doNotUseEccTempBiasCurveGI = isa(obj.ecctempbiascurve,'griddedInterpolant') && all(obj.ecctempbiascurve.Values == 0);
            obj.doNotUseAtmoPressCurveGI = (obj.atmohgt > 0 && isempty(obj.atmopresscurve));

            obj.parentBodyInfoNeedsUpdate = true;
            obj.getParBodyInfo(obj.celBodyData);

            obj.childrenBodyInfoNeedsUpdate = true;
            obj.getChildrenBodyInfo(obj.celBodyData);

            obj.setBaseBodySurfaceTexture();

            if(isempty(obj.sunDotNormalCache))
                obj.sunDotNormalCache = SunDotNormalDataCache(obj);
            end

            if(isempty(obj.atmoTempCache))
                obj.atmoTempCache = AtmoTempDataCache(obj);
            end

            if(isempty(obj.numIntStateCache))
                obj.numIntStateCache = CelestialBodySunRelStateDataCache(obj);
            end

            obj.setPropTypeEnum();

            if(isempty(obj.bodyInertialFrameCache))
                obj.bodyInertialFrameCache = BodyCenteredInertialFrame(obj, obj.celBodyData);
            end

            if(isempty(obj.bodyFixedFrameCache))
                obj.bodyFixedFrameCache = BodyFixedFrame(obj, obj.celBodyData);
            end

            if(obj.propTypeEnum == BodyPropagationTypeEnum.TwoBody)
                obj.propTypeIsTwoBody = true;
                obj.propTypeIsNumerical = false;
            elseif(obj.propTypeEnum == BodyPropagationTypeEnum.Numerical)
                obj.propTypeIsTwoBody = false;
                obj.propTypeIsNumerical = true;
            end
        end

        function bodyObj = getObjFromBodyInfoStruct(bodyInfo)
            bodyObj = KSPTOT_BodyInfo();

            bodyObj.epoch = bodyInfo.epoch;
            bodyObj.sma = bodyInfo.sma;
            bodyObj.ecc = bodyInfo.ecc;
            bodyObj.inc = bodyInfo.inc;
            bodyObj.raan = bodyInfo.raan;
            bodyObj.arg = bodyInfo.arg;
            bodyObj.mean = bodyInfo.mean;
            bodyObj.gm = bodyInfo.gm;
            bodyObj.radius = bodyInfo.radius;
            bodyObj.atmohgt = bodyInfo.atmohgt;
            bodyObj.atmopresscurve = bodyInfo.atmopresscurve;
            bodyObj.atmotempcurve = bodyInfo.atmotempcurve;
            bodyObj.atmotempsunmultcurve = bodyInfo.atmotempsunmultcurve;
            bodyObj.lattempbiascurve = bodyInfo.lattempbiascurve;
            bodyObj.lattempsunmultcurve = bodyInfo.lattempsunmultcurve;
            bodyObj.atmomolarmass = bodyInfo.atmomolarmass;
            bodyObj.rotperiod = bodyInfo.rotperiod;
            bodyObj.rotini = bodyInfo.rotini;
            bodyObj.bodycolor = bodyInfo.bodycolor;
            bodyObj.canbecentral = bodyInfo.canbecentral;
            bodyObj.canbearrivedepart = bodyInfo.canbearrivedepart;
            bodyObj.parent = bodyInfo.parent;
            bodyObj.parentid = bodyInfo.parentid;
            bodyObj.name = bodyInfo.name;
            bodyObj.id = bodyInfo.id;
            bodyObj.celBodyData = bodyInfo.celBodyData;

            bodyObj.bodyZAxis = bodyInfo.bodyZAxis;
            bodyObj.bodyXAxis = bodyInfo.bodyXAxis;
            bodyObj.bodyRotMatFromGlobalInertialToBodyInertial = bodyInfo.bodyRotMatFromGlobalInertialToBodyInertial;
        end
    end
end


text:
% Test file for bndfun/sum.m

function pass = test_sum(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Set a domain
dom = [-2 7];

% Create seed for random number generator
seedRNG(6178);

% Generate a few random points to use as test values.
x = diff(dom) * rand(1000, 1) + dom(1);

%%
% Spot-check integrals for a couple of functions.
f = bndfun(@(x) exp(x) - 1, struct('domain', dom), pref);
pass(1) = (abs(sum(f) -  1.087497823145222e3) < ...
    10*get(f, 'vscale')*eps);

f = bndfun(@(x) 1./(1 + x.^2), struct('domain', dom), pref);
pass(2) = (abs(sum(f) - (atan(-dom(1))+atan(dom(2)))) < ...
    10*get(f, 'vscale')*eps);

f = bndfun(@(x) cos(1e4*x), struct('domain', dom), pref);
exact = (sin(1e4*dom(2))-sin(1e4*dom(1)))/1e4;
pass(3) = (abs(sum(f) - exact) < 100*get(f, 'vscale')*eps);

z = exp(2*pi*1i/6);
f = bndfun(@(t) sinh(t*z), struct('domain', dom), pref);
exact = ((cos(sqrt(3) - 1i) - cos((7*sqrt(3))/2 - 7i/2))*(sqrt(3) + 1i)*1i)/2;
pass(4) = (abs(sum(f)-exact) < 10*get(f, 'vscale')*eps);

%%
% Check a few basic properties.
a = 2;
b = -1i;
f = bndfun(@(x) x.*sin(x.^2) - 1, struct('domain', dom), pref);
df = diff(f);
g = bndfun(@(x) exp(-(x/10).^2), struct('domain', dom), pref);
dg = diff(g);
fg = f.*g;
gdf = g.*df;
fdg = f.*dg;

tol_f = 10*max(get(f, 'vscale')*eps);
tol_g = 10*max(get(g, 'vscale')*eps);
tol_fg = 10*max(get(fg, 'vscale')*eps);
tol_df = 10*max(get(df, 'vscale')*eps);
tol_dg = 10*max(get(dg, 'vscale')*eps);
tol_gdf = 10*max(get(gdf, 'vscale')*eps);
tol_fdg = 10*max(get(fdg, 'vscale')*eps);

% Linearity.
pass(5) = (abs(sum(a*f + b*g) - (a*sum(f) + b*sum(g))) < max(tol_f, tol_g));

% Integration-by-parts.
pass(6) = (abs(sum(fdg) - (feval(fg, dom(2)) - ...
    feval(fg, dom(1)) - sum(gdf))) < max([tol_fdg ; tol_fg ; tol_gdf]));

% Fundamental Theorem of Calculus.
pass(7) = (abs(sum(df) - (feval(f, dom(2)) - feval(f, dom(1)))) < ...
    max(tol_df, tol_f));
pass(8) = (abs(sum(dg) - (feval(g, dom(2)) - feval(g, dom(1)))) < ...
    max(tol_dg, tol_g));

%%
% Check operation for array-valued bndfun objects.
f = bndfun(@(x) [sin(x) x.^2 exp(1i*x)], struct('domain', dom), pref);
I = sum(f);
I_exact = [(cos(dom(1)) - cos(dom(2))) (dom(2)^3 - dom(1)^3)/3 ...
    1i*(exp(1i*dom(1)) - exp(1i*dom(2)))];
pass(9) = (max(abs(I - I_exact)) < 10*max(get(f, 'vscale')*eps));

% DIM option with array-valued input.
g = sum(f, 2);
h = @(x) sin(x) + x.^2 + exp(1i*x);
pass(10) = (norm(feval(g, x) - h(x), inf) < ...
    10*max(get(g, 'vscale')*eps));

% DIM option with non-array-valued input should leave everything as it was.
h = bndfun(@(x) cos(x), struct('domain', dom), pref);
sumh2 = sum(h, 2);
pass(11) = all((feval(h, x) == feval(sumh2, x)));

%% Test on singular function:

pow = -0.5;
op = @(x) (x - dom(1)).^pow.*sin(x);
pref.blowup = true;
data.domain = dom;
data.exponents = [pow 0];
f = bndfun(op, data, pref);
I = sum(f);
I_exact = -1.92205524578386613;
pass(12) = ( abs(I-I_exact) < 200*eps*abs(I_exact) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function show_r3_point_set(points_x,varargin)
%SHOW_R3_POINT_SET 3D illustration of a point set
%
%Syntax
% show_r3_point_set(POINTS_X,options);
%
%Description
% SHOW_R3_POINT_SET(POINTS_X) uses a 3d plot to illustrate a point set in relation to
% the unit sphere S^2.
%
% The argument POINTS_X must be an array of real numbers of size (3 by N), where N is a
% positive integer, representing N points of R^3.
%
% SHOW_R3_POINT_SET(POINTS_X,options) also recognizes a number of illustration
% options, which are specified as name, value pairs.
% Any number of pairs can be used, in any order.
%
% The following illustration options are used.
%
% SHOW_R3_POINT_SET(POINTS_X,'fontsize',size)
% Font size used in titles (numeric, default 16).
%
% SHOW_R3_POINT_SET(POINTS_X,'title','show')
% SHOW_R3_POINT_SET(POINTS_X,'title','hide')
% Show or hide title (default 'hide').
%
% SHOW_R3_POINT_SET(POINTS_X,'sphere','show')
% SHOW_R3_POINT_SET(POINTS_X,'sphere','hide')
% Show or hide the unit sphere S^2 (default 'hide').
%
% For more detail on illustration options, see HELP ILLUSTRATION_OPTIONS.
%
%Note
% This function is primarily for use with the point set POINTS_X as a subset of the
% unit sphere S^2, but this is not assumed and not checked.
% If you show the unit sphere S^2 and POINTS_X contains points closer than radius 1
% from the origin, the sphere will hide these points.
%
%Examples
% > points_x
% points_x =
%          0    0.0000   -0.0000    0.0000
%          0    1.0000   -1.0000         0
%     1.0000    0.0000    0.0000   -1.0000
%
% > show_r3_point_set(points_x,'sphere','hide')
% > show_r3_point_set(points_x,'sphere','show')
%
%See also
% ILLUSTRATION_OPTIONS, SHOW_S2_PARTITION, PROJECT_POINT_SET

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

gdefault.fontsize = 16;
gdefault.show_title  = false;
gdefault.show_sphere = false;

gopt = illustration_options(gdefault, varargin{:});

N = size(points_x,2);

surf_jet;

if gopt.show_title
    titlestr = sprintf(...
        '\nPoint set containing %d points.',N);
    title(titlestr,'FontWeight','bold','FontUnits','normalized',...
        'FontSize',gopt.fontsize/512);
end

if gopt.show_sphere
    show_s2_sphere;
    hold on
end

[X,Y,Z] = sphere;

r = min(0.05,N^(-1/2)/2);
rX = r*X; rY = r*Y; rZ = r*Z;
for n = 1:N
   surf(points_x(1,n)+rX,points_x(2,n)+rY,points_x(3,n)+rZ,ones(size(rZ)),...
   'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
end
%
axis equal
axis off
grid off
hold off

repo: gspbox


text:
function print_progress(n, N)
%print_progress(n, N)
%
% print work progress of loops (tic, for n=1:N)

%% idea by stefan harmeling

cent= round(100*n/N);
togo= toc/n*(N-n);
fprintf('\r%3d%%, %ds to go    ', cent, floor(togo));
if n==N,
  fprintf('\r%40s\r','');
end

author: bbci


text:
clear();

CONVOLUTION_SHAPE_FULL  = 1;
CONVOLUTION_SHAPE_SAME  = 2;
CONVOLUTION_SHAPE_VALID = 3;

maxThr = 1e-9;

tic();
for numRowsImage = 28:32
    for numColsImage = 28:32

        mI = rand(numRowsImage, numColsImage);

        for numRowsKernel = 3:7
            for numColsKernel = 3:7

                mH = rand(numRowsKernel, numColsKernel);

                for convShape = 1:3

                    switch(convShape)
                        case(CONVOLUTION_SHAPE_FULL)
                            numRowsOut = numRowsImage + numRowsKernel - 1;
                            numColsOut = numColsImage + numColsKernel - 1;

                            convShapeString = 'full';
                        case(CONVOLUTION_SHAPE_SAME)
                            numRowsOut = numRowsImage;
                            numColsOut = numColsImage;

                            convShapeString = 'same';
                        case(CONVOLUTION_SHAPE_VALID)
                            numRowsOut = numRowsImage - numRowsKernel + 1;
                            numColsOut = numColsImage - numColsKernel + 1;

                            convShapeString = 'valid';
                    end

                    mORef   = conv2(mI, mH, convShapeString);
                    % mK      = CreateConvMtx2D(mH, numRowsImage, numColsImage, convShape);
                    mK      = CreateConvMtx2DSparse(mH, numRowsImage, numColsImage, convShape);
                    mO      = reshape(mK * mI(:), numRowsOut, numColsOut);

                    disp([' ']);
                    disp(['Validating solution for the following parameters:']);
                    disp(['Image Size - [', num2str(numRowsImage), ' x ', num2str(numColsImage), ']']);
                    disp(['Kernel Size - [', num2str(numRowsKernel), ' x ', num2str(numColsKernel), ']']);
                    disp(['Convolution Shape - ', convShapeString]);

                    mE = mO - mORef;
                    maxAbsDev = max(abs(mE(:)));
                    if(maxAbsDev >= maxThr)
                        disp([' ']);
                        disp(['Validation Failed']);
                        disp([' ']);
                    end
                    assert(maxAbsDev < maxThr);

                end
            end
        end
    end
end

toc();

author: RoyiAvital

repo: StackExchangeCodes