// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2010 - INRIA - Serge Steer
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function champ3d1(x,y,z,fx,fy,fz,arfact)
//  x,  y,  z:  vectors ( x(i), y(i), z(i))  define the coordinates of the ith vector base
// fx, fy, fz:  vectors (fx(i),fy(i),fz(i)) define the  components of the ith vector
  if argn(2)==0 then //demonstration
    //SIR differential equation
    ns=10;ni=10;nr=10;
    S=linspace(0,1,ns);
    I=linspace(0,1,ni);
    R=linspace(0,1,nr);

    a=1;b=0.3;
    
    dS=matrix(-a*S.*.I.*.ones(R),ns,ni,nr);
    dI=matrix((a*S-b).*.I.*.ones(R),ns,ni,nr);;
    dR=matrix(b*ones(S).*.I.*.ones(R),ns,ni,nr);;
    f=gcf();f.color_map=hotcolormap(128);
    clf;champ3d1(S,I,R,dS,dI,dR,0.2)
    a=gca()
    a.x_label.text='S';a.y_label.text='I';a.z_label.text='R';
    function ydot=rhs(t,y)
       a=1;b=0.3;
       ydot=[-a*y(1)*y(2);(a*y(1)-b)*y(2);b*y(2)]
    endfunction
    for s=[0.6 0.8 0.9 0.999]
      y=ode([s;1-s;0],0,0:1000,rhs)
      xpoly(0,0);e=gce();e.data=y';e.thickness=3;
    end
    return
  end
  
  if exists('arfact','local')==0 then arfact=1,end
  if type(x)<>1|~isreal(x) then
    error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"champ3d1",1))
  end
  if type(y)<>1|~isreal(y) then
    error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"champ3d1",2))
  end
  if type(z)<>1|~isreal(z) then
    error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"champ3d1",3))
  end
  if typeof(fx)<>"hypermat"|~isreal(fx)|ndims(fx)<>3 then
    error(msprintf(_("%s: Wrong type for input argument #%d: Real 3D hypermatrix expected.\n"),"champ3d1",4))
  end
  if typeof(fy)<>"hypermat"|~isreal(fy)|ndims(fy)<>3 then
    error(msprintf(_("%s: Wrong type for input argument #%d: Real 3D hypermatrix expected.\n"),"champ3d1",5))
  end
  if typeof(fz)<>"hypermat"|~isreal(fz)|ndims(fz)<>3 then
    error(msprintf(_("%s: Wrong type for input argument #%d: Real 3D hypermatrix expected.\n"),"champ3d1",6))
  end
  nx=size(x,'*');
  ny=size(y,'*');
  nz=size(z,'*');
  
  if or(size(fx)<>[nx ny nz]) then
    error(msprintf(_("%s: Incompatible dimensions for input argument #%d.\n"),"champ3d1",4))
  end
  if or(size(fy)<>[nx ny nz]) then
    error(msprintf(_("%s: Incompatible length for input arguments #%d and #%d.\n"),"champ3d1",5))
  end
  if  or(size(fz)<>[nx ny nz]) then
    error(msprintf(_("%s: Incompatible length for input arguments #%d and #%d.\n"),"champ3d1",6))
  end


  x = matrix(x ,1,-1);  y=matrix( y,1,-1); z=matrix( z,1,-1);
  n1=sqrt((max(x)-min(x))^2+(max(y)-min(y))^2+(max(z)-min(z))^2)/30;
  fx= matrix(fx,1,-1); fy=matrix(fy,1,-1);fz=matrix(fz,1,-1);
  n=sqrt(fx.^2+fy.^2+fz.^2)/n1;
  n(n==0)=%eps;
  
  //create the mesh 
  X=x.*.ones(y).*.ones(z)
  Y=ones(x).*.y.*.ones(z)
  Z=ones(x).*.ones(y).*.z
  
  f=gcf();dm=f.immediate_drawing;f.immediate_drawing='off'; //delay drawings
  a=gca();
  a.view='3d';
  if a.children==[] then
    a.data_bounds=[ min(X),min(Y),min(Z);
                    max(X),max(Y),max(Z)];
    a.axes_visible='on';

  else
    a.data_bounds=[min(a.data_bounds(1,:), [min(X),min(Y),min(Z)]);
                   max(a.data_bounds(2,:), [max(X),max(Y),max(Z)])];
  end
  xpoly(X,Y,'marks');
  e=gce();e.data(:,3)=Z(:);
  X=[X;X+fx./n]; 
  Y=[Y;Y+fy./n]; 
  Z=[Z;Z+fz./n];

  //builds the segments
  xsegs([0,0],[0,0]);e=gce();e.arrow_size=arfact;

  e.data=[X(:),Y(:),Z(:)];
  nc=size(f.color_map,1);
  c=n-min(n);c=1+((nc-1)*(c./max(c)));
  e.segs_color=c;

  f.immediate_drawing=dm;
endfunction