Commit 664b5408 authored by jclamber's avatar jclamber
Browse files

more work on rectify project

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@100 f264a43e-d52d-4b82-913b-c2bd5215a18a
parent 09a34876
......@@ -199,5 +199,6 @@ INSTALL(FILES ${PROJECT_SOURCE_DIR}/man/uns2uns.1 DESTINATION man/man1)
INSTALL(FILES ${PROJECT_SOURCE_DIR}/man/uns_projects.1 DESTINATION man/man1)
INSTALL(FILES ${PROJECT_SOURCE_DIR}/lib/utils/nemodep/cfalcon.h DESTINATION inc/uns)
INSTALL(FILES ${PROJECT_SOURCE_DIR}/lib/utils/csnaptools.h DESTINATION inc/uns)
INSTALL(FILES ${PROJECT_SOURCE_DIR}/lib/projects/nemodep/crectify.h DESTINATION inc/uns)
INSTALL(FILES ${PROJECT_SOURCE_DIR}/lib/utils/ctimer.h DESTINATION inc/uns)
#
......@@ -17,10 +17,14 @@
#include <cstdio>
#include "crectify.h"
extern "C" {
#include "nrutil.h"
//#include "nrutil.h"
float **fmatrix(long nrl, long nrh, long ncl, long nch);
float *fvector(long nl, long nh);
#include "nr.h"
}
#include <cvecutils.h>
#include <iomanip>
using namespace uns_proj;
using namespace jclut;
......@@ -36,13 +40,27 @@ CRectify::CRectify()
}
//------------------------------------------------------------------------------
// Rectify
// this method rectify position and vitesse for tge current time, according to
// the followings rules :
// _use_rho = true, use density as weighting factor, if _rho array is empty then
// local density is computed
// = false, use mass as weighting factor
// _threshold=when compute density, program will keep only particles above threshold
// which must be a value between 0.0% and 100.0%
// _rectify = true, positions and velocities will be rectified
// = false, vectors only are computed
// _cod_file= is a filename contening COD of snapshots
// if you don't provide this file, you must set _use_tho=true
// _rect_file=output file with results of the rectification, with the following format
// time 6xCOD 9xVectors
//
bool CRectify::rectify(const int _nbody,const float _time,
float * _pos, float * _vel,float * _mass, float * _rho,
const bool _use_rho, const bool _rectify,
std::string _cod_file,std::string _rect_file, const float _threshold)
{
bool ok = true;
nbody= _nbody;
pos = _pos;
vel = _vel;
......@@ -59,9 +77,9 @@ bool CRectify::rectify(const int _nbody,const float _time,
process();
return ok;
}
// PUBLIC CDataIndex
//------------------------------------------------------------------------------
//
......@@ -81,7 +99,8 @@ void CRectify::init()
}
//------------------------------------------------------------------------------
//
// process
// run all the process on snapshot
void CRectify::process()
{
......@@ -98,14 +117,16 @@ void CRectify::process()
if (rect) {
snapTransform();
}
saveRectVector();
}
//------------------------------------------------------------------------------
//
// findCenter
// center snapshots
void CRectify::findCenter()
{
bool is_codfile=false;
if (cod_file.length()!=0) { // read cod from file
is_codfile=CSnaptools::getTimeDataFile(cod_file,time,6,cod);
is_codfile=CSnaptools::getTimeDataFile<double>(cod_file,(double)time,6,cod);
if (! is_codfile ) {
std::cerr << "From file ["<<cod_file<<"] unable to find corresponding time ["<<time<<"]"
<< "aborting program....\n";
......@@ -133,7 +154,8 @@ void CRectify::findCenter()
}
//------------------------------------------------------------------------------
//
// processRho
// compute rho if requested
void CRectify::processRho()
{
......@@ -196,7 +218,8 @@ void CRectify::processRho()
}
//------------------------------------------------------------------------------
//
// findMoment
// compute moment
void CRectify::findMoment()
{
CLRM(w_qpole);
......@@ -217,13 +240,13 @@ void CRectify::findMoment()
ADDM(w_qpole, w_qpole, tmpm);
}
} else { // data are stored in arrays
for (unsigned int i=0; i<nbody; i++) {
for ( int i=0; i<nbody; i++) {
w_sum += mass[i];
}
float tmpv[3], pos_b[3];
matrix tmpm;
for (unsigned int i=0; i<nbody; i++) {
for ( int i=0; i<nbody; i++) {
float w_b = mass[i];
vectutils::subv(pos_b, &pos[i*3], codf); // shift to cod
MULVS(tmpv, pos_b, w_b);
......@@ -235,10 +258,10 @@ void CRectify::findMoment()
std::cerr << " find moment, done\n";
}
//------------------------------------------------------------------------------
//
// computeVectors
// compute Vectors
void CRectify::computeVectors()
{
//initOldFrame(); // !!!!!
// computing Eigen vectors according to moment w_qpole
eigenFrame(frame, w_qpole);
float tmp;
......@@ -258,12 +281,53 @@ void CRectify::computeVectors()
vectutils::setv(oldframe[i], frame[i]);
}
//------------------------------------------------------------------------------
//
// snapTransform STATIC
// rectify positions and velocities using rect_file
// this function returns :
// true if time value has been found in rect file
// false otherwise
bool CRectify::snapTransform(const int nbody,const float time,
float * pos, float * vel,
std::string rect_file, int &status)
{
bool ok=false;
if (!CSnaptools::isFileExist(rect_file)) {
std::cerr << "\nRequested rect file["<<rect_file<<"] does not exist, abort\n\n";
std::exit(1);
}
float data[15];
ok=CSnaptools::getTimeDataFile<float>(rect_file,time,15,data);
if (ok) {
float pos_b[3], vel_b[3];
// LOOP on all particles not only those from threshold !!!
for ( int i=0; i<nbody; i++) {
// center according to COD
vectutils::subv(&pos[i*3], &pos[i*3], data); // cod x data[0-2]
vectutils::subv(&vel[i*3], &vel[i*3], data+3); // cod v data[3-5]
// transfor according eigens vectors
for (int ii = 0; ii < NDIM; ii++) {
vectutils::dotvp(pos_b[ii],&pos[i*3], data+6+ii*3); // data[6-8,9-11,12-14]
vectutils::dotvp(vel_b[ii],&vel[i*3], data+6+ii*3); // data[6-8,9-11,12-14]
//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);
}
}
return ok;
}
//------------------------------------------------------------------------------
//snapTransform
// rectify positions and velocities using vectors
void CRectify::snapTransform()
{
float pos_b[3], vel_b[3];
// LOOP on all particles not only those from threshold !!!
for (unsigned int i=0; i<nbody; i++) {
for ( 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++) {
......@@ -278,7 +342,8 @@ void CRectify::snapTransform()
}
//------------------------------------------------------------------------------
//
// eigenFrame
// compute eigen values
void CRectify::eigenFrame(float frame[3][3], matrix mat)
{
float **q, *d, **v;
......@@ -297,7 +362,8 @@ void CRectify::eigenFrame(float frame[3][3], matrix mat)
frame[i-1][j-1] = v[j][i];
}
//------------------------------------------------------------------------------
//
// printVec
// print vectors
void CRectify::printvec(std::string name, float vec[3])
{
double PI=acos(-1.);
......@@ -319,3 +385,38 @@ void CRectify::xyz2rtp(float xyz[3], float rtp[3])
rtp[2] = atan2(xyz[1], xyz[0]); /* phi: in range -PI .. PI */
rtp[0] = sqrt(z*z+w*w);
}
//------------------------------------------------------------------------------
// saveRectVector
// save informations in rect file with the following format"
// time codx cody codz codvx codvy codvz v11 v12 v12 v21 v22 v23 v31 v32 v33
void CRectify::saveRectVector()
{
std::fstream fd; // file descriptor
if (rect_file.length() != 0) {
fd.open(rect_file.c_str(),std::ios::in | std::ios::out | std::ios::app); // open file for appending
if (fd.is_open()) {
//fd << scientific << left << time;
std::ostringstream ss;
// save time
ss << scientific << left << time << " ";
// save cod
for (int i=0; i<6;i++ ){
ss << codf[i] << " ";
}
// save vectors
for (int i=0; i<3;i++ ){
ss << frame[i][0] << " " << frame[i][1] << " " << frame[i][2] << " ";
}
ss << "\n";
fd << ss.str(); // save one line
fd.close();
} else {
std::cerr << "Unable to open file ["<< rect_file.c_str()<<"] in appending mode, abort...\n";
std::exit(1);
}
}
}
......@@ -14,6 +14,11 @@
@author Jean-Charles Lambert <Jean-Charles.Lambert@lam.fr>
*/
/*
This class aims to rectify a galaxy in XY plan. It's based on snaprect program
from the NEMO project, see:
http://carma.astro.umd.edu/nemo/man_html/snaprect.1.html
*/
#ifndef CRECTIFY_H
#define CRECTIFY_H
#include <string>
......@@ -23,6 +28,7 @@
extern "C" {
#include <vectmath.h>
}
#include "csnaptools.h"
namespace uns_proj {
......@@ -49,13 +55,36 @@ public:
class CRectify {
public:
CRectify();
//------------------------------------------------------------------------------
// Rectify
// this method rectify position and vitesse for tge current time, according to
// the followings rules :
// _use_rho = true, use density as weighting factor, if _rho array is empty then
// local density is computed
// = false, use mass as weighting factor
// _threshold=when compute density, program will keep only particles above threshold
// which must be a value between 0.0% and 100.0%
// _rectify = true, positions and velocities will be rectified
// = false, vectors only are computed
// _cod_file= is a filename contening COD of snapshots
// if you don't provide this file, you must set _use_tho=true
// _rect_file=output file with results of the rectification, with the following format
// time 6xCOD 9xVectors
//
bool rectify(const int nbody,const float time,
float * pos, float * vel,float * mass, float * weight,
const bool _use_rho, const bool rectify=false,
std::string cod_file="",std::string rect_file="", const float threshold=0.0);
void process();
//------------------------------------------------------------------------------
// snapTransform STATIC
// rectify positions and velocities using rect_file
// this function returns :
// true if time value has been found in rect file
// false otherwise
static bool snapTransform(const int nbody,const float time,
float * pos, float * vel,std::string rect_file, int &status);
void initOldFrame() {
float init_frame[][3] = {
{ 1.0, 0.0, 0.0, },
......@@ -82,6 +111,7 @@ private:
void eigenFrame(float frame[3][3], matrix mat);
void xyz2rtp(float xyz[3], float rtp[3]);
void printvec(std::string name, float vec[3]);
void saveRectVector();
matrix w_qpole;
double cod[6];
......@@ -91,4 +121,15 @@ private:
};
}
extern "C" {
bool rectify_snap_(const int * nbody,const float * time,
float * pos, float * vel,
const char * rect_file, int *status, const int lenstring) {
std::string filename=jclut::CSnaptools::fixFortran(rect_file,lenstring);
return uns_proj::CRectify::snapTransform(*nbody,*time,pos,vel,filename,*status);
}
}
#endif // CRECTIFY_H
......@@ -41,7 +41,7 @@ bool CSnaptools::isFileExist(std::string test_file)
//
// ----------------------------------------------------------------------------
// getTimeDataFile
bool CSnaptools::getTimeDataFile(std::string input_file,const float time,const int n,double data[],const double offset,const bool verbose)
template <class T> bool CSnaptools::getTimeDataFile(std::string input_file,const T time,const int n,T data[],const T offset,const bool verbose)
{
bool ok=false;
std::ifstream fd; // file descriptor
......@@ -75,6 +75,8 @@ bool CSnaptools::getTimeDataFile(std::string input_file,const float time,const i
return ok;
}
template bool CSnaptools::getTimeDataFile<float>(std::string input_file,const float time,const int n,float data[],const float offset,const bool verbose);
template bool CSnaptools::getTimeDataFile<double>(std::string input_file,const double time,const int n,double data[],const double offset,const bool verbose);
//
// ----------------------------------------------------------------------------
// fixFortran
......
......@@ -28,7 +28,8 @@ 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 getTimeDataFile(std::string _file,const float time,const int n,double data[],const double offset=0.0001,const bool verbose=false);
template <class T> static bool getTimeDataFile(std::string _file,const T time,const int n,
T data[],const T offset=0.0001,const bool verbose=false);
static bool isFileExist(std::string);
static std::string basename(const std::string);
static std::string dirname(const std::string);
......
......@@ -35,7 +35,7 @@ const char * defv[] = {
"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",
"threshold=0.0\n percentage of density kept >= threshold",
"times=all\n selected time",
"first=f\n add a trailing numbering to the first output file",
"verbose=f\n verbose on/off",
......
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