/**************************************************************************/
/******************  Craig solid angle calculations  **********************/
/**************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "sterdefn.h"
#include "integrat.h"
#include "stermem.h"
#include "steraid.h"
#include "stercalc.h"
#include "stertext.h"
#include "craig.h"

/**************************************************************************/
/******************  ellipse memory management  ***************************/
/**************************************************************************/

void Initialize_Ellipse(Elps *elps, Atms *at)
{
  elps->a=0.0;                     /* long axis                           */
  elps->b=0.0;                     /* short axis                          */
  elps->c=0.0;                     /* shift off origin                    */
  elps->delta=0.0;                 /* atom angle from vector G            */
  elps->A=at;                      /* pointer to atom memory              */
  elps->I[0]=Vequal(0.0,0.0,0.0);
  elps->I[1]=Vequal(0.0,0.0,0.0);
  elps->phi[0]=PI/2;
  elps->phi[1]=PI/2;               /* angles of intersepts from main axis */
  elps->mode=0;                    /* for special ellipse modes           */
}

/**************************************************************************/

void Initialize_Over(Over *over, Atms *A[MAX_OVER], int order)
{
  int i=0;
  if(order>MAX_OVER) order=MAX_OVER;
  if(over!=NULL)
  {
    for(i=0;i<MAX_OVER;i++) Initialize_Ellipse(&over->El[i],A[i]);
    over->G=Vequal(0.0,0.0,0.0);     /* vector to center of overlap       */
    over->order=order;               /* overlap order (eg. 3=triple       */
    over->nell=order;                /* number of ellipses bounding over  */
    over->solid=0.0;                 /* solid angle calculated            */
  }
/* next and prev pointers have been set so don't touch                    */
}

/**************************************************************************/

void Copy_Ellipse(Elps *new, Elps *old)
{
  new->a=old->a;
  new->b=old->b;
  new->c=old->c;
  new->delta=old->delta;
  new->I[0]=VequalV(old->I[0]);
  new->I[1]=VequalV(old->I[1]);
  new->phi[0]=old->phi[0];
  new->phi[1]=old->phi[1];
}

/**************************************************************************/

void Copy_Over(Over *new, Over *old)
{
  int i=0;

  if((new!=NULL)&&(old!=NULL))
  {
    for(i=0;i<MAX_OVER;i++) Copy_Ellipse(&new->El[i],&old->El[i]);
    new->G=VequalV(old->G);
    new->order=old->order;
    new->nell=old->nell;
    new->solid=old->solid;
  }
}

/**************************************************************************/
/**************************************************************************/

Over *New_Over(Over *old)
{
  Over *new=NULL;
  if((new=(Over *)malloc(sizeof(Over)))==NULL)
  {
	Error_Message(E_NOMEM,"New Over");
	return(NULL);
  }
  new->next=NULL;
  new->prev=NULL;
  if(old==NULL) return(new);
  new->next=old->next;
  new->prev=old;
  if(old->next) old->next->prev=new;
  old->next=new;
  return(new);
}

/**************************************************************************/
/**************************************************************************/

Over *First_Over(Over *current)
{
  if(current==NULL) return(NULL);
  while(current->prev!=NULL) current=current->prev;
  return(current);
}

/**************************************************************************/

Over *Last_Over(Over *current)
{
  if(current==NULL) return(NULL);
  while(current->next!=NULL) current=current->next;
  return(current);
}

/**************************************************************************/
/**************************************************************************/

int Get_Over_Number(Over *current)
{
  int count=0;
  if(current==NULL) return(0);
  while(count++,current->prev!=NULL) current=current->prev;
  return(count);
}

/**************************************************************************/
/**************************************************************************/

Over *Goto_Over(Over *current, int num)
{
  int count=0;
  current=First_Over(current);
  if(current==NULL) return(NULL);
  while(count++,(current->next!=NULL)&&(count!=num)) current=current->next;
  return(current);
}

/**************************************************************************/
/**************************************************************************/

Over *Next_Over(Over *current)
{
  if(current==NULL) return(NULL);
  if(current->next==NULL) return(current);
  return(current->next);
}

/**************************************************************************/

Over *Previous_Over(Over *current)
{
  if(current==NULL) return(NULL);
  if(current->prev==NULL) return(current);
  return(current->prev);
}

/**************************************************************************/
/**************************************************************************/

Over *Close_Over(Over *current)
{
  Over *old;
  if(current==NULL) return(NULL);
  if(current->prev==NULL)
  {
    if(current->next==NULL)
    {
      free(current);
      return(NULL);
    }
    current=current->next;
    free(current->prev);
    current->prev=NULL;
    return(current);
  }
  if(current->next==NULL)
  {
    current=current->prev;
    free(current->next);
    current->next=NULL;
    return(current);
  }
  old=current;
  current=old->prev;
  old->prev->next=old->next;
  old->next->prev=old->prev;
  free(old);
  return(current);
}

/**************************************************************************/

void Close_All_Overlaps(Over *over, char *title, unsigned short mode)
{
  Over *o=NULL;
  char line[LINELEN];
  int count=0;
  for(o=First_Over(over),count=0;o!=NULL;count++,o=o->next);
  if((count>0)&&(mode&VIS_BIT))
  {
    sprintf(line,"Closing temporary overlap calculation memory for %d %s overlaps",count,title);
    Out_Message(line,O_NEWLN);
  }
  for(over=First_Over(over);over!=NULL;over=Close_Over(over));
}

/**************************************************************************/
/******************  ellipse calculations  ********************************/
/**************************************************************************/

void Find_2atom_Intersepts(double alpha, double beta
                          ,Vector V[2], Vector I[2])
{
  int i=0;
  double D=0.0,T=1.0;
  double p=0.0,q=0.0,l=0.0,m=0.0;
  double a=0.0,b=0.0,c=0.0;

  while((i<3)&&(AlmostZero(T= V[0].x * V[1].y    - V[1].x * V[0].y)))
  {
    Rotate_Indices_Right(&V[0]);
    Rotate_Indices_Right(&V[1]);
    i++;
  }

  p=  ( V[0].x * cos(beta) - V[1].x * cos(alpha))/T;

  q=  ( V[1].x * V[0].z - V[0].x * V[1].z)/T;

  l=  cos(alpha) / V[0].x  -  p * V[0].y / V[0].x ;

  m=  V[0].z / V[0].x  +  q * V[0].y / V[0].x ;

  a= m*m + q*q + 1 ;
  b= 2*p*q - 2*l*m ;
  c= l*l + p*p - 1 ;

  D= b*b -4*a*c;

  if(D<0.0)
  {
    I[0]=Vequal(0.0,0.0,0.0);
    I[1]=Vequal(0.0,0.0,0.0);
  }
  else
  {
    I[0].z = (-b-sqrt(D))/(2*a);
    I[0].x = l-m*I[0].z;
    I[0].y = p+q*I[0].z;

    I[1].z = (-b+sqrt(D))/(2*a);
    I[1].x = l-m*I[1].z;
    I[1].y = p+q*I[1].z;

    I[0]=unit_vector(I[0]);
    I[1]=unit_vector(I[1]);
  }
  for(;i>0;i--)
  {
    Rotate_Indices_Left(&V[0]);
    Rotate_Indices_Left(&V[1]);
    Rotate_Indices_Left(&I[0]);
    Rotate_Indices_Left(&I[1]);
  }

}

/**************************************************************************/

void Find_2atom_G(Over *over, Vector V[2])
{
  double gap=0.0;  
  Find_2atom_Intersepts(over->El[0].A->SVangle,over->El[1].A->SVangle
                       ,V,over->El[0].I);
  over->El[1].I[0]=VequalV(over->El[0].I[0]);
  over->El[1].I[1]=VequalV(over->El[0].I[1]);
  over->G=unit_vector(Vsum(over->El[0].I[0],over->El[0].I[1]));
  over->El[0].delta=VangleV(over->G,over->El[0].A->v);
  over->El[1].delta=VangleV(over->G,over->El[1].A->v);
  gap=over->El[0].delta+over->El[1].delta
     -VangleV(over->El[0].A->v,over->El[1].A->v);
  if((gap>0.0)&&!(Small(gap)))   /* fix to solve inadequacies in ellipse  */
  {                                             /* parameter calculations */
    if(over->El[0].delta<over->El[1].delta) over->El[0].mode|=S_NEGS;
    else over->El[1].mode|=S_NEGS;
  }
}

/**************************************************************************/

void Get_Ellipse_Mode(Elps *Elly, Vector G)
{
  Vector A={0.0,0.0,0.0};                  /* vector to intersept average */
  double angle=0.0;                        /* angle between G and A       */

/* work out whether ellipse is +ve or -ve                                 */

  A=SxV(0.5,Vsum(Elly->I[0],Elly->I[1]));
  angle=VangleV(A,G);
  if(!AlmostZero(angle))
  {
    if(VangleV(Elly->A->v,A)<angle) Elly->mode|=S_NEGS;
  }
}

/**************************************************************************/

void Calculate_Ellipse_Intersept_Angles(Elps *Elly, double theta, double phi)
{
  double angle=0.0;
  Vector X[2];                             /* rotated intersept vectors   */
  Vector A={0.0,0.0,0.0};                  /* vector to atom rotated      */
  int i=0;

/* calculate intersept vectors in plane normal to vector G                */

  A=Rotate_About_Z(Elly->A->v,-1.0*phi);
  A=Rotate_About_Y(A,-1.0*theta);
  A.z=0.0;
  for(i=0;i<2;i++)
  {
    X[i]=Rotate_About_Z(Elly->I[i],-1.0*phi);
    X[i]=Rotate_About_Y(X[i],-1.0*theta);
    X[i].z=0.0;
  }

/* calculate intersept angles phi[]                                       */

  angle=VangleV(X[0],X[1]);
  Elly->phi[0]=VangleV(X[0],A);
  Elly->phi[1]=VangleV(X[1],A);
  if(AlmostZero(2*PI-angle-Elly->phi[0]-Elly->phi[1]))
  {                          /* normal situation with +ve and -ve phi's */
    Elly->phi[0]=PI-Elly->phi[0];
    Elly->phi[1]=PI-Elly->phi[1];
  }
  else
  {                                          /* both phi's on same side */
    if(Elly->phi[0]<Elly->phi[1])
    {
      Elly->phi[0]=PI-Elly->phi[0];
      Elly->phi[1]-=PI;
    }
    else
    {
      Elly->phi[1]=PI-Elly->phi[1];
      Elly->phi[0]-=PI;
    }
  }
}

/**************************************************************************/

void Calculate_Ellipse_Parameters(Elps *Elly)
{
  double s=0.0;

/* calculate ellipse coefficients using the Leach equations (White et al.)*/

  s=cos(2*Elly->A->SVangle)+cos(2*Elly->delta);
  Elly->a=sin(2*Elly->A->SVangle)/s;
  Elly->b=msqrt(2*pow(sin(Elly->A->SVangle),2)/s,"msqrt -> Calculate Ellipse Parameters (b)");
  Elly->c=sin(2*Elly->delta)/s;
  
}

/**************************************************************************/

void Setup_Two_Ellipses(Over *over, Atms *A[MAX_OVER], char mode)
{
  Vector V[2]={{0.0,0.0,0.0},{0.0,0.0,0.0}};
  double theta=0.0,phi=0.0;

  switch(mode)
  {
    case FIRST  : over->G=unit_vector(A[0]->v);
                  over->nell=1;
                  over->El[0].A=A[0];
                  break;
    case SECOND : over->G=unit_vector(A[1]->v);
                  over->nell=1;
                  over->El[0].A=A[1];
                  break;
    case BOTH   : V[0]=unit_vector(over->El[0].A->v);
                  V[1]=unit_vector(over->El[1].A->v);
                  Find_2atom_G(over,V);
                  theta=Get_Vector_Theta(over->G);
                  phi=Get_Vector_Phi(over->G);
                  over->El[0].delta=VangleV(over->El[0].A->v,over->G);
                  over->El[1].delta=VangleV(over->El[1].A->v,over->G);
                  Calculate_Ellipse_Intersept_Angles(&over->El[0],theta,phi);
                  Calculate_Ellipse_Intersept_Angles(&over->El[1],theta,phi);
                  Calculate_Ellipse_Parameters(&over->El[0]);
                  Calculate_Ellipse_Parameters(&over->El[1]);
                  break;
    default     : break;
  }
}

/**************************************************************************/
/******************  integration calculations  ****************************/
/**************************************************************************/

static double iA;
static double iB;
static double iC;
static int isign;

void Prepare_Ellipse(Elps *Elly)
{
  iA=Elly->a;
  iB=Elly->b;
  iC=Elly->c;
  if(Elly->mode&S_NEGS)
  {
    isign=-1;
    printf("n");
  }
  else isign=1;
}
/**************************************************************************/

double Ellipse(double phi)      /* this is the general ellipse equation */
{
  double p,q;
  double s;
  double eta;

  p=cos(phi)*iC/(iA*iA);
  q=pow(cos(phi)/iA,2)+pow(sin(phi)/iB,2);
  s=p*p+q*(1-(iC/iA)*(iC/iA));
  if (s<0.0) s=0.0;
  eta=(msqrt(s,"msqrt -> Ellipse (eta)")-isign*p)/q;
  return(1/msqrt(1+eta*eta,"msqrt -> Ellipse"));
}

/**************************************************************************/
/**************************************************************************/
/******************  solid angle calculations  ****************************/
/**************************************************************************/
/**************************************************************************/

void Add_Phi(Over *over, int order, double excess)
{
  int count,n;
  double phi=0.0;
  while(excess>0.0)
  {
    for(count=0,n=0;count<order;count++)
    {                                                     /* find max phi */
      if(over->El[count].phi[1]<PI/2.0)
      {
        n++;
        if(over->El[count].phi[1]>phi) phi=over->El[count].phi[1];
      }
      if(over->El[count].phi[0]<PI/2.0)
      {
        n++;
        if(over->El[count].phi[0]>phi) phi=over->El[count].phi[0];
      }
    }
    phi=PI/2.0-phi;
    if(excess/(double)n<=phi) phi=excess/(double)n;
    for(count=0;count<order;count++)
    {
      if(over->El[count].phi[0]>0.0) over->El[count].phi[0]-=phi;
      if(over->El[count].phi[1]>0.0) over->El[count].phi[1]-=phi;
    }
    excess-=phi*(double)n;
    if(AlmostZero(excess)) break;
  }
}

/**************************************************************************/

void Cut_Phi(Over *over, int order, double excess)
{
  int count,n;
  double phi=100.0;
  while(excess>0.0)
  {
    for(count=0,n=0;count<order;count++)
    {                                                     /* find min phi */
      if(over->El[count].phi[1]>0.0)
      {
        n++;
        if(over->El[count].phi[1]<phi) phi=over->El[count].phi[1];
      }
      if(over->El[count].phi[0]>0.0)
      {
        n++;
        if(over->El[count].phi[0]<phi) phi=over->El[count].phi[0];
      }
    }
    if(excess/(double)n<=phi) phi=excess/(double)n;
    for(count=0;count<order;count++)
    {
      if(over->El[count].phi[0]>0.0) over->El[count].phi[0]-=phi;
      if(over->El[count].phi[1]>0.0) over->El[count].phi[1]-=phi;
    }
    excess-=phi*(double)n;
    if(AlmostZero(excess)) break;
  }
}

/**************************************************************************/
/******************  single atom calculations  ****************************/
/**************************************************************************/

double Get_All_Single_Solids(Mol *M, unsigned mode)
{
  int count=0;
  double solid=0.0;
  Atms *A=NULL;
  char line[LINELEN];

  for(A=M->atoms,count=0;A!=NULL;A=A->next)
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Get All Single Solids");
    if(A->stat&CALC_BIT)
    {
      if(mode&VIS_BIT) Out_Message("+ ",O_BLANK);
      solid+=Single_Atom_Solid_Angle(A,mode);
      if(mode&VIS_BIT)
      {
        sprintf(line,"= %f",solid);
        Out_Message(line,O_NEWLN);
      }
      A->count=1;
    }
    else A->count=0;
  }
  return(solid);
}

/**************************************************************************/
/****************  two atom overlap calculation  **************************/
/**************************************************************************/

double Double_Overlap_Solid(Over *over, double eps, unsigned mode)
{
  double totphi=0.0;
  double integral=0.0;
  int count=0;
  char line[150];

/* check for excessively large phi range                                  */

  for(count=0;count<2;count++)
    totphi+=over->El[count].phi[0]+over->El[count].phi[1];
  totphi-=2.0*PI;
  if(!AlmostZero(totphi))
  {
    if(mode&VIS_BIT)
    {
      sprintf(line,"Double Overlap Solid [%7.4f]",totphi);
      Error_Message(E_PNOTPI,line);
    }
    if(totphi>0.0) Cut_Phi(over,2,totphi);
    if(totphi<0.0) Add_Phi(over,2,-1.0*totphi);
  }

/* calculate right ellipse                                                */

  Prepare_Ellipse(&over->El[0]);
  integral+=qsimp(Ellipse,0.0,over->El[0].phi[1],eps);

/* calculate left ellipse                                                 */

  Prepare_Ellipse(&over->El[1]);
  integral+=qsimp(Ellipse,0.0,over->El[1].phi[1],eps);

/* calculate solid angle                                                  */

  integral=2*PI-2*integral;

/* display answers if in normal total calculation                         */

  if(mode&VIS_BIT)
  {
    sprintf(line,"double(%s[%d]-%s[%d]):[%f] "
                ,over->El[0].A->name,Get_Atom_Number(over->El[0].A)
                ,over->El[1].A->name,Get_Atom_Number(over->El[1].A),integral);
    Out_Message(line,O_BLANK);
  }

  if(integral<0.0)
  {
    Error_Message(E_NEGSLD,"Double Overlap Solid");
    return(0.0);
  }
  return(integral);
}

/**************************************************************************/

double Remove_All_Double_Overlaps(Mol *M, Over **current, double eps
                               ,unsigned mode, double solid)
{
  int count,i;
  Over *new=NULL;
  Over *over=NULL;
  double chi=0.0;
  Atms *A[MAX_OVER];
  char line[LINELEN];

  for(i=0;i<MAX_OVER;i++) A[i]=NULL;
                                          /* count A[0] though whole list */
  for(A[0]=M->atoms,count=0;A[0]!=NULL;A[0]=A[0]->next)
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Remove All Double Overlaps");
    if((A[0]->SVangle!=0.0)&&(A[0]->stat&CALC_BIT))
    {                                     /* count A[1] though whole list */
      for(A[1]=A[0]->next,i=0;A[1]!=NULL;A[1]=A[1]->next)
      {
        i++; if(i>M->num_atoms) Error_Message(E_ATMNUM,"Remove All Double Overlaps");
        if ((A[0]!=A[1])&&(A[1]->SVangle!=0.0)&&(A[1]->stat&CALC_BIT))
        {
          chi=VangleV(A[0]->v,A[1]->v);
          if (chi<(A[0]->SVangle+A[1]->SVangle))
          {                           /* A[0]-A[1] steric overlap found   */
            if((new=New_Over(over))!=NULL)
            {
              over=new;          /* create the memory for counted doubles */
              Initialize_Over(over,A,2);
              if(mode&VIS_BIT) Out_Message("- ",O_BLANK);
              if (chi<=fabs(A[0]->SVangle-A[1]->SVangle))
              {                                   /* total overlap found  */
                Error_Message(E_UNSING,"Remove All Double Overlaps");
                if(mode&VIS_BIT) Out_Message("- ",O_BLANK);
                if (A[0]->SVangle>=A[1]->SVangle)
                {                                   /* A[0] overlaps A[1] */
                  Setup_Two_Ellipses(over,A,SECOND);
                  solid-=Single_Atom_Solid_Angle(A[1],mode);
                }
                else
                {                                   /* A[1] overlaps A[0] */
                  Setup_Two_Ellipses(over,A,FIRST);
                  solid-=Single_Atom_Solid_Angle(A[0],mode);
                }
                if(mode&VIS_BIT)
                {
                  sprintf(line,"= %f",solid);
                  Out_Message(line,O_NEWLN);
                }
              }
              else                               /* partial overlap found */
              {
                Setup_Two_Ellipses(over,A,BOTH);
                solid-=Double_Overlap_Solid(over,eps,mode);
                if(mode&VIS_BIT)
                {
                  sprintf(line,"= %f",solid);
                  Out_Message(line,O_NEWLN);
                }
              }
            }
          }
        }
      }
    }
  }
  *current=over;
  return(solid);
}

/**************************************************************************/
/**************************************************************************/
/****************  multiple atom overlap calculation  *********************/
/**************************************************************************/
/**************************************************************************/

/**************************************************************************/
/** test for case where one atom is enveloped by another atom            **/

int Inside_Atom(Atms *A, Atms *B)
{
  if(A==B) return(1);
  if(A->SVangle+VangleV(A->v,B->v)<=B->SVangle) return(1);
  return(0);
}

/**************************************************************************/
/** test for that atom is not enveloped by any other atom                **/

int Outside_All_Atoms(Atms *A)
{
  Atms *atoms;
  if(A==NULL) return(0);
  for(atoms=First_Atom(A);atoms!=NULL;atoms=atoms->next)
  {
    if(atoms==A) atoms=atoms->next;
    if(atoms==NULL) break;
    if(Inside_Atom(A,atoms)) return(0);
  }
  return(1);
}

/**************************************************************************/
/** test for case where overlap region in entirely within an atom        **/

int Atom_Covers_Overlap(Atms *A, Over *O)
{
  int i;
  for(i=0;i<O->nell;i++)
  {
    if(VangleV(O->El[i].I[0],A->v)>A->SVangle) return(0);
    if(VangleV(O->El[i].I[1],A->v)>A->SVangle) return(0);
  }
  return(1);
}

/**************************************************************************/
/* find overlap region that doesn't involve specified atom               **/

int Find_Over_Without_Atom(Over *po[MAX_OVER], Atms *at, int o[2]
                          ,int order, int start)
{
  int i,k,a;
  
  for(i=start;i<order;i++)  
  {
    a=0;
    for(k=0;k<order-1;k++)
    {                                                  /* find extra atom */
      if(at==po[i]->El[k].A) a++;
    }
    if(a>1) return(0);                       /* found same atom twice ??  */
    if(a==0)
    {
      if(at==NULL) return(0);
      else
      {
        o[1]=i;                                             /* found it ! */
        return(1);
      }
    }
  }
  return(0);
}
/**************************************************************************/
/* swap positions of two pointers to overs in the array po                */

void Swap_Overs(Over *po[MAX_OVER], int o[2])
{
  Over *temp;
  temp=po[o[1]];
  po[o[1]]=po[o[0]];
  po[o[0]]=temp;
}

/**************************************************************************/
/* order the n* "(n-1) overlaps" to facilitate discovery of nth overlap   */

int Order_Previous_Overlaps(Over *po[MAX_OVER], Atms *at[MAX_OVER], int order)
{
  Atms *atom=NULL;
  int i,k,o[2],a=0;
  
  if(order<3) return(0);             /* only triple or greater considered */
  for(i=0;i<order;i++)
  {                               /* n * "(n-1) overlaps" must be present */
    if(po[i]==NULL) return(0);
  }
  for(i=1;i<order;i++)
  {                                      /* map po[0] to last (n-1) atoms */
    at[i]=po[0]->El[i-1].A;
    if(at[i]==NULL) return(0);
  }                                 /* find appropriate overlap for at[1] */
  if(!Find_Over_Without_Atom(po,at[1],o,order,1)) return(0);
  o[0]=1;
  Swap_Overs(po,o);                             /* swap over 1 and o[1]   */
  for(k=0;k<order-1;k++)
  {                                                    /* find extra atom */
    a=0;
    atom=po[1]->El[k].A;
    if(atom==NULL) return(0);
    if(at[1]==atom) return(0);         /* Find_Over_Without_Atom() failed */
    for(i=2;i<order;i++)      
    {
      if(at[i]==atom) a++;
    }
    if(a>1) return(0);                       /* found same atom twice ??  */
    if(a==0)
    {
      if(at[0]==NULL) at[0]=atom;
      else return(0);                               /* two unknowns found */
    }
  }
                                        /* order the rest of the overlaps */  
  for(i=2;i<order;i++)
  {
    if(!Find_Over_Without_Atom(po,at[i],o,order,i)) return(0);
    o[0]=i;
    Swap_Overs(po,o);
  }

  for(i=0;i<order;i++)     /* final check for the correct number of atoms */
  {
    for(k=0;k<order-1;k++)
    {
      for(a=0;a<order;a++)
      {
        if(po[i]->El[k].A==at[a]) break;
      }
      if(a==order) return(0);
    }
  }
  
  return(1);
}

/**************************************************************************/
/* find the correct atomic intersepts for particular pair of atoms        */

int Find_Atomic_Intersepts(Over *dover, Atms *A[MAX_OVER]
                           ,int a, int b
                           ,Vector I[2], int inter[2])
{
  int i;
  double angle;
  I[0]=Vequal(0.0,0.0,0.0);
  I[1]=Vequal(0.0,0.0,0.0);
  inter[0]=0;
  inter[1]=0;
  for(dover=First_Over(dover);dover!=NULL;dover=dover->next)
  {                                                    /* find ab overlap */
    if(((dover->El[0].A==A[a])&&(dover->El[1].A==A[b]))
     ||((dover->El[1].A==A[a])&&(dover->El[0].A==A[b])))
    {
      inter[0]=1;      /* check that both intersepts are within all atoms */
      inter[1]=1;
      for(i=0;((i<MAX_OVER)&&(A[i]!=NULL));i++)
      {                                            /* look for exceptions */
        angle=VangleV(dover->El[0].I[0],A[i]->v)-A[i]->SVangle;
        if((angle>0.0)&&(!AlmostZero(angle))) inter[0]=0;
        angle=VangleV(dover->El[0].I[1],A[i]->v)-A[i]->SVangle;
        if((angle>0.0)&&(!AlmostZero(angle))) inter[1]=0;
      }
      if(inter[0]) I[0]=VequalV(dover->El[0].I[0]);
      if(inter[1]) I[1]=VequalV(dover->El[0].I[1]);
      return(1);
    }
  }
  return(0);
}

/**************************************************************************/
/* determine the nature and parameters of the new nth order overlap       */

Over *New_Multiple_Overlap(Over *over, Over *dover, int order
                          ,Over *O[MAX_OVER], Atms *A[MAX_OVER])
{
  int i=0,l=0,k=0,K=0,n=0,p=0;
  Vector I[MAX_OVER-1][2];                  /* vectors for phi[] calc.    */
  int inter[MAX_OVER-1][2];
  char mode=S_NORMAL;             /* normal overlap combination expected  */
  int multi_atom=0;
  double theta=0.0,phi=0.0;                       /* G's angular position */

  for(i=0;i<order;i++)
  {
    if(O[i]==NULL) return(NULL);
    if(A[i]==NULL) return(NULL);
  }
  if((over=New_Over(over))==NULL) return(NULL);
  Initialize_Over(over,A,order);

  for(i=0;i<order-1;i++)
  {
    I[i][0]=Vequal(0.0,0.0,0.0);
    I[i][1]=Vequal(0.0,0.0,0.0);
  }
  
  for(i=0;i<order;i++)
  {
    if(Atom_Covers_Overlap(A[i],O[i]))
    {                               /* region already done for (n-1) case */
      Copy_Over(over,O[i]);
      over->order=order;
      return(over);
    }
  }

  for(i=0;i<order;i++)                 /* loop through all relevant atoms */
  {
    over->El[i].A=A[i];                               /* set ellipse atom */
    for(l=0;l<order-1;l++) 
    {
      inter[l][0]=0;
      inter[l][1]=0;
    }
    for(k=0,n=0,p=0;k<order;k++,n++)      /* loop through all other atoms */
    {
      if(k==i) k++;
      if(k>=order) break;
      Find_Atomic_Intersepts(dover,A,i,k,I[n],inter[n]);
      if(inter[n][0]) p++;
      if(inter[n][1]) p++;
    }
    if(p<2) return(Close_Over(over));        /* no multiple overlap found */
    if((p!=2)&&(p!=2*n)) return(Close_Over(over));       /* unknown combo */
    else if(p==2)
    {                                         /* normal type of intersept */
      for(K=0,n=0;K<order-1;K++)
      {
        if(inter[K][0])
        {
          over->El[i].I[n]=VequalV(I[K][0]);
          n++;
        }
        if(n>2) return(Close_Over(over));            /* unexpected error */
        if(inter[K][1])
        {
          over->El[i].I[n]=VequalV(I[K][1]);
          n++;
        }
        if(n>2) return(Close_Over(over));             /* unexpected error */
      }
      if(n!=2) return(Close_Over(over));              /* unexpected error */
    }
    else if(K==2*n)
    {                                        /* special type of intersept */
      mode=S_GCENT;
      multi_atom=i;                /* atom i used for partial single atom */
      over->El[i].mode|=S_GCENT;
    }
  }
  
/* calculate the vector G within the n-overlap region                     */

  if(mode&S_GCENT)
  {
    over->G=VequalV(unit_vector(over->El[multi_atom].A->v));
  }
  else
  {
    over->G=Vequal(0.0,0.0,0.0);
    for(i=0;i<order;i++)
    {
      over->G=Vsum(over->G,over->El[i].I[0]);
      over->G=Vsum(over->G,over->El[i].I[1]);
    }
    over->G=unit_vector(over->G);
  }
  
/* check that G is not too near any atom centres                          */

  for(i=0;i<order;i++)
  {
    over->El[i].delta=VangleV(over->G,over->El[i].A->v);
    if((!(mode&S_GCENT))
       &&(Similar(over->El[i].A->SVangle
                 ,over->El[i].A->SVangle-over->El[i].delta
                 ,S_PERC)))
    {
      over->El[i].mode|=S_GCENT;
      over->G=unit_vector(A[i]->v);
      mode|=S_GCENT;
    }
  }

/* make absolutely sure that this vector is indeed within all atoms       */

  for(i=0;i<order;i++)
  {
    over->El[i].delta=VangleV(over->G,over->El[i].A->v);
    if(over->El[i].delta>A[i]->SVangle)
    {
      Error_Message(E_BAD_G,"New Multiple Overlap");
      return(Close_Over(over));
    }
  }
  
/* calculate the ellipse parameters for all necessary ellipses            */

  for(i=0;i<order;i++)
  {
    Get_Ellipse_Mode(&over->El[i],over->G);
    if(!(over->El[i].mode&S_GCENT))
    {
      theta=Get_Vector_Theta(over->G);
      phi=Get_Vector_Phi(over->G);
      Calculate_Ellipse_Intersept_Angles(&over->El[i],theta,phi);
      Calculate_Ellipse_Parameters(&over->El[i]);
    }
  }
  
  if(mode&S_GCENT)
  {                            /* centred ellipses dealt with differently */
    for(i=0;i<order;i++)
    {
      if(over->El[i].mode&S_GCENT)
      {
        over->El[i].phi[0]=PI;
        over->El[i].phi[1]=PI;
        for(n=0;n<order;n++)
        {
          if(n!=i)
          {
            over->El[i].phi[0]-=over->El[n].phi[0];
            over->El[i].phi[1]-=over->El[n].phi[1];
          }
        }
      }
    }
  }
  return(over);
}

/**************************************************************************/
/* perform the necessary integration on the new nth order overlap         */

double Multiple_Overlap_Solid(Over *over, int order
                             ,double eps, unsigned mode)
{
  double totphi=0.0;
  double integral=0.0;
  double multi_solid=0.0;           /* extra solid angle for S_GCENT mode */
  int count;
  char line[150];

/* check for excessively large phi range                                  */

  for(count=0;count<order;count++)
    totphi+=over->El[count].phi[0]+over->El[count].phi[1];
  totphi-=2.0*PI;
  if(!AlmostZero(totphi))
  {
    if(mode&VIS_BIT)
    {
      sprintf(line,"Multiple Overlap Solid [%7.4f]",totphi);
      Error_Message(E_PNOTPI,line);
    }
    if(totphi>0.0) Cut_Phi(over,order,totphi);
    if(totphi<0.0) Add_Phi(over,order,-1.0*totphi);
    for(count=0;count<order;count++)
      totphi+=over->El[count].phi[0]+over->El[count].phi[1];
    totphi-=2.0*PI;
  }

/* integrate over each of the ellipses present                            */

  if(AlmostZero(totphi))
  {
    for(count=0;count<order;count++)
    {
      if(over->El[count].mode&S_GCENT)
      {
        printf("m");
	totphi=over->El[count].phi[0]+over->El[count].phi[1];
        multi_solid=totphi/(2*PI)
                   *Single_Atom_Solid_Angle(over->El[count].A,0);
      }
      else
      {
        Prepare_Ellipse(&over->El[count]);
        if(over->El[count].phi[0]<0.0)
        {
          integral+=qsimp(Ellipse,-over->El[count].phi[0]
                                  ,over->El[count].phi[1],eps);
        }
        else if(over->El[count].phi[1]<0.0)
        {
          integral+=qsimp(Ellipse,-over->El[count].phi[1]
                                  ,over->El[count].phi[0],eps);
        }
        else
        {
          integral+=qsimp(Ellipse,0.0,over->El[count].phi[0],eps);
          integral+=qsimp(Ellipse,0.0,over->El[count].phi[1],eps);
        }
      }
    }
  
/* calculate solid angle from total integration results                   */

    integral=2*PI-totphi-integral+multi_solid;
  }

/* display answers if in normal total calculation                         */

  if(mode&VIS_BIT)
  {
    sprintf(line,"[%f] ",integral);
    Out_Message(line,O_BLANK);
  }

  if(integral<0.0)
  {
    Error_Message(E_NEGSLD,"Multiple Overlap Solid");
    return(0.0);
  }
  return(integral);
}

/**************************************************************************/
/* find nth order overlap and calculate the relevant solid angle          */

double Calc_Multiple_Overlap(Over *po[MAX_OVER], Over **over
                         ,Over *dover, int order
                         ,double eps, unsigned mode
                         ,char *multiple, double solid)
{
  Over *new=NULL;
  char line[LINELEN];
  char atmstr[25];
  Atms *A[MAX_OVER];
  Over *O[MAX_OVER];
  int i;
  short sign=1;
  
  /* initialize */
  
  if(order/2*2==order) sign=1; else sign=-1;
  for(i=0;i<MAX_OVER;i++)
  {
    A[i]=NULL;
    O[i]=po[i];
  }
  
  /* test for correct number of atoms, and correct possible overlaps      */
  
  if(!Order_Previous_Overlaps(O,A,order+1)) return(solid); 
                                            /* not a possible combination */
  
  if((new=New_Multiple_Overlap(*over,dover,order+1,O,A))!=NULL)
  {
    if((new!=NULL)&&(new!=*over))              /* multiple overlap found  */
    {
      *over=new;
      if(mode&VIS_BIT)
      {
        if(sign>0) sprintf(line,"+ %s[%d](",multiple,new->nell);
        else sprintf(line,"- %s[%d](",multiple,new->nell);
        for(i=0;(i<MAX_OVER)&&(new->El[i].A!=NULL);i++)
        {
          sprintf(atmstr,"%s[%d]-",A[i]->name,Get_Atom_Number(A[i]));
          strcat(line,atmstr);
        }
        line[strlen(line)-1]=0;
        strcat(line,") ");
        Out_Message(line,O_BLANK);
      }
      if(new->solid==0.0) 
        solid+=sign*Multiple_Overlap_Solid(new,new->nell,eps,mode);
      else solid+=sign*new->solid;
      if(mode&VIS_BIT)
      {
        sprintf(line,"= %f",solid);
        Out_Message(line,O_NEWLN);
      }
    }
  }
  return(solid);
}

/**************************************************************************/
/******* recursive routine to search possible (n-1) combinations **********/
/**************************************************************************/

double Multiple_Overlap(Over *pover, Over **over, Over *dover
                       ,int order, int level, double eps
                       ,unsigned mode, char *multiple, double solid)
{
  static Over *po[MAX_OVER];
  int i;

  if(pover==NULL) return(solid);               /* miss call               */
  if(level==0)                                 /* initialize              */
  {
    while(pover->prev!=NULL) pover=pover->prev;
    for(i=0;i<MAX_OVER;i++) po[i]=NULL;
  }
  if(level>order)
  {
    Error_Message(E_LEVHI,"Multiple Overlap");
    return(solid);                             /* miss call               */
  }
  if(pover->order!=order)
  {
    Error_Message(E_BADORD,"Multiple Overlap");
    return(solid);                             /* miss call               */
  }
  
  for(po[level]=pover;po[level]!=NULL;pover=pover->next,po[level]=pover)
  {
    if(level<order)                            /* recurse further         */
    {
      solid=Multiple_Overlap(pover->next,over,dover,order,1+level,eps,mode,multiple,solid);
    }
    else                                       /* final depth, calculate  */
    {
      solid=Calc_Multiple_Overlap(po,over,dover,order,eps,mode,multiple,solid);
    }
  }
  return(solid);
}

/**************************************************************************/
/**************************************************************************/
/******  Main Counting algorithm to look for all steric overlaps  *********/
/**************************************************************************/
/**************************************************************************/

double New_Craig_Counting(Mol *M, double eps, unsigned mode)
{
  Over *over[MAX_OVER];
  char multiple[10][MAX_OVER];
  double solid=0.0;
  char line[LINELEN];
  Atms *at=NULL;
  int i=0;
  
  if(M->multi>MAX_OVER) M->multi=MAX_OVER;
  for(i=0;i<M->multi;i++)
  {
    over[i]=NULL;
    switch(i+1)
    {
      case 2  : strcpy(multiple[i],"double"); break;
      case 3  : strcpy(multiple[i],"triple"); break;
      case 4  : strcpy(multiple[i],"quadruple"); break;
      case 5  : strcpy(multiple[i],"quintuple"); break;
      case 6  : strcpy(multiple[i],"sextuple"); break;
      case 7  : strcpy(multiple[i],"septuple"); break;
      case 8  : strcpy(multiple[i],"octuple"); break;
      case 9  : strcpy(multiple[i],"nonuple"); break;
      case 10 : strcpy(multiple[i],"dectuple"); break;
      default : sprintf(multiple[i],"order=%d",i+1); break;
    }
  }
  M->atoms=First_Atom(M->atoms);        /* start main loops at first atom */

/****************** ignore all fully overlapped atoms *********************/
   
  for(at=M->atoms;at!=NULL;at=at->next)
  {
    if((at->stat&MAIN_BIT)&&(Outside_All_Atoms(at)))
    at->stat=at->stat|CALC_BIT; else at->stat=at->stat&(FULL_BIT^CALC_BIT);
  }

/****************** add all single atom solid angles **********************/

  if(mode&VIS_BIT) Out_Message("Calculating all single atom solid angles",O_NEWLN);
  solid=Get_All_Single_Solids(M,mode);

/****************** remove all double atom overlaps ***********************/

  if(mode&VIS_BIT) Out_Message("Removing all double overlaps",O_NEWLN);
  solid=Remove_All_Double_Overlaps(M,&over[1],eps,mode,solid);

/****************** add/subtract all higher order overlaps ****************/

  for(i=2;i<M->multi;i++)
  {
    if(mode&VIS_BIT)
    { 
      if(i/2*2==i) sprintf(line,"Adding all %s overlaps",multiple[i]);
      else sprintf(line,"Removing all %s overlaps",multiple[i]);
      Out_Message(line,O_NEWLN);
    }
    solid=Multiple_Overlap(over[i-1],&over[i],over[1]
                          ,i,0,eps,mode,multiple[i],solid);
    if(over[i]==NULL) break;
  }

/**************************** finish off **********************************/

  for(i=1;i<M->multi;i++) Close_All_Overlaps(over[i],multiple[i],mode);
   
  if (solid>4*PI) { Error_Message(E_STOBIG,"New Craig Counting"); solid=4*PI; }
  if (solid<0.0)  { Error_Message(E_STOSML,"New Craig Counting"); solid=0.0; }
  return(solid);
}

/**************************************************************************/
/******************  The End ...  *****************************************/
/**************************************************************************/