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;
}
|