yorick banner

Home

Manual

Packages

Global Index

Keywords

Quick Reference


/*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
gadget.i defines
gadget structure
CHalo structure
is_real(x) 
read_halo_tree(fname) 
grid_read(fname,big_endian=) 
gadget_read(fname) 
_gadget_read(file, &address, what, ..) 
gadget_write(s,fname) 
_gadget_write(file,what,&address) 
gadgetcenter(q,pos,inplace=) 
showxyz(x,y,z,txt=,nsqr=,r200=,center=,spin=,axis=,noerase=,color=) 
showsim(fname,imin=,imax=,pers=,vel=,dens=,delay=,jpg=,isovel=,
xrange=,yrange=,zrange=,sub=,npix=,col=,iso=,dir=,center=,notxt=,cut=,nodisp=,
marker=,comp=,rot=,flag=,color=,noerase=)
checksim(dir) 
readSnapTimings(dummy,dir=) 
readTimings(fname,nf=) 
showhaloenv(q,center,cut=,range=) 
showgal(q,rot=,comp=) 
gadget_merge(q1,q2,rot=,trans=,vtrans=,scale=,flip=) 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/


func is_real (x) { s= structof(x); return s==double || s==float; }

struct GadgetSnapshot {
  /* Header part (note: total header size is 256 bytes). */
  int npart(6);
  double massarr(6);
  double time;
  double redshift;
  int flag_sfr;
  int flag_feedback;
  int nparttotal(6);
  int flag_cool
  int nfiles
  double boxsize
  double omegam
  double omegav
  double hubble
  /* Data part (note: data part starts at offset 256). */
  pointer pos, vel, id, mass, u, rho;
}



struct CHalo {
  long   level  ;
  long   mother;
  long   firstchild;
  long   nsisters;
  long   sister;
  double rho_saddle;
  double density;
  double densmax;
  double radius ;
  double mass   ;
  double position(3);
// level : integer number corresponding to the level in the tree.
//         The root, node(1), corresponds to level=0. The root corresponds
//         to all the selected particles (above the density threshold)
//         in the simulation box: it is not physically relevant.
// mother : the id of the mother of this node, i.e. the id of substructure
//         which contains this substructure (for the root, there is no
//         mother, i.e. mother=0).
// firstchild : the id of the first child of this node, i.e. the first
//         substructure it contains. Other substructures will be found
//         by running over the list their sisters. If we are sitting
//         at the end of a branch, firstchild=0 (i.e. this node correspond
//         to an elementary substructure) .
// nsisters : the number of sisters of this substructure + 1.
// sister : the id of the sister ``at the right'' of this structure. A structure
//         that has not sister at its ``right'' has sister=0.
// rho_saddle : the density of the saddle(s) point(s) connecting this
//         structure to its sisters. In principle this could be calculated
//         uniquely so that a structure would always have at most one
//         sister, but as explained in the important note above,
//         rho_saddle is computed approximately, as a result, it can
//         correspond to more than a single saddle point.
//         For a given density threshold, rho_th, if rho_th > rho_saddle,
//         the substructure is connected to its sisters, whereas if
//         rho_th < rho_saddle, the substructure is disconnected from its
//         sisters.
// density : average density of the substructure, calculated from the
//         average SPH density of all particles with SPH density larger than
//         rho_saddle and belonging to the substructure.
// densmax : maximum density within the substructure : it correspond
//         to the largest SPH density of all the particles within the
//         substructure.
// radius : size in Mpc of the substructure calculated by computing the mean
//         square distance of each particle belonging to it with respect
//         to its center of mass.
// mass :  number of particles within the substructure.
// position(3) : coordinates in Mpc of the center of mass of the
//         substructure.
}


func read_halo_tree (fname,cut=)
/* DOCUMENT 
     read_halo_tree(fname) reads a 
     anahop file and returns a CHalo structure.
     EXAMPLE
     h1=read_halo_tree("seed064853998_128_012_000.nodes");
     ws1; PL,h1.density,h1.mass;
   SEE ALSO:
 */
{

 ff=open(fname);

  vutils=array(long,2);
  read,ff,vutils;

  nnode=vutils(1);

  Node=array(double,15,nnode);
  read,ff,Node;
  close,ff;

   hh=array(CHalo,nnode); 
   hh.level=Node(2,);
   hh.mother=Node(3,);
   hh.firstchild=Node(4,);
   hh.nsisters=Node(5,);
   hh.sister=Node(6,);
   hh.rho_saddle=Node(7,);
   hh.density=Node(8,);
   hh.densmax=Node(9,);
   hh.radius=Node(10,);
   hh.mass=Node(11,);
   hh.position=Node(13:15,);
  return hh;

  //    ff=open(fname);
//    a=0; nn=0; 
//    read,ff,nn,a;
//    if(!is_void(cut)) nn =min(cut,nn);
//    res=[];
//    for(i=1;i<=nn;i++)  {
//      hh=CHalo(); 
//      a=array(0,[1,6]);read,ff, a;
//      hh.level=a(2);
//      hh.mother=a(3);
//      hh.firstchild=a(4);
//      hh.nsisters=a(5);
//      hh.sister=a(6);
//      a=array(0.,[1,9]);read,ff, a;
//      hh.rho_saddle=a(1);
//      hh.density=a(2);
//      hh.densmax=a(3);
//      hh.radius=a(4);
//      hh.mass=a(5)   ;
//      hh.position=a(7:9);
//      grow,res,hh;
//    }
//   close,ff; 
//   return res;

}



/*
  address   size description
  0         4    dummy
  4         256  header
  260       4    dummy
  264       4    dummy
  268       ?    data   pos
  ?         4    dummy
  ?         4    dummy
  ?         ?    data   vel
  ?         4    dummy
  ?         4    dummy
  ?         ?    data   id
  ?         4    dummy
  if (ntot_withmasses>0) skip;
  if (ntot_withmasses>0) skip;
*/









func grid_read (fname,big_endian=)
/* DOCUMENT    grid_read(fname)  return a grid structure
 */      
{
  local stream, address;
  
  /* Open file. */
  if (is_void(fname)) error," grid_read needs an argument";
  if (!is_void(big_endian)) encoding = "sun";
  stream = _grid_open(fname, encoding=encoding);
  address = 0;

  /* Read the header and the data. */
  header = _grid_read(int, 4);
  ngridx = header(2);
  ngridy = header(3);
  ngridz = header(4);
  write, format="reading a %d%d%d array\n", ngridx, ngridy, ngridz;
  tab = _grid_read(float, [3,ngridx,ngridy,ngridz]);
  dummy = _grid_read(int);
  close, stream;
  return tab;
}

func grid_write (tab,fname,big_endian=)
/* DOCUMENT    grid_write(tab,fname)
*/
{
  local stream, address;

  /* Open file. */
  if (is_void(fname)) error," grid_write needs an argument";
  if (!is_void(big_endian)) encoding = "sun";
  stream = _grid_new(fname, encoding=encoding, overwrite=1);
  address = 0;
  
  /* Write the header and the data. */
  tab = float(tab);
  dd = dimsof(tab); 
  header = array(int, 4);
  header(2) = dd(2);
  header(3) = dd(3);
  header(4) = dd(4);
  _grid_write, header;
  _grid_write, float(tab);
  _grid_write, int(0);
  close, stream;
}


func _grid_open (filename, encoding=)
{
  stream = open(filename, "rb");
  if (is_void(encoding)) encoding= "native";
  if (encoding != "native") {
    set_encoding = symbol_def(encoding+"_primitives");
    if (is_func(set_encoding) != 1) error, "bad encoding \""+encoding+"\"";
    set_encoding, stream;
  }
  save, stream, complex; /* make stream aware of the definition of a complex */
  return stream;
}

func _grid_new (filename, overwrite=, encoding=)
{
  if (! overwrite && open(filename, "r", 1))
    error, "file \""+filename+"\" already exists";
  logfile = filename + "L";
  if (open(logfile, "r", 1)) logfile = string(0);
  stream = open(filename, "wb");
  if (logfile) remove, logfile;
  if (is_void(encoding)) encoding= "native";
  if (encoding != "native") {
    set_encoding = symbol_def(encoding+"_primitives");
    if (is_func(set_encoding) != 1) error, "bad encoding \""+encoding+"\"";
    set_encoding, stream;
  }
  save, stream, complex; /* make stream aware of the definition of a complex */
  return stream;
}

func _grid_read (type, dimlist)
{
  extern stream, address;
  data = array(type, dimlist);
  _read, stream, address, data;
  address += sizeof(data);
  return data;
}

func _grid_write (..)
{
  extern stream, address;
  local data;
  while (more_args()) {
    eq_nocopy, data, next_arg();
    _write, stream, address, data;
    address += sizeof(data);
  }
}




func gadget_read (fname,big_endian=)
/* DOCUMENT    gadget_read(fname)  return a gadget structure
    q= gadget_read("./snapshot_000")
    for instance  *q.pos holds the  positions;
    modified to be gadget 1.1 compatible
*/      
{
  /* Open file. */
  if (is_void(fname)) error," gadget_read needs an argument";
  file= open(fname, "rb");
 if (!is_void(big_endian)) sun_primitives,file; 

  s= GadgetSnapshot();

  /* Read the header. */
  address=4 ;
  s.npart=         _gadget_read(file,address,int,6);
  s.massarr=       _gadget_read(file,address,double,6);
  s.time=          _gadget_read(file,address,double);
  s.redshift=      _gadget_read(file,address,double);
  s.flag_sfr=      _gadget_read(file,address,int);
  s.flag_feedback= _gadget_read(file,address,int);
  s.nparttotal=    _gadget_read(file,address,int,6);
  s.flag_cool= _gadget_read(file,address,int);
  s.nfiles= _gadget_read(file,address,int);
  s.boxsize=_gadget_read(file,address,double);
  s.omegam=_gadget_read(file,address,double);
  s.omegav=_gadget_read(file,address,double);
  s.hubble=_gadget_read(file,address,double);
  //bytesleft= 256 - (6*4 + 6*8 + 8 + 8 + 4 + 4 + 6*4 + 4 + 4 + +8 +8 +8 +8);
  //la=            _gadget_read(file,address,char,bytesleft);
  address= 260;
  

  /* Read the data part. */
  ind= where((s.npart > 0) * (s.massarr == 0)) ;
  if (is_array(ind)) {
    nwithmass= sum(s.npart(ind));
  } else {
    nwithmass= 0;
  }
  n= sum(s.npart);  address+= 8;
  s.pos= &_gadget_read(file,address,float,3,n);  address+= 8;
  s.vel= &_gadget_read(file,address,float,3,n);  address+= 8;
  s.id=  &_gadget_read(file,address,int,n);
  if (nwithmass > 0) {    address+= 8;
    s.mass= &_gadget_read(file,address,float,nwithmass);
  }
  
  ngas= s.npart(1);
  if (ngas > 0) {    address+= 8;
    s.u=   &_gadget_read(file,address,float,ngas);    address+= 8;
    s.rho= &_gadget_read(file,address,float,ngas);
  }
  //  info,s;  write,s.npart;
  return s;
}




func hydra_read (fname,big_endian=)
/* DOCUMENT    hydra_read(fname)  return a hydra slice
    q= hydra_read("./snapshot_000")
*/      
{
  /* Open file. */
  if (is_void(fname)) error," hydra_read needs an argument";
  file= open(fname, "rb");
 if (!is_void(big_endian)) sun_primitives,file; 

  address= 0;
  n=long(2028);
  slice= _gadget_read(file,address,float,n,n); 
  return slice;
}


func _gadget_read (file, &address, what, ..)
/* DOCUMENT _gadget_read(file, address, expr)
       -or- _gadget_read, file, address, variable;

     Read unformated, binary data from FILE at offset ADDRESS (must be a
     scalar long integer).  On return, ADDRESS is incremented by the size of
     the data read.  The third argument gives the data type and dimension list
     of the array to read.  If called as a function, the argument EXPR may
     either be an expression or a variable which is filled with the data read
     and is returned as the result of the call.  If called as a subroutine,
     the third argument should be a predefined variable.

  SEE ALSO gadget_read,_gadget_write, open, _read. */
{
  if (! is_array(what)) {
    if (typeof(what)=="struct_definition") {
      type= what;
      dimslist= [];
      nargs= 0;
      while (more_args()) {
        ++nargs;
        grow, dimslist, next_arg();
      }
      what= [];
      if (nargs==0) {
        what= array(type);
      } else if (is_array(dimslist)
                 && (structof(dimslist)==long || structof(dimslist)==int)
                 && noneof(dimslist <= 0)) {
        if (nargs==1 && dimslist(1)==numberof(dimslist)-1) {
          what= array(type, dimslist);
        } else {
          dims= dimslist;
          ndims= numberof(dims);
          dimslist= array(long, ndims+1);
          dimslist(1)= ndims;
          dimslist(2:)= dims;
          what= array(type, dimslist);
        }
      }
    }
    if (! is_array(what)) error, "bad array specification";
  }
  if (_read(file, address, what) != numberof(what)) error, "short file";
  address+= sizeof(what);
  return what;
}




func gadget_write (s,fname)
 /* DOCUMENT
        writes a gadget structure into the file
        er= gadget_write(s,"test.dat");
     */      
{
  /* Open file. */
  file= open(fname, "wb");
  /* Writes the header. */
  address= 0; dummy=int(256);
  er=_gadget_write(file,dummy,&address);
  er=_gadget_write(file,s.npart,&address);
  er=_gadget_write(file,s.massarr,&address);
  er=_gadget_write(file,s.time,&address);
  er=_gadget_write(file,s.redshift,&address);
  er=_gadget_write(file,s.flag_sfr,&address);
  er=_gadget_write(file,s.flag_feedback,&address);
  er=_gadget_write(file,s.nparttotal,&address);
  er=_gadget_write(file,s.flag_cool,&address);
  er=_gadget_write(file,s.nfiles,&address);
  er=_gadget_write(file,s.boxsize,&address);
  er=_gadget_write(file,s.omegam,&address);
  er=_gadget_write(file,s.omegav,&address);
  er=_gadget_write(file,s.hubble,&address);
  for(i=address;i<260;i+=4)
  er=_gadget_write(file,dummy,&address);

  // Writes the data part. 
     ind= where((s.npart > 0) * (s.massarr == 0)) ;
  if (is_array(ind)) {
    nwithmass= sum(s.npart(ind));
  } else {
    nwithmass= 0;
  }
  n= sum(s.npart);

  er=_gadget_write(file,dummy,&address); 
  er=_gadget_write(file,dummy,&address);
  er=_gadget_write(file,*s.pos,&address);
  er=_gadget_write(file,dummy,&address);
  er=_gadget_write(file,dummy,&address);
  er=_gadget_write(file,*s.vel,&address);
  er=_gadget_write(file,dummy,&address);
  er=_gadget_write(file,dummy,&address);
  er=_gadget_write(file,*s.id,&address);
  er=_gadget_write(file,dummy,&address);
  
  if (nwithmass > 0) {   
    er=_gadget_write(file,dummy,&address);
    er=_gadget_write(file,*s.mass,&address);
    er=_gadget_write(file,dummy,&address);
  }

  ngas= s.npart(1);
  if (ngas > 0) {   
    er=_gadget_write(file,dummy,&address);
    er=_gadget_write(file,*s.u,&address);
    er=_gadget_write(file,dummy,&address);
    er=_gadget_write(file,dummy,&address);
    er=_gadget_write(file,*s.rho,&address);
    er=_gadget_write(file,dummy,&address);
  }
  close,file;
  return info(GadgetSnapshot);
}



func _gadget_write (file,what,&address)
/* DOCUMENT _gadget_write(file, address, expr)

     Write unformated, binary data from FILE.  
     On return, ADDRESS is incremented by the size of
     the data written.  The third argument gives the data to be written.
     
  SEE ALSO gadget_write,_gadget_read, open, _read. */
{
er=  _write(file,*address,what);
*address+=sizeof(what);
return er;
}



func gadgetcenter (q,pos,inplace=)
/* DOCUMENT
   recenters a given gadget snapshot on a given coordinate
   SEE ALSO:
 */
{ 
  x=(*q.pos)(1,);
  y=(*q.pos)(2,);
  z=(*q.pos)(3,);
  xmax=max(x);ymax=max(y);zmax=max(z);
  xmin=min(x);ymin=min(y);zmin=min(z);
  xmax-=xmin;
  ymax-=ymin;
  zmax-=zmin;

  zpos=pos-[xmax,ymax,zmax]/2.;
  x= (2*xmax+x-zpos(1)) % long(xmax);
  y= (2*ymax+y-zpos(2)) % long(ymax);
  z= (2*zmax+z-zpos(3)) % long(zmax);
  x += xmin;
  y += ymin;
  z += zmin;
  xyz= transpose([x,y,z]);
  x=y=z=[];
if(is_void(inplace)){
  q2=q;
  q2.pos=&xyz;
  return  q2;}
   else
     {
  q.pos=&xyz;

     }
}


func showxyz (x,y,z,txt=,nsqr=,r200=,center=,spin=,axis=,noerase=,color=)
/* DOCUMENT 
   display a cluster of points in xyz projections
   SEE ALSO:
 */
{
  if (typeof(x)=="struct_instance")
    {
       y= (*x.pos)(2,);
       z= (*x.pos)(3,);
       x= (*x.pos)(1,);
    } else
  if (dimsof(x)(1)==2)
    {
       y= (x)(2,);
       z= (x)(3,);
       x= (x)(1,);
     
    }

  require,"jgraph.i";
if(is_void(noerase)) window,style="win22.gs";
  plsys,1; pl,y,x; if(is_void(nsqr)) limits,square=1;
  if(!is_void(r200)) er=JDrawCircle(center(1),center(2),r200,color="red",width=4,number=30);
  if(!is_void(spin)) { plg,[center(2),center(2)+spin(2)],[center(1),center(1)+spin(1)],color="green",width=7;}
  if(!is_void(axis)) { plg,[center(2),center(2)+axis(2)],[center(1),center(1)+axis(1)],color="blue",width=7;}
 plsys,2; pl,z,x; if(is_void(nsqr)) limits,square=1;
  if(!is_void(r200)) er=JDrawCircle(center(1),center(3),r200,color="red",width=4,number=30);
  if(!is_void(spin)) { plg,[center(3),center(3)+spin(3)],[center(1),center(1)+spin(1)],color="green",width=7;}
  if(!is_void(axis)) { plg,[center(3),center(3)+axis(3)],[center(1),center(1)+axis(1)],color="blue",width=7;}
  plsys,3; pl,z,y;if(is_void(nsqr)) limits,square=1;
  if(!is_void(r200)) er=JDrawCircle(center(2),center(3),r200,color="red",width=4,number=30);
  if(!is_void(spin)) { plg,[center(3),center(3)+spin(3)],[center(2),center(2)+spin(2)],color="green",width=7;}
  if(!is_void(axis)) { plg,[center(3),center(3)+axis(3)],[center(2),center(2)+axis(2)],color="blue",width=7;}
  if (!is_void(txt)) {plsys,4; for(i=1;i<=numberof(txt);i++){ plt,txt(i),0.05,1 -(i+1)/float(2+numberof(txt)),tosys=1,height=8;}
   limits,0,1.5,0,1.1; }
else {plsys,4; plt,"x-y,x-z,\n y-z",0.375,0.85,tosys=1;}
}



func showsim (fname,imin=,imax=,pers=,vel=,dens=,delay=,jpg=,isovel=,
             xrange=,yrange=,zrange=,sub=,npix=,col=,iso=,dir=,center=,notxt=,cut=,nodisp=,
             marker=,comp=,rot=,flag=,color=,noerase=,ramses=)
/* DOCUMENT
   loads and   displays gadget binary format simulations
   
   imin= minimum indice
     imax= maximum indice
     pers= perspective plot
     vel= velocity plot (dodgy)
     dens= isodensity plot
     delay= delay between shots
     jpg= saves picture
     xrange= either e.g. 1e4 or [1e4,2e4]
     yrange= ""         (by default equal xrange)
     zrange= ""
     sub = if non zero plots one particle in ten
     npix = number of pixel for density maps
     col = color plot w.r.t. depth
     iso = 3D isocontour
     isovel = 3D isocontour of vx 1 vy 2 vz componant
     dir = directory for output
     center= position of new frame center
     notxt= if non void remove file name
     cut = if non void returns subsim
     nodisp = if non void doesn't plot anything
     comp=  component filter (if void is equal to [1,1,1,1,1] (all parts displayed);
            uses black blue red green cyan magenta)
     rot= rotation (0=[x,y], 1=[y,z], 2=[z,x] )
     flag= should contain a index range of particles to select
     color= obvious
     noerase= obvious
     SEE ALSO:
 */
{
  require,"Eric/system.i";
  require,"Eric/regul.i";

  if(!is_void(nodisp)) notxt=1;
  if(is_void(comp)) comp=indgen(6);
  comp0=(indgen(6)*0); comp0(comp)=1; comp=comp0;

  if (is_void(rot)) rot=0;
  xpos = 1+(rot %3);
  ypos = 1+((rot+1) %3);
  
  if((structof(fname)==string))   
    {  files = exec("ls "+pr1(fname)+"*");
    len= dimsof(files)(0);
    if (is_void(dir)) dir ="./";
    if (is_void(imin)) imin=1;
    if (is_void(imax)) imax=len;
    } else {imin=imax=1;}
  for(i=imin;i<=imax;i++) 
    {  if(structof(fname)==string){   q=gadget_read(files(i));}
    else {q=fname; files=[""];}

    if(!is_void(flag)) {
      a=(*q.pos)(,flag);      q.pos=&a;
      a=(*q.vel)(,flag);      q.vel=&a;
      //   a=(*q.mass)(flag);      q.mass=&a;
      q.npart=[0,numberof(flag),0,0,0,0];q.nparttotal=q.npart;
    }
    
    if (comp(sum)<6)
      { 
        a=[]; av=[];  am=[];
        for(k=1;k<=6;k++)
          {
            if ((q.npart*comp)(k))
              {
                grow,a,(*q.pos)(,(q.npart(1:((k-1)?k-1:1))(sum)+1):(q.npart(1:k)(sum)));
                grow,av,(*q.vel)(,(q.npart(1:((k-1)?k-1:1))(sum)+1):(q.npart(1:k)(sum)));
                grow,am,(*q.mass)((q.npart(1:((k-1)?k-1:1))(sum)+1):(q.npart(1:k)(sum)));
              }
          }
        q.pos=&a;
        q.vel=&av; 
        q.mass=&am;
        q.npart=q.npart*(comp>0);
        q.nparttotal=q.npart;
        q.id=&(indgen(q.npart(sum)));
      }
    ff =indgen((q.npart)(sum))*0+1;
    if (!is_void(center)) {gadgetcenter(q,center,inplace=1);}
    if(is_void(nodisp) && is_void(noerase) )    fma;
    nn=long(((q.npart*comp)(sum))^(1/3.))+1;
    if (!is_void(yrange) && is_void(xrange) )     xrange=yrange;
    if (!is_void(xrange) || !is_void(yrange) || !is_void(zrange))
      { if (!is_void(xrange) )
        if (dimsof(xrange)(1)==1)
          ffx=((*q.pos)(1,)<xrange(2))*((*q.pos)(1,)>xrange(1));
        else ffx=((*q.pos)(1,)<xrange);
      if (!is_void(xrange) && is_void(yrange) )     yrange=xrange;
      if (!is_void(yrange) )
        if (dimsof(yrange)(1)==1)
          ffy=((*q.pos)(2,)<yrange(2))*((*q.pos)(2,)>yrange(1));
        else ffy=((*q.pos)(2,)<yrange);
      if (!is_void(zrange) )
        if (dimsof(zrange)(1)==1)
          ffz=((*q.pos)(3,)<zrange(2))*((*q.pos)(3,)>zrange(1));
        else ffz=((*q.pos)(3,)<zrange);
      }
       
    if (!is_void(ffx)) ff =ff*ffx;
    if (!is_void(ffy)) ff = ff*ffy;
    if (!is_void(ffz)) ff = ff*ffz;
    w= where(ff);
       
    if (is_void(pers) && is_void(vel) && is_void(dens) && is_void(iso) &is_void(isovel))
      {
        if (is_void(nodisp)) if (is_void(col) && is_void(marker) )
          {
            if(comp(sum)==6){ pl,(*q.pos)(ypos,w),(*q.pos)(xpos,w),color=color;} else
              { n=q.npart;
              if (n(1) ) pl,(*q.pos)(ypos,1:n(1)),(*q.pos)(xpos,1:n(1)),
                           marker=marker,color="black";
              if (n(2))  pl,(*q.pos)(ypos,n(1)+1:n(1:2)(sum)),(*q.pos)(xpos,n(1)+1:n(1:2)(sum)),
                           marker=marker,color="blue";
              if (n(3))  pl,(*q.pos)(ypos,n(1:2)(sum)+1:n(1:3)(sum)),
                           (*q.pos)(xpos,n(1:2)(sum)+1:n(1:3)(sum)),
                           marker=marker,color="red";
              if (n(4))  pl,(*q.pos)(ypos,n(1:3)(sum)+1:n(1:4)(sum)),
                           (*q.pos)(xpos,n(1:3)(sum)+1:n(1:4)(sum)),
                           marker=marker,color="green";
              if (n(5))  pl,(*q.pos)(ypos,n(1:4)(sum)+1:n(1:5)(sum)),
                           (*q.pos)(xpos,n(1:4)(sum)+1:n(1:5)(sum)),
                           marker=marker,color="cyan";
              if (n(6))  pl,(*q.pos)(ypos,n(1:5)(sum)+1:n(1:6)(sum)),
                           (*q.pos)(xpos,n(1:5)(sum)+1:n(1:6)(sum)),
                           marker=marker,color="magenta";
              }
          }
 
        else if (!is_void(col)) {plcol,(*q.pos)(3,w),(*q.pos)(2,w),(*q.pos)(1,w);}
        else if (!is_void(marker)) { PL,(*q.pos)(2,w),(*q.pos)(1,w),msize=marker;}
        LL=max((*q.pos)(,w));
        if (!is_void(xrange))  if (dimsof(xrange)(1)==1)
          limits,xrange(1),xrange(2) ,yrange(1),yrange(2) 
            if(is_void(notxt)&& is_void(nodisp)) plt,files(i),LL/5.,LL/3.,color="blue",height=14,tosys=1;
        if (!is_void(jpg)) jpeg,dir+files(i)+".jpg",nodisp=1;
      }
    else
      if (!is_void(pers))
        { LL=max((*q.pos)(,w));
        plp3,(*q.pos)(3,w)(::4),(*q.pos)(1,w)(::4),(*q.pos)(2,w)(::4);
        if(is_void(notxt))    plt,files(i),-LL/2.,LL/2.,color="red",height=14,tosys=1;
        } 
      else             if (!is_void(vel)) 
        {
          vx2=reform((*q.vel)(1,),[3,nn,nn,nn]);
          vy2=reform((*q.vel)(2,),[3,nn,nn,nn]);
          vx2 /= (vx2)(*)(rms);
          vy2 /= (vy2)(*)(rms);
          //   if (!is_void(range)) {
          //     z=reform((*q.pos)(3,),[3,nn,nn,nn]);
          //     w =where(z<range);
          //  }
          xmesh=span(0,1,nn)(,-:1:nn);
          ymesh=span(0,1,nn)(-:1:nn,);
          ncouche=5; // coupe suivant Z
          plv,vel*vy2(,,1:5)(,,avg),vel*vx2(,,1:5)(,,avg),ymesh,xmesh,scale=0.035; 
          if(is_void(notxt)) plt,files(i),0.5,0.5,color="red",height=14,tosys=1;
        }
      else             if (!is_void(dens)) 
        {
          if (!is_void(w) )
            {
              x=(*q.pos)(xpos,w);
              y=(*q.pos)(ypos,w);
              if(!is_void(ramses)) rhos=(*q.rho)(w);
            } else
              {
                x=(*q.pos)(xpos,);
                y=(*q.pos)(ypos,);
                if(!is_void(ramses)) rhos=(*q.rho);
              }
           if(is_void(xrange)) x = (x-min(x))/(max(x)-min(x)); else  x = (x-xrange(1))/(xrange(dif));
           if(is_void(yrange)) y = (y-min(y))/(max(y)-min(y)); else  y = (y-yrange(1))/(yrange(dif));
          if (is_void(npix)) npix =50;
          if(is_void(ramses))    u=histo2d([x,y],span(-0.01,1.01,npix),span(-0.01,1.01,npix)); else
            u=histo2d([x,y],span(-0.01,1.01,npix),span(-0.01,1.01,npix),wght=rhos);
       if(is_void(nodisp)){   if(!is_void(xrange))
            {
              limits,xrange(1),xrange(2),yrange(1),yrange(2);
              plk,log(u+1),span(yrange(1),yrange(2),npix-1)(-:1:npix-1,),
                span(xrange(1),xrange(2),npix-1)(,-:1:npix-1);
            } else {plk,log(u+1);}
       
          if (!is_void(jpg)) jpeg,dir+files(i)+".jpg",nodisp=1;
          if(is_void(notxt)) plt,files(i),5,5,color="blue",height=14,tosys=1;
       }
       }
      else    if (!is_void(isovel)) 
        { 
          if (!is_void(w) )
            {
              x=(*q.pos)(1,w);
              y=(*q.pos)(2,w);
              z=(*q.pos)(3,w);
              vx=(*q.vel)(1,w);
              vy=(*q.vel)(2,w);
              vz=(*q.vel)(3,w);
            } else
              {
                x=(*q.pos)(1,);
                y=(*q.pos)(2,);
                z=(*q.pos)(3,);
                vx=(*q.vel)(1,);
                vy=(*q.vel)(2,);
                vz=(*q.vel)(3,);
              }
 if(is_void(xrange)) x = (x-min(x))/(max(x)-min(x)); else  x = (x-xrange(1))/(xrange(dif));
 if(is_void(yrange)) y = (y-min(y))/(max(y)-min(y)); else  y = (y-yrange(1))/(yrange(dif));
  if(is_void(yrange)) z = (z-min(z))/(max(z)-min(z)); else  z = (z-zrange(1))/(zrange(dif));
          if (is_void(npix)) npix =31;
          niso=numberof(isovel);
          if(niso>1 ) {
            vv1=[vx,vy,vz];// required since possible masking
            wght=vx*0.+1;
            for(i=1;i<=niso;i++)
              wght*=vv1(,isovel(i));
            uv=histo3d([x,y,z],span(-0.01,1.01,npix),span(-0.01,1.01,npix),span(-0.01,1.01,npix),wght=wght);
          }
          else{
            if (isovel==1 )  uv=histo3d([x,y,z],span(-0.01,1.01,npix),span(-0.01,1.01,npix),span(-0.01,1.01,npix),wght=vx);
            if (isovel==2)    uv=histo3d([x,y,z],span(-0.01,1.01,npix),span(-0.01,1.01,npix),span(-0.01,1.01,npix),wght=vy);
            if (isovel==3)    uv=histo3d([x,y,z],span(-0.01,1.01,npix),span(-0.01,1.01,npix),span(-0.01,1.01,npix),wght=vz);
          }
 
          u=histo3d([x,y,z],span(-0.01,1.01,npix),span(-0.01,1.01,npix),span(-0.01,1.01,npix));
 
          u=uv/(u+ (u==0));
          u1=smooth(log(abs(u)+1),2);
          if (!is_real(isovel))    vv=centile(u1,ucut=0.1)(2); else vv=isovel; 
          if(is_void(nodisp))   slice,u1,val=vv; 
          if (!is_void(jpg)) jpeg,dir+files(i)+".jpg",nodisp=1;
          LL=max((*q.pos)(,w));
          if(is_void(notxt))   plt,files(i),-LL/2.,LL/2.,color="red",height=18,tosys=1;
        }

      else    if (!is_void(iso)) 
        {
          if (!is_void(w) )
            {
              x=(*q.pos)(1,w);
              y=(*q.pos)(2,w);
              z=(*q.pos)(3,w);
            } else
              {
                x=(*q.pos)(1,);
                y=(*q.pos)(2,);
                z=(*q.pos)(3,);
              }
 if(is_void(xrange)) x = (x-min(x))/(max(x)-min(x)); else  x = (x-xrange(1))/(xrange(dif));
 if(is_void(yrange)) y = (y-min(y))/(max(y)-min(y)); else  y = (y-yrange(1))/(yrange(dif));
  if(is_void(yrange)) z = (z-min(z))/(max(z)-min(z)); else  z = (z-zrange(1))/(zrange(dif));
          if (is_void(npix)) npix =30;
          u=histo3d([x,y,z],span(-0.01,1.01,npix),span(-0.01,1.01,npix),span(-0.01,1.01,npix));
          u1=smooth(log(u+1),2);
          if (!is_real(iso))    vv=centile(u1,ucut=0.1)(2); else vv=iso; 
          if(is_void(nodisp))    slice,u1,val=vv; 
          if (!is_void(jpg)) jpeg,dir+files(i)+".jpg",nodisp=1;
          LL=max((*q.pos)(,w));
          if(is_void(notxt) )   plt,files(i),-LL/2.,LL/2.,color="red",height=18,tosys=1;
        }



    if (!is_void(delay)) pause,long(delay*1000);
    } // loop on snapshots
  if(is_void(nodisp)) animate,0;
  if (!is_void(dens) || !is_void(iso) || !is_void(isovel)) return u;
  if (!is_void(w) && !is_void(cut) )
    {
      x=(*q.pos)(1,w);
      y=(*q.pos)(2,w);
      z=(*q.pos)(3,w);
      vx=(*q.vel)(1,w);
      vy=(*q.vel)(2,w);
      vz=(*q.vel)(3,w);
      xyz=transpose([x,y,z]);
      vxvyvz=transpose([vx,vy,vz]);
      (q.pos)=&xyz;
      (q.vel)=&vxvyvz;
      (q.npart)(2) =numberof(w);
      return q;
    }

  return q;
}



func checksim (dir)
/* DOCUMENT 
     checks the length of files
   SEE ALSO:
 */
{ require,"Eric/system.i";
cd,dir;
 dir=get_cwd();
  ll=exec("ls -tdr G*/"); 
for(i=1;i<=numberof(ll);i++) {
  cd,ll(i);
rr= exec("\ls -l seed* | cut -c34-43 | sort | uniq | wc");
 if (strtok(rr)(1)=="1") {diag= "ok";} else {diag="bad!";} 

 rr= exec("\ls -l seed* | cut -c34-43 | sort | wc");
 if (strtok(rr)(1)=="37") {completed= "completed";} else {completed="in progress "+strtok(rr)(1);} 
 write,format= "%s %s %s \n ",ll(i),diag,completed;
  cd,dir;}
}



func readSnapTimings (dummy,dir=)
{
  if (is_void(dir)) dir="./";
  ll=exec( "\ls -l "+dir+ "seed* | cut -c49-55 | tr ':'  ' '");
  print,ll;
}



func readTimings (fname,nf=)
/* DOCUMENT  reads gadget Timing file 
   EXEMPLE  plg,readTimings("~/gadget/test/timings.txt")(dif);
   SEE ALSO asciiRead, 
*/
{
ff=open(fname,"r");
  xc=[];
 do
   {
     ll=rdline(ff);
     l1 = strtok(ll);

     if ((l1(1)=="Step=")*(is_void(nf))) {
       l1 = strtok(l1(2));
       l1 = strtok(l1(2));
       l1 = strtok(l1(2));
       x=0.; 
       sread,format="%f\n",l1(1),x;
           grow,xc,x;
     }
else 
  if ((l1(1)=="Nf=") *(!is_void(nf))) {
       l1 = strtok(l1(2));
       x=0.; 
       sread,format="%f\n",l1(1),x;
       grow,xc,x;
} // looks at neighbours 
   }
     while (ll)
   close,ff;
 return xc;
}

func showhaloenv (q,center,cut=,range=)
{
  if (is_void(range))  range=0.25e4;
  if (!is_void(cut)) q1=showsim(q,zrange=[center(3)-range,center(3)+range],
                                xrange=[center(1)-range,center(1)+range],
                                yrange=[center(2)-range,center(2)+range],cut=1);
else showsim,q,zrange=[center(3)-range,center(3)+range];
  ww=where(abs(center(3,)-center(3))<range);
  /*  plmk2,center(2,ww)(1:9),center(1,ww)(1:9),incolor=-5,width=10,msize=0.7,marker=4;
  plmk2,center(2,ww)(9:),center(1,ww)(9:),incolor=-6,width=10,msize=0.5,marker=4;
  */
  plmk2,center(2),center(1),incolor=-10,width=10,msize=1,marker=4;
  // limits,center(1)-1e4,center(1)+1e4,center(2)-1e4,center(2)+1e4;
  limits,center(1)-range,center(1)+range,center(2)-range,center(2)+range;
   return q1;
}


func showgal (q,rot=,comp=)
/* DOCUMENT
 showgal(q,rot=,comp=)  
 displays galaxy with componants using different colors
*/
{ ws;
 if (is_void(comp)) comp=[2,3,4]; 
  q=showsim(q,comp=comp,rot=rot);
  return q;
}


func _rotate (xyz0,rot)
/* DOCUMENT 
     
   SEE ALSO:
 */
{
  local x,y,z,xyz,ca,sa;
  rot=double(rot);

  xyz=xyz0;
  for(i=1;i<=numberof(rot);i++)
    { 
      x=xyz(((i-1) %3)+1,);
      y=xyz(((i-0) %3)+1,);
          z=xyz(((i+1) %3)+1,);
          ca= cos(rot(i));
          sa= sin(rot(i));
         if(abs(rot(i))>1e-5){
           a= x;
           x=  ca*a + sa*y;
           y= -sa*a + ca*y;
         }
         xyz=transpose([x,y,z]);
         
    }
  return xyz;
}

func gadget_merge (q1,q2,rot=,trans=,vtrans=,scale=,flip=)
/* DOCUMENT 
   returns the merge gadget snapshot with
   a possible rotation  rot=[anglexy,[angleyz],[anglezx]]
   translation trans=[ , ,]
   velocity translation vtrans=[ , ,]
   velocity flip = 1/0 or void
   SEE ALSO:
 */
{
  require,"pl3d.i";
  
  q=GadgetSnapshot();
  a=*q1.mass;
  grow,a,*q2.mass;
  q.mass=&a;
  a=*q1.pos;
  b=*q2.pos;
  av=*q1.vel;
  bv=*q2.vel;
  am=*q1.mass;
  bm=*q2.mass;
  
  q.npart=[0,
          q1.npart(2)+q2.npart(2),
          q1.npart(3)+q2.npart(3),
          q1.npart(4)+q2.npart(4),0,0];
  q.nparttotal=q.npart;

  ab=[];
  abv=[];
  abm=[];
  
  for(k=1;k<=6;k++){
    if((q1.npart(k)>0)||(q2.npart(k)>0))
      {
        if(q1.npart(k)>0)
          { a1=a(,(q1.npart(1:((k-1)?k-1:1))(sum)+1):(q1.npart(1:k)(sum)));
            av1=av(,(q1.npart(1:((k-1)?k-1:1))(sum)+1):(q1.npart(1:k)(sum)));
            am1=am((q1.npart(1:((k-1)?k-1:1))(sum)+1):(q1.npart(1:k)(sum)));
          }
        else {a1=[]; av1=[]; am1=[];}
        if(q2.npart(k)>0)
          {
            b1=b(,(q2.npart(1:((k-1)?k-1:1))(sum)+1):(q2.npart(1:k)(sum)));
             bv1=bv(,(q1.npart(1:((k-1)?k-1:1))(sum)+1):(q1.npart(1:k)(sum)));
             bm1=bm((q1.npart(1:((k-1)?k-1:1))(sum)+1):(q1.npart(1:k)(sum)));
          }
        else {b1=[]; bv1=[]; bm1=[];}                     

        grow,ab,a1;grow,abv,av1;grow,abm,am1;

        if (!is_void(flip)) bv1=-bv1;
        if (is_void(trans)) trans=[0,0,0];
        if (is_void(vtrans)) vtrans=[0,0,0];

        if (is_void(rot)) { grow,ab,b1+trans; grow,abv,bv1+vtrans; } else
             {
              grow,ab,_rotate(b1,rot)+trans; grow,abv,_rotate(bv1,rot)+vtrans;
             }
        
        
        grow,abm,bm1;
      }
  }
  q.pos=&ab;
  q.vel=&abv;
  q.mass=&abm;
  q.id=&(indgen(long(numberof(*q1.id)+numberof(*q2.id))));
    return q;
  }