text:
function medit_to_fem ( prefix )

%*****************************************************************************80
%
%% MAIN is the main program for TRIANGLE_TO_FEM.
%
%  Discussion:
%
%    MEDIT_TO_FEM converts mesh data from MEDIT to FEM format.
%
%  Usage:
%
%    medit_to_fem prefix
%
%    where 'prefix' is the common filename prefix:
%
%    * 'prefix'.mesh is the MEDIT mesh file.
%    * 'prefix'_nodes.txt will contain the node coordinates.
%    * 'prefix'_elements.txt will contain the element node connectivity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 November 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MEDIT_TO_FEM\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Read a mesh description created by MEDIT\n' );
  fprintf ( 1, '  * "prefix".mesh.\n' );
  fprintf ( 1, '  Write two simple FEM files listing nodes and elements.\n' );
  fprintf ( 1, '  * "prefix"_nodes.txt, node coordinates.\n' );
  fprintf ( 1, '  * "prefix"_elements.txt, element connectivity.\n' );
%
%  Get the filename prefix.
%
  if ( nargin < 1 )

    prefix = input ( 'Enter the filename prefix:  ' );

  end
%
%  Create the filenames.
%
  medit_filename = strcat ( prefix, '.mesh' );
  fem_node_filename = strcat ( prefix, '_nodes.txt' );
  fem_element_filename = strcat ( prefix, '_elements.txt' );
%
%  Read MEDIT sizes.
%
  [ dim, vertices, edges, triangles, quadrilaterals, tetrahedrons, ...
    hexahedrons ] = mesh_size_read ( medit_filename );
%
%  Report sizes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of dimensions = %d\n', dim );
  fprintf ( 1, '  Number of vertices = %d\n', vertices );
  fprintf ( 1, '  Number of edges = %d\n', edges );
  fprintf ( 1, '  Number of triangles = %d\n', triangles );
  fprintf ( 1, '  Number of quadrilaterals = %d\n', quadrilaterals );
  fprintf ( 1, '  Number of tetrahedrons = %d\n', tetrahedrons );
  fprintf ( 1, '  Number of hexahedrons = %d\n', hexahedrons );
%
%  Read MEDIT data.
%
  [ vertex_coordinate, vertex_label, edge_vertex, edge_label, ...
    triangle_vertex, triangle_label, quadrilateral_vertex, ...
    quadrilateral_label, tetrahedron_vertex, tetrahedron_label, ...
    hexahedron_vertex, hexahedron_label ] = mesh_data_read ( medit_filename, ...
    dim, vertices, edges, triangles, quadrilaterals, tetrahedrons, ...
    hexahedrons );
%
%  Choose the FEM data.
%
%  We need to assume that there is only one element type.
%  If there are elements of multiple dimension, take the highest.
%
  m = dim;
  node_num = vertices;
  r8mat_write ( fem_node_filename, dim, vertices, vertex_coordinate );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created node coordinate file "%s"\n', fem_node_filename );

  if ( 0 < hexahedrons && dim == 3 )

    element_order = 8;
    element_num = hexahedrons;
    i4mat_write ( fem_element_filename, element_order, element_num, ...
      hexahedron_vertex );

  elseif ( 0 < tetrahedrons && dim == 3 )

    element_order = 4;
    element_num = tetrahedrons;
    i4mat_write ( fem_element_filename, element_order, element_num, ...
      tetrahedron_vertex );

  elseif ( 0 < quadrilaterals && dim == 2 )

    element_order = 4;
    element_num = quadrilaterals;
    i4mat_write ( fem_element_filename, element_order, element_num, ...
      quadrilateral_vertex );

  elseif ( 0 < triangles && dim == 2 )

    element_order = 3;
    element_num = triangles;
    i4mat_write ( fem_element_filename, element_order, element_num, ...
      triangle_vertex );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'MEDIT_TO_FEM - Fatal error!\n' );
    fprintf ( 1, '  Unexpected combination of spatial dimension\n' );
    fprintf ( 1, '  and number of nonzero objects.\n' );
    error ( 'MEDIT_TO_FEM - Fatal error!' );

  end

  fprintf ( 1, '  Created element connectivity file "%s"\n', ...
    fem_element_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MEDIT_TO_FEM:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function c = ch_cap ( c )

%*****************************************************************************80
%
%% CH_CAP capitalizes a single character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to capitalize.
%
%    Output, character C, the capitalized character.
%
  if ( 'a' <= c && c <= 'z' )
    c = c + 'A' - 'a';
  end

  return
end
function truefalse = ch_eqi ( c1, c2 )

%*****************************************************************************80
%
%% CH_EQI is a case insensitive comparison of two characters for equality.
%
%  Example:
%
%    CH_EQI ( 'A', 'a' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C1, C2, the characters to compare.
%
%    Output, logical TRUEFALSE, is TRUE (1) if the characters are equal.
%
  FALSE = 0;
  TRUE = 1;

  if ( ch_cap ( c1 ) == ch_cap ( c2 ) )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

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

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  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, integer TABLE(M,N), the points.
%

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ vertex_coordinate, vertex_label, edge_vertex, edge_label, ...
  triangle_vertex, triangle_label, quadrilateral_vertex, ...
  quadrilateral_label, tetrahedron_vertex, tetrahedron_label, ...
  hexahedron_vertex, hexahedron_label ] = mesh_data_read ( filename, dim, ...
  vertices, edges, triangles, quadrilaterals, tetrahedrons, hexahedrons )

%*****************************************************************************80
%
%% MESH_READ reads a MESH file defining a mesh and returns the data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Pascal Frey,
%    MEDIT: An interactive mesh visualization software,
%    Technical Report RT-0253,
%    Institut National de Recherche en Informatique et en Automatique,
%    03 December 2001.
%
%  Parameters:
%
%    Input, string FILENAME, the name of the MESH file.
%
%    Input, integer DIM, the spatial dimension, which should be 2 or 3.
%
%    Input, integer VERTICES, the number of vertices.
%
%    Input, integer EDGES, the number of edges (may be 0).
%
%    Input, integer TRIANGLES, the number of triangles (may be 0).
%
%    Input, integer QUADRILATERALS, the number of quadrilaterals (may be 0).
%
%    Input, integer TETRAHEDRONS, the number of tetrahedrons (may be 0).
%
%    Input, integer HEXAHEDRONS, the number of hexahedrons (may be 0).
%
%    Output, real VERTEX_COORDINATE(DIM,VERTICES), the coordinates
%    of each vertex.
%
%    Output, integer VERTEX_LABEL(VERTICES), a label for each vertex.
%
%    Output, integer EDGE_VERTEX(2,EDGES), the vertices that form each edge.
%
%    Output, integer EDGE_LABEL(EDGES), a label for each edge.
%
%    Output, integer TRIANGLE_VERTEX(3,TRIANGLES), the vertices that form
%    each triangle.
%
%    Output, integer TRIANGLE_LABEL(TRIANGLES), a label for each triangle.
%
%    Output, integer QUADRILATERAL_VERTEX(4,QUADRILATERALS), the vertices that
%    form each quadrilateral.
%
%    Output, integer QUADRILATERAL_LABEL(QUADRILATERALS), a label for
%    each quadrilateral.
%
%    Output, integer TETRAHEDRON_VERTEX(4,TETRAHEDRONS), the vertices that
%    form each tetrahedron.
%
%    Output, integer TETRAHEDRON_LABEL(TETRAHEDRONS), a label for
%    each tetrahedron.
%
%    Output, integer HEXAHEDRON_VERTEX(8,HEXAHEDRONS), the vertices that form
%    each hexahedron.
%
%    Output, integer HEXAHEDRON_LABEL(HEXAHEDRONS), a label for each hexahedron.
%

%
%  Initialize everything to nothing.
%
  vertex_coordinate = [];
  vertex_label = [];
  edge_vertex = [];
  edge_label = [];
  triangle_vertex = [];
  triangle_label = [];
  quadrilateral_vertex = [];
  quadrilateral_label = [];
  tetrahedron_vertex = [];
  tetrahedron_label = [];
  hexahedron_vertex = [];
  hexahedron_label = [];
%
%  Open the file.
%
  unit = fopen ( filename, 'rt' );

  if ( unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file "%s".\n', filename );
    error ( 'MESH_DATA_READ - Error!' );
    return
  end
%
%  Read lines til you get alphanumerics and determine a "mode"
%
  line_num = 0;
  keyword = 'NONE';

  while ( 1 )

    text = fgetl ( unit );

    if ( text == -1 )
      break
    end

    line_num = line_num + 1;

    if ( s_len_trim ( text ) == 0 )
      keyword = 'NONE';
      continue
    end

    if ( text(1) == '#' )
      continue
    end
%
%  Remove initial blanks.
%

%
%  Expecting a keyword.
%
        if ( s_eqi ( text, 'CORNERS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'DIMENSION' ) )

      keyword = 'DIMENSION';

    elseif ( s_eqi ( text, 'EDGES' ) )

      keyword = 'EDGES';

    elseif ( s_eqi ( text, 'END' ) )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  END statement encountered.\n' );
      break

    elseif ( s_eqi ( text, 'HEXAHEDRA' ) || s_eqi ( text, 'HEXAHEDRONS' ) )

      keyword = 'HEXAHEDRONS';

    elseif ( s_begin ( text, 'MESHVERSIONFORMATTED' ) )

    elseif ( s_eqi ( text, 'NORMALATQUADRILATERALVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALATTRIANGLEVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALATVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'QUADRILATERALS' ) )

      keyword = 'QUADRILATERALS';

    elseif ( s_eqi ( text, 'REQUIREDEDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'REQUIREDVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'RIDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TANGENTATEDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TANGENTS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TETRAHEDRA' ) || s_eqi ( text, 'TETRAHEDRONS' ) )

      keyword = 'TETRAHEDRONS';

    elseif ( s_eqi ( text, 'TRIANGLES' ) )

      keyword = 'TRIANGLES';

    elseif ( s_eqi ( text, 'VERTICES' ) )

      keyword = 'VERTICES';
%
%  Presumably, numeric data to be processed by keyword.
%
    elseif ( s_eqi ( keyword, 'DIMENSION' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       dim = value;
      end

      keyword = 'NONE';

    elseif ( s_eqi ( keyword, 'EDGES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       edges = value;
      end

      keyword = 'EDGE_VERTEX';
      edge = 0;
      edge_vertex = zeros ( 2, edges );
      edge_label = zeros ( 1, edges );

    elseif ( s_eqi ( keyword, 'EDGE_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d' );

      edge = edge + 1;
      edge_vertex(1:2,edge) = value(1:2);
      edge_label(edge) = value(3);

    elseif ( s_eqi ( keyword, 'HEXAHEDRONS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       hexahedrons = value;
      end

      keyword = 'HEXAHEDRON_VERTEX';
      hexahedron = 0;
      hexahedron_vertex = zeros ( 8, hexahedrons );
      hexahedron_label = zeros ( 1, hexahedrons );

    elseif ( s_eqi ( keyword, 'HEXAHEDRON_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d  %d  %d  %d  %d' );

      hexahedron = hexahedron + 1;
      hexahedron_vertex(1:8,hexahedron) = value(1:8);
      hexahedron_label(hexahedron) = value(9);

    elseif ( s_eqi ( keyword, 'QUADRILATERALS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       quadrilaterals = value;
      end

      keyword = 'QUADRILATERAL_VERTEX';
      quadrilateral = 0;
      quadrilateral_vertex = zeros ( 4, quadrilateral );
      quadrilateral_label = zeros ( 1, quadrilateral );

    elseif ( s_eqi ( keyword, 'QUADRILATERAL_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d' );

      quadrilateral = quadrilateral + 1;
      quadrilateral_vertex(1:4,quadrilateral) = value(1:4);
      quadrilateral_label(quadrilateral) = value(5);

    elseif ( s_eqi ( keyword, 'TETRAHEDRONS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       tetrahedrons = value;
      end

      keyword = 'TETRAHEDRON_VERTEX';
      tetrahedron = 0;
      tetrahedron_vertex = zeros ( 4, tetrahedron );
      tetrahedron_label = zeros ( 1, tetrahedron );

    elseif ( s_eqi ( keyword, 'TETRAHEDRON_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d' );

      tetrahedron = tetrahedron + 1;
      tetrahedron_vertex(1:4,tetrahedron) = value(1:4);
      tetrahedron_label(tetrahedron) = value(5);

    elseif ( s_eqi ( keyword, 'TRIANGLES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       triangles = value;
      end

      keyword = 'TRIANGLE_VERTEX';
      triangle = 0;
      triangle_vertex = zeros ( 3, triangle );
      triangle_label = zeros ( 1, triangle );

    elseif ( s_eqi ( keyword, 'TRIANGLE_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d' );

      triangle = triangle + 1;
      triangle_vertex(1:3,triangle) = value(1:3);
      triangle_label(triangle) = value(4);

    elseif ( s_eqi ( keyword, 'VERTICES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
%       vertices = value;
      end

      keyword = 'VERTEX_COORDINATE';
      vertex = 0;
      vertex_coordinate = zeros ( dim, vertices );
      vertex_label = zeros ( 1, vertices );

    elseif ( s_eqi ( keyword, 'VERTEX_COORDINATE' ) )

      if ( dim == 2 )
        [ value, count ] = sscanf ( text, '%f  %f  %d' );
      elseif ( dim == 3 )
        [ value, count ] = sscanf ( text, '%f  %f  %f  %d' );
      end

      vertex = vertex + 1;
      vertex_coordinate(1:dim,vertex) = value(1:dim);
      vertex_label(vertex) = value(dim+1);

    elseif ( s_eqi ( keyword, 'SKIP' ) )

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'MESH_DATA_READ - Fatal error!\n' );
      fprintf ( 1, '  Could not find keyword while reading line %d:\n', line_num );
      fprintf ( 1, '"%s"\n', text );
      error ( 'MESH_DATA_READ - Fatal error!\n' );

    end
  end
%
%  Close the file.
%
  fclose ( unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read %d lines from "%s".\n', line_num, filename );

  return
end
function [ dim, vertices, edges, triangles, quadrilaterals, tetrahedrons, ...
  hexahedrons ] = mesh_size_read ( filename )

%*****************************************************************************80
%
%% MESH_SIZE_READ reads a MESH file defining a mesh and returns sizes only.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Pascal Frey,
%    MEDIT: An interactive mesh visualization software,
%    Technical Report RT-0253,
%    Institut National de Recherche en Informatique et en Automatique,
%    03 December 2001.
%
%  Parameters:
%
%    Input, string FILENAME, the name of the MESH file.
%
%    Output, integer DIM, the spatial dimension, which should be 2 or 3.
%
%    Output, integer VERTICES, the number of vertices.
%
%    Output, integer EDGES, the number of edges (may be 0).
%
%    Output, integer TRIANGLES, the number of triangles (may be 0).
%
%    Output, integer QUADRILATERALS, the number of quadrilaterals (may be 0).
%
%    Output, integer TETRAHEDRONS, the number of tetrahedrons (may be 0).
%
%    Output, integer HEXAHEDRONS, the number of hexahedrons (may be 0).
%

%
%  Initialize everything to nothing.
%
  dim = 0;
  vertices = 0;
  edges = 0;
  triangles = 0;
  quadrilaterals = 0;
  tetrahedrons = 0;
  hexahedrons = 0;
%
%  Open the file.
%
  unit = fopen ( filename, 'rt' );

  if ( unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_SIZE_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file "%s".\n', filename );
    error ( 'MESH_SIZE_READ - Error!' );
    return
  end
%
%  Read lines til you get alphanumerics and determine a "mode"
%
  line_num = 0;
  keyword = 'NONE';

  while ( 1 )

    text = fgetl ( unit );

    if ( text == -1 )
      break
    end

    line_num = line_num + 1;

    if ( s_len_trim ( text ) == 0 )
      keyword = 'NONE';
      continue
    end

    if ( text(1) == '#' )
      continue
    end
%
%  Remove initial blanks.
%

%
%  Expecting a keyword.
%
        if ( s_eqi ( text, 'CORNERS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'DIMENSION' ) )

      keyword = 'DIMENSION';

    elseif ( s_eqi ( text, 'EDGES' ) )

      keyword = 'EDGES';

    elseif ( s_eqi ( text, 'END' ) )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  END statement encountered.\n' );
      break

    elseif ( s_eqi ( text, 'HEXAHEDRA' ) || s_eqi ( text, 'HEXAHEDRONS' ) )

      keyword = 'HEXAHEDRONS';

    elseif ( s_begin ( text, 'MESHVERSIONFORMATTED' ) )

    elseif ( s_eqi ( text, 'NORMALATQUADRILATERALVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALATTRIANGLEVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALATVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'NORMALS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'QUADRILATERALS' ) )

      keyword = 'QUADRILATERALS';

    elseif ( s_eqi ( text, 'REQUIREDEDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'REQUIREDVERTICES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'RIDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TANGENTATEDGES' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TANGENTS' ) )

      keyword = 'SKIP';

    elseif ( s_eqi ( text, 'TETRAHEDRA' ) || s_eqi ( text, 'TETRAHEDRONS' ) )

      keyword = 'TETRAHEDRONS';

    elseif ( s_eqi ( text, 'TRIANGLES' ) )

      keyword = 'TRIANGLES';

    elseif ( s_eqi ( text, 'VERTICES' ) )

      keyword = 'VERTICES';
%
%  Presumably, numeric data to be processed by keyword.
%
    elseif ( s_eqi ( keyword, 'DIMENSION' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        dim = value;
      end

      keyword = 'NONE';

    elseif ( s_eqi ( keyword, 'EDGES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        edges = value;
      end

      keyword = 'EDGE_VERTEX';

    elseif ( s_eqi ( keyword, 'EDGE_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d' );

    elseif ( s_eqi ( keyword, 'HEXAHEDRONS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        hexahedrons = value;
      end

      keyword = 'HEXAHEDRON_VERTEX';

    elseif ( s_eqi ( keyword, 'HEXAHEDRON_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d  %d  %d  %d  %d' );

    elseif ( s_eqi ( keyword, 'QUADRILATERALS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        quadrilaterals = value;
      end

      keyword = 'QUADRILATERAL_VERTEX';

    elseif ( s_eqi ( keyword, 'QUADRILATERAL_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d' );

    elseif ( s_eqi ( keyword, 'TETRAHEDRONS' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        tetrahedrons = value;
      end

      keyword = 'TETRAHEDRON_VERTEX';

    elseif ( s_eqi ( keyword, 'TETRAHEDRON_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d  %d' );

    elseif ( s_eqi ( keyword, 'TRIANGLES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        triangles = value;
      end

      keyword = 'TRIANGLE_VERTEX';

    elseif ( s_eqi ( keyword, 'TRIANGLE_VERTEX' ) )

      [ value, count ] = sscanf ( text, '%d  %d  %d  %d' );

    elseif ( s_eqi ( keyword, 'VERTICES' ) )

      [ value, count ] = sscanf ( text, '%d' );

      if ( count == 1 )
        vertices = value;
      end

      keyword = 'VERTEX_COORDINATE';

    elseif ( s_eqi ( keyword, 'VERTEX_COORDINATE' ) )

      if ( dim == 2 )
        [ value, count ] = sscanf ( text, '%f  %f  %d' );
      elseif ( dim == 3 )
        [ value, count ] = sscanf ( text, '%f  %f  %f  %d' );
      end

    elseif ( s_eqi ( keyword, 'SKIP' ) )

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'MESH_SIZE_READ - Fatal error!\n' );
      fprintf ( 1, '  Could not find keyword while reading line %d:\n', line_num );
      fprintf ( 1, '"%s"\n', text );
      error ( 'MESH_SIZE_READ - Fatal error!\n' );

    end
  end
%
%  Close the file.
%
  fclose ( unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read %d lines from "%s".\n', line_num, filename );

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

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  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.
%
%  Alternative print statements include:
%
%     fprintf ( output_unit, '  %24.16e', table(i,j) );
%     fprintf ( output_unit, '  %14.6e', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %g', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function value = s_begin ( s1, s2 )

%*****************************************************************************80
%
%% S_BEGIN is TRUE if one string matches the beginning of the other.
%
%  Discussion:
%
%    The strings are compared, ignoring blanks and capitalization.
%
%  Example:
%
%     S1              S2      S_BEGIN
%
%    'Bob'          'BOB'     TRUE
%    '  B  o b '    ' bo b'   TRUE
%    'Bob'          'Bobby'   TRUE
%    'Bobo'         'Bobb'    FALSE
%    ' '            'Bob'     FALSE    (Do not allow a blank to match
%                                       anything but another blank string.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character S1(*), S2(*), the strings to be compared.
%
%    Output, logical S_BEGIN, is TRUE if the strings match up to
%    the end of the shorter string, ignoring case.
%
  len1 = s_len_trim ( s1 );
  len2 = s_len_trim ( s2 );
%
%  If either string is blank, then both must be blank to match.
%  Otherwise, a blank string matches anything, which is not
%  what most people want.
%
  if ( len1 == 0 || len2 == 0 )

    if ( len1 == 0 && len2 == 0 )
      value = 1;
    else
      value = 0;
    end

    return

  end

  i1 = 0;
  i2 = 0;
%
%  Find the next nonblank in S1.
%
  while ( 1 )

    while ( 1 )

      i1 = i1 + 1;

      if ( len1 < i1 )
        value = 1;
        return
      end

      if ( s1(i1) ~= ' ' )
        break
      end

    end
%
%  Find the next nonblank in S2.
%
    while ( 1 )

      i2 = i2 + 1;

      if ( len2 < i2 )
        value = 1;
        return
      end

      if ( s2(i2) ~= ' ' )
        break
      end

    end
%
%  If the characters match, get the next pair.
%
    if ( ~ch_eqi ( s1(i1), s2(i2) ) )
      break
    end

  end

  value = 0;

  return
end
function value = s_eqi ( s1, s2 )

%*****************************************************************************80
%
%% S_EQI is a case insensitive comparison of two strings for equality.
%
%  Example:
%
%    S_EQI ( 'Anjana', 'ANJANA' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, S2, the strings to compare.
%
%    Output, logical VALUE, is TRUE if the strings are equal.
%
  len1 = length ( s1 );
  len2 = length ( s2 );
  lenc = min ( len1, len2 );

  for i = 1 : lenc

    c1 = ch_cap ( s1(i) );
    c2 = ch_cap ( s2(i) );

    if ( c1 ~= c2 )
      value = 0;
      return
    end

  end

  for i = lenc + 1 : len1
    if ( s1(i) ~= ' ' )
      value = 0;
      return
    end
  end

  for i = lenc + 1 : len2
    if ( s2(i) ~= ' ' )
      value = 0;
      return
    end
  end

  value = 1;

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  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 anim = nrsfmXiaoKanade( W, nBasis )
% Compute orthographic non-rigid structure from motion using Xiao-Kanade
%
% Just check the Xiao's CVPR04 paper
%
% USAGE
%  anim = nrsfmXiaoKanade( W, nBasis )
%
% INPUTS
%  W             - [ 2 x nPoint x nFrame ] set of 2D points
%  nBasis        - {Computed } number of bases to use
%
% OUTPUTS
%  anim          - Animation object (help Animation for details)
%
% EXAMPLE
%
% See also COMPUTESMFROMW
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2008-2010 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

nFrame = size( W, 3 ); nPoint = size( W, 2 );
animBest=Animation; animBest.W = W;

W = reshape( permute( W, [ 1 3 2 ] ), [], nPoint );
animBest.t = reshape( mean(W,2), 2, nFrame ); animBest.t(3,:) = 0;
W = bsxfun(@minus,W,mean(W,2));

if nargin<2 || isempty(nBasis)
  [ U S V ] = svd(W,'econ'); S = diag(S); sumVal = S(1); Kd = 1;
  while ( Kd<=length(S) ) && ( sumVal < 0.995*sum(abs(S)) )
    if S(Kd)<0; break; end
    Kd = Kd + 1; sumVal = sumVal + S(Kd);
  end
  W = U(:,1:Kd)*diag(S(1:Kd))*V(:,1:Kd)';
else
  Kd = 3*nBasis;
  % Check if Kd is not smaller (degenerate deformations)
  [ U S V ] = svd(W,'econ');
  S=cumsum(diag(S));
  ind = find(S>=0.995*S(end));
  Kd=min(Kd,ind(1));
end

%%%%%%%%%%%%%%%%%%%%%% Determine K3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K3 = floor(Kd/3);
% Basis constraints
WTildeOri = [ W(1:2:end,:) W(2:2:end,:) ];

while 1
  % Perform random search for 5 seconds
  t0 = clock; mini = Inf;
  while etime( clock, t0 )<5
    indWTmp = randSample(nFrame,K3);
    tmp = cond( WTildeOri( indWTmp, : ) );
    if tmp < mini; mini = tmp; indW = indWTmp; end
  end

  if mini>1e2 && K3>0; K3=K3-1; else break; end
end

%%%%%%%%%%%%%%%%%%%%%% Deal with the K3 basis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute initial variables
[ MTilde S disc ] = svd( W );
MTilde = MTilde( :, 1:Kd )*S(1:Kd,1:Kd); gj3=zeros(Kd,3,K3);

% Rotation constraints
Qi = sdpvar( Kd, Kd );
F = set( Qi >= 0 );
% vect(AXB)=kron(B',A)*vect(X)
kron1=zeros(nFrame,Kd^2); kron2=zeros(nFrame,Kd^2);
for i=1:nFrame
  kron1(i,:)=kron(MTilde(2*i-1,:),MTilde(2*i-1,:)) - ...
    kron(MTilde(2*i,:),MTilde(2*i,:));
  kron2(i,:)=kron(MTilde(2*i,:),MTilde(2*i-1,:));
end
objRotation=abs(kron1*Qi(:))+abs(kron2*Qi(:));

% Solve for Qk using all the constraints
kk = 1;
for i = indW
  nc = 2*nFrame + 1;
  objBasis = 0;
  disp('Computing the basis constraints');
  for m = indW
    if m==i; interN = m; else interN = 1 : nFrame; end
    for n = interN % could be sampled
      if m==n && m==i
        objBasis = objBasis +abs(MTilde(2*m-1,:)*Qi*MTilde(2*n-1,:)'-1) ...
          + abs( MTilde(2*m,:)*Qi*MTilde(2*n,:)' - 1);
      else
        objBasis = objBasis +abs(MTilde(2*m-1,:)*Qi*MTilde(2*n-1,:)') ...
          + abs( MTilde(2*m,:)*Qi*MTilde(2*n,:)' );
      end

      objBasis = objBasis + abs( MTilde(2*m-1,:)*Qi*MTilde(2*n,:)' ) + ...
        abs( MTilde(2*m,:)*Qi*MTilde(2*n-1,:)' );
    end
  end

  % Sample the constraints
  disp('Solving the SDP');
  diagno = solvesdp( F,objRotation+objBasis,sdpsettings('solver',...
    'sdpa,csdp,sedumi,*','dualize',1,'verbose',0));

  QiTmp = double(Qi);

  [ U S disc ] = svd( QiTmp ); gj3(:,:,kk) = U(:,1:3)*sqrt(S(1:3,1:3));
  kk = kk + 1;
end

%  [U S V]=svd(constr,'econ');
%  rank(constr)
%  size(null(constr(1:2*nFrame,:)))
%  size(constr)
%  tmp=diag(S)';
%  tmp(end-3:end)
%  size(null(constr),2)
%  return
%%%%%%%%%%%%%%%%%%%%%% Get K2 and K1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The explanations given by Xiao04 seems wrong to determine K1 and K2
%% So, let's do it the lazy way: let's try all possible combinations of K1
%% and K2
errBest = Inf;
fprintf('\nSo far, Kd=%d, K3=%d\n', Kd, K3 );
for K2 = floor((Kd-3*K3)/2) : -1 : 0
  K1 = Kd - 3*K3 - 2*K2;
  K = K3+K2+K1;

  fprintf('\nTrying, Kd=%d, K1=%d, K2=%d, K3=%d\n', Kd, K1, K2, K3 );

  %%%%%%%%%%%%%%%%%%%%%% Get the sets of rotations %%%%%%%%%%%%%%%%%%%%%%%%
  G = reshape( gj3, Kd, [] );
  M = MTilde*G; M = M/max(abs(M(:)))*K3; % Just for numerical stability
  R = zeros( 2*nFrame, 3*K3 ); l = zeros( K, nFrame ); thres = 0.001;
  rotIsBad = zeros( K3, nFrame );
  for i = 1 : nFrame
    for k = 1 : K3
      tmp = M(2*i-1:2*i,3*k-2:3*k);
      l(k,i) = norm(tmp);

      if abs( l(k,i) )>thres; R(2*i-1:2*i,3*k-2:3*k) = tmp/l(k,i);
      else rotIsBad(k,i)=1;
      end
    end
  end

  % Rectify R (not really mentioned by Xiao)
  for k=2:K3
    for i = 2:nFrame
      goodFrame = find( sum( rotIsBad([1 k], 1:i), 1 )==0);
      goodFrame = sort( [ 2*goodFrame 2*goodFrame-1 ] );
      if length(goodFrame)<4; continue; end
      R1 = R( goodFrame, 1:3 );
      R2 = R( goodFrame, 3*k-2:3*k );
      R3 = R2; R3(end-1:end,:) = -R3(end-1:end,:);

      if norm( R2*(R2\R1)-R1, 'fro' )>norm( R3*(R3\R1)-R1, 'fro' )
        R( goodFrame(end-1:end), 3*k-2:3*k ) = ...
          -R( goodFrame(end-1:end), 3*k-2:3*k );
      end
    end

    goodFrame = find( sum( rotIsBad([1 k], 1:nFrame), 1 )==0);
    goodFrame = sort( [ 2*goodFrame 2*goodFrame-1 ] );
    R1 = R( goodFrame, 1:3 );
    R2 = R( goodFrame, 3*k-2:3*k );
    tmp=R2\R1;

    R( :, 3*k-2:3*k ) = R( :, 3*k-2:3*k )*tmp;
    M( :, 3*k-2:3*k ) = M( :, 3*k-2:3*k )*tmp;
  end

  % Get all the rotations and coefficients in K3 (if non-degenerate)
  RTot = zeros(3,3,nFrame);
  for i = 1 : nFrame
    % Get the average rotation matrix
    tmp=find(~rotIsBad(:,i)');
    if isempty(tmp); continue; end
    RTmp=zeros(3,3,length(tmp));

    qTot=cell(1,2);
    for k=1:2
      kk=1;
      for j=tmp
        RTmp(:,:,kk)=rotationMatrix( R(2*i-1:2*i,3*j-2:3*j)*(2*k-3) );
        kk=kk+1;
      end

      q=quaternion(RTmp);
      for j=2:size(q,2)
        if norm(-q(:,j)-q(:,1))<norm(q(:,j)-q(:,1)); q(:,j)=-q(:,j); end
      end
      qTot{k}=mean(q,2);
    end

    q=qTot{1};
    if i>=2
      q0 = quaternion( RTot(:,:,i-1) );
      if min( norm(qTot{2}-q0), norm(-qTot{2}-q0) )<...
          min( norm(qTot{1}-q0), norm(-qTot{1}-q0) )
        q=qTot{2};
      end
    end

    RTot(:,:,i) = quaternion(q);

    % Recover the coefficients
    l(1:K3,i)=(reshape(RTot(1:2,:,i),[],1)\reshape(M(2*i-1:2*i,1:3*K3),6,K3))';

    % Compute MTilde
    MTilde(2*i-1:2*i,1:3*K3)=kron(l(1:K3,i)',RTot(1:2,:,i));
  end
  for k = 1 : K3
    gj3(:,:,k) = 0; gj3(3*k-2:3*k,:,k) = eye(3);
  end
  R = RTot;
  rotIsUnknown = sum(abs(l(1:K3,:))<thres,1)==K3;
  if sum(~rotIsUnknown)==0
    [disc indMax]=max(sum(l(1:K3,:).^2,1)); rotIsUnknown(indMax)=1;
  end
  goodFrame = find( ~rotIsUnknown );
  tmp=sort([2*goodFrame,2*goodFrame-1]);

  % Get the MTilde where l(1:K3,i) is close to 0
  M=MTilde(tmp,:); B=M\W(tmp,:); M2=W/B;
  for i=find(rotIsUnknown)
    MTilde(2*i-1:2*i,:)=M2(2*i-1:2*i,:);
  end

  %   find(rotIsUnknown)
  %   plot(l')
  %     B=MTilde\W;
  %   norm(MTilde*B-W,'fro')
  %   pause

  % The first 3*K3 columns of MTilde and MHat are now identical
  % We now need to deal with the last columns

  %%%%%%%%%%%%%%%%%%%%%% Find gj and rj %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % We here perform a full optimization and not the hacky alternate
  % optimization from the paper
  MNull = zeros( 2*K3, Kd );
  for i = 1 : K3
    MNull(2*i-1:2*i,:) = MTilde( 2*indW(i)-1:2*indW(i), : );
  end
  [ U S V ] = svd( MNull, 'econ' );
  [ V disc ] =rq(V);
  gj2 = reshape( V(:,end-(2*K2+K1)+1:end-K1), [], 2, K2);rj2=zeros(3,2,K2);
  gj1 = V(:,end-K1+1:end); rj1=zeros(3,K1);

  % Perform iterations for K2 bases
  opt = optimset( 'Display', 'on', 'GradObj', 'off', 'Hessian', 'off',...
    'LargeScale', 'off' );
  for j=1:K2
    % Get the best rj2
    ATmp=zeros(4*length(goodFrame),3);
    for kk=1:2
      n = 1;
      for m=goodFrame
        for k=1:2
          ATmp(n,:)=gj2(:,k,j)'*MTilde(2*m-1,:)'*R(2,:,m) - ...
            gj2(:,k,j)'*MTilde(2*m,:)'*R(1,:,m);
          n = n + 1;
        end
      end

      % Force the rj2 to be independent
      zeroStart=2+kk;
      [ U S V ] = svd( ATmp(:,1:zeroStart-1), 'econ' );
      rj2(1:zeroStart-1,kk,j) = V(:,end)/V(end,end);
      rj2(zeroStart:end,kk,j) = 0;
    end

    % Perform the full optimization
    onePos=3*K3+1+2*(j-1);
    x=[gj2(1:onePos-1,1,j); gj2(1:onePos,2,j); rj2(1:2,1,j); rj2(1:2,2,j)];
    MTildeTmp= MTilde(:,1:onePos+1); MNullTmp=MNull(:,1:onePos+1);
    x = fminunc( @(X)( optimizeGj2( X,goodFrame, MTildeTmp, ...
      R, MNullTmp ) ), x,opt);

    % Save the best results
    rj2(:,:,j)=[x(end-3:end-2) x(end-1:end); 0 1 ];
    gj2(:,:,j)=0;
    gj2(1:onePos,1,j) = [ x(1:onePos-1); 1 ];
    gj2(1:onePos+1,2,j) = [ x(1:onePos); 1 ];
  end

  % Perform iterations for K1 bases
  for j=1:K1
    % Get the best rj1
    ATmp=zeros(length(goodFrame),3);
    for m=goodFrame
      ATmp(m,:) = gj1(:,j)'*( MTilde( 2*m-1 , : )'*R(2,:,m) - ...
        MTilde( 2*m , : )'*R(1,:,m) );
    end

    [ U S V ] = svd( ATmp, 'econ' );
    rj1(:,j) = V(:,end);

    % Perform the full optimization
    onePos=3*K3+2*K2+j;
    x=[gj1(1:onePos-1,j); rj1(:,j) ];
    MTildeTmp= MTilde(:,1:onePos); MNullTmp=MNull(:,1:onePos);
    x = fminunc( @(X)( optimizeGj1( X,goodFrame, MTildeTmp, ...
      R, MNullTmp ) ), x,opt);

    % Save the best results
    rj1(:,j)=x(end-2:end);
    gj1(:,j)=0; gj1(1:onePos,j) = [ x(1:onePos-1); 1 ];
  end

  % % quality
  % err=zeros(1,length(goodFrame));
  % for j = 1:K2
  %   for m=goodFrame
  %     err(m)=err(m)+abs(gj2(:,1,j)'*( MTilde( 2*m-1 , : )'*R(2,:,m) - ...
  %       MTilde( 2*m , : )'*R(1,:,m) )*rj2(:,1,j)) ...
  %        + abs(gj2(:,2,j)'*( MTilde( 2*m-1 , : )'*R(2,:,m) - ...
  %       MTilde( 2*m , : )'*R(1,:,m) )*rj2(:,2,j));
  %   end
  % end
  %       norm(err)
  % %       norm(MNull*gj1(:,j))
  %       plot(err)
  %
  % pause

  %%%%%%%%%%%%%%%%%%%%%% Recover the coefficients for K2 K1 %%%%%%%%%%%%%%%
  % Rectify the right columns of MTilde
  for m=1:nFrame
    for i = 1 : K2
      MTilde(2*m-1:2*m,3*K3+2*i-1:3*K3+2*i)=MTilde(2*m-1:2*m,:)*gj2(:,:,i);
    end
    for i = 1 : K1
      MTilde(2*m-1:2*m,3*K3+2*K2+i)=MTilde(2*m-1:2*m,:)*gj1(:,i);
    end
  end
  %   l(K3+i,m)*R(1:2,:,m)*rj2(:,:,i)=
  %                   MTilde(2*m-1:2*m,3*K3+2*i-1:3*K3+2*i)*RAmbiguityi
  %   l(K3+i,m)*R(1:2,:,m)*rj1(:,i)=
  %                   MTilde(2*m-1:2*m,3*K3+2*K2+i)
  for m=goodFrame
    for i = 1 : K2
      tmp1=R(1:2,:,m)*rj2(:,:,i); tmp1=tmp1*tmp1';
      tmp2=MTilde(2*m-1:2*m,3*K3+2*i-1:3*K3+2*i); tmp2=tmp2*tmp2';
      l(K3+i,m) = sqrt( tmp1(:)\tmp2(:) );
    end
    for i = 1 : K1
      tmp1=R(1:2,:,m)*rj1(:,i);
      tmp2=MTilde(2*m-1:2*m,3*K3+2*K2+i);
      l(K3+K2+i,m) = tmp1\tmp2;
    end
  end

  %   B=MTilde\W;
  %   norm(MTilde*B-W,'fro')
  %   plot(l')
  %   pause

  % Change the signs of the coefficients
  for i = 1 : K2
    A=zeros(2*nFrame,2);
    A(1:2,:)=l(K3+i,goodFrame(1))*R(1:2,:,goodFrame(1))*rj2(:,:,i);

    for j=2:length(goodFrame)
      m=goodFrame(j); A(2*j-1:2*j,:)=l(K3+i,m)*R(1:2,:,m)*rj2(:,:,i);
      tmp=goodFrame(1:j); tmp=sort( [ 2*tmp-1, 2*tmp ] );
      MTildeLoc=MTilde(tmp,3*K3+2*i-1:3*K3+2*i);

      RAmb1=MTildeLoc\A(1:2*j,:);
      err1=norm(MTildeLoc*RAmb1-A(1:2*j,:),'fro');

      A(2*j-1:2*j,:)=-A(2*j-1:2*j,:);
      RAmb2=MTildeLoc\A(1:2*j,:);
      err2=norm(MTildeLoc*RAmb2-A(1:2*j,:),'fro');

      if err2<err1; l(K3+i,m)=-l(K3+i,m);
      else A(2*j-1:2*j,:)=-A(2*j-1:2*j,:);
      end
    end
  end

  B=MTilde\W;
  %   norm(MTilde*B-W,'fro')
  %   plot(l')
  %   pause

  %%%%%%%%%%%%%%%%%%%%%% Recover the shape basis %%%%%%%%%%%%%%%%%%%%%%%%%%
  M=zeros(2*length(goodFrame),3*K);
  for i=1:length(goodFrame)
    M(2*i-1:2*i,:)=kron(l(:,goodFrame(i))',R(1:2,:,goodFrame(i)));
  end

  tmp=sort([2*goodFrame,2*goodFrame-1]);
  B=M\W(tmp,:);
  %   norm(M*B-W(tmp,:),'fro')
  %   tmp=sum((M*B-W(tmp,:)).^2,2);
  %   plot( tmp(1:2:end)+tmp(2:2:end) )
  %   'coiin'
  %   pause

  % Recover the full basis
  SBasis = zeros( 3, nPoint, K );
  for i = 1 : K; SBasis(:,:,i) = B(3*i-2:3*i,:); end

  %%%%%%%%%%%%%%%%%%%%%% Recover the missing coefficients %%%%%%%%%%%%%%%%%
  opt = optimset( 'Display', 'off', 'GradObj', 'off', 'Hessian', 'off',...
    'LargeScale', 'off' );

  while 1
    for m=1:nFrame
      if rotIsUnknown(m)
        if m==1 || rotIsUnknown(m-1)
          if m==nFrame; continue
          else if ~rotIsUnknown(m+1); iIni=m+1; else continue; end
          end
        else iIni=m-1;
        end

        % Optimize over the shape coefficients/rotation
        coeff=[ l(:,iIni); quaternion( R(:,:,iIni) ) ];
        Wm=W(2*m-1:2*m,:);

        [ c res ] = fminunc( @(X)( optimRl( X, K,B,Wm ) ), coeff, opt );

        l(:,m)=c(1:K); R(:,:,m)=quaternion(c(K+1:end));
        rotIsUnknown(m)=0;
      else
        if m>=2
          %           plot(quaternion(R)')
          %           pause
          q0 = quaternion( R(:,:,m-1) );
          q1 = quaternion( R(:,:,m) );
          q2 = quaternion( rotationMatrix( -R(1:2,:,m) ) );
          if min( norm(q2-q0), norm(-q2-q0) )<min(norm(q1-q0),norm(-q1-q0))
            l(:,m)=-l(:,m); R(:,:,m)=rotationMatrix(-R(1:2,:,m));
          end
        end
      end
    end
    if ~rotIsUnknown(1); break; end
  end

  %%%%%%%%%%%%%%%%%%%%%% Finalize the anim object %%%%%%%%%%%%%%%%%%%%%%%%%
  anim=Animation('SBasis', SBasis, 'l', l, 'R', R, 'isProj', false, ...
    't', animBest.t, 'W', animBest.W );
  err = anim.computeError(); err=err(1);
  fprintf( 'Reprojection error: %f \n', err );

  if err<errBest; errBest=err; animBest=anim; end
end
anim=animBest;

fprintf( 'Best reprojection error: %f \n', errBest );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res = optimRl(X,K,B,Wm)
lm=X(1:K);
Rm=quaternion(X(K+1:end)/norm(X(K+1:end)));
res = norm(kron(lm',Rm(1:2,:))*B-Wm,'fro')^2;
end

function val = optimizeGj2( x, goodFrame, MTilde, R, MNull )
% gjk' MTilde_{2m-1}' R_m2 rjkk - gjk' MTilde_2m' R_m1 rjkk = 0,
% k and kk in {1 2}
% rj can be multiplied by any 2x2
rj2=[x(end-3:end-2) x(end-1:end); 0 1 ]; rj2(:,1)=rj2(:,1)/norm(rj2(:,1));
n=(length(x)-4-1)/2;
gj2=[ [ x(1:n); 1; 0 ] [ x(n+1:2*n+1); 1 ] ];
n=n+2;

n=1;
val=zeros(1,length(goodFrame));
A=gj2'*MTilde';
for m=goodFrame
  val(n)=norm( ( A(:,2*m-1,:)*R(2,:,m) - A(:,2*m)*R(1,:,m) )*rj2, 'fro' );
  n=n+1;
end
val=norm([ val reshape(MNull*gj2,[],1)' ]);
end

function val = optimizeGj1( x, goodFrame, MTilde, R, MNull )
% gj' ( MTilde_{2m-1} R_m2 - MTilde_2m R_m1 ) rj = 0
gj1=[ x(1:end-3); 1 ];
rj1=x(end-2:end); rj1=rj1/norm(rj1);

n=1;
A=gj1'*MTilde';
B=zeros(length(goodFrame),3);
for m=goodFrame
  B(n,:)=A(:,2*m-1)*R(2,:,m) - A(:,2*m)*R(1,:,m);
  n=n+1;
end
val=norm([ (B*rj1); reshape(MNull*gj1,[],1) ]);

end

author: vrabaud


text:
function [ output_meta ] = meta2sicd_acft( acft_struct )
%META2SICD_ACFT Converts metadata stored in the ACFTA/B NITF TRE into a
% SICD metadata format
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

M2F=3.28083989501312; % Set meters-to-feet conversion constant

if ~isempty(acft_struct.SENSOR_ID)
    output_meta.CollectionInfo.CollectorName=acft_struct.SENSOR_ID;
end

% This NITF value is actually the ground plane spacing.  Should actually be
% multiplied by cos of graze angle for SICD slant plane spacing.
output_meta.Grid.Col.SS=acft_struct.ROW_SPACING;
% This NITF value is actually the ground plane spacing.  Should actually be
% multiplied by cos of tilt angle for SICD slant plane spacing.
output_meta.Grid.Row.SS=acft_struct.COL_SPACING;

if ~isfield(acft_struct,'ROW_SPACING_UNITS')||acft_struct.ROW_SPACING_UNITS=='f'
    output_meta.Grid.Col.SS=output_meta.Grid.Col.SS/M2F;
end
if ~isfield(acft_struct,'COL_SPACING_UNITS')||acft_struct.COL_SPACING_UNITS=='f'
    output_meta.Grid.Row.SS=output_meta.Grid.Row.SS/M2F;
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function net = initModel(param)
caffe.reset_all();
if exist(param.modelFile, 'file') == 0
  fprintf('%s does not exist. Start downloading ...\n', param.modelFile);
  downloadModel(param.modelName);
end
if ~exist(param.protoFile,'file')
  error('%s does not exist.', param.protoFile);
end

if param.useGPU
  fprintf('Using GPU Mode\n');
  caffe.set_mode_gpu();
  caffe.set_device(param.GPUID);
else
  fprintf('Using CPU Mode\n');
  caffe.set_mode_cpu;
end

net = caffe.Net(param.protoFile, param.modelFile, 'test');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [B,Theta_diag,C_diag,Alpha_diag,Fit_diag,...
    C_offdiag,Theta_offdiag,Alpha_offdiag,Fit_offdiag phi_t]=smooth_data(C,K,Q,nloop_diag,nloop_offdiag)

%% This function is a wrapper for the smooth_mcmc function, which does the actual smoothing
%% This was written to smoot the diagonal and off-diagonal conectivity time series separately
%% and to save the coefficents and basis functions from each

N=size(C,1);
M_i=zeros(N,1);
for i=1:N
    M_i(i)=size(C{i},1);
end

%% Get smoothed time-varying connectivity coefficients

% smooth diagonal connectivities
C_diag=cell(sum(M_i),1);
ind=0;
for i=1:N
    for j=1:M_i(i)
        ind=ind+1;
        C_diag{ind}=squeeze(C{i}(j,j,:));
    end
end
[ALPHA ALPHA_BAR THETA FIT phi_t]=smooth_mcmc(C_diag,K,Q,nloop_diag);
first=round(nloop_diag/2); last=nloop_diag+1;
Alpha_diag=mean(ALPHA(:,:,first:last),3);
Theta_diag=mean(THETA(:,:,first:last),3);
Fit_diag=phi_t'*Theta_diag*Alpha_diag;

B=cell(N,1);
ind=0;
for i=1:N
    B{i}=cell(M_i(i));
    for j=1:M_i(i)
        ind=ind+1;
        B{i}{j,j}=Alpha_diag(:,ind)';
    end
end

% smooth off-diagonal connectivities
M_offdiag=M_i(1)^2-M_i(1);
for i=2:N
    M_offdiag=M_offdiag+M_i(i)^2-M_i(i);
end

C_offdiag=cell(M_offdiag,1);
ind=0;
for i=1:N
    for j1=1:M_i(i)
        for j2=1:M_i(i)
            if(not(j1==j2))
                ind=ind+1;
                C_offdiag{ind}=squeeze(C{i}(j1,j2,:));
            end
        end
    end
end
[ALPHA ALPHA_BAR THETA FIT phi_t]=smooth_mcmc(C_offdiag,K,Q,nloop_offdiag);
first=round(nloop_offdiag/2); last=nloop_offdiag+1;
Alpha_offdiag=mean(ALPHA(:,:,first:last),3);
Theta_offdiag=mean(THETA(:,:,first:last),3);
Fit_offdiag=phi_t'*Theta_offdiag*Alpha_offdiag;

ind=0;
for i=1:N
    for j1=1:M_i(i)
        for j2=1:M_i(i)
            if not(j1==j2)
                ind=ind+1;
                B{i}{j1,j2}=Alpha_offdiag(:,ind)';
            end
        end
    end
end

repo: REST


text:
% ESTIMATION OF ZERO COUPON YIELDS FROM COUPON BOND PRICES USING THE
% NELSON-SIEGEL AND SVENSSON MODEL
%
% Model = 'NS' (Nelson-Siegel) | 'Svensson'
%
% Optimization Structure
%
% .Method = 'price' | 'ytm'
%  'price' fits (weighted)) cupon bond prices
%  'ytm' fits yields-to-maturity
%
% .Weights  = 'MD' | 'LA' | empty (Weights are irrelevant for the YTM method.)
%  'MD' weights use modified duration of bonds
%  'LA' weights work as the correct linear approaximation of ytm;
%  LA weights equal MD weights divideded by observed cupon bond prices.
%  Note that weights are irrelevant for the YTM method.
%  Further note that short-rates (LIBOR rates) are underweighted, see the
%  NSobjP.m for details and "The Czech Treasury Yield Curve from 1999 to the
%  Present" for a discussion.
%
% .Algorithm = 'lsqnonlin' | 'fminsearch'
%  Supports two standard Matlab optimization routines. The nonlinear least
%  squares, lsqnonlin, which is tha part of Optimization Tbx is recommended.
%
% .Multistart = 'yes' | 'no'
%  The estimation algorothm creates a grid of starting values for optimization, which
%  gives better hope for finding global optima. Recommended.
%
%
% Kamil Kladivko
% email: kladivko@gmail.com
%
% Cite as:
% Kladivko Kamil (2010). The Czech Treasury Yield Curve from 1999 to the Present,
% Czech Journal of Economics and Finance, 60(4): 307-335
% =========================================================================

clc
close all
clear variables

% 1) Set model
Model = 'NS';
Optimization.Method = 'price';
Optimization.Weights = 'LA';
Optimization.Algorithm = 'lsqnonlin';
Optimization.MultiStart = 'yes';

% 2) Read data from EXCEL
% The EXCEL file contains Czech Governemnt bond prices and PRIBOR rates (short rates) from March 2, 2007
% There were 13 Czech Government bonds available on March 2, 2007. Furter, 4 PRIBOR rates are used
% to fix the short-end of the yield curve. See Section 5 of "The Czech Treasury Yield Curve from 1999 to the Present"

[B Btext] = xlsread('bonddata.xls', 'data', 'B6:F18');
[junk Settle] = xlsread('bonddata.xls', 'data', 'C2:C2');
[SR SRtext] = xlsread('bonddata.xls', 'data', 'J6:L9');

% The DATE format in EXCEL depends on the "national settings". Needs to be
% set accordingly!
Bonds.Issue = datenum(Btext(:,1), 'dd.mm.yyyy');
Bonds.Maturity = datenum(Btext(:,2), 'dd.mm.yyyy');
Bonds.Settle = datenum(Settle, 'dd.mm.yyyy');
Bonds.Prices = B(:,2);
% Make sure that bond ID (isseu number) was read as text! If not adjust
% either this code or EXCEL.
Bonds.Emission = Btext(:,5);
Bonds.Coupon = B(:,1)./100;
% Assumes that coupon is paid once a year. The first cupon can be payed
% later than after one year. See czbondcf.m and czbondfuturecf.m for
% details
Bonds.Notional = 100;
Bonds.Basis = 1;
% Bonds.Basis = '1' | '2'
% Bond.Basis is Day-Count basis for calculation coupon cash flows
% 1: 30E/360 (default), 2: Actual/360
% Go to czbondcf.m and czbondfuturecf.m to adjust for your basis

ShortRates.IR = SR(:,1)'./100; % This must be as row vector opposed to the bond data
if ~isempty(ShortRates), ShorRates.IR = log(1+ShortRates.IR); end % Transform SR to continuous compounding
ShortRates.TimeFactor = SR(:,2)';
ShortRates.Emission = SRtext';

% If Short Rates are not to be used then supply "ShotRates = []" into the
% estimation routine below.
%ShortRates = [];

% 3) Nelson-Siegel Parameters Estimation (calls the main estimation routine)
[Params Fval MaxTime2Mat ParamsMS FvalMS] = NSest(Bonds, ShortRates, Model, Optimization);

% 4) Estimation Error;
IncludeSRinError = 'no';
[PriceE YTME] = NSerror(Bonds, ShortRates, Model, Params, IncludeSRinError);

% 5) Plot
Plots1D(Params, Model, YTME, ShortRates);

author: Sable

repo: mcbench-benchmarks


text:
%ISOLATE_AXES Isolate the specified axes in a figure on their own
%
% Examples:
%   fh = isolate_axes(ah)
%   fh = isolate_axes(ah, vis)
%
% This function will create a new figure containing the axes specified, and
% also their associated legends and colorbars. The axes specified must all
% be in the same figure, but they will generally only be a subset of the
% axes in the figure.
%
% IN:
%    ah - An array of axes handles, which must come from the same figure.
%    vis - A boolean indicating whether the new figure should be visible.
%          Default: false.
%
% OUT:
%    fh - The handle of the created figure.

% Copyright (C) Oliver Woodford 2011

% Thank you to Rosella Blatt for reporting a bug to do with axes in GUIs

function fh = isolate_axes(ah, vis)
% Make sure we have an array of handles
if ~all(ishandle(ah))
    error('ah must be an array of handles');
end
% Check that the handles are all for axes, and are all in the same figure
fh = ancestor(ah(1), 'figure');
nAx = numel(ah);
for a = 1:nAx
    if ~strcmp(get(ah(a), 'Type'), 'axes')
        error('All handles must be axes handles.');
    end
    if ~isequal(ancestor(ah(a), 'figure'), fh)
        error('Axes must all come from the same figure.');
    end
end
% Tag the axes so we can find them in the copy
old_tag = get(ah, 'Tag');
if nAx == 1
    old_tag = {old_tag};
end
set(ah, 'Tag', 'ObjectToCopy');
% Create a new figure exactly the same as the old one
fh = copyfig(fh); %copyobj(fh, 0);
if nargin < 2 || ~vis
    set(fh, 'Visible', 'off');
end
% Reset the axes tags
for a = 1:nAx
    set(ah(a), 'Tag', old_tag{a});
end
% Find the objects to save
ah = findall(fh, 'Tag', 'ObjectToCopy');
if numel(ah) ~= nAx
    close(fh);
    error('Incorrect number of axes found.');
end
% Set the axes tags to what they should be
for a = 1:nAx
    set(ah(a), 'Tag', old_tag{a});
end
% Keep any legends and colorbars which overlap the subplots
lh = findall(fh, 'Type', 'axes', '-and', {'Tag', 'legend', '-or', 'Tag', 'Colorbar'});
nLeg = numel(lh);
if nLeg > 0
    ax_pos = get(ah, 'OuterPosition');
    if nAx > 1
        ax_pos = cell2mat(ax_pos(:));
    end
    ax_pos(:,3:4) = ax_pos(:,3:4) + ax_pos(:,1:2);
    leg_pos = get(lh, 'OuterPosition');
    if nLeg > 1;
        leg_pos = cell2mat(leg_pos);
    end
    leg_pos(:,3:4) = leg_pos(:,3:4) + leg_pos(:,1:2);
    for a = 1:nAx
            % Overlap test
            ah = [ah; lh(leg_pos(:,1) < ax_pos(a,3) & leg_pos(:,2) < ax_pos(a,4) &...
                         leg_pos(:,3) > ax_pos(a,1) & leg_pos(:,4) > ax_pos(a,2))];
    end
end
% Get all the objects in the figure
axs = findall(fh);
% Delete everything except for the input axes and associated items
delete(axs(~ismember(axs, [ah; allchildren(ah); allancestors(ah)])));
return

function ah = allchildren(ah)
ah = allchild(ah);
if iscell(ah)
    ah = cell2mat(ah);
end
ah = ah(:);
return

function ph = allancestors(ah)
ph = [];
for a = 1:numel(ah)
    h = get(ah(a), 'parent');
    while h ~= 0
        ph = [ph; h];
        h = get(h, 'parent');
    end
end
return

function fh = copyfig(fh)
% Is there a legend?
if isempty(findobj(fh, 'Type', 'axes', 'Tag', 'legend'))
    % Safe to copy using copyobj
    fh = copyobj(fh, 0);
else
    % copyobj will change the figure, so save and then load it instead
    tmp_nam = [tempname '.fig'];
    hgsave(fh, tmp_nam);
    fh = hgload(tmp_nam);
    delete(tmp_nam);
end
return

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function xyz = rtp_to_xyz ( r, theta, phi )

%*****************************************************************************80
%
%% RTP_TO_XYZ converts (R,Theta,Phi) to (X,Y,Z) coordinates.
%
%  Discussion:
%
%    R measures the distance of the point to the origin.
%
%    Theta measures the "longitude" of the point, between 0 and 2 PI.
%
%    PHI measures the angle from the "north pole", between 0 and PI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, THETA, PHI, the radius, longitude, and
%    declination of a point.
%
%    Output, real XYZ(3,1), the corresponding Cartesian coordinates.
%
  xyz = zeros ( 3, 1 );

  xyz(1,1) = r * cos ( theta ) * sin ( phi );
  xyz(2,1) = r * sin ( theta ) * sin ( phi );
  xyz(3,1) = r *                 cos ( phi );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% EKF/UKF toolbox for Matlab 7.x
% Version 1.3, August 12, 2011
%
% Copyright (C) 2005-2011 Simo S�rkk�, <simo.sarkka@hut.fi>
%               2007-2011 Jouni Hartikainen <jmjharti@cc.hut.fi>
%               2010-2011 Arno Solin <arno.solin@tkk.fi>
% History:
%   12.08.2011 JH & AS & SS Updated to version 1.3
%   04.09.2007 JH & SS Updated for version 1.1
%   06.08.2007 JH Updated for version 1.0
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
%
%
% Kalman filtering
%   KF_PREDICT    Perform Kalman Filter prediction step
%   KF_UPDATE     Kalman Filter update step
%   KF_LHOOD      Kalman Filter measurement likelihood
%   RTS_SMOOTH    Rauch-Tung-Striebel Smoother
%   TF_SMOOTH     Smoother based on combination of two Kalman filters
%
% Extended Kalman filtering
%   EKF_PREDICT1  1st order Extended Kalman Filter prediction step
%   EKF_UPDATE1   1st order Extended Kalman Filter update step
%   EKF_PREDICT2  2nd order Extended Kalman Filter prediction step
%   EKF_UPDATE2   2nd order Extended Kalman Filter update step
%   ERTS_SMOOTH1  1st order Extended RTS Smoother
%   ETF_SMOOTH1   Smoother based on two 1. order extended Kalman filters
%
% Nonlinear transform based filtering
%   UT_WEIGHTS    Generate weights for sigma points using the summation form
%   UT_MWEIGTS    Generate weights for sigma points using the matrix form
%   UT_SIGMAS     Generate Sigma Points for Unscented Transformation
%   UT_TRANSFORM  Makes the Unscented Transformation of x and y
%   UKF_PREDICT1  Nonaugmented UKF prediction step
%   UKF_UPDATE1   Nonaugmented UKF update step
%   UKF_PREDICT2  Augmented (state and process noise) UKF prediction step
%   UKF_UPDATE2   Augmented (state and measurement noise) UKF update step
%   UKF_PREDICT3  Augmented (state, process and measurement noise) UKF prediction step
%   UKF_UPDATE3   Augmented (state, process and measurement noise) UKF update step
%   URTS_SMOOTH1  Nonaugmented unscented RTS-smoother
%   URTS_SMOOTH2  Augmented unscented RTS-smoother
%   UTF_SMOOTH    Smoother based on combination of two unscented Kalman filters
%   GH_TRANSFORM  Gauss-Hermite transform of random variables
%   GHKF_PREDICT  Gauss-Hermite Kalman filter prediction step
%   GHKF_UPDATE   Gauss-Hermite Kalman filter update step
%   GHRTS_SMOOTH  Additive form Gauss-Hermite Rauch-Tung-Striebel smoother
%   CKF_TRANSFORM Cubature Kalman filter transform of random variables
%   CKF_PREDICT   Cubature Kalman filter prediction step
%   CKF_UPDATE    Cubature Kalman filter update step
%   CRTS_SMOOTH - Additive form cubature Rauch-Tung-Striebel smoother
%
% Multiple Model Filtering
%   IMM_PREDICT   IMM filter prediction step
%   IMM_UPDATE    IMM filter update step
%   IMM_SMOOTH    IMM smoothing
%   EIMM_PREDICT  IMM-EKF filter prediction step
%   EIMM_UPDATE   IMM-EKF filter update step
%   EIMM_SMOOTH   IMM-EKF smoothing
%   UIMM_PREDICT  IMM-UKF filter prediction step
%   UIMM_UPDATE   IMM-UKF filter update step
%   UIMM_SMOOTH   IMM-UKF smoothing
%
%
% Misc.
%   GAUSS_PDF     Multivariate Gaussian PDF
%   GAUSS_RND     Multivariate Gaussian random variables
%   LTI_INT       Integrate LTI ODE with Gaussian Noise
%   LTI_DISC      Discretize LTI ODE with Gaussian Noise
%   RK4           Runge-Kutta integration
%   DER_CHECK     Check derivatives using finite differences
%   SCHOL         Positive semidefinite matrix Cholesky factorization
%   RESAMPSTR     Stratified resampling
%
% /DEMOS/
%
%   /KF_CWPA_DEMO/
%      KF_CWPA_DEMO       CWPA model demonstration with Kalman filter
%
%   /EKF_SINE_DEMO/
%      EKF_SINE_F         Dynamic model function (needed by the augmented UKF)
%      EKF_SINE_H         Measurement model function
%      EKF_SINE_DH_DX     1st order derivative of the measurement model
%      EKF_SINE_D2H_DX2   2nd order derivative of the measurement model
%      EKF_SINE_DEMO      Random Sine Signal demonstration
%
%   /UNGM_DEMO/
%      UNGM_F             Dynamic model function
%      UNGM_DF_DX         1st order derivative of the dynamic model
%      UNGM_D2F_DX2       2nd order derivative of the dynamic model (not used)
%      UNGM_H             Measurement model function
%      UNGM_DH_DX         1st order derivative of the measurement model
%      UNGM_D2H_DX2       2nd order derivative of the measurement model (not used)
%      UNGM_DEMO          UNGM model demonstration
%
%   /BOT_DEMO/
%      BOT_H              Measurement model function
%      BOT_DH_DX          1st order derivative of the measurement model
%      BOT_D2H_DX2        2nd order derivative of the measurement model
%      BOT_DEMO_ALL       BOT demo with EKF and UKF
%      EKFS_BOT_DEMO      BOT demo with EKF
%      UKFS_BOT_DEMO      BOT demo with UKF
%      GHKFS_BOT_DEMO     BOT demo with GHKF
%      CKFS_BOT_DEMO      BOT demo with CKF
%
%   /REENTRY_DEMO/
%      REENTRY_F          Dynamic model function
%      REENTRY_DF         Derivative of the dynamic model
%      REENTRY_H          Measurement model function
%      REENTRY_DH         Derivative of the measurement model
%      REENTRY_IF         Inverse prediction of the dynamic model
%      REENTRY_COND       Generates condition numbers for simulation data
%      MAKE_REENTRY_DATA  Generates the simulation data for reentry dynamics
%      REENTRY_DEMO       Reentry Vehicle Tracking demonstration
%
%   /IMM_DEMO/
%      IMM_DEMO           Tracking a Target with Simple Manouvers demonstration
%
%   /EIMM_DEMO/
%      F_TURN             Dynamic model function for the coordinated turn model
%      F_TURN_DX          Jacobian of the coordinated turn model's dynamic model
%      F_TURN_INV         Inverse dynamics of the coordinated turn model
%      CT_DEMO            Coordinated Turn Model demonstration
%      BOT_H              Measurement model function
%      BOT_DH_DX          1st order derivative of the measurement model
%      BOT_D2H_DX2        2nd order derivative of the measurement model
%      BOTM_DEMO          Bearings Only Tracking of a Manouvering Target Demonstration
%
% Demos currently included in the toolbox, but not documented:
%
% /KF_SINE_DEMO/
%      KF_SINE_DEMO       Sine signal demonstration with Kalman filter

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function vol = ft_headmodel_infinite(varargin)

% FT_HEADMODEL_INFINITE returns an infinitely large homogenous
% volume conduction model. For EEG the volume conductor can be used
% to compute the leadfield of electric current dipoles, for MEG it
% can be used for computing the leadfield of magnmetic dipoles.
%
% Use as
%   vol = ft_headmodel_infinite;
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD

% this is an easy one
vol = [];
vol.type = 'infinite';

author: PatternRecognition

repo: OpenBMI


text:
function neqn = p20_neqn ( )

%*****************************************************************************80
%
%% P20_NEQN returns the number of equations for problem P20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 4;

  return
end


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN THIS SCRIPT TO REMOVE LOCAL PATHS ADDED WHEN RUNNING THE
% CONTROLLER.
%
% In the Simulink model, this script is run every time the user presses
% the 'terminate' button.

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rmpath(genpath('./src/'))

% Create a folder for collecting data
if Config.SAVE_WORKSPACE

    if (~exist(['experiments',date],'dir'))

        mkdir(['experiments',date]);
    end

   matFileList = dir(['./experiments',date,'/*.mat']);
   c           = clock;

   save(['./experiments',date,'/exp_',num2str(c(4)),'-',num2str(c(5)),'.mat'])
end

% Compare the Yarp time with the Simulink time to catch if the required
% integration time step is respected during the simulation
if Config.CHECK_INTEGRATION_TIME && exist('yarp_time','var')

    sim_time   = yarp_time.time;

    % normalize the yarp time over the first value (at t_sim = 0)
    yarp_time0 = yarp_time.signals.values - yarp_time.signals.values(1);

    % fast check of yarp time vs. Simulink time
    if Config.PLOT_INTEGRATION_TIME

        figure
        hold on
        plot(yarp_time0,0:length(yarp_time0)-1,'o')
        plot(sim_time,0:length(sim_time)-1,'o-')
        xlabel('Time [s]')
        ylabel('Iterations')
        grid on
        legend('Yarp Time','Simulink Time')
        title('Yarp time vs. Simulink time')
    end

    % number of times the real time step was bigger than twice the
    % desired time step value
    numOfTimeStepViolations = sum(diff(yarp_time0) > 2*Config.Ts);

    if numOfTimeStepViolations > 1 && numOfTimeStepViolations <= 50

        warning(['The time step tolerance of ', num2str(Config.Ts), '[s] has been violated at least once.'])

    elseif numOfTimeStepViolations > 50

        warning(['The time step tolerance of ', num2str(Config.Ts), '[s] has been violated  more than 50 times.'])
    end
end

author: robotology

repo: whole-body-controllers


text:
function pf =loadPoleFigure_nja(fname,varargin)

assertExtension(fname,'.nja');

comment = '';
fid = efopen(fname);
try
  % read header
  comment = textscan(fid,'%s',45,'delimiter','&','whitespace','');
  fclose(fid);
catch
  interfaceError(fname,fid);
end

try
  % read data
  d = dlmread(fname,'',21,0);
  th = d(:,1)*degree;
  rh = d(:,2)*degree;
  bg  = d(:,4);
  d  = d(:,3);

  h = string2Miller(regexprep([comment{1}{19:21}], '[A-Z=\s]', '', 'ignorecase'));
  r = vector3d.byPolar(th,rh,'antipodal');

  pf = PoleFigure(h,r,d,'BACKGROUND',bg,varargin{:});
catch
  interfaceError(fname);
end

author: mtex-toolbox

repo: mtex

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


text:
function [B, tau] = mex_dgeqrf(A)
% B = mex_dgeqrf(A)
%
% Interface to LAPACK's DGEQRF function.
% Upper part of B is the triangular factor.
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.


text:
function [ output_meta ] = meta2sicd_nisar( HDF5_fid )
%META2SICD_NISAR Converts NISAR SLC HDF5 into a SICD-style metadata structure
%
% Takes as input a file identifier to an open HDF5 file-- as is returned by
% H5F.open
%
% Note that MATLAB's h5disp is generally a good tool for manually browsing
% through HDF5 metadata like that found in NISAR format.
%
% Known issues with simulated NISAR data (from UAVSAR) to date:
%
% 1) ISCEVersion not populated
% 2) Doppler rate provided is positive.  JPL says will it fix in future
% update.
% 3) processedAzimuthBandwidth is acknowledged by JPL to be wrong and will
% be fixed in a future update.
% 4) Weighting description is currently a placeholder, not an accurate
% description of UAVSAR data used for simulated data foramt.
%
% Written by: Wade Schwartzkopf, NGA/Research
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

SECONDS_IN_A_DAY = 24*60*60;

% Query for frequencies and polarization (each of which would create its
% own SICD)
freqs_id = H5D.open(HDF5_fid,'/science/LSAR/identification/listOfFrequencies');
freqs = H5D.read(freqs_id);
freqs = freqs(1,:);
H5D.close(freqs_id);
for i=1:numel(freqs)
    pol_id = H5D.open(HDF5_fid,['/science/LSAR/SLC/swaths/frequency' freqs(i) '/listOfPolarizations']);
    pols{i} = H5D.read(pol_id)';
    H5D.close(pol_id);
end

%% CollectionInfo
% CollectorName could also be deblank(get_hdf_data(HDF5_fid,'/science/LSAR/identification','missionId')')
output_meta.CollectionInfo.CollectorName=get_hdf_attribute(HDF5_fid,'/','mission_name');
% JPL suggested this as best way to form unique string for each collection
output_meta.CollectionInfo.CoreName=[...
    num2str(get_hdf_data(HDF5_fid,'/science/LSAR/identification','absoluteOrbitNumber'),'%07u') ...
    get_hdf_data(HDF5_fid,'/science/LSAR/identification','trackNumber')'];
output_meta.CollectionInfo.CollectType='MONOSTATIC';
output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
output_meta.CollectionInfo.Classification='UNCLASSIFIED';

%% ImageCreation
try
    % TODO: UNTESTED
    % Sample data did not have valid values for this
    output_meta.ImageCreation.Application=['ISCE ' ...
        deblank(get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/processingInformation/algorithms','ISCEVersion')')];
end
output_meta.ImageCreation.Profile='Prototype';

%% ImageData
output_meta.ImageData=struct(); % Just a placeholder
% Most subfields added below in "per band" section
% Used for computing SCP later

%% GeoData
output_meta.GeoData.EarthModel='WGS_84';
% Initially, we just seed this with a very rough value.  Later we will
% put in something more precise.
poly_str = get_hdf_data(HDF5_fid,'/science/LSAR/identification','boundingPolygon')';
poly_str(isletter(poly_str)|poly_str=='('|poly_str==')')='';
poly_str(poly_str==',')=';';
bound_poly=str2num(poly_str);
bound_poly=bound_poly(1:(end-1),:);
rough_center = mean(bound_poly);
output_meta.GeoData.SCP.LLH.Lat=rough_center(2);
output_meta.GeoData.SCP.LLH.Lon=rough_center(1);
output_meta.GeoData.SCP.LLH.HAE=...
    mean(get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/processingInformation/parameters/','referenceTerrainHeight'));
% NISAR stores this as single, which is sufficient for storage, but we need
% double precision for some of the computations later on, to include
% geopositioning of exact SCP.
output_meta.GeoData.SCP.LLH.HAE=double(output_meta.GeoData.SCP.LLH.HAE);
ecf=geodetic_to_ecf([output_meta.GeoData.SCP.LLH.Lat output_meta.GeoData.SCP.LLH.Lon output_meta.GeoData.SCP.LLH.HAE]);
output_meta.GeoData.SCP.ECF.X=ecf(1);
output_meta.GeoData.SCP.ECF.Y=ecf(2);
output_meta.GeoData.SCP.ECF.Z=ecf(3);
% Calling derived_sicd_fields at the end will populate these fields
% with the sensor model, so we don't need to do it here.
% latlon=get_hdf_attribute(dset_id(i),'Top Left Geodetic Coordinates');
% output_meta.GeoData.ImageCorners.ICP.FRFC.Lat=bound_poly(1,1);
% output_meta.GeoData.ImageCorners.ICP.FRFC.Lon=bound_poly(1,2);
% latlon=get_hdf_attribute(dset_id(i),'Bottom Left Geodetic Coordinates');
% output_meta.GeoData.ImageCorners.ICP.FRLC.Lat=bound_poly(2,1);
% output_meta.GeoData.ImageCorners.ICP.FRLC.Lon=bound_poly(2,2);
% latlon=get_hdf_attribute(dset_id(i),'Bottom Right Geodetic Coordinates');
% output_meta.GeoData.ImageCorners.ICP.LRLC.Lat=bound_poly(3,1);
% output_meta.GeoData.ImageCorners.ICP.LRLC.Lon=bound_poly(3,2);
% latlon=get_hdf_attribute(dset_id(i),'Top Right Geodetic Coordinates');
% output_meta.GeoData.ImageCorners.ICP.LRFC.Lat=bound_poly(4,1);
% output_meta.GeoData.ImageCorners.ICP.LRFC.Lon=bound_poly(4,2);

%% Grid
output_meta.Grid.ImagePlane='SLANT';
output_meta.Grid.Type='RGZERO';
output_meta.Grid.Row.Sgn=-1; % Always true for NISAR
output_meta.Grid.Col.Sgn=-1; % Always true for NISAR
output_meta.Grid.Col.KCtr=0;
output_meta.Grid.Row.DeltaKCOAPoly=0;
% TODO: JPL states that uniform weighting in data simulated from UAVSAR is
% placeholder, not an accurate description of the data.  At this point, it
% is not clear what the final weighting description for NISAR will be.
output_meta.Grid.Row.WgtFunct = get_hdf_data(HDF5_fid,...
    '/science/LSAR/SLC/metadata/processingInformation/parameters','rangeChirpWeighting');
if all(output_meta.Grid.Row.WgtFunct(1)==output_meta.Grid.Row.WgtFunct)
    output_meta.Grid.Row.WgtType.WindowName = 'UNIFORM'; % In the sample datasets
else
    output_meta.Grid.Row.WgtType.WindowName = 'UNKNOWN';
end
output_meta.Grid.Col.WgtFunct = get_hdf_data(HDF5_fid,...
    '/science/LSAR/SLC/metadata/processingInformation/parameters','azimuthChirpWeighting');
if all(output_meta.Grid.Col.WgtFunct(1)==output_meta.Grid.Col.WgtFunct)
    output_meta.Grid.Col.WgtType.WindowName = 'UNIFORM'; % In the sample datasets
else
    output_meta.Grid.Col.WgtType.WindowName = 'UNKNOWN';
end
% More subfields added below in "per band" section

%% Timeline
% Using zero Doppler times.  Pulse transmission likely started before this
[collectStart, collectStartFrac]=datenum_w_frac(...
    deblank(get_hdf_data(HDF5_fid, '/science/LSAR/identification', ...
    'zeroDopplerStartTime')'));
[collectEnd, collectEndFrac]=datenum_w_frac(...
    deblank(get_hdf_data(HDF5_fid, '/science/LSAR/identification', ...
    'zeroDopplerEndTime')'));
output_meta.Timeline.CollectStart=collectStart + (collectStartFrac/SECONDS_IN_A_DAY);
output_meta.Timeline.CollectDuration=...
    round((collectEnd-collectStart)*SECONDS_IN_A_DAY) + ... % Convert days to seconds
    (collectEndFrac-collectStartFrac); % Handle fractional seconds
output_meta.Timeline.IPP.Set.TStart=0;
output_meta.Timeline.IPP.Set.TEnd=0; % Apply real value later.  Just a placeholder.
output_meta.Timeline.IPP.Set.IPPStart=uint32(0);
% More subfields added below in "per band" section

%% Position
% Compute polynomial from state vectors
% Get reference time
ref_time = get_ref_time(HDF5_fid,'/science/LSAR/SLC/metadata/orbit/time');
% Times in SICD are with respect to time from start of collect.
ref_time_offset = round((ref_time-collectStart)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
    collectStartFrac; % Handle fractional seconds
state_vector_T  = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/orbit','time')';
state_vector_T = state_vector_T + ref_time_offset; % Make with respect to Timeline.CollectStart
state_vector_pos  = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/orbit','position');
% sv2poly.m shows ways to determine best polynomial order, but we just use 6th here
polyorder=min(6, numel(state_vector_T) - 1);
old_state = warning('off','MATLAB:polyfit:RepeatedPointsOrRescale');
P_x  = polyfit(state_vector_T, state_vector_pos(1,:), polyorder);
P_y  = polyfit(state_vector_T, state_vector_pos(2,:), polyorder);
P_z  = polyfit(state_vector_T, state_vector_pos(3,:), polyorder);
% We don't use these since they are derivable from the position polynomial
% state_vector_vel = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/orbit','velocity');
% state_vector_acc = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/orbit','acceleration');
% P_vx = polyfit(state_vector_T, state_vector_vel(1,:), polyorder);
% P_vy = polyfit(state_vector_T, state_vector_vel(2,:), polyorder);
% P_vz = polyfit(state_vector_T, state_vector_vel(3,:), polyorder);
% P_ax = polyfit(state_vector_T, state_vector_acc(1,:), polyorder);
% P_ay = polyfit(state_vector_T, state_vector_acc(2,:), polyorder);
% P_az = polyfit(state_vector_T, state_vector_acc(3,:), polyorder);
warning(old_state);
% Store position polynomial
output_meta.Position.ARPPoly.X  = P_x(end:-1:1).';
output_meta.Position.ARPPoly.Y  = P_y(end:-1:1).';
output_meta.Position.ARPPoly.Z  = P_z(end:-1:1).';

%% SCPCOA
% This should almost always be left looking for NISAR.
output_meta.SCPCOA.SideOfTrack=get_hdf_data(HDF5_fid,'/science/LSAR/identification','lookDirection');
output_meta.SCPCOA.SideOfTrack=upper(output_meta.SCPCOA.SideOfTrack(1));
% Most subfields added below in "per band" section.

%% ImageFormation
output_meta.ImageFormation.RcvChanProc=struct('NumChanProc',1,'PRFScaleFactor',1);
output_meta.ImageFormation.ImageFormAlgo='RMA';
output_meta.ImageFormation.TStartProc=0;
output_meta.ImageFormation.TEndProc=output_meta.Timeline.CollectDuration;
output_meta.ImageFormation.STBeamComp='NO';
output_meta.ImageFormation.ImageBeamComp='SV';
output_meta.ImageFormation.AzAutofocus='NO';
output_meta.ImageFormation.RgAutofocus='NO';
% More subfields added below in "per band" section
output_meta.RMA.RMAlgoType='OMEGA_K';
output_meta.RMA.ImageType='INCA';
output_meta.RMA.INCA.DopCentroidCOA=true;
% Get reference time
ref_time = get_ref_time(HDF5_fid,'/science/LSAR/SLC/swaths/zeroDopplerTime');
% Times in SICD are with respect to time from start of collect.
ref_time_offset = round((ref_time-collectStart)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
    collectStartFrac; % Handle fractional seconds
zd_t = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/swaths','zeroDopplerTime') + ref_time_offset;
ss_az_s = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/swaths','zeroDopplerTimeSpacing');
if output_meta.SCPCOA.SideOfTrack(1)=='L'
    zd_t = zd_t(end:-1:1);
    ss_az_s = -ss_az_s;
end
grid_r = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/processingInformation/parameters','slantRange');
% Get reference time
ref_time = get_ref_time(HDF5_fid,'/science/LSAR/SLC/metadata/processingInformation/parameters/zeroDopplerTime');
% Times in SICD are with respect to time from start of collect.
ref_time_offset = round((ref_time-collectStart)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
    collectStartFrac; % Handle fractional seconds
grid_zd = get_hdf_data(HDF5_fid,'/science/LSAR/SLC/metadata/processingInformation/parameters','zeroDopplerTime') + ref_time_offset;

%% Radiometric
cal_fields = {'beta0','gamma0','sigma0'};
sicd_cal_fields = {'BetaZeroSFPoly','GammaZeroSFPoly','SigmaZeroSFPoly'};
cal_data = cell(3,1);  % beta0, gamma0, sigma0
for i = 1:numel(cal_fields)
    cal_data{i} = get_hdf_data(HDF5_fid, '/science/LSAR/SLC/metadata/calibrationInformation/geometry',cal_fields{i}).';
    % Get fill value
    gid = H5G.open(HDF5_fid, '/science/LSAR/SLC/metadata/calibrationInformation/geometry/');
    data = H5D.open(gid, cal_fields{i});
    H5G.close(gid);
    attribute = H5A.open(data, '_FillValue');
    H5D.close(data);
    fill = H5A.read(attribute,H5A.get_type(attribute));
    H5A.close(attribute);
    fill = cast(fill,class(cal_data{i}));
    cal_data{i}(cal_data{i}==fill) = NaN;
    if any(diff(cal_data{i}(isfinite(cal_data{i}))))
        % We expect NISAR to be constant beta, although simulated data from
        % UAVSAR is constant sigma
        if strcmp(cal_fields{i},'beta0')
            warning('META2SICD_NISAR:unexpectedCalibrationValues',...
                'Beta-0 values expected to be constant.');
        end
    else
        output_meta.Radiometric.(sicd_cal_fields{i}) = mean(cal_data{i}(:),'omitnan');
    end
    % We will derive non-constant Radiometric fields from the constant one
    % (likely beta).
end
output_meta.Radiometric.NoiseLevel.NoiseLevelType = 'ABSOLUTE';

%% Process fields specific to each polarimetric band
band_independent_meta=output_meta; % Values that are consistent across all bands
grouped_meta=cell(sum(cellfun(@(x) size(x,1),pols)),1);
meta_i=0;
for i=1:numel(freqs)
    freq_meta = band_independent_meta;

    freq_id=H5G.open(HDF5_fid,['/science/LSAR/SLC/swaths/frequency' freqs(i)]);

    %% Grid
    row_ss_id=H5D.open(freq_id,'slantRangeSpacing');
    freq_meta.Grid.Row.SS=H5D.read(row_ss_id);
    H5D.close(row_ss_id);
    % Col.SS is derived after DRateSFPoly below, rather than used from this
    % given field, so that SICD metadata can be internally consistent:
    % col_ss_id=H5D.open(freq_id,'sceneCenterAlongTrackSpacing');
    % output_meta.Grid.Col.SS=H5D.read(col_ss_id);
    % H5D.close(col_ss_id);
    proc_bw_id=H5D.open(freq_id,'processedRangeBandwidth');
    freq_meta.Grid.Row.ImpRespBW=2*H5D.read(proc_bw_id)/SPEED_OF_LIGHT;
    H5D.close(proc_bw_id);
    freq_meta.Grid.Row.DeltaK1=-freq_meta.Grid.Row.ImpRespBW/2;
    freq_meta.Grid.Row.DeltaK2=-freq_meta.Grid.Row.DeltaK1;

    %% Timeline
    prf_id=H5D.open(freq_id,'nominalAcquisitionPRF');
    prf=H5D.read(prf_id);
    H5D.close(prf_id);
    freq_meta.Timeline.IPP.Set.IPPEnd=uint32(floor(prf*band_independent_meta.Timeline.CollectDuration));
    freq_meta.Timeline.IPP.Set.IPPPoly=[0; prf];
    freq_meta.Timeline.IPP.Set.TEnd=freq_meta.Timeline.CollectDuration;

    %% RadarCollection
    fc_id=H5D.open(freq_id,'acquiredCenterFrequency');
    fc=H5D.read(fc_id);
    H5D.close(fc_id);
    bw_id=H5D.open(freq_id,'acquiredRangeBandwidth');
    bw=H5D.read(bw_id);
    H5D.close(bw_id);
    freq_meta.RadarCollection.TxFrequency.Min=fc-(bw/2);
    freq_meta.RadarCollection.TxFrequency.Max=fc+(bw/2);
    % No waveform parameters provided
    for j=1:size(pols{i},1)
        freq_meta.RadarCollection.RcvChannels.ChanParameters(j).TxRcvPolarization=[pols{i}(j,1) ':' pols{i}(j,2)];
    end
    tx_pol = unique(pols{i}(:,1));
    if numel(tx_pol)==1
        freq_meta.RadarCollection.TxPolarization = tx_pol;
    else
        freq_meta.RadarCollection.TxPolarization = 'SEQUENCE';
        for j = 1:numel(tx_pol)
            freq_meta.RadarCollection.TxSequence.TxStep(j).WFIndex = j;
            freq_meta.RadarCollection.TxSequence.TxStep(j).TxPolarization = tx_pol(j);
        end
    end

    %% ImageFormation
    fc_id=H5D.open(freq_id,'processedCenterFrequency');
    fc=H5D.read(fc_id);
    H5D.close(fc_id);
    bw_id=H5D.open(freq_id,'processedRangeBandwidth');
    bw=H5D.read(bw_id);
    H5D.close(bw_id);
    freq_meta.ImageFormation.TxFrequencyProc.MinProc=fc-(bw/2);
    freq_meta.ImageFormation.TxFrequencyProc.MaxProc=fc+(bw/2);
    range_id=H5D.open(freq_id,'slantRange');
    r_ca_sampled=H5D.read(range_id);
    H5D.close(range_id);
    % TODO: processedAzimuthBandwidth acknowledged by JPL to be wrong in
    % simulated datasets.
    dop_bw_id=H5D.open(freq_id,'processedAzimuthBandwidth');
    dop_bw=H5D.read(dop_bw_id);
    H5D.close(dop_bw_id);
    H5G.close(freq_id);
    dopcentroid_sampled = get_hdf_data(HDF5_fid,[...
        '/science/LSAR/SLC/metadata/processingInformation/parameters/frequency' freqs(i)], ...
        'dopplerCentroid').';
    doprate_sampled = get_hdf_data(HDF5_fid,[...
        '/science/LSAR/SLC/metadata/processingInformation/parameters/frequency' freqs(i)], ...
        'azimuthFMRate').';

    for j=1:size(pols{i},1)
        pol_meta=freq_meta;

        %% ImageData
        dset_id=H5D.open(HDF5_fid,['/science/LSAR/SLC/swaths/frequency' freqs(i) '/' pols{i}(j,:)]);
        dspace_id=H5D.get_space(dset_id);
        [~, datasize] = H5S.get_simple_extent_dims(dspace_id);
        pol_meta.ImageData.NumCols=uint32(datasize(1));
        pol_meta.ImageData.NumRows=uint32(datasize(2));
        pol_meta.ImageData.FullImage=pol_meta.ImageData;
        if (H5T.get_class(H5D.get_type(dset_id)) == H5ML.get_constant_value('H5T_COMPOUND'))
            if (H5T.get_member_class(H5D.get_type(dset_id),0) == H5ML.get_constant_value('H5T_FLOAT'))
                % Could be either 16- or 32-bit float.  Either way it will
                % cast to this SICD type.
                pol_meta.ImageData.PixelType='RE32F_IM32F';
            elseif (H5T.get_member_class(H5D.get_type(dset_id),0) == H5ML.get_constant_value('H5T_INTEGER')) && ...
                    H5T.get_size(H5D.get_type(dset_id)) == 4
                pol_meta.ImageData.PixelType='RE16I_IM16I';
            else
                error('META2SICD_NISAR:unrecognizedDataType','Must be 16- or 32-bit complex float.');
            end
        else
            error('META2SICD_NISAR:unrecognizedDataType','Must be 16- or 32-bit complex float.');
        end
        H5S.close(dspace_id);
        H5D.close(dset_id);
        pol_meta.ImageData.FirstRow=uint32(0); pol_meta.ImageData.FirstCol=uint32(0);
        % SCP pick is arbitrary.  Just pick something near center.
        pol_meta.ImageData.SCPPixel.Col = floor(datasize(1)/2);
        pol_meta.ImageData.SCPPixel.Row = floor(datasize(2)/2);

        %% ImageFormation
        pol_meta.ImageFormation.RcvChanProc.ChanIndex=meta_i;
        pol_meta.ImageFormation.TxRcvPolarizationProc=...
            freq_meta.RadarCollection.RcvChannels.ChanParameters(j).TxRcvPolarization;

        %% RMA
        pol_meta.RMA.INCA.R_CA_SCP = r_ca_sampled(pol_meta.ImageData.SCPPixel.Row+1);
        scp_ca_time = zd_t(pol_meta.ImageData.SCPPixel.Col+1);
        % Compute DRateSFPoly
        % For the purposes of the DRateSFPoly computation, we ignore any
        % changes in velocity or doppler rate over the azimuth dimension.
        pos_coefs = [P_x(:) P_y(:) P_z(:)];
        % Velocity is derivate of position.
        vel_coefs=pos_coefs(1:end-1,:).*repmat(((size(pos_coefs,1)-1):-1:1)',[1 3]);
        vel_x = polyval(vel_coefs(:,1), scp_ca_time);
        vel_y = polyval(vel_coefs(:,2), scp_ca_time);
        vel_z = polyval(vel_coefs(:,3), scp_ca_time);
        vm_ca_sq = vel_x.^2 + vel_y.^2 + vel_z.^2; % Magnitude of the velocity squared
        r_ca_poly = [pol_meta.RMA.INCA.R_CA_SCP; 1]; % Polynomial representing range as a function of range distance from SCP
        [~, min_ind] = min(abs(grid_zd - scp_ca_time));  % Closest Doppler rate polynomial to SCP
        coords_rg_m = grid_r - pol_meta.RMA.INCA.R_CA_SCP;
        % TODO: Is it the NISAR convention to give the absolute value of
        % Doppler rate? (It should be negative.) Sample datasets appear
        % this way.
        old_warning_state=warning('off','MATLAB:polyfit:RepeatedPointsOrRescale');
        dop_rate_poly=polyfit(coords_rg_m, -doprate_sampled(min_ind,:).', 4);
        warning(old_warning_state);
        % figure; plot(r,dop_rate(:,min_ind),r,polyval(dop_rate_poly,r));  % Show closeness of fit
        dop_rate_poly = dop_rate_poly(end:-1:1);
        pol_meta.RMA.INCA.DRateSFPoly = - conv(dop_rate_poly,r_ca_poly) * ... % Multiplication of two polynomials is just a convolution of their coefficients
            SPEED_OF_LIGHT / (2 * fc * vm_ca_sq(1)); % Assumes a SGN of -1
        % Fields dependent on Doppler rate
        pol_meta.Grid.Col.SS = sqrt(vm_ca_sq(1)) * abs(ss_az_s) * pol_meta.RMA.INCA.DRateSFPoly(1,1);
        % Should be close to this:
        % colss_id=H5D.open(freq_id,'sceneCenterAlongTrackSpacing');
        % pol_meta.Grid.Col.SS = H5D.read(colss_id)
        % H5D.close(colss_id);
        pol_meta.Grid.Col.ImpRespBW = ... % Convert to azimuth spatial bandwidth (cycles per meter)
            min(dop_bw*abs(ss_az_s),1)/pol_meta.Grid.Col.SS; % Can't have more bandwidth in data than sample spacing
        pol_meta.RMA.INCA.TimeCAPoly = [scp_ca_time; ... % With respect to start of collect
            ss_az_s/pol_meta.Grid.Col.SS]; % Convert zero doppler spacing from sec/pixels to sec/meters
        % TimeCOAPoly/DopCentroidPoly/DeltaKCOAPoly
        POLY_ORDER = 3; % Order of polynomial which we want to compute
        coords_az_m = (grid_zd - scp_ca_time) * pol_meta.Grid.Col.SS / ss_az_s;
        timeca_sampled = repmat(grid_zd,[1 numel(grid_r)]);
        % TimeCOAPoly=TimeCA+(DopCentroid/dop_rate)
        timecoa_sampled = timeca_sampled+(dopcentroid_sampled./doprate_sampled);
        pol_meta.RMA.INCA.DopCentroidPoly=polyfit2d(dopcentroid_sampled, coords_az_m, coords_rg_m, POLY_ORDER);
        pol_meta.Grid.Col.DeltaKCOAPoly=...
            pol_meta.RMA.INCA.DopCentroidPoly*ss_az_s/pol_meta.Grid.Col.SS;
        pol_meta.Grid.TimeCOAPoly=polyfit2d(timecoa_sampled, coords_az_m, coords_rg_m, POLY_ORDER);

        %% Radiometric
        % In the simulated data (from UAVSAR), we expect these noise values
        % to be only valid for the same incidence angle as NISAR (small
        % part of swath).  We also expect the values to be in dB, not
        % linear.
        nesz = get_hdf_data(HDF5_fid,['/science/LSAR/SLC/metadata/calibrationInformation/frequency' freqs(i) '/' pols{i}(j,:)], 'nes0').';
        sigma0sf = cal_data{3};
        noise_samples = nesz - (10*log10(sigma0sf));
        pol_meta.Radiometric.NoiseLevel.NoisePoly = ...
            polyfit2d(noise_samples, coords_az_m, coords_rg_m, POLY_ORDER);
        % The conversion between radiometric fields in simulated data don't
        % seem to be valid, so we will skip using the non-constant terms
        % for now.  But this is how one would do it:
        % Fit polynomial only to valid data in image
        % az_bounds1 = [-double(pol_meta.ImageData.SCPPixel.Col) ...
        %     double(pol_meta.ImageData.NumCols - pol_meta.ImageData.SCPPixel.Col - 1)] *...
        %     pol_meta.Grid.Col.SS;
        % rg_bounds1 = [-double(pol_meta.ImageData.SCPPixel.Row) ...
        %     double(pol_meta.ImageData.NumRows - pol_meta.ImageData.SCPPixel.Row - 1)] *...
        %     pol_meta.Grid.Row.SS;
        % for k = 1:numel(cal_fields)
        %     % Only fit valid values
        %     rg_val = all(isfinite(cal_data{k}),1)&coords_rg_m'>rg_bounds1(1)&coords_rg_m'<rg_bounds1(2);
        %     if ~any(diff(cal_data{k}(isfinite(cal_data{k}))))  % Constant across image
        %         pol_meta.Radiometric.(sicd_cal_fields{k}) = mean(cal_data{k}(:),'omitnan');
        %     elseif all(all(diff(cal_data{k},1,1)==0) | all(~isfinite(cal_data{k})))% Only varies in range
        %         rad_poly = polyfit(coords_rg_m(rg_val), cal_data{k}(1,rg_val), 7);
        %         pol_meta.Radiometric.(sicd_cal_fields{k}) = rad_poly(end:-1:1).';
        %     else  % Varies 2-dimensionally.  Not tested.  POLY_ORDER may not be reasonable.
        %         pol_meta.Radiometric.(sicd_cal_fields{k}) = ...
        %             polyfit2d(cal_data{k}(:,rg_val), coords_az_m, coords_rg_m(rg_val), POLY_ORDER);
        %     end
        % end

        %% GeoData
        % Now that sensor model fields have been populated, we can populate
        % GeoData.SCP more precisely.
        ecf = point_image_to_ground([pol_meta.ImageData.SCPPixel.Row;pol_meta.ImageData.SCPPixel.Col],pol_meta);
        pol_meta.GeoData.SCP.ECF.X=ecf(1);
        pol_meta.GeoData.SCP.ECF.Y=ecf(2);
        pol_meta.GeoData.SCP.ECF.Z=ecf(3);
        llh=ecf_to_geodetic([pol_meta.GeoData.SCP.ECF.X pol_meta.GeoData.SCP.ECF.Y pol_meta.GeoData.SCP.ECF.Z]);
        pol_meta.GeoData.SCP.LLH.Lat=llh(1);
        pol_meta.GeoData.SCP.LLH.Lon=llh(2);
        pol_meta.GeoData.SCP.LLH.HAE=llh(3);

        pol_meta = derived_sicd_fields(pol_meta);
        meta_i=meta_i+1;
        grouped_meta{meta_i}=pol_meta;
    end
end
output_meta = grouped_meta;

end

% Reads HDF5 attribute
function value = get_hdf_attribute(hid_t,path,attribute_name)
gid = H5G.open(hid_t,path);
attribute = H5A.open(gid,attribute_name);
value = H5A.read(attribute,H5A.get_type(attribute));
% The behavior of whether HDF5 attributes are returned in cells seems to
% have changed in 2020a. The behavior also changed in h5readatt.
if iscell(value)
    value = value{1};
end
H5G.close(gid);
H5A.close(attribute);
end

% Reads HDF5 data
function value = get_hdf_data(hid_t,path,data_name)
gid = H5G.open(hid_t,path);
data = H5D.open(gid,data_name);
value = H5D.read(data);
H5G.close(gid);
H5D.close(data);
end

function ref_time = get_ref_time(hid_t,data_path)
    data = H5D.open(hid_t,data_path);
    attribute = H5A.open(data,'units');
    ref_time_string = H5A.read(attribute,H5A.get_type(attribute))';
    H5A.close(attribute);
    H5D.close(data);
    [ref_time_string,~] = regexp(ref_time_string,'seconds since (\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2})\s*$','tokens','match');
    if numel(ref_time_string)==1  % Should always be the case according to JPL
        ref_time = datenum(ref_time_string{1}{1},'yyyy-mm-dd HH:MM:SS');
    else
        error('META2SICD_NISAR:unrecognizedTimeUnit','Time units must be seconds since a reference UTC time.');
    end
end

% MATLAB's datenum function won't handle precise times down under a
% millisecond, because 1) It won't accept a format with more than 3 .FFF in
% the string description of the date format and 2) the resulting serial
% date number is stored in days from 00-JAN-0000 and just doesn't have
% enough bits to handle fractional seconds to the level we want.  Here we
% handle the fractional seconds separately so we can read date with the
% precision we need.
function [datenum_s, datenum_frac] = datenum_w_frac(datestring)
    datenum_s = datenum(datestring,'yyyy-mm-ddTHH:MM:SS');
    datenum_frac = str2double(regexp(datestring,'\.\d*','match'));
    if any(isnan(datenum_frac))||isempty(datenum_frac), datenum_frac = 0; end
end

function p = polyfit2d(samples, dim1, dim2, poly_order)
    % Least squares fit for 2D polynomial
    % A*x = b
    [dim1_2d, dim2_2d] = ndgrid(dim1, dim2);
    a = zeros(numel(samples), (poly_order+1)^2);
    for k = 0:poly_order
        for m = 0:poly_order
            a(:,k*(poly_order+1)+m+1) = (dim2_2d(:).^m).*(dim1_2d(:).^k);
        end
    end
    b = zeros((poly_order+1)^2,1);
    for k=1:((poly_order+1)^2)
        b(k)=sum(samples(:).*a(:,k)); % center of aperture
    end
    A=zeros((poly_order+1)^2);
    for k=1:((poly_order+1)^2)
        for m=1:((poly_order+1)^2)
            A(k,m)=sum(a(:,k).*a(:,m));
        end
    end
    old_warning_state=warning('off','MATLAB:nearlySingularMatrix');
    x_coa=A\b; % MATLAB often flags this as badly scaled, but results still appear valid
    warning(old_warning_state);
    p=reshape(x_coa, poly_order+1, poly_order+1);
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////


text:
function [statusBoxMsg] = setINIStatusBoxMsg()
%setINIStatusBoxMsg Summary of this function goes here
%   Detailed explanation goes here

    statusBoxMsg = {['KSP Trajectory Optimization Tool v', getKSPTOTVersionNumStr(), sprintf(' (R%s)', version('-release'))], 'Written By Arrowstar (C) 2022', statusBoxHR()};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function pde = elli3DcircIntf2(am,ap,bm,bp,r,x0,y0,z0,a11,a12,a)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,...
    'B',@B,'Bm',@Bm,'Bp',@Bp);

pde.am = am;
pde.ap = ap;
pde.bm = bm;
pde.bp = bp;
%% interface function
    function u = intf(x,y,z)
        u = ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(1/2)/r-1;
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    function u = um1(x,y,z)
        u1 = uker1(x,y,z).*(x-x0);
        u2 = uker2(x,y,z);
        u = (u1 + u2)/bm;
    end
    function u = um2(x,y,z)
        u1 = uker1(x,y,z).*(y-y0);
        u2 = uker2(x,y,z);
        u = (u1 + u2)/bm;
    end
    function u = um3(x,y,z)
        u1 = uker1(x,y,z).*(z-z0);
        u2 = uker2(x,y,z);
        u = (u1 + u2)/bm;
    end
    function u = up1(x,y,z)
        u1 = uker1(x,y,z).*(x-x0);
        u2 = uker2(x,y,z);
        u = (u1 + u2)/bp;
    end
    function u = up2(x,y,z)
        u1 = uker1(x,y,z).*(y-y0);
        u2 = uker2(x,y,z);
        u = (u1 + u2)/bp;
    end
    function u = up3(x,y,z)
        u1 = uker1(x,y,z).*(z-z0);
        u2 = uker2(x,y,z);
        u = (u1 + u2)/bp;
    end
    function u = uker1(x,y,z)
        u = nthroot((x-x0).^2 + (y-y0).^2 + (z-z0).^2 - r^2,a11);
        u = u.^a12;
    end
    function u = uker2(x,y,z)
        u = ((x-x0).^2 + (y-y0).^2 + (z-z0).^2 - r^2).^a;
    end
%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = (Duker(x,y,z).*(y-y0)*2 - Duker(x,y,z).*(z-z0)*2)/bm;
    end
    function u = Dyum(x,y,z)
        u = (Duker(x,y,z).*(z-z0)*2 - Duker(x,y,z).*(x-x0)*2)/bm;
    end
    function u = Dzum(x,y,z)
        u = (Duker(x,y,z).*(x-x0)*2 - Duker(x,y,z).*(y-y0)*2)/bm;
    end
    function u = Dxup(x,y,z)
        u = (Duker(x,y,z).*(y-y0)*2 - Duker(x,y,z).*(z-z0)*2)/bp;
    end
    function u = Dyup(x,y,z)
        u = (Duker(x,y,z).*(z-z0)*2 - Duker(x,y,z).*(x-x0)*2)/bp;
    end
    function u = Dzup(x,y,z)
        u = (Duker(x,y,z).*(x-x0)*2 - Duker(x,y,z).*(y-y0)*2)/bp;
    end

    function u = Duker(x,y,z)
        u = a*((x-x0).^2 + (y-y0).^2 + (z-z0).^2 - r^2).^(a-1);
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        xh = x-x0; yh = y-y0; zh = z-z0;
        u = (-4*Duker(x,y,z) + DDuker(x,y,z).*(xh.*(yh+zh)-yh.^2-zh.^2)) + bm*um1(x,y,z);
    end
    function u = fm2(x,y,z)
        xh = x-x0; yh = y-y0; zh = z-z0;
        u = (-4*Duker(x,y,z) + DDuker(x,y,z).*(yh.*(xh+zh)-xh.^2-zh.^2)) + bm*um2(x,y,z);
    end
    function u = fm3(x,y,z)
        xh = x-x0; yh = y-y0; zh = z-z0;
        u = (-4*Duker(x,y,z) + DDuker(x,y,z).*(zh.*(xh+yh)-xh.^2-yh.^2)) + bm*um3(x,y,z);
    end
    function u = fp1(x,y,z)
        xh = x-x0; yh = y-y0; zh = z-z0;
        u = (-4*Duker(x,y,z) + DDuker(x,y,z).*(xh.*(yh+zh)-yh.^2-zh.^2)) + bp*up1(x,y,z);
    end
    function u = fp2(x,y,z)
        xh = x-x0; yh = y-y0; zh = z-z0;
        u = (-4*Duker(x,y,z) + DDuker(x,y,z).*(yh.*(xh+zh)-xh.^2-zh.^2)) + bp*up2(x,y,z);
    end
    function u = fp3(x,y,z)
        xh = x-x0; yh = y-y0; zh = z-z0;
        u = (-4*Duker(x,y,z) + DDuker(x,y,z).*(zh.*(xh+yh)-xh.^2-yh.^2)) + bp*up3(x,y,z);
    end

    function u = DDuker(x,y,z)
        u = a*(a-1)*((x-x0).^2 + (y-y0).^2 + (z-z0).^2 - r^2).^(a-2);
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = am*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = ap*ones(size(x));
    end

%% Mass coefficient function
    function u = B(x,y,z)
        u = Bm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Bp(x(id),y(id),z(id));
    end
    function u = Bm(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Bp(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function success=testJAABAWindowdata()

% This tests to make sure that loading cached windowdata and computing
% windowdata again learn the same classifier.

% AL20150602. This test doesn't work when I run it on the current codebase;
% I think it's because window feature computation is no longer
% deterministic, even with the setting of RNG state below.
% ComputeWindowFeatures() is called within a parfor and has a call to
% randsample.

%%
success = false;
jabFileName='/groups/branson/home/kabram/bransonlab/projects/JAABA/test_data/test_windowdata.jab';
gtMode = false;
data=JLabelData('setstatusfn',@(str)(fprintf('%s\n',str)), ...
                'clearstatusfn',@()(nop()));
data.openJabFile(jabFileName,gtMode);

if matlabpool('size')>0, matlabpool('close'); end
oldcl = data.classifier;

s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);

data.Train;
newcl1 = data.classifier;

data.setWindowFeaturesParams(data.windowfeaturesparams);

s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);

data.Train;
newcl2 = data.classifier;

if ~isequal(oldcl,newcl1) || ~isequal(newcl1,newcl2),
  error('Classifiers trained are not same');

end
success = true;

author: kristinbranson

repo: JAABA