text:
function [features, labels, labelprobs, weights] = ...
    mcmcGenerateGoodSegments(im, imsegs, vclassifierSP, hclassifierSP,eclassifier, ...
    segclassifier, spdata, adjlist, edata)
% Computes the marginals of the geometry for the given input image
% spdata, adjlist, edata are optional inputs

niter = 2500;

features = zeros(imsegs.nseg*6, 102);
labels = zeros(imsegs.nseg*6, 1);
labelprobs = zeros(imsegs.nseg*6, 7);
weights = zeros(imsegs.nseg*6, 1);

count = 0;

[pvSP, phSP, pE, smap] = mcmcInitialize(spdata, edata, ...
    adjlist, imsegs, vclassifierSP, hclassifierSP, eclassifier, 'labels');

imdata = mcmcComputeImageData(im, imsegs);
imdata.pvSP = pvSP;
imdata.phSP = phSP;

nseg = max(smap);

segProb = computeSegProb(segclassifier, pvSP, phSP, pE, adjlist, smap, 1:nseg);

[features, labels, labelprobs, weights, count] = ...
    addFeatures(features, labels, labelprobs, weights, count, ...
    imsegs, spdata, imdata, smap);

% do mcmc
for iter = 2:niter

    [smap, newseg, origseg, segadjmat] = mcmcGenerateProposals(pE, adjlist, smap);
   % disp(num2str(iter))

    nseg2 = max(smap);

    tmpSegProb = segProb;
    if nseg2 > nseg
        tmpSegProb(newseg) = computeSegProb(segclassifier, pvSP, phSP, pE, adjlist, smap, newseg);
        if origseg~=newseg
            tmpSegProb(origseg) = computeSegProb(segclassifier, pvSP, phSP, pE, adjlist, smap, origseg);
        end
    end

    % get all possible segmentation labelings for newseg
    neighborseg = [newseg find(segadjmat(newseg, :))];
    nn = numel(neighborseg);

    % compute transition probability to each neighbor
    transprob = zeros(nn, 1);

    sind = find(smap==newseg);
    adj1 = find((smap(adjlist(:, 1))==newseg) | (smap(adjlist(:, 2))==newseg));
    for ni = 1:nn

        s = neighborseg(ni);

        if s==newseg
            transprob(ni) = 1;
        else
            eind = adj1(find((smap(adjlist(adj1, 1))==s) | (smap(adjlist(adj1, 2))==s)));
            smap2 = smap;
            smap2(sind) = s;

            sind2 = find(smap2 > newseg);
            smap2(sind2) = smap2(sind2) - 1;
            s2 = s - (s > newseg);

            probni = computeSegProb(segclassifier, pvSP, phSP, pE, adjlist, smap2, s2);
            transprob(ni) = probni / tmpSegProb(newseg) / tmpSegProb(s);
            transprob(ni) = transprob(ni)*prod(1-pE(eind));
        end

    end

    %disp(num2str(transprob'))

    % randomly select neighbor according to transprob
    transprob = cumsum(transprob / sum(transprob));
    s = neighborseg(find(rand(1) < transprob));
    s = s(1);

    sind = find(smap==newseg);
    segProb = tmpSegProb;
    if origseg==newseg && s==newseg  % nothing changed
        % do nothing
    elseif origseg~=newseg && s==newseg  % newseg split from origseg
    else % s~=newseg - newseg merged into s

        smap(sind) = s;
        sind2 = find(smap > newseg);
        smap(sind2) = smap(sind2) - 1;
        s = s - (s > newseg);

        segProb(1:end-1) = segProb([1:newseg-1 newseg+1:end]);
        segProb(s) = computeSegProb(segclassifier, pvSP, phSP, pE, adjlist, smap, s);

    end

    if mod(iter, 500)==0
        disp(num2str(iter))
        [features, labels, labelprobs, weights, count] = ...
            addFeatures(features, labels, labelprobs, weights, count, ...
            imsegs, spdata, imdata, smap);
    end

end

features = features(1:count, :);
labels = labels(1:count);
labelprobs = labelprobs(1:count, :);
weights = weights(1:count);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function segProb = computeSegProb(segclassifier, pvSP, phSP, pE, adjlist, smap, sind)

nseg = numel(sind);
segProb = zeros(nseg, 1);

for i = 1:nseg

    s = sind(i);

    data = mcmcGetSegmentationFeatures(pvSP, phSP, pE, adjlist, smap, s);
    conf = test_boosted_dt_mc(segclassifier, data);
    segProb(i) = 1 / (1+exp(-conf));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [features, labels, labelprobs, weights, count] = ...
    addFeatures(features, labels, labelprobs, weights, count, ...
    imsegs, spdata, imdata, smap)

nseg = max(smap);

for s = 1:nseg

    features(count+s, :) = mcmcGetSegmentFeatures(imsegs, spdata, imdata, smap, s);

    spind = find(smap==s);
    labelcount = zeros(1, 7);
    for sp = spind
        if imsegs.labels(sp)~=0
            labelcount(imsegs.labels(sp)) = ...
                labelcount(imsegs.labels(sp)) + imsegs.npixels(sp);
        end
    end
    labels(count+s) = 0;
    if sum(labelcount)>0
        labelprobs(count+s, :) = labelcount / sum(labelcount);
        [maxval, maxind] = max(labelprobs(count+s, :));
        if maxval >= 0.99
            labels(count+s) = maxind;
        end
    else
        labelprobs(count+s, :) = 0;
    end
    weights(count+s) = features(count+s, 54);

end

count = count + nseg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pixind = getPixels(pixlist, smap, npixels, s)

ind = find(smap==s);
npix = sum(npixels(ind));
pixind = zeros(npix, 1);

pixlist = pixlist(ind);

c = 0;
for k = 1:numel(pixlist)
    nkpix = numel(pixlist{k});
    pixind(c+1:c+nkpix) = pixlist{k};
    c = c + nkpix;
end

author: zouchuhang

repo: LayoutNet


text:
function [Q, qMap] = mef_ssim(imgSeq, fI,  K, window)

if (nargin < 2 || nargin > 4)
   Q = -Inf;
   qMap = Inf;
   return;
end

if (~exist('K', 'var'))
   K = 0.03;
end

if (~exist('window', 'var'))
   window = fspecial('gaussian', 11, 1.5);
end

imgSeq = double(imgSeq);
fI = double(fI);
[s1, s2, s3] = size(imgSeq);
wSize = size(window,1);
sWindow = ones(wSize) / wSize^2; % square window used to calculate the distance
bd = floor(wSize/2);
mu = zeros(s1-2*bd, s2-2*bd, s3);
ed = zeros(s1-2*bd, s2-2*bd, s3);
for i = 1:s3
    img = squeeze(imgSeq(:,:,i));
    mu(:,:,i) = filter2(sWindow, img, 'valid');
    muSq = mu(:,:,i) .* mu(:,:,i);
    sigmaSq = filter2(sWindow, img.*img, 'valid') - muSq;
    ed(:,:,i) =  sqrt( max( wSize^2 * sigmaSq, 0 ) ) + 0.001; % add a small constant to avoid instability
end

R = zeros(s1-2*bd,s2-2*bd); % consistency map which could be used as an output if necessary
for i = bd+1:s1-bd
    for j = bd+1:s2-bd
        vecs = reshape(imgSeq(i-bd:i+bd,j-bd:j+bd,:),[wSize*wSize, s3]);
        denominator = 0;
        for k = 1:s3
            denominator = denominator + norm(vecs(:,k) - mu(i-bd,j-bd,k));
        end
        numerator = norm(sum(vecs,2) - mean(sum(vecs,2)));
        R(i-bd,j-bd) = (numerator + eps) / (denominator + eps);
    end
end

R(R > 1) = 1 - eps; % get rid of numerical instability
R(R < 0) = 0 + eps;

p = tan(pi/2 * R);
p( p >  10 ) = 10; % to avoid blow up (large number such as 10 is equivalent to taking maximum)
p = repmat(p,[1,1,s3]);

wMap = (ed / wSize).^p + eps; % to avoid blowing up
normalizer = sum(wMap,3);
wMap = wMap ./ repmat(normalizer,[1,1,s3]);
maxEd = max(ed,[],3);

C = (K * 255)^2;
qMap = zeros(s1-2*bd, s2-2*bd);
for i = bd+1:s1-bd
    for j = bd+1:s2-bd
        blocks = imgSeq(i-bd:i+bd,j-bd:j+bd,:);
        rBlock = zeros(wSize,wSize);
        for k = 1 : s3
            rBlock = rBlock  + wMap(i-bd,j-bd,k) * ( blocks(:,:,k) - mu(i-bd,j-bd,k) ) / ed(i-bd,j-bd,k);
        end
        if norm(rBlock(:)) > 0
            rBlock = rBlock / norm(rBlock(:)) * maxEd(i-bd,j-bd);
        end
        fBlock = fI(i-bd:i+bd,j-bd:j+bd);
        rVec = rBlock(:);
        fVec = fBlock(:);
        mu1 = sum( window(:) .* rVec );
        mu2 = sum( window(:) .* fVec );
        sigma1Sq = sum( window(:) .* (rVec - mu1).^2 );
        sigma2Sq = sum( window(:) .* (fVec - mu2).^2 );
        sigma12 = sum(  window(:) .* (rVec - mu1) .* (fVec - mu2)  );
        qMap(i-bd,j-bd) = ( 2 * sigma12 + C ) ./ ( sigma1Sq + sigma2Sq + C );
    end
end

Q = mean2(qMap);

author: Linfeng-Tang

repo: Image-Fusion


text:
function rgb = skin_medium

% This returns a predefined color as [red green blue] values
%   red               = [255   0   0]/255;
%   green             = [  0 192   0]/255;
%   blue              = [  0   0 255]/255;
%   magenta           = [255 255   0]/255;
%   cyan              = [  0 255 255]/255;
%   yellow            = [255 255   0]/255;
%   white             = [255 255 255]/255;
%   black             = [  0   0   0]/255;
%
%   skull             = [140  85  85]/255
%   cortex            = [255 213 119]/255;
%   cortex_light      = [199 194 169]/255;
%   cortex_dark       = [100  97  85]/255;
%   skin              = [249 223 192]/255;
%   skin_light        = [249 223 192]/255;
%   skin_medium_light = [225 194 158]/255;
%   skin_medium       = [188 142 106]/255;
%   skin_medium_dark  = [155 102	65]/255;
%   skin_dark         = [ 91  71  61]/255;
%
% The different skin-based colors follow the Fitzpatrick scale with type I and II
% combined, and return RGB values that approximate those used by Apple in the emoji
% skin tones. See also https://emojipedia.org/emoji-modifier-sequence/
%
% If no specific skin tone is specified, this function returns a light skin color.
% This corresponds with that of one of the developers who approximated his own skin
% color more than 15 years ago upon the first implementation of this function.

% returns a predefined color as [red green blue] values
rgb = [188	142	106]/255;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Script demo_set_geometry
% Illustrates the different components that define an image geometry.
%
%  demo_set_geometry
%
%
%   See also definition_of_geometry
%
% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-11-05
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get example data
dataPath = tapas_uniqc_get_path('data');
niftiFile4D = fullfile(dataPath, 'nifti', 'rest', 'fmri_short.nii');
dataLoad = MrImage(niftiFile4D);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. Create MrImage object from matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% no geometry information is supplied, defaults are used
dataRaw = MrImage(dataLoad.data); % as reference, remains unchanged
data = MrImage(dataLoad.data); % geometry and dimInfo will be adapted
dataRaw = dataRaw.select('t', 1);
data = data.select('t', 1);
data.parameters.save.fileName = 'orig.nii';
disp_centre_and_origin(data);
data.plot('plotType', 'spmi');
%   Note 1: no dimInfo/geometry information are supplied, defaults are used:
%   nSamples is derived from the data matrix
%   resolutions is assumed to be 1
%   the origin (voxel at position [0 0 0] mm) is in the center of the image

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. Add resolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot data using the classical way, but adding sampling points
% (this option is only available with the image processing toolbox)
iptsetpref('ImshowAxesVisible', 'on');
f = data.plot;
a = f.CurrentAxes;
nX = round(a.XLim(2)/data.dimInfo.nSamples(1));
xAxis = repmat(data.dimInfo.samplingPoints{1}, [1,nX]);
a.XTickLabel = xAxis(a.XTick);
nY = round(a.YLim(2)/data.dimInfo.nSamples(2));
yAxis = repmat(data.dimInfo.samplingPoints{2}, [1,nY]);
a.YTickLabel = yAxis(a.YTick);

% add additional resolution information
data.dimInfo.resolutions = dataLoad.dimInfo.resolutions;

f = data.plot;
a = f.CurrentAxes;
nX = round(a.XLim(2)/data.dimInfo.nSamples(1));
xAxis = repmat(data.dimInfo.samplingPoints{1}, [1,nX]);
a.XTickLabel = xAxis(a.XTick);
nY = round(a.YLim(2)/data.dimInfo.nSamples(2));
yAxis = repmat(data.dimInfo.samplingPoints{2}, [1,nY]);
a.YTickLabel = yAxis(a.YTick);

disp_centre_and_origin(data);
data.plot('plotType', 'spmi', 'overlayImages', dataRaw);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. Add Shear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% none of these options will affect matrix plot
data.affineTransformation.shear = [0.5 0 0];
data.plot('plotType', 'spmi', 'overlayImages', dataRaw);
disp_centre_and_origin(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4. Add Rotation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data.affineTransformation.shear = [0 0 0];
data.affineTransformation.rotation_deg = [0 30 0];
data.plot('plotType', 'spmi', 'overlayImages', dataRaw);
disp_centre_and_origin(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4a. Add Translation (offcentre_mm) in the affineTrafo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data.affineTransformation.offcenter_mm(3) = 10;
data.plot('plotType', 'spmi', 'overlayImages', dataRaw);
% world space origin is changed (but note, that, since the transformation is
% applied last, the origin changes in the two dimension which are affected by the rotation)
disp_centre_and_origin(data);
% but voxel space origin is maintained
disp(data.dimInfo.get_origin());

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4b. Change translation in dimInfo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now the translation is applied first
data2 = data.select('t', 1).copyobj();
data2.affineTransformation.offcenter_mm(3) = -10;
data2.dimInfo.set_dims(3, 'firstSamplingPoint', data2.dimInfo.samplingPoints{3}(1) + 20);
data2.plot('plotType', 'spmi', 'overlayImages', data.select('t', 1));
disp_centre_and_origin(data2);

author: translationalneuromodeling

repo: tapas


text:
function [c,b,c1] = cvx_foopsi(y,b,c1,sn,b_lb,g,w,keep)

% implementation of constrained foopsi in CVX
% Written by Eftychios Pnevmatikakis

    if isempty(b)
        bas_est = 1;
    else
        bas_est = 0;
    end
    if isempty(c1)
        c1_est = 1;
    else
        c1_est = 0;
    end
    gd = max(roots([1,-g(:)']));
    T = length(y);
    G = spdiags(ones(T,1)*[-g(end:-1:1)',1],-length(g):0,T,T);
    gd_vec = gd.^((0:T-1)');
    cvx_begin quiet
        variable c2(T)
        if bas_est; variable b; end
        if c1_est; variable c1; end
        minimize(w'*(G*c2))
        subject to
            G*c2>=0;
            norm(y(keep)-c2(keep)-b-c1*gd_vec(keep))<=sqrt(sum(keep))*sn;
            if bas_est; b>=b_lb; end
            if c1_est; c1>=0; end
    cvx_end
    if strcmpi(cvx_status,'Infeasible');
        %disp('Problem is infeasible, adjusting noise value.');
        cvx_begin quiet
            variable c2(T)
            if bas_est; variable b; end
            if c1_est; variable c1; end
            minimize(norm(y(keep)-c2(keep)-b-c1*gd_vec(keep)))
            subject to
                G*c2>=0;
                if bas_est; b>=b_lb; end
                if c1_est; c1>=0; end
        cvx_end
        sn = cvx_optval/sqrt(sum(keep));
    end
    c = c2;
end

author: flatironinstitute

repo: CaImAn-MATLAB


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

function difference = H_diff(H1, H2)
    difference = H1 \ H2;
    logR = logm(difference(1:3,1:3));
    v = [-logR(1,2) logR(1,3) -logR(2,3) difference(1:3,4)'];
    difference = norm(v);
end

author: UMich-BipedLab


text:
function[varargout]=makefigs_gulftransfer(str)
%MAKEFIGS_TRANSFER Makes all figures for Lilly and Elipot (2021).
%
%   This function makes all figures for
%
%        Lilly, J. M. and S. Elipot (2021). A unifying perspective on
%            transfer function solutions to the unsteady Ekman problem.
%            Fluids, 6 (2): 85, 1--36.
%
%   This may take a while as some of the plots are computationally
%   intensive.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2021 J.M. Lilly --- type 'help jlab_license' for details

if nargin==0
  str='noprint';
end

if strcmp(str,'--f')
     makefigs_gulfdrifters('noprint');return
end

%This is where you want things to be printed to
dirname='/Users/lilly/Desktop/Dropbox/Projects/transfer/theory/figures';

%/*************************************************************************
%basic figure of transfer function
rho=1027;
fc=1e-4;
omega=[-8.5:.01:8.5]*fc;
z=[0:.1:50]';
zo=20;
delta=20;zo=20;h=50;mu=delta.^2./zo;
G1=rho.*fc.*windtrans(omega*24*3600,z,fc*24*3600,delta,mu,inf);
G2=rho.*fc.*windtrans(omega*24*3600,z,fc*24*3600,delta,mu,50);
G1(~isfinite(G1))=1e6;
[ommat,zmat]=meshgrid(omega,z);
zeta=2*sqrt(2)*frac(zo,delta).*sqrt((1+zmat./zo).*abs(1+ommat./fc));

%str='Transfer Function Magnitude (m^2s/kg)';
str='Log10 Transfer Function Magnitude (m^{-1})';

figure,
clf
subplot(1,2,1),contourf(omega./fc,z,log10(abs(G1))',-10:.1:6),nocontours,flipy,hold on
text(-8,48,'(a)','color','w')
subplot(1,2,2),contourf(omega./fc,z,log10(abs(G2))',-10:.1:6),nocontours,flipy,hold on
text(-8,48,'(b)','color','w')

for i=1:2
    subplot(1,2,i)
    caxis([-4 -1]),hlines(20,'0.5k:'),%vlines(-1,'w'),
    %hlines([10 40],'2w')
    if i==1
        hlines([10 40],'0.5D')
    else
        hlines([40],'0.5D')
    end
    xtick([-10:2:10])
    ylabel('Depth (m)')
    xlabel('Nondimensional Frequency $\omega/|f|$','interpreter','latex')
    %contour(omega./fc,z,zeta,[4 7],'w','linewidth',2)
    contour(omega./fc,z,zeta,[4 7.5],'w')
end
subplot(1,2,1)
text(-1,5,'I','color',[1 1 1],'HorizontalAlignment','center')
text(-1,25,'IV','color',[1 1 1],'HorizontalAlignment','center')
text(-1,45,'VII','color',[1 1 1],'HorizontalAlignment','center')
text(2,5,'II','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-4,5,'II','color',[1 1 1]*0,'HorizontalAlignment','center')
text(1,25,'V','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-3,25,'V','color',[1 1 1]*0,'HorizontalAlignment','center')
text(0.4,45,'VIII','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-2.4,45,'VIII','color',[1 1 1]*0,'HorizontalAlignment','center')
text(5.5,5,'III','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-7.5,5,'III','color',[1 1 1]*0,'HorizontalAlignment','center')
text(5.5,25,'VI','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-7.5,25,'VI','color',[1 1 1]*0,'HorizontalAlignment','center')
text(5.5,45,'IX','color',[1 1 1],'HorizontalAlignment','center')
text(-7.5,45,'IX','color',[1 1 1],'HorizontalAlignment','center')

subplot(1,2,2)
%text(-1,5,'I-$h$','color',[1 1 1],'HorizontalAlignment','center')
text(-1,25,'IV-$h$','color',[1 1 1],'HorizontalAlignment','center')
text(-1,45,'VII-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
%text(2,5,'II-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
%text(-4,5,'II-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(1,25,'V-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-3,25,'V-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(0.4,45,'VIII-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-2.4,45,'VIII-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
%text(5.5,5,'III-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
%text(-7.5,5,'III-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(5.5,25,'VI-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(-7.5,25,'VI-$h$','color',[1 1 1]*0,'HorizontalAlignment','center')
text(5.5,45,'IX-$h$','color',[1 1 1],'HorizontalAlignment','center')
text(-7.5,45,'IX-$h$','color',[1 1 1],'HorizontalAlignment','center')

ha=packfig(1,2,'col');
axes(ha(1)),%ylim([0 62]),ytick([0:10:50])
hc=colorbar('SouthOutside');
hc.Label.String=str;
pos1=get(ha(1),'position');
pos2=get(ha(2),'position');
pos=get(hc,'position');%pos(4)=pos(4)/2;
%set(hc,'position',[pos(1)+0.195 pos(2) pos(3) pos(4)])
set(hc,'position',[pos(1)+0.195 pos(2)+.02 pos(3) pos(4)/2])
set(ha(2),'position',[pos2(1) pos1(2) pos2(3) pos1(4)])
set(ha(1),'position',pos1)
fontsize 10 10 9 8
set(gcf,'paperposition',[1 1 10 6])
if strcmp(str,'print')
    jprint(dirname,'transferfunctionschematic')
end
% %--------------------------------------------------------------------------
% %Green's function
% rho=1027;
% fc=1e-4;
% %omega=[-100:.01:100]*fc;
% omega=[-1000:.01:1000]*fc;
% z=[0:.1:49.9]';
% zo=20;
% delta=20;zo=20;h=50;mu=delta.^2./zo;
% %G1=rho.*fc.*windtrans(omega*24*3600,z,fc*24*3600,delta,mu,inf);
% G2=rho.*fc.*windtrans(omega*24*3600,z,fc*24*3600,delta,mu,50);
% g1=ifft(ifftshift(G1,1));
% g2=ifft(ifftshift(G2,1));
%
% figure,
% jpcolor(1:1e4,z,log10(abs(g2(1:1e4,:)))'),caxis([-7 -3.5]),flipy,hold on
% contour(1:1e4,z,real(g2(1:1e4,:))',[0 0],'k')
%\*************************************************************************

%/*************************************************************************
%basic figure of transfer function, flower version
rho=1027;
fc=1e-4;
omega=[-8.5:.01:8.5]*fc;
z=[0:5:45]';
zo=20;
delta=20;zo=20;h=50;mu=delta.^2./zo;
lnomega=logspace(2,100,10000);
omega=[-fliplr(lnomega) -100:.1:-8.6 -8.5:.001:-1-0.005 -1-0.005:1e-7:-1+0.005  -1+0.005:0.001 :8.5 8.6:.1:100 lnomega]*fc;

G1=rho.*abs(fc).*windtrans(omega*24*3600,z,fc*24*3600,delta,mu,inf);
G2=rho.*abs(fc).*windtrans(omega*24*3600,z,fc*24*3600,delta,mu,50);
figure
clf,
subplot(1,2,1)
hl=plot(G1);hold on,linestyle(hl,'0.5D'),
hl=plot(G1(:,1)); linestyle(hl,'2k')
hl=plot(G1(:,end));linestyle(hl,'2k--')
subplot(1,2,2)
hl=plot(G2);hold on,linestyle(hl,'0.5D'),
hl=plot(G2(:,1)); linestyle(hl,'2k')
hl=plot(G2(:,end));linestyle(hl,'2k--')
for i=1:2
    subplot(1,2,i)
    axis square, axis equal
    if i==1
        axis([-.1 1.5 -.8 .8]),
        xtick([-0.8:.2:1.6])
        text(1.35,0.7,['(' setstr(96+i) ')'])
    elseif i==2
        axis([-.01 .13 -.07 .07])
        text(0.117,0.06,['(' setstr(96+i) ')'])
    end

    vlines(0,'0.5k:'),hlines(0,'0.5k:')
    xlabel('Real Part of $\widetilde G(\omega,z)$ (m$^{-1}$)','interpreter','latex')%(Dimensionless)')
    ylabel('Imaginary Part of $\widetilde G(\omega,z)$ (m$^{-1}$)','interpreter','latex')% (Dimensionless)')
    %xtick([-.4:.2:1.2]),ytick([-.6:.2:.6])
end
set(gcf,'paperposition',[1 1 10 5])
if strcmp(str,'print')
    jprint(dirname,'schematicflowers')
end
%\*************************************************************************

%/*************************************************************************
%green's functions for schematic
rho=1027;
fc=1e-4;
%omega=[-100:.0001:100]*fc;
omega=[-10000:.01:10000]*fc;
z=[0:1/2:49.5]';
zo=20;
delta=20;zo=20;h=50;mu=delta.^2./zo;

[G1,G2]=vzeros(length(omega),length(z));

%z=[0:5:45]';
parfor i=1:length(z)
    i
%   G1(:,i)=rho.*abs(fc).*windtrans(omega*24*3600,z(i),fc*24*3600,delta,mu,inf);
    G2(:,i)=rho.*abs(fc).*windtrans(omega*24*3600,z(i),fc*24*3600,delta,mu,50);
end
G2(:,end+1)=1e-16+0*G2(:,end);
z=[z;50];
%g1=ifft(ifftshift(G1,1));
g2=ifft(ifftshift(G2,1));
g2=g2./maxmax(abs(g2));

%fc=2 pi rad/ 10^4 seconds
%fN= 2 pi rad
%dt= 1/2 seconds
%Tf = 2 pi / fc = 10^4 seconds
t=[1:size(g2,1)]'*frac(1,2)./10^4;

figure
ii=1:40000;
subplot(3,1,1)
contourf(t(ii),z,log10(abs(g2(ii,:)))',[-10:.1:0]);nocontours,hold on,flipy
contour(t(ii),z,log10(abs(g2(ii,:)))',[-12:1:0],'w');
hlines(z(1:10:end),'0.2D')
hlines(z(1),'2k')
hlines(z(end-10),'2k--')
%contour(t(ii),z,real(g2(ii,:))',[0 0],'k:');
caxis([-5.5 -1.5]),ylim([-.01 49.999]),xlim([0 2]),xtick([0:.25:2])
vlines(0.25005+[0:1/2:2],'0.25k:')
ylabel('Depth (m)')
text(1.875,2.5,'(a)','color','k')
subplot(3,1,2)
h=plot(t(ii),real(g2(ii,1:10:end))*1000);ylim(1000*[-.009 .009]+1000*0.005)
linestyle(h,'0.5D'),linestyle(h(1),'2k'),linestyle(h(end-1),'2k--')
ylabel('Real Part of $g(t,z)$ ($\times 1000 / \mathrm{max}\{g\}$)','interpreter','latex')
text(1.875,12.7,'(b)','color','k')
hold on,vlines(0.25005+[0:1/2:2],'0.25k:'),xtick([0:.25:2])
%h=plot(t(ii),abs(g2(ii,1:10:end))*1000);
subplot(3,1,3)
h=plot(t(ii),imag(g2(ii,1:10:end))*1000);ylim(1000*[-.009 .009])
linestyle(h,'0.5D'),linestyle(h(1),'2k'),linestyle(h(end-1),'2k--')
ylabel('Imaginary Part of $g(t,z)$ ($\times 1000 / \mathrm{max}\{g\}$)','interpreter','latex')
xlabel('Time (Inertial Periods $2\pi / f$)','interpreter','latex')
text(1.875,8,'(c)','color','k')
hold on,vlines(0.25005+[0:1/2:2],'0.25k:'),xtick([0:.25:2])
packfig(3,1)

set(gcf,'paperposition',[1 1 5 10])
if strcmp(str,'print')
    jprint(dirname,'greensfunction')
end
%\*************************************************************************

if 0
%/*************************************************************************
%green's functions for schematic
rho=1027;
fc=1e-4;
%omega=[-100:.0001:100]*fc;
omega=[-10000:.01:10000]*fc;
z=[0:1/2:49.5]';
delta=20;zo=20;h=50;mu=delta.^2./zo;
delta=10;zo=10;h=50;mu=delta.^2./zo;

[G1,G2]=vzeros(length(omega),length(z));

%z=[0:5:45]';
parfor i=1:length(z)
    i
   % G1(:,i)=rho.*abs(fc).*windtrans(omega*24*3600,z(i),fc*24*3600,delta,mu,inf);
    G2(:,i)=rho.*abs(fc).*windtrans(omega*24*3600,z(i),fc*24*3600,delta,mu,50);
end
%G1(~isfinite(G1))=0;
G2(:,end+1)=1e-16+0*G2(:,end);
z=[z;50];
%g1=ifft(ifftshift(G1,1));
g2=ifft(ifftshift(G2,1));
g2=g2./maxmax(abs(g2));

%fc=2 pi rad/ 10^4 seconds
%fN= 2 pi rad
%dt= 1/2 seconds
%Tf = 2 pi / fc = 10^4 seconds
t=[1:size(g2,1)]'*frac(1,2)./10^4;

clf
ii=1:40000;
subplot(3,1,1)
contourf(t(ii),z,log10(abs(g2(ii,:)))',[-12:.1:0]);nocontours,hold on,flipy
contour(t(ii),z,log10(abs(g2(ii,:)))',[-12:1/2:0],'w');
hlines(z(1:10:end),'0.2D'),hlines(z(1),'2k'),hlines(z(end-10),'2k--')
%contour(t(ii),z,real(g2(ii,:))',[0 0],'k:');
caxis([-5.5 -1.5]),ylim([-.01 49.999]),xlim([0 2]),xtick([0:.25:2])
vlines(0.25005+[0:1/2:2],'0.25k:')
ylabel('Depth (m)')
text(1.875,2.5,'(a)','color','k')
subplot(3,1,2)
h=plot(t(ii),real(g2(ii,1:10:end))*1000);ylim(1000*[-.009 .009]+1000*0.005-1.5)
linestyle(h,'0.5D'),linestyle(h(1),'2k'),linestyle(h(end-1),'2k--')
ylabel('Real Part of $g(t,z)$ ($\times 10^3 / \mathrm{max}\{g\}$)','interpreter','latex')
text(1.875,12.7-1.5,'(b)','color','k')
hold on,vlines(0.25005+[0:1/2:2],'0.25k:'),xtick([0:.25:2])
%h=plot(t(ii),abs(g2(ii,1:10:end))*1000);
subplot(3,1,3)
h=plot(t(ii),imag(g2(ii,1:10:end))*1000);ylim(1000*[-.009 .009])
linestyle(h,'0.5D'),linestyle(h(1),'2k'),linestyle(h(end-1),'2k--')
ylabel('Imaginary Part of $g(t,z)$ ($\times 10^3 / \mathrm{max}\{g\}$)','interpreter','latex')
text(1.875,8,'(c)','color','k')
hold on,vlines(0.25005+[0:1/2:2],'0.25k:'),xtick([0:.25:2])
xlabel('Time (Inertial Periods $2\pi / f$)','interpreter','latex')
% subplot(4,1,4)
% %gke=5*fliplr(log10(cumsum(squared(fliplr(g2)),2)));
% gke=5*log10(cumsum(squared(g2),2));
% contourf(t(ii),z,gke(ii,:)',[-50:.5:0]);nocontours,hold on,flipy
% contour(t(ii),z,gke(ii,:)',[-50:5:0],'w');
% hlines(z(1:10:end),'0.2D'),hlines(z(1),'2k'),hlines(z(end-10),'2k--')
% caxis([-40 -5]),ylim([-.01 49.999]),xlim([0 2]),xtick([0:.25:2])
% vlines(0.25005+[0:1/2:2],'0.25k:')
% text(1.875,2.5,'(d)','color','k')
% %ylabel('Kinetic Energy $\int_{-h}^z |g(t,x)|^2 \mathrm{d}x$')
% ylabel('Kinetic Energy $\int_{0}^z |g(t,x)|^2 \mathrm{d}x$')
% xlabel('Time (Inertial Periods $2\pi / f$)','interpreter','latex')
% packfig(4,1)
%set(gcf,'paperposition',[1 1 5 13])

packfig(3,1)
set(gcf,'paperposition',[1 1 5 10])
if strcmp(str,'print')
    jprint(dirname,'greensfunction_broader')
end
%\*************************************************************************
end

%/*************************************************************************
%illustration of problems with overflow
log10delta=linspace(-1,4,100);
log10mu=linspace(-7,4,100);
h=100;
%h=1000;
%h=1e6;
%h=1e8;
fc=1e-4;

Gi1=vzeros(length(log10mu),length(log10delta));
Gi2=vzeros(length(log10mu),length(log10delta));
Gi3=vzeros(length(log10mu),length(log10delta));
Gi4=vzeros(length(log10mu),length(log10delta));
%f=0;
f=2.*fc*24*3600;
%f=fc*24*3600*0.99;
%f=fc*24*3600*8;
tic
for i=1:length(log10mu)
    i
    for j=1:length(log10delta)
     Gi1(i,j)=windtrans(f,15,fc*24*3600,10.^log10delta(j),10.^log10mu(i),h,'two');
     Gi2(i,j)=windtrans(f,15,fc*24*3600,10.^log10delta(j),10.^log10mu(i),h,'far');
     Gi3(i,j)=windtrans(f,15,fc*24*3600,10.^log10delta(j),10.^log10mu(i),h,'general');
     Gi4(i,j)=windtrans(f,15,fc*24*3600,10.^log10delta(j),10.^log10mu(i),h,'expansion');
    end
end
%etime1=toc
%etime2=toc

figure,
clf
subplot(1,3,1),contourf(log10delta,log10mu,log10(abs(Gi3)),[-200:.1:0]),nocontours,hold on
contour(log10delta,log10mu,log10(abs(Gi3)),[-200:10:0],'w'),caxis([-12 maxmax(log10(abs(Gi1)))])
hc=colorbar('South');
hc.Label.String='Log10 Magnitude at $\omega = -2f$';
hc.Label.Interpreter='latex';
pos=get(hc,'position');
set(hc,'position',[pos(1)+pos(1)/10 pos(2) pos(3)+pos(1)/10 pos(4)/2])
subplot(1,3,2),contourf(log10delta,log10mu,log10(abs(Gi1)),[-200:.1:0]),nocontours,hold on
contour(log10delta,log10mu,log10(abs(Gi1)),[-200:10:0],'w'),caxis([-12 maxmax(log10(abs(Gi1)))])
subplot(1,3,3),contourf(log10delta,log10mu,log10(abs(Gi2-Gi4)./abs(Gi4)),[-14:.1:0]),nocontours,hold on
caxis([-14 -3])
hc2=colorbar('North');
hc2.Label.String='Log10 Fractional Error at $\omega = -2f$';
hc2.Label.Interpreter='latex';
pos2=get(hc2,'position');
set(hc2,'position',[pos2(1)-pos2(1)/28 pos2(2)-0.03 pos2(3)+pos(1)/10 pos2(4)/2])

[demat,mumat]=meshgrid(10.^log10delta,10.^log10mu);
for i=1:3
    subplot(1,3,i),%contour(log10delta,log10mu,log10(abs(Gi1-Gi2)./abs(Gi1)),[-3.4 -3.4],'k','linewidth',2)
    text(3.5,3.7,['(' setstr(96+i) ')'],'color','k')
    xlim([-1 3.99])
    axis equal,xlim([-1 3.99]),ylim([-5 4])
    ylabel('Log10 Madsen Depth $\mu$ (m)','interpreter','latex')
    xlabel('Log10 Ekman Depth $\delta$ (m)','interpreter','latex')
    %contour(log10delta,log10mu,log10(demat./mumat*2*sqrt(2)*sqrt(3)),[1 1]*3,'color',[1 1 1]*0.6,'linewidth',2)
    contour(log10delta,log10mu,log10(demat./mumat*2*sqrt(2)*sqrt(3)),[1 1]*2.9,'color',[1 1 1]*0.6,'linewidth',2)
end
packfig(1,3,'columns')
set(hc,'position',[pos(1)+pos(1)/10 pos(2)+.1 pos(3)+pos(1)/10 pos(4)/2])
set(hc2,'position',[pos2(1)-pos2(1)/28 pos2(2)-0.03-.07 pos2(3)+pos(1)/10 pos2(4)/2])
%set(gcf,'paperposition',[1 1 11 5])
set(gcf,'paperposition',[1 1 11 8])
if strcmp(str,'print')
    jprint(dirname,'overflowproblem')
end
%\*************************************************************************

%/*************************************************************************
%symmetry plots
rho=1027;
%log10delta=linspace(-1,4,101);
%log10mu=linspace(-7,4,100);
log10delta=linspace(-4,4,201);
log10mu=linspace(-8,4,200);
%h=10000;
h=[16 1e2 1e3 1e4 1e5 1e6];
%h=[15.01 15.1 16 26 1e2 1e3 1e4 1e5 1e6];
fc=1e-4;

[demat,mumat]=meshgrid(10.^log10delta,10.^log10mu);
zomat=demat.^2./mumat;

A=zeros(size(demat,1),size(demat,2),length(h));
A0=zeros(size(demat,1),size(demat,2),length(h));
for i=1:length(h)
    A(:,:,i)=frac(2,mumat).*log(frac(1+h(i)./zomat,1+15./zomat));
    A0(:,:,i)=frac(2,mumat).*log(1+h(i)./zomat);
end

%Gf=vzeros(length(log10mu),length(log10delta),length(h));
%Gf0=vzeros(length(log10mu),length(log10delta),length(h));
% for i=1:length(log10mu)
%     i
%     parfor j=1:length(log10delta)
%         for k=1:6
%           Gf(i,j,k)=rho.*fc.*windtrans(-fc*24*3600,15,fc*24*3600,10.^log10delta(j),10.^log10mu(i),h(k));
%          Gf0(i,j,k)=rho.*fc.*windtrans(-fc*24*3600,0,fc*24*3600,10.^log10delta(j),10.^log10mu(i),h(k));
%         end
%     end
% end
% maxmax(abs(A-Gf))
% maxmax(abs(A0-Gf0))
% jpcolor(log10(abs(A(:,:,1)-Gf(:,:,1))./A(:,:,1)))
% %these agree to the numerical noise level

omega=[-3:.001:1]*fc;
log10zo=[-12:.1:12]';
ii=1:10:length(log10zo);

clear dei mui dei0 mui0 dei1 mui1
zo=10.^log10zo;
for j=1:length(h)
    dei{j,1}=sqrt(frac(2*zo,1).*log(frac(1+h(j)./zo,1+15./zo)));
    mui{j,1}=dei{j}.^2./zo;
    dei1{j,1}=sqrt(frac(2*zo,10).*log(frac(1+h(j)./zo,1+15./zo)));
    mui1{j,1}=dei1{j}.^2./zo;
    dei0{j,1}=sqrt(frac(2*zo,1).*log(frac(1+h(j)./zo,1)));
    mui0{j,1}=dei0{j}.^2./zo;
end
%--------------------------------------------------------------------------
figure,
clf
contourf(log10delta,log10mu,log10(A(:,:,3)),[-10:.1:10]),nocontours,hold on
contour(log10delta,log10mu,log10(A(:,:,3)),[0 0 ],'k','linewidth',1.5)
contour(log10delta,log10mu,log10(A(:,:,3)),[-10:1:10],'k','linewidth',0.5)
contour(log10delta,log10mu,log10(zomat),[-20:20],'w','linewidth',0.5)
contour(log10delta,log10mu,log10(zomat),[0 0],'w','linewidth',1.5)
plot(log10(dei{3}(1:10:end)),log10(mui{3}(1:10:end)),'ko','markerfacecolor','k','markersize',4)
caxis([-4.7 7.3]),axis equal,ylim([-7 4]),xlim([-4 4]),
ylabel('Log10 Madsen Depth $\mu$ (m)','interpreter','latex')
xlabel('Log10 Ekman Depth $\delta$ (m)','interpreter','latex')
hc=colorbar('SouthOutside');
hc.Label.String='Log10 Inertial Amplitude A (m$^{-1}$)';
hc.Label.Interpreter='latex';
hc.Ticks=[-4:1:8];
set(gcf,'paperposition',[1 1 4*1.5 5*1.5])
%text(-0.92,3.7,'M'),text(3.8,-6.8,'E')
text(-3.925,3.7,'M'),text(3.7,-6.8,'E','color','w')
if strcmp(str,'print')
    jprint(dirname,'deltamuplane')
end
%--------------------------------------------------------------------------
%figure,cellplot(dec,muc)
%figure,cellplot(dec0,muc0)
omega=[-100:.1:-8.6 -8.5:.001:-1-0.005 -1-0.005:1e-7:-1+0.005  -1+0.005:0.001:8.5 8.6:.1:100]*fc;
%omega=[-1e4:10:-1000 -1000:1:-100 -100:.1:-8.6 -8.5:.001:-1-0.005 -1-0.005:1e-7:-1+0.005  -1+0.005:0.001 :8.5 8.6:.1:100 100:1:1000 1000:10:1e4]*fc;
G=nan*zeros(length(omega),length(ii),length(h));
[Gm,Ge,Gminf,Geinf]=vzeros(length(omega),length(h));
%G0=nan*zeros(length(omega),length(ii),length(h));
%G1=nan*zeros(length(omega),length(ii),length(h));
for j=1:size(G,3)
    j
   Gm(:,j)=rho.*abs(fc).*windtrans(omega*24*3600,15,fc*24*3600,0,mui{j}(1),h(j));
   Ge(:,j)=rho.*abs(fc).*windtrans(omega*24*3600,15,fc*24*3600,dei{j}(end),0,h(j));
   parfor k=1:length(ii)
        k
        G(:,k,j)=rho.*abs(fc).*windtrans(omega*24*3600,15,fc*24*3600,dei{j}(ii(k)),mui{j}(ii(k)),h(j));
         %G1(:,k,j)=rho.*abs(fc).*windtrans(omega*24*3600,15,fc*24*3600,dei1{j}(ii(k)),mui1{j}(ii(k)),h(j));
        %G0(:,k,j)=rho.*abs(fc).*windtrans(omega*24*3600,0,fc*24*3600,dei0{j}(ii(k)),mui0{j}(ii(k)),h(j));
    end
end
%--------------------------------------------------------------------------
figure
clf
%[~,jj]=min(abs(omega));
for j=1:size(G,3)
    subplot(2,3,j)
    hl=plot(G(:,:,j));axis equal, axis square
    axis([-.3 1 -.65 .65]*1.05),vlines(0,'0.5k:'),hlines(0,'0.5k:')
    linestyle(hl,'0.5D')
    hl=plot(G(:,1,j));linestyle(hl,'2k--')
    hl=plot(G(:,end,j));linestyle(hl,'2k')
    hl=plot(Ge(:,j));linestyle(hl,'0.5w')
    hl=plot(Gm(:,j));linestyle(hl,'0.5k')
    text(0.9,0.59,['(' setstr(96+j) ')'])
%    xlabel('Real Part of Transfer Function (m^2s / kg)')%(Dimensionless)')
%    ylabel('Imaginary Part of Transfer Function (m^2s / kg)')% (Dimensionless)')
    xlabel('Real Part of $\widetilde G(\omega,15$ m$)$ (m$^{-1}$)','interpreter','latex')%(Dimensionless)')
    ylabel('Imaginary Part of $\widetilde G(\omega,15$ m$)$ (m$^{-1}$)','interpreter','latex')% (Dimensionless)')
    xtick([-.4:.2:1.2]),ytick([-.6:.2:.6])
    %plot(G(jj,:,j),'ko','markerfacecolor','k','markersize',4)
end
packfig(2,3)
set(gcf,'paperposition',[1 1 9.4 6])
fontsize 10 10 9 9
if strcmp(str,'print')
    jprint(dirname,'similarflowers')
end
%%check invariance
%figure,plot(G(:,:,1)),hold on,plot(G1(:,:,1)/10,'r')
%figure,plot(G(:,:,end)),hold on,plot(G1(:,:,end)/10,'r')
% %--------------------------------------------------------------------------
% figure
% clf
% for j=1:size(Gf,3)
%     subplot(2,3,j)
%     hl=plot(omega./fc,abs(G(:,2:end-1,j)));hold on
%     axis([-1.99 0 1.001*1e-4 1.25]),ylog
%     linestyle(hl,'0.5D')
%     hl=plot(omega./fc,abs(G(:,1,j)));
%     linestyle(hl,'2k')
%     hl=plot(omega./fc,abs(G(:,end,j)));
%     linestyle(hl,'2k--')
%     ylabel('Log10 Magnitude of $G(\omega,15$ m$)$','interpreter','latex')
%     xlabel('Nondimensional Frequency $\omega/f$','interpreter','latex')
%     text(-.17,0.8,['(' setstr(96+j) ')'])
%     xtick([-2:.5:0])
%     vlines(-1,'0.5k:')
% end
% packfig(2,3)
% set(gcf,'paperposition',[1 1 9*1.25 5*1.25])
% jprint(dirname,'similarmagnitudes')
% %--------------------------------------------------------------------------
% figure
% clf
% for j=1:size(Gf,3)
%     subplot(2,3,j)
%     hl=plot(omega./fc,(360/2/pi)*unwrap(angle(G(:,2:end-1,j))));hold on
%     axis([-1.99 0 -179.9 180]),
%     linestyle(hl,'0.5D')
%     hl=plot(omega./fc,(360/2/pi)*unwrap(angle(G(:,1,j))));
%     linestyle(hl,'2k')
%     hl=plot(omega./fc,(360/2/pi)*unwrap(angle(G(:,end,j))));
%     linestyle(hl,'2k--')
%     ylabel('Angle of $G(\omega,15$ m$) (degrees) $','interpreter','latex')
%     xlabel('Nondimensional Frequency $\omega/f$','interpreter','latex')
%     %text(-.17,0.8,['(' setstr(96+j) ')'])
%     xtick([-2:.5:0])
%     vlines(-1,'0.5k:')
% end
% packfig(2,3)
% set(gcf,'paperposition',[1 1 9*1.25 5*1.25])
% jprint(dirname,'similarphases')
%--------------------------------------------------------------------------
% figure
% clf
% for j=1:size(Gf0,3)
%     subplot(2,3,j)
%     hl=plot(G0(:,2:end-1,j));axis equal, axis square
%     axis([-.05 1.05 -.505 .505]),vlines(0,'0.5k:'),hlines(0,'0.5k:')
%     linestyle(hl,'0.5D'),
%     hl=plot(G0(:,end,j)); linestyle(hl,'2k')
%     hl=plot(G0(:,1,j));linestyle(hl,'2k--')
%     text(0.9,0.44,['(' setstr(96+j) ')'])
%     xlabel('Real Part of $G(\omega,0$ m$)$','interpreter','latex')%(Dimensionless)')
%     ylabel('Imaginary Part of $G(\omega,0$ m$)$','interpreter','latex')% (Dimensionless)')
%     xtick([-.4:.2:1.2]),ytick([-.6:.2:.6])
% end
% packfig(2,3)
% set(gcf,'paperposition',[1 1 9.4 6])
% jprint(dirname,'similarflowers0')
lnomega=logspace(2,100,10000);
omega=[-fliplr(lnomega) -100:.1:-8.6 -8.5:.001:-1-0.005 -1-0.005:1e-7:-1+0.005  -1+0.005:0.001 :8.5 8.6:.1:100 lnomega]*fc;
G0=vzeros(length(omega),length(ii));
Gm0=rho.*abs(fc).*windtrans(omega*24*3600,0,fc*24*3600,0,mui0{3}(1),h(3));
Ge0=rho.*abs(fc).*windtrans(omega*24*3600,0,fc*24*3600,dei0{3}(end),0,h(3));
%Gm0=rho.*abs(fc).*windtrans(omega*24*3600,1e-10,fc*24*3600,0,mui0{3}(ii(end)),h(3));
%Ge0=rho.*abs(fc).*windtrans(omega*24*3600,0,fc*24*3600,dei0{3}(ii(1)),0,h(3));
%zeta=2*sqrt(2)*frac(zo,delta).*sqrt(abs(1+ommat./fc));

for k=1:length(ii)
    k
    G0(:,k)=rho.*abs(fc).*windtrans(omega*24*3600,0,fc*24*3600,dei0{3}(ii(k)),mui0{3}(ii(k)),h(3));
 %   G1(:,k)=rho.*abs(fc).*windtrans(omega*24*3600,0,fc*24*3600,dei0{3}(ii(k)),mui0{3}(ii(k)),h(1));
end
figure
clf
hl=plot(G0);axis equal, axis square
axis([-.05 1.05 -.505 .505]),vlines(0,'0.5k:'),hlines(0,'0.5k:')
linestyle(hl,'0.5D'),
hl=plot(G0(:,1)); linestyle(hl,'2k--')
hl=plot(G0(:,end)); linestyle(hl,'2k')
hl=plot(Ge0); linestyle(hl,'0.5w')
hl=plot(Gm0);linestyle(hl,'0.5k')
%text(0.9,0.44,['(' setstr(96+j) ')'])
xlabel('Real Part of $\widetilde G(\omega,0$ m$)$ (m$^{-1}$)','interpreter','latex')%(Dimensionless)')
ylabel('Imaginary Part of $\widetilde G(\omega,0$ m$)$ (m$^{-1}$)','interpreter','latex')% (Dimensionless)')
xtick([-.4:.2:1.2]),ytick([-.6:.2:.6])
set(gcf,'paperposition',[1 1 4 4])
plot(2*h(3)/squared(dei0{3}(end)),0,'wo','markerfacecolor','k')
fontsize 10 10 9 9
if strcmp(str,'print')
    jprint(dirname,'similarflowers0')
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [Y, X] = minandmax2( f )
%MINANDMAX2   Find global minimum and maximum of a SEPARABLEAPPROX.
%   Y = minandmax2(F) returns the minimum and maximum value of a SEPARABLEAPPROX over
%   its domain. Y is a vector of length 2 such that Y(1) = min(f(x,y)) and Y(2)
%   = max(f(x,y)).
%
%   [Y, X] = minandmax2(F) also returns the position of the minimum and
%   maximum. For example,
%
%       F(X(1,1),X(1,2)) = Y(1)     and      F(X(2,1),X(2,2)) = Y(2)
%
% See also MAX2, MIN2, NORM.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% check for empty SEPARABLEAPPROX.
if ( isempty( f ) )
    Y = [];
    X = [];
    return
end

% This is a default implementation for separableApprox children that do not
% have know how to properly implement this method.

% Maximum possible sample matrix size:
maxsize = 4e3;

% Is the function the zero function?
if ( iszero( f )  )
    dom = f.domain;
    X = [ (dom(2) + dom(1))/2 (dom(4) + dom(3))/2 ];
    X = [ X ; X ];
    Y = [0 ; 0];
    return;
end

% Extract low rank representation:
frows = f.rows;
fcols = f.cols;
piv = f.pivotValues;
dom = f.domain;

% Convert rows and columns to chebfuns.
frows = simplify(chebfun(frows, dom(1:2), maxsize));
fcols = simplify(chebfun(fcols, dom(3:4), maxsize));

% Share out scaling:
sgn = sign( piv ).';
sq = 1 ./ sqrt( abs( piv ) );
frows = frows * diag( sq.'.*sgn );
fcols = fcols * diag( sq );

if ( length(f) == 1 ) % Rank-1 is easy:
    % We can find it from taking maximum and minimum in x and y direction.

    % Find minandmax of rows and columns:
    [yr, xr] = minandmax( frows );
    [yc, xc] = minandmax( fcols );
    % All possible combinations:
    vv = [yr(1)*yc(1), yr(1)*yc(2), yr(2)*yc(1), yr(2)*yc(2)];

    [Y(2), indmx] = max( vv );
    [Y(1), indmn] = min( vv );

    % Work out the location of the maximum.
    X = zeros(2);
    X(1,1) = xr(2);
    X(1,2) = xc(2);

    if ( indmn <= 2 )
        X(1,1) = xr(1);
    end
    if ( mod(indmn,2) == 1 )
        X(1,2) = xc(1);
    end
    X(2,1) = xr(2);
    X(2,2) = xc(2);

    if ( indmx <= 2 )
        X(2,1) = xr(1);
    end

    if ( mod(indmx,2) == 1 ),
        X(2,2) = xc(1);
    end

elseif ( length(f) <= maxsize )

    % We seek a fast initial guess. So we first truncate the SEPARABLEAPPROX.
    ypts = chebpts(length(fcols), fcols.domain);
    xpts = chebpts(length(frows), frows.domain);
    cvals = feval(fcols, ypts);
    rvals = feval(frows, xpts);

    A = cvals*rvals.';
    % Maximum entry in discretisation.
    [ignored, ind] = min( A(:) ); %#ok<ASGLU>
    [row, col] = ind2sub(size(A), ind);
    X(1,1) = xpts( col );
    X(1,2) = ypts( row );
    Y(1) = feval( f, X(1,1), X(1,2) );
    % Minimum entry in discretisation.
    [ignored, ind] = max(A(:)); %#ok<ASGLU>
    [row, col] = ind2sub(size(A), ind);
    X(2,1) = xpts(col);
    X(2,2) = ypts(row);
    Y(2) = feval(f, X(2,1), X(2,2));

    % Get more digits with optimisation algorithms.
    lb = [ dom(1) ; dom(3) ];
    ub = [ dom(2) ; dom(4) ];

    try

        % If the optimization toolbox is available then use it to get a better maximum.

        warnstate = warning;
        warning('off'); %#ok<WNOFF> % Disable verbose warnings from fmincon.
        options = optimset('Display', 'none', 'TolFun', eps, 'TolX', eps, ...
            'algorithm', 'active-set');
        [mn, Y(1)] = fmincon(@(x,y) feval(f, x(1), x(2)), X(1, :), ...
            [], [], [], [], lb, ub, [], options);
        [mx, Y(2)] = fmincon(@(x) -feval(f, x(1), x(2)), X(2,:), ...
            [], [], [], [], lb, ub, [], options);
        Y(2) = -Y(2);
        X(1,:) = mn;
        X(2,:) = mx;
        warning(warnstate);

    catch

        try
            % Try converting to an unconstrained problem and using built-in solver.

            % Maps from [-1, 1] to [dom(1:2)] and [dom(3:4)], respectively.
            map1 = bndfun.createMap(dom(1:2));
            map2 = bndfun.createMap(dom(3:4));
            % Unconstrained initial guesses:
            Z(:,1) = asin(map1.Inv(X(:,1)));
            Z(:,2) = asin(map2.Inv(X(:,2)));
            % Maps from R to [dom(1), dom(2)] and [dom(3), dom(4)], respectively.
            map1 = @(x) map1.For(sin(x));
            map2 = @(x) map2.For(sin(x));
            % Set options:
            options = optimset('Display', 'off', 'TolFun', eps, 'TolX', eps);
            warnstate = warning;
            warning('off'); %#ok<WNOFF> % Disable verbose warnings from fminsearch.
            f_mapped = @(x) feval(f, map1(x(1)), map2(x(2)));
            [mn, Y(1)] = fminsearch(@(x) f_mapped(x), Z(1, :), options);
            [mx, Y(2)] = fminsearch(@(x) -f_mapped(x), Z(2, :), options);
            Y(2) = -Y(2);
            X(1:2,1) = map1([mn(1) ; mx(1)]);
            X(1:2,2) = map2([mn(2) ; mx(2)]);
            warning(warnstate);

        catch

            % Nothing is going to work so initial guesses will have to do.

        end

    end

elseif ( length(f) >= maxsize )

    error('CHEBFUN:SEPARABLEAPPROX:minandmax2:length', 'Rank is too large.');

end

end

%%%
% Use the approach below when bivariate rootfinding is fully implemented.
%%%
% % Use bivariate rootfinding to find all the local extrema:
% F = gradient( f );
% r = roots( F );
% if ( ~isempty( r ) )
%     [inMax, idInMax] = max( feval(f, r(:,1), r(:,2) ) );
%     [inMin, idInMin] = min( feval(f, r(:,2), r(:,2) ) );
% else
%     inMax = inf;   % max and min must occur on boundary.
%     inMin = inf;
% end
%
% % Search along boundary:
% dom = f.domain;
% left = feval(f, dom(1), ':');
% right = feval(f, dom(2), ':');
% down = feval(f, ':', dom(3));
% up = feval(f, ':', dom(4));
% [Yleft, Xleft] = minandmax( left );
% [Yright, Xright] = minandmax( right );
% [Yup, Xup] = minandmax( up );
% [Ydown, Xdown] = minandmax( down );
%
% % Store Min/Max location for later:
% BcMinLocations = [ Xleft(1,:) ; Xright(1,:) ; Xup(1,:) ; Xdown(1,:) ];
% BcMaxLocations = [ Xleft(2,:) ; Xright(2,:) ; Xup(2,:) ; Xdown(2,:) ];
%
% [ BcMax, idBcMax ] = max( [ Yleft(2), Yright(2), Yup(2), Ydown(2) ].' );
% [ BcMin, idBcMin ] = min( [ Yleft(1), Yright(1), Yup(1), Ydown(1) ].' );
%
% % What is the global min and max?:
% [Ymax, inOrOutMax] = max( [inMax, BcMax].' );
% [Ymin, inOrOutMin] = min( [inMin, BcMin].' );
% Y = [Ymin, Ymax];
% X = zeros(2);
%
% % Unravel to find locations:
% if ( inOrOutMin == 1 )
%     X(1,:) = [ r(idInMin,1), r(idInMin,2) ];
% else
%     X(1,:) = BcMinLocations(idBcMin, :);
% end
% if ( inOrOutMax == 1 )
%     X(2,:) = [ r(idInMax,1), r(idInMax,2) ];
% else
%     X(2,:) = BcMaxLocations(idBcMax, :);
% end
%
% end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function calpak_test336 ( )

%*****************************************************************************80
%
%% CALPAK_TEST336 tests MONTH_LENGTH_EG_CIVIL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2013
%
%  Author:
%
%    John Burkardt
%
  n_test = 2

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST336\n' );
  fprintf ( 1, '  For the Egyptian Civil calendar,\n' );
  fprintf ( 1, '  MONTH_LENGTH_EG_CIVIL returns month lengths.\n' );

  y_test(1) = 3;
  y_test(2) = 4;

  for i_test = 1 : n_test

    y = y_test(i_test);
    sy = y_to_s_eg_civil ( y );
    months = year_length_months_eg_civil ( y );
    days = year_length_eg_civil ( y );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %d\n', y );
    fprintf ( 1, '  %s\n', sy );
    fprintf ( 1, '  Year length in months = %d\n', months );
    fprintf ( 1, '  Year length in days = %d\n', days );
    fprintf ( 1, '\n' );

    for m = 1 : months
      month_name = month_to_month_name_eg_civil ( m );
      fprintf ( 1, '  %10s  %2d\n', month_name, month_length_eg_civil ( y, m ) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fname] = gtopo30s(latlim,lonlim)

    %GTOPO30S (30-arc-sec resolution) DEM file names
    %
    % fname = GTOPO30S(LATLIM,LONLIM) returns a cellarray of the file
    % names covering the geographic region for GTOPO30 digital elevation maps
    % (also referred to as "30-arc second" DEMs).  The region is specified by
    % scalar latitude and longitude points, or two element vectors of latitude
    % and longitude limits in units of degrees.
    %
    % The data is available over the Internet via anonymous ftp from
    % <ftp://edcftp.cr.usgs.gov/pub/data/gtopo30/global>. The data and
    % some documentation is also available over the World-Wide-Web from
    % <http://edcwww.cr.usgs.gov/landdaac/gtopo30/gtopo30.html> and
    % <http://edcwww.cr.usgs.gov/landdaac/gtopo30/README.html>.
    %
    % See also: GTOPO30

    %  Written by:  A. Kim, W. Stumpf, L. Job
    %  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
    %  $Revision: 1399 $ $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    % ensure row vectors
    latlim = latlim(:)';
    lonlim = lonlim(:)';

    if nargin~=2
        error('Incorrect number of arguments')
    end

    if  isequal(size(latlim),[1 1])
        latlim = latlim*[1 1];
    elseif ~isequal(size(latlim),[1 2])
        error('Latitude limit input must be a scalar or 2 element vector')
    end

    if isequal(sort(size(lonlim)),[1 1])
        lonlim = lonlim*[1 1];
    elseif ~isequal(sort(size(lonlim)),[1 2])
        error('Longitude limit input must be a scalar or 2 element vector')
    end

    fid = fopen('gtopo30s.dat','r');
    if fid==-1
        error('Couldn''t open gtopo30s.dat')
    end

    % preallocate bounding rectangle data for speed

    YMIN = zeros(1,33); YMAX = YMIN;
    XMIN = YMIN; XMAX = YMIN;

    % read names and bounding rectangle limits

    for n=1:33
        fnames{n,1} = fscanf(fid,'%s',1);
        YMIN(n) = fscanf(fid,'%d',1);
        YMAX(n) = fscanf(fid,'%d',1);
        XMIN(n) = fscanf(fid,'%d',1);
        XMAX(n) = fscanf(fid,'%d',1);
    end
    fclose(fid);

    % case where dateline is not crossed
    if lonlim(1) <= lonlim(2)
        do = ...
            find( ...
            (...
            (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
            (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
            (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
            (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
            ) ...
            &...
            (...
            (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
            (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
            (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
            (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
            )...
            );
    end

    % case where the dateline is crossed
    if lonlim(1) > lonlim(2)
        lmin = lonlim(1); lmax = lonlim(2);
        lonlim(2) = 180;
        % do eastern side of the dateline first
        doEAST = ...
            find( ...
            (...
            (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
            (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
            (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
            (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
            ) ...
            &...
            (...
            (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
            (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
            (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
            (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
            )...
            );
        % do western side of the dateline second
        lonlim(1) = -180; lonlim(2) = lmax;
        doWEST = ...
            find( ...
            (...
            (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
            (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
            (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
            (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
            ) ...
            &...
            (...
            (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
            (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
            (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
            (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
            )...
            );
        % concatenate indices
        do = [doEAST doWEST];
    end

    if ~isempty(do)
        fname = fnames(do);
    else
        fname = [];
    end

author: CelsoReyes


text:
function C = cumsummat(N, dom, disc)
%CUMSUMMAT   Indefinite integration matrix.
%   C = CUMSUMMAT(N) returns the NxN indefinite integration matrix associated
%   with the Chebyshev spectral collocation method at second-kind Chebyshev
%   points. By convection, the arbitrary constant is chosen so that the result
%   is zero at -1. See CHEBCOLLOC2.CUMSUMMAT for further details.
%
%   D = CUMSUMMAT(N, DOM) scales the indefinite integration matrix D to the
%   domain DOM. DOM should be a 1x2 vector.
%
%   D = CUMSUMMAT(N, DOM, DISC) or CUMSUMMAT(N, DISC) returns the indefinite
%   integration matrix associated with the OPDISCRETIZATION DISC.
%
% See also CUMSUM, CHEBCOLLOC2.DIFFMAT, DIFFMAT.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%% Parse the inputs:
if ( (nargin == 2) )
    if ( isa(dom, 'function_handle') || ischar(dom) )
        disc = dom;
        dom = cheboppref().domain;
    else
        disc = chebcolloc2();
    end
elseif ( nargin == 1 )
    disc = chebcolloc2();
    dom = cheboppref().domain;
end
% Ensure DISC is a discretization:
if ( ischar(disc) )
    disc = str2func(disc);
end
if ( isa(disc, 'function_handle') )
    disc = disc();
end
% No breakpoints allowed:
if ( numel(dom) > 2 )
    dom = dom([1 end]);
    warning('CHEBFUN:cumsummat:noBreaks', ...
        'CUMSUMMAT does not support domains with breakpoints.');
end

%% Call DISC.DIFFMAT(N) and scale appropriately.
scl = .5*(dom(end) - dom(1));
C = scl*disc.cumsummat(N);

end


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = pcev(t,f,f0,sigma, beta)
%Probability Distribution from CEV dynamics
%   dS(t) = sigma S(t)^beta dW(t)
%   t = maturity
%   f0 = spot
%   f = forward
%   sigma = volatility
%   beta = CEV exponent

nu = sigma^2*t;
k = 2 ./(nu * (2-beta)^2);
x = k.*f0^(2-beta);
w = k.*f.^(2-beta);

q = 1 / ( 2 - beta);

wbess = (2-beta) .* (x .* w.^(1-2.*beta)).^(1/(4-2*beta)) ...
    .* exp(-x-w) .*k .^(1 ./(2-beta));

Z = 2 * sqrt(x.*w);

bess = besseli(q,Z);

y= wbess.* bess;

end

author: Sable

repo: mcbench-benchmarks


text:
function timer_tictoc_test04 ( )

%*****************************************************************************80
%
%% TEST04 times the 2D nearest neighbor problem.
%
%  Discussion:
%
%    For the MATLAB implementation of this test, the unvectorized
%    test limits had to be DRASTICALLY reduced.  But then I noticed
%    that some kind of MATLAB cleverness causes all the loops
%    after the first cycle to be computed in almost no time.
%    This test will have to be rethought.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2006
%
%  Author:
%
%    John Burkardt
%
  n_log_min = 10;
  n_log_max = 18;
  n_min = 2^n_log_min;
  n_max = 2^n_log_max;
  n_rep = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  Time the 2D nearest neighbor problem.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Given X(2,N) and Y(2),\n' );
  fprintf ( 1, '    find X(2,*) closest to Y(2).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    for i = 1 : n\n' );
  fprintf ( 1, '      if distance ( x(2,i), y ) < minimum so far\n' );
  fprintf ( 1, '        x_min = x(2,i)\n' );
  fprintf ( 1, '    end\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data vectors will be of minimum size %d\n', n_min );
  fprintf ( 1, '  Data vectors will be of maximum size %d\n', n_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', n_rep );

  x = rand (2,n_max);
  y = rand(1,2);

  for i_rep = 1 : n_rep

    for n_log = n_log_min : n_log_max

      n = 2^n_log;

      tic;

      dist_min = r8_huge ( );
      i_min = 0;
      for i = 1 : n
        dist_i = sum ( ( x(1:2,i) - y(1:2)' ).^2 );
        if ( dist_i < dist_min )
          dist_min = dist_i;
          i_min = i;
        end
      end

      delta(n_log,i_rep) = toc;

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04 Results:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Vector Size  Rep #1        Rep #2        Rep #3        ' );
  fprintf ( 1, 'Rep #4        Rep #5\n' );
  fprintf ( 1, '\n' );
  for n_log = n_log_min : n_log_max
    n = 2^n_log;
    fprintf ( 1, '%10d', n );
    for j = 1 : n_rep
      fprintf ( 1, '%14f', delta(n_log,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
% Test gmux.
% The following model, where Y is a gmux node,
% and M is set to 1, should be equivalent to X1 -> Y
%
% X1 Xn M
% \ |  /
%   Y

n = 3;
N = n+2;
Xs = 1:n;
M = n+1;
Y = n+2;
dag = zeros(N,N);
dag([Xs M], Y)=1;

dnodes = M;
ns = zeros(1, N);
sz = 2;
ns(Xs) = sz;
ns(M) = n;
ns(Y) = sz;

bnet = mk_bnet(dag, ns, 'discrete', M, 'observed', [M Y]);

psz = ns(Xs(1));
selfsz = ns(Y);

W = randn(selfsz, psz);
mu = randn(selfsz, 1);
Sigma = eye(selfsz, selfsz);

bnet.CPD{M} = root_CPD(bnet, M);
for i=Xs(:)'
  bnet.CPD{i} = gaussian_CPD(bnet, i, 'mean', zeros(psz, 1), 'cov', eye(psz, psz));
end
bnet.CPD{Y} = gmux_CPD(bnet, Y, 'mean', mu, 'weights', W, 'cov', Sigma);

evidence = cell(1,N);
yval = randn(selfsz, 1);
evidence{Y} = yval;
m = 2;
%notm = not(m-1)+1; % only valid for n=2
notm = mysetdiff(1:n, m);
evidence{M} = m;

engines = {};
engines{end+1} = jtree_inf_engine(bnet);
engines{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');

for e=1:length(engines)
  engines{e} = enter_evidence(engines{e}, evidence);
  mXm{e} = marginal_nodes(engines{e}, Xs(m));

  % Since M=m, only Xm was updated.
  % Hence the posterior on Xnotm should equal the prior.
  for i=notm(:)'
    mXnotm = marginal_nodes(engines{e}, Xs(i));
    assert(approxeq(mXnotm.mu, zeros(psz,1)))
    assert(approxeq(mXnotm.Sigma, eye(psz, psz)))
  end
end

% Check that all engines give the same posterior
for e=2:length(engines)
  assert(approxeq(mXm{e}.mu, mXm{1}.mu))
  assert(approxeq(mXm{e}.Sigma, mXm{1}.Sigma))
end

% Compute the correct posterior by building Xm -> Y

N = 2;
dag = zeros(N,N);
dag(1, 2)=1;
ns = [psz selfsz];
bnet = mk_bnet(dag, ns, 'discrete', [], 'observed', 2);

bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', zeros(psz, 1), 'cov', eye(psz, psz));
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', mu, 'cov', Sigma, 'weights', W);

jengine  = jtree_inf_engine(bnet);
evidence = {[], yval};
jengine = enter_evidence(jengine, evidence); % apply Bayes rule to invert the arc
mX = marginal_nodes(jengine, 1);

for e=1:length(engines)
  assert(approxeq(mX.mu, mXm{e}.mu))
  assert(approxeq(mX.Sigma, mXm{e}.Sigma))
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function r8vec_sorted_split_test ( )

%*****************************************************************************80
%
%% R8VEC_SORTED_SPLIT_TEST tests R8VEC_SORTED_SPLIT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 25;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_SORTED_SPLIT_TEST\n' );
  fprintf ( 1, '  R8VEC_SORTED_SPLIT splits a sorted vector into\n' );
  fprintf ( 1, '  entries less than and greater than a\n' );
  fprintf ( 1, '  splitting value.\n' );

  b = 0.0;
  c = 10.0;
  seed = 123456789;

  [ a, seed ] = r8vec_uniform_ab ( n, b, c, seed );

  a(1:n) = round ( a(1:n) ) / 2.0;

  a = r8vec_sort_heap_a ( n, a );

  split = 0.5 * ( a(1) + a(n) );

  r8vec_print ( n, a, '  The sorted array:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Splitting value is %f\n', split );

  [ i_lt, i_gt ] = r8vec_sorted_split ( n, a, split );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Lower index I_LT = %d\n', i_lt );
  fprintf ( 1, '  Upper index I_GT = %d\n', i_gt );

  return
end


text:
function cloudInterface( src, event )
% interface function for several animation display functions
%
% USAGE
%  cloudInterface( src, event )
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also CLOUDBOUNDARY
%
% Vincent's Structure From Motion Toolbox      Version NEW
% Copyright (C) 2008-2010 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

userData=get(gcf,'UserData');

anim = userData.anim; cam = userData.cam; hCam = userData.hCam;
hPoint = userData.hPoint; hConn = userData.hConn; hGT = userData.hGT;

if userData.isIndep
  animRange = 1;
else
  animRange = 1:numel(anim);
end

if ~iscell(hGT); hGT = {hGT}; end

% Deal with a pressed key to change the view or quit the animation
diffView=[0 0];
if ~isfield(event,'Key') || ~any(strcmp(event.Key,userData.possibleKey))
  event.Key='';
end
switch event.Key
  case 'c', % Show the connectivity
    userData.showConn = ~userData.showConn;
  case 'f', % Show the first point
    userData.showFirst = ~userData.showFirst;
  case 'g', % Show the ground truth
    userData.showGT = ~userData.showGT;
  case 'i', % make a view independent from the other ones
    userData.isIndep=~userData.isIndep;
  case 'p', % toggle between pretty axes
    userData.showPrettyAxes=~userData.showPrettyAxes;
  case 'q',
    userData.doReturn=1;
  case 'r', % rotate the axes
    for i=animRange
      if ~isempty(anim{i}.S)
        anim{i}.S = anim{i}.S([2 3 1],:,:);
        anim{i}.R = anim{i}.R(:,[2 3 1],:);
        anim{i}.t = anim{i}.t([2 3 1],:);
        for j=1:2
          tmp = [ get(hPoint{i}(j),'XData'); get(hPoint{i}(j),'YData'); ...
            get(hPoint{i}(j),'ZData') ];
          set(hPoint{i}(j),'XData',tmp(2,:));
          set(hPoint{i}(j),'YData',tmp(3,:));
          set(hPoint{i}(j),'ZData',tmp(1,:));
          tmp = [ get(hGT{i}(j),'XData'); get(hGT{i}(j),'YData'); ...
            get(hGT{i}(j),'ZData') ];
          set(hGT{i}(j),'XData',tmp(2,:));
          set(hGT{i}(j),'YData',tmp(3,:));
          set(hGT{i}(j),'ZData',tmp(1,:));
        end
      end
    end
  case 't', % Show the title
    userData.showTitle = ~userData.showTitle;
  case {'1','2','3'} % switch 2D/3D view
    if userData.camMode==str2double(event.Key) && src>=0 && ...
        any(strcmp( 'g', userData.possibleKey ))
      userData.alignGT = ~userData.alignGT;
    else
      for i=animRange; cam{i}(userData.camMode).axis = axis; end
      userData.camMode = str2double(event.Key);
      for i=animRange
        switch userData.camMode
          case 1,
            axis( cam{i}(userData.camMode).gca, 'on', 'vis3d' );
            set(gcf,'Color',0.8*[1 1 1]);
            if ~isempty(hCam{i})
              set(reshape(hCam{i}(:,1:8),[],1),'Visible','on');
            end
          case 2,
            axis( cam{i}(userData.camMode).gca, 'off', 'equal' );
            set(gcf,'Color',[1 1 1]);
            if ~isempty(hCam{i})
              set(reshape(hCam{i}(:,1:8),[],1),'Visible','off');
            end
          case 3,
            axis( cam{i}(userData.camMode).gca, 'on', 'vis3d' );
            set(gcf,'Color',[1 1 1]);
            if ~isempty(hCam{i})
              set(reshape(hCam{i}(:,1:8),[],1),'Visible','on');
            end
        end
      end

      for i=animRange
        axis( cam{i}(userData.camMode).gca, cam{i}(userData.camMode).axis);
        if userData.camMode~=2
          set( cam{i}(userData.camMode).gca, 'CameraPositionMode', 'auto');
        else
          set( cam{i}(userData.camMode).gca, 'CameraPosition', [0 -5 0] );
        end
        set( cam{i}(userData.camMode).gca, 'CameraTarget', ...
          cam{i}(userData.camMode).target );
      end
    end
  case 'numpad0' % re-center the camera
    for i=animRange; cam{i}(userData.camMode).diff=0; end
  case 'subtract', % slow down the speed
    userData.fps = userData.fps-5;
    if userData.fps<0; userData.fps=1; end
  case 'add', % increase the speed
    userData.fps = userData.fps+5;
  case 'space', % pause
    userData.doPause = ~userData.doPause;
  case 'leftarrow',
    diffView = [ 10 0 ];
  case 'rightarrow',
    diffView = [ -10 0];
  case 'uparrow',
    diffView = [ 0 -10];
  case 'downarrow',
    diffView = [ 0 10 ];
end

% Update the look of the data
offon = {'off' 'on'};
for i=animRange
  cam{i}(userData.camMode).diff = cam{i}(userData.camMode).diff + diffView;

  % Set some stuff visible or not
  set(hPoint{i}(2),'Visible',offon{1+userData.showFirst});
  set(get(cam{i}(userData.camMode).gca,'Title'),'Visible',...
    offon{1+userData.showTitle});

  if (~isempty(hGT)) && (~isempty(hGT{i}))
    if hGT{i}(1)>=0
      set(hGT{i}(1),'Visible',offon{1+userData.showGT});
    end
    if (length(hGT{i}) >=2 ) && (hGT{i}(2)>=0)
      set(hGT{i}(2),'Visible',offon{1+userData.showGT*userData.showFirst});
    end
  end
  if (~isempty(hConn))
    set(hConn{i},'Visible',offon{1+userData.showConn});
  end
end

if userData.camMode~=2 || isempty(event.Key)
  for i=animRange
    set( cam{i}(userData.camMode).gca, 'View', ...
      cam{i}(userData.camMode).view + cam{i}(userData.camMode).diff);

    if userData.showPrettyAxes
      set(cam{i}(userData.camMode).gca,'XTick',[],'YTick',[],'ZTick',[],...
        'LineWidth',3);
    else
      set(cam{i}(userData.camMode).gca,'XTickMode','auto','YTickMode',...
        'auto','ZTickMode','auto','LineWidth',1);
    end
  end
end

% Save the data back
userData.anim = anim; userData.cam = cam; userData.hPoint = hPoint;
userData.hGT = hGT;

set(gcf,'UserData',userData);

author: vrabaud


text:
function origin = get_origin(this)
% Computes the origin of this image geometry.
%
%   Y = MrImageGeometry()
%   Y.get_origin()
%
% This is a method of class MrImageGeometry.
%
% OUT   Voxel indices x, y, z of the origin, i.e. those voxel who are
%       at location [0 0 0].
%
% EXAMPLE
%   origin = Y.get_origin()
%
%   See also MrImageGeometry demo/MrImageGeometry/definition_of_geometry

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-11-06
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

invA = inv(this.get_affine_matrix());
origin = invA(1:3,4);
end


text:
function oHDEVArray = hDevCalculation(fracFreq, sPeriod, tau)
%function for doing ovelapping Hadamard calculations for an array tau
%values. input arguments include array of fractional frequency readings,
%sPeriod of readings for overlap factor, and array of tau values

 phaseError = calculatePhaseError((1/sPeriod),fracFreq); %turn freq readings to frac freq values

 tCount = numel(tau); %get the size of the tau array
 oHDEVArray = zeros(1,tCount); %allocate array for Allan Dev values

 %add a wait bar so user knows status
j = 1/tCount;
h = waitbar(0,'Performing HDEV calculations...','CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
 setappdata(h,'canceling',0)
 %loop through tau values and calculate HDEV at each value
 for i = 1:tCount
     oHDEVArray(i) = calculateHDEV(tau(i),sPeriod,phaseError);
     if getappdata(h,'canceling')
        oHDEVArray = [];
        oHDEVArray = -42881;
        delete(h);
        return;
    end
     waitbar((i*j),h,'Performing HDEV calculations...');
 end

 waitbar(1.0,h,'HDEV calculations done');
 delete(h);


text:
function result = intersect2(cell)

% PURPOSE:
%   This function finds the common elements of several 1D numeric arrays
%   of the same or different sizes and returns an array consisting of only
%   those elements.
%
% INPUT:
%   The input variable "cell" has to be a cell array, with each cell
%   occupied by an 1D numeric array. For example, if you want to find the
%   intersection (common elements) of the following three arrays
%       a = [ 1 3 4 6 8 9 ];
%       b = [ 3 1 0 8 6 4 ];
%       c = [ 7 8 1 9 3 4 ];,
%   you need to first put all of them into a cell array, i.e.
%       cell = {a, b, c};
%   Then you could use "cell" as the input variable to this function, i.e.
%       result = intersect2(cell);
%
% OUTPUT:
%   The output of this function is simply an array consisting of elements
%   that are common to all the arrays stored in "cell". For the particular
%   case in the example above, the output would be
%       result = [ 1 3 4 8 ];
%   Also note that the numbers in the result are sorted in ascending order.

n = max(size(cell));

if n == 2
    result = intersect(cell{1, 1}, cell{1, 2});
elseif (n > 2 && rem(n,2) == 0)
    for i = 1:(n/2)
        intersections{1, i} = intersect(cell{1, i}, cell{1, n-i+1});
    end
elseif (n > 2 && rem(n,2) == 1)
    for i = 1:((n-1)/2)
        intersections{1, i} = intersect(cell{1, i}, cell{1, n-i+1});
    end
        intersections{1, i+1} = cell{1, (n+1)/2};
else
    error('You need at least two arrays as inputs');
end

if n > 2
    result = intersect2(intersections);
end


text:
function [label, model, llh] = mixLinReg(X, y, k, lambda)
% Mixture of linear regression
% input:
%   X: d x n data matrix
%   y: 1 x n responding vector
%   k: number of mixture component
%   lambda: regularization parameter
% output:
%   label: 1 x n cluster label
%   model: trained model structure
%   llh: loglikelihood
% Written by Mo Chen (sth4nth@gmail.com).
if nargin < 4
    lambda = 1;
end
n = size(X,2);
X = [X;ones(1,n)]; % adding the bias term
d = size(X,1);
label = ceil(k*rand(1,n));  % random initialization
R = full(sparse(label,1:n,1,k,n,n));
tol = 1e-6;
maxiter = 500;
llh = -inf(1,maxiter);
Lambda = lambda*eye(d);
W = zeros(d,k);
Xy = bsxfun(@times,X,y);
beta = 1;
for iter = 2:maxiter
    % maximization
    nk = sum(R,2);
    alpha = nk/n;
    for j = 1:k
        Xw = bsxfun(@times,X,sqrt(R(j,:)));
        U = chol(Xw*Xw'+Lambda);
        W(:,j) = U\(U'\(Xy*R(j,:)'));  % 3.15 & 3.28
    end
    D = bsxfun(@minus,W'*X,y).^2;
    % expectation
    logRho = (-0.5)*beta*D;
    logRho = bsxfun(@plus,logRho,log(alpha));
    T = logsumexp(logRho,1);
    logR = bsxfun(@minus,logRho,T);
    R = exp(logR);
    llh(iter) = sum(T)/n;
    if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter)); break; end
end
llh = llh(2:iter);
model.alpha = alpha; % mixing coefficient
model.beta = beta; % mixture component precision
model.W = W;  % linear model coefficent
[~,label] = max(R,[],1);
model.label = label;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
% THis function will search...
% INputs:
% lines: cell array of entries to be searched on
% pattern: query
% start:

function found = searchInCellBinary(lines, pattern, start, stop)

found = 0;

if stop-start <= 10
    for i=start:stop
        if strcmp(lines{i}, pattern)
            found = i;
            return;
        end
    end

else

    middle = floor( (start+stop)/2 );
    notEqual = 0;
    middleWord = lines{middle};
    min_length = min( length(pattern), length(middleWord) );
    for i=1:min_length
        if int16(pattern(i)) > int16(middleWord(i))
            notEqual = 1;
            found = searchInCellBinary(lines, pattern, middle, stop);
            break;
        elseif int16(pattern(i)) < int16(middleWord(i))
            notEqual = 1;
            found = searchInCellBinary(lines, pattern, start, middle);
            break;
        end
    end

    if notEqual == 0
        if length(pattern) == length(middleWord)
            found = middle;
        elseif length(pattern) > length(middleWord)
            found = searchInCellBinary(lines, pattern, middle+1, stop);
        elseif length(pattern) < length(middleWord)
            found = searchInCellBinary(lines, pattern, start, middle-1);
        end
    end
end

author: singaxiong

repo: SignalGraph


text:
% rovXYZ 3 x N vectors, inillh 3 x N vectors
% output 3 x N vectors
function dNED=posdiff_v001(rovXYZ, inillh)
% for blh2xyz, input N x 3, output N x 3;
startXYZ= blh2xyz(inillh');
startCen=llh2dcm_v000(inillh(1:2),[0,1]);
dXYZ = rovXYZ-startXYZ';
dNED = startCen*dXYZ;

repo: ekfmonoslam


text:
function A = adddiag(A,v,k)
%ADDDIAG  Modify the diagonal(s) of a matrix.
% ADDDIAG(A,V)
% ADDDIAG(A,V,K)
%
% See also DIAG, FINDDIAG.

if nargin < 3
  k = 0;
end
i = finddiag(A,k);
A(i) = A(i) + v;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function f = p01_f ( x )

%*****************************************************************************80
%
%% P01_F evaluates the objective function for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = ( x - 2.0 ) * ( x - 2.0 ) + 1.0;

  return
end


text:
% Copyright (C) 2016, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function jacobian_symbolic_spherical_v
% link lengths
L3 = 0.38;
L6 = 0.065;

syms q1 q2 q3
% matrices DH
A01 = dh_sym(q1, L3, 0, pi/2);
A12 = dh_sym(q2, 0, 0, -pi/2);
A23 = dh_sym(q3, L6, 0, 0);

A02 = A01*A12;
A03 = A01*A12*A23;

z0 = [0 0 1]';
z1 = A01(1:3,3);
z2 = A02(1:3,3);

p03=A03(1:3,4);
p13=A03(1:3,4)-A01(1:3,4);
p23=A03(1:3,4)-A02(1:3,4);

Jv = [cross(z0, p03) cross(z1, p13) cross(z2, p23)];

singularities = det(Jv)

function A = dh_sym(theta, d, a, alpha)
syms q1 q2 q3
% avoid almost zero elements in cos(alpha) and sin(alpha)
ca = cos(alpha);
sa = sin(alpha);
if abs(ca) < 1e-6
    ca = 0;
end
if abs(sa) < 1e-6
    sa = 0;
end

A=[cos(theta)  -ca*sin(theta)   sa*sin(theta)   a*cos(theta);
   sin(theta)   ca*cos(theta)  -sa*cos(theta)   a*sin(theta);
            0              sa             ca             d;
            0         0                     0              1];

repo: ARTE


text:
function n = p33_n ( )

%*****************************************************************************80
%
%% P33_N returns the number of variables for problem 33.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 4;

  return
end


text:
function r = or(p,q)
% this function compares the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    r=p;
    r.xvalue = p.xvalue|q.xvalue;
    r.yvalue = p.yvalue|q.yvalue;
    r.zvalue = p.zvalue|q.zvalue;
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue|q;
    r.yvalue = p.yvalue|q;
    r.zvalue = p.zvalue|q;
else
    r=q;
    r.xvalue = p|q.xvalue;
    r.yvalue = p|q.yvalue;
    r.zvalue = p|q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [k, sk, n2] = rbfperiodic2KernCompute(kern, x, x2)

% RBFPERIODIC2KERNCOMPUTE Compute the RBFPERIODIC2 kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the RBF periodic covariance with variying period
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the RBF periodic covariance with variying period
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : rbfperiodic2KernParamInit, kernCompute, kernCreate, rbfperiodic2KernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2009
%
% MODIFICATIONS : Andreas C. Damianou, 2011
%
% MODIFICATIONS : Michalis K. Titsias, 2011

% KERN

factor = kern.factor; % Default (if period is fixed: 2*pi/kern.period)
if nargin < 3
  n2 = sin(0.5*factor*(repmat(x, 1, size(x, 1)) - repmat(x', size(x, 1), 1)));
  n2 = n2.*n2;
  wi2 = (2 .* kern.inverseWidth);
  sk = exp(-n2*wi2);
else
  n2 = sin(0.5*factor*(repmat(x, 1, size(x2, 1)) - repmat(x2', size(x, 1), 1)));
  n2 = n2.*n2;
  wi2 = (2 .* kern.inverseWidth);
  sk = exp(-n2*wi2);
end
k = kern.variance*sk;
% Test kernel with: kernTest('rbfperiodic2',1)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [O]=titgen(S)
%these function give a string for example from poly2strs.m and if power of
%if each element >=10 then place it between { }
%for example if input: s^11+s^10+s^9+1  then output is s^{11}+s^{10}+s^9+1
%these is used for title.m for power >=10

if isempty(S)
    display('Error: Input String is empty');
    return;
end

len=length(S);
pointer=0;
k=1;
for i=1:len
    if S(i)=='^'
        pointer(1,k)=i;
        k=k+1;
    end
    if (S(i)=='+' || S(i)=='-') && i~=1
        pointer(1,k)=i;
        k=k+1;
    end
end
lp=length(pointer);
str{1}=S(1:pointer(1));
for j=1:lp
    if j~=lp
        str{j+1}=S(pointer(j)+1:pointer(j+1));
    end
end
str{j+1}=S(pointer(lp)+1:len);

for m=2:2:length(str)
    temp=str{m};
    if temp(length(temp))=='+'
        temp(length(temp))='';
        str{m}=temp;
        str{m+1}=strcat('+',str{m+1});
    elseif temp(length(temp))=='-'
        temp(length(temp))='';
        str{m}=temp;
        str{m+1}=strcat('-',str{m+1});
    end

    num_pow=str2double(str{m});
    if num_pow>=10
        str{m}=sprintf('{%g}',num_pow);
    end
end

O=strcat(str{:});


text:
function [ know, x ] = p05_sol ( m, know )

%*****************************************************************************80
%
%% P05_SOL returns the solution for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution, if known.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = 0.0;
  else
    know = 0;
    x(1:m,1) = 0.0;
  end

  return
end


text:
%% gpatch
% Below is a demonstration of the features of the |gpatch| function

%%
clear; close all; clc;

%% Syntax
% |[hp]=gpatch(F,V,C,CE,A,L)|;

%% Description
% This function is a short-hand version of the |patch| command. The inputs
% for |gpatch| are the faces (F), the vertices (V), the color description
% (C), the edge color description CE, the transparancy (A), and the edge
% width (L).
% The color data descriptions C (or equivalently CE for edges) can be:
% 1) A string such as 'g' for green
% 2) A triplet of RGD values e.g. [1 0 0] is blue
% 3) A nx1 or a mx1 array of colormapped colors (where n=size(F,1) or m=size(V,1))
% 4) (simiarl to 3) A nx3 or a mx3 RGB color value array for the faces or vertices respectively.

%% Examples

%%
% Create example mesh data

[F,V,~]=geoSphere(2,1); %Faces and vertices
CV=V(:,1); %Color information for vertices
CF=vertexToFaceMeasure(F,CV); %Color information for faces
CF_rgb=abs(vertexToFaceMeasure(F,V)); %Color information for faces

%% Example: Introduction to using |gpatch| for mesh visualization
% The below visualization show the syntax require using |patch| and
% |gpatch|. Essentially |gpatch| is just a short-hand version of |patch|
% allowing for quick and easy visualization using patch graphics.

%%
% Using patch graphics in MATLAB see documentation on |patch| for more information

cFigure;
subplot(2,3,1);
title('Single face and edge color');
patch('Faces',F,'Vertices',V,'FaceColor','r','EdgeColor','g','FaceAlpha',0.5);
axisGeom;

subplot(2,3,2);
title('Colormapped face colors');
patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',CF,'EdgeColor','k','FaceAlpha',1);
colormap gjet; colorbar;
axisGeom;

subplot(2,3,3);
title('RGB driven face colors');
patch('Faces',F,'Vertices',V,'FaceColor','flat','FaceVertexCData',CF_rgb,'EdgeColor','k','FaceAlpha',1,'LineWidth',3);
axisGeom;

subplot(2,3,4);
title('Colormapped edge colors');
patch('Faces',F,'Vertices',V,'FaceColor',grayColor(0.5),'EdgeColor','flat','CData',CV,'LineWidth',3);
colormap gjet; colorbar;
axisGeom;

subplot(2,3,5);
title('alpha mapping');
patch('Faces',F,'Vertices',V,'FaceColor','g','FaceVertexAlphaData',CV,'EdgeColor','none','FaceAlpha','flat');
axisGeom;

drawnow;

% Using |gpatch| shorthand alternative to |patch|

cFigure;
subplot(2,3,1);
title('Single face and edge color');
gpatch(F,V,'r','g',0.5);
axisGeom;

subplot(2,3,2);
title('Colormapped face colors');
gpatch(F,V,CF);
colormap gjet; colorbar;
axisGeom;

subplot(2,3,3);
title('RGB driven face colors');
gpatch(F,V,CF_rgb,'k',1,3);
axisGeom;

subplot(2,3,4);
title('Colormapped edge colors');
gpatch(F,V,'kw',CV,1,3); %kw -> grayColor(0.5)
colormap gjet; colorbar;
axisGeom;

subplot(2,3,5);
title('alpha mapping');
gpatch(F,V,'g','none',CV);
axisGeom;

drawnow;

%% Using |gpatch| with cell arrays containing patch data

%%
% Use of |gpatch| when both the faces and vertices are stored in a cell
% array

% Create example cell arrays
[F1,V1]=graphicsModels(1); V1=V1-mean(V1,1); V1=V1./max(V1(:)); V1=V1-1;
[F2,V2]=graphicsModels(2); V2=V2-mean(V2,1); V2=V2./max(V2(:));
[F3,V3]=graphicsModels(3); V3=V3-mean(V3,1); V3=V3./max(V3(:)); V3=V3+1;
F={F1,F2,F3}; %Cell array containing faces
V={V1,V2,V3}; %Cell array containing vertices
C={ones(size(F1,1),1),2*ones(size(F2,1),1),3*ones(size(F3,1),1)}; %Cell array containing color labels

cFigure;
title('Cell array of patch data');
hp=gpatch(F,V,C,'none',1);
axisGeom; camlight headlight;
colormap spectral; icolorbar;
drawnow;

%%
% Use of |gpatch| when only the faces are stored in a cell array

% Create example cell array

[F1,V1]=graphicsModels(1); V1=V1-mean(V1,1); V1=V1./max(V1(:)); V1=V1-1;
[F2,V2]=graphicsModels(2); V2=V2-mean(V2,1); V2=V2./max(V2(:));
[F3,V3]=graphicsModels(3); V3=V3-mean(V3,1); V3=V3./max(V3(:)); V3=V3+1;
F={F1,F2+size(V1,1),F3+size(V1,1)+size(V2,1)}; %Cell array containing faces
V=[V1;V2;V3]; %Normal array containing all vertices
C={ones(size(F1,1),1),2*ones(size(F2,1),1),3*ones(size(F3,1),1)}; %Cell array containing color labels

cFigure;
title('Cell array of patch data');
hp=gpatch(F,V,C,'none',1);
axisGeom; camlight headlight;
colormap spectral; icolorbar;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%[1995]-"Particle Swarm Optimization"
%[1998]-"A modified particle swarm optimizer"

% (9/12/2020)

function PSO = jParticleSwarmOptimization(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;
c1    = 2;              % cognitive factor
c2    = 2;              % social factor
w     = 0.9;            % inertia weight
Vmax  = (ub - lb) / 2;  % Maximum velocity

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'c1'), c1 = opts.c1; end
if isfield(opts,'c2'), c2 = opts.c2; end
if isfield(opts,'w'), w = opts.w; end
if isfield(opts,'Vmax'), Vmax = opts.Vmax; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
V   = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
  end
end
% Fitness
fit  = zeros(1,N);
fitG = inf;
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  % Gbest update
  if fit(i) < fitG
    Xgb  = X(i,:);
    fitG = fit(i);
  end
end
% PBest
Xpb  = X;
fitP = fit;
% Pre
curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
% Iterations
while t <= max_Iter
  for i = 1:N
    for d = 1:dim
      r1 = rand();
      r2 = rand();
      % Velocity update (2a)
      VB = w * V(i,d) + c1 * r1 * (Xpb(i,d) - X(i,d)) + ...
        c2 * r2 * (Xgb(d) - X(i,d));
      % Velocity limit
      VB(VB > Vmax) = Vmax;  VB(VB < -Vmax) = -Vmax;
      V(i,d) = VB;
      % Position update (2b)
      X(i,d) = X(i,d) + V(i,d);
    end
    % Boundary
    XB = X(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    X(i,:) = XB;
    % Fitness
    fit(i) = fun(feat,label,(X(i,:) > thres),opts);
    % Pbest update
    if fit(i) < fitP(i)
      Xpb(i,:) = X(i,:);
      fitP(i)  = fit(i);
    end
    % Gbest update
    if fitP(i) < fitG
      Xgb  = Xpb(i,:);
      fitG = fitP(i);
    end
  end
  curve(t) = fitG;
  fprintf('\nIteration %d Best (PSO)= %f',t,curve(t))
  t = t + 1;
end
% Select features based on selected index
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
PSO.sf = Sf;
PSO.ff = sFeat;
PSO.nf = length(Sf);
PSO.c  = curve;
PSO.f  = feat;
PSO.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function M = cross2Matrix(x)
% CROSS2MATRIX  Antisymmetric matrix corresponding to a 3-vector
%
% Computes the antisymmetric matrix M corresponding to a 3-vector x such
% that M*y = cross(x,y) for all 3-vectors y.
%
% Input:
%   - x(3,1) : vector
%
% Output:
%   - M(3,3) : antisymmetric matrix
%

M = [0,    -x(3),  x(2); ...
     x(3),   0,   -x(1); ...
    -x(2),  x(1),   0  ];

end

repo: Stereo-Odometry-SOFT


text:
%compute dxcentral_mm

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

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dxcentral_mm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dxcentral_mm{1,i}=(trx(larva).xcentral_mm(2:end)-trx(larva).xcentral_mm(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dxcentral_mm;

author: kristinbranson

repo: JAABA


text:
function cs_demo (do_pause, matrixpath)
%CS_DEMO run all CSparse demos.
%   cs_demo(0) will run all demos without pausing.
%
% Example:
%   cs_demo
% See also: cs_demo1, cs_demo2, cs_demo3

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

help cs_demo
if (nargin < 1)
    do_pause = 1 ;
end
if (nargin < 2)
    matrixpath = [] ;
end

figure (1)
clf

fprintf ('\n\n-------------------------------------------------------\n') ;
help cs_demo1 ;
cs_demo1 (matrixpath) ;

fprintf ('\n\n-------------------------------------------------------\n') ;
help cs_demo2
cs_demo2 (do_pause, matrixpath) ;

fprintf ('\n\n-------------------------------------------------------\n') ;
help cs_demo3
cs_demo3 (do_pause, matrixpath) ;

fprintf ('\n\n-------------------------------------------------------\n') ;
help private/ex_1
ex_1

fprintf ('\n\n-------------------------------------------------------\n') ;
help private/ex2
ex2

fprintf ('\n\n-------------------------------------------------------\n') ;
help private/ex3
ex3

fprintf ('\nAll CSparse demos finished.\n') ;

author: fieldtrip

repo: fieldtrip


text:
function upsilon = lfmjpComputeUpsilonVector(gamma, sigma2, t)

% LFMJPCOMPUTEUPSILONVECTOR Upsilon vector jolt. pos. with t1, t2 limits
% FORMAT
% DESC computes a portion of the LFM kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

sigma = sqrt(sigma2);

upsilon = gamma^2*lfmvpComputeUpsilonVector(gamma, sigma2, t) ...
    + (4/(sqrt(pi)*sigma^3))*exp(-(t.^2)./sigma2).*(t.*(gamma + 2*t/sigma2) - 1);


text:
function [hdr, nameFlag] = read_buffer_offline_header(headerfile)

% function [hdr, nameFlag] = read_buffer_offline_header(headerfile)
%
% This function reads a FCDC buffer header from a binary file or text file
%
% On return, nameFlag has one of the following values:
%   0 = No labels were generated (fMRI etc.)
%   1 = Fake labels were generated
%   2 = Got channel labels from chunk information

% (C) 2010 S. Klanke

type = {
  'char'
  'uint8'
  'uint16'
  'uint32'
  'uint64'
  'int8'
  'int16'
  'int32'
  'int64'
  'single'
  'double'
};

wordsize = {
  1 % 'char'
  1 % 'uint8'
  2 % 'uint16'
  4 % 'uint32'
  8 % 'uint64'
  1 % 'int8'
  2 % 'int16'
  4 % 'int32'
  8 % 'int64'
  4 % 'single'
  8 % 'double'
};

hdr = [];
hdr.nSamples    = 0;
%hdr.nEvents     = 0;
hdr.nSamplesPre = 0;
hdr.nTrials     = 1;
hdr.orig =[];

nameFlag = 0;

txt=[];
finfo=dir([headerfile '.txt']);
if ( isempty(finfo) || finfo.bytes==0 )
  ft_warning(sprintf('Could not open text header file : %s',headerfile));
else
  txt.nSamples=0;
  FA = fopen_or_error([headerfile '.txt'], 'r');
  while ~feof(FA)
    s = fgetl(FA);
    indEq = find(s=='=');
    if numel(indEq)==1
      key = s(1:(indEq-1));
      val = s((indEq+1):end);
      switch key
       case 'fSample'
        txt.Fs = str2num(val);
       case 'nSamples'
        txt.nSamples = str2num(val);
       case 'nEvents'
        %        txt.nEvents = str2num(val);
       case 'nChans'
        txt.nChans = str2num(val);
        txt.label = cell(txt.nChans,1);
       case 'dataType'
        if strcmp(val,'float32')
          txt.orig.data_type = 'single';
        elseif strcmp(val,'float64')
          txt.orig.data_type = 'double';
        else
          txt.orig.data_type = val;
        end
        txt.orig.wordsize = wordsize{strmatch(txt.orig.data_type,type)};
        txt.data_type = strmatch(txt.orig.data_type,type)-1;
        txt.wordsize  = txt.orig.wordsize;
       case 'version'
        txt.orig.version = str2num(val);
       case 'endian'
        if strcmp(val,'big')
          txt.orig.endianness='ieee-be';
        elseif strcmp(val,'little')
          txt.orig.endianness='ieee-le';
        else
          ft_error('Invalid value for key ''endian''.');
        end
      end
      continue;
    end
    indCol = find(s==':');
    if numel(indCol)>=1
      ind = str2num(s(1:(indCol-1)));
      name = s((indCol+1):end);
      if ~isempty(ind) && ind>=1 && ind<=txt.nChans
        txt.label{ind} = name;
        nameFlag = 2;
      else
        ft_error('Invalid channel name definition - broken header file?');
      end
    end
  end

  fclose(FA);
end

bin=[];
finfo=dir(headerfile);
if ( isempty(finfo) || finfo.bytes==0 )
  ft_warning(sprintf('Couldnt open binary header file : %s',headerfile));
  hdr = copyfields(txt, hdr, fieldnames(txt)); % ensure that the predefined stuff still exists
else
  if ( ~isfield(hdr.orig,'endianness') )
    hdr.orig.endianness='native';
  end
  F = fopen_or_error(headerfile,'rb',hdr.orig.endianness);
  hdr.nChans    = double(fread(F,1,'uint32'));
  hdr.nSamples  = double(fread(F,1,'uint32'));
  hdr.nEvents   = double(fread(F,1,'uint32'));
  hdr.Fs        = double(fread(F,1,'single'));
  hdr.data_type = fread(F,1,'uint32');
  hdr.bufsize   = fread(F,1,'uint32');

  if isfield(txt,'nChans') && ~isequal(hdr.nChans,txt.nChans)
    ft_error('Number of channels in binary header does not match ASCII definition');
  end
  if isfield(txt,'nChans') && ~isequal(hdr.nSamples,txt.nSamples)
    ft_warning('Number of samples in binary header does not match ASCII definition');
  end
  if isfield(txt,'nEvents') && ~isequal(hdr.nEvents,txt.nEvents)
    ft_error('Number of events in binary header does not match ASCII definition');
  end
  if isempty(hdr.data_type)
    hdr.data_type = strmatch(txt.orig.data_type,type)-1;
  end
  %if strcmp(hdr.data_type, type{hdr.data_type+1})
  hdr.orig.data_type = type{hdr.data_type+1};
  hdr.orig.wordsize = wordsize{hdr.data_type+1};
  %else
  %  ft_error('Data type in binary header does not match ASCII definition');
  %end

  % TODO: add chunk handling
  while ~feof(F)
    typeAndSize = fread(F, 2, 'uint32');
    if numel(typeAndSize) < 2
      break
    end
    switch typeAndSize(1)
     case 1 % channel names
            % already dealt with, TODO: maybe check consistency with ASCII stuff
      dummy = fread(F, typeAndSize(2), 'uint8=>uint8');
     case 2 % FT_CHUNK_CHANNEL_FLAGS
            % FIXME: ignored for now
      dummy = fread(F, typeAndSize(2), 'uint8=>uint8');
     case 3 % FT_CHUNK_RESOLUTIONS
      if typeAndSize(2) == 8*hdr.nChans
        hdr.resolutions = fread(F, [hdr.nChans, 1], 'double');
      else
        ft_warning('Invalid size of RESOLUTIONS chunk - skipping');
        dummy = fread(F, typeAndSize(2), 'uint8=>uint8');
      end
     case 4 % FT_CHUNK_ASCII_KEYVAL
      dummy = fread(F, typeAndSize(2), 'uint8');
      % FIXME: ignore for now
     case 5 % FT_CHUNK_NIFTI1
      hdr.orig.nifti_1 = fread(F, [1, typeAndSize(2)], 'uint8=>uint8');
      if typeAndSize(2) == 348
        hdr.nifti_1 = decode_nifti1(hdr.orig.nifti_1);
      else
        ft_warning('Invalid size of NIFTI_1 chunk - skipping');
      end
     case 6 % FT_CHUNK_SIEMENS_AP = 6
      hdr.orig.siemensap = fread(F, typeAndSize(2), 'uint8=>uint8');
      if exist('sap2matlab')==3
        hdr.siemensap = sap2matlab(hdr.orig.siemensap);
      end
     case 7 % FT_CHUNK_CTF_RES4 = 7
      hdr.orig.ctf_res4 = fread(F, typeAndSize(2), 'uint8=>uint8');
      tmp_name = tempname;
      FT = fopen_or_error(tmp_name, 'wb');
      fwrite(FT, orig.ctf_res4, 'uint8');
      fclose(FT);
      R4F = read_ctf_res4(tmp_name);
      delete(tmp_name);
      % copy over the labels
      hdr.label = R4F.label;
      % copy over the 'original' header
      hdr.orig = R4F;
    end
  end

  fclose(F);
end

% binary has invalid samples info (or isn't there), get from samples file instead..
if ~isfield(hdr,'nSamples') || hdr.nSamples == 0
  datafile = [headerfile(1:end-6) 'samples'];
  FD = fopen_or_error(datafile, 'rb');
  fseek(FD,0,'eof');
  fileSize = ftell(FD);
  fclose(FD);

  n = fileSize / hdr.orig.wordsize / hdr.nChans;
  hdr.nSamples = floor(n);
  if hdr.nSamples ~= n
    ft_warning('Size of ''samples'' is not a multiple of the size of one sample');
  end
end

if nameFlag < 2 && hdr.nChans < 2000
  nameFlag = 1; % fake labels generated - these are unique
  for i=1:hdr.nChans
    hdr.label{i} = sprintf('%d', i);
  end
end