source: Sophya/trunk/AddOn/TAcq/mfits2spec.cc@ 3559

// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"

// include standard c/c++
#include <math.h>
#include <stdio.h>

#include <iostream>
#include <string>

#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
#include "tarrinit.h" 
#include "timestamp.h"
#include "fftpserver.h"
#include "fftwserver.h"

#include "FFTW/fftw3.h"

// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"


// include mini-fits lib 
#include "minifits.h"


inline r_4 Zmod2(complex<r_4> z) 
{ return (z.real()*z.real()+z.imag()*z.imag()); }

int main(int narg, char* arg[])
{
  if (narg < 2) {
    cout << " Usage: mfits2spec file1.fits [file2.fits ...] " << endl;
    return 1;
  }
  cout << " ---------- mfits2spec.cc Start ------------- " << endl;

  TArrayInitiator  _inia;

  int rc = 0;
  try {
        ResourceUsage resu;
        TVector<float> spectre;
        sa_size_t nzm = 0;   // Nb de spectres moyennes
    for(int ifile=1; ifile<narg; ifile++) {
      string ffname = arg[ifile];
// -------------- Lecture de bytes      
      cout << ifile <<"-Ouverture/lecture fichier " << ffname << endl;
      MiniFITSFile mff(ffname, MF_Read);
      cout << "... Type=" << mff.DataTypeToString() << " NAxis1=" << mff.NAxis1() 
           << " NAxis2=" << mff.NAxis2() << endl;
      if (mff.DataType() != MF_Byte) {
        cout << " PB : DataType!=MF_Byte --> skipping " << endl;
      }
      size_t sx = mff.NAxis1(); 
      size_t sy = mff.NAxis2();

      Byte* data = new Byte[sx];
      TVector<r_4> vx(sx);
      TVector< complex<r_4> > cfour;
      FFTPackServer ffts;
      
      for(int j=0; j<sy; j++) {
            mff.ReadB(data, sx, j*sx);
        // On convertit en float 
        for(sa_size_t ix=0; ix<vx.Size(); ix++)  vx(ix) = (r_4)(data[ix]);
        // On fait le FFT
        ffts.FFTForward(vx, cfour);
        if (!spectre.IsAllocated())  spectre.SetSize(cfour.Size());
        
        // On cumule les modules carres
        for(sa_size_t jf=1; jf<spectre.Size(); jf++) 
          spectre(jf) += Zmod2(cfour(jf)); 
        nzm++;

//      cout << " j=" << j << " data[0...3]=" << (int)data[0] << " , " 
//           << (int)data[1] << " , " <<  (int)data[2] << endl;
      }
      cout << "---- FIN lecture " << ffname << endl;      
    }   
    cout << "---- Ecriture spectre moyenne ds mspectre.ppf " << endl;      
    spectre /= (r_4)(nzm);
    POutPersist po("mspectre.ppf");
    po << spectre;
    cout << resu ; 
  }
  catch (MiniFITSException& exc) {
    cerr << " mfits2spec.cc catched MiniFITSException " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n mfits2spec.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
    rc = 78;
  }
  catch (...) {
    cerr << " mfits2spec.cc catched unknown (...) exception  " << endl; 
    rc = 79; 
  } 

  cout << ">>>> mfits2spec.cc ------- FIN ----------- RC=" << rc << endl;
  return rc;

}