EasySpin forum

Code: Select all

function varargout=d01read(varargin)

%  d01read Read data from .d01/.exp SpecMan data files
%
%   data = d01read(filename,...)
%   [ax,data] = d01read(filename,...);
%   [ax,data,dsc] = d01read(filename,...);
%
%   Reads data and axis information from Bruker
%   files  DSC/DTA  and  PAR/SPC.   Aditional 
%   arguments  are  always  coming  in pairs 
%   field, value. Fields are:
%     'Dataset', [real dataset, imag dataset] 
%   ax contains fields x,y,z  depends on data 
%   dimension and may contain other fields as 
%   well. dsc is a cell array  of description 
%   strings.

% Original code by Boris Epel & Alexey Silakov
% MPI of Bioinorganic Chemistry, Muelhaim an der Ruhr, 2003
% Free for non-commercial use. Use the program at your own risk. 
% The authors retains all rights. 
% Contact: epel@mpi-muelheim.mpg.de


if ~nargin
    [name,path] = uigetfile('*.d01','Load Specman File');
    if ~ischar(name)
        varargout(1:nargout) = {[]};
        return;
    end
    fname = fullfile(path,name);
    dscname=strrep(fname, 'd01', 'exp');
else
   [path,name] = fileparts(varargin{1}); 
   fname = fullfile(path,[name,'.d01']);
   dscname = fullfile(path,[name,'.exp']);

end


fid=fopen(char(fname),'r', 'ieee-le');
if fid<1, error(['File ''',fname,''' can not be open for read.']);end

ndim1=fread(fid, 1,'uint32');       % number of headers, re/im etc.
dformat=fread(fid,1,'uint32');     % format:0-double,1-float
if dformat==1
  sformat='float32';
else
  sformat='double';
end

dstrms = {};
ntotal = 1;

for k=1:ndim1
  ndim2       = fread(fid,1,'int32');
  dstrms{end+1}.dim = fread(fid,4,'int32');
  dstrms{end}.dim(ndim2+1:end) = 1;
  dstrms{end}.first = ntotal;
  dstrms{end}.total = fread(fid,1,'int32');
  ntotal = ntotal + dstrms{end}.total;
end
tmpdat=fread(fid,ntotal,sformat);
fclose(fid);

switch(ndim1)
    case 0,
        error('No data present');
    case 2,
        spec=tmpdat(dstrms{1}.first:(dstrms{1}.first+dstrms{1}.total-1))+...
            1i*tmpdat((dstrms{2}.first:dstrms{2}.first+dstrms{2}.total-1));
        spec=reshape(spec,dstrms{1}.dim');
    case 1,
        spec=reshape(tmpdat,dstrms{1}.dim');
    otherwise,
        %% find if all data have the same dimensions  
        dim = dstrms{1}.dim;
        isthesame = 1;
        is1D      = (sum(dim~=1) == 1);
        for k=2:ndim1
            if sum(dim == dstrms{k}.dim) ~= 4,
                isthesame = false;
                break;
            elseif is1D && sum(dstrms{k}.dim~=1)~=1 
                is1D = false;    
            end
        end
        if isthesame && is1D,
            % read all as columns
            xdim = dim(dim~=1);
            spec=reshape(tmpdat, xdim, ndim1);
        else
            spec=tmpdat;
        end
end

dsc = SpecMandsc(dscname);
ax = SpecManpar(dsc);

if ~isfield(ax, 'x') || size(ax.x, 1)~=size(spec, 1)
  ax.x = [1:size(spec, 1)];
end
ax.type = 'data';
if isfield(dsc,'general_freq1')
    ax.freq1 = kvgetvalue(dsc.general_freq1);
end

% assign output depending on number of output arguments
% requested
switch nargout
 case 1,
   varargout = {spec};
 case 2,
   varargout = {ax, spec};
 case 3,
   varargout = {ax, spec, dsc};
end
return

function par = SpecMandsc(filename)
h = fopen(filename);
if h<0
  disp('Description file was not found.');
  par = [];
  return;
end
olda = '';
par = [];
forbidden = '~!@#$%^&*()./\';
section = '';
text = 1; prg = 1;
while feof(h)<1
  s = strtrim(fgetl(h));
  if isempty(s), continue; end;
  sect = find(s=='[' | s==']');
  %   this is a section header   
  if size(sect, 2)==2 && sect(1)==1
    section = s(sect(1)+1:sect(2)-1);
%     par = setfield(par, section, 'section');
  else
    switch section
      case 'text'
        par = setfield(par, ['text', num2str(text)], s);
        text = text + 1;
      case 'program'
        par = setfield(par, ['prg', num2str(prg)], s);
        prg = prg + 1;
      otherwise
        [a,s]=strtok(s, '=');
        a = strtrim(a);
        a(a=='/' | a=='\' | a==' ')='_';
        par = setfield(par, [section,'_',a], s(2:end));
    end  
  end
end
fclose(h);
return

function res = SpecManpar(par)

prefix = ['n', 'u', 'm', 'k', 'M', 'G'];
koeff  = [1E-9, 1E-6, 1E-3, 1E3, 1E6, 1E9];
res.title = safeget(par, 'general_name', '?');
sweepax = {};
key = 'transient';
fullfield = ['sweep_', key];
idx = 0;
triggers = str2num(safeget(par, 'streams_triggers', '1'));
while isfield(par, fullfield)
    [ax.t, str] = strtok(getfield(par, fullfield), ',');
    ax.t = strtrim(ax.t);
    [ax.size, str] = strtok(str(2:end), ',');
    if ax.t=='S' || ax.t=='I' || ax.t=='A' || ax.t=='R', ax.size = 1; else ax.size = str2num(ax.size); end
    [ax.reps, str] = strtok(str(2:end), ',');
    ax.reps = str2num(ax.reps);
    ax.var = {};
    while ~isempty(str)
        [ax.var{end+1}, str] = strtok(str(2:end), ',');
    end
    sweepax{end+1,1} = ax;
    fullfield = ['sweep_sweep', num2str(idx)];
    idx = idx +1;
end
sweepax{1}.size=sweepax{1}.size*triggers;
res.sweepax = sweepax;

axislabel = ['xyz'];
counter = 1;
for k = 1:size(sweepax, 1)
    arr = [];
    asize = sweepax{k}.size;
    if asize > 1
        switch sweepax{k}.t
            case {'I', 'A'}
                tempparam = 'trans';
                par.params_trans = '1sl step 1sl;';
            case 'T'
                tempparam = 'trans';
                dwell_time_str = safeget(par, 'streams_dwelltime', '1 ns');
                dwell_time_str = strtrim(gettoken(dwell_time_str, ','));
                par.params_trans = ['0 ns step ', dwell_time_str,';'];
            otherwise
                tempparam = sweepax{k}.var{1};
                tempparam(strfind(tempparam, ' ')) = '_';
        end
        % check if this is a parameter
        if isfield(par, ['params_', tempparam])
            str = getfield(par, ['params_', tempparam]);
            if ~isempty(strfind(str,'step'))
                [tk1, str1] = gettoken(str, 'step');
                % string of the type 10ns step 6 ns
                tk2 = strtrim(gettoken(str1, ';'));
                [minval, unit] = kvgetvalue(tk1);
                step = kvgetvalue(tk2);
                arr = [0:asize-1]*step+minval;
            elseif ~isempty(strfind(str,'logto'))
                [tk1, str1] = gettoken(str, 'logto');
                % string of the type 10ns logto 60 ns
                tk2 = strtrim(gettoken(str1, ';'));
                [minval, unit] = kvgetvalue(tk1);
                maxval = kvgetvalue(tk2);
                arr = logspace(log10(minval), log10(maxval),asize);
            elseif ~isempty(strfind(str,'to'))
                [tk1, str1] = gettoken(str, 'to');
                % string of the type 10ns to 60 ns
                tk2 = strtrim(gettoken(str1, ';'));
                [minval, unit] = kvgetvalue(tk1);
                maxval = kvgetvalue(tk2);
                arr = [0:1/(asize-1):1]*(maxval-minval)+minval;
            else
                % string of the type 10ns, 20ns, 30ns;
                [str1] = gettoken(str, ';');
                [tk1, str1] = gettoken(str1, ',');
                while ~isempty(tk1)
                    [arr(end+1),unit] = kvgetvalue(tk1);
                    [tk1, str1] = gettoken(str1, ',');
                    if isempty(tk1) && ~isempty(str1)
                        tk1 = str1; str1 = [];
                    end
                end
            end
        else
            str = getfield(par, ['aquisition_', tempparam]);
            arr = [0:asize-1]';
            unit = 's';
        end
        
        % Unit normalization
        switch unit
        case 'G',
        case 'K',
        case 's',
        otherwise
        umax = max(abs(arr));
        for kk = length(koeff):-1:1
            if umax > koeff(kk)
                uk = koeff(kk);
                unit = [prefix(kk), unit];
                arr = arr./uk;
                break;
            end
        end
        end

        res = setfield(res, axislabel(counter), arr');
        res = setfield(res, [axislabel(counter), 'label'], ...
            [sweepax{k}.var{1}, ', ',unit]);
        counter = counter + 1;
    end
end
return

function [tk,rstr] = gettoken(istr, tok)

pos = strfind(istr, tok);

if isempty(pos)
    tk=istr;
    rstr='';
else
    tk=strtrim(istr(1:pos-1));
    rstr=strtrim(istr(pos+length(tok):end));
end
return

function res = safeget(strct, fld, deflt)
% function res = safeget(strct, fld, deflt)
% Returns the field of the structure or default 
% value if field is absent.
% if default value is array then return value
% has not less elements than in this array 
% (not in char case)

if isfield(strct, fld)
    res = getfield(strct, fld);
    sd = size(deflt, 2);
    sr = size(res, 2);
    if sr < sd
        if iscell(deflt)
            [res{sr+1:sd}] = deal(deflt{sr+1:sd});
        elseif ~ischar(deflt)
            res(sr+1:sd) = deflt(sr+1:sd);
        end
        
    end
else
    res = deflt;
end
return

function [val, unit, pref, pref_val] = kvgetvalue(str)
% KVGETVALUE read string value in ci-standard units

% [val, str_unit, str_koefficient] = kvgetvalue(str)


prefix = ['p','n', 'u', 'm', 'k', 'M', 'G', 'T'];
koeff  = [1E-12, 1E-9, 1E-6, 1E-3, 1E3, 1E6, 1E9, 1E12];
idx = (str >= '0' & str <='9') | str == '.' | ...
    upper(str) == 'E' | str == '+' | str == '-';
pref = '';
pref_val = 1;
val = str2num(str(idx));
unit = str(~idx);
unit = unit(unit~=' ');
if length(unit) > 1
    if ~isempty(unit)
%         kk = findstr(prefix, unit(1));
          kk = strfind(prefix, unit(1));
        if ~isempty(kk)
            val = val * koeff(kk);
            unit = unit(2:end);
            pref = prefix(kk);
            pref_val = koeff(kk);
        end
    end
end
return