/**************************************************************************/
/**************************************************************************/
/**************************   "steric"   **********************************/
/**************************************************************************/
/*************     Program to calculate ligand cone    ********************/
/*************     angles as a measure of steric size  ********************/
/**************************************************************************/
/**************************************************************************/

/**************************************************************************/
/****************** overlap calculation routines **************************/
/**************************************************************************/
/******************        This module is        **************************/
/******************      system independant      **************************/
/**************************************************************************/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "sterdefn.h"      /* declaration of structures and globals       */
#include "stercomm.h"      /* definitions of all menu command options     */
#include "stercalc.h"      /* main calculation routines                   */
#include "craig.h"         /* solid angle calculation routines            */
#include "sterover.h"      /* steric overlap (SA and VA) calculations     */
#include "stermem.h"       /* dynamic memory management module            */
#include "stertext.h"      /* all text functions for text mode            */
#include "steraid.h"       /* additional functions needed                 */

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

int bond_allowed(Atms *A, Atms *B, char mode)
{
  Bond *a=NULL;
  if(mode)
  {
    for(a=First_Bond(A->bond);a!=NULL;a=a->next)
    {
      if((a->type)&&(a->atom==B)) return(0);
    }    
  }
  return(1);
}

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

int group_allowed(Atms *A, Atms *B, char mode)
{
  if((mode)&&(A->group!=0)&&(B->group!=0))
  {
    if(A->group!=B->group) return(0);
  }
  return(1);
}

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

double Total_Overlap(Atms *A, char mode)
{
  if(mode) return(Single_Atom_Solid_Angle(A,0));
  else return(A->SVangle);
}

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

double Two_Atom_Overlap(Atms *A[2], char mode, double eps, double chi)
{
  Over *over=NULL;
  double Overlap=0.0;
  if(mode)
  {
    if((over=New_Over(over))!=NULL)
    {                            /* create the memory for counted doubles */
      Initialize_Over(over,A,2);
      Setup_Two_Ellipses(over,A,BOTH);
      Overlap=Double_Overlap_Solid(over,eps,0);
      Close_Over(over);
    }
    return(Overlap);
  }
  else return(A[0]->SVangle+A[1]->SVangle-chi);
}

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

double Steric_Overlap(Mol *M, Set *set, unsigned short mode)
{
  char line[LINELEN];
  double Overlap=0.0,sval;
  double chi;
  Atms *A[2]={NULL,NULL};

  for (A[0]=First_Atom(M->atoms);A[0]!=NULL;A[0]=A[0]->next)
  {
    if ((A[0]->SVangle!=0.0)&&(A[0]->stat&MAIN_BIT))
    {
      for (A[1]=A[0]->next;A[1]!=NULL;A[1]=A[1]->next)
      {
        if ((A[0]!=A[1])&&(A[1]->SVangle!=0.0)&&(A[1]->stat&MAIN_BIT))
        {
          if((bond_allowed(A[0],A[1],set->mode&NBO_BIT))
           &&(group_allowed(A[0],A[1],set->mode&NGR_BIT)))
          {
            chi=VangleV(A[0]->v,A[1]->v);
            if (chi<(A[0]->SVangle+A[1]->SVangle))                   /* overlap */
            {
              if (chi<=fabs(A[0]->SVangle-A[1]->SVangle))
              {
                if (A[0]->SVangle>=A[1]->SVangle)
                sval=Total_Overlap(A[0],mode&SA_OV);
                else sval=Total_Overlap(A[1],mode&SA_OV);
                Overlap+=sval;                            /* total overlap */
                if (mode&VIS_BIT)
                {
                  sprintf(line,"+ totalOL[%s%d-%s%d]:(%f) = %f"
                              ,A[0]->name,Get_Atom_Number(A[0])
                              ,A[1]->name,Get_Atom_Number(A[1])
                              ,sval,Overlap);
                  Out_Message(line,O_NEWLN);
                }
              }
              else
              {
                sval=Two_Atom_Overlap(A,mode&SA_OV,set->eps,chi);
                Overlap+=sval;
                if (mode&VIS_BIT)
                {
                  sprintf(line,"+ overlap[%s%d-%s%d]:(%f) = %f"
                              ,A[0]->name,Get_Atom_Number(A[0])
                              ,A[1]->name,Get_Atom_Number(A[1])
                              ,sval,Overlap);
                  Out_Message(line,O_NEWLN);
                }
              }
            }
          }
        }
      }
    }
  }
  return(Overlap);
}

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