Commit 5f2a943e authored by jclamber's avatar jclamber
Browse files

new rectify based snaprec from NEMO

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@93 f264a43e-d52d-4b82-913b-c2bd5215a18a
parent 3e9eafe9
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2014
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 6110
// ============================================================================
/*
@author Jean-Charles Lambert <Jean-Charles.Lambert@lam.fr>
*/
#include "csnaptools.h"
#include <cstdio>
#include "crectify.h"
extern "C" {
#include "nrutil.h"
#include "nr.h"
}
#include <cvecutils.h>
using namespace uns_proj;
using namespace jclut;
using namespace std;
using namespace uns;
// PUBLIC CRectify
CRectify::CRectify(uns::CunsIn * _uns)
{
unsin = _uns;
init();
}
CRectify::CRectify()
{
init();
}
bool CRectify::rectify(int _nbody, float * _pos, float * _vel,float * _mass, float * _rho,
float _threshold, std::string _cod_file)
{
nbody= _nbody;
pos = _pos;
vel = _vel;
mass = _mass;
rho = _rho;
threshold = _threshold;
cod_file = _cod_file;
process();
}
bool CRectify::rectify(std::string component,float threshold, std::string cod_file)
{
}
// PUBLIC CDataIndex
bool CDataIndex::sortData(const CDataIndex &a, const CDataIndex &b)
{
return a.data < b.data;
}
// PRIVATE
void CRectify::init()
{
density=NULL;
nbody=0;
}
void CRectify::process()
{
if (density) delete density;
findCenter();
findMoment();
snapTransform();
}
void CRectify::findCenter()
{
if (!rho) {
// Instantiate a density object
density = new CDensity(nbody,&pos[0],&mass[0]);
density->compute(0,32,1,8); // estimate density
}
rho_di.clear();
// store rho in structure
for (int i=0; i<nbody; i++) {
CDataIndex p;
p.foo();
if (rho) {
p.setDI(rho[i],i);
} else {
p.setDI(density->getRho()[i],i);
}
rho_di.push_back(p);
}
// sort by rho
std::sort(rho_di.begin(),rho_di.end(),CDataIndex::sortData);
double minrho=log(rho_di[0].data);
double maxrho=log(rho_di[rho_di.size()-1].data);
vpos.clear();
vvel.clear();
vmass.clear();
vrho.clear();
// keep only particles > threshold
int ii=0;
for (std::vector<CDataIndex>::iterator it=rho_di.begin(); it != rho_di.end(); it++) {
float logrho = log((*it).data);
if (((logrho - minrho)*100./(maxrho-minrho))>=threshold) { // we keep particles
int sindex=(*it).index;
vpos.push_back(pos[sindex*3+0]);
vpos.push_back(pos[sindex*3+1]);
vpos.push_back(pos[sindex*3+2]);
if (vel) {
vvel.push_back(vel[sindex*3+0]);
vvel.push_back(vel[sindex*3+1]);
vvel.push_back(vel[sindex*3+2]);
}
vmass.push_back(mass[sindex]);
vrho.push_back(rho[sindex]);
ii++;
}
}
std::cerr << "Nbody = "<< nbody << " sel=" << ii << "\n";
// compute COD without shifting
CSnaptools::moveToCod<float>(vmass.size(),&vpos[0],&vvel[0],&vmass[0],&vrho[0],cod,false); // compute cod don't move
for (int i=0; i<6; i++) {
codf[i] = (float) cod[i];
}
}
void CRectify::findMoment()
{
CLRM(w_qpole);
float w_sum=0.0;
for (unsigned int i=0; i<vmass.size(); i++) {
w_sum += vrho[i] * vmass[i];
}
float tmpv[3], pos_b[3];
matrix tmpm;
for (unsigned int i=0; i<vmass.size(); 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);
OUTVP(tmpm, tmpv, pos_b);
ADDM(w_qpole, w_qpole, tmpm);
}
DIVMS(w_qpole, w_qpole, w_sum);
std::cerr << " find moment, done\n";
}
void CRectify::snapTransform()
{
float oldframe[][3] = {
{ 1.0, 0.0, 0.0, },
{ 0.0, 1.0, 0.0, },
{ 0.0, 0.0, 1.0, }};
float pos_b[3], vel_b[3];
eigenFrame(frame, w_qpole);
float tmp;
vectutils::dotvp(tmp,oldframe[0], frame[0]);
if (tmp < 0.0)
vectutils::mulvs(frame[0], frame[0], (float) -1.0);
vectutils::dotvp(tmp,oldframe[2], frame[2]);
if (tmp < 0.0)
vectutils::mulvs(frame[2], frame[2], (float )-1.0);
CROSSVP(frame[1], frame[2], frame[0]);
printvec("e_x:", frame[0]);
printvec("e_y:", frame[1]);
printvec("e_z:", frame[2]);
//for (b = btab; b < btab+nbody; b++) {
// !!!!
// LOOP on all particles not only those from threshold !!!
for (unsigned int i=0; i<nbody; i++) {
vectutils::subv(&pos[i*3], &pos[i*3], codf);
vectutils::subv(&vel[i*3], &vel[i*3], codf+3);
for (int ii = 0; ii < NDIM; ii++) {
vectutils::dotvp(pos_b[ii],&pos[i*3], frame[ii]);
vectutils::dotvp(vel_b[ii],&vel[i*3], frame[ii]);
//acc_b[i] = dotvp(Acc(b), frame[i]);
}
vectutils::setv(&pos[i*3], pos_b);
vectutils::setv(&vel[i*3], vel_b);
//SETV(Acc(b), acc_b);
}
for (int i = 0; i < NDIM; i++)
vectutils::setv(oldframe[i], frame[i]);
}
void CRectify::eigenFrame(float frame[3][3], matrix mat)
{
float **q, *d, **v;
int i, j, nrot;
q = fmatrix(1, 3, 1, 3);
for (i = 1; i <= 3; i++)
for (j = 1; j <= 3; j++)
q[i][j] = mat[i-1][j-1];
d = fvector(1, 3);
v = fmatrix(1, 3, 1, 3);
jacobi(q, 3, d, v, &nrot);
eigsrt(d, v, 3);
for (i = 1; i <= 3; i++)
for (j = 1; j <= 3; j++)
frame[i-1][j-1] = v[j][i];
}
void CRectify::printvec(std::string name, float vec[3])
{
double PI=acos(-1.);
float rtp[3]; /* radius - theta - phi */
xyz2rtp(vec,rtp);
fprintf(stderr,"%12s %10.5f %10.5f %10.5f %10.5f %5.1f %6.1f\n",
name.c_str(), rtp[0], vec[0], vec[1], vec[2],
rtp[1]*180.0/PI, rtp[2]*180.0/PI);
}
void CRectify::xyz2rtp(float xyz[3], float rtp[3])
{
double PI=acos(-1.);
float z = xyz[2];
float w = sqrt((xyz[0]*xyz[0])+(xyz[1]*xyz[1]));
rtp[1] = atan(w/z); /* theta: in range 0 .. PI */
if (z<0) rtp[1] += PI;
rtp[2] = atan2(xyz[1], xyz[0]); /* phi: in range -PI .. PI */
rtp[0] = sqrt(z*z+w*w);
}
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2014
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 6110
// ============================================================================
/*
@author Jean-Charles Lambert <Jean-Charles.Lambert@lam.fr>
*/
#ifndef CRECTIFY_H
#define CRECTIFY_H
#include <string>
#include "uns.h"
#include "cfalcon.h"
#define SINGLEPREC
#define real float
extern "C" {
#include <vectmath.h>
}
namespace uns_proj {
class CDataIndex {
public:
CDataIndex() {}
CDataIndex(const float _data, const int _index){
setDI(_data, _index);
}
void foo() {
}
void setDI(float _data, int _index ) {
data = _data;
index = _index;
}
static bool sortData(const CDataIndex& a, const CDataIndex& b);
float data;
int index;
};
class CRectify {
public:
CRectify(uns::CunsIn * uns);
CRectify();
bool rectify(std::string component,float theshold=70.0, std::string cod_file="");
bool rectify(int nbody,float * pos, float * vel,float * mass, float * rho,
float threshold=70.0, std::string cod_file="");
void process();
private:
uns::CunsIn * unsin;
jclut::CDensity * density;
float * pos, * vel, * mass, * rho;
float threshold;
std::string cod_file;
int nbody;
std::vector<CDataIndex> rho_di;
void init();
void findCenter();
void findMoment();
void snapTransform();
void eigenFrame(float frame[3][3], matrix mat);
void xyz2rtp(float xyz[3], float rtp[3]);
void printvec(std::string name, float vec[3]);
matrix w_qpole;
double cod[6];
float codf[6];
float frame[3][3];
std::vector <float> vpos,vvel,vmass,vrho;
};
}
#endif // CRECTIFY_H
......@@ -10,6 +10,7 @@
// ============================================================================
#include <iostream>
#include <fstream> // C++ file I/O
#include <cmath>
#include <sstream>
#include <assert.h>
......@@ -19,6 +20,41 @@
using namespace jclut;
//
// ----------------------------------------------------------------------------
// getCodFromFile
bool CSnaptools::getCodFromFile(std::string cod_file,const float time,double cod[6], double offset,bool verbose)
{
bool ok=false;
std::ifstream fd; // file descriptor
fd.open(cod_file.c_str(), std::ios::in); // open file
if (fd.is_open()) {
bool stop=false;
std::string line;
while (!stop && ! fd.eof() ) {
// Read line
getline(fd,line);
if (! fd.eof()) {
double r_time;
std::istringstream ss(line);
ss >> r_time;
if (r_time-offset <= time && time <= r_time+offset) { // found good time
stop = true;
ok = true;
for (int i=0; i<6; i++) { // read cod
ss >> cod[i];
}
}
}
}
}
return ok;
}
//
// ----------------------------------------------------------------------------
// fixFortran
......
......@@ -28,6 +28,7 @@ namespace jclut {
template <class T> static void moveToCod(const int nbody,T * pos, T * Vel, T * mass, T * rho, double cod[6], bool move, bool verbose=false);
template <class T> static void moveToCom(const int nbody,T * pos, T * mass, bool verbose=false);
static bool getCodFromFile(std::string cod_file,const float time,double cod[6],double offset=0.0000001, bool verbose=false);
static std::string basename(const std::string);
static std::string dirname(const std::string);
static std::string parseString(std::string & next_string, const std::string sep=",");
......
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2010-2014
// Centre de donneeS Astrophysiques de Marseille (CeSAM)
// e-mail: Jean-Charles.Lambert@lam.fr
// address: Aix Marseille Universite, CNRS, LAM
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS UMR 7326
// ============================================================================
#include <iostream> // C++ I/O
#include <fstream> // C++ file I/O
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <assert.h>
#include <cmath>
#include <nemo.h>
#include <iomanip>
// ------------------------------------------------------------
// Include file
#include "uns.h"
#include "crectify.h"
// ------------------------------------------------------------
// Nemo variable
const char * defv[] = {
"in=???\n UNS input snapshot",
"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",
"codfile=\n use cod from file to recenter",
"threshold=70.0\n percentage of density kept >= threshold",
"times=all\n selected time",
"verbose=f\n verbose on/off",
"VERSION=1.0\n compiled on <"__DATE__"> JCL ",
NULL,
};
const char * usage="Rectify an UNS snapshot";
using namespace uns_proj;
//------------------------------------------------------------------------------
// M A I N
//------------------------------------------------------------------------------
// main
int main(int argc, char ** argv )
{
// start NEMO
initparam(const_cast<char**>(argv),const_cast<char**>(defv));
if (argc) {;} // remove compiler warning :)
// Get input parameters
std::string simname (getparam ((char *) "in" ));
std::string outname (getparam ((char *) "out" ));
std::string codf (getparam ((char *) "codfile" ));
std::string rectf (getparam ((char *) "rectf" ));
bool rho (getbparam((char *) "rho" ));
std::string select_c (getparam ((char *) "select" ));
std::string select_t (getparam ((char *) "times" ));
float threshold(getdparam ((char *) "threshold"));
bool verbose = getbparam((char *) "verbose" );
std::cerr << "codfile["<<codf<<"] length="<< codf.length()<<"\n";
std::exit(1);
// -----------------------------------------------
// instantiate a new UNS input object (for reading)
uns::CunsIn * uns = new uns::CunsIn(simname,select_c,select_t,verbose);
if (uns->isValid()) { // input file is known by UNS lib
while(uns->snapshot->nextFrame()) { // there is a new frame
std::cerr << "Input file is of type :"<<uns->snapshot->getInterfaceType()<<"\n";
bool ok;
int cnbody,nbody;
float * pos=NULL, * vel=NULL, * mass=NULL, * rho=NULL, * hsml=NULL,time;
// get the input number of bodies according to the selection
ok=uns->snapshot->getData("nsel",&nbody);
assert(ok==true);
// get the simulation time
ok=uns->snapshot->getData("time",&time);
// get POS from input snapshot
ok=uns->snapshot->getData("pos" ,&cnbody,&pos);
assert(ok==true);
// get VEL from input snapshot
ok=uns->snapshot->getData("vel" ,&cnbody,&vel);
// get MASS from input snapshot
ok=uns->snapshot->getData("mass",&cnbody,&mass);
assert(ok==true);
// get RHO if exist
ok=uns->snapshot->getData("rho",&cnbody,&rho);
// get HSML if exist
ok=uns->snapshot->getData("hsml",&cnbody,&hsml);
std::cerr << "nbody=" << nbody << " time="<<time <<"\n";
CRectify * crectify = new CRectify();
crectify->rectify(nbody,pos,vel,mass,rho, threshold);
uns::CunsOut * unsout = new uns::CunsOut(outname,"nemo",verbose);
unsout->snapshot->setData("time",time);
unsout->snapshot->setData("pos" ,nbody,pos,false);
unsout->snapshot->setData("mass" ,nbody,mass,false);
if (rho)
unsout->snapshot->setData("rho" ,nbody,rho ,false);
if (hsml)
unsout->snapshot->setData("hsml",nbody,hsml,false);
unsout->snapshot->save();
delete unsout;
}
}
}
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