text:
function varargout = maindalta(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @untitled1_OpeningFcn, ...
                   'gui_OutputFcn',  @untitled1_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before untitled1 is made visible.
function untitled1_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to untitled1 (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = untitled1_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
global a;
handles.text7.String = 'Frame t';
handles.text6.String = 'Threshold c';
handles.text2.String = 'image f';
set(handles.edit1,'string',num2str(get(hObject,'value')));

dalta = get(handles.slider1,'Value');

[I,blur,et,newEt] = generalclean(dalta);
imshow(blur,'Parent',handles.axes1)
imshow(mat2gray(newEt),'Parent',handles.axes2)
imshow((I),'Parent',handles.axes3)

ts = get(handles.slider3,'Value');
[It,et] = reconstructAvideo(dalta,ts);
%
imshow((It),'Parent',handles.axes4)
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)

if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function edit1_Callback(hObject, eventdata, handles)

global a;
handles.text7.String = 'Frame t';
handles.text6.String = 'Threshold c';
handles.text2.String = 'image f';
set(handles.slider1,'value',str2double(get(hObject,'string')));
sliderValue = get(handles.edit1,'String');
ang = str2double(sliderValue);
set(handles.slider1,'Value',ang);

dalta = get(handles.slider1,'Value');

[I,blur,et,newEt] = generalclean(dalta);
imshow(blur,'Parent',handles.axes1)
imshow(mat2gray(newEt),'Parent',handles.axes2)
imshow((I),'Parent',handles.axes3)

ts = get(handles.slider3,'Value');
[It,et] = reconstructAvideo(dalta,ts);
%
imshow((It),'Parent',handles.axes4)
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
%
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on slider movement.
function slider3_Callback(hObject, eventdata, handles)
global b;
handles.text7.String = 'Frame t';
handles.text6.String = 'Threshold c';
handles.text2.String = 'image f';
set(handles.edit2,'string',num2str(get(hObject,'value')));
dalta = get(handles.slider1,'Value');
ts = get(handles.slider3,'Value');
[It,Et] = reconstructAvideo(dalta,ts);%_old

imshow(mat2gray((Et)),'Parent',handles.axes2)
imshow(((It)),'Parent',handles.axes4)

% --- Executes during object creation, after setting all properties.
function slider3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function edit2_Callback(hObject, eventdata, handles)
global b;
handles.text7.String = 'Frame t';
handles.text6.String = 'Threshold c';
handles.text2.String = 'image f';
set(handles.slider3,'value',str2double(get(hObject,'string')));
sliderValue = get(handles.edit2,'String');
ang = str2double(sliderValue);
set(handles.slider3,'Value',ang);

dalta = get(handles.slider1,'Value');
ts = get(handles.slider3,'Value');
[It,Et] = reconstructAvideo(dalta,ts);%_old

imshow(mat2gray((Et)),'Parent',handles.axes2)
imshow(((It)),'Parent',handles.axes4)

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: panpanfei


text:
% Test file for ADCHEBFUN trigonometric and related functions, part 3/4

function pass = test_trig3

% List of trigonometric functions to test.
funcList = {@atand, @atanh, @cos, @cosd, @cosh, @cot, @cotd, @coth, ...
    @csc, @cscd};

% Call the ADCHEBFUN testUnary() method to do the tests.
pass = adchebfun.testUnary(funcList);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%  Figure 3.35      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.35
clf;
t=0:.05:10;
k=1;
ki=0;
numCL=[k ki];
denCL=[1 3 2+k ki];
y1=step(numCL,denCL,t);
ki=1;
numCL=[k ki];
denCL=[1 3 2+k ki];
y2=step(numCL,denCL,t);
k=10;
ki=5;
numCL=[k ki];
denCL=[1 3 2+k ki];
y3=step(numCL,denCL,t);
axis([0 10 0 1.2]);
plot(t,y1,t,y2,t,y3);
grid;
xlabel('Time (sec)');
ylabel('y(t)');
title('Fig. 3.35 Transient responses');
axis('normal');
text(2,.25,'K=1, K_I=0');
text(3,0.75,'K=1, K_I=1');
text(1.5,1.1,'K=10, K_I=5');

author: Sable

repo: mcbench-benchmarks


text:
%  Figure 10.23      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_23.m is a script to generate Fig. 10.23,
%  the root locus of the LQR symmetric root locus compensator of the
%  satellite position control, non-collocated case WITH ESTIMATOR

% Parameter values
m=[1, .1]; k=[0, .091] ; d=[0, .0036]; k1=[0, 0.4];

% call model
[f,g,h,j] = twomass(m,k,d);

% compute feedforward values
s=[f, g;h, 0];
r=[0*g;1];
n=s\r;nx=n(1:4);
nu=n(5);

% call model
[f1,g,h,j] = twomass(m,k1,d);

% form G(s)G(-s) in state-space
a=[f, 0*f;
-h'*h, -f'];
b=[g;0*g];
c=[0*h, g'];
d=[0];
hold off; clf
P=eig(a-b*c*0.1621);
pc=P(real(P<0)==1);
K=place(f,g,pc);
nbar=nu+K*nx;
% eig(f-g*K)
P=eig(a-b*c*3.056e7);
pe=P(real(P<0)==1);
L=place(f',h',pe)';
ac=f-g*K-L*h ;bc=L;cc=K;dc=0;
[Aol,Bol,Col,Dol]=series(ac,bc,cc,dc,f,g,h,j);
[acl,bcl,ccl,dcl]=feedback(f,g,h,j,ac,bc,cc,dc);
[acl1,bcl,ccl,dcl]=feedback(f1,g,h,j,ac,bc,cc,dc);
bcl= nbar*[g;g];
rlocus(Aol,Bol,Col,Dol);
grid;
v=[-2, 2, -1.5, 1.5]
axis(v);
title('Fig. 10.23 Root locus for D_4(s)G(s).')


text:
function [f, df] = func_split_spcp(x,Y,params,errFcn)
% [f, df] = func_split_spcp(x,Y,params,errFunc)
% [errHist] = func_split_spcp();
% [S] = func_split_spcp(x,Y,params,'S');
%
% Compute function and gradient of split-SPCP objective
%
%   lambda_L/2 (||U||_F^2 + ||V||_F^2) + phi(U,V)
%
% where
%
%   phi(U,V) = min_S .5|| U*V' + S - Y ||_F^2 + lambda_S ||S||_1.

persistent errHist

if nargin==0
   f = errHist;
   errHist = [];
   return;
end

if nargin<3, params=[]; end
m   = params.m;
n   = params.n;
k   = params.k;
lambdaL = params.lambdaL;
lambdaS = params.lambdaS;
useGPU  = params.gpu;

U = reshape(x(1:m*k),m,k);
V = reshape(x(m*k+1:m*k+n*k),n,k);

L   = U*V';
LY = vec(Y-L);

soft_thresh  = @(LY,lambdaS) sign(LY).*max(abs(LY) - lambdaS,0);
S = soft_thresh(LY,lambdaS);
if nargout==1
    f = S;
    return;
end
SLY = reshape(S-LY,m,n);

fS = norm(SLY,'fro')^2/2 + lambdaS*norm(S(:),1);

f = lambdaL/2*(norm(U,'fro')^2 + norm(V,'fro')^2) + fS;

if nargout > 1

    if useGPU
        df = gpuArray.zeros(m*k + n*k,1);
    else
        df = zeros(m*k + n*k,1);
    end

    df(1:m*k) = vec(lambdaL*U) + vec((SLY)*V);
    df(m*k+1:m*k+n*k) = vec(lambdaL*V) + vec((SLY)'*U);

end

errHist(end+1,1) = toc;
errHist(end,2) = gather(f);

if ~isempty(errFcn)
    errHist(end,3) = errFcn(x);
end

tic;

end

function out = setOpts( params, field, default )
    if ~isfield( params, field )
        params.(field)    = default;
    end
    out = params.(field);
    params    = rmfield( params, field ); % so we can do a check later
end

author: stephenbeckr

repo: fastRPCA


text:
function [ x, w ] = ncoh_rule ( n )

%*****************************************************************************80
%
%% NCOH_RULE: Newton Cotes Open Half quadrature rule.
%
%  Discussion:
%
%    Newton Cotes Open Half rule on [-1,+1];
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), W(N), the points and weights.
%
  x = ( 1 : 2 : 2*n-1 )';
  x = x / ( 2 * n );
%
%  Linear transformation from [0,1] to [-1,+1].
%
  x = 2.0 * x - 1.0;
%
%  Compute the weights.
%
  x_min = -1.0;
  x_max = 1.0;

  w = nc_compute ( n, x_min, x_max, x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [b,t,p,hh] = ancova(groups,x,y,varargin)
% :Usage:
% ::
%
%     b,t,p,pthandles] = ancova(groups,x,y,[plot],[covs of no interest])
%
% :Outputs:
%
%   **Elements of b, t, p:**
%        1st = intercept, 2nd = group effect, 3rd = slope, 4th = grp x slope
%        interaction
%
% recursive -- call ancova repeately if y is a matrix
% get pairwise standardized slopes (corrs)
% and group diffs and slope interaction for all
% pairs of y vectors
%

hh = [];
doplot = 0; if length(varargin) > 0, doplot = varargin{1};,end

% covariates of no interest
if length(varargin) > 1,
    covs = varargin{2};,

    [x,y,r,p,rrob,prob] = partialcor([x covs],y,1);
    % doesn't adjust groups now...

end

if size(y,2) > 1 & isempty(x)

    for i = 1:size(y,2)-1
        for j = i+1:size(y,2)

            [b,t,p] = ancova(groups,y(:,i),y(:,j));
            bb{1}(i,j) = b(2); tt{1}(i,j) = t(2); pp{1}(i,j) = p(2);  % group effect = cell 1
            bb{2}(i,j) = b(3); tt{2}(i,j) = t(3); pp{2}(i,j) = p(3);  % slope = cell 2
            bb{3}(i,j) = b(4); tt{3}(i,j) = t(4); pp{3}(i,j) = p(4);  % grp x slope = cell 3

        end
    end

    % clean up last row
    for i = 1:length(bb)
        bb{i}(end+1,:) = 0;  bb{i} = bb{i} + bb{i}';
        tt{i}(end+1,:) = 0;  tt{i} = tt{i} + tt{i}';
        pp{i}(end+1,:) = 0;  pp{i} = pp{i} + pp{i}';
    end

    b = bb;
    t = tt;
    p = pp;
    return

end

% -----------------------------------------
% basic ancova program
% -----------------------------------------

% center and format
x = [scale(groups,1) scale(x,1)];
%y = scale(y);

% add interaction in slopes (diff. slopes model)

x(:,3) = scale(x(:,1).*x(:,2));

[b,dev,stats]=glmfit(x,y);
bi = b(end); ti = stats.t(end); ppi = stats.p(end);

if ppi > .05

    % no sig. interaction, drop interaction and use parallel slopes

    [b,dev,stats]=glmfit(x(:,1:2),y);
    b(end+1) = bi; stats.t(end+1) = ti; stats.p(end+1) = ppi;

end

t = stats.t;
p = stats.p;

% -----------------------------------------
%  ancova plot
% -----------------------------------------

if doplot
    hh = [];
    colors = {'yo' 'm^'};
    [uni,b,grps] = unique(groups);
    legstr = {['Group 1: ' num2str(uni(1))] ['Group 2: ' num2str(uni(2))]};
        uni = 1:max(grps);
    % median split, if continuous, high then low
    if length(uni) > 2,
        uni = [1 2]; grps = grps.*0;
        grps(groups>median(groups)) = 1;
        grps(groups<median(groups)) = 2;
        legstr = {'High' 'Low'};
    end
    %tor_fig;
    for i = 1:length(uni)   % for each group, make a plot

        [tmp,tmp,tmp,hh(i)] = plot_correlation_samefig(x(find(grps==uni(i)),2),y(find(grps==uni(i))),[],colors{i},0,1);
    end
    xlabel('x'); ylabel('y');
    legend(hh,legstr)
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function f = ldivide(f, c, varargin)
%.\   Left array divide for a TRIGTECH.
%   C .\ F is equivalent to F ./ C.
%
% See also MLDIVIDE, RDIVIDE, TIMES.

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

f = rdivide(c, f, varargin{:});

end


text:
clear
close all
clc;

%time
dt=0.005;

%coeffs
Cg=200.;
Cd=0;
minSpeed=10;

name='greencone.tiff'
sanitizedname=[strrep(name,'.','') '_static_shortstep']
mkdir([sanitizedname 'outputdir'])
%green
%greenheight=double(imread('testgreen.tiff'))/256*1.7;
%greenheight=double(imread('greencone.tiff'))/256*.3;
%greenheight=double(imread('greenconeinv.tiff'))/256*.3;
greenheight=-fliplr(rot90(rot90(double(imread(name))/256*.3))).^-1 * 1000;
greenheightPLOT=max(-fliplr(rot90(rot90(double(imread(name))/256*.3))).^-1 * 500,-100);
[gradX,gradY]=gradient(greenheight(:,:,1));
forceFun=@(R,RL) (RL-R)/dt*Cd + [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))]*Cg;
%forceFun=@(R,RL) [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))];
l=size(greenheight,1);

%pin
holeLoc=size(greenheight)'/2;
holeRadius=3;

%graphics
%sph=sphere(30)
alt=30.85+20;
az=-39+20;
cmapgreen=flipud([181,228,138; 160,220,104; 124,216,87; 81,212,76; 52,188,67; 37,167,60; 32,154,61; 1,114,56; 0,86,19]/255);
stickX=[holeLoc(1),holeLoc(1)];
stickY=[holeLoc(2),holeLoc(2)];
stickZ=[greenheightPLOT(floor(holeLoc(1)),floor(holeLoc(1))),greenheightPLOT(floor(holeLoc(1)),floor(holeLoc(1)))+100];
flagX=[stickX;stickX+40];
flagY=[stickY;stickY];
flagZ=[stickZ(2),stickZ(2)-30;stickZ(2),stickZ(2)-30];
flagC=zeros(2,2,3);
flagC(:,:,1)=.9;

% balls
startLoc=[64;64]*3;

%angles=(45)*pi/180;
%speeds=400;
%angles=linspace((45-15)*pi/180,(45+15)*pi/180,10);%-pi/4;
angles=linspace((-45-20)*pi/180,(-45+20)*pi/180,5);%-pi/4;
speeds=linspace(60,80,20);

numParticles=length(speeds)*length(angles);

speeds=reshape(speeds,1,1,length(speeds));
angleVectors=[cos(angles);sin(angles)];
startConditions=bsxfun(@times,speeds,angleVectors);
startConditions=reshape(startConditions,2,length(speeds)*length(angles));

r=repmat(startLoc,1,numParticles);
rl=r-startConditions*dt;

rstart=r;
rlstart=rl;

%sequentialhits staging
%each ball laucnhes after 1/4 second (300 iters/4) = 80 ticks
launchTicks=80

i=0;
while sum(sum(r~=rl))
    %energy verlet + drag
    i=i+1;
    rn=2*r-rl+(forceFun(r,rl))*dt^2;
    rl=r;
    r=rn;

    %%%LAUNCH CONDITIONS
    if i<50
        r=rstart;
        rl=rlstart;
    end

    %%%STOP CONDITIONS
    %static friction
    dr=rl-r;
    s=sqrt(dr(1,:).^2+dr(2,:).^2)/dt;
    %walls
    haltBallsEdge=(r>l)|(r<1);
    %in the hole
    distToHole=bsxfun(@minus,holeLoc,rl);%use rl not r to check if the ball is in the hole so it can't escape
    distToHole=sqrt(distToHole(1,:).^2+distToHole(2,:).^2);
    %stop those in stop conditions
    holeBalls=(distToHole<holeRadius);
    rl(:,holeBalls)=repmat(holeLoc,1,sum(holeBalls));
    haltBalls=(haltBallsEdge(1,:)|haltBallsEdge(2,:))|(s<minSpeed)|(distToHole<holeRadius);
    r(:,haltBalls)=rl(:,haltBalls);
        if mod(i,20)==0
        figure(1)
        surf(greenheightPLOT,'edgecolor','none')
        colormap(cmapgreen)
        set(gca,'DataAspectRatio',[1 1 1])
        hold on
        scatter3(r(1,:),r(2,:),3+greenheightPLOT(sub2ind(size(greenheightPLOT), floor(r(2,:)),floor(r(1,:)))),'wo','filled')
        plot3(stickX,stickY,stickZ,'k-')
        surf(flagX,flagY,flagZ,flagC,'edgecolor','none')
        hold off
        %az=az+.1;
        view([az alt])
%         imshow(greenheight,[min(min(greenheight)),max(max(greenheight))],'colormap',colormap('parula'))
        hold on
        %rectangle('Position',[holeLoc(1)-holeRadius holeLoc(2)-holeRadius holeRadius*2 holeRadius*2],'Curvature',[1,1],'FaceColor','black');
        %scatter(r(1,:),r(2,:),'w.')
        hold off
        xlim([0,512])
        ylim([0,512])
        zlim([-100,0])
        %title(num2str(s(1)))
        set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
        saveas(gcf,[[sanitizedname 'outputdir'] '\' num2str(i/20,'%05i') '.png'])
    end
end

%phase space

figure(2)
maxDist=256;
maxColor=reshape([1,0,0],1,1,3);
minColor=reshape([1,1,1],1,1,3);
dists=fliplr(flipud(reshape(distToHole,length(angles),length(speeds))'));
imshow(dists,[holeRadius,maxDist],'colormap',colormap('jet'))
hold on
h=imshow(repmat(ones(size(dists)),1,1,3));
hold off
set(h,'AlphaData',dists<holeRadius);
ylabel('Putt speed')
xlabel('putt angle')

author: BrianHaidet

repo: AlphaPhoenix


text:
% reshape() - reshape of memory mapped underlying array
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function obj = reshape(obj,d1,d2,d3)

    % decode length
    % -------------
    if nargin > 3
        d1 = [ d1 d2 d3 ];
    elseif nargin > 2
        d1 = [ d1 d2 ];
    end;

    if prod(size(obj)) ~= prod(d1)
        error('Wrong dimensions for reshaping');
    end;

    if obj.transposed
        d1 = [d1(2:end) d1(1)];
    end;

    obj.dimensions = d1;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
filename='RVE_Square_Triangle';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'enforceCh_CCstar_L2','perimeter'};%enforceCh_CCstar_L2
weights=[1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;%mma CON 1E-4
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.6;
Perimeter_target=3;
optimality_final =1e-4;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
selectiveC_Cstar = 'Vfrac06';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Xco Tau MT] = HC_explicit_nh (Tin,TC,DE,HC,SE,DI,TS,TT,TBC,BC)
%
% DESCRIPTION:
%
% Explicit numerical method for one dimensional unsteady state
% heat transfer by conduction for non-homogenous material
%
% Basic PDE equation :
%
%     dT         d^2 T
%   ------ = TD ------             for 0 < x < DI  and  tau > 0
%    dtau        dx^2
%
% Principle of method :
%       ^
%  time |
%       |      m-1     m     m+1
%       |       |      |      |
%       |    ---+------x------+-----   p+1
%       |  dtau |      |      |
%       |    ---o------o------o-----   p
%       |  dtau |      |      |
%       |    ---+------+------+-----   p-1
%       |       |  dx  |  dx  |
%       |
%       +-------------------------------->
%                                        x
% Space description:
%                  DI
%   |<------------------------------>|
%   +   +   +   +   +  . . .     +   +
%   1   2   3   4   5           n-1  n
%   |<->|
%     SE
%
% Difference equation :
%
%   T(m,p+1) - T(m,p)      TD
%  ------------------- = ------- (T(m-1,p) - 2.T(m,p) + T(m+1,p))
%          TS             SE.SE
%
%       TD.TS
%  M = -------- <= 0,5
%       SE.SE
%
%  T(m,p+1) = M.T(m-1,p) + (1 - 2.M)T(m,p) + M.T(m+1,p)
%
%  T(m,p+1) = M(m-1).T(m-1,p) + (1 - (M(m-1)+M(m)))T(m,p) + M(m).T(m+1,p)
%  for m <2;n-1>
%
%  First-type boundary condition:
%    - temperatures T(1,p) and T(n,p) are given
%    - temperatures T(m,p) for m = 2, n-1 :
%      T(m,p+1) = M (T(m-1,p) + T(m+1,p)) + T(m,p).(1 - 2.M)
%
%  Second-type boundary condition:
%    - heat flux iQ(1,p) and iQ(n,p) are given
%    - temperature T(1,p) :
%                                   SE
%      T(1,p) = T(2,p-1) - iQ(1,p) ----
%                                   TC
%    - temperature T(n,p) :
%                                       SE
%      T(end,p) = T(n-1,p-1) - iQ(n,p) ----
%                                       TC
%    - temperatures T(m,p) for m = 2, n-1 :
%      T(m,p+1) = M (T(m-1,p) + T(m+1,p)) + T(m,p).(1 - 2.M)
%
%
% INPUTS:
%
% Tin  - initialization temperatures                        [ K           ]
% TC   - thermal conductivity                               [ W.m^-1.K^-1 ]
% DE   - density                                            [ kg.m^-3     ]
% HC   - specific heat capacity                             [ kg.m^-3     ]
% SE   - size of element                                    [ m           ]
% DI   - distance                                           [ m           ]
% TS   - time step                                          [ s           ]
% TT   - total time of simulation                           [ s           ]
% TBC  - type of boundary condition
%        (1 = first-type, 2 = second-type)                  [ -           ]
% BC   - boundary condition
%         - first-type : BC(1) = T(1),  BC(2) = T(n)        [ K           ]
%         - second-type: BC(1) = iQ(1), BC(2) = iQ(n)       [ W.m^-2      ]
%
% OUTPUTS:
%
% Xco  - vector of X coordinate                             [ m           ]
% Tau  - vector of time                                     [ s           ]
% MT   - matrix of temperatures                             [ K           ]
%
% AUXILIARY VARIABLE:
%
% r    - time  step  number                                 [ -           ]
% n    - space point number                                 [ -           ]
% TD   - thermal diffusivity                                [ m^2.m^-1    ]
% M    - modul (Fourier number)                             [ -           ]
%
% Copyright (C) 2013, Technical University of Kosice
% Author  : Zecova Monika, Terpak Jan
% Revision: 27.08.2013
%
  r = round(TT/TS) + 1;
  n = round(DI/SE) + 1;

  Xco = 0.0:SE:DI;
  Tau = 0.0:TS:TT;
  MT  = zeros(r,n);
  MT(1,:) = Tin;

  for j=2:r
    [T]     = explicit_nh (Tin,TC,DE,HC,SE,TS,TBC,BC);
    Tin     = T;
    MT(j,:) = Tin;
  end

end
function [T] = explicit_nh (Tin,TC,DE,HC,SE,TS,TBC,BC)

  n  = length(Tin);
  TD = zeros(1,n-1);
  for i=1:n-1
    TD(i)= TC(i)/(DE(i)*HC(i));
  end
  M  = TD*TS/(SE*SE);

  T  = Tin;

  if (M <= 0.5)

     if TBC == 1
        T(1)   = BC(1);
        T(end) = BC(2);
     end

     if TBC == 2
        T(1)   = T(2)   - BC(1)*SE/TC(1);
        T(end) = T(end-1) - BC(2)*SE/TC(end);
     end

     for m=2:length(Tin)-1
       T(m) = M(m-1)*Tin(m-1) + (1 - (M(m-1)+M(m)))*Tin(m) + M(m)*Tin(m+1);
     end
  else
     disp(['error: unstable solution, M must be <= 0,5; M = ' num2str(M)]);
  end
end


text:
function D = diffmat(N, varargin)
%DIFFMAT   Spectral differentiation matrix.
%   D = DIFFMAT(N) returns the N x N differentiation matrix associated with the
%   Chebyshev spectral collocation method at second-kind Chebyshev points.
%
%   D = DIFFMAT(N, P) returns the N x N differentiation matrix of order P.
%
%   D = DIFFMAT(N, P, DOM) scales the differentiation matrix D to the domain
%   DOM. DOM should be a 1x2 vector.
%
%   D = DIFFMAT(N, P, DOM, GRID) returns the square differentiation matrix on
%   grid specified by GRID. The specifier GRID can be 'chebkind1' (first-kind
%   Chebyshev grid) or 'chebkind2' (second-kind Chebyshev grid) or 'leg'
%   (Legendre grid).
%
%   D = DIFFMAT(N, P, DOM, GRID, LBC, RBC) returns the square differentiation
%   matrix on grid GRID with the corresponding row(s) replaced by boundary
%   conditions specified by string specifiers LBC and RBC for the left and right
%   boundary respectively. The specifier LBC and RBC can be any of 'dirichlet',
%   'neumann', and 'sum' with 'dirichlet' and 'neumann' indicating Dirichlet and
%   Neumann boundary conditions respectively. The string 'sum' indicates a side
%   condition of definite integral i.e., sum(U), where U is the implied solution
%   to the system formed by D. For a differentiation matrix of order higher than
%   1, if there are multiple boundary conditions at one boundary, the
%   corresponding specifiers need to be grouped in a cell using curly brackets.
%   If there is no boundary condition at a boundary, then it needs to be
%   indicated by either an empty array, i.e. [], or an empty cell, i.e. {}. Note
%   that the number of the boundary conditions given by LBC and RBC must total
%   P, i.e. the order of differentiation, otherwise an error is thrown. Also
%   note that RBC can be omitted if there are only left boundary conditions.
%
%   Example 1: D = DIFFMAT(N, 1, DOM, GRID, [], 'dirichlet') replaces the last
%   row of an N x N differentiation matrix by a Dirichlet boundary condition.
%
%   Example 2: D = DIFFMAT(N, 1, DOM, GRID, 'dirichlet') replaces the first
%   row of an N x N differentiation matrix by a Dirichlet boundary condition.
%
%   Example 3: D = DIFFMAT(N, 2, DOM, GRID, 'dirichlet', 'neumann') replaces
%   the first and last rows of an N x N differentiation matrix by a Dirichlet
%   and a Neumann boundary conditions respectively.
%
%   Example 4: D = DIFFMAT(N, 3, DOM, GRID, 'dirichlet', {'dirichlet' 'neumann'})
%   replaces the first row of an N x N differentiation matrix by a Dirichlet
%   boundary condition and the last two rows by Dirichlet and Neumann boundary
%   conditions.
%
%   D = DIFFMAT(N, 'periodic') returns the N x N first-order Fourier
%   differentiation matrix on the default interval [-1 1]. The tag
%   'periodic' can be replaced by 'trig'.
%
%   D = DIFFMAT(N, P, 'periodic') returns the N x N Fourier differentiation
%   matrix of order P  on the default interval [-1 1].
%
%   D = DIFFMAT(N, P, 'periodic', DOM) scales the Pth-order Fourier
%   differetiation matrix to domain DOM.
%
%   The remaining options concern rectangular spectral differentiation matrices,
%   as introduced in [Driscoll & Hale 2014]. Chebfun uses these to solve ODE
%   boundary value problems in the mode GRID1 = 'chebkind2', GRID2 =
%   'chebkind1'. See also [Xu & Hale 2014].
%
%   D = DIFFMAT([M N]) returns the M x N first-order rectangular differentiation
%   matrix which maps from an N-point Chebyshev grid of the second kind to an
%   M-point Chebyshev grid of the same kind.
%
%   D = DIFFMAT([M N], P) returns an M x N rectangular differentiation matrix of
%   order P which maps from an N-point to an M-point Chebyshev grid, both of
%   second kind.
%
%   D = DIFFMAT([M N], P, DOM) returns the same D but scaled to the domain DOM.
%
%   D = DIFFMAT([M N], P, DOM, GRID) returns an M x N first-order rectangular
%   differentiation matrix which maps from an N-point grid of type GRID to an
%   M-point grid of the same type. The specifier GRID can be 'chebkind1',
%   'chebkind2', or 'leg'.
%
%   D = DIFFMAT([M N], P, DOM, GRID1, GRID2) returns an M x N first-order
%   rectangular differentiation matrix which maps from an N-point grid of type
%   GRID1 to an M-point grid of type GRID2. The specifier GRID1 and GRID2 can be
%   any of 'chebkind1', 'chebkind2', and 'leg'.
%
%   D = DIFFMAT(N, P, DOM, GRID1, GRID2, LBC, RBC) returns a square
%   differentiation matrix which is squared up and deflated by appending
%   boundary conditions specified by string specifiers LBC and RBC for the left
%   and the right boundary respectively.
%
%   Example: D = DIFFMAT([M N], P, DOM, GRID1, GRID2, {}, {'neumann' 'sum'})
%   appends two rows to an M x N rectangular differentiation matrix to square it
%   up. The first appended row corresponds to a Neumann boundary condition at
%   the right endpoint while the second appended row corresponds to a side
%   condition sum(U) = I for a scalar I, where U is the solution to the
%   resulting system.
%
% See also DIFF, CHEBCOLLOC2.DIFFMAT, CUMSUMMAT, DIFFROW, INTMAT, INTROW.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References for rectangular differentiation matrices:
%
%   Driscoll, T. and Hale, N., Rectangular spectral collocation, submitted 2014.
%
%   Xu, K. and Hale, N., Explicit construction of rectangular differentiation
%   matrices, submitted 2014.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parse the inputs:
[m, n, p, dom, bc, nlbc, nrbc, mapFrom, mapTo] = parseInputs(N, varargin{:});

%% Different cases:
if ( strcmpi(mapFrom, mapTo) && ( m == n ) ) % Square case:
    if ( strcmpi(mapFrom, 'chebkind1') )
        D = chebcolloc1.diffmat(n, p);
    elseif ( strcmpi(mapFrom, 'chebkind2') )
        D = chebcolloc2.diffmat(n, p);
    elseif ( strcmpi(mapFrom, 'periodic') )
        D = trigcolloc.diffmat(n, p);
    else
        [x, ignored, v] = legpts(n); %#ok<ASGLU>
        [y, ignored, w] = chebpts(n); %#ok<ASGLU>
        z = legpts(n);
        P1 = barymat(y, x, v);
        P2 = barymat(z, y, w);
        D = chebcolloc2.diffmat(n, p);
        D = P2*D*P1;

        % Flipping trick for symmetry:
        d = diag(rot90(D));
        d = sign(d).*(abs(d) + abs(flipud(d)))/2;
        D(logical(flipud(eye(N)))) = d;
        DRot = rot90(D, 2);
        idxTo = rot90(~triu(ones(N)));
        D(idxTo) = (-1)^p*DRot(idxTo);
        if ( mod(N, 2) == 1 )
            D((N+1)/2,(N+1)/2) = 0;
        end
    end

elseif ( strcmpi(mapTo, 'chebkind1') )

    if ( strcmpi(mapFrom, 'chebkind1') )
        if ( p == 1 )
            D = rectdiff1(m, n);
        else
            D = rectdiff_rec(m, n, p, 1);
        end
    elseif ( strcmpi(mapFrom, 'chebkind2') )
        if ( p == 1 )
            D = rectdiff2(m, n);
        else
            D = rectdiff_rec(m, n, p, 2);
        end
    else
        [x, ignored, v] = legpts(n);  %#ok<ASGLU>
        [y, ignored, w] = chebpts(n);  %#ok<ASGLU>
        z = chebpts(m, 1);
        P1 = barymat(y, x, v);
        P2 = barymat(z, y, w);
        D = chebcolloc2.diffmat(n, p);
        D = P2*D*P1;
    end

elseif ( strcmpi(mapTo, 'chebkind2') )
    [z, ignored, ignored, s] = chebpts(m);  %#ok<ASGLU>
    if ( strcmpi(mapFrom, 'chebkind1') )
        D = chebcolloc1.diffmat(n, p);
        [x, ignored, v, r] = chebpts(n, 1);  %#ok<ASGLU>
        P = barymat(z, x, v, s, r, 1);
        D = P*D;
    elseif ( strcmpi(mapFrom, 'chebkind2') )
        D = chebcolloc2.diffmat(n, p);
        [x, ignored, v, r] = chebpts(n);  %#ok<ASGLU>
        P = barymat(z, x, v, s, r, 1);
        D = P*D;
    else
        [x, ignored, v] = legpts(n);  %#ok<ASGLU>
        [y, ignored, w] = chebpts(n);  %#ok<ASGLU>
        P1 = barymat(y, x, v);
        P2 = barymat(z, y, w);
        D = chebcolloc2.diffmat(n, p);
        D = P2*D*P1;
    end

elseif ( strcmpi(mapTo, 'leg') )
    z = legpts(m);
    if ( strcmpi(mapFrom, 'chebkind1') )
        D = chebcolloc1.diffmat(n, p);
        [x, ignored, v] = chebpts(n, 1);  %#ok<ASGLU>
        P = barymat(z, x, v);
        D = P*D;
    elseif ( strcmpi(mapFrom, 'chebkind2') )
        D = chebcolloc2.diffmat(n, p);
        [x, ignored, v] = chebpts(n);  %#ok<ASGLU>
        P = barymat(z, x, v);
        D = P*D;
    else
        [x, ignored, v] = legpts(n);  %#ok<ASGLU>
        [y, ignored, w] = chebpts(n);  %#ok<ASGLU>
        P1 = barymat(y, x, v);
        P2 = barymat(z, y, w);
        D = chebcolloc2.diffmat(n, p);
        D = P2*D*P1;
    end

end

%% Rescaling:
scl = (2/(dom(end) - dom(1)))^p;
D = scl*D;

%% Boundary conditions:
if ( ~isempty(bc) )
    nbc = nlbc + nrbc;
    BC = zeros(nbc, n);
    for j = 1:nbc
        if ( j <= nlbc )
            z = -1;
            r = pi;
            idx = 1;
        else
            z = 1;
            r = 0;
            idx = n;
        end

        switch bc{j}
            case 'dirichlet'
                if ( strcmpi(mapFrom, 'chebkind1') )
                    [x, ignored, v, t] = chebpts(n, 1);  %#ok<ASGLU>
                    BC(j,:) = barymat(z, x, v, r, t);
                elseif ( strcmpi(mapFrom, 'chebkind2') )
                    I = eye(n);
                    BC(j,:) = I(idx,:);
                else
                    [x, ignored, v] = legpts(n);  %#ok<ASGLU>
                    BC(j,:) = barymat(z, x, v);
                end

            case 'neumann'

                if ( strcmpi(mapFrom, 'chebkind1') )
                    DD = diffmat(n, 1, dom, 'chebkind1');
                    [x, ignored, v, t] = chebpts(n, 1);  %#ok<ASGLU>
                    P = barymat(z, x, v, r, t);
                    DD = P*DD;
                    BC(j,:) = DD;
                elseif ( strcmpi(mapFrom, 'chebkind2') )
                    DD = diffmat(n, 1, dom);
                    BC(j,:) = DD(idx,:);
                else
                    DD = diffmat(n, 1, dom);
                    [x, ignored, v] = legpts(n);  %#ok<ASGLU>
                    y = chebpts(n);
                    P = barymat(y, x, v);
                    DD = DD*P;
                    BC(j,:) = DD(idx, :);
                end

            case 'sum'

                if ( strcmpi(mapFrom, 'chebkind1') )
                    [ignored, w] = chebpts(n, dom, 1);  %#ok<ASGLU>
                elseif ( strcmpi(mapFrom, 'chebkind2') )
                    [ignored, w] = chebpts(n, dom);  %#ok<ASGLU>
                else
                    [ignored, w] = legpts(n, dom);  %#ok<ASGLU>
                end

                BC(j,:) = w;

            otherwise
                error('CHEBFUN:diffmat:wrongBC', ...
                    'Unknown type of boundary conditions.');
        end
    end
end

%% Replace or append the boundary conditions:
if ( ~isempty(bc) && ( m == n ) )
    % Replacement for square case:
    D(1:nlbc, :) = BC(1:nlbc, :);
    D(end-nrbc+1:end,:) = BC(nlbc+1:end,:);
elseif ( ~isempty(bc) )
    D = [BC(1:nlbc,:); D; BC(nlbc+1:end,:)];
end

end

function D = rectdiff1(m, n)
%RECTDIFF1  Explicit constrcution of 1st-order rectangular differentiation
%matrix mapping from 1st-kind grid.

% mapping-from grid (angles):
T = chebtech1.angles(n).';

% difference between dimensions:
c = n-m;

% mapping-to grid (angles):
TAU = chebtech1.angles(m);

% Denominator:
denom = bsxfun(@(u,v) 2*sin((v+u)/2).*sin((v-u)/2), T, TAU);

% Sign:
sgn = ones(m, n);
if ( mod(c, 2) )
    sgn(1:2:end,2:2:end) = -1;
    sgn(2:2:end,1:2:end) = -1;
else
    sgn(1:2:end,1:2:end) = -1;
    sgn(2:2:end,2:2:end) = -1;
end

D = sgn.*((cos(c*TAU)./sin(TAU))*sin(T)-sin(c*TAU)/n*sin(T)./denom)./denom;

% indices for applying negative-sum trick:
[~, idx] = min(abs(denom), [], 2);
idx = sub2ind([m n], 1:m, idx.');
D(idx) = 0;
D(idx) = -sum(D, 2);

if ( c == 1 )
    % Flipping trick:
    ii = logical(rot90(tril(ones(m, n)), 2));
    rot90D = rot90(D,2);
    D(ii) = -rot90D(ii);
end

end

function D = rectdiff2(m, n)
%RECTDIFF2  Explicit constrcution of 1st-order rectangular differentiation
%matrix mapping from a 2nd-kind grid.

nm1 = n - 1;                    % For convenience.
cm1 = nm1 - m;                  % Difference between dimensions:
t = chebpts(n).';               % Second-kind grid.
tau = chebpts(m, 1);            % First-kind grid.
T = chebtech2.angles(n).';          % Second-kind grid (angles).
TAU = chebtech1.angles(m);   % First-kind grid (angles).

% Explicit expression:
denom = 2*bsxfun( @(u,v) sin((v+u)/2) .* sin((v-u)/2), T, TAU );
numer = bsxfun( @times, 1 - tau*t, cos(cm1*TAU)./sin(TAU) );

sgn = (-1)^cm1;

if ( cm1 == 0 )
    D = numer ./ denom.^2 / nm1;
else
    D = repmat(sin(cm1*TAU), 1, n)./denom + numer./denom.^2 / nm1;
    D = sgn*D;
end
D(:,[1,n]) = .5*D(:,[1,n]);     % Scaling for first and last columns.

if ( cm1 == 0 )
    % Flipping trick:
    ii = logical(rot90(tril(ones(m, n)), 2));
    rot90D = rot90(D,2);
    D(ii) = sgn*rot90D(ii);
end

% Sign:
D(1:2:end,1:2:end) = -D(1:2:end,1:2:end);
D(2:2:end,2:2:end) = -D(2:2:end,2:2:end);

% Negative sum trick:
[~, idx] = min(abs(denom), [], 2);
idx = sub2ind([m n], 1:m, idx.');
D(idx) = 0; D(idx) = -sum(D, 2);

if ( cm1 == 0 )
    % Fix corner values:
    D(1) = -.25/(nm1*sin(pi/(2*m))*sin(pi/(4*m))^2); D(end) = -D(1);
    % Negative sum trick for corner entries:
    D(1,2) = -sum(D(1,[1 3:end])); D(end,end-1) = -D(1,2);
end

end

function D = rectdiff_rec(m, n, p, kind)
%Recursive construction for high-order rectangular differentiation matrices

% Sign and scaling:
sgn = ones(1, n);
sgn(1:2:end) = -1;

if ( kind == 1 )
    % mapped-from angles (1st-kind):
    T = chebtech1.angles(n).';

    % Compute the first order diff matrix:
    D = rectdiff1(m, n);

    % Preparation for higher order (p>1):
    a = [zeros(n,1); 1];

    % Signs:
    sgn = (-1)^(n-1)*sgn.*sin(T)/n;

else
    % mapped-from grid (2nd-kind):
    T = chebtech2.angles(n).';

    % Compute the first order diff matrix:
    D = rectdiff2(m, n);

    % Preparation for higher order (p>1):
    a = [zeros(n-2,1); -1; 0; 1];

    % Signs:
    sgn = (-1)^(n-1)*sgn/(2*(n-1));
    sgn([1 n]) = sgn([1 n])/2;

end

% mapped-to grid (1st-kind):
tau = chebpts(m, 1);
TAU = chebtech1.angles(m);

a = computeDerCoeffs(a);

% Denominator:
denom = bsxfun(@(u,v) 2*sin((v+u)/2).*sin((v-u)/2), T, TAU);

% indices for applying negative-sum trick:
[~, idx] = min(abs(denom), [], 2);
idx = sub2ind([m n], 1:m, idx.');

% Recursion for higher-order matrices:
for l = 2:p

    % Compute coefficients of the derivative of T_n:
    a = computeDerCoeffs(a);

    % Evaluating at tau by Clenshaw method:
    Tt = chebtech.clenshaw(tau, a);
    D = (Tt*sgn + l*D)./denom;

end

% negative-sum trick:
D(idx) = 0;
D(idx) = -sum(D, 2);

end

function cout = computeDerCoeffs(c)
%COMPUTEDERCOEFFS   Recurrence relation for coefficients of derivative.
%   C is the matrix of Chebyshev coefficients of a (possibly array-valued)
%   CHEBTECH object.  COUT is the matrix of coefficients for a CHEBTECH object
%   whose columns are the derivatives of those of the original.

    [n, m] = size(c);
    cout = zeros(n-1, m);                       % Initialize vector {c_r}
    w = repmat(2*(1:n-1)', 1, m);
    v = w.*c(2:end,:);                          % Temporal vector
    cout(n-1:-2:1,:) = cumsum(v(n-1:-2:1,:));   % Compute c_{n-2}, c_{n-4},...
    cout(n-2:-2:1,:) = cumsum(v(n-2:-2:1,:));   % Compute c_{n-3}, c_{n-5},...
    cout(1,:) = .5*cout(1,:);                   % Adjust the value for c_0
end

function [m, n, p, dom, bc, nlbc, nrbc, mapFrom, mapTo] = parseInputs(N, varargin)
% Parse the inputs to DIFFMAT.

p = 1;
dom = [-1 1];
mapFrom = [];
mapTo = [];
lbc = {};
nlbc = 0;
rbc = {};
nrbc = 0;
isLbcGiven = 0;
isRbcGiven = 0;

if ( isscalar(N) )
    n = N;
    m = n;
elseif ( N(2) < 0 )
    n = N(1);
    m = n + N(2);
else
    m = N(1);
    n = N(2);
end

for j = 1:numel(varargin)
    v = varargin{j};
    if ( isnumeric(v) )
        if ( isempty(v) )
            if ( ~isLbcGiven )
                lbc = {};
                isLbcGiven = 1;
            elseif ( ~isRbcGiven )
                rbc = {};
                isRbcGiven = 1;
            end
        elseif ( isscalar(v) )
            p = v;
        else
            dom = v;
        end
    elseif ( ischar(v) )
        switch v
            case 'rect'

                if ( strcmpi(mapFrom, 'periodic') )
                    error('CHEBFUN:diffmat:wrongInput', ...
                        ['Rectangular Fourier differentiation matrices are '...
                        'not supported.']);
                end

                if ( isscalar(N) )
                    m = n - p;
                end

            case {'periodic', 'trig'}
                mapFrom = 'periodic';
                mapTo = 'periodic';
                if ( m ~= n )
                    error('CHEBFUN:diffmat:wrongInput', ...
                        ['Rectangular Fourier differentiation matrices are '...
                        'not supported.']);
                end

            case {'chebkind1', 'chebkind2', 'leg'}
                if ( isempty(mapFrom) )
                    mapFrom = v;
                elseif ( isempty(mapTo) )
                    mapTo = v;
                else
                    error('CHEBFUN:diffmat:unknown', ...
                        'Too many inputs for grid type.');
                end

            case {'dirichlet', 'neumann', 'sum', []}
                if ( ~isLbcGiven )
                    lbc = {v};
                    isLbcGiven = 1;
                elseif ( ~isRbcGiven )
                    rbc = {v};
                    isRbcGiven = 1;
                else
                    error('CHEBFUN:diffmat:unknown', ...
                        ['Too many inputs for boundary condition. ' ...
                        'Use curly brackets to group left and right ' ...
                        'boundary conditions, if multiple boundary ' ...
                        'conditions are considered at one boundary.']);
                end
            otherwise
                error('CHEBFUN:diffmat:unknown', ['Unknown input ', v]);
        end

    elseif ( iscell(v) )
        if ( ~isLbcGiven )
            lbc = v;
            isLbcGiven = 1;
        elseif ( ~isRbcGiven )
            rbc = v;
            isRbcGiven = 1;
        else
            error('CHEBFUN:diffmat:unknown', ...
                'Unrecognized boundary condition.');
        end
    else
        error('CHEBFUN:diffmat:unknown', ...
            'Unknown input of type %s.', class(v));
    end
end

% If only one grid is given, then let
if ( ~isempty(mapFrom) && isempty(mapTo) )
    mapTo = mapFrom;
end

if ( isempty(mapFrom) )
    mapFrom = 'chebkind2';
    mapTo = 'chebkind2';
end

if ( p < 0 )
    error('CHEBFUN:diffmat:wrongInput', ...
            'The order of differentiation matrix must be non-negative.');
end

% No breakpoints allowed:
if ( numel(dom) > 2 )
    dom = dom([1 end]);
    warning('CHEBFUN:diffmat:noBreaks', ...
        'DIFFMAT does not support domains with breakpoints.');
end

% Boundary conditions:
bc = [lbc rbc];
if ( ~isempty(bc) )
    if ( isa(mapTo, 'periodic')  )
        error('CHEBFUN:diffmat:wrongBC', ...
            ['For periodic functions, there is no need to specify boundary ' ...
             'conditions.']);
    end
    nlbc = numel(lbc);
    nrbc = numel(rbc);
    nbc = nlbc + nrbc;
    if ( nbc ~= p )
        error('CHEBFUN:diffmat:wrongBC', ...
            ['The number of boundary conditions must match differentiation ' ...
             'order p.']);
    end
end

end


text:
function topLevelBodyInfo = getTopLevelCentralBody(celBodyData)
    %getTopLevelCentralBody Summary of this function goes here
    %   Detailed explanation goes here
    try
        topLevelBodyInfo = celBodyData.emptyBodyInfo;
    catch
        topLevelBodyInfo = KSPTOT_BodyInfo.empty(1,0);
    end
    celBodyFields = fields(celBodyData);
    for(i=1:length(celBodyFields))
        bInfo = celBodyData.(celBodyFields{i});

        if(isempty(bInfo.parentBodyInfo))
            topLevelBodyInfo = bInfo;
            break;
        end
    end

    if(isempty(topLevelBodyInfo))
        topLevelBodyInfo = celBodyData.(celBodyFields{1});
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function separatePausedNSx(varargin)

% separatePausedNsx    Obsoleted. Use splitNSxPauses instead.

disp('This function is obsoleted.')
disp('Use splitNSxPauses instead.');

splitNSxPauses;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cc] = cm32cc(cm3)
% Convert volume from cubic centimeters to cubic centimeters.
% Chad Greene 2012
cc = cm3;


text:
function [f,J,D] = tapas_ceode_fx_cmc(x, u, P, M, xtau)
% [f,J,D] = tapas_ceode_fx_cmc(x, u, P, M, xtau)
%
% Computation of the dynamical (state) equation for a neural mass model
% (canonical microcircuit). Compatibile with continuous extension for ODE
% methods. Adapted from spm_fx_cmc.m (original function license below).
% Function output signature kept the same to ensure compatibility.
%
%
% INPUT
%   x           mat             vector of current state activity
%   u           mat             vector of current driving input
%   P           struct          parameter structure
%   M           struct          model specification
%   xtau        mat             vector of delayed state activity
%
% OUTPUT
%   f           mat             matrix of dx/dt
%   J           mat             Jacobian of the system df/dx. Fixed to 0
%                               since irrelevant for this integration
%                               scheme.
%   D           mat             matrix of delays. Empty since irrelevant
%                               for this integration scheme.

% -------------------------------------------------------------------------
%
% Author: Dario Schöbi
% Created: 2020-08-10
% Copyright (C) 2020 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS ceode Toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
%
% ORIGINAL FUNCTION LICENSE
%
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging
%
% Karl Friston
% $Id: spm_fx_cmc.m 6720 2016-02-15 21:06:55Z karl $
%--------------------------------------------------------------------------

% get dimensions and configure state variables
%--------------------------------------------------------------------------
if nargin < 5
    xtau = x;
end

% get dimensions and configure state variables
%--------------------------------------------------------------------------
x  = spm_unvec(x,M.x);            % neuronal states
xtau = spm_unvec(xtau, M.x);
n  = size(x,1);                   % number of sources

% [default] fixed parameters
%--------------------------------------------------------------------------
E  = [1 1/2 1 1/2]*200;           % extrinsic (forward and backward)
G  = [4 4 8 4 4 2 4 4 2 1]*200;   % intrinsic connections
T  = [2 2 16 28];                 % synaptic time constants
R  = 2/3;                         % slope of sigmoid activation function
B = 0;                            % bias or background (sigmoid)

% [specified] fixed parameters
%--------------------------------------------------------------------------
if isfield(M,'pF')
    try, E = M.pF.E; end
    try, G = M.pF.G; end
    try, T = M.pF.T; end
    try, R = M.pF.R; end
end

% Extrinsic connections
%--------------------------------------------------------------------------
% ss = spiny stellate
% sp = superficial pyramidal
% dp = deep pyramidal
% ii = inhibitory interneurons
%--------------------------------------------------------------------------
if n > 1
    A{1} = exp(P.A{1})*E(1);      % forward  connections (sp -> ss)
    A{2} = exp(P.A{2})*E(2);      % forward  connections (sp -> dp)
    A{3} = exp(P.A{3})*E(3);      % backward connections (dp -> sp)
    A{4} = exp(P.A{4})*E(4);      % backward connections (dp -> ii)
else
    A    = {0,0,0,0};
end

% detect and reduce the strength of reciprocal (lateral) connections
%--------------------------------------------------------------------------
for i = 1:length(A)
    L    = (A{i} > exp(-8)) & (A{i}' > exp(-8));
    A{i} = A{i}./(1 + 4*L);
end

% input connections
%--------------------------------------------------------------------------
C    = exp(P.C);

% pre-synaptic inputs: s(V)
%--------------------------------------------------------------------------
R    = R.*exp(P.S);              % gain of activation function
F    = 1./(1 + exp(-R*xtau + B));   % firing rate
S    = F - 1/(1 + exp(B));       % deviation from baseline firing

% input
%==========================================================================
if isfield(M,'u')

    % endogenous input
    %----------------------------------------------------------------------
    U = u(:)*512;

else
    % exogenous input
    %----------------------------------------------------------------------
    U = C*u(:)*32;

end

% time constants and intrinsic connections
%==========================================================================
T    = ones(n,1)*T/1000;
G    = ones(n,1)*G;

% extrinsic connections
%--------------------------------------------------------------------------
% forward  (i)   2  sp -> ss (+ve)
% forward  (ii)  1  sp -> dp (+ve)
% backward (i)   2  dp -> sp (-ve)
% backward (ii)  1  dp -> ii (-ve)
%--------------------------------------------------------------------------
% free parameters on time constants and intrinsic connections
%--------------------------------------------------------------------------
% G(:,1)  ss -> ss (-ve self)  4
% G(:,2)  sp -> ss (-ve rec )  4
% G(:,3)  ii -> ss (-ve rec )  4
% G(:,4)  ii -> ii (-ve self)  4
% G(:,5)  ss -> ii (+ve rec )  4
% G(:,6)  dp -> ii (+ve rec )  2
% G(:,7)  sp -> sp (-ve self)  4
% G(:,8)  ss -> sp (+ve rec )  4
% G(:,9)  ii -> dp (-ve rec )  2
% G(:,10) dp -> dp (-ve self)  1
%--------------------------------------------------------------------------
% Neuronal states (deviations from baseline firing)
%--------------------------------------------------------------------------
%   S(:,1) - voltage     (spiny stellate cells)
%   S(:,2) - conductance (spiny stellate cells)
%   S(:,3) - voltage     (superficial pyramidal cells)
%   S(:,4) - conductance (superficial pyramidal cells)
%   S(:,5) - current     (inhibitory interneurons)
%   S(:,6) - conductance (inhibitory interneurons)
%   S(:,7) - voltage     (deep pyramidal cells)
%   S(:,8) - conductance (deep pyramidal cells)
%--------------------------------------------------------------------------
j     = [1 2 3 4];
for i = 1:size(P.T,2)
    T(:,j(i)) = T(:,j(i)).*exp(P.T(:,i));
end
j     = [7 2 3 4];
for i = 1:size(P.G,2)
    G(:,j(i)) = G(:,j(i)).*exp(P.G(:,i));
end

% Modulatory effects of dp depolarisation on intrinsic connection j(1)
%--------------------------------------------------------------------------
if isfield(P,'M')
    G(:,j(1)) = G(:,j(1)).*exp(-P.M*32*S(:,7));
end

% Motion of states: f(x)
%--------------------------------------------------------------------------

% Conductance
%==========================================================================

% Granular layer (excitatory interneurons): spiny stellate: Hidden causes
%--------------------------------------------------------------------------
u      =   A{1}*S(:,3) + U;
u      = - G(:,1).*S(:,1) - G(:,3).*S(:,5) - G(:,2).*S(:,3) + u;
f(:,2) =  (u - 2*x(:,2) - x(:,1)./T(:,1))./T(:,1);

% Supra-granular layer (superficial pyramidal cells): Hidden causes - error
%--------------------------------------------------------------------------
u      = - A{3}*S(:,7);
u      =   G(:,8).*S(:,1) - G(:,7).*S(:,3) + u;
f(:,4) =  (u - 2*x(:,4) - x(:,3)./T(:,2))./T(:,2);

% Supra-granular layer (inhibitory interneurons): Hidden states - error
%--------------------------------------------------------------------------
u      = - A{4}*S(:,7);
u      =   G(:,5).*S(:,1) + G(:,6).*S(:,7) - G(:,4).*S(:,5) + u;
f(:,6) =  (u - 2*x(:,6) - x(:,5)./T(:,3))./T(:,3);

% Infra-granular layer (deep pyramidal cells): Hidden states
%--------------------------------------------------------------------------
u      =   A{2}*S(:,3);
u      = - G(:,10).*S(:,7) - G(:,9).*S(:,5) + u;
f(:,8) =  (u - 2*x(:,8) - x(:,7)./T(:,4))./T(:,4);

% Voltage
%==========================================================================
f(:,1) = x(:,2);
f(:,3) = x(:,4);
f(:,5) = x(:,6);
f(:,7) = x(:,8);

if nargout < 2; return, end

% Jacobian and delay matrix (will be properly computed in
% tnudcm_compute_xtau)
%==========================================================================
J  = 0;
D  = [];

author: translationalneuromodeling

repo: tapas


text:
function student_cdf_values_test ( )

%*****************************************************************************80
%
%% STUDENT_CDF_VALUES_TEST demonstrates the use of STUDENT_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'STUDENT_CDF_VALUES_TEST:\n' );
  fprintf ( 1, '  STUDENT_CDF_VALUES returns values of \n' );
  fprintf ( 1, '  the Student T Cumulative Density Function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      C     X       STUDENT_CDF(A,X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, c, x, fx ] = student_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %9f  %9f  %24.16f\n', c, x, fx );

  end

  return
end


text:
function [R,subgrad,data] = risk_svm(data,W)
% RISK_SVM evaluates SVM risk term given by L1-hinge loss.
%
% Synopsis:
%  [R,subgrad,data] = risk_svm(data)
%  [R,subgrad,data] = risk_svm(data,W)
%
% Description:
%  Let the risk term be defined as
%
%  R(W) = 1/m sum_{i=1}^m max(0, 1-data.y(i)*(W'*[data.X(:,i);data.X0]))
%
%  This function returns value R and subgradient SUBGRAD of the
%  risk R(W) at W.
%

if nargin < 2
    % evaluates the risk and the subgradient of at W = 0

    nExamples = size(data.X,2);
    R = 1;
    if data.X0 == 0
        subgrad = -data.X*data.y(:)/nExamples;
    else
        subgrad = [-data.X*data.y(:)/nExamples; ...
                   -data.X0*sum(data.y)/nExamples];
    end

else

    nExamples = size(data.X,2);

    if data.X0 == 0
        dfce = 1 - (W'*data.X).*reshape(data.y,1,nExamples);

        idx = find(dfce>0);
        R = sum(dfce(idx))/nExamples;

        subgrad = -data.X(:,idx)*reshape(data.y(idx),length(idx),1)/nExamples;
    else

        dfce = 1 - (W(1:end-1)'*data.X + W(end)*data.X0).*(data.y(:)');
        idx = find(dfce>0);

        R = sum(dfce(idx))/nExamples;
        subgrad = [-data.X(:,idx)*reshape(data.y(idx),length(idx),1)/nExamples; ...
                   -data.X0*sum(data.y(idx))/nExamples];
    end

end

return;
% EOF

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function demoBallisticScenario()
%%DEMOBALLISTICSCENARIO Demonstrate how the Joint Probabilistic Data
%             Association Filter (JPDAF) and the Global Nearest Neighbor
%             (GNN)-JPDAF can be used when tracking  two ballistic targets
%             with the possibility of missed detections and false alarms.
%             Track initiation and termination are not considered. This
%             function demonstrates how nonlinear filtering algorithms can
%             be used with the JPDAF/ GNN-JPDAF. Clistering is not
%             performed.
%
%A simple scenario of a ballistic target going from the Sea of Japan to
%an area in the middle of the Pacific (ignoring atmospheric drag, using a
%simple J2 model for the Earth's gravitational field is considered
%(ignoring light pressure and the gravitational attraction of other
%celestial bodies). Measurements are obtained by a radar in the Aleutian
%islands using a spherical measurement model. Refraction and other
%atmospheric effects are not simulated; Doppler is not used. Tracking and
%measurements are given in an Earth-centered Earth-fixed (ECEF) coordinate,
%system, requiring corrections for the Coriolis effect. Measurement-warping
%effects due to the Earth's rotation are not considered.
%
%The true trajectories for simulation are created using the
%createBallisticTrajectories subroutine in this file. Measurements are
%assumed to be available when the target is above a certain elevation in
%the receiver's coordinate system. The elevation in the receiver's
%coordinate system is measured with respect to the local tangent plane to
%the WGS84 reference ellipsoid. The trajectories go from the Sea of Japan
%to the middle of the Pacific. The function orbVelDet2PtMinEng is used to
%get an approximate initial velocity for a minimum energy trajectory
%between the points. However, the function is only an approximation as it
%is meant for use in an Earth-centered-inertial coordinate system (thus
%trajectories are minimum energy ignoring the energy imparted by the
%rotation of the Earth. Additionally, the endpoints are specified in an
%ECEF coordinate system, so the initial conditions outputted by the
%function are only approximate (a ballistic projectile would miss its
%target due to the Earth's rotation). Additionally, atmospheric drag is
%not taken into account and the function assumes a Keplerian dynamic model.
%Thus, the results might be useful as initial estimates for a
%higher-fidelity trajectory estimation routine, but would be unsuitable for
%use on their own. For the purpose of this simulation, the trajectory is
%sufficiently realistic.
%
%To make two targets sufficiently close that they will contest
%measurements (making a JPDAF relevant), the second target's trajectory was
%set to be that of the first target, but with a 50m offset in range. The
%resulting trajectory is no longer truly ballistic, but is close enough for
%the simulation.
%
%The tracks are started using the spher2CartOrbitCubature function (which
%is for a Keplerian model, not a J2 model). The correct measurements for
%two scans are used to start the tracks. After that, data association
%uncertainty can exist. This method of starting the tracks demonstrates how
%cubature integration can be used with an arbitrary dynamic model and
%measurement function to start a track. Note that additional steps would
%need to be taken if one were tracking satellites and multiple orbits
%between signtings were allowed.
%
%The JPDAF and GNN-JPDAF algorithms require a likelihood matrix to
%make their (soft) associations of measurements to targets. The function
%makeStandardCartOnlyLRMatHyps creates the likelihood matrix based on the
%exponent of the dimensionless score function of [1]. This requires that
%the clutter density in measurement space and the detection probability be
%known (or estimated). A nominal detection probability of 80% is used.
%
%The makeStandardCartOnlyLRMatHyps function also updates the tracks for
%each of the measurements. The sqrtKalmanUpdate function is used with
%Cartesian-converted measurements. The measurement conversion is performed
%using the spher2CartCubature function via fifth-order cubature
%integration.
%
%The target state is predicted using the momentMatchPred function, which
%can handle the nonlinear dynamic model. The momentMatchPred function
%utilizes non-stochastic Runge-Kutta methods to propagate the mean and
%square root covariance matrix.
%
%The trajectories are plotted along with the measurements. Also plotted is
%the offset from the range estimates of the target and the location of the
%first target. One could plot the mean squared error of the trajectories as
%well, but this is less informative regarding how well targets are
%separated as the cross-range errors for a monostatic tracking scenario
%such as this tend to be very high.
%
%REFERENCES:
%[1] Y. Bar-Shalom, S. S. Blackman, and R. J. Fitzgerald, "Dimensionless
%    score function for multiple hypothesis tracking," IEEE Transactions on
%    Aerospace and Electronic Systems, vol. 43, no. 1, pp. 392-400, Jan.
%    2007.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

disp('An example of simple nonlinear filtering and data association for a ballistic model')

%%%SETUP THE MODEL
disp('Setting parameters and creating the simulated trajectories.')
PD=0.8;%Detection probability --same for all targets.
%The number of observations in the simulation.
numObs=100;

%Assumed clutter parameter in terms of false alarms per meter-radian^2
%(spherical coordinates). False alarms are generated in the measurement
%domain of the radar. The tracker can still work if the value for lambda
%does not perfectly match the value simulated/ the real value. However, one
%should note that lambda in the tracker cannot be arbitrarily large lest
%the tracker eventually ignore all measurements (false alarms become always
%more likely than true tracks).
lambda=1e-2;
numMethods=2;%2 filters, a JPDAF and an GNN-JPDAF.

%The AlgSel parameters are inputs to the singleScanUpdate function.
%Method 1 uses the JPDAF
algSel1(1)=1;
algSel2(1)=0;
%Method 2 uses the GNN-JPDAF
algSel1(2)=0;
algSel2(2)=0;

%The location of the simulated radar in the Aleutian islands.
posRx=[51.86310; -177.98669]*(pi/180);
rRx=ellips2Cart([posRx;0]);

%The continuous-time dynamic model.
rDot=@(rState,t)aJ2Gravity(rState);

%Get the local coordinate axes of the receiver.
uENU=getENUAxes(posRx);

%Rotation matrix from global Cartesian coordinates to ENU
RGlob2ENU=findTransParam(eye(3,3),uENU);

%The nonlinear measurement function
h=@(x)Cart2Sphere(x(1:3,:),0,true,rRx,rRx,RGlob2ENU);

%Next, create the ballistic trajectories
[tObs,xObs,zObsTrue]=createBallisticTrajectories(rDot,h,numObs);
numTargets=size(zObsTrue,2);
PD=repmat(PD,[numTargets,1]);
zDim=size(zObsTrue,1);

%The time between samples is constant in this simulation, though missed
%detections are possible.
deltaT=tObs(2)-tObs(1);

sigmaR=10;
sigmaAzel=0.1*(pi/180);
%The square root measurement covariance matrix, local spherical
%coordinates. Here, it is assumed diagonal.
SR=diag([sigmaR;sigmaAzel;sigmaAzel]);

%Region in range around the targets in which to generate false alarms.
clutR=50e3;%50km
%Angular region in which to generate random false alarms around the
%targets.
clutAngs=2*(pi/180);
%The volume used for generating false alarms based times the given false
%alarm density is the expected number of false alarms per scan.
lambdaV=lambda*(2*clutR)*(2*clutAngs)^2;

%For measurement conversion
[xi3,w3]=fifthOrderCubPoints(3);

%For state propagation
[xi6,w6]=fifthOrderCubPoints(6);

%Initializations
q0=processNoiseSuggest('PolyKal-ROT',1e-2,deltaT);

a=rDot;
D=@(x,t)DPoly(1,q0,1,3);

disp('Initializing the tracks.')
%Generate the first two observations explicitly to start the tracks using an
%information filter. All methods get the same initialization.
xStates=zeros(6,numTargets,numObs,numMethods);
SStates=zeros(6,6,numTargets,numObs,numMethods);
zScans=cell(numObs,1);

zScans{1}=[];
zScans{2}=[];
zScansSpher{1}=[];
zScansSpher{2}=[];
for curTrack=1:numTargets
    z1=zObsTrue(:,curTrack,1)+SR*randn(3,1);
    [zCart1,RCart1]=spher2CartCubature(z1,SR,0,true,rRx,rRx,RGlob2ENU,xi3,w3);
    z2=zObsTrue(:,curTrack,2)+SR*randn(3,1);
    [zCart2,RCart2]=spher2CartCubature(z2,SR,0,true,rRx,rRx,RGlob2ENU,xi3,w3);

    zScansSpher{1}=z1;
    zScansSpher{2}=z2;

    [xInit,PInit]=spher2CartOrbitCubature([z1,z2],SR,0,true,deltaT,rRx,rRx,RGlob2ENU,true,xi6,w6);

    SInit=chol(PInit,'lower');

    for curMethod=1:numMethods
        xStates(:,curTrack,2,curMethod)=xInit;
        SStates(:,:,curTrack,2,curMethod)=SInit;
    end
    zScans{1}=[zScans{1},zCart1];
    zScans{2}=[zScans{2},zCart2];
end

disp('Running the tracker simulation (this will take a little while).')
%Run the tracker starting with the third measurement.
for curObs=3:numObs
    %Determine the number of false alarms to generate.
    numFalse=PoissonD.rand(1,lambdaV);

    %Determine which, if any, targets should be detected.
    isDet=rand(numTargets,1)<PD;

    %Allocate space for the detections.
    numMeas=numFalse+sum(isDet);
    zCur=zeros(3,numMeas);
    curDet=1;

    %Generate the detection from the targets, if any.
    for curTar=1:numTargets
        if(isDet(curTar)||curObs<5)
            zCur(:,curDet)=zObsTrue(:,curTar,curObs)+SR*randn(3,1);
            curDet=curDet+1;
        end
    end

    %Generate the false alarm detections, if any. They are centered around
    %the first target.
    rClutBounds=zObsTrue(1,1,curObs)+[-clutR;clutR];
    azBounds=zObsTrue(2,1,curObs)+[-clutAngs;clutAngs];
    elBounds=zObsTrue(3,1,curObs)+[-clutAngs;clutAngs];
    for curFalse=1:numFalse
        r=UniformD.rand(1,rClutBounds);
        az=UniformD.rand(1,azBounds);
        el=UniformD.rand(1,elBounds);

        zCur(:,curDet)=[r;az;el]+SR*randn(3,1);
        curDet=curDet+1;
    end

    %Convert all of the measurements to Cartesian coordinates using
    %cubature integration.
    zCart=zeros(zDim,numMeas);
    SRCart=zeros(zDim,zDim,numMeas);
    measJacobDet=zeros(numMeas,1);
    for curMeas=1:numMeas
        [zCart(:,curMeas),RCart]=spher2CartCubature(zCur(:,curMeas),SR,0,true,rRx,rRx,RGlob2ENU,xi3,w3);
        SRCart(:,:,curMeas)=chol(RCart,'lower');

        %The measurement Jacobian determinants are used for transforming
        %the clutter density.
        measJacobDet(curMeas)=det(calcSpherConvJacob(zCur(:,curMeas),0,true,rRx,rRx,RGlob2ENU));
    end
    zScans{curObs}=zCart;
    zScansSpher{curObs}=zCur;

    tPrev=tObs(curObs-1);
    tCur=tObs(curObs);

    %Update all of the tracks for the methods.
    for curMethod=1:numMethods
        %Predict the tracker to the current time for each target.
        x=xStates(:,:,curObs-1,curMethod);
        S=SStates(:,:,:,curObs-1,curMethod);
        for curTar=1:numTargets
            [x(:,curTar),S(:,:,curTar)]=momentMatchPred(x(:,curTar),S(:,:,curTar),a,D,3,tPrev,tCur,xi6,w6);
        end

        %Brute-force hypothesis formation and likelihood matrix
        %computation. A is the likelihood matrix.
        [A,xHyp,PHyp]=makeStandardCartOnlyLRMatHyps(x,S,zCart,SRCart,[],PD,lambda,[],[],measJacobDet);

        %Perform the single scan assignment
        [xPost,PPost]=singleScanUpdate(xHyp,PHyp,A,algSel1(curMethod),algSel2(curMethod));

        xStates(:,:,curObs,curMethod)=xPost;
        for curTar=1:numTargets
            SStates(:,:,curTar,curObs,curMethod)=chol(PPost(:,:,curTar),'lower');
        end
    end
end

disp('Plotting the measurements and tracks over the Earth.')
%Indicates how the lines for the method are drawn.
methodStyle{1}='r';
methodStyle{2}='b';

%%%Plot the estimated trajectories on the true trajectory
figure(1)
clf
hold on
clf
plotMapOnEllipsoid();
hold on;

%Plot all of the measurements for all of the scans
for curObs=1:numObs
    zCur=zScans{curObs};
    scatter3(zCur(1,:),zCur(2,:),zCur(3,:),'.b')
end

%Plot the tracks
plot3(xObs(1,:),xObs(2,:),xObs(3,:),'-g','linewidth',4)
for curMethod=1:numMethods
    %Track 1 is drawn in red; track 2 is drawn in blue.
    for curTrack=1:numTargets
        if(curTrack==1)
            lineOpts=['-',methodStyle{curMethod}];
        else
            lineOpts=['--',methodStyle{curMethod}];
        end

        x=xStates(:,curTrack,:,curMethod);
        plot3(x(1,2:end),x(2,2:end),x(3,2:end),lineOpts,'LineWidth',2)
    end
end

disp('Plotting the offset in range of the track estimates from the first target.')
disp('The lines for track 1 are solid; those for track 2 are dashed.')
disp('JPDAF estimates are red; GNN estimates are blue.')
disp('The observations are shown in black. The estimates are relatively close to the observations.')
%%%Compute and then plot the offset in range from the primary target.
rErrTargets=zeros(numObs,numTargets,numMethods);
for curObs=2:numObs
    for curMethod=1:2
        for curTrack=1:numTargets
            xState=xStates(:,curTrack,curObs,curMethod);
            SState=SStates(:,:,curTrack,curObs,curMethod);
            r=unbiasedR(xState,SState,RGlob2ENU,rRx,xi6,w6);

            rErrTargets(curObs,curTrack,curMethod)=r-zObsTrue(1,1,curObs);
        end
    end
end

figure(3)
clf
hold on
for curMethod=1:numMethods
    %Track 1 is drawn in red; track 2 is drawn in blue.
    for curTrack=1:numTargets
        if(curTrack==1)
            lineOpts=['-',methodStyle{curMethod}];
        else
            lineOpts=['--',methodStyle{curMethod}];
        end

        r=rErrTargets(:,curTrack,curMethod);
        plot(tObs(2:end),r(2:end),lineOpts,'LineWidth',2)
    end
end

%Plot the range of all of the measurements for all of the scans too
for curObs=2:numObs
    zCur=zScansSpher{curObs};
    numObs=size(zCur,2);
    if(numObs>0)
        scatter(repmat(tObs(curObs),[numObs,1]),bsxfun(@minus,zCur(1,:),zObsTrue(1,1,curObs)),'.k');
    end
end

h1=xlabel('Time (seconds)');
h2=ylabel('Offset in Range (meters)');

set(gca,'FontSize',18,'FontWeight','bold','FontName','Times')
set(h1,'FontSize',20,'FontWeight','bold','FontName','Times')
set(h2,'FontSize',20,'FontWeight','bold','FontName','Times')
axis([tObs(2), tObs(end), -300, 300])

end

function [tObs,xObs,zObsTrue]=createBallisticTrajectories(rDot,h,numObs)
    posStart=([39.61693; 134.54975])*pi/180;%The Sea of Japan
    %The middle of the Pacific.
    posEnd=([40.56751; -159.53953])*pi/180;
    rStart=ellips2Cart([posStart;0]);
    rEnd=ellips2Cart([posEnd;0]);

    minEL=10*(pi/180);%Minimum elevation for tracking.

    %Get starting velocity
    [vStart,tMinEllip]=orbVelDet2PtMinEng(rStart,rEnd);
    tSpan=[0; tMinEllip*1.25];

    %The function RKAdaptiveOverRange could be used to integrate over the
    %trajectory to a higher polynomial accuracy, but it is simpler to use
    %ode45, because the interpolation routine deval for the solution is
    %easy to use to get observations at desired times.
    options=odeset('RelTol',1e-8,'AbsTol',1e-10,'Refine',20);
    [tFull,xFull]=ode45(@(t,rState)rDot(rState,t),tSpan,[rStart;vStart],options);
    xFull=xFull';

    %Noise-free observations in the local coordinate system of the receiver.
    zTrue=h(xFull(1:3,:));
    elIdx=find(zTrue(3,:)>minEL);

    %Allocate space depending on the chosen scenario.
    zObsTrue=zeros(3,2,numObs);%The main target and one other return.

    %Get the primary trajectory at desired times
    tObs=linspace(tFull(elIdx(1)),tFull(elIdx(end)),numObs);
    options=odeset('RelTol',1e-8,'AbsTol',1e-10);
    sol=ode45(@(t,rState) rDot(rState,t),tSpan,[rStart;vStart],options);
    xObs=deval(sol,tObs);
    zObsTrue(:,1,:)=h(xObs(1:3,:));

    %The second target is offset by a constant amount in range.
    zObsTrue(2:end,2,:)=zObsTrue(2:end,1,:);
    zObsTrue(1,2,:)=zObsTrue(1,1,:)+50;
end

function r=unbiasedR(xState,SState,RGlob2ENU,rRx,xi,w)
    %Get an unbiased conversion of a Cartesian state into a range to a
    %sensor using cubature conversion (ignoring all propagation effects).

    %Transform the cubature points.
    xi=bsxfun(@plus,SState*xi,xState);

    %Subtract the receiver location and rotate into ENU coordinates.
    xiLocal=RGlob2ENU*bsxfun(@minus,xi(1:3,:),rRx);

    %Compute ranges
    rCub=sqrt(sum(xiLocal(1:3,:).^2,1));
    r=sum(bsxfun(@times,rCub(:),w(:)));
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