cosmos2m( fNameDir, fNameRoot , meshonly , meshindeces , whichModes)

function [ modedata , theMesh ]= cosmos2m( fNameDir, fNameRoot , meshonly , meshindeces , whichModes)
% COSMOS2M  Create eigenvectors/eigenvalue matrices from CosmosWorks FEA
%    
% Processes all the files beginning with fNameRoot, where each file contains the
% entire node data for a given mode. Returns vectors of eigenvalues and eigenmodes
% 
% Assumes modeshapes are in files names [ fnameRoot '_mode' foo ]
% and frequencies in [ fnameRoot '_freq' ]
%
% The argument meshindeces shrinks the mesh to contain only the specified rows
% of themesh obtained from the concatenation of mesh files. Note that
% meshindeces does NOT contain the node ids as returned by CosmosWorks in the
% first column; it only contains the indeces into the mesh as would have been
% returned by this function with a meshonly==true call.
%
% NOTE: if only N mode shape files are found, then only the first N frequencies
% are returned.

% Copyright 2006, The MathWorks, Inc.
   
    p = pwd;

    if nargin < 3
        meshonly = [];
    end
    
    if nargin < 4
        meshindeces = [];
    end
    
    if nargin < 5
        whichModes = [];
    end
    
    if isempty( meshonly )
        meshonly = false;
    end
    
    try
        cd( fNameDir );
        

        
        %
        % Mesh
        % 
        fileNames = dir( [ fNameRoot '_mesh*' ] );
        numberOfModesRead = length( fileNames );
        
        theMesh = [];
        for thisFileIndex = 1:numberOfModesRead;
            
            thisFileName = fileNames( thisFileIndex ).name;
            
            fprintf( 1 , 'Processing %s\n' , thisFileName );
                 
            theMesh = [theMesh ; load( thisFileName ) ];
        end
        theMesh = unique( theMesh , 'rows' );


        if meshonly
            modedata = [];           
        else
        
        meshNodeIds = theMesh( meshindeces , 1 );
        if ~isempty( meshindeces )
            theMesh = theMesh( meshindeces , : );
        end
                    
        %        meshindeces = 1:length(theMesh);
        
        fileNames = dir( [ fNameRoot '_mode*' ] );
        numberOfModesRead = length( fileNames );
        
        for thisFileIndex = 1:numberOfModesRead;
            
            thisFileName = fileNames( thisFileIndex ).name;
            
            fprintf( 1 , 'Processing %s\n' , thisFileName );
            modedata.mode( thisFileIndex ).filename = thisFileName;
            
            thisFile = fopen( thisFileName );
            if thisFile < 0 
                error( 'Could not open file %s' , thisFileName )
            end
            
            modedata.mode( thisFileIndex ).scanned.header = textscan( thisFile , '%*[^\n]' , 5); 
            
            modedata.mode( thisFileIndex ).scanned.info = textscan( thisFile , '%d,%n,%n,%n,%n' );
            nodal_displacement_X =  modedata.mode( thisFileIndex ).scanned.info( 2 );
            nodal_displacement_Y =  modedata.mode( thisFileIndex ).scanned.info( 3 );
            nodal_displacement_Z =  modedata.mode( thisFileIndex ).scanned.info( 4 );
            nodal_displacement_R =  modedata.mode( thisFileIndex ).scanned.info( 5 );
            
            
            modedata.mode( thisFileIndex ).nodal_displacement_X =  nodal_displacement_X{:}(meshNodeIds);
            modedata.mode( thisFileIndex ).nodal_displacement_Y =  nodal_displacement_Y{:}(meshNodeIds);
            modedata.mode( thisFileIndex ).nodal_displacement_Z =  nodal_displacement_Z{:}(meshNodeIds);
            modedata.mode( thisFileIndex ).nodal_displacement_R =  nodal_displacement_R{:}(meshNodeIds);
            
            fclose( thisFile );
            
        end
        
        
        %
        % Frequencies
        %
        
        fileNames = ls( [ fNameRoot '_freq*' ] );
        thisFileName = fileNames( 1 : (end) );

        
        thisFile = fopen( thisFileName );
        if thisFile < 0 
            error( 'Could not open file %s' , thisFileName )
        end   
        modedata.freq.scanned.header = textscan( thisFile , '%*[^\n]' , 3); 
        modedata.freq.scanned.info   = textscan( thisFile , '%d,%n,%n,%n%*[^\n]' );    
        fclose( thisFile );
        
        frequency_angular      = modedata.freq.scanned.info{ 2 };
        frequency_linear       = modedata.freq.scanned.info{ 3 };
        frequency_period       = modedata.freq.scanned.info{ 4 };
        
        modedata.frequency_angular =  frequency_angular( 1 : numberOfModesRead ) ;
        modedata.frequency_linear  =  frequency_linear ( 1 : numberOfModesRead ) ;
        modedata.frequency_period  =  frequency_period ( 1 : numberOfModesRead ) ;
        
        
        
        
        % 
        % Compact and condition the data
        %    
        modedata.nodal_displacements_X = [ modedata.mode.nodal_displacement_X ];
        modedata.nodal_displacements_Y = [ modedata.mode.nodal_displacement_Y ];
        modedata.nodal_displacements_Z = [ modedata.mode.nodal_displacement_Z ];
        modedata.nodal_displacements_R = [ modedata.mode.nodal_displacement_R ];
       
        %
        % Select modes
        %
        if ~isempty( whichModes )
            modedata.frequency_angular =  modedata.frequency_angular( whichModes );
            modedata.frequency_linear  =  modedata.frequency_linear(  whichModes );
            modedata.frequency_period  =  modedata.frequency_period(  whichModes );
            modedata.nodal_displacements_X = modedata.nodal_displacements_X( : , whichModes );
            modedata.nodal_displacements_Y = modedata.nodal_displacements_Y( : , whichModes );
            modedata.nodal_displacements_Z = modedata.nodal_displacements_Z( : , whichModes );
            modedata.nodal_displacements_R = modedata.nodal_displacements_R( : , whichModes );
            modedata.num_modes = length( whichModes );

        end
        
        [ modedata.num_nodes , modedata.num_modes ] = size( modedata.nodal_displacements_X );
        
        end % ~meshOnly
        
        cd( p );
    
    catch
        cd( p );
        rethrow( lasterror );
    end