text:
function [ x, y, z, w ] = ld5294 ( )

%*****************************************************************************80;
%
%% LD5294 computes the 5294 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(5294,1);
  y = zeros(5294,1);
  z = zeros(5294,1);
  w = zeros(5294,1);
  a = 0.0;
  b = 0.0;
  v = 0.9080510764308163E-04;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.2084824361987793E-03;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.2303261686261450E-01;
  v = 0.5011105657239616E-04;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3757208620162394E-01;
  v = 0.5942520409683854E-04;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5821912033821852E-01;
  v = 0.9564394826109721E-04;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.8403127529194872E-01;
  v = 0.1185530657126338E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1122927798060578;
  v = 0.1364510114230331E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1420125319192987;
  v = 0.1505828825605415E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1726396437341978;
  v = 0.1619298749867023E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2038170058115696;
  v = 0.1712450504267789E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2352849892876508;
  v = 0.1789891098164999E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2668363354312461;
  v = 0.1854474955629795E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2982941279900452;
  v = 0.1908148636673661E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3295002922087076;
  v = 0.1952377405281833E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3603094918363593;
  v = 0.1988349254282232E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3905857895173920;
  v = 0.2017079807160050E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4202005758160837;
  v = 0.2039473082709094E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4490310061597227;
  v = 0.2056360279288953E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4769586160311491;
  v = 0.2068525823066865E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5038679887049750;
  v = 0.2076724877534488E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5296454286519961;
  v = 0.2081694278237885E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5541776207164850;
  v = 0.2084157631219326E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5990467321921213;
  v = 0.2084381531128593E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6191467096294587;
  v = 0.2083476277129307E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6375251212901849;
  v = 0.2082686194459732E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6540514381131168;
  v = 0.2082475686112415E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6685899064391510;
  v = 0.2083139860289915E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6810013009681648;
  v = 0.2084745561831237E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6911469578730340;
  v = 0.2087091313375890E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6988956915141736;
  v = 0.2089718413297697E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.7041335794868720;
  v = 0.2092003303479793E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.7067754398018567;
  v = 0.2093336148263241E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3840368707853623E-01;
  v = 0.7591708117365267E-04;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.9835485954117399E-01;
  v = 0.1083383968169186E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.1665774947612998;
  v = 0.1403019395292510E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.2405702335362910;
  v = 0.1615970179286436E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.3165270770189046;
  v = 0.1771144187504911E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.3927386145645443;
  v = 0.1887760022988168E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.4678825918374656;
  v = 0.1973474670768214E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.5408022024266935;
  v = 0.2033787661234659E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.6104967445752438;
  v = 0.2072343626517331E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.6760910702685738;
  v = 0.2091177834226918E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.6655644120217392E-01;
  b = 0.1936508874588424E-01;
  v = 0.9316684484675566E-04;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.9446246161270182E-01;
  b = 0.4252442002115869E-01;
  v = 0.1116193688682976E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1242651925452509;
  b = 0.6806529315354374E-01;
  v = 0.1298623551559414E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1553438064846751;
  b = 0.9560957491205369E-01;
  v = 0.1450236832456426E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1871137110542670;
  b = 0.1245931657452888;
  v = 0.1572719958149914E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2192612628836257;
  b = 0.1545385828778978;
  v = 0.1673234785867195E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2515682807206955;
  b = 0.1851004249723368;
  v = 0.1756860118725188E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2838535866287290;
  b = 0.2160182608272384;
  v = 0.1826776290439367E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3159578817528521;
  b = 0.2470799012277111;
  v = 0.1885116347992865E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3477370882791392;
  b = 0.2781014208986402;
  v = 0.1933457860170574E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3790576960890540;
  b = 0.3089172523515731;
  v = 0.1973060671902064E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4097938317810200;
  b = 0.3393750055472244;
  v = 0.2004987099616311E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4398256572859637;
  b = 0.3693322470987730;
  v = 0.2030170909281499E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4690384114718480;
  b = 0.3986541005609877;
  v = 0.2049461460119080E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4973216048301053;
  b = 0.4272112491408562;
  v = 0.2063653565200186E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5245681526132446;
  b = 0.4548781735309936;
  v = 0.2073507927381027E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5506733911803888;
  b = 0.4815315355023251;
  v = 0.2079764593256122E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5755339829522475;
  b = 0.5070486445801855;
  v = 0.2083150534968778E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1305472386056362;
  b = 0.2284970375722366E-01;
  v = 0.1262715121590664E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1637327908216477;
  b = 0.4812254338288384E-01;
  v = 0.1414386128545972E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1972734634149637;
  b = 0.7531734457511935E-01;
  v = 0.1538740401313898E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2308694653110130;
  b = 0.1039043639882017;
  v = 0.1642434942331432E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2643899218338160;
  b = 0.1334526587117626;
  v = 0.1729790609237496E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2977171599622171;
  b = 0.1636414868936382;
  v = 0.1803505190260828E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3307293903032310;
  b = 0.1942195406166568;
  v = 0.1865475350079657E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3633069198219073;
  b = 0.2249752879943753;
  v = 0.1917182669679069E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3953346955922727;
  b = 0.2557218821820032;
  v = 0.1959851709034382E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4267018394184914;
  b = 0.2862897925213193;
  v = 0.1994529548117882E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4573009622571704;
  b = 0.3165224536636518;
  v = 0.2022138911146548E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4870279559856109;
  b = 0.3462730221636496;
  v = 0.2043518024208592E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5157819581450322;
  b = 0.3754016870282835;
  v = 0.2059450313018110E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5434651666465393;
  b = 0.4037733784993613;
  v = 0.2070685715318472E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5699823887764627;
  b = 0.4312557784139123;
  v = 0.2077955310694373E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5952403350947741;
  b = 0.4577175367122110;
  v = 0.2081980387824712E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2025152599210369;
  b = 0.2520253617719557E-01;
  v = 0.1521318610377956E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2381066653274425;
  b = 0.5223254506119000E-01;
  v = 0.1622772720185755E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2732823383651612;
  b = 0.8060669688588620E-01;
  v = 0.1710498139420709E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3080137692611118;
  b = 0.1099335754081255;
  v = 0.1785911149448736E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3422405614587601;
  b = 0.1399120955959857;
  v = 0.1850125313687736E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3758808773890420;
  b = 0.1702977801651705;
  v = 0.1904229703933298E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4088458383438932;
  b = 0.2008799256601680;
  v = 0.1949259956121987E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4410450550841152;
  b = 0.2314703052180836;
  v = 0.1986161545363960E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4723879420561312;
  b = 0.2618972111375892;
  v = 0.2015790585641370E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5027843561874343;
  b = 0.2920013195600270;
  v = 0.2038934198707418E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5321453674452458;
  b = 0.3216322555190551;
  v = 0.2056334060538251E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5603839113834030;
  b = 0.3506456615934198;
  v = 0.2068705959462289E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5874150706875146;
  b = 0.3789007181306267;
  v = 0.2076753906106002E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6131559381660038;
  b = 0.4062580170572782;
  v = 0.2081179391734803E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2778497016394506;
  b = 0.2696271276876226E-01;
  v = 0.1700345216228943E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3143733562261912;
  b = 0.5523469316960465E-01;
  v = 0.1774906779990410E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3501485810261827;
  b = 0.8445193201626464E-01;
  v = 0.1839659377002642E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3851430322303653;
  b = 0.1143263119336083;
  v = 0.1894987462975169E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4193013979470415;
  b = 0.1446177898344475;
  v = 0.1941548809452595E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4525585960458567;
  b = 0.1751165438438091;
  v = 0.1980078427252384E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4848447779622947;
  b = 0.2056338306745660;
  v = 0.2011296284744488E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5160871208276894;
  b = 0.2359965487229226;
  v = 0.2035888456966776E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5462112185696926;
  b = 0.2660430223139146;
  v = 0.2054516325352142E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5751425068101757;
  b = 0.2956193664498032;
  v = 0.2067831033092635E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6028073872853596;
  b = 0.3245763905312779;
  v = 0.2076485320284876E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6291338275278409;
  b = 0.3527670026206972;
  v = 0.2081141439525255E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3541797528439391;
  b = 0.2823853479435550E-01;
  v = 0.1834383015469222E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3908234972074657;
  b = 0.5741296374713106E-01;
  v = 0.1889540591777677E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4264408450107590;
  b = 0.8724646633650199E-01;
  v = 0.1936677023597375E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4609949666553286;
  b = 0.1175034422915616;
  v = 0.1976176495066504E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4944389496536006;
  b = 0.1479755652628428;
  v = 0.2008536004560983E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5267194884346086;
  b = 0.1784740659484352;
  v = 0.2034280351712291E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5577787810220990;
  b = 0.2088245700431244;
  v = 0.2053944466027758E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5875563763536670;
  b = 0.2388628136570763;
  v = 0.2068077642882360E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6159910016391269;
  b = 0.2684308928769185;
  v = 0.2077250949661599E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6430219602956268;
  b = 0.2973740761960252;
  v = 0.2082062440705320E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4300647036213646;
  b = 0.2916399920493977E-01;
  v = 0.1934374486546626E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4661486308935531;
  b = 0.5898803024755659E-01;
  v = 0.1974107010484300E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5009658555287261;
  b = 0.8924162698525409E-01;
  v = 0.2007129290388658E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5344824270447704;
  b = 0.1197185199637321;
  v = 0.2033736947471293E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5666575997416371;
  b = 0.1502300756161382;
  v = 0.2054287125902493E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5974457471404752;
  b = 0.1806004191913564;
  v = 0.2069184936818894E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6267984444116886;
  b = 0.2106621764786252;
  v = 0.2078883689808782E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6546664713575417;
  b = 0.2402526932671914;
  v = 0.2083886366116359E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5042711004437253;
  b = 0.2982529203607657E-01;
  v = 0.2006593275470817E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5392127456774380;
  b = 0.6008728062339922E-01;
  v = 0.2033728426135397E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5726819437668618;
  b = 0.9058227674571398E-01;
  v = 0.2055008781377608E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6046469254207278;
  b = 0.1211219235803400;
  v = 0.2070651783518502E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6350716157434952;
  b = 0.1515286404791580;
  v = 0.2080953335094320E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6639177679185454;
  b = 0.1816314681255552;
  v = 0.2086284998988521E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5757276040972253;
  b = 0.3026991752575440E-01;
  v = 0.2055549387644668E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6090265823139755;
  b = 0.6078402297870770E-01;
  v = 0.2071871850267654E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6406735344387661;
  b = 0.9135459984176636E-01;
  v = 0.2082856600431965E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6706397927793709;
  b = 0.1218024155966590;
  v = 0.2088705858819358E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6435019674426665;
  b = 0.3052608357660639E-01;
  v = 0.2083995867536322E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6747218676375681;
  b = 0.6112185773983089E-01;
  v = 0.2090509712889637E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function significant = spectsignificance(tests,alpha)
%
% From a estiamtion of p-values through the specttest function,
% It returns 0 or 1 for each power of coherence value, if the p-value is
% below the level of significance 'alpha'.
%
% INPUTS:
%
% tests                 The output of specttest, i.e. a cell
%                       where each element corresponds to one subject
% alpha                 Significance level
%
% OUTPUT:
%
% significant           Struct with 4 fields: lower, higher, lower_corr, higher_corr.
%                       ('_corr' reflects to multiple comparisons correction)
%                       Each field has a field 'states', just like each
%                       element of tests, with fields state(k).psd and
%                       state(k).coh containing 1 if the value is
%                       significant and 0 otherwise
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2017)

if nargin < 2, alpha = 0.01; end

significant = struct();
significant.higher = struct();
significant.higher.state = struct();
significant.higher_corr = struct();
significant.higher_corr.state = struct();
significant.lower = struct();
significant.lower.state = struct();
significant.lower_corr = struct();
significant.lower_corr.state = struct();

K = length(tests.higher.state);
for k = 1:K
    significant.higher.state(k).psd = tests.higher.state(k).psd < alpha;
    significant.higher_corr.state(k).psd = tests.higher_corr.state(k).psd < alpha;
    significant.lower.state(k).psd = tests.lower.state(k).psd < alpha;
    significant.lower_corr.state(k).psd = tests.lower_corr.state(k).psd < alpha;
    significant.higher.state(k).coh = tests.higher.state(k).coh < alpha;
    significant.higher_corr.state(k).coh = tests.higher_corr.state(k).coh < alpha;
    significant.lower.state(k).coh = tests.lower.state(k).coh < alpha;
    significant.lower_corr.state(k).coh = tests.lower_corr.state(k).coh < alpha;
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function pb = plot_scatter_patch(X,r,col,a)
%PLOT_SCATTER_PATCH Summary of this function goes here
%   Detailed explanation goes here

t= 0:pi/100:2*pi;

pb = zeros(size(X,1),1);

for i=1:size(X,1)
    pb(i)=patch( (r(i)*sin(t)+ X(i,1)) , (r(i)*cos(t)+X(i,2)) ,col,'MarkerEdgeColor','none');
    alpha(pb(i),a);
end

end

author: epfl-lasa


text:
%SYMEXPR2SLBLOCK Create symbolic embedded MATLAB Function block
%
% symexpr2slblock(VARARGIN) creates an Embedded MATLAB Function block
% from a symbolic expression. The input arguments are just as used with
% the functions emlBlock or matlabFunctionBlock.
%
% Notes::
% - In Symbolic Toolbox versions prior to V5.7 (2011b) the function to
%   create Embedded Matlab Function blocks from symbolic expressions is
%   'emlBlock'.
% - Since V5.7 (2011b) there is another function named
%   'matlabFunctionBlock' which replaces the old function.
% - symexpr2slblock is a wrapper around both functions, which
%   checks for the installed Symbolic Toolbox version and calls the
%   required function accordingly.
%
% Authors::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also emlBlock, matlabFunctionBlock.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function [] = symexpr2slblock(varargin)

% V5.8 (R2012a)
if verLessThan('symbolic','5.7')
    emlBlock(varargin{:});
elseif verLessThan('symbolic','5.11')
    matlabFunctionBlock(varargin{:});
else
    % Work around a bug in matlabFunctionBlock.m

    %% Read orignal file
	fid = fopen(which('matlabFunctionBlock.m'),'r');
    funStr = fscanf(fid, '%c',inf);
    fclose(fid);

    %% Create Temporary Workaround File
    tmpFName = fullfile(pwd,'tmp_workaround_matlabFunctionBlock.m');

    % perform necessary modifications
    funStr = strrep(funStr,...
        'if isa(b,''Stateflow.EMChart'')',...   Buggy expression to be replaced
        'if ~strcmp(class(b),''Stateflow.EMChart'')'); % Working expression from previous version.
    funStr = strrep(funStr,...
        'matlabFunctionBlock',...   Buggy expression to be replaced
        'tmp_workaround_matlabFunctionBlock'); % Working expression from previous version.

    % write modified function to temporary file
    fid = fopen(tmpFName, 'w');
    fprintf(fid,'%s\n','% This is a temporary file used as workaround for a bug in matlabFunctionBlock.m (R2013b).');
    fprintf(fid,'%s\n','% It should have been deleted automatically by the Robotics Toolbox. If not, feel free to do it manually.');
    fprintf(fid,'%s',funStr);
    fclose(fid);

    % update function database
    rehash

    % perform the actual block generation
    tmp_workaround_matlabFunctionBlock(varargin{:});

    % clean up
    delete(tmpFName);
end

author: petercorke

repo: robotics-toolbox-matlab


text:
% Procedure to get nodal stresses using Extrapolation method
function [sigma] = Extrapolation(stressGP)

%--------------------------------------------------------------------------
% Purpose : Extrapolating stress at Gaussian points to get nodal stress
%           values
%
% Synopsis : sigma = extrapolation(stressGP)
%
% Variable Description :
%       sigma - Stress at nodal points
%       stressGP - Stressat gaussian points
%--------------------------------------------------------------------------
% Extrapolation matrix obtained from Shape Functions
% Stress points are at (-1,-1),(-1,1),(1,1),1,-1)
explmt = [1+sqrt(3)/2    -1/2          1-sqrt(3)/2    -1/2;
           -1/2       1-sqrt(3)/2          -1/2    1+sqrt(3)/2 ;
          1-sqrt(3)/2    -1/2          1+sqrt(3)/2    -1/2;
           -1/2       1+sqrt(3)/2         -1/2     1-sqrt(3)/2 ] ;

sigmax = explmt*stressGP(:,1) ;
sigmay = explmt*stressGP(:,2) ;
sigmaxy = explmt*stressGP(:,3) ;

sigma = [sigmax sigmay sigmaxy];

author: Sable

repo: mcbench-benchmarks


text:
function out = path_manager_dubins(in,p_drone,start_of_simulation)

% PATH_MANAGER_DUBINS - follow Dubins paths between waypoint configurations
%
% Inputs:
%   nb_wpts      - number of waypoint configurations
%   wpt_list     - an array of dimension 5 by p_drone.size_wpt_array.
%                - the first nb_wpts rows define waypoint
%                  configurations
%                - format for each waypoint configuration:
%                  [wn, we, wd, chi_d, Va_d]
%                  where the (wn, we, wd) is the NED position of the
%                  waypoint, chi_d is the desired course at the waypoint,
%                  and Va_d is the desired airspeed along the path.
%
% Outputs:
%   flag   - if flag==1, follow waypoint path
%            if flag==2, follow orbit
%   Va_d   - desired airspeed
%   r      - inertial position of start of waypoint path
%   q      - unit vector that defines inertial direction of waypoint path
%   c      - center of orbit
%   rho    - radius of orbit
%   lambda - direction of orbit (+1 for CW, -1 for CCW)
%

  NN = 0;
  nb_wpts = in(1+NN);
  wpt_list = reshape(in(2+NN:5*p_drone.size_waypoint_array+1+NN),5,p_drone.size_waypoint_array);
  NN = NN + 1 + 5*p_drone.size_waypoint_array;
  pn        = in(1+NN);
  pe        = in(2+NN);
  h         = in(3+NN);
  % Va      = in(4+NN);
  % alpha   = in(5+NN);
  % beta    = in(6+NN);
  % phi     = in(7+NN);
  % theta   = in(8+NN);
  chi     = in(9+NN);
  % p       = in(10+NN);
  % q       = in(11+NN);
  % r       = in(12+NN);
  % Vg      = in(13+NN);
  % wn      = in(14+NN);
  % we      = in(15+NN);
  % psi     = in(16+NN);
  state     = in(1+NN:16+NN);
  NN = NN + 16;
  t         = in(1+NN);

  p = [pn; pe; -h];
  R = p_drone.radius;                 % radius of circular orbits

  persistent wpt_list_prev      % stored copy of old waypoints
  persistent idx_a              % waypoint pointer
  persistent flag_transition    % state of transition state machine
  persistent dubinspath
  persistent flag_need_new_wpts % flag that request new waypoints from path planner
  persistent flag_first_time_in_state

  if start_of_simulation
      wpt_list_prev = zeros(5,p_drone.size_wpt_array);
      flag_need_new_wpts = 0;
      flag_transition = 0;
      flag_first_time_in_state = 1;
  end

  % If waypoints have changed, update the waypoint pointer and plan new
  % Dubin's path
  if min(min(wpt_list==wpt_list_prev))==0
      wpt_list_prev = wpt_list;
      if start_of_simulation
        flag_transition = 0;
      end
      idx_a = 1;
      idx_b = 2;
      start_node = [wpt_list(1:4,idx_a)', 0, 0];
      end_node   = [wpt_list(1:4,idx_b)', 0, 0];
      dubinspath = compute_dubins_param(start_node, end_node, R);
      flag_need_new_wpts = 0;
      flag_first_time_in_state = 1;
  end

  % Define transition state machine
  switch flag_transition
      case 0             % beginning of simulation
          flag   = 1;
          Va_d   = p_drone.Va0;
          r      = dubinspath.ps;
          q      = dubinspath.q1;
          c      = dubinspath.cs;       % not used for waypoint path
          rho    = dubinspath.R;        % not used for waypoint path
          lambda = dubinspath.lams;     % not used for waypoint path
          if flag_first_time_in_state
              flag_first_time_in_state =0;
          end
          flag_transition = 1;

      case 1	% follow first orbit on Dubins path until intersect H1
          flag   = 2;
          Va_d   = p_drone.Va0;
          r      = dubinspath.w1;
          q      = dubinspath.q1;
          c      = dubinspath.cs;       % not used for waypoint path
          rho    = dubinspath.R;        % not used for waypoint path
          lambda = dubinspath.lams;     % not used for waypoint path

          if ((p-dubinspath.w1)'*dubinspath.q1 >= 0)&&(flag_first_time_in_state==1) % start in H1
              flag_transition = 2;
              flag_first_time_in_state = 1;
           elseif (p-dubinspath.w1)'*dubinspath.q1 >= 0 % entering H1
              flag_transition = 3;
              flag_first_time_in_state = 1;
          else
              flag_first_time_in_state = 0;
          end

      case 2    % follow first orbit on Dubins path until on right side of H1
          flag   = 2;
          Va_d   = p_drone.Va0;
          r      = dubinspath.w1;
          q      = dubinspath.q1;
          c      = dubinspath.cs;       % not used for waypoint path
          rho    = dubinspath.R;        % not used for waypoint path
          lambda = dubinspath.lams;     % not used for waypoint path

          if (p-dubinspath.w1)'*dubinspath.q1 < 0 % get to right side H1
              flag_transition = 1;
              flag_first_time_in_state = 1;
          else
              flag_first_time_in_state = 0;
          end

      case 3    % follow straight line on Dubins path until intersect H2
          flag   = 1;
          Va_d   = p_drone.Va0;
          r      = dubinspath.w1;
          q      = dubinspath.q1;
          c      = dubinspath.cs;       % not used for waypoint path
          rho    = dubinspath.R;        % not used for waypoint path
          lambda = dubinspath.lams;     % not used for waypoint path
          flag_first_time_in_state = 0;

          if (p-dubinspath.w2)'*dubinspath.q1 >= 0 % entering H2
              flag_transition = 4;
              flag_first_time_in_state = 1;
          end

      case 4    % follow second orbit on Dubins path until intersect H3
          flag   = 2;
          Va_d   = p_drone.Va0;
          r      = dubinspath.w1;
          q      = dubinspath.q1;
          c      = dubinspath.ce;       % not used for waypoint path
          rho    = dubinspath.R;        % not used for waypoint path
          lambda = dubinspath.lame;     % not used for waypoint path
          flag_first_time_in_state = 0;

          if ((p-dubinspath.w3)'*dubinspath.q3 >= 0)&&(flag_first_time_in_state==1) % start in H3
              flag_transition = 5;
              flag_first_time_in_state=1;
          elseif (p-dubinspath.w3)'*dubinspath.q3 >= 0 % entering H3
              % increase the waypoint pointer
              if idx_a==nb_wpts-1
                  flag_need_new_wpts = 1;
                  idx_b = idx_a+1;
              else
                  idx_a = idx_a+1;
                  idx_b = idx_a+1;
                  flag_transition = 1;
                  flag_first_time_in_state = 1;
              end
              % Plan new Dubin's path to next waypoint configuration
              start_node = [wpt_list(1:4,idx_a)', 0, 0];
              end_node   = [wpt_list(1:4,idx_b)', 0, 0];
              dubinspath = compute_dubins_param(start_node, end_node, R);
          else
              flag_first_time_in_state = 0;
          end

      case 5    % follow first orbit on Dubins path until on right side of H3
          flag   = 2;
          Va_d   = p_drone.Va0;
          r      = dubinspath.w1;
          q      = dubinspath.q1;
          c      = dubinspath.cs;       % not used for waypoint path
          rho    = dubinspath.R;        % not used for waypoint path
          lambda = dubinspath.lams;     % not used for waypoint path
          flag_first_time_in_state = 0;

          if (p-dubinspath.w1)'*dubinspath.q1 < 0 % get to right side of H3
              flag_transition = 4;
              flag_first_time_in_state = 1;
          else
              flag_first_time_in_state = 0;
          end

  end

  out = [flag; Va_d; r; q; c; rho; lambda; state; flag_need_new_wpts];

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function [mu,vecs,vals] = onlineImagePCA_radon(files,batchSize,scale,pixels,thetas,numPerFile)
%onlineImagePCA_radon finds postural eigenmodes based upon a set of
%aligned images (called by findPosturalEigenmodes.m).
%
%   Input variables:
%
%       files -> cell array of VideoReader objects
%       batchSize -> # of images to process at once
%       scale -> image scaling factor
%       pixels -> radon-transform space pixels to use (Lx1 or 1xL array)
%       thetas -> angles used in Radon transform
%       numPerFile -> # of images to use per file
%
%
%   Output variables:
%
%       mu -> mean value for each of the pixels
%       vecs -> postural eignmodes (LxL array).  Each column (vecs(:,i)) is
%                   an eigenmode corresponding to the eigenvalue vals(i)
%       vals -> eigenvalues of the covariance matrix
%
%
% (C) Gordon J. Berman, 2014
%     Princeton University

    if nargin < 3 || isempty(scale);
        scale = 10/7;
    end

    if nargin < 6 || isempty(numPerFile)
        numPerFile = -1;
    end

    Nf = length(files);
    lengths = zeros(Nf,1);
    for i=1:Nf
        lengths(i) = files{i}.NumberOfFrames;
    end

    L = length(pixels);

    testImage = read(files{1},1);
    testImage = testImage(:,:,1);
    s = size(testImage);
    nX = round(s(1)/scale);
    nY = round(s(2)/scale);
    s = [nX nY];

    firstBatch = true;
    tempMu = zeros(1,L);
    totalImages = 0;
    for t=1:Nf

        fprintf(1,'Processing File #%5i out of %5i\n',t,Nf);

        M = lengths(t);
        if numPerFile == -1
            currentNumPerFile = M;
        else
            currentNumPerFile = numPerFile;
        end

        if M < currentNumPerFile
            currentIdx = 1:M;
        else
            currentIdx = randperm(M,currentNumPerFile);
        end
        M = min([lengths(t) currentNumPerFile]);

        if M < batchSize
            currentBatchSize = M;
        else
            currentBatchSize = batchSize;
        end
        num = ceil(M/currentBatchSize);

        currentVideoReader = files{t};

        currentImage = 0;
        X = zeros(currentBatchSize,L);
        for j=1:num

            fprintf(1,'\t Batch #%5i out of %5i\n',j,num);

            if firstBatch

                firstBatch = false;

                parfor i=1:currentBatchSize

                    a = read(currentVideoReader,currentIdx(i));
                    a = double(imresize(a(:,:,1),s));
                    lowVal = min(a(a>0));
                    highVal = max(a(a>0));
                    a = (a - lowVal) / (highVal - lowVal);

                    R = radon(a,thetas);
                    X(i,:) = R(pixels);

                end
                currentImage = currentBatchSize;

                mu = sum(X);
                C = cov(X).*currentBatchSize + (mu'*mu)./ currentBatchSize;

            else

                if j == num
                    maxJ = M - currentImage;
                else
                    maxJ = currentBatchSize;
                end

                tempMu(:) = 0;
                iterationIdx = currentIdx((1:maxJ) + currentImage);
                parfor i=1:maxJ

                    a = read(currentVideoReader,iterationIdx(i));
                    a = double(imresize(a(:,:,1),s));

                    lowVal = min(a(a>0));
                    highVal = max(a(a>0));
                    a = (a - lowVal) / (highVal - lowVal);

                    R = radon(a,thetas);
                    y = R(pixels);
                    X(i,:) = y';
                    tempMu = tempMu + y';

                end

                mu = mu + tempMu;
                C = C + cov(X(1:maxJ,:)).*maxJ + (tempMu'*tempMu)./maxJ;
                currentImage = currentImage + maxJ;

            end

        end

        totalImages = totalImages + currentImage;

    end

    mu = mu ./ totalImages;
    C = C ./ totalImages - mu'*mu;

    fprintf(1,'Finding Principal Components\n');
    [vecs,vals] = eig(C);

    vals = flipud(diag(vals));
    vecs = fliplr(vecs);

author: gordonberman

repo: MotionMapper


text:
function varargout = gui_layout(varargin)
% GUI_LAYOUT: Automated management of the visualization figures (size and position).
%
% USAGE: gui_layout('Update')                   : Apply current layout
%        gui_layout('GetDecorationSize',jFrame) : Get figure decorations
%        gui_layout('GetScreenClientArea')      : Get the screen areas that are used by Brainstorm
%        gui_layout('GetScreenBrainstormAreas') : Get figure area(s)
%        gui_layout('GetFigureGroups')          : Get all the registered Brainstorm figures, grouped by dataset
%        gui_layout('TileWindows', UseWeights)  : Use 'TileWindows' window layout
%        gui_layout('PositionFigure', hFigure, figArea, decorationSize)
%        gui_layout('ShowAllWindows')           : Show all figures

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2019

eval(macro_method);
end

%% ===== UPDATE LAYOUT =====
function Update() %#ok<DEFNU>
    % Cancel call in headless mode
    if (bst_get('GuiLevel') == -1)
        return;
    end
    % Call the appropriate function
    switch bst_get('Layout', 'WindowManager')
        case 'WeightWindows'
            % Call the tile windows layout
            TileWindows(1);
        case 'TileWindows'
            % Call the tile windows layout WITHOUT WEIGHTS CORRECTION
            TileWindows(0);
        case 'FullArea'
            [jArea, figArea] = GetScreenBrainstormAreas();
            FixedSizeWindows(figArea);
        otherwise
            % Nothing to do
    end
end

%% ===== GET DECORATION SIZE =====
% Get decorations size for a window: [left, top, right, bottom, menubarHeight, toolbarHeight]
function decorationSize = GetDecorationSize(jBstWindow)
    % Do not try to get it when running in server or nogui mode
    if (bst_get('GuiLevel') <= 0)
        decorationSize = [0 0 0 0 0 0];
        return;
    end
    % Get brainstorm window if needed
    if (nargin < 1) || isempty(jBstWindow)
        jBstWindow = bst_get('BstFrame');
    end
    % Get decorations size (left, top, right, bottom, menubarHeight, toolbarHeight) for a window
    decorationSize = [...
        jBstWindow.getRootPane.getBounds.getX(), ...
        jBstWindow.getRootPane.getBounds.getY(), ...
        jBstWindow.getBounds.getWidth() - jBstWindow.getRootPane.getBounds.getWidth() - jBstWindow.getRootPane.getBounds.getX(), ...
        jBstWindow.getBounds.getHeight() - jBstWindow.getRootPane.getBounds.getHeight() - jBstWindow.getRootPane.getBounds.getY(), ...
        20, ...% jBstWindow.getJMenuBar.getSize.getHeight(), ...
        28]; % TOOLBAR HEIGHT
    % For windows 10 and macos, remove the borders of the figures (they are transparent)
    if ispc && ~isempty(strfind(system_dependent('getos'), '10'))
        decorationSize(1) = 0;
        decorationSize(2) = 31;
        decorationSize(3) = 2;
        decorationSize(4) = 1;
    end
end

%% ===== GET MAXIMUM WINDOW =====
% Size and position of the rectangle usable by client applications on each screen.
%    - jMaxWindow : Position as seen by Java/Swing (real number of pixels on the screen)
%    - MaxWindow  : Position as seen by Matlab (scaled if there is a scale factor applied on high-DPI screens)
%
% Possible ways to get these information:
%    tk = java.awt.Toolkit.getDefaultToolkit();
%    ge = java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment();
%
%    1) get(0,'ScreenSize')
%       tk.getScreenSize()
%       ge.getDefaultScreenDevice().getDefaultConfiguration().getBounds()
%        - WINDOWS: current screen size (including toolbars...)
%        - LINUX  : total screen size (all the screens together, including toolbars...)
%
%    2) get(0,'MonitorPositions')
%     ** - WINDOWS: one entry per screen (including toolbars...)
%        - LINUX  : identical to get(0,'ScreenSize')
%
%    4) ge.getMaximumWindowBounds()
%     ** - WINDOWS: current screen size (EXCLUDING toolbars...)
%        - LINUX  : identical to get(0,'ScreenSize')
%        - MACOS  : Available area (corrects for position of dock and top menu)
%
%    6) ge.getScreenDevices()[].getDefaultConfiguration().getBounds()
%     ** - WINDOWS: size of each screen
%        - LINUX  : identical to get(0,'ScreenSize')
%
%    7) jFrame.setExtendedState(JFrame.MAXIMIZED_BOTH)
%        - WINDOWS: works well
%        - LINUX  : does not work on simple terminals (but works on CentOS/KDE workstations)
%
%    8) java.awt.Toolkit.getDefaultToolkit().getScreenInsets(jScreens(2).getDefaultConfiguration())
%        => Works for each screen separately
% function [jMaxWindow, MaxWindow, ZoomFactor] = GetMaximumWindow()
%     % Get graphic environment
%     ge = java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment();
%     % Get the default size
%     jMaxWindow = ge.getMaximumWindowBounds();
%     MaxWindow = [jMaxWindow.getX(), jMaxWindow.getY(), jMaxWindow.getWidth(), jMaxWindow.getHeight()];
%     % Try to get better for X11 systems
%     try
%         % Check if X11 classes exists
%         isX11 = strcmpi(class(ge), 'sun.awt.X11GraphicsEnvironment');
%         % In case of X11: rebuild information
%         if isX11
%             % Run xprop to get information about the desktop
%             panelWorkaround = java.lang.ProcessBuilder({'xprop', '-root', '-notype', '_NET_WORKAREA'});
%             proc = panelWorkaround.start();
%             % Get the result
%             br = java.io.BufferedReader(java.io.InputStreamReader(proc.getInputStream()));
%             res = br.readLine();
%             % Split to get the values
%             xprop = res.split('=');
%             xprop = str2num(char(xprop(2)));
%             % Read values
%             MaxWindow = xprop(1:4);
%         end
%     catch
%         % Nothing to do, just keep the regular values
%     end
%     % Return values in a Rectangle object
%     jMaxWindow = java.awt.Rectangle(MaxWindow(1), MaxWindow(2), MaxWindow(3), MaxWindow(4));
%     % Adjust with the OS scaling factor (for high-DPI screens)
%     javaScreenSize   = ge.getDefaultScreenDevice().getDisplayMode();
%     matlabScreenSize = get(0, 'ScreenSize');
%     ZoomFactor = double(javaScreenSize.getWidth()) ./ double(matlabScreenSize(3));
%     ZoomFactor = round(ZoomFactor * 100) / 100;
%     % For Matlab R2015b and above: the maximum window is not scaled coordinates, need to fix this
%     if (ZoomFactor >= 1.05) && (bst_get('MatlabVersion') >= 806)
%         MaxWindow = floor(MaxWindow ./ ZoomFactor);
%     end
% end

%% ===== GET CLIENT SCREEN SIZE =====
% Process information from GetMaximumWindow
function ScreenDef = GetScreenClientArea()
    % If running in headless mode, return a fake configuration
    if (bst_get('GuiLevel') == -1)
        ScreenDef.screenInsets = java.awt.Insets(0,0,0,0);
        ScreenDef.javaPos      = java.awt.Rectangle(1,1,1024,768);
        ScreenDef.matlabPos    = [0 0 1024 768];
        ScreenDef.zoomFactor   = 1;
        return;
    end
    % Get Java GraphicsEnvironment
    ge = java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment();
    tk = java.awt.Toolkit.getDefaultToolkit();
    % Get list of screens
    jScreens = ge.getScreenDevices();
    % Matlab monitor positions
    MonitorPositions = get(0, 'MonitorPositions');
    isOldPositions = (bst_get('MatlabVersion') < 804);
    % Fix discrepancies (reported in: https://neuroimage.usc.edu/forums/t/using-brainstorm-on-two-screens-under-linux/28418)
    if length(jScreens) > size(MonitorPositions,1)
        jScreens = jScreens(1:size(MonitorPositions,1));
    end
    % Find default screen
    % iDefaultScreen = ge.getDefaultScreenDevice().getScreen() + 1;   %%% CRASHES ON JAVA 7/Matlab2013 ON MACOSX
    jDefScreen = ge.getDefaultScreenDevice();
    iDefaultScreen = 1;
    for i = 2:length(jScreens)
        if (jScreens(i) == jDefScreen)
            iDefaultScreen = i;
        end
    end
    % Put the default screen in first position
    if (iDefaultScreen > 1)
        jScreens = [jScreens(iDefaultScreen), jScreens(1)];
        if (size(MonitorPositions,1) > 1)
            MonitorPositions = [MonitorPositions(iDefaultScreen,:); MonitorPositions(1,:)];
        end
    end
    % If there is a second screen but it is disabled: keep only the first one
    isDoubleScreen = bst_get('Layout', 'DoubleScreen');
    if (length(jScreens) > 1) && (isempty(isDoubleScreen) || ~isDoubleScreen)
        jScreens = jScreens(1);
        MonitorPositions = MonitorPositions(1,:);
    end

    % Get values for each screen
    for i = 1:length(jScreens)
        % === SCREEN SIZE ===
        % Get screen size
        jBounds = jScreens(i).getDefaultConfiguration().getBounds();
        % Get the screen insets
        jInsets = tk.getScreenInsets(jScreens(i).getDefaultConfiguration());
        ScreenDef(i).screenInsets = jInsets;
        % Maximum window size
        MaxWin = [jBounds.getX() + 1 + jInsets.left, ...
                  jBounds.getY() + 1 + jInsets.bottom, ...
                  jBounds.getWidth() - jInsets.left - jInsets.right, ...
                  jBounds.getHeight() - jInsets.top - jInsets.bottom];
        % Convert to a Java rectangle
        ScreenDef(i).javaPos = java.awt.Rectangle(MaxWin(1), MaxWin(2), MaxWin(3), MaxWin(4));
        % For Matlab < 2014b:  Adjust position of secondary monitors based on the size of the first one
        if isOldPositions
            % Secondary screens
            if (i > 1) && (i <= size(MonitorPositions,1))
                isVertical = (MonitorPositions(i,1) == MonitorPositions(1,1));
                % Vertical extension
                if isVertical && (MonitorPositions(i,2) ~= MonitorPositions(1,2))
                    MaxWin(1) = MaxWin(1) - MaxWin(3) + MonitorPositions(1,3);
                % Horizontal extension
                elseif (MonitorPositions(i,4) ~= MonitorPositions(1,4))
                    MaxWin(2) = MaxWin(2) - MaxWin(4) + MonitorPositions(1,4);
                end
            end
        % Newer Matlab: Rely on MonitorPositions directly
        else
            MaxWin = [MonitorPositions(i,1) + jInsets.left, ...
                      MonitorPositions(i,2) + jInsets.bottom, ...
                      MonitorPositions(i,3) - jInsets.left - jInsets.right, ...
                      MonitorPositions(i,4) - jInsets.top - jInsets.bottom];
        end

        % === SCALING ===
        try
            % Adjust with the OS scaling factor (for high-DPI screens)
            javaScreenSize = jScreens(i).getDefaultConfiguration().getBounds().getWidth();
            if (i == 1) || (i > size(MonitorPositions,1))
                matlabScreenSize = MonitorPositions(1,3) - MonitorPositions(1,1) + 1;
            else
                % For Matlab < 2014b:
                if isOldPositions
                    matlabScreenSize = MonitorPositions(i,3) - (MonitorPositions(i,1) - MonitorPositions(i-1,3)) + 1;
                % Newer Matlab
                else
                    matlabScreenSize = MonitorPositions(i,3);
                end
            end
            ZoomFactor = double(javaScreenSize) ./ double(matlabScreenSize);
            ZoomFactor = round(ZoomFactor * 100) / 100;
%             % For Matlab R2015b and above: the screen size not in scaled coordinates, need to fix this
%             if (ZoomFactor >= 1.05) && (bst_get('MatlabVersion') >= 806)
%                 MaxWin = floor(MaxWin ./ ZoomFactor);
%             end
        catch
            ZoomFactor = 1;
        end
        % Save maximum window position and zoom factor
        ScreenDef(i).matlabPos = MaxWin;
        ScreenDef(i).zoomFactor = ZoomFactor;
    end
end

%% ===== GET BRAINSTORM SCREEN AREAS =====
% Two workspace configurations are available :
%      .---------------.           .---------------.
%   1) | BST | Figures |   OR   2) | Figures | BST |
%      '---------------'           '---------------'
function [jBstArea, FigArea, nbScreens, jFigArea, jInsets] = GetScreenBrainstormAreas(jBstWindow)
    % Jave window not provided
    if (nargin < 1) || isempty(jBstWindow)
        jBstWindow = bst_get('BstFrame');
    end
    % Get client area in whole screen
    ScreenDef = bst_get('ScreenDef');
    nbScreens = length(ScreenDef);

    % ===== NO GUI OR FULLSCREEN =====
    if ~bst_get('isGUI') || bst_get('Layout', 'FullScreen')
        jBstArea = java.awt.Rectangle(0,0,0,0);
        FigArea  = ScreenDef(1).matlabPos;
        jFigArea = ScreenDef(1).javaPos;
        jInsets  = ScreenDef(1).screenInsets;

    % ===== ONE SCREEN =====
    elseif (nbScreens == 1)
        javaArea   = ScreenDef.javaPos;
        matlabArea = ScreenDef.matlabPos;
        jInsets    = ScreenDef.screenInsets;
        ZoomFactor = ScreenDef.zoomFactor;
        % Check that Brainstorm window is completely inside the client area
        fixedBounds = jBstWindow.getBounds.intersection(javaArea);
        % If the window is outside of the principal screen: reset its default positions
        if ((fixedBounds.getWidth() < 200) || (fixedBounds.getHeight() < 300))
            fixedBounds.setRect(...
                javaArea.getX(), ...
                javaArea.getY(), ...
                round(450 * bst_get('InterfaceScaling') / 100), ...
                javaArea.getHeight() .* .9);
        end
        jBstWindow.setBounds(fixedBounds);
        % Compute distance from main BST window to borders
        distanceToLeft  = jBstWindow.getBounds.getX() - javaArea.getX();
        distanceToRight = javaArea.getWidth() ...
                          - (jBstWindow.getBounds.getX() + jBstWindow.getBounds.getWidth());
        % === CONFIG 1 ===
        % If BST window is more on the left : configuration 1)
        if (distanceToLeft < distanceToRight)
            % Get the area covered by the Brainstorm window
            jBstArea = java.awt.Rectangle(...
                            javaArea.getX(), ...
                            javaArea.getY(), ...
                            jBstWindow.getBounds.getX() + jBstWindow.getBounds.getWidth() - javaArea.getX(), ...
                            javaArea.getHeight());
            % Get the limit of the Brainstorm window
            bstWinX = (jBstArea.getX() + jBstArea.getWidth() + 1);
            % For Matlab R2015b and above: the screen size not in scaled coordinates, need to fix this
            if (ZoomFactor >= 1.05) && (bst_get('MatlabVersion') >= 806)
                bstWinX = floor(bstWinX ./ ZoomFactor);
            end
            FigArea = [bstWinX, ...
                       matlabArea(2), ...
                       matlabArea(1) + matlabArea(3) - bstWinX, ...
                       matlabArea(4)];
        % === CONFIG 2 ===
        else
            jBstArea = java.awt.Rectangle(...
                            jBstWindow.getBounds.getX(), ...
                            javaArea.getY(), ...
                            javaArea.getX() + javaArea.getWidth() - jBstWindow.getBounds.getX(), ...
                            javaArea.getHeight());
            % Get the limit of the Brainstorm window
            bstWinX = jBstArea.getX() - javaArea.getX();
            % For Matlab R2015b and above: the screen size not in scaled coordinates, need to fix this
            if (ZoomFactor >= 1.05) && (bst_get('MatlabVersion') >= 806)
                bstWinX = floor(bstWinX ./ ZoomFactor);
            end
            FigArea = [matlabArea(1), ...
                       matlabArea(2), ...
                       bstWinX, ...
                       matlabArea(4)];
        end
        jFigArea = java.awt.Rectangle(floor((FigArea(1)-1) .* ZoomFactor), ...
                                      javaArea.getY(), ...
                                      floor(FigArea(3) .* ZoomFactor), ...
                                      floor(FigArea(4) .* ZoomFactor));

    % ===== TWO SCREENS (OR MORE) =====
    elseif (nbScreens >= 2)
        % If Brainstorm window is on screen 2
        tol = 30;
        if (jBstWindow.getBounds.getX() >= ScreenDef(2).javaPos.getX() - tol) && (jBstWindow.getBounds.getX() <= ScreenDef(2).javaPos.getX() + ScreenDef(2).javaPos.getWidth() + tol) && ...
           (jBstWindow.getBounds.getY() >= ScreenDef(2).javaPos.getY() - tol) && (jBstWindow.getBounds.getY() <= ScreenDef(2).javaPos.getY() + ScreenDef(2).javaPos.getHeight() + tol)
            jBstArea = ScreenDef(2).javaPos;
            FigArea  = ScreenDef(1).matlabPos;
            jFigArea = ScreenDef(1).javaPos;
            jInsets  = ScreenDef(1).screenInsets;
        % If Brainstorm window is on screen 1
        else
            jBstArea = ScreenDef(1).javaPos;
            FigArea  = ScreenDef(2).matlabPos;
            jFigArea = ScreenDef(2).javaPos;
            jInsets  = ScreenDef(2).screenInsets;
        end
    end
end

%% ===== GET FIGURE GROUPS =====
function Figures = GetFigureGroups(isSkipMriViewer)
    global GlobalData;
    % Parse inputs
    if (nargin < 1) || isempty(isSkipMriViewer)
        isSkipMriViewer = 1;
    end

    % Get visualization figures list
    Figures = repmat(struct('fDataTimeSeries',    [], ...
                            'fRawViewer',         [], ...
                            'fTopography',        [], ...
                            'f3DViz',             [], ...
                            'fResultsTimeSeries', [], ...
                            'fTimefreq',          [], ...
                            'fSpectrum',          [], ...
                            'fOther',             [], ...
                            'fMriViewer',         [], ...
                            'fConnect',           [], ...
                            'fPac',               [], ...
                            'fImage',             [], ...
                            'fVideo',             [], ...
                            'iDataSet',           [], ...
                            'VerticalWeight',     [], ...
                            'nFigures',           [], ...
                            'nRawChannels',       []), 0);
    % Search in all the DataSets
    for iDS = 1:length(GlobalData.DataSet)
        dsFigs = GlobalData.DataSet(iDS).Figure;
        isRaw = strcmpi(GlobalData.DataSet(iDS).Measures.DataType, 'raw');
        % If there are visible figures registered in this DataSet
        if ~isempty(dsFigs) && any(~strcmpi(get([dsFigs.hFigure], 'WindowStyle'), 'docked'))
            % Detect raw viewer figures: Create one block of figure PER RAW VIEWER
            iRawViewer = [];
            if isRaw
                for iFig = 1:length(dsFigs)
                    if strcmpi(dsFigs(iFig).Id.Type, 'DataTimeSeries') && ~isempty(dsFigs(iFig).Id.Modality) && (dsFigs(iFig).Id.Modality(1) ~= '$')
                        iRawViewer = [iRawViewer, iFig];
                        % Create a new group
                        if isempty(Figures) || isempty(Figures(iFigureGroups).fRawViewer) || (Figures(iFigureGroups).iDataSet ~= iDS)
                            iFigureGroups = length(Figures) + 1;
                            Figures(iFigureGroups).iDataSet       = iDS;
                            Figures(iFigureGroups).VerticalWeight = 4;
                            Figures(iFigureGroups).fRawViewer     = dsFigs(iFig).hFigure;
                            Figures(iFigureGroups).nFigures = 1;
                        % Add figure to existing group
                        else
                            Figures(iFigureGroups).fRawViewer(end+1) = dsFigs(iFig).hFigure;
                            Figures(iFigureGroups).nFigures = Figures(iFigureGroups).nFigures + 1;
                        end
                        % EEG, MEG: Ration 3:1 with respect with isolated channels
                        Figures(iFigureGroups).nRawChannels(end+1) = length(dsFigs(iFig).SelectedChannels);
                    end
                end
                % Remove raw viewer figures from the list of figures
                if ~isempty(iRawViewer)
                    dsFigs(iRawViewer) = [];
                end
                % No more figures in this dataset: skip to the next one
                if isempty(dsFigs)
                    continue;
                end
            end

            % If MriViewer figure only in this dataset: put in in the first figure group always (does not count as a figure group)
            if isSkipMriViewer && (length(dsFigs) == 1) && strcmpi(dsFigs(1).Id.Type, 'MriViewer')
                iFigureGroups = 1;
            else
                iFigureGroups = length(Figures) + 1;
                Figures(iFigureGroups).nFigures = 0;
            end
            % Get all figures for this figure group
            Figures(iFigureGroups).iDataSet = iDS;
            for iFig = 1:length(dsFigs)
                % Figure must not be docked
                if strcmpi(get(dsFigs(iFig).hFigure,'WindowStyle'), 'docked')
                    continue;
                end
                % Add figure to the list
                Figures(iFigureGroups).(['f' dsFigs(iFig).Id.Type]) = [Figures(iFigureGroups).(['f' dsFigs(iFig).Id.Type]), dsFigs(iFig).hFigure];
                Figures(iFigureGroups).VerticalWeight = 1;
                Figures(iFigureGroups).nFigures = Figures(iFigureGroups).nFigures + 1;
            end
        end
    end

    % ===== GROUP ISOLATED FIGURES =====
    % Combine together lists of isolated figure groups (example: many results from different datasets in memory)
    % => group them so that they can be displayed as a mosaique instead of flat vertical list
    if (length(Figures) > 3) && all([Figures.nFigures] == 1) && ...
            (~isempty(Figures(1).fMriViewer) || ~isempty(Figures(1).fTopography) || ~isempty(Figures(1).f3DViz)    || ~isempty(Figures(1).fResultsTimeSeries) || ...
             ~isempty(Figures(1).fTimefreq)   || ~isempty(Figures(1).fSpectrum) || ~isempty(Figures(1).fConnect) || ~isempty(Figures(1).fPac) || ~isempty(Figures(1).fImage) || ~isempty(Figures(1).fVideo))
        uniqueFigures = Figures(1);
        uniqueFigures.fMriViewer         = [Figures.fMriViewer];
        uniqueFigures.fTopography        = [Figures.fTopography];
        uniqueFigures.f3DViz             = [Figures.f3DViz];
        uniqueFigures.fRawViewer         = [Figures.fRawViewer];
        uniqueFigures.fDataTimeSeries    = [Figures.fDataTimeSeries];
        uniqueFigures.fResultsTimeSeries = [Figures.fResultsTimeSeries];
        uniqueFigures.fTimefreq          = [Figures.fTimefreq];
        uniqueFigures.fSpectrum          = [Figures.fSpectrum];
        uniqueFigures.fConnect           = [Figures.fConnect];
        uniqueFigures.fPac               = [Figures.fPac];
        uniqueFigures.fImage             = [Figures.fImage];
        uniqueFigures.fVideo             = [Figures.fVideo];
        uniqueFigures.fOther             = [Figures.fOther];
        uniqueFigures.nFigures           = length(Figures);
        Figures = uniqueFigures;
    end
end

%% ===== TILE WINDOWS =====
% Divide vertically the available area for figures between the DataSets
% For each DataSet with figures :
%     - One line for DataTimeSeries
%     - One line for TopoGraphy and 3DViz figure
%     - One line for ResultsTimeSeries
function TileWindows(UseWeights)
    % Get Brainstorm window
    jBstWindow = bst_get('BstFrame');
    % Get the figures
    isSkipMriViewer = 0;
    Figures = GetFigureGroups(isSkipMriViewer);
    % Get Brainstorm Main window and figures areas
    [jBstArea, FigArea, nbScreens] = GetScreenBrainstormAreas(jBstWindow);
    % Check that the space to display figures is sufficient, else use the whole screen
    if (FigArea(3) < 250)
        % Use all whole screen
        FigArea = [1, FigArea(2), jBstArea.getWidth() + FigArea(3), FigArea(4)];
    end
    % Get decoration dimensions on the current operating system (left, top, right, bottom, menubarHeight)
    decorationSize = bst_get('DecorationSize');

    % For each Figures block
    nbBlocks = length(Figures);
    % Define maximum size for each figure type (in pixels)
    if UseWeights || (nbBlocks > 1)
        maxSizeOtherFigures = [4096, 4096];
    elseif (nbScreens == 1)
        maxSizeOtherFigures = [700, 600];
    else
        maxSizeOtherFigures = [2048, 600];
    end

    % Total weights
    totalWeight = sum([Figures.VerticalWeight]);
    % Process each block
    for iBlock = 1:nbBlocks
        if UseWeights
            blockWeight = Figures(iBlock).VerticalWeight ./ totalWeight;
            prevWeight = sum([Figures(iBlock+1:end).VerticalWeight]) ./ totalWeight;
        else
            blockWeight = 1/nbBlocks;
            prevWeight = 1 - iBlock/nbBlocks;
        end
        % Define the area reserved for this figures block
        blockArea = [FigArea(1), ...
                     FigArea(2) + round(FigArea(4) * prevWeight), ...
                     FigArea(3), ...
                     round(FigArea(4) * blockWeight)];

        % ===== 3DViz, Topography and ResultsTimeSeries figures =====
        OtherFigures = [Figures(iBlock).fDataTimeSeries, Figures(iBlock).fTopography, ...
                        Figures(iBlock).f3DViz, Figures(iBlock).fResultsTimeSeries, ...
                        Figures(iBlock).fTimefreq, Figures(iBlock).fSpectrum, Figures(iBlock).fConnect,...
                        Figures(iBlock).fPac, Figures(iBlock).fImage, Figures(iBlock).fVideo, Figures(iBlock).fOther];
        % Add MRI Viewer
        if ~isSkipMriViewer
            OtherFigures = [OtherFigures, Figures(iBlock).fMriViewer];
        end
        nbOtherFigures = length(OtherFigures);
        % Always display full screen MRI viewer if it is the only figure
        isFullScreenMri = isSkipMriViewer || ((nbBlocks == 1) && (nbOtherFigures == 1));
        % Position regular figures
        if (nbOtherFigures >= 1)
            % Divide list of figures in many rows
            if (nbBlocks == 1) && (nbOtherFigures == 3)
                nbCols = 2;
                nbRows = 2;
            else
                if (blockArea(3) <= blockArea(4))
                    nbCols = floor(sqrt(nbOtherFigures));
                    nbRows = ceil(nbOtherFigures / nbCols);
                elseif (blockArea(3) > blockArea(4))
                    nbRows = floor(sqrt(nbOtherFigures));
                    nbCols = ceil(nbOtherFigures / nbRows);
                end
            end

            % Divide horizontally the available space between all DataTimeSeries figures
            for iOtherFigures = 1:nbOtherFigures
                % If this figure is a full screen MRI figure: skip
                if isFullScreenMri && ~isempty(Figures(iBlock).fMriViewer) && ismember(OtherFigures(iOtherFigures), Figures(iBlock).fMriViewer)
                    continue;
                end
                % Get figure row and column
                iRow = floor((iOtherFigures-1) / nbCols);
                iCol = mod((iOtherFigures-1), nbCols);

                % Compute new position for this DataTimeSeries figure
                figArea = [blockArea(1) + iCol * round(blockArea(3) / nbCols), ...
                           blockArea(2) + (nbRows - iRow - 1) * round(blockArea(4) / nbRows), ...
                           round(blockArea(3) / nbCols), ...
                           round(blockArea(4) / nbRows)];

                % Limit figure size
                if (nbOtherFigures == 1) && ~ismember(OtherFigures, Figures(iBlock).fDataTimeSeries) && ~ismember(OtherFigures, Figures(iBlock).fResultsTimeSeries)
                    figArea(3) = min(figArea(3), maxSizeOtherFigures(1));
                end
                figArea(2) = figArea(2) + (figArea(4) - min(figArea(4), maxSizeOtherFigures(2)));
                figArea(4) = min(figArea(4), maxSizeOtherFigures(2));

                % Apply this new position
                PositionFigure(OtherFigures(iOtherFigures), figArea, decorationSize);
            end
        end

        % === RAW VIEWER ===
        % Get the display weight of each view
        if UseWeights && ~isempty(Figures(iBlock).nRawChannels)
            rawWeight = [Figures(iBlock).nRawChannels] ./ max([Figures(iBlock).nRawChannels]);
            rawWeight(rawWeight == 1)  = 7;
            rawWeight(rawWeight < .1) = 1;
            rawWeight(rawWeight < .2) = 4;
            rawWeight(rawWeight < .4) = 5;
            rawWeight(rawWeight <  1)  = 6;
        else
            rawWeight = ones(1, length(Figures(iBlock).fRawViewer));
        end
        % Position each raw viewer
        for iRawFig = 1:length(Figures(iBlock).fRawViewer)
            % If only one block: fill vertically AND horizontally
            if (nbBlocks == 1)
                rawArea = FigArea;
            % If many blocks: fill only horizontally (full block area)
            else
                rawArea = blockArea;
            end
            % Assign a weight to each figure based on the number of sensors on each view
            curHeight = round(rawWeight(iRawFig) ./ sum(rawWeight) .* rawArea(4));
            prevHeight = round(sum(rawWeight(iRawFig+1:end)) ./ sum(rawWeight) .* rawArea(4));
            % Calculate position
            figRawArea = [rawArea(1), ...
                          rawArea(2) + prevHeight, ...
                          rawArea(3), ...
                          curHeight];
            % Set figure position
            PositionFigure(Figures(iBlock).fRawViewer(iRawFig), figRawArea, decorationSize);
        end
        % === MRI VIEWER FIGURES ===
        % Set to the whole avaiable space
        if isSkipMriViewer || ((nbBlocks == 1) && (nbOtherFigures == 1))
            for iMriFig = 1:length(Figures(iBlock).fMriViewer)
                PositionFigure(Figures(iBlock).fMriViewer(iMriFig), FigArea, decorationSize);
            end
        end
    end

    % === MASK EDITOR ===
    % Set to the entire avaiable space
    global gEditMaskData;
    if ~isempty(gEditMaskData) && isfield(gEditMaskData, 'Handles') && ishandle(gEditMaskData.Handles.mri_editMask)
        PositionFigure(gEditMaskData.Handles.mri_editMask, FigArea, decorationSize);
    end
end

%% ===== FIXED LAYOUT =====
% Set the same size/position for all the figures
function FixedSizeWindows(figArea)
    % Get the figures
    Figures = GetFigureGroups();
    % Get decoration dimensions on the current operating system (left, top, right, bottom, menubarHeight)
    decorationSize = bst_get('DecorationSize');
    % Process each block of figures
    for iBlock = 1:length(Figures)
        % Get all the figures
        hAllFig = [Figures(iBlock).fDataTimeSeries, Figures(iBlock).fTopography, ...
                   Figures(iBlock).f3DViz, Figures(iBlock).fResultsTimeSeries, ...
                   Figures(iBlock).fTimefreq, Figures(iBlock).fSpectrum, Figures(iBlock).fConnect, Figures(iBlock).fPac, Figures(iBlock).fImage, Figures(iBlock).fVideo, ...
                   Figures(iBlock).fOther, Figures(iBlock).fRawViewer, Figures(iBlock).fMriViewer];
        % Set the position
        for iFig = 1:length(hAllFig)
            PositionFigure(hAllFig(iFig), figArea, decorationSize);
        end
    end
end

%% ===== POSITION FIGURES =====
function PositionFigure(hFigure, figArea, decorationSize)
    drawnow;
    if ~ishandle(hFigure)
        return
    end
    % If figure has a menu bar
    if ~strcmpi(get(hFigure, 'MenuBar'), 'none')
        menubarHeight = decorationSize(5);
    else
        menubarHeight = 0;
    end
    % If figure has a figure toolbar
    toolbarHeight = 0;
    if ~strcmpi(get(hFigure, 'ToolBar'), 'none')
        toolbarHeight = toolbarHeight + decorationSize(6);
    end
    if ~isempty(findobj(hFigure, '-depth', 1, 'Tag', 'AlignToolbar'))
        toolbarHeight = toolbarHeight + decorationSize(6);
    end
    % If figure has a plot edit toolbar
    isPlotEditToolbar = getappdata(hFigure, 'isPlotEditToolbar');
    if ~isempty(isPlotEditToolbar) && isPlotEditToolbar
        toolbarHeight = toolbarHeight + decorationSize(6);
    end
    % Remove decorations dimensions
    figDim = figArea + [decorationSize(1), ...
                        decorationSize(4), ...
                        - decorationSize(1) - decorationSize(3), ...
                        - decorationSize(2) - decorationSize(4) - menubarHeight - toolbarHeight];
    % Older versions of Matlab/MacOSX have random 22px positionning problem (fixed for versions >= R2011b)
    if strncmp(computer, 'MAC', 3) && (bst_get('MatlabVersion') <= 712)
        figDim(2) = figDim(2) - 22;
    end
    % Force all figures to have positive dimensions
    figDim(3:4) = max(figDim(3:4), [10 10]);
%     % Get client area in whole screen
%     ScreenDef = bst_get('ScreenDef');
%     % If the figure is on a secondary screen and the two screens do not have the same size, need to fix the "Y=0" of the screen
%     if (length(ScreenDef) > 1) && (figDim(1) >= ScreenDef(1).matlabPos(1) + ScreenDef(1).matlabPos(3)) && (ScreenDef(1).matlabPos(4) ~= ScreenDef(2).matlabPos(4))
%         figDim(2) = figDim(2) - ScreenDef(2).matlabPos(4) + ScreenDef(1).matlabPos(4);
%     end
    % Apply position to figure
    set(hFigure, 'Position', figDim);
end

%% ===== SHOW ALL WINDOWS =====
function ShowAllWindows() %#ok<DEFNU>
    % Get Brainstorm window
    jBstWindow = bst_get('BstFrame');
    % Get the figures
    Figures = GetFigureGroups();
    % Set focus to each figure sequentially
    if ~isempty(Figures)
        fieldNames = {'fDataTimeSeries', 'fRawViewer', 'fTopography', 'f3DViz', 'fResultsTimeSeries', 'fTimefreq', 'fSpectrum', 'fMriViewer', 'fConnect', 'fPac', 'fImage', 'fVideo', 'fOther'};
        nbBlocks = length(Figures);
        for iBlock = 1:nbBlocks
            for iType = 1:length(fieldNames)
                if ~isempty(Figures(iBlock).(fieldNames{iType}))
                    for hFig = Figures(iBlock).(fieldNames{iType})
                        if strcmpi(get(hFig, 'Visible'), 'on')
                            figure(hFig);
                        end
                    end
                end
            end
        end
    end
    % Put focus on the main brainstorm window again
    drawnow;
    jBstWindow.setVisible(1);
end

%% ===================================================================================================================================
%  ===== DISPLAY SETUPS ==============================================================================================================
%  ===================================================================================================================================

%% ===== CREATE SETUP MENU =====
function SetupMenu(jMenu) %#ok<DEFNU>
    import org.brainstorm.icon.*;
%     % Get all the open figures
%     hAllFig = bst_figures('GetAllFigures');
%     if isempty(hAllFig)
%         return;
%     end
    fontSize = [];
    % Get current setups
    UserSetups = bst_get('Layout', 'UserSetups');
    % List all the pipelines
    for iSetup = 1:length(UserSetups)
        gui_component('MenuItem', jMenu, [], UserSetups(iSetup).Name, IconLoader.ICON_LAYOUT_CASCADE, [], @(h,ev)LoadSetup(iSetup), fontSize);
    end
    % Separator
    if ~isempty(UserSetups)
        jMenu.addSeparator();
    end
    % Create new setup
    gui_component('MenuItem', jMenu, [], 'New setup', IconLoader.ICON_SAVE, [], @(h,ev)CreateNewSetup(), fontSize);
    % Delete entries
    if ~isempty(UserSetups)
        jMenuDel = gui_component('Menu', jMenu, [], 'Delete setup', IconLoader.ICON_DELETE, [], [], fontSize);
        % List all the pipelines
        for iSetup = 1:length(UserSetups)
            gui_component('MenuItem', jMenuDel, [], UserSetups(iSetup).Name, IconLoader.ICON_DELETE, [], @(h,ev)DeleteSetup(iSetup), fontSize);
        end
    end
end

%% ===== CREATE NEW SETUP =====
function CreateNewSetup()
    global GlobalData;
    % Get layout structure
    UserSetups = GlobalData.Preferences.Layout.UserSetups;
    % Create new structure
    sSetup.Name = '';
    sSetup.Figures = repmat(struct(...
        'FigureId', [], ...
        'AppData',  [], ...
        'Position', [], ...
        'Color',    [], ...
        'Camera',   []), 0);
    DataFile = [];
    nWarningSkip = 0;
    % Loop on all the figures
    for iDS = 1:length(GlobalData.DataSet)
        for iFig = 1:length(GlobalData.DataSet(iDS).Figure)
            % Get figure infor
            Figure = GlobalData.DataSet(iDS).Figure(iFig);
            AppData = getappdata(Figure.hFigure);
            % Keep only simple figures with a data file
            if ~ismember(Figure.Id.Type, {'DataTimeSeries','3DViz','Topography'}) || isempty(AppData.DataFile) || (isfield(AppData, 'Timefreq') && ~isempty(AppData.Timefreq)) || strcmpi(file_gettype(AppData.DataFile), 'pdata')
                nWarningSkip = nWarningSkip + 1;
                continue;
            end
            % Remove unwanted objects from the AppData structure
            for field = fieldnames(AppData)'
                if ~isempty(strfind(field{1}, 'uitools')) || isjava(AppData.(field{1})) || ismember(field{1}, {'SubplotDefaultAxesLocation', 'SubplotDirty'})
                    AppData = rmfield(AppData, field{1});
                    continue;
                end
            end
            % Keep only the figures with the same DataFile as the first one
            if ~isempty(DataFile)
                if ~file_compare(DataFile, AppData.DataFile)
                    nWarningSkip = nWarningSkip + 1;
                    continue;
                end
            else
                DataFile = AppData.DataFile;
            end
            % Save figure
            iSaveFig = length(sSetup.Figures) + 1;
            sSetup.Figures(iSaveFig).FigureId = Figure.Id;
            sSetup.Figures(iSaveFig).AppData  = AppData;
            sSetup.Figures(iSaveFig).Position = get(Figure.hFigure, 'Position');
            sSetup.Figures(iSaveFig).Color    = get(Figure.hFigure, 'Color');
            % 3D figures: Get camera
            hAxes = findobj(Figure.hFigure, '-depth', 1, 'Tag', 'Axes3D');
            if strcmpi(Figure.Id.Type, '3DViz') && ~isempty(hAxes)
                % Copy view angle and camup
                [cam.az, cam.el] = view(hAxes);
                cam.up = camup(hAxes);
                sSetup.Figures(iSaveFig).Camera = cam;
            else
                sSetup.Figures(iSaveFig).Camera = [];
            end
        end
    end
    % Check the number of figures
    if (length(sSetup.Figures) < 1)
        bst_error('Not enough figures related with the same dataset to create a screen setup.', 'Screen setup', 0);
        return;
    elseif (nWarningSkip > 0)
        disp(['BST> Warning: ' num2str(nWarningSkip) ' figure(s) could not be included in the saved setup.']);
    end
    % Ask user the name for the new setup
    newName = java_dialog('input', 'Enter a name for the new screen setup:', 'Screen setup');
    if isempty(newName)
        return;
    end
    % Check if setup already exists
    if ~isempty(UserSetups) && any(strcmpi({UserSetups.Name}, newName))
        bst_error('This name already exists.', 'Screen setup', 0);
        return
    end
    % Save new setup name
    sSetup.Name = newName;
    % Add setup to list
    if isempty(UserSetups)
        UserSetups = sSetup;
    else
        UserSetups(end+1) = sSetup;
    end
    % Save modifications
    GlobalData.Preferences.Layout.UserSetups = UserSetups;
end

%% ===== DELETE SETUP =====
function DeleteSetup(iSetup)
    global GlobalData;
    % Ask confirmation
    if ~java_dialog('confirm', ['Delete setup "' GlobalData.Preferences.Layout.UserSetups(iSetup).Name '"?'], 'Screen setup');
        return;
    end
    % Select first item in the pipeline
    GlobalData.Preferences.Layout.UserSetups(iSetup) = [];
end

%% ===== LOAD SETUP =====
function LoadSetup(iSetup)
    global GlobalData;
    % Get all the loaded data files
    iDS = [];
    DataFile = [];
    for i = 1:length(GlobalData.DataSet)
        if ~isempty(GlobalData.DataSet(i).DataFile) && ~isempty(DataFile) && ~file_compare(GlobalData.DataSet(i).DataFile, DataFile)
            disp(['BST> Warning: The screen setup is applied only to the first loaded dataset, skipping data file: ' GlobalData.DataSet(i).DataFile]);
        elseif ~isempty(GlobalData.DataSet(i).DataFile)
            iDS = i;
            DataFile = GlobalData.DataSet(iDS).DataFile;
        end
    end
    % Check for loaded data
    if isempty(DataFile)
        bst_error('You have to load one dataset before applying a screen setup.', 'Screen setup', 0);
        return
    end
    % Get all the old figures
    hFigOldAll = bst_figures('GetAllFigures');
    % Get screen setup to load
    sSetup = GlobalData.Preferences.Layout.UserSetups(iSetup);
    % Disable temporarily the layout manager
    curWM = bst_get('Layout', 'WindowManager');
    if ~isempty(curWM)
        bst_set('Layout', 'WindowManager', sSetup.Name);
    end
    % Reload all the figures
    hNewFig = [];
    for i = 1:length(sSetup.Figures)
        % Progress bar
        bst_progress('start', 'Screen setup', 'Loading figures...');
        % Re-use existing figures: Get existing figures with the same FigureId
        [hFigOld, iFigOld] = bst_figures('GetFigure', iDS, sSetup.Figures(i).FigureId);
        hFigOld = setdiff(hFigOld, hNewFig);
        if ~isempty(hFigOld)
            hFig = hFigOld(1);
            iFig = iFigOld(1);
        % Else: Force creation
        else
            [hFig, iFig] = bst_figures('CreateFigure', iDS, sSetup.Figures(i).FigureId, 'AlwaysCreate');
        end
        if isempty(hFig)
            disp(['BST> Error: Could not create figure #' num2str(i)]);
            continue;
        end
        hNewFig(end+1) = hFig;
        % Update AppData structure
        AppData = sSetup.Figures(i).AppData;
        if isfield(AppData, 'DataFile') && ~isempty(AppData.DataFile)
            AppData.DataFile = DataFile;
        end
        if isfield(AppData, 'TsInfo') && isfield(AppData.TsInfo, 'FileName') && ~isempty(AppData.TsInfo.FileName)
            AppData.TsInfo.FileName = DataFile;
            AppData.TsInfo = struct_copy_fields(AppData.TsInfo, db_template('TsInfo'), 0);
        end
        if isfield(AppData, 'TopoInfo') && isfield(AppData.TopoInfo, 'FileName') && ~isempty(AppData.TopoInfo.FileName)
            AppData.TopoInfo.FileName = DataFile;
            AppData.TopoInfo = struct_copy_fields(AppData.TopoInfo, db_template('TopoInfo'), 0);
        end
        if isfield(AppData, 'Surface') && ~isempty(AppData.Surface)
            % Save old surfaces
            Surface = AppData.Surface;
            % Remove all existing surfaces
            TessInfo = getappdata(hFig, 'Surface');
            for iDel = 1:length(TessInfo)
                panel_surface('RemoveSurface', hFig, iDel);
            end
            AppData.Surface = repmat(db_template('TessInfo'),0);
        else
            Surface = [];
        end
        if isfield(AppData, 'StudyFile') && ~isempty(AppData.StudyFile)
            sStudy = bst_get('AnyFile', DataFile);
            AppData.StudyFile = sStudy.FileName;
            sSubject = bst_get('Subject', sStudy.BrainStormSubject);
            AppData.SubjectFile = sSubject.FileName;
        end
        % Set AppData
        for field = fieldnames(AppData)'
            setappdata(hFig, field{1}, AppData.(field{1}));
        end
        % Position figure
        set(hFig, 'Position', sSetup.Figures(i).Position);
        % Switch according to figure type
        switch (sSetup.Figures(i).FigureId.Type)
            case '3DViz'
                % Loop on the surface files to add
                for iSurf = 1:length(Surface)
                    % Get source file for the data
                    if ~isempty(Surface(iSurf).DataSource) && strcmpi(Surface(iSurf).DataSource.Type, 'Source')
                        % Get first results file for data file
                        [sStudy, iStudy, iResults] = bst_get('ResultsForDataFile', DataFile);
                        if isempty(iResults)
                            disp('BST> Warning: Data file does not have any attached source file.');
                            close(hFig);
                            continue;
                        end
                        ResultsFile = sStudy.Result(iResults(1)).FileName;
                        % Get surface file
                        ResultsMat = in_bst_results(ResultsFile, 0, 'SurfaceFile');
                        % Add surface
                        panel_surface('AddSurface', hFig, ResultsMat.SurfaceFile);
                        % Load source file
                        bst_memory('LoadResultsFile', ResultsFile);
                        % Set surface data
                        panel_surface('SetSurfaceData', hFig, iSurf, 'Source', ResultsFile, 0);
                        % Set surface properties
                        panel_surface('SetSurfaceSmooth',       hFig, iSurf, Surface(iSurf).SurfSmoothValue, 0);
                        panel_surface('SetSurfaceTransparency', hFig, iSurf, Surface(iSurf).SurfAlpha);
                        panel_surface('SetShowSulci',           hFig, iSurf, Surface(iSurf).SurfShowSulci);
                        panel_surface('SetDataThreshold',       hFig, iSurf, Surface(iSurf).DataThreshold);
                        panel_surface('SetSizeThreshold',       hFig, iSurf, Surface(iSurf).SizeThreshold);
                    else
                        % Get surface type
                        SurfaceType = Surface(iSurf).Name;
                        % Get subject structure
                        sStudy = bst_get('DataFile', DataFile);
                        sSubject = bst_get('Subject', sStudy.BrainStormSubject);
                        if ismember(SurfaceType, {'Cortex', 'Scalp', 'InnerSkull', 'InnerSkull'})
                            SurfaceFile = sSubject.Surface(sSubject.(['i' SurfaceType])).FileName;
                        elseif strcmpi(SurfaceType, 'Anatomy')
                            SurfaceFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
                        else
                            SurfaceFile = [];
                        end
                        % Add default surface
                        if ~isempty(SurfaceFile)
                            panel_surface('AddSurface', hFig, SurfaceFile);
                        else
                            disp(['BST> Warning: Cannot reload surface #' num2str(iSurf) '.']);
                        end
                    end
                end
                % Find 3D axes
                hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');
                if ~isempty(hAxes)
                    % Set camera
                    view(hAxes, sSetup.Figures(i).Camera.az, sSetup.Figures(i).Camera.el);
                    camup(hAxes, sSetup.Figures(i).Camera.up);
                    % Update head light position
                    camlight(findobj(hAxes, '-depth', 1, 'Tag', 'FrontLight'), 'headlight');
                end
                % Make it visible
                set(hFig, 'Visible', 'on');

            case 'Topography'
                % Plot figure
                figure_topo('PlotFigure', iDS, iFig, 1);
                % Make it visible
                set(hFig, 'Visible', 'on');

            case 'DataTimeSeries'
                view_timeseries(DataFile, sSetup.Figures(i).FigureId.Modality, AppData.TsInfo.RowNames, hFig);

            otherwise
                bst_figures('ReloadFigures', hFig);
        end
        % Set background
        if isfield(sSetup.Figures(i), 'Color') && ~isempty(sSetup.Figures(i).Color)
            set(hFig, 'Color', sSetup.Figures(i).Color);
        end
    end
    % Close all the figures that haven't been updated
    hFigClose = setdiff(hFigOldAll, hNewFig);
    if ~isempty(hFigClose)
        close(hFigClose);
    end
    % Close progress bar
    bst_progress('stop');
end

%% ===== UPDATE MAXIMUM BST WINDOW SIZE =====
function UpdateMaxBstSize() %#ok<DEFNU>
    % Get screen definition
    jBstFrame = bst_get('BstFrame');
    ScreenDef = bst_get('ScreenDef');
    sLayout   = bst_get('Layout');
    nbScreens = length(ScreenDef);

    % Max size for Brainstorm window
    if ~isempty(jBstFrame) && ~isempty(ScreenDef)
        % Detect on which screen was Brainstorm window at the previous session
        tol = 30;
        % If Brainstorm window is on screen 2
        if (length(ScreenDef) > 1) && ...
           (sLayout.MainWindowPos(1) >= ScreenDef(2).javaPos.getX() - tol) && (sLayout.MainWindowPos(1) <= ScreenDef(2).javaPos.getX() + ScreenDef(2).javaPos.getWidth() + tol) && ...
           (sLayout.MainWindowPos(2) >= ScreenDef(2).javaPos.getY() - tol) && (sLayout.MainWindowPos(2) <= ScreenDef(2).javaPos.getY() + ScreenDef(2).javaPos.getHeight() + tol)
            javaMax = ScreenDef(2).javaPos;
        % If Brainstorm window is on screen 1
        else
            javaMax = ScreenDef(1).javaPos;
        end
        % One screen: Half size
        if (nbScreens == 1)
            jBstFrame.setMaximumSize(java.awt.Dimension(javaMax.getWidth() / 2, javaMax.getHeight()));
            jBstFrame.setMaximizedBounds(java.awt.Rectangle(0,0,javaMax.getWidth() / 2, javaMax.getHeight()));
        % More screens: No limit
        else
            jBstFrame.setMaximumSize([]);
            jBstFrame.setMaximizedBounds([]);
        end
    end
end

author: brainstorm-tools


text:
function [s, h] = getReplicationSettings(repName)
% GETREPLICATIONSETTINGS Get settings object corresponding to previous
% replication classes. Aditionally, return a second output with the help
% text of the replication containing the information of the authors of that
% study.

switch repName

    case "bear_settings_"
        s = bear_settings_;
        fn = 'bear_settings_';
    case "61"
        s = bear_settings_61;
        fn = 'bear_settings_61';
    case "AAU2009"
        s = bear_settings_AAU2009;
        fn = 'bear_settings_AAU2009';
    case "BBE2005"
        s = bear_settings_BBE2005;
        fn = 'bear_settings_BBE2005';
    case "BvV2018"
        s = bear_settings_BvV2018;
        fn = 'bear_settings_BvV2018';
    case "CH2019"
        s = bear_settings_CH2019;
        fn = 'bear_settings_CH2019';
    case "WGP2016"
        s = bear_settings_WGP2016;
        fn = 'bear_settings_WGP2016';
    otherwise
        error('getReplicationSettings:ReplicationNotFound','Unable to find the specified replication')

end
% Change results folder to current folder
s.results_path = pwd();

% Get model definition on non-deployed setups
if ~isdeployed
    h = help(fn);
else
    h = '';
end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function sofa_pw = extrapolate_farfield_hrtfset(sofa,conf)
%EXTRAPOLATE_FARFIELD_HRTFSET far-field extrapolation of a given HRTF dataset
%
%   Usage: sofa = extrapolate_farfield_hrtfset(sofa,conf)
%
%   Input parameters:
%       sofa    - IR data set for the virtual secondary sources
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       sofa    - IR data set extra polated to conation plane wave IRs
%
%   EXTRAPOLATE_FARFIELD_HRTFSET(SOFA,conf) generates a far-field extrapolated
%   set of impulse responses, using the given irs set. Far-field means that the
%   resulting impulse responses are plane waves. The extrapolation is done via
%   WFS.
%
%   See also: get_ir, driving_function_imp_wfs
%
%   References:
%       Spors and Ahrens (2011) - "Generation of far-field head-related
%       transfer functions using sound field synthesis", 37th German Annual
%       Conference on Acoustics (DAGA), pp. 673-674,
%       http://pub.dega-akustik.de/DAGA_2011/data/articles/000370.pdf

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);
isargstruct(sofa,conf);

%% ===== Configuration ===================================================
fs = conf.fs;
dimension = conf.dimension;
showprogress = conf.showprogress;
conf.ir.usehcomp = false;

%% ===== Variables ======================================================
[nls,~,N] = size(sofa.Data.IR);
APV = SOFAcalculateAPV(sofa);
phi = rad(APV(:,1));
theta = rad(APV(:,2));
R = APV(:,3);
conf.secondary_sources.number = nls;
[conf.secondary_sources.x0(:,1), ...
 conf.secondary_sources.x0(:,2), ...
 conf.secondary_sources.x0(:,3)] = sph2cart(phi,theta,R);
conf.secondary_sources.x0(:,4:6) = ...
    direction_vector(conf.secondary_sources.x0,repmat(conf.xref,nls,1));
% Weights
if strcmp('3D',dimension)
    % Use rectangular grid to get a first approximation of the grid weights
    % R^2 * cos(theta) is the integrational weight for integration on a sphere
    conf.secondary_sources.x0(:,7) = ...
        weights_for_points_on_a_sphere_rectangle(phi,theta) .* ...
        R.^2 .* cos(theta);
else
    conf.secondary_sources.x0(:,7) = ones(nls,1);
end
% Check if we have a 2D or 3D secondary source setup
if any(phi-phi(1)>eps('single')) && any(theta-theta(1)>eps('single'))
    % 3D case
    conf.usetapwin = false;
    if ~strcmp('3D',dimension)
        warning(['You have a 3D HRTF data set, but are not using ', ...
            'conf.dimension="3D".']);
    end
else
    if strcmp('3D',dimension)
        warning(['You are using a 2D HRTF data set, but are using ', ...
            'conf.dimension="3D".']);
    end
end
if strcmp('2.5D',dimension)
    % Apply a amplitude correction, due to 2.5D. This will result in a correct
    % reproduced ILD in the resulting impulse responses, see Spors and
    % Ahrens (2011).
    amplitude_correction = -1.7 * sin(phi);
else
    amplitude_correction = zeros(nls,1);
end

%% ===== Computation =====================================================
% Get virtual secondary source positions
x0_all = secondary_source_positions(conf);
conf.wfs.hpreflow = 50;
conf.wfs.hprefhigh = aliasing_frequency(x0_all,conf);

% Initialize new irs set
sofa_pw = sofa;
sofa_pw.GLOBAL_Comment = 'Extrapolated HRTF set containing plane waves';
sofa_pw.Data.IR = zeros(nls,2,N);
sofa_pw.SourcePosition = [Inf 0 0];

% Get all HRTFs for the secondary source positions
ir_all = zeros(nls,2,N);
X = SOFAconvertCoordinates(sofa.ListenerPosition, ...
                           sofa.ListenerPosition_Type, ...
                           'cartesian');
head_orientation = SOFAconvertCoordinates(sofa.ListenerView, ...
                                          sofa.ListenerView_Type, ...
                                          'spherical');
head_orientation = rad(head_orientation(1,1:2));
warning('off','SFS:get_ir'); % Disable warning for short N
for ii=1:nls
    ir_all(ii,:,:) = get_ir(sofa,X,head_orientation, ...
                            x0_all(ii,1:3),'cartesian',conf)';
end
warning('on','SFS:get_ir');

% Generate a impulse response set for all given angles
for ii=1:nls

    % Show progress
    if showprogress, progress_bar(ii,nls); end;

    % Direction of plane wave
    [xs(1),xs(2),xs(3)] = sph2cart(phi(ii),theta(ii),R(ii));
    xs = -xs;

    % Calculate active virtual speakers
    [x0,idx] = secondary_source_selection(x0_all,xs,'pw');
    ir = ir_all(idx,:,:);
    % Apply tapering window
    x0 = secondary_source_tapering(x0,conf);

    % Get driving signals, temporarely deactivate WFS pre-filter, because it
    % will be applied once at the end
    tmp_usehpre = conf.wfs.usehpre;
    conf.wfs.usehpre = false;
    [~,delay,weight] = driving_function_imp_wfs(x0,xs,'pw',conf);
    conf.wfs.usehpre = tmp_usehpre;
    % Sum up contributions from individual virtual speakers
    for jj=1:size(x0,1)
        % Delay and weight HRTFs
        sofa_pw.Data.IR(ii,:,:) = squeeze(sofa_pw.Data.IR(ii,:,:)) + ...
            delayline(squeeze(ir(jj,:,:))',delay(jj),weight(jj),conf)';
    end
    sofa_pw.Data.IR(ii,1,:) = sofa_pw.Data.IR(ii,1,:)/10^(amplitude_correction(ii)/20);
    sofa_pw.Data.IR(ii,2,:) = sofa_pw.Data.IR(ii,2,:)/10^(-amplitude_correction(ii)/20);

end

%% ===== Pre-equalization ===============================================
sofa_pw.Data.IR(:,1,:) = wfs_preequalization(squeeze(sofa_pw.Data.IR(:,1,:))',conf)';
sofa_pw.Data.IR(:,2,:) = wfs_preequalization(squeeze(sofa_pw.Data.IR(:,2,:))',conf)';

author: sfstoolbox

repo: sfs-matlab