Commit a915efa9 authored by LAMBERT Jean-charles's avatar LAMBERT Jean-charles

cleaning code

parent 5b180b4d
#!/usr/bin/env python
"""
This script is called by external python script through execfile() call
"""
import matplotlib
matplotlib.use('Agg') # for offscreen rendering
from simulations.cinert import *
from simulations.ccod import *
from simulations.ccom import *
from simulations.crectify import *
from simulations.cmovie import *
from simulations.c2dplot import *
run_cod=1
run_com=0
run_rectify=0
run_movie=0
run_movie2=0
run_inert=0
print ("Matplotlib backend Using:",matplotlib.get_backend(),file=sys.stderr)
##
## COD
##
if run_cod:
data.lock_id=0 # lock file ID for COD analysis
data.cod_dir=data.sim_info['dir']+"/ANALYSIS/xCODvCOM_2"
#data.cod_select="halo: disk: stars: gas: bndry: bulge: halo_1: halo_2: halo,disk: halo,stars: disk,stars: halo,disk,stars"
data.cod_select="halo_2"
if data.first:
cod=CCod("",analysis=data,threshold=-20,verbose_debug=True) # initialyze COD analysis
cod.smartAnalysis(analysis=data)
##
## COM
##
if run_com:
data.lock_id=1
data.com_select="halo: disk: stars: gas: disk,stars: halo,disk,stars"
if data.first:
com=CCom(analysis=data,verbose_debug=True) # initialyze COM analysis
com.smartAnalysis()
##
## RECTIFY
##
if run_rectify:
data.lock_id=2
data.rectify_dir=data.sim_info['dir']+"/ANALYSIS/rectify2"
data.rectify_select="stars#halo,disk,stars: stars#stars: gas#stars@10.0 " # syntax: "component#cod@rcut"
if data.first:
rectify=CRectify(analysis=data,verbose_debug=True) # initialyze RECTIFY analysis
rectify.smartAnalysis()
##
## MOVIE
##
if run_movie:
data.lock_id=3
#data.movie_dir="ANALYSIS/"+data.simname+"/movie"
#data.movie_select syntax: component # propertie # percen # newradius # extdir
data.movie_select="gas#none#99#1#whole : stars#none#99#1#whole : disk#none#99#1#whole:"+ \
"gas#none#-30#1#center : stars#none#-30#1#center : disk#none#-30#1#center:"+ \
"stars#age#99#1#age "
if data.first:
movie=CMovie(analysis=data,verbose_debug=True) # initialyze MOVIE analysis
movie.smartAnalysis()
##
## MOVIE2
##
if run_movie2:
data.lock_id=4
data.movie_dir=data.sim_info['dir']+"/ANALYSIS/movie2"
data.cod_dir=data.sim_info['dir']+"/ANALYSIS/xCODvCOM"
#data.movie_select syntax: component # propertie # percen # newradius # extdir
data.movie_select="gas,stars,halo"
mergers="@"+data.sim_info['name']
center_cod="@"+data.sim_info['name']
cmap="gas:jet,stars:Paired,halo:Accent"
if data.first:
movie=C2dplot(analysis=data,verbose_debug=True) # initialyze MOVIE analysis
movie.smartAnalysis(component=data.movie_select,xrange=20.,mergers=mergers,center_cod=center_cod,cmap=cmap,contour=False,noxz=False)
movie.smartAnalysis(component=data.movie_select,xrange=20.,mergers=mergers,center_cod=center_cod,cmap=cmap,contour=True,noxz=False)
##
## INERT
##
if run_inert:
data.lock_id=5 # lock file ID for INERT analysis
# common data
ncut=20
codf="@"+data.simname
comp="halo"
fract_min=0.1
fract_max=0.9
tree_density=True
data.cod_dir=data.sim_info['dir']+"/ANALYSIS/xCODvCOM"
#data.uns_snap.debugOn()
print("DATA SIM INFO dir [%s]"%(data.sim_info['dir']))
# 10% density threshold
tree_threshold=-10.0
data.inert_dir=data.sim_info['dir']+"/ANALYSIS/inert/inertia_10"
if data.first:
inert10=CInert(analysis=data,verbose_debug=True,ncut=ncut,component="halo",center_file=codf,
fract_min=fract_min, fract_max=fract_max, tree_threshold=tree_threshold,
tree_density=tree_density)
inert10.smartAnalysis(analysis=data)
# 20% density threshold
tree_threshold=-20.0
data.inert_dir=data.sim_info['dir']+"/ANALYSIS/inert/inertia_20"
if data.first:
inert20=CInert(analysis=data,verbose_debug=True,ncut=ncut,component="halo",center_file=codf,
fract_min=fract_min, fract_max=fract_max, tree_threshold=tree_threshold,
tree_density=tree_density)
inert20.smartAnalysis(analysis=data)
data.first=False #MANDATORY
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2019
// e-mail: Jean-Charles.Lambert@lam.fr
// address: CeSAM
// 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> // C++ I/O
#include <fstream> // C++ file I/O
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <assert.h>
#include <cmath>
namespace nos { // use namespace to avoid conflict between
// typedef char * string from stdinc.h
// and using std::string from ccfits.h
typedef char * string;
#include "getparam.h"
}
// --------------------------------------------------------------
#include "uns.h"
#include "csnaptools.h"
#include "ctree.h"
#include "cneibors.h"
#include "ctimer.h"
using namespace jclut;
using namespace jcltree;
// --------------------------------------------------------------
const char * defv[] = {
"in=???\n UNS input snapshot",
"intp=???\n UNS input test particles snapshot",
"out=???\n output snapshot",
"select=???\n select particles (range, or component name)",
"neib=6\n minimum #neibours",
"direct=f\n use direct method (very slow)",
//"smr=f\n stop max radius?",
"times=all\n selected time",
"verbose=f\n verbose on/off"
"VERSION=1.O\n compiled on <" __DATE__ "> JCL ",
NULL,
};
const char * usage="compute density on test particles";
const long long one=1;
// --------------------------------------------------------------
// main
int main(int argc, char ** argv )
{
if (argc) {;}
// start NEMO
nos::initparam(const_cast<char**>(argv),const_cast<char**>(defv));
/* recuperation des parametres en entree */
std::string simname(nos::getparam((char *) "in" ));
std::string simnametp(nos::getparam((char *) "intp" ));
std::string outname(nos::getparam((char *) "out" ));
std::string select_c(nos::getparam((char *) "select" ));
std::string select_t(nos::getparam((char *) "times" ));
int neib=nos::getiparam((char *) "neib" );
bool direct = (nos::getbparam((char *) "direct" ));
//bool smr = (nos::getbparam((char *) "smr" ));
bool verbose = (nos::getbparam((char *) "verbose"));
// instantiate a new uns object
//s::Cuns * uns = new uns::Cuns(simname,select_c,select_t);
float * postp=NULL, *masstp=NULL, *veltp=NULL;
int *idtp=NULL;
int nbodytest=0;
uns::CunsIn * unstp = new uns::CunsIn(simnametp,"all","all",verbose);
if (unstp->isValid()) {
if (unstp->snapshot->nextFrame("mxvI")) {
bool ok;
// get the input number of bodies for tests particles
ok=unstp->snapshot->getData("nsel",&nbodytest);
assert(ok);
ok = unstp->snapshot->getData("pos",&nbodytest,&postp);
assert(ok);
}
}
// instantiate a new uns object
//s::Cuns * uns = new uns::Cuns(simname,select_c,select_t);
uns::CunsIn * uns = new uns::CunsIn(simname,select_c,select_t,verbose);
if (uns->isValid()) {
while(uns->snapshot->nextFrame("mxvI")) {
int nbody;
float time;
// get the input number of bodies according to the selection
uns->snapshot->getData("nsel",&nbody);
// get the simulation time
uns->snapshot->getData("time",&time);
std::cerr << "nbody=" << nbody << " time="<< time <<"\n";
std::cerr << "filename = " << uns->snapshot->getFileName() << "\n";
bool ok;
float * pos=NULL, * mass=NULL, *vel=NULL;
int * id=NULL;
ok = uns->snapshot->getData("id",&nbody,&id);
ok = uns->snapshot->getData("mass",&nbody,&mass);
ok = uns->snapshot->getData("vel",&nbody,&vel);
ok = uns->snapshot->getData("pos",&nbody,&pos);
CTimer timing;
double time_maketree;
if (ok) {
CTree<float> * tree = new CTree<float>(nbody,pos,(float *) NULL);
time_maketree = timing.cpu();
std::cerr << "Tree done in :" << time_maketree << " secondes\n";
timing.restart();
// create neibors object
CNeibors<float> * neibors = new CNeibors<float>(tree);
std::vector<CDistanceId> tabneib;
float * rho = new float[nbodytest];
float * hsml= new float[nbodytest];
double rmax=tree->getRsize()/ double(one<<(tree->getLevelMax()/2));
double pi=acos(-1);
double vsphere=4*pi/3.0;
int count_search=0;
double radius=tree->getRsize()/ double(one<<(tree->getLevelMax()));
std::cerr << "First Radius = "<<radius<<"\n";
std::cerr << "Rmax = "<<rmax<<"\n";
if (direct) {
std::cerr << "Using direct method...";
} else {
std::cerr << "Using tree method...\n";
}
float onepercent=0.;
for (int i=0; i<nbodytest;i++) {
tabneib.clear();
//radius /= 2.0;
//double radius=tree->distanceBodyToMesh(i)*1.5;
assert(radius>=0.);
//neibors->setStopAtMaxRadius(smr);
//neibors->setMaxRadius((radius)*1.5);
float * pp=postp+(i*3);
if (!direct) {
neibors->process(pp,neib,&tabneib);
} else {
neibors->direct(pp,neib,&tabneib);
}
onepercent++;
if (onepercent>=0.1*nbodytest/100.) {
onepercent=0.;
fprintf(stderr,"\rcompleted : %.2f %%",i*100./nbodytest);
}
int nneib=1;
radius = tabneib[tabneib.size()-1].getDistance();
nneib= tabneib.size();
count_search++;
//std::cerr << "#neibs = "<<nneib<<"\n";
assert(radius>=0.);
assert(nneib>=neib);
//rho[i]=(Level(tree->getBodyData()+i))/radius/radius/radius;
rho[i]=nneib/(vsphere*radius*radius*radius);
hsml[i]=radius;
//std::cerr << "i="<< i<< " radius ="<<radius << " rho="<<rho<< " level=" << Level(tree->getBodyData()+i)<<"\n";
}
std::cerr << "\n";
std::cerr << "Search neighbours :"<<timing.cpu() << " secondes\n";
std::cerr << "Level max="<<tree->getLevelMax() << "\n";
std::cerr << "TREE Level min="<<tree->getLevelMin() << "\n";
uns::CunsOut * unsout = new uns::CunsOut(outname,"nemo",verbose);
unsout->snapshot->setData("time",time);
unsout->snapshot->setData("all","pos",nbodytest,postp,false);
if (masstp)
unsout->snapshot->setData("all","mass",nbodytest,masstp,false);
if (veltp)
unsout->snapshot->setData("all","vel",nbodytest,veltp,false);
if (rho)
unsout->snapshot->setData("all","rho",nbodytest,rho,false);
if (hsml)
unsout->snapshot->setData("all","hsml",nbodytest,hsml,false);
if (idtp)
unsout->snapshot->setData("all","id",nbodytest,idtp,false);
// save snapshot
unsout->snapshot->save();
delete tree;
delete neibors;
/*
if (pos) delete [] pos;
if (vel) delete [] vel;
if (mass) delete [] mass;
if (id) delete [] id;
*/
if (rho) delete [] rho;
if (hsml) delete [] hsml;
delete unsout;
}
}
}
// finish NEMO
nos::finiparam();
}
//
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