#include "source.h"

/* Parametri del circuito */
double SourceFrequency=1.E4;
double SourceR= 2.;
double SourceC= 1.E-6;
double SourceL= 1.E-3;
double SourceV= 1.;
double (*Vin)(double t);

/* Forzante sinusoidale */
double SinPulser(double t) {
  return SourceV*sin(2*M_PI*SourceFrequency*t);
}

/* Forzante quadra */
double SqrPulser(double t) {
  return SourceV*(double)(((int)(t*SourceFrequency))%2);
}

/* Equazione differenziale per circuito RC */
void RCdiffeq(double t,double vout[], double dvout[]) {
  dvout[0] = ((*Vin)(t)-vout[0])/(SourceR*SourceC);
  return;
}

/* Equazione differenziale per circuito LC */
void LCdiffeq(double t,double vout[], double dvout[]) {
  dvout[0] = vout[1];
  dvout[1] = ((*Vin)(t)-vout[0])/(SourceL*SourceC);
  return;
}

/* Equazione differenziale per circuito RLC */
void RLCdiffeq(double t,double vout[], double dvout[]) {
  dvout[0] = vout[1];
  dvout[1] = ((*Vin)(t)-vout[0])/(SourceL*SourceC)-vout[1]*(SourceR/SourceL);
  return;
}

/*
 * Seguono le funzioni che soddisfano le interfacce definite per
 * TDS1002Wave.
 * Queste definiscono con Vin il tipo di forzante, usano le variabili 
 * statich yold ed ynew per aver memoria dello stato dell'integrazione
 * e poi chiamano la funzione rkad di RungeKutta.c con la corretta equazione
 * differenziale per far avanzare l'integrazione di un passo temporale.
 */

void SimpleSinusoid(double* t0, double dummy1, double dummy2,
                    double  dt, double* sinusoid, double* cosinusoid) {
  (*t0)+=dt;
  *sinusoid   = sin(2*M_PI*SourceFrequency*(*t0));
  *cosinusoid = cos(2*M_PI*SourceFrequency*(*t0));
  return;
}

void RCSinusoid(    double* t0, double dummy1, double oldvcapacitance,
                    double  dt, double* pulser, double* vcapacitance) {
  double yold[1], ynew[1];
  Vin=SinPulser;
  yold[0]=oldvcapacitance;
  rkad(*t0,(*t0)+dt,1,yold,ynew,RKPRECISION,RCdiffeq);
  (*t0)+=dt;
  *pulser=(*Vin)(*t0);
  *vcapacitance=ynew[0];
  return;
}

void RCSquare(      double* t0, double dummy1, double oldvcapacitance,
                    double  dt, double* pulser, double* vcapacitance) {
  double yold[1], ynew[1];
  Vin=SqrPulser;
  yold[0]=oldvcapacitance;
  rkad(*t0,(*t0)+dt,1,yold,ynew,RKPRECISION,RCdiffeq);
  (*t0)+=dt;
  *pulser=(*Vin)(*t0);
  *vcapacitance=ynew[0];
  return;
}

void LCSinusoid(    double* t0, double dummy1, double oldvcapacitance,
                    double  dt, double* pulser, double* vcapacitance) {
  static double yold[2]={0.,0.};
  static double ynew[2]={0.,0.};
  Vin=SinPulser;
  yold[0]=oldvcapacitance;
  yold[1]=ynew[1];
  rkad(*t0,(*t0)+dt,2,yold,ynew,RKPRECISION,LCdiffeq);
  (*t0)+=dt;
  *pulser=(*Vin)(*t0);
  *vcapacitance=ynew[0];
  return;
}

void LCSquare(    double* t0, double dummy1, double oldvcapacitance,
                    double  dt, double* pulser, double* vcapacitance) {
  static double yold[2]={0.,0.};
  static double ynew[2]={0.,0.};
  Vin=SqrPulser;
  yold[0]=oldvcapacitance;
  yold[1]=ynew[1];
  rkad(*t0,(*t0)+dt,2,yold,ynew,RKPRECISION,LCdiffeq);
  (*t0)+=dt;
  *pulser=(*Vin)(*t0);
  *vcapacitance=ynew[0];
  return;
}

void RLCSinusoid(    double* t0, double dummy1, double oldvcapacitance,
                    double  dt, double* pulser, double* vcapacitance) {
  static double yold[2]={0.,0.};
  static double ynew[2]={0.,0.};
  Vin=SinPulser;
  yold[0]=oldvcapacitance;
  yold[1]=ynew[1];
  rkad(*t0,(*t0)+dt,2,yold,ynew,RKPRECISION,RLCdiffeq);
  (*t0)+=dt;
  *pulser=(*Vin)(*t0);
  *vcapacitance=ynew[0];
  return;
}

void RLCSquare(    double* t0, double dummy1, double oldvcapacitance,
                    double  dt, double* pulser, double* vcapacitance) {
  static double yold[2]={0.,0.};
  static double ynew[2]={0.,0.};
  Vin=SqrPulser;
  yold[0]=oldvcapacitance;
  yold[1]=ynew[1];
  rkad(*t0,(*t0)+dt,2,yold,ynew,RKPRECISION,RLCdiffeq);
  (*t0)+=dt;
  *pulser=(*Vin)(*t0);
  *vcapacitance=ynew[0];
  return;
}


