Back to home page

MITgcm

 
 

    


Warning, /utils/matlab/rdmds.m is written in an unsupported language. File is not indexed.

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
37662e6918 Jean*0001 function [AA,itrs,MM] = rdmds(fnamearg,varargin)
3323410f65 Alis*0002 % RDMDS  Read MITgcmUV meta/data files
a15d675b90 Alis*0003 %
00398243b4 Alis*0004 % A = RDMDS(FNAME)
                0005 % A = RDMDS(FNAME,ITER)
                0006 % A = RDMDS(FNAME,[ITER1 ITER2 ...])
52e6460920 Alis*0007 % A = RDMDS(FNAME,NaN)
                0008 % A = RDMDS(FNAME,Inf)
e6cfd13877 Jean*0009 % [A,ITS,M] = RDMDS(FNAME,[...])
96629ec1c2 Alis*0010 % A = RDMDS(FNAME,[...],'rec',RECNUM)
00398243b4 Alis*0011 %
3323410f65 Alis*0012 %   A = RDMDS(FNAME) reads data described by meta/data file format.
                0013 %   FNAME is a string containing the "head" of the file names.
e6cfd13877 Jean*0014 %
3323410f65 Alis*0015 %   eg. To load the meta-data files
                0016 %       T.0000002880.000.000.meta, T.0000002880.000.000.data
                0017 %       T.0000002880.001.000.meta, T.0000002880.001.000.data
                0018 %       T.0000002880.002.000.meta, T.0000002880.002.000.data
                0019 %       T.0000002880.003.000.meta, T.0000002880.003.000.data
                0020 %   use
                0021 %      >> A=rdmds('T.0000002880');
                0022 %      >> size(A)
                0023 %   ans =
                0024 %      64    32     5
                0025 %   eg. To load a multiple record file
                0026 %      >> A=rdmds('pickup.0000002880');
                0027 %      >> size(A)
                0028 %   ans =
                0029 %      64    32     5    61
e6cfd13877 Jean*0030 %
                0031 %
3323410f65 Alis*0032 %   A = RDMDS(FNAME,ITER) reads data described by meta/data file format.
                0033 %   FNAME is a string containing the "head" of the file name excluding the
                0034 %   10-digit iterartion number.
                0035 %   ITER is a vector of positive integers that will expand to the 10-digit
                0036 %   number in the file name.
52e6460920 Alis*0037 %   If ITER=NaN, all iterations will be read.
                0038 %   If ITER=Inf, the last (highest) iteration will be read.
96629ec1c2 Alis*0039 %
3323410f65 Alis*0040 %   eg. To repeat above operation
                0041 %      >> A=rdmds('T',2880);
                0042 %   eg. To read multiple time steps
                0043 %      >> A=rdmds('T',[0 1440 2880]);
52e6460920 Alis*0044 %   eg. To read all time steps
                0045 %      >> [A,ITS]=rdmds('T',NaN);
                0046 %   eg. To read the last time step
                0047 %      >> [A,IT]=rdmds('T',Inf);
3323410f65 Alis*0048 %   Note: this form can not read files with no iteration count in file name.
96629ec1c2 Alis*0049 %
                0050 %
                0051 %   A = RDMDS(FNAME,[...],'rec',RECNUM) reads individual records from
                0052 %   multiple record files.
e6cfd13877 Jean*0053 %
96629ec1c2 Alis*0054 %   eg. To read a single record from a multi-record file
                0055 %      >> [A,IT]=rdmds('pickup.ckptA',11);
                0056 %   eg. To read several records from a multi-record file
                0057 %      >> [A,IT]=rdmds('pickup',Inf,'rec',[1:5 8 12]);
                0058 %
e6cfd13877 Jean*0059 %
3323410f65 Alis*0060 %   A = RDMDS(FNAME,ITER,MACHINEFORMAT) allows the machine format to be
96629ec1c2 Alis*0061 %   A = RDMDS(FNAME,MACHINEFORMAT)
3323410f65 Alis*0062 %   specified which MACHINEFORMAT is on of the following strings:
                0063 %     'n' 'l' 'b' 'd' 'g' 'c' 'a' 's'  - see FOPEN for more details
96629ec1c2 Alis*0064 %
e6cfd13877 Jean*0065 
52e6460920 Alis*0066 AA=[];
37662e6918 Jean*0067 itrs=[];
e6cfd13877 Jean*0068 MM=[];
a15d675b90 Alis*0069 
                0070 % Default options
                0071 ieee='b';
00398243b4 Alis*0072 fname=fnamearg;
96629ec1c2 Alis*0073 userecords=0;
                0074 recnum=[];
a15d675b90 Alis*0075 
                0076 % Check optional arguments
00398243b4 Alis*0077 for ind=1:size(varargin,2);
                0078  arg=varargin{ind};
                0079  if ischar(arg)
                0080   if strcmp(arg,'n') | strcmp(arg,'native')
                0081    ieee='n';
                0082   elseif strcmp(arg,'l') | strcmp(arg,'ieee-le')
                0083    ieee='l';
                0084   elseif strcmp(arg,'b') | strcmp(arg,'ieee-be')
                0085    ieee='b';
                0086   elseif strcmp(arg,'c') | strcmp(arg,'cray')
                0087    ieee='c';
                0088   elseif strcmp(arg,'a') | strcmp(arg,'ieee-le.l64')
                0089    ieee='a';
                0090   elseif strcmp(arg,'s') | strcmp(arg,'ieee-be.l64')
                0091    ieee='s';
96629ec1c2 Alis*0092   elseif strcmp(arg,'rec')
                0093    userecords=1;
00398243b4 Alis*0094   else
                0095    error(['Optional argument ' arg ' is unknown'])
                0096   end
a15d675b90 Alis*0097  else
96629ec1c2 Alis*0098   if userecords==1
                0099    recnum=arg;
37662e6918 Jean*0100   elseif isempty(itrs)
356a4fc521 Alis*0101   if isnan(arg)
37662e6918 Jean*0102    itrs=scanforfiles(fname);
                0103    disp([ sprintf('Reading %i time levels:',size(itrs,2)) sprintf(' %i',itrs) ]);
29a315c92f Alis*0104   elseif isinf(arg)
37662e6918 Jean*0105    itrs=scanforfiles(fname);
                0106    if isempty(itrs)
                0107     AA=[];itrs=[];return;
2f1302ba5d Alis*0108    end
37662e6918 Jean*0109    disp([ sprintf('Found %i time levels, reading %i',size(itrs,2),itrs(end)) ]);
                0110    itrs=itrs(end);
22586ff4d2 Alis*0111 % elseif prod(double(arg>=0)) & prod(double(round(arg)==arg))
60f503df53 Jean*0112 % elseif prod(arg>=0) & prod(round(arg)==arg)
                0113   elseif min(arg)>=0 & isempty(find(round(arg)~=arg))
29a315c92f Alis*0114    if arg>=9999999999
                0115     error(sprintf('Argument %i > 9999999999',arg))
                0116    end
37662e6918 Jean*0117    itrs=arg;
aa4af8ca63 Jean*0118   elseif length(arg) == 1 & arg == -1
                0119    itrs=arg;
00398243b4 Alis*0120   else
                0121    error(sprintf('Argument %i must be a positive integer',arg))
                0122   end
96629ec1c2 Alis*0123   else
                0124    error('Multiple iterations should be specified as a vector')
                0125   end
a15d675b90 Alis*0126  end
00398243b4 Alis*0127 end
                0128 
37662e6918 Jean*0129 if isempty(itrs)
                0130  itrs=-1;
96629ec1c2 Alis*0131 end
                0132 
00398243b4 Alis*0133 % Loop over each iteration
37662e6918 Jean*0134 for iter=1:size(itrs,2);
                0135  if itrs(iter)>=0
                0136   fname=sprintf('%s.%10.10i',fnamearg,itrs(iter));
                0137  end
a15d675b90 Alis*0138 
bf3d93be62 Alis*0139 % Figure out if there is a path in the filename
37662e6918 Jean*0140  NS=findstr('/',fname);
                0141  if size(NS)>0
                0142   Dir=fname(1:NS(end));
                0143  else
                0144   Dir='./';
                0145  end
bf3d93be62 Alis*0146 
a15d675b90 Alis*0147 % Match name of all meta-files
aa4af8ca63 Jean*0148   %fprintf(' search for file "%s".*meta\n',fname);
                0149  allfiles=dir( sprintf('%s.*meta',fname) );
a15d675b90 Alis*0150 
37662e6918 Jean*0151  if size(allfiles,1)==0
aa4af8ca63 Jean*0152   disp(sprintf('No files match the search: %s.*meta',fname));
37662e6918 Jean*0153  %allow partial reads%  error('No files found.')
                0154  end
a15d675b90 Alis*0155 
                0156 % Loop through allfiles
37662e6918 Jean*0157  for j=1:size(allfiles,1);
aa4af8ca63 Jean*0158   %fprintf(' file # %3i : %s\n',j,allfiles(j).name);
a15d675b90 Alis*0159 
                0160 % Read meta- and data-file
37662e6918 Jean*0161   [A,N,M,mG] = localrdmds([Dir allfiles(j).name],ieee,recnum);
e6cfd13877 Jean*0162 
                0163 %- Merge local Meta file content (M) to MM string:
37662e6918 Jean*0164   if j > 0, %- to comment out this block: "if j < 0" (since j is always > 0)
                0165    ii=findstr(M,' timeStepNumber');
                0166    if isempty(ii), ii1=0; ii2=0;
                0167    else ii1=ii; ii2=ii+min(findstr(M(1+ii:end),'];')); end
                0168    ii=findstr(M,' timeInterval');
                0169    if isempty(ii), jj1=0; jj2=0;
                0170    else jj1=ii; jj2=ii+min(findstr(M(1+ii:end),'];')); end
                0171    if iter==1 & j==1,
aa4af8ca63 Jean*0172     MM=M; ind1=0; ind2=0; is1=ii1; js1=jj1; M3='';
                0173     if ii1*jj1 > 0,
37662e6918 Jean*0174      %ind1=min(ii1,jj1); ind2=max(ii2,jj2);
                0175      %if ii1 < jj1, ii3=ii2+1; jj3=jj1-1;
                0176      %else  ii3=jj2+1; jj3=ii1-1; end
                0177       order=sort([ii1 ii2 jj1 jj2]);
                0178       ind1=order(1); ii3=order(2)+1; jj3=order(3)-1; ind2=order(4);
                0179       M2=M(ii1:ii2); M4=M(jj1:jj2); M3=M(ii3:jj3);
aa4af8ca63 Jean*0180     elseif ii1 > 0,
37662e6918 Jean*0181       ind1=ii1; ind2=ii2;
                0182       M2=M(ii1:ii2); M4='';
aa4af8ca63 Jean*0183     elseif jj1 > 0,
37662e6918 Jean*0184       ind1=jj1; ind2=jj2;
                0185       M4=M(jj1:jj2); M2='';
                0186     end
                0187     M5=M(1+ind2:end);
                0188     %fprintf(' ii1,ii2 = %i %i ; jj1,jj2= %i %i ;', ii1,ii2, jj1,jj2);
                0189     %fprintf(' ii3,jj3= %i %i ; ind1,ind2= %i %i\n',ii3,jj3,ind1,ind2);
                0190     %fprintf('M(1:ind1)=%s<\n',M(1:ind1));
                0191     %fprintf(' M2=%s<\n',M2);
                0192     %fprintf(' M3=%s<\n',M3);
                0193     %fprintf(' M4=%s<\n',M4);
                0194     %fprintf(' M5=%s<\n',M5);
                0195    else
aa4af8ca63 Jean*0196     if ii1*jj1 > 0,
37662e6918 Jean*0197          order=sort([ii1 ii2 jj1 jj2]);
                0198          ind=order(1); ii3=order(2)+1; jj3=order(3)-1; ind2=order(4);
                0199     else ind=max(ii1,jj1); ind2=ii2+jj2; end
                0200     compar=(ind == ind1);    ii=0;
                0201     if compar & ind1 == 0,   ii=1; compar=strcmp (MM,M); end
                0202     if compar & ind1 > 0,    ii=2; compar=strncmp(MM,M,ind1) ; end
                0203     if compar & ind1 > 0,    ii=3; compar=strcmp(M5,M(1+ind2:end)); end
                0204     if compar & ii1*jj1 > 0, ii=4; compar=strcmp(M3,M(ii3:jj3)); end
                0205     if ~compar,
                0206      fprintf('WARNING: Meta file (%s) is different (%i) from 1rst one:\n', ...
                0207               allfiles(j).name,ii);
                0208      fprintf(' it=%i :MM:%s\n',itrs(1),MM);
                0209      fprintf(' it=%i :M :%s\n\n',itrs(iter),M);
                0210     elseif ind1 > 0,
                0211      if ii1 > 0,
                0212       Mj=M(ii1:ii2); ii=findstr(Mj,'['); Mj=Mj(1+ii:end);
e6cfd13877 Jean*0213 %   add it-number from Mj to M2 (if different):
37662e6918 Jean*0214       if isempty(findstr(M2,Mj)), M2=[deblank(M2(1:end-1)),Mj]; end
                0215      end
                0216      if jj1 > 0,
                0217       Mj=M(jj1:jj2); ii=findstr(Mj,'['); Mj=Mj(1+ii:end);
                0218 %   add time interval from Mj to M4 (if different):
                0219       if isempty(findstr(M4,Mj)), M4=[deblank(M4(1:end-1)),' ;',Mj]; end
                0220      end
                0221     end
                0222    end
e6cfd13877 Jean*0223 %  save modifications:
aa4af8ca63 Jean*0224    if ind1>0 & j==size(allfiles,1) & iter==size(itrs,2),
                0225      if ii1 < jj1, MM=[MM(1:ind1-1),M2,M3,M4,M5];
                0226      else          MM=[MM(1:ind1-1),M4,M3,M2,M5]; end
37662e6918 Jean*0227    end
                0228   end
a15d675b90 Alis*0229 
e6cfd13877 Jean*0230 %- put local data file content in global array AA:
37662e6918 Jean*0231   bdims=N(1,:);
                0232   r0=N(2,:);
                0233   rN=N(3,:);
                0234   ndims=prod(size(bdims));
                0235   if j==1 & iter==1, AA=zeros([bdims size(itrs,2)]); end
                0236   if mG(1)==0 & mG(2)==1,
                0237     if     (ndims == 1)
                0238      AA(r0(1):rN(1),iter)=A;
                0239     elseif (ndims == 2)
                0240      AA(r0(1):rN(1),r0(2):rN(2),iter)=A;
                0241     elseif (ndims == 3)
                0242      AA(r0(1):rN(1),r0(2):rN(2),r0(3):rN(3),iter)=A;
                0243     elseif (ndims == 4)
                0244      AA(r0(1):rN(1),r0(2):rN(2),r0(3):rN(3),r0(4):rN(4),iter)=A;
                0245     elseif (ndims == 5)
                0246      AA(r0(1):rN(1),r0(2):rN(2),r0(3):rN(3),r0(4):rN(4),r0(5):rN(5),iter)=A;
                0247     else
                0248      error('Dimension of data set is larger than currently coded. Sorry!')
                0249     end
                0250   elseif     (ndims == 1)
                0251      AA(r0(1):rN(1),iter)=A;
128322677f Jean*0252   else
                0253 %- to debug: do simple stransfert (with a loop on 2nd index);
                0254 %  will need to change this later, to improve efficiency:
37662e6918 Jean*0255    for i=0:rN(2)-r0(2),
                0256     if (ndims == 2)
                0257      AA(r0(1)+i*mG(1):rN(1)+i*mG(1),r0(2)+i*mG(2),iter)=A(:,1+i);
                0258     elseif (ndims == 3)
                0259      AA(r0(1)+i*mG(1):rN(1)+i*mG(1),r0(2)+i*mG(2), ...
128322677f Jean*0260                                     r0(3):rN(3),iter)=A(:,1+i,:);
37662e6918 Jean*0261     elseif (ndims == 4)
                0262      AA(r0(1)+i*mG(1):rN(1)+i*mG(1),r0(2)+i*mG(2), ...
128322677f Jean*0263                         r0(3):rN(3),r0(4):rN(4),iter)=A(:,1+i,:,:);
37662e6918 Jean*0264     elseif (ndims == 5)
                0265      AA(r0(1)+i*mG(1):rN(1)+i*mG(1),r0(2)+i*mG(2), ...
128322677f Jean*0266             r0(3):rN(3),r0(4):rN(4),r0(5):rN(5),iter)=A(:,1+i,:,:,:);
37662e6918 Jean*0267     else
                0268      error('Dimension of data set is larger than currently coded. Sorry!')
                0269     end
                0270    end
128322677f Jean*0271   end
a15d675b90 Alis*0272 
37662e6918 Jean*0273  end % files
00398243b4 Alis*0274 end % iterations
a15d675b90 Alis*0275 
                0276 %-------------------------------------------------------------------------------
                0277 
128322677f Jean*0278 function [A,N,M,map2glob] = localrdmds(fname,ieee,recnum)
a15d675b90 Alis*0279 
                0280 mname=strrep(fname,' ','');
                0281 dname=strrep(mname,'.meta','.data');
                0282 
128322677f Jean*0283 %- set default mapping from tile to global file:
                0284 map2glob=[0 1];
                0285 
a15d675b90 Alis*0286 % Read and interpret Meta file
                0287 fid = fopen(mname,'r');
                0288 if (fid == -1)
                0289  error(['File' mname ' could not be opened'])
                0290 end
                0291 
                0292 % Scan each line of the Meta file
                0293 allstr=' ';
                0294 keepgoing = 1;
                0295 while keepgoing > 0,
                0296  line = fgetl(fid);
                0297  if (line == -1)
                0298   keepgoing=-1;
                0299  else
                0300 % Strip out "(PID.TID *.*)" by finding first ")"
                0301 %old  ind=findstr([line ')'],')'); line=line(ind(1)+1:end);
                0302   ind=findstr(line,')');
                0303   if size(ind) ~= 0
                0304     line=line(ind(1)+1:end);
                0305   end
                0306 % Remove comments of form //
fe3c47efa0 Jean*0307   line=[line,' //']; ind=findstr(line,'//'); line=line(1:ind(1)-1);
                0308 % Add to total string (without starting & ending blanks)
                0309   while line(1:1) == ' ', line=line(2:end); end
128322677f Jean*0310   if strncmp(line,'map2glob',8), eval(line);
                0311   else allstr=[allstr,deblank(line),' '];
                0312   end
a15d675b90 Alis*0313  end
                0314 end
                0315 
                0316 % Close meta file
                0317 fclose(fid);
                0318 
                0319 % Strip out comments of form /* ... */
                0320 ind1=findstr(allstr,'/*'); ind2=findstr(allstr,'*/');
                0321 if size(ind1) ~= size(ind2)
                0322  error('The /* ... */ comments are not properly paired')
                0323 end
                0324 while size(ind1,2) > 0
fe3c47efa0 Jean*0325  allstr=[deblank(allstr(1:ind1(1)-1)) allstr(ind2(1)+2:end)];
                0326 %allstr=[allstr(1:ind1(1)-1) allstr(ind2(1)+3:end)];
a15d675b90 Alis*0327  ind1=findstr(allstr,'/*'); ind2=findstr(allstr,'*/');
                0328 end
                0329 
                0330 % This is a kludge to catch whether the meta-file is of the
                0331 % old or new type. nrecords does not exist in the old type.
52e6460920 Alis*0332 nrecords = NaN;
a15d675b90 Alis*0333 
e6cfd13877 Jean*0334 %- store the full string for output:
                0335 M=strrep(allstr,'format','dataprec');
                0336 
a15d675b90 Alis*0337 % Everything in lower case
                0338 allstr=lower(allstr);
                0339 
                0340 % Fix the unfortunate choice of 'format'
                0341 allstr=strrep(allstr,'format','dataprec');
                0342 
                0343 % Evaluate meta information
                0344 eval(allstr);
                0345 
                0346 N=reshape( dimlist , 3 , prod(size(dimlist))/3 );
e6cfd13877 Jean*0347 rep=[' dimList = [ ',sprintf('%i ',N(1,:)),']'];
96629ec1c2 Alis*0348 if ~isnan(nrecords) & nrecords > 1 & isempty(recnum)
3323410f65 Alis*0349  N=[N,[nrecords 1 nrecords]'];
96629ec1c2 Alis*0350 elseif ~isempty(recnum) & recnum>nrecords
                0351  error('Requested record number is higher than the number of available records')
                0352 end
                0353 
e6cfd13877 Jean*0354 %- make "dimList" shorter (& fit output array size) in output "M":
                0355  pat=' dimList = \[(\s*\d+\,?)*\s*\]';
                0356  M=regexprep(M,pat,rep);
                0357 %  and remove double space within sq.brakets:
                0358 ind1=findstr(M,'['); ind2=findstr(M,']');
                0359 if length(ind1) == length(ind2),
                0360  for i=length(ind1):-1:1, if ind1(i) < ind2(i),
                0361   M=[M(1:ind1(i)),regexprep(M(ind1(i)+1:ind2(i)-1),'(\s+)',' '),M(ind2(i):end)];
                0362  end; end
                0363 else error('The [ ... ] brakets are not properly paired')
                0364 end
                0365 
96629ec1c2 Alis*0366 if isempty(recnum)
                0367  recnum=1;
3323410f65 Alis*0368 end
a15d675b90 Alis*0369 
52e6460920 Alis*0370 if isnan(nrecords)
a15d675b90 Alis*0371 % This is the old 'meta' method that used sequential access
                0372 
                0373 A=allstr;
                0374 % Open data file
                0375 fid=fopen(dname,'r',ieee);
                0376 
                0377 % Read record size in bytes
                0378 recsz=fread(fid,1,'uint32');
                0379 ldims=N(3,:)-N(2,:)+1;
                0380 numels=prod(ldims);
                0381 
                0382 rat=recsz/numels;
                0383 if rat == 4
                0384  A=fread(fid,numels,'real*4');
                0385 elseif rat == 8
                0386  A=fread(fid,numels,'real*8');
                0387 else
                0388  sprintf(' Implied size in meta-file = %d', numels )
                0389  sprintf(' Record size in data-file = %d', recsz )
                0390  error('Ratio between record size and size in meta-file inconsistent')
                0391 end
                0392 
                0393 erecsz=fread(fid,1,'uint32');
                0394 if erecsz ~= recsz
                0395  sprintf('WARNING: Record sizes at beginning and end of file are inconsistent')
                0396 end
                0397 
                0398 fclose(fid);
                0399 
                0400 A=reshape(A,ldims);
                0401 
                0402 else
                0403 % This is the new MDS format that uses direct access
                0404 
                0405  ldims=N(3,:)-N(2,:)+1;
96629ec1c2 Alis*0406  for r=1:size(recnum(:),1);
a15d675b90 Alis*0407  if dataprec == 'float32'
96629ec1c2 Alis*0408   A(:,r)=myrdda(dname,ldims,recnum(r),'real*4',ieee);
a15d675b90 Alis*0409  elseif dataprec == 'float64'
96629ec1c2 Alis*0410   A(:,r)=myrdda(dname,ldims,recnum(r),'real*8',ieee);
a15d675b90 Alis*0411  else
                0412   error(['Unrecognized format in meta-file = ' format]);
                0413  end
96629ec1c2 Alis*0414  end
                0415 
                0416  A=reshape(A,[ldims size(recnum(:),1)]);
                0417  if size(recnum(:),1)>1
                0418   N(1,end+1)=size(recnum(:),1);
                0419   N(2,end)=1;
                0420   N(3,end)=size(recnum(:),1);
                0421  end
a15d675b90 Alis*0422 
                0423 end
                0424 
                0425 %-------------------------------------------------------------------------------
                0426 
                0427 % result = RDDA( file, dim, irec [options] )
                0428 %
                0429 % This routine reads the irec'th record of shape 'dim' from the
                0430 % direct-access binary file (float or double precision) 'file'.
                0431 %
                0432 % Required arguments:
                0433 %
                0434 %   file  - string  - name of file to read from
                0435 %   dim   - vector  - dimensions of the file records and the resulting array
                0436 %   irec  - integer - record number in file in which to write data
                0437 %
                0438 % Optional arguments (must appear after the required arguments):
                0439 %   prec  - string  - precision of storage in file. Default = 'real*8'.
                0440 %   ieee  - string  - IEEE bit-wise representation in file. Default = 'b'.
                0441 %
                0442 % 'prec' may take the values:
                0443 %       'real*4' - floating point, 32 bits.
                0444 %       'real*8' - floating point, 64 bits - the efault.
                0445 %
                0446 % 'ieee' may take values:
                0447 %    'ieee-be'     or 'b' - IEEE floating point with big-endian
                0448 %                           byte ordering - the default
                0449 %    'ieee-le'     or 'l' - IEEE floating point with little-endian
                0450 %                           byte ordering
                0451 %    'native'      or 'n' - local machine format
                0452 %    'cray'        or 'c' - Cray floating point with big-endian
                0453 %                           byte ordering
                0454 %    'ieee-le.l64' or 'a' - IEEE floating point with little-endian
                0455 %                           byte ordering and 64 bit long data type
                0456 %    'ieee-be.l64' or 's' - IEEE floating point with big-endian byte
                0457 %                           ordering and 64 bit long data type.
                0458 %
                0459 % eg.   T=rdda('t.data',[64 64 32],1);
                0460 %       T=rdda('t.data',[256],4,'real*4');
                0461 %       T=rdda('t.data',[128 64],2,'real*4','b');
                0462 function [arr] = myrdda(file,N,k,varargin)
                0463 
                0464 % Defaults
                0465 WORDLENGTH=8;
                0466 rtype='real*8';
                0467 ieee='b';
                0468 
                0469 % Check optional arguments
                0470 args=char(varargin);
                0471 while (size(args,1) > 0)
                0472  if deblank(args(1,:)) == 'real*4'
                0473   WORDLENGTH=4;
                0474   rtype='real*4';
                0475  elseif deblank(args(1,:)) == 'real*8'
                0476   WORDLENGTH=8;
                0477   rtype='real*8';
                0478  elseif deblank(args(1,:)) == 'n' | deblank(args(1,:)) == 'native'
                0479   ieee='n';
                0480  elseif deblank(args(1,:)) == 'l' | deblank(args(1,:)) == 'ieee-le'
                0481   ieee='l';
                0482  elseif deblank(args(1,:)) == 'b' | deblank(args(1,:)) == 'ieee-be'
                0483   ieee='b';
                0484  elseif deblank(args(1,:)) == 'c' | deblank(args(1,:)) == 'cray'
                0485   ieee='c';
                0486  elseif deblank(args(1,:)) == 'a' | deblank(args(1,:)) == 'ieee-le.l64'
                0487   ieee='a';
                0488  elseif deblank(args(1,:)) == 's' | deblank(args(1,:)) == 'ieee-be.l64'
                0489   ieee='s';
                0490  else
                0491   error(['Optional argument ' args(1,:) ' is unknown'])
                0492  end
                0493  args=args(2:end,:);
                0494 end
                0495 
                0496 nnn=prod(N);
                0497 
                0498 [fid mess]=fopen(file,'r',ieee);
                0499 if fid == -1
                0500  error('Error while opening file:\n%s',mess)
                0501 end
                0502 st=fseek(fid,nnn*(k-1)*WORDLENGTH,'bof');
                0503 if st ~= 0
                0504  mess=ferror(fid);
                0505  error('There was an error while positioning the file pointer:\n%s',mess)
                0506 end
                0507 [arr count]=fread(fid,nnn,rtype);
                0508 if count ~= nnn
                0509  error('Not enough data was available to be read: off EOF?')
                0510 end
                0511 st=fclose(fid);
96629ec1c2 Alis*0512 %arr=reshape(arr,N);
356a4fc521 Alis*0513 
                0514 %
37662e6918 Jean*0515 function [itrs] = scanforfiles(fname)
356a4fc521 Alis*0516 
37662e6918 Jean*0517 itrs=[];
fe06077190 Alis*0518 allfiles=dir([fname '.*.001.001.meta']);
29a315c92f Alis*0519 if isempty(allfiles)
fe06077190 Alis*0520  allfiles=dir([fname '.*.meta']);
29a315c92f Alis*0521  ioff=0;
f4f3cf4445 Alis*0522 else
                0523  ioff=8;
29a315c92f Alis*0524 end
356a4fc521 Alis*0525 for k=1:size(allfiles,1);
                0526  hh=allfiles(k).name;
37662e6918 Jean*0527  itrs(k)=str2num( hh(end-14-ioff:end-5-ioff) );
356a4fc521 Alis*0528 end
37662e6918 Jean*0529 itrs=sort(itrs);