Commit 32c7143a authored by jclamber's avatar jclamber
Browse files

fixes

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@114 f264a43e-d52d-4b82-913b-c2bd5215a18a
parent 6d50c15c
......@@ -112,7 +112,6 @@ void CRectify::process()
}
if (density) delete density;
findCenter();
findMoment();
computeVectors();
......@@ -211,7 +210,7 @@ void CRectify::processRho()
if (rho) {
vrho.push_back(rho[sindex]);
} else {
vrho.push_back(density->getRho()[sindex]);
vrho.push_back(density->getRho()[ii]); //[sindex]);
}
ii++;
}
......@@ -230,16 +229,10 @@ void CRectify::processRho()
void CRectify::findMoment()
{
CLRM(w_qpole);
float w_sum=0.0;
double w_sum=0.0;
float distance=0.0;
if ( use_rho) { // data are stored in Vectors
for (unsigned int i=0; i<vmass.size(); i++) {
ABSV(distance,&vpos[i*3]);
if (distance <=radius) {
w_sum += vrho[i] * vmass[i];
}
}
float tmpv[3], pos_b[3];
matrix tmpm;
......@@ -247,6 +240,7 @@ void CRectify::findMoment()
vectutils::subv(pos_b, &vpos[i*3], codf); // shift to cod
ABSV(distance,pos_b);
if (distance <=radius) {
w_sum += (vrho[i] * vmass[i]);
float w_b = vrho[i] * vmass[i];
//vectutils::subv(pos_b, &vpos[i*3], codf); // shift to cod
MULVS(tmpv, pos_b, w_b);
......@@ -255,12 +249,6 @@ void CRectify::findMoment()
}
}
} else { // data are stored in arrays
for ( int i=0; i<nbody; i++) {
ABSV(distance,&pos[i*3]);
if (distance <=radius) {
w_sum += mass[i];
}
}
float tmpv[3], pos_b[3];
matrix tmpm;
......@@ -268,6 +256,7 @@ void CRectify::findMoment()
vectutils::subv(pos_b, &pos[i*3], codf); // shift to cod
ABSV(distance,pos_b);
if (distance <=radius) {
w_sum += mass[i];
float w_b = mass[i];
//vectutils::subv(pos_b, &pos[i*3], codf); // shift to cod
MULVS(tmpv, pos_b, w_b);
......
......@@ -134,6 +134,8 @@ int main(int argc, char ** argv )
unsout->snapshot->setData("time",time);
unsout->snapshot->setData("pos" ,nbody,pos,false);
if (vel)
unsout->snapshot->setData("vel" ,nbody,vel,false);
unsout->snapshot->setData("mass" ,nbody,mass,false);
if (rho)
unsout->snapshot->setData("rho" ,nbody,rho ,false);
......
......@@ -33,9 +33,9 @@ const char * defv[] = {
"out=\n NEMO output file ",
"select=???\n component selected (disk,stars,halo,gas,range,all",
"rectf=\n rectified info file",
"rho=t\n use rho as weight factor, otherwise mass",
"rho=f\n use rho as weight factor, otherwise mass",
"codfile=\n use cod from file to recenter",
"radius=50.0\n cutting radius",
"rcut=50.0\n cutting radius",
"dmin=0.0\n keep particles if %log(density) >= dmin",
"dmax=100.0\n keep particles if %log(density) <= dmax",
"times=all\n selected time",
......@@ -65,7 +65,7 @@ int main(int argc, char ** argv )
bool use_rho (getbparam((char *) "rho" ));
std::string select_c (getparam ((char *) "select" ));
std::string select_t (getparam ((char *) "times" ));
float radius(getdparam ((char *) "radius"));
float radius(getdparam ((char *) "rcut"));
float dmin(getdparam ((char *) "dmin"));
float dmax(getdparam ((char *) "dmax"));
bool first (getbparam((char *) "first" ));
......
......@@ -49,7 +49,7 @@ void readtipsystd_(const char * filename,std::string outname,float scale)
fp = fopen(filename,"r");
xdrstdio_create(&xdr,fp,XDR_DECODE);
/* read header */
xdr_header(&xdr,&head);
std::cerr << "Scale = "<<scale <<"\n";
......@@ -71,7 +71,9 @@ void readtipsystd_(const char * filename,std::string outname,float scale)
bool ok;
// set time
ok=unsout->setData("time",time);
ok=unsout->snapshot->setData("time",time);
if (ok) ; // remove compiler warning
/* reads data for gas particles */
......@@ -104,11 +106,11 @@ void readtipsystd_(const char * filename,std::string outname,float scale)
}
if (head.nsph) {
ok=unsout->setData("gas","pos" ,head.nsph,xx ,false);
ok=unsout->setData("gas","vel" ,head.nsph,vv ,false);
ok=unsout->setData("gas","mass",head.nsph,mass,false);
ok=unsout->setData("gas","rho" ,head.nsph,rho ,false);
ok=unsout->setData("gas","hsml",head.nsph,hsml,false);
ok=unsout->snapshot->setData("gas","pos" ,head.nsph,xx ,false);
ok=unsout->snapshot->setData("gas","vel" ,head.nsph,vv ,false);
ok=unsout->snapshot->setData("gas","mass",head.nsph,mass,false);
ok=unsout->snapshot->setData("gas","rho" ,head.nsph,rho ,false);
ok=unsout->snapshot->setData("gas","hsml",head.nsph,hsml,false);
delete [] mass;
delete [] xx;
delete [] vv;
......@@ -141,9 +143,9 @@ void readtipsystd_(const char * filename,std::string outname,float scale)
}
if (head.ndark) {
ok=unsout->setData("halo","pos" ,head.ndark,xx ,false);
ok=unsout->setData("halo","vel" ,head.ndark,vv ,false);
ok=unsout->setData("halo","mass",head.ndark,mass,false);
ok=unsout->snapshot->setData("halo","pos" ,head.ndark,xx ,false);
ok=unsout->snapshot->setData("halo","vel" ,head.ndark,vv ,false);
ok=unsout->snapshot->setData("halo","mass",head.ndark,mass,false);
delete [] mass;
delete [] xx;
delete [] vv;
......@@ -173,9 +175,9 @@ void readtipsystd_(const char * filename,std::string outname,float scale)
}
if (head.nstar) {
ok=unsout->setData("stars","pos" ,head.nstar,xx ,false);
ok=unsout->setData("stars","vel" ,head.nstar,vv ,false);
ok=unsout->setData("stars","mass",head.nstar,mass,false);
ok=unsout->snapshot->setData("stars","pos" ,head.nstar,xx ,false);
ok=unsout->snapshot->setData("stars","vel" ,head.nstar,vv ,false);
ok=unsout->snapshot->setData("stars","mass",head.nstar,mass,false);
delete [] mass;
delete [] xx;
delete [] vv;
......@@ -185,7 +187,7 @@ void readtipsystd_(const char * filename,std::string outname,float scale)
//printf("Data for %d particles read from %s \n",i,filename);
unsout->save();
unsout->snapshot->save();
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment