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

#include "a7p2-0-base.c"

#define NN (640)

int main(int argc,char *(argv[])) {
  int Ix,Th1Ix,Th2Ix,Th3Ix,NumPts,PtIx;
  double Vol,VolCpt,Th1,Th2,Th3,r1,r2,hh0,hh1,hh2,hh3,
    sMax,ss,tt,VolFact,fTmp,DeltaVol,DeltaVolCpt,
    Q,JacFact;
  Pt_t Pt0;
  int PtsCnt[MaxNumPts];
  double FaceVol[MaxNumPts],FaceVolCpt[MaxNumPts];

  MakeGpK();
  // PrintMatrices(GpK,GpKSiz,"GpK");
  for(Ix=0;Ix<GpKSiz;Ix++)assert(CheckUnitary(GpK[Ix],FF));

  PtsFile=fopen(argv[1],"r");
  assert(PtsFile != NULL);
  NumPts=ReadPts(Pts,MaxNumPts);

  // PrintPts(Pts,NumPts);

  memset(PtsCnt,0,sizeof(PtsCnt));
  memset(FaceVol,0,sizeof(FaceVol));
  memset(FaceVolCpt,0,sizeof(FaceVol));
  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)=%f sMax=%f\n",Th3Ix,
	     Th3/(M_PI/2.0),sMax);
    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.z[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.z[1]=r2*(cos(Th2)+I*sin(Th2));
	ss=sVal(Pt0,NumPts,&PtIx);
	sMax=MAX(ss,sMax);
	tt=ss/R0;
	fTmp=tt*tt/(1.0-tt*tt);
	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=%.15f FaceVolCpt=%.15f\n"
	     "FaceVolNorm=%.15f FaceVolCptNorm=%.15f 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));
      PrintPt(Pts[Ix]);
    }
  }
    
  printf("VolCpt=%.15f Normalized VolCpt=%.15f\n",
  	 VolCpt,2.0*VolCpt/(M_PI*M_PI),sMax);
  printf("Vol=%.15f Normalized Vol=%.15f sMax=%.15f\n",
  	 Vol,2.0*Vol/(M_PI*M_PI),sMax);
}