function varargout = ftraj(z,x,t,flag)
  % index 2 - trajectory control
  % Parameters:
 if isempty(flag)
     [varargout{1},varargout{2}]= f(z,x,t);
  else
    switch flag
     case ''
      [varargout{1},varargout{2}] =  f(z,x,t);
     case 'dfy'
      varargout{1} = dfy(z,x,t);
     case 'dfx'
      varargout{1} = dfx(z,x,t);
     case 'dft'
      varargout{1} = dft(z,x,t);
     case 'd'
      varargout{1} = d(x,t);
     case 'dx'
      varargout{1} = dx(x,t);
     otherwise
      error(['Unknown flag ''' flag '''.']);
    end
  end

  function [fxt,ok] = f(z,x,t)
  ok=1;
  ae=0.209029e+8;
  mue=0.1407653916e+17;
  ome= 0.72921159e-4;
  m=0.2890532728e+1;
  s=1;
  fxt = zeros(size(x,1),1);
  %   if(abs(x(4)).lt.uround) return 1 
  r=x(1)+ae;
  %    if(abs(r).lt.uround) return 1
  g=mue/r/r;
  hpi=pi/18d1;
  hsing=sin(x(5)*hpi);
  hsina=sin(x(6)*hpi);
  hcosg=cos(x(5)*hpi);
  hcosa=cos(x(6)*hpi);
  hsinl=sin(x(3)*hpi);
  hcosl=cos(x(3)*hpi);
  hn1=r*hcosl;
  %     if(abs(hn1).lt.uround) return 1
  hn2=x(4)*hcosg;
  %     if(abs(hn2).lt.uround) return 1       
  vr2=x(4)*x(4);
  hcl=cl(x(7));
  % check the value of x(1)!
  rho=0.002378e0*exp(-x(1)/238e2);
  l=0.5*rho*hcl*s*vr2;
  cd=0.04e0+0.1*hcl*hcl;
  d=0.5*rho*cd*s*vr2;
  ome2=ome*ome;
  sc=hsinl*hcosa;
  fxt(1)=z(1)-x(4)*hsing;
  fxt(2)=z(2)-x(4)*hcosg*hsina/hn1;
  fxt(3)=z(3)-x(4)*hcosg*hcosa/r;
  fxt(4)=z(4)+(d/m+g*hsing+ome2*r*hcosl*(sc*hcosg-hcosl*hsing));
  fxt(5)=z(5)+(-l*cos(x(8)*hpi)/m/x(4)-hcosg/x(4)*(vr2/r-g)-2d0*...
      ome*hcosl*hsina-ome2*r*hcosl/x(4)*(sc*hsing+hcosl*hcosg));
  fxt(6)=z(6)+(-l*sin(x(8)*hpi)/m/hn2-x(4)/r*hcosg*hsina*tan(x(3)*...
      hpi)+2*ome*(hcosl*hcosa*tan(x(5)*hpi)-hsinl)-ome2*r*hcosl*...
      hsinl*hsina/x(4)/hcosg);
  fxt(7)=x(5)+1+9*t*t/9e4;
  fxt(8)=x(6)-45-90*t*t/9e4 ; 
 if (sum(isnan(fxt)) > 0 | sum(isinf(fxt)) > 0)
      ok=0;  
      error('funktion f_traj not calculable');
  end;

function [dfyzxt]=dfy(z,x,t,varargin);
dfyzxt=[speye(size(z,1));sparse(size(x,1)-size(z,1),size(z,1))];

function [dfxzxt]=dfx(z,x,t,varargin);
dfyzxt=sparse(size(x,1),size(x,1));
error('not yet realized');
dfyzxt(1,4)=-hsing;

function [dftzxt]=dft(z,x,t,varargin);
dftzxt=[sparse(6,1);2*t/1e4;-2*t/1e3];


function [dxt]=d(x,t);
dxt=x(1:6);

function [dxxt]=dx(x,t);
dxxt=[speye(6),sparse(6,2)];

function clp=cl(p)
clp=1e-2*p;









