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

13 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
namespace uns {

// ----------------------------------------------------------------------------
30
// READING constructor
31 32
template <class T>
CSnapshotGadgetIn<T>::CSnapshotGadgetIn(const std::string _name,
33
                                  const std::string _comp,
jcl's avatar
jcl committed
34
                                  const std::string _time,
35
                  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;
75

jcl's avatar
jcl committed
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
    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;
103 104
    if (rho)      delete [] rho;
    if (hsml)     delete [] hsml;
jcl's avatar
jcl committed
105 106 107
    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
  }
  crv.clear();
}
// ============================================================================
122
// 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
    }
  }
  return &crv;
}
//// ============================================================================
141
//// nextFrame
jcl's avatar
jcl committed
142 143 144
//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
      read(user_select); // read Gadget
      status = 1;
    }
  }
  return status;
}

// ============================================================================
166 167 168 169 170
// 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
  }
239
  for(int k=0;k<6;k++) {
jcl's avatar
jcl committed
240
    if (header.npart[k]>0) { // there are particles for the component
241
      int idx=index2[npartOffset[k]];
jcl's avatar
jcl committed
242
      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
  }
jcl's avatar
jcl committed
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
{
  const uns::t_indexes_tab * index=user_select.getIndexesTab();
  const int nsel=user_select.getNSel();
385

jcl's avatar
jcl committed
386 387 388 389
  int npartOffset[6], // store absolute component offset (between components)
      compOffset[6];  // store relative component offset
  if (! is_read ) {
    //checkCompBits(index,nsel);
390

jcl's avatar
jcl committed
391
    // check component bits`
392
    this->comp_bits=user_select.compBits();
jcl's avatar
jcl committed
393 394 395

    //ComponentRange::list(&crv);
    //ComponentRange::list(user_select.getCrvFromSelection());
396

jcl's avatar
jcl committed
397 398 399 400 401 402 403 404
    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
          << "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";
410

jcl's avatar
jcl committed
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
    // 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
        if ((block_name=="SSL" || block_name=="ssl") && this->req_bits&SSL_BIT && this->comp_bits&STARS_BIT) {
604
          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
          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
649

jcl's avatar
jcl committed
650
      // 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
  }
  return 1;
}
// ============================================================================
705
// 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
  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 */
728
#ifdef GRAPE_UNIT
jcl's avatar
jcl committed
729 730 731 732 733 734 735
  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

  Xh= 0.76;  /* mass fraction of hydrogen */
742

jcl's avatar
jcl committed
743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
  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;
    }
763

jcl's avatar
jcl committed
764 765 766 767
#if 0
#ifndef NOAGE
    if(P[i].Type==4)
      {
768 769 770 771 772 773 774 775
    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);
      }
jcl's avatar
jcl committed
776 777 778 779 780 781 782 783

      }
#endif
#endif
  }

}
// ============================================================================
784 785 786
// 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
    std::cerr << "header.flag_cooling = " << header.flag_cooling << "\n";
  }
814 815 816 817 818
  t_header.BoxSize = (T) header.BoxSize;
  t_header.redshift = (T) header.redshift;
  t_header.Omega0  = (T) header.Omega0;
  t_header.OmegaLambda  = (T) header.OmegaLambda;
  t_header.HubbleParam  = (T) header.HubbleParam;
jcl's avatar
jcl committed
819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836
  //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];
      }
837
      if (this->verbose) {
jcl's avatar
jcl committed
838
        std::cerr << "mass["<<k<<"]="<<header.mass[k]<<"\n";
839
      }
jcl's avatar
jcl committed
840 841 842
    }
    storeComponents();
  }
jclamber's avatar
jclamber committed
843 844
  in_start_block = in.tellg(); // save poistion of input file

jcl's avatar
jcl committed
845 846 847
  return 0;
}
// ============================================================================
848
// readBlockName : read Gadget2 file format block
849 850
template <class T>
bool CSnapshotGadgetIn<T>::readBlockName()
jcl's avatar
jcl committed
851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
{
  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();
871
    if (status && block_name!= "HEAD" && this->verbose)
jcl's avatar
jcl committed
872 873 874 875 876
      std::cerr << "Reading Block Name : <" << block_name << ">\n";
  }
  return status;
}
// ============================================================================
877 878 879
// guessVersion()
// detect Gadget  file format version ( 1 or 2 )
// return true if successfully detected otherwise false
880 881
template <class T>
bool CSnapshotGadgetIn<T>::guessVersion()
jcl's avatar
jcl committed
882 883 884 885 886 887 888 889 890 891 892 893 894
{
  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
895
     status = false;                             // not a gadget file
jcl's avatar
jcl committed
896 897 898 899 900
  }
  if (status) {
    if (dummy==256) version=1; // gadget1
    else            version=2; // gadget2
    in.seekg(0,std::ios::beg); // rewind
901
    //std::cerr << "gadget Version : " <<version << "\n";
jcl's avatar
jcl committed
902 903 904 905 906
  }

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

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

jcl's avatar
jcl committed
920 921 922 923 924 925 926 927 928
    // 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);
929
    // data bigger on file (double) than array(float)
jcl's avatar
jcl committed
930 931 932 933 934 935 936 937 938 939 940
    // 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

    }
941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965
  }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
966 967 968 969 970

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

jcl's avatar
jcl committed
980
  int nbody,first,last;
981 982
  bool status=false;

983
  if (comp!="STREAM" && comp!="HEADER") {
984
    status=this->getRangeSelect(comp.c_str(),&nbody,&first,&last,false); // find components ranges
985 986 987
    if (!status && comp=="all") { // retreive all particles selected by the user
      status=1;
      first=0;
988
      nbody=this->getNSel();
989
    }
jcl's avatar
jcl committed
990
  }
991

992
  switch(CunsOut2<T>::s_mapStringValues[name]) {
jcl's avatar
jcl committed
993 994 995 996 997 998 999
  case uns::Nbody :
    if (status) {
      *data = NULL;
      *n = nbody;
    } else {
      ok = false;
    }
1000
    break;
jcl's avatar
jcl committed
1001 1002 1003 1004 1005
  case uns::Nsel   :
    if (status) {
      *n    = nbody;
    } else {
      ok=false;
1006
    }
jcl's avatar
jcl committed
1007 1008 1009
  case uns::Pos   :
    if (status && getPos()) {
      *data = &getPos()[first*3];
1010
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1011 1012 1013 1014 1015 1016 1017
    } else {
      ok=false;
    }
    break;
  case uns::Vel  :
    if (status && getVel()) {
      *data = &getVel()[first*3];
1018
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1019 1020 1021 1022 1023 1024 1025
    } else {
      ok=false;
    }
    break;
  case uns::Acc  :
    if (status && getAcc()) {
      *data = &getAcc()[first*3];
1026
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1027 1028 1029 1030 1031 1032 1033
    } else {
      ok=false;
    }
    break;
  case uns::Pot  :
    if (status && getPot()) {
      *data = &getPot()[first];
1034
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1035 1036 1037 1038 1039 1040 1041
    } else {
      ok=false;
    }
    break;
  case uns::Mass  :
    if (status && getMass()) {
      *data = &getMass()[first];
1042
      *n    = nbody;//this->getNSel();
jcl's avatar
jcl committed
1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058
    } 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;
1059
    }
jcl's avatar
jcl committed
1060 1061 1062 1063 1064 1065
    break;
  case uns::Hsml :
    if (status && comp=="gas" && getHsml(*n)) {
      *data = getHsml(*n);
    } else {
      ok=false;
1066
    }
jcl's avatar
jcl committed
1067 1068 1069 1070 1071 1072
    break;
  case uns::Temp :
    if (status && comp=="gas" && getTemp(*n)) {
      *data = getTemp(*n);
    } else {
      ok=false;
1073
    }
jcl's avatar
jcl committed
1074 1075 1076 1077 1078 1079
    break;
  case uns::Age :
    if (status && comp=="stars" && getAge(*n)) {
      *data = getAge(*n);
    } else {
      ok=false;
1080
    }
jcl's avatar
jcl committed
1081 1082
    break;
  case uns::Metal :
1083
    if (status && comp=="gas" && this->ckloadBit(METAL_BIT)) {
jcl's avatar
jcl committed
1084 1085
      *data = getMetalGas(*n);
    } else
1086
      if (status && comp=="stars" && this->ckloadBit(METAL_BIT)) {
jcl's avatar
jcl committed
1087 1088 1089
      *data = getMetalStars(*n);
    } else {
      ok=false;
1090
    }
jcl's avatar
jcl committed
1091
    break;
1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104
  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
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 1132 1133 1134 1135 1136 1137
  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
1138 1139 1140 1141 1142 1143 1144
  case uns::Ssl :
    if (status && comp=="stars" && getSsl(*n)) {
      *data = getSsl(*n);
    } else {
      ok=false;
    }
    break;
jclamber's avatar
jclamber committed
1145 1146
  default: // unkown name
      if (comp=="STREAM") { // data stream reading
1147
        if (data_vector[name].size()==0) { // map key does not exist
1148
          std::vector<T> * p = &(data_vector[name]);
1149 1150 1151 1152 1153 1154
          int status=readStreamBlock(name,*p);
          if (status>0) {
            ok=true;
          }
        } else {
          ok = true;
jclamber's avatar
jclamber committed
1155
        }
1156 1157 1158 1159
        if (ok) {
          *n = data_vector[name].size();
          *data = &data_vector[name][0];
        }
jclamber's avatar
jclamber committed
1160
      }
1161

jcl's avatar
jcl committed
1162
  }
1163
  if (comp != "HEADER" && ok && !*data &&
1164 1165 1166
      (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
1167
    if (ok) {
1168
      std::cerr << "CSnapshotGadgetIn::getData name["<<name<<"]=" << CunsOut2<T>::s_mapStringValues[name] << "\n";
jcl's avatar
jcl committed
1169 1170 1171 1172 1173 1174 1175
    } else {
      std::cerr << "**WARNING** CSnapshotGadgetIn::getData Value ["<<name<<"] for component <"<<comp<<"> does not exist...\n";
    }
  }
  return ok;
}
// ============================================================================
1176 1177
// getData
// return requested array according 'name' selection
1178 1179
template <class T>
bool CSnapshotGadgetIn<T>::getData(const std::string name,int *n,T **data)
jcl's avatar
jcl committed
1180 1181 1182 1183
{
  bool ok=true;
  *data=NULL;
  *n = 0;
1184

1185
  switch(CunsOut2<T>::s_mapStringValues[name]) {
jcl's avatar
jcl committed
1186 1187
  case uns::Pos   :
    *data = getPos();