function [out,out1,out2,out3]=fe_curve(varargin)

%FE_CURVE Curve handling and signal processing utilities 
%
%Supported curve handling  commands are
%
%           out=fe_curve('DataTypes');
%           out=fe_curve('DataTypes',DesiredType);
%           out=fe_curve('GetCurve',model,curve_name);
%           out=fe_curve('GetId #Id',model);
%           model=fe_curve(model,'Set',name,curve_def);
%           out=fe_curve('read /path/filename.cyt');
%           Y = fe_curve('returny',curve,X)
%
%           curve=fe_curve('getcurve',model,curve_name)
%           list=fe_curve('getcurve',model)
%           model=fe_curve(model,'set',curve_name,data_structure)
%           model=fe_curve(model,'set',curve_name,string)
%
%        Plot :
%           fe_curve('plot',curve);
%           fe_curve('plot',handle,curve);
%           fe_curve('plot',model,curve_name);
%           fe_curve('plot',handle,model,curve_name);
%
%
%-------------------
%Signal processing  commands are
%-------------------
%
%           out=fe_curve('FFT',frame); % continuous fourier transform
%           out=fe_curve('TimeFreq',Input,Transfert);
%
%           out=fe_curve('BandPass Hz f_min f_max',frames);
%           out=fe_curve('BandPass Rd w_min w_max',frames);
%
%           out=fe_curve('window window_name',Nb_pts);
%
%           out=fe_curve('H1H2 input_channels',frame,window,weighing);
%                  window and weighing can be omitted
%
%           out=fe_curve('ResSpectrum Type',frames,freq,damp);
%             
%           out=fe_curve('Noise',Nw_pt,sf);
%           out=fe_curve('Noise',Nw_pt,sf,f_max);
%           out=fe_curve('Noise',Nw_pt,sf,Spec);
% 
%           out=fe_curve('TestFrame');
%
%           out=fe_curve('TestFunc');% Fun=cos, sin, etc.
%           out=fe_curve('TestFunc',t);%t:time vector, ex: t=linspace(0,1,199)

%	Etienne Balmes, Mathieu Corus
%       Copyright (c) 2001-2005 by INRIA and SDTools, All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license

if nargin==0; CAM=''; 
else CAM=varargin{1};  carg=2; 
end
if ~ischar(CAM)&isfield(CAM,'Stack');   model=CAM;  CAM=varargin{2}; carg=3; 
else;    model=[]; 
end

[CAM,Cam]=comstr(CAM,1);

RunOpt.ExtrapBy0=0;
if ~isempty(findstr(Cam,'-extrapby0')); RunOpt.ExtrapBy0=1; end


%---------------------------------------------------------------
if comstr(Cam,'types')
  
  
  out={'1 (t,f(t))','2 f(t+kdt)','3 Constant','4 Sinusoid','5 Staircase', ...
      '6','7'}; %xxx 
  %-- Please fill in with comments --%
  
  
%---------------------------------------------------------------
% Fail safe procedure to get the x axis values
elseif comstr(Cam,'getx')

 val=varargin{carg};carg=carg+1;
 if ~isfield(val,'X')|~isfield(val,'Y')
    error('Not a valid curve'); 
 end

 t=val.X;if iscell(t);t=t{1};end; if isstruct(t);t=t.data;end
 if size(t,1)==1&size(t,2)>1; t=t(:);end
 N=length(t); 
 if N~=size(val.Y,1) 
    error(sprintf('x-axis (%i) is not consistent with .Y rows (%i)', ...
    N,size(val.Y,1))); 
 end
 out=t;

 if nargout>1&comstr(Cam,'getxtime') % Build the frequency vector
  r1=diff(t);if std(r1)/mean(r1)>1e-10; 
   error('Can''t compute FFT for uneven time steps')
  end
  f=1/mean(r1)*[0:length(t)-1]/length(t);
  out1=f; out2=length(t);
 elseif nargout>1&comstr(Cam,'getxfreq') % Build the time vector
  f=t; r1=diff(f); if std(r1)/mean(r1)>1e-10; 
   error('Can''t compute FFT for uneven frequency steps')
  end
  if f(1)~=0; error('First frequency should be zero');end
  N=length(f); 
  if abs(val.Y(2,1)./conj(val.Y(N,1))-1)>1e-10; N=2*N-2;end
  t=[0:N-1]/N*f(end); out=t;out1=f;out2=length(t);

 end

% indices for spectrum halves and real values
elseif comstr(Cam,'getia')

 N=varargin{carg};carg=carg+1;
 if round(rem(N,2))==1; 
    ia=N:-1:(N+3)/2;ib=2:(N+1)/2;ind=1:(N+1)/2;ic=[1];
 else; 
   ia=N:-1:N/2+2;ib=2:N/2;ind=1:N/2+1;ic=[1 N/2+1];
 end
 out=struct('ia',ia,'ib',ib,'ic',ic,'ind',ind);


%---------------------------------------------------------------
elseif comstr(Cam,'get')

  st=Cam; [CAM,Cam]=comstr(CAM,4);
  if isempty(model)&nargin>=carg ; model=varargin{carg};carg=carg+1; end

  % get all curves (by name)
  if nargin<carg & comstr(Cam,'curve') 
   out={};
   out=stack_get(model,'curve'); out=out(:,3)';

   Case=fe_case('getcase',model);
   st={'DOFLoad','FVol','FSurf'};
   for j3=1:length(st)
    C1=stack_get(Case,st{j3});
    for j1=1:size(C1,1)
     if isfield(C1{j1,3},'curve')
        i2=[];
        for j2=1:length(C1{j1,3}.curve); 
         if ~ischar(C1{j1,3}.curve{j2});i2(end+1)=j2;end
        end
        out(end+[1:length(i2)])=C1{j1,3}.curve(i2);
     end
    end % j1
   end % j3
   return;
  end % return all curves


  if comstr(Cam,'curve') % 'GetCurve'
    
    curve_name=varargin{carg};carg=carg+1;
    if ischar(curve_name)
      if isempty(curve_name); out={}; return; end
      [out,i1]=stack_get(model,'curve',curve_name);
      if isempty(i1) ;
       try
        Case=fe_case('getcase',model); st={'DOFLoad','FVol','FSurf'};
        for j3=1:length(st)
          C1=stack_get(Case,st{j3});
          for j1=1:size(C1,1)
           if isfield(C1{j1,3},'curve')
            i2=[];
            for j2=1:length(C1{j1,3}.curve); 
             if ~ischar(C1{j1,3}.curve{j2});i2(end+1)=j2;end
            end
            out=C1{j1,3}.curve(i2);return;
           end
          end % j1
        end % j3
       catch
        error('Unknown curve name'); 
       end
      end 
      out1=out{1,2}; out=out{1,3};
    elseif isa(curve_name,'cell')
      out={}; out1={}; for j1=1:length(curve_name)
        if isempty(curve_name{j1}) r1=[];
        else [r1,st1]=fe_curve('getcurve',model,curve_name{j1});
        end
        if ~isempty(r1) out{end+1}=r1; out1{end+1}=st1;return;end
      end
    elseif isstruct(curve_name) & isfield(curve_name,'name')
      out=fe_curve('getcurve',model,curve_name.name);
    else
      out=fe_curve(model,sprintf('GetId %i',curve_name)); 
    end
    
  elseif comstr(Cam,'id') % 'GetId'
    
    i1=comstr(CAM(3:end),[-1]);
    if nargin>=carg  
    elseif isempty(model) error('You must give a model'); end %xxx
    if length(i1)==0 error('ID not provided'); end
    [out,i2]=stack_get(model,'curve');
    if isempty(out) out=[];
    else
      i2=[];
      for j1=1:size(out,1)
        if isfield(out{j1,3},'ID') & any(out{j1,3}.ID==i1) 
          i2(end+1)=j1;
        end
      end
      out=out(i2,3); if length(out)==1 out=out{1};end
    end
  else sdtw('''Get%s'' unknown',CAM);
  end
  
%---------------------------------------------------------------
% Check format consistency and reformat if needed
elseif comstr(Cam,'check');

curve=varargin{carg};carg=carg+1;

% SDT xf format 
if isfield(curve,'w')&isfield(curve,'xf');
 out=feutil('rmfield',curve,{'w','xf','dof'});
 out.X=iigui(curve.w); out.Y=iigui(curve.xf); 
 if isfield(curve,'dof'); out.data=iigui(curve.dof);
  if size(out.data,2)>=8%Universal file 58 : Z axis 
    out.Z=out.data(:,8);
  end
 end
else; out=curve;
end


%---------------------------------------------------------------
elseif comstr(Cam,'set'); [CAM,Cam]=comstr(CAM,4);

  curve_name=varargin{carg};carg=carg+1;
  curve=varargin{carg};carg=carg+1;

  if comstr(Cam,'loadcurve') % add curve in a specified load case
   % model=fe_curve(model,'set LoadCurve','Point load 1',2,'step 1e-4*10');
   data=varargin{carg};carg=carg+1;
   Case=fe_case('getcase',model);
   [r1,i1]=stack_get(Case,'DOFLoad',curve_name);
   r1{1,3}.curve{curve}=data;
   out = fe_case(model,r1{1,1},r1{1,2},r1{1,3});

  else % in add curve in model.Stack

   if ischar(curve)
     curve=struct('ID',1,'X',[],'Y',curve,'xunit',[],'yunit',[], ...
       'name',curve_name,'Interp','','PlotFcn','');
   elseif isstruct(curve) & isfield(curve,'name') ...
                          & ~comstr(curve.name,curve_name)
     sdtw('_nb',sprintf('Curve name changed from %s to %s',...
                         curve.name,curve_name));
     curve.name=curve_name;
   elseif isstruct(curve) & ~isfield(curve,'name')
     sdtw('_nb',sprintf('Curve name set to %s',curve_name));
     curve.name=curve_name;
   end
   out=stack_set(model,'curve',curve_name,curve);

  end 
%---------------------------------------------------------------
% These are datatypes defined in universal file format
elseif comstr(Cam,'datatype');[CAM,Cam]=comstr(CAM,9);
  
 if comstr(Cam,'fix') % go from old field format to cell array
  r1=varargin{carg};carg=carg+1;
  if isfield(r1,'label');
   out={r1.label r1.unit r1.type};
  else; out=r1;
  end
 else
  out=DataType(varargin{carg:end});
  if ~isempty(strfind(Cam,'cell')); out={out.label out.unit out.type};end
 end
% Analysis types are also defined in UFF 55
elseif comstr(Cam,'analysistype'); out=DataType(varargin{carg:end});
  
  out={0','Unknown';1,'Static';2,'Normal Mode'
    3,'Cplx eig. order 1'
    4,'Transient'
    5,'Frequency Response'
    6,'Buckling'
    7,'Cplx eig. order 2'};
  
% DOF Units - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'doflab');[CAM,Cam]=comstr(CAM,7);

if comstr(Cam,'load');
 st1={'Force','N',[13 0 1 0 0];
      'Moment','N/m',[100 -1 1 0 0]
      'd*S','m^3',[100 3 0 0 0]};
else
 st1={DataType('Displacement -cell');
      DataType('Rotation -cell')
      DataType('Pressure -cell')};
 for j1=1:size(st1,1);st1(j1,1:3)=st1{j1,1};end
end

out=num2cell([1:12 19]');
i1=[1 2 3 7 8 9];out(i1,2)=st1(1,1);out(i1,3)=st1(1,2);out(i1,4)=st1(1,3);
i1=[4 5 6 10 11 12];out(i1,2)=st1(2,1);out(i1,3)=st1(2,2);out(i1,4)=st1(2,3);
i1=13;out(i1,2)=st1(3,1);out(i1,3)=st1(3,2);out(i1,4)=st1(3,3);

if nargin>1; % Return dof labels if needed
 DOF=varargin{carg};carg=carg+1;
 if size(DOF,2)==1; 
  out1=fe_c(DOF);out1=out1(:);
  i2=[out{:,1}];nind=sparse(i2,1,1:length(i2),99,1);
  i2=nind(round(rem(DOF,1)*100));i3=find(i2);
  out1(i3,2:3)=out(i2(i3),3:4);
 elseif size(DOF,2)==5; % 5 column sensor format [Tag NodeId Dir[xyz]]
  try;
   out1=sprintf('%20.2f_(%20i)',DOF(:,1:2)');
   out1=strrep(out1,'.00','   ');
   out1=reshape(out1,43,length(out1)/43)';
   out1(:,find(all(out1==' ')))=''; out1(find(out1=='_'))=' ';
  catch;out1=num2str(DOF(:,2));
  end
  out1=cellstr(out1);  out1(:,2)={out{1,3}}; out1(:,3)={out{1,4}};
 end
 out=out1;
end
 
%---------------------------------------------------------------
% Basic plotting capabilities
elseif comstr(Cam,'plot')

 if sp_util('issdt');
   eval('out=ii_plp(varargin{:});');
   return;
 end  
 if ishandle(varargin{carg})
   ca=varargin{carg};carg=carg+1;
   if ~isempty(get(ca,'children'))
     delete(get(ca,'children'));
     axes(ca);set(ca,'visible','on','fontsize',7);
   end;
 else 
   figure;
   ca=gca;
 end
  
  curve=varargin{carg};carg=carg+1;
  ind_curve=[];
  
  if isfield(curve,'Stack')
    curve_name=varargin{carg};carg=carg+1;
    for i1=1:size(curve.Stack)
      if ischar(curve.Stack{i1,1}) & ischar(curve.Stack{i1,2})
        if comstr(curve.Stack{i1,1},'curve') & ...
            comstr(curve.Stack{i1,2},curve_name)
          ind_curve=[ind_curve i1];
        end;
      end;
    end;
    
    if isempty(ind_curve) error('Unknown curve name'); end
    
    if length(ind_curve)>1
      warning(['More than one curve named ' curve_name]);
      disp('First curve retained');
      ind_curve=ind_curve(1);
    end;
    curve=curve.Stack{ind_curve,3};
  end
  
  if isfield(curve,'PlotFcn')
    if ischar(curve.PlotFcn)
      eval(curve.PlotFcn);
    else error('Not expected');
    end
  else
    line(curve.X,curve.Y);
  end
  
 curve=fe_curve('datatype fix',curve);
 try;if isfield(curve,'xunit');
    st=curve.xunit{1,1}; 
    if ~strcmpi(curve.xunit{1,2},'none');st=[st ' ' curve.xunit{1,2}];end
    xlabel(st);
 end;end
  
 try; if isfield(curve,'yunit');
    st=curve.yunit{1,1}; 
    if ~strcmpi(curve.yunit{1,2},'none');st=[st ' ' curve.yunit{1,2}];end
    ylabel(st);
 end;end
  
 if isfield(curve,'name'); title(curve.name); end
  
  %---------------------------------------------------------------
  %---------------------------------------------------------------
  
  %-- Compute Time response to a signal given the transfer    
elseif comstr(Cam,'timefreq') [CAM,Cam]=comstr(CAM,5);
  
  val=varargin{carg};carg=carg+1;
  TF=varargin{carg};carg=carg+1;
  
  %% Fourier transform of input signal
  if isempty(val); 
   w=fe_curve('getx',TF); RunOpt=fe_curve('getia',length(w));  
   if abs(TF.Y(RunOpt.ia(1),1)./TF.Y(RunOpt.ib(1),1)-1)>1e-10
    N=2*length(w); w=[0:N-1]*mean(diff(w));
   end
   val=struct('X',w,'Y',ones(size(w)));
  else;w=fe_curve('getx',val);   N=length(w);   
  end
  RunOpt=fe_curve('getia',N);  
  if size(val.Y,1)==N; ffta=fft(val.Y,[],1);
  else;  ffta=fft(val.Y.',[],1);
  end;
  
  if isfield(TF,'xf');r1=r1.w; xf=TF.xf;
  elseif isfield(TF,'Y');r1=fe_curve('getx',TF);xf=TF.Y;
  else; r1=w; 
  end

  if size(xf,1)~=N;xf=xf.';end
  if ~any(size(xf)==N);error('Inconsistent Time Signal and Transfer');end
 
  step=mean(diff(w));  w=2*pi/step/N*[0:N-1];time=[0:(N-1)]'*step;
  
  
  Xt=zeros(size(xf));
  
  for j1=1:size(xf,2)
    X=ffta.*xf(:,j1);
    %-- Before IFFT of X padd with conjugate f(i+N/2+1)=conj(f(N/2-i+1))    
    X(RunOpt.ia)=conj(X(RunOpt.ib));  X(RunOpt.ic)=real(X(RunOpt.ic));
    Xt(:,j1)=real(ifft(X));
  end
  
  out=struct('X',time,'Y',Xt, ...
    'name','Time-Freq response');
  if sp_util('issdt');out=fe_def('fixcurvetime',out,TF);end  
  if nargout==0; fe_curve('plot',out);end
  
  
%---------------------------------------------------------------
%-- Response spectrum : see "Engineering Vibration", Daniel J. Inman
%--                         Prentice Hall, 1994           
elseif comstr(Cam,'resspectrum'); [CAM,Cam]=comstr(CAM,12);
  
  Stack={'Flexibility','H= 1./(s2+2*zeta*omj*s+omj^2)'
         'Admittance','H= s./(s2+2*zeta*omj*s+omj^2)'
         'Accelerance','H=s2./(s2+2*zeta*omj*s+omj^2)'
         'Disp transmis.','H=[2*zeta*omj*s+omj^2]./(s2+2*zeta*omj*s+omj^2)';
         'RDisp transmis.','H=s2./(s2+2*zeta*omj*s+omj^2)'};
  % xxx stiffness on disp input
  if nargin<4
    Stack=Stack';fprintf('%20s : %s\n',Stack{:}); return;
  end;
  
  val=varargin{carg};carg=carg+1;
  freq=varargin{carg};carg=carg+1;
  if carg<=nargin; damp=varargin{carg};carg=carg+1;
  else; damp=[1e-3 1-2 1e-1];
  end
  damp=damp(:); freq=freq(:);
 
  
  % use closest actual frequency to avoid leakage
  freq=sort(freq); 
  [t,F,N]=fe_curve('getxtime',val);
  for j1=1:length(freq); 
   [r1,i1]=min(abs(F-freq(j1)));freq(j1)=F(i1);
  end
  
  RunOpt=fe_curve('getia',N);
  %F(find(abs(F)<eps))=eps;
  s=i*2*pi*F(ind);s2=s.^2; U=fft(val.Y);U=U(ind);

  % derivatives of output signal
  RunOpt.Deriv=[]; RunOpt.tplot=0;
  i1=strfind(Cam,'-v'); 
  if ~isempty(i1);CAM(i1+[0:1])=''; RunOpt.Deriv=s; end
  Cam=comstr(CAM,-27);i1=strfind(Cam,'-a'); 
  if ~isempty(i1);CAM(i1+[0:1])=''; RunOpt.Deriv=s2; end
  Cam=comstr(CAM,-27);i1=strfind(Cam,'+v'); 
  if ~isempty(i1);CAM(i1+[0:1])=''; RunOpt.Deriv=1./s; end
  Cam=comstr(CAM,-27);i1=strfind(Cam,'+a'); 
  if ~isempty(i1);CAM(i1+[0:1])=''; RunOpt.Deriv=1./s2; end
  Cam=comstr(CAM,-27);i1=strfind(Cam,'-tplot');
  if ~isempty(i1);CAM(i1+[0:5])=''; RunOpt.tplot=1;end
  [CAM,Cam]=comstr(CAM,1);


  j1=strmatch(Cam(1:min(3,length(Cam))),comstr(Stack(:,1),-27));
  if isempty(j1);j1=1;end
  RunOpt.HFcn=horzcat(Stack{j1,2},';');

  out=struct('X',freq(:)*ones(1,length(damp)), ...
             'Z',ones(size(freq(:)))*damp(:)', ...
             'name','Stack{j1,1}', ...
             'xunit',{fe_curve('datatypecell','Frequency')});
  z=[]; 
  for j1=1:length(freq); for j2=1:length(damp)
    omj=2*pi*freq(j1); zeta=damp(j2);
    eval(RunOpt.HFcn); Y=H.*U; 
    if ~isempty(RunOpt.Deriv);Y=Y.*RunOpt.Deriv;end 
    Y(RunOpt.ia)=conj(Y(RunOpt.ib)); 
    Y(RunOpt.ic)=real(Y(RunOpt.ic)); y=real(ifft(Y));
    if RunOpt.tplot; z(:,end+1)=y;end
    out.Y(j1,j2)=max(abs(y));
  end;end
  % Debug for tplot
  if RunOpt.tplot; figure(100);clf; plot(val.X,z);iimouse; end

  if ~isempty(strfind(Cam,'pseudov'))  % Pseudo-velocity
    out.Y=out.Y.*freq(:,ones(1,size(out.Y,2)));
  elseif ~isempty(strfind(Cam,'pseudoa')) % Pseudo-acceleration
    f2=freq(:).^2; out.Y=out.Y.*f2(:,ones(1,size(out.Y,2)));
  end
   
%---------------------------------------------------------------
%  out =fe_curve('IFFT',frames); continuous inverse transform
elseif comstr(Cam,'ifft');[CAM,Cam]=comstr(CAM,5);
  
  frames=varargin{carg};carg=carg+1;
  [t,f,N]=fe_curve('getxfreq',frames);  
  RunOpt=fe_curve('getia',N);  
  if comstr(Cam,'shock'); RunOpt.Coef=diff(t(1:2));[CAM,Cam]=comstr(CAM,6);
  else; RunOpt.Coef=1/N;
  end

  out=frames;
  if ~iscell(out.X);
    out.X={t 1:size(frames.Y,2)}; out.Xlab={'Time','channel'};
  else; out.X{1}=t;
  end
  if ~isfield(out,'Xlab')|~isa(out.Xlab,'cell')
   out.Xlab={'Time','channel'};
  end
  r2=frames.Y;out.Y(length(t),1)=0;
  for j1=1:size(r2,3); 
   X=squeeze(frames.Y(:,:,j1));
   X(RunOpt.ia)=conj(X(RunOpt.ib,:));  X(RunOpt.ic,:)=real(X(RunOpt.ic,:));
   out.Y(:,:,j1)=real(ifft(X))/RunOpt.Coef;% compute continuous IFT
  end
 
  
  %---------------------------------------------------------------
  %  out =fe_curve('FFT',frames);
elseif comstr(Cam,'fft');[CAM,Cam]=comstr(CAM,4);
  
  frames=varargin{carg};carg=carg+1;
  [t,f,N]=fe_curve('getxtime',frames);  
  RunOpt=fe_curve('getia',N);  
  if comstr(Cam,'shock'); RunOpt.Coef=diff(t(1:2));[CAM,Cam]=comstr(CAM,6);
  else; RunOpt.Coef=1/N;
  end
  if ~isempty(Cam) fc=comstr(Cam,[-1 1]); 
  else; fc=f(end);
  end

  ind=ceil(length(f)*fc/f(end));  if ind>=length(f); ind=length(f) ;end
  ind=1:ind;

  out=frames;i1=size(frames.Y);i1(1)=length(ind);
  if ~iscell(out.X);
    out.X={f(ind) 1:size(frames.Y,2)};
    out.Xlab={'Frequency','channel'};
  else; out.X{1}=f(ind);
  end
  if ~isfield(out,'Xlab')|~isa(out.Xlab,'cell')
   out.Xlab={'Frequency','channel'};
  end
  r2=frames.Y;
  for j1=1:size(r2,3); 
   SPEC=fft(squeeze(frames.Y(:,j1)))*RunOpt.Coef; % DFT or ShockDFT
   out.Y(:,:,j1)=SPEC(ind);
  end
  if isempty(strfind(Cam,'full')); % if not full return the half spectrum
   i1=unique([RunOpt.ib RunOpt.ic]);
   out.X{1}=out.X{1}(i1);  out.Y=out.Y(i1,:);
  end
  
  %---------------------------------------------------------------
  % [H1,H2,Coh]=fe_curve('h1h2',frames,window);
elseif comstr(Cam,'h1h2') [CAM,Cam]=comstr(CAM,5);
  
  frames=varargin{carg};carg=carg+1;
  if carg<=nargin win=varargin{carg};carg=carg+1;
  else win=[];
  end
  if carg<=nargin; pond=varargin{carg};carg=carg+1;
  else; pond=[];
  end
  
  i1=comstr(Cam,[-1 1]); 
  
  i2=1:size(frames{1}.Y,2);i2(i1)=0;i2=find(i2);    
  
  t=frames{1}.X; N=length(t);
  if isempty(win) win='None'; end;
  if ischar(win) win=fe_curve(['window' win],N);end
  win=win(:);
  
  f=1/diff(frames{1}.X(1:2))*[0:length(frames{1}.X)-1]/length(frames{1}.X);
  iw=1:round(450/1024*length(f));
  
  %-- check validity for multiple frames case --%
  for j1=2:length(frames) 
    if length(frames{1}.X)~=length(frames{j1}.X);
      error('Inconsistent frames length');
    end;
  end;
  
  if length(i1)==1  %-- SIMO case
    k1=size(frames{1}.Y);
    k2=1; iout=2:k1(2); % input channel
    
    
    Gyu=zeros(length(iw),k1(2)); Gyy=zeros(length(iw),length(iout));Guu=zeros(length(iw),1);
    SPEC=zeros(size(frames{1}.Y));
    
    %-- Using PSD averaging
    %-- instead of linear spectra
    
    for j2=1:length(frames) % fill in 
      r2=frames{j2}.Y;
      r2 = r2.*win(:,ones(1,size(r2,2)));     % apply time window here
      SPEC=fft(r2);                      % compute FFT
      Gyu=Gyu+SPEC(iw,:).*conj(SPEC(iw,k2*ones(1,k1(2))));
      Gyy=Gyy+SPEC(iw,iout).*conj(SPEC(iw,iout));
      Guu=Guu+SPEC(iw,1).*conj(SPEC(iw,1));
    end
    
    
    H1=Gyu(:,iout)./(Gyu(:,k2*ones(size(iout))));
    H2=conj(Gyy./Gyu(:,iout));
    COH=real(H1./H2);
    out=struct('X',f(iw),'H1',H1,'H2',H2,'COH',COH,'Gyy',Gyy/length(frames)/length(iw),...
      'Guu',Guu/length(frames)/length(iw),'Gyu',Gyu/length(frames)/length(iw));
    
  else  %-- MIMO case
    
    %-- set indices to extract matrix element from vector
    %-- compute only the lower triangular part of cross correlation matrix
    
    ind_comp=[i1 i2];
    ind_1=ind_comp(ones(length(ind_comp),1),:);
    ind_1=ind_1(:);
    ind_2=ind_comp(ones(length(ind_comp),1),:)';
    ind_2=ind_2(:);
    
    if isempty(pond); % default
        pond=norm(frames{1}.Y,inf);
        for j1 = 2:length(frames); pond = pond + norm(j1.Y,inf);end
        pond=diag(pond/length(frames));
    elseif pond==0; pond=ones(size(frames{1}.Y,2),1); 
    else; pond=1./(pond+eps);  %-- to avoid singularities 
    end
    
    if size(pond,2)<2 | size(pond,1)<2 pond=diag(pond(:)); end;
    
    if size(pond,1)<size(frames{1}.Y,2) 
      error('Ponderation error - check for weighing vector');
    end;
    
    
    %-- Using PSD averaging
    %-- instead of linear spectra
    
    Sy=zeros(size(frames{1}.Y));Gaa=zeros(length(iw),length(ind_comp).^2);
    for j1=1:length(frames)
      Y=(frames{j1}.Y*diag(1./diag(pond))).*...
        win(:,ones(size(frames{j1}.Y,2),1));
      Sy=fft(Y,[],1);
      Gaa=Gaa+Sy(iw,ind_1).*conj(Sy(iw,ind_2)); % SM 18/01/2005
    end;
    
    %-- compute estimators
    li1li2=length(i1)*length(i2);
    H1=zeros(length(iw),li1li2);
    Hv=zeros(length(iw),length(i1)*length(i2));
    co=[];
    for j1=1:length(iw)
      gaa=zeros(length(ind_comp));
      gaa=reshape(Gaa(j1,:),length(ind_comp),length(ind_comp)); 
      gyu=gaa(i2,i1);
      guu=gaa(i1,i1);
      gyy=gaa(i2,i2);
      %-- H1 estimator --%
      H1_=gyu*(guu\eye(length(i1)));
      H1(j1,:)=H1_(:)';
      %-- Hv estimator 
      %-- Using pseudo inverse solution for least square problem
       Hv_=([gyu gyy]*[guu gyu']')*...
         (([guu gyu']*[guu gyu']')\eye(length(i1)));
      Hv(j1,:)=Hv_(:)';
      [u,s,v]=svd([guu gyu']*[guu gyu']');
      co(j1,:)=diag(s)';
    end;
    
    if size(pond,2)>=2 & size(pond,1)>=2 pond=diag(pond); end;
    
    Pond(1,:,:)=(pond(:,ones(size(pond,1),1))./pond(:,ones(size(pond,1),1))').^2;
    pond_vec=sqrt(Pond(1,i2,i1));
    pond_vec=pond_vec(:)';
    H1=H1.*pond_vec(ones(length(iw),1),:);
    Hv=Hv.*pond_vec(ones(length(iw),1),:);
    
    out=struct('X',f(iw),'H1',H1,'Hv',Hv,'Gyy',gyy/length(frames)/length(iw),...
      'Guu',guu/length(frames)/length(iw),'Gyu',gyu/length(frames)/length(iw),'Pond',pond_vec,'SingVal',co);
  end;
  
  %---------------------------------------------------------------
  % [Spec]=fe_curve('Spectra ref_channels',frames,nb_pts,window,overlap);
  
elseif comstr(Cam,'spectra') [CAM,Cam]=comstr(CAM,8);
  
  clear frames
  frames=varargin{carg};carg=carg+1;
  nb_pts=varargin{carg};carg=carg+1;
  
  if carg<=nargin win=varargin{carg};carg=carg+1; else win=[]; end;
  if carg<=nargin ovlp=varargin{carg};carg=carg+1; else ovlp=1; end;
  if ovlp>=1 & ovlp<=100 ovlp=ovlp/100;
  elseif ovlp > 100 error('Overlap can not exceed 100%');
  end; 
  if isstruct(frames) frames{1}=frames; end;
  
  i1=comstr(Cam,[-1 1]);
  i1=i1(1); 
  i2=1:size(frames{1}.Y,2);i2(i1)=0;i2=find(i2);
  if isempty(i1) error('No reference channel(s) specified'); end;
  
  
  li1=length(i1);li2=length(i2);li12=size(frames{1}.Y,2);

  Time=frames{1}.X; N=length(Time);
  if isempty(win) win='None'; end;
  if ischar(win) win=fe_curve(['window' win],nb_pts);end;
  win=win(:);
  f=1/diff(frames{1}.X(1:2))*[0:nb_pts-1]/nb_pts;
  iw=1:round(450/1024*length(f));
  win=win(:,ones(li12,1));
  Sy=zeros(li12,nb_pts);
  
  for j1=1:length(frames)
    Gyy{j1}=zeros(length(iw),li12^2);
    
    for k1=1:floor((length(Time)-nb_pts)/((1-ovlp)*nb_pts))+1   
      Y=frames{j1}.Y((1:nb_pts)+round((1-ovlp)*nb_pts*(k1-1)),:).*...
        win;
      Sy=Sy+fft(Y',[],2);
    end;
    
    for l1=1:length(iw)
      Gt=Sy(:,l1)*Sy(:,l1)';
      Gyy{j1}(l1,:)=Gt(:).'./nb_pts./k1;
    end;
  end;
  
  out=struct('X',f(iw),'Y',Gyy,'xlabel',...
    fe_curve('datatype',18),'ylabel',fe_curve('datatype',0),...
    'name','Power spectra');
  %---------------------------------------------------------------
elseif comstr(Cam,'returny')

  if nargin>=carg  curve=varargin{carg};carg=carg+1; end
  if nargin<carg; error('Bad number of input argument'); end
  x=varargin{carg};carg=carg+1;
  if nargin>=carg; model=varargin{carg}; end
  
  if isa(curve,'cell') % array of curves 
    out={};for j1=1:length(curve)
      out{end+1}=fe_curve(Cam,curve{j1},x,model);
    end
    out=reshape([out{:}],length(x),length(out));
    return; 
  end

  if isempty(curve); out=ones(size(x)); return; end

  if ischar(curve); curve=fe_curve('getcurve',model,curve); end

  if ~isfield(curve,'X')|~isfield(curve,'Y')
    curve_name=varargin{carg};carg=carg+1;
    [row,i1]=stack_get(model,'curve',curve_name);
    if isempty(i1); 
     sdtw('_nb','Unknown curve name'); out=ones(size(x));return;
    end 
    curve=row{3};
  end
  
  if isempty(curve.X) curve.X=x; end
  if ischar(curve.Y)
    if comstr(comstr(curve.Y,-27),'step') % Step Tmax
      r1=comstr(curve.Y(5:length(curve.Y)),-1);
      out=zeros(size(x)); out(find(x<r1))=1; 
    else eval(curve.Y);
    end
  else
   if RunOpt.ExtrapBy0
    % extrapolate if necessary 
    i1=find(x>max(curve.X));
    if ~isempty(i1)
     curve.X=[curve.X(:)' x(i1)]; curve.Y=[curve.Y(:)' zeros(1,length(i1))];
     disp('extrapolation')
    end
   end
   eval('out=feutilb(''computer'',curve.Y,curve.X,x,2);');
  end
  
  
  %---------------------------------------------------------------
  % Read
  % val=fe_curve('read h:/sdtdata/gefdyn/cas_test/input.cyt');
  
elseif comstr(Cam,'read'); [CAM,Cam]=comstr(CAM,5);
  
  fin=0;
  if carg<=nargin; fin=varargin{carg};CAM=fopen(fin);Cam=comstr(CAM,-27);end
  if isempty(Cam)
    [fname,wd]=uigetfile(CAM,'Pick a curve containing a file');
    if ~ischar(name) return; end
    [wd,fname,ext]=fileparts(fullfile(wd,fname));
  else
    [wd,fname,ext]=fileparts(CAM);
  end
  if isempty(wd) wd=pwd;end
  FileName=fullfile(wd,[fname ext]);
  
  if comstr(ext,'.cyt') % Cyberquake text format
    
    [fin,mes]=fopen(FileName); if ~isempty(mes) error(mes);end
    
    st=fgetl(fin); if comstr(comstr(st,-27),'#cyber')
      st=fgetl(fin);
      %[i1,st]=comstr(st,'#','%i');
      %st=st(find(st=='=')+1:end);i1(2)=comstr(st,-1);
      r1=fscanf(fin,'%g');
      val=struct('name',fname,'type','fe_curve', ...
        'X',r1(1)+[0:r1(3)-1]'*r1(2),'Y',r1(5:end));
      %[1 2 r1(3) 0 0 0],'data',r1([1:2 5:end]));
      if length(val.Y)==2*r1(3)
        val.X=val.Y(1:2:end); val.Y=val.Y(2:2:end);
      end
    end
    
  elseif comstr(ext,'.cyb') % Cyberquake binary format
    
    [fin,mes]=fopen(FileName,'r','b'); if ~isempty(mes) error(mes);end
    fseek(fin,0,-1);  i1=fread(fin,1,'int8');
    val.name=char(fread(fin,i1,'char')')
    
    i1=fread(fin,2,'int32')
    i1=fread(fin,6,'float32')
    i1=fread(fin,2,'int32')
    
    fclose(fin);
 % Dactron sig format very initial reading - - - - - - - - - - - - -
  elseif comstr(ext,'.sig') 

   if fin<2; [fin,mes]=fopen(FileName,'r'); end
   fseek(fin,0,1);filelen=ftell(fin); fseek(fin,0,-1);
   val=struct('version',{{fread(fin,1,'int32') char(fread(fin,8,'char')')}});
   % ? ? Blocksize ? ?
   val.i1=fread(fin,8,'int32')';
   val.header=reshape(char(fread(fin,30,'char')),10,3)';
   st1=char(fread(fin,70))';
   ind=unique([1 find(st1<26)]);for j1=1:length(ind)-1;
     val.st1{1,j1}=st1(ind(j1)+1:ind(j1+1)-1);
   end
   i2=filelen-val.i1(3)*4-4-ftell(fin);i2=i2/4;
   r1=fread(fin,i2,'int32')';fseek(fin,-i2*4,0);
   r1(2,:)=fread(fin,i2,'float')';
   val.r1=r1;
   
   fseek(fin,filelen-val.i1(3)*4-4,-1);
   val.Y=fread(fin,val.i1(3),'float');
   i1=fread(fin,'int32'); if i1~=0; warning('unexpected end');end
   %figure(1);plot(val.Y);set(gca,'xlim',[1 1e3]);

  else sdtw('''Read%s'' is not a valid call',ext);val=[]; end
  
  if isempty(val)
  elseif carg<=nargin mdl=varargin{carg};
    out=stack_set(mdl,'curve',fname,val);
  else 
    out=val;
  end
  
  % -----------------------------------------------------------------
elseif comstr(Cam,'window'); [N,CAM,Cam]=comstr(CAM,'window','%i');
  
  % Symmetric windows
  if isempty(N)&carg<=nargin N=varargin{carg};carg=carg+1;end
  
  if isempty(Cam)
    out={'None','Exponential','Hanning','Hamming'};
  elseif comstr(Cam,'none')
    out=ones(N,1);
  elseif comstr(Cam,'exponential')
    [CAM,Cam]=comstr(CAM,12); opt=comstr(CAM,[-1 0 10 10]);
    opt(1:2)=round(opt(1:2)); 
    
    N=N-opt(1)-opt(2)-1;  
    if N<0 error('Improper number of zero and flat top points');end
    out=[zeros(opt(1),1);ones(opt(2),1); exp(-[0:N]'/N*opt(3))];
    
  elseif comstr(Cam,'hanning') % Hanning window
    out = (1 - cos(2*pi/(N-1)*(0:N-1)'))/2; 
  elseif comstr(Cam,'hamming') % Hamming window
    out = .54 - .46*cos(2*pi/(N-1)*(0:N-1)');
  end

% -----------------------------------------------------------------
elseif comstr(Cam,'butter'); [opt,CAM,Cam]=comstr(CAM,'butter','%g');
  
if rem(opt(1),2)==1 opt(1)=opt(1)+1; end % only even orders are supported

% opt(2) = freq_filt/(F_ech/2);
%continuous poles
r1 = exp(i*([1:2:opt(1)-1]'/opt(1) + 1)*pi/2);r1=[r1;conj(r1)];
%discrete poles
r1= (1+r1/opt(2))./(1+r1/opt(2));
r1=real(poly(r1));
poly(-ones(1,opt(1)))

  % -----------------------------------------------------------------
elseif comstr(Cam,'cvs')
  
  out='$Revision: 1.69 $  $Date: 2006/05/17 09:50:56 $';
  
%---------------------------------------------------------------
%-- Perfect band pass filter --%  
elseif comstr(Cam,'bandpass'); [CAM,Cam]=comstr(CAM,10);
  
  if comstr(Cam,'hz') [CAM,Cam]=comstr(CAM,3); f_band=str2num(Cam);
  elseif comstr(Cam,'rd') [CAM,Cam]=comstr(CAM,3); f_band=str2num(Cam)/2/pi;
  else sdtw('_nb','Assuming fmin and fmax are in Hz !!');f_band=str2num(Cam);
  end;
  
  
  if ~isa(varargin{carg},'cell') frames={varargin{carg}};carg=carg+1;
  else frames=varargin{carg};
  end;    
  
  for i1=1:length(frames)
    
    [t,f,N]=fe_curve('getxtime',frames{i1});
    RunOpt=fe_curve('getia',N);  
    out{i1}.X=t;  
    
    if find(size(frames{i1}.Y)==N)==1 Y=fft(frames{i1}.Y);
    else Y=fft((frames{i1}.Y).');
    end;
    Y(find( f < min(f_band) | f > max(f_band) ),:)=0;
    Y(RunOpt.ia,:)=conj(Y(RunOpt.ib,:)); 
    Y(RunOpt.ic,:)=real(Y(RunOpt.ic,:));
    out{i1}.Y=real(ifft(Y)); 
  end
  
  if i1==1 out=out{1}; end;

%---------------------------------------------------------------
%-- Zoom FFT 
elseif comstr(Cam,'zoomfft'); [CAM,Cam]=comstr(CAM,8);

  if ~isa(varargin{carg},'cell') frames={varargin{carg}};carg=carg+1;
  else frames=varargin{carg};
  end;    

  ftarg=comstr(CAM,-1); if length(ftarg)>2; 
   error('You must provide center frequency and zoom factor');
  end

  for i1=1:length(frames)
    

    t=frames{i1}.X(:);
    r1=frames{i1}.Y; if find(size(r1)==length(t))==2;r1=r1';end
    r1=fft(r1.*exp(-sqrt(-1)*2*pi*ftarg(1)*t)); % shift and transform

     N=length(t)/ftarg(2); r1(N:size(r1,1)-N)=0;  % filter  
     r1(end:-1:end-N+1)=conj(r1(2:N+1));                  % make symmetric
     r1(1)=real(r1(1));
 
     x=real(ifft(r1));x=x(1:ftarg(2):end);tx=t(1:ftarg(2):end);
     fx=[1/diff(tx(1:2))*[0:length(tx)-1]'/length(tx)]+ftarg(1);
     r1=fft(x)*ftarg(2);%/length(fx);
 
     ind=round(length(fx)*[.05 .45]);ind=ind(1):ind(2);  
     out=struct('X',fx(ind),'Y',r1(ind,:), ...
        'xunit',{fe_curve('datatypecell','Freq')});
    
  end
 
  
%---------------------------------------------------------------
%-- Unitary & given PSD spectrum noise generation from randomly 
%-- dephased cosines addition
  
elseif comstr(Cam,'noise') [CAM,Cam]=comstr(CAM,6);
  
  M=varargin{carg};carg=carg+1;
  
  if mod(M,2)==0 M=M/2; is_even=1;
  else M=(M+1)/2; is_even=0;
  end;
  
  fe=varargin{carg};carg=carg+1;
  f_max=[];
  Spec=[];
  
  if nargin==4
    if length(varargin{4})==1
      f_max=varargin{carg};carg=carg+1;
      if f_max>fe/2
        error('Select a max frequency <= sampling rate / 2');
      end;
    elseif comstr(class(varargin{4}),'double')
      Spec=varargin{4};
      if length(Spec)~=M
        warning('resampling PSD spectrum...');
        if length(Spec)<M
          Spec=interp1((0:length(Spec)-1)/(length(Spec)-1),...
            Spec,(0:M-1)/(M-1));
        else
          M_ech=lcm(M,length(Spec));
          Spec=interp1(round(linspace(1,M_ech,length(Spec))),...
            Spec,1:M_ech);
          Spec=Spec(round(linspace(1,M_ech,M)));
        end;
      end;
    else
      error('Bad argument type');
    end;
  end;
  
  if isempty(f_max) & isempty(Spec)
    Ck=ones(M,1).*exp(j*rand(M,1)*2*pi);
  elseif ~isempty(f_max) & isempty(Spec)
    Ck=[ones(ceil(M*(2*f_max/fe)),1) ; ...
        zeros(M-length(ones(ceil(M*(2*f_max/fe)),1)),1)].*...
      exp(j*rand(M,1)*2*pi);
  elseif isempty(f_max) & ~isempty(Spec)
    Ck=sqrt(Spec(:)).*exp(j*rand(M,1)*2*pi);
  end;
  
  
  Ak=[Ck ; zeros(M,1)];
  
  out.X=0:1/fe:2*M/fe-1/fe;
  out.Y=2*(fe)^(1/2)*real(ifft(Ak));
  
  if is_even==0
    out.Y=out.Y(1:end-1)-mean(out.Y(1:end-1));
    out.X=out.X(1:end-1);
  else 
    out.Y=out.Y-mean(out.Y);
  end;
  
  %---------------------------------------------------------------
  %---------------------------------------------------------------
  
elseif comstr(Cam,'test'); [CAM,Cam]=comstr(CAM,5);

% Find ID if given
ID=[];i1=strfind(Cam,'-id'); if ~isempty(i1);
 [ID,i2,i3,i4]=sscanf(CAM(i1+3:end),'%i');
 CAM(i1:i1+1+i4)=''; [CAM,Cam]=comstr(CAM,1);
end
if isempty(ID); ID=1;end

% 'TestFrame'
if comstr(Cam,'frame');
  
  
  fs=512;  %-- sampling frequency
  ech_length=4;  %-- sample length (s)
  
  noise=fe_curve('Noise',fs*ech_length,fs);  %-- computes noise
  
  %-- build the curve associated to the time signal of noise
  
  out{1}=struct('X',noise.X,'Y',noise.Y,'xunit',...
    fe_curve('DataType','Time'),'yunit',...
    fe_curve('DataType','Excit. force'),'name','Input at DOF 2');
  
  
  
  %-- set up an oscillator with 3 DOF --%
  
  Puls = [30 80 150]'*2*pi;  %-- natural frequencies
  Damp = [.02 .015 .01]';  %-- damping
  Amp = [1 2 -1;2 -1 1;-1 1 2];  %-- pseudo "mode shapes"
  Amp=Amp./det(Amp);
  
  %-- set up command and observation matrices --%
  
  C=[1 0 0];  %-- Observation matrix
  B=[0 1 0]';  %-- Command matrix
  
  %-- Compute frequencies vector --%
  
  freq=([0:length(noise.X)-1]/length(noise.X))*fs*2*pi;
  
  %-- Eliminating frequencies corresponding to the aliased part of the noise
  %-- spectrum
  
  freq=freq(1:length(noise.X)/2);
  
  %-- Compute tranfert function 
  
  FRF=nor2xf(Puls,Damp,Amp*B,C*Amp,freq);
  
  %-- Compute the time response to input noise
  
  Resp=fe_curve('TimeFreq',noise,[FRF ; zeros(length(FRF),1)].');
  
  %-- build the curve associated to the time signal of response
  
  out{2}=struct('X',Resp.X,'Y',Resp.Y,'xunit',fe_curve('DataType','Time'),...
    'yunit',fe_curve('DataType','Displacement'),'name','Output at DOF 1');
  
  if nargout==0
    h1=figure(1);
    fe_curve('Plot',h1,out{1});    
    figure(2);
    semilogy(freq/2/pi,abs(FRF));
    title('FRF DOF 2 / DOF 1')
    xlabel('Frequency Hz');
    ylabel('Transfert m/N');
    h2=figure(3);
    fe_curve('Plot',h2,out{2});
  end;


%TestRamp NStep FinalValue
elseif comstr(Cam,'ramp');[CAM,Cam]=comstr(CAM,5);

 opt=comstr(CAM,[-1 5 1]);
 r1=1:opt(1);r1=r1(:);
 out=struct('X',r1,'Y',r1/r1(end)*opt(2),'name','ramp');

%TestRicker SignalDuration Nstep Amplitude TotalTime
elseif comstr(Cam,'ricker');[CAM,Cam]=comstr(CAM,7);

 opt=comstr(CAM,[-1 5 1]);

 t=linspace(0.,opt(4),opt(2));
 ts=opt(1)/2; tp=opt(1)/2;
 a = (t-ts).*(t-ts)/(tp*tp);
 y = opt(3)*(1 - exp(1)*(a .* exp(-a)));
 i1=find(t>2*ts); y(i1)=0.;

 out=struct('ID',ID,'X',t,'Y',y,'name','ricker');

%'Test[Sin, Cos, ...] SignalDuration Nstep Amplitude TotalTime
%                     standard curve generation 
else % if ~isempty(strmatch(Cam,{'sin','cos','tan','exp'}));

 if nargin>1
  r1=varargin{carg};
  if min(size(r1)==1); r1=r1(:);end
  out=struct('X',r1,'Y',feval(Cam,r1),'ID',ID, ...
      'name',sprintf('%s(X)',Cam));
 else
  f1=Cam(1:3);
  [CAM,Cam]=comstr(CAM,4);
  opt=comstr(CAM,[-1 5 1]);
  t=linspace(0.,opt(4),opt(2));
  if any(exist(f1)==[2 3 5 6])
    ft=zeros(length(t));   i1=find(t<=opt(1));
    ft(i1)=feval(f1,pi*t(i1)/opt(1));
    out=struct('X',t,'Y',ft,'ID',ID, ...
     'name',sprintf('%s(X)',f1));
  else sdtw('_nb','''Test%s'' is not a known command',CAM);    
  end
 end % nargin
end;
if nargout==0; fe_curve('plot',out);end

  
%---------------------------------------------------------------
% This validates the structure of a given curve structure and makes any
% appropriate modification
elseif comstr(Cam,'fixcurve')
  
  curve=varargin{carg}; carg=carg+1;
  
  st=fieldnames(curve);
  
  % This curve is a SDT <5.1 frequency/time response curve
  if length(intersect(st,{'w','xf'}))==2
    
    out = struct('X',curve.w,'Y',curve.xf,'xunit',[],'yunit',[], ...
      'zunit',[],'name',[],'type','fe_curve','unit','USER', ...
      'Interp','','PlotFcn',[]);
    st=setdiff(st,{'w','xf','x','yd','yn'});
    for j1=1:length(st); out=setfield(out,st{j1},getfield(curve,st{j1}));end
  else;  out=curve;
  end
  if size(out.X,1)==1&size(out.X,2)>1;out.X=out.X(:);end
  if size(out.Y,1)==1&size(out.Y,2)>1;out.Y=out.Y(:);end

  if ~isfield(out,'name'); out.name=''; end
  if isfield(curve,'name'); out.name=curve.name;end
  if isfield(curve,'x'); out.xunit=curve.x; end
  try;if isfield(curve,'yn')&isfield(curve,'yd'); out.yunit={};
   for j1=1:size(curve.yn,1)
    r2=curve.yn{j1,2}-curve.yd{j1,2};r2(1)=0;
    out.yunit(j1,1:3)={sprintf('%s/%s',curve.yn{j1,1},curve.yd{j1,1}), ...
     sprintf('%s/%s',curve.yn{j1,2},curve.yd{j1,2}),r2};
  end;end;end

  curve=fe_curve('datatype fix',curve);

  %---------------------------------------------------------------
else sdtw('_nb','''%s'' is not a known command',CAM);    
end;

%---------------------------------------------------------------
%---------------------------------------------------------------

%---------------------------------------------------------------------
%---------------------------------------------------------------------
%---------------------------------------------------------------------
function r2 = DataType(r1,def,unit);

if nargin<2|isempty(def); def={0,'None'}; end
if nargin>0&ischar(r1)&~isempty(strfind(r1,'-cell'));
 OutTyp='cell';r1(strfind(r1,'-cell')+[0:4])='';r1=comstr(r1,1);
else;OutTyp='struct';
end

% Typ [Length Force Temp Time exponents]
DataTypes={[0 0 0 0 0],'Unknown','NONE';
  [1 0 0 0 0],'General','NONE';
  [2 -2 1 0 0],'Stress','Pa';
  [3  0 0 0 0],'Strain','';
  [5 0 0 1 0],'Temperature','^o';
  [6 1 1 0 0],'Heat flux','NONE';
  [8  1 0 0 0],'Displacement','m';
  [9  0 1 0 0],'React. force','N';
  [11 1 0 0 -1],'Velocity','m/s';
  [12 1 0 0 -2],'Acceleration','ms^{-2}';
  [13 0 1 0 0],'Excit. force','N';
  [15 -2 1 0 0],'Pressure','Pa';
  [16 -1 1 0 0],'Mass','Kg';
  [17 0 0 0 1],'Time','s';
  [18 0 0 0 -1],'Frequency','Hz';
  [19 0 0 0 -1],'RPM','Hz'; 
  [20 0 0 0 0],'Order','?';    % ABOVE HERE UFF SPECIFICATION
  [100 -1 1 0 0],'Moment','N'; % BELOW HERE OPENFEM Specific
  [101 0 0 0 0],'Rotation','rad';
  };

i1=reshape([DataTypes{:,1}],5,size(DataTypes,1))';
if nargin==0
  r2=struct('label',[],'value',i1);r2.label=DataTypes(:,2);
  fprintf('Available data types\n\n')
  out=[num2cell(i1(:,1)) r2.label(:)]';
  fprintf('%4i : %s\n',out{:});
  return;
elseif isnumeric(r1) % match against number
  i2=find(i1(:,1)==r1(1)); 
  if ~isempty(i2) r1=DataTypes(i2,[2 3 1]); else; r1=def;end
elseif isstruct(r1)
elseif iscell(r1)&size(r1,2)==3  % already cell format do nothing
elseif comstr(comstr(r1,-27),'_datatype')
  r2=struct('label',[],'value',i1);r2.label=DataTypes(:,2); return; 
elseif ischar(r1)
  i2=strmatch(comstr(r1,-27),comstr(DataTypes(:,2),-27));
  if isempty(i2) sdtw('_nb',sprintf('''%s'' not matched',r1));r1=def;
  else; 
   r1=DataTypes(i2,[2 3 1]);
   ;r1=struct('label',r1{1,1},'unit',r1{1,2},'type',r1{1,3});
  end
  
else; r1=def;
end
if isa(r1,'cell') % Newer cell array format
  if size(r1,2)<3; r1(:,3)={[0 0 0 0 0]};end
  for j1=1:size(r1,1);
   if isempty(r1{j1,3})|(r1{j1,3}(1)==0 &  def{1}); r1{j1,3}(1)=def{1}; end
   if  r1{j1,3}(1)~=0; r1{j1,3}(2:5)=i1(find(i1(:,1)==r1{j1,3}(1)),2:5); end
   % MISSING label checks
  end
elseif isstruct(r1)
 % checks on type
 if r1.type(1)==0 &  def{1}; r1.type(1)=def{1}; end
 if ~any(r1.type(1)==i1(:,1)); r1.type(1)=0; end
 if r1.type(1)~=0;  r1.type(2:5)=i1(find(i1(:,1)==r1.type(1)),2:5); end

 if strcmpi(r1.label,'none')|strcmpi(r1.label,'unknown')|isempty(r1.label)
   r1.label=DataTypes{find(i1(:,1)==def{1}),2};
 end
 if strcmpi(r1.label,'none')|isempty(r1.label) 
   r1.label=DataTypes{find(i1(:,1)==r1.type(1)),2};
 end
 if strcmpi(r1.unit,'none')|isempty(r1.unit) r1.unit=def{2}; end
end

if strcmp(OutTyp,'cell')&isstruct(r1);
 r2={r1.label,r1.unit,r1.type};
else; r2=r1;
end




%---------------------------------------------------------------
%---------------------------------------------------------------

function [r1,list,cur]=get_curve(GEF,cur);

[r1,ind]=stack_get(GEF,'curve');
list=r1(:,2);r1=r1(:,3);

if isempty(r1) i1=0; list={'New ...'};r1=cell(1,1);
else list{end+1,1}='New ...'; r1{end+1,1}=[];end

if nargin==1 cur=[];elseif cur<=length(ind) cur=ind(cur); else cur=[];end
