function [m, n, S, names] = readmfasta(filename);
%READFASTA read set of sequence from multiple fasta format.

% Use:
% [m, n, S, names] = readfasta(filename), 
% 'filename' is a name of a multiple sequence fasta file;
% 'S' is an m x n matrix ('m' sequences of lenght 'n')
% 'names' is the array of strings of sequence names   
% Latest modified by Steinar Thorvaldsen 4. Mars 2005

MAXLENGTH = 10000;
MAXNAME = 200;

file = fopen(filename, 'r');
disp(['Reading ',filename]);

% ----------------------------------------------------------
% First determine the maximum length of the sequence.
m=0;    % the number of sequences
n=0;    % the maximum length
cn = 0; % current sequence length

while 1
	[x,nr] = fscanf(file,'%c',1);
   if nr == 0 break; end;
   if x =='>'  % new sequence started
		%if cn > n n=cn; end; %moved!
		cn = 0;
		fgets(file);
		m=m+1;
   else
        if isstrprop(x,'alphanum') | x=='-' | x=='~'
		%if isletter(x) | x=='-' | x=='~'
			cn=cn+1;
		end;
        if cn > n n=cn; end;
	end;
end

% ----------------------------------------------------------
% Go throught the file
Ss = char(n); S = [];
str = zeros(1,MAXNAME);
frewind(file);
names=[];
i=0;j=1;
while 1
	[x,nr] = fscanf(file,'%c',1);
   if nr == 0 break; end;
	if x =='>'  % new sequence started
		if i~= 0 % store the sequence
			S=strvcat(S,Ss);
			Ss = []; Ss = char(n);
		end;
		str=fgetl(file); % read the name, we remove the '>' symbol
		names=strvcat(names,str);
		disp(str);
		i=i+1;
		j=1;
    else
        if isstrprop(x,'alphanum') | x=='-' | x=='~'
		%if isletter(x) | x=='-' | x=='~' OLD ver
		% processing the sequence symbol
	   	Ss(j) = upper(x); 
			j=j+1;
		end;
	end;
end

S=strvcat(S,Ss);
fclose(file);

