#include <iostream>
#include "TArrayD.h"
#include "TMath.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TVirtualFFT.h"
#include "TWignerd.h"

void wig_points(TArrayD& arr,Int_t l,Int_t m,Int_t n);

void wig_graph(Int_t l,Int_t m,Int_t n)
{
   // Number of points that will be used to represent the function.
   Int_t sizet = 2*l+1;
   if(sizet<100) sizet = 100;
   const Int_t size = sizet;

   TArrayD arr(size);
   wig_points(arr,l,m,n);

   Double_t vecx[size];
   for(Int_t j=0;j<size;++j) 
      vecx[j] = 2*j*TMath::Pi()/(size);

   TGraph *gr = new TGraph(size,vecx,arr.GetArray());

   gr->SetTitle("Wigner d function");
   gr->SetLineColor(kSpring+3);
   gr->SetFillColor(kSpring+3);
   gr->SetLineWidth(3);
   gr->Draw("ac");

}

void wig_points(TArrayD& arr,Int_t l,Int_t m,Int_t n)
{
    // Fills arr with d^l_{mn}. Number of points sampled
    // will be equal to arr.GetSize(), therefore a big enough 
    // array must be provided. The performance of this macro 
    // could be enhanced using even/odd symmetries of the 
    // d-functions and sine and cosine transforms.
    
   Int_t size = arr.GetSize();
   if(size < (2*l+1)){
      std::cout << "Need bigger array." << std::endl;
      return;
   }

   if(TMath::Abs(m)>l || TMath::Abs(n)>l){
      std::cout << "Absolute value of m and n must be <= l" << std::endl;
      return;
   }

   TWignerd wig(l);
   wig.Advance(l);

   TVirtualFFT *B = TVirtualFFT::FFT(1,&size,"C2R M");
   Int_t half = size/2 + 1;

   // Set only positive frequencies

   for(Int_t u=0;u<half;++u){
      TComplex c = wig.FourierCoeff(m,n,u);
      B->SetPointComplex(u,c);
   }

   B->Transform();

   arr.Set(size,B->GetPointsReal());
   delete B;
}