text:
function chebyshev1_rule ( order, a, b, filename )

%*****************************************************************************80
%
%% CHEBYSHEV1_RULE generates a Gauss-Chebyshev type 1 quadrature rule.
%
%  Discussion:
%
%    This program computes a standard Gauss-Chebyshev type 1 quadrature rule
%    and writes it to a file.
%
%    The user specifies:
%    * the ORDER (number of points) in the rule;
%    * A, the left endpoint;
%    * B, the right endpoint;
%    * FILENAME, the root name of the output files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV1_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a Gauss-Chebyshev type 1 rule for approximating\n' );
  fprintf ( 1, '    Integral ( A <= x <= B ) f(x) / sqrt ( ( x - A ) * ( B - x ) ) dx\n' );
  fprintf ( 1, '  of order ORDER.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies ORDER, A, B, and FILENAME.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER is the number of points;\n' );
  fprintf ( 1, '  A is the left endpoint;\n' );
  fprintf ( 1, '  B is the right endpoint;\n' );
  fprintf ( 1, '  FILENAME is used to generate 3 files:\n' );
  fprintf ( 1, '    filename_w.txt - the weight file\n' );
  fprintf ( 1, '    filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '    filename_r.txt - the region file.\n' );
%
%  Initialize the parameters.
%
  alpha = 0.0;
  beta = 0.0;
%
%  Get ORDER.
%
  if ( nargin < 1 )
    order = input ( '  Enter the rule order ORDER:  ' );
  elseif ( ischar ( order ) )
    order = str2num ( order );
  end
%
%  Get A.
%
  if ( nargin < 2 )
    a = input ( '  Enter the left endpoint A:  ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end
%
%  Get B.
%
  if ( nargin < 3 )
    b = input ( '  Enter the right endpoint B:  ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Get FILENAME.
%
  if ( nargin < 4 )
    fprintf ( 1,  '\n' );
    fprintf ( 1,  '  FILENAME specifies the ''root name'' of the quadrature files).\n' );
    filename = input ( '  Enter the value of FILENAME as a quoted string:  ' );
  end
%
%  Input summary.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER = %d\n', order );
  fprintf ( 1, '  A = %f\n', a );
  fprintf ( 1, '  B = %f\n', b );
  fprintf ( 1, '  FILENAME = "%s".\n', filename );
%
%  Construct the rule.
%
  kind = 2;
  alpha = 0.0;
  beta = 0.0;
  [ x, w ] = cgqf ( order, kind, alpha, beta, a, b );
%
%  Write the rule.
%
  r = [ a, b ]';
  rule_write ( order, filename, x, w, r );
%
%  Terminate.
%
  fprintf ( 1,  '\n' );
  fprintf ( 1,  'CHEBYSHEV1_RULE:\n' );
  fprintf ( 1,  '  Normal end of execution.\n' );
  fprintf ( 1,  '\n' );
  timestamp ( );

  return
end
function [ t, wts ] = cdgqf ( nt, kind, alpha, beta )

%*****************************************************************************80
%
%% CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with a classical weight function with default values for A and B,
%    and only simple knots.
%
%    There are no moments checks and no printing is done.
%
%    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%
  parchk ( kind, 2 * nt, alpha, beta );
%
%  Get the Jacobi matrix and zero-th moment.
%
  [ aj, bj, zemu ] = class_matrix ( kind, nt, alpha, beta );
%
%  Compute the knots and weights.
%
  [ t, wts ] = sgqf ( nt, aj, bj, zemu );

  return
end
function [ t, wts ] = cgqf ( nt, kind, alpha, beta, a, b )

%*****************************************************************************80
%
%% CGQF computes knots and weights of a Gauss quadrature formula.
%
%  Discussion:
%
%    The user may specify the interval (A,B).
%
%    Only simple knots are produced.
%
%    The user may request that the routine print the knots and weights,
%    and perform a moment check.
%
%    Use routine EIQFS to evaluate this quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Compute the Gauss quadrature formula for default values of A and B.
%
  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );
%
%  All knots have multiplicity = 1.
%
  mlt = zeros(nt,1);
  mlt(1:nt) = 1;
%
%  NDX(I) = I.
%
  ndx = ( 1 : nt );
%
%  Scale the quadrature rule.
%
  [ t, wts ] = scqf ( nt, t, mlt, wts, nt, ndx, kind, alpha, beta, a, b );

  return
end
function [ aj, bj, zemu ] = class_matrix ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
%
%  Discussion:
%
%    This routine computes the diagonal AJ and subdiagonal BJ
%    elements of the order M tridiagonal symmetric Jacobi matrix
%    associated with the polynomials orthogonal with respect to
%    the weight function specified by KIND.
%
%    For weight functions 1-7, M elements are defined in BJ even
%    though only M-1 are needed.  For weight function 8, BJ(M) is
%    set to zero.
%
%    The zero-th moment of the weight function is returned in ZEMU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the Jacobi matrix.
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real AJ(M), BJ(M), the diagonal and subdiagonal
%    of the Jacobi matrix.
%
%    Output, real ZEMU, the zero-th moment.
%
  temp = eps;

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < abs ( ( gamma ( temp2 ) )^2 - pi ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLASS - Fatal error!\n' );
    fprintf ( 1, '  Gamma function does not match machine parameters.\n' );
    error ( 'CLASS - Fatal error!' );
  end

  bj = zeros(m,1);
  aj = zeros(m,1);

  if ( kind == 1 )

    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod ( i, 2 );
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 2 )

    zemu = pi;

    aj(1:m) = 0.0;

    bj(1) =  sqrt ( 0.5 );
    bj(2:m) = 0.5;

  elseif ( kind == 3 )

    ab = alpha * 2.0;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 )^2 ...
      / gamma ( ab + 2.0 );

    aj(1:m) = 0.0;
    bj(1) = 1.0 / ( 2.0 * alpha + 3.0 );
    for i = 2 : m
      bj(i) = i * ( i + ab ) / ( 4.0 * ( i + alpha )^2 - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 4 )

    ab = alpha + beta;
    abi = 2.0 + ab;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
      * gamma ( beta + 1.0 ) / gamma ( abi );
    aj(1) = ( beta - alpha ) / abi;
    bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
      / ( ( abi + 1.0 ) * abi * abi );
    a2b2 = beta * beta - alpha * alpha;

    for i = 2 : m
      abi = 2.0 * i + ab;
      aj(i) = a2b2 / ( ( abi - 2.0 ) * abi );
      abi = abi^2;
      bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
        / ( ( abi - 1.0 ) * abi );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 5 )

    zemu = gamma ( alpha + 1.0 );

    for i = 1 : m
      aj(i) = 2.0 * i - 1.0 + alpha;
      bj(i) = i * ( i + alpha );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 6 )

    zemu = gamma ( ( alpha + 1.0 ) / 2.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0;
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 7 )

    ab = alpha;
    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod(i,2);
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 8 )

    ab = alpha + beta;
    zemu = gamma ( alpha + 1.0 ) * gamma ( - ( ab + 1.0 ) ) ...
      / gamma ( - beta );
    apone = alpha + 1.0;
    aba = ab * apone;
    aj(1) = - apone / ( ab + 2.0 );
    bj(1) = - aj(1) * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
    for i = 2 : m
      abti = ab + 2.0 * i;
      aj(i) = aba + 2.0 * ( ab + i ) * ( i - 1 );
      aj(i) = - aj(i) / abti / ( abti - 2.0 );
    end

    for i = 2 : m - 1
      abti = ab + 2.0 * i;
      bj(i) = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i ) ...
        / ( abti^2 ) * ( ab + i ) / ( abti + 1.0 );
    end

    bj(m) = 0.0;
    bj(1:m) =  sqrt ( bj(1:m) );

  end

  return
end
function [ d, z ] = imtqlx ( n, d, e, z )

%*****************************************************************************80
%
%% IMTQLX diagonalizes a symmetric tridiagonal matrix.
%
%  Discussion:
%
%    This routine is a slightly modified version of the EISPACK routine to
%    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
%
%    The authors thank the authors of EISPACK for permission to use this
%    routine.
%
%    It has been modified to produce the product Q' * Z, where Z is an input
%    vector and Q is the orthogonal matrix diagonalizing the input matrix.
%    The changes consist (essentialy) of applying the orthogonal transformations
%    directly to Z as they are generated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%    Roger Martin, James Wilkinson,
%    The Implicit QL Algorithm,
%    Numerische Mathematik,
%    Volume 12, Number 5, December 1968, pages 377-383.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the diagonal entries of the matrix.
%
%    Input, real E(N), the subdiagonal entries of the
%    matrix, in entries E(1) through E(N-1).
%
%    Input, real Z(N), a vector to be operated on.
%
%    Output, real D(N), the diagonal entries of the diagonalized matrix.
%
%    Output, real Z(N), the value of Q' * Z, where Q is the matrix that
%    diagonalizes the input symmetric tridiagonal matrix.
%
  itn = 30;

  prec = eps;

  if ( n == 1 )
    return
  end

  e(n) = 0.0;

  for l = 1 : n

    j = 0;

    while ( 1 )

      for m = l : n

        if ( m == n )
          break
        end

        if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) )
          break
        end

      end

      p = d(l);

      if ( m == l )
        break
      end

      if ( j == itn )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'IMTQLX - Fatal error!\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        error ( 'IMTQLX - Fatal error!' );
      end

      j = j + 1;
      g = ( d(l+1) - p ) / ( 2.0 * e(l) );
      r =  sqrt ( g * g + 1.0 );
      g = d(m) - p + e(l) / ( g + r8_sign ( g ) * abs ( r ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ii = 1 : mml

        i = m - ii;
        f = s * e(i);
        b = c * e(i);

        if ( abs ( f ) >= abs ( g ) )
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e(i+1) = f * r;
          s = 1.0 / r;
          c = c * s;
        else
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e(i+1) = g * r;
          c = 1.0 / r;
          s = s * c;
        end

        g = d(i+1) - p;
        r = ( d(i) - g ) * s + 2.0 * c * b;
        p = s * r;
        d(i+1) = g + p;
        g = c * r - b;
        f = z(i+1);
        z(i+1) = s * z(i) + c * f;
        z(i) = c * z(i) - s * f;

      end

      d(l) = d(l) - p;
      e(l) = g;
      e(m) = 0.0;

    end

  end

  for ii = 2 : n

     i = ii - 1;
     k = i;
     p = d(i);

     for j = ii : n
       if ( d(j) < p )
         k = j;
         p = d(j);
       end
     end

     if ( k ~= i )
       d(k) = d(i);
       d(i) = p;
       p = z(i);
       z(i) = z(k);
       z(k) = p;
     end

  end

  return
end
function parchk ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% PARCHK checks parameters ALPHA and BETA for classical weight functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the highest moment to
%    be calculated.  This value is only needed when KIND = 8.
%
%    Input, real ALPHA, BETA, the parameters, if required
%    by the value of KIND.
%
  if ( kind <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND <= 0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
%
  if ( 3 <= kind && alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  3 <= KIND and ALPHA <= -1.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check BETA for Jacobi.
%
  if ( kind == 4 && beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND == 4 and BETA <= -1.0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA and BETA for rational.
%
  if ( kind == 8 )
    tmp = alpha + beta + m + 1.0;
    if ( 0.0 <= tmp || tmp <= beta )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PARCHK - Fatal error!\n' );
      fprintf ( 1, '  KIND == 8 but condition on ALPHA and BETA fails.\n' );
      error ( 'PARCHK - Fatal error!' );
    end
  end

  return
end
function value = r8_sign ( x )

%*****************************************************************************80
%
%% R8_SIGN returns the sign of an R8.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( 0 <= x )
    value = +1.0;
  else
    value = -1.0;
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function rule_write ( order, filename, x, w, r )

%*****************************************************************************80
%
%% RULE_WRITE writes a quadrature rule to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, string FILENAME, specifies the output files.
%    write files 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining
%    weights, abscissas, and region.
%
%    Input, real X(ORDER), the abscissas.
%
%    Input, real W(ORDER), the weights.
%
%    Input, real R(2), the region.
%
  filename_x = strcat ( filename, '_x.txt' );
  filename_w = strcat ( filename, '_w.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1, order, w' );
  r8mat_write ( filename_x, 1, order, x' );
  r8mat_write ( filename_r, 1, 2,     r' );

  return
end
function [ t, wts ] = scqf ( nt, t, mlt, wts, nwts, ndx, kind, alpha, ...
  beta, a, b )

%*****************************************************************************80
%
%% SCQF scales a quadrature formula to a nonstandard interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real T(NT), the original knots.
%
%    Input, integer MLT(NT), the multiplicity of the knots.
%
%    Input, real WTS(NWTS), the weights.
%
%    Input, integer NWTS, the number of weights.
%
%    Input, integer NDX(NT), used to index the array WTS.
%    For more details see the comments in CAWIQ.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the scaled knots.
%
%    Output, real WTS(NWTS), the scaled weights.
%
  temp = eps;

  parchk ( kind, 1, alpha, beta )

  if ( kind == 1 )

    al = 0.0;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 2 )

    al = -0.5;
    be = -0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 3 )

    al = alpha;
    be = alpha;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 4 )

    al = alpha;
    be = beta;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 5 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / b;
    al = alpha;
    be = 0.0;

  elseif ( kind == 6 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / sqrt ( b );
    al = alpha;
    be = 0.0;

  elseif ( kind == 7 )

    al = alpha;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 8 )

    if ( a + b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  A + B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = a + b;
    al = alpha;
    be = beta;

  elseif ( kind == 9 )

    al = 0.5;
    be = 0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  end

  p = slp^( al + be + 1.0 );

  for k = 1 : nt

    t(k) = shft + slp * t(k);
    l = abs ( ndx(k) );

    if ( l ~= 0 )
      tmp = p;
      for i = l : l + mlt(k) - 1
        wts(i) = wts(i) * tmp;
        tmp = tmp * slp;
      end
    end

  end

  return
end
function [ t, wts ] = sgqf ( nt, aj, bj, zemu )

%*****************************************************************************80
%
%% SGQF computes knots and weights of a Gauss Quadrature formula.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with simple knots from the Jacobi matrix and the zero-th
%    moment of the weight function, using the Golub-Welsch technique.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real AJ(NT), the diagonal of the Jacobi matrix.
%
%    Input, real BJ(NT), the subdiagonal of the Jacobi
%    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
%
%    Input, real ZEMU, the zero-th moment of the weight function.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Exit if the zero-th moment is not positive.
%
  if ( zemu <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SGQF - Fatal error!\n' );
    fprintf ( 1, '  ZEMU <= 0.\n' );
    error ( 'SGQF - Fatal error!' );
  end
%
%  Set up vectors for IMTQLX.
%
  wts = zeros ( nt, 1 );

  wts(1) = sqrt ( zemu );
  wts(2:nt) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt) = wts(1:nt).^2;

  return
end
function timestamp ( )

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

  return
end


text:
function [ v, more ] = simplex_lattice_layer_point_next ( n, c, v, more )

%*****************************************************************************80
%
%% SIMPLEX_LATTICE_LAYER_POINT_NEXT: next simplex lattice layer point.
%
%  Discussion:
%
%    The simplex lattice layer L is bounded by the lines
%
%      0 <= X(1:N),
%      L - 1 < sum X(1:N) / C(1:N)  <= L.
%
%    In particular, layer L = 0 always contains just the origin.
%
%    This function returns, one at a time, the points that lie within
%    a given simplex lattice layer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, integer C(N+1), coefficients defining the
%    lattice layer in entries 1 to N, and the laver index in C(N+1).
%    The coefficients should be positive, and C(N+1) must be nonnegative.
%
%    Input/output, integer V(N).  On first call for a given layer,
%    the input value of V is not important.  On a repeated call for the same
%    layer, the input value of V should be the output value from the previous
%    call.  On output, V contains the next lattice layer point.
%
%    Input/output, logical MORE.  On input, set MORE to FALSE to indicate
%    that this is the first call for a given layer.  Thereafter, the input
%    value should be the output value from the previous call.  On output,
%    MORE is TRUE if the returned value V is a new point.
%    If the output value is FALSE, then no more points were found,
%    and V was reset to 0, and the lattice layer has been exhausted.
%

%
%  Treat layer C(N+1) = 0 specially.
%
  if ( c(n+1) == 0 )
    if ( ~more )
      v(1:n) = 0;
      more = 1;
    else
      more = 0;
    end
    return
  end
%
%  Compute the first point.
%
  if ( ~more )

    v(1) = ( c(n+1) - 1 ) * c(1) + 1;
    v(2:n) = 0;
    more = 1;

  else

    c1n = i4vec_lcm ( n, c );

    rhs1 = c1n * ( c(n+1) - 1 );
    rhs2 = c1n *   c(n+1);
%
%  Try to increment component I.
%
    for i = 1 : n

      v(i) = v(i) + 1;

      v(1:i-1) = 0;

      if ( 1 < i )
        v(1) = rhs1;
        for j = 2 : n
          v(1) = v(1) - ( c1n / c(j) ) * v(j);
        end
        v(1) = floor ( ( c(1) * v(1) ) / c1n );
        v(1) = max ( v(1), 0 );
      end

      lhs = 0;
      for j = 1 : n
        lhs = lhs + ( c1n / c(j) ) * v(j);
      end

      if ( lhs <= rhs1 )
        v(1) = v(1) + 1;
        lhs = lhs + c1n / c(1);
      end

      if ( lhs <= rhs2 )
        return
      end

    end

    v(1:n) = 0;
    more = 0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function filum_test14 ( )

%*****************************************************************************80
%
%% TEST14 tests FILENAME_INC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
  ntest = 4;

  string(1,1:11) = 'file???.dat';
  string(2,1:11) = 'file072.dat';
  string(3,1:11) = '2cat9.dat  ';
  string(4,1:11) = 'fred98.txt ';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST14\n' );
  fprintf ( 1, '  FILENAME_INC increments a string\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     Input             Output\n' );

  for i = 1 : ntest
    fprintf ( 1, '\n' );
    filename = string(i,:);
    for j = 1 : 4
      filename_old = filename;
      filename = filename_inc ( filename );
      fprintf ( 1, '  %20s  %20s\n', filename_old, filename );
      if ( s_len_trim ( filename ) == 0 )
        fprintf ( 1, '  (Empty output string.  Quit loop!)\n' );
        break
      end
    end
  end

  return
end


text:
function y = pmean(x,p)
%Calculates generalized mean
% x = vector subject to generalized mean calculations. Can be also matrix.
% p = parameter for generalized mean.
% y = generalized mean value.
y = (sum(x.^p,1)/size(x,1)).^(1/p);

author: Sable

repo: mcbench-benchmarks


text:
function [gps_week, gps_sow, gps_dow] = date2gps(date)

% SYNTAX:
%   [gps_week, gps_sow, gps_dow] = date2gps(date);
%
% INPUT:
%   date = date [year, month, day, hour, min, sec]
%
% OUTPUT:
%   gps_week = GPS week
%   gps_sow  = GPS seconds of week
%   gps_dow  = GPS day of week
%
% DESCRIPTION:
%   Conversion from calendar date to GPS time.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

gps_start_datenum = 723186; %This is datenum([1980,1,6,0,0,0])

%number of days since the beginning of GPS time
%deltat   = (datenum([date(:,1), date(:,2), date(:,3)]) - gps_start_datenum);
% hack: datenummmx is faster cause it does not check argins
deltat   = (datenummx(date(:,1:3)) - gps_start_datenum);

gps_week = floor(deltat/7);            %GPS week
gps_dow  = floor(deltat - gps_week*7); %GPS day of week
gps_sow  = (deltat - gps_week*7)*86400;
gps_sow = gps_sow + date(:,4)*3600 + date(:,5)*60 + date(:,6); %GPS seconds of week

% %alternative way, using the Julian day
% jd = date2jd(date);
% [gps_week, gps_dow, gps_sow] = jd2gps(jd);

author: goGPS-Project


text:
function [MB] = bit2MB(bit)
% Convert computery things from bits to megabytes.
% Chad A. Greene 2012
MB = bit*2^-23;


text:
%Ray3D Ray in 3D space
%
% This object represents a ray in 3D space, defined by a point on the ray
% and a direction unit-vector.
%
% Methods::
% intersect    Intersection of ray with plane or ray
% closest      Closest distance between point and ray
% char         Ray parameters as human readable string
% display      Display ray parameters in human readable form
%
% Properties::
% P0    A point on the ray (3x1)
% d     Direction of the ray, unit vector (3x1)
%
% Notes::
% - Ray3D objects can be used in vectors and arrays

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

    methods
        function r = Ray3D(P0, d)
        %Ray3D.Ray3D Ray constructor
        %
        % R = Ray3D(P0, D) is a new Ray3D object defined by a point on the ray P0
        % and a direction vector D.
            r.P0 = P0(:);
            r.d = unit(d(:));
        end

        function [x,e] = intersect(r1, r2)
        %Ray3D.intersect Intersetion of ray with line or plane
        %
        % X = R.intersect(R2) is the point on R that is closest to the ray R2.
        % If R is a vector then then X has multiple columns, corresponding to
        % the intersection of R(i) with R2.
        %
        % [X,E] = R.intersect(R2) as above but also returns the closest distance
        % between the rays.
        %
        % X = R.intersect(P) returns the point of intersection between the
        % ray R and the plane P=(a,b,c,d) where aX + bY + cZ + d = 0.
        % If R is a vector then X has multiple columns, corresponding to
        % the intersection of R(i) with P.

            if isa(r2, 'Ray3D')
                % ray intersect ray case
                if length(r1) ~= length(r2)
                    error('can only intersect rays pairwise');
                end

                for i=1:length(r1)
                    alpha = [-r1(i).d r2(i).d cross(r1(i).d,r2(i).d)] \ ...
                        (r1(i).P0-r2(i).P0);
                    x(:,i) = r1(i).P0 + alpha(1)*r1(i).d;
                    e(i) = norm(r1(i).P0 - r2(i).P0 + ...
                        alpha(1)*r1(i).d -alpha(2)*r2(i).d);
                end
            else
                % ray intersect plane case
                % plane is P = (a,b,c,d), P.x = 0

                for i=1:length(r1)
                    n = r2(1:3); d = r2(4);
                    alpha = -(d + dot(n, r1.P0)) / ( dot(n, r1.d) );
                    x(:,i) = r1.P0 + alpha*r1.d;
                end
            end
        end

        % closest distance between point and line
        function [x,e] = closest(r1, P)
        %Ray3D.closest Closest distance between point and ray
        %
        % X = R.closest(P) is the point on the ray R closest to the point P.
        %
        % [X,E] = R.closest(P) as above but also returns the distance E between
        % X and P.
            alpha = dot(P - r.P0, r.d);
            x = r1.P0 + alpha*r1.d;
            e = norm(x-P);
        end

        function display(rays)
        %Ray3D.display Display value
        %
        % R.display() displays a compact human-readable representation of the Ray3D's
        % value.  If R is a vector then the elements are printed one per line.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a Ray3D object and the command has no trailing
        %   semicolon.
        %
        % See also Ray3D.char.
            disp(' ');
            disp([inputname(1), ' = '])
            disp(' ');
            if length(rays) > 20
                fprintf('%d corresponding points (listing suppressed)\n', length(rays));
            else
                disp( char(rays) );
            end
        end % display()

        function s = char(rays)
        %Ray3D.char Convert to string
        %
        % S = R.char() is a compact string representation of the Ray3D's value.
        % If R is a vector then the string has multiple lines, one per element.

            s = '';
            for r=rays
                ss = sprintf('d=(%g, %g, %g), P0=(%g, %g, %g)\n', ...
                    r.d, r.P0);
                s = strvcat(s, ss);
            end
        end

    end
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function ch=rayleigh(n_samp,fd,fs)

% ch=rayleigh(n_samp,fd,fs)
%
% Create n_samp samples of a Rayleigh channel whose doppler spread is fd
% and which is sampled at frequency fs.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

error(nargchk(3,3,nargin));
error(chk_param(n_samp,'n_samp','scalar','real','integer','>=',1));
error(chk_param(fd,'fd','scalar','real','>=',0));
error(chk_param(fs,'fs','scalar','real','>',0));
if (fd>fs/2)
  error('Doppler spread must be smaller than fs/2');
end

% Special case when requested doppler spread is zero.
if (fd==0)
  ch=rayleigh(1,0.5,2);
  ch=repmat(ch,1,n_samp);
  return;
end

% First, a rayleigh channel is created with a doppler frequency of 0.5
% and a sampling frequency of 2. Then, this channel is interpolated
% up to the desired sampling frequency.

% Calculate the interpolation factor.
intp=fs/fd/4;

% Rayleigh transfer function
rtf=@(f)(abs(f)<.5).*sqrt(1./((pi*0.5)*sqrt(1-(f/(0.5+eps(0.5))).^2)));

% Create the source Rayleigh signal
n_samp_source=max(2048,ceil(n_samp/intp));
bb=sqrt(2)*cyc_filt(blnoise(n_samp_source),rtf);

% Now, interpolate up to the desired sampling rate.
ch=interpft(bb,round(n_samp_source*intp));
ch=ch(1:n_samp);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function geometry_test198 ( )

%*****************************************************************************80
%
%% TEST198 tests SHAPE_POINT_NEAR_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  nside = 6;
  ntest = 8;
  center(1:2,1) = [ 3.0; 0.0 ];
  p1(1:2,1) = [ 5.0; 0.0 ];
  ptest(1:2,1:ntest) = [ ...
     3.0, 0.0; ...
     5.0, 0.0; ...
     4.0, 0.0; ...
    10.0, 0.0; ...
     4.0, 1.7320508; ...
     5.0, 2.0 * 1.7320508; ...
     3.0, 1.7320508; ...
     3.0, 1.7320508 / 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST198\n' );
  fprintf ( 1, '  For a shape in 2D,\n' );
  fprintf ( 1, '  SHAPE_POINT_NEAR_2D computes the nearest\n' );
  fprintf ( 1, '    point to a point;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sides:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %d\n', nside );

  r8vec_print ( 2, center, '  Hexagon center:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Hexagon vertex #1:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %8f  %8f\n', p1(1:2,1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I       X            Y              PN     Dist\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    p(1:2,1) = ptest(1:2,i);

    [ pn, dist ] = shape_point_near_2d ( center, p1, nside, p ) ;

    fprintf ( 1, '  %6d  %10f  %10f  %10f  %10f  %10f\n', ...
      i, p(1:2,1), pn(1:2,1), dist );

  end

  return
end


text:
function test_failed=test_idft
Lr=[1, 19, 20];

test_failed=0;

disp(' ===============  TEST_IDFT ==============');

for jj=1:length(Lr)
  L=Lr(jj);
    for n = 1:2

    if (n==1)
       type = 'complex';
       c=tester_crand(L,1);
    elseif (n==2)
       type = 'real';
       c=tester_rand(L,1);
    end

    f1=idft(c);
    f2=ref_idft(c);

    res=norm(f1-f2);
    [test_failed,fail]=ltfatdiditfail(res,test_failed);
    s=sprintf('IDFT %6s  L:%3i %0.5g %s',type,L,res,fail);
    disp(s);
    end
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%-----------------------------------------------------------------------
% Job saved on 01-Apr-2016 11:17:08 by cfg_util (rev $Rev: 6134 $)
% spm SPM - SPM12 (6225)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{2}.spm.stats.fmri_est.spmmat = cfg_dep('fMRI model specification: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
matlabbatch{2}.spm.stats.fmri_est.write_residuals = 0;
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.space.volume.block_type = 'Slices';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.signal = 'Uninformative';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.ARP = 3;
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.noise.UGL = 1;
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.LogEv = 'Yes';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.anova.first = 'No';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.anova.second = 'Yes';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.gcon(1).name = 'main positive';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.gcon(1).convec = 1;
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.gcon(2).name = 'main negative';
matlabbatch{2}.spm.stats.fmri_est.method.Bayesian.gcon(2).convec = -1;

author: translationalneuromodeling

repo: tapas


text:
% hierarchical bayesian mixture model
function [MCMC_State] = grp_hbmm_mcmc(varargin)

% inputs:
%   B: {subj}(M_i x M_i x Q) matrix of (possibly time-varying)  (obtain from EEG(i).CAT.Conn)
%      connectivities or basis coefficients
%   S: {subj}(M_i x 3) matrix of dipole locations (obtain from EEG(i).dipfit)
%   M_i: [N x 1] vector of number of sources for each subject
%   M:  number of clusters required
% automatically determined:
%   N:  number of subjects (length(B))
%   Q:  number of basis coefficients (size(B{1},3))
%   k:  clustering mode (1 = kmeans, 2=gmm)

arg_define([0 3], varargin, ...
    arg_norep({'B','Connectivity'},mandatory,[],'Connectivity matrices. B{i} is an (M_i x M_i x Q) matrix of (possibly time-varying) connectivity values or basis coefficients for the ith subject. M_i is the number of components for the ith subject.'), ...
    arg_norep({'S','DipoleLocations'},mandatory,[],'Dipole locations. S{i} is an (M_i x 3) matrix of [X Y Z] dipole locations for the ith subject. M_i is the number of components for the ith subject.'), ...
    arg_norep({'MCMC_InitState'},struct([]),[],sprintf('Initial state of MCMC sampler. This can be computed from grp_hbmm_initMCMC or can be the output of a previous call to grp_hmbb_mcmc (i.e. the last state of the MCMC iterator). This structure must contain the following fields:\nZ\t: Mi x M matrices of group indicators\nS_BAR\t: cluster centroid locations\nSIGMA_S\t: cluster centroid variances\nB_BAR\t: group level connectivities\nSIGMA_B\t: variances of connectivities\nN_k\t: number of components that belong to each cluster\nN_k1k2\t: pairwise counts of group level clusters'),'type','struct'), ...
    arg_sub({'hyperparams','Hyperparameters'},[],...
    {...
        arg({'c','ConnMeanPriorVar'},10000,[eps Inf],'Between-cluster mean connectivity prior variance. Variance of gaussian prior for between-cluster mean connectivity. Larger --> more uncertainty'), ...
        arg({'eta','DipoleLocVarPriorShape'},[],[],'Within-cluster dipole location prior variance D.O.F. Degrees of freedom for inverse-wishart prior distribution for within-cluster dipole location covariance matrices. Leave empty to compute from initial clustering.'), ...
        arg({'SS','DipoleLocVarPriorScale'},[],[],'Within-cluster dipole location prior variance scale. Scale matrix of inverse-wishart prior distribution for within-cluster dipole location covariance matrices. Leave empty to compute from initial clustering.'), ...
        arg({'p1','ConnVarPriorShape','a'},1,[],'Between-cluster mean connectivity variance shape.'), ...
        arg({'p2','ConnVarPriorScale','b'},1,[],'Between-cluster mean connectivity variance scale.'), ...
    },'MCMC hyperparameters'), ...
    arg({'nMCMCiters','NumMCMCIters','niter','niters'},1000, [1 Inf], 'Number of MCMC iterations for spline fitting'), ...
    arg({'burnInFraction','BurnInFractionForMCMC'},0.5,[0 0.99],'Fraction of initial MCMC samples to discard (burn in period). The number of MCMC samples is taken to be the smaller of MCMCitersOffDiag and MCMCitersDiag.'), ...
    arg({'thinFactor','ThinningFactor'},1,[1 Inf],'Thinning factor for MCMC. We keep every kth MCMC sample, where k=ThinningFactor. This is useful when we have limited available memory to store MCMC results since successive MCMC estimates are more likely to be correlated.'), ...
    arg({'normlog','NormalizeAndLogTransform'},true,[],'Transform data before smoothinFactorg. Normalize across last dim (time), add 1 and take logarithm. Inverse transform is applied after smoothinFactorg'), ...
    arg({'verb','VerbosityLevel'},2,{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical'), ...
    arg({'appendLastState','AppendLastState'},false,[],'Append new state to initial MCMC state') ...
    );

%     arg({'basisCoeffVarPrior'},1000,[eps Inf],'Variance of basis coefficient gaussian prior. Larger --> more wiggling allowed'), ...
%     arg({'noiseVarPriorShape'},0.01,[eps Inf],'Shape (D.O.F) of noise variance prior. This is the "alpha" parameter of the inverse gamma prior distribution. Increasing noiseVarPriorShape --> decreased variance of noise variance distribution.'), ...
%     arg({'noiseVarPriorScale'},0.01,[eps Inf],'Scale parameter of noise variance prior. This is the "theta" (1/beta) parameter of inverse gamma prior distribution. Increasing noiseVarPriorScale --> right-shift of distribution --> (increase in expected noise variance). In general MEAN(noiseVariance) = noiseVarPriorScale/noiseVarPriorShape and MODE(noiseVariance) = noiseVarPriorScale/(noiseVarPriorShape-1) for noiseVarPriorShape>=1.'), ...
%     arg({'initNoiseVariance'},0.1,[eps Inf],'Initial noise variance'), ...

% set up initial state of MCMC
if isempty(MCMC_InitState)
    error('You must provide an initial state for the MCMC iterator (MCMC_InitState)');
end

% initialize vars based on last state of MCMC_InitState
% this will initialize the following variables:
% 'Z','S_BAR','SIGMA_S','B_BAR','SIGMA_B','N_k','N_k1k2'
varnames = setdiff_bc(fieldnames(MCMC_InitState),'initstate');
if MCMC_InitState.initstate
    % initialize to initial state of MCMC iterator
    for i = 1:length(varnames)
        eval(sprintf('%s=MCMC_InitState.(''%s'');',varnames{i},varnames{i}));
    end
else
    % initialize to last state of MCMC iterator
    for i = 1:length(varnames)
        eval(sprintf('%s=MCMC_InitState.(''%s''){end};',varnames{i},varnames{i}));
    end
end
if ~appendLastState
    clear MCMC_InitState;
end

% define some vars
M   = length(N_k);                      % number of clusters
N   = length(B);                        % number of subjects
M_i = cellfun(@(B_i) size(B_i,1),B);    % number of sources for each subject
Q   = size(B{1},3);                     % connectivity time-series dimension

% initialize prior probabilities of cluster membership
if ~exist('MU','var')
    MU  = ones(M,1)/M; end
% initialize hyperparams
if isempty(hyperparams.eta)
    hyperparams.eta=median(N_k); end
if isempty(hyperparams.SS)
    hyperparams.SS = hyperparams.eta*mean(SIGMA_S,3); end

% compute number of burn-in samples
numBurnInSamples = floor(burnInFraction*nMCMCiters);
niterToKeep      = round((nMCMCiters-numBurnInSamples)/thinFactor);
if verb,
    fprintf(['I will discard %d burn-in samples.\n' ...
             'I will thin the distribution by a factor of %d samples\n', ...
             'The distribution of the estimator will have %d samples\n'], ...
            numBurnInSamples,thinFactor,niterToKeep);
end

%% initialize arrays of posterior draws
Z_array         = cell(niterToKeep,1);
S_BAR_array     = cell(niterToKeep,1);
SIGMA_S_array   = cell(niterToKeep,1);
B_BAR_array     = cell(niterToKeep,1);
SIGMA_B_array   = cell(niterToKeep,1);
MU_array        = cell(niterToKeep,1);
N_k_array       = cell(niterToKeep,1);
N_k1k2_array    = cell(niterToKeep,1);
iter_array      = 0;

%% run MCMC
for iter=1:nMCMCiters

    % Draw S_BAR (group centroid locations)
    SSinv = double(inverse(hyperparams.SS));
    for k=1:M
        mu_s_bar_k=zeros(3,1);
        Sigma_s_bar_k=double(inverse(N_k(k)*double(inverse(SIGMA_S(:,:,k)))+SSinv));
        for i=1:N
            if max(Z{i}(:,k))==1
                j_k=find(Z{i}(:,k)==1);
                for j=1:length(j_k)
                    mu_s_bar_k=mu_s_bar_k+S{i}(j_k(j),:)';
                end
            end
        end
        mu_s_bar_k=Sigma_s_bar_k/SIGMA_S(:,:,k)*mu_s_bar_k;
        S_BAR(k,:)=mvnrnd(mu_s_bar_k,Sigma_s_bar_k);
    end

    % Draw SIGMA_S (group centroid covariance matrices)
    for k=1:M
        eta_k=hyperparams.eta+.5*N_k(k);
        SS_k=hyperparams.SS;
        for i=1:N
            if max(Z{i}(:,k))==1
                j_k=find(Z{i}(:,k)==1);
                for j=1:length(j_k)
                    SS_k=SS_k+.5*(S{i}(j_k(j),:)-S_BAR(k,:))'*(S{i}(j_k(j),:)-S_BAR(k,:));
                end
            end
        end
        SS_k=.5*(SS_k+SS_k');
        inv_SS_k=double(inverse(SS_k));
        inv_SS_k=.5*(inv_SS_k+inv_SS_k');
        SIGMA_S(:,:,k)=double(inverse(wishrnd(inv_SS_k,eta_k)));
    end

    % Draw B_BAR (group mean connectivites)
    for k1=1:M
        for k2=1:M
            mu_b_bar_k1k2=zeros(Q,1);
            Sigma_sq_b_k1k2=1/(1/hyperparams.c+N_k1k2(k1,k2)/SIGMA_B(k1,k2))*eye(Q);
            for i=1:N
                n_k1=sum(Z{i}(:,k1));
                n_k2=sum(Z{i}(:,k2));
                if(and(n_k1>=1,n_k2>=1))
                    j_k1=find(Z{i}(:,k1)==1);
                    j_k2=find(Z{i}(:,k2)==1);
                    for j1=1:n_k1
                        for j2=1:n_k2
                            mu_b_bar_k1k2=mu_b_bar_k1k2+squish(B{i}(j_k1(j1),j_k2(j2),:));
                        end
                    end
                end
            end
            mu_b_bar_k1k2=Sigma_sq_b_k1k2*mu_b_bar_k1k2/SIGMA_B(k1,k2);
            B_BAR(k1,k2,:)=mvnrnd(mu_b_bar_k1k2,Sigma_sq_b_k1k2);
        end
    end

    % Draw SIGMA_B (group connectivity variances)
    for k1=1:M
        for k2=1:M
            p1_k1k2=hyperparams.p1+.5*Q*N_k1k2(k1,k2);
            p2_k1k2=hyperparams.p2;
            for i=1:N
                n_k1=sum(Z{i}(:,k1));
                n_k2=sum(Z{i}(:,k2));
                if(and(n_k1>=1,n_k2>=1))
                    j_k1=find(Z{i}(:,k1)==1);
                    j_k2=find(Z{i}(:,k2)==1);
                    for j1=1:n_k1
                        for j2=1:n_k2
                            p2_k1k2=p2_k1k2+.5*sum((B{i}(j_k1(j1),j_k2(j2),:)-B_BAR(k1,k2,:)).^2);
                        end
                    end
                end
            end
            SIGMA_B(k1,k2)=1/gamrnd(p1_k1k2,1/p2_k1k2);
        end
    end

    % Draw Z (indicators of group membership)
    ind=randperm(N);
    for i=ind
        ind_i=randperm(M_i(i));
        for j=ind_i
            Log_p_ij=zeros(M,1);
            for k=1:M
                Sigma_s_ijk=SIGMA_S(:,:,k);
                s_ijk=S{i}(j,:)-S_BAR(k,:);
                log_p_s_ijk=-.5*log(det(Sigma_s_ijk))-.5*s_ijk/Sigma_s_ijk*s_ijk';
                log_p_b_ijk=0;
                for j1=1:M_i(i)
                    if j==j1
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(k,k)))...
                        -.5*sum((B{i}(j,j,:)-B_BAR(k,k,:)).^2)/SIGMA_B(k,k);
                    end
                    if and(not(j==j1),find(Z{i}(j1,:)==1)==k)
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(k,k)))...
                        -.5*sum((B{i}(j1,j1,:)-B_BAR(Z{i}(j1,:)==1,Z{i}(j1,:)==1,:)).^2)/...
                            SIGMA_B(Z{i}(j,:)==1,Z{i}(j,:)==1);
                    end
                    if not(find(Z{i}(j,:)==1)==find(Z{i}(j1,:)==1))
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(k,Z{i}(j1,:)==1)))...
                        -.5*sum((B{i}(j,j1,:)-B_BAR(k,Z{i}(j1,:)==1,:)).^2)/...
                              SIGMA_B(k,Z{i}(j1,:)==1);
                    end
                end
                for j2=1:M_i(i)
                    if not(find(Z{i}(j2,:)==1)==find(Z{i}(j,:)==1))
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(Z{i}(j2,:)==1,k)))...
                        -.5*sum((B{i}(j2,j,:)-B_BAR(Z{i}(j2,:)==1,k,:)).^2)/...
                                SIGMA_B(Z{i}(j2,:)==1,k);
                    end
                end
                Log_p_ij(k)=MU(k)+log_p_s_ijk+log_p_b_ijk;
            end
            Log_p_ij=Log_p_ij-max(Log_p_ij);
            P_ij=exp(Log_p_ij)/sum(exp(Log_p_ij));
            v=rand;
            Z{i}(j,:)=0*Z{i}(j,:);
            if v<=P_ij(1)
                   Z{i}(j,1)=1;
                   log_p_ij=Log_p_ij(1);
            end
            for k=1:(M-1)
               if(and(v>sum(P_ij(1:k)),v<=sum(P_ij(1:(k+1)))))
                   Z{i}(j,k+1)=1;
                   log_p_ij=Log_p_ij(k);
               end
            end
        end
    end

    % Draw MU (probabilities of clusters)
    N_k1k2=zeros(M);
    for k1=1:M
        for k2=1:M
            for i=1:N
                n_k1=sum(Z{i}(:,k1));
                n_k2=sum(Z{i}(:,k2));
                if(and(n_k1>=1,n_k2>=1))
                    if not(k1==k2)
                        N_k1k2(k1,k2)=N_k1k2(k1,k2)+n_k1*n_k2;
                    end
                    if k1==k2
                        N_k1k2(k1,k2)=N_k1k2(k1,k2)+n_k1;
                    end
                end
            end
        end
    end
    N_k=diag(N_k1k2);
    MU=drchrnd(1/M*ones(1,M)+N_k',1);

    % Posterior Draw arrays
    if and(iter>=numBurnInSamples,thinFactor*round(iter/thinFactor)==iter)
        iter_array=iter_array+1;
        Z_array{iter_array}=Z;
        S_BAR_array{iter_array}=S_BAR;
        SIGMA_S_array{iter_array}=SIGMA_S;
        B_BAR_array{iter_array}=B_BAR;
        SIGMA_B_array{iter_array}=SIGMA_B;
        N_k_array{iter_array}=N_k;
        N_k1k2_array{iter_array}= N_k1k2;
        MU_array{iter_array}=MU;

        if verb
            disp(['iter: ' num2str(iter)]);
%             % compute and print out intermediate estimates of source locations
%             S_BAR_mean=S_BAR_array{1};
%             for j=2:iter_array
%                 S_BAR_mean=S_BAR_mean+S_BAR_array{j};
%             end
%             S_BAR_mean=S_BAR_mean/iter_array;
%             disp([S_BAR zeros(M,1) S_BAR_mean])
        end

        %[sqrt(sum(min(pdist2(S_BAR_true,S_BAR_st)').^2))/M sqrt(sum(min(pdist2(S_BAR_true,S_BAR_mean)').^2))/M]

        if verb
%             %compute and print intermediate indicator probs for random subject
%             i=discretesample(ones(N,1)/N,1);
%             Z_i_mean=Z_array{1}{i};
%             for j=2:iter_array
%                Z_i_mean=Z_i_mean+Z_array{j}{i};
%             end
%             Z_i_mean=Z_i_mean/iter_array;
% %         Z_st{i}
%            disp(Z_i_mean);
        end
    end
end

% return MCMC state
varnames = {'Z_array','S_BAR_array','SIGMA_S_array','B_BAR_array','SIGMA_B_array','MU_array','N_k1k2_array','N_k_array'};
for i = 1:length(varnames)
    % state fieldnames have trailing '_array' stripped
    fname = strrep(varnames{i},'_array','');
    vname = varnames{i};
    if appendLastState
        % append current state estimates to last state estimate
        MCMC_State.(fname) = [MCMC_State.(fname); eval(vname)];
    else
        % replace last state estimate with current state estimates
        MCMC_State.(fname) = eval(vname);
    end
end
MCMC_State.initstate = false;

function r = drchrnd(a,n)
    p = length(a);
    r = gamrnd(repmat(a,n,1),1,n,p);
    r = r ./ repmat(sum(r,2),1,p);

repo: REST


text:
function f = dejongsfcn(x)

if strcmpi(x,'init')
    f.Vectorized = 'on' ;
    f.PopInitRange = [-5; 5] ;
    f.KnownMin = [0 0] ; % For plotting only
else
    f = sum(x.*x,2) ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function cmp = cstrcmp(a, b)
%CSTRCMP  Compare two character arrays according to ASCII dictionary order in C
%style.
%
%   INPUTS:
%
%   -|a|, |b|: character arrays, not necessarily of the same length
%
%   OUTPUT:
%
%   -|cmp|: integer holding the result of the comparison. Equal to 0 when the
%    two inputs are equal, positive when |a| is greater than |b|, and negative
%    when |a| is less than |b|.

% Force the arrays to equal length.
x = char({a;b});

% Subtract one from the other.
d = x(1,:) - x(2,:);

if any(d)
    cmp = d(find(d, 1));
else
    cmp = 0;
end

end

author: sakaridis


text:
function Gnew = gsp_kron_reduction( G,ind )
%GSP_KRON_REDUCTION See gsp_kron_reduce

Gnew = gsp_kron_reduce(G,ind);

end

repo: gspbox


text:
function [path, fileList, matches] = DPMpcNEW(leak, planC, stateS, nhist, OutputError)
% written by KZ
% Calculate IM struct, call DPM, get
%
% Usage: IIMCalc = calculateBeamIM(0.018, 'DPM_10beams_planC', 1, 1, 5)
%   leak = 0.018
%   planC_File = file name contained planC, eg. 'DPM_10beams_planC'
%   indexBeam = Beam Index: 1, 2, ...
%   imin = Beamlet index
%   imax = Beamlet index
%
%   eval(['load Beam',num2str(indexBeam),'_w']);
%   where, Beam1_w.mat (Beam2_w.mat) need to be in the current directory. It's the weight of
%   the beamlets for each beam.
%
% JC. Aug 10, 2005, Need to make planC read from the disk, since in the
% stand-alone mode, no matlab run is expected.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% open planC_File
% In matlab release 14, "eval" can be compiled.
Error = []; %Initialize Error to NULL.
% load(planC_File);
%eval(['load ',planC_File]);
%leak = 0.018; %leakage in per *100
% leak = str2num(leak)
currentDir = cd;

[FileName,path] = uigetfile('*.*','Select DICOM Dose File, RD');

if path == 0
    errordlg('DICOM Dose File Should exist');
    error('DICOM Dose File Should exist');
end

cd(path);

fileList = dir;

filesToRun = {};

matches = strmatch('RD', {fileList.name});

cd(currentDir);

for i = 1:length(matches)
    [maxDose, dA] = getDICOMMaxDose([path fileList(matches(i)).name]);
    dM(i) = maxDose;
end

for indexBeam = 1 : planC{7}.FractionGroupSequence.Item_1.NumberOfBeams

    bs = planC{7}.BeamSequence.(['Item_' num2str(indexBeam)]);

    LS = getDICOMLeafPositions(bs);

    [inflMap, xV, yV, colDividerXCoord, rowDividerYCoord, rowLeafPositions] = getLSInfluenceMapFactor(LS,leak,bs.BeamNumber);

    gA = bs.ControlPointSequence.Item_1.GantryAngle;
    iC = bs.ControlPointSequence.Item_1.IsocenterPosition;

    IMBase.beams(1).gantryAngle = gA;
    IMBase.beams(1).isocenter.x = iC(1)/10;
    IMBase.beams(1).isocenter.y = -iC(2)/10;
    IMBase.beams(1).isocenter.z = -iC(3)/10;

    IMBase.beams(1).isodistance = bs.SourceAxisDistance/10;

    IMBase.beams(1).xRel = IMBase.beams(1).isodistance * sindeg(IMBase.beams(1).gantryAngle);

    IMBase.beams(1).yRel = IMBase.beams(1).isodistance * cosdeg(IMBase.beams(1).gantryAngle);

    IMBase.beams(1).beamEnergy = bs.ControlPointSequence.Item_1.NominalBeamEnergy;

    [PBX] = getPBXCoorFromFluence(inflMap, rowLeafPositions, 2, 50, colDividerXCoord, rowDividerYCoord, 10);

    [xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, w_field] = getPBVectors(PBX,rowLeafPositions,colDividerXCoord,rowDividerYCoord,inflMap);

    %Visualization
%         figure;hAxis1 = axes;imagesc([xV(1)+.05 xV(end)-.05], [yV(1)+.05, yV(end)-.05], inflMap);
%         set(gcf, 'renderer', 'zbuffer')
%         figure;hAxis2 = axes;hold on;
%         set(gcf, 'renderer', 'zbuffer')
%         xL = get(hAxis1, 'xlim');
%         yL = get(hAxis1, 'ylim');
%         set(hAxis2, 'xlim', xL);
%         set(hAxis2, 'ylim', yL);
%         axis(hAxis2, 'manual');
%     %     w_colors = floor((w_field ./ max(w_field))*255)+1;
%         set(gcf, 'doublebuffer', 'on');
%         for i=1:length(xPosV)
%             patch([xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2], [yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2], w_field(i));
%         end
%         axis([hAxis1 hAxis2], 'ij');
%         kids = get(hAxis2, 'children');
%     %     set(kids, 'edgecolor', 'none');
%         cMap = colormap('jet');
%         set(hAxis2, 'color', cMap(1,:));

    xPosV = xPosV/10;
    yPosV = yPosV/10;
    beamlet_delta_x = beamlet_delta_x/10;
    beamlet_delta_y = beamlet_delta_y/10;

    IMBase.goals.PBMargin = 0.5;
    IMBase.goals.structNum = 2;
    IMBase.params.xyDownsampleIndex = 1;
    IMBase.goals.isTarget(1) = 'y';
    IMBase.goals.xySampleRate = 2;
    IMBase.params.numCTSamplePts = 300;

    IMBase.beams(1).zRel = 0;
    IMBase.beams(1).xRel =  IMBase.beams(1).isodistance * sindeg(IMBase.beams(1).gantryAngle);
    IMBase.beams(1).yRel =  IMBase.beams(1).isodistance * cosdeg(IMBase.beams(1).gantryAngle);

    %IMBase.params.algorithm = 'VMC++';
    %IMCalc = IMRTP_MapCheck(IMBase);

    IMBase.params.algorithm = 'DPM';
    IMCalc = IMRTP_MapCheck(IMBase, planC, stateS, xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
    % JC. 11 Aug, 05 Added more arguments to IMRTP_MapCheck,

    %disp('OK_BIG1'); pause(5);
    % JC. Aug. 3, 2005
    % Add fields to IMCalc, necessary for generateDPMInfluence.m
    IMCalc.beams.couchAngle = 0;
    IMCalc.beams.collimatorAngle = 0;
    IMCalc.params.ScatterMethod = 'exponential';
    IMCalc.params.Scatter.Threshold = 0.01;
    IMCalc.params.Scatter.RandomStep = 30;
    IMCalc.beamlets = struct;
    %eval(['load Beam',num2str(indexBeam),'_w']);
    IMCalc = DPMInfluenceNEW(IMCalc,'photon.spectrum',nhist, OutputError,planC,stateS)
    %     filename = ['IMCalc_Beam',num2str(indexBeam)];
    %     save(filename, 'IMCalc');

    IMCalc = DoseScale(IMCalc);

    dose3D = getIMDose(IMCalc, w_field, 1);
    %dose3D = (dose3D/max(dose3D(:)));

%     if ~exist('dose')
%         dose = dose3D;
%     else
%         dose = dose + dose3D;
%     end

    filename = ['dose3D_',num2str(indexBeam)];
    currentDir = pwd;
    mkdir([getCERRPath,'planCheck\MC\doseMC']);
    cd([getCERRPath,'planCheck\MC\doseMC']);
    save(filename, 'dose3D');
    save dM dM;
    cd(currentDir);
    %    eval(['save IMCalc_Beam',num2str(indexBeam),' IMCalc']);

   clear  dose3D;

end

% dose3D = dose;

%showIMDosePlanCheck;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%IMMOMENTS PRTools routine for computing central moments of object images
%
%	  M = IMMOMENTS(A,TYPE,MOMENTS)
%
% INPUT
%   A        Dataset with object images dataset
%   TYPE     Desired type of moments
%   MOMENTS  Desired moments
%
% OUTPUT
%   M        Dataset with moments replacing images
%
% DESCRIPTION
% Computes for the image A a (1*N) vector M moments as defined by TYPE
% and MOMENTS. The following types are supported:
%
% TYPE = 'none'     Standard moments as specified in the Nx2 array MOMENTS.
%                   Moments are computed with respect to the image center.
%                   This is the default for TYPE.
%                   Default MOMENTS = [1 0; 0 1];
% TYPE = 'central'  Central moments as specified in the Nx2 array MOMENTS.
%                   Moments are computed with respect to the image mean
%                   Default MOMENTS = [2 0; 1 1; 0 2], which computes
%                   the variance in the x-direction (horizontal), the
%                   covariance between x and y and the variance in the
%                   y-direction (vertical).
% TYPE = 'scaled'   Scale-invariant moments as specified in the Nx2 array
%                   MOMENTS. Default MOMENTS = [2 0; 1 1; 0 2].
%                   After: M. Sonka et al.,
%                   Image processing, analysis and machine vision.
% TYPE = 'hu'       Calculates 7 moments of Hu, invariant to translation,
%                   rotation and scale.
%                   After: M. Sonka et al.,
%                   Image processing, analysis and machine vision.
% TYPE = 'zer'      Calculates the Zernike moments up to the order as
%                   specified in the scalar MOMENTS (1 <= MOMENTS <= 12).
%                   MOMENTS = 12 generates in total 47 moments.
%                   After: A. Khotanzad and Y.H. Hong, Invariant image
%                   recognition by Zernike moments, IEEE-PAMI, vol. 12,
%                   no. 5, 1990, 489-497.
%
% See DATASETS

% Copyright: D. de Ridder, R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = immoments(a,type,mom)

	if nargin < 3, mom = []; end
	if nargin < 2 | isempty(type), type = 'none'; end
	isobjim(a);
	[m,k] = size(a);
	im = data2im(a);
	for i = 1:m
		imi = im(:,:,i);
		switch type
		case {'none'}
			if isempty(mom)
				mom = [1 0; 0 1];
			end
			mout = moments(imi,mom(:,1),mom(:,2),0,0);
		case {'central'}
			if isempty(mom)
				mom = [2 0; 1 1; 0 2];
			end
			mout = moments(imi,mom(:,1),mom(:,2),1,0);
		case {'scaled'}
			if isempty(mom)
				mom = [2 0; 1 1; 0 2];
			end
			mout = moments(imi,mom(:,1)',mom(:,2)',1,1);
		case {'hu' 'Hu'}
			mout = hu_moments(imi);
		case {'zer' 'zernike' 'Zernike'}
			if isempty(mom)
				mom = 12;
			end
			mout = zernike_moments(imi,mom);
		otherwise
			error('Moments should be of type none, central, scaled, hu or zer')
		end
		if i==1
			b = zeros(m,length(mout));
		end
		b(i,:) = mout;
	end
	b = setdata(a,b);

% M = MOMENTS (IM, P, Q, CENTRAL, SCALED)
%
% Calculates moments of order (P+Q) (can be arrays of indentical length)
% on image IM. If CENTRAL is set to 1 (default: 0), returns translation-
% invariant moments; if SCALED is set to 1 (default: 0), returns scale-
% invariant moments.
%
% After: M. Sonka et al., Image processing, analysis and machine vision.

function m = moments (im,p,q,central,scaled)

	if (nargin < 5), scaled = 0; 	end;
	if (nargin < 4), central = 0; end;
  if (nargin < 3)
  	error ('Insufficient number of parameters.');
  end;

	if (length(p) ~= length(q))
  	error ('Arrays P and Q should have equal length.');
  end;

  if (scaled & ~central)
  	error ('Scale-invariant moments should always be central.');
  end;

  % xx, yy are grids with co-ordinates
  [xs,ys] = size(im);
  [xx,yy] = meshgrid(-(ys-1)/2:1:(ys-1)/2,-(xs-1)/2:1:(xs-1)/2);

	if (central)

  	% Calculate zeroth and first order moments
	  m00 = sum(sum(im));
	  m10 = sum(sum(im.*xx));
	  m01 = sum(sum(im.*yy));

    % This gives the center of gravity
    xc  = m10/m00;
    yc  = m01/m00;

    % Subtract this from the grids to center the object
    xx  = xx - xc;
    yy  = yy - yc;

  end;

  % Calculate moment(s) (p,q).
  for i = 1:length(p)
		m(i) = sum(sum((xx.^p(i)).*(yy.^q(i)).*im));
  end;

  if (scaled)

  	c = 1 + (p+q)/2;

    % m00 should be known, as scaled moments are always central
    m = m ./ (m00.^c);

	end;

return;

% M = HU_MOMENTS (IM)
%
% Calculates 7 moments of Hu on image IM, invariant to translation,
% rotation and scale.
%
% After: M. Sonka et al., Image processing, analysis and machine vision.

function m = hu_moments (im)

	p = [ 1 0 2 1 2 0 3 ];
	q = [ 1 2 0 2 1 3 0 ];

  n = moments(im,p,q,1,1);

  m(1) = n(2) + n(3);
  m(2) = (n(3) - n(2))^2   + 4*n(1)^2;
  m(3) = (n(7) - 3*n(4))^2 + (3*n(5) - n(6))^2;
  m(4) = (n(7) +   n(4))^2 + (  n(5) + n(6))^2;
  m(5) = (  n(7) - 3*n(4)) * (n(7) + n(4)) * ...
           (  (n(7) + n(4))^2 - 3*(n(5) + n(6))^2) + ...
         (3*n(5) -   n(6)) * (n(5) + n(6)) * ...
           (3*(n(7) + n(4))^2 -   (n(5) + n(6))^2);
  m(6) = (n(3) - n(2)) * ((n(7) + n(4))^2 - (n(5) + n(6))^2) + ...
          4*n(1) * (n(7)+n(4)) * (n(5)+n(6));
  m(7) = (3*n(5) -   n(6)) * (n(7) + n(4)) * ...
           (  (n(7) + n(4))^2 - 3*(n(5) + n(6))^2) - ...
         (  n(7) - 3*n(4)) * (n(5) + n(6)) * ...
           (3*(n(7) + n(4))^2 -   (n(5) + n(6))^2);

return;

% M = ZERNIKE_MOMENTS (IM, ORDER)
%
% Calculates Zernike moments up to and including ORDER (<= 12) on image IM.
% Default: ORDER = 12.

function m = zernike_moments (im, order)

  if (nargin < 2),             order = 12;                      end;
  if (order < 1 | order > 12), error ('order should be 1..12'); end;

  % xx, yy are grids with co-ordinates

  [xs,ys] = size(im);
  [xx,yy] = meshgrid(-(ys-1)/2:1:(ys-1)/2,-(xs-1)/2:1:(xs-1)/2);

  % Calculate center of mass and distance of any pixel to it

  m  = moments (im,[0 1 0],[0 0 1],0,0);
  xc = m(2)/m(1); yc = m(3)/m(1);
  xx = xx - xc; yy = yy - yc;

  len     = sqrt(xx.^2+yy.^2);
  max_len = max(max(len));

  % Map pixels to unit circle; prevent divide by zero.

  rho        = len/max_len;
  rho_tmp    = rho; rho_tmp(find(rho==0)) = 1;
  theta      = acos((xx/max_len)./rho_tmp);

  % Flip angle for pixels above center of mass

  yneg            = length(find(yy(:,1)<0));
  theta(:,1:yneg) = 2*pi - theta(:,1:yneg);

  % Calculate coefficients

  c = zeros(order,order);
  s = zeros(order,order);

  i = 1;
  for n = 2:order
    for l = n:-2:0
      r    = polynomial (n,l,rho);
      c    = sum(sum(r.*cos(l*theta)))*((n+1)/(pi*max_len^2));
      s    = sum(sum(r.*sin(l*theta)))*((n+1)/(pi*max_len^2));
      m(i) = sqrt(c^2+s^2);
      i    = i + 1;
    end;
  end;

return

function p = polynomial (n,l,rho)

  switch (n)
    case 2, switch (l)
        case 0, p = 2*(rho.^2)-1;
        case 2, p =   (rho.^2);
      end;
    case 3, switch (l)
        case 1, p = 3*(rho.^3)-2*rho;
        case 3, p =   (rho.^3);
      end;
    case 4, switch (l)
        case 0, p = 6*(rho.^4)-6*(rho.^2)+1;
        case 2, p = 4*(rho.^4)-3*(rho.^2);
        case 4, p =   (rho.^4);
      end;
    case 5, switch (l)
        case 1, p = 10*(rho.^5)-12*(rho.^3)+3*rho;
        case 3, p =  5*(rho.^5)- 4*(rho.^3);
        case 5, p =    (rho.^5);
      end;
    case 6, switch (l)
        case 0, p = 20*(rho.^6)-30*(rho.^4)+12*(rho.^2)-1;
        case 2, p = 15*(rho.^6)-20*(rho.^4)+ 6*(rho.^2);
        case 4, p =  6*(rho.^6)- 5*(rho.^4);
        case 6, p =    (rho.^6);
      end;
    case 7, switch (l)
        case 1, p = 35*(rho.^7)-60*(rho.^5)+30*(rho.^3)-4*rho;
        case 3, p = 21*(rho.^7)-30*(rho.^5)+10*(rho.^3);
        case 5, p =  7*(rho.^7)- 6*(rho.^5);
        case 7, p =    (rho.^7);
      end;
    case 8, switch (l)
        case 0, p = 70*(rho.^8)-140*(rho.^6)+90*(rho.^4)-20*(rho.^2)+1;
        case 2, p = 56*(rho.^8)-105*(rho.^6)+60*(rho.^4)-10*(rho.^2);
        case 4, p = 28*(rho.^8)- 42*(rho.^6)+15*(rho.^4);
        case 6, p =  8*(rho.^8)-  7*(rho.^6);
        case 8, p =    (rho.^8);
      end;
    case 9, switch (l)
        case 1, p = 126*(rho.^9)-280*(rho.^7)+210*(rho.^5)-60*(rho.^3)+5*rho;
        case 3, p =  84*(rho.^9)-168*(rho.^7)+105*(rho.^5)-20*(rho.^3);
        case 5, p =  36*(rho.^9)- 56*(rho.^7)+ 21*(rho.^5);
        case 7, p =   9*(rho.^9)-  8*(rho.^7);
        case 9, p =     (rho.^9);
      end;
    case 10, switch (l)
        case  0, p = 252*(rho.^10)-630*(rho.^8)+560*(rho.^6)-210*(rho.^4)+30*(rho.^2)-1;
        case  2, p = 210*(rho.^10)-504*(rho.^8)+420*(rho.^6)-140*(rho.^4)+15*(rho.^2);
        case  4, p = 129*(rho.^10)-252*(rho.^8)+168*(rho.^6)- 35*(rho.^4);
        case  6, p =  45*(rho.^10)- 72*(rho.^8)+ 28*(rho.^6);
        case  8, p =  10*(rho.^10)-  9*(rho.^8);
        case 10, p =     (rho.^10);
      end;
    case 11, switch (l)
        case  1, p = 462*(rho.^11)-1260*(rho.^9)+1260*(rho.^7)-560*(rho.^5)+105*(rho.^3)-6*rho;
        case  3, p = 330*(rho.^11)- 840*(rho.^9)+ 756*(rho.^7)-280*(rho.^5)+ 35*(rho.^3);
        case  5, p = 165*(rho.^11)- 360*(rho.^9)+ 252*(rho.^7)- 56*(rho.^5);
        case  7, p =  55*(rho.^11)-  90*(rho.^9)+  36*(rho.^7);
        case  9, p =  11*(rho.^11)-  10*(rho.^9);
        case 11, p =     (rho.^11);
      end;
    case 12, switch (l)
        case  0, p = 924*(rho.^12)-2772*(rho.^10)+3150*(rho.^8)-1680*(rho.^6)+420*(rho.^4)-42*(rho.^2)+1;
        case  2, p = 792*(rho.^12)-2310*(rho.^10)+2520*(rho.^8)-1260*(rho.^6)+280*(rho.^4)-21*(rho.^2);
        case  4, p = 495*(rho.^12)-1320*(rho.^10)+1260*(rho.^8)- 504*(rho.^6)+ 70*(rho.^4);
        case  6, p = 220*(rho.^12)- 495*(rho.^10)+ 360*(rho.^8)-  84*(rho.^6);
        case  8, p =  66*(rho.^12)- 110*(rho.^10)+  45*(rho.^8);
        case 10, p =  12*(rho.^12)-  11*(rho.^10);
        case 12, p =     (rho.^12);
      end;
  end;

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function result = Main_Process(Img, n)
if ndims(Img) == 3
    I = rgb2gray(Img);
else
    I = Img;
end
g1 = [0 1 0
    0 1 0
    0 1 0];
g2 = [0 0 0
    1 1 1
    0 0 0];
g3 = [0 0 1
    0 1 0
    1 0 0];
g4 = [1 0 0
    0 1 0
    0 0 1];
g5 = [0 1 0
    1 1 1
    0 1 0];
Gi1 = Multi_Process(I, g1, n);
Gi2 = Multi_Process(I, g2, n);
Gi3 = Multi_Process(I, g3, n);
Gi4 = Multi_Process(I, g4, n);
Gi5 = Multi_Process(I, g5, n);
G{1} = Gi1;
G{2} = Gi2;
G{3} = Gi3;
G{4} = Gi4;
G{5} = Gi5;
ua1 = Coef(Gi1, G);
ua2 = Coef(Gi2, G);
ua3 = Coef(Gi3, G);
ua4 = Coef(Gi4, G);
ua5 = Coef(Gi5, G);
u = [ua1, ua2, ua3, ua4, ua5];
u = u/sum(u);
Gf1 = Edge_One(G, u);
result = Gf1;

author: decouples


text:
function asc = read_eyelink_asc(filename)

% READ_EYELINK_ASC reads the header information, input triggers, messages
% and all data points from an Eyelink *.asc file. The output events are
% represented as matlab tables (after Aug 2022)
%
% Use as
%   asc = read_eyelink_asc(filename)

% Copyright (C) 2010-2015, Robert Oostenveld
% Copyright (C) 2022, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% read the whole file at once
fid = fopen_or_error(filename, 'rt');
aline = fread(fid, inf, 'char=>char');          % returns a single long string
fclose(fid);

aline(aline==uint8(sprintf('\r'))) = [];        % remove cariage return
aline = tokenize(aline, uint8(newline));        % split on newline

sel = startsWith(aline, '0') | ...
      startsWith(aline, '1') | ...
      startsWith(aline, '2') | ...
      startsWith(aline, '3') | ...
      startsWith(aline, '4') | ...
      startsWith(aline, '5') | ...
      startsWith(aline, '6') | ...
      startsWith(aline, '7') | ...
      startsWith(aline, '8') | ...
      startsWith(aline, '9');

datline = aline(sel);
aline   = aline(~sel);

% according to the eyelink documentation, the sample lines have different
% flavours, depending on how the data was collected (monocular or
% binocular), and possibly the conversion settings from edf2asc
% Monocular,                              <time> <xp>  <yp>  <ps>
% Monocular, with velocity                <time> <xp>  <yp>  <ps>  <xv>  <yv>
% Monocular, with resolution              <time> <xp>  <yp>  <ps>  <xr>  <yr>
% Monocular, with velocity and resolution <time> <xp>  <yp>  <ps>  <xv>  <yv>  <xr>  <yr>
% Binocular,                              <time> <xpl> <ypl> <psl> <xpr> <ypr> <psr>
% Binocular, with velocity                <time> <xpl> <ypl> <psl> <xpr> <ypr> <psr> <xvl> <yvl> <xvr> <yvr>
% Binocular, with resolution              <time> <xpl> <ypl> <psl> <xpr> <ypr> <psr> <xr> <yr>
% Binocular, with velocity and resolution <time> <xpl> <ypl> <psl> <xpr> <ypr> <psr> <xvl> <yvl> <xvr> <yvr><xr> <yr>
%
% Then, there may be additional columns, if recorded with cornea reflection
% mode:
%
% MONOCULAR Corneal Reflection (CR) Samples
% "..." if no warning for sample
% first character is "I" if sample was interpolated
% second character is "C" if CR missing
% third character is "R" if CR recovery in progress
%
% BINOCULAR Corneal Reflection (CR) Samples?
% "....." if no warning for sample
% first character is "I" if sample was interpolated
% second character is "C" if LEFT CR missing
% third character is "R" if LEFT CR recovery in progress
% fourth character is "C" if RIGHT CR missing
% fifth character is "R" if RIGHT CR recovery in progress
%
% Then, there may be additional columns, if data collection was done in
% remote mode:
%
% Data files recorded using the Remote Mode have extra columns to encode the
% target distance, position, and eye/target status information. The first three
% columns are:
% <target x>: X position of the target in camera coordinate (a value from 0 to 10000).
% Returns "MISSING_DATA" (-32768) if target is missing.
% <target y>: Y position of the target in camera coordinate (a value from 0 to 10000).
% Returns "MISSING_DATA" (-32768) if target is missing.
% <target distance>: Distance between the target and camera (in millimeters).
% Returns "MISSING_DATA" (-32768) if target is missing.
% The next thirteen fields represent warning messages for that sample relating to
% the target and eye image processing."............." if no warning for target and eye image
% first character is "M" if target is missing
% second character is "A" if extreme target angle occurs
% third character is "N" if target is near eye so that the target window and eye window overlap
% fourth character is "C" if target is too close
% fifth character is "F" if target is too far
% sixth character is "T" if target is near top edge of the camera image
% seventh character is "B" if target is near bottom edge of the camera image
% eighth character is "L" if target is near left edge of the camera image
% ninth character is "R" if target is near right edge of the camera image
% tenth character is "T" if eye is near top edge of the camera image
% eleventh character is "B" if eye is near bottom edge of the camera image
% twelfth character is "L" if eye is near left edge of the camera image
% thirteenth character is "R" if eye is near right edge of the camera image
% For a binocular recording, there will be seventeen target/eye status columns,
% with the last eight columns reporting the warning messages for the left and
% right eyes separately.

% So, anecdotally, parsing the file is a bit tedious, and not very robust.
% Reading and converting per line is super slow, so here I chose to pipe
% the datalines through a checker that removes all '...' '.....'
% '.............' and '.................' (or containing warning) stuff
% so that we end up with more or less well behaved data.

% check if the formatting for all datlines is similar
seltab = char(datline)==sprintf('\t');
ntab   = sum(seltab,2);

% check whether all lines have the same number of columns
assert(all(ntab==ntab(1)));

% identify the chunks of consecutive '...' and comments
str = char(datline)';
str(end+1, :) = sprintf('\t');
siz = size(str);
str = uint8(str(:)');

boolval = str==46 | (str>=65&str<=90); % dots or capital letters
begsmp  = find(boolval==1 & [0 boolval(1:end-1)]==0);
endsmp  = find(boolval==1 & [boolval(2:end) 1]==0);
sel     = ~ismember(endsmp-begsmp+1, [3 5 13 17]);
boolval(begsmp(sel)) = false;
str(boolval) = 32;

str = reshape(char(str), siz); % now all the warnings etc should have been removed
datline = cellstr(str(1:end-1,:)');

% check again if the formatting for all datlines is similar
nline  = numel(datline);
seltab = char(datline)==sprintf('\t');
ntab   = sum(seltab,2);

% check whether all lines have the same number of columns
assert(all(ntab==ntab(1)));

fmt = repmat('%s ', [1 ntab(1)+1]);
fmt = fmt(1:end-1);

% convert datline in a long string again
datline = char(datline)';
datline(end+1, :) = sprintf('\t');
datline = datline(:)';

C   = textscan(datline, fmt, nline, 'Delimiter', sprintf('\t'));
dat = zeros(nline, numel(C)) + nan;
for i=1:numel(C)
  notok = strcmp(C{i}, '.')|strcmp(C{i}, '  ')|strcmp(C{i}, '     ')|cellfun('isempty', C{i})|startsWith(C{i}, '. ');
  ok    = ~notok;
  tmp = reshape([char(C{i}(ok))';repmat(' ',1,sum(ok))],[],1)';

  % convert to floats, much faster than str2double
  if sum(ok)
    tmp = textscan(tmp, '%f');
    dat(ok, i) = tmp{1};
  end
end
asc.dat = dat';

% remove any all-nan rows
asc.dat(sum(~isfinite(asc.dat),2)==size(asc.dat,2), :) = [];

selsfix    = startsWith(aline, 'SFIX');
selefix    = startsWith(aline, 'EFIX');
selssacc   = startsWith(aline, 'SSACC');
selesacc   = startsWith(aline, 'ESACC');
selsblink  = startsWith(aline, 'SBLINK');
seleblink  = startsWith(aline, 'EBLINK');
selmsg     = startsWith(aline, 'MSG');
selhdr     = startsWith(aline, '**');
selinput   = startsWith(aline, 'INPUT');

asc.sfix   = tocell(aline(selsfix)); % convert cell-vector of lines into a cell array
asc.efix   = tocell(aline(selefix));
asc.ssacc  = tocell(aline(selssacc));
asc.esacc  = tocell(aline(selesacc));
asc.sblink = tocell(aline(selsblink));
asc.eblink = tocell(aline(seleblink));
asc.msg    = tocell(aline(selmsg), 1);
asc.header = char(aline(selhdr));
asc.input  = tocell(aline(selinput));

% convert the cell-arrays into tables
if ~isempty(asc.sfix),   asc.sfix   = totable(asc.sfix(:, 2:end),   {'eye' 'stime'}); end
if ~isempty(asc.sblink), asc.sblink = totable(asc.sblink(:, 2:end), {'eye' 'stime'}); end
if ~isempty(asc.ssacc),  asc.ssacc  = totable(asc.ssacc(:, 2:end),  {'eye' 'stime'}); end
if ~isempty(asc.input),  asc.input  = totable(asc.input(:, 2:end),  {'timestamp', 'value'}, [1 2]); end

% these lines may have been formatted in different ways
if ~isempty(asc.efix)
  try
    asc.efix = totable(asc.efix(:,2:end), {'eye', 'stime', 'etime', 'dur', 'axp', 'ayp', 'aps'});
  catch
    asc.efix = totable(asc.efix(:,2:end), {'eye', 'stime', 'etime', 'dur', 'axp', 'ayp', 'aps', 'xr', 'yr'});
  end
end
if ~isempty(asc.eblink)
  asc.eblink = totable(asc.eblink(:,2:end), {'eye', 'stime', 'etime', 'dur'});
end
if ~isempty(asc.esacc)
  try
    asc.esacc = totable(asc.esacc(:,2:end), {'eye', 'stime', 'etime', 'dur', 'sxp', 'syp', 'exp', 'eyp', 'ampl', 'pv'});
  catch
    asc.esacc = totable(asc.esacc(:,2:end), {'eye', 'stime', 'etime', 'dur', 'sxp', 'syp', 'exp', 'eyp', 'ampl', 'pv', 'xr', 'yr'});
  end
end

function lineout = convertline(linein, msgflag)

% split the line by horizontal tabs, and the first output element
% thereof by spaces.

%tmp     = strrep(linein, sprintf('\t'), ' ');
lineout = strsplit(linein, '\t');
lineout = [strsplit(lineout{1}, ' ') lineout(2:end)];
if msgflag
  tmp = strsplit(lineout{2}, ' ');
  lineout = [lineout(1) tmp(1) {sprintf('%s ', tmp{2:end})}];
end

function tableout = totable(cellin, fnames, convert)

if nargin<3
  convert = 2:numel(fnames);
end

% convert into table, and do a str2double conversion for the numeric data,
% assuming only the first variable to be kept as a string

tableout = cell2table(cellin, 'VariableNames', fnames);
for k = convert
  tableout.(fnames{k}) = str2double(tableout.(fnames{k}));
end

function cellout = tocell(textin, msgflag)

if nargin<2
  msgflag = false;
end

cellout = {};
for k = 1:numel(textin)
  tline   = convertline(textin{k}, msgflag);
  if k==1
    cellout = cell(numel(textin), numel(tline));
  end
  cellout(k,:) = tline;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function tIncCh = tIncCh_Initialize(tIncCh, d, ml)
if isempty(tIncCh)
    tIncCh = [ones(size(d,1), size(ml,1)); ml'];
end

author: BUNPC


text:
function a = subsasgn(a,S,val)
% A.FIELD = val; %Subscripted assignment
%
% SEE ALSO, SET

if isempty(a)
    % Ex. a(1) = chanhandle(filename,channel_num);
    if isequal(S.type, '()')&isa(val,'chanhandle')&isequal(S.subs{:},1)
        a = val;
        return;
    else
        newline= sprintf('\n');
        str = ['Could not initialize chanhandle object with given parameters',newline,...
        'Please note that:',newline,...
        '  1. Right-hand side should be a valid chanhandle object',newline,...
        '  2. Only the first channel should be set',newline,...
        '  3. Property of chanhandle cannot be set until the object itself is initialzed'];
        error(str);
    end;
end;
if nargin == 3
    switch length(S)
    case 1
        switch S.type
        case '()'
            if isa(val,'chanhandle')
                if isequal(S.subs{1}, 1:length(val)) % a(i) = chanhandle_obj;
                    a.handle(S.subs{1}) = getfield(struct(val),'handle');
                elseif isequal(length(S.subs{:}),length(getfield(struct(val),'handle')))
                    % Ex. a(i1,i2,...) = [chanhandle_obj1,chanhandle_obj2,...];
                    for i = 1:length(S.subs)
                        a.handle(S.subs{i}) = getfield(struct(val(i)),'handle');
                    end;
                else
                    error('Input chanhandle structure and Channel number do not match');
                end;
            elseif isstruct(val)  %Ex. a(i) = chanhandle_obj.handle;
                if isequal(length(S.subs{:}),length(val))
                    for i = 1:length(S.subs)
                        a.handle(S.subs{i}) = val(i);
                    end;
                else
                    error('Input chanhandle structure and Channel number do not match');
                end;
            else
                error('RHS must be either a valid structure or a chanhandle object');
            end;
        case '.'
            a = set(a,S.subs,val);      % a.field = val;
        end;
    case 2
        if strcmp(S(1).type,'()')&strcmpi(S(2).type,'.')    % a(i).field = val;
            a = set(a,S(1).subs{:},S(2).subs,val);
        elseif strcmp(S(1).type,'.')&strcmpi(S(2).type,'()')% a.field(i) = val;
            a = set(a,S(1).subs,S(2).subs{:},val);
        elseif strcmp(S(1).type,'.')&strcmpi(S(2).type,'.') % a.field1.field2 = val;
            a = set(a,[S(1).subs,'.',S(2).subs],val);
        else
            error('Sorry - output structure not supported');
        end;
    otherwise
        error('Sorry - output structure not supported');
    end;
else
    error('incorrect number of arguments');
end;


text:
function dist = distancePoints(p1, p2, varargin)
%DISTANCEPOINTS Compute distance between two points.
%
%   D = distancePoints(P1, P2)
%   Return the Euclidean distance between points P1 and P2.
%
%   If P1 and P2 are two arrays of points, result is a N1-by-N2 array
%   containing distance between each point of P1 and each point of P2.
%
%   D = distancePoints(P1, P2, NORM)
%   Compute distance using the specified norm. NORM=2 corresponds to usual
%   euclidean distance, NORM=1 corresponds to Manhattan distance, NORM=inf
%   is assumed to correspond to maximum difference in coordinate. Other
%   values (>0) can be specified.
%
%   D = distancePoints(..., 'diag')
%   compute only distances between P1(i,:) and P2(i,:).
%
%   See also:
%   points2d, minDistancePoints, nndist, hausdorffDistance
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Copyright 2009 INRA - Cepia Software Platform.
% created the 24/02/2004.
%

%   HISTORY :
%   25/05/2004: manage 2 array of points
%   07/04/2004: add option for computing only diagonal.
%   30/10/2006: generalize to any dimension, and manage different norms
%   03/01/2007: bug for arbitrary norm, and update doc
%   28/08/2007: fix bug for norms 2 and infinite, in diagonal case

%% Setup options

% default values
diag = false;
norm = 2;

% check first argument: norm or diag
if ~isempty(varargin)
    var = varargin{1};
    if isnumeric(var)
        norm = var;
    elseif strncmp('diag', var, 4)
        diag = true;
    end
    varargin(1) = [];
end

% check last argument: diag
if ~isempty(varargin)
    var = varargin{1};
    if strncmp('diag', var, 4)
        diag = true;
    end
end

% number of points in each array and their dimension
n1  = size(p1, 1);
n2  = size(p2, 1);
d   = size(p1, 2);

if diag
    % compute distance only for apparied couples of pixels
    dist = zeros(n1, 1);

    if norm == 2
        % Compute euclidian distance. this is the default case
        % Compute difference of coordinate for each pair of point
        % and for each dimension. -> dist is a [n1*n2] array.
        for i = 1:d
            dist = dist + (p2(:,i)-p1(:,i)).^2;
        end
        dist = sqrt(dist);

    elseif norm == inf
        % infinite norm corresponds to maximal difference of coordinate
        for i = 1:d
            dist = max(dist, abs(p2(:,i)-p1(:,i)));
        end

    else
        % compute distance using the specified norm.
        for i = 1:d
            dist = dist + power((abs(p2(:,i)-p1(:,i))), norm);
        end
        dist = power(dist, 1/norm);
    end
else
    % compute distance for all couples of pixels
    dist = zeros(n1, n2);

    if norm == 2
        % Compute euclidian distance. This is the default case.
        % Compute difference of coordinate for each pair of point
        % and for each dimension. -> dist is a [n1*n2] array.
        for i = 1:d
            % equivalent to:
            % dist = dist + ...
            %   (repmat(p1(:,i), [1 n2])-repmat(p2(:,i)', [n1 1])).^2;
            dist = dist + bsxfun (@minus, p1(:,i), p2(:, i)').^2;
        end
        dist = sqrt(dist);

    elseif norm == inf
        % infinite norm corresponds to maximal difference of coordinate
        for i = 1:d
            dist = max(dist, abs(bsxfun (@minus, p1(:,i), p2(:, i)')));
        end

    else
        % compute distance using the specified norm.
        for i = 1:d
            % equivalent to:
            % dist = dist + power((abs(repmat(p1(:,i), [1 n2]) - ...
            %     repmat(p2(:,i)', [n1 1]))), norm);
            dist = dist + power(abs(bsxfun(@minus, p1(:,i), p2(:, i)')), norm);
        end
        dist = power(dist, 1/norm);
    end
end


text:
classdef InteriorSubCellsConnecComputer < handle

    properties (GetAccess = public, SetAccess = private)
        connec
        xCoordsIso
        cellContainingSubcell
    end

    properties (Access = private)
        isSubCellInterior
        subCellCases
        allSubCellsConnec
        xNodesInSubCells
    end

    properties (Access = private)
        localCells
        allSubCellsConnecParams
        cutElems
        nSubCellsByElem
    end

    methods (Access = public)

        function obj = InteriorSubCellsConnecComputer(cParams)
            obj.init(cParams);
            obj.compute();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.allSubCellsConnecParams = cParams.allSubCellsConnecParams;
            obj.cutElems                = cParams.cutElems;
            obj.nSubCellsByElem         = cParams.nSubCellsByElem;
            obj.isSubCellInterior       = cParams.isSubCellInterior;
        end

        function compute(obj)
            obj.computeAllSubCellsConnec();
            obj.computeLocalCellContainingSubCell();
            obj.computeConnec();
            obj.computeXcoordsIso();
            obj.computeCellContainingSubCell();
        end

        function computeAllSubCellsConnec(obj)
            s = obj.allSubCellsConnecParams;
            a = AllSubCellsConnecComputer(s);
            a.compute();
            obj.allSubCellsConnec = a.allSubCellsConnec;
            obj.xNodesInSubCells  = a.xNodesInSubCells;
        end

        function cells = computeLocalCellContainingSubCell(obj)
            isInterior = obj.isSubCellInterior;
            nCutElem   = size(isInterior,2);
            localSubCellsInCell = repmat(1:nCutElem,obj.nSubCellsByElem,1);
            cells = localSubCellsInCell(isInterior);
            obj.localCells = cells;
        end

        function computeConnec(obj)
            allConnec  = obj.allSubCellsConnec;
            isInterior = obj.isSubCellInterior(:);
            obj.connec = allConnec(isInterior,:);
        end

        function computeXcoordsIso(obj)
            isInterior = obj.isSubCellInterior(:);
            x = obj.xNodesInSubCells;
            nDim = size(x,3);
            nElem = sum(isInterior);
            nNode = size(x,1);
            obj.xCoordsIso = zeros(nDim,nNode,nElem);
            for idim = 1:nDim
                obj.xCoordsIso(idim,:,:) = x(:,isInterior,idim);
            end
        end

        function computeCellContainingSubCell(obj)
            localToGlobalCut = obj.cutElems;
            lCells  = obj.localCells;
            globalCell = localToGlobalCut(lCells);
            obj.cellContainingSubcell = globalCell;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan