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

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

#include "sterdefn.h"
#include "stermem.h"
#include "stertext.h"
#include "ryan.h"
#include "ryan_nsa.c"

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

double Convert_To_NSA(Mol *M, double eps, unsigned short mode)
{
/* craigs variables  */
  double solid=0.0;
  Atms *at=NULL;

/* ryans variables   */
  atomType atoms[maxAtoms] ;   /* Global instance of atoms. */
  int numAtoms  ;              /* Contains the number of atoms in atoms[] */
  int parity , i , j , index1 , index2;
  atomType ovAtom[maxOv] ;   /* Contains atoms being checked for overlap. */
  double SA;
  double tmp;

/**************** copy data from craigs format to ryans *******************/
  j = 0;

  for(i=0,at=First_Atom(M->atoms);at!=NULL;i++,at=at->next)
  {                                          /* loop through craigs atoms */
     
    if( fabs(at->SVangle) >= 0.000001 )
    {
     j++;
     atoms[j].coords.x = at->v.x ;
     atoms[j].coords.y = at->v.y ;
     atoms[j].coords.z = at->v.z ;
     atoms[j].vertex = at->SVangle;
    } 
}
  numAtoms=j; 
 
/******************** prepare for ryan calculation  ***********************/

  obs.x = 0.0;  /* repair 110195 */
  obs.y = 0.0;
  obs.z = 0.0;
 
  InitAtoms( atoms , numAtoms );
  parity=1;
  SA = 0;

/************************* do ryan calculation  ***************************/

  for(i=1;i<=numAtoms;i++)
  {
   tmp=SA;
   SA = SA + single_cone( atoms[i] );
   tmp=SA-tmp;
   printf("\n%d : %lf", i , tmp);
  }

  for(j=2;j<=numAtoms;j++)
  {

    parity=parity*(-1);

    GenPerms( numAtoms , j ) ;

    for(index1=1;index1<=permCount;index1++)
    {

      printf("\n");
      for(index2=1;index2<=j;index2++)
      {
       printf("%d ",perm[index1][index2]);
       ovAtom[index2]=atoms[ perm[index1][index2] ] ;
      }

      tmp=SA;
      SA = SA + nOverlap ( ovAtom , j )*parity  ;
      tmp=tmp - SA;
      printf(": %lf",fabs(tmp));
    }

  }

  printf("\nThe total solid angle is %lf.\n", SA);


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

  solid=SA;
  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 ...  *****************************************/
/**************************************************************************/