source: Sophya/trunk/Poubelle/archTOI.old/formepulse.cc@ 3955

#include <iostream.h>
#include <fstream.h>

// D. Yvon, CE Saclay, DAPNIA/SPP, 08/99

#include <math.h>
#include "mayer_fft.h"
#include "dyffttools.h"


#include "transfelec.h"
#include "formepulse.h"

#include "ssthardware.h"

#define LTab (2048)
#define LSave (256)
#define Tpas (5.e-4)
//float PulseOptOrigin=LTab/16.*Tpas;
float PulseOptOrigin=0.;

//#define DYTestRun

/*
FormePulse::FormePulse() {
        gainElec=0.;
        FullwidthHalfMax=0;
        ShapeArr=new double[LSave];
        for(int i=0; i<LSave; i++)  ShapeArr[i]=0.;
}
*/

FormePulse::FormePulse(TransFuncElec TF) {

        Dcomplex j(0.,1.);
        double Fpas=1./(Tpas*LTab);
        
        double *PulseInReel=new double[LTab];
        double *PulseInImg=new double[LTab];
        double *PulseOut=new double[LTab];

        Dcomplex* PulseFFT=new Dcomplex[LTab];
        

        for(long i=0; i<LTab; i++) PulseInReel[i]=pulseOptique(i*Tpas);
        for(long i=0; i<LTab; i++) PulseInImg[i]=0.;
        
/*
//      On fenetre le signal.
        for(long i=0; i<LTab; i++) PulseInReel[i]=PulseInReel[i]*DYWindow(i,LTab);
*/

        ofstream dlout; // sert a sortir les print test
        
#ifdef DYTestRun
                //On ecrit le signal entrant
                char s[32]="WaveIn.txt";
                dlout.open(s,ios::out | ios::trunc);
                if (!dlout.is_open())
                {       cerr<<"erreur a l'ouverture du fichier de resultats: "<<s<<endl;
                        exit(-1);
                }
                dlout<<"temps"<<"\tPartieReelle"<<"\tPartieImag"<<endl;
                
                for (long k=0; k<LTab; k++) {
                        dlout<<k*Tpas<<"\t"<<PulseInReel[k]<<"\t"<<PulseInImg[k]<<endl;
                }
                dlout.close();

#endif

// C'est la fft!
        fft(LTab, PulseInReel, PulseInImg);


// On Filtre dans l'espace des frequences
        for (int i=0; i<LTab; i++) {
                PulseFFT[i]=PulseInReel[i]+PulseInImg[i]*j;
        //On applique le filtre Iici
                PulseFFT[i]=TF.TsFunc(2*3.1416*i*Fpas)*PulseFFT[i];
                PulseInReel[i]=real(PulseFFT[i]);
                PulseInImg[i]=imag(PulseFFT[i]);
        }
        
#ifdef DYTestRun        
        // On écrit la FFT      
                char s1[32]="PulseFFT.txt";
                dlout.open(s1,ios::out | ios::trunc);
                if (!dlout.is_open())
                {       cerr<<"erreur a l'ouverture du fichier de resultats: "<<s1<<endl;
                        exit(-1);
                }
                dlout<<"Frequence"<<"\tPartieReelle"<<"\tPartieImag"<<endl;
                
                for (long k=0; k<LTab; k++) {
                        dlout<<k*Fpas<<"\t"<<PulseInReel[k]<<"\t"<<PulseInImg[k]<<endl;
                }
                dlout.close();
#endif

// FFT inverse
        ifft(LTab,PulseInReel,PulseInImg);
        for (int i=0; i<LTab; i++) {
                PulseInReel[i]/=LTab;
                PulseInImg[i]/=LTab;
        }
        

// On recherche le maximum d'amplitude
        double Max=0.;          //En volt
        int MaxIndex=0;
        
        for(int i=0; i<LSave; i++) 
                if(PulseInReel[i]>Max) {
                        Max=PulseInReel[i];
                        MaxIndex=i;
                }
                
        gainElec=Max/SSTFPulseCourant;
                //En Volt/ampere        
                
// On calcule la LMH

        int HalfMaxInf;
                for(HalfMaxInf=MaxIndex;HalfMaxInf>0;HalfMaxInf--)
                        if(fabs(PulseInReel[HalfMaxInf])<(Max/2.)) break;
        int HalfMaxSup; 
                for(HalfMaxSup=MaxIndex;HalfMaxSup<LTab;HalfMaxSup++)
                        if(fabs(PulseInReel[HalfMaxSup])<(Max/2.)) break;
        
        
        FullwidthHalfMax=(HalfMaxSup-HalfMaxInf)*Tpas;
                
// On normalise à un Volt. Et on ecrit en memoire le resultat
        ShapeArr=new double[LSave];
        for(int i=0; i<LSave; i++)  ShapeArr[i]=PulseInReel[i]/Max;

//#ifdef DYTestRun      
                //On ecrit le resultat dans un fichier 
                char s2[32]="Resultatfiltrage.txt";
                dlout.open(s2,ios::out | ios::trunc);
        //      if (!dlout.is_open()){  
        //              cerr<<"erreur a l'ouverture du fichier de resultats: "<<s2<<endl;
        //              exit(-1);
        //      }

                dlout<<"Impulsion courant de: "<<SSTFPulseCourant<<" Ampères"<<endl;
                dlout<<"Maximum d'amplitude: "<<Max<< " Volts"<<endl;
                dlout<<" gain chaine preampli: "<<gainElec<<endl;
                dlout<<"largueur a mi hauteur: "<<FullwidthHalfMax<<endl<<endl;
                dlout<<"Temps"<<"\tPartie Reelle"<<endl;

                for (int k=0; k<LSave; k++) {
                        dlout<<k*Tpas<<"\t"<<ShapeArr[k]<<endl;
                }       
//#endif


        
//On tue les tableaux de passage....
        delete[] PulseInReel;
        delete[] PulseInImg;
        delete[] PulseOut;
        delete[] PulseFFT;

}

FormePulse& FormePulse::operator =(const FormePulse& FPlse) {
        ShapeArr=new double[LSave];
        for(int i=0; i<LSave; i++) ShapeArr[i]=FPlse.ShapeArr[i];
        return *this;
}

FormePulse::~FormePulse() {
        delete[] ShapeArr;
}

double FormePulse::TMax() {
        double returnVal=LSave*Tpas;
        return returnVal;
}

double FormePulse::PulseShape(double temps) {
        double returnValue;
        int ientier=(int)(temps/Tpas);
        if (ientier<0) returnValue=0.;
        else if (ientier>=LSave) returnValue= ShapeArr[LSave-1];
        else {
                double Min=ShapeArr[ientier];
                double Max=ShapeArr[ientier+1];
//              double test=temps/Tpas-ientier;
                returnValue=Min+(Max-Min)*(temps/Tpas-ientier);
        }
        return returnValue;
}


double FormePulse::LongPulseOptique() {
        double returnVal;
        returnVal=SSTLongPhotDiode/SSTLongFocale*PeriodRotTeles/2/3.1415/sin(SSTTelesElevation);
        //en secondes
        return returnVal;
//      return returnVal=30.e-3;                        //Gros BUGG!!!!!
}
        
double FormePulse::pulseOptique(double temps) {
        double value;
        if((temps>PulseOptOrigin)&&(temps<=PulseOptOrigin+LongPulseOptique()))  
                value=SSTFPulseCourant;                  //     Amplitude: 1 nanoampère

        else value=0.;
//      for(long i=0; i<LTab; i++) value=sin(2*3.1416*temps/5.e-3); //pour test FFT
        return value;
}