Commit e3aac64e authored by jclamber's avatar jclamber
Browse files

fits2uns first working implementation

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@118 f264a43e-d52d-4b82-913b-c2bd5215a18a
parent 488aed29
......@@ -33,15 +33,18 @@ using namespace CCfits;
const char * defv[] = { // use `::'string because of 'using namespace std'
"in=???\n tipsy input file ",
"out=???\n output file)",
"type=nemo\n nemo|gadget2",
"type=gadget2\n nemo|gadget2",
"zmin=-1\n #z plane min",
"zmax=-1\n #z plane max",
"verbose=f\n verbose mode ",
"VERSION=1.0\n compiled on <"__DATE__"> JCL ",
NULL
};
const char * usage="Simple converter fits to uns compatible output format";
int readImage(std::string in, std::string out, std::string type);
int readImage(std::string in, std::string out, std::string type,int zmin, int zmax);
int main(int argc, char ** argv )
{
......@@ -51,6 +54,8 @@ int main(int argc, char ** argv )
// Get input parameters
std::string simname (nos::getparam ((char *) "in" ));
std::string outname (nos::getparam ((char *) "out" ));
int zmin (nos::getiparam ((char *) "zmin" ));
int zmax (nos::getiparam ((char *) "zmax" ));
std::string typeout = (nos::getparam ((char *) "type" ));
bool verbose = (nos::getbparam((char *) "verbose" ));
......@@ -58,7 +63,7 @@ int main(int argc, char ** argv )
try {
if (!readImage(simname,outname,typeout)) std::cerr << " readImage() \n";
if (!readImage(simname,outname,typeout,zmin,zmax)) std::cerr << " readImage() \n";
}
catch (FitsException&)
// will catch all exceptions thrown by CCfits, including errors
......@@ -73,7 +78,7 @@ int main(int argc, char ** argv )
nos::finiparam();
}
int readImage(std::string in, std::string out, std::string type)
int readImage(std::string in, std::string out, std::string type, int zmin, int zmax)
{
// std::auto_ptr<FITS> pInfile(new FITS("/windows7/JDL/ngc6503.fits",Read,true));
std::auto_ptr<FITS> pInfile(new FITS(in,Read,false));
......@@ -83,7 +88,7 @@ int readImage(std::string in, std::string out, std::string type)
HDU * image;// = &pInfile->extension(1);
std::cerr << " Begin \n";
std::valarray<unsigned long> contents;
std::valarray<float> contents;
try {
image = &pInfile->extension(1);
......@@ -100,7 +105,6 @@ int readImage(std::string in, std::string out, std::string type)
((PHDU *) image)->read(contents);
}
std::cerr << " yo yo \n";
#if 1
......@@ -120,15 +124,52 @@ int readImage(std::string in, std::string out, std::string type)
// this doesn't print the data, just header info.
//std::cout << *image << std::endl;
long ax(image->axis(0));
long ay(image->axis(1));
long az;
if (image->axes()>=3) {
az=(image->axis(2));
} else {
az=0;
}
if (zmin==-1) zmin=0;
if (zmax==-1) zmax=az;
std::vector <float> pos,hsml,rho;
long nan=0;
std::cout << "value = " << contents.size() << "\n";
for (unsigned long i=0; i<contents.size(); i++) {
if (isnan(contents[i])) {
std::cerr << "Nan ["<<i<<"]\n";
if (std::isfinite(contents[i])) {
long z_i=int(i/(ax*ay));
assert(z_i<=az);
if (z_i>=zmin && z_i<=zmax) {
long nxy=i-z_i*(ax*ay);
long y_i=int(nxy/ax);
assert(y_i<ay);
long x_i=nxy-(y_i*ax);
assert(z_i<az);
pos.push_back(x_i*1.0);
pos.push_back(y_i*1.0);
pos.push_back(z_i*1.0);
rho.push_back(contents[i]);
hsml.push_back(0.8);
}
if (isinf(contents[i])) {
std::cerr << "Inf ["<<i<<"]\n";
//std::cerr << x_i << " " << y_i << " " << z_i << " " << contents[i] << "\n";
} else { // NAN or INF
nan++;
}
}
std::cerr << " #NAN values = "<< nan << "\n";
uns::CunsOut * unsout = new uns::CunsOut(out,type,false);
unsout->snapshot->setData("time",0.0);
unsout->snapshot->setData("gas","pos",pos.size()/3,&pos[0],false);
unsout->snapshot->setData("gas","rho",rho.size(),&rho[0],false);
unsout->snapshot->setData("gas","hsml",hsml.size(),&hsml[0],false);
// save snapshot
unsout->snapshot->save();
return 0;
long ax1(image->axis(0));
......
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