#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>

#include "../numeric/SU21-numeric.h"
#include "../numeric/SU21-numeric-base.c"
#include "../numeric/SU21-fundamental-domain.c"

int main(int argc,char *(argv[])) {

  // Input files
  FILE *FewEltsFile;

  char *Prefix_p,*NN_p,FileNam[MAX_FILENAME_LEN]="../";

  Elt_t *FewElts;
  int NumFewElts,NumPts;

  int NN,Ix,Th1Ix,Th2Ix,Th3Ix,PtIx;
  FLOAT_t Vol,VolCpt,Th1,Th2,Th3,r1,r2,hh0,hh1,hh2,hh3;
  FLOAT_t sMax,ss,VolFact,fTmp,DeltaVol,DeltaVolCpt;
  FLOAT_t Q,JacFact;
  Pt_t *Pts,Pt0;
  int *PtsCnt;
  FLOAT_t *FaceVol,*FaceVolCpt;

  // Get the prefix
  Prefix_p=getenv("PREFIX");
  assert(Prefix_p != NULL);

  // Open the files
  {
    size_t len=strlen(Prefix_p);
    assert(3+len+4<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix_p);
    // Input files
    strcpy(&FileNam[3]+len,"-FewElts");
    printf("%s\n",FileNam);
    FewEltsFile=fopen(FileNam,"r");
    assert(FewEltsFile != NULL);
    // Output files
  }

  // Get the elements
  GetEltsString(FewEltsFile,"NumFewElts",&NumFewElts,&FewElts);
  // copy the points into the points vector
  SetupPts(&Pts,&NumPts,FewElts,NumFewElts);
  printf("NumFewElts=%i NumPts=%i\n",NumFewElts,NumPts);
  fflush(stdout);
  free(FewElts);

  // Allocate memory
  PtsCnt=malloc(NumPts*sizeof(int));
  assert(PtsCnt != NULL);
  FaceVol=malloc(NumPts*sizeof(FLOAT_t));
  assert(FaceVol != NULL);
  FaceVolCpt=malloc(NumPts*sizeof(FLOAT_t));
  assert(FaceVolCpt != NULL);

  // Get NN
  NN_p=getenv("NN");
  assert(NN_p != NULL);
  assert(sscanf(NN_p,"%i",&NN) == 1);

  memset(PtsCnt,0,NumPts*sizeof(int));
  memset(FaceVol,0,NumPts*sizeof(FLOAT_t));
  memset(FaceVolCpt,0,NumPts*sizeof(FLOAT_t));
  sMax=0.0;
  Vol=0.0;
  VolCpt=0.0;
  hh0=((M_PI/2.0)/NN)*((2.0*M_PI)/NN)*((2.0*M_PI)/NN);

  /* Run through all possible directions on the 3-sphere. */
  Th3=0.0;
  for(Th3Ix=0;Th3Ix<=NN;Th3Ix++) {
    if(Th3Ix % 3 == 0) {
      printf("A: Th3Ix=%i Th3/(pi/2)="FOShortCnv" sMax="FOShortCnv"\n",
	     Th3Ix,Th3/(M_PI/2.0),sMax);
      fflush(stdout);
    }
    if(Th3Ix==0 || Th3Ix==NN)
      hh1=hh0*0.5;
    else
      hh1=hh0;
    r1=COS(Th3);
    r2=SIN(Th3);
    JacFact=r1*r2;
//Th1=0.0;  // Here we use symmetry under complex conjugation
    Th1=(-M_PI);
    for(Th1Ix=0;Th1Ix<=NN;Th1Ix++) {
      // printf("B: Th3Ix=%i Th3=%f Th1Ix=%i Th1=%f sMax=%f\n",
      //     Th3Ix,Th3,Th1Ix,Th1,sMax);
      if(Th1Ix==0 || Th1Ix==NN)
	hh2=hh1*0.5;
      else
	hh2=hh1;
      Pt0[0]=r1*(COS(Th1)+I*SIN(Th1));
      Th2=(-M_PI);
      for(Th2Ix=0;Th2Ix<=NN;Th2Ix++) {
	if(Th2Ix==0 || Th2Ix==NN)
	  hh3=hh2*0.5;
	else
	  hh3=hh2;
	Pt0[1]=r2*(COS(Th2)+I*SIN(Th2));
	ss=sVal(Pt0,Pts,NumPts,&PtIx);
	sMax=MAX(ss,sMax);
	fTmp=ss*ss/(1.0-ss*ss);
	DeltaVol=hh3*JacFact*fTmp*fTmp/4.0;
	DeltaVolCpt=hh3*JacFact/4.0;
	Vol+=DeltaVol;
	VolCpt+=DeltaVolCpt;
        PtsCnt[PtIx]++;
	FaceVol[PtIx]+=DeltaVol;
	FaceVolCpt[PtIx]+=DeltaVolCpt;

	Th2+=(2.0*M_PI)/NN;
      }
      Th1+=(2.0*M_PI)/NN;
    }
    Th3+=(M_PI/2.0)/NN;
  }

  // Vol*=2.0; // Because we've divided out by complex conjugation
  // VolCpt*=2.0;

  for(Ix=0;Ix<NumPts;Ix++) {
    if(PtsCnt[Ix]>0) {
      printf("Ix=%5i Cnt=%10i FaceVol="FOCnv" FaceVolCpt="FOCnv"\n"
	     "FaceVolNorm="FOCnv" FaceVolCptNorm="FOCnv" Pt=\n",
	     Ix,PtsCnt[Ix],FaceVol[Ix],FaceVolCpt[Ix],
	     2.0*FaceVol[Ix]/(M_PI*M_PI),
	     2.0*FaceVolCpt[Ix]/(M_PI*M_PI));
      OutputPt(stdout,Pts[Ix]);
    }
  }
    
  printf("\nNormalized VolCpt="FOCnv" sMax="FOCnv"\n",
  	 2.0*VolCpt/(M_PI*M_PI),sMax);
  printf("NN=%i Vol="FOCnv" Normalized Vol="FOCnv"\n",
  	 NN,Vol,2.0*Vol/(M_PI*M_PI));
}