Commit c1299ea5 authored by jclamber's avatar jclamber

remove projects from the trunk

git-svn-id: http://svn.oamp.fr/repos/unsio/trunk@105 ce2cc22f-6084-46ce-a062-084b172ee5dc
parent 678963fa
......@@ -324,8 +324,10 @@ int CAmr::loadData(uns::CParticles * particles,
if (take) {
particles->indexes.push_back(0); // GAS particles
particles->ngas++; // One more gas particles
particles->ntot++; // one more total particles
}
particles->ntot++; // one more total particles
// this variable count all particles
#endif
}
nbody++;
......
......@@ -152,14 +152,13 @@ int CPart::loadData(uns::CParticles * particles,
}
if (take && agetmp[k]!=0) {
particles->indexes.push_back(4); // save star positions
particles->nstars++;
particles->ntot++; // one more total particles
particles->nstars++;
}
if (take && agetmp[k]==0) {
particles->indexes.push_back(1); // save DM positions
particles->ndm++;
particles->ntot++; // one more total particles
}
particles->ntot++; // one more total particles
offset++;
//assert(offset<=nselect);
......@@ -192,9 +191,9 @@ int CPart::loadData(uns::CParticles * particles,
}
if (take) {
particles->indexes.push_back(1); // save DM positions
particles->ndm++;
particles->ntot++; // one more total particles
particles->ndm++;
}
particles->ntot++; // one more total particles
offset++;
//assert(offset<=nselect);
}
......
......@@ -122,6 +122,7 @@ int CSnapshotRamsesIn::nextFrame(uns::UserSelection &user_select)
// close operation
int CSnapshotRamsesIn::close()
{
return 1;
}
// ============================================================================
// reorderParticles
......@@ -133,7 +134,7 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
std::cerr <<"Nbody particles loaded="<<particles->ntot<<"\n";
std::vector <int> offset_comp(6,-1); // init 6 offsets with value =-1
std::vector <int> npart_comp(6,0); // initialise #part per component to ZERO
char * comp[] = { "gas","halo","disk","bulge","stars","bndry","all" };
char * comp[] = { (char*) "gas",(char*) "halo",(char*) "disk",(char*) "bulge",(char*) "stars",(char*) "bndry",(char*) "all" };
// get ordering
std::vector <int> select_order = user_select.selectOrder();
......@@ -150,7 +151,7 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
bool first=true;
// according to user selection
// we reformat offset of each components
for (int i=0;i<select_order.size(); i++) {
for (unsigned int i=0;i<select_order.size(); i++) {
assert(select_order[i]<7);
if (first) { // first in the order
if (npart_comp[select_order[i]] > 0 ) { // part exist for the component
......@@ -185,8 +186,8 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
// insert "all" at the beginning of componentRangeVector
crv.insert(it,cr);
for (int i=0;i<offset_comp.size(); i++) {
std::cerr << "i="<<i<<" comp offset="<<offset_comp[i]<<"\n";
for (unsigned int i=0;i<offset_comp.size(); i++) {
std::cerr << "i="<<i<<" ***comp offset="<<offset_comp[i]<<"\n";
}
if (verbose)
uns::ComponentRange::list(&crv);
......@@ -211,7 +212,7 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
// positions
if (particles->pos.size()>0) {
assert((istart*3)+2<particles->pos.size());
assert((istart*3)+2<(int)particles->pos.size());
found=true;
pos[istart*3+0] = particles->pos[i*3+0]; // x
pos[istart*3+1] = particles->pos[i*3+1]; // y
......@@ -220,7 +221,7 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
// velocities
if (particles->vel.size()>0) {
assert((istart*3)+2<particles->vel.size());
assert((istart*3)+2<(int)particles->vel.size());
found=true;
vel[istart*3+0] = particles->vel[i*3+0]; // x
vel[istart*3+1] = particles->vel[i*3+1]; // y
......@@ -229,7 +230,7 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
// masses
if (particles->mass.size()>0) {
assert(istart<particles->mass.size());
assert(istart<(int)particles->mass.size());
found=true;
mass[istart] = particles->mass[i];
}
......@@ -243,7 +244,7 @@ int CSnapshotRamsesIn::reorderParticles(uns::UserSelection &user_select)
particles->vel=vel;
particles->mass=mass;
}
return 1;
}
// ============================================================================
// getData
......@@ -357,6 +358,91 @@ bool CSnapshotRamsesIn::getData(const std::string comp, std::string name, int *n
bool CSnapshotRamsesIn::getData(const std::string name,int *n,float **data)
{
bool ok=false;
*data = NULL;
*n = 0;
int first=0;
bool status=true;
switch(CunsOut::s_mapStringValues[name]) {
case uns::Nsel :
if (status) {
*n = particles->ntot;;
} else {
ok=false;
}
case uns::Pos :
if (particles->pos.size()>0) {
*data = &particles->pos[first*3];
*n = particles->pos.size()/3;
} else {
ok=false;
}
break;
case uns::Vel :
if (particles->vel.size()>0) {
*data = &particles->vel[first*3];
*n = particles->vel.size()/3;
} else {
ok=false;
}
break;
case uns::Mass :
if (particles->mass.size()) {
*data = &particles->mass[first];
*n = particles->mass.size();
} else {
ok=false;
}
break;
case uns::Rho :
if (particles->rho.size()>0) {
*data = &particles->rho[0];
*n=particles->rho.size();
} else {
ok=false;
}
break;
case uns::U :
ok=false;
break;
case uns::Hsml :
if (particles->hsml.size()>0) {
*data = &particles->hsml[0];
*n=particles->hsml.size();
} else {
ok=false;
}
break;
case uns::Temp :
if (particles->temp.size()>0) {
*data = &particles->temp[0];
*n=particles->temp.size();
} else {
ok=false;
}
break;
case uns::Age :
if (particles->age.size()>0) {
*data = &particles->age[0];
*n=particles->age.size();
} else {
ok=false;
}
break;
default: ok=false;
}
if (ok && !*data &&
(CunsOut::s_mapStringValues[name]!=uns::Nbody &&
CunsOut::s_mapStringValues[name]!=uns::Nsel)) ok = false; // not ok because array is NULL
if (verbose) {
if (ok) {
std::cerr << "CSnapshotGadgetIn::getData name["<<name<<"]=" << CunsOut::s_mapStringValues[name] << "\n";
} else {
std::cerr << "**WARNING** CSnapshotGadgetIn::getData Value ["<<name<<"] does not exist...\n";
}
}
return ok;
}
// ============================================================================
......
# ============================================================================
# Copyright Jean-Charles LAMBERT - 2010
# 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
# ============================================================================
# CMakeListst.txt file for UNS library
# ============================================================================
cmake_minimum_required(VERSION 2.6)
# project name
project (uns_projects)
# set CMAKE Module PATH (NEMO)
SET(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_BINARY_DIR})
MESSAGE( STATUS "CMAKE_MODULE_PATH=" ${CMAKE_MODULE_PATH} )
# require NEMO
FIND_PACKAGE(NEMO REQUIRED)
# Sqlite3
FIND_PACKAGE(SQLITE3)
# detect MacOS
SET(OSX FALSE)
if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
SET(OSX TRUE)
endif (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
# require boost
find_package( Boost 1.36.0)
set(LIBBOOST boost_thread-mt.a boost_system-mt.a)
set(DNOBOOST "")
set(Boost_FOUND FALSE)
if(NOT Boost_FOUND)
MESSAGE (STATUS " Boost not found, uns_2dplot will run slowly.....")
set(DNOBOOST "-DNOBOOST")
set(LIBBOOST "")
endif(NOT Boost_FOUND)
# Some necessary variable
set(DEHNEN $ENV{NEMO}/usr/dehnen)
set(NEMOLIB $ENV{NEMOLIB})
#FIND_PACKAGE(NEMO REQUIRED)
MESSAGE( STATUS "NEMOLIB : " ${NEMOLIB} )
# contains the full path to the top level directory of your build tree
MESSAGE( STATUS "PROJECT_BINARY_DIR: " ${PROJECT_BINARY_DIR} )
#MESSAGE ( STATUS "Boost_LIBRARY_DIRS : " ${Boost_LIBRARY_DIRS})
# contains the full path to the root of your project source directory,
# i.e. to the nearest directory where CMakeLists.txt contains the PROJECT() command
MESSAGE( STATUS "PROJECT_SOURCE_DIR: " ${PROJECT_SOURCE_DIR} )
# Find all the sources for the utils LIB
FILE(GLOB LIBUTILS ../lib/utils/*.cc)
# create the library "JCLutils"
add_library (JCLutils SHARED ${LIBUTILS})
if(OSX)
set_target_properties(JCLutils PROPERTIES LINK_FLAGS "-undefined suppress -flat_namespace")
endif(OSX)
# Find all the sources for the projects LIB
FILE(GLOB LIBPROJECTS ../lib/projects/*.cc)
# create the library "JCLutils"
add_library (JCLprojects SHARED ${LIBPROJECTS})
if(OSX)
set_target_properties(JCLprojects PROPERTIES LINK_FLAGS "-undefined suppress -flat_namespace")
endif(OSX)
# Destination path for the lib
SET(LIBRARY_OUTPUT_PATH ../lib)
# compilation options
add_definitions( -DNO_CUDA -O3 -DfalcON_NEMO -DfalcON_SINGLE -funroll-loops -Wall ${DNOBOOST})
# Executbale output path
set(EXECUTABLE_OUTPUT_PATH ../bin)
MESSAGE( STATUS "EXECUTABLE_OUTPUT_PATH: " ${EXECUTABLE_OUTPUT_PATH} )
# Make sure the compiler can find include files from our Hello library.
include_directories (../lib/utils ../lib/projects ${DEHNEN}/falcON/inc ${DEHNEN}/falcON/utils/inc $ENV{NEMOINC} ${NEMOLIB} $ENV{NEMOINC}/uns)
# Make sure the linker can find the Hello library once it is built.
link_directories (${NEMOLIB} ${DEHNEN}/falcON/lib ${DEHNEN}/utils/lib $ENV{PGPLOT_DIR} ${G2C_DIR} ${PROJECT_BINARY_DIR}/../lib /usr/lib64 /usr/X11/lib ${FC_GFORT_PATH} ${FC_G77_PATH}
)
# ----------------------------------------------------------
# Make sure the compiler can find include files from our Hello library.
#include_directories (${UNS_SOURCE_DIR}/src $ENV{NEMOINC} ${NEMOLIB})
# Find all LIB's sources files
FILE(GLOB execpp_sources ../src/*.cc ../src/*.c)
# build cpp executables according to the source
FOREACH(exe_cpp ${execpp_sources})
get_filename_component(exe ${exe_cpp} NAME_WE)
MESSAGE( STATUS "New executable ==> " ${exe})
add_executable (${exe} ${exe_cpp})
# Link the executable to the Hello library.
# target_link_libraries (${exe} JCLutils JCLprojects cpgplot pgplot X11 g2c unsio nemo sqlite3 ${Boost_LIBRARIES} boost_thread-mt.a boost_system-mt.a)
target_link_libraries (${exe} JCLutils JCLprojects cpgplot pgplot.a ${FC_GFORT_LIB} ${FC_G77_LIB} X11 unsio WDutils falcON nemo ${SQLITE3_LIB} gomp ${LIBBOOST} pthread dl)
ENDFOREACH(exe_cpp ${execpp_sources})
# ----------------------------------------------------------
# install target
set (CMAKE_INSTALL_PREFIX $ENV{NEMO})
IF (OSX) # Mac OSX
SET(SOEXT "dylib")
ELSE (OSX) # Linux
SET(SOEXT "so")
ENDIF(OSX)
#
# install bin targets
FOREACH(exe_cpp ${execpp_sources})
get_filename_component(exe ${exe_cpp} NAME_WE)
INSTALL(PROGRAMS ${PROJECT_BINARY_DIR}/../bin/${exe} DESTINATION bin)
ENDFOREACH(exe_cpp ${execpp_sources})
INSTALL(FILES ${PROJECT_BINARY_DIR}/../lib/libJCLutils.${SOEXT} DESTINATION lib)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../lib/libJCLprojects.${SOEXT} DESTINATION lib)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../man/uns_2dplot.1 DESTINATION man/man1)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../man/uns_stack.1 DESTINATION man/man1)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../man/uns2uns.1 DESTINATION man/man1)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../man/uns_projects.1 DESTINATION man/man1)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../lib/utils/cfalcon.h DESTINATION inc/uns)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../lib/utils/csnaptools.h DESTINATION inc/uns)
INSTALL(FILES ${PROJECT_BINARY_DIR}/../lib/utils/ctimer.h DESTINATION inc/uns)
#
SET(NEMO_INSTALLED FALSE)
FILE(GLOB GLOB_TEMP_VAR $ENV{NEMO})
IF(GLOB_TEMP_VAR)
SET(NEMO_INSTALLED TRUE)
ENDIF(GLOB_TEMP_VAR)
if (NOT NEMO_INSTALLED)
MESSAGE(SEND_ERROR "NEMO environement not loaded, cannot continue....")
endif (NOT NEMO_INSTALLED)
SET(NEMO_INSTALLED FALSE)
FILE(GLOB GLOB_TEMP_VAR $ENV{NEMO})
IF(GLOB_TEMP_VAR)
SET(NEMO_INSTALLED TRUE)
# Detect fortran compiler from ${NEMOLIB}/makedefs
# files (variable FC = )
FILE(READ $ENV{NEMOLIB}/makedefs MAKEDEFS)
STRING(REGEX MATCH "FC = gfortran" GFORTRAN_IS_SET ${MAKEDEFS})
IF ("FC = gfortran" STREQUAL "${GFORTRAN_IS_SET}")
MESSAGE("GFortran compiler Detected......")
SET (FORT_COMP gfortran)
ELSE("FC = gfortran" STREQUAL "${GFORTRAN_IS_SET}")
# g77 stuffs
MESSAGE("Assuming G77 compiler")
SET (FORT_COMP g77)
ENDIF ("FC = gfortran" STREQUAL "${GFORTRAN_IS_SET}")
# gfortran stuffs
SET (FC_GFORT_LIB "")
SET (FC_GFORT_PATH "")
set (GFORTPATH "")
execute_process(COMMAND gfortran -print-file-name=libgfortran.a
OUTPUT_VARIABLE GFORTPATH OUTPUT_STRIP_TRAILING_WHITESPACE)
IF (NOT ${GFORTPATH} STREQUAL "")
MESSAGE("GFORTRAN compiler exist, adding libgfortran......")
execute_process(COMMAND dirname ${GFORTPATH}
OUTPUT_VARIABLE FC_GFORT_PATH OUTPUT_STRIP_TRAILING_WHITESPACE)
SET (FC_GFORT_LIB gfortran)
ENDIF(NOT ${GFORTPATH} STREQUAL "")
# g77 stuffs
SET (FC_G77_LIB "")
SET (FC_G77_PATH "")
set (G2CFPATH "")
execute_process(COMMAND g77 -print-file-name=libg2c.a
OUTPUT_VARIABLE G2CFPATH OUTPUT_STRIP_TRAILING_WHITESPACE)
IF (NOT ${G2CFPATH} STREQUAL "")
MESSAGE("G77 compiler exist, adding libg2c......")
execute_process(COMMAND dirname ${G2CFPATH}
OUTPUT_VARIABLE FC_G77_PATH OUTPUT_STRIP_TRAILING_WHITESPACE)
SET (FC_G77_LIB g2c)
ENDIF(NOT ${G2CFPATH} STREQUAL "")
ENDIF(GLOB_TEMP_VAR)
if (NOT NEMO_INSTALLED)
MESSAGE(STATUS "NEMO environement not loaded.... using nemo_light")
endif (NOT NEMO_INSTALLED)
# --------------------------------------------------------
# CMake module to detect SQLITE3 library
# --------------------------------------------------------
SET(SQLITE3_FOUND FALSE)
SET(DNOSQLITE3 "-DNOSQLITE3") # if sqlite3 does not exist
SET(SQLITE3_LIB_PATH "")
SET(SQLITE3_H_PATH "")
SET(SQLITE3_LIB "")
find_file(SQLITE3_H NAMES "sqlite3.h" PATH /usr/local/include)
IF (SQLITE3_H)
get_filename_component(SQLITE3_H_PATH ${SQLITE3_H} PATH)
MESSAGE(STATUS "Found sqlite3.h:" ${SQLITE3_H})
MESSAGE(STATUS "Found sqlite3.h path:" ${SQLITE3_H_PATH})
find_library(SQLITE3 NAMES sqlite3 PATH /usr/local/lib)
IF (SQLITE3)
SET(SQLITE3_FOUND TRUE)
SET(DNOSQLITE3 "") # SQLITE3 exist
SET(SQLITE3_LIB sqlite3)
MESSAGE(STATUS "Found library here :" ${SQLITE3})
get_filename_component(SQLITE3_LIB_PATH ${SQLITE3} PATH)
MESSAGE(STATUS "Found library PATH :" ${SQLITE3_LIB_PATH})
ENDIF(SQLITE3)
ENDIF(SQLITE3_H)
This diff is collapsed.
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2009-2011
// 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@oamp.fr>
*/
#ifndef C2DPLOT_H
#define C2DPLOT_H
#include <string>
#include <vector>
#include <iostream>
#ifndef NOBOOST
#include <boost/thread/thread.hpp>
#include <boost/thread/mutex.hpp>
#endif
#include "cgaussian.h"
#define NTHREAD_MAX 256
#ifndef NO_CUDA
#define NB_BLOCK 2
#endif
namespace uns_proj {
typedef struct {
int devid;
int nblock, offset;
} t_cuda_grid;
template <class T> class C2dplot {
public:
C2dplot(const int,const int,const int, const int, const T);
void compute(std::string pic, const int _no_frame,const int _nbody, T * _pos ,
float _range[3][2], std::string _title,
std::string _sel_comp, std::string _filename, float _time,
bool _xy, bool _xz, bool _zy, bool _sview);
private:
int nthreads, dimx, dimy, pixel;
jclut::CGaussian<float> * gaussian;
T g;
float * tab[NTHREAD_MAX];
std::vector <int> indexes;
// drawImage
void drawImage(const bool disp,const int xaxis, const int yaxis, const int nbview, int &showtext);
// computeImage
void computeImage(const int xaxis, const int yaxis);
// buildFrameName
std::string buildFrameName(std::string label,const int idx);
// display text
void displayText(bool sview);
void findIndexes(const int xaxis, const int yaxis);
void displayIndexes() {
for (unsigned int i=0; i<indexes.size(); i++) {
std::cerr << indexes[i] << "\n";
}
}
// Parallel tasks
void startWorkers(const int nbody, T * data, const int xaxis, const int yaxis,float& zmin,float& zmax);
void worker(const int ithread, const int offset, const int nbody, T * data, const int xaxis, const int yaxis);
#ifndef NOBOOST
std::vector<boost::shared_ptr<boost::thread> > v_threads;
boost::mutex io_mutex;
#endif
float xmin,xmax,ymin,ymax;
std::string dev, title, sel_comp, filename;
int no_frame, nbody;
float time,range[3][2];
T * pos;
bool xy, xz, zy, sview;
#ifndef NO_CUDA
// CUDA stuffs
int GPU_N;
t_cuda_grid * device_data;
float
* cuda_pos; // particles's positions
int
c_nbody; // #cuda bodies particles
void initCuda();
void startCuda(const int nbody, T * data, const int xaxis, const int yaxis,float& zmin,float& zmax);
void workerCudaThread(const int id, T * data, const int xaxis, const int yaxis);
#endif
};
}
#endif
//
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2011
// 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
// ============================================================================
#include <iostream>
#include <cmath>
#include <sstream>
#include "cbar.h"
#include "uns.h"
#include "csnaptools.h"
using namespace uns_proj;
using namespace std;
#define PI 3.141592653589793238462643
// ----------------------------------------------------------------------------
// contructor
CBar::CBar(const int _nbody, float * _pos, float * _vel, float * _mass, float *_rho, float *_hsml)
{
// link arguments
nbody = _nbody;
pos = _pos;
vel = _vel;
mass= _mass;
density = NULL;
rho = _rho;
hsml=_hsml;
// sort by rho
sortRho();
}
// ----------------------------------------------------------------------------
// contructor
CBar::~CBar()
{
if (density) delete density;
}
// ----------------------------------------------------------------------------
// computeAngle();
// Based on anglebar_impr.F
// with new method using density instead of mass
float CBar::computeAngle(const float dmin, const float dmax, const bool mvcod)
{
// shift to COD
double cod[6]={0.,0.,0.,0.,0.,0.};
if (mvcod)
CSnaptools::moveToCod<float>(nbody,&pos[0],&vel[0],&mass[0],rho,cod,false);
float alpha=0.0;
float beta =0.0;
float minrho=log(rho[vec_rho.at(0).index]);
float maxrho=log(rho[vec_rho.at(nbody-1).index]);
float binf=minrho+dmin*(maxrho-minrho);
float bsup=minrho+dmax*(maxrho-minrho);
std::cerr << "binf/bsup :"<< binf<<"/"<<bsup<<"\n";
int cpt=0;
//for (int i=dmin*nbody; i<dmax*nbody; i++) {
for (int i=0; i<nbody;i++) {
int ii= vec_rho.at(i).index;
float rhoii = log(rho[ii]);
if ((rhoii)>=binf && (rhoii)<=bsup ) {
cpt++;
float x = pos[ii*3+0]-cod[0];
float y = pos[ii*3+1]-cod[1];
float x2 = x*x;
float y2 = y*y;
float r2 = x2+y2;
float cos2th = (x2 - y2) / r2;
float sin2th = 2. * x * y / r2;
#if 1
alpha += rho[ii]*sin2th;
beta += rho[ii]*cos2th;
#else
alpha += mass[ii]*sin2th;
beta += mass[ii]*cos2th;
#endif
}
}
std::cerr << "Found ["<<cpt<<"] particles into the range.\n";
assert(cpt>0);
float bar_angle = 0.5 * atan2(alpha,beta);
return bar_angle;
}
// ----------------------------------------------------------------------------
// computeAngle();
// Based on anglebar_impr.F
// with new method using density instead of mass
// select max #particles per density shells
// guess rotation angle from selected particles
float CBar::computeAngle(const bool mvcod)
{
float minrho=log(rho[vec_rho.at(0).index]);
float maxrho=log(rho[vec_rho.at(nbody-1).index]);
// reset array
for (int i=0; i<100; i++) {
data_histo[i] = 0;
}
// compute #particles per density shell
for (int i=0; i<nbody;i++) {
int ii= vec_rho.at(i).index;
int index=((log(rho[ii])-minrho)*99.)/(maxrho-minrho);
//std::cerr << "CBar::computeAngle index="<<index<<"\n";
assert(index<100);
data_histo[index]++;
}
// find max shell
int maxshell=data_histo[0];
int ishell=0;
for (int i=1;i<100;i++) {
if (data_histo[i]>maxshell) {
maxshell = data_histo[i];
ishell=i;
}
}
float dmax=std::max(ishell+5,ishell);
float dmin=std::max(0.,(ishell-20.));
std::cerr << "CBar::computeAngle dmin="<<dmin<<"/ dmax="<<dmax<<"\n";
return computeAngle(dmin/100.,dmax/100.,mvcod);
}
// ----------------------------------------------------------------------------
// rotateOnX
void CBar::rotateOnX(const float bar_angle)
{
float theta = -bar_angle;
rotate(theta);
}
// ----------------------------------------------------------------------------
// rotateOnY
void CBar::rotateOnY(const float bar_angle)
{
float theta = 0.5*M_PI-bar_angle;
rotate(theta);
}
// ----------------------------------------------------------------------------
// rotate
void CBar::rotate(const float theta)
{
for (int i=0; i<nbody; i++) {
float rx = pos[i*3+0] * cos(theta) - pos[i*3+1] * sin(theta);
float ry = pos[i*3+0] * sin(theta) + pos[i*3+1] * cos(theta);
pos[i*3+0] = rx;
pos[i*3+1] = ry;
}
if (vel) { // velocities exist
for (int i=0; i<nbody; i++) {
float rx = vel[i*3+0] * cos(theta) - vel[i*3+1] * sin(theta);
float ry = vel[i*3+0] * sin(theta) + vel[i*3+1] * cos(theta);
vel[i*3+0] = rx;
vel[i*3+1] = ry;
}
}
}
// ----------------------------------------------------------------------------
// save
void CBar::save(std::string out, const float timu, const bool mvcod)
{
// shift to COD
double cod[6]={0.,0.,0.,0.,0.,0.};