text:
function peakmask = getpeaks(w)
   %GETPEAKS return mask for peak values for a waveform
   %     PEAKMASK = GETPEAKS(WAVEFORM)
   %     returns a mask of all maximums in the waveform.
   %
   %   Input Arguments
   %       WAVEFORM: waveform object  (SINGLE WAVEFORM)
   %
   %   Output
   %       PEAKMASK: a mask of all local maximums within the waveform.
   %       This mask can then be used to reference the peak points in a
   %       waveform
   %
   %   a MASK is a logical array the same size as the waveform's data,
   %   where the value is TRUE when that data point is a peak, and is zero
   %   otherwise.
   %
   %   in some cases, it may be preferable to get the peakmask for the
   %   absolute value of the data, that way both highs and lows are marked.
   %
   % See also WAVEFORM/HILBERT

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   if ~isscalar(w)
      error('Waveform:getpeaks:tooManyWaveforms',...
         'getpeaks can only be used with individual waveforms');
   end
   BiggerLeft = [(w.data(1:end-1) >= w.data(2:end)); true];
   BiggerRight = [true ; (w.data(1:end-1) < w.data(2:end))];
   peakmask = (BiggerLeft & BiggerRight);
end

author: geoscience-community-codes

repo: GISMO


text:
function x= rectwin(N)
% Function RECTWIN(N) creates the N-point Rectangular window.

x(1 : N)=ones(1, N);

author: Sable

repo: mcbench-benchmarks


text:
function out = proj_simplex(x,r,eq_flag)
%PROJ_SIMPLEX computes the orthogonal projection onto the r-simplex {x:sum(x)=r,x>=0}
%                                                or r-full simplex {x:sum(x)<=r,x>=0}
%
%  Usage:
%  out = PROJ_SIMPLEX(x,[r],[eq_flag])
%  =============================================================
%  Input:
%  x - point to be projected (vector/matrix)
%  r - a positive scalar [default: 1]
%  eq_flag - a flag that determines whether the projection is onto the r-simplex ('eq') or
%                 the r-full simplex ('ineq') [defualt: 'eq']
%  ==============================================================
%  Assumptions:
%  r > 0
%  ==============================================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

%reading the user x and setting defalut values when reqired.
if (nargin < 1)
    error ('usage: proj_simplex(x,r,[eq_flag])') ;
end

if (( nargin < 3) || (isempty(eq_flag)) )
    %eq_flag is not given, setting to default value: true
    eq_flag = 'eq' ;
end

if ((nargin < 2) || (isempty(r)))
    %r is not given, setting to default value: 1
    r = 1 ;
end

if (r < 0)
    error('Set is infeasible') ;
end

if (strcmp(eq_flag,'eq') == 1)
    %call proj_hyperplane_box
   out = proj_hyperplane_box(x,ones(size(x)),r,0) ;
else
    if (strcmp(eq_flag,'ineq') == 1)
        %call proj_halfspace_box
        out = proj_halfspace_box(x,ones(size(x)),r,0) ;
    else
        error('usage: proj_simplex(x,r,[eq_flag]) - eq_flag should be either eq or ineq') ;
    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
%P = vgg_P_from_F(F)  Compute cameras from fundamental matrix.
%   F has size (3,3), P has size (3,4).
%
%   If x2'*F*x1 = 0 for any pair of image points x1 and x2,
%   then the camera matrices of the image pair are
%   P1 = eye(3,4) and P2 = vgg_P_from_F(F), up to a scene homography.

% Tomas Werner, Oct 2001

function P = vgg_P_from_F(F)

[U,S,V] = svd(F);
e = U(:,3);
P = [-vgg_contreps(e)*F e];

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function hhh=tapas_uniqc_hline(y,in1,in2)
% function h=hline(y, linetype, label)
%
% Draws a horizontal line on the current axes at the location specified by 'y'.  Optional arguments are
% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line.  The
% label appears in the same color as the line.
%
% The line is held on the current axes, and after plotting the line, the function returns the axes to
% its prior hold state.
%
% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
% legends, but it is not findable by using findobj.  Specifying an output argument causes the function to
% return a handle to the line, so it can be manipulated or deleted.  Also, the HandleVisibility can be
% overridden by setting the root's ShowHiddenHandles property to on.
%
% h = tapas_uniqc_hline(42,'g','The Answer')
%
% returns a handle to a green horizontal line on the current axes at y=42, and creates a text object on
% the current axes, close to the line, which reads "The Answer".
%
% tapas_uniqc_hline also supports vector inputs to draw multiple lines at once.  For example,
%
% tapas_uniqc_hline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
%
% draws three lines with the appropriate labels and colors.
%
% By Brandon Kuczenski for Kensington Labs.
% brandon_kuczenski@kensingtonlabs.com
% 8 November 2001

if length(y)>1  % vector input
    for I=1:length(y)
        switch nargin
        case 1
            linetype='r:';
            label='';
        case 2
            if ~iscell(in1)
                in1={in1};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            label='';
        case 3
            if ~iscell(in1)
                in1={in1};
            end
            if ~iscell(in2)
                in2={in2};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            if I>length(in2)
                label=in2{end};
            else
                label=in2{I};
            end
        end
        h(I)=tapas_uniqc_hline(y(I),linetype,label);
    end
else
    switch nargin
    case 1
        linetype='r:';
        label='';
    case 2
        linetype=in1;
        label='';
    case 3
        linetype=in1;
        label=in2;
    end

    g=ishold(gca);
    hold on

    x=get(gca,'xlim');
    h=plot(x,[y y],linetype);
    if ~isempty(label)
        yy=get(gca,'ylim');
        yrange=yy(2)-yy(1);
        yunit=(y-yy(1))/yrange;
        if yunit<0.2
            text(x(1)+0.02*(x(2)-x(1)),y+0.02*yrange,label,'color',get(h,'color'))
        else
            text(x(1)+0.02*(x(2)-x(1)),y-0.02*yrange,label,'color',get(h,'color'))
        end
    end

    if g==0
    hold off
    end
    set(h,'tag','hline','handlevisibility','off') % this last part is so that it doesn't show up on legends
end % else

if nargout
    hhh=h;
end

author: translationalneuromodeling

repo: tapas


text:
function cartPoints=ellips2Cart(points,a,f)
%%ELLIPS2CART Convert ellipsoidal coordinates to ECEF Cartesian
%             coordinates.
%
%INPUTS: points One or more points given in geodetic latitude and
%               longitude, in radians, and height, in meters that are to be
%               converted to Cartesian coordinates. To convert N points,
%               points is a 3XN matrix with each column having the format
%               [latitude;longitude; height].
%             a The semi-major axis of the reference ellipsoid. If this
%               argument is omitted or an empty matrix is passed, the value
%               in Constants.WGS84SemiMajorAxis is used.
%             f The flattening factor of the reference ellipsoid. If this
%               argument is omitted or an empty matrix is passed, the value
%               in Constants.WGS84Flattening is used.
%
%OUTPUTS: cartPoints For N points, cartPoints is a 3XN matrix of the
%               converted points with each column having the format
%               [x;y;z].
%
%The conversions are mentioned in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "Simulating aerial targets in 3D accounting for the
%    Earth's curvature," Journal of Advances in Information Fusion, vol.
%    10, no. 1, Jun. 2015.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(f))
    f=Constants.WGS84Flattening;
end

if(nargin<2||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

%The geodetic latitudes
phi=points(1,:);
%The longitudes
lambda=points(2,:);

if(size(points,1)==3)
    %The altitudes
    h=points(3,:);
else
    h=0;
end

sinP=sin(phi);
cosP=cos(phi);
sinL=sin(lambda);
cosL=cos(lambda);

%The square of the first numerical eccentricity
e2=2*f-f^2;
%The normal radii of curvature.
Ne=a./sqrt(1-e2*sinP.^2);

x=(Ne+h).*cosP.*cosL;
y=(Ne+h).*cosP.*sinL;
z=(Ne*(1-e2)+h).*sinP;

cartPoints=[x;y;z];
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function f = sparsnessCol2(Hkt_T, varargin)
k=varargin{1};
t=varargin{2};
Hkt=exp(reshape(Hkt_T, k, t));
sumH_t = sum(Hkt,2);
% figure
% bar(sumH_t)
sumH_t_sq = sum(sumH_t.^2);
sumH = sum(sumH_t);
nh = length(sumH_t);
f = (sqrt(nh)-sumH/(sqrt(sumH_t_sq)))/(sqrt(nh)-1);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% Spontaneous imbibition in a core, IMPES, works fine
% not as fast as expected, can be improved
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
%% define the geometry
Nx = 20; % number of cells in x direction
Ny = 50; % number of cells in y direction
W = 0.02; % [m] length of the domain in x direction
H = 0.07; % [m] length of the domain in y direction
% m = createMesh1D(Nx, W);
m = createMeshCylindrical2D(Nx, Ny, W, H); % creates a 2D mesh
%% define the physical parametrs
krw0_ww = 0.3;
krw0_ow = 1.0;
kro0_ww = 0.6;
kro0_ow = 0.76;
nw = 2.4;
no = 2.0;
sor_ww=0.1;
sor_ow=0.12;
swc_ww=0.09;
swc_ow=0.09;
SF=createFaceVariable(m, 0.0); % 1 is water wet, 0 is oil wet
krw0=krw0_ww*SF+krw0_ow*(1-SF);
kro0=kro0_ww*SF+kro0_ow*(1-SF);
sor=sor_ww*SF+sor_ow*(1-SF);
swc=swc_ww*SF+swc_ow*(1-SF);
sws=@(sw, sor, swc)((sw>swc).*(sw<1-sor).*(sw-swc)./(1-sor-swc)+(sw>=1-sor));
kro=@(sw, kro0, sor, swc)((sw>=swc).*kro0.*(1-sws(sw, sor, swc)).^no+(sw<swc).*(1+(kro0-1)./swc.*sw));
krw=@(sw, krw0, sor, swc)((sw<=1-sor).*krw0.*sws(sw, sor, swc).^nw+(sw>1-sor).*(-(1-krw0)./sor.*(1.0-sw)+1.0));
dkrwdsw=@(sw, krw0, sor, swc)((sw<=1-sor).*nw.*krw0.*(1./(1-sor-swc)).*sws(sw, sor, swc).^(nw-1)+(sw>1-sor).*((1-krw0)./sor));
dkrodsw=@(sw, kro0, sor, swc)((sw>=swc).*(-kro0.*no.*(1-sws(sw, sor, swc)).^(no-1))./(-swc-sor+1)+(sw<swc).*((kro0-1)./swc));
p0 = 100e5; % [bar] pressure
pin = 150e5; % [bar] injection pressure at the left boundary
u_in= 1.0/(24*3600); % [m/s] equal to 1 m/day
sw0 = swc_ww+0.01;
% sw0(10:end-10, 10:end-10)=swc+0.2;
% sw0 = swc+0.1; % initial water saturation
sw_in = 1;
mu_oil = 2e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 0.005e-12; % [m^2] average reservoir permeability
phi0 = 0.45; % average porosity
teta_ow=deg2rad(30);
gama_ow=0.03; % N/m
labda=2.4;
eps1=1e-6;
clx=0.05;
cly=0.05;
V_dp=0.01; % Dykstra-Parsons coef.
perm_val= k0; %field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
pce=0.1*gama_ow*cos(teta_ow)*(phi./k).^0.5;
pce_face=0.05*gama_ow*cos(teta_ow)*(arithmeticMean(phi)./geometricMean(k)).^0.5;
sco=swc+0.01;
pc=@(sw, sor, swc)(pce.*(sws(sw, sor, swc)+eps).^(-1/labda)); % it can also be defined for each block
dpc=@(sw, sor, swc)((-1/labda)*(1./(1-sor-swc)).*pce_face.*(sws(sw, sor, swc)+eps).^(-1/labda-1));
dpcdk=@(sw, sor, swc)(0.5*gama_ow*cos(teta_ow)*(geometricMean(k)./arithmeticMean(phi)).^0.5.*(sws(sw, sor, swc)+eps).^(-1/labda));
grad_phik=gradientTerm(phi./k);
% sw_plot=linspace(0,1, 10000);
% plot(sw_plot, pc(sw_plot))
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
%% Define the boundaries: all fixed Sw=1, fixed pressure everywhere(?)
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% left boundary pressure gradient
% BCp.left.a(:)=(krw(sw_in)*lw.xvalue(1,:)+kro(sw_in)*lo.xvalue(1,:)); BCp.left.b(:)=0; BCp.left.c(:)=-u_in;
% change the right boandary to constant pressure (Dirichlet)
% BCp.left.a(:)=0; BCp.left.b(:)=1; BCp.left.c(:)=p0;
BCp.right.a(:)=0; BCp.right.b(:)=1; BCp.right.c(:)=p0;
BCp.top.a(:)=0; BCp.top.b(:)=1; BCp.top.c(:)=p0;
BCp.bottom.a(:)=0; BCp.bottom.b(:)=1; BCp.bottom.c(:)=p0;
% change the left boundary to constant saturation (Dirichlet)
% BCs.left.a(:)=0; BCs.left.b(:)=1; BCs.left.c(:)=1.0-sor;
BCs.right.a(:)=0; BCs.right.b(:)=1; BCs.right.c(:)=1.0;
BCs.top.a(:)=0; BCs.top.b(:)=1; BCs.top.c(:)=1.0;
BCs.bottom.a(:)=0; BCs.bottom.b(:)=1; BCs.bottom.c(:)=1.0;
%% define the time step and solver properties
% dt = 1000; % [s] time step
% dt=(W/Nx)/u_in/20; % [s]
dt=1;
t_end = 10*3600*24; % [s] final time
eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-4; % saturation accuracy
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
oil_init=domainInt(1-sw_old);
p = p_old;
uw = -gradientTerm(p_old); % an estimation of the water velocity
%% start the main loop
% generate intial pressure profile (necessary to initialize the fully
% implicit solver)
rec_fact=0;
t_day=0;
t = 0;
dt0=dt;
dsw_alwd= 0.01;
dp_alwd= 100; % Pa
while (t<t_end)
% for i=1:5
    error_p = 1e5;
    error_sw = 1e5;
    % Implicit loop
%     while ((error_p>eps_p) || (error_sw>eps_sw))
    while(1)
        % calculate parameters
        pgrad = gradientTerm(p);
%         pcgrad=gradientTerm(pc(sw));
        sw_face = upwindMean(sw, -pgrad); % average value of water saturation
        sw_grad=gradientTerm(sw);
        sw_ave=arithmeticMean(sw);
        pcgrad=dpc(sw_ave, sor, swc).*sw_grad+dpcdk(sw_ave, sor, swc).*grad_phik;
        % solve for pressure at known Sw
        labdao = lo.*funceval(kro, sw_face, kro0, sor, swc);
        labdaw = lw.*funceval(krw, sw_face, krw0, sor, swc);
        labda = labdao+labdaw;
        % compute [Jacobian] matrices
        Mdiffp1 = diffusionTerm(-labda);
        RHSpc1=divergenceTerm(labdao.*pcgrad);
        [Mbcp, RHSbcp] = boundaryCondition(BCp);
        RHS1 = RHSpc1+RHSbcp; % with capillary
        p_new=solvePDE(m, Mdiffp1+Mbcp, RHS1);

        % solve for Sw
        pgrad = gradientTerm(p_new);
        uw=-labdaw.*pgrad;
        [Mbcsw, RHSbcsw] = boundaryCondition(BCs);
        RHS_sw=-divergenceTerm(uw);
        sw_new=solveExplicitPDE(sw_old, dt, RHS_sw, BCs, phi);

        error_p = max(abs((p_new.value(:)-p.value(:))./p_new.value(:)))
        error_sw = max(abs(sw_new.value(:)-sw.value(:)))
        dt_new=dt*min(dp_alwd/error_p, dsw_alwd/error_sw);
        % assign new values of p and sw
        if error_sw>dsw_alwd
            dt=dt*(dsw_alwd/error_sw)
        else
            t=t+dt;
            p = p_new;
            sw = sw_new;
            p_old = p;
            sw_old = sw;
            dt=min(dt*(dsw_alwd/error_sw), 10*dt);
            break;
        end
    end

    rec_fact=[rec_fact (oil_init-domainInt(1-sw))/oil_init];
    t_day=[t_day t];
    figure(1);visualizeCells(1-sw); drawnow;
    figure(2); plot(t_day/3600/24, rec_fact)
    xlabel('time [day]');
    ylabel('recovery factor');
    title([num2str(t/3600/24) ' day']); drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function pass = test_fevalm( pref )
% Test chebfun2/fevalm

if ( nargin == 0)
    pref = chebfunpref;
end

rng(2016);
tol = 100*pref.cheb2Prefs.chebfun2eps;

% Check empty chebfun2:
f = chebfun2;
s = 2*rand(5,1) - 1;
t = 2*rand(5,1) - 1;
B = fevalm(f, s, t);
pass(1) = isempty( B );

% Check symmetric function:
f = chebfun2(@(x,y) cos(x.*y));
s = 2*rand(5,1) - 1;
t = 2*rand(5,1) - 1;
[ss, tt] = meshgrid( s, t);
A = feval(f, ss, tt);
B = fevalm(f, s, t);
pass(2) = norm( A - B ) < tol;

% Check essentially one dimensional function:
f = chebfun2(@(x,y) cos(y-.1));
s = 2*rand(5,1) - 1;
t = 2*rand(5,1) - 1;
[ss, tt] = meshgrid( s, t);
A = feval(f, ss, tt);
B = fevalm(f, s, t);
pass(3) = norm( A - B ) < tol;

% Check complex-valued function:
f = chebfun2(@(z) cos(z));
s = 2*rand(6,1) - 1;
t = 2*rand(6,1) - 1;
[ss, tt] = meshgrid( s, t);
A = feval(f, ss, tt);
B = fevalm(f, s, t);
pass(4) = norm( A - B ) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function V = spm_write_vol(V,Y)
% Write an image volume to disk, setting scales and offsets as appropriate
% FORMAT V = spm_write_vol(V,Y)
% V (input)  - a structure containing image volume information (see spm_vol)
% Y          - a one, two or three dimensional matrix containing the image voxels
% V (output) - data structure after modification for writing.
%
% Note that if there is no 'pinfo' field, then SPM will figure out the
% max and min values from the data and use these to automatically determine
% scalefactors.  If 'pinfo' exists, then the scalefactor in this is used.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

use_offset = false;

if ndims(Y)>3, error('Can only handle a maximum of 3 dimensions.'), end
dim = [size(Y) 1 1 1];
if ~all(dim(1:3) == V.dim(1:3)),
    error('Incompatible dimensions.');
end

if ~isfield(V,'pinfo'),
    V.pinfo = [1;0;0];
    rescal  = true;
elseif ~all(isfinite(V.pinfo(1:2))) || V.pinfo(1) == 0,
    V.pinfo(1:2) = [1;0];
    rescal  = true;
else
    rescal = false;
end

if rescal,
    % Set scalefactors and offsets
    %-----------------------------------------------------------------------
    dt           = V.dt(1);
    s            = find(dt == [2 4 8 256 512 768]);
    if isempty(s),
        V.pinfo(1:2) = [1;0];
    else
        dmnmx        = [0 -2^15 -2^31 -2^7 0 0 ; 2^8-1 2^15-1 2^31-1 2^7-1 2^16-1 2^32-1];
        dmnmx        = dmnmx(:,s);
        mxs          = zeros(dim(3),1)+NaN;
        mns          = zeros(dim(3),1)+NaN;

        for p=1:dim(3),
            tmp    = double(Y(:,:,p));
            tmp    = tmp(isfinite(tmp));
            if ~isempty(tmp),
                mxs(p) = max(tmp);
                mns(p) = min(tmp);
            end
        end

        mx = max(mxs(isfinite(mxs)));
        mn = min(mns(isfinite(mns)));
        if isempty(mx), mx = 0; end;
        if isempty(mn), mn = 0; end;
        if mx~=mn,
            if use_offset,
                V.pinfo(1,1) = (mx-mn)/(dmnmx(2)-dmnmx(1));
                V.pinfo(2,1) = (dmnmx(2)*mn-dmnmx(1)*mx)/(dmnmx(2)-dmnmx(1));
            else
                if dmnmx(1)<0,
                    V.pinfo(1) = max(mx/dmnmx(2),mn/dmnmx(1));
                else
                    V.pinfo(1) = mx/dmnmx(2);
                end
                V.pinfo(2) = 0;
            end
        else
            V.pinfo(1,1) = mx/dmnmx(2);
            V.pinfo(2,1) = 0;
        end
    end
end

%-Create and write image
%-----------------------------------------------------------------------
V = spm_create_vol(V);
V = spm_write_plane(V,Y,':');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Check rate of convergence of FEM for 3D elliptic interface problem
%     -div(A grad u)  = f,    x\in \Omega
%      where A is a piecewise constant on Omega^+ and Omega^-.
%
% Domain: Rectangular domain: [xmin,xmax] X [ymin,ymax] X [zmin,zmax]
% Mesh: Structured rectangular mesh.
% Method: Conforming Lagrange Type FEMs (Q1: Tri-Linear FEM).
%
% Last modified by Xu Zhang 08/06/2020.

%% Geometry and Boundary Conditions
clear
close all
%clc

domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1]; % Dirichelet BC

%% Finite Element Type
femtype = 'P1';
disp(['FEM Type =  Conforming ', femtype]);

%% Initial Partition
nx0 = 10;
ny0 = nx0;
nz0 = nx0;

%% Task
showErr = 0;
showMesh = 0;
computErr = 1;

%% PDE
test = 2;
switch test
    case 0
        pde = poissonNonPoly3D;
    case 1 % circular interface
        r = pi/4; bm = 1; bp = pi/(2*r^2);
        x0 = 0; y0 = 0; z0 = 0; rx = pi/4; ry = pi/4; rz = pi/4;
        pde = elli3DcircIntf(bm,bp,r,x0,y0,z0,rx,ry,rz);
    case 2 % orthotorus interface
        domain = [-1.2,1.2,-1.2,1.2,-1.2,1.2];
        bm = 1; bp = 1;
        rx = 1; ry = 0.075; rz = 3;
        pde = elli3DorthocircIntf(bm,bp,rx,ry,rz);
    case 3 % line interface
        bm = 1; bp = 100;
        rx = 1; ry = 0; rz = 1; cx = pi/10; cy = 0; cz = 0;
        cm = 1; cp = 1; a = 1;
        pde = elli3DlinIntf(bm,bp,cx,cy,cz,rx,ry,rz,a,cm,cp);
end

%% Max Iteration
maxIt = 16;
time = zeros(7,maxIt);

for i = 1:maxIt

    %% 1. Generate Mesh
    tic
    nx = nx0 + 10*(i-1); h = (domain(2) - domain(1))/nx;
    ny = ny0 + 10*(i-1);
    nz = nz0 + 10*(i-1);
    time(1,i) = nx;
    disp(' ')
    disp('*******************************************************************************')
    disp(['Partition =  ',int2str(nx),' X ',int2str(ny),' X ',int2str(nz)]);
    disp(' ')

    mesh = genMesh3D(domain, nx, ny, nz);
    disp(['number of element =  ', int2str(length(mesh.t))]);
%     mesh = genIntfMesh3D(mesh,pde.intf);
%     disp(['number of interface element =  ', int2str(-min(mesh.tLoc)),...
%         ', is ', num2str(100*-min(mesh.tLoc)/length(mesh.t)), '% of all elements']);
    if showMesh == 1
        tetramesh(mesh.t,mesh.p,ones(size(mesh.t)));
    end
    time(2,i) = toc;

    %% 2. Generate FEM DoF
    tic
    fem = genFEM3D(mesh,femtype,bc);
    disp(['number of DoF =  ', int2str(length(fem.p))]);
    time(3,i) = toc;

    %% 3. Assemble Matrix
    tic
    disp(' '); disp('Start Assembling Matrix');
    matrix = genMatEll3D(pde,mesh,fem);
    time(4,i) = toc;

    %% 3. Solve the linear system Au = f
    tic
    % L = ichol(matrix.A,struct('michol','on'));
    L = ichol(matrix.A);
    [u,flag,relres,iter,resvec] = pcg(matrix.A,matrix.f,1e-8,300,L,L');
    time(5,i) = toc;
    tu = pde.exactu(fem.p(:,1),fem.p(:,2),fem.p(:,3));
    uh = tu; uh(fem.mapper) = u;

    %% 4. Postprocess: Calculating Errors
    tic
    errND = max(abs(uh - tu)); % Error on nodes
    err.nd = errND; err.inf = 0; err.l2 = 0; err.h1 = 0;
    if computErr == 1
        disp(' ');  disp('Start computing error in Inf norm');
        errInf = getErrInf3D(uh,pde.exactu,fem,[0,0,0]);
        eNorm = 'L2'; disp(['Start computing error in ',eNorm,'  norm']);
        [errL2,errL2K] = getErr3D(uh, pde.exactu, fem, eNorm);
        eNorm = 'H1x'; disp(['Start computing error in ',eNorm,' norm']);
        [errH1x,errH1xK] = getErr3D(uh, pde.Dxu, fem, eNorm);
        eNorm = 'H1y'; disp(['Start computing error in ',eNorm,' norm']);
        [errH1y,errH1yK] = getErr3D(uh, pde.Dyu, fem, eNorm);
        eNorm = 'H1z'; disp(['Start computing error in ',eNorm,' norm']);
        [errH1z,errH1zK] = getErr3D(uh, pde.Dzu, fem, eNorm);
        err.inf = max([errND,errInf]);
        err.l2 = errL2;
        err.h1 = sqrt(errH1x^2+errH1y^2+errH1z^2);
    end
    time(6,i) = toc;
    time(7,i) = 1e6*sum(time(2:6,i))/length(mesh.t);

    %% 5: Output

    disp(' ')
    disp('Errors')
    disp('Node        Inf norm    L2 norm     H1 norm')
    formatSpec = '%6.4e  %6.4e  %6.4e  %6.4e\n';
    fprintf(formatSpec, err.nd, err.inf, err.l2, err.h1)

    if i > 1
        format short
        rNd = log(err0.nd/err.nd)./log(h0/h);
        rInf = log(err0.inf/err.inf)./log(h0/h);
        rL2 = log(err0.l2/err.l2)./log(h0/h);
        rH1 = log(err0.h1/err.h1)./log(h0/h);
        disp(' ')
        disp('Convergence Rate')
        disp('Node        Inf norm    L2 norm     H1 norm')
        formatSpec = '%6.4f      %6.4f      %6.4f      %6.4f\n';
        fprintf(formatSpec, rNd, rInf, rL2, rH1)
    end
    err0 = err; h0 = h;

    disp(' '); disp('CPU Time')
    disp('   N     Mesh     FEM      Matrix   Solve    Error    Time/1M cell')
    formatSpec = '%4i  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f   %7.2f\n';
    fprintf(formatSpec, time(:,1:i))

    %% 6. Plot Solution and Error
    if showErr == 1
    end
end

repo: ifem


text:
function [ c_est, r_est, c_err, r_err ] = ellipsoid_random_centralize ( ...
  m, n, c, r )

%*****************************************************************************80
%
%% ELLIPSOID_RANDOM_CENTRALIZE estimates the center and radius of an ellipsoid.
%
%  Discussion:
%
%    We generate N points on the M dimensional ellipsoid of center C,
%    and radiuses R.  We seek to estimate the values of C and R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of random points to generate.
%
%    Input, real C(M,1), the center.
%
%    Input, real R(M,1), the radiuses.
%
%    Output, real C_EST(M,1), R_EST, the estimated center and radius.
%
%    Output, real C_ERR, R_ERR, the errors in the estimates.
%
  c = c(:);
  r = r(:);
%
%  Construct a diagonal matrix from R.
%
  a = diag ( r );
%
%  Multiply it by a random orthogonal matrix.
%
  q = r8mat_orth_uniform ( m );
  a = q * a * q';
%
%  Now generate points on the ellipsoid surface (X-C)' * A * (X-C) = R * R.
%
  rr = 1.0;
  x = ellipsoid_surface_sample ( m, n, a, c, rr );

  c_est = sum ( x, 2 ) / n;
  c_err = norm ( c - c_est );

  r_mat = x - repmat ( c_est, 1, n );
  r_mat = r_mat.^2;
  r_mat = sum ( r_mat, 1 );
  r_mat = sqrt ( r_mat );
  r_est = sum ( r_mat ) / n;
% r_err = abs ( r - r_est );
  r_err = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function m_grid(varargin)
% M_GRID make a grid on a map.
%        M_GRID('parameter','value',...) with any number (or no)
%        optional parameters is used to draw a lat/long grid for a
%        previously initialized map projection.
%
%        The optional parameters allow the user
%        to control the look of the grid. These parameters are listed
%        by M_GRID('get'), with default parameters in M_GRID('set');
%
%        see also M_PROJ

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
%  19/6/97 - set visibility of titles and so forth to 'on' (they
%            default to 'off' when axes visibility is turned off)
%  2/11/97 - for M5.1, the old way of making the patch at the bottom (i.e.
%            by rearranging the axes children) instead causes matlab to lose
%            track of titles. Try a different fix.
% 11/01/98 - Added way of making longitude lines cut off to prevent crowding near poles (you have
%            to specify a vector for allowabale latitudes for this to work).
% 16/02/98 - Made a little fudge to allow the user to fully specify grid location
%            without getting the edge points. It doesn't quite work if only *one* edge
%            point is desired....but I hope it will be OK that way.
% 19/02/98 - PC-users complain about layers getting out of order! Their fault for using
%            such an awful OS...however (with help from Eric Firing) I think I have
%            a fix.
%  7/04/98 - Another fix to grid locations to not automatically add edge points
%            (as requested by EF)
%  7/05/98 - Added 'fancy' outline box.
% 14/11/98 - Changed tag names from m_* to m_grid_*.
% 11/07/99 - Apparently fontname changing didn't work (thanks to Dave McCollum)
% 28/04/04 - Changed m_grid_color code so it works right under unix; old
%            way retained for windows (ugh).
% 16/10/05 - Kirk Ireson discovered that the way to fix those annoying 'cut-throughs'
%            in fancy_box was to add a 'large' zdata...so I've adapted his fix in
%            fancybox and fancybox2.
% 21/11/06 - added 'backcolor'
% 16/4/07  - sorted ticklabels when user-specified (prevents an odd problem near in
%            azimuthal projections).
% 4/DEc/11 - isstr to ischar
% 7/Dec/11 - cd /ocean/rich/home/matlab 3.2.3 compatibility
% 8/Sep/13 - added 'tickstyle' parameter
% 27/Sep/13 - matlab 2013b out, includes graphic bug. Workaround provided by
%             Corinne Bassin.
% 10/Jul/14 - in 2014a BITMAX starts not to be used, changed to FLINTMAX...
% 13/Nov/14 - 2014b graphics changes are biting; a number of version-dependent
%              fixes implemented.
% 17/Nov/17 - All kinds of "new graphics" changes in the background patch
% 27/Nov/17 - made a specialized fix for xaxis tick directions in Robinson and
%             Sinusoidal projections. (problem noted by J. Tsoa, but I didn't like
%             the way he fixed it. Also improved consistency of "lake"
%             handling with non-white background colours, and more Octave
%             compatibility.
%  8/Feb/18  - changed ^o to char(176) in labels. (thanks Brian Dushaw)
%  21/06/2018 added axes option to specify the axes to draw the grid on (instead of
%              gca) - thanks to Yoann Ladroit.
%  23/12/2018 - aesthetic improvement to mollweide projection.

%
% Note that much of the work in generating line data
% is done by calls to the individual projections -
% most of M_GRID is concerned with the mechanics of plotting

% These structures are initialized by m_proj()

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

% Otherwise we are drawing a grid!

% Set current projection to geographic
Currentmap=m_coord('set');
m_coord(MAP_PROJECTION.coordsystem.name);

% Default parameters for grid

ax_idx=find(strcmpi('axes',varargin));
if isempty(ax_idx)
    ax=gca;
else
    ax=varargin{ax_idx+1};
end

xtick=6;
ytick=6;
xlabels=NaN;
ylabels=NaN;
gcolor='k';
ggridcolor=[];
gbackcolor='w'; %%get(gcf,'color');
glinestyle=':';
glinewidth=get(ax,'linewidth');
gbox='on';
gfontsize=get(ax,'fontsize');
gfontname=get(ax,'fontname');
gxaxisloc=get(ax,'xaxislocation');
gyaxisloc=get(ax,'yaxislocation');
gtickdir=get(ax,'tickdir');
gticklen=get(ax,'ticklength'); gticklen=gticklen(1);
gxticklabeldir='middle';
gyticklabeldir='end';
gtickstyle='dm';

dpatch=5; % interpolation factor for fancy grids

% Parse parameter list for options. I really should do some
% error checking here, but...

k=1;
while k<=length(varargin)
  switch lower(varargin{k}(1:3))
    case 'box'
      gbox=varargin{k+1};
    case 'xti'
      if length(varargin{k})==5
        xtick=sort(varargin{k+1});   % Added 'sort' here for people who put things in
      else                           % a random order near poles
        xlabels=varargin{k+1};
      end
    case 'yti'
      if length(varargin{k})==5
        ytick=sort(varargin{k+1});
      else
        ylabels=varargin{k+1};
      end
    case 'xla'
      gxticklabeldir=varargin{k+1};
    case 'yla'
      gyticklabeldir=varargin{k+1};
    case 'col'
      gcolor=varargin{k+1};
    case 'gri'
      ggridcolor=varargin{k+1};
    case 'bac'
      gbackcolor=varargin{k+1};
    case 'lin'
      switch lower(varargin{k}(1:5))
         case 'linew'
           glinewidth=varargin{k+1};
         case 'lines'
           glinestyle=varargin{k+1};
      end
    case 'fon'
       switch lower(varargin{k}(1:5))
         case 'fonts'
           gfontsize=varargin{k+1};
         case 'fontn'
           gfontname=varargin{k+1};
       end
    case 'xax'
      gxaxisloc=varargin{k+1};
    case 'yax'
      gyaxisloc=varargin{k+1};
    case 'tic'
      switch lower(varargin{k}(1:5))
        case 'tickl'
           gticklen=varargin{k+1};
        case 'tickd'
           gtickdir=varargin{k+1};
	    case 'ticks'
	       gtickstyle=varargin{k+1};
      end
    case {'get','usa'}
      disp('      ''axes'',( gca | axis handle)');
      disp('      ''box'',( ''on'' | ''fancy'' | ''off'' )');
      disp('      ''xtick'',( num | [value1 value2 ...])');
      disp('      ''ytick'',( num | [value1 value2 ...])');
      disp('      ''xticklabels'',[label1;label2 ...]');
      disp('      ''yticklabels'',[label1;label2 ...]');
      disp('      ''xlabeldir'', ( ''middle'' | ''end'' )');
      disp('      ''ylabeldir'', ( ''end'' | ''middle'' )');
      disp('      ''ticklength'',value');
      disp('      ''tickdir'',( ''in'' | ''out'' )');
      disp('      ''tickstyle'',(''dm'' | ''da'' | ''dd''  )');  % deg-min or abbreviated deg-min or decimal-deg
      disp('      ''color'',colorspec');
      disp('      ''gridcolor'',colorspec');
      disp('      ''backgroundcolor'',colorspec');
      disp('      ''linewidth'', value');
      disp('      ''linestyle'', ( linespec | ''none'' )');
      disp('      ''fontsize'',value');
      disp('      ''fontname'',name');
      disp('      ''XaxisLocation'',( ''bottom'' | ''middle'' | ''top'' ) ');
      disp('      ''YaxisLocation'',( ''left'' | ''middle'' | ''right'' ) ');
      return;
    case 'set'
      disp(['      axes = ' ax]);
      disp(['      box = ' gbox]);
      disp(['      xtick = ' num2str(xtick)]);
      disp(['      ytick = ' num2str(ytick)]);
      disp(['      ticklength = ' num2str(gticklen)]);
      disp(['      tickdir = ' gtickdir]);
      disp(['      tickstyle = ' gtickstyle]);
      disp(['      xlabeldir = ' gxticklabeldir]);
      disp(['      ylabeldir = ' gyticklabeldir]);
      disp(['      color = ' gcolor]);
      disp(['      gridcolor = ' gcolor]);
      disp(['      backgroundcolor - ' gbackcolor]);
      disp(['      linewidth = ' num2str(glinewidth)]);
      disp(['      linestyle = ' glinestyle]);
      disp(['      fontsize = ' num2str(gfontsize)]);
      disp(['      fontname = ' gfontname]);
      disp(['      XaxisLocation = ' gxaxisloc]);
      disp(['      YaxisLocation = ' gyaxisloc]);
      return;
  end
  k=k+2;
end

if isempty(ggridcolor)
    ggridcolor=gcolor;
end

if strcmp(gbox,'fancy')
  if strcmp(MAP_VAR_LIST.rectbox,'on') || strcmp(MAP_VAR_LIST.rectbox,'circle')
   gbox='on';
   warning([' No fancy outline with ''rectbox'' set to ''' MAP_VAR_LIST.rectbox '''']);
  end
end

[X,Y]=feval(MAP_PROJECTION.routine,'box');

% These color parts make Octave fail so I moved them here, otherwise
% they should be before this if block

% If there are lakes around, make them the backgroundcolor..
hh=get(ax,'children');
hh_tags=get(hh,'tag');
if length(hh)==1, hh_tags={hh_tags}; end
for i=1:length(hh_tags)
   if ~isempty(strfind(hh_tags{i},'_lake')) && strcmp(get(hh(i),'type'),'patch')
     set(hh(i),'facecolor',gbackcolor);
   end
end

% Set the axes colour so lakes get coloured right (if a coastline added
% next)

set(ax,'color',gbackcolor');

if MAP_PROJECTION.IsOctave

  %  Octave handling here all provided by Jamie Tsao

  % Note that both gnuplot and fltk doesn't really work well with the z-axis,
  % so here we simply ignore it.

  % The tag is important for M_COAST to catch the correct underlying color.
  % As a result, though, M_COAST should be called after M_GRID.

  % IN 2018, this workaround no longer needed.
  if 0,%strcmp(graphics_toolkit(), 'gnuplot')
    % Unfortunately, gnuplot also has a bug with using patch, where a linestyle
    % of 'none' would cause it to barf. This bug has been reported since 2013,
    % but apparently it is still not fixed. It may have something to do with
    % how Octave interacts with the underlying gnuplot backend.

    % On the other hand, this means I can set the linewidth to 0 and it should
    % work exactly as intended. If there is no problems, then I can refactor
    % this part, assuming ftlk doesn't barf on a zero linewidth.
    patch('xdata', X(:), 'ydata', Y(:), 'facecolor', gbackcolor, ...
          'edgecolor', 'k', 'linewidth', 0, 'tag', 'm_grid_color','parent',ax);
  else

    patch('xdata', X(:), 'ydata', Y(:), 'facecolor', gbackcolor, ...
          'edgecolor', 'k', 'linestyle', 'none', 'tag', 'm_grid_color','parent',ax);
  end

  % Now I set it at the bottom of the children list so it gets drawn first
  % (i.e. doesn't cover anything)

  % A bug in Octave 3.8.1 makes set(ax, 'children', -) insert the new
  % permutation before the old one instead of replacing it. This only occurs
  % when the ShowHiddenHandles flag is set to on instead of off. Although fltk
  % is not affected, gnuplot still is. Combined with fltk's lack of support for
  % unicode and hence the degree sign, this is still important. Furthermore,
  % this bug causes delete() to function incorrectly.

  show = get(0, 'ShowHiddenHandles');
  set(0, 'ShowHiddenHandles', 'on');

  hh = get(ax,'children');
  htags = get(hh,'tag');
  k = strcmp('m_grid_color', htags);
  hht = hh;
  hh(k) = [];
  hh = [hh;hht(k)];

  % Continuing on the bug, there is no particular need for why this needs to be
  % placed under even hidden handles. Even so, as a workaround, we can grab all
  % hidden handles, set them to be visible, turn off flag, reorder, then set
  % them back as hidden, as noted below.

  hht_inv = strcmp('off', get(hht, 'HandleVisibility'));
  set(hht(hht_inv), 'HandleVisibility', 'on');
  set(0, 'ShowHiddenHandles', 'off');
  set(ax,'children',hh);
  set(hht(hht_inv), 'HandleVisibility', 'off');

  set(0, 'ShowHiddenHandles', show);

  % Patch clipping does not work correctly in Octave, so I have to
  % manually set the Z-limits to that fancy boxes (with 3D patches)
  % works correctly.
  zlim([-1 MAP_PROJECTION.LARGVAL]);
else   % For MATLAB

   % This is a very problematic part of the code. It turns out the the interaction between
   % PATCH objects and CONTOURF objects does not work correctly in the Painters renderer -
   % this is true in all versions up to 7.7 at least. Patches with large negative Z just
   % don't get drawn under contourgroup patches.
   %
   % There are several possible workarounds:
   %
   %  1) Make sure you use the 'v6' option in contourf calls (see m_contourf.m to see
   %     how I have tried to do that for some versions of matlab)
   %      - problem: the 'v6' option is going away soon, also you may want the
   %                contourgroup object that the v6 option destroys.
   %
   %  2) Change the renderer to something else:
   %        set(gcf,'renderer','opengl') or
   %        set(gcf,'renderer','zbuffer')
   %      - problem: These other renderers are not available on all systems, they may also
   %                give less-precise results.
   %
   %  3) Use the painters renderer, but reorder the children so that the patch is at the
   %     bottom of the list (painters cares about child order)
   %      - problem: sometimes the child order is rearranged if you click on the figure,
   %                 also (at least in some versions of matlab) this causes labels to
   %                 disappear.
   %
   % With version 7.4 onwards I have discovered that reordering the children apparently
   % is Mathworks-blessed (c.f. the UISTACK function).  So I am going to try to implement
   % the latter as a default.
   %
   % Now, putting in a white background works under linux (at least) and
   % NOT under windows...I don't know about macs.
   %
   % This broke with matlab 2014b - but eventually (I sent in a "bug
   % report" and found it was a "feature") I learned of the 'sortmethod'
   % tag, which made it work again - as long as I didn't use negative z
   % values in the patch command. I don't know WHY that is...

   if ~MAP_PROJECTION.newgraphics
        patch('xdata',X(:),'ydata',Y(:),'zdata',-MAP_PROJECTION.LARGVAL*ones(size(X(:))),'facecolor',gbackcolor,...
	   'edgecolor','k','linestyle','none','tag','m_grid_color','Parent',ax);
	  % Now I set it at the bottom of the children list so it gets drawn first (i.e. doesn't
	  % cover anything)
	 show=get(0, 'ShowHiddenHandles');
	 set(0, 'ShowHiddenHandles', 'on');
	 hh=get(ax,'children');
	 htags = get(hh,'tag');
	 k = strcmp('m_grid_color',htags);
	 hht = hh;
	 hh(k) = [];
	 hh = [hh;hht(k)];
	 set(ax,'children',hh);
	 set(0, 'ShowHiddenHandles', show);

   else  % after 2014b

     patch('xdata',X(:),'ydata',Y(:),'zdata',-MAP_PROJECTION.LARGVAL*ones(size(X(:))),'facecolor',gbackcolor,...
             'edgecolor','k','linestyle','none','tag','m_grid_color','parent',ax);
%     patch('xdata',X(:),'ydata',Y(:),'facecolor',gbackcolor,...
%             'edgecolor','k','linestyle','none','tag','m_grid_color');
	 set(ax,'sortmethod','childorder');
     uistack(findobj(ax,'tag','m_grid_color'),'bottom');

   end

end

% X-axis labels and grid

if ~isempty(xtick)

 % Tricky thing - if we are drawing a map with the poles, its nasty when the lines get too close
 % together. So we can sort of fudge this by altering MAP_VAR_LIST.lats to be slightly smaller,
 % and then changing it back again later.
 fudge_north='n';fudge_south='n';
 if ~isempty(ytick) && length(ytick)>1
  if MAP_VAR_LIST.lats(2)==90
    fudge_north='y';
    MAP_VAR_LIST.lats(2)=ytick(end);
  end
  if MAP_VAR_LIST.lats(1)==-90
    fudge_south='y';
    MAP_VAR_LIST.lats(1)=ytick(1);
  end
 end

 [X,Y,lg,lgI]=feval(MAP_PROJECTION.routine,'xgrid',xtick,gxaxisloc,gtickstyle);
 [labs,scl]=m_labels('lon',lg,xlabels,gtickstyle);

 % Draw the grid. Every time we draw something, I first reshape the matrices into a long
 % row so that a) it plots faster, and b) all lines are given the same handle (which cuts
 % down on the number of children hanging onto the axes).

 [n,m]=size(X);
 line(reshape([X;NaN+ones(1,m)],(n+1)*m,1),reshape([Y;NaN+ones(1,m)],(n+1)*m,1),...
      'linestyle',glinestyle,'color',ggridcolor,'linewidth',0.1,'tag','m_grid_xgrid','parent',ax);

 % Get the tick data
 [ltx,lty,utx,uty]=maketicks(X,Y,gticklen,gtickdir);

 % Draw ticks if labels are on top or bottom (not if they are in the middle)

 if strcmp(gxticklabeldir,'middle')
  if lgI==size(X,1) && strcmp(gxaxisloc,'top')  % Check to see if the projection supports this option.
   vert='bottom';horiz='center';drawticks=1;
   xx=utx(1,:);yy=uty(1,:);rotang=atan2(diff(uty),diff(utx))*180/pi+90;
  elseif lgI==1 && strcmp(gxaxisloc,'bottom')
   vert='top';horiz='center';drawticks=1;
   xx=ltx(1,:);yy=lty(1,:);rotang=atan2(diff(lty),diff(ltx))*180/pi-90;
  else
   vert='middle';horiz='center';lgIp1=lgI+1;drawticks=0;
   xx=X(lgI,:); yy=Y(lgI,:);rotang=atan2(Y(lgIp1,:)-Y(lgI,:),X(lgIp1,:)-X(lgI,:))*180/pi-90;
  end
 else
  if lgI==size(X,1) && strcmp(gxaxisloc,'top')  % Check to see if the projection supports this option.
   vert='middle';horiz='left';drawticks=1;
   xx=utx(1,:);yy=uty(1,:);rotang=atan2(diff(uty),diff(utx))*180/pi+180;
  elseif lgI==1 && strcmp(gxaxisloc,'bottom')
   vert='middle';horiz='right';drawticks=1;
   xx=ltx(1,:);yy=lty(1,:);rotang=atan2(diff(lty),diff(ltx))*180/pi;
  else
   vert='top';horiz='center';lgIp1=lgI+1;drawticks=0;
   xx=X(lgI,:); yy=Y(lgI,:);rotang=atan2(Y(lgIp1,:)-Y(lgI,:),X(lgIp1,:)-X(lgI,:))*180/pi;
  end
 end

 % For some projections it is better if the x-grid labels are horizontal or
 % vertical

 if any(strcmp(MAP_PROJECTION.name,{'Robinson','Sinusoidal'}) )
     if strcmp(gxaxisloc,'bottom') || strcmp(gxaxisloc,'top')
         if strcmp(gxticklabeldir,'middle')
            rotang = zeros(size(rotang));    % flat
            yy=repmat(mean(yy),size(yy));    % Put them all on the same y-value.
         else
            rotang = zeros(size(rotang))+90;  % upright
         end
     end
 end

 if strcmp(gbox,'fancy')
    if gtickdir(1)=='i'
      fancybox(lg,MAP_VAR_LIST.longs,'xgrid','bottom',dpatch,gticklen,gtickstyle,ax);
      drawticks=0;
    else
      fancybox2(lg,MAP_VAR_LIST.longs,'xgrid','bottom',dpatch,gticklen,gtickstyle,ax);
    end
 end
 if drawticks
   [n,m]=size(ltx);
   line(reshape([ltx;NaN+ones(1,m)],(n+1)*m,1),reshape([lty;NaN+ones(1,m)],(n+1)*m,1),...
        'linestyle','-','color',gcolor,'linewidth',glinewidth,'tag','m_grid_xticks-lower','clipping','off','parent',ax);
   line(reshape([utx;NaN+ones(1,m)],(n+1)*m,1),reshape([uty;NaN+ones(1,m)],(n+1)*m,1),...
        'linestyle','-','color',gcolor,'linewidth',glinewidth,'tag','m_grid_xticks-upper','clipping','off','parent',ax);
 end

 % Add the labels! (whew)

 ik=1:size(X,2);

 for k=ik
   [rotang(k), horizk, vertk] = upright(rotang(k), horiz, vert);
   text(xx(k),yy(k),labs{k},'horizontalalignment',horizk,'verticalalignment',vertk, ...
        'rotation',rotang(k),'fontsize',gfontsize*scl(k),'color',gcolor,...
        'tag','m_grid_xticklabel','fontname',gfontname,'parent',ax);
 end

 if fudge_north=='y'
   MAP_VAR_LIST.lats(2)=90;
 end
 if fudge_south=='y'
   MAP_VAR_LIST.lats(1)=-90;
 end

end

if ~isempty(ytick)
 % Y-axis labels and grid

 [X,Y,lt,ltI]=feval(MAP_PROJECTION.routine,'ygrid',ytick,gyaxisloc,gtickstyle);
 [labs,scl]=m_labels('lat',lt,ylabels,gtickstyle);

 % Draw the grid
 [n,m]=size(X);
 line(reshape([X;NaN+ones(1,m)],(n+1)*m,1),reshape([Y;NaN+ones(1,m)],(n+1)*m,1),...
      'linestyle',glinestyle,'color',ggridcolor,'linewidth',0.1,'tag','m_grid_ygrid','parent',ax);

 % Get the tick data
 [ltx,lty,utx,uty]=maketicks(X,Y,gticklen,gtickdir);

 % Draw ticks if labels are on left or right (not if they are in the middle)
 if strcmp(gyticklabeldir,'end')
  if ltI==size(X,1) && strcmp(gyaxisloc,'right')  % Check to see if the projection supports this option.
   horiz='left';vert='middle';drawticks=1;
   xx=utx(1,:);yy=uty(1,:);rotang=atan2(diff(uty),diff(utx))*180/pi+180;
  elseif ltI==1 && strcmp(gyaxisloc,'left')
   horiz='right';vert='middle';drawticks=1;
   xx=ltx(1,:);yy=lty(1,:);rotang=atan2(diff(lty),diff(ltx))*180/pi;
  else
   horiz='center';vert='top';ltIp1=ltI+1;drawticks=0;
   xx=X(ltI,:); yy=Y(ltI,:);rotang=atan2(Y(ltIp1,:)-Y(ltI,:),X(ltIp1,:)-X(ltI,:))*180/pi;
  end
 else
  if ltI==size(X,1) && strcmp(gyaxisloc,'right')  % Check to see if the projection supports this option.
   horiz='center';vert='top';drawticks=1;
   xx=utx(1,:);yy=uty(1,:);rotang=atan2(diff(uty),diff(utx))*180/pi+270;
  elseif ltI==1 && strcmp(gyaxisloc,'left')
   horiz='center';vert='bottom';drawticks=1;
   xx=ltx(1,:);yy=lty(1,:);rotang=atan2(diff(lty),diff(ltx))*180/pi+90;
  else
   horiz='left';vert='middle';ltIp1=ltI+1;drawticks=0;
   xx=X(ltI,:); yy=Y(ltI,:);rotang=atan2(Y(ltIp1,:)-Y(ltI,:),X(ltIp1,:)-X(ltI,:))*180/pi+90;
  end
 end

 % Added special case for latitudes. Dec/2018
 if any(strcmp(MAP_PROJECTION.name,{'Mollweide'}) )
     if strcmp(gyaxisloc,'left')
         rotang=-lt;
     elseif strcmp(gyaxisloc,'right')
         rotang=lt;
     end
 end

 if strcmp(gbox,'fancy')
    if gtickdir(1)=='i'
      fancybox(lt,MAP_VAR_LIST.lats,'ygrid','left',dpatch,gticklen,gtickstyle,ax);
      drawticks=0;
    else
      fancybox2(lt,MAP_VAR_LIST.lats,'ygrid','left',dpatch,gticklen,gtickstyle,ax);
    end
 end
 if drawticks
   [n,m]=size(ltx);
   line(reshape([ltx;NaN+ones(1,m)],(n+1)*m,1),reshape([lty;NaN+ones(1,m)],(n+1)*m,1),...
        'linestyle','-','color',gcolor,'linewidth',glinewidth,'tag','m_grid_yticks-left','clipping','off','parent',ax);
   line(reshape([utx;NaN+ones(1,m)],(n+1)*m,1),reshape([uty;NaN+ones(1,m)],(n+1)*m,1),...
        'linestyle','-','color',gcolor,'linewidth',glinewidth,'tag','m_grid_yticks-right','clipping','off','parent',ax);
 end

 % Finally - the labels!
 ik=1:size(X,2);

 for k=ik
   [rotang(k), horizk, vertk] = upright(rotang(k), horiz, vert);
   text(xx(k),yy(k),labs{k},'horizontalalignment',horizk,'verticalalignment',vertk,...
        'rotation',rotang(k),'fontsize',gfontsize*scl(k),'color',gcolor,...
	'tag','m_grid_yticklabel','fontname',gfontname,'parent',ax);
 end

end

% Draw the plot box

if strcmp(gbox,'on')
  [X,Y]=feval(MAP_PROJECTION.routine,'box');
  line(X(:),Y(:),'linestyle','-','linewidth',glinewidth,'color',gcolor,'tag','m_grid_box','clipping','off','parent',ax);
end

% Give a 1-1 aspect ratio and get rid of the matlab-provided axes stuff.

if isempty(strfind(version,'R2013b'))  % 27/Sept/13 - Handling for 2013b provided by CB.
    set(ax,'visible','off',...
            'dataaspectratio',[1 1 1],...
            'xlim',MAP_VAR_LIST.xlims,...
            'ylim',MAP_VAR_LIST.ylims);
else
    set(ax,'visible','off',...
            'dataaspectratio',[1 1 1e16],...
            'xlim',MAP_VAR_LIST.xlims,...
            'ylim',MAP_VAR_LIST.ylims);
end

set(ax,'tag',['m_grid_' MAP_PROJECTION.name]);

set(get(ax,'title'),'visible','on');
set(get(ax,'xlabel'),'visible','on');
set(get(ax,'ylabel'),'visible','on');

% Set coordinate system back

m_coord(Currentmap.name);

end

%-------------------------------------------------------------
% upright simply turns tick labels right-side up while leaving
% their positions unchanged.
% Sat  98/02/21 Eric Firing
%
function   [rotang, horiz, vert] = upright(rotang, horiz, vert);
if isnan(rotang), rotang=0; end   % Added in 2014!
if rotang > 180, rotang = rotang - 360; end
if rotang < -180, rotang = rotang + 360; end
if rotang > 90
   rotang = rotang - 180;
elseif rotang < -90
   rotang = 180 + rotang;
else
   return    % no change needed.
end
switch horiz(1)
   case 'l'
      horiz = 'right';
   case 'r'
      horiz = 'left';
end
switch vert(1)
   case 't'
      vert = 'bottom';
   case 'b'
      vert = 'top';
end

end

%--------------------------------------------------------------------------
function [L,fs]=m_labels(dir,vals,uservals,tickstyle)
% M_LONLABEL creates longitude labels
%         Default values are calculated automatically when the grid is
%         generated. However, the user may wish to specify the labels
%         as either numeric values or as strings (in the usual way
%         for axes).
%
%         If auto-labelling occurs, minutes are labelled in a different
%         (smaller) fontsize than even degrees.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997

% If the user has specified [] (i.e. no labels), we return blanks.

if isempty(uservals)
  L=cellstr(char(' '*ones(length(vals),1)));
  fs=1.0*ones(length(L),1);
  return;
end

% If the user has specified strings, we merely need to make
% sure that there are enough to cover all ticks.

if any(ischar(uservals))
  L=cellstr( uservals((rem([0:length(vals)-1],length(uservals))+1),:) );
  fs=1.0*ones(length(L),1);
  return;
end

% Otherwise we are going to have to generate labels from numeric
% data.

if length(uservals)==1 && isnan(uservals)  % use default values
  vals=vals(:)'; % make it a row.
else                                       % or ones provided
  lv=length(vals);
  vals=uservals(:)';
  while length(vals)<lv
    vals=[vals uservals(:)'];
  end
end

% longitudes and latitudes have some differences....
% and keep a label for 0 or 180 where we don't want a letter
if strfind(dir,'lat')
  labname=['S';'N';' '];
else
  labname=['W';'E';' '];
  vals=rem(vals+540,360)-180;
end

%i=[vals<0;vals>0;vals==0];  % get the 'names' (i.e. N/S or E/W)
%modified by Changyong He the label for 0 or +/-180 will not show the
%an E or W.
i=[-180~=vals & vals<0; 180~=vals & vals>0; vals==0 | abs(vals) == 180];

L=cell(length(vals),1);
fs=ones(length(vals),1);

if strcmp(tickstyle,'dm')
   vals=abs(vals);             % Convert to +ve values

   % For each label we have different options:
   %  1 - even degrees are just labelled as such.
   %  2 - ticks that fall on even minutes are just labelled as even minutes
   %      in a smaller fontsize.
   %  3 - fractional minutes are labelled to 2 decimal places in the
   %      smaller fontsize.
   for k=1:length(vals)
     if rem(vals(k),1)==0
       nam=find(i(:,k));
       L{k}=sprintf([' %3.0f' char(176) labname(nam) ' '],vals(k));
     elseif abs(vals*60-round(vals*60))<.01
       L{k}=sprintf([' %2.0f'' '],rem(vals(k),1)*60);
       fs(k)=0.75;
     else
       L{k}=sprintf([' %2.2f'' '],rem(vals(k),1)*60);
       fs(k)=0.75;
     end
   end

   % In most cases, the map will have at least one tick with an even degree label,
   % but for very small regions (<1 degree in size) this won't happen so we
   % want to force one label to show degrees *and* minutes.

   if ~any(fs==1)
    k=round(length(vals)/2);
    nam=find(i(:,k));
    L{k}={sprintf([' %3.0f' char(176) labname(nam) ' '],fix(vals(k))),...
	  sprintf([' %2.2f'' '],rem(vals(k),1)*60)};
    fs(k)=1;
   end

elseif strcmp(tickstyle,'da')
   % Abbreviated option suggested by Brooke Jones - Dec/2018
   % For each label we have different options:
   %  1 - even degrees are just labelled as such.
   %  2 - 2 decimal place intervals use 2 decimal places
   %  3 - the rest fo to 4
   for k=1:length(vals)
     if rem(vals(k),1)==0
       nam=find(i(:,k));
       L{k}=sprintf([' %4.0f '],vals(k));
     elseif abs(vals*100-round(vals*100))<0.01
       L{k}=sprintf([' %2.0f'' '],rem(abs(vals(k)),1)*60);
       fs(k)=0.6;
     else
       L{k}=sprintf([' %2.2f'' '],rem(abs(vals(k)),1)*60);
       fs(k)=0.6;
     end
   end

   % In most cases, the map will have at least one tick with an even degree label,
   % but for very small regions (<1 degree in size) this won't happen so we
   % want to force one label to show degrees *and* minutes.

   if ~any(fs==1)
    k=round(length(vals)/2);
    nam=find(i(:,k));
    L{k}={sprintf([' %4.0f    '],fix(vals(k))),...
	  sprintf([' %2.2f'' '],abs(rem(vals(k),1))*60)};
    fs(k)=1;
   end

elseif strcmp(tickstyle,'dd')
   vals=abs(vals);             % Convert to +ve values
   % For each label we have different options:
   %  1 - even degrees are just labelled as such.
   %  2 - 2 decimal place intervals use 2 decimal places
   %  3 - the rest fo to 4
   for k=1:length(vals)
     if rem(vals(k),1)==0
       nam=find(i(:,k));
       L{k}=sprintf([' %3.0f' char(176) labname(nam) ' '],vals(k));
     elseif abs(vals*100-round(vals*100))<0.01
       L{k}=sprintf([' %2.2f'],vals(k));
       fs(k)=0.75;
     else
       L{k}=sprintf([' %6.4f'],vals(k));
       fs(k)=0.75;
     end
   end

  % write code.
end
end

%---------------------------------------------------------
function [ltx,lty,utx,uty]=maketicks(X,Y,gticklen,gtickdir)
% MAKETICKS makes the axis ticks.
%           AXes ticks are based on making short lines at
%           the end of the grid lines X,Y.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997

global MAP_VAR_LIST

tlen=gticklen*max( diff(MAP_VAR_LIST.xlims),diff(MAP_VAR_LIST.ylims));

lx=sqrt((X(2,:)-X(1,:)).^2+(Y(2,:)-Y(1,:)).^2);

if strcmp(gtickdir,'out')
  ltx=[X(1,:)-tlen*(X(2,:)-X(1,:))./lx;X(1,:)];
  lty=[Y(1,:)-tlen*(Y(2,:)-Y(1,:))./lx;Y(1,:)];
else
  ltx=[X(1,:);X(1,:)+tlen*(X(2,:)-X(1,:))./lx];
  lty=[Y(1,:);Y(1,:)+tlen*(Y(2,:)-Y(1,:))./lx];
end

lx=sqrt((X(end,:)-X(end-1,:)).^2+(Y(end,:)-Y(end-1,:)).^2);

if strcmp(gtickdir,'out')
  utx=[X(end,:)-tlen*(X(end-1,:)-X(end,:))./lx;X(end,:)];
  uty=[Y(end,:)-tlen*(Y(end-1,:)-Y(end,:))./lx;Y(end,:)];
else
  utx=[X(end,:);X(end,:)+tlen*(X(end-1,:)-X(end,:))./lx];
  uty=[Y(end,:);Y(end,:)+tlen*(Y(end-1,:)-Y(end,:))./lx];
end
end

%---------------------------------------------------------
function fancybox(vals,lims,gridarg1,gridarg2,dpatch,gticklen,gridarg3,ax);
%
%  FANCYBOX  - draws fancy outlines for either top/bottom or left/right sides,
%              depending on calling parameters.

global MAP_PROJECTION

% Get xlocations including endpoints
xval=sort([lims(1) vals(vals>lims(1) & vals<lims(2)) lims(2)]);
% Add all half-way points as well.
xval=sort([xval,xval(1:end-1)+diff(xval)/2]);

% Interpolate extra points to handle curved boundary conditions.

xval=(xval(1:end-1)'*ones(1,dpatch)+diff(xval)'*[0:dpatch-1]/dpatch)';
xval=[xval(:);lims(2)];

% Get lat/long positions for everything

[X2,Y2,lg2,lgI2]=feval(MAP_PROJECTION.routine,gridarg1,xval,gridarg2,gridarg3);
[l2x,l2y,u2x,u2y]=maketicks(X2,Y2,gticklen,'in');

if gridarg1(1)=='x', sig=1; else sig=-1; end

id=[1:dpatch size(l2x,2)+[-dpatch+1:0]];
dx=diff(l2x(:,id));
l2x(2,id)=l2x(2,id)+diff(l2y(:,id)).*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;
l2y(2,id)=l2y(2,id)-dx.*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;
dx=diff(u2x(:,id));
u2x(2,id)=u2x(2,id)-diff(u2y(:,id)).*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;
u2y(2,id)=u2y(2,id)+dx.*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;

% Now actually draw the patches.

% Added the z-values 16/Oct/05.

px=prod(size(l2x));
kk=[0:(dpatch*4):px-3]'*ones(1,dpatch*2+2);
kk=kk+ones(size(kk,1),1)*[1 2:2:(dpatch*2+2) (dpatch*2+1):-2:3];
patch(reshape(u2x(kk),size(kk,1),size(kk,2))',...
      reshape(u2y(kk),size(kk,1),size(kk,2))',...
      repmat(MAP_PROJECTION.LARGVAL  ,size(kk,2),size(kk,1)),'w','edgecolor','k','clipping','off','tag','m_grid_fancybox1','parent',ax);
patch(reshape(l2x(kk),size(kk,1),size(kk,2))',...
      reshape(l2y(kk),size(kk,1),size(kk,2))',...
      repmat(MAP_PROJECTION.LARGVAL-1,size(kk,2),size(kk,1)),'k','clipping','off','tag','m_grid_fancybox1','parent',ax);

kk=[dpatch*2:(dpatch*4):px-3]'*ones(1,dpatch*2+2);
kk=kk+ones(size(kk,1),1)*[1 2:2:(dpatch*2+2) (dpatch*2+1):-2:3];
patch(reshape(l2x(kk),size(kk,1),size(kk,2))',...
      reshape(l2y(kk),size(kk,1),size(kk,2))',...
      repmat(MAP_PROJECTION.LARGVAL  ,size(kk,2),size(kk,1)),'w','edgecolor','k','clipping','off','tag','m_grid_fancybox1','parent',ax);
patch(reshape(u2x(kk),size(kk,1),size(kk,2))',...
      reshape(u2y(kk),size(kk,1),size(kk,2))',...
      repmat(MAP_PROJECTION.LARGVAL-1,size(kk,2),size(kk,1)),'k','clipping','off','tag','m_grid_fancybox1','parent',ax);
end

%---------------------------------------------------------
function fancybox2(vals,lims,gridarg1,gridarg2,dpatch,gticklen,gridarg3,ax)
%
%  FANCYBOX  - draws fancy outlines for either top/bottom or left/right sides,
%              depending on calling parameters.

global MAP_PROJECTION

% Get xlocations including endpoints
xval=sort([lims(1) vals(vals>lims(1) & vals<lims(2)) lims(2)]);
% Add all half-way points as well.
xval=sort([xval,xval(1:end-1)+diff(xval)/2]);

% Interpolate extra points to handle curved boundary conditions.

xval=(xval(1:end-1)'*ones(1,dpatch)+diff(xval)'*[0:dpatch-1]/dpatch)';
xval=[xval(:);lims(2)];

% Get lat/long positions for everything

[X2,Y2,lg2,lgI2]=feval(MAP_PROJECTION.routine,gridarg1,xval,gridarg2,gridarg3);
[l2x,l2y,u2x,u2y]=maketicks(X2,Y2,gticklen,'in');

if gridarg1(1)=='x'
    sig=1;
else
    sig=-1;
end

id=[1:dpatch size(l2x,2)+[-dpatch+1:0]];
dx=diff(l2x(:,id));
l2x(2,id)=l2x(2,id)+diff(l2y(:,id)).*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;
l2y(2,id)=l2y(2,id)-dx.*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;
dx=diff(u2x(:,id));
u2x(2,id)=u2x(2,id)-diff(u2y(:,id)).*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;
u2y(2,id)=u2y(2,id)+dx.*([dpatch:-1:1 -1:-1:-dpatch]/(dpatch))*sig;

% Now actually draw the patches.

% Added large z-values 16/Oct/05

px=prod(size(l2x));
kk=[0:(dpatch*2):px-3]'*ones(1,dpatch*2+2);
kk=kk+ones(size(kk,1),1)*[1 2:2:(dpatch*2+2) (dpatch*2+1):-2:3];
 patch(reshape(l2x(kk),size(kk,1),size(kk,2))',...
       reshape(l2y(kk),size(kk,1),size(kk,2))',...
       repmat(MAP_PROJECTION.LARGVAL-1,size(kk,2),size(kk,1)),...
       'w','edgecolor','k','clipping','off','linewidth',.2,'tag','m_grid_fancybox2','parent',ax);
 patch(reshape(u2x(kk),size(kk,1),size(kk,2))',...
       reshape(u2y(kk),size(kk,1),size(kk,2))',...
       repmat(MAP_PROJECTION.LARGVAL-1,size(kk,2),size(kk,1)),...
       'w','edgecolor','k','clipping','off','linewidth',.2,'tag','m_grid_fancybox2','parent',ax);
kk=[0:(dpatch*2):size(l2x,2)-dpatch-1]'*ones(1,dpatch+1);
kk=(kk+ones(size(kk,1),1)*[1:dpatch+1])';
[k1,k2]=size(kk);
line(reshape(mean(l2x(:,kk)),k1,k2),reshape(mean(l2y(:,kk))',k1,k2),...
     repmat(MAP_PROJECTION.LARGVAL,k1,k2),'color','k','clipping','off','tag','m_grid_fancybox2','parent',ax);

kk=[dpatch:(dpatch*2):size(l2x,2)-dpatch-1]'*ones(1,dpatch+1);
kk=(kk+ones(size(kk,1),1)*[1:dpatch+1])';
[k1,k2]=size(kk);
line(reshape(mean(u2x(:,kk))',k1,k2),reshape(mean(u2y(:,kk))',k1,k2),...
     repmat(MAP_PROJECTION.LARGVAL,k1,k2),'color','k','clipping','off','tag','m_grid_fancybox2','parent',ax);

end

author: goGPS-Project


text:
%%%  PosControl_Sim
clear
path('./icon/',path);
Init;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;
%% Initial condition
ModelInit_PosE = [0, 0, -100];
ModelInit_VelB = [0, 0, 0];
ModelInit_AngEuler = [0, 0, 0];
ModelInit_RateB = [0, 0, 0];
ModelInit_Rads = 0;
%% control parameter
%Attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
%Position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
PosControl_Sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function Wop = normalize_standardize_threshold_wavelet_factory_1d(N,filter_options,scat_options,epsilon)

if nargin <4
    epsilon =2^(-20);
end

filters = filter_bank(N, filter_options);

ren_op=@(X)(func_output(...
            @renorm_wavelet_layer_1d,[1,2],X,epsilon));

std_op= @(X,m)(func_output(...
            @standardize_wavelet_layer_1d,[1,2],X,scat_options.sigmas{m+1}));

	for m = 0:scat_options.M
		filt_ind = min(numel(filters), m+1);
        if m<scat_options.M
		Wop{m+1} = @(X)(threshold_wavelet_layer_1d(...
            std_op(...
            ren_op(func_output(@wavelet_layer_1d,[1,2],X,...
            filters{filt_ind},scat_options)),m),scat_options.threshold));
        else
        Wop{m+1} = @(X)(threshold_wavelet_layer_1d(...
                      ren_op(func_output(@wavelet_layer_1d,[1,2],X,...
            filters{filt_ind},scat_options)),scat_options.threshold));
	end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function shim = cvx_glpk( shim )

% CVX_SOLVER_SHIM	GLPK interface for CVX.
%   This procedure returns a 'shim': a structure containing the necessary
%   information CVX needs to use this solver in its modeling framework.

if ~isempty( shim.solve ),
    return
end
if isempty( shim.name ),
    fname = 'glpk.m';
    ps = pathsep;
    shim.name = 'GLPK';
    shim.config = struct( 'dualize', 1, 'nonnegative', 1, ...
    	'integer', 1, 'binary', 1 );
    flen = length(fname);
    fpaths = which( fname, '-all' );
    if ~iscell(fpaths),
      fpaths = { fpaths };
    end
    old_dir = pwd;
    oshim = shim;
    shim = [];
    for k = 1 : length(fpaths),
        fpath = fpaths{k};
        if ~exist( fpath, 'file' ) || any( strcmp( fpath, fpaths(1:k-1) ) ),
            continue
        end
        new_dir = fpath(1:end-flen-1);
        cd( new_dir );
        tshim = oshim;
        tshim.fullpath = fpath;
        tshim.version = 'unknown';
        tshim.location = new_dir;
        if isempty( tshim.error ),
            tshim.solve = @solve;
            if k ~= 1,
                tshim.path = [ new_dir, ps ];
            end
        end
        shim = [ shim, tshim ]; %#ok
    end
    cd( old_dir );
    if isempty( shim ),
        shim = oshim;
        shim.error = 'http://glpkmex.sourceforge.net/';
    end
else
    shim.solve = @solve;
end

function [ x, status, tol, iters, y, z ] = solve( At, b, c, nonls, quiet, prec, settings )

n  = length( c );
m  = length( b );
lb = -Inf(n,1);
ub = +Inf(n,1);
vtype = 'C';
vtype = vtype(ones(n,1));
ctype = 'S';
ctype = ctype(ones(m,1));
rr = zeros(0,1);
cc = rr;
vv = rr;
zinv = rr;
is_ip = false;
for k = 1 : length( nonls ),
    temp = nonls( k ).indices;
    tt = nonls( k ).type;
    switch tt,
    case 'nonnegative',
      lb(temp) = 0;
    case 'integer',
      is_ip = true;
      vtype(temp) = 'I';
    case 'binary',
      is_ip = true;
      vtype(temp) = 'I';
      lb(temp) = 0;
      ub(temp) = 1;
    otherwise,
      cvx_throw( 'GLPK does not support nonlinear constraints.' );
    end
end
if ~isempty(rr),
  znorm = [1:n]';
  znorm(zinv) = [];
  rr = [ rr ; znorm ];
  cc = [ cc ; znorm ];
  vv = [ vv ; ones(size(znorm)) ];
  reord = sparse( rr, cc, vv, n, n );
  At = reord' * At;
  c  = reord' * c;
end
if quiet,
  param.msglev = 0;
else
  param.msglev = 2;
end
param.scale = 128;
param.tolbnd = prec(1);
param.toldj = prec(1);
param.tolobj = prec(1);
[ xx, fmin, errnum, extra ] = cvx_run_solver( @glpk, c, At', b, lb, ub, ctype, vtype, 1, param, 'xx', 'fmin', 'errnum', 'extra', settings, 9 ); %#ok
iters = [];
x = full( xx );
y = full( extra.lambda );
z = full( extra.redcosts );
if ~isempty( rr ),
  x = reord * x;
  z = reord * z;
  z(zinv) = z(zinv) * 0.5;
end
status = 'Failed';
switch errnum,
case 0,
  switch extra.status,
  case 2,
    if is_ip,
      status = 'Suboptimal';
    elseif errnumn == 0,
      status = 'Solved';
    else
      status = 'Inaccurate/Solved';
    end
  case {3,4},
    status = 'Infeasible';
  case 5,
    status = 'Solved';
  case 6,
    status = 'Unbounded';
  end
case 10,
  status = 'Infeasible';
case 11,
  status = 'Unbounded';
case {5,17,15,19}
  status = 'Failed';
case {6,7,8,9,13,14},
  switch extra.status,
  case {2,5}
    if is_ip,
      status = 'Suboptimal';
    else
      status = 'Inaccurate/Solved';
    end
  case { 3,4 },
    status = 'Inaccurate/Infeasible';
  case 6,
    status = 'Inaccurate/Unbounded';
  end
end
if strcmp(status,'Failed'),
  tol = Inf;
elseif strncmp(status,'Inaccurate/',11),
  tol = prec(3);
else
  tol = prec(2);
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function [au] = km2au(km)
% Convert length from kilometers to astronomical units.
% Chad A. Greene 2012
au = km*6.684587122671e-9;


text:
function Plugin_ReverseK
global VVar

VVar.Kz=-VVar.Kz;
VVar.Ky=-VVar.Ky;
VVar.Kx=-VVar.Kx;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [S] = L1QP_FeatureSign_Set(X, B, Sigma, beta, gamma)

[dFea, nSmp] = size(X);
nBases = size(B, 2);

% sparse codes of the features
S = sparse(nBases, nSmp);

A = B'*B + 2*beta*Sigma;

for ii = 1:nSmp,
    b = -B'*X(:, ii);
%     [net] = L1QP_FeatureSign(gamma, A, b);
    S(:, ii) = L1QP_FeatureSign_yang(gamma, A, b);
end

author: thomas-koehler

repo: SupER


text:
%The basic idea behind the false position method is similar to the bisection
%method in that we continuously shrink the interval the root lies on until
%the algorithm converges on the root. Unlike the bisection method, the false
%position method does not halve the interval with each iteration. Instead of
%using the midpoint of a and b to create the new interval, the false position
%method uses the x-intercept of the line connecting f(a) and f(b). This
%algorithm converges faster than the bisection method.

%INPUTS:
%Function handle f
%endpoint a
%endpoint b
%maximum tolerated error

%OUTPUTS:
%An approximated value for the root of f within the defined interval.

%Written by MatteoRaso

function y = false_position(f, a, b, error)
  if ~(f(a) < 0)
    disp("f(a) must be less than 0")
  elseif ~(f(b) > 0)
    disp("f(b) must be greater than zero")
  else
    c = 100000;
    while abs(f(c)) > error
      %Formula for the x-intercept
      c = -f(b) * (b - a) / (f(b) - f(a)) + b;
      if f(c) < 0
        a = c;
      else
        b = c;
      endif
      disp(f(c))
    endwhile
    x = ["The root is approximately located at ", num2str(c)];
    disp(x)
    y = c;
  endif
endfunction

author: TheAlgorithms

repo: MATLAB-Octave


text:
function AddedPath = requireFEXpackage(FEXSubmissionID)
%Function requireFEXpackage -
%installs Matlab Central File Exchange (FEX) submission
%with given ID into the directory chosen by the user.
%A new FEX submissions may use previous FEX submissions as its part.
%The function 'requireFEXpackage' helps in adding those previous
%submissions to the user's MATLAB installation.
%
%This function is a part of File Exchange submission 31069.
%Download the entire submission:
%http://www.mathworks.com/matlabcentral/fileexchange/31069
%
% SYNTAX:
%    AddedPath = requireFEXpackage(FEXSubmissionID)
%
% INPUT:
%    ID of the required submission to File Exchange
%
% OUTPUT:
%    the path to that submission added to the user's MATLAB path.
%
% HOW TO CALL:
%    The command
%           P = requireFEXpackage(8277)
%    will download and install the package with ID 8277
%    (namely, nice 'fminsearchbnd' by John D'Errico)
%
% EXAMPLES -- HOW TO USE:
%
% EXAMPLE 1 (using 'exist' command):
%
%     % first, somewhere in the very beginning of your code,
%     % check if the function 'fminsearchbnd' from the FEX package 8277
%     % is on your MATLAB path, and if it is not there,
%     % require the FEX package 8277:
%     if ~(exist('fminsearchbnd', 'file') == 2)
%         P = requireFEXpackage(8277);  % fminsearchbnd is part of 8277
%     end
%
%     % Then just use 'fminsearchbnd' where you need it:
%     syms x
%     RosenbrockBananaFunction = @(x) (1-x(1)).^2 + 100*(x(2)-x(1).^2).^2;
%     x = fminsearchbnd(RosenbrockBananaFunction,[3 3])

% EXAMPLE 2 (using 'try-catch' command):
%
%     syms x
%     RosenbrockBananaFunction = @(x) (1-x(1)).^2+100*(x(2)-x(1).^2).^2;
%     try
%        % if function 'fminsearchbnd' already exists in your MATLAB
%        % installation, just use it:
%        x = fminsearchbnd(RosenbrockBananaFunction,[3 3])
%     catch
%        % if function 'fminsearchbnd' is not present in your MATLAB
%        % installation, first get the package 8277 (to which it belongs)
%        % from the MATLAB Central File Exchange (FEX)
%        P = requireFEXpackage(8277);  % fminsearchbnd is part of 8277
%        % and then use that function:
%        x = fminsearchbnd(RosenbrockBananaFunction,[3 3])
%     end
%
%
% NOTE: on Mac platform, the title of the dialog box for
% choosing the directory for installing the required FEX package
% is not shown; this is not a bug, this is how UIGETDIR works on Macs --
% see the documentation for UIGETDIR
% http://www.mathworks.com/help/techdoc/ref/uigetdir.html
%
% (C) Igor Podlubny, 2011

ID = num2str(FEXSubmissionID);

% Ask user for the confirmation of the installation
% of the required FEX package
yes = ['YES, Install package ' ID];
no = 'NO, do not install';
userchoice = questdlg(['The Matlab function/toolbox, which you are running, ' ...
    'requires the presence of the package ' ID ...
    ' from Matlab Central File Exchange.' ...
    'Would you like to install the FEX package ' ID ' now?'] , ...
	['Required package ' ID], ...
	yes, no, yes);

% Handle response
switch userchoice
    case yes,
        install = 1;
    case no,
        install = 0;
    otherwise,
        install = 0;
end

if install == 1
    baseURL = 'http://www.mathworks.com/matlabcentral/fileexchange/';
    query = '?download=true';

    location = uigetdir(pwd, ['Select the directory for installing the required FEX package' ID ]);
    if location ~= 0
        % download package 'ID' from Matlab Central File Exchange
        filetosave = [location filesep ID '.zip'];
        FEXpackage = [baseURL ID query];
        [f,status] = urlwrite(FEXpackage,filetosave);
        if status==0
            warndlg(['No connection to Matlab Central File Exchange,' ' or package ' ID ' does not exist.' ...
                ' Package ' ID ' has not been installed. ' ...
                ' Check you internet settings and the ID of the required package, and try again. '] , ...
                ['No connection to Matlab Central File Exchange' ' or package ' ID ' does not exist'], ...
                'modal');
            AddedPath = '';
            return
        end
        % unzip the downloaded file to the subdirectory 'ID'
        todir = [location filesep ID];
        % if the directory 'ID' doesn't exist at given location, create it
        if ~(exist([location filesep ID], 'dir') == 7)
            mkdir(location, ID);
        end
        try
            unzip(filetosave, todir);
            % after unzipping, delete the downloaded ZIP file
            delete(filetosave);
            % prepend the paths to the downloaded package to the MATLAB path
            P = genpath([location filesep ID]);
            path(P,path);
        catch
            % if the FEX package is not ZIP, then it is a single m-file
            % just move the file to the ID directory
            [pathstr, name, ext] = fileparts(filetosave);
            movefile(filetosave, [todir filesep name '.m']);
            P = genpath([location filesep ID]);
            path(P,path);
        end
    else
        P = '';
    end

    AddedPath = P;
else
    AddedPath = '';
end

if install == 1,
    % Ask user about reviewing and saving the modified MATLAB path,
    % and take him to PATHTOOL, if the user wants to save the modified path
    yes = 'YES, I want to review and save the MATLAB path';
    no = 'NO, I don''t want to save the path permanently';
    userchoice = questdlg(['After adding the package ' ID ...
        ' from Matlab Central File Exchange to your MATLAB installation,' ...
        ' the MATLAB path has been modified accordingly. ', ...
        'Would you like to review and save the modified MATLAB path?'] , ...
        'Review and save the modified MATLAB path for future use?', ...
        yes, no, yes);

    % Handle response
    switch userchoice
        case yes,
            pathtool;
        case no,
        otherwise,
    end
end

author: Sable

repo: mcbench-benchmarks


text:
% Replace zeros f0 values by linear interpolation
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function f0sout = fillf0(f0sin)

    f0sout = f0sin;

    idx = find(f0sin(:,2)~=0);
    if length(idx)==1
        f0sout(:,2) = f0sin(idx,2);
    elseif length(idx)>1
        f0sout(:,2) = 2.^(interp1(f0sin(idx,1), log2(f0sin(idx,2)), f0sin(:,1), 'linear', NaN));
    end

    f0sout(:,2) = interp1_fixnan(f0sout(:,2));

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function y = newvar( prob, name, siz, str, geo )
error( nargchk( 2, 5, nargin ) );

%
% Check problem
%

if ~isa( prob, 'cvxprob' ),
    error( 'First argument must be a cvxprob object.' );
end
p = index( prob );
global cvx___

%
% Check name
%

if ischar( name ),
    if ~isempty( name ),
        if size( name, 1 ) ~= 1,
            error( 'Second argument must be a string or a subscript structure array.' );
        elseif ~isvarname( name ),
            error( 'Invalid variable name: %s', name );
        elseif isfield( cvx___.problems( p ).variables, name ),
            error( 'Variable already exists: %s', name );
        end
        nstr = struct( 'type', '.', 'subs', name );
    else
        nstr = [];
    end
elseif ~isstruct( name ),
    error( 'Second argument must be a string or a subscript structure array.' );
else
    nstr = name;
    name = cvx_subs2str( name, [ 1, 0, 1 ], 1 );
    name(1) = [];
    if ~isequal( nstr(1).type, '.' ),
        error( 'Invalid subscript structure: first element must be a field reference.' );
    end
end

%
% Retrieve an existing variable, and check for conflicts
%

vars = cvx___.problems( p ).variables;
if ~isempty( nstr ),
    try
        y = builtin( 'subsref', vars, nstr );
    catch
        y = [];
    end
    if nargin == 2,
        if ~isempty( y ), return; end
        error( [ 'Unknown variable: ', name ] );
    elseif ~isempty( y ),
        error( [ 'Duplicate variable name: ', name ] );
    end
elseif nargin == 2,
    error( 'Second argument must be a non-empty string or subscript structure array.' );
else
    y = []; %#ok
end

%
% Check for conflict with dual variable
%

if ~isempty( nstr ) && isfield( cvx___.problems( p ).duals, nstr(1).subs ),
    error( 'Primal/dual variable name conflict: %s', nstr(1).subs );
end

%
% Quick exit for retrieval mode
%

if nargin == 2,
    y = cvx___.problems( p ).variables;
    try
        y = subsref( y, name );
        return
    catch
        error( [ 'Unknown variable: ', nstr ] );
    end
end

%
% Create the variable
%

if isa( siz, 'cvx' ),

    %
    % Out of an existing object
    %

    y = newsubst( prob, siz );

else

    %
    % Check size
    %

    [ temp, siz ] = cvx_check_dimlist( siz, true );
    if ~temp,
        error( 'Invalid size vector.' );
    end

    %
    % Check structure
    %

    len = prod( siz );
    if nargin < 4 || isempty( str ),
        dof = len;
        str = [];
    elseif ~isnumeric( str ) || ndims( str ) > 2 || size( str, 2 ) ~= len,
        error( 'Fourth argument must be a valid structure matrix.' );
    elseif nnz( str ) == 0,
        error( 'Structure matrix cannot be identically zero.' );
    else
        temp = any( str, 2 );
        dof = full( sum( temp ) );
        if dof ~= length( temp ),
            str = str( temp, : );
        end
    end

    %
    % Geometric flag
    %

    if nargin < 5 || isempty( geo ),
        geo = false;
    elseif ~isnumeric( geo ) && ~islogical( geo ) || length( geo ) ~= 1,
        error( 'Fifth argument must be true or false.' );
    end

    %
    % Allocate the raw variable data
    %

    geo  = any( geo( : ) );
    odim = length( cvx___.reserved );
    ndim = odim + 1 : odim + ( geo + 1 ) * dof;
    cvx___.reserved(    ndim, 1 ) = 0;
    cvx___.vexity(      ndim, 1 ) = 0;
    cvx___.canslack(    ndim, 1 ) = true;
    cvx___.readonly(    ndim, 1 ) = p;
    cvx___.geometric(   ndim, 1 ) = 0;
    cvx___.logarithm(   ndim, 1 ) = 0;
    cvx___.exponential( ndim, 1 ) = 0;
    if geo,
        ndim2 = ndim( 1 : dof );
        ndim  = ndim( dof + 1 : end );
        cvx___.vexity(      ndim,  1 ) = +1;
        cvx___.canslack(    ndim,  1 ) = true;
        cvx___.geometric(   ndim,  1 ) = +1;
        cvx___.exponential( ndim2, 1 ) = ndim( : );
        cvx___.logarithm(   ndim,  1 ) = ndim2( : );
    end
    cvx___.x = [];
    cvx___.y = [];

    %
    % Create the variable object
    %

    if dof == 0,
        str2 = sparse( 1, 0 );
    else
        str2 = sparse( ndim, 1 : dof, 1 );
        if ~isempty( str ),
            str2 = str2 * str;
        end
    end
    y = cvx( siz, str2, dof * ( 1 - 2 * geo ), false );

end

%
% If the variable is named, save it in the problem structure
%

if ~isempty( nstr ),
    try
        cvx___.problems( p ).variables = builtin( 'subsasgn', vars, nstr, y );
    catch
        error( [ 'Invalid variable name: ', name ] );
    end
end

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

repo: REST


text:
% LOGDET_CHOL  -  Computes the log-determinant of a matrix using its
% Cholesky factor.
%
% Y = LOGDET_CHOL(U)
%
% Computes the log-determinant of the matrix X whose Cholesky factor is
% the triangular matrix U. U can be lower or upper triangular.
%
% More generally, X could be a product of two instances of U and/or U'. For
% instance, X=U'*U' or X=U'*U.

% Last modified 2010-06-04
% Copyright (c) Jaakko Luttinen

function y = logdet_chol(U)

y = 2 * logdet_tri(U);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%TRAINSEG	Interactively train a color segmentation
%
%	map = trainseg(im)
%
%	Two windows are displayed, one the bivariant histogram in
%	normalized (r,g) coordinates, the other the original image.
%
%	For each pixel selected and clicked in the original image a point
%	is marked in the bivariant histogram.  By selecting numerous points
%	in the color region of interest, its extent in the (r,g) plane
%	develops.
%
%	This map can be smoothed, expanded and filled in using morphological
%	operations.
%
% SEE ALSO: colorseg imorph

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.
function map = colorseg(im)

	% convert image to (r,g) coordinates
	y = sum(im, 3);
	r = round( im(:,:,1) ./ y * 255);
	g = round( im(:,:,2) ./ y * 255);

	% display the original image

	% create and display the map
	map = zeros(256, 256);
	hm = figure
	set(gcf, 'Units', 'normalized', 'Position', [0.1 0.5 0.8 0.4])
	subplot(121)
	image(im)
	axis('equal')
	title('input image')

	subplot(122)
	image(map)
	axis('equal')
	xlabel('r');
	ylabel('g');
	colormap(gray(2));
	title('segmentation map')

	while 1,
		subplot(121)
		[y,x] = ginput(1);
		if isempty(y),
			break;
		end
		x = round(x);
		y = round(y);
		map(r(x,y), g(x,y)) = 256;
		subplot(122)
		image(map)
		drawnow
	end

	map = map';

author: petercorke

repo: machinevision-toolbox-matlab


text:
function stlButtonCallback(blockHandle)
% Callback function for STLbutton
% Allows user to specify what .stl file to use for body

    [filename,filepath] = uigetfile('.stl');               % Return filename, filepath of user-chosen file
    % Don't set value if no file is chosen, or prompt canceled.
    if ~isequal(filename,0) && ~isequal(filepath,0)
        mask = Simulink.Mask.get(blockHandle);
        fileParam = mask.getParameter('geometryFile');
        fileParam.Value = [filepath,filename];             % Update geometryFile with new filename
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function comp_enum_test ( )

%*****************************************************************************80
%
%% COMP_ENUM_TEST tests COMP_ENUM;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMP_ENUM_TEST\n' );
  fprintf ( 1, '  COMP_ENUM counts compositions;\n' );
  fprintf ( 1, '\n' );
  for n = 0 : 10
    for k = 1 : 10
      num = comp_enum ( n, k );
      fprintf ( 1, '  %6d', num );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function [iSel, iEcg, chTags, chInd] = ImaGIN_select_channels(chNames, isSEEG)
% IMAGIN_SELECT_CHANNELS Keep only channels of interest.
%
% USAGE:  [iSel, iEcg, chTags, chInd] = ImaGIN_select_channels(chNames, isSEEG=1)
%
% INPUT:
%    - chNames : Cell-array of strings
%    - isSEEG  : 1 if the data is SEEG, 0 if regular EEG
%
% OUTPUT:
%    - iSel   : Array of indices of the channels that are considered as valid
%    - iEcg   : Array of indices of the channels that are identified as ECG
%    - chTags : Array of electrode names for the selected sEEG channels
%    - chInd  : Array of indices of the selected contact on the corresponding sEEG electrode

% -=============================================================================
% This function is part of the ImaGIN software:
% https://f-tract.eu/
%
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE AUTHORS
% DO NOT ASSUME ANY LIABILITY OR RESPONSIBILITY FOR ITS USE IN ANY CONTEXT.
%
% Copyright (c) 2000-2019 Inserm U1216
% =============================================================================-
%
% Authors: Francois Tadel, 2017-2019

% Parse inputs
if (nargin < 2) || isempty(isSEEG)
    isSEEG = 1;
end
% Make sure the names are in one row
chNames = chNames(:)';

% Get all names: remove special characters
AllNames = cellfun(@(c)c(~ismember(c, ' .,?!-_@#$%^&*+*=()[]{}|/')), chNames, 'UniformOutput', 0);
AllTags  = cell(size(AllNames));
AllInd   = cell(size(AllNames));
isNoInd  = zeros(size(AllNames));
% Separate characters and numbers in the names
for i = 1:length(AllNames)
    % Find the last letter in the name
    iLastLetter = find(~ismember(AllNames{i}, '0123456789'), 1, 'last');
    if isempty(iLastLetter)
        iLastLetter = 0;
    end
    AllTags{i} = AllNames{i}(1:iLastLetter);
    % For EEG: The separate name/index does not make sense
    if ~isSEEG
        AllInd{i} = sprintf('%d',i);
        isNoInd(i) = 0;
    % If there are digits at the end of the name: use them as the index of the contact
    elseif (iLastLetter < length(AllNames{i}))
        AllInd{i} = AllNames{i}(iLastLetter+1:end);
    else
        isNoInd(i) = 1;
        AllInd{i} = '0';
    end
end
% Convert indices to double values
AllInd = cellfun(@str2num, AllInd);

% Remove all channels with indices > 18 (this is not SEEG)
% Tolerance for indices 101-118 and 201-218 (examples: electrode "T1" contact "11" => "T11",  electrode "T1" contact "2" => "T12" or "T02")
MAX_CONTACTS = 18;
if isSEEG
    uniqueTags = unique(upper(AllTags));
    for iTag = 1:length(uniqueTags)
        % Get channels of this tag
        iChTag = find(strcmpi(uniqueTags{iTag}, AllTags));
        contactInd = AllInd(iChTag);
        % Reject channels "E" without rejecting "E1"
        if any(contactInd == 0) && any(contactInd > 0)
            % Remove unwanted channels
            iChTagDel = iChTag(contactInd == 0);
            AllNames(iChTagDel) = {'XXXXX'};
            AllTags(iChTagDel) = {'XXXXX'};
            % Remove from the list, so we can process the rest of the channels
            iChTag(contactInd == 0) = [];
            contactInd(contactInd == 0) = [];
        end
        % If the tag is "X": Used by Nihon Kohden for technical tracks, but can be a SEEG electrode too. Keep only indices <= 18
        if isequal(uniqueTags{iTag}, 'X')
            iChTagExtra = iChTag(contactInd > MAX_CONTACTS);
            AllNames(iChTagExtra) = {'XXXXX'};
            AllTags(iChTagExtra) = {'XXXXX'};
        % Accept if all indices < 18
        elseif all(((contactInd >= 1) & (contactInd <= MAX_CONTACTS)))
            % OK
        % Electrode names ending in "1" or "2" (two possible nameing conventions, with or without zeros)
        elseif all(((contactInd >= 11)  & (contactInd <= 19))  | ...                  % Ending in 1, without zeros:  T11..T19, T110..T118
                   ((contactInd >= 110) & (contactInd <= 100 + MAX_CONTACTS)) | ...
                   ((contactInd >= 21)  & (contactInd <= 29))  | ...                  % Ending in 2, without zeros:  T21..T29, T210..T218
                   ((contactInd >= 210) & (contactInd <= 200 + MAX_CONTACTS))) || ...
               all(((contactInd >= 101)  & (contactInd <= 100 + MAX_CONTACTS)) | ...  % Ending in 1, with zeros:  T101..T118
                   ((contactInd >= 201)  & (contactInd <= 200 + MAX_CONTACTS)))       % Ending in 2, with zeros:  T201..T218
            % All channel indices starting with 1 => Add "1" to electrode name
            iElec1 = iChTag(((contactInd >= 11) & (contactInd <= 19))  |  (contactInd >= 101) & (contactInd <= 100 + MAX_CONTACTS));
            for i = 1:length(iElec1)
                AllTags{iElec1(i)} = [AllTags{iElec1(i)}, '1'];
                if (AllInd(iElec1(i)) < 100)
                    AllInd(iElec1(i)) = AllInd(iElec1(i)) - 10;
                elseif (AllInd(iElec1(i)) > 100)
                    AllInd(iElec1(i)) = AllInd(iElec1(i)) - 100;
                end
            end
            % All channel indices starting with 2 => Add "2" to electrode name
            iElec2 = iChTag(((contactInd >= 21) & (contactInd <= 29))  |  (contactInd >= 201) & (contactInd <= 200 + MAX_CONTACTS));
            for i = 1:length(iElec2)
                AllTags{iElec2(i)} = [AllTags{iElec2(i)}, '2'];
                if (AllInd(iElec2(i)) < 200)
                    AllInd(iElec2(i)) = AllInd(iElec2(i)) - 20;
                elseif (AllInd(iElec2(i)) > 200)
                    AllInd(iElec2(i)) = AllInd(iElec2(i)) - 200;
                end
            end
        % Mark as bad all the other channels with Label+Index (but not the channels without labels, because we want to keep the ECG)
        elseif all(contactInd > 0)
            AllNames(iChTag) = {'XXXXX'};
            AllTags(iChTag) = {'XXXXX'};
        end
        % Removed: ((length(iChTag) < 2) && ~any(chNames{iChTag(1)} == '-'))   % less than 2 (but only if not designating a bipolar montage, with '-' in the name)
        % Keep electrodes with one contact only
    end
end

% Process all the channels
iSel = [];
iEcg = [];
for i = 1:length(AllNames)
    % ECG: Accept (should be labelled as such)
    if ismember(lower(AllTags{i}), {'ecg', 'ekg'}) || ~isempty(strfind(lower(AllNames{i}), 'ecg')) || ~isempty(strfind(lower(AllNames{i}), 'ekg'))
        iEcg(end+1) = i;
    % No index or does not end with a digit
    elseif isSEEG && (isNoInd(i) || ~ismember(AllNames{i}(end), '0123456789'))
        continue;
    % Does not contain at least a letter
    elseif isSEEG && ~any(ismember(lower(AllNames{i}), 'abcdefghijklmnopqrstuvwxyz'))
        continue;
    % Contains a /: exclude
    elseif isSEEG && any(ismember(lower(chNames{i}), '/'))
        continue;
    % Unwanted labels (strict)
    elseif ismember(lower(AllTags{i}), {'xxxxx', 'mark', 'dc', 'veo', 'heo', 'dd', 'dg', 'el', 'pulse', 'spo2', 'lpar', 'rpar', 'tib'})
        % MAYBE ADD 'oc' ??
        continue;
    % Unwanted labels (including variants): EOG, EMG, MYO, MAST, REF
    elseif ~isempty(strfind(lower(AllTags{i}), 'eog')) ...
        || ~isempty(strfind(lower(AllTags{i}), 'emg')) ...
        || ~isempty(strfind(lower(AllTags{i}), 'myo')) ...
        || ~isempty(strfind(lower(AllTags{i}), 'mast')) ...
        || ~isempty(strfind(lower(AllTags{i}), 'ref'))
        continue;
    % Unwanted EEG labels
    elseif isSEEG && ismember(lower(AllTags{i}), {'cz', 'fz', 'pz', 'oz', 'nz', 'fpz'})
        continue;
    % Unwanted electrodes that are explicitely scalp electrodes (case sensitive)
    elseif isSEEG && ismember(AllNames{i}, {'sFp1','sFp2','sF4','sF3','sC3','sC4','sP4','sP3','sO2','sO1','sF8','sF7','sT8/T4','sT7/T3','sP8/T6','sP7/T5','sPz','sFz','sIO1','sIO2','sAF9','sAF10','sF9','sF10','sCB1','sCB2','sTP7','sTP9','sTP10','sTP8','sOz','sIz','sPO4','sPO3','sCP5','sCP6','sCP1','sCP2','sFT9','sFT10','sFC2','sFC1','sAF3','sAF4','sFC6','sFC5','sCPz','sP1','sPOz','sP2','sP6','sC6','sP5','sC1','sC2','sC5','sF2','sF6','sF1','sAF8','sF5','sAF7','sFpz','sFCz','sCz'})
        continue;
    % Unwanted FPx/Fx/Cx/Tx/Pz/Ox labels
    elseif isSEEG && ismember(lower(AllNames{i}), {'fp1','fp2'}) && ~any(ismember({'fp3','fp4','fp5','fp6','fp9','fp10','fp11','fp12','fp13','fp14','fp15','fp16','fp17','fp18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'f3','f4','f7','f8'}) && ~any(ismember({'f1','f2','f5','f6','f9','f10','f11','f12','f13','f14','f15','f16','f17','f18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'ft1','ft2','ft9','ft10'}) && ~any(ismember({'ft5','ft6','ft11','ft12','ft13','ft14','ft15','ft16','ft17','ft18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'tp9','tp10'}) && ~any(ismember({'tp1','tp2','tp5','tp6','tp11','tp12','tp13','tp14','tp15','tp16','tp17','tp18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'c3','c4'}) && ~any(ismember({'c1','c2','c5','c6','c7','c8','c9','c10','c11','c12','c13','c14','c15','c16','c17','c18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'t3','t4','t5','t6','t9','t10'}) && ~any(ismember({'t1','t2','t7','t8','t11','t12','t13','t14','t15','t16','t17','t18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'p3','p4'}) && ~any(ismember({'p1','p2','p5','p6','p7','p8','p9','p10','p11','p12','p13','p14','p15','p16','p17','p18'}, lower(AllNames)))
        continue;
    elseif isSEEG && ismember(lower(AllNames{i}), {'o1','o2'}) && ~any(ismember({'o3','o4','o4','o5','o6','o7','o8','o9','o10','o11','o12','o13','o14','o15','o16','o17','o18'}, lower(AllNames)))
        continue;
    % Otherwise: accept
    else
        iSel(end+1) = i;
    end
end

% Sort channels by tag and index
AllTags = AllTags(iSel);
AllInd = AllInd(iSel);
uniqueTags = unique(upper(AllTags));
iOrderSel = [];
for iTag = 1:length(uniqueTags)
    iChTag = find(strcmpi(AllTags, uniqueTags{iTag}));
    [tmp, iOrderInd] = sort(AllInd(iChTag));
    iOrderSel = [iOrderSel, iChTag(iOrderInd)];
end

% Add ECG channels at the end of the file
iSel = [iSel(iOrderSel), iEcg];

% Return lists of selected channels and electrode names
chTags = AllTags;
chInd = AllInd;

author: brainstorm-tools


text:
function [t, u] = swept_sine(f,cyc,res,plt)
% generate a swept sine curve
% function [t, u] = swept_sine(f,cyc,res,plt)
%
% inputs 4 - 3 optional
% f      frequency interval ([low high] (Hz)       class real
% cyc    number of cycles to cover (total)         class real      - optional
% res    number of points per cycle (resolution)   class integer   - optional
% plt    plot function (0 / 1)                     class integer   - optional
%
% outputs 2
% t      time increment vector (non-linear)        class real
% u      steering input content                    class real
%
% michael arant - Aug 7, 2009
%
% ex
%	[t u] = swept_sine([.1 5],25,25,1,0)

if nargin < 1; help swept_sine; error('Need inputs'); end
if ~exist('cyc','var'); cyc = ceil((f(2) / f(1)) / 5); end
if ~exist('res','var'); res = 50; end
if ~exist('plt','var'); plt = 0; end

% define sine cycles
a = linspace(0,pi*cyc*2,res*cyc)';
% generate sine
u = sin(a);

% time step vector (time between points)
ti = linspace(1/(res * f(1)),1/(res * f(2)),res*cyc)';
% time vector (point location in time) - scale by number of points and cycles
t = (cumsum(ti) - ti(1)) * cyc / res;

% proof of linear time reduction
% diff(diff(t(1:res:end)))

% debug plot
if plt
	figure; set(gcf,'color','w'); plot(u); grid on;
	title('Steering Input'); ylabel('Amplitude (deg)'); xlabel('Increment')
	figure; set(gcf,'color','w'); plot(t,u); grid on;
	title('Steering Input'); ylabel('Amplitude (deg)'); xlabel('Time (s)')
end