text:
function [Jul1,Jul2]=TT2TDB(Jul1,Jul2,deltaTTUT1,clockLoc)
%%TT2TDB Convert from terrestrial time (TT) to barycentric dynamical time
%        (TDB) to an accuracy of nanoseconds (if deltaT is accurate)
%        using the routines from the International Astronomical Union's
%        library that do not require external ephemeris data.
%
%INPUTS: Jul1,Jul2 Two parts of a Julian date given in TT. The units of
%                  the date are days. The full date is the sum of both
%                  terms. The date is broken into two parts to provide
%                  more bits of precision. It does not matter how the date
%                  is split.
%       deltaTTUT1 An optional parameter specifying the offset between TT
%                  and UT1 in seconds. If this parameter is omitted or an
%                  empty matrix is passed, then the value of the function
%                  getEOP will be used.
%         clockLoc An optional 3X1 vector specifying the location of the
%                  clock in the Terrestrial Intermediate Reference System
%                  (TIRS), though it would not make much of a difference
%                  if the International Terrestrial Reference System
%                  (ITRS) were used. The units are meters. Due to
%                  relativistic effects, clocks that are synchronized with
%                  respect to TT are not synchronized with respect to TDB.
%                  If this parameter is omitted, then a clock at the
%                  center of the Earth is used.
%
%OUTPUTS:Jul1,Jul2 Two parts of a Julian date given in TDB.
%
%This function relies on a number of functions in the International
%Astronomical Union's Standards of Fundamental Astronomy library. The
%implementation is essentially the same as TT2TCB except a final
%conversion step is omitted.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[Jul1,Jul2]=TT2TDB(Jul1,Jul2,deltaTTUT1,clockLoc);
%or
%[Jul1,Jul2]=TT2TDB(Jul1,Jul2);
%
%Many temporal coordinate systems standards are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report,
%    Naval Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016,
%    173 pages.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%{
Author: Deepti Ghadiyaram

Description: Given an input RGB image, this method transforms the input
into HSI color space.
%}

function hsi= convertRGBToHSI(rgb)
    hsi = (colorspace('RGB->HSI',rgb));
end

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function A=storageclass(lwl,hwl,sto_class)
%STORAGECLASS berfungsi untuk membagi kapasitas tampungan dalam beberapa
%kelas
%lwl adalah low water level operation (m3)
%hwl adalah high water level operation (m3)
%sto_class adalah jumlah storage class
%berbeda dengan inflowclass, nilai kelas disini adalah nilai batas
incr=(hwl-lwl)/(sto_class-1);
%penentuan nilai maksimum per kelas
for a=0:(sto_class-1)
    A(a+1,:)=lwl+incr*a;
end
A;

author: Sable

repo: mcbench-benchmarks


text:
function inds = fernsInds( data, fids, thrs )
% Compute indices for each input by each fern.
%
% USAGE
%  inds = fernsInds( data, fids, thrs )
%
% INPUTS
%  data     - [NxF] N length F binary feature vectors
%  fids     - [MxS] feature ids for each fern for each depth
%  thrs     - [MxS] threshold corresponding to each fid
%
% OUTPUTS
%  inds     - [NxM] computed indices for each input by each fern
%
% EXAMPLE
%
% See also fernsClfTrain, fernsClfApply
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

inds = fernsInds1( data, fids, thrs );

%%% OLD MATLAB CODE -- NOW IN MEX
% [M,S]=size(fids); N=size(data,1);
% inds = zeros(N,M,'uint32');
% for n=1:N
%   for m=1:M
%     for s=1:S
%       inds(n,m)=inds(n,m)*2;
%       if( data(n,fids(m,s))<thrs(m,s) )
%         inds(n,m)=inds(n,m)+1;
%       end
%     end
%   end
% end
% inds=inds+1;
% end

author: zouchuhang

repo: LayoutNet


text:
function [Bx, By, Bz] = igrf(time, latitude, longitude, altitude, coord)

% IGRF Earth's magnetic field from IGRF model.
%
% Usage: [BX, BY, BZ] = IGRF(TIME, LATITUDE, LONGITUDE, ALTITUDE, COORD)
%     or [BX, BY, BZ] = IGRF(COEFS, LATITUDE, LONGITUDE, ALTITUDE, COORD)
%     or B = IGRF(TIME, LATITUDE, LONGITUDE, ALTITUDE, COORD)
%     or B = IGRF(COEFS, LATITUDE, LONGITUDE, ALTITUDE, COORD)
%
% Calculates the components of the Earth's magnetic field using the
% International Geomagnetic Reference Field (IGRF) model. The inputs for
% the position can be scalars or vectors (in the latter case each should
% have the same number of elements or be a scalar), but TIME must be a
% scalar.
%
% When all the coordinate inputs are scalars, the function can be run more
% efficiently by providing the proper IGRF coefficient vector for a given
% time rather than the time itself. This mode is useful when making
% multiple calls to the function while keeping the time the same (meaning
% the coefficients will be the same for each run) as loading the
% coefficients can be the most time-consuming part of the function. The
% coefficient vector can be easily loaded using the function LOADIGRFCOEFS.
% This mode is assumed when all the coordinate inputs are scalars and the
% first input is a vector. In this case, the coefficient vector should be
% formatted as (LOADIGRFCOEFS provides this):
%
%   [g(n=1,m=0) g(n=1,m=1) h(n=1,m=1) g(n=2,m=0) g(n=2,m=1) h(n=2,m=1) ...]
%
% Regardless of the size of the inputs, the outputs will be column vectors.
% If only one output is requested, B = [BX(:), BY(:), BZ(:)] is output.
% Note that the other parameters the IGRF gives can be computed from BX,
% BY, and BZ as:
%
%   Horizonal intensity: hypot(BX, BY) (i.e., sqrt(BX.^2 + BY.^2) )
%   Total intensity: hypot(BX, hypot(BY, BZ))
%   Declination: atan2(BY, BX)
%   Inclination: atan(BZ./hypot(BX, BY))
%
% This function relies on having the file igrfcoefs.mat in the MATLAB
% path to function properly when a time is input. If this file cannot be
% found, this function will try to create it by calling GETIGRFCOEFS.
%
% The IGRF is a spherical harmonic expansion of the Earth's internal
% magnetic field. Currently, the IGRF model is valid between the years 1900
% and 2015. See the health warning for the IGRF model here:
% http://www.ngdc.noaa.gov/IAGA/vmod/igrfhw.html
%
% Reference:
% International Association of Geomagnetism and Aeronomy, Working Group
% V-MOD (2010), International Geomagnetic Reference Field: the eleventh
% generation, _Geophys. J. Int._, _183_(3), 1216-1230,
% doi:10.1111/j.1365-246X.2010.04804.x.
%
% Inputs:
%   -TIME: Time to get the magnetic field values either in MATLAB serial
%   date number format or a string that can be converted into MATLAB serial
%   date number format using DATENUM with no format specified (see
%   documentation of DATENUM for more information).
%   -COEFS: Instead of inputting a time, you can simply specify the proper
%   coefficients for the time you want by inputting in the first argument
%   the proper coefficient vector from igrfcoefs.mat.
%   -LATITUDE: Geocentric or geodetic latitude in degrees.
%   -LONGITUDE: Geocentric or geodetic longitude in degrees.
%   -ALTITUDE: For geodetic coordiates, the height in km above the Earth's
%   surface. For geocentric coordiates, the radius in km from the center of
%   the Earth.
%   -COORD: String specifying the coordinate system to use. Either
%   'geocentric' or 'geodetic' (optional, default is geodetic). Note that
%   only geodetic coordinates have been verified.
%
% Outputs:
%   -BX: Northward component of the magnetic field in nanoteslas (nT).
%   -BY: Eastward component of the magnetic field in nT.
%   -BZ: Downward component of the magnetic field in nT.
%   -B: [BX(:), BY(:), BZ(:)].
%
% See also: LOADIGRFCOEFS, GETIGRFCOEFS, IGRFLINE, DATENUM, IGRF11MAGM.

% Run IGRFS if all position inputs are scalars.
if isscalar(latitude) && isscalar(longitude) && isscalar(altitude)
    if nargin < 5
        [Bx, By, Bz] = igrfs(time, latitude, longitude, altitude);
    else
        [Bx, By, Bz] = igrfs(time, latitude, longitude, altitude, coord);
    end
% Otherwise run IGRFV.
else
    if nargin < 5
        [Bx, By, Bz] = igrfv(time, latitude, longitude, altitude);
    else
        [Bx, By, Bz] = igrfv(time, latitude, longitude, altitude, coord);
    end
end

if nargout <= 1
    Bx = [Bx(:), By(:), Bz(:)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%        IGRF vector function.        %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Bx, By, Bz] = igrfv(time, latitude, longitude, altitude, coord)

% Fundamental constant.
Rearth_km = 6371.2;

%%% CHECK INPUT VALIDITY %%%
% Convert time to a datenumber if it is a string.
if ischar(time)
    time = datenum(time);
end

% Make sure time has only one element.
if numel(time) > 1
    error('igrf:timeInputInvalid', ['The input TIME can only have one ' ...
        'element.']);
end

% Check that the inputs all have either one or the same number of elements.
numlat = numel(latitude);
numlon = numel(longitude);
numalt = numel(altitude);
if numlat > 1
    if numlon == 1
        longitude = repmat(longitude, size(latitude));
    end
    if numalt == 1
        altitude = repmat(altitude, size(latitude));
    end
elseif numlon > 1
    latitude = repmat(latitude, size(longitude));
    if numalt == 1
        altitude = repmat(altitude, size(latitude));
    end
elseif numalt > 1
    latitude = repmat(latitude, size(altitude));
    longitude = repmat(longitude, size(altitude));
end
numlat = numel(latitude);
numlon = numel(longitude);
numalt = numel(altitude);
if numlat ~= numlon || numlat ~= numalt || numlon ~= numalt
    error('igrf:inputNotSameSize', ['The input coordinates must have ' ...
        'the same number of elements.']);
end

%%% SPHERICAL COORDINATE CONVERSION %%%
% Convert the latitude, longitude, and altitude coordinates input into
% spherical coordinates r (radius), theta (inclination angle from +z axis),
% and phi (azimuth angle from +x axis). Also, make the coordinates go down
% the rows.
% We want cos(theta) and sin(theta) rather than theta itself.
costheta = cos((90 - latitude(:))*pi/180);
sintheta = sin((90 - latitude(:))*pi/180);

% Convert from geodetic coordinates to geocentric coordinates if necessary.
% This method was adapted from igrf11syn, which was a conversion of the
% IGRF subroutine written in FORTRAN.
if nargin < 5 || isempty(coord) || strcmpi(coord, 'geodetic') || ...
        strcmpi(coord, 'geod') || strcmpi(coord, 'gd')
    a = 6378.137; f = 1/298.257223563; b = a*(1 - f);
    rho = hypot(a*sintheta, b*costheta);
    r = sqrt( altitude(:).^2 + 2*altitude(:).*rho + ...
        (a^4*sintheta.^2 + b^4*costheta.^2) ./ rho.^2 );
    cd = (altitude(:) + rho) ./ r;
    sd = (a^2 - b^2) ./ rho .* costheta.*sintheta./r;
    oldcos = costheta;
    costheta = costheta.*cd - sintheta.*sd;
    sintheta = sintheta.*cd + oldcos.*sd;
elseif strcmpi(coord, 'geocentric') || strcmpi(coord, 'geoc') || ...
        strcmpi(coord, 'gc')
    r = altitude(:);
    cd = 1;
    sd = 0;
else
    error('igrf:coordInputInvalid', ['Unrecognized command ' coord ...
        ' for COORD input.']);
end

% Special case when sin(theta) = 0.
sintheta0 = sintheta == 0;
anysintheta0 = any(sintheta0);
anysinthetanot0 = any(~sintheta0);

% Convert longitude to radians.
phi = longitude(:)*pi/180;

%%% GET PROPER IGRF COEFFICIENTS %%%
[g, h] = loadigrfcoefs(time);
nmax = size(g, 1);

% We need cos(m*phi) and sin(m*phi) multiple times, so precalculate into a
% matrix here:
cosphi = cos(bsxfun(@times, 0:nmax, phi));
sinphi = sin(bsxfun(@times, 0:nmax, phi));

%%% BEGIN MAGNETIC FIELD CALCULATION %%%
% Initialize variables used in for loop below.
Br = zeros(size(r));
Bt = zeros(size(r));
Bp = zeros(size(r));
lastP = 1;
lastdP_1 = 0;
lastdP_2 = 0;

% Sum for each n value.
for n = 1 : nmax

    m = 0 : n;

    % Calculate legendre values. The output of the function has each m
    % value going down the rows, but since m goes along the columns
    % (coordinates go down the rows, remember?), permute it.
    P = legendre(n, costheta, 'sch').';

    % We also need the derivative of the legendre with respect to theta. It
    % is given by a recursive function of both the previous legendre values
    % as well as the previous derivatives. Functionally, it is:
    % dP(0, 0) = 0, dP(1, 1) = cos(theta)
    % dP(n, n) = sqrt(1 - 1/(2n))*(sin(theta)*dP(n-1, n-1) +
    %     cos(theta)*P(n-1, n-1))
    % dP(n, m) = (2n - 1)/sqrt(n^2 - m^2)*(cos(theta)*dP(n-1, m) -
    %     sin(theta)*P(n-1, m)) - sqrt(((n-1)^2 - m^2)/(n^2 - m^2))*
    %     dP(n-2, m)
    dP = [bsxfun(@minus, bsxfun(@times, ...
        (2*n - 1)./sqrt(n^2 - m(1:end-1).^2), ...
        bsxfun(@times, costheta, lastdP_1) - bsxfun(@times, sintheta, ...
        lastP)), bsxfun(@times, sqrt(((n - 1)^2 - m(1:end-1).^2)./...
        (n^2 - m(1:end-1).^2)), lastdP_2)), zeros(size(costheta))];
    if n > 1
        dP(:, end) = sqrt(1 - 1/(2*n))*...
            (sintheta*lastdP_1(end) + costheta*lastP(end));
        lastdP_2 = [lastdP_1 zeros(size(costheta))];
    else
        dP(:, end) = costheta;
        lastdP_2 = lastdP_1;
    end
    lastP = P;
    lastdP_1 = dP;

    % Multiply coefficients by proper longitude trigonemetric term.
    gcos = bsxfun(@times, g(n, m + 1), cosphi(:, m + 1));
    gsin = bsxfun(@times, g(n, m + 1), sinphi(:, m + 1));
    hcos = bsxfun(@times, h(n, m + 1), cosphi(:, m + 1));
    hsin = bsxfun(@times, h(n, m + 1), sinphi(:, m + 1));

    % Calculate the magnetic field components as a running sum. Find
    % explicit expressions for these in Global Earth Physics: a Handbook of
    % Physical Constants by Thomas J. Aherns (1995), pg. 49. Link:
    % http://books.google.com/books?id=aqjU_NHyre4C&lpg=PP1&dq=Global%20
    % earth%20physics%3A%20a%20handbook%20of%20physical%20constants&pg=PA49
    % #v=onepage&q&f=false
    % (except equation 6 is missing a required 1/sin(theta) and m; correct
    % equations on page 5 (equations 3a-3c) of:
    % http://hanspeterschaub.info/Papers/UnderGradStudents/
    % MagneticField.pdf)
    a_r = (Rearth_km./r).^(n + 2);
    Br = Br + a_r.*(n+1).*sum((gcos + hsin).*P, 2);
    Bt = Bt - a_r.*sum((gcos + hsin).*dP, 2);
    % Different case when sin(theta) == 0 for phi component.
    if anysinthetanot0
        Bp(~sintheta0) = Bp(~sintheta0) - 1./sintheta(~sintheta0).*...
            a_r(~sintheta0).*sum(bsxfun(@times, m, ...
            (-gsin(~sintheta0, :) + hcos(~sintheta0, :)).*...
            P(~sintheta0, :)), 2);
    end
    if anysintheta0
        Bp(sintheta0) = Bp(sintheta0) - costheta(sintheta0).*...
            a_r(sintheta0).*sum((-gsin(sintheta0, :) ...
            + hcos(sintheta0, :)).*dP(sintheta0, :), 2);
    end

end

% Convert from spherical to (x,y,z) = (North,East,Down).
Bx = -Bt;
By = Bp;
Bz = -Br;

% Convert back to geodetic coordinates if necessary.
Bx_old = Bx;
Bx = Bx.*cd + Bz.*sd;
Bz = Bz.*cd - Bx_old.*sd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%        IGRF scalar function.        %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Bx, By, Bz] = igrfs(time, latitude, longitude, altitude, coord)

% Fundamental constant.
Rearth_km = 6371.2;

%%% CHECK INPUT VALIDITY %%%
% Convert time to a datenumber if it is a string.
if ischar(time)
    time = datenum(time);
end

% Check that the input coordinates are scalars.
if ~isscalar(latitude) || ~isscalar(longitude) || ~isscalar(altitude)
    error('igrf1:inputNotScalar', ...
        'The input coordinates must be scalars.');
end

%%% SPHERICAL COORDINATE CONVERSION %%%
% Convert the latitude, longitude, and altitude coordinates input into
% spherical coordinates r (radius), theta (inclination angle from +z axis),
% and phi (azimuth angle from +x axis).
% We want cos(theta) and sin(theta) rather than theta itself.
costheta = cos((90 - latitude)*pi/180);
sintheta = sin((90 - latitude)*pi/180);

% Convert from geodetic coordinates to geocentric coordinates if necessary.
% This method was adapted from igrf11syn, which was a conversion of the
% IGRF subroutine written in FORTRAN.
if nargin < 5 || isempty(coord) || strcmpi(coord, 'geodetic') || ...
        strcmpi(coord, 'geod') || strcmpi(coord, 'gd')
    a = 6378.137; f = 1/298.257223563; b = a*(1 - f);
    rho = hypot(a*sintheta, b*costheta);
    r = sqrt( altitude.^2 + 2*altitude.*rho + ...
        (a^4*sintheta.^2 + b^4*costheta.^2) ./ rho.^2 );
    cd = (altitude + rho) ./ r;
    sd = (a^2 - b^2) ./ rho .* costheta.*sintheta./r;
    oldcos = costheta;
    costheta = costheta.*cd - sintheta.*sd;
    sintheta = sintheta.*cd + oldcos.*sd;
elseif strcmpi(coord, 'geocentric') || strcmpi(coord, 'geoc') || ...
        strcmpi(coord, 'gc')
    r = altitude;
    cd = 1;
    sd = 0;
else
    error('igrf:coordInputInvalid', ['Unrecognized command ' coord ...
        ' for COORD input.']);
end

% Convert longitude to radians.
phi = longitude*pi/180;

%%% GET PROPER IGRF COEFFICIENTS %%%
if isscalar(time)
    gh = loadigrfcoefs(time);
    nmax = sqrt(numel(gh) + 1) - 1;
% Assume a vector input means the coefficients are the input.
else
    gh = time;
    nmax = sqrt(numel(gh) + 1) - 1;
    % nmax should be an integer.
    if nmax - round(nmax) ~= 0
        error('igrf:timeInputInvalid', ['TIME input should either be ' ...
            'a single date or a valid coefficient vector.']);
    end
end

% We need cos(m*phi) and sin(m*phi) multiple times, so precalculate into a
% vector here:
cosphi = cos((1:nmax)*phi);
sinphi = sin((1:nmax)*phi);

Pmax = (nmax+1)*(nmax+2)/2;

%%% BEGIN MAGNETIC FIELD CALCULATION %%%
% Initialize variables used in for loop below.
Br = 0; Bt = 0; Bp = 0;
 P = zeros(1, Pmax);  P(1) = 1;  P(3) = sintheta;
dP = zeros(1, Pmax); dP(1) = 0; dP(3) = costheta;

% For this initial condition, the first if will result in n = 1, m = 0.
m = 1; n = 0; coefindex = 1;

a_r = (Rearth_km/r)^2;

% Increment through all the n's and m's. gh will be a vector with g
% followed by h for incrementing through n and m except when h would be
% redundant (i.e., when m = 0).
for Pindex = 2:Pmax

    % Increment to the next n when m becomes larger than n.
    if n < m
        m = 0;
        n = n + 1;
        a_r = a_r*(Rearth_km/r); % We need (Rearth_km./r)^(n+2)
    end

    % Calculate P and dP. They are given recursively according to:
    %
    % P(0, 0) = 1, P(1, 1) = sin(theta) <- Specified above
    % P(n, n) = sqrt(1 - 1/(2n))*sin(theta)*P(n-1, n-1)
    % P(n, m) = (2n - 1)/sqrt(n^2 - m^2)*cos(theta)*P(n-1, m) -
    %     sqrt(((n-1)^2 - m^2) / (n^2 - m^2)) * P(n-2, m)
    %
    % dP(0, 0) = 0, dP(1, 1) = cos(theta) <- Specified above
    % dP(n, n) = sqrt(1 - 1/(2n))*(sin(theta)*dP(n-1, n-1) +
    %     cos(theta)*P(n-1, n-1))
    % dP(n, m) = (2n - 1)/sqrt(n^2 - m^2)*(cos(theta)*dP(n-1, m) -
    %     sin(theta)*P(n-1, m)) - sqrt(((n-1)^2 - m^2)/(n^2 - m^2))*
    %     dP(n-2, m)
    if m < n && Pindex ~= 3 % (Pindex=3 is n=1, m=1, initial cond. above)
        last1n = Pindex - n;
        last2n = Pindex - 2*n + 1;
        P(Pindex) = (2*n - 1)/sqrt(n^2 - m^2)*costheta*P(last1n) - ...
            sqrt(((n-1)^2 - m^2) / (n^2 - m^2)) * P(last2n);
        dP(Pindex) = (2*n - 1)/sqrt(n^2 - m^2)*(costheta*dP(last1n) - ...
            sintheta*P(last1n)) - sqrt(((n-1)^2 - m^2) / (n^2 - m^2)) * ...
            dP(last2n);
    elseif Pindex ~= 3
        lastn = Pindex - n - 1;
        P(Pindex) = sqrt(1 - 1/(2*m))*sintheta*P(lastn);
        dP(Pindex) = sqrt(1 - 1/(2*m))*(sintheta*dP(lastn) + ...
            costheta*P(lastn));
    end

    % Calculate the magnetic field components as a running sum. Find
    % explicit expressions for these in Global Earth Physics: a Handbook of
    % Physical Constants by Thomas J. Aherns (1995), pg. 49. Link:
    % http://books.google.com/books?id=aqjU_NHyre4C&lpg=PP1&dq=Global%20
    % earth%20physics%3A%20a%20handbook%20of%20physical%20constants&pg=PA49
    % #v=onepage&q&f=false
    % (except equation 6 is missing a required 1/sin(theta) and m; correct
    % equations on page 5 (equations 3a-3c) of:
    % http://hanspeterschaub.info/Papers/UnderGradStudents/
    % MagneticField.pdf)
    if m == 0 % Implies h = 0, so only coefficient in gh is g
        coef = a_r*gh(coefindex); %*cos(0*phi) = 1
        Br = Br + (n+1)*coef*P(Pindex);
        Bt = Bt - coef*dP(Pindex);
        % Bp is 0 for m = 0.
        coefindex = coefindex + 1; % Only need to skip over g this time.
    else
        coef = a_r*(gh(coefindex)*cosphi(m) + gh(coefindex+1)*sinphi(m));
        Br = Br + (n+1)*coef*P(Pindex);
        Bt = Bt - coef*dP(Pindex);
        if sintheta == 0 % Use different formula when dividing by 0.
            Bp = Bp - costheta*a_r*(-gh(coefindex)*sinphi(m) + ...
                gh(coefindex+1)*cosphi(m))*dP(Pindex);
        else
            Bp = Bp - 1/sintheta*a_r*m*(-gh(coefindex)*sinphi(m) + ...
                gh(coefindex+1)*cosphi(m))*P(Pindex);
        end
        coefindex = coefindex + 2; % Skip over g and h this time.
    end

    % Increment m.
    m = m + 1;

end

% Convert from spherical to (x,y,z) = (North,East,Down).
Bx = -Bt;
By = Bp;
Bz = -Br;

% Convert back to geodetic coordinates if necessary.
Bx_old = Bx;
Bx = Bx.*cd + Bz.*sd;
Bz = Bz.*cd - Bx_old.*sd;


text:
classdef geodata
    %   GEODATA: Geographical data class
    %   Handles geographical data describing coastlines or other features in
    %   the form of a shapefile and topobathy in the form of a DEM
    %   Copyright (C) 2018  Keith Roberts & William Pringle
    %
    %   The following inputs (with default values) are available for the
    %   geodata method:
    %     defval = 0; % placeholder value if arg is not passed.
    %     addOptional(p,'bbox',defval);
    %     addOptional(p,'shp',defval);
    %     addOptional(p,'h0',defval);
    %     addOptional(p,'dem',defval);
    %     addOptional(p,'backupdem',defval);
    %     addOptional(p,'fp',defval);
    %     addOptional(p,'weirs',defval);
    %     addOptional(p,'pslg',defval);
    %     addOptional(p,'boubox',defval);
    %     addOptional(p,'window',defval);
    %
    %   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/>.
    properties
        bbox  % bounding coordinates
        boubox % bbox coordinates as a CCW polygon
        mainland % mainland boundary.
        outer % outer boundary.
        inner % island boundary(ies)
        mainlandb % height of mainland
        innerb    % height of inner
        mainlandb_type % type of mainland (lake or river or ocean)
        innerb_type    % type of inner (lake or river or ocean)
        weirs % weir crestlines
        weirPfix % boundaries of weir
        weirEgfix % edges of weir
        ibconn_pts % cell-array of paired weir nodes
        inpoly_flip % reverse the notion of "in"
        contourfile %  cell-array of shapefile
        demfile % filename of dem
        BACKUPdemfile % filename of dem
        h0 % min. edgelength (meters)
        window  % smoothing window on boundary (default 5 points)
        fp % deprecated but kept for backwards capability
        Fb % linear gridded interpolant of DEM
        Fb2 % linear gridded interpolant of backup DEM
        x0y0 % bottom left of structure grid or position (0,0)
        pslg % piecewise liner straight line graph
        spacing = 2.0 ; %Relative spacing along polygon, large effect on computational efficiency of signed distance.
        gridspace
        shapefile_3d % if the shapefile has a height attribute
    end

    methods

        function obj = geodata(varargin)
            % Class constructor to parse NetCDF DEM data, NaN-delimited vector,
            % or shapefile that defines polygonal boundary of meshing
            % domain.
            % options
            %addOptional(p,'bbox',defval);
            %addOptional(p,'shp',defval);
            %addOptional(p,'h0',defval);
            %addOptional(p,'dem',defval);
            %addOptional(p,'backupdem',defval);
            %addOptional(p,'fp',defval);
            %addOptional(p,'weirs',defval);
            %addOptional(p,'pslg',defval);
            %addOptional(p,'boubox',defval);

            % Check for m_map dir
            M_MAP_EXISTS=0 ;
            if exist('m_proj','file')==2
              M_MAP_EXISTS=1 ;
            end
            if M_MAP_EXISTS~=1
              error('Where''s m_map? Please read the user guide')
            end

            % Check for utilties dir
            UTIL_DIR_EXISTS=0  ;
            if exist('inpoly.m','file')
              UTIL_DIR_EXISTS=1;
            end
            if UTIL_DIR_EXISTS~=1
              error('Where''s the utilities directory? Please read the user guide')
            end

            % Check for dataset dir
            DATASET_DIR_EXISTS=0 ;
            if exist('datasets','dir')==7
                DATASET_DIR_EXISTS=1 ;
            end
            if DATASET_DIR_EXISTS~=1
                warning('We suggest you to place your files in a directory called "datasets". Please read the user guide')
            end

            p = inputParser;

            defval = 0; % placeholder value if arg is not passed.
            % add name/value pairs
            addOptional(p,'bbox',defval);
            addOptional(p,'shp',defval);
            addOptional(p,'h0',defval);
            addOptional(p,'dem',defval);
            addOptional(p,'backupdem',defval);
            addOptional(p,'fp',defval);
            addOptional(p,'weirs',defval);
            addOptional(p,'pslg',defval);
            addOptional(p,'boubox',defval);
            addOptional(p,'window',defval);
            addOptional(p,'shapefile_3d',defval);

            % parse the inputs
            parse(p,varargin{:});
            % store the inputs as a struct
            inp=p.Results;
            % get the fieldnames of the edge functions
            fields = fieldnames(inp);
            % loop through and determine which args were passed.
            % also, assign reasonable default values if some options were
            % not assigned.
            obj.boubox = [] ;
            for i = 1 : numel(fields)
                type = fields{i};
                switch type
                    % parse aux options first
                    case('bbox')
                        obj.bbox = inp.(fields{i});
                        if any(obj.bbox ~=0)
                            obj.bbox = inp.(fields{i});
                            %elseif ~ischar(inp.(fields{contains(fields,'dem')}))
                            %    error('No bbox specified!');
                        end
                    case('h0')
                        obj.h0= inp.(fields{i});
                        if obj.h0 ~=0
                            obj.h0 = inp.(fields{i});
                        else
                            error('No h0 specified!');
                        end
                    case('shp')
                        obj.contourfile = inp.(fields{i});
                        if ~iscell(obj.contourfile) && ~ischar(obj.contourfile)
                            obj.contourfile = [];
                        end
                    case('dem')
                        obj.demfile= inp.(fields{i});
                        if obj.demfile ~=0
                            obj.demfile = inp.(fields{i});
                        else
                            obj.demfile = [];
                        end
                    case('backupdem')
                        obj.BACKUPdemfile= inp.(fields{i});
                        if obj.BACKUPdemfile ~=0
                            obj.BACKUPdemfile = inp.(fields{i});
                        else
                            obj.BACKUPdemfile = [];
                        end
                    case('fp')
                        obj.fp= inp.(fields{i});
                        if obj.fp ~=0
                            obj.fp = inp.(fields{i});
                        end
                    case('outer')
                        obj.outer = inp.(fields{i});
                        if obj.outer(1) ~=0
                            obj.outer = inp.(fields{i});
                        end
                    case('pslg')
                        obj.pslg = inp.(fields{i});
                        if obj.pslg(1) ~=0
                            obj.pslg = inp.(fields{i});
                        else
                            obj.pslg = [];
                        end
                    case('mainland')
                        obj.mainland = inp.(fields{i});
                        if obj.mainland(1) ~=0
                            obj.mainland = inp.(fields{i});
                        end
                    case('inner')
                        obj.inner = inp.(fields{i});
                        if obj.inner(1) ~=0
                            obj.inner = inp.(fields{i});
                        end
                    case('boubox')
                        obj.boubox = inp.(fields{i}) ;
                    case('window')
                        obj.window = inp.(fields{i}) ;
                        if obj.window == 0
                            % Default value
                            obj.window = 5;
                        end
                    case('shapefile_3d')
                         obj.shapefile_3d = inp.(fields{i}) ;
                    case('weirs')
                        if ~iscell(inp.(fields{i})) && ~isstruct(inp.(fields{i})) && inp.(fields{i})==0, continue; end
                        if ~iscell(inp.(fields{i})) && ~isstruct(inp.(fields{i}))
                            error('Data for weirs must be in a cell-array or struct. Please see the user guide.');
                        end
                        obj.weirs = inp.(fields{i}) ;
                        noWeirs   = length(obj.weirs) ;
                        disp(['INFO: User has passed ',num2str(noWeirs),' weir crestlines.']) ;
                        obj.weirPfix = [] ; obj.weirEgfix = [] ;
                        if iscell(obj.weirs)
                            ss = obj.weirs{1}(1);
                        else
                            ss = obj.weirs(1).X(1);
                        end
                        if ss ~= 0
                            for ii = 1 : noWeirs
                                if iscell(obj.weirs)
                                    crestlines = obj.weirs{ii}(:,1:2) ;
                                    width      = obj.weirs{ii}(1,3)/111e3 ;
                                    weir_min_ele = obj.weirs{ii}(1,4)/111e3;
                                else
                                    crestlines = [obj.weirs(ii).X obj.weirs(ii).Y];
                                    width      = obj.weirs(ii).width/111e3 ;
                                    weir_min_ele = obj.weirs(ii).min_ele/111e3;
                                end
                                if width == 0
                                  error('Please specify non-zero width of weir in meters!');
                                end
                                % user-defined spacing along face of weir
                                if weir_min_ele ~= 0
                                    [tempPfix,tmpEgfix,obj.ibconn_pts{ii}] = GenerateWeirGeometry(crestlines,width,...
                                          weir_min_ele,0) ;
                                else
                                % by default spacing along face of weir is
                                % MIN_EL
                                    [tempPfix,tmpEgfix,obj.ibconn_pts{ii}] = GenerateWeirGeometry(crestlines,width,...
                                        2*(obj.h0/111e3),0) ;
                                end
                                [tmpEgfix]=renumberEdges(tmpEgfix) ;
                                weirLength(ii) = length(tempPfix) ;
                                if ii~=1
                                    tmpEgfix = tmpEgfix + length(obj.weirPfix) ;
                                end
                                obj.weirPfix = [obj.weirPfix; tempPfix] ;
                                obj.weirEgfix = [obj.weirEgfix; tmpEgfix] ;
                            end
                        end
                end
            end
            % Basic error handling should go here.
            % if not bbox and no dem, you're outta luck
            if size(obj.bbox,1) == 1 && isempty(obj.demfile) && isempty(obj.pslg)
                error('No DEM supplied and no bbox supplied, sorry pal');
            end
            obj.gridspace    = abs(obj.h0)/111e3; %point spacing along polygon,
            % Get bbox information from demfile if not supplied
            if size(obj.bbox,1) == 1 && isempty(obj.pslg)
                obj = ParseDEM(obj,'bbox');
            end

            if obj.bbox == 0
                error('No bbox supplied. If you are using the pslg option then you need to supply a bbox.')
            elseif size(obj.bbox,1) == 2
                % Typical square bbox type
                % Make the bounding box 5 x 2 matrix in clockwise order
                obj.boubox = bbox_to_bou(obj.bbox);
            else
                % Handle non-square bbox type regions
                obj.boubox = obj.bbox;
                if ~isnan(obj.boubox(end,1))
                    obj.boubox(end+1,:) = [NaN NaN];
                end
                obj.bbox = [min(obj.boubox(:,1)) max(obj.boubox(:,1))
                    min(obj.boubox(:,2)) max(obj.boubox(:,2))] ;
            end

            obj = ParseShoreline(obj) ;

            % kjr Add the weir faux islands to the inner geometry
            if ~isempty(obj.weirPfix)
                idx = [0; cumsum(weirLength)']+1 ;
                tmp = [] ;
                for ii = 1 : noWeirs
                    tmp =  [tmp;
                        [obj.weirPfix(idx(ii):idx(ii+1)-1,:)
                        obj.weirPfix(idx(ii),:)]
                        NaN NaN] ;
                end
                obj.inner = [obj.inner ; NaN NaN ;  tmp ] ;
            end

            % Ensure inpoly flip is correct
            obj = check_connectedness_inpoly(obj);

            disp(['Read in meshing boundary: ',obj.contourfile]);

            if obj.BACKUPdemfile~=0
                obj = ParseDEM(obj,'BackUp') ;
            end

            obj = ParseDEM(obj) ;

        end

        function obj = ParseShoreline(obj)
            % Read in the geometric meshing boundary information from a
            % ESRI-shapefile
            if ~isempty(obj.contourfile)
                if ~iscell(obj.contourfile)
                    obj.contourfile = {obj.contourfile};
                end

                polygon_struct = Read_shapefile( obj.contourfile, [], ...
                    obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d);

                % Unpack data from function Read_Shapefile()s
                obj.outer     = polygon_struct.outer;
                obj.mainland  = polygon_struct.mainland;
                obj.inner     = polygon_struct.inner;
                obj.mainlandb = polygon_struct.mainlandb;
                obj.innerb    = polygon_struct.innerb;
                obj.mainlandb_type = polygon_struct.mainlandb_type;
                obj.innerb_type    = polygon_struct.innerb_type;

                % Read in the geometric meshing boundary information from a
                % NaN-delimited vector.
            elseif ~isempty(obj.pslg)

                % Handle the case for user defined mesh boundary information
                polygon_struct = Read_shapefile( [], obj.pslg, ...
                    obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d);

                % Unpack data from function Read_Shapefile()s
                obj.outer     = polygon_struct.outer;
                obj.mainland  = polygon_struct.mainland;
                obj.inner     = polygon_struct.inner;
                obj.mainlandb = polygon_struct.mainlandb;
                obj.innerb    = polygon_struct.innerb;
                obj.mainlandb_type = polygon_struct.mainlandb_type;
                obj.innerb_type    = polygon_struct.innerb_type;

            else
                % set outer to the boubox
                obj.outer = obj.boubox;
            end
            obj = ClassifyShoreline(obj) ;
        end

        function obj = ClassifyShoreline(obj)
            % Helper function to...
            % 1) Read the data from supplied file
            % 2) Classify segments as either inner, outer, or mainland.
            % 3) Ensure point spacing is adequate to support mesh sizes
            %    and ensure computational efficiency when calculating the
            %    signed distance function.

            % Check if no mainland segments, set outer
            % to the boubox.
            if isempty(obj.mainland)
                obj.outer = [ ];
                obj.outer = obj.boubox;
            end

            % Make sure the shoreline components have spacing of
            % gridspace/spacing
            [la,lo] = my_interpm(obj.outer(:,2),obj.outer(:,1),...
                obj.gridspace/obj.spacing);
            obj.outer = [];  obj.outer(:,1) = lo; obj.outer(:,2) = la;
            outerbox = obj.outer(1:find(isnan(obj.outer(:,1)),1,'first'),:);

            if ~isempty(obj.mainland)
                [la,lo] = my_interpm(obj.mainland(:,2),obj.mainland(:,1),...
                    obj.gridspace/obj.spacing);
                obj.mainland = []; obj.mainland(:,1) = lo; obj.mainland(:,2) = la;
            end

            if ~isempty(obj.inner)
                [la,lo]=my_interpm(obj.inner(:,2),obj.inner(:,1),...
                    obj.gridspace/obj.spacing);
                obj.inner = []; obj.inner(:,1) = lo; obj.inner(:,2) = la;
            end
            clearvars lo la

            % Smooth the coastline (apply moving average filter).
            if obj.window > 1
                disp(['Smoothing coastline with ' ...
                    num2str(obj.window) ' point window'])
                if ~isempty(obj.outer)
                    obj.outer = smooth_coastline(obj.outer,obj.window,0);
                end
                if ~isempty(obj.mainland)
                    obj.mainland = smooth_coastline(obj.mainland,obj.window,0);
                end
                if ~isempty(obj.inner)
                    obj.inner = smooth_coastline(obj.inner,obj.window,0);
                end
            else
                disp('No smoothing of coastline enabled')
            end

            % KJR: Coarsen portions of outer, mainland
            % and inner outside bbox.
            iboubox = bbox_to_bou(obj.bbox);
            iboubox(:,1) = 1.10*iboubox(:,1)+(1-1.10)*mean(iboubox(1:end-1,1));
            iboubox(:,2) = 1.10*iboubox(:,2)+(1-1.10)*mean(iboubox(1:end-1,2));

            % Coarsen outer
            obj.outer = coarsen_polygon(obj.outer,iboubox);

            % Coarsen inner and move parts that overlap with bounding box
            % to mainland
            if ~isempty(obj.inner)
                obj.inner = coarsen_polygon(obj.inner,iboubox);
            end

            % Coarsen mainland and remove parts that overlap with bounding
            % box (note this doesn't change the polygon used for inpoly,
            % only changes the distance function used for edgefx)
            if ~isempty(obj.mainland)
                obj.mainland = coarsen_polygon(obj.mainland,iboubox);
            end

        end

        function obj = ParseDEM(obj,varargin)
            % obj = ParseDEM(obj,varargin)
            % set 'BackUp' for reading backupdem
            % set 'bbox' for reading only bbox
            fname = obj.demfile ; backup = 0;
            if any(strcmp(varargin,'BackUp'))
                backup = 1;
                fname = obj.BACKUPdemfile ;
            end

            % Process the DEM for the meshing region.
            if ~isempty(fname)

                % Find name of x, y (one that has 1 dimensions)
                % z values (use one that has 2 dimensions)
                [xvn, yvn, zvn] = getdemvarnames(fname);

                % Read x and y
                x = double(ncread(fname,xvn));
                y = double(ncread(fname,yvn));

                if any(strcmp(varargin,'bbox'))
                    obj.bbox = [min(x) max(x); min(y) max(y)];
                    return;
                end

                modbox = [0 0];
                if obj.bbox(1,2) > 180 && obj.bbox(1,1) < 180 && min(x) < 0
                    % bbox straddles 180/-180 line
                    loop = 2;
                else
                    loop = 1;
                    if max(x) > 180
                        if obj.bbox(1,1) < 0; modbox(1) = 1; end
                        if obj.bbox(1,2) < 0; modbox(2) = 1; end
                    elseif min(x) < 0
                        if obj.bbox(1,1) > 180; modbox(1) = -1; end
                        if obj.bbox(1,2) > 180; modbox(2) = -1; end
                    end
                end
                J = find(y >= obj.bbox(2,1) & y <= obj.bbox(2,2));
                I = []; demz = [];
                for nn = 1:loop
                    bboxt = obj.bbox;
                    bboxt(1,:) =  bboxt(1,:) + modbox.*360;
                    if loop == 2
                        if nn == 1
                            bboxt(1,2) = 180;
                        else
                            bboxt(1,1) = -180;
                            bboxt(1,2) = bboxt(1,2) - 360;
                        end
                    end
                    It = find(x >= bboxt(1,1) & x <= bboxt(1,2));
                    I = [I; It];

                    % At this point, detect the memory footprint of the
                    % DEM subset and compute the required stride necessary
                    % to satisfy the memory requirements
                    if nn == 1
                        AVAILABLE_MEMORY=4; % 4 gb;
                        mult = (obj.bbox(1,2) - obj.bbox(1,1))/...
                                   (bboxt(1,2) - bboxt(1,1));
                        peak_mem = mult*length(I)*length(J)*4/1e9; % in GB assuming single
                        STRIDE_MEM = ceil(sqrt(peak_mem/AVAILABLE_MEMORY));
                        DEM_GRIDSPACE = (x(2)-x(1))*111e3; % in meters
                        STRIDE_H0 = ceil(obj.h0/DEM_GRIDSPACE); % skip # of DEM entires
                        STRIDE = max(STRIDE_H0, STRIDE_MEM);
                        if STRIDE > 1
                            if STRIDE_MEM > STRIDE_H0
                                warning(['DEM would occupy ',num2str(peak_mem),...
                                         'GB of RAM. DEM will be downsampled ' ...
                                         'by a stride of ' num2str(STRIDE)])
                            else
                                warning(['DEM will be downsampled by a stride ' ...
                                         'of ' num2str(STRIDE) ' to match h0'])
                            end
                        end

                    end
                    % grab only the portion that was requested
                    if STRIDE > 1 && STRIDE_MEM > STRIDE_H0
                        % with a stride to save memory
                        LX = length(It(1:STRIDE:end));
                        LY = length(J(1:STRIDE:end));
                        demzt = single(ncread(fname,zvn,[It(1) J(1)],...
                                       [LX LY],[STRIDE,STRIDE]));
                    else
                        % faster to read in with no stride if not memory
                        % bound
                        demzt = single(ncread(fname,zvn,[It(1) J(1)],...
                                       [length(It) length(J)]));
                        if STRIDE > 1
                            demzt = demzt(1:STRIDE:end,1:STRIDE:end);
                        end
                    end
                    if isempty(demz)
                        demz = demzt;
                    else
                        demz = cat(1,demz,demzt);
                    end
                end
                x = x(I(1:STRIDE:end)); y = y(J(1:STRIDE:end));
                if obj.bbox(1,2) > 180
                    x(x < 0) = x(x < 0) + 360;
                    [x1,IA] = unique(x);
                    if length(x) > length(x1)
                        x = x1; demz = demz(IA,:);
                    end
                end
                % handle DEMS from packed starting from the bottom left
                if y(2) < y(1)
                    y = flipud(y) ;
                    demz = fliplr(demz) ;
                end
                % Determine bottom left corner of DEM
                % (after possible flipping of DEM packing)
                obj.x0y0 = [x(1),y(1)];

                % check for any invalid values
                bad = isnan(demz);
                if sum(bad(:)) > 0 && ~backup
                    warning('ALERT: Invalid and/or missing DEM values detected..check DEM');
                    if ~isempty(obj.BACKUPdemfile)
                        disp('Replacing invalid values with back-up DEMfile');
                        [demx,demy] = ndgrid(x,y) ;
                        demz(bad) = obj.Fb2(demx(bad),demy(bad)) ;
                    else
                        % just replace it with NaNs
                        demz(bad) = NaN ;
                    end
                end

                % creating back up interpolant
                if backup
                    obj.Fb2   = griddedInterpolant({x,y},demz,...
                        'linear','nearest');
                    % clear data from memory
                    clear x y demz
                else
                    % main interpolant
                    if obj.BACKUPdemfile~=0
                        obj.Fb   = griddedInterpolant({x,y},demz,...
                            'linear','none'); % no extrapolation (so use Fb2)
                    else
                        obj.Fb   = griddedInterpolant({x,y},demz,...
                            'linear','nearest');
                    end
                    % clear data from memory
                    clear x y demz
                end

                disp(['Read in demfile ',fname]);
            end

            % Handle the case of no dem
            if isempty(obj.x0y0)
                obj.x0y0 = [obj.bbox(1,1), obj.bbox(2,1)];
            end

            function [xvn, yvn, zvn] = getdemvarnames(fname)
                % Define well-known variables for longitude and latitude
                % coordinates in Digital Elevation Model NetCDF file (CF
                % compliant).
                xvn = []; yvn = []; zvn = [];
                wkv_x = {'x','Longitude','longitude','lon'} ;
                wkv_y = {'y','Latitude', 'latitude','lat'} ;
                finfo = ncinfo(fname);
                for ii = 1:length(finfo.Variables)
                    if ~isempty(xvn) && ~isempty(yvn) && ~isempty(zvn); break; end
                    if length(finfo.Variables(ii).Size) == 1
                        if isempty(xvn) && ...
                                any(strcmp(finfo.Variables(ii).Name,wkv_x))
                            xvn = finfo.Variables(ii).Name;
                        end
                        if isempty(yvn) && ...
                                any(strcmp(finfo.Variables(ii).Name,wkv_y))
                            yvn = finfo.Variables(ii).Name;
                        end
                    elseif length(finfo.Variables(ii).Size) == 2
                        if isempty(zvn)
                            zvn = finfo.Variables(ii).Name;
                        end
                    end
                end
                if isempty(xvn)
                    error('Could not locate x coordinate in DEM') ;
                end
                if isempty(yvn)
                    error('Could not locate y coordinate in DEM') ;
                end
                if isempty(zvn)
                    error('Could not locate z coordinate in DEM') ;
                end
            end

        end

        function obj = check_connectedness_inpoly(obj)
            % Check for connected polygons and if not connected,
            %  whether to flip the inpoly result to ensure meshing
            %  on the coastal side of the polygon.

            obj.inpoly_flip = 0;
            % return if outer polygon is connected
            shpEnd = find(isnan(obj.outer(:,1)));
            [~,loc] = max(diff(shpEnd));
            shpEnd = vertcat(0,shpEnd); loc = loc+1;
            if abs(sum(obj.outer(shpEnd(loc)+1,:)) - ...
                    sum(obj.outer(shpEnd(loc+1)-1,:))) > eps
                disp('Warning: Shapefile is unconnected... continuing anyway')
            end

            % check for inpoly goodness read the GSHHS checker
            ps = Read_shapefile( {'GSHHS_l_L1'}, [], ...
                obj.bbox, obj.gridspace, obj.boubox, 0, 0 );

            % make a "fake" tester grid
            x = linspace(obj.bbox(1,1),obj.bbox(1,2),100);
            y = linspace(obj.bbox(2,1),obj.bbox(2,2),100);
            edges = Get_poly_edges( [ps.outer; ps.inner] );
            in_Test = inpoly([x',y'],[ps.outer; ps.inner],edges);
            polytester = [obj.outer; obj.inner];
            edges = Get_poly_edges( polytester );
            in_Shpf = inpoly([x',y'],polytester,edges);
            % if more than half of thepoints disagree between Test and Shpf
            % lets flip the inpoly
            if length(find(xor(in_Shpf,in_Test))) > 50
                obj.inpoly_flip = 1;
                disp(['Shapefile inpoly is inconsistent ' ...
                    'with GHSSS test file, flipping the inpoly test'])
            end

            % if flooplaind meshing, flip the inpoly test
            if obj.fp
                obj.inpoly_flip = mod(1,obj.inpoly_flip);
            end
        end

        function obj = close(obj)
            % Clips the mainland segment with the boubox.
            % Performs a breadth-first search given a seed position
            % of the piecewise-straight line graph (PSLG) that is used to define the meshing boundary.
            % This returns back an updated geodata class instance with the outer boundary clipped with the boubox.
            % kjr,und,chl 2018

            if isempty(obj.Fb)
                warning('ALERT: Meshing boundary is not a polygon!')
                warning('ALERT: DEM is required to clip line segment with polygon')
                seed=input('Enter coordinate of seed location to clip: ');
            else
                % Guess seed location: submerged portion of domain within
                % bbox?
                [demx,demy] = ndgrid(obj.x0y0(1):obj.h0/111e3:obj.bbox(1,2), ...
                    obj.x0y0(2):obj.h0/111e3:obj.bbox(2,2));
                demz = obj.Fb(demx,demy);
                [idx]=find(demz(:) < -10) ; % find points deeper than -10 m below sea level
                allpts=[demx(:),demy(:)] ;
                cands = allpts(idx,:);
                if(isnan(obj.boubox(end,1)))
                    edges=Get_poly_edges(obj.boubox) ;
                    [in]=inpoly(cands,obj.boubox,edges) ; %  determine which of these points is in the domain
                else
                    [in]=inpoly(cands,obj.boubox) ;       %  determine which of these points is in the domain
                end
                cands2=cands(in,:) ;
                seed = cands2(50,:) ;
            end

            geom = [obj.mainland; obj.boubox] ;

            [NODE,PSLG]=getnan2(geom) ;

            [NODE2,PSLG2,PART2] = bfsgeo2(NODE,PSLG,seed) ;

            POLY = extdom_polygon(PSLG2(PART2{1},:),NODE2,-1) ;

            new_outer = cell2mat(POLY') ;

            [la,lo] = my_interpm(new_outer(:,2),new_outer(:,1),((obj.h0/2)/111e3))  ;

            new_outer = [lo la] ;

            obj.outer = new_outer ;

            % reset this to default
            obj.inpoly_flip = 0 ;
        end

        function obj = extractContour(obj,ilev)
            % Extract a geometric contour from the DEM at elevation ilev.
            % obj = extractContour(obj,ilev)
            %
            % Can use to get the mean sea level contour, e.g.;
            % gdat = geodata('pslg',0,'h0',min_el,'dem',dem); % make the dummy gdat for the dem extents;
            % lmsl = extractContour(gdat,0); %using the dummy gdat with dem info to get the 'lmsl' gdat with the 0-m contour.

            [node,edge] = ...
                getiso2( obj.Fb.GridVectors{1},obj.Fb.GridVectors{2},...
                double(obj.Fb.Values'),ilev) ;

            polyline = cell2mat(extdom_polygon(edge,node,-1,1,10)') ;

            obj = geodata('pslg',polyline,'bbox',obj.bbox,...
                          'h0',obj.h0,'dem',obj.demfile) ;

        end

        function plot(obj,type,projection,holdon)
            % plot(obj,type,projection,holdon)
            % Plot geodata class info
            %
            % Inputs:
            % obj  : geodata class object [required input]
            %
            % optional inputs =>....
            % type : i) 'shp' [default] - plots the shapelines (shoreline) only
            %       ii) 'dem' - plots the dem bathy in addition to shapelines
            %      iii) 'omega' - hatches the meshing domain in addition to plotting shapelines
            % projection : choose the projection type from m_map options
            %              (default is Mercator)
            % holdon : plot on a new (= 0 [default]) or existing (= 1) figure?

            if nargin == 1 || isempty(type)
                type = 'shp';
            end
            if nargin < 3 || isempty(projection)
                projection = 'Mercator';
            end
            if nargin < 4 || isempty(holdon)
                holdon = 0;
            end

            % plotting on new or existing figure?
            if ~holdon
                % setup the projection
                bufx = 0.2*(obj.bbox(1,2) - obj.bbox(1,1));
                bufy = 0.2*(obj.bbox(2,2) - obj.bbox(2,1));
                if ~isempty(regexp(projection,'ste'))
                    m_proj(projection,'lat',min(obj.bbox(2,:)),...
                        'long',mean(obj.bbox(1,:)),'radius',...
                        min(179.9,1.20*max(diff(obj.bbox(2,:)))));
                else
                    lmin = -180; lmax = +180;
                    if obj.bbox(1,2) > 180; lmax = 360; lmin = 0; end
                    lon1 = max(lmin,obj.bbox(1,1) - bufx);
                    lon2 = min(lmax,obj.bbox(1,2) + bufx);
                    lat1 = max(- 90,obj.bbox(2,1) - bufy);
                    lat2 = min(+ 90,obj.bbox(2,2) + bufy);
                    m_proj(projection,...
                           'long',[lon1, lon2],'lat',[lat1, lat2]);
                end
                % plot on new figure
                figure;
                colori = 'g-'; % set island color
                colorm = 'r-'; % set mainland color
            else
                colori = 'b-'; % set island color
                colorm = 'm-'; % set mainland color
            end
            hold on
            % select optional types
            switch type
                case('dem')
                    % interpolate DEM's bathy linearly onto our
                    % edgefunction grid (or a coarsened version of it for
                    % memory considerations)
                    mem = inf; stride = obj.h0/111e3;
                    while mem > 1
                        xx = obj.x0y0(1):stride:obj.bbox(1,2);
                        yy = obj.x0y0(2):stride:obj.bbox(2,2);
                        xs = whos('xx'); ys = whos('yy');
                        mem = xs.bytes*ys.bytes/1e9;
                        stride = stride*2;
                    end
                    [demx,demy] = ndgrid(xx,yy);
                    demz = obj.Fb(demx,demy);
                    m_fastscatter(demx(:),demy(:),demz(:));
                    cb = colorbar; ylabel(cb,'topo-bathy depth [m]')
                case('omega')
                    % hatch the meshing domain, Omega
                    [demx,demy] = ndgrid(obj.x0y0(1):obj.h0/111e3:obj.bbox(1,2), ...
                        obj.x0y0(2):obj.h0/111e3:obj.bbox(2,2));
                    edges = Get_poly_edges( [obj.outer; obj.inner] );
                    in = inpoly([demx(:),demy(:)],[obj.outer; obj.inner], edges);
                    long = demx(~in); lati = demy(~in);
                    m_hatch(obj.boubox(1:end-1,1),...
                        obj.boubox(1:end-1,2),'cross',45,0.05);
                    m_plot(long,lati,'.','Color','white')
            end
            if ~isempty(obj.mainland) && obj.mainland(1) ~= 0
                h1 = m_plot(obj.mainland(:,1),obj.mainland(:,2),...
                    colorm,'linewi',1); hold on;
            end
            if ~isempty(obj.inner) && obj.inner(1) ~= 0
                h2 = m_plot(obj.inner(:,1),obj.inner(:,2),...
                    colori,'linewi',1); hold on;
            end
            if ~isempty(obj.weirs)
                if isstruct(obj.weirs) ~= 0
                    for ii =1 : length(obj.weirs)
                        h3 = m_plot(obj.weirs.X,obj.weirs.Y,...
                            'm-','linewi',1); hold on;
                    end
                else
                    for ii =1 : length(obj.weirs)
                        h3 = m_plot(obj.weirs{ii}(:,1),obj.weirs{ii}(:,2),...
                            'm-','linewi',1); hold on;
                    end
                end
            end
            [la,lo] = my_interpm(obj.boubox(:,2),obj.boubox(:,1),...
                0.5*obj.h0/111e3);
            m_plot(lo,la,'k--','linewi',2);
            if ~holdon
                m_grid('xtick',10,'tickdir','out','yaxislocation','left','fontsize',10);
            end
            if exist('h1','var') && exist('h2','var') && exist('h3','var')
                legend([h1 h2,h3],{'mainland' 'inner' 'weirs'},'Location','NorthWest')
            elseif exist('h1','var') && exist('h2','var')
                legend([h1 h2],{'mainland' 'inner'},'Location','NorthWest')
            end
        end

    end

end

author: CHLNDDEV


text:
function ou_euler_maruyama ( theta, mu, sigma, x0, tmax, n, r, seed )

%*****************************************************************************80
%
%% OU_EULER_MARUYAMA applies Euler-Maruyama to the Ornstein-Uhlenbeck SDE.
%
%  Discussion:
%
%    The stochastic differential equation (SDE) is:
%
%      dx = theta * ( mu - x(t) ) dt + sigma dW,
%      x(0) = x0,
%
%    The discretized Brownian path uses a constant stepsize.
%
%    A "large" time step DT_LARGE is used for the smooth portion of the
%    equation, while a smaller time step DT_SMALL is used for the
%    discretized Brownian path.  We take R small steps to equal one
%    large step, so that:
%
%      dt_large = r * dt_small = tmax / n
%
%    For an SDE of the form:
%
%      dx = f(x(t)) dt + g(x(t)) dW(t)
%
%    the Euler-Maruyama method has the form:
%
%      x(j) = x(j-1) + f(X(j-1)) * dt_large + g(X(j-1)) * dW(j-1)
%
%    where dW(j-1) is approximated by the sum of R normal random values
%    multiplied by the square root of DT_SMALL.
%
%    Note that if SIGMA is zero, the problem becomes deterministic,
%    with solution
%
%      x(t) = mu + ( x0 - mu ) * exp ( - theta * t )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Desmond Higham,
%    An Algorithmic Introduction to Numerical Simulation of
%    Stochastic Differential Equations,
%    SIAM Review,
%    Volume 43, Number 3, September 2001, pages 525-546
%
%  Parameters:
%
%    Input, real THETA, MU, SIGMA, the value of problem parameters.
%
%    Input, real X0, the initial condition.  When studying many
%    realizations of this problem, it is usual for X0 to be chosen
%    from a normal distribution.
%
%    Input, real TMAX, the final time.
%
%    Input, integer N, the number of large scale time steps.
%
%    Input, integer R, the number of small scale time steps per single
%    large scale time step.
%
%    Input, integer SEED, a seed for the random number generator.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'OU_EULER_MARUYAMA:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use an Euler-Maruyama method to approximate the solution of\n' );
  fprintf ( 1, '  the Ornstein-Uhlenbeck stochastic differential equation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    d x(t) = theta * ( mu - x(t) ) dt + sigma dW\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  with initial condition x(0) = x0.\n' );

  if ( nargin < 1 )
    theta = 2.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default decay rate THETA = %g\n', theta );
  end

  if ( nargin < 2 )
    mu = 1.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default mean MU = %g\n', mu );
  end

  if ( nargin < 3 )
    sigma = 0.15;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default variance SIGMA = %g\n', sigma );
  end

  if ( nargin < 4 )
    x0 = 2.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default initial value X0 = %g\n', x0 );
  end

  if ( nargin < 5 )
    tmax = 3.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default final time TMAX = %g\n', tmax );
  end

  if ( nargin < 6 )
    n = 10000;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default number of large timesteps N = %d\n', n );
  end

  if ( nargin < 7 )
    r = 8;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default ratio of small to large timesteps R = %d\n', r );
  end

  if ( nargin < 8 )
    seed = 123456789;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using default value of random SEED = %d\n', seed );
  end
%
%  Initialize the random number generator.
%  The following way to initialize the random number generator
%  may not be available in older versions of MATLAB.
%
  rng ( seed )
%
%  Set time steps.
%
  dt_large = tmax / n;
  dt_small = tmax / n / r;
%
%  Carry out the Euler-Maruyama approximate integration process.
%
  t = linspace ( 0, tmax, n + 1 );
  x = zeros ( 1, n + 1 );

  x(1) = x0;
  for j = 1 : n
    dw = sqrt ( dt_small ) * randn ( 1, r );
    x(j+1) = x(j) + dt_large * theta * ( mu - x(j) ) + sigma * sum ( dw(1:r) );
  end
%
%  Plot the approximate solution.
%
  plot ( t, x, 'r-', 'LineWidth', 3 )
  xlabel ( 't', 'FontSize', 16 )
  ylabel ( 'X(t)', 'FontSize', 16, 'Rotation', 0, 'HorizontalAlignment', 'right' )
  title ( 'Euler-Maruyama solution of Ornstein-Uhlenbeck SDE', 'FontSize', 16 )
  grid on

  filename = 'ou_euler_maruyama.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function SaveImToTiffStack(IM,tiffdir)
% IM: 1-d cell array containing all image planes

% display each plane and save as tif
h = figure;
for i_plane = 1:length(IM)
    im = IM{i_plane};
    image(im);axis equal; axis off
    drawnow;
    % save tiff
    if (i_plane == 1)
        imwrite(im, tiffdir, 'compression','none','writemode','overwrite')
    else
        imwrite(im, tiffdir, 'compression','none','writemode','append')
    end
    %     pause(0.2)
end
close(h)

author: xiuyechen

repo: FishExplorer


text:
classdef KeyVal < handle
%KEYVAL A class that encapsulates a key and a value, which could be an
%       instance of a handle class. It is assumed that standard comparison
%       operations can be performed on the keys. This class is useful for
%       inserting elements into a BinaryHeap or other data structure that
%       requires a key and a value.
%
%DEPENDENCIES: NONE
%
%Instances of the KeyVal class can be compared to each other by key or to
%other objects that have the same type of key. For example,
%keyVal1<keyVal2 will compare the two objects by their keys.
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    properties
       key
       value
    end
    methods
        function newKeyVal=KeyVal(key,val)
        %KEYVAL The constructor method.
            if(nargin<2)
                val=[];
            end
            if(nargin<1)
                key=[];
            end

            newKeyVal.key=key;
            newKeyVal.value=val;
        end

        function val=lt(obj1,obj2)
        %%LT Less than comparison

            if(isa(obj1,'KeyVal')==true&&isa(obj2,'KeyVal')==true)
                val=obj1.key<obj2.key;
            elseif(isa(obj1,'KeyVal')==true)
                val=obj1.key<obj2;
            else
                val=obj1<obj2.key;
            end
        end

        function val=gt(obj1,obj2)
        %%GT Greater than comparison

            if(isa(obj1,'KeyVal')==true&&isa(obj2,'KeyVal')==true)
                val=obj1.key>obj2.key;
            elseif(isa(obj1,'KeyVal')==true)
                val=obj1.key>obj2;
            else
                val=obj1>obj2.key;
            end
        end

        function val=le(obj1,obj2)
        %%LE Less than or equal to comparison

            if(isa(obj1,'KeyVal')==true&&isa(obj2,'KeyVal')==true)
                val=obj1.key<=obj2.key;
            elseif(isa(obj1,'KeyVal')==true)
                val=obj1.key<=obj2;
            else
                val=obj1<=obj2.key;
            end
        end

        function val=ge(obj1,obj2)
        %%GE Greater than or equal to comparison

            if(isa(obj1,'KeyVal')==true&&isa(obj2,'KeyVal')==true)
                val=obj1.key>=obj2.key;
            elseif(isa(obj1,'KeyVal')==true)
                val=obj1.key>=obj2;
            else
                val=obj1>=obj2.key;
            end
        end

        function val=ne(obj1,obj2)
        %%NE Not equal to comparison

            if(isa(obj1,'KeyVal')==true&&isa(obj2,'KeyVal')==true)
                val=obj1.key~=obj2.key;
            elseif(isa(obj1,'KeyVal')==true)
                val=obj1.key~=obj2;
            else
                val=obj1~=obj2.key;
            end
        end

        function val=eq(obj1,obj2)
        %%EQ Equal to comparison

            if(isa(obj1,'KeyVal')==true&&isa(obj2,'KeyVal')==true)
                val=obj1.key==obj2.key;
            elseif(isa(obj1,'KeyVal')==true)
                val=obj1.key==obj2;
            else
                val=obj1==obj2.key;
            end
        end

        function objCopy=copy(obj)
        %%COPY Create a new KeyVal class where the keys and values are
        %      duplicates of those in this class. This is a shallow copy.
        %      We do not call the copy function on handle class values to
        %      make deep copies.

            objCopy=KeyVal(obj.key,obj.value);
        end
    end
end

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


text:
function [sys,x0,str,ts]=ADRC_Two_order(t,x,u,flag,r0,h,B01,B02,B03,D,b0,c,r1,h1)

switch flag
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes(h);
    case 2
        sys=mdlUpdate(x,u,r0,h,B01,B02,B03,b0,D);
    case 3
        sys=mdlOutputs(x,c,r1,h1,b0);
    case 4,
        sys=mdlGetTimeOfNextVarHit(t,h);
    case {1,9}
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes(h)
    sizes=simsizes;
    sizes.NumContStates=0;
    sizes.NumDiscStates=5;
    sizes.NumOutputs=1;
    sizes.NumInputs=3;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0;0;0;0];
    str=[];
    ts=[h 0];
function sys=mdlUpdate(x,u,r0,h,B01,B02,B03,b0,D)
    e1=x(1)-u(1);
    fh=fhan(e1,x(2),r0,h);
    sys(1)=x(1)+h*x(2);
    sys(2)=x(2)+h*fh;
    e2=x(3)-u(2);
    sys(3)=x(3)+h*(x(4)-B01*e2);
    sys(4)=x(4)+h*(x(5)-B02*fal(e2,0.5,D)+b0*u(3));
    sys(5)=x(5)+h*(-B03*fal(e2,0.25,D));

function sys=mdlOutputs(x,c,r1,h1,b0)
    e3=x(1)-x(3);
    e4=x(2)-x(4);
    sys=-fhan(e3,c*e4,r1,h1)-x(5)/b0;
function sys=mdlGetTimeOfNextVarHit(t,h)
    sys=t+h;

 function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end

function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function [Knots_n] = refine_linear_grid_list(Knots, old_spacing, ds, volsz, varargin)
    k_down = 0.5;
%     new_vol_sz = [];
    if nargin >= 5
        k_down = varargin{1};
    end

    nd = numel(ds);
    if nd == 2
        Knots_sz = ceil(volsz ./ old_spacing) + 1;
        Knots_n = zeros([Knots_sz, 2, size(Knots, 4)]);
        for i = 1 : size(Knots, 4)
            Knots_n(:,:,:, i) = refine_linear_grid_2d(Knots(:,:,:, i), old_spacing, ds, volsz, k_down);
        end
    elseif nd == 3
        Knots_sz = ceil(volsz ./ old_spacing) + 1;
        Knots_n = zeros([Knots_sz, 3, size(Knots, 5)]);
        for i = 1 : size(Knots, 5)
            Knots_n(:,:,:,:, i) = refine_linear_grid_3d(Knots(:,:,:,:, i), old_spacing, ds, volsz, k_down);
        end
    end
end

repo: pTVreg


text:
% Copyright (C) 2000-2002 Paul Kienzle <pkienzle@users.sf.net>
%
% 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/>.
%
% usage:  w = triang (L)
%
% Returns the filter coefficients of a triangular window of length L.
% Unlike the bartlett window, triang does not go to zero at the edges
% of the window.  For odd L, triang(L) is equal to bartlett(L+2) except
% for the zeros at the edges of the window.

function w = triang(L)
  if nargin ~= 1
    help(mfilename);
  elseif ~isscalar(L) || L ~= fix (L) || L < 1
    error('triang: L has to be an integer > 0');
  end % if
  w = 1 - abs ((-(L-1):2:(L-1))' / (L+rem(L,2)));
end

%!error triang
%!error triang(1,2)
%!error triang([1,2]);
%!assert (triang(1), 1)
%!assert (triang(2), [1; 1]/2)
%!assert (triang(3), [1; 2; 1]/2);
%!assert (triang(4), [1; 3; 3; 1]/4);
%!test
%! x = bartlett(5);
%! assert (triang(3), x(2:4));

%!demo
%! subplot(221); axis([-1, 1, 0, 1.3]); grid('on');
%! title('comparison with continuous for odd n');
%! n=7; k=(n-1)/2; t=[-k:0.1:k]/(k+1);
%! plot(t,1-abs(t),';continuous;',[-k:k]/(k+1),triang(n),'g*;discrete;');
%!
%! subplot(222); axis([-1, 1, 0, 1.3]); grid('on');
%! n=8; k=(n-1)/2; t=[-k:0.1:k]/(k+1/2);
%! title('note the higher peak for even n');
%! plot(t,1+1/n-abs(t),';continuous;',[-k:k]/(k+1/2),triang(n),'g*;discrete;');
%!
%! subplot(223); axis; grid('off');
%! title('n odd, triang(n)==bartlett(n+2)');
%! n=7;
%! plot(0:n+1,bartlett(n+2),'g-*;bartlett;',triang(n),'r-+;triang;');
%!
%! subplot(224); axis; grid('off');
%! title('n even, triang(n)!=bartlett(n+2)');
%! n=8;
%! plot(0:n+1,bartlett(n+2),'g-*;bartlett;',triang(n),'r-+;triang;');
%!
%! subplot(111); title('');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function xyz_write_test ( file_out_name )

%*****************************************************************************80
%
%% XYZ_WRITE_TEST tests the XYZ write routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_OUT_NAME, the name of the output file.
%
  point_num = xyz_example_size ( );
%
%  Set the data.
%
  xyz = xyz_example ( point_num );
%
%  Write the data to the file.
%
  xyz_write ( file_out_name, point_num, xyz );

  return
end


text:
%  Matlab (MatConvNet) implementation of our paper
%  DFM: A Performance Baseline for Deep Feature Matching
%  at CVPR 2021 Image Matching Workshop.
%
%  See details at https://github.com/ufukefe/DFM
%
% @authors: ufukefe, kutalmisince
% Created on March 23, 2021
% @Middle East Technical University, Center for Image Analysis
%  Last Edited on July 1, 2021

function [pointsA, pointsB] = DFM_VGG_stage1(imgA, imgB, model, ratios)

    resize_coeff_A = 1;
    resize_coeff_B = 1;
    % Check the image size
    if size(imgA,1)>1600 || size(imgA,2)>1600
        resize_coeff_A = 1600/max(size(imgA,1),size(imgA,2));
        imgA = imresize(imgA,resize_coeff_A);
    end

    if size(imgB,1)>1600 || size(imgB,2)>1600
        resize_coeff_B = 1600/max(size(imgB,1),size(imgB,2));
        imgB = imresize(imgB,resize_coeff_B);
    end

    size_org_A = size(imgA);
    size_org_B = size(imgB);

    % Normalize Images
    imgA = single(imgA);
    imgA = bsxfun(@minus, imgA, model.meta.normalization.averageImage);

    imgB = single(imgB);
    imgB = bsxfun(@minus, imgB, model.meta.normalization.averageImage);

    % zero padding for vgg (canvas should be a multiple of 16)
    imgA = ZeroPadding4VGG(imgA);
    imgB = ZeroPadding4VGG(imgB);

    % Assign layers to be used
    layers_to_use_A = {'conv5_2';'conv4_2';'conv3_2';'conv2_2';'conv1_2'};
    layers_to_use_B = {'conv5_2';'conv4_2';'conv3_2';'conv2_2';'conv1_2'};
    model = ArrangeNetwork(model,layers_to_use_A);

    % Assign ratio tests for layers
    ratios_s1 = ratios(5:-1:1);

    % get activations
    activationsA = GetActivations(imgA, model, layers_to_use_A);
    activationsB = GetActivations(imgB, model, layers_to_use_B);

    % initiate matches
    [pointsA, pointsB] = DenseFeatureMatching(activationsA{1,2}, activationsB{1,2},ratios_s1(1));
    [pointsA, pointsB] = insideImage(pointsA,pointsB,size_org_A,size_org_B,16);

    for k = 2:5
        [pointsA, pointsB] = RefinePoints(pointsA, pointsB, activationsA{k,2}, activationsB{k,2}, ratios_s1(k));
    end

    % Reject matches at the side of the images
    [pointsA,pointsB] = rejectSideMatches(pointsA,pointsB,size_org_A,size_org_B,16);

    if resize_coeff_A ~= 1
       pointsA = (1/resize_coeff_A)*(pointsA-0.5)+0.5;
    end

    if resize_coeff_B ~= 1
       pointsB = (1/resize_coeff_B)*(pointsB-0.5)+0.5;
    end

end

author: ufukefe

repo: DFM

save_path: github-repos/MATLAB/ufukefe-DFM


text:
function [hl,ht,hp,ha] = visualize_rotations(rotorder,rotangle_deg,axhandle)

% FUNCTION
%   [hl,ht,hp,ha] = visualize_rotations(rotorder,rotangle_deg,axhandle)
%
% DESCRIPTION
%   Visualize Euler rotation sequences by displaying rotated axes systems {E_r}
%   and rotation planes, where {E_r} is the orthonormal set of vectors [i_r; j_r; k_r]
%
% INPUT ARGUMENTS
%   rotorder = string containing order of rotations, e.g. 'zxy' or 'xyx'
%   rotangle_deg = vector of rotation angles corresponding to rotorder [degrees]
%   axhandle = optional handle to existing axes
%
% OUTPUT ARGUMENTS
%   hl = 3-by-n+1 matrix of line handles for the x,y,z-axes for E0 to En
%   ht = same for axis text labels
%   hp = 3-by-n matrix of handles to the colored patches
%   ha = handle to annotation
%
% EXAMPLE:
%   visualize_rotations('xzxyz',[25 25 -25 45 10])
%
%   or call the function without input arguments to show an example
%
%
% Author: D.J. van Gerwen
% Created: 22-Feb-2011
% Revised: 13-Jun-2013

% Process input arguments
if nargin<2
    rotangle_deg = [25 35 -35]; % [deg]
end
if nargin<1
    rotorder = 'yzx';
end

% Parameters
lw = 1;
fs = 9;
bgc = 'none';
xyztxt = 'xyz';
xlim = [-1 1];
ylim = xlim;
zlim = xlim;
azel = [135,30]; % Azimuth and elevation of view

% Convert angles
rotangle_rad = rotangle_deg/180*pi; % [-]

% Define rotation matrices for {E_i+1}=[R]{E_i} as inline functions (based on right-hand rule)
Rx = inline('[1 0 0; 0 cos(x) sin(x); 0 -sin(x) cos(x)]'); % Rotation about x-axis by angle x
Ry = inline('[cos(y) 0 -sin(y); 0 1 0; sin(y) 0 cos(y)]'); % Rotation about y-axis by angle y
Rz = inline('[cos(z) sin(z) 0; -sin(z) cos(z) 0; 0 0 1]'); % Rotation about z-axis by angle z

% Define global reference system
E = [1 0 0;  % i0
     0 1 0;  % j0
     0 0 1]; % k0 (could use eye(3), but that is less informative)

% Apply transformations: {E1}=[R1]{E0}, {E2}=[R2]{E1}=[R2][R1]{E0}, etc
for r = 1:length(rotorder)
    switch rotorder(r)
        case 'x'
            % Evaluate rotation matrix
            Rr(:,:,r) = Rx(rotangle_rad(r));
            % Define patch
            pv(:,:,r) = [[0 1 0].*Rr(2,:,r); Rr(2,:,r); [0 0 1].*Rr(2,:,r);
                [0 0 1].*Rr(3,:,r); Rr(3,:,r); [0 1 0].*Rr(3,:,r)]; % NOTE: to show triangles, set the third and sixth row to zero
        case 'y'
            % Evaluate rotation matrix
            Rr(:,:,r) = Ry(rotangle_rad(r));
            % Define patch
            pv(:,:,r) = [[1 0 0].*Rr(1,:,r); Rr(1,:,r); [0 0 1].*Rr(1,:,r);
                [0 0 1].*Rr(3,:,r); Rr(3,:,r); [1 0 0].*Rr(3,:,r)]; % NOTE: to show triangles, set the third and sixth row to zero
        case 'z'
            % Evaluate rotation matrix
            Rr(:,:,r) = Rz(rotangle_rad(r));
            % Define patch
            pv(:,:,r) = [[1 0 0].*Rr(1,:,r); Rr(1,:,r); [0 1 0].*Rr(1,:,r);
                [0 1 0].*Rr(2,:,r); Rr(2,:,r); [1 0 0].*Rr(2,:,r)]; % NOTE: to show triangles, set the third and sixth row to zero
    end
    % Transformation to E0
    E(:,:,r+1) = Rr(:,:,r)*E(:,:,r);
end

% Initialize figure
if ~exist('axhandle','var')
    figure('color','white','numbertitle','off','name','Visualize Euler rotation sequence')
    set(gca,'xlim',xlim,'ylim',ylim,'zlim',zlim,'projection','orthographic') % NOTE: projection is either 'orthographic' or 'perspective'
    view(azel)
    axis off
    title(sprintf(['Rotation sequence: %s (%3.0f\\circ' repmat(',%3.0f\\circ',1,length(rotorder)-1) ')'],rotorder,rotangle_deg))
else
    axes(axhandle)
end
%cmp = colormap(lines);
cmp = [1 0 0; 0 0.5 0; 0 0 1; 1 1 0; 0 1 1; 1 0 1];

% Display global reference system {E0} and rotated reference systems {E.}
for r = 1:length(rotorder)+1
    for i = 1:3
        % Axis and label
        hl(i,r) = line([0 E(i,1,r)],[0 E(i,2,r)],[0 E(i,3,r)],'color','k','linewidth',lw);
        ht(i,r) = text('position',E(i,:,r),'string',sprintf('%s_%u',xyztxt(i),r-1),'fontsize',fs,'backgroundcolor',bgc);
        % Create patches
        if r<=length(rotorder)
            pvr = pv(:,:,r);
            for rr = r-1:-1:1
                pvr = pvr*Rr(:,:,rr); % Transform patch vector from local to global reference
            end
            hp(i,r) = patch(pvr(:,1),pvr(:,2),pvr(:,3),cmp(r,:),'edgecolor','none','facealpha',0.1);
        end
    end
end

% Correct axis labels for the rotation axes
for i = 1:3
    for r = strfind(rotorder,xyztxt(i))
        str1 = get(ht(i,r),'string');
        str2 = get(ht(i,r+1),'string');
        set(ht(i,r+1),'string',[str1 ',' str2])
    end
end

% Text
ha = annotation('textbox',[0.4 0 0.6 0.14],'string',{'Colored patches represent components of the local rotation matrix in the plane of rotation. Definitions are according to right-hand rule.'},'edgecolor','none');


text:
dir = 'C:\Users\xiuye\Dropbox\ZBrainOutput\light-sheet data (11 fish)\cell locations';

Collection = cell(1,11);
for i_fish = 1:11,
    i_fish
    fishdir = fullfile(dir,['F',num2str(i_fish),'_cell_info.mat']);
    load(fishdir);
    XY = reshape([cell_info.center],2,[]);
    Z = [cell_info.slice];
    m = vertcat(XY,Z)';
    text = cell(length(m),1);
    for i = 1:length(m),
        text{i} = [num2str(m(i,1)),' ',num2str(m(i,2)),' ',num2str(m(i,3))];
    end
    Collection{i_fish} = text;
end

save(fullfile(dir,'Allfish.mat'),'Collection');

%%

for i_fish = 1:11,
    i_fish
    %%
    fileID = fopen(['F',num2str(i_fish),'_XYZ.txt'],'w');

    fishdir = fullfile(dir,['F',num2str(i_fish),'_cell_info.mat']);
    load(fishdir);
    XY = reshape([cell_info.center],2,[]);
    Z = [cell_info.slice];
    m = vertcat(XY,Z)';
    for i = 1:length(m),
       s = [num2str(m(i,2)),' ',num2str(m(i,1)),' ',num2str(m(i,3))];
       fprintf(fileID,[s,'\r\n']);
    end

    fclose(fileID);
    %%
end

%%

for i_fish = 1:11,
    fileID = fopen(['F',num2str(i_fish),'_XYZ.txt'],'w');
    text =Collection{i_fish};
    for i = 1:length(text);
        fprintf(fileID,[text{i},'\r\n']);
    end
    fclose(fileID);
    %%
end


text:
function [tSeries, tSerr] = meanTSeries(vw, scanNum, roi, getRawData)
% Computes the mean tSeries, averaged over ROIcoords, for scanNum.
%
%   tSeries = meanTSeries(vw, scanNum, [roi], [getRawData])
%
% scanNum: scan number
% roi:     roi specification (ROI coords, name of loaded ROI, or ROI
%           struct). It can also be a cell array, in which case the
%           returned tSeries will also be a cell array, one entry for each
%           ROI.
% getRawData: if true, then do not detrend and do not convert to % signal
%
% djh, 7/98
% djh, 2/2001
% - updated to mrLoadRet-3.0
% - updated to work for gray and flat views
% bw/ab 6/2003
% - wrote getSlicesROI, ROIcoords2cellArray and inserted here so
%   we could use those routines in other places, too.
% Ress, 4/05 One fix and a new feature:
% Fix: Added code to properly ignore bad-flagged (NaN) values when taking
%       the mean. Feature: Now returns spatial SEM as an optional second
%       output argument; note that this estimate needs to be scaled up by a
%       voxel-size dependent spatial correlation factor
if notDefined('roi'),        roi         = vw.selectedROI; end
if notDefined('getRawData'), getRawData  = false;            end

% make sure the ROI is properly specified
roi = tc_roiStruct(vw, roi);

% if the ROI is empty, don't bother with the rest of this function
if length(roi) == 1 && isempty(roi.coords)
    tSeries = [];
    tSerr = [];
    return
end

% get the set of ROI coords as a cell array
for n = 1:length(roi)
    ROIcoords{n} = roi(n).coords;
end

% Find the slice indices for this collection of ROIs
sliceInds = getSlicesROI(vw, ROIcoords);

nROIs = length(ROIcoords);

tSeries = cell(1,nROIs);
tSerr   = cell(1, nROIs);

nFrames               = viewGet(vw, 'numFrames', scanNum);
detrend               = detrendFlag(vw,scanNum);
inhomoCorrection      = inhomoCorrectionFlag(vw, scanNum);
temporalNormalization = temporalNormalizationFlag(vw, scanNum);
smoothFrames          = detrendFrames(vw,scanNum);

if getRawData
    detrend = false;
    inhomoCorrection      = false;
    temporalNormalization = false;
    %smoothFrames = 0;
end

% Take first pass through ROIs to see which slices to load
switch vw.viewType
    case {'Inplane' 'Flat'}
        sliceInds = [];
        for r=1:nROIs
            if isempty(ROIcoords{r})
                disp('MeanTSeries: ignoring empty ROI in this slice')
            else
                sliceInds = [sliceInds, ROIcoords{r}(3,:)];
            end
        end
        sliceInds = unique(sliceInds);
    case {'Gray' 'Volume'}
        sliceInds = 1;
    otherwise
        myErrorDlg('meanTSeries: Only for Inplane, Gray, or Flat views.');
end

% Load tSeries & divide by mean, but don't detrend yet.
% Otherwise, detrending the entire tSeries is much slower. DJH
detrendNow    = false;
noMeanRemoval = false;
if getRawData,	noMeanRemoval = true; end % just get the raw data!

vw = percentTSeries(vw,scanNum,sliceInds,detrendNow,inhomoCorrection,temporalNormalization,noMeanRemoval);

for r=1:nROIs
    % Extract time-series
    subtSeries = getTSeriesROI(vw,ROIcoords{r});
    if ~isempty(subtSeries)

        % Detrend now (faster to do it now after extracting subtSeries for a small subset of the voxels)
        subtSeries = detrendTSeries(subtSeries, detrend, smoothFrames);
        % Add 'em up
        if isempty(tSeries{r})
            numPts{r} = sum(isfinite(subtSeries), 2);
            subtSeries(~isfinite(subtSeries)) = 0;
            tSerr{r} = sum(subtSeries.^2, 2);
            tSeries{r} = sum(subtSeries, 2);
        else
            numPts{r} = numPts{r} + sum(isfinite(subtSeries), 2);
            subtSeries(~isfinite(subtSeries)) = 0;
            tSerr{r} = tSerr{r} + sum(subtSeries.^2, 2);
            tSeries{r} = tSeries{r} + sum(subtSeries,2);
        end
    end
end

% Final pass through ROIs to turn sum into mean
for r=1:nROIs
    if isempty(numPts{r})
        tSeries{r} = zeros(nFrames,1);
        tSerr{r} = zeros(nFrames, 1);
    else
        tSeries{r} = tSeries{r} ./ numPts{r};
        tSerr{r} = sqrt((tSerr{r} ./ numPts{r}.^2) - tSeries{r}.^2./numPts{r});
    end
end

if (nROIs==1)
    tSeries = tSeries{1};
    tSerr = tSerr{1};
end

% Clean up (because didn't detrend the whole tSeries properly)
vw.tSeries=[];
vw.tSeriesScan=NaN;
vw.tSeriesSlice=NaN;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the first column of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first column of W corresponds to the parameters for the bias unit
%

% =========================================================================

end

author: vkosuri

repo: CourseraMachineLearning


text:
function varargout = waterfall(varargin)
%WATERFALL   Waterfall plot of a SPHEREFUN.
%   WATERFALL(F) displays the waterfall plot of F.
%
%   WATERFALL(F, S) displays the column and row chebfuns of F that are used for
%   its approximation.  This is a 3D version of plot(f,S), where S is a string
%   (see PLOT).
%
%   WATERFALL(F, S, 'nslices', N) displays the first min(N,length(f)) column
%   and rows.
%
%   WATERFALL supports passing options to the plot, similar to standard Matlab
%   plot commands. The options supported are:
%       'color':      Color of lines and markers plotted.
%       'marker':     Marker for pivot points.
%       'markersize': Size of markers plotted at pivot points.
%
%   H = WATERFALL(...) returns a handle to a waterfall plot object.
%
% See also SPHEREFUN/PLOT.

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

[varargout{1:nargout}] = waterfall@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [px, py, gx, gy] = position_7_link(z,P)
%[px,py,gx,gy] = POSITION_7_LINK(Z,P)
%
%FUNCTION:  This function computes the cartesian positions
%    of the CoM and tip of each link
%INPUTS:
%
%
%OUTPUTS:
%    px = [nLink X nTime] position of the end of each link
%    py = [nLink X nTime] position of the end of each link
%    gx = [nLink X nTime] position of the CoM of each link
%    gy = [nLink X nTime] position of the CoM of each link
%
%NOTES:
%    This file was automatically generated by writePosition.m

g  = P.g ; %gravity
m1 = P.m(1); % Link 1 mass
m2 = P.m(2); % Link 2 mass
m3 = P.m(3); % Link 3 mass
m4 = P.m(4); % Link 4 mass
m5 = P.m(5); % Link 5 mass
m6 = P.m(6); % Link 6 mass
m7 = P.m(7); % Link 7 mass
l1 = P.l(1); % Link 1 length
l2 = P.l(2); % Link 2 length
l3 = P.l(3); % Link 3 length
l4 = P.l(4); % Link 4 length
l5 = P.l(5); % Link 5 length
l6 = P.l(6); % Link 6 length
l7 = P.l(7); % Link 7 length
I1 = P.I(1); % Link 1 moment of inertia about its center of mass
I2 = P.I(2); % Link 2 moment of inertia about its center of mass
I3 = P.I(3); % Link 3 moment of inertia about its center of mass
I4 = P.I(4); % Link 4 moment of inertia about its center of mass
I5 = P.I(5); % Link 5 moment of inertia about its center of mass
I6 = P.I(6); % Link 6 moment of inertia about its center of mass
I7 = P.I(7); % Link 7 moment of inertia about its center of mass
d1 = P.d(1); % Link 1 distance between center of mass and parent joint
d2 = P.d(2); % Link 2 distance between center of mass and parent joint
d3 = P.d(3); % Link 3 distance between center of mass and parent joint
d4 = P.d(4); % Link 4 distance between center of mass and parent joint
d5 = P.d(5); % Link 5 distance between center of mass and parent joint
d6 = P.d(6); % Link 6 distance between center of mass and parent joint
d7 = P.d(7); % Link 7 distance between center of mass and parent joint

th1 = z(1,:);
th2 = z(2,:);
th3 = z(3,:);
th4 = z(4,:);
th5 = z(5,:);
th6 = z(6,:);
th7 = z(7,:);

nTime = length(th1);
px = zeros(7,nTime);
py = zeros(7,nTime);
gx = zeros(7,nTime);
gy = zeros(7,nTime);

px(1,:) = l1.*cos(th1);
py(1,:) = l1.*sin(th1);
gx(1,:) = d1.*cos(th1);
gy(1,:) = d1.*sin(th1);

px(2,:) = l1.*cos(th1) + l2.*cos(th2);
py(2,:) = l1.*sin(th1) + l2.*sin(th2);
gx(2,:) = d2.*cos(th2) + l1.*cos(th1);
gy(2,:) = d2.*sin(th2) + l1.*sin(th1);

px(3,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3);
py(3,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3);
gx(3,:) = d3.*cos(th3) + l1.*cos(th1) + l2.*cos(th2);
gy(3,:) = d3.*sin(th3) + l1.*sin(th1) + l2.*sin(th2);

px(4,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4);
py(4,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4);
gx(4,:) = d4.*cos(th4) + l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3);
gy(4,:) = d4.*sin(th4) + l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3);

px(5,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4) + l5.*cos(th5);
py(5,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4) + l5.*sin(th5);
gx(5,:) = d5.*cos(th5) + l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4);
gy(5,:) = d5.*sin(th5) + l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4);

px(6,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4) + l5.*cos(th5) + l6.*cos(th6);
py(6,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4) + l5.*sin(th5) + l6.*sin(th6);
gx(6,:) = d6.*cos(th6) + l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4) + l5.*cos(th5);
gy(6,:) = d6.*sin(th6) + l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4) + l5.*sin(th5);

px(7,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4) + l5.*cos(th5) + l6.*cos(th6) + l7.*cos(th7);
py(7,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4) + l5.*sin(th5) + l6.*sin(th6) + l7.*sin(th7);
gx(7,:) = d7.*cos(th7) + l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4) + l5.*cos(th5) + l6.*cos(th6);
gy(7,:) = d7.*sin(th7) + l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4) + l5.*sin(th5) + l6.*sin(th6);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%%% void convert_video4d_to_avi(V,filename)
%
function convert_video4d_to_avi(video,filename)
  movobj(1:size(video,4)) = struct('cdata',[], 'colormap',[]);

  for i = 1:size(video,4)
    frame = video(:,:,:,i);

    if(size(frame,3) == 1)
      frame = repmat(frame,[1 1 3]);
    end

    movobj(i).cdata = frame;
  end

  movie2avi(movobj, filename, 'compression', 'None');
end

author: andrewssobral

repo: lrslibrary