[C] Problema calcolo coefficienti serie di Fourier

Ciao a tutti... devo realizzare questo programmino per il calcolo del valore dei vari coefficienti della serie di Fourier.. con la scelta fra varie funzioni da studiare.. è ancora incompleto, proprio perché mi trovo davanti a questo problema: usando un metodo di integrazione (avevo impostato Cavalieri simpson ma ho ridotto ai rettangoli per vedere se c'era qualcosa che non tornava) i coefficienti vengono calcolati con il valore che mi aspetto, ma con il segno sbagliato.. ovvero: sono tutti positivi quando invece dovrebbero avere anche segno negativo!! I valori in modulo però sembrano essere corretti... non capisco dove sia l'errore...
Se qualcuno ha tempo da perdere... beh, ecco il codice

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/********************************/
/* ###  PROTOTIPI FUNZIONI  ### */
/********************************/

/* PROTOTIPO ONDA QUADRA */
double ondaq(double t);
/* PROTOTIPO ONDA TRIANGOLARE */
double ondat(double t);
/* PROTOTIPO DENTE DI SEGA */
double dentes(double t);
/* PROTOTIPO FUNZIONE POLIGONALE */
double polig(double t);
/* PROTOTIPO COEFFICIENTE An*/
double An(double t, int n, double T, double (*f)(double));
/* PROTOTIPO COEFFICIENTE Bn*/
double Bn(double t, int n, double T, double (*f)(double));
/* PROTOTIPO METODO CAVALIERI-SIMPSON */
double fCS(double t, double IR, int N, int n, double T, double (*f)(double), double (*integranda)(double,int,double,double(*f)(double)));


/**********************************/ 
/* ###  PROGRAMMA PRINCIPALE  ### */
/**********************************/
int main()
{
  /* Inizializzo le variabili */
  int pari; /* parametro indicatore della parità o meno della funzione scelta */
  int fscelta; /* parametro indicatore della funzione scelta */
  int n=0, j=0; /* contatori dei cicli */
  int N=1000; /* numero di passi da effettuare nel calcolo degli integrali */
  int ordine; /* ordine a cui fermare lo sviluppo in serie */
  double T; /* periodo della funzione */
  double t=0; /* variabile indipendente (tempo) */
  double *a; /* puntatori ai vettori degli n coefficienti */
  double *b; /* della serie di Fourier */
  double funzione; /* variabile che contiene il valore della funzione f(t) */
  double IR=0; /* variabile che tiene conto del risultato dell'integrazione */
  double serie=0; /* variabile che tiene conto della somma della serie di Fourier */
  double (*f)(double); /* puntatore a funzione che punta ad una funzione specificata in seguito, */
  /* a seconda della scelta effettuata dall'utente sulla funzione da studiare */
  double (*integranda)(double,int,double,double);
  
  /* Chiedo quale funzione studiare */
  printf("\n\nChe funzione studiare? Digitare:\n1 = ONDA QUADRA\n2 = OND TRIANGOLARE\n3 = DENTE DI SEGA\n4 = FUNZIONE POLINOMIALE\n");
  scanf("%d", &fscelta);
  printf("\nFUNZIONE SCELTA: ");
  if (fscelta==1)
  {
    pari=1;
    printf("Onda quadra\n");
    f=ondaq;
    T=2.*M_PI;
    t=0.;
  } else if (fscelta==2) {
    pari=0;
    printf("Onda triangolare\n");
    f=ondat;
  } else if (fscelta==3) {
    pari=1;
    printf("Dente di sega\n");
    f=dentes;
  } else if (fscelta==4) {
    pari=2;
    printf("Funzione polinomiale\n");
    f=polig;
  }
  
  printf("\nA che ordine fermare la serie?\n");
  scanf("%d", &ordine);
  a = (double*) malloc(ordine*sizeof(double));
  b = (double*) malloc(ordine*sizeof(double));

  
  printf("\t------------------------------------------------------------------------\n");
  /*scanf("%lf", &t);*/
  funzione=(*f)(t);
  printf("%lf, %lf, %lf", funzione, (*f)(t), t);
  printf("\n\nARMONICA     An       Bn\n");
  for (n=0;n<ordine;n++)
  {
    a[n]=(2./T)*fCS(t,IR,N,n,T,f,An);
    b[n]=(2./T)*fCS(t,IR,N,n,T,f,Bn);
    printf("%d           %+1.4lf     %1.2lf\n", n+1, a[n],b[n]);
  }
  
  /*for (t=0;t<T;t+h)
  {
    printf
  }*/
  
  
  printf("\n\n");
  return 0;
}


/********************************************************/ 
/* ### FUNZIONI UTILIZZATE NEL PROGRAMMA PRINCIPALE ### */
/********************************************************/

/* IMPLEMENTAZIONE DELLA FUNZIONE INTEGRANDA PER IL CALCOLO DEL COEFFICIENTE An */
double An(double t, int n, double T, double (*f)(double))
{
  return ((*f)(t)*sin((n*M_PI*t)/(T/2.)));
}

/* IMPLEMENTAZIONE DELLA FUNZIONE INTEGRANDA PER IL CALCOLO DEL COEFFICIENTE Bn */
double Bn(double t, int n, double T, double (*f)(double))
{
  return ((*f)(t)*cos((n*M_PI*t)/(T/2.)));
}

/* IMPLEMENTAZIONE METODO CAVALIERI-SIMPSON 1/3 */
double fCS(double t, double IR, int N, int n, double T, double (*f)(double), double (*integranda)(double,int,double,double(*f)(double)))
{
  /*int u=0,k=0;
  double h=(double)(T/N); passo di integrazione
  double f0=(*integranda)(t,n,T,f);
  double fN=(*integranda)(T/2.,n,T,f);
  IR=(h/3.)*(f0+fN);
  for (u=1;u<(N/2.);u++)
  {    
    t=t+2*h;
    IR=IR+(2.*h/3.)*(*integranda)(t,n,T,f);
  }   
  t=h; CONTROLLA!!!
  for (k=1;k<(N/2.-1);k++)
  {
    t=t+2*h;
    IR=IR+(4.*h/3.)*(*integranda)(t,n,T,f);
  }
  return (IR);*/
  int u=0;    
  double h=(double)(T/N);
  t=0;
  for (u=0;u<N;u++)
  {
    IR=IR+h*(*integranda)(t,n,T,f);
    t=t+h;
  }
  return (IR);
}

/* ONDA QUADRA */
double ondaq(double t)
{
  double oq, T=2.*M_PI;
  if ((t>=0.) && (t<(T/2.))) {
    oq=1.;
  } else {
    oq=0.;
  }
  return oq;
}

/* ONDA TRIANGOLARE */
double ondat(double t)
{
  double oq, T;
  if (t>=0 && t<T/2.) {
    oq=1;
  } else {
    oq=0;
  }
  return oq;
}

/* DENTE DI SEGA */
double dentes(double t)
{
  double oq, T;
  if (t>=0 && t<T/2.) {
    oq=1;
  } else {
    oq=0;
  }
  return oq;
}

/* FUNZIONE POLIGONALE */
double polig(double t)
{
  double oq, T;
  if (t>=0 && t<T/2.) {
    oq=1;
  } else {
    oq=0;
  }
  return oq;
}