snapshotgadget.cc 81.6 KB
Newer Older
jcl's avatar
jcl committed
1
// ============================================================================
2
// Copyright Jean-Charles LAMBERT - 2008-2016
jclamber's avatar
jclamber committed
3 4 5
//           Centre de donneeS Astrophysiques de Marseille (CeSAM)              
// e-mail:   Jean-Charles.Lambert@lam.fr                                      
// address:  Aix Marseille Universite, CNRS, LAM 
jcl's avatar
jcl committed
6 7 8 9
//           Laboratoire d'Astrophysique de Marseille                          
//           Pole de l'Etoile, site de Chateau-Gombert                         
//           38, rue Frederic Joliot-Curie                                     
//           13388 Marseille cedex 13 France                                   
jclamber's avatar
jclamber committed
10
//           CNRS UMR 7326                                       
jcl's avatar
jcl committed
11 12 13
// ============================================================================

/* 
jclamber's avatar
jclamber committed
14
	@author Jean-Charles Lambert <Jean-Charles.Lambert@lam.fr>
jcl's avatar
jcl committed
15 16 17 18 19 20 21 22 23 24 25
 */


#include <sstream>
#include <cstdlib>
#include "snapshotgadget.h"
#include "uns.h"
#include "ctools.h"
#define DEBUG 0
#include "unsdebug.h"

26

jcl's avatar
jcl committed
27 28 29 30
namespace uns {

// ----------------------------------------------------------------------------
// READING constructor                                                                 
31 32
template <class T>
CSnapshotGadgetIn<T>::CSnapshotGadgetIn(const std::string _name,
jcl's avatar
jcl committed
33 34 35
                                  const std::string _comp, 
                                  const std::string _time,
				  const bool verb)
36
   :CSnapshotInterfaceIn<T>(_name, _comp, _time, verb)
jcl's avatar
jcl committed
37 38
{
  filename = _name;
39
  //this->valid=false;
jcl's avatar
jcl committed
40 41 42 43 44
  first_loc = true;
  status=false;
  is_open=false;
  is_read=false;
  swap   = false;
45

jcl's avatar
jcl committed
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  mass   = NULL;
  pos    = NULL;
  vel    = NULL;
  acc    = NULL;
  pot    = NULL;
  id     = NULL;
  age    = NULL;
  metal  = NULL;
  intenerg=NULL;
  temp   = NULL;
  rho    = NULL;
  hsml   = NULL;
  zs     = NULL;
  zsmt   = NULL;
  im     = NULL;
jclamber's avatar
jclamber committed
61
  ssl    = NULL;
62
  cm     = NULL;
jcl's avatar
jcl committed
63
  bits   = 0;
64
  this->load_bits = 0;
jcl's avatar
jcl committed
65 66 67 68 69 70 71 72
  tframe = 0.;
  redshift=0.;
  frecord_offset = 4;
  czs = 0;
  czsmt= 0;
  bytes_counter = 0;
  multiplefiles = 0;
  lonely_file   = true;
73
  ntotmasses = 0;
74
  this->verbose = verb;
jcl's avatar
jcl committed
75 76 77
  
  int fail = open(filename);
  if (!fail) {
78
    this->valid=true;
jcl's avatar
jcl committed
79 80
    std::ostringstream stm;
    stm << getVersion();
81 82 83
    this->interface_type = "Gadget" + stm.str();
    this->interface_index=1;
    this->file_structure = "component";
jcl's avatar
jcl committed
84
  }
85

jcl's avatar
jcl committed
86 87 88
}
// ----------------------------------------------------------------------------
// destructor                                                                 
89 90
template <class T>
CSnapshotGadgetIn<T>::~CSnapshotGadgetIn()
jcl's avatar
jcl committed
91
{
92
  if (this->valid) {
jcl's avatar
jcl committed
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
    if (mass)     delete [] mass;
    if (pos)      delete [] pos;
    if (vel)      delete [] vel;
    if (acc)      delete [] acc;
    if (pot)      delete [] pot;
    if (id)       delete [] id;
    if (age)      delete [] age;
    if (metal)    delete [] metal;
    if (intenerg) delete [] intenerg;
    if (temp)     delete [] temp;
    if (rho)      delete [] rho;  
    if (hsml)     delete [] hsml;  
    if (zs)       delete [] zs;
    if (zsmt)     delete [] zsmt;
    if (im)       delete [] im;
jclamber's avatar
jclamber committed
108
    if (ssl)      delete [] ssl;
109
    if (cm)       delete [] cm;
110 111

    typename std::map<std::string,std::vector<T>  >::const_iterator it;
112 113 114 115 116 117
    for (it =  data_vector.begin(); it != data_vector.end(); it++) {
          std::cerr << "~CSnapshotGadgetIn(): Delete key = " <<  it->first
                    << "  size=" << it->second.size() << "\n";
          data_vector[it->first].clear();

      }
jcl's avatar
jcl committed
118 119 120 121 122
  }
  crv.clear();
}
// ============================================================================
// getSnapshotRange                                                            
123 124
template <class T>
ComponentRangeVector * CSnapshotGadgetIn<T>::getSnapshotRange()
jcl's avatar
jcl committed
125 126
{
  //crv.clear();
127
  if (this->valid && crv.size()) {
jcl's avatar
jcl committed
128 129
    //crv = getCRV();
    //ComponentRange::list(&crv);
130 131 132 133
    if (this->first) {
      this->first       = false;
      this->crv_first   = crv;
      this->nbody_first = getNtotal();
jcl's avatar
jcl committed
134
      //std::cerr << "CSnapshotGadgetIn::getSnapshotRange() = " << nbody_first << "\n";
135
      this->time_first  = getTime();
jcl's avatar
jcl committed
136 137 138 139 140 141 142 143 144
    }
  }
  return &crv;
}
//// ============================================================================
//// nextFrame 
//int CSnapshotGadgetIn::nextFrame()
//{
//  user_select.setSelection(getSelectPart(),getSnapshotRange());
145
//  return nextFrame(user_select.getIndexesTab(),user_select.this->getNSel());
jcl's avatar
jcl committed
146 147 148 149
//}

// ============================================================================
// nextFrame 
150 151
template <class T>
int CSnapshotGadgetIn<T>::nextFrame(uns::UserSelection &user_select)
jcl's avatar
jcl committed
152 153
{
  int status=0;
154
  assert(this->valid==true);
jcl's avatar
jcl committed
155 156
  if (first_loc) {
    first_loc = false;
157
    if (this->checkRangeTime(getTime())) {
jcl's avatar
jcl committed
158 159 160 161 162 163 164 165 166 167 168 169 170
      read(user_select); // read Gadget
      status = 1;
    }
  }
  return status;
}

// ============================================================================
// open() :                                                                    
// open file and return :                                                      
// 0 : success                                                                 
// 1 : unable to open                                                          
// 2 : not a GADGET file                                                       
171 172
template <class T>
int CSnapshotGadgetIn<T>::open(const std::string myfile)
jcl's avatar
jcl committed
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
{
  int fail=0;
  in.clear();
  in.open(myfile.c_str(),std::ios::in | std::ios::binary);
  if ( ! in.is_open()) {
   in.close();
   in.clear(); // mandatory under win32 !
   // try to open a multiple gadget file
   file0 = myfile+".0";
   PRINT("GadgetIn::open In gadget open file0 : " << file0 << "\n";)
   in.open(file0.c_str(),std::ios::in | std::ios::binary);
   if (in.is_open()) {
//     assert(in.good());
     lonely_file=false;
     PRINT("It's a multiple gadget file.\n";)
   }
  }
  if ( ! in.is_open()) {
    fail = 1;               // unable to open
    PRINT("In gadget, failed to open...\n";)
  }
  else {
    is_open=true;
    if (guessVersion()) {
      fail = readHeader(0);    // try to read header
      if (!fail) status=true; // valid header
      else close();           // not valid header
    }
    else {
      fail=1;
      close();             // not a valid gadgetfile
    }
  }
  return fail;
}
// ============================================================================
209 210
template <class T>
int CSnapshotGadgetIn<T>::close()
jcl's avatar
jcl committed
211 212 213 214 215 216
{
  if (is_open) in.close();
  is_open = false;
  return 1;
}
// ============================================================================
217 218 219
template <class T>
template <class U> int CSnapshotGadgetIn<T>::readCompData(U ** data, const int * index2,
                                                          const int * npartOffset,const int dim, const int nsel)
jcl's avatar
jcl committed
220 221 222
{
  bytes_counter=0;
  int len1 = readFRecord();
223 224
  //std::cerr << " SizeofT="<<sizeof(T) << " sizeofU="<<sizeof(U)<<"\n";
  checkFileVsArray(len1,sizeof(U),npart_total_local*dim);
jcl's avatar
jcl committed
225 226
  // allocate memory for the array
  if (*data==NULL) { // the first time
227
    *data = new U[dim*nsel];
jcl's avatar
jcl committed
228
  }
229
  U * ptr = *data;
230
  int fskip; // file skip factor
231
  if (array_vs_file_size<=1) {
232 233
    // 0 keep
    // 1 X 2
234
    fskip = sizeof(U) << array_vs_file_size;
235
  } else {
236
    // 2 divide by 2
237
    fskip= sizeof(U) >> 1;
238
  }
jcl's avatar
jcl committed
239 240 241 242
  for(int k=0;k<6;k++) {   
    if (header.npart[k]>0) { // there are particles for the component
      int idx=index2[npartOffset[k]];                
      if (idx != -1) {
243
        readData((char *) &ptr[dim*idx], sizeof(U), dim*header.npart[k]);
jcl's avatar
jcl committed
244
      } else {
245
        skipData(fskip*dim*header.npart[k]);
jcl's avatar
jcl committed
246 247 248 249
      }
    }
  }  
  int len2 = readFRecord();
250 251
  if (len2==len1) ; // remove warning....
  assert(len2==len1 && in.good() && len1==bytes_counter);
jcl's avatar
jcl committed
252 253 254
  return 1;
}
// ============================================================================
255 256 257
template <class T>
template <class U>
int CSnapshotGadgetIn<T>::readGasStarsUnknownArray(U ** data, int * nguess,const int * compOffset)
jcl's avatar
jcl committed
258 259 260
{
  bytes_counter=0;
  int len1 = readFRecord();
261
  *nguess = len1/sizeof(U)/(header.npart[0]+header.npart[4]); // geuss #records
jcl's avatar
jcl committed
262
  //std::cerr << "NGuess constant = " << *nguess << "\n";
263
  checkFileVsArray(len1,sizeof(U),(*nguess)*(header.npart[0]+header.npart[4]));
jcl's avatar
jcl committed
264 265
  // allocate memory for the array
  if (*data==NULL) { // the first time
266
    *data = new U[(*nguess)*(header.npartTotal[0]+header.npartTotal[4])];
jcl's avatar
jcl committed
267
  }
268
  U * ptr = *data;
jcl's avatar
jcl committed
269 270 271
  // read gas DATA
  int idx=compOffset[0]*(*nguess); // set index to the correct offset
  assert((idx+(*nguess)*header.npart[0])<=(*nguess)*(header.npartTotal[0]+header.npartTotal[4]));
272
  readData((char *) &ptr[idx], sizeof(U),(*nguess)*header.npart[0]);
jcl's avatar
jcl committed
273 274 275
  // read stars DATA
  idx=header.npartTotal[0]*(*nguess)+compOffset[4]*(*nguess);
  assert((idx+(*nguess)*header.npart[4])<=(*nguess)*(header.npartTotal[0]+header.npartTotal[4]));
276
  readData((char *) &ptr[idx], sizeof(U),(*nguess)*header.npart[4]);
jcl's avatar
jcl committed
277
  int len2 = readFRecord();
278 279
  if (len2==len1) ; // remove warning....
  assert(in.good() && len2==len1 && len1==bytes_counter);
jcl's avatar
jcl committed
280 281 282
  return 1;
}
// ============================================================================
283 284 285
template <class T>
template <class U>
int CSnapshotGadgetIn<T>::readOneArray(U ** data, const int compid,const int * compOffset)
jcl's avatar
jcl committed
286 287 288
{
  bytes_counter=0;
  int len1 = readFRecord();
289
  checkFileVsArray(len1,sizeof(U),header.npart[compid]);
jcl's avatar
jcl committed
290 291
  // allocate memory for the array
  if (*data==NULL) { // the first time
292
    *data = new U[header.npartTotal[compid]];
jcl's avatar
jcl committed
293
  }
294
  U * ptr = *data;
jcl's avatar
jcl committed
295 296
  int idx=compOffset[0]; // set index to the correct offset
  assert((idx+header.npart[compid])<=header.npartTotal[compid]);
297
  readData((char *) &ptr[idx], sizeof(U),header.npart[compid] );
jcl's avatar
jcl committed
298
  int len2 = readFRecord();
299
  if (len2==len1) ; // remove warning....
jcl's avatar
jcl committed
300 301 302
  assert(in.good() && len1==len2 && len1==bytes_counter);
  return 1;
}
jclamber's avatar
jclamber committed
303
// ============================================================================
304 305 306
template <class T>
template <class U>
int CSnapshotGadgetIn<T>::readStreamBlock(const std::string req_blockname, std::vector<U> & data)
jclamber's avatar
jclamber committed
307 308 309 310 311 312 313
{
  if (is_read) { // file has already been read
    close();
  }
  int fail = open(filename);
  if (fail) assert(0);     // fail is true abort

314
  //assert(sizeof(U)==sizeof(float)); // works only for float NOW
jclamber's avatar
jclamber committed
315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339

  int last_len=0;
  // loop on all the files
  for (int i=0; i<header.num_files || (i==0 && header.num_files==0);i++) {
    std::string infile;
    if (header.num_files > 0 ) { // more than one file
      std::ostringstream stm;
      stm << "." << i;               // add ".XX" extension
      infile = filename + stm.str(); // new filename
      if (i>0) {
        close();                 // close previous file
        int fail=open(infile);   // open new file,read header
        if (fail) assert(0);     // fail is true abort
      }
    }
    else infile=filename;        // lonely file

    bool stop=false;
    bool ok=false;

    while(readBlockName() && !stop) {
      if (req_blockname == block_name) { // found requested block name
        bytes_counter=0;
        int len1 = readFRecord();

340
        if (data.size()==0) { //  first time
341
          data.resize(len1/sizeof(U)); // !! TO CHECK
jclamber's avatar
jclamber committed
342 343 344
        }
        else {
          if (i>0) { // more than one file
345
            data.resize((last_len+len1)/sizeof(U));
jclamber's avatar
jclamber committed
346 347 348 349 350
          } else {
            std::cerr << "WARNING, pointer not allocated...\n";
          }
        }

351 352
        U * ptr = &data[0];
        readData((char *) &ptr[last_len], sizeof(U), len1/sizeof(U) );
jclamber's avatar
jclamber committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372
        int len2 = readFRecord();
        if (len2==len1) ; // remove warning....
        assert(in.good() && len1==len2 && len1==bytes_counter);
        last_len += len1;
        stop=true; // we must stop
        ok=true;
      } else {
        if (!ok) {
          if (in.eof()) { // it's end of file, must stop loop
            stop=true;
          }
          else {
            skipBlock();
          }
        }

      }
    }
  }

373 374 375
//  if (last_len>0) { // found something
//    *nbody=last_len/sizeof(T);
//  }
jclamber's avatar
jclamber committed
376 377
  return 1;
}
jcl's avatar
jcl committed
378 379

// ============================================================================
380 381
template <class T>
int CSnapshotGadgetIn<T>::read(uns::UserSelection &user_select)
jcl's avatar
jcl committed
382 383 384 385 386 387 388 389 390 391
{
  const uns::t_indexes_tab * index=user_select.getIndexesTab();
  const int nsel=user_select.getNSel();
  
  int npartOffset[6], // store absolute component offset (between components)
      compOffset[6];  // store relative component offset
  if (! is_read ) {
    //checkCompBits(index,nsel);
    
    // check component bits`
392
    this->comp_bits=user_select.compBits();
jcl's avatar
jcl committed
393 394 395 396 397 398 399 400 401 402 403 404

    //ComponentRange::list(&crv);
    //ComponentRange::list(user_select.getCrvFromSelection());
    
    is_read=true;
    assert(nsel<=npartTotal);
    // compute offset in case of multiple gadget file
    npartOffset[0] = 0;
    compOffset[0]  = 0;
    for (int i=1;i<=5;i++) {
      npartOffset[i] = npartOffset[i-1]+header.npartTotal[i-1];
      compOffset[i] = 0;
405
      if (this->verbose) std::cerr
jcl's avatar
jcl committed
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437
          << "npartOffset["<<i<<"]="<<npartOffset[i]<<" npartOffset["<<i-1<<"]="
          <<npartOffset[i-1]<<" header.npartTotal["<<i-1<<"]="<<header.npartTotal[i-1]<<"\n";
    }
    //for (int i=0;i<6;i++)   std::cerr << "npartOffset["<<i<<"]="<<npartOffset[i] <<"\n";
    
    // allocate array to store indexes
    int * index2  = new int[npartTotal];
    for (int i=0; i<npartTotal; i++)
      index2[i] = -1; // index2 init

    for (int i=0, cpt=0; i<npartTotal; i++) {
      int idx=index[i].i;
      assert(idx<npartTotal);
      //std::cerr << i<< " " << index[i].i << " "<< index[i].idx << " " << nsel << " " <<npartTotal << "\n";
      if (idx != -1 ) {
        index2[idx] = cpt++;
      }
    }

    int z_offset  =0; // metalicity offset
    int age_offset=0; // age offset
    int im_offset=0; // age offset
    // loop on all the files
    for (int i=0; i<header.num_files || (i==0 && header.num_files==0);i++) {
      std::string infile;
      if (header.num_files > 0 ) { // more than one file
        std::ostringstream stm;
        stm << "." << i;               // add ".XX" extension
        infile = filename + stm.str(); // new filename
        if (i>0) {
          close();                 // close previous file
          int fail=open(infile);   // open new file,read header
438
          if (fail) assert(0);     // fail is true abort
jcl's avatar
jcl committed
439 440 441 442 443 444 445 446 447 448
        }
      }
      else infile=filename;        // lonely file
      int len1=0,len2=0;
      bool stop=false;
      std::string next_block_name;
      if (version==1) next_block_name="POS";

      // Read the whole file till a valid block exist
      // or stop = true in case of Gadget1 file
449
      while (this->req_bits!=0 && readBlockName() && !stop) {
jcl's avatar
jcl committed
450 451 452
        if (version==1) block_name=next_block_name;
        bool ok=false;
        // --> Postions block
453 454
        if (block_name=="POS" && this->req_bits&POS_BIT) {
          this->load_bits |= POS_BIT;
jcl's avatar
jcl committed
455 456 457 458 459
          ok=true;
          readCompData(&pos,index2,npartOffset,3, nsel);
          if (version==1) next_block_name="VEL";
        }
        // --> Velocities block
460 461
        if (block_name=="VEL" && this->req_bits&VEL_BIT) {
          this->load_bits |= VEL_BIT;
jcl's avatar
jcl committed
462 463 464 465 466
          ok=true;
          readCompData(&vel,index2,npartOffset,3,nsel);
          if (version==1) next_block_name="ID";
        }
        // --> IDs block
467 468
        if (block_name=="ID" && this->req_bits&ID_BIT) {
          this->load_bits |= ID_BIT;
jcl's avatar
jcl committed
469 470 471 472 473
          ok=true;
          readCompData(&id,index2,npartOffset,1,nsel);
          if (version==1) next_block_name="MASS";
        }
        // --> MASS block
474 475
        if (block_name=="MASS" && this->req_bits&MASS_BIT ) {
          this->load_bits |= MASS_BIT;
jcl's avatar
jcl committed
476 477
          ok=true;
          bytes_counter=0;
478
          if (ntotmasses>0) {    // different masses
jcl's avatar
jcl committed
479
            len1 = readFRecord(); // we must read from disk
480
            checkFileVsArray(len1,sizeof(T),ntotmasses);
jcl's avatar
jcl committed
481 482
          }
          // allocate memory if NULL pointer
483
          if (! mass )  mass = new T[nsel  ];
jcl's avatar
jcl committed
484 485 486 487 488
          for(int k=0;k<6;k++) {
            for(int n=0;n<header.npart[k];n++){
              int idx=index2[npartOffset[k]+n];
              assert(idx<nsel);
              if (idx != -1 && (header.mass[k] == 0)) {          // variable mass selected
489
                readData((char *) &mass[idx], sizeof(T), 1); // read from disk
jcl's avatar
jcl committed
490 491
              } else {
                if (idx == -1 && (header.mass[k] == 0)) {        // variable mass **not** selected
492 493
                  T tmp;
                  readData((char *) &tmp, sizeof(T), 1);     // read from disk (skip it)
jcl's avatar
jcl committed
494 495 496 497 498 499 500 501
                } else {
                  if (idx != -1 && (header.mass[k] != 0)) {      // constant mass selected
                    mass[idx] = header.mass[k];                  // read from header
                  }
                }
              }
            } // for n
          } // for k
502
          if (ntotmasses > 0) {  // different masses
jcl's avatar
jcl committed
503
            len2 = readFRecord(); // we must read from disk
504 505
            if (len2==len1) ; // remove warning....
            assert(in.good() && len2==len1 && len1==bytes_counter);
jcl's avatar
jcl committed
506 507 508 509
          }
          if (version==1) stop=true; // we stop reading for gadget1
        }
        // --> POT block
510 511
        if (block_name=="POT" && this->req_bits&POT_BIT) {
          this->load_bits |= POT_BIT;
jcl's avatar
jcl committed
512 513 514 515
          ok=true;
          readCompData(&pot,index2,npartOffset,1,nsel);
        }
        // --> Acceleration block
516 517
        if (block_name=="ACCE" && this->req_bits&ACC_BIT) {
          this->load_bits |= ACC_BIT;
jcl's avatar
jcl committed
518 519 520 521
          ok=true;
          readCompData(&acc,index2,npartOffset,3,nsel);
        }
        // --> U block (Internal energy)
522 523
        if (block_name=="U" && this->req_bits&U_BIT && this->comp_bits&GAS_BIT) {
          this->load_bits |= U_BIT;
jcl's avatar
jcl committed
524 525 526 527 528
          assert(header.npart[0]>0); // (gas only)
          ok=true;
          readOneArray(&intenerg,0,compOffset);
        }
        // --> Temperature block
529 530
        if (block_name=="NE" && this->req_bits&TEMP_BIT && this->comp_bits&GAS_BIT) {
          this->load_bits |= TEMP_BIT;
jcl's avatar
jcl committed
531 532 533 534 535
          assert(header.npart[0]>0); // (gas only)
          ok=true;
          readOneArray(&temp,0,compOffset);
        }
        // --> RHO block (density)
536 537
        if (block_name=="RHO" && this->req_bits&RHO_BIT && this->comp_bits&GAS_BIT) {
          this->load_bits |= RHO_BIT;
jcl's avatar
jcl committed
538 539 540 541 542
          assert(header.npart[0]>0); // (gas only)
          ok=true;
          readOneArray(&rho,0,compOffset);
        }
        // --> HSML block (neighbours size)
543 544
        if (block_name=="HSML" && this->req_bits&HSML_BIT && this->comp_bits&GAS_BIT) {
          this->load_bits |= HSML_BIT;
jcl's avatar
jcl committed
545 546 547 548 549
          assert(header.npart[0]>0); // (gas only)
          ok=true;
          readOneArray(&hsml,0,compOffset);
        }
        // --> Z block (Metalicity)
550 551
        if (block_name=="Z" && this->req_bits&METAL_BIT && (this->comp_bits&GAS_BIT || this->comp_bits&STARS_BIT)) {
          this->load_bits |= METAL_BIT;
jcl's avatar
jcl committed
552 553 554 555
          assert((header.npart[0]+header.npart[4])>0); // gas+stars
          ok=true;
          bytes_counter=0;
          len1 = readFRecord();
556
          checkFileVsArray(len1,sizeof(T),header.npart[0]+header.npart[4]);
jcl's avatar
jcl committed
557 558
          // allocate memory for array only if NULL pointer
          if (!metal) {
559
            metal = new T[header.npartTotal[0]+header.npartTotal[4]];
jcl's avatar
jcl committed
560 561 562 563
          }
          // read gas metal
          int idx=compOffset[0];
          assert((idx+header.npart[0])<=(header.npartTotal[0]+header.npartTotal[4]));
564
          readData((char *) &metal[idx], sizeof(T),header.npart[0]);
jcl's avatar
jcl committed
565 566 567
          // read stars metal
          idx=header.npartTotal[0]+compOffset[4];
          assert((idx+header.npart[4])<=(header.npartTotal[0]+header.npartTotal[4]));
568
          readData((char *) &metal[idx], sizeof(T),header.npart[4]);
jcl's avatar
jcl committed
569 570 571 572
          len2 = readFRecord();
          assert(in.good() && len1==len2 && len1==bytes_counter);
        }
        // --> AGE block
573 574
        if (block_name=="AGE" && this->req_bits&AGE_BIT && this->comp_bits&STARS_BIT) {
          this->load_bits |= AGE_BIT;
jcl's avatar
jcl committed
575 576 577
          ok=true;
          bytes_counter=0;
          len1 = readFRecord();
578
          checkFileVsArray(len1,sizeof(T),header.npart[4]);
jcl's avatar
jcl committed
579 580
          // allocate memory for array only if NULL pointer
          if (!age) {
581
            age = new T[header.npartTotal[4]];
jcl's avatar
jcl committed
582
          }
583
          if ((int)(len1/sizeRealOnFile()) != header.npart[4]) {
jcl's avatar
jcl committed
584 585 586 587 588 589
            std::cerr << "\nWARNING: Wang's AGE bug detected.......skipping age\n";
            in.seekg(len1,std::ios::cur);
          } else {
            assert(header.npart[4]>0); // stars only
            int idx=age_offset;
            assert((idx+header.npart[4])<=(header.npartTotal[4]));
590
            readData((char *) &age[idx], sizeof(T),header.npart[4]);
jcl's avatar
jcl committed
591 592 593 594 595
          }
          len2 = readFRecord();
          assert(in.good() && len1==len2);// && len1==bytes_counter);
        }
        // --> IM block
596 597
        if (block_name=="iM" && this->req_bits&IM_BIT && this->comp_bits&STARS_BIT) {
          this->load_bits |= IM_BIT;
jcl's avatar
jcl committed
598 599 600 601
          assert((header.npart[4])>0); // stars
          ok=true;
          readOneArray(&im,4,compOffset);
        }
jclamber's avatar
jclamber committed
602
        // --> SSL block
603 604
        if (block_name=="SSL" && this->req_bits&SSL_BIT && this->comp_bits&STARS_BIT) {
          this->load_bits |= SSL_BIT;
jclamber's avatar
jclamber committed
605 606 607 608
          assert((header.npart[4])>0); // stars
          ok=true;
          readOneArray(&ssl,4,compOffset);
        }
609
        // --> CM block
610 611 612
        if (block_name=="cM" && this->req_bits&CM_BIT && (this->comp_bits&GAS_BIT || this->comp_bits&STARS_BIT)) {
          assert((this->load_bits&CM_BIT)==0); // if failed means that multiple file not supported for this block
          this->load_bits |= CM_BIT;
613 614 615 616 617 618
          assert((header.npart[0]+header.npart[4])>0); // gas+stars
          ok=true;
          int ccm;
          readGasStarsUnknownArray(&cm,&ccm,compOffset);
          assert(ccm==1); // there is only one value per particle
        }
jcl's avatar
jcl committed
619
        // --> Zs block
620 621 622
        if (block_name=="Zs" && this->req_bits&ZS_BIT && (this->comp_bits&GAS_BIT || this->comp_bits&STARS_BIT)) {
          assert((this->load_bits&ZS_BIT)==0); // if failed means that multiple file not supported for this block
          this->load_bits |= ZS_BIT;
jcl's avatar
jcl committed
623 624 625 626 627
          assert((header.npart[0]+header.npart[4])>0); // gas+stars
          ok=true;
          readGasStarsUnknownArray(&zs,&czs,compOffset);
        }
        // --> ZSMT block
628 629 630
        if (block_name=="ZSMT" && this->req_bits&ZSMT_BIT && (this->comp_bits&GAS_BIT || this->comp_bits&STARS_BIT)) {
          assert((this->load_bits&ZSMT_BIT)==0); // if failed means that multiple file not supported for this block
          this->load_bits |= ZSMT_BIT;
jcl's avatar
jcl committed
631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650
          assert((header.npart[0]+header.npart[4])>0); // gas+stars
          ok=true;
          readGasStarsUnknownArray(&zsmt,&czsmt,compOffset);
        }
        if (!ok) {
          if (in.eof()) {
            stop=true;
          }
          else {
            skipBlock();
            if (version==1) { // we set the nextblockname
              if (block_name=="POS") next_block_name="VEL";
              if (block_name=="VEL") next_block_name="ID";
              if (block_name=="ID" ) next_block_name="MASS";
            }
          }
        }
      } // end of while readBlock
      
      // add masses if no BLOCK_MASS present (mass inside the header)
651 652
      if (ntotmasses == 0 && this->req_bits&MASS_BIT) { // no BLOCK_MASS present (mass inside the header)
        this->load_bits |= MASS_BIT;
jcl's avatar
jcl committed
653
        // allocate memory if NULL pointer
654
        if (! mass )  mass = new T[nsel  ];
jcl's avatar
jcl committed
655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684
        for(int k=0;k<6;k++) {
          for(int n=0;n<header.npart[k];n++){
            int idx=index2[npartOffset[k]+n];
            assert(idx<nsel);
            if (idx != -1 && (header.mass[k] != 0)) {      // constant mass selected
              mass[idx] = header.mass[k];                  // read from header
            }
          } // for n
        } // for k
      } // if

      // correct the offset of the particles which have been read
      for (int i=0;i<6;i++) {
        npartOffset[i] = npartOffset[i]+header.npart[i];
        compOffset[i] = compOffset[i]+header.npart[i];
      }
      // correction for metalicity and age
      z_offset   += (header.npart[0]+header.npart[4]);
      age_offset += header.npart[4];
      im_offset  += header.npart[4];
    } // end of loop on numfiles

    // garbage collecting
    delete [] index2;

    // convert to temperature units
    if (header.npartTotal[0] > 0) {
      //unitConversion();
    }
    // check bits
685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700
    this->freeNotLoadedData(&mass     ,MASS_BIT);
    this->freeNotLoadedData(&pos      ,POS_BIT);
    this->freeNotLoadedData(&vel      ,VEL_BIT);
    this->freeNotLoadedData(&acc      ,ACC_BIT);
    this->freeNotLoadedData(&pot      ,POT_BIT);
    this->freeNotLoadedData(&age      ,AGE_BIT);
    this->freeNotLoadedData(&metal    ,METAL_BIT);
    this->freeNotLoadedData(&temp     ,TEMP_BIT);
    this->freeNotLoadedData(&intenerg ,U_BIT);
    this->freeNotLoadedData(&rho      ,RHO_BIT);
    this->freeNotLoadedData(&hsml     ,HSML_BIT);
    this->freeNotLoadedData(&zs       ,ZS_BIT);
    this->freeNotLoadedData(&zsmt     ,ZSMT_BIT);
    this->freeNotLoadedData(&im       ,IM_BIT);
    this->freeNotLoadedData(&ssl      ,SSL_BIT);
    this->freeNotLoadedData(&cm       ,CM_BIT);
jcl's avatar
jcl committed
701 702 703 704 705
  }
  return 1;
}
// ============================================================================
// unitConversion()                                                            
706 707
template <class T>
void CSnapshotGadgetIn<T>::unitConversion()
jcl's avatar
jcl committed
708
{
709
  double BOLTZMANN, PROTONMASS;
jcl's avatar
jcl committed
710
  double UnitLength_in_cm, UnitMass_in_g, UnitVelocity_in_cm_per_s;
711 712
  double UnitTime_in_s, UnitDensity_in_cgs, UnitEnergy_in_cgs;
  double Xh;
jcl's avatar
jcl committed
713

714 715 716 717
#if 0
  double UnitPressure_in_cgs, HubbleParam=0.65m G;
  double GRAVITY;
#endif
jcl's avatar
jcl committed
718 719 720
  double MeanWeight, u, gamma;
  double RhoUniverse_omegabar;
  /* physical constants in cgs units */
721
  //GRAVITY   = 6.672e-8;
jcl's avatar
jcl committed
722 723 724 725 726 727 728 729 730 731 732 733 734 735
  BOLTZMANN = 1.3806e-16;
  PROTONMASS = 1.6726e-24;

  /* internal unit system of the code */
  UnitLength_in_cm= 3.085678e21;   /*  code length unit in cm/h */
  UnitMass_in_g= 1.989e43;         /*  code mass unit in g/h */
#ifdef GRAPE_UNIT 
  UnitVelocity_in_cm_per_s= 207.9*1.0e5;
#else
  UnitVelocity_in_cm_per_s= 1.0e5;
#endif

  UnitTime_in_s= UnitLength_in_cm / UnitVelocity_in_cm_per_s;
  UnitDensity_in_cgs= UnitMass_in_g/ pow(UnitLength_in_cm,3);
736
  //UnitPressure_in_cgs= UnitMass_in_g/ UnitLength_in_cm/ pow(UnitTime_in_s,2);
jcl's avatar
jcl committed
737 738
  UnitEnergy_in_cgs= UnitMass_in_g * pow(UnitLength_in_cm,2) / pow(UnitTime_in_s,2);

739
  //G=GRAVITY/ pow(UnitLength_in_cm,3) * UnitMass_in_g * pow(UnitTime_in_s,2);
jcl's avatar
jcl committed
740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786

  Xh= 0.76;  /* mass fraction of hydrogen */
  
  RhoUniverse_omegabar=1.9e-29*0.04;

  assert(intenerg != NULL);

  for(int i=0;i<header.npart[0];i++) {

    MeanWeight= 4.0/(3*Xh+1+4*Xh*temp[i]) * PROTONMASS;

    /* convert internal energy to cgs units */

    u  = intenerg[i] * UnitEnergy_in_cgs/ UnitMass_in_g;

    gamma= 5.0/3;

    /* get temperature in Kelvin */

    temp[i] = MeanWeight/BOLTZMANN * (gamma-1) * u;
    if (rho) {
      rho[i] *= UnitDensity_in_cgs/RhoUniverse_omegabar;
    }
	  
#if 0
#ifndef NOAGE
    if(P[i].Type==4)
      {
	if(P[i].Age!=0)
	  P[i].Age=1.0/P[i].Age-1;
	else
	  {
	    //cout<<"Error Age"<<endl;
	    //cout<<"N="<<i<<P[i].Age<< endl;
	    //exit(0);
	  }

      }
#endif
#endif
  }

}
// ============================================================================
// readHeader():                                                               
// read gadget header structure                                                
// return values : 0 success, >0 not a Gadget header data                      
787 788
template <class T>
int CSnapshotGadgetIn<T>::readHeader(const int id)
jcl's avatar
jcl committed
789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810
{
  int len1,len2;

  // Header block
  readBlockName();
  bytes_counter=0;
  len1 = readFRecord();
  readData((char *)  header.npart         , sizeof(int   ),  6);
  readData((char *)  header.mass          , sizeof(double),  6);
  readData((char *) &header.time          , sizeof(double),  1);
  readData((char *) &header.redshift      , sizeof(double),  1);
  readData((char *) &header.flag_sfr      , sizeof(int   ),  1);
  readData((char *) &header.flag_feedback , sizeof(int   ),  1);
  readData((char *)  header.npartTotal    , sizeof(int   ),  6);
  readData((char *) &header.flag_cooling  , sizeof(int   ),  1);
  readData((char *) &header.num_files     , sizeof(int   ),  1);
  readData((char *) &header.BoxSize       , sizeof(double),  1);
  readData((char *) &header.Omega0        , sizeof(double),  1);
  readData((char *) &header.OmegaLambda   , sizeof(double),  1);
  readData((char *) &header.HubbleParam   , sizeof(double),  1);
  readData((char *)  header.fill          , sizeof(char  ), 96);
  len2 = readFRecord();
811
  if (this->verbose) {
jcl's avatar
jcl committed
812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831
    std::cerr << "header.flag_cooling = " << header.flag_cooling << "\n";
  }
  //std::cerr << "len1="<<len1<< " len2="<<len2<< " bytes_counter="<<bytes_counter<<"\n";
  if (in.bad() || len1!=len2 || len1!=bytes_counter)
    return 2;
  if (id==0) {                         // first time
    tframe = header.time;
    redshift = header.redshift;
    npartTotal = 0;
    npart_total_local = 0;
    ntotmasses  = 0 ; // !!!!!!!!!!! NEW
    for(int k=0; k<6; k++)  {
      npartTotal += header.npartTotal[k];    // count global total  particles
      npart_total_local += header.npart[k];  // count local total particles
    }
    for(int k=0; k<6; k++) {
      if (header.mass[k] == 0 ) {
        //ntotmasses += header.npartTotal[k];
        ntotmasses += header.npart[k];
      }
832
      if (this->verbose) {
jcl's avatar
jcl committed
833 834 835 836 837
        std::cerr << "mass["<<k<<"]="<<header.mass[k]<<"\n";
      }      
    }
    storeComponents();
  }
jclamber's avatar
jclamber committed
838 839
  in_start_block = in.tellg(); // save poistion of input file

jcl's avatar
jcl committed
840 841 842 843
  return 0;
}
// ============================================================================
// readBlockName : read Gadget2 file format block                              
844 845
template <class T>
bool CSnapshotGadgetIn<T>::readBlockName()
jcl's avatar
jcl committed
846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865
{
  bool status=true;
  if (version == 2 ) { // gadget2 file format
    int dummy,nextblock;
    char name[5];
    array_vs_file_size = 0;
    readData((char *) &dummy    , sizeof(int) , 1); // read
    readData((char *) name      , sizeof(char), 4); // read label
    readData((char *) &nextblock, sizeof(int) , 1); // read nextblock
    readData((char *) &dummy    , sizeof(int) , 1); // read
#if 0
    int i=0; while (isupper(name[i])&& i<4) i++;
    name[i]='\0';
#else
    int i=0; while (name[i]!=' ' && i<4) i++;
    name[i]='\0';
#endif

    block_name=name;
    status = in.good();
866
    if (status && block_name!= "HEAD" && this->verbose)
jcl's avatar
jcl committed
867 868 869 870 871 872 873 874
      std::cerr << "Reading Block Name : <" << block_name << ">\n";
  }
  return status;
}
// ============================================================================
// guessVersion()                                                              
// detect Gadget  file format version ( 1 or 2 )                               
// return true if successfully detected otherwise false                        
875 876
template <class T>
bool CSnapshotGadgetIn<T>::guessVersion()
jcl's avatar
jcl committed
877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903
{
  bool status=true;
  // try to read 32bits length fortran record
  int dummy;
  swap = false; // no swapping
  array_vs_file_size = 0;
  readData((char *) &dummy, sizeof(int), 1); // read
  //std::cerr << "GadgetIn::guessVersion dummy = "<<dummy<<"\n";
  if( dummy != 256  && dummy != 8  ) {           // unknow number
    swap = true;                                 // need bytes swapping
    swapBytes((int *) &dummy, sizeof(int));      // swap
    //std::cerr << "GadgetIn::guessVersion dummy swapped= "<<dummy<<"\n";
    if( dummy != 256  && dummy != 8  )           // unknow swapped number
     status = false;                             // not a gadget file    
  }
  if (status) {
    if (dummy==256) version=1; // gadget1
    else            version=2; // gadget2
    in.seekg(0,std::ios::beg); // rewind
    //std::cerr << "gadget Version : " <<version << "\n"; 
  }

  return status;
}
// ============================================================================
// readData:                                                                     
// perform IO (Read,Write) operation on Data                                   
904 905
template <class T>
int CSnapshotGadgetIn<T>::readData(char * ptr,const size_t size_bytes,const int items)
jcl's avatar
jcl committed
906
{
907
  if (array_vs_file_size==0) { // no conversion, same floating format
jcl's avatar
jcl committed
908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923
    bytes_counter += (size_bytes*items);

    // get data from file
    in.read(ptr,size_bytes*items);
    //assert(in.good());
    if (! in.good()) return 0;
    
    // We SWAP data
    if (swap && (size_bytes != CHAR)) { // swapping requested
      for (int i=0; i<items; i++) {
        swapBytes(ptr,size_bytes);
        ptr += size_bytes;
      }
    }
  } else if (array_vs_file_size==1) {
    bytes_counter += (size_bytes*2*items);
924
    // data bigger on file (double) than array(float)
jcl's avatar
jcl committed
925 926 927 928 929 930 931 932 933 934 935
    // we must use a temporary variable
    double tmp;
    for (int i=0; i<items; i++) {
      in.read((char *) &tmp,sizeof(double)); // read one double from file
      if (swap && (size_bytes != CHAR)) { // swapping requested
        swapBytes((char *) &tmp,sizeof(double));
      }
      float tofloat=(float) tmp; // cast to float
      memcpy(ptr+sizeof(float)*i,(char*)&tofloat,sizeof(float)); // copy back to array

    }
936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960
  }else {
    assert(array_vs_file_size==2);
    long int size_bytes_file = (size_bytes>>1);
    // data bigger on array (double) than file (float)
    // we use the destination array to store all read data
    bytes_counter += (size_bytes_file*items);

    // get data from file
    in.read(ptr+(size_bytes_file*items),size_bytes_file*items);
    //assert(in.good());
    if (! in.good()) return 0;

    // convert float read to double array

    for (int i=0; i<items; i++) {
      float  * value = (float * ) (ptr+(size_bytes_file*items)+i*sizeof(float));
      if (swap && (size_bytes != CHAR)) { // swapping requested
        swapBytes((char *) value,sizeof(float));
      }
      double tmp=(double) *value;
      char * p=(char *) &tmp;
      for (unsigned int j=0; j<size_bytes; j++) {
        ptr[size_bytes*i+j]=p[j];
      }
    }
jcl's avatar
jcl committed
961 962 963 964 965 966 967

  }
  return 1;
}
// ============================================================================
// getData                                                               
// return requested array according 'name' selection                                               
968 969
template <class T>
bool CSnapshotGadgetIn<T>::getData(const std::string comp, std::string name, int *n,T **data)
jcl's avatar
jcl committed
970 971 972 973 974 975
{
  bool ok=true;
  *data=NULL;
  *n = 0;
  
  int nbody,first,last;
976 977 978
  bool status=false;

  if (comp!="STREAM") {
979
    status=this->getRangeSelect(comp.c_str(),&nbody,&first,&last,false); // find components ranges
980 981 982
    if (!status && comp=="all") { // retreive all particles selected by the user
      status=1;
      first=0;
983
      nbody=this->getNSel();
984
    }
jcl's avatar
jcl committed
985
  }
986
  switch(CunsOut2<T>::s_mapStringValues[name]) {
jcl's avatar
jcl committed
987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003
  case uns::Nbody :
    if (status) {
      *data = NULL;
      *n = nbody;
    } else {
      ok = false;
    }
    break;  
  case uns::Nsel   :
    if (status) {
      *n    = nbody;
    } else {
      ok=false;
    }  
  case uns::Pos   :
    if (status && getPos()) {
      *data = &getPos()[first*3];
1004
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1005 1006 1007 1008 1009 1010 1011
    } else {
      ok=false;
    }
    break;
  case uns::Vel  :
    if (status && getVel()) {
      *data = &getVel()[first*3];
1012
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1013 1014 1015 1016 1017 1018 1019
    } else {
      ok=false;
    }
    break;
  case uns::Acc  :
    if (status && getAcc()) {
      *data = &getAcc()[first*3];
1020
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1021 1022 1023 1024 1025 1026 1027
    } else {
      ok=false;
    }
    break;
  case uns::Pot  :
    if (status && getPot()) {
      *data = &getPot()[first];
1028
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1029 1030 1031 1032 1033 1034 1035
    } else {
      ok=false;
    }
    break;
  case uns::Mass  :
    if (status && getMass()) {
      *data = &getMass()[first];
1036
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066
    } else {
      ok=false;
    }
    break;
  case uns::Rho :
    if (status && comp=="gas" && getRho(*n)) {
      *data = getRho(*n);
    } else {
      ok=false;
    }
    break;
  case uns::U :
    if (status && comp=="gas" && getU(*n)) {
      *data = getU(*n);
    } else {
      ok=false;
    } 
    break;
  case uns::Hsml :
    if (status && comp=="gas" && getHsml(*n)) {
      *data = getHsml(*n);
    } else {
      ok=false;
    }  
    break;
  case uns::Temp :
    if (status && comp=="gas" && getTemp(*n)) {
      *data = getTemp(*n);
    } else {
      ok=false;
1067
    }
jcl's avatar
jcl committed
1068 1069 1070 1071 1072 1073 1074 1075 1076
    break;
  case uns::Age :
    if (status && comp=="stars" && getAge(*n)) {
      *data = getAge(*n);
    } else {
      ok=false;
    } 
    break;
  case uns::Metal :
1077
    if (status && comp=="gas" && this->ckloadBit(METAL_BIT)) {
jcl's avatar
jcl committed
1078 1079
      *data = getMetalGas(*n);
    } else
1080
      if (status && comp=="stars" && this->ckloadBit(METAL_BIT)) {
jcl's avatar
jcl committed
1081 1082 1083 1084 1085
      *data = getMetalStars(*n);
    } else {
      ok=false;
    }    
    break;
1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098
  case uns::Cm :
    if (status && comp=="gas") {
      *data = getCmGas(*n);
    } else
      if (status && comp=="stars") {
        *data = getCmStars(*n);
      } else
        if (status && comp=="all") {
          *data = getCm(*n);
        } else {
          ok=false;
        }
    break;
jcl's avatar
jcl committed
1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131
  case uns::Zs :
    if (status && comp=="gas") {
      *data = getZsGas(*n);
    } else
      if (status && comp=="stars") {
        *data = getZsStars(*n);
      } else
        if (status && comp=="all") {
          *data = getZs(*n);
        } else {
          ok=false;
        }
    break;
  case uns::ZSMT :
    if (status && comp=="gas") {
      *data = getZsmtGas(*n);
    } else
      if (status && comp=="stars") {
        *data = getZsmtStars(*n);
      } else
        if (status && comp=="all") {
          *data = getZsmt(*n);
        } else {
          ok=false;
        }
    break;
  case uns::Im :
    if (status && comp=="stars" && getIm(*n)) {
      *data = getIm(*n);
    } else {
      ok=false;
    }
    break;
jclamber's avatar
jclamber committed
1132 1133 1134 1135 1136 1137 1138
  case uns::Ssl :
    if (status && comp=="stars" && getSsl(*n)) {
      *data = getSsl(*n);
    } else {
      ok=false;
    }
    break;
jclamber's avatar
jclamber committed
1139 1140
  default: // unkown name
      if (comp=="STREAM") { // data stream reading
1141
        if (data_vector[name].size()==0) { // map key does not exist
1142
          std::vector<T> * p = &(data_vector[name]);
1143 1144 1145 1146 1147 1148
          int status=readStreamBlock(name,*p);
          if (status>0) {
            ok=true;
          }
        } else {
          ok = true;
jclamber's avatar
jclamber committed
1149
        }
1150 1151 1152 1153 1154
        if (ok) {
          *n = data_vector[name].size();
          *data = &data_vector[name][0];
        }

jclamber's avatar
jclamber committed
1155 1156 1157
      } else {
        ok=false;
      }
jcl's avatar
jcl committed
1158 1159
  }
  if (ok && !*data &&
1160 1161 1162
      (CunsOut2<T>::s_mapStringValues[name]!=uns::Nbody &&
       CunsOut2<T>::s_mapStringValues[name]!=uns::Nsel)) ok = false; // not ok because array is NULL
  if (this->verbose) {
jcl's avatar
jcl committed
1163
    if (ok) {
1164
      std::cerr << "CSnapshotGadgetIn::getData name["<<name<<"]=" << CunsOut2<T>::s_mapStringValues[name] << "\n";
jcl's avatar
jcl committed
1165 1166 1167 1168 1169 1170 1171 1172 1173
    } else {
      std::cerr << "**WARNING** CSnapshotGadgetIn::getData Value ["<<name<<"] for component <"<<comp<<"> does not exist...\n";
    }
  }
  return ok;
}
// ============================================================================
// getData                                                               
// return requested array according 'name' selection                                               
1174 1175
template <class T>
bool CSnapshotGadgetIn<T>::getData(const std::string name,int *n,T **data)
jcl's avatar
jcl committed
1176 1177 1178 1179 1180
{
  bool ok=true;
  *data=NULL;
  *n = 0;
  
1181
  switch(CunsOut2<T>::s_mapStringValues[name]) {
jcl's avatar
jcl committed
1182 1183
  case uns::Pos   :
    *data = getPos();
1184
    *n    = this->getNSel();
jcl's avatar
jcl committed
1185 1186 1187
    break;
  case uns::Vel  :
    *data = getVel();
1188
    *n    = this->getNSel();
jcl's avatar
jcl committed
1189 1190 1191
    break;
  case uns::Mass  :
    *data = getMass();
1192
    *n    = this->getNSel();
jcl's avatar
jcl committed
1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208
    break;
  case uns::Rho :
    *data = getRho(*n);
    break;
  case uns::U :
    *data = getU(*n);
    break;
  case uns::Hsml :
    *data = getHsml(*n);
    break;
  case uns::Temp :
    *data = getTemp(*n);
    break;
  case uns::Age :
    *data = getAge(*n);
    break;
jclamber's avatar
jclamber committed
1209 1210
  case uns::Acc :
    *data = getAcc();
1211
    *n    = this->getNSel();
jclamber's avatar
jclamber committed
1212
    break;
jcl's avatar
jcl committed
1213
  case uns::Metal :
1214
    if (this->comp_bits&STARS_BIT && this->comp_bits&GAS_BIT)      // gas and stars requested
jcl's avatar
jcl committed
1215 1216
      *data = getMetal(*n);
    else
1217
      if (this->comp_bits&STARS_BIT)                         // only stars requested
jcl's avatar
jcl committed
1218 1219
        *data = getMetalStars(*n);
      else
1220
        if (this->comp_bits&GAS_BIT)                         // only gas requested
jcl's avatar
jcl committed
1221 1222 1223
          *data = getMetalGas(*n);
    break;
  case uns::GasMetal :
1224
    if (this->ckloadBit(METAL_BIT))
jcl's avatar
jcl committed
1225 1226 1227
      *data = getMetalGas(*n);
    break;
  case uns::StarsMetal :
1228
    if (this->ckloadBit(METAL_BIT))
jcl's avatar
jcl committed
1229 1230 1231
      *data = getMetalStars(*n);
    break;
  case uns::Zs :
1232
    if (this->comp_bits&STARS_BIT && this->comp_bits&GAS_BIT)      // gas and stars requested
jcl's avatar
jcl committed
1233 1234
      *data = getZs(*n);
    else
1235
      if (this->comp_bits&STARS_BIT)                         // only stars requested
jcl's avatar
jcl committed
1236 1237
        *data = getZsStars(*n);
      else
1238
        if (this->comp_bits&GAS_BIT)                         // only gas requested
jcl's avatar
jcl committed
1239 1240
          *data = getZsGas(*n);
    break;
1241
  case uns::Cm :
1242
    if (this->comp_bits&STARS_BIT && this->comp_bits&GAS_BIT)      // gas and stars requested
1243 1244
      *data = getCm(*n);
    else
1245
      if (this->comp_bits&STARS_BIT)                         // only stars requested
1246 1247
        *data = getCmStars(*n);
      else
1248
        if (this->comp_bits&GAS_BIT)                         // only gas requested
1249 1250
          *data = getCmGas(*n);
    break;
jcl's avatar
jcl committed
1251
  case uns::ZSMT :
1252
    if (this->comp_bits&STARS_BIT && this->comp_bits&GAS_BIT)      // gas and stars requested
jcl's avatar
jcl committed
1253 1254
      *data = getZsmt(*n);
    else
1255
      if (this->comp_bits&STARS_BIT)                         // only stars requested
jcl's avatar
jcl committed
1256 1257
        *data = getZsmtStars(*n);
      else
1258
        if (this->comp_bits&GAS_BIT)                         // only gas requested
jcl's avatar
jcl committed
1259 1260 1261 1262 1263
          *data = getZsmtGas(*n);
    break;
  case uns::Im :
    *data = getIm(*n);
    break;
jclamber's avatar
jclamber committed
1264 1265 1266
  case uns::Ssl :
    *data = getSsl(*n);
    break;