/**************************************************************************/
/********************  Projected Area calculations  ***********************/
/**************************************************************************/

#include 
#include 
#include 
#include 
#include "sterdefn.h"
#include "integrat.h"
#include "stermem.h"
#include "steraid.h"
#include "stercalc.h"
#include "stertext.h"
#include "craig.h"
#include "proja.h"

/**************************************************************************/
/******************  circle memory management  ***************************/
/**************************************************************************/

Pover *New_Pover(Pover *old)
{
  Pover *new=NULL;
  if((new=(Pover *)malloc(sizeof(Pover)))==NULL)
  {
    Error_Message(E_NOMEM,"New Pover");
    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);
}

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

Circ *New_Circ(Circ *old)
{
  Circ *new=NULL;
  if((new=(Circ *)malloc(sizeof(Circ)))==NULL)
  {
    Error_Message(E_NOMEM,"New Circ");
    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);
}

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

Mint *New_Mint(Mint *old)
{
  Mint *new=NULL;
  if((new=(Mint *)malloc(sizeof(Mint)))==NULL)
  {
    Error_Message(E_NOMEM,"New Mint");
    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);
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Pover *Goto_Pover(Pover *current, int num)
{
  int count=0;
  current=First_Pover(current);
  if(current==NULL) return(NULL);
  while(count++,(current->next!=NULL)&&(count!=num)) current=current->next;
  return(current);
}

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

Circ *Goto_Circ(Circ *current, int num)
{
  int count=0;
  current=First_Circ(current);
  if(current==NULL) return(NULL);
  while(count++,(current->next!=NULL)&&(count!=num)) current=current->next;
  return(current);
}

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

Mint *Goto_Mint(Mint *current, int num)
{
  int count=0;
  current=First_Mint(current);
  if(current==NULL) return(NULL);
  while(count++,(current->next!=NULL)&&(count!=num)) current=current->next;
  return(current);
}

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

Pover *Next_Pover(Pover *current)
{
  if(current==NULL) return(NULL);
  if(current->next==NULL) return(current);
  return(current->next);
}

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

Pover *Previous_Pover(Pover *current)
{
  if(current==NULL) return(NULL);
  if(current->prev==NULL) return(current);
  return(current->prev);
}

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

Circ *Next_Circ(Circ *current)
{
  if(current==NULL) return(NULL);
  if(current->next==NULL) return(current);
  return(current->next);
}

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

Circ *Previous_Circ(Circ *current)
{
  if(current==NULL) return(NULL);
  if(current->prev==NULL) return(current);
  return(current->prev);
}

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

Mint *Next_Mint(Mint *current)
{
  if(current==NULL) return(NULL);
  if(current->next==NULL) return(current);
  return(current->next);
}

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

Mint *Previous_Mint(Mint *current)
{
  if(current==NULL) return(NULL);
  if(current->prev==NULL) return(current);
  return(current->prev);
}

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

Pover *Close_Pover(Pover *current)
{
  Pover *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);
}

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

Circ *Close_Circ(Circ *current)
{
  Circ *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);
}

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

Mint *Close_Mint(Mint *current)
{
  Mint *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);
}

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

Circ *Close_All_Circles(Circ *circ)
{
  for(circ=First_Circ(circ);circ!=NULL;circ=Close_Circ(circ));
  return(circ);
}

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

Mint *Close_All_Mints(Mint *mint)
{
  for(mint=First_Mint(mint);mint!=NULL;mint=Close_Mint(mint));
  return(mint);
}

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

Pover *Close_Current_Pover(Pover *over)
{
  over->Cr=Close_All_Circles(over->Cr);
  over->Cr=Close_All_Circles(over->Cr);
  return(Close_Pover(over));
}

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

void Close_All_Poverlaps(Pover *over, char *title, unsigned short mode)
{
  Pover *o=NULL;
  char line[LINELEN];
  int count=0;
  for(o=First_Pover(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_Pover(over);over!=NULL;over=Close_Current_Pover(over));
}

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

void Initialize_Mint(Mint *mint, Atms *A, Atms *B)
{
  mint->I=Vequal(0.0,0.0,0.0);
  mint->A[0]=A;
  mint->A[1]=B;
#ifdef DEBUG
  mint->num=0;
#endif
}

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

void Initialize_Circle(Circ *circ, Atms *at)
{
  circ->radius=0.0;                /* circle radius                       */
  if(at!=NULL) circ->radius=at->tradius;
  circ->delta=0.0;                 /* atom distance from position G       */
  circ->A=at;                      /* pointer to atom memory              */
  circ->I[0]=Vequal(0.0,0.0,0.0);
  circ->I[1]=Vequal(0.0,0.0,0.0);
  circ->theta=0.0;                 /* angle between intersepts            */
  circ->phi=0.0;                   /* angle between i[0] and g            */
  circ->gamma=0.0;                 /* angle between gi[0] and gi[1]       */
  circ->mode=0;                    /* for special circle modes            */
#ifdef DEBUG
  circ->num=0;
#endif
}

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

int Add_Circle(Pover *over, Atms *A)
{
  Circ *cr=NULL;
  if(over==NULL) return(0);
  if((cr=New_Circ(over->Cr))==NULL) return(0);
  over->Cr=cr;
  Initialize_Circle(over->Cr,A);
  return(1);
}

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

int Make_Circles(Pover *over, Atms *A[], int order)
{
  int i=0;
  if(order>MAX_OVER) order=MAX_OVER;
  if(over==NULL) return(0);
  over->Cr=Close_All_Circles(over->Cr);
  for(i=0;(incir=order;
  return(1);
}

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

int Initialize_Pover(Pover *over, int order)
{
  if(order>MAX_OVER) order=MAX_OVER;
  if(over==NULL) return(0);
  over->Cr=NULL;
  over->Mi=NULL;
  over->G=Vequal(0.0,0.0,0.0);       /* vector to center of overlap       */
  over->order=order;                 /* overlap order (eg. 3=triple       */
  over->ncir=order;                  /* number of circles bounding over   */
  over->area=0.0;                    /* area calculated                   */
#ifdef DEBUG
  over->num=0;
#endif
  return(1);
/* next and prev pointers have been set so don't touch                    */
}

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

void Copy_Mint(Mint *new, Mint *old)
{
  new->I=VequalV(old->I);
  new->A[0]=old->A[0];
  new->A[1]=old->A[1];
}

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

void Copy_Circle(Circ *new, Circ *old)
{
  new->radius=old->radius;
  new->delta=old->delta;
  new->A=old->A;
  new->I[0]=VequalV(old->I[0]);
  new->I[1]=VequalV(old->I[1]);
  new->theta=old->theta;
  new->phi=old->phi;
  new->gamma=old->gamma;
  new->mode=old->mode;
}

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

void Copy_Pover(Pover *new, Pover *old)
{
  Circ *cr=NULL;
  Circ *ncr=NULL;
  Mint *mi=NULL;
  Mint *nmi=NULL;

  if((new!=NULL)&&(old!=NULL))
  {
    new->Cr=Close_All_Circles(new->Cr);
    for(cr=First_Circ(old->Cr);cr!=NULL;cr=cr->next)
    {
      if((ncr=New_Circ(new->Cr))==NULL) break;
      new->Cr=ncr;
      Copy_Circle(new->Cr,cr);
    }
    new->Mi=Close_All_Mints(new->Mi);
    for(mi=First_Mint(old->Mi);mi!=NULL;mi=mi->next)
    {
      if((nmi=New_Mint(new->Mi))==NULL) break;
      new->Mi=nmi;
      Copy_Mint(new->Mi,mi);
    }
    new->G=VequalV(old->G);
    new->order=old->order;
    new->ncir=old->ncir;
    new->area=old->area;
  }
}

/**************************************************************************/
/*******************  circle calculations  ********************************/
/**************************************************************************/

double deltaAB(Atms *A, Atms *B)
{
  double d=0.0, dx, dy;
  dx=A->tv.x-B->tv.x;
  dy=A->tv.y-B->tv.y;
  d = msqrt(dx*dx+dy*dy,"deltaAB");
  return(d);
}

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

int Find_2atom_Intersepts_P(double rA, double rB, Vector V[2], Vector I[2])
{
  double dAB, dx, dy;
  double cosP, sinP, cosT, sinT;
  int q=0;
  double a, l, p, t, s, m;
  static double Qx[4]={ 1,-1,-1, 1};
  static double Qy[4]={ 1, 1,-1,-1};

/* test for unusual cases */

  if((AlmostZero(rA))||(AlmostZero(rB)))
  {
    if(rA>=rB) return(1);
    else return(2);
  }
  dx=fabs(V[0].x-V[1].x);
  dy=fabs(V[0].y-V[1].y);
  dAB = msqrt(dx*dx+dy*dy,"dAB in Find 2atom Intersepts P");
  if(dAB>=(rA+rB))
  {
    Error_Message(E_NOOVER,"Find 2atom Intersepts P");
    return(-1);
  }
  if(AlmostZero(dAB))
  {
    if(rA>=rB) return(1);
    else return(2);
  }

/* find quadrant for B coords centred on xA, yA */

  if(V[1].y>=V[0].y) { if(V[1].x>=V[0].x) q=0; else q=1; }
  else               { if(V[1].x>=V[0].x) q=3; else q=2; }

/* now start actual intersept calculation */

  a = (rA*rA - rB*rB + dAB*dAB) / (2*dAB);
  l = msqrt(rA*rA - a*a,"Find 2atom Intersepts P");
                           /* above tests should ensure no error */

  cosP = a/rA;
  sinP = l/rA;
  cosT = dx/dAB;
  sinT = dy/dAB;
  
  p = cosT * cosP;
  t = sinT * sinP;
  s = sinT * cosP;
  m = cosT * sinP;
  
  I[0].z = 0.0;
  I[0].x = V[0].x + Qx[q] * rA * (p-t);
  I[0].y = V[0].y + Qy[q] * rA * (s+m);

  I[1].z = 0.0;
  I[1].x = V[0].x + Qx[q] * rA * (p+t);
  I[1].y = V[0].y + Qy[q] * rA * (s-m);

  return(0);
}

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

int Calculate_Circle_Intersept_Angles(Circ *circ, Vector *G, Vector *A)
{
  Vector i[2], g;
  double dAG, g2;


/**/
/* first to calculate gamma for checking purposes */
/**/
/* calculate i[] at G origin */

  i[0].x=circ->I[0].x-G->x;
  i[0].y=circ->I[0].y-G->y;
  i[0].z=0.0;
  i[1].x=circ->I[1].x-G->x;
  i[1].y=circ->I[1].y-G->y;
  i[1].z=0.0;

/* calculate angle between two intersept vectors */

  circ->gamma=VangleV(i[0],i[1]);

/**/
/* now to calculate the usable values, theta and phi */
/**/
/* calculate g at A origin */

  g.x=G->x-A->x;
  g.y=G->y-A->y;
  g.z=0.0;

  dAG = msqrt(g.x*g.x+g.y*g.y,"dAG in Calculate Circle Intersept Angles");
  if(dAG>=circ->radius)
  {
    Error_Message(E_BAD_G,"Calculate Circle Intersept Angles");
    return(0);
  }

/* calculate i[] at A origin */

  i[0].x=circ->I[0].x-A->x;
  i[0].y=circ->I[0].y-A->y;
  i[0].z=0.0;
  i[1].x=circ->I[1].x-A->x;
  i[1].y=circ->I[1].y-A->y;
  i[1].z=0.0;

/* calculate angle between two intersept vectors */

  circ->theta=VangleV(i[0],i[1]);

  if(circ->mode&A_NEGS) circ->theta=2*PI-circ->theta;

  if(AlmostZero(dAG))
  {
    circ->phi=0.0;
    circ->mode|=A_GCENT;     /* G centred on A */
    return(1);
  }

/* if G not at A calculate G */

  g2 = dot_product(g,i[1]);

  circ->phi=VangleV(i[0],g);
  
  if(((circ->mode&A_NEGS)&&(g2>=0.0))
   ||((!(circ->mode&A_NEGS))&&(g2<0))) circ->phi=2*PI-circ->phi;

/**/
/* first to calculate gamma for checking purposes */
/**/
/* calculate i[] at G origin */

  i[0].x=circ->I[0].x-G->x;
  i[0].y=circ->I[0].y-G->y;
  i[0].z=0.0;
  i[1].x=circ->I[1].x-G->x;
  i[1].y=circ->I[1].y-G->y;
  i[1].z=0.0;

/* calculate angle between two intersept vectors */

  circ->gamma=VangleV(i[0],i[1]);

  return(1);
}

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

int Find_2atom_G_P(Pover *over)
{
  Vector V[2]={{0.0,0.0,0.0},{0.0,0.0,0.0}};
  double gap=0.0;
  int n=0;
  Circ *Cr[2]={NULL,NULL};
  Mint *mint=NULL;
  
  if(((Cr[0]=Goto_Circ(over->Cr,1))==NULL)
   ||((Cr[1]=Goto_Circ(over->Cr,2))==NULL)
   ||(Cr[0]->A==NULL)||(Cr[1]->A==NULL))
    return(Error_Message(E_BDCIRC,"Find 2atom G P"));
  V[0]=VequalV(Cr[0]->A->tv);
  V[1]=VequalV(Cr[1]->A->tv);
  if((n=Find_2atom_Intersepts_P(Cr[0]->radius,Cr[1]->radius,V,Cr[0]->I))!=0) return(n);
  Cr[1]->I[0]=VequalV(Cr[0]->I[0]);
  Cr[1]->I[1]=VequalV(Cr[0]->I[1]);

/* add extra intersept memory here */

  over->Mi=Close_All_Mints(over->Mi);
  if((mint=New_Mint(over->Mi))==NULL) return(0);
  Initialize_Mint(mint,Cr[0]->A,Cr[1]->A);
  mint->I=VequalV(Cr[0]->I[0]);
  over->Mi=mint;
  if((mint=New_Mint(over->Mi))==NULL) return(0);
  Initialize_Mint(mint,Cr[0]->A,Cr[1]->A);
  mint->I=VequalV(Cr[0]->I[1]);
  over->Mi=mint;

/* end of new addition 6/12/95 */

  over->G=SxV(0.5L,Vsum(Cr[0]->I[0],Cr[0]->I[1]));
  Cr[0]->delta=seperate(over->G,V[0]);
  Cr[1]->delta=seperate(over->G,V[1]);
  Cr[0]->mode=A_NORMAL;
  Cr[1]->mode=A_NORMAL;
  gap=Cr[0]->delta+Cr[1]->delta-seperate(V[0],V[1]);
  if((gap>0.0)&&!(Small(gap)))  /* case where theta is greater than PI */
  {
    if(Cr[0]->deltadelta) Cr[0]->mode|=A_NEGS;
    else Cr[1]->mode|=A_NEGS;
  }
  Calculate_Circle_Intersept_Angles(Cr[0],&over->G,&V[0]);
  Calculate_Circle_Intersept_Angles(Cr[1],&over->G,&V[1]);
  return(0);
}

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

void Get_Circle_Mode(Circ *circ, Vector G)
{
  Vector V={0.0,0.0,0.0};                  /* vector to intersept average */
  double dAG=0.0;                          /* distance between A and G    */
  double dAV=0.0;                          /* distance between A and V    */
  double aVG=0.0;                          /* angle between AG and AV     */

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

  V=SxV(0.5,Vsum(circ->I[0],circ->I[1]));
  dAG=seperate(circ->A->tv,G);
  dAV=seperate(circ->A->tv,V);
  if(dAG>dAV)
  {
    aVG=VangleV(Vdif(V,circ->A->tv),Vdif(G,circ->A->tv));
                                               /* angle between AV and AG */
    aVG=cos(aVG);
    if((!AlmostZero(aVG))&&(aVG>0.0)&&(dAG>dAV/aVG)) circ->mode|=A_NEGS;
  }
}

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

void Setup_Two_Circles(Pover *over, Atms *A[MAX_OVER], char mode)
{
  Atms *at[2]={NULL,NULL};

  while(mode)
  {
    switch(mode)
    {
      case FIRST :
        at[0]=A[0];
        Make_Circles(over,at,1);
        over->G=VequalV(A[0]->tv);
        mode=0;
        break;
      case SECOND : 
        at[0]=A[1];
        Make_Circles(over,at,1);
        over->G=VequalV(A[1]->tv);
        mode=0;
        break;
      case BOTH : 
        Make_Circles(over,A,2);
        switch(Find_2atom_G_P(over))
        {
          case 0 : mode=0; break;        /* normal */
          case 1 : mode=FIRST;  break;   /* FIRST  */
          case 2 : mode=SECOND; break;   /* SECOND */
          default: mode=0; break;
        }
        break;
      default : mode=0; break;
    }
  }
}

/**************************************************************************/
/**************************************************************************/
/******************  area calculations  ****************************/
/**************************************************************************/
/**************************************************************************/

double Single_Atom_Area(Atms *A, unsigned mode)
{
  double area;
  char line[LINELEN];
  area=PI*A->tradius*A->tradius;
  if(mode&VIS_BIT)
  {
    sprintf(line,"single(%s[%d]): [%f]",A->name,Get_Atom_Number(A),area);
    Out_Message(line,O_BLANK);
  }
  return(area);
}

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

double Get_All_Single_Areas(Mol *M, unsigned mode)
{
  int count=0;
  double area=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 Areas");
    if(A->stat&CALC_BIT)
    {
      if(mode&VIS_BIT) Out_Message("+ ",O_BLANK);
      area+=Single_Atom_Area(A,mode);
      if(mode&VIS_BIT)
      {
        sprintf(line,"= %f",area);
        Out_Message(line,O_NEWLN);
      }
      A->count=1;
    }
    else A->count=0;
  }
  return(area);
}

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

double Double_Poverlap_Area(Pover *over, unsigned mode)
{
  double totgamma=0.0;
  double area=0.0;
  char line[150];
  Circ *cr=NULL;

/* check for excessively large theta range                                  */

  for(cr=First_Circ(over->Cr);cr!=NULL;cr=cr->next) totgamma+=cr->gamma;
  totgamma-=2.0*PI;
  if(!AlmostZero(totgamma))
  {
    if(mode&VIS_BIT)
    {
      sprintf(line,"Double Poverlap Area [%7.4f]",totgamma);
      Error_Message(E_GNOTPI,line);
    }
  }

/* calculate areas                                                        */

  for(cr=First_Circ(over->Cr);cr!=NULL;cr=cr->next)
  {
    if(cr->mode&A_IGNORE) continue;
    area+= 0.5 * cr->radius * (cr->theta * cr->radius
               - cr->delta * (sin(cr->phi) + sin(cr->theta-cr->phi)));
  }

/* display answers if in normal total calculation                         */

  if(mode&VIS_BIT)
  {
    if((cr=Goto_Circ(over->Cr,1))==NULL) Error_Message(E_BDCIRC,"Double Poverlap Area");
    else
    {
      sprintf(line,"double(%s[%d]",cr->A->name,Get_Atom_Number(cr->A));
      Out_Message(line,O_BLANK);
    }
    if((cr=Goto_Circ(over->Cr,2))==NULL) Error_Message(E_BDCIRC,"Double Poverlap Area");
    else
    {
      sprintf(line,"-%s[%d]):[%f] ",cr->A->name,Get_Atom_Number(cr->A),area);
      Out_Message(line,O_BLANK);
    }
  }

  if(area<0.0)
  {
    Error_Message(E_NEGARE,"Double Poverlap Area");
    return(0.0);
  }
  return(area);
}

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

double Remove_All_Double_Poverlaps(Mol *M, Pover **current
                               ,unsigned mode, double area)
{
  int count,i;
  Pover *new=NULL;
  Pover *over=NULL;
  double delta=0.0;
  Atms *A[MAX_OVER];
  char line[LINELEN];

  for(i=0;iatoms,count=0;A[0]!=NULL;A[0]=A[0]->next)
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Remove All Double Poverlaps");
    if((!AlmostZero(A[0]->tradius))&&(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 Poverlaps");
        if ((A[0]!=A[1])&&(!AlmostZero(A[1]->tradius))&&(A[1]->stat&CALC_BIT))
        {
          delta=deltaAB(A[0],A[1]);
          if (delta<(A[0]->tradius+A[1]->tradius))
          {                           /* A[0]-A[1] steric overlap found   */
            if((new=New_Pover(over))!=NULL)
            {
              over=new;          /* create the memory for counted doubles */
              Initialize_Pover(over,2);
              if(mode&VIS_BIT) Out_Message("- ",O_BLANK);
              if (delta<=fabs(A[0]->tradius-A[1]->tradius))
              {                                   /* total overlap found  */
                Error_Message(E_UNSING,"Remove All Double Poverlaps");
                if(mode&VIS_BIT) Out_Message("- ",O_BLANK);
                if (A[0]->tradius>=A[1]->tradius)
                {                                   /* A[0] overlaps A[1] */
                  Setup_Two_Circles(over,A,SECOND);
                  area-=Single_Atom_Area(A[1],mode);
                }
                else
                {                                   /* A[1] overlaps A[0] */
                  Setup_Two_Circles(over,A,FIRST);
                  area-=Single_Atom_Area(A[0],mode);
                }
                if(mode&VIS_BIT)
                {
                  sprintf(line,"= %f",area);
                  Out_Message(line,O_NEWLN);
                }
              }
              else                               /* partial overlap found */
              {
                Setup_Two_Circles(over,A,BOTH);
                area-=Double_Poverlap_Area(over,mode);
                if(mode&VIS_BIT)
                {
                  sprintf(line,"= %f",area);
                  Out_Message(line,O_NEWLN);
                }
              }
            }
          }
        }
      }
    }
  }
  *current=over;
  return(area);
}

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

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

int Inside_Atom_P(Atms *A, Atms *B)
{
  double dAB, dR;
  if(A==B) return(1);
  dAB=deltaAB(A,B);
  dR=B->tradius-A->tradius;
  if(dAB<=dR) return(1);
  if((AlmostZero(dAB))&&(AlmostZero(dR))) return(1);
  return(0);
}

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

int Outside_All_Atoms_P(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(!atoms->stat&CALC_BIT) continue;
    if(Inside_Atom_P(A,atoms)) return(0);
  }
  return(1);
}

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

int Atom_Covers_Poverlap(Atms *A, Pover *O)
{
  Circ *cr=NULL;
  for(cr=First_Circ(O->Cr);cr!=NULL;cr=cr->next)
  {
    if(seperate(cr->I[0],A->tv)>A->tradius) return(0);
    if(seperate(cr->I[1],A->tv)>A->tradius) return(0);
  }
  return(1);
}

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

int Find_Pover_Without_Atom(Pover *po[MAX_OVER], Atms *at, int o[2]
                          ,int order, int start)
{
  int i,a;
  Circ *cr=NULL;
  
  for(i=start;iCr);cr!=NULL;cr=cr->next)
    {                                                  /* find extra atom */
      if(at==cr->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_Povers(Pover *po[MAX_OVER], int o[2])
{
  Pover *temp;
  temp=po[o[1]];
  po[o[1]]=po[o[0]];
  po[o[0]]=temp;
}

/**************************************************************************/
/* Remove redundant intersepts                                            */

Mint *Clean_Mints(Mint *mint, Atms *A[MAX_OVER], int order)
{
  Mint *mi=NULL;
  Mint *mp=NULL;
  int i;
  double sep=0.0;

/* check that all intersepts are within every atom */
  
  for(mi=First_Mint(mint);mi!=NULL;)
  {
    for(i=0;itv,mi->I)-A[i]->tradius)>0.0)
       &&(!AlmostZero(sep))) break;
    if(iprev==NULL) mi=mint;
      else mi=mint->next;
    }
    else mi=mi->next;
  }

/* check for redundant intersepts */

  for(mi=First_Mint(mint);mi!=NULL;mi=mi->next)
  {
    for(mp=mi->next;mp!=NULL;)
    {
      if((((mp->A[0]==mi->A[0])&&(mp->A[1]==mi->A[1]))
        ||((mp->A[0]==mi->A[1])&&(mp->A[1]==mi->A[0])))
       &&(AlmostZero(seperate(mp->I,mi->I))))     /* identical intersepts */
      {
        mint=Close_Mint(mp);
        if(mint==NULL) break;
        else mp=mint->next;
      }
      else mp=mp->next;
    }
    if(mint==NULL) break;
  }
  return(mint);
}

/**************************************************************************/
/* Sort the mints to facilitate multi atom feature                        */

Mint *Sort_Mints(Mint *mint)
{
  Mint *m=NULL;
  Mint *new=NULL;
  Mint *smallest=NULL;
  double mag, cosphi, sinphi;

  for(m=First_Mint(mint);m!=NULL;m=m->next)
  {
    mag=Vmag(m->I);
    cosphi=m->I.x/mag;
    sinphi=m->I.y/mag;
    m->phi=arcCos(cosphi);
    if(sinphi<0.0) m->phi=2*PI-m->phi;
  }
  while(mint!=NULL)
  {
    for(m=First_Mint(mint);m!=NULL;m=m->next)
    {
      if(smallest==NULL) smallest=m;
      else if(smallest->phi>m->phi) smallest=m;
    }
    if(smallest==NULL) break;
    if(mint==smallest)
    {
      if(smallest->prev!=NULL) mint=smallest->prev;
      else if(smallest->next!=NULL) mint=smallest->next;
      else mint=NULL;
    }
    if(smallest->next!=NULL) smallest->next->prev=smallest->prev;
    if(smallest->prev!=NULL) smallest->prev->next=smallest->next;
    smallest->next=NULL;
    smallest->prev=NULL;
    if(new==NULL) new=smallest;
    else
    {
      new->next=smallest;
      smallest->prev=new;
      new=smallest;
    }
    smallest=NULL;
  }
  return(new);
}

/**************************************************************************/
/* Sort the circles, determining their intersepts and there modes         */

int Sort_Circles_Using_Mints(Pover *over)
{
  int i,n,a,numseg=0;
  Circ *cr=NULL, *multi=NULL, *new=NULL;
  Mint *m[2]={NULL,NULL};
  Atms *com[2]={NULL,NULL};
  
  for(cr=First_Circ(over->Cr);cr!=NULL;cr=cr->next)
  {
    a=0;
    cr->mode=A_NORMAL;
    for(m[0]=First_Mint(over->Mi);m[0]!=NULL;m[0]=m[0]->next)
    {
      i=-1;
      if(m[0]->A[0]==cr->A) { i=0; n=1; }
      if(m[0]->A[1]==cr->A) { i=1; n=0; }
      if(i<0) continue;
      if(a<2) cr->I[a]=VequalV(m[0]->I);
      a++;
    }
    switch(a)
    {
      case 0 :                    /* multi atom situation on another atom */
        cr->mode=A_IGNORE;
        break;
      case 1 : 
        Error_Message(E_BDCIRC,"Sort Circles Using Mints");
        cr->mode=A_IGNORE;
        break;
      case 2 : break;    /* normal situation each atom has two intersepts */
      default:
        if((a/2)*2!=a)        /* should be multiple of two for multi atom */
        {
          Error_Message(E_BDCIRC,"Sort Circles Using Mints");
          cr->mode=A_IGNORE;
          break;
        }
        if(multi!=NULL)                  /* should be only one multi atom */
        {
          Error_Message(E_BDCIRC,"Sort Circles Using Mints");
          cr->mode=A_IGNORE;
          break;
        }
        numseg=a;
        multi=cr;
        cr->mode|=A_MULTI;
        break;
    }
    Get_Circle_Mode(cr,over->G);
  }

  /* deal with multi atom circle */
  
  if(multi!=NULL)
  {
    over->Mi=Sort_Mints(over->Mi);
    for(m[0]=First_Mint(over->Mi);m[0]!=NULL;m[0]=m[0]->next)
    {
      com[0]=NULL;
      com[1]=NULL;
      if(m[0]->next!=NULL) m[1]=m[0]->next;
      else m[1]=First_Mint(over->Mi);
      if((m[0]->A[0]==m[1]->A[0])
       ||(m[0]->A[0]==m[1]->A[1])) com[0]=m[0]->A[0];
      if((m[0]->A[1]==m[1]->A[0])
       ||(m[0]->A[1]==m[1]->A[1])) com[1]=m[0]->A[1];

      if((com[0]==NULL)&&(com[1]==NULL))
      {
        Error_Message(E_BDMINT,"Sort Circles Using Mints");
        continue;
      }
      if((com[0]!=NULL)&&(com[1]!=NULL)) continue;
      if((com[0]!=NULL)&&(com[1]==NULL)) i=0;
      if((com[0]==NULL)&&(com[1]!=NULL)) i=1;
      if(com[i]!=multi->A) Error_Message(E_BDMINT,"Sort Circles Using Mints");

      if((new=New_Circ(over->Cr))==NULL) break;
      Initialize_Circle(new,multi->A);
      new->I[0]=VequalV(m[0]->I);
      new->I[1]=VequalV(m[1]->I);
      new->mode=A_NORMAL|A_TEMP;
      over->Cr=new;
    }
  }
  
  /* count number of usable circles and calculate their info */

  for(i=0,cr=First_Circ(over->Cr);cr!=NULL;cr=cr->next)
  {
    if(cr->mode==A_IGNORE) continue;
    i++;
    cr->delta=seperate(cr->A->tv,over->G);
    if(!(cr->mode&A_MULTI)) Calculate_Circle_Intersept_Angles(cr,&over->G,&cr->A->tv);
  }
  over->ncir=i;
  return(over->ncir);
}

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

int Order_Previous_Poverlaps(Pover *po[MAX_OVER], Atms *at[MAX_OVER], int order)
{
  Atms *atom=NULL;
  Circ *cr=NULL;
  int i,o[2],a=0;
  
  if(order<3) return(0);             /* only triple or greater considered */
  for(i=0;iCr);(inext,i++)
  {                                      /* map po[0] to last (n-1) atoms */
    at[i]=cr->A;
    if(at[i]==NULL) return(0);
  }                                 /* find appropriate overlap for at[1] */
  if(!Find_Pover_Without_Atom(po,at[1],o,order,1)) return(0);
  o[0]=1;
  Swap_Povers(po,o);                             /* swap over 1 and o[1]  */
  for(cr=First_Circ(po[1]->Cr);cr!=NULL;cr=cr->next)
  {                                                    /* find extra atom */
    a=0;
    atom=cr->A;
    if(atom==NULL) return(0);
    if(at[1]==atom) return(0);        /* Find_Pover_Without_Atom() failed */
    for(i=2;i1) 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;iCr);cr!=NULL;cr=cr->next)
    {
      for(a=0;aA==at[a]) break;
      if(a==order) return(0);
    }
  }
  
  return(1);
}

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

int Find_Atomic_Intersepts_P(Pover *dover, Atms *A[MAX_OVER],
                             int a, int b,
                             Vector I[2], int inter[2])
{
  int i;
  double delta;
  Circ *dCr[2]={NULL,NULL};

  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_Pover(dover);dover!=NULL;dover=dover->next)
  {                                                    /* find ab overlap */
    if(((dCr[0]=Goto_Circ(dover->Cr,1))==NULL)
     ||((dCr[1]=Goto_Circ(dover->Cr,2))==NULL)
     ||(dCr[0]->A==NULL)||(dCr[1]->A==NULL))
      return(Error_Message(E_BDCIRC,"Find Atomic Intersepts P"));
    if(((dCr[0]->A==A[a])&&(dCr[1]->A==A[b]))
     ||((dCr[1]->A==A[a])&&(dCr[0]->A==A[b])))
    {
      inter[0]=1;      /* check that both intersepts are within all atoms */
      inter[1]=1;
      for(i=0;((iI[0],A[i]->tv)-A[i]->tradius;
        if((delta>0.0)&&(!AlmostZero(delta))) inter[0]=0;
        delta=seperate(dCr[0]->I[1],A[i]->tv)-A[i]->tradius;
        if((delta>0.0)&&(!AlmostZero(delta))) inter[1]=0;
      }
      if(inter[0]) I[0]=VequalV(dCr[0]->I[0]);
      if(inter[1]) I[1]=VequalV(dCr[0]->I[1]);
      return(1);
    }
  }
  return(0);
}

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

int Get_All_Relevant_Intersept_Vectors(Pover *over, int order,
                                       Pover *O[MAX_OVER], Atms *A[MAX_OVER])
{
  int i=0;
  Mint *mint=NULL;
  Mint *new=NULL;

  for(i=0;iMi);mint!=NULL;mint=mint->next)
    {
      if((new=New_Mint(over->Mi))==NULL) return(0);
      Copy_Mint(new,mint);
      over->Mi=new;
    }
  }

  if((over->Mi=Clean_Mints(over->Mi,A,order))==NULL) return(0);

/* calculate G as average of intersept vectors */

  over->G=Vequal(0.0,0.0,0.0);
  for(i=0,mint=First_Mint(over->Mi);mint!=NULL;mint=mint->next,i++)
    over->G=Vsum(over->G,mint->I);
  over->G=SxV(1.0/(double)i,over->G);

  if(!Sort_Circles_Using_Mints(over)) return(0);

  return(1);
}

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

Pover *New_Multiple_Poverlap(Pover *over, int order
                          ,Pover *O[MAX_OVER], Atms *A[MAX_OVER])
{
  int i=0;
  Circ *cr=NULL;

  for(i=0;iCr))==NULL) return(Close_Current_Pover(over));
    over->Cr=cr;
    Initialize_Circle(cr,A[i]);
  }

  if(!Get_All_Relevant_Intersept_Vectors(over,order,O,A))
    return(Close_Current_Pover(over));
  
  return(over);
}

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

double Multiple_Poverlap_Area(Pover *over, int order, unsigned mode)
{
  double totgamma=0.0;
  double area=0.0;
  double multi_area=0.0;                   /* extra area for A_GCENT mode */
  char line[LINELEN];
  Circ *cr=NULL;
#ifdef DEBUG
  char type[LINELEN];
  type[0]=0;
#endif

/* check for excessively large theta range                                  */

  for(cr=First_Circ(over->Cr);cr!=NULL;cr=cr->next)
    totgamma+=cr->gamma;
  totgamma-=2.0*PI;
  if(!AlmostZero(totgamma))
  {
    if(mode&VIS_BIT)
    {
      sprintf(line,"Multiple Poverlap Area [%7.4f]",totgamma);
      Error_Message(E_GNOTPI,line);
    }
  }

/* calculate area from each of the circles present                        */

  for(cr=First_Circ(over->Cr);cr!=NULL;cr=cr->next)
  {
    if(cr->mode&A_IGNORE) continue;
    if(cr->mode&A_MULTI) continue;
#ifdef DEBUG
    if(cr->mode&A_TEMP) strcat(type,"|m|");
    else if(cr->mode&A_GCENT) strcat(type,"|c|");
    else strcat(type,"|n|");
#endif
    if(cr->mode&A_GCENT)
      multi_area=cr->theta/(2*PI)*Single_Atom_Area(cr->A,0);
    else
      area+= 0.5 * cr->radius * (cr->theta * cr->radius
                 - cr->delta * (sin(cr->phi) + sin(cr->theta-cr->phi)));
  }

/* display answers if in normal total calculation                         */

  if(mode&VIS_BIT)
  {
    sprintf(line,"[%f] ",area);
#ifdef DEBUG
    strcat(line,type);
#endif
    Out_Message(line,O_BLANK);
  }

  if(area<0.0)
  {
    Error_Message(E_NEGARE,"Multiple Poverlap Area");
    return(0.0);
  }
  return(area);
}

/**************************************************************************/
/* find nth order overlap and calculate the relevant area                 */

double Calc_Multiple_Poverlap(Pover *po[MAX_OVER], Pover **over,
                              int order, unsigned mode,
                              char *multiple, double area)
{
  Pover *new=NULL;
  char line[LINELEN];
  char atmstr[25];
  char atomsline[LINELEN];
  Atms *A[MAX_OVER];
  Pover *O[MAX_OVER];
  int i;
  short sign=1;
  
  /* initialize */
  
  if(order/2*2==order) sign=1; else sign=-1;
  for(i=0;iname,Get_Atom_Number(A[i]));
          strcat(atomsline,atmstr);
        }
        atomsline[strlen(atomsline)-1]=0;
        if(sign>0) sprintf(line,"+ %s[%d](%s) ",multiple,new->ncir,atomsline);
        else sprintf(line,"- %s[%d](%s) ",multiple,new->ncir,atomsline);
        Out_Message(line,O_BLANK);
      }
      if(new->area==0.0) 
        area+=sign*Multiple_Poverlap_Area(new,new->ncir,mode);
      else area+=sign*new->area;
      if(mode&VIS_BIT)
      {
        sprintf(line,"= %f",area);
        Out_Message(line,O_NEWLN);
      }
    }
  }
  return(area);
}

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

double Multiple_Poverlap(Pover *pover, Pover **over,
                         int order, int level, unsigned mode,
                         char *multiple, double area)
{
  static Pover *po[MAX_OVER];
  int i;

  if(pover==NULL) return(area);               /* miss call               */
  if(level==0)                                 /* initialize              */
  {
    while(pover->prev!=NULL) pover=pover->prev;
    for(i=0;iorder)
  {
    Error_Message(E_LEVHI,"Multiple Poverlap");
    return(area);                             /* miss call               */
  }
  if(pover->order!=order)
  {
    Error_Message(E_BADORD,"Multiple Poverlap");
    return(area);                             /* miss call               */
  }
  
  for(po[level]=pover;po[level]!=NULL;pover=pover->next,po[level]=pover)
  {
    if(levelnext,over,order,1+level,mode,multiple,area);
    }
    else                                       /* final depth, calculate  */
    {
      area=Calc_Multiple_Poverlap(po,over,order,mode,multiple,area);
    }
  }
  return(area);
}

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

double New_Craig_Counting_P(Mol *M, unsigned mode)
{
  Pover *over[MAX_OVER];
  char multiple[10][MAX_OVER];
  double area=0.0;
  char line[LINELEN];
  Atms *at=NULL;
  int i=0;
#ifdef DEBUG
  Pover *o=NULL;
  Circ *c=NULL;
  int numo, numc;
#endif
  
  if(M->multi>MAX_OVER) M->multi=MAX_OVER;
  for(i=0;imulti;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) at->stat=at->stat|CALC_BIT;
  for(at=M->atoms;at!=NULL;at=at->next)
  {
    if((at->stat&MAIN_BIT)&&(Outside_All_Atoms_P(at)))
    at->stat=at->stat|CALC_BIT; else at->stat=at->stat&(FULL_BIT^CALC_BIT);
  }

/****************** add all single atom areas *****************************/

  if(mode&VIS_BIT) Out_Message("Calculating all single atom projected areas",O_NEWLN);
  area=Get_All_Single_Areas(M,mode);

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

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

#ifdef DEBUG
  for(numo=1,o=First_Pover(over[1]);o!=NULL;o=o->next,numo++)
  {
    for(numc=1,c=First_Circ(o->Cr);c!=NULL;c=c->next,numc++)
    {
      c->num=numc;
    }
    o->num=numo;
  }
#endif

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

  for(i=2;imulti;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);
    }
    area=Multiple_Poverlap(over[i-1],&over[i],i,0,mode,multiple[i],area);
    if(over[i]==NULL) break;
#ifdef DEBUG
    for(numo=1,o=First_Pover(over[i]);o!=NULL;o=o->next,numo++)
    {
      for(numc=1,c=First_Circ(o->Cr);c!=NULL;c=c->next,numc++)
      {
        c->num=numc;
      }
      o->num=numo;
    }
#endif
  }

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

  for(i=1;imulti;i++) Close_All_Poverlaps(over[i],multiple[i],mode);
   
  if (area<0.0)  { Error_Message(E_ATOSML,"New Craig Counting P"); area=0.0; }
  return(area);
}

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