//******************************************************************************
// Programme de simulation du comportement d'un dipole RC excité par un signal
// en créneaux
//
// Dominique Lefebvre - TangenteX.com  Décembre 2010
//******************************************************************************

#include <cstdlib>
#include <iostream>
#include <math.h>
#include <dislin_d.h>

using namespace std;

//******************************************************************************
// Définition des constantes
//******************************************************************************
#define N        1         // nombre d'équations du système différentiel
#define FEch     1000      // fréquence d'échantillonnage


//******************************************************************************
// Définition des variables globales
//******************************************************************************
double F,R,C,V0;
int K;

//******************************************************************************
// Définition des subroutines de service
//******************************************************************************
void InitGraphDislin(void)
{
    metafl("XWIN");
    disini();
    pagera();
    pagfll(255);
    color("black");
    winfnt("arial");
    
    return;     
}

//******************************************************************************
// Définition du système différentiel
// du/dt = (1/RC)*(-u + V)
//******************************************************************************
void Derivee(double Y[], double dY[])
{
    dY[0] = 1/(R*C)*(-Y[0] + K*V0);   
    return;
}

//******************************************************************************
// Routine RK4
//******************************************************************************
void RK4(double yin[], double yout[], double h) 
{ 
    int i; 
    double dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N]; 
  
    Derivee(yin, dydx);
    for (i = 0; i < N; i++) 
     { 
      k1[i] = h*dydx[i]; 
      yt[i] = yin[i] + 0.5*k1[i]; 
     }

    Derivee(yt, dydxt); 
    for (i = 0; i < N; i++) 
     { 
      k2[i] = h*dydxt[i]; 
      yt[i] = yin[i] + 0.5*k2[i]; 
     }

    Derivee(yt, dydxt); 
    for (i = 0; i < N; i++) 
     { 
      k3[i] = h*dydxt[i]; 
      yt[i] = yin[i] + k3[i]; 
     }  

    Derivee(yt, dydxt); 
    for (i = 0; i < N; i++) 
     { 
      k4[i] = h*dydxt[i]; 
      yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.; 
     }  
    return; 
}

//******************************************************************************
// Routines de traitement 
//******************************************************************************
void BaseTemps(const double P,const int NbP,const double dt,double T[])
{
    int i;
    double ti;
   
    T[0] = 0.0;
    for(i=1,ti=0.0;ti<P*NbP;ti+=dt,i++) T[i] = T[i-1]+dt;
    return;     
}


void SignalCarre(const double P, const int NbP,const double dt,double Vs[])
{
    int i, i0, j, k, NbPoints;
    double ti,tp;
       
    i0 = 0;
    for(k=0;k<NbP;k++)
      {
        for (ti=0.0+k*P,i=i0;ti<=(P/2)+k*P;ti+=dt,i++)Vs[i]=V0; 
        tp = ti;
        for (ti=tp,j=i;ti<=(k+1)*P;ti+=dt,j++)Vs[j]=-V0;
        i0 = j;         
      }
    return;     
}


//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
    double t0,ti, dt, Periode;
    double Yin[N], Yout[N], ytemp[N];
    double xmin, xmax, xpas,ymin, ymax, ypas;
    double *T, *Vs, *U;
    int i, j, i0, nbstep, NbPeriode;
    
    // En tête
    cout << "Reponse d'un dipole RC a un signal carre\n\n";
    
    // Saisie des paramètres
    cout << "Frequence du signal carre (en Hz): " << flush; cin >> F;   
    cout << "Amplitude max du signal carre (en volt): " << flush; cin >> V0;
    cout << "Resistance (en ohm) : " << flush; cin >> R;
    cout << "Capacite (en farad): " << flush; cin >> C; 
    cout << "Nombre de periodes a afficher: " << flush; cin >> NbPeriode;
    cout << "\n\n";
    cout << "La constante de temps du dipole est: " << R*C << " s\n";
    cout << "La periode du signal de charge est: " << 1/F << " s\n";   

    
    // Autres paramètres de simulation
    Periode = 1/F;
    dt = Periode*NbPeriode/FEch;
    nbstep = FEch;
       
    // Initialisation des variables de calcul
    T = new double [2*nbstep];
    Vs = new double [2*nbstep];
    U = new double [2*nbstep];
   
    // Conditions initiales de la simulation
    U[0] = -V0; 
    
    // Calcul du vecteur temps
    BaseTemps(Periode,NbPeriode,dt,T);
        
    // Calcul du signal carré entre 0 et NbPeriode 
    SignalCarre(Periode,NbPeriode,dt,Vs);   
    
    // boucle de calcul par RK4 pour une période
    i0 = 0;i = 0;  
    for(j=1;j<=NbPeriode;j++)
      {   
        ytemp[0] = U[i0];
        K = 1;  // pour la 1er demi periode - charge  
        for (ti=0.0;ti<=Periode/2;ti+=dt) 
         {
           Yin[0] = ytemp[0];
           RK4(Yin, Yout, dt);
           ytemp[0] = Yout[0];
           U[i++] = ytemp[0];
          } 
       ytemp[0] = U[i-1];
       K = -1;  // pour la seconde demi periode - décharge  
       for (ti=Periode/2;ti<=Periode;ti+=dt) 
         {
           Yin[0] = ytemp[0]; 
           RK4(Yin, Yout, dt);
           ytemp[0] = Yout[0];
           U[i++] = ytemp[0];
         }
       i0 = i-1;         
      }
        
    // Initialisation de la page graphique
    InitGraphDislin();
    wintit("Dipole RC");
    titlin("Reponse d'un dipole RC a un signal creneaux",1);
    name("Temps (s)","X");
    name("Tension (V)","Y");
    labdig(4,"X");
    xmin = T[0]; xmax = Periode*NbPeriode; xpas = xmax/5;
    ymin = -V0-0.1; ymax = V0+0.1; ypas = 1.0;
    graf(xmin, xmax, xmin,xpas,ymin, ymax, -V0, ypas);

    title();
    rline(0,0,xmax,0);
    color("blue");
    curve(T, Vs,nbstep);
    color("red");
    curve(T,U,nbstep);
    
    // Fin du programme
    disfin();
    delete [] T;
    delete [] Vs;
    delete [] U;
      
    return EXIT_SUCCESS;
}