source: Sophya/trunk/Poubelle/archTOI.old/toiiter.cc@ 426

// toiiter.cc
// Eric Aubourg         CEA/DAPNIA/SPP   juillet 1999

#ifndef nb_max_bolo
#define _archeops               //  Acquisition Archeops (avec transputer) 
#define programme
extern "C" {
#include "archeops.h"
}
#endif


#include "toiiter.h"
#include "toiinterpolator.h"
#include "archparam.h"
#include <dirent.h>
#include <iostream.h>
#include "templocator.h"

// Format bloc GPS
// $GPGGA,hhmmss.ss,ddmm.mmmm,n,dddmm.mmmm,e,q,ss,y.y,a.a,z,

TOIIter::TOIIter() {
 // Toutes les valeurs initialisees par defaut. Le TOISvr les positionnera,
 // puis lancera l'initialisation
 
 file = NULL;
 directory = "";
 
 files.clear();

 isOnBoardRecorder = false;
   
 imes=0;

 tStart = -9.e99;
 tEnd   = 9.e99;
 
 //tBlock0 = -1;
 // perEch = -1;
 // Values for Trapani 99 = default values
 //tBlock0 = 1376.8358818;
 //perEch = 0.005836818076;
 
 trigMask = 0;
 rawIter = NULL;
 interp  = NULL;
 lastSample = -1;
 maxLookAhead = 10000;
 
 auxGPS = NULL;
 
 initDone = false;
 
 underSample = 1;
}

TOIIter::TOIIter(TOIIter const& x) {
  directory = x.directory;
  files = x.files;
  if (x.initDone && x.curFile != x.files.end()) curFile = files.find(*(x.curFile));
  isOnBoardRecorder = x.isOnBoardRecorder;
  imes = x.imes;
  tStart = x.tStart;
  tEnd = x.tEnd;
  trigMask = x.trigMask;
  infos = x.infos;   

  if (x.file)
    file = new ArcheopsFile(*x.file);
  else 
    file = NULL;
    
  if (x.rawIter) {
    rawIter = new TOIIter(*x.rawIter);
    interp = new TOIInterpolator[infos.size()];
    for (int i=0; i<infos.size(); i++)
      interp[i] = x.interp[i];
  } else {
    rawIter = NULL;
    interp  = NULL;
  }
  
  lastSample = x.lastSample;
  maxLookAhead = x.maxLookAhead;
  
  auxGPS = x.auxGPS; 
  if (auxGPS) auxGPS = auxGPS->clone();
  
  initDone = x.initDone;
  underSample = x.underSample;
}

TOIIter::~TOIIter() {
  delete file;
  delete rawIter;
  delete[] interp;
  delete auxGPS;
}

#ifdef __MWERKS__
#define filesep ':'
#else
#define filesep '/'
#endif

void TOIIter::Init() {
 if (initDone) return;
 initDone = true;
 // On a soit un repertoire, soit une liste de fichiers....
  if (directory == "") {
    if (files.empty()) {   // Ni repertoire, ni fichiers
      cerr << "toiiter : pas de repertoire, pas de fichiers" << endl;
      exit(-1);
    } else {
      // On a deja une liste de fichiers
    }
  } else { // On a un repertoire a explorer
    // On cherche soit les fichiers dans le repertoire donne, soit des fichiers
    // dans un sous-repertoire "arch-YY_MM_DD". Les fichiers ont un nom en
    // "hYY_MM_DD-hh_mm_ss".
    // Pour l'enregistreur de vol, les fichiers ont un nom en ARKxxxxxx.DAT
    if (directory[directory.length()-1] != filesep) 
      directory += filesep;
    DIR* dir = opendir(directory.c_str());
    struct dirent* ent;
    while ((ent = readdir(dir)) != NULL) {
      // si c'est un repertoire, avec un nom de jour, il faut l'explorer...
      if (!strncmp(ent->d_name, "arch-", 5)) {
        double mjd = ArcheopsFile::decodeMJD(ent->d_name+5) - 2./24.; // ENTIER + .5 en temps local!
        if (mjd >= tStart - 1. && mjd <= tEnd) {
          string direc2 = directory + ent->d_name + filesep;
          DIR* dir2 = opendir(direc2.c_str());
          struct dirent* ent2;
          while ((ent2 = readdir(dir2)) != NULL) {
            if (*ent2->d_name == 'h') {
              double mjd2 = ArcheopsFile::decodeMJD(ent->d_name+1) - 2./24.;
              if (mjd2 >= tStart - 1./24. && mjd2 <= tEnd) {
                files.insert(direc2 + ent2->d_name);
              }
            }
          }
        }
      } else {
        if (!isOnBoardRecorder && *ent->d_name == 'h') {
          double mjd = ArcheopsFile::decodeMJD(ent->d_name+1) - 2./24.; // $CHECK$ UTCOffset
          if (mjd >= tStart - 1./24. && mjd <= tEnd) {
            files.insert(directory + ent->d_name);
          }
        } else if (isOnBoardRecorder) {
          if (strncmp(ent->d_name, "ARK", 3) && strncmp(ent->d_name, "ark", 3)) continue;
          char * sfx = ent->d_name + strlen(ent->d_name) - 4;
          if (strcmp(sfx, ".DAT") && strcmp(sfx, ".dat")) continue;
          files.insert(directory + ent->d_name);
        }
      }
    }
    closedir(dir);
  }
    
  ScanFiles();
  
  curFile = files.begin();
  file = new ArcheopsFile((*curFile).c_str());
  cout << "opening file " << (*curFile).c_str() << endl;

  // On avance jusqu'a avoir au moins un bloc param et un bloc reglage,
  // car on ne peut rien faire sans...
  // Si on a des donnees de l'enregistreur de vol, pas de bloc param, et
  // on en simule un
  double oldTStart = tStart;
  tStart = -9.e99; // pour init, on accepte des blocs avant tstart....
  
  if (!file->lastParam()) {
    if (isOnBoardRecorder) {
      extern param_bolo  parametr;
      block_type_param block;
      block.param = parametr;
      valide_block((block_type_modele*)&block, block_param, 0);
      file->forceBlock((block_type_modele*)&block);
    } else {
      file->nextBlock(block_param_mask);
    }
  }
  if (!file->lastReglage()) file->nextBlock(block_reglage_mask);
  
  // On cherche un bloc GPS pour avoir la correspondance timestamp/UTC.
  // Pour le moment, on se fonde sur le premier bloc GPS. On pourra faire
  // mieux en prenant le min de tous les delta_T, a condition d'avoir un 
  // peu plus de details sur la facon dont le GPS est lu.
  
  if (archParam.acq.tBlock0 < 0) {
    archParam.acq.tBlock0 = file->getStartMJD();

    file->pushMark();
    if (file->lastGPS() || file->nextBlock(block_gps_mask)) {
      // le temps du bloc courant, en secondes
      double dt = file->blockNum() * file->perBlock();
      archParam.acq.tBlock0 = file->getGPSMJD() - dt/86400.;
    } else { // pas de bloc GPS...
      archParam.acq.tBlock0 = file->getStartMJD();
    }
    file->popMark();
  }
  tStart = oldTStart;  // on restaure
  
  if (archParam.acq.perEch < 0)
    archParam.acq.perEch = file->perEchant();

  bool hasInterp = false;
  
  trigMask = 0;
  for (vector<info>::iterator i = infos.begin(); i != infos.end(); i++) {
    if ((*i).interpolated) hasInterp = true;
    if ((*i).triggering) {
      switch ((*i).kind) {
      case boloTens:
      case boloTens2:
      case boloRaw:
      case boloRawCN:
      case boloTensCN:
        trigMask |= block_bolo_mask;
        break;
      case boloTens2T:
      case boloRes:
      case boloTemp:
         trigMask |= block_bolo_mask | block_reglage_mask;
        break;
     case dilDAC:
      case dilSwitch:
        trigMask |= block_dilution_mask;
      case boloGainAmpli:
      case boloDACV:
      case boloDACI:
        trigMask |= block_reglage_mask;
        break;
      case gpsTime:
      case longitude:
      case latitude:
      case altitude:
        trigMask |= block_gps_mask;
        break;
      case azimutBolo:
        trigMask |= block_bolo_mask;
        //file->needSSTProcessMask(SSTHandler::findPeriod);
        //trigMask |= block_sst_mask;
        break;
      case sstStarCnt:
      case sstStarZ:
      case sstStarF:
      case sstStarT:
        file->needSSTProcessMask(SSTHandler::findStars);
        trigMask |= block_sst_mask;
        break;
      case sstDiode:
      case sstDiodeCN:
        file->needSSTProcessMask(SSTHandler::permDiode);
        trigMask |= block_sst_mask;
        break;
      case sstChannel:
      case sstChannelCN:
        trigMask |= block_sst_mask;
        break;
      case gyroRaw:
      case gyroTens:
        trigMask |= block_gyro_mask;
        break;
      case gyroSpeed:
        trigMask |= block_gyro_mask; // $CHECK$ + info to calibrate gyros
        break;
      case alphaAxis:
      case deltaAxis:
      case alphaSst:
      case deltaSst:
      case alphaBolo:
      case deltaBolo:
        trigMask |= block_bolo_mask;
        //file->needSSTProcessMask(SSTHandler::findAxis);
        //trigMask |= block_sst_mask;
        break;
      }
    }
  }
  
  if (trigMask & (block_bolo_mask | block_sst_mask)) {
    imes = 9999;
  } else {
    imes = 0;
  }
  
  if (hasInterp) {
    rawIter = new TOIIter(*this);
    interp = new TOIInterpolator[infos.size()];
    for (int i=0; i<infos.size(); i++) {
      rawIter->infos[i].interpolated = false;
    }
    delete file; file = NULL; // on ne travaille plus sur le fichier directement...
  }
}

void TOIIter::ScanFiles() {
  file1stSamp.clear();
  cout << "Scanning all files" << endl;
  // Petite astuce pour les STL non conformes comme celles de digital
  // qui ne supportent pas files.erase(i) suivi de i++....
  set<string> copy = files;
  for (set<string>::iterator i = copy.begin(); i != copy.end(); i++) {
    ArcheopsFile fich((*i).c_str());
    if (fich.nextBlock()) {
      file1stSamp[*i] = fich.blockNum()*72; // premier numsample
      cout << "File " << *i << " 1st sample = " << fich.blockNum()*72 << endl;
    } else {
      cout << "File " << *i << " unrecoverable, skipping" << endl;
      files.erase(*i);
    }
  }
  cout << "Scan done" << endl;
  
  // Et maintenant, on ne garde que ceux qui tombent dans l'intervalle...
  copy = files;
  string prev="";
  for (set<string>::iterator i = copy.begin(); i != copy.end(); i++) {
    double smp = file1stSamp[*i];
    double t   = archParam.acq.tBlock0 + smp * archParam.acq.perEch/86400.;
    if (t>tEnd) {  // premier echantillon apres tEnd
     files.erase(*i);
     prev = "";
     continue;
    }
    if (t<tStart) { // premier echantillon avant tStart -> on vire le precedent si existe
      if (prev != "") {
        files.erase(prev);
      }
    }
    prev = *i;
  }
}

bool TOIIter::NextFile() {
  if (rawIter)
    return rawIter->NextFile();
    
  if (files.empty()) {
    return false;
  } else {
    if (curFile == files.end()) return false;
    curFile++;
    if (curFile == files.end()) return false;
    cout << "opening file " << (*curFile).c_str() << endl;
    ArcheopsFile* newfile = new ArcheopsFile((*curFile).c_str());
    newfile->grabLastBlocs(*file);
    delete file;
    file = newfile;
    return true;
  }
}

bool TOIIter::Next() {
  for (int ii=0; ii<underSample; ii++)
    if (!Next1()) return false;
  return true;
}

bool TOIIter::Next1() {
  if (!initDone) Init();
  while (1) {
    if (!NextSample()) return false; // end of files
    double t = getMJD();
    if (t < tStart) continue;
    if (t > tEnd) return false;
    return true;
  }
}

bool TOIIter::NextSample() {
  if (rawIter) {  // Delegation pour interpolation
    // Trouve prochain sample disponible
    for (int k=0; k<2; k++) {
      long smp = 2147483647L;
      for (int i=0; i<infos.size(); i++) {
        long ss = interp[i].nextSample(lastSample+1);
        if (ss > 0 && ss < smp) smp=ss;
      }
      if (smp != 2147483647L) {
        lastSample = smp;
        break;
      }
      if (!fetchAhead())  // tout le monde etait en bout de course,
        return false;     // on lit un echantillon, ca suffit, d'ou le k<2
    }
    // Verifie que tous les interpolateurs ont assez de donnees pour
    // trouver la valeur correspondante
    for (int i=0; i<infos.size(); i++) {
      //rif (infos[i].interpolated) 
        while (interp[i].needMoreDataFor(lastSample) && 
               rawIter->getSampleIndex() - lastSample < maxLookAhead)
          if (!fetchAhead()) return false;
    }
    
    // On est pret...
    return true;
  }

  // trigger sur info indexee dans bloc bolo, bloc gyro ou bloc sst ?
  if (trigMask & (block_bolo_mask | block_sst_mask | block_gyro_mask )) {
    imes++;
    if (imes < file->nEchBlock()) return true;
    imes = 0;
  }
  
  // soit pas d'info indexee, soit fin bloc courant...
  while (1) {
    if (file->nextBlock(trigMask)) {
      while (file->sameBlockNumAhead()) {  // tant que meme numero de bloc, on lit
        if (!file->nextBlock()) {          // fin de fichier ?
          if (NextFile()) file->nextBlock(); // fichier suivant
          else return false;    // tout fini
        }
      }
      return true;
    }
    if (!NextFile()) return false;
  }
}

/* double TOIIter::getTime() {                // MJD 
   // le temps du bloc courant, en secondes
   double dt = file->blockNum() * file->perBlock();
   return tBlock0 + dt/86400. + imes*file->perEchant()/86400.;
 }
 */
 
bool TOIIter::canGetValue(int column) {
  if (!initDone) Init();
   if (column < 0 || column >= infos.size()) return false;
   TOIKind kind = infos[column].kind;
   if (auxGPS && 
       (kind == longitude || kind == latitude || kind == altitude || kind == tsid)) {
         double dummy;
         return auxGPS->getLocation(getMJD(), dummy, dummy, dummy) == 0;
   }
   if (rawIter) {
       return interp[column].canGet(lastSample);
   }
   int index = infos[column].index;
   switch (kind) {
   case sampleNum:
   case internalTime:
   case mjd:
     return true;
   case boloTens:
   case boloTens2:
   case boloTens2T:
   case boloRes:
   case boloRawCN:
     if (imes==0 && file->llastBolo()==NULL) return false;
     return file->lastBolo() != NULL;
   case boloTensCN:
   case boloRaw:
     return file->lastBolo() != NULL;
   case boloGainAmpli:
   case boloDACV:
   case boloDACI:
     return file->lastReglage() != NULL;
   case dilDAC:
   case dilSwitch:
     return file->lastDilution() != NULL;
   case sstDiode:
   case sstChannel:
   case sstDiodeCN:
   case sstChannelCN:
     return file->lastSST() != NULL;
   case sstStarCnt:
   case sstStarZ:
   case sstStarF:
   case sstStarT:{
     if (file->lastSST() == NULL) return false;
     int n = file->getNumbStar(imes);
     return (n > 0 && index < n);
   }
   case gyroRaw:
   case gyroTens:
   case gyroSpeed:
     return (file->lastGyro() != NULL);
   case gpsTime:
     return file->hasGPSTime();
   case longitude:
   case latitude:
     return file->hasGPSPos();
   case altitude:
     return file->hasGPSAlt();
   case tsid:
     return file->hasGPSPos();
   case azimutBolo:
   case alphaAxis:
   case deltaAxis:
   case alphaBolo:
   case deltaBolo:
     return true;
     //return (file->lastGPS() != NULL && file->lastSST() != NULL);
   case alphaSst:
   case deltaSst:
     return false;
     
   case boloTemp:
     return false;
   }
   return false;
}


double TOIIter::getValue(int column) {
  if (!initDone) Init();
   if (column < 0 || column >= infos.size()) return -1;
   TOIKind kind = infos[column].kind;
   if (auxGPS && 
       (kind == longitude || kind == latitude || kind == altitude || kind == tsid)) {
         double lat,lon,alt;
         if (auxGPS->getLocation(getMJD(), lat, lon, alt)) return -99999;
         if (kind == longitude) return lon;
         if (kind == latitude)  return lat;
         if (kind == altitude)  return alt;
         if (kind == tsid) {
           tSid.setLongitude(lon);
           return tSid.getLST(getMJD());
         }
   }
   if (rawIter) {
     if (infos[column].interpolated)
       return interp[column].getIValue(lastSample);
     else
       return interp[column].getEValue(lastSample);
   }
   int index = infos[column].index;
   switch (kind) {
   case sampleNum:
     return getSampleIndex();
   case internalTime:
     return getSampleIndex() * archParam.acq.perEch;
   case mjd:
     return getMJD();
   case boloTens:
     return file->getMuVBolo(index, imes);
   case boloTens2:
     return file->getMuVBolo2(index, imes);
   case boloRaw:
     return file->getRawBolo(index, imes);
   case boloRawCN:
     return file->getRawBoloCN(index, imes);
   case boloTensCN:
     return file->getMuVBoloCN(index, imes);
   case boloGainAmpli:
     return file->getGainAmpli(index);
   case boloDACV:
     return file->getDACV(index);
   case boloDACI:
     return file->getDACI(index);
   case boloTens2T:
     return file->getMuVBolo2T(index, imes);
   case boloRes:
     return file->getBoloRes(index, imes);
   case dilDAC:
     return file->getADCDil(index);
   case dilSwitch:
     return file->getSwitchDil();
   case sstDiode:
     return file->getSSTSignal(index, imes);
   case sstChannel:
     return file->getSSTRawSignal(index, imes);
   case sstDiodeCN:
     return file->getSSTSignalCN(index, imes);
   case sstChannelCN:
     return file->getSSTRawSignalCN(index, imes);
   case sstStarCnt:
     return file->getNumbStar(imes);
   case sstStarZ:
     return file->getSSTStarZ(index, imes);
   case sstStarF:
     return file->getSSTStarF(index, imes);
   case sstStarT:
     return file->getSSTStarT(index, imes);
   case gyroRaw:
     return file->getGyroRaw(index, imes);
   case gyroTens:
     return file->getGyroTens(index, imes);
   case gyroSpeed:
     return file->getGyroSpeed(index, imes);
   case gpsTime:
     return file->getGPSUTC();
   case longitude:
     return file->getGPSLong();
   case latitude:
     return file->getGPSLat();
   case altitude:
     return file->getGPSAlt();
   case tsid:
     tSid.setLongitude(file->getGPSLong());
     return tSid.getLST(getMJD());
   case azimutBolo:
   case alphaAxis:
   case deltaAxis: 
   case alphaBolo:
   case deltaBolo: {
        double lat, lon;
        if (auxGPS) {
          double alti;
          auxGPS->getLocation(getMJD(), lat, lon, alti);
        } else {
          lon = file->getGPSLong();
          lat = file->getGPSLat();
        }
        tSid.setLongitude(lon);
        double ts = tSid.getLST(getMJD());
        tempLocator.setEarthPos(lon, lat);
        tempLocator.setTSid(ts);
        if ((kind) == alphaAxis) return tempLocator.getAlphaZenith();
        if ((kind) == deltaAxis) return tempLocator.getDeltaZenith();
        if ((kind) == azimutBolo) return tempLocator.getAzimutBolo(getSampleIndex(),index);
        if ((kind) == alphaBolo) return tempLocator.getAlphaBolo(getSampleIndex(),index);
        if ((kind) == deltaBolo) return tempLocator.getDeltaBolo(getSampleIndex(),index);
        return 0;
     }
   }
   return -1;
}
 
bool   TOIIter::newValue(int column) {
  if (!initDone) Init();
   if (column < 0 || column >= infos.size()) return false;
   TOIKind kind = infos[column].kind;
   if (auxGPS && 
       (kind == longitude || kind == latitude || kind == altitude)) {
      return true;
   }
   if (rawIter) {
       return interp[column].isNewValue(lastSample);
   }
   switch (kind) {
     case sampleNum:
     case internalTime:
     case mjd:
     case tsid:
       return true;
     case boloTens:
     case boloTens2:
     case boloTens2T:
     case boloRes:
       return file->blockNum() == file->getBoloBlockNum();
     case boloRaw:
     case boloRawCN:
     case boloTensCN:
       return file->blockNum() == file->getBoloBlockNum();
     case boloGainAmpli:
       return file->blockNum() == file->getReglageBlockNum() && imes==0;
     case boloDACI:
     case boloDACV:
       return file->blockNum() == file->getReglageBlockNum() && imes==0;
     case dilDAC:
     case dilSwitch:
       return file->blockNum() == file->getDilutionBlockNum() && imes==0;
     case sstChannel:
     case sstDiode:
     case sstStarCnt:
     case sstStarZ:
     case sstStarF:
     case sstStarT:
       return file->blockNum() == file->getSSTBlockNum();
     case sstChannelCN:
     case sstDiodeCN:
       return file->blockNum() == file->getSSTCompBlockNum();
     case gyroRaw:
     case gyroTens:
     case gyroSpeed:
       return file->blockNum() == file->getGyroBlockNum();
     case gpsTime:
       return file->blockNum() == file->getGPSBlockNum() && imes==0;
     case longitude:
       return file->blockNum() == file->getGPSBlockNum() && imes==0;
     case latitude:
       return file->blockNum() == file->getGPSBlockNum() && imes==0;
     case altitude:
       return file->blockNum() == file->getGPSBlockNum() && imes==0;
     case azimutBolo:
     case alphaBolo:
     case deltaBolo:
       return true; // $CHECK$ with SSTHandler
     case alphaAxis:
       return true; // $CHECK$ with SSTHandler
     case deltaAxis:
       return true; // $CHECK$ with SSTHandler
   }
   return false;
}
 
bool   TOIIter::extendValue(int column) {
   return (!infos[column].interpolated && !newValue(column));
}
 
bool   TOIIter::interpValue(int column) {
   return (infos[column].interpolated && !newValue(column));
}

bool TOIIter::isTrig(int column) {
   if (column < 0 || column >= infos.size()) return false;
   return infos[column].triggering;
}

 
TOIKind TOIIter::getKind(int column) {
   if (column < 0 || column >= infos.size()) return (TOIKind)-1;
   return infos[column].kind;
}
 
int TOIIter::getIndex(int column) {
   if (column < 0 || column >= infos.size()) return (TOIKind)-1;
   return infos[column].index;
}

int TOIIter::getColTOI(TOIKind kind, int index) {
   for (int i=0; i<infos.size(); i++)
     if (infos[i].kind == kind && infos[i].index == index)
       return i;
   return -1;
}

bool TOIIter::canGetTOI(TOIKind kind, int index) {
   int col = getColTOI(kind, index);
   if (col<0) return false;
   return canGetValue(col);
}

double TOIIter::getTOI(TOIKind kind, int index) {
   int col = getColTOI(kind, index);
   if (col<0) return -9.e99;
   return getValue(col);
}

   
int TOIIter::getBlockSampleIndex() {
   return imes;
}
 
int TOIIter::getUnderSampling() {
  return underSample;
}

int TOIIter::getSampleIndex() {
  if (!initDone) Init();
   if (file) {
     return file->blockNum() * file->nEchBlock() + imes;
   } else {
     return lastSample;
   }
}

double TOIIter::getMJD() {
  if (!initDone) Init();
  int sample = getSampleIndex();
  return archParam.acq.tBlock0 + sample*archParam.acq.perEch/86400.;
}
 
bool TOIIter::fetchAhead() { // Seulement si delegation
  if (!rawIter) return false;
  if (!rawIter->Next()) return false;
  long sample = rawIter->getSampleIndex();
  for (int i=0; i<infos.size(); i++) {
    if (rawIter->canGetValue(i) && rawIter->newValue(i)) {
      interp[i].enterValue(rawIter->getValue(i), sample);
    }
  }
  return true;
}

block_type_param* TOIIter::lastParam() {
  if (!initDone) Init();
  if (file) return file->lastParam(); 
  else return rawIter->lastParam();
}