source: Sophya/trunk/Poubelle/archTOI.old/archeopsfile.cc@ 315

#define utilitaires_de_block_archeops
#include "archeopsfile.h"
extern "C" {
#include "compress.h"
}
#include <iostream.h>

// BlockSet is the current "state" of the file reader automate.
// It contains the last blocks of each kind, the current block and
// the look-ahead block.

class BlockSet {
public:
  BlockSet();
  BlockSet(BlockSet const&);
  ~BlockSet();
  void setBloc(block_type_modele const& blk);
  void setRawBloc(block_type_modele const& blk);

  block_type_param*         lastParam;
  block_type_journal*       lastJournal;
  block_type_reglage*       lastReglage;
  block_type_dilution*      lastDilution;
  block_type_gps*           lastGPS;
  block_type_une_periode*   lastUnePeriode;
  block_type_synchro_sol*   lastSynchroSol;
  block_type_pointage_sol*  lastPointageSol;
  block_type_bolo*          lastBolo;   // can be uncompressed bolo_comp
  block_type_bolo*          llastBolo;   
  block_type_gyro*          lastGyro;
  block_type_sst*           lastSST;
  block_type_bolo_comprime* lastBoloComp;   // can be uncompressed bolo_comp
  block_type_gyro_comprime* lastGyroComp;
  block_type_sst_comprime*  lastSSTComp;
    
  block_type_modele         curBlock;
  block_type_modele         peekBlock;
  
  SSTHandler sstHandler;
};

BlockSet::BlockSet() {
  lastParam      = NULL;
  lastJournal    = NULL;
  lastReglage    = NULL;
  lastDilution   = NULL;
  lastGPS        = NULL;
  lastUnePeriode = NULL;
  lastSynchroSol = NULL;
  lastPointageSol= NULL;
  lastBolo       = NULL;   
  llastBolo      = NULL;   
  lastGyro       = NULL;
  lastSST        = NULL;
  lastBoloComp   = NULL;   
  lastGyroComp   = NULL;
  lastSSTComp    = NULL;
  
  memset(&curBlock,0,sizeof(block_type_modele));
  memset(&peekBlock,0,sizeof(block_type_modele));
}

BlockSet::~BlockSet() {
  delete lastParam;
  delete lastJournal;
  delete lastReglage;
  delete lastDilution;
  delete lastGPS;
  delete lastUnePeriode;
  delete lastSynchroSol;
  delete lastPointageSol;
  delete lastBolo;   
  delete llastBolo;   
  delete lastGyro;
  delete lastSST;
  delete lastBoloComp;   
  delete lastGyroComp;
  delete lastSSTComp;
}

BlockSet::BlockSet(BlockSet const& x) 
: sstHandler(x.sstHandler)
{
  lastParam      = NULL;
  lastJournal    = NULL;
  lastReglage    = NULL;
  lastDilution   = NULL;
  lastGPS        = NULL;
  lastUnePeriode = NULL;
  lastSynchroSol = NULL;
  lastPointageSol= NULL;
  lastBolo       = NULL;   
  llastBolo      = NULL;   
  lastGyro       = NULL;
  lastSST        = NULL;
  lastBoloComp   = NULL;   
  lastGyroComp   = NULL;
  lastSSTComp    = NULL;
  
  if (x.lastParam) {
    lastParam = new block_type_param;
    *lastParam = *x.lastParam;
  }
  if (x.lastJournal) {
    lastJournal = new block_type_journal;
    *lastJournal = *x.lastJournal;
  }
  if (x.lastReglage) {
    lastReglage = new block_type_reglage;
    *lastReglage = *x.lastReglage;
  }
  if (x.lastDilution) {
    lastDilution = new block_type_dilution;
    *lastDilution = *x.lastDilution;
  }
  if (x.lastGPS) {
    lastGPS = new block_type_gps;
    *lastGPS = *x.lastGPS;
  }
  if (x.lastUnePeriode) {
    lastUnePeriode = new block_type_une_periode;
    *lastUnePeriode = *x.lastUnePeriode;
  }
  if (x.lastSynchroSol) {
    lastSynchroSol = new block_type_synchro_sol;
    *lastSynchroSol = *x.lastSynchroSol;
  }
  if (x.lastPointageSol) {
    lastPointageSol = new block_type_pointage_sol;
    *lastPointageSol = *x.lastPointageSol;
  }
  if (x.lastBolo) {
    lastBolo = new block_type_bolo;
    *lastBolo = *x.lastBolo;
  }
  if (x.llastBolo) {
    llastBolo = new block_type_bolo;
    *llastBolo = *x.llastBolo;
  }
  if (x.lastGyro) {
    lastGyro = new block_type_gyro;
    *lastGyro = *x.lastGyro;
  }
  if (x.lastSST) {
    lastSST = new block_type_sst;
    *lastSST = *x.lastSST;
  }
  if (x.lastBoloComp) {
    lastBoloComp = new block_type_bolo_comprime;
    *lastBoloComp = *x.lastBoloComp;
  }
  if (x.lastGyroComp) {
    lastGyroComp = new block_type_gyro_comprime;
    *lastGyroComp = *x.lastGyroComp;
  }
  if (x.lastSSTComp) {
    lastSSTComp = new block_type_sst_comprime;
    *lastSSTComp = *x.lastSSTComp;
  }
  
  memcpy(&curBlock,&x.curBlock,sizeof(block_type_modele));
  memcpy(&peekBlock,&x.peekBlock,sizeof(block_type_modele));

}

void BlockSet::setBloc(block_type_modele const& blk)
{
  int kind = type_block(&blk);
  switch (kind) {
    case block_param:
      if (!lastParam) lastParam = new block_type_param;
      memcpy(lastParam, &blk, sizeof(block_type_param));
      // Les fichiers fournis contiennent des valeurs debiles...
      if (lastParam->param.n_max_bolo < 0)
        lastParam->param.n_max_bolo = nb_max_bolo;
      if (lastParam->param.n_per_block < 0)
        lastParam->param.n_per_block = nb_per_block;
      if (lastParam->param.n_max_mes_per < 0)
        lastParam->param.n_max_mes_per = nb_max_mes_per;
      break;
    case block_journal:
      if (!lastJournal) lastJournal = new block_type_journal;
      memcpy(lastJournal, &blk, sizeof(block_type_journal));
      break;
    case block_reglage:
      if (!lastReglage) lastReglage = new block_type_reglage;
      memcpy(lastReglage, &blk, sizeof(block_type_reglage));
      break;
    case block_dilution:
      if (!lastDilution) lastDilution = new block_type_dilution;
      memcpy(lastDilution, &blk, sizeof(block_type_dilution));
      break;
    case block_gps:
      if (!lastGPS) lastGPS = new block_type_gps;
      memcpy(lastGPS, &blk, sizeof(block_type_gps));
      break;
    case block_une_periode:
      if (!lastUnePeriode) lastUnePeriode = new block_type_une_periode;
      memcpy(lastUnePeriode, &blk, sizeof(block_type_une_periode));
      break;
    case block_synchro_sol:
      if (!lastSynchroSol) lastSynchroSol = new block_type_synchro_sol;
      memcpy(lastSynchroSol, &blk, sizeof(block_type_synchro_sol));
      break;
    case block_pointage_sol:
      if (!lastPointageSol) lastPointageSol = new block_type_pointage_sol;
      memcpy(lastPointageSol, &blk, sizeof(block_type_pointage_sol));
      break;
    case block_bolo:
      if (!lastBolo) lastBolo = new block_type_bolo;
      else if (!llastBolo) {
        llastBolo = new block_type_bolo;
        memcpy(llastBolo, lastBolo, sizeof(block_type_bolo));
      }
      memcpy(lastBolo, &blk, sizeof(block_type_bolo));
      break;
    case block_gyro:
      if (!lastGyro) lastGyro = new block_type_gyro;
      memcpy(lastGyro, &blk, sizeof(block_type_gyro));
      break;
    case block_sst:
      if (!lastSST) lastSST = new block_type_sst;
      memcpy(lastSST, &blk, sizeof(block_type_sst));
      break;
  }
}

void BlockSet::setRawBloc(block_type_modele const& blk)
{
  int kind = type_block(&blk);
  switch (kind) {
    case block_bolo_comprime:
      if (!lastBoloComp) lastBoloComp = new block_type_bolo_comprime;
      memcpy(lastBoloComp, &blk, sizeof(block_type_bolo_comprime));
      break;
    case block_gyro_comprime:
      if (!lastGyroComp) lastGyroComp = new block_type_gyro_comprime;
      memcpy(lastGyroComp, &blk, sizeof(block_type_gyro_comprime));
      break;
    case block_sst_comprime:
      if (!lastSSTComp) lastSSTComp = new block_type_sst_comprime;
      memcpy(lastSSTComp, &blk, sizeof(block_type_sst_comprime));
      break;
  }
}



ArcheopsFile::ArcheopsFile(string const& fname) {
  f = fopen(fname.c_str(), "rb");
  fn = fname;
  if (!f) throw ArchExc("file not found");
  fseek(f,0,SEEK_END);
  fLen = ftell(f);
  curPos = -1;
  peekPos = -1;
  curKind = -1;
  curRawKind = -1;
  blockSet = new BlockSet;
  utcOffset=2;
  computeMJD(fname);
}

ArcheopsFile::ArcheopsFile(ArcheopsFile const& x) {
  blockSet = x.blockSet ? new BlockSet(*x.blockSet) : NULL;
  fposStack = x.fposStack;
  
  stack<BlockSet*> tmp;
  ArcheopsFile& y = (ArcheopsFile&) x;
  while (!y.blockStack.empty()) {
    tmp.push(y.blockStack.top());
    y.blockStack.pop();
  }
  while (!tmp.empty()) {
    y.blockStack.push(tmp.top());
    blockStack.push(new BlockSet(*tmp.top()));
    tmp.pop();
  }

  //curBlock   = x.curBlock;
  curKind    = x.curKind;
  curRawKind = x.curRawKind;
  curPos     = x.curPos;
  peekPos    = x.peekPos;
  fLen       = x.fLen;
  fn         = x.fn;
  f = fopen(fn.c_str(), "rb");
}

ArcheopsFile::~ArcheopsFile() {
  if (f) fclose(f);
  delete blockSet;
  while (!blockStack.empty())
    { delete blockStack.top(); blockStack.pop();}
}

void ArcheopsFile::grabLastBlocs(ArcheopsFile const& oldFile) {
  delete blockSet;
  blockSet = new BlockSet(*oldFile.blockSet);
  //curBlock = oldFile.curBlock;
  //peekBlock = oldFile.peekBlock;
  peekPos = 0;
  setUTCOffset(oldFile.utcOffset);
}


def_nom_block

bool ArcheopsFile::nextBlock() {
  // si pas de peek, alors en lire un pour commencer...
  if (peekPos<0) {
    peekPos=0; if (!nextBlock()) return false;
  }
  block_type_modele* curBlock=currentBlock();
  block_type_modele* peekBlock=aheadBlock();
  
  memcpy(curBlock, peekBlock, sizeof(block_type_modele));
  
  if (peekPos+12 > fLen) return false;
  fseek(f,peekPos,SEEK_SET);
  size_t read = fread(peekBlock,1,sizeof(block_type_modele),f);
  swapEntete(peekBlock);
  if (read < longueur_block(peekBlock)) return false;
  swapContent(peekBlock);
  if (verifie_block(peekBlock) != block_correct) {
  printf("block invalide\n"); throw ArchExc("invalid block");
  }
  peekPos += longueur_block(peekBlock);
  if (curPos < 0) {
    curPos = 0; return true;
  }
  curPos += longueur_block(curBlock);
  curRawKind = curKind = type_block(curBlock);
  //printf("block %d  :  %s \n",numero_block(&curBlock),nom_block[curKind]);
  postProcessBlock();
  saveCurBlock();
  return true;
}

bool ArcheopsFile::nextBlock(long mask) {
  if (!mask) return false;
  while (1) {
    if (!nextBlock()) return false;
    if (( 1 << curRawKind) & mask) return true;
  }
}

int  ArcheopsFile::blockKind() {
  return curKind;
}

int  ArcheopsFile::blockRawKind() {
  return curRawKind;
}

int  ArcheopsFile::blockNum() {
  return numero_block(currentBlock());
}

#define  bitmot                 24              // nb de bit horloge dans un mot ADC

double ArcheopsFile::perEchant() { // periode d'echantillonage pour le dernier bloc reglage
  double p,f1,f2,f3;
  int pp;
  pp=lastReglage()->reglage.horloge.periode;
  p=pp/5.;
  f1=1000/p;f2=f1/bitmot;f3=f2*1000./(double)(lastReglage()->reglage.horloge.nb_mesures);
  return 0.5/f3;                                        //  2 fois la frequence de modulation
}

double ArcheopsFile::perBlock() {  // duree (en secondes) correspondant a un bloc bolo
  return perEchant() * (double)lastParam()->param.n_per_block*2.;
}

int ArcheopsFile::nEchBlock() {  // Nb d'echantillons dans un bloc
  return lastParam()->param.n_per_block*2.;
}

string  ArcheopsFile::blockKdName() {
  return nom_block[curKind];
}

string  ArcheopsFile::blockRawKdName() {
  return nom_block[curRawKind];
}


block_type_modele*  ArcheopsFile::currentBlock() {
  //if (curPos<0) return NULL;
  return &blockSet->curBlock;
}

block_type_modele*  ArcheopsFile::aheadBlock() {
  //if (peekPos<0) return NULL;
  return &blockSet->peekBlock;
}

bool ArcheopsFile::sameBlockNumAhead() {
  if (curPos<0) return false;
  return (numero_block(&blockSet->curBlock) == numero_block(&blockSet->peekBlock));
}


block_type_param* ArcheopsFile::lastParam() {
  return blockSet->lastParam;
}
block_type_journal* ArcheopsFile::lastJournal() {
  return blockSet->lastJournal;
}
block_type_reglage* ArcheopsFile::lastReglage() {
  return blockSet->lastReglage;
}
block_type_dilution* ArcheopsFile::lastDilution() {
  return blockSet->lastDilution;
}
block_type_gps* ArcheopsFile::lastGPS() {
  return blockSet->lastGPS;
}
block_type_une_periode* ArcheopsFile::lastUnePeriode() {
  return blockSet->lastUnePeriode;
}
block_type_synchro_sol* ArcheopsFile::lastSynchroSol() {
  return blockSet->lastSynchroSol;
}
block_type_pointage_sol* ArcheopsFile::lastPointageSol() {
  return blockSet->lastPointageSol;
}
block_type_bolo* ArcheopsFile::lastBolo() {
  return blockSet->lastBolo;
} 
block_type_bolo* ArcheopsFile::llastBolo() {
  return blockSet->llastBolo;
} 
block_type_gyro* ArcheopsFile::lastGyro() {
  return blockSet->lastGyro;
}
block_type_sst* ArcheopsFile::lastSST() {
  return blockSet->lastSST;
}
block_type_bolo_comprime* ArcheopsFile::lastBoloComp() {
  return blockSet->lastBoloComp;
} 
block_type_gyro_comprime* ArcheopsFile::lastGyroComp() {
  return blockSet->lastGyroComp;
}
block_type_sst_comprime* ArcheopsFile::lastSSTComp() {
  return blockSet->lastSSTComp;
}

void ArcheopsFile::postProcessBlock() {
  switch (curKind) {
     case block_bolo_comprime: {
       blockSet->setRawBloc(blockSet->curBlock);
       block_type_bolo blk2;
       block_type_bolo_comprime* blk = (block_type_bolo_comprime*) &blockSet->curBlock;
       for(int j=0;j<nb_bolo_util;j++)  
                {
                   decompress_7_2(blk->data_bolo[j],blk2.data_bolo[j],nb_per_block*2);
                }
       valide_block((block_type_modele*)&blk2,block_bolo,numero_block(blk));    
       memcpy(&blockSet->curBlock,&blk2,sizeof(blk2));
       curKind = block_bolo;
       break;
     }
     case block_sst_comprime: {  // $CHECK$ TBD
     }
     case block_sst : {
       blockSet->sstHandler.ProcessBlock((block_type_sst*)&blockSet->curBlock);
     }
  }
}

void ArcheopsFile::saveCurBlock() {
  blockSet->setBloc(blockSet->curBlock);
}

void ArcheopsFile::pushMark() {    // push current file position, and "last" blocks`
  fposStack.push(curPos);
  fposStack.push(peekPos);
  blockStack.push(new BlockSet(*blockSet));
}

void ArcheopsFile::popMark() {    // pops last file position and "last" blocks
  if (! blockStack.empty()) {
    delete blockSet;
    blockSet = blockStack.top();
    blockStack.pop();
    peekPos = fposStack.top();
    fposStack.pop();
    curPos = fposStack.top();
    fposStack.pop();
  }
}

typedef unsigned int4 uint_4;
typedef unsigned short uint_2;

static inline void bswap4(void* p)
{
  uint_4 tmp = *(uint_4*)p;
  *(uint_4*)p = ((tmp >> 24) & 0x000000FF) | 
                ((tmp >> 8)  & 0x0000FF00) |
                ((tmp & 0x0000FF00) << 8)  |
                ((tmp & 0x000000FF) << 24);
}

static inline void bswap2(void* p)
{
  uint_2 tmp = *(uint_2*)p;
  *(uint_2*)p = ((tmp >> 8) & 0x00FF) | 
                ((tmp & 0x00FF) << 8);
}

#ifdef __DECCXX
#define SWAP
#endif
#if defined(Linux)  || defined(linux)
#define SWAP
#endif

void ArcheopsFile::swapEntete(block_type_modele* blk) {
#ifdef SWAP
  bswap4(&(blk->debut));
  bswap4(&(blk->code1));
  bswap4(&(blk->code2));
  long lg = longueur_block(blk);
  bswap4(((char*)blk) + lg - 4);
#endif
}

void ArcheopsFile::swapContent(block_type_modele* blk) {
#ifdef SWAP
  int typ = type_block(blk);

  if (typ == block_gps) return; // char only, no swap
  if (typ == block_une_periode) { // une_periode, des shorts
    block_type_une_periode* b = (block_type_une_periode*) blk;
    for (int i=0; i<nb_max_bolo; i++)
      for (int j=0; j<nb_max_mes_per; j++)
        bswap2(&b->bol_per[i][j]);
    return;
  }

  for (int i=0; i<longueur_block(blk)/4-4; i++)
    bswap4(blk->mot + i);
  // On deswappe ce qui ne devait pas etre swappe...
  switch (typ) {
  case block_param: {
    block_type_param* b = (block_type_param*) blk;
    for (int i=0; i<nb_max_bolo; i++)
      for (int j=0; j<8; j++)
        bswap4(&b->param.bolo[i].bolo_nom[4*j]);
  }
  
  }
#endif
}
  
  


static int mlen[] = {31,28,31,30,31,30,31,31,30,31,30,31};

void ArcheopsFile::computeMJD(string const& fname) {
  //h99_04_29-15h36m22
  short y,m,d,hh,mm,ss;
  char const* p = fname.c_str();
  char const* p2 = p  + fname.length()-1;
  while (*p2 != ':' && *p2 != '/' && p2 != p) p2--;
  if (*p2 == ':' || *p2 == '/') p2++;
  if (*p2 == 'h') p2++;
  y  = atoi(p2); p2+=3;
  m  = atoi(p2); p2+=3;
  d  = atoi(p2); p2+=3;
  hh = atoi(p2); p2+=3;
  mm = atoi(p2); p2+=3;
  ss = atoi(p2);
  
  if (y<50) y += 100;
  // 1. depuis 0/1/97 minuit
  rawMJD = (int) (365.25 * (y-97));
  for (int i=0; i<m-1; i++)
    rawMJD += mlen[i];
  if (y%4 == 0 && m > 2) rawMJD++;
  rawMJD += d;
  rawMJD += hh/24. + mm/24./60. + ss/24./60./60;
  
  rawMJD += 448.5; // 0/1/97 minuit
  startMJD = rawMJD - utcOffset/24.;
}

void ArcheopsFile::setUTCOffset(int UTCOffset) {
  utcOffset = UTCOffset;
  startMJD = rawMJD - utcOffset/24.;
}

double ArcheopsFile::getStartMJD() {
  return startMJD;
}

int ArcheopsFile::getBoloBlockNum() {
  if (!lastBolo()) return -1;
  return numero_block(lastBolo());
}

int ArcheopsFile::getSSTBlockNum() {
  if (!lastSST()) return -1;
  return numero_block(lastSST());
}


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


int ArcheopsFile::getGPSBlockNum() {
  if (!lastGPS()) return -1;
  return numero_block(lastGPS());
}

double ArcheopsFile::getGPSUTC() {  // en secondes depuis minuit UTC jour courant
  if (!lastGPS()) return -9.e99;
  char* p = lastGPS()->gps;
  if (strncmp(p, "$GPGGA,", 7)) return -9.e99;
  p += 7;
  double t;
  sscanf(p, "%lg", &t);
  int h = int(t/10000); t -= h*10000;
  int m = int(t/100);   t -= m*100;
  return h*3600. + m*60. + t;
}

double ArcheopsFile::getGPSMJD() {  // modified julian day du dernier bloc GPS, JD - 2450000
  double t = getGPSUTC()/86400.;
  if (t<0) return t;
  if (t > (startMJD - int(startMJD))) {
    return int(startMJD) + t;
  } else {
    return int(startMJD) + 1. + t;
  }
}

double ArcheopsFile::getGPSLat() {  // degres, +  = NORD
  if (!lastGPS()) return -9.e99;
  char* p = lastGPS()->gps;
  if (strncmp(p, "$GPGGA,", 7)) return -9.e99;
  char* fence = p+80;
  p += 7;
  while (*p != ',' && p<fence) p++;
  if (*p != ',') return -9.e99; 
  p++;
  double t;
  sscanf(p, "%lg", &t);
  int d = int(t/100); t -= d*100;
  t = d + t/60;
  while (*p != ',' && p<fence) p++;
  if (*p != ',') return -9.e99; 
  p++;
  if (*p == 'S') t = -t;
  return t;
}

double ArcheopsFile::getGPSLong() { // degres, +  = EST
  if (!lastGPS()) return -9.e99;
  char* p = lastGPS()->gps;
  if (strncmp(p, "$GPGGA,", 7)) return -9.e99;
  char* fence = p+80;
  p += 7;
  for (int i=0; i<3; i++) {
    while (*p != ',' && p<fence) p++;
    if (*p != ',') return -9.e99; 
    p++;
  }
  double t;
  sscanf(p, "%lg", &t);
  int d = int(t/100); t -= d*100;
  t = d + t/60;
  while (*p != ',' && p<fence) p++;
  if (*p != ',') return -9.e99; 
  p++;
  if (*p == 'W') t = -t;
  return t;
}



// Bolo

long ArcheopsFile::getRawBolo(int ibolo, int imesure) {    // donnee brute, avec seulement soustraction offset
  int nb_coups,aa;  
  block_type_bolo* blk = imesure >= 0 ? lastBolo() : llastBolo();
  if (imesure < 0) imesure += nEchBlock();
  if (!blk) return 0;
  block_type_reglage* reglage = lastReglage();
  block_type_param* param = lastParam();
  
  nb_coups= reglage->reglage.horloge.nb_mesures/2 - reglage->reglage.horloge.temp_mort;
  aa = (nb_coups<<14) + (nb_coups*190) ;
  
  int s = imesure % 2 ? 1 : -1;

  return s*(((blk->data_bolo[ibolo][imesure]-aa)<<1)/nb_coups);
}

def_gains 

double ArcheopsFile::getMuVBolo(int ibolo, int imesure) { // microvolts, filtre avec filtre carre
  double y = (getRawBolo(ibolo, imesure-1) - getRawBolo(ibolo, imesure))/2.;
  if (imesure%2) y=-y;
  block_type_reglage* reglage = lastReglage();
  return ((1e5*y)/(65536.*gain_ampli(reglage->reglage.bolo[ibolo])));
}


// SST, gyros...

void ArcheopsFile::needSSTProcessMask(int mask) {
  blockSet->sstHandler.NeedProcess(mask);
}


long ArcheopsFile::getSSTSignal(int idiode, int imesure) {
  return blockSet->sstHandler.getSignal(imesure, idiode);
}
// $CHECK$ TBD


double ArcheopsFile::getAzimut(int imesure) {return imesure*360/nEchBlock();}
double ArcheopsFile::getPendDirect(int /*imesure*/) {return 0;}
double ArcheopsFile::getPendOrth(int /*imesure*/) {return 0;}
double ArcheopsFile::getAlpha(int /*imesure*/) {return 0;}
double ArcheopsFile::getDelta(int /*imesure*/) {return 0;}