text:
classdef TestBilateralFilter
    %TestBilateralFilter

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_simple
            img = cv.imread(TestBilateralFilter.im, 'ReduceScale',2);
            result = cv.bilateralFilter(img);
            validateattributes(result, {class(img)}, {'size',size(img)});
        end

        function test_options
            img = cv.imread(TestBilateralFilter.im, 'ReduceScale',2);
            result = cv.bilateralFilter(img, ...
                'Diameter',7, 'SigmaColor',50, 'SigmaSpace',50, ...
                'BorderType','Default');
            validateattributes(result, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.bilateralFilter();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end

        function test_error_non_existant_option
            img = imread(TestBilateralFilter.im);
            try
                cv.bilateralFilter(img, 'foo','bar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end

        function test_error_invalid_option_value
            img = imread(TestBilateralFilter.im);
            try
                cv.bilateralFilter(img, 'BorderType','foo');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv


text:
function [ l, u ] = vand2_lu ( n, x )

%*****************************************************************************80
%
%% VAND2_LU returns the LU factors of the VAND2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Halil Oruc, George Phillips,
%    Explicit factorization of the Vandermonde matrix,
%    Linear Algebra and its Applications,
%    Volume 315, Number 1-3, 15 August 2000, pages 113-123.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real X(N), the values that define the matrix.
%
%    Output, real L(N,N), U(N,N), the LU factors of the matrix.
%
  l = zeros ( n, n );

  for i = 1 : n
    for j = 1 : i
      l(i,j) = 1.0;
      for k = 1 : j - 1
        l(i,j) = l(i,j) * ( x(i) - x(k) ) / ( x(j) - x(k) );
      end
    end
  end

  u = zeros ( n, n );

  for i = 1 : n
    for j = i : n
      u(i,j) = complete_symmetric_poly ( i, j - i, x );
      for k = 1 : i - 1
        u(i,j) = u(i,j) * ( x(i) - x(k) );
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function alpha = db2neper(alpha, y)
%DB2NEPER   Convert decibels to nepers.
%
% DESCRIPTION:
%       db2neper converts an attenuation coefficient in units of
%       dB/(MHz^y cm) to units of Nepers/((rad/s)^y m).
%
% USAGE:
%       alpha = db2neper(alpha)
%       alpha = db2neper(alpha, y)
%
% INPUTS:
%       alpha       - attenuation in dB/(MHz^y cm)
%
% OPTIONAL INPUTS:
%       y           - power law exponent (default = 1)
%
% OUTPUTS:
%       alpha       - attenuation in Nepers/((rad/s)^y m)
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 27th March 2009
%       last update - 5th December 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also neper2db

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

if nargin == 1
    y = 1;
end
alpha = 100*alpha.*(1e-6/(2*pi)).^y./(20*log10(exp(1)));

repo: Aero-matlab


text:
function res = cputime_res
    % clock resolution:
    res = 0; while (res == 0), res = abs(cputime - cputime); end
end

%!test
%! res = cputime_res;

author: Sable

repo: mcbench-benchmarks


text:
%% bfs
load_gaimc_graph('bfs_example.mat') % use the dfs example from Boost
d = bfs(A,1)

%% bipartite_matching
A = rand(10,8); % bipartite matching between random data
[val mi mj] = bipartite_matching(A);
val

%% clustercoeffs
load_gaimc_graph('clique-10');
cc = clustercoeffs(A) % they are all equal! as we expect in a clique

%% dfs
load_gaimc_graph('dfs_example.mat') % use the dfs example from Boost
d = dfs(A,1)

%% dijkstra
% Find the minimum travel time between Los Angeles (LAX) and
% Rochester Minnesota (RST).
load_gaimc_graph('airports')
A = -A; % fix funny encoding of airport data
lax=247; rst=355;
[d pred] = dijkstra(A,lax);
fprintf('Minimum time: %g\n',d(rst));
% Print the path
fprintf('Path:\n');
path =[]; u = rst; while (u ~= lax) path=[u path]; u=pred(u); end
fprintf('%s',labels{lax});
for i=path; fprintf(' --> %s', labels{i}); end, fprintf('\n');

%% dirclustercoeffs
load_gaimc_graph('celegans'); % load the C elegans nervous system network
cc=dirclustercoeffs(A);
[maxval maxind]=max(cc)
labels(maxind) % most clustered vertex in the nervous system

%% graph_draw
load_gaimc_graph('dfs_example');
graph_draw(A,xy);

%% mst_prim
load_gaimc_graph('airports'); % A(i,j) = negative travel time
A = -A; % convert to travel time.
A = max(A,A'); % make the travel times symmetric
T = mst_prim(A);
gplot(T,xy); % look at the minimum travel time tree in the US

%% scomponents
% scomponents
load_gaimc_graph('cores_example'); % the graph A has three components
ci = scomponents(A)
ncomp = max(ci)               % should be 3
R = sparse(1:size(A,1),ci,1,size(A,1),ncomp); % create a restriction matrix
CG = R'*A*R;                  % create the graph with each component
                              % collapsed into a single node.

%% load_gaimc_graph
% equivalent to load('graphs/airports.mat') run from the gaimc directory
load_gaimc_graph('airports')
% equivalent to P=load('graphs/kt-7-2.mat') run from the gaimc directory
P=load_gaimc_graph('kt-7-2.mat')
% so you don't have to put the path in for examples!

%% largest_component
load_gaimc_graph('dfs_example')
[Acc p] = largest_component(A); % compute the largest component
xy2 = xy(p,:); labels2 = labels(p); % get component metadata
% draw original graph
subplot(1,2,1); graph_draw(A,xy,'labels',labels); title('Original');
% draw component
subplot(1,2,2); graph_draw(Acc,xy2,'labels',labels2); title('Component');

%% corenums
load_gaimc_graph('cores_example'); % the graph A has three components
corenums(A)

%% sparse_to_csr
A=sparse(6,6); A(1,1)=5; A(1,5)=2; A(2,3)=-1; A(4,1)=1; A(5,6)=1;
[rp ci ai]=sparse_to_csr(A);

%% csr_to_sparse
A=sparse(6,6); A(1,1)=5; A(1,5)=2; A(2,3)=-1; A(4,1)=1; A(5,6)=1;
[rp ci ai]=sparse_to_csr(A);
A2 = csr_to_sparse(rp,ci,ai);

author: ckczzj

repo: CHAN

save_path: github-repos/MATLAB/ckczzj-CHAN


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function save_nii_hdr(hdr, fid)

   if ~exist('hdr','var') || ~exist('fid','var')
      error('Usage: save_nii_hdr(hdr, fid)');
   end

   if ~isequal(hdr.hk.sizeof_hdr,348),
      error('hdr.hk.sizeof_hdr must be 348.');
   end

   if hdr.hist.qform_code == 0 && hdr.hist.sform_code == 0
      hdr.hist.sform_code = 1;
      hdr.hist.srow_x(1) = hdr.dime.pixdim(2);
      hdr.hist.srow_x(2) = 0;
      hdr.hist.srow_x(3) = 0;
      hdr.hist.srow_y(1) = 0;
      hdr.hist.srow_y(2) = hdr.dime.pixdim(3);
      hdr.hist.srow_y(3) = 0;
      hdr.hist.srow_z(1) = 0;
      hdr.hist.srow_z(2) = 0;
      hdr.hist.srow_z(3) = hdr.dime.pixdim(4);
      hdr.hist.srow_x(4) = (1-hdr.hist.originator(1))*hdr.dime.pixdim(2);
      hdr.hist.srow_y(4) = (1-hdr.hist.originator(2))*hdr.dime.pixdim(3);
      hdr.hist.srow_z(4) = (1-hdr.hist.originator(3))*hdr.dime.pixdim(4);
   end

   write_header(hdr, fid);

   return;					% save_nii_hdr

%---------------------------------------------------------------------
function write_header(hdr, fid)

        %  Original header structures
	%  struct dsr				/* dsr = hdr */
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

   header_key(fid, hdr.hk);
   image_dimension(fid, hdr.dime);
   data_history(fid, hdr.hist);

   %  check the file size is 348 bytes
   %
   fbytes = ftell(fid);

   if ~isequal(fbytes,348),
      msg = sprintf('Header size is not 348 bytes.');
      warning(msg);
   end

   return;					% write_header

%---------------------------------------------------------------------
function header_key(fid, hk)

   fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                      /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */

   fwrite(fid, hk.sizeof_hdr(1),    'int32');	% must be 348.

   % data_type = sprintf('%-10s',hk.data_type);	% ensure it is 10 chars from left
   % fwrite(fid, data_type(1:10), 'uchar');
   pad = zeros(1, 10-length(hk.data_type));
   hk.data_type = [hk.data_type  char(pad)];
   fwrite(fid, hk.data_type(1:10), 'uchar');

   % db_name   = sprintf('%-18s', hk.db_name);	% ensure it is 18 chars from left
   % fwrite(fid, db_name(1:18), 'uchar');
   pad = zeros(1, 18-length(hk.db_name));
   hk.db_name = [hk.db_name  char(pad)];
   fwrite(fid, hk.db_name(1:18), 'uchar');

   fwrite(fid, hk.extents(1),       'int32');
   fwrite(fid, hk.session_error(1), 'int16');
   fwrite(fid, hk.regular(1),       'uchar');	% might be uint8

   % fwrite(fid, hk.hkey_un0(1),    'uchar');
   % fwrite(fid, hk.hkey_un0(1),    'uint8');
   fwrite(fid, hk.dim_info(1),      'uchar');

   return;					% header_key

%---------------------------------------------------------------------
function image_dimension(fid, dime)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%			/*
	%				pixdim[] specifies the voxel dimensions:
	%				pixdim[1] - voxel width
	%				pixdim[2] - voxel height
	%				pixdim[3] - interslice distance
	%				pixdim[4] - volume timing, in msec
	%					..etc
	%			*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

   fwrite(fid, dime.dim(1:8),        'int16');
   fwrite(fid, dime.intent_p1(1),  'float32');
   fwrite(fid, dime.intent_p2(1),  'float32');
   fwrite(fid, dime.intent_p3(1),  'float32');
   fwrite(fid, dime.intent_code(1),  'int16');
   fwrite(fid, dime.datatype(1),     'int16');
   fwrite(fid, dime.bitpix(1),       'int16');
   fwrite(fid, dime.slice_start(1),  'int16');
   fwrite(fid, dime.pixdim(1:8),   'float32');
   fwrite(fid, dime.vox_offset(1), 'float32');
   fwrite(fid, dime.scl_slope(1),  'float32');
   fwrite(fid, dime.scl_inter(1),  'float32');
   fwrite(fid, dime.slice_end(1),    'int16');
   fwrite(fid, dime.slice_code(1),   'uchar');
   fwrite(fid, dime.xyzt_units(1),   'uchar');
   fwrite(fid, dime.cal_max(1),    'float32');
   fwrite(fid, dime.cal_min(1),    'float32');
   fwrite(fid, dime.slice_duration(1), 'float32');
   fwrite(fid, dime.toffset(1),    'float32');
   fwrite(fid, dime.glmax(1),        'int32');
   fwrite(fid, dime.glmin(1),        'int32');

   return;					% image_dimension

%---------------------------------------------------------------------
function data_history(fid, hist)

	% Original header structures
	%struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

   % descrip     = sprintf('%-80s', hist.descrip);     % 80 chars from left
   % fwrite(fid, descrip(1:80),    'uchar');
   pad = zeros(1, 80-length(hist.descrip));
   hist.descrip = [hist.descrip  char(pad)];
   fwrite(fid, hist.descrip(1:80), 'uchar');

   % aux_file    = sprintf('%-24s', hist.aux_file);    % 24 chars from left
   % fwrite(fid, aux_file(1:24),   'uchar');
   pad = zeros(1, 24-length(hist.aux_file));
   hist.aux_file = [hist.aux_file  char(pad)];
   fwrite(fid, hist.aux_file(1:24), 'uchar');

   fwrite(fid, hist.qform_code,    'int16');
   fwrite(fid, hist.sform_code,    'int16');
   fwrite(fid, hist.quatern_b,   'float32');
   fwrite(fid, hist.quatern_c,   'float32');
   fwrite(fid, hist.quatern_d,   'float32');
   fwrite(fid, hist.qoffset_x,   'float32');
   fwrite(fid, hist.qoffset_y,   'float32');
   fwrite(fid, hist.qoffset_z,   'float32');
   fwrite(fid, hist.srow_x(1:4), 'float32');
   fwrite(fid, hist.srow_y(1:4), 'float32');
   fwrite(fid, hist.srow_z(1:4), 'float32');

   % intent_name = sprintf('%-16s', hist.intent_name);	% 16 chars from left
   % fwrite(fid, intent_name(1:16),    'uchar');
   pad = zeros(1, 16-length(hist.intent_name));
   hist.intent_name = [hist.intent_name  char(pad)];
   fwrite(fid, hist.intent_name(1:16), 'uchar');

   % magic	= sprintf('%-4s', hist.magic);		% 4 chars from left
   % fwrite(fid, magic(1:4),           'uchar');
   pad = zeros(1, 4-length(hist.magic));
   hist.magic = [hist.magic  char(pad)];
   fwrite(fid, hist.magic(1:4),        'uchar');

   return;					% data_history

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function odf = calcKernelODF(ori,varargin)
% calculate ODF from individuel orientations via kernel density estimation
%
% *calcKernelODF* is one of the core function of the MTEX toolbox.
% It estimates an ODF from a set of individual crystal orientations by
% <EBSD2ODF.html kernel density estimation>.
%
% The function *calcKernelODF* has several options to control the halfwidth
% of the kernel functions, the resolution, etc. Most important the
% estimated ODF is affected by the *halfwidth* of the kernel function.
%
% If the halfwidth is large the estimated ODF is smooth whereas a small
% halfwidth results in a sharp ODF. It depends on your prior information
% about the ODF to choose this parameter right. Look at this
% <EBSD2ODF.html description> for an exhausive discussion.
%
% Syntax
%   odf = calcKernelODF(ori)
%   odf = calcKernelODF(grains.meanOrientation,'weigths',grains.area)
%   odf = calcKernelODF(ebsd.orientations,'halfwidth',5*degree)
%
% Input
%  ori - @orientation
%
% Output
%  odf - @SO3Fun
%
% Options
%  halfwidth  - halfwidth of the kernel function
%  resolution - resolution of the grid where the ODF is approximated
%  kernel     - kernel function (default -- de la Valee Poussin kernel)
%  weights    - list of weights for the orientations
%
% Flags
%  exact      - no approximation to a corser grid
%
% See also
% ebsd_demo EBSD2odf EBSDSimulation_demo EBSD/load EBSD/calcKernel kernel/kernel

% maybe there is nothing to do
if isempty(ori), odf = ODF; return, end

% extract weights
if check_option(varargin,'weights')
  weights = get_option(varargin,'weights');
  varargin = delete_option(varargin,'weights',1);
else
  weights = ones(1,length(ori));
end

% remove nan orientations and weights
weights = weights(~isnan(ori));
ori = subSet(ori,~isnan(ori));

% normalize weights
weights = weights ./ sum(weights(:));

% extract kernel function
psi = SO3DeLaValleePoussinKernel('halfwidth',10*degree,varargin{:});
psi = get_option(varargin,'kernel',psi);
hw = psi.halfwidth;

% if we have to many orientation approximate them on a grid
if length(ori) > 1000 && ~check_option(varargin,'exact')

  [ori,weights] = gridify(ori,'weights',weights,...
    'resolution',max(0.75*degree,hw / 2), varargin{:});

end

% set up exact ODF
odf = unimodalODF(ori,psi,ori.CS,ori.SS,'weights',weights);

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function result = sphere_05_nd ( func, n, xc, r )

%*****************************************************************************80
%
%% SPHERE_05_ND approximates surface integrals on a sphere in ND.
%
%  Integration region:
%
%    R1**2 <= sum ( X(1:N) - XC(1:N) )**2 <= R2**2.
%
%  Discussion:
%
%    A 2*N+2**N points 5-th degree formula is used, Stroud number UN:5-2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F(X), at the N dimensional point
%    X, of the form
%      function value = func ( n, x )
%
%    Input, integer N, the dimension of the space.
%
%    Input, real XC(N), the center of the sphere.
%
%    Input, real R, the radius of the sphere.
%
%    Output, real RESULT, the approximate integral of the function.
%
  x1 = 1.0;
  x2 = 1.0 / sqrt ( n );

  w1 = 1.0 / ( n * ( n + 2 ) );
  w2 = n / ( ( n + 2 ) * 2^n );

  x(1:n) = xc(1:n);

  quad = 0.0;

  for i = 1 : n
    x(i) = xc(i) + r * x1;
    quad = quad + w1 * feval ( func, n, x );
    x(i) = xc(i) - r * x1;
    quad = quad + w1 * feval ( func, n, x );
    x(i) = xc(i);
  end
%
%  For this computation, we need to keep R8VEC_MIRROR_NEXT happy
%  by storing in X the "perturbation" from XC, for which we
%  interested in computing all possible sign variations.
%
  x(1:n) = abs ( r ) * x2;

  while ( 1 )

    quad = quad + w2 * feval ( func, n, xc(1:n) + x(1:n) );

    [ x, done ] = r8vec_mirror_next ( n, x );

    if ( done )
      break
    end

  end

  volume = sphere_area_nd ( n, r );
  result = quad * volume;

  return
end


text:
function A=posmatr(matriks)
%fungsi ini untuk menghitung pada baris ke berapa suatu matriks mencapai
%nilai maksimum
[m,n,o,p]=size(matriks);
for j=1:p
        A(1,j)=posmax(matriks(:,j));
end
A;


text:
%% Header
% Written by Abraham Asfaw
% asfaw@princeton.edu
% Inspired by regular plot version here:
% http://www.mathworks.com/matlabcentral/fileexchange/10656-data-space-to-figure-units-conversion

%% Function
% x_point = x value that you want to be normalized
% y_point = y value that corresponds to the x_point
% axes = axes of the current plot, usually found using get(gca, 'Position')
% xlims_ = x limits of the current plot, usually found using get(gca, 'xlim')
% ylims_ = y limits of the current plot, usually found using get(gca, 'ylim')
% xlog_bool = 1 (any nonzero, really) if the plot is semilogx or loglog, zero otherwise
% ylog_bool = 1 (any nonzero, really) if the plot is semilogy or loglog, zero otherwise

function [nx, ny] = normalize_coordinate(x_point_, y_point_, axes, xlims_, ylims_, xlog_bool, ylog_bool)
    if xlog_bool && ylog_bool, % loglog plots
        x_point = log(x_point_);
        y_point = log(y_point_);
        xlims = log(xlims_);
        ylims = log(ylims_);
    elseif xlog_bool, % semilogx plots
        x_point = log(x_point_);
        y_point = y_point_;
        xlims = log(xlims_);
        ylims = ylims_;
    elseif ylog_bool, % semilogy plots
        x_point = x_point_;
        y_point = log(y_point_);
        xlims = xlims_;
        ylims = log(ylims_);
    else % plot plots
        x_point = x_point_;
        y_point = y_point_;
        xlims = xlims_;
        ylims = ylims_;
    end

    nx = ((x_point-xlims(1))/(xlims(2) - xlims(1)))*axes(3);
    ny = ((y_point-ylims(1))/(ylims(2) - ylims(1)))*axes(4);
    nx = axes(1) + nx;
    ny = axes(2) + ny;
end


text:
function [VI,L]=TriScatteredInterp_nat_near(X,V,XI)

%Natural neighbour interpolation
F=scatteredInterpolant(X,V,'natural');
VI=F(XI);

%Use nearest neighbour for points outside of convex-hull
L=any(isnan(VI),2);
if any(L)
    F=scatteredInterpolant(X,V,'nearest');
    VI(L,:)=F(XI(L,:));
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ d ] = dist2(x0,Y)
% Evaluate the distance between a point X in R^2 and each point Y(1,:)

d = sqrt((x0(1)-Y(:,1)).^2+(x0(2)-Y(:,2)).^2);

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [g]=comp_iwfac(gf,L,a,M)
%COMP_IWFAC  Compute inverse window factorization
%   Usage: g=comp_iwfac(gf,a,M);
%
%   Input parameters:
%         gf    : Factored Window
%         a     : Length of time shift.
%         M     : Number of frequency bands.
%   Output parameters:
%         g     : Window function.
%
%   References: so07-2 st98-8

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: OK
%   REFERENCE: OK

% Calculate the parameters that was not specified
R=prod(size(gf))/L;

N=L/a;
b=L/M;

% The four factorization parameters.
c=gcd(a,M);
p=a/c;
q=M/c;
d=N/q;

gf=reshape(gf,p,q*R,c,d);

% Scale by the sqrt(M) comming from Walnuts representation
gf=gf/sqrt(M);

% fft them
if d>1
  gf=ifft(gf,[],4);
end;

g=zeros(L,R,assert_classname(gf));

% Set up the small matrices
for w=0:R-1
  for s=0:d-1
    for l=0:q-1
      for k=0:p-1
	g((1:c)+mod(k*M-l*a+s*p*M,L),w+1)=reshape(gf(k+1,l+1+q*w,:,s+1),c,1);
      end;
    end;
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Find handle graphics object with user data check.
% Retrieves those handle graphics objects (HGOs) that have the specified
% Tag property and whose UserData property satisfies the given predicate.
%
% Input arguments:
% fcn:
%    a predicate (a function that returns a logical value) to test against
%    the HGO's UserData property
% tag (optional):
%    a string tag to restrict the set of controls to investigate
%
% See also: findobj

% Copyright 2010 Levente Hunyadi
function h = findobjuser(fcn, tag)

validateattributes(fcn, {'function_handle'}, {'scalar'});
if nargin < 2 || isempty(tag)
    tag = '';
else
    validateattributes(tag, {'char'}, {'row'});
end

%hh = get(0, 'ShowHiddenHandles');
%cleanup = onCleanup(@() set(0, 'ShowHiddenHandles', hh));  % restore visibility on exit or exception

if ~isempty(tag)
    % look among all handles (incl. hidden handles) to help findobj locate the object it seeks
    h = findobj(findall(0), '-property', 'UserData', '-and', 'Tag', tag);  % more results if multiple matching HGOs exist
else
    h = findobj(findall(0), '-property', 'UserData');
end
h = unique(h);
try
    for k=1:length(h)
        pred = fcn(get(h(k), 'UserData'));
        if isempty(pred)
            pred = false; end
        f(k) = pred;
    end
    %f = arrayfun(@(handle) fcn(get(handle, 'UserData')), h, 'UniformOutput',false);
catch
    1
end
h = h(f);

repo: REST


text:
function tests = SO2Test
  tests = functiontests(localfunctions);
end

% we will assume that the primitives rotx,trotx, etc. all work

function constructor_test(tc)

    tc.verifyEqual(SO2().double, eye(2,2));

    %% from angle

    tc.verifyEqual(SO2(0).double, eye(2,2), 'AbsTol', 1e-10  );
    tc.verifyEqual(SO2(pi/2).double, rot2(pi/2), 'AbsTol', 1e-10  );
    tc.verifyEqual(SO2(90, 'deg').double, rot2(pi/2), 'AbsTol', 1e-10  );

    %% from R

    tc.verifyEqual(SO2( eye(2,2) ).double, eye(2,2), 'AbsTol', 1e-10  );

    tc.verifyEqual(SO2( rot2(pi/2) ).double, rot2(pi/2), 'AbsTol', 1e-10  );
    tc.verifyEqual(SO2( rot2(pi) ).double, rot2(pi), 'AbsTol', 1e-10  );

    %% from T
    tc.verifyEqual(SO2( trot2(pi/2) ).double, rot2(pi/2), 'AbsTol', 1e-10  );
    tc.verifyEqual(SO2( trot2(pi) ).double, rot2(pi), 'AbsTol', 1e-10  );

    %% R,T
    tc.verifyEqual(SO2( eye(2,2) ).R, eye(2,2), 'AbsTol', 1e-10  );

    tc.verifyEqual(SO2( rot2(pi/2) ).R, rot2(pi/2), 'AbsTol', 1e-10  );

    %% vectorised forms of R
    R = [];
    for theta = [-pi/2 0 pi/2 pi]
        R = cat(3, R, rot2(theta));
    end
    tc.verifyEqual(SO2(R).R, R, 'AbsTol', 1e-10);

    %% copy constructor
    r = SO2(rot2(0.3));
    tc.verifyEqual(SO2(r), r, 'AbsTol', 1e-10);

end

function concat_test(tc)
    x = SO2();
    xx = [x x x x];

    tc.verifyClass(xx, 'SO2');
    tc.verifySize(xx, [1 4]);
end

function primitive_convert_test(tc)
    % char

    s = char( SO2() );

end

function staticconstructors_test(tc)

    %% exponential
    tc.verifyEqual(SO2.exp( skew(0.3) ).R, rot2(0.3), 'AbsTol', 1e-10  );

end

function isa_test(tc)

    verifyTrue(tc, SO2.isa(rot2(0)) );
    verifyTrue(tc, SO2.isa(rot2(0), 'valid') );

    verifyFalse(tc, SO2.isa(1) )
end

function resulttype_test(tc)

    r = SO2();
    verifyClass(tc, r, 'SO2');

    verifyClass(tc, r*r, 'SO2');

    verifyClass(tc, r/r, 'SO2');

    verifyClass(tc, inv(r), 'SO2');

end

function multiply_test(tc)

    vx = [1 0]'; vy = [0 1]';
    r0 = SO2(0);
    r1 = SO2(pi/2);
    r2 = SO2(pi);
    u = SO2();

    %% SO2-SO2 product
    % scalar x scalar

    tc.verifyEqual(r0*u, r0);
    tc.verifyEqual(u*r0, r0);

    % vector x vector
    tc.verifyEqual([r0 r1 r2] * [r2 r0 r1], [r0*r2 r1*r0 r2*r1]);

    % scalar x vector
    tc.verifyEqual(r1 * [r0 r1 r2], [r1*r0 r1*r1 r1*r2]);

    % vector x scalar
    tc.verifyEqual([r0 r1 r2] * r2, [r0*r2 r1*r2 r2*r2]);

    %% SO2-vector product
    % scalar x scalar

    tc.verifyEqual(r1*vx, vy, 'AbsTol', 1e-10);

    % vector x vector
    tc.verifyEqual([r0 r1 r0] * [vy vx vx], [vy vy vx], 'AbsTol', 1e-10);

    % scalar x vector
    tc.verifyEqual(r1 * [vx vy -vx], [vy -vx -vy], 'AbsTol', 1e-10);

    % vector x scalar
    tc.verifyEqual([r0 r1 r2] * vy, [vy -vx -vy], 'AbsTol', 1e-10);
end

function divide_test(tc)

    r0 = SO2(0);
    r1 = SO2(pi/2);
    r2 = SO2(pi);
    u = SO2();

    % scalar / scalar
    % implicity tests inv

    tc.verifyEqual(r1/u, r1);
    tc.verifyEqual(r1/r1, u);

    % vector / vector
    tc.verifyEqual([r0 r1 r2] / [r2 r1 r0], [r0/r2 r1/r1 r2/r0]);

    % vector / scalar
    tc.verifyEqual([r0 r1 r2] / r1, [r0/r1 r1/r1 r2/r1]);
end

function conversions_test(tc)

    T = SO2(pi/2).SE2;
    verifyClass(tc, T, 'SE2');
    tc.verifyEqual(T.T, trot2(pi/2));

    %% Lie stuff
    th = 0.3;
    RR = SO2(th);
    tc.verifyEqual(RR.log, skew(th), 'AbsTol', 1e-10 );

end

function miscellany_test(tc)

    r = SO2( 0.3 );
    tc.verifyEqual(det(r), 1, 'AbsTol', 1e-10  );

    tc.verifyEqual(dim(r), 2);

    verifyEqual( tc, eig(r), eig(r.double) );

    verifyEqual( tc, isSE(r), false );

    verifyClass(tc, r.new, 'SO2');

    verifyClass(tc, SO2.convert(r), 'SO2');
    verifyClass(tc, SO2.convert( rot2(0.3) ), 'SO2');
    z = SO2.convert(r);
    tc.verifyEqual(double(z), double(r));

    z = SO2.convert(rot2(0.3));
    tc.verifyEqual(double(z), rot2(0.3));

    T = r.SE2;
    verifyClass(tc, SO2.convert(T), 'SO2');
end

function display_test(tc)

    R = SO2( 0.3 );

    R.print
    trprint2(R)   % old style syntax

    R.plot
    trplot2(R)   % old style syntax

    R2 = SO2(0.6);
    R.animate
    R.animate(R2)
    tranimate(R2)   % old style syntax
    tranimate(R, R2)   % old style syntax
    tranimate2(R2)   % old style syntax
    tranimate2(R, R2)   % old style syntax
end

author: petercorke

repo: spatialmath-matlab


text:
function TFM = createRotationAboutPoint3d(ROT, point)
%CREATEROTATIONABOUTPOINT3D Rotate about a point using a rotation matrix.
%
%   TFM = createRotationAboutPoint3d(ROT, POINT) Returns the transformation
%   matrix corresponding to a translation(-POINT), rotation with ROT and
%   translation(POINT). Ignores a possible translation in ROT(1:3,4).
%
%   See also
%   transforms3d, transformPoint3d, createRotationOx, createRotationOy,
%   createRotationOz, createRotation3dLineAngle, createRotationVector3d,
%   createRotationVectorPoint3d

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2021-01-31
% Copyright 2021-2022

% Extract only the rotation
ROT = [ROT(1:3,1:3), [0 0 0]'; [0 0 0 1]];

TFM = createTranslation3d(point) * ROT * createTranslation3d(-point);

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function ar=lpczz2ar(zz)
%LPCZZ2AR Convert z-place poles to ar coefficients AR=(ZZ)
% The complex poles must occur in complex conjugate pairs
% but the order is unimportant.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpczz2ar.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p]=size(zz);
ar=zeros(nf,p+1);
for k=1:nf
  ar(k,:)=real(poly(zz(k,:)));
end

author: decouples


text:
function b = r8gb_mxv ( m, n, ml, mu, a, x )

%*****************************************************************************80
%
%% R8GB_MXV multiplies a R8GB matrix times a vector.
%
%  Discussion:
%
%    An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU
%    is assumed to be entirely zero, except for the main diagonal, and
%    entries in the ML nearest subdiagonals, and MU nearest superdiagonals.
%
%    LINPACK and LAPACK "R8GB" storage for such a matrix generally includes
%    room for ML extra superdiagonals, which may be required to store
%    nonzero entries generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%    LINPACK and LAPACK storage of general band matrices requires
%    an extra ML upper diagonals for possible fill in entries during
%    Gauss elimination.  This routine does not access any entries
%    in the fill in diagonals, because it assumes that the matrix
%    has NOT had Gauss elimination applied to it.  If the matrix
%    has been Gauss eliminated, then the routine R8GB_MU must be
%    used instead.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than min(M,N)-1.
%
%    Input, real A(2*ML+MU+1,N), the R8GB matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(M), the product A * x.
%
  b(1:m) = 0.0;

  for i = 1 : m
    jlo = max ( 1, i - ml );
    jhi = min ( n, i + mu );
    for j = jlo : jhi
      b(i) = b(i) + a(i-j+ml+mu+1,j) * x(j);
    end
  end

  return
end


text:
function y = saturateInput(u, umin, umax)

    % SATURATEINPUT saturates the input value.
    %
    % FORMAT: y = saturateInput(u, umin, umax)
    %
    % INPUT:  - umin = [n * 1] min values;
    %         - umax = [n * 1] max values;
    %         - u    = [n * 1] current values;
    %
    % OUTPUT: - y = [n * 1] saturated values;
    %
    % Author(s): Gabriele Nava
    %
    %            all authors are with the Italian Istitute of Technology (IIT)
    %            email: name.surname@iit.it
    %
    % Genoa, Jul 2020
    %

    %% --- Initialization ---
    assert(isequal(size(umin), size(umax)), 'umin and umax must be same size')

    if length(umin) == 1

        y           = u;
        y(y > umax) = umax;
        y(y < umin) = umin;

    else
        assert(length(umin) == length(u), 'input and saturation must have same size');
        y = u;

        for i = 1:length(umin)

            if y(i) > umax(i)

                y(i) = umax(i);

            elseif y(i) < umin(i)

                y(i) = umin(i);
            end
        end
    end
end

author: robotology

repo: whole-body-controllers


text:
function [ x, seed ] = square01_sample ( n, seed )

%*****************************************************************************80
%
%% SQUARE01_SAMPLE samples points in the unit square in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real X(2,N), the points.
%
  m = 2;

  [ x, seed ] = r8mat_uniform_01 ( m, n, seed );

  return
end


text:
function fem1d_adaptive ( )

%*****************************************************************************80
%
%% MAIN is the main program for FEM1D_ADAPTIVE.
%
%  Discussion:
%
%    FEM1D_ADAPTIVE solves a 1D problem using an adaptive finite element method.
%
%    The equation to be treated is:
%
%      -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x)  =  F(x)
%
%    by the finite-element method using piecewise linear basis
%    functions.
%
%    An adaptive method is used to try to reduce the maximum
%    error by refining the mesh in certain places.
%
%    Here U is an unknown scalar function of X defined on the
%    interval [XL,XR], and P, Q and F are given functions of X.
%
%    The values of U at XL and XR are also specified.
%
%    The interval [XL,XR] is "meshed" with N+1 points,
%
%      XN(0) = XL, XN(1) = XL+H, XN(2) = XL+2*H, ..., XN(N) = XR.
%
%    This creates N subintervals, with interval I having endpoints
%    XN(I-1) and XN(I).
%
%    The algorithm tries to guarantee a certain amount
%    of accuracy by examining the current solution, estimating the error
%    in each subinterval, and, if necessary, subdividing one or more
%    subintervals and repeating the calculation.
%
%    We can think of the adaptive part of the algorithm as a refined
%    problem.  The program re-solves the problem on the pair of
%    intervals J and J+1, which extend from node J-1 to node J+1.
%    The values of U that were just computed at nodes J-1 and J+1
%    will be used as the boundary values for this refined problem.
%    The intervals J and J+1 will each be evenly divided into NY
%    smaller subintervals.  This boundary value problem is solved,
%    and the derivatives of the original and refined solutions are
%    then compared to get an estimate of the error.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    real ETA(N).
%    ETA(I) is the error estimate for interval I.  It is computed
%    as the sum of two quantities, one associated with the left
%    and one with the right node of the interval.
%
%    real F(NU).
%    ASSEMBLE stores into F the right hand side of the linear
%    equations.
%    SOLVE replaces those values of F by the solution of the
%    linear equations.
%
%    real FY(M).
%    FY is the right hand side of the linear system of the refined
%    problem.
%
%    real H(N)
%    H(I) is the length of subinterval I.  This code uses
%    equal spacing for all the subintervals.
%
%    real HY(M).
%    HY(I) is the length of subinterval I in the refined problem.
%
%    integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    integer IBCY.
%    IBCY declares the boundary conditions for the refined problem
%    which should always be that the value of U is specified at
%    both the left and right endpoints.  This corresponds to a
%    value of IBCY = 3.
%
%    integer INDX(0:N).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    integer INDY(0:M).
%    INDY(I) records the index of the unknown associated with
%    node I for the refined problem.
%
%    integer JADD(N).
%    JADD(I) is 1 if the error estimates show that interval I
%    should be subdivided.
%
%    integer KOUNT, the number of adaptive steps that have been taken.
%
%    integer M.
%    M is the number of subintervals used in the refined problem.
%    M is equal to NY for computations centered at node 0 or node N,
%    and otherwise, M is equal to 2*NY.
%
%    integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    integer NMAX, the maximum number of unknowns that can be handled.
%
%    integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    integer NODEY(NL,M).
%    NODEY performs the same function for the refined problem that
%    NODE performs for the full problem, recording the node numbers
%    associated with a particular subinterval.
%
%    integer NQUAD
%    The number of quadrature points used in a subinterval.
%
%    integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    integer NUY.
%    The number of unknowns in the refined problem.
%
%    integer NY.
%    NY is the number of subintervals into which a given interval
%    will be subdivided, before solving the refined probelm.
%
%    integer PROBLEM, chooses the problem to be solved.
%    The user must choose this value by setting it in routine GETPRB.
%    * 1, u = x, p = 1, q = 0, f = 0, ibc = 3, ul = 0, ur = 1.
%    The program should find the solution exactly, and the
%    adaptive code should find that there is no reason to
%    subdivide any interval.
%    * 2, u = x*x, p = 1, q = 0, f = -2, ibc = 3, ul = 0, ur = 1.
%    This problem should find the solution exactly, and
%    the adaptive code should again find there is nothing
%    to do.
%    *3, u = sin(pi*x/2), p = 1, q = 0, ibc = 3, f = 0.25*pi*pi*sin(pi*x/2),
%    ul = 0, ur = 1.
%    *4, u = cos(pi*x/2), p = 1, q = 0, ibc = 3, f = 0.25*pi*pi*cos(pi*x/2),
%    ul = 1, ur = 0.
%    *5: u = x**(beta+2)/((beta+2)*(beta+1)), p = 1, q = 1, ibc = 3,
%    f = -x**beta + (x**(beta+2))/((beta+2)*(beta+1)),
%    ul = 0, ur = 1/((beta+2)*(beta+1))
%    (beta must be greater than -2, and not equal to -1)
%    *6: u = atan((x-0.5)/alpha), p = 1, q = 0, ibc = 3,
%    f =  2*alpha*(x-0.5) / (alpha^2 + (x-0.5)^2) ^2,
%    ul = u(0), ur = u(1)
%
%    integer STATUS, reports status of subdivision.
%    0, a new subdivision was carried out.
%    1, no more subdivisions are needed.
%    -1, no more subdivisions can be carried out.
%
%    real TOL.
%    A tolerance that is used to determine whether the estimated
%    error in an interval is so large that it should be subdivided
%    and the problem solved again.
%
%    real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    real WQUAD(NQUAD).
%    WQUAD(I) is the weight associated with the I-th point
%    of an NQUAD point Gaussian quadrature rule.
%
%    real XL.
%    XL is the left endpoint of the interval over which the
%    differential equation is being solved.
%
%    real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    real XQUAD(NQUAD,NMAX), the I-th quadrature point
%    in interval J.
%
%    real XQUADY(NQUAD,NMAY ), the I-th quadrature point
%    in subinterval J of the refined problem.
%
%    real XR.
%    XR is the right endpoint of the interval over which the
%    differential equation is being solved.
%
%    Workspace, double precision XT(0:NMAX), used to compute a new
%    set of nodes.
%
%    real YN(0:M).
%    YN(I) is the location of the I-th node in the refined
%    problem.
%
  nl = 2;
  nmax = 30;
  nquad = 2;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_ADAPTIVE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Solve the two-point boundary value problem:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x)  =  F(x)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'on the interval [0,1], specifying the value\n' );
  fprintf ( 1, 'of U at each endpoint.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'The number of basis functions per element is %d\n', nl );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'The number of quadrature points per element is %d\n', nquad );

  problem = get_problem ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Problem index = %d\n', problem );
  fprintf ( 1, '\n' );

  if ( problem == 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  "Linear" problem:\n' );
    fprintf ( 1, '  (No refinement needed)\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(X) =  X\n' );
    fprintf ( 1, '  P(X) =  1.0\n' );
    fprintf ( 1, '  Q(X) =  0.0\n' );
    fprintf ( 1, '  F(X) =  0.0\n' );
    fprintf ( 1, '  IBC  =  3\n' );
    fprintf ( 1, '  UL   =  0.0\n' );
    fprintf ( 1, '  UR   =  1.0\n' );

  elseif ( problem == 2 )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  "Quadratic" problem:\n' );
    fprintf ( 1, '  (No refinement needed)\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(X) =  X*X\n' );
    fprintf ( 1, '  P(X) =  1.0\n' );
    fprintf ( 1, '  Q(X) =  0.0\n' );
    fprintf ( 1, '  F(X) = -2.0\n' );
    fprintf ( 1, '  IBC  =  3\n' );
    fprintf ( 1, '  UL   =  0.0\n' );
    fprintf ( 1, '  UR   =  1.0\n' );

  elseif ( problem == 3 )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  "SINE" problem:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(X) =  SIN(PI*X/2)\n' );
    fprintf ( 1, '  P(X) =  1.0\n' );
    fprintf ( 1, '  Q(X) =  0.0\n' );
    fprintf ( 1, '  F(X) =  PI*PI*SIN(PI*X/2)/4\n' );
    fprintf ( 1, '  IBC  =  3\n' );
    fprintf ( 1, '  UL   =  0.0\n' );
    fprintf ( 1, '  UR   =  1.0\n' );

  elseif ( problem == 4 )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  "COSINE" problem:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(X) =  COS(PI*X/2)\n' );
    fprintf ( 1, '  P(X) =  1.0\n' );
    fprintf ( 1, '  Q(X) =  0.0\n' );
    fprintf ( 1, '  F(X) =  PI*PI*COS(PI*X/2)/4\n' );
    fprintf ( 1, '  IBC  =  3\n' );
    fprintf ( 1, '  UL   =  0.0\n' );
    fprintf ( 1, '  UR   =  1.0\n' );

  elseif ( problem == 5 )

    beta = get_beta ( );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  "RHEINBOLDT" problem:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(X) =  X**(B+2)/((B+2)*(B+1))\n' );
    fprintf ( 1, '  P(X) =  1.0\n' );
    fprintf ( 1, '  Q(X) =  1.0\n' );
    fprintf ( 1, '  F(X) =  -X**B+(X**B+2))/((B+2)*(B+1))\n' );
    fprintf ( 1, '  IBC  =  3\n' );
    fprintf ( 1, '  UL   =  0.0\n' );
    fprintf ( 1, '  UR   =  1/((B+2)*(B+1))\n' );
    fprintf ( 1, '  B    = %f\n', beta );

  elseif ( problem == 6 )

    alpha = get_alpha ( );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  "ARCTAN" problem:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(X) =  ATAN((X-0.5)/A)\n' );
    fprintf ( 1, '  P(X) =  1.0\n' );
    fprintf ( 1, '  Q(X) =  0.0\n' );
    fprintf ( 1, '  F(X) =  2*A*(X-0.5)/(A^2+(X-0.5)^2)^2\n' );
    fprintf ( 1, '  IBC  =  3\n' );
    fprintf ( 1, '  UL   =  ATAN(-0.5/A)\n' );
    fprintf ( 1, '  UR   =  ATAN( 0.5/A)\n' );
    fprintf ( 1, '  A    = %f\n', alpha );

  end
%
%  Start out with just 4 subintervals.
%
  n = 4;
%
%  Initialize values that define the problem.
%
  [ ibc, tol, ul, ur, xl, xn, xr ] = init ( n );
%
%  Start the iteration counter off at 0.
%
  kount = 0;
%
%  Begin the next iteration.
%
  while ( 1 )

    kount = kount + 1;

    fprintf ( 1, '\n' );
    fprintf ( 1, 'Begin new iteration with %d nodes.\n', n );
    fprintf ( 1, '\n' );
%
%  Solve the regular problem.
%
    [ u, h, nu ] = solvex ( ibc, kount, n, nl, nmax, ul, ur, xn );
%
%  Solve N subproblems to get the error estimators.
%
    eta = solvey ( u, h, n, nu, ul, ur, xn );
%
%  Examine the error estimators, and see how many intervals should
%  be subdivided.
%
    [ n, xn, status ] = subdiv ( eta, kount, n, nmax, tol, xn );

    if ( status ~= 0 )
      break
    end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_ADAPTIVE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ adiag, aleft, arite, f ] = assemble ( h, n, indx, node, nu, nl, ...
  nquad, nmax, ul, ur, wquad, xn, xquad )

%*****************************************************************************80
%
%% ASSEMBLE assembles the global matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real H(N)
%    H(I) is the length of subinterval I.  This code uses
%    equal spacing for all the subintervals.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer INDX(0:N).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Input, integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    Input, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Input, integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    Input, integer NQUAD
%    The number of quadrature points used in a subinterval.
%
%    Input, integer NMAX, the maximum number of unknowns that can be handled.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real WQUAD(NQUAD).
%    WQUAD(I) is the weight associated with the I-th point
%    of an NQUAD point Gaussian quadrature rule.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Input, real XQUAD(NQUAD,NMAX), the I-th quadrature point
%    in interval J.
%
%    Output, real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    Output, real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    Output, real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    Output, real F(NU).
%    ASSEMBLE stores into F the right hand side of the linear
%    equations.
%    SOLVE replaces those values of F by the solution of the
%    linear equations.
%

%
%  Zero out the entries.
%
  f(1:nu) = 0.0;
  aleft(1:nu) = 0.0;
  arite(1:nu) = 0.0;
  adiag(1:nu) = 0.0;
%
%  For each interval,
%
  for ie = 1 : n

    he = h(ie);
    xleft = xn(node(1,ie)+1);
    xrite = xn(node(2,ie)+1);
%
%  For each quadrature point in the interval,
%
    for iq = 1 : nquad

      xquade = xquad(iq,ie);
      wquade = wquad(iq);
%
%  Pick a basis function which defines the equation,
%
      for il = 1 : nl

        ig = node(il,ie);
        iu = indx(ig+1);

        if ( 0 < iu )

          [ phii, phiix ] = phi ( il, xquade, xleft, xrite );
          f(iu) = f(iu) + he * wquade * ff ( xquade ) * phii;
%
%  Take care of boundary conditions specifying the value of U'.
%
          if ( ig == 0 )
            x = 0.0;
            f(iu) = f(iu) - pp ( x ) * ul;
          elseif ( ig == n )
            x = 1.0;
            f(iu) = f(iu) + pp ( x ) * ur;
          end
%
%  Pick a basis function which defines the coefficient
%  being computed.
%
          for jl = 1 : nl

            jg = node(jl,ie);
            ju = indx(jg+1);
            [ phij, phijx ] = phi ( jl, xquade, xleft, xrite );

            aij = he * wquade * ...
               ( pp ( xquade ) * phiix * phijx ...
               + qq ( xquade ) * phii * phij );
%
%  Decide where the coefficient is to be added.
%
            if ( ju <= 0 )

              if ( jg == 0 )
                f(iu) = f(iu) - aij * ul;
              elseif ( jg == n )
                f(iu) = f(iu) - aij * ur;
              end

            elseif ( iu == ju )
              adiag(iu) = adiag(iu) + aij;
            elseif ( ju < iu )
              aleft(iu) = aleft(iu) + aij;
            else
              arite(iu) = arite(iu) + aij;
            end

          end
        end
      end
    end
  end

  return
end
function value = ff ( x )

%*****************************************************************************80
%
%% FF evaluates the function F in the differential equation.
%
%  Discussion:
%
%    This is the function F(X) that appears on the right hand
%    side of the equation:
%
%      -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x)  =  F(x)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%

%
%  Find out which problem we're working on.
%
  problem = get_problem ( );

  if ( problem == 1 )

    value = 0.0;

  elseif ( problem == 2 )

    value = -2.0 * x;

  elseif ( problem == 3 )

    value = 0.25 * pi^2 * sin ( 0.5 * pi * x );

  elseif ( problem == 4 )

    value = 0.25 * pi^2 * cos ( 0.5 * pi * x );

  elseif ( problem == 5 )

    beta = get_beta ( );

    value = - ( x^beta ) + ( x^( beta + 2.0 ) ) ...
      / ( ( beta + 2.0 ) * ( beta + 1.0 ) );

  elseif ( problem == 6 )

    alpha = get_alpha ( );

    value = 2.0 * alpha * ( x - 0.5 ) ...
      / ( alpha^2 + ( x - 0.5 )^2 )^2;

  end

  return
end
function [ h, indx, node, nu, wquad, xquad ] = geometry ( ibc, n, nl, nmax, ...
  nquad, xn )

%*****************************************************************************80
%
%% GEOMETRY sets up some of the geometric information for the problem.
%
%  Discussion:
%
%    Note, however, that the location of the nodes
%    is done outside of this routine, and, in fact, before this
%    routine is called.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    Input, integer NMAX, the maximum number of unknowns that can be handled.
%
%    Input, integer NQUAD
%    The number of quadrature points used in a subinterval.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Output, real H(N)
%    H(I) is the length of subinterval I.  This code uses
%    equal spacing for all the subintervals.
%
%    Output, integer INDX(0:N).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Output, integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    Output, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Output, real WQUAD(NQUAD).
%    WQUAD(I) is the weight associated with the I-th point
%    of an NQUAD point Gaussian quadrature rule.
%
%    Output, real XQUAD(NQUAD,NMAX), the I-th quadrature point
%    in interval J.
%

%
%  Store in NODE the fact that interval I has node I-1
%  as its left endpoint, and node I as its right endpoint.
%
  for i = 1 : n
    node(1,i) = i-1;
    node(2,i) = i;
  end
%
%  For every node that is associated with an unknown, we
%  record the number of the unknown in INDX.
%
  nu = 0;

  for i = 0 : n

    if ( i == 0 & ( ibc == 1 | ibc == 3 ) )
      indx(i+1) = -1;
    elseif ( i == n & ( ibc == 2 | ibc == 3 ) )
      indx(i+1) = -1;
    else
      nu = nu + 1;
      indx(i+1) = nu;
    end

  end
%
%  We compute the width of each interval.
%
  for i = 1 : n
    igl = node(1,i);
    igr = node(2,i);
    h(i) = xn(igr+1) - xn(igl+1);
  end
%
%  We compute the location of the quadrature points in each
%  interval.
%
  for i = 1 : n

    xl = xn(node(1,i)+1);
    xr = xn(node(2,i)+1);

    if ( nquad == 1 )
      xquad(1,i) = 0.5 * ( xl + xr );
    elseif ( nquad == 2 )
      alfa = -0.577350;
      xquad(1,i) = ( ( 1.0 - alfa ) * xl   ...
                   + ( 1.0 + alfa ) * xr ) ...
                   /   2.0;
      alfa = +0.577350;
      xquad(2,i) = ( ( 1.0 - alfa ) * xl   ...
                   + ( 1.0 + alfa ) * xr ) ...
                   /   2.0;
    elseif ( nquad == 3 )
      alfa = -0.774597;
      xquad(1,i) = ( ( 1.0 - alfa ) * xl   ...
                   + ( 1.0 + alfa ) * xr ) ...
                   /   2.0;
      xquad(2,i) = 0.5 * ( xl + xr );
      alfa = +0.774597;
      xquad(3,i) = ( ( 1.0 - alfa ) * xl   ...
                   + ( 1.0 + alfa ) * xr ) ...
                   /   2.0;
    end

  end
%
%  Store the weights for the quadrature rule.
%
  if ( nquad == 1 )
    wquad(1) = 1.0;
  elseif ( nquad == 2 )
    wquad(1) = 0.5;
    wquad(2) = 0.5;
  elseif ( nquad == 3 )
    wquad(1) = 4.0 / 9.0;
    wquad(2) = 5.0 / 18.0;
    wquad(3) = 4.0 / 9.0;
  end

  return
end
function alpha = get_alpha ( )

%*****************************************************************************80
%
%% GET_ALPHA returns the value of ALPHA, for use by problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Output, real ALPHA, the value of ALPHA.
%
  alpha = 0.01;

  return
end
function beta = get_beta ( )

%*****************************************************************************80
%
%% GET_BETA returns the value of BETA, for use by problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2014
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Output, real BETA, the value of BETA.
%
  beta = -0.9;

  return
end
function problem = get_problem ( )

%*****************************************************************************80
%
%% GETPRB returns the value of the current problem number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Output, integer PROBLEM, the index of the problem.
%
  problem = 6;

  return
end
function [ ibc, tol, ul, ur, xl, xn, xr ] = init ( n )

%*****************************************************************************80
%
%% INIT initializes some parameters that define the problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Output, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Output, real TOL.
%    A tolerance that is used to determine whether the estimated
%    error in an interval is so large that it should be subdivided
%    and the problem solved again.
%
%    Output, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Output, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Output, real XL.
%    XL is the left endpoint of the interval over which the
%    differential equation is being solved.
%
%    Output, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Output, real XR.
%    XR is the right endpoint of the interval over which the
%    differential equation is being solved.
%
  tol = 0.01;
%
%  Find out which problem we're working on.
%
  problem = get_problem ( );
%
%  Set the boundary conditions for the problem, and
%  print out its title.
%
  if ( problem == 1 )

    ibc = 3;
    ul = 0.0;
    ur = 1.0;
    xl = 0.0;
    xr = 1.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Exact solution is U = X\n' );

  elseif ( problem == 2 )

    ibc = 3;
    ul = 0.0;
    ur = 1.0;
    xl = 0.0;
    xr = 1.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Exact solution is U = X*X\n' );

  elseif ( problem == 3 )

    ibc = 3;
    ul = 0.0;
    ur = 1.0;
    xl = 0.0;
    xr = 1.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Exact solution is U = SIN(PI*X/2)\n' );

  elseif ( problem == 4 )

    ibc = 3;
    ul = 1.0;
    ur = 0.0;
    xl = 0.0;
    xr = 1.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Exact solution is U = COS(PI*X/2)\n' );

  elseif ( problem == 5 )

    ibc = 3;
    beta = get_beta ( 'DUMMY' );
    ul = 0.0;
    ur = 1.0 / ( ( beta + 2.0 ) * ( beta + 1.0 ) );
    xl = 0.0;
    xr = 1.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Rheinboldt problem\n' );

  elseif ( problem == 6 )

    ibc = 3;
    alpha = get_alpha ( );
    xl = 0.0;
    xr = 1.0;
    ul = uexact ( xl );
    ur = uexact ( xr );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Arctangent problem\n' );

  end
%
%  The nodes are defined here, and not in the geometry routine.
%  This is because each new iteration chooses the location
%  of the new nodes in a special way.
%
  for i = 0 : n
    xn(i+1) = ( ( n - i ) * xl   ...
              + (     i ) * xr ) ...
              / ( n     );
  end

  fprintf ( 1, 'The equation is to be solved for\n' );
  fprintf ( 1, 'X greater than %f\n', xl );
  fprintf ( 1, ' and less than %f\n', xr );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'The boundary conditions are:\n' );
  fprintf ( 1, '\n' );

  if ( ibc == 1 | ibc == 3 )
    fprintf ( 1, '  At X = XL, U = %f\n', ul );
  else
    fprintf ( 1, '  At X = XL, U'' = %f\n', ul );
  end

  if ( ibc == 2 | ibc == 3 )
    fprintf ( 1, '  At X = XR, U = %f\n', ur );
  else
    fprintf ( 1, '  At X = XR, U'' %f\n= ', ur );
  end

  return
end
function output ( f, ibc, indx, n, nu, ul, ur, xn )

%*****************************************************************************80
%
%% OUTPUT prints out the computed solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real F(NU).
%    ASSEMBLE stores into F the right hand side of the linear
%    equations.
%    SOLVE replaces those values of F by the solution of the
%    linear equations.
%
%    Input, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Input, integer INDX(0:N).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Node    X(I)        U(X(I))        Uexact        Error\n' );
  fprintf ( 1, '\n' );

  for i = 0 : n

    if ( i == 0 )

      if ( ibc == 1 | ibc == 3 )
        u = ul;
      else
        u = f(indx(i+1));
      end

    elseif ( i == n )

      if ( ibc == 2 | ibc == 3 )
        u = ur;
      else
        u = f(indx(i+1));
      end

    else

      u = f(indx(i+1));

    end

    uex = uexact ( xn(i+1) );
    error = u - uex;

    fprintf ( 1, '  %4d  %12f  %12f  %12f  %12f\n', i, xn(i+1), u, uex, error );

  end

  return
end
function [ phii, phiix ] = phi ( il, x, xleft, xrite )

%*****************************************************************************80
%
%% PHI evaluates a linear basis function and its derivative.
%
%  Discussion:
%
%    The functions are evaluated at a point X in an interval.  In any
%    interval, there are just two basis functions.  The first
%    basis function is a line which is 1 at the left endpoint
%    and 0 at the right.  The second basis function is 0 at
%    the left endpoint and 1 at the right.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer IL, the local index of the basis function.
%
%    Input, real X, the evaluation point.
%
%    Input, real XLEFT, XRITE, the endpoints of the interval.
%
%    Output, real PHII, PHIIX, the value of the basis function
%    and its derivative.
%
  if ( xleft <= x & x <= xrite )

    if ( il == 1 )
      phii = ( xrite - x ) / ( xrite - xleft );
      phiix = -1.0 / ( xrite - xleft );
    else
      phii = ( x - xleft ) / ( xrite - xleft );
      phiix = 1.0 / ( xrite - xleft );
    end
%
%  If X is outside of the interval, then the basis function
%  is always zero.
%
  else
    phii = 0.0;
    phiix = 0.0;
  end

  return
end
function value = pp ( x )

%*****************************************************************************80
%
%% PP evaluates the function P in the differential equation.
%
%  Discussion:
%
%    The function P(X) occurs in the differential equation:
%
%      -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x)  =  F(x)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of P(X).
%

%
%  Find out which problem we're working on.
%
  problem = get_problem ( );

  if ( problem == 1 )
    value = 1.0;
  elseif ( problem == 2 )
    value = 1.0;
  elseif ( problem == 3 )
    value = 1.0;
  elseif ( problem == 4 )
    value = 1.0;
  elseif ( problem == 5 )
    value = 1.0;
  elseif ( problem == 6 )
    value = 1.0;
  end

  return
end
function prsys ( adiag, aleft, arite, f, nu )

%*****************************************************************************80
%
%% PRSYS prints out the tridiagonal linear system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    Input, real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    Input, real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    Input, real F(NU).
%    ASSEMBLE stores into F the right hand side of the linear
%    equations.
%    SOLVE replaces those values of F by the solution of the
%    linear equations.
%
%    integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Printout of tridiagonal linear system:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Equation  A-Left  A-Diag  A-Rite  RHS\n' );
  fprintf ( 1, '\n' );

  for i = 1 : nu
    if ( i == 1 )
      fprintf ( 1, '%3d                %12f  %12f  %12f\n', ...
        i, adiag(i), arite(i), f(i) );
    elseif ( i < nu )
      fprintf ( 1, '%3d  %12f  %12f  %12f  %12f\n', ...
        i, aleft(i), adiag(i), arite(i), f(i) );
    else
      fprintf ( 1, '%3d  %12f  %12f                %12f\n', ...
        i, aleft(i), adiag(i), f(i) );
    end
  end

  return
end
function value = qq ( x )

%*****************************************************************************80
%
%% QQ evaluates the function Q in the differential equation.
%
%  Discussion:
%
%    The function Q(X) occurs in the differential equation:
%
%      -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x)  =  F(x)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of Q(X).
%

%
%  Find out which problem we're working on.
%
  problem = get_problem ( );

  if ( problem == 1 )
    value = 0.0;
  elseif ( problem == 2 )
    value = 0.0;
  elseif ( problem == 3 )
    value = 0.0;
  elseif ( problem == 4 )
    value = 0.0;
  elseif ( problem == 5 )
    value = 1.0;
  elseif ( problem == 6 )
    value = 0.0;
  end

  return
end
function u = solve ( adiag, aleft, arite, f, nu )

%*****************************************************************************80
%
%% SOLVE solves a tridiagonal matrix system of the form A*x = b.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real ADIAG(NU), ALEFT(NU), ARITE(NU).
%    the diagonal, left and right entries of the equations.
%
%    Input, real F(NU), the right hand side of the linear system.
%
%    Input, INTEGER NU, the number of equations to be solved.
%
%    Output, real U(NU), the solution of the linear system.
%

%
%  Handle the special case of a single equation.
%
  if ( nu == 1 )

    u(1) = f(1) / adiag(1);
%
%  The general case, when NU is greater than 1.
%
  else

    arite(1) = arite(1) / adiag(1);
    for i = 2 : nu-1
      adiag(i) = adiag(i) - aleft(i) * arite(i-1);
      arite(i) = arite(i) / adiag(i);
    end
    adiag(nu) = adiag(nu) - aleft(nu) * arite(nu-1);

    u(1) = f(1) / adiag(1);
    for i = 2 : nu
      u(i) = ( f(i) - aleft(i) * u(i-1) ) / adiag(i);
    end

    for i = nu-1 : -1 : 1
      u(i) = u(i) - arite(i) * u(i+1);
    end

  end

  return
end
function [ u, h, nu ] = solvex ( ibc, kount, n, nl, nmax, ul, ur, xn )

%*****************************************************************************80
%
%% SOLVEX discretizes and solves a differential equation given the nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    Input, integer KOUNT, the number of adaptive steps that have been taken.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    Input, integer NMAX, the maximum number of unknowns that can be handled.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Output, real U(NU), the finite element coefficients of the
%    solution of the differential equation.
%
%    Output, real H(N)
%    H(I) is the length of subinterval I.  This code uses
%    equal spacing for all the subintervals.
%
%    Output, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%  Local parameters:
%
%    Local, real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    Local, real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    Local, real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    Local, integer INDX(0:N).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    Local, integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    Local, integer NQUAD
%    The number of quadrature points used in a subinterval.
%
%    Local, real WQUAD(NQUAD).
%    WQUAD(I) is the weight associated with the I-th point
%    of an NQUAD point Gaussian quadrature rule.
%
%    Local, real XQUAD(NQUAD,NMAX), the I-th quadrature point
%    in interval J.
%
  nquad = 2;
%
%  Given a set of N nodes (where N increases on each iteration),
%  compute the other geometric information.
%
  [ h, indx, node, nu, wquad, xquad ] = geometry ( ibc, n, nl, nmax, nquad, xn );
%
%  Assemble the linear system.
%
  [ adiag, aleft, arite, f ] = assemble ( h, n, indx, node, nu, nl, ...
    nquad, nmax, ul, ur, wquad, xn, xquad );
%
%  Print out the linear system, just once.
%
  if ( kount == 1 )
    prsys ( adiag, aleft, arite, f, nu );
  end
%
%  Solve the linear system.
%
  u = solve ( adiag, aleft, arite, f, nu );
%
%  Print out the solution.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Basic solution\n' );

  output ( u, ibc, indx, n, nu, ul, ur, xn );

  return
end
function eta = solvey ( u, h, n, nu, ul, ur, xn )

%*****************************************************************************80
%
%% SOLVEY computes error estimators for a finite element solution.
%
%  Discussion:
%
%    SOLVEY accepts information about the solution of a finite element
%    problem on a grid of nodes with coordinates XN.  It then starts
%    at node 0, and for each node, computes two "error estimators",
%    one for the left, and one for the right interval associated with the
%    node.  These estimators are found by solving a finite element problem
%    over the two intervals, using the known values of the original
%    solution as boundary data, and using a mesh that is "slightly"
%    refined over the original one.
%
%    Note that the computations at the 0-th and N-th nodes only involve
%    a single interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real U(NU), the finite element representation of the
%    solution of the current problem.
%
%    Input, real H(N)
%    H(I) is the length of subinterval I.  This code uses
%    equal spacing for all the subintervals.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    Input, real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    Input, real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Output, real ETA(N).
%    ETA(I) is the error estimate for interval I.  It is computed
%    as the sum of two quantities, one associated with the left
%    and one with the right node of the interval.
%
  nl = 2;
  ny = 2;
  nquad = 2;
  nmay = 2 * ny;
%
%  Initialize the error estimators to zero.
%
  eta(1:n) = 0.0;
%
%  Set the boundary conditions for each subproblem to be
%  known values of U at the left and right.
%
%
%  For each node, subdivide its left and right hand intervals
%  into NY subintervals.
%
%  Set up and solve the differential equation again on this
%  smaller region.
%
%  The 0-th and N-th nodes are special cases.
%
  ibcy = 3;

  for j = 0 : n

    if ( j == 0 )
      m = ny;
      jlo = j;
      jmid = j + 1;
      jhi = j + 1;
    elseif ( j == n )
      m = ny;
      jlo = j - 1;
      jmid = j;
      jhi = j;
    else
      m = 2 * ny;
      jlo = j - 1;
      jmid = j;
      jhi = j + 1;
    end
%
%  Set the location of the nodes in the subintervals.
%
    yl = xn(jlo+1);
    ym = xn(jmid+1);
    yr = xn(jhi+1);

    for i = 0 : ny
      yn(i+1) = ( ( ny - i ) * yl   ...
                + (      i ) * ym ) ...
                / ( ny     );
    end

    for i = ny+1 : m
      yn(i+1) = ( ( m - i      ) * ym   ...
                + (     i - ny ) * yr ) ...
                / ( m -     ny );
    end
%
%  Set up the geometry of the sub-problem.
%
    [ hy, indy, nodey, nuy, wquad, xquady ] = geometry ( ibcy, m, nl, ...
      nmay, nquad, yn );
%
%  Set the boundary values for the sub-problem.
%
    if ( j <= 1 )
      uly = ul;
    else
      uly = u(j-1);
    end

    if ( n - 1 <= j )
      ury = ur;
    else
      ury = u(j+1);
    end
%
%  Assemble the matrix for the sub-problem.
%
    [ adiag, aleft, arite, fy ] =  assemble ( hy, m, indy, nodey, nuy, nl, ...
      nquad, nmay, uly, ury, wquad, yn, xquady );
%
%  Solve the system.
%
    uy = solve ( adiag, aleft, arite, fy, nuy );
%
%  Compute the weighted sum of the squares of the differences
%  of the original computed slope and the refined computed slopes.
%
%  Calculation for left interval.
%
    if ( 1 <= j )

      if ( j <= 1 )
        uleft = ul;
        urite = u(1);
      elseif ( j == n )
        uleft = u(j-1);
        urite = ur;
      else
        uleft = u(j-1);
        urite = u(j);
      end

      uprime = ( urite - uleft ) / h(j);

      total = 0.0;
      for i = 1 : ny

        yl = yn(i-1+1);
        yr = yn(i+1);

        if ( i == 1 )
          vlval = uly;
          vrval = uy(i);
        elseif ( i == m )
          vlval = uy(i-1);
          vrval = ury;
        else
          vlval = uy(i-1);
          vrval = uy(i);
        end

        vprime = ( vrval - vlval ) / hy(i);

        ulval = ( real ( ny - i + 1 ) * uleft   ...
                + real (      i - 1 ) * urite ) ...
                / real ( ny         );

        urval = ( real ( ny - i ) * uleft   ...
                + real (      i ) * urite ) ...
                / real ( ny     );
%
%  Compute the integral of
%
%    p(x)*(u'(x)-v'(x))^2 + q(x)*(u(x)-v(x))^2
%
        for k = 1 : nquad

          y  =  xquady(k,i);

          uval = ( ( yl - y      ) * urval   ...
                 + (      y - yr ) * ulval ) ...
                 / ( yl     - yr );

          vval = ( ( yl - y      ) * vrval   ...
                 + (      y - yr ) * vlval ) ...
                 / ( yl     - yr );

          total = total + 0.5 * wquad(k) * hy(i) * ...
            ( pp ( y ) * ( uprime - vprime )^2 ...
            + qq ( y ) * ( uval - vval )^2 );

        end

      end

      eta(j) = eta(j) + 0.5 * sqrt ( total );

    end
%
%  Calculation for right interval.
%
    if ( j <= n - 1 )

      if ( j == 0 )
        uleft = ul;
        urite = u(j+1);
      elseif ( n - 1 <= j )
        uleft = u(j);
        urite = ur;
      else
        uleft = u(j);
        urite = u(j+1);
      end

      uprime = ( urite - uleft ) / h(j+1);

      total = 0.0;
      for i = m+1-ny : m

        yl = yn(i);
        yr = yn(i+1);

        if ( i == 1 )
          vlval = uly;
          vrval = uy(i);
        elseif ( i == m )
          vlval = uy(i-1);
          vrval = ury;
        else
          vlval = uy(i-1);
          vrval = uy(i);
        end

        vprime = ( vrval - vlval ) / hy(i);

        ulval = ( (      m - i + 1 ) * uleft   ...
                + ( ny - m + i - 1 ) * urite ) ...
                / ( ny             );

        urval = ( (      m - i ) * uleft   ...
                + ( ny - m + i ) * urite ) ...
                / ( ny         );
%
%  Compute the integral of
%
%    p(x)*(u'(x)-v'(x))^2 + q(x)*(u(x)-v(x))^2
%
        for k = 1 : nquad

          y  =  xquady(k,i);

          uval = ( ( yl - y      ) * urval   ...
                 + (      y - yr ) * ulval ) ...
                 / ( yl     - yr );

          vval = ( ( yl - y      ) * vrval   ...
                 + (      y - yr ) * vlval ) ...
                 / ( yl     - yr );

          total = total + 0.5 * wquad(k) * hy(i) * ...
            ( pp ( y ) * ( uprime - vprime )^2 ...
            + qq ( y ) * ( uval - vval )^2 );

        end

      end

      eta(j+1) = eta(j+1) + 0.5 * sqrt ( total );

    end

  end
%
%  Print out the error estimators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ETA\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %12f\n', eta(j) );
  end

  return
end
function [ n, xn, status ] = subdiv ( eta, kount, n, nmax, tol, xn )

%*****************************************************************************80
%
%% SUBDIV decides which intervals should be subdivided.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real ETA(N).
%    ETA(I) is the error estimate for interval I.  It is computed
%    as the sum of two quantities, one associated with the left
%    and one with the right node of the interval.
%
%    Input, integer KOUNT, the number of adaptive steps that have been taken.
%
%    Input, integer N
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    Input, integer NMAX, the maximum number of unknowns that can be handled.
%
%    Input, real TOL.
%    A tolerance that is used to determine whether the estimated
%    error in an interval is so large that it should be subdivided
%    and the problem solved again.
%
%    Input, real XN(0:N).
%    XN(I) is the location of the I-th node.  XN(0) is XL,
%    and XN(N) is XR.
%
%    Output, integer N, the updated number of nodes.
%
%    Output, real XN(0:N), the updated nodes.
%
%    Output, integer STATUS, reports status of subdivision.
%    0, a new subdivision was carried out.
%    1, no more subdivisions are needed.
%    -1, no more subdivisions can be carried out.
%
%  Local Parameters:
%
%    Local, integer JADD(N).
%    JADD(I) is 1 if the error estimates show that interval I
%    should be subdivided.
%
  status = 0;
%
%  Add up the ETA's, and get their average.
%
  ave = sum ( eta(1:n) ) / n;
%
%  Look for intervals whose ETA value is relatively large,
%  and note in JADD that these intervals should be subdivided.
%
  k = 0;
  temp = max ( 1.2 * ave + 0.00001, tol^2 / n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Tolerance = %f\n', temp );
  fprintf ( 1, '\n' );

  for j = 1 : n

    if ( temp < eta(j) )
      k = k + 1;
      jadd(j) = 1;
      fprintf ( 1, 'Subdivide interval %d\n', j );
    else
      jadd(j) = 0;
    end

  end
%
%  If no subdivisions needed, we're done.
%
  if ( k <= 0 )
    fprintf ( 1, 'Success on step %d\n', kount );
    status = 1;
    return
  end
%
%  See if we're about to go over our limit.
%
  if ( nmax < n + k )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'The iterations did not reach their goal.\n' );
    fprintf ( 1, 'The next value of N is %d\n', n + k );
    fprintf ( 1, 'which exceeds NMAX = %d\n', nmax );
    status = -1;
    return
  end
%
%  Insert new nodes where needed.
%
  k = 0;
  xt(1) = xn(1);
  for j = 1 : n

    if ( 0 < jadd(j) )
      xt(j+k+1) = 0.5 * ( xn(j+1) + xn(j) );
      k = k + 1;
    end

    xt(j+k+1) = xn(j+1);

  end
%
%  Update the value of N, and copy the new nodes into XN.
%
  n = n + k;

  xn(1:n+1) = xt(1:n+1);

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function value = uexact ( x )

%*****************************************************************************80
%
%% UEXACT returns the value of the exact solution at any point X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of the exact solution at X.
%

%
%  Find out which problem we're working on.
%
  problem = get_problem ( );

  if ( problem == 1 )
    value = x;
  elseif ( problem == 2 )
    value = x^2;
  elseif ( problem == 3 )
    value = sin ( pi * x / 2.0 );
  elseif ( problem == 4 )
    value = cos ( pi * x / 2.0 );
  elseif ( problem == 5 )
    beta = get_beta ( );
    value = ( x^( beta + 2.0 ) ) ...
      / ( ( beta + 2.0 ) * ( beta + 1.0 ) );
  elseif ( problem == 6 )
    alpha = get_alpha ( );
    value = atan ( ( x - 0.5 ) / alpha );
  else
    value = 0.0;
  end

  return
end


text:
function [signalMatrix, signalSpikes] = viewSpikeSummary(signalMatrix,signalSpikes,varargin)
	% gives several summary plots for a set of signals
	% biafra ahanonu
	% started: 2013.10.28
	% inputs
		% signalMatrix - 2D matrix of [signals frames] that contains the raw analog signal
		% signalSpikes - 2D matrix of [signals frames] that contains 0/1 indicating whether a peak occurred at that point
	% outputs
		% signalMatrix - 2D matrix of [signals frames] that contains the raw analog signal
		% signalSpikes - 2D matrix of [signals frames] that contains 0/1 indicating whether a peak occurred at that point

	[signalMatrix, signalSpikes] = ciapkg.view.viewSpikeSummary(signalMatrix,signalSpikes,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah