#ifndef FUNDAMENTAL_DOMAIN_C
#define FUNDAMENTAL_DOMAIN_C

#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"

FLOAT_t sVal(Pt_t Pt0,Pt_t *Pts,int NumPts,int *PtIx_p);
FLOAT_t sValOnePt(Pt_t Pt0,Pt_t z);
int InDirichlet(Pt_t P,Pt_t *Pts,int NumPts);
inline void CheckGeom1Simp(Geom1Simp_t GSimp);
inline void CheckGeom2Simp(Geom2Simp_t GSimp);
inline void CheckGeom3Simp(Geom3Simp_t GSimp);
inline FLOAT_t InRadius1(Geom1Simp_t GSimp);
FLOAT_t InRadius2(Geom2Simp_t GSimp);
FLOAT_t InRadius3(Geom3Simp_t GSimp);
FLOAT_t Eccentricity(Geom3Simp_t G3Simp);
FLOAT_t MaxEdgeLen(Geom3Simp_t G3Simp);
void ReduceGeom3Simp(Geom3Simp_t *NewG3Simp_p,Geom3Simp_t G3Simp);
void SplitGeom3Simp(Geom3Simp_t *NewG3Simps_p,Geom3Simp_t G3Simp);

// Given:
//
//   a unit vector, Pt0, in $\CC_2$ and
//   a table, Pts, of points in the unit ball of $\CC_2$
//
// this routine identifies the largest positive real multiple of Pt0, which
// is hyperbolically as close to the origin as to any of the points of Pts.
//
// The returned value, sVal, is the necessary positive real multiplier of
// Pt0.  If there is no upper bound, that is if all positive real multiples
// of Pt0 lying in the unit ball are hyperbolically closer to the origin
// than to any of the points of Pts, then the returned value is HUGE_FLOAT.
//
// In addition *PtIx_p is set to the index of the point in Pts which
// fixes the bound for the positive real multiplier.  If no such point
// has been found, *PtIx_p is set to -1.
//
// Normally the origin itself is one of the elements of Pts, and the routine
// deals with this gracefully.
//
// Recall that QQ(Pt1,Pt2) gives a quantity which depends monotonically on
// the hyperbolic distance, and that QQ1(Pt1,Pt2) is the same, but missing
// a factor depending only on Pt2.

FLOAT_t sVal(Pt_t Pt0,Pt_t *Pts,int NumPts,int *PtIx_p) {
  int Ix,Ix0;
  FLOAT_t complex u1;
  FLOAT_t R,limit,real1,a1,a2,zAbsSq,test,s,s0,QQSav,QQ0,QQz;
  Pt_t LPt,ZeroPt;


  // Make sure that Pt0 is a unit vector
  assert(CABS(Scal(Pt0,Pt0)-1.0)<epsBig);

  ZeroPt[0]=ZeroPt[1]=0.0;   // Zero vector
  LPt[0]=LPt[1]=0.0; // Positive multiple of Pt0
  s0=HUGE_FLOAT;    // Multiplier to be returned
  Ix0=(-1);   // Index of the point of Pts which fixes the multiplier
  QQSav=HUGE_FLOAT; // Value to compare to QQ(Pts[Ix],LPt).
  limit=HUGE_FLOAT; // Value to compare to QQ1(Pts[Ix],LPt)

  for(Ix=0;Ix<NumPts;Ix++) {
    // printf("sVal 100: Pt0="); PrintPt(Pt0);
    // printf("          z="); PrintPt(Pts[Ix]);

    // Is LPt closer to Pts[Ix] than to the origin?
    if(QQ1(Pts[Ix],LPt)<limit) {
      assert(QQ(Pts[Ix],LPt)-QQSav < epsBig*QQSav);
      u1=Scal(Pt0,Pts[Ix]);
      // printf("sVal 200: u1=% .15f%+.15f\n",u1);
      real1=CREAL(u1);
      // printf("sVal 300: real1=% .15f\n",real1);

      // Do there exist positive multiples of Pt0 which are closer to
      // Pts[Ix] than to the origin?  This can fail only if _limit_
      // still has its initial value of HUGE_FLOAT.
      if(real1>0) {
	a1=CABS(u1);
	// printf("sVal 400:  a1=% .15f\n",a1);
	zAbsSq=PtAbsSq(Pts[Ix]);
	// printf("sVal 500: zAbsSq=% .15f\n",zAbsSq);

	test=real1*real1-(a1*a1*zAbsSq);
	// printf("sVal 600: test=% .15f\n",test);

	// If test<=0, then there aren't any multiples of Pt0 which
	// are closer to Pts[Ix] than to the origin?  This can happen
	// only if _limit_ still has its initial value of HUGE_FLOAT.
	if(test>0) {
	  // Calculate the positive multiplier s0 such that s0*Pt0
	  // is equidistant from 0 and Pts[Ix]
	  s=zAbsSq/(real1+SQRT(test));
	  // printf("sVal 700: Pt0="); PrintPt(Pt0);
	  // printf("              z="); PrintPt(Pts[Ix]);
	  // printf("              s=% .15f\n",s);

	  // Make sure the point found is in the open unit ball.  This
	  // test can fail only if _limit_ still has its initial value
	  // of HUGE_FLOAT.
	  if(s<1.0) {
	    // Here we have found a new value for s0.
	    Ix0=Ix;
	    // We should have s<s0, but round-off errors in the above
	    // calculations may introduce a tiny error.
	    assert(s<s0+epsSmall);
	    s0=s;
	    // Calculate s*Pt0
	    LPt[0]=s*Pt0[0];
	    LPt[1]=s*Pt0[1];
	    // Verify that d(0,s*Pt0) = d(Pts[Ix],s*Pt0)
	    QQ0=QQ(ZeroPt,LPt);
            QQz=QQ(Pts[Ix],LPt);
	    assert(fabs((QQz-QQ0)/QQ0)<epsBig);
	    // printf("QQSav=%15.10f QQ=%15.10f Ix0=%4i\n",
	    // QQSav,QQ(ZeroPt,LPt),Ix0);
	    // To compare with new values of QQ(Pts[Ix],LPt)
	    QQSav=QQ0;
	    // To compare with new values of QQ1(Pts[Ix],LPt)
	    limit=QQ1(Pts[Ix],LPt);
	  }  
	}
      }
    }
  }
  // assert(Ix0 != -1);
  
  // Diagnostics for the case that no point Pts[Ix] has been found
  // which sets an upper bound to the multiple of Pt0 which is as at
  // least as close to the origin as it is to Pts[Ix].
  if(Ix0 == -1) {
    printf("Bad Pt0: ");
    OutputPt(stdout,Pt0);
    printf("  args=(" FOSpCnv "," FOSpCnv ") "
	   "abs=(" FOCnv "," FOCnv ")\n",
   	   CARG(Pt0[0])/M_PI,CARG(Pt0[1])/M_PI,
   	   CABS(Pt0[0]),CABS(Pt0[1]));
  }
  *PtIx_p=Ix0;
  // printf("Ix0=%i s0=%f\n",Ix0,s0);
  return s0;
}

// Given:
//
//   a unit vector, Pt0, in $\CC_2$ and
//   a point z in the unit ball of $\CC_2$
//
// this routine identifies a postivie real multiple of Pt0
// which is hyperbolically equidistant between the origin and z.
//
// The returned value, sVal, is the necessary positive real multiplier
// of Pt0.  If there is no such positive real multiplier which works,
// then the returned value is HUGE_FLOAT.

FLOAT_t sValOnePt(Pt_t Pt0,Pt_t z) {
  FLOAT_t real1,s,a1,zAbsSq,test;
  COMPLEX_t u1;
  
  u1=Scal(Pt0,z);
  // printf("sVal 200: u1=% .15f%+.15f\n",u1);
  real1=CREAL(u1);
  if(real1 <= 0)return HUGE_FLOAT;
  a1=CABS(u1);
  zAbsSq=PtAbsSq(z);
  test=real1*real1-(a1*a1*zAbsSq);
  if(test <= 0)return HUGE_FLOAT;
  s=zAbsSq/(real1+SQRT(test));
  if(s>=1)return HUGE_FLOAT;
  return s;
}

// Checks whether a given point of the complex unit ball lies in the
// Dirichlet fundamental domain.
int InDirichlet(Pt_t P,Pt_t *Pts,int NumPts) {
  Pt_t P0;
  double Norm;
  int dummy;
  
  Norm=SQRT(PtAbsSq(P));
  if(Norm<epsBig)return TRUE;

  // Divide by the norm and obtain a unit vector
  P0[0]=P[0]/Norm;
  P0[1]=P[1]/Norm;
  assert(FABS(PtAbsSq(P0)-1.0)<epsSmall);
  
  return Norm < sVal(P0,Pts,NumPts,&dummy)+epsBig;
}

// Checks that the 2 vertices of a 1-simplex in~$\CC^2$ both lie on the
// unit sphere.
inline void CheckGeom1Simp(Geom1Simp_t GSimp) {
  Pt_t Vertex;
  int VIx;
  assert(FABS(PtAbsSq(GSimp.Vertex1) - 1.0) < epsSmall);
  AddPts(Vertex,GSimp.Vertex1,GSimp.Diff0);
  assert(FABS(PtAbsSq(Vertex) - 1.0) < epsSmall);
  assert(FABS(RealScal(GSimp.Diff0,GSimp.Vertex1) 
	      / PtAbsSq(GSimp.Diff0)
	      + 0.5 < epsBig));

}

// Checks that the 3 vertices of a 2-simplex in~$\CC^2$ all lie on the
// unit sphere.
inline void CheckGeom2Simp(Geom2Simp_t GSimp) {
  Pt_t Vertex;
  int VIx;
  assert(FABS(PtAbsSq(GSimp.Vertex2) - 1.0) < epsSmall);
  for(VIx=0;VIx<2;VIx++) {
    AddPts(Vertex,GSimp.Vertex2,GSimp.Diff[VIx]);
    assert(FABS(PtAbsSq(Vertex) - 1.0) < epsSmall);
    assert(FABS(RealScal(GSimp.Diff[VIx],GSimp.Vertex2) 
		/ PtAbsSq(GSimp.Diff[VIx])
		+ 0.5 < epsBig));
  }
}

// Checks that the 4 vertices of a 3-simplex in~$\CC^2$ all lie on the
// unit sphere.
inline void CheckGeom3Simp(Geom3Simp_t GSimp) {
  Pt_t Vertex;
  int VIx;
  FLOAT_t LenDiff,LenDiff2;
  assert(FABS(PtAbsSq(GSimp.Vertex3) - 1.0) < epsSmall);
  for(VIx=0;VIx<3;VIx++) {
    AddPts(Vertex,GSimp.Vertex3,GSimp.Diff[VIx]);
    assert(FABS(PtAbsSq(Vertex) - 1.0) < epsSmall);
    LenDiff2=PtAbsSq(GSimp.Diff[VIx]);
    LenDiff=SQRT(LenDiff2);
    if(LenDiff2>epsSmall)
      assert((FABS(RealScal(GSimp.Diff[VIx],GSimp.Vertex3)/LenDiff2
		   + 0.5 < epsBig/LenDiff)));
  }
}

// Calculates the distance of a 1-simplex in~$\CC^2$ to the origin.
// This calculation works under the hypothesis that both vertices
// lie on the unit sphere.
FLOAT_t InRadius1(Geom1Simp_t GSimp) {
  Pt_t Nearest;
  Geom1Simp_t OneSimp;
  FLOAT_t Sum;

  CheckGeom1Simp(GSimp);

  // Given that the two vertices lie on the unit sphere, the nearest
  // point to the origin is their midpoint.
  Nearest[0]=GSimp.Vertex1[0]+0.5*GSimp.Diff0[0];
  Nearest[1]=GSimp.Vertex1[1]+0.5*GSimp.Diff0[1];
  Sum=1.0-0.25*PtAbsSq(GSimp.Diff0);
  assert(FABS(Sum - PtAbsSq(Nearest)) < epsSmall);
  return SQRT(PtAbsSq(Nearest));
}

// Calculates the distance of a 2-simplex in~$\CC^2$ to the origin.
// This calculation works under the hypothesis that all the vertices
// lie on the unit sphere.
FLOAT_t InRadius2(Geom2Simp_t GSimp) {
  FLOAT_t Coeff[2],M[2][2],C[2],DetM;
  Pt_t Nearest;
  Geom1Simp_t OneSimp;

  CheckGeom2Simp(GSimp);

  // We find the point of the linear 2-space containing the 2-simplex
  // which is closest to the origin.  This will be
  //
  //   Nearest = Vertex2 + \sum_j Coeff[j]*Diff[j]
  //
  // Since this must be perpendicular to Diff[k] for each k
  // we have 
  //
  //   (Vertex2,Diff[k]) + \sum_j Coeff[j] (Diff[j],Diff[k]) = 0
  //
  // From |Vertex2|^2 = |Vertex2+Diff[k]|^2 one sees:
  //
  //   2 (Vertex3,Diff[k]) + (Diff[k],Diff[k]) = 0 
  //
  // whence
  //
  //    \sum_j Coeff[j] (Diff[j],Diff[k]) = 1/2 (Diff[k],Diff[k])

  M[0][0]=PtAbsSq(GSimp.Diff[0]);
  M[1][1]=PtAbsSq(GSimp.Diff[1]);
  M[0][1]=M[1][0]=RealScal(GSimp.Diff[0],GSimp.Diff[1]);

  DetM=M[0][0]*M[1][1]-M[0][1]*M[1][0];
  assert(DetM>0);

  C[0]=0.5*PtAbsSq(GSimp.Diff[0]);
  C[1]=0.5*PtAbsSq(GSimp.Diff[1]);

  Coeff[0]=(M[1][1]*C[0]-M[0][1]*C[1])/DetM;
  Coeff[1]=(-M[1][0]*C[0]+M[0][0]*C[1])/DetM;

  memcpy(Nearest,GSimp.Vertex2,sizeof(Pt_t));
  Nearest[0]+=
    (Coeff[0]*GSimp.Diff[0][0]+
     Coeff[1]*GSimp.Diff[1][0]);
  Nearest[1]+=
    (Coeff[0]*GSimp.Diff[0][1]+
     Coeff[1]*GSimp.Diff[1][1]);
     
  assert(FABS(RealScal(Nearest,GSimp.Diff[0]))
	 < epsBig*SQRT(PtAbsSq(GSimp.Diff[0])));
  assert(FABS(RealScal(Nearest,GSimp.Diff[1]))
	 < epsBig*SQRT(PtAbsSq(GSimp.Diff[1])));
     
  // This nearest point is in fact the circumcenter of the simplex.
  // Points in the linear 2-space containing the simplex are closer to
  // the origin when they are closer to this circumcenter.  If the
  // circumcenter lies within the simplex we are done.  Otherwise we
  // have to find the point of the simplex which lies nearest the
  // circumcenter.  If, for any of the bounding lines of the simplex,
  // the circumcenter and the simplex lie on opposite sides, then the
  // point we are looking for must lie on that  bounding line.
     
  if(Coeff[0]<0) {  // Circumcenter opposite the line through
                    // vertices 1 and 2
    memcpy(OneSimp.Vertex1,GSimp.Vertex2,sizeof(Pt_t));
    memcpy(OneSimp.Diff0,GSimp.Diff[1],sizeof(Pt_t));
  }
  else if(Coeff[1]<0) { // Circumcenter opposite the line
                        // vertices 0, 2
    memcpy(OneSimp.Vertex1,GSimp.Vertex2,sizeof(Pt_t));
    memcpy(OneSimp.Diff0,GSimp.Diff[0],sizeof(Pt_t));
  }
  else if(Coeff[0]+Coeff[1]>1.0) {
    // Circumcenter opposite the line through
    // vertices 0, 1
    AddPts(OneSimp.Vertex1,GSimp.Vertex2,GSimp.Diff[1]);
    SubPts(OneSimp.Diff0,GSimp.Diff[0],GSimp.Diff[1]);
  }
  else
    return SQRT(PtAbsSq(Nearest)); // Circumcenter inside the simplex

  return InRadius1(OneSimp);
}     

// Calculates the distance of a 3-simplex in~$\CC^2$ to the origin.
// This calculation works under the hypothesis that all the vertices
// lie on the unit sphere.
FLOAT_t InRadius3(Geom3Simp_t GSimp) {
  FLOAT_t Coeff[3],M[3][3],C[3];
  Pt_t Nearest;
  Geom2Simp_t TwoSimp;

  CheckGeom3Simp(GSimp);

  // We find the point of the linear 3-space containing the 3-simplex
  // which is closest to the origin.  This will be
  //
  //   Nearest = Vertex3 + \sum_j Coeff[j]*Diff[j]
  //
  // Since this must be perpendicular to Diff[k] for each k
  // we have 
  //
  //   (Vertex3,Diff[k]) + \sum_j Coeff[j] (Diff[j],Diff[k]) = 0
  //
  // From |Vertex3|^2 = |Vertex3+Diff[k]|^2 one sees:
  //
  //   2 (Vertex3,Diff[k]) + (Diff[k],Diff[k]) = 0 
  //
  // whence
  //
  //    \sum_j Coeff[j] (Diff[j],Diff[k]) = 1/2 (Diff[k],Diff[k])

  M[0][0]=PtAbsSq(GSimp.Diff[0]);
  M[1][1]=PtAbsSq(GSimp.Diff[1]);
  M[2][2]=PtAbsSq(GSimp.Diff[2]);
  M[0][1]=M[1][0]=RealScal(GSimp.Diff[0],GSimp.Diff[1]);
  M[1][2]=M[2][1]=RealScal(GSimp.Diff[1],GSimp.Diff[2]);
  M[2][0]=M[0][2]=RealScal(GSimp.Diff[2],GSimp.Diff[0]);

  assert(Det3By3Real(M)+epsBig >= 0);

  C[0]=0.5*PtAbsSq(GSimp.Diff[0]);
  C[1]=0.5*PtAbsSq(GSimp.Diff[1]);
  C[2]=0.5*PtAbsSq(GSimp.Diff[2]);

  Solve3By3Real(Coeff,M,C);

  memcpy(Nearest,GSimp.Vertex3,sizeof(Pt_t));
  Nearest[0]+=
    (Coeff[0]*GSimp.Diff[0][0]+
     Coeff[1]*GSimp.Diff[1][0]+
     Coeff[2]*GSimp.Diff[2][0]);
  Nearest[1]+=
    (Coeff[0]*GSimp.Diff[0][1]+
     Coeff[1]*GSimp.Diff[1][1]+
     Coeff[2]*GSimp.Diff[2][1]);

  assert(FABS(RealScal(Nearest,GSimp.Diff[0]))
	 < epsBig*SQRT(PtAbsSq(GSimp.Diff[0])));
  assert(FABS(RealScal(Nearest,GSimp.Diff[1]))
	 < epsBig*SQRT(PtAbsSq(GSimp.Diff[1])));
  assert(FABS(RealScal(Nearest,GSimp.Diff[2]))
	 < epsBig*SQRT(PtAbsSq(GSimp.Diff[2])));

  // This nearest point is in fact the circumcenter of the simplex.
  // Points in the linear 3-space containing the simplex are closer to
  // the origin when they are closer to this circumcenter.  If the
  // circumcenter lies within the simplex we are done.  Otherwise we
  // have to find the point of the simplex which lies nearest the
  // circumcenter.  If, for any of the bounding planes of the simplex,
  // the circumcenter and the simplex lie on opposite sides, then the
  // point we are looking for must lie on that  bounding plane.

  if(Coeff[0]<0) {  // Circumcenter opposite the plane through
                    // vertices 1, 2, and 3
    memcpy(TwoSimp.Vertex2,GSimp.Vertex3,sizeof(Pt_t));
    memcpy(TwoSimp.Diff[0],GSimp.Diff[1],sizeof(Pt_t));
    memcpy(TwoSimp.Diff[1],GSimp.Diff[2],sizeof(Pt_t));
  }
  else if(Coeff[1]<0) { // Circumcenter opposite the plane through
                        // vertices 0, 2, and 3
    memcpy(TwoSimp.Vertex2,GSimp.Vertex3,sizeof(Pt_t));
    memcpy(TwoSimp.Diff[0],GSimp.Diff[0],sizeof(Pt_t));
    memcpy(TwoSimp.Diff[1],GSimp.Diff[2],sizeof(Pt_t));
  }
  else if(Coeff[2]<0) { // Circumcenter opposite the plane through
                        // vertices 0, 1, and 3
    memcpy(TwoSimp.Vertex2,GSimp.Vertex3,sizeof(Pt_t));
    memcpy(TwoSimp.Diff[0],GSimp.Diff[0],sizeof(Pt_t));
    memcpy(TwoSimp.Diff[1],GSimp.Diff[1],sizeof(Pt_t));
  }
  else if(Coeff[0]+Coeff[1]+Coeff[2]>1.0) {
    // Circumcenter opposite the plane through
    // vertices 0, 1, and 2
    AddPts(TwoSimp.Vertex2,GSimp.Vertex3,GSimp.Diff[2]);
    SubPts(TwoSimp.Diff[0],GSimp.Diff[0],GSimp.Diff[2]);
    SubPts(TwoSimp.Diff[1],GSimp.Diff[1],GSimp.Diff[2]);
  }
  else
    return SQRT(PtAbsSq(Nearest));

  return InRadius2(TwoSimp);
}

// Calculate the eccentricity of a 3-simplex, defined as the ratio of
// squares of the longest and shortest edges.
FLOAT_t Eccentricity(Geom3Simp_t G3Simp) {
  FLOAT_t MaxSq,MinSq,SideSq;
  int VIx1,VIx2;
  
  MaxSq=0.0;
  MinSq=HUGE_FLOAT;
  for(VIx1=0;VIx1<3;VIx1++) {
    SideSq=PtAbsSq(G3Simp.Diff[VIx1]);
    if(SideSq>MaxSq)MaxSq=SideSq;
    if(SideSq<MinSq)MinSq=SideSq;

    for(VIx2=VIx1+1;VIx2<3;VIx2++){
      SideSq=EucDist2(G3Simp.Diff[VIx1],G3Simp.Diff[VIx2]);
      if(SideSq>MaxSq)MaxSq=SideSq;
      if(SideSq<MinSq)MinSq=SideSq;
    }
  }
  return MaxSq/MinSq;
}

// Calculates the maximum edge length of a 3-simplex
FLOAT_t MaxEdgeLen(Geom3Simp_t G3Simp) {
  FLOAT_t SideSq,MaxSq;
  int VIx1,VIx2;
  
  MaxSq=0.0;
  for(VIx1=0;VIx1<3;VIx1++) {
    SideSq=PtAbsSq(G3Simp.Diff[VIx1]);
    if(SideSq>MaxSq)MaxSq=SideSq;

    for(VIx2=VIx1+1;VIx2<3;VIx2++){
      SideSq=EucDist2(G3Simp.Diff[VIx1],G3Simp.Diff[VIx2]);
      if(SideSq>MaxSq)MaxSq=SideSq;
    }
  }
  return SQRT(MaxSq);
}

// Load sixteen simplexes, each of whose vertices is a unit vector,
// which exactly cover all directions to infinity
// in~$B(\CC^2)$ into a vector of geometric 3-simplexes
void Load16Geom3Simps(Geom3Simp_t *NewG3Simps_p) {
  int z1ReSgn,z1ImSgn,z2ReSgn,z2ImSgn,CIx;
  Pt_t V0,V1,V2,V3;

  for(z1ReSgn=1;z1ReSgn>=-1;z1ReSgn-=2)
    for(z1ImSgn=1;z1ImSgn>=-1;z1ImSgn-=2)
      for(z2ReSgn=1;z2ReSgn>=-1;z2ReSgn-=2)
	for(z2ImSgn=1;z2ImSgn>=-1;z2ImSgn-=2) {
	  V0[0] = (COMPLEX_t) z1ReSgn;
	  V0[1] = 0;
	  V1[0] = z1ImSgn*I;
	  V1[1] = 0;
	  V2[0] = 0;
	  V2[1] = (COMPLEX_t) z2ReSgn;
	  V3[0] = 0;
	  V3[1] = z2ImSgn*I;

	  memcpy(NewG3Simps_p->Vertex3,V3,sizeof(Pt_t));
	  SubPts(NewG3Simps_p->Diff[0],V0,V3);
	  SubPts(NewG3Simps_p->Diff[1],V1,V3);
	  SubPts(NewG3Simps_p->Diff[2],V2,V3);

	  NewG3Simps_p++;
	}
}

// Make a copy of a geometric three-simplex centered where the
// original was and somewhat smaller.
void ReduceGeom3Simp(Geom3Simp_t *NewG3Simp_p,Geom3Simp_t G3Simp) {
  const FLOAT_t Factor=0.75;
  Pt_t Displacement,Displacement3,MidPtDisp,MidPtDiff;
  int VIx;

  // Calculate the difference vector for the midpoint of 
  // the simplex;
  MidPtDiff[0]=
    0.25*((G3Simp.Diff)[0][0]+
	  (G3Simp.Diff)[1][0]+
	  (G3Simp.Diff)[2][0]);
  MidPtDiff[1]=
    0.25*((G3Simp.Diff)[0][1]+
	  (G3Simp.Diff)[1][1]+
	  (G3Simp.Diff)[2][1]);
  // normalize it so that it corresponds to a unit vector;
  NormalizeDiff(MidPtDiff,G3Simp.Vertex3);
  // multiply it by (1-Factor);
  memcpy(Displacement3,MidPtDiff,sizeof(Pt_t));
  Displacement3[0]*=(1.0-Factor);
  Displacement3[1]*=(1.0-Factor);
  // normalize the product so that it corresponds to a unit vector;
  NormalizeDiff(Displacement3,G3Simp.Vertex3);
  // and add the result to the old Vertex3 to get the new Vertex3
  AddPts(NewG3Simp_p->Vertex3,G3Simp.Vertex3,Displacement3);
  
  // For each of the first 3 vertices
  for(VIx=0;VIx<3;VIx++) {
    // adjust its difference vector to be with respect to the midpoint
    // of the simplex;
    SubPts(MidPtDisp,(G3Simp.Diff)[VIx],MidPtDiff);
    // multiply that adjustment by Factor;
    MidPtDisp[0]*=Factor;
    MidPtDisp[1]*=Factor;
    // unadjust the resulting product so that to be based at the 
    // old Vertex3 again;
    AddPts(Displacement,MidPtDisp,MidPtDiff);
    // normalize the result
    NormalizeDiff(Displacement,G3Simp.Vertex3);
    // and adjust it one more time to be with respect to the new
    // Vertex3.
    SubPts((NewG3Simp_p->Diff)[VIx],Displacement,Displacement3);
  }
}

// Slits a geometric three-simplex into eight smaller simplexes, and
// put these onto a vector list of geometic three-simplexes
void SplitGeom3Simp(Geom3Simp_t *NewG3Simps_p,Geom3Simp_t G3Simp) {
  Pt_t DMidPt01,DMidPt02,DMidPt03,DMidPt12,DMidPt13,DMidPt23;

  // Construct the midpoints of the edges: first in the form
  // of doubled-difference vectors;
  AddPts(DMidPt01,G3Simp.Diff[0],G3Simp.Diff[1]); // V0 and V1
  AddPts(DMidPt02,G3Simp.Diff[0],G3Simp.Diff[2]); // V0 and V2
  memcpy(DMidPt03,G3Simp.Diff[0],sizeof(Pt_t));   // V0 and V3
  AddPts(DMidPt12,G3Simp.Diff[1],G3Simp.Diff[2]); // V1 and V2
  memcpy(DMidPt13,G3Simp.Diff[1],sizeof(Pt_t));   // V1 and V3
  memcpy(DMidPt23,G3Simp.Diff[2],sizeof(Pt_t));   // V2 and V3
  // and then in the form of difference vectors;
  HalvePt(DMidPt01); HalvePt(DMidPt02); HalvePt(DMidPt03);
  HalvePt(DMidPt12); HalvePt(DMidPt13); HalvePt(DMidPt23);
  // Normalize the midpoint difference vectors so that the points
  // vectors are unit vectors in~$\CC^2$.
  NormalizeDiff(DMidPt01,G3Simp.Vertex3);
  NormalizeDiff(DMidPt02,G3Simp.Vertex3);
  NormalizeDiff(DMidPt03,G3Simp.Vertex3);
  NormalizeDiff(DMidPt12,G3Simp.Vertex3);
  NormalizeDiff(DMidPt13,G3Simp.Vertex3);
  NormalizeDiff(DMidPt23,G3Simp.Vertex3); 
  
  // Construct the eight simplexes (tetrahedra) making up the original
  // simplex.

  // V0, MidPt01, MidPt02, and Vertex3=MidPt03
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt03);
  SubPts((NewG3Simps_p->Diff)[0],G3Simp.Diff[0],DMidPt03);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt01,DMidPt03);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt02,DMidPt03);

  NewG3Simps_p++;

  // V1, MidPt01, MidPt12, and Vertex3=MidPt13
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt13);
  SubPts((NewG3Simps_p->Diff)[0],G3Simp.Diff[1],DMidPt13);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt01,DMidPt13);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt12,DMidPt13);

  NewG3Simps_p++;

  // V2, MidPt23, MidPt12, and Vertex3=MidPt02
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt02);
  SubPts((NewG3Simps_p->Diff)[0],G3Simp.Diff[2],DMidPt02);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt23,DMidPt02);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt12,DMidPt02);

  NewG3Simps_p++;

  // V3, MidPt23, MidPt13, and Vertex3=MidPt03
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt03);
  NegPt((NewG3Simps_p->Diff)[0],DMidPt03);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt23,DMidPt03);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt13,DMidPt03);

  NewG3Simps_p++;

  // MidPt01, MidPt23, MidPt02, and Vertex3=MidPt12
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt12);
  SubPts((NewG3Simps_p->Diff)[0],DMidPt01,DMidPt12);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt23,DMidPt12);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt02,DMidPt12);

  NewG3Simps_p++;
  
  // MidPt01, MidPt23, MidPt12, and Vertex3=MidPt13
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt13);
  SubPts((NewG3Simps_p->Diff)[0],DMidPt01,DMidPt13);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt23,DMidPt13);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt12,DMidPt13);

  NewG3Simps_p++;
  
  // MidPt01, MidPt23, MidPt13, and Vertex3=MidPt03
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt03);
  SubPts((NewG3Simps_p->Diff)[0],DMidPt01,DMidPt03);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt23,DMidPt03);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt13,DMidPt03);

  NewG3Simps_p++;
  
  // MidPt01, MidPt23, MidPt03, and Vertex3=MidPt02
  AddPts(NewG3Simps_p->Vertex3,G3Simp.Vertex3,DMidPt02);
  SubPts((NewG3Simps_p->Diff)[0],DMidPt01,DMidPt02);
  SubPts((NewG3Simps_p->Diff)[1],DMidPt23,DMidPt02);
  SubPts((NewG3Simps_p->Diff)[2],DMidPt03,DMidPt02);
}

#endif