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

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

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

#include "sterdefn.h"      /* declaration of structures and globals       */
#include "stercomm.h"      /* definitions of all menu command options     */
#include "sterfile.h"      /* file input and output functions             */
#include "stergrp.h"       /* file with group dependant functions         */
#include "stercalc.h"      /* main calculation routines                   */
#include "stermem.h"       /* dynamic memory management module            */
#include "stertext.h"      /* all text functions for text mode            */
#include "steraid.h"       /* additional functions needed                 */
#include "craig.h"         /* functions for craig solid angle             */
#include "proja.h"         /* functions for projected area                */
#include "leach.h"         /* functions for old leach solid angle         */
#include "sterover.h"      /* steric overlap (VAO and SAO) calculations   */

#ifdef _RYAN_
#include "ryan.h"          /* functions for ryan solid angle              */
#endif

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

void Initialize_Steric_Parameter(Ster *ster, char *name, char type
                                 ,double min, double max)
{
  strcpy(ster->name,name); /* steric parameter name                       */
  ster->type=type;         /* type of calculation performed               */
  ster->tot_val=0.0;       /* total value of steric parameter             */
  ster->err_val=0.0;       /* calculated maximum error, if any            */
  ster->tot_con=0.0;       /* conformer average steric parameter          */
  ster->err_con=0.0;       /* calculated maximum error of conf. average   */
  ster->max_val=0.0;       /* maximum value in profile                    */
  ster->pr_area=0.0;       /* area under radial profile                   */
  ster->peak_R=0.0;        /* radius at profile peak                      */
  ster->min=min;
  ster->max=max;           /* min and max profile range                   */
  ster->size=0;            /* array size                                  */
  ster->val=NULL;          /* pointer to profile array                    */
  ster->conf=NULL;         /* pointer to individual conformer memory      */
/* next and prev pointers have already been assigned, do not fiddle !!!   */
}

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

Ster *Get_Steric_Type(Ster *ster, char *name, char type, Set *set)
{
  Ster *current=NULL;
  if((ster!=NULL)&&(ster->type==type)) return(ster);
  for(current=First_Steric(ster);current!=NULL;current=current->next)
  {
    if(current->type==type) return(current);
  }
  if(current==NULL)
  {
    if((current=New_Steric(ster))!=NULL)
    {
      if(type==S_ANGU)
        Initialize_Steric_Parameter(current,name,type,set->a_min,set->a_max);
      else
        Initialize_Steric_Parameter(current,name,type,set->min,set->max);
    }
    else return(ster);
  }
  return(current);
}

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

Ster *Which_Ster(char arg, Ster *ster, Set *set)
{
  Ster *S=NULL;
  switch(arg)
  {
    case ANGU : S=Get_Steric_Type(ster,S_ANGUS,S_ANGU,set); break;
    case CONE : S=Get_Steric_Type(ster,S_CONES,S_CONE,set); break;
    case TOLM : S=Get_Steric_Type(ster,S_TOLMS,S_TOLM,set); break;
    case OLDL : S=Get_Steric_Type(ster,S_OLDLS,S_OLDL,set); break;
    case RYAN : S=Get_Steric_Type(ster,S_RYANS,S_RYAN,set); break;
    case CRAI : S=Get_Steric_Type(ster,S_CRAIS,S_CRAI,set); break;
    case NUME : S=Get_Steric_Type(ster,S_NUMES,S_NUME,set); break;
    case VAOV : S=Get_Steric_Type(ster,S_VAOVS,S_VAOV,set); break;
    case SAOV : S=Get_Steric_Type(ster,S_SAOVS,S_SAOV,set); break;
    case VOLM : S=Get_Steric_Type(ster,S_VOLMS,S_VOLM,set); break;
    case PROJ : S=Get_Steric_Type(ster,S_PROJS,S_PROJ,set); break;
    default   : return(NULL);
  }
  return(S);
}

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

void Set_Projected_XY(Mol *M, double theta, double phi)
{
  int n=0;
  Atms *at=NULL;
  Vector Y;

  Y=unit_vector(cross_product(M->plane,M->plane_x));
  for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Set Projected XY");
    if(at->stat&MAIN_BIT)
    {
      at->tv.x=cos(theta)*dot_product(at->v,M->plane_x);
      at->tv.y=cos(phi)*dot_product(at->v,Y);
      at->tv.z=0.0;
      at->tradius=at->radius;
    }
  }
}

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

void Set_Total_SVAngles(Mol *M)
{
  int n=0;
  Atms *at=NULL;
  double dis=0.0, rad=0.0;

  for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Set Total SVAngles");
    if(at->stat&MAIN_BIT)
    {
      dis=at->distance;
      rad=at->radius;
      if (rad==0) at->SVangle=0.0;
      else if (rad>dis) at->SVangle=PI;
      else if (rad==dis) at->SVangle=PI_2;
      else at->SVangle=asin(rad/dis);
    }
    else at->SVangle=0.0;
    if(at->SVangle>=PI) Error_Message(E_SVTOBG,"Set Total SVAngles");
  }
}

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

void Set_Radial_SVAngles(Mol *M, double Prad)
{
  Atms *at=NULL;
  int n=0;
  double dis=0.0, rad=0.0;

  for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Set Radial SVAngles");
    if(at->stat&MAIN_BIT)
    {
      dis=at->distance;
      rad=at->radius;
      if (rad==0) at->SVangle=0.0;
      else if ((Prad>(dis-rad))&&(Prad<(dis+rad)))
      at->SVangle=cosrule_angle(Prad,dis,rad);
      else at->SVangle=0.0;
    }
    else at->SVangle=0.0;
    if(at->SVangle>=PI) Error_Message(E_SVTOBG,"Set Radial SVAngles");
  }
}

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

int Get_Theta_Positions(Mol *M)
{
  Vector V={0.0,0.0,0.0};
  double D=0.0;
  Atms *at=NULL;
  int n=0;

  if((M==NULL)||(M->atoms==NULL)) return(0);
  if((M->main_atom!=NULL)
   &&(Vcmp(M->main_atom->v,Vequal(0.0,0.0,0.0))))
  {
    V=unit_vector(M->main_atom->v);
  }
  else 
  {
    for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
    {
      n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Get Theta Positions");
      if(!AlmostZero(at->SVangle))
      {
        V=Vsum(V,at->v);
        D+=at->distance;
      }
    }
    D/=(double)n;
    V=SxV(1.0/(double)n,V);
    if(Vmag(V)<D/3.0) V=First_Atom(M->atoms)->v;
    V=unit_vector(V);
  }
  for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Get Theta Positions");
    if(!AlmostZero(at->SVangle)) at->theta=VangleV(V,at->v);
  }
  M->basis_z=VequalV(V);
  return(1);
}

/**************************************************************************/
/*******************  cone angle aids  ************************************/
/**************************************************************************/

double Find_Total_Cone(Mol *M)
{
  Atms *at=NULL;
  double T=0.0;
  int n=0;

  for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Find Total Cone");
    if(!AlmostZero(at->SVangle)) T=fmax(T,at->SVangle+at->theta);
  }
  return(T);
}

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

double calc_cone(double alpha, double theta, double phi, double Pang)
{
  double T,D,a,b,c;
  double sin_cone,cos_cone;

  if(alpha<=0.0) return(0.0);
/*  if(phi>=PI) phi=phi-2*PI;  */
  if(AlmostZero(theta)) return(alpha);

  T=sin(theta)*(cos(phi)*cos(Pang)+sin(phi)*sin(Pang));
  a=T*T+cos(theta)*cos(theta);
  b=-2*T*cos(alpha);
  c=cos(alpha)*cos(alpha)-cos(theta)*cos(theta);
  D=b*b-4*a*c;
  if(D<0.0)
  {
    if(!AlmostZero(D)) return(0.0);
    else D=0.0;
  }
  else D=sqrt(D);

/* check that correct root is chosen out of two possible intersepts       */ 
  if(phi>Pang+PI) phi-=2*PI;
  if(phi<Pang-PI) phi+=2*PI;
  if(fabs(Pang-phi)>PI/2.0) D*=-1.0;

  sin_cone=(D-b)/(2*a);
  if(sin_cone<0.0) return(0.0);

  b=-2*cos(theta)*cos(alpha);
  c=cos(alpha)*cos(alpha)-T*T;
  D=b*b-4*a*c;
  if(D<0.0)
  {
    if(!AlmostZero(D)) return(0.0);
    else D=0.0;
  }
  else D=sqrt(D);
 
/* check that correct root is chosen out of two possible intersepts       */ 
  if(phi>Pang+PI) phi-=2*PI;
  if(phi<Pang-PI) phi+=2*PI;
  if(fabs(Pang-phi)>PI/2.0) D*=-1.0;

  cos_cone=(-D-b)/(2*a);

  return(acos(cos_cone));
}

/**************************************************************************/
/*******************  steric value calculations  **************************/
/**************************************************************************/

double Find_Min_Cone(Mol *M, double Pang)
{
  double cone=0.0;
  Atms *at=NULL;
  int n=0;

  for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Find Min Cone");
    cone=fmax(cone,calc_cone(at->SVangle,at->theta,at->phi,Pang));
  }
  return(cone);
}

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

double Single_Atom_Solid_Angle(Atms *atoms, unsigned mode)
{
  double solid;
  char line[100];
  if (atoms->SVangle==PI) solid=2*PI;
  else if (atoms->SVangle>PI) solid=4*PI;
  else solid=2*(PI)*(1-cos(atoms->SVangle));  /*<- single atom solid angle*/
  if(mode&VIS_BIT)
  {
    sprintf(line,"single(%s[%d]): [%f]",atoms->name,Get_Atom_Number(atoms),solid);
    Out_Message(line,O_BLANK);
  }
  return(solid);
}

/**************************************************************************/
/*******************  main menu options  **********************************/
/**************************************************************************/

double Cone(Mol *M, Ster *ster)
{
  ster->tot_val=2.0*Find_Total_Cone(M);
  ster->err_val=0.0;
  return(ster->tot_val);
}

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

double Tolman(Mol *M, Ster *ster)
{
  ster->tot_val=Find_Max_Group_SVangles(M);
  ster->err_val=0.0;
  return(ster->tot_val);
}

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

double Solid_Leach(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  ster->tot_val=Craig_Counting(M,set->eps,mode);
  ster->err_val=set->eps;
  return(ster->tot_val);
}

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

double Solid_Ryan(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  ster->tot_val=0.0;
#ifdef _RYAN_
  ster->tot_val=Convert_To_NSA(M,set->eps,mode);
#endif
  ster->err_val=set->eps;
  return(ster->tot_val);
}

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

double Solid_Craig(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  ster->tot_val=New_Craig_Counting(M,set->eps,mode);
  ster->err_val=set->eps;
  return(ster->tot_val);
}

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

double Molecular_Projection(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  ster->tot_val=New_Craig_Counting_P(M,mode);
  ster->err_val=set->eps;
  return(ster->tot_val);
}

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

double VA_Overlap(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  ster->tot_val=Steric_Overlap(M,set,mode);
  ster->err_val=0.0;
  return(ster->tot_val);
}

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

double SA_Overlap(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  ster->tot_val=Steric_Overlap(M,set,mode|SA_OV);
  ster->err_val=set->eps;
  return(ster->tot_val);
}

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

int Find_Encompassing_Box(Mol *M, int gnum, Vector *min, Vector *max, Vector *d)
{
  Vector average={0.0,0.0,0.0};
  Atms *A=NULL;
  double V=0.0;
  char line[LINELEN];
  int count;
  Grps *group=NULL;

  if(M==NULL) return(0);
  if((gnum)&&((group=Goto_Group(M->groups,gnum))==NULL)) gnum=0;
  if((gnum)&&(group->main==NULL)) gnum=0;

/* find centre of molecule */

  for(A=First_Atom(M->atoms),count=0;A!=NULL;A=A->next)
  {
    if(((A->stat&MAIN_BIT)&&(gnum==0))||(gnum==A->group))
    {
      average=Vsum(average,A->v);
      count++;
    }
  }
  if(count==0) return(0);
  average=SxV(1.0/(double)count,average);
  min->x=average.x;  max->x=average.x;
  min->y=average.y;  max->y=average.y;
  min->z=average.z;  max->z=average.z;

/* find encompassing box */

  for(A=First_Atom(M->atoms),count=0;A!=NULL;count++,A=A->next)
  {
    if(((A->stat&MAIN_BIT)&&(gnum==0))||(gnum==A->group))
    {
      if(min->x>A->v.x-A->radius) min->x=A->v.x-A->radius;
      if(max->x<A->v.x+A->radius) max->x=A->v.x+A->radius;

      if(min->y>A->v.y-A->radius) min->y=A->v.y-A->radius;
      if(max->y<A->v.y+A->radius) max->y=A->v.y+A->radius;

      if(min->z>A->v.z-A->radius) min->z=A->v.z-A->radius;
      if(max->z<A->v.z+A->radius) max->z=A->v.z+A->radius;
    }
  }
  if(gnum)                            /* make sure box encompasses sphere */
  {
    if(min->x>group->main->v.x-group->volrad) min->x=group->main->v.x-group->volrad;
    if(max->x<group->main->v.x+group->volrad) max->x=group->main->v.x+group->volrad;

    if(min->y>group->main->v.y-group->volrad) min->y=group->main->v.y-group->volrad;
    if(max->y<group->main->v.y+group->volrad) max->y=group->main->v.y+group->volrad;

    if(min->z>group->main->v.z-group->volrad) min->z=group->main->v.z-group->volrad;
    if(max->z<group->main->v.z+group->volrad) max->z=group->main->v.z+group->volrad;
  }

  d->x=max->x-min->x; d->y=max->y-min->y; d->z=max->z-min->z;
  V=d->x*d->y*d->z;

  sprintf(line,"average x=%f y=%f z=%f",average.x,average.y,average.z);
  Out_Message(line,O_NEWLN);
  sprintf(line,"min     x=%f y=%f z=%f",min->x,min->y,min->z);
  Out_Message(line,O_NEWLN);
  sprintf(line,"max     x=%f y=%f z=%f",max->x,max->y,max->z);
  Out_Message(line,O_NEWLN);
  sprintf(line,"D       x=%f y=%f z=%f",d->x,d->y,d->z);
  Out_Message(line,O_NEWLN);
  sprintf(line,"V=%f",V);
  Out_Message(line,O_NEWLN);
  
  return(1);
}

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

int Vector_In_Cell(Mol *M, Vector *p)
{
  Vector h={0.0,0.0,0.0};
  Atms *A=NULL;
  double v2=0.0;
  
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
  {
    if(strcmp(A->type,"A")==0)
    {
      v2=Vmag(A->v);
      v2*=v2;
      switch(A->name[1])
      {
        case '1' : h.x=dot_product(A->v,*p)/v2; break;
        case '2' : h.y=dot_product(A->v,*p)/v2; break;
        case '3' : h.z=dot_product(A->v,*p)/v2; break;
        default  : break;
      }
    }
  }
  if((h.x>=0.0)&&(h.x<1.0)
   &&(h.y>=0.0)&&(h.y<1.0)
   &&(h.z>=0.0)&&(h.z<1.0)) return(1);
  return(0);
}

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

int Vector_In_Sphere(Mol *M, Vector *p, int gnum)
{
  Grps *group=NULL;
  if(gnum<0) gnum*=-1;
  if(gnum==0) return(1);
  if((group=Goto_Group(M->groups,gnum))==NULL) return(1);
  if(AlmostZero(group->volrad)) return(1);
  if(group->main==NULL) return(1);
  if(Vmag(Vdif(group->main->v,*p))>group->volrad) return(0);
  return(1);
}

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

int Test_Atoms_Volume(Mol *M, int gnum, Vector *p, int count)
{
  Atms *A=NULL;

  if(gnum<0)
  {
    count++;
    for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
    {
      if(!A->stat&MAIN_BIT) continue;
      if(-1*gnum==A->group) continue;
      if(Vmag(Vdif(A->v,*p))<A->radius)
      {
        count--;
        break;
      }
    }
  }
  else if(gnum==0)
  {
    for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
    {
      if((A->stat&MAIN_BIT)
       &&(Vmag(Vdif(A->v,*p))<A->radius))
      {
        count++;
        break;
      }
    }
  }
  else
  {
    for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
    {
      if((A->stat&MAIN_BIT)
       &&(gnum==A->group)
       &&(Vmag(Vdif(A->v,*p))<A->radius))
      {
        count++;
        break;
      }
    }
  }
  return(count);
}

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

double Monte_Carlo_Volume(Mol *M, Set *set, int gnum
                         ,Vector *min, Vector *d, unsigned mode)
{
  Vector p={0.0,0.0,0.0};
  double V,v=0.0;
  long count=0,i=0,n=0,o=0,l=0;
  time_t t,tp;
  char line[LINELEN];
  double T_val=0.0;
  double tot_val=0.0;
  double dat_val=0.0;
  double cel_val=0.0;
  double sph_val=0.0;

  V=d->x*d->y*d->z;
  t=time(&t);
  tp=t;
  srand(t);
  while(i<set->maxmcv)
  {
    i++;
    p.x=min->x+((double)rand()/(double)RAND_MAX)*d->x;
    p.y=min->y+((double)rand()/(double)RAND_MAX)*d->y;
    p.z=min->z+((double)rand()/(double)RAND_MAX)*d->z;
    if((mode==VOL_UC)&&(!Vector_In_Cell(M,&p))) continue;
    o++;
    if((mode==VOL_RD)&&(!Vector_In_Sphere(M,&p,gnum))) continue;
    n++;
    count=Test_Atoms_Volume(M,gnum,&p,count);
    dat_val=V*(double)count/(double)i;
    t=time(&t);
    if(t-tp>set->tgap)
    {
      tp=t;
      sprintf(line,"n= %ld in N= %ld -> V= %f",count,i,dat_val);
      Out_Message(line,O_NEWLN);
    }
    if(i>set->avmmcv)
    {
      l++;
      T_val+=dat_val;
      v=tot_val;
      tot_val=T_val/(double)(l);
      if((i>set->minmcv)
       &&(!AlmostZero(tot_val))
       &&(fabs(v-tot_val)/tot_val<set->veps)) break;
    }
  }
  sprintf(line,"Data   n= %ld in N= %ld -> V= %f",count,i,tot_val);
  Out_Message(line,O_NEWLN);
  cel_val=V*(double)o/(double)i;
  sprintf(line,"Cell   n= %ld in N= %ld -> V= %f",o,i,cel_val);
  Out_Message(line,O_NEWLN);
  sph_val=V*(double)n/(double)i;
  sprintf(line,"Sphere n= %ld in N= %ld -> V= %f",n,i,sph_val);
  Out_Message(line,O_NEWLN);
  return(tot_val);
}

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

double Fixed_Grid_Volume(Mol *M, Set *set, int gnum, Vector *min, Vector *max
                        ,Vector *d, long *total, unsigned mode)
{
  Vector p={0.0,0.0,0.0};
  double V,v=0.0;
  long count=0,i=0,n=0,o=0;
  char line[LINELEN];
  double tot_val=0.0;
  double cel_val=0.0;
  double sph_val=0.0;

  V=d->x*d->y*d->z;
  d->x/=set->vx; d->y/=set->vy; d->z/=set->vz;
  v=d->x*d->y*d->z;
  
  count=0;
  for(i=0,p.x=min->x;p.x<=max->x;p.x+=d->x)
  {
    for(p.y=min->y;p.y<=max->y;p.y+=d->y)
    {
      for(p.z=min->z;p.z<=max->z;p.z+=d->z,i++)
      {
        if((mode==VOL_UC)&&(!Vector_In_Cell(M,&p))) continue;
        o++;
        if((mode==VOL_RD)&&(!Vector_In_Sphere(M,&p,gnum))) continue;
        n++;
        count=Test_Atoms_Volume(M,gnum,&p,count);
      }
    }
  }
  tot_val=V*(double)count/(double)i;
  *total=i;
  sprintf(line,"Data   n= %ld in N= %ld -> V= %f",count,i,tot_val);
  Out_Message(line,O_NEWLN);
  cel_val=V*(double)o/(double)i;
  sprintf(line,"Cell   n= %ld in N= %ld -> V= %f",o,i,cel_val);
  Out_Message(line,O_NEWLN);
  sph_val=V*(double)n/(double)i;
  sprintf(line,"Sphere n= %ld in N= %ld -> V= %f",n,i,sph_val);
  Out_Message(line,O_NEWLN);
  return(tot_val);
}

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

double Molecular_Volume(Mol *M, Ster *ster, Set *set, Grps *group)
{
  Vector min={0.0,0.0,0.0}, max={0.0,0.0,0.0}, d={0.0,0.0,0.0};
  double V;
  long total;
  double tot_val=0.0;
  double err_val=0.0;
  int gnum=0;
  unsigned mode=0;

  if(group!=NULL) gnum=Get_Group_Number(group);
  if((gnum>0)&&(group->volrad>0.0)) mode=VOL_RD;

/* find encompassing box */

  Find_Encompassing_Box(M,gnum,&min,&max,&d);
  if((group!=NULL)&&(group->mode&CAVV_BIT)) gnum*=-1;

/* choose two ways of calulating volume */

  if(set->mode&MTC_BIT)   /* Monte Carlo approach */
  {
    tot_val=Monte_Carlo_Volume(M,set,gnum,&min,&d,mode);
    err_val=set->veps;
  }
  else                    /* systemmatic approach */
  {
    tot_val=Fixed_Grid_Volume(M,set,gnum,&min,&max,&d,&total,mode);
    err_val=V/(double)total;
  }

  if(ster!=NULL)
  {
    ster->tot_val=tot_val;
    ster->err_val=err_val;
  }
  if(group!=NULL)
  {
    if(group->mode&CAVV_BIT)
    {
      if(group->mode&RADV_BIT) group->radcav=tot_val;
      else group->cavity=tot_val;
    }
    else group->volume=tot_val;
  }
  return(tot_val);
}

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

double Solid_Numerical(Mol *M, Ster *ster, Set *set, double Prad)
{
  int n=0;
  double cone=0.0;
  double Pang=0.0;
  double total[2]={0.0,0.0};
  double gap;

  gap=2*PI/(double)(set->n_size);
  for(n=0;n<set->n_size;n++)
  {
    Pang=gap*n;
    cone=Find_Min_Cone(M,Pang);
    if(set->contour!=NULL)
      fprintf(set->contour,"%10.6f %10.6f %10.6f\n",Prad,Pang,cone);
    total[0]+=(1.0-cos(cone))*gap;
    if((n/2)*2==n) total[1]+=(1.0-cos(cone))*2*gap;
  }
  if((set->contour!=NULL)&&(set->mode&PERS_BIT)) fprintf(set->contour,"\n");
  ster->tot_val=total[0];
  ster->err_val=fabs(total[1]-total[0]);
  return(ster->tot_val);
}

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

int Calc_Total_Steric(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  char perc='%';
  char line[LINELEN];

  sprintf(line,"Calculating total %s for %s ligand",ster->name,M->name);
  Out_Message(line,O_NEWLN);
  Set_Projected_XY(M,M->plane_T,M->plane_P);
  Set_Total_SVAngles(M);
  switch(mode)
  {
    case CONE : Cone(M,ster); break;
    case TOLM : Tolman(M,ster); break;
    case OLDL : Solid_Leach(M,ster,set,VIS_BIT|((MOVG_BIT|SCOR_BIT)&set->mode)); break;
    case RYAN : Solid_Ryan(M,ster,set,VIS_BIT); break;
    case CRAI : Solid_Craig(M,ster,set,VIS_BIT); break;
    case NUME : Solid_Numerical(M,ster,set,0.0); break;
    case VAOV : VA_Overlap(M,ster,set,VIS_BIT); break;
    case SAOV : SA_Overlap(M,ster,set,VIS_BIT); break;
    case VOLM : Molecular_Volume(M,ster,set,NULL); break;
    case PROJ : Molecular_Projection(M,ster,set,VIS_BIT); break;
    default   : return(0);
  }
  if((mode==CONE)||(mode==TOLM)||(mode==VAOV))
    sprintf(line,"Total %s for %s: %8.5f(%7.5f)rad %6.2fdeg %4.0f%c circle"
                ,ster->name,M->name
                ,ster->tot_val,ster->err_val,RtoD(ster->tot_val)
                ,100.0*ster->tot_val/(PI*2),perc);
  else if(mode==VOLM)
    sprintf(line,"Total %s for %s: %8.5f(%7.5f)cubic angstroms"
                ,ster->name,M->name
                ,ster->tot_val,ster->err_val);
  else if(mode==PROJ)
    sprintf(line,"Total %s for %s: %8.5f(%7.5f)square angstroms"
                ,ster->name,M->name
                ,ster->tot_val,ster->err_val);
  else
    sprintf(line,"Total %s for %s: %8.5f(%7.5f)sr %6.2fdeg %4.0f%c sphere"
                ,ster->name,M->name
                ,ster->tot_val,ster->err_val,StoD(ster->tot_val)
                ,100.0*ster->tot_val/(PI*4),perc);
  Out_Message(line,O_NEWLN);
  return(1);
}

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

int Calc_Radial_Profile(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  char perc='%';
  int n=0;
  double value=0.0;
  double area=0.0;
  double tot_old=0.0;
  double err_old=0.0;
  double Prad=0.0;
  double gap=0.0;
  char line[LINELEN];

  New_Profile(ster,set->size);
  if(!ster->size) return(0);
  if(set->mode&RSET_BIT)
  {
    ster->min=set->min;
    ster->max=set->max;
  }
  else if(set->mode&RDEF_BIT)
  {
    ster->min=M->minR;
    ster->max=M->maxR;
  }
  ster->max_val=0.0;

  sprintf(line,"Calculating the %s radial profile for %s ligand",ster->name,M->name);
  Out_Message(line,O_NEWLN);

  tot_old=ster->tot_val;
  err_old=ster->err_val;
  gap=(ster->max-ster->min)/(ster->size-1);
  for(n=0;n<ster->size;n++)
  {
    Prad=ster->min+gap*n;
    Set_Radial_SVAngles(M,Prad);
    switch(mode)
    {
      case CONE : value=Cone(M,ster); break;
      case TOLM : value=Tolman(M,ster); break;
      case OLDL : value=Solid_Leach(M,ster,set,(MOVG_BIT|SCOR_BIT)&set->mode); break;
      case RYAN : value=Solid_Ryan(M,ster,set,0); break;
      case CRAI : value=Solid_Craig(M,ster,set,0); break;
      case NUME : value=Solid_Numerical(M,ster,set,Prad); break;
      case VAOV : value=VA_Overlap(M,ster,set,0); break;
      case SAOV : value=SA_Overlap(M,ster,set,0); break;
      default   : return(0);
    }
    area+=value*gap;
    ster->val[n]=value;
    if(ster->max_val<value)
    {
      ster->max_val=value;
      ster->peak_R=Prad;
    }
    ster->max_val=fmax(ster->max_val,value);
    if((mode==CONE)||(mode==TOLM)||(mode==VAOV))
      sprintf(line,"\r%-3d, %-8s: %s at %6.3f Angstroms is %6.3f radians  ",
        ster->size-n-1,M->name,ster->name,Prad,value);
    else
      sprintf(line,"\r%-3d, %-8s: %s at %6.3f Angstroms is %6.3f steradians  ",
        ster->size-n-1,M->name,ster->name,Prad,value);
    Out_Message(line,O_CARET);
  }
  ster->tot_val=tot_old;
  ster->err_val=err_old;
  ster->pr_area=area;
  if((mode==CONE)||(mode==TOLM)||(mode==VAOV))
  {
    sprintf(line,"Maximum %s for %s was %8.5frad %6.2fdeg %4.0f%c circle"
                ,ster->name,M->name,ster->max_val,RtoD(ster->max_val)
                ,100.0*ster->max_val/(PI*2),perc);
    Out_Message(line,O_NEWLN);
    sprintf(line,"Area under %s profile for %s was %8.5frad.angstroms"
                ,ster->name,M->name,ster->pr_area);
    Out_Message(line,O_NEWLN);
  }
  else
  {
    sprintf(line,"Maximum %s for %s was %8.5fsr %6.2fdeg %4.0f%c sphere"
                ,ster->name,M->name,ster->max_val,StoD(ster->max_val)
                ,100.0*ster->max_val/(PI*4),perc);
    Out_Message(line,O_NEWLN);
    sprintf(line,"Area under %s profile for %s was %8.5fsr.angstroms"
                ,ster->name,M->name,ster->pr_area);
    Out_Message(line,O_NEWLN);
  }
  sprintf(line,"Radius at %s profile peak for %s was %8.5fangstroms"
              ,ster->name,M->name,ster->peak_R);
  Out_Message(line,O_NEWLN);
  if((!AlmostZero(ster->tot_val))&&(ster->max_val>ster->tot_val))
  {
    Error_Message(E_MTOBIG,"Calc Radial Profile");
    if(mode==OLDL) Error_Message(E_MTBOLD,"Calc Radial Profile");
  }
  if(ster->max_val<0.0) Error_Message(E_MTOSML,"Calc Radial Profile");
  return(1);
}

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

int Calc_Angular_Profile(Mol *M, Ster *ster, Set *set)
{
  char perc='%';
  int n=0;
  double cone=0.0;
  double Pang=0.0;
  double total[2]={0.0,0.0};
  double error=0.0;
  double gap=0.0;
  char line[LINELEN];

  New_Profile(ster,set->size);
  if(!ster->size) return(0);
  if(set->mode&RSET_BIT)
  {
    ster->min=set->a_min;
    ster->max=set->a_max;
  }
  else if(set->mode&RDEF_BIT)
  {
    ster->min=0.0;
    ster->max=2*PI;
  }
  ster->max_val=0.0;

  sprintf(line,"Calculating the %s angular profile for %s ligand",ster->name,M->name);
  Out_Message(line,O_NEWLN);

  Set_Total_SVAngles(M);

  gap=(ster->max-ster->min)/((double)ster->size);
  for(n=0;n<ster->size;n++)
  {
    Pang=ster->min+n*gap;
    cone=Find_Min_Cone(M,Pang);
    ster->val[n]=cone;
    ster->max_val=fmax(ster->max_val,cone);
    total[0]+=(1.0-cos(cone))*gap;
    if((n/2)*2==n) total[1]+=(1.0-cos(cone))*2*gap;
    sprintf(line,"\r%-3d, %-8s: Cone-angle at %6.3f Radians is %6.3f radians ",
	ster->size-n-1,M->name,Pang,cone);
    Out_Message(line,O_CARET);
  }
  error=fabs(total[1]-total[0]);
  sprintf(line,"Maximum %s for %s was %8.5frad %6.2fdeg %4.0f%c circle"
              ,ster->name,M->name,ster->max_val,RtoD(ster->max_val)
              ,100.0*ster->max_val/(PI*2),perc);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Total Solid Angle for %s: %8.5f(%7.5f)sr %4.0f%c sphere"
		,M->name,total[0],error,100.0*total[0]/(PI*4),perc);
  Out_Message(line,O_NEWLN);
  if((!AlmostZero(ster->tot_val))&&(ster->max_val>ster->tot_val))
    Error_Message(E_MTOBIG,"Calc Angular Profile");
  if(ster->max_val<0.0) Error_Message(E_MTOSML,"Calc Angular Profile");
  return(1);
}

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

int Calc_Projection_Map(Mol *M, Ster *ster, Set *set)
{
  char perc='%';
  int n=0;
  double cone=0.0;
  double Pang=0.0;
  double total[2]={0.0,0.0};
  double error=0.0;
  double gap=0.0;
  char line[LINELEN];

  New_Profile(ster,set->size);
  if(!ster->size) return(0);
  if(set->mode&RSET_BIT)
  {
    ster->min=set->a_min;
    ster->max=set->a_max;
  }
  else if(set->mode&RDEF_BIT)
  {
    ster->min=0.0;
    ster->max=2*PI;
  }
  ster->max_val=0.0;

  sprintf(line,"Calculating the %s angular profile for %s ligand",ster->name,M->name);
  Out_Message(line,O_NEWLN);

  Set_Total_SVAngles(M);

  gap=(ster->max-ster->min)/((double)ster->size);
  for(n=0;n<ster->size;n++)
  {
    Pang=ster->min+n*gap;
    cone=Find_Min_Cone(M,Pang);
    ster->val[n]=cone;
    ster->max_val=fmax(ster->max_val,cone);
    total[0]+=(1.0-cos(cone))*gap;
    if((n/2)*2==n) total[1]+=(1.0-cos(cone))*2*gap;
    sprintf(line,"\r%-3d, %-8s: Cone-angle at %6.3f Radians is %6.3f radians ",
	ster->size-n-1,M->name,Pang,cone);
    Out_Message(line,O_CARET);
  }
  error=fabs(total[1]-total[0]);
  sprintf(line,"Maximum %s for %s was %8.5frad %6.2fdeg %4.0f%c circle"
              ,ster->name,M->name,ster->max_val,RtoD(ster->max_val)
              ,100.0*ster->max_val/(PI*2),perc);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Total Solid Angle for %s: %8.5f(%7.5f)sr %4.0f%c sphere"
		,M->name,total[0],error,100.0*total[0]/(PI*4),perc);
  Out_Message(line,O_NEWLN);
  if((!AlmostZero(ster->tot_val))&&(ster->max_val>ster->tot_val))
    Error_Message(E_MTOBIG,"Calc Angular Profile");
  if(ster->max_val<0.0) Error_Message(E_MTOSML,"Calc Angular Profile");
  return(1);
}

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

int Calculate_Areas_OneD(Mol *M, Set *set, char mode)
{
  int n=0;
  double area=0.0;
  double angle=0.0;
  double maxang=0.0;
  double gap=0.0;
  char line[LINELEN];
  char name[LINELEN];
  Ster *ster=NULL;

  if((ster=Which_Ster(PROJ,M->ster,set))==NULL) return(0);
  M->ster=ster;
  if(mode==PHIP) New_Profile(ster,set->p_size);
  else           New_Profile(ster,set->t_size);
  if(!ster->size) return(0);
  ster->max_val=0.0;

  switch(mode)
  {
    case PHIP :
      ster->min=M->plane_minP;
      ster->max=M->plane_maxP;
      M->plane_P=M->plane_maxP;
      strcpy(name,"Phi");
      break;
    default :
      ster->min=M->plane_minT;
      ster->max=M->plane_maxT;
      M->plane_T=M->plane_maxT;
      strcpy(name,"Theta");
      break;
  }

  sprintf(line,"Calculating the %s vs %s 1D profile for %s ligand",ster->name,name,M->name);
  Out_Message(line,O_NEWLN);

  gap=(ster->max-ster->min)/((double)ster->size);
  for(n=0;n<ster->size;n++)
  {
    angle=ster->min+n*gap;
    if(mode==PHIP) Set_Projected_XY(M,M->plane_T,angle);
    else           Set_Projected_XY(M,angle,M->plane_P);
    area=Molecular_Projection(M,ster,set,0);
    ster->val[n]=area;
    if(area>ster->max_val)
    {
      ster->max_val=area;
      maxang=angle;
    }
    sprintf(line,"%-3d, %-8s: Projected area at %6.3f radians is %6.3f square angstroms ",
	ster->size-n-1,M->name,angle,area);
    Out_Message(line,O_CARET);
  }
  sprintf(line,"Maximum %s for %s was %8.5f square angstroms at %6.3f radians"
              ,ster->name,M->name,ster->max_val,maxang);
  Out_Message(line,O_NEWLN);
  if(ster->max_val<0.0) Error_Message(E_MTOSML,"Calculate Areas OneD");
  return(1);
}

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

int Calculate_Areas_TwoD(Mol *M, Set *set)
{
  int n=0,l=0;
  double area=0.0;
  double theta=0.0, phi=0.0;
  double maxtheta=0.0, maxphi=0.0;
  double gapT=0.0;
  double gapP=0.0;
  char line[LINELEN];
  Ster *ster=NULL;

  if((ster=Which_Ster(PROJ,M->ster,set))==NULL) return(0);
  M->ster=ster;
  New_Profile(ster,set->p_size*set->t_size);
  if(!ster->size) return(0);
  ster->min=M->plane_minT*M->plane_minP;
  ster->max=M->plane_maxT*M->plane_maxP;
  ster->max_val=0.0;

  sprintf(line,"Calculating the %s vs theta and phi 2D profile for %s ligand",ster->name,M->name);
  Out_Message(line,O_NEWLN);

  gapT=(M->plane_maxT-M->plane_minT)/((double)set->t_size);
  gapP=(M->plane_maxP-M->plane_minP)/((double)set->p_size);
  for(n=0;n<set->t_size;n++)
  {
    theta=M->plane_minT+n*gapT;
    for(l=0;l<set->p_size;l++)
    {
      phi=M->plane_minP+l*gapP;
      if(n*l>=ster->size) Error_Message(E_EXRANG,"Calculate Areas TwoD");
      Set_Projected_XY(M,theta,phi);
      area=Molecular_Projection(M,ster,set,0);
      ster->val[n]=area;
      if(area>ster->max_val)
      {
        ster->max_val=area;
        maxtheta=theta;
        maxphi=phi;
      }
      sprintf(line,"%-6d, %-8s: Projected area at (%6.3f,%6.3f) is %6.3f square angstroms ",
	ster->size-n*set->p_size-l-1,M->name,theta,phi,area);
      Out_Message(line,O_CARET);
    }
  }
  M->plane_T=theta;
  M->plane_P=phi;
  sprintf(line,"Maximum %s for %s was %8.5f square angstroms at %6.3f x %6.3f"
              ,ster->name,M->name,ster->max_val,maxtheta,maxphi);
  Out_Message(line,O_NEWLN);
  if(ster->max_val<0.0) Error_Message(E_MTOSML,"Calculate Areas TwoD");
  return(1);
}

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

int Calc_Conformer_Average(Mol *M, Ster *ster, Set *set, unsigned mode)
{
  Conf *first=NULL;
  Conf *current=NULL;
  Conf *prev=NULL;
  Atms *at=NULL;
  unsigned n,i,num;
  double Ec=0.0, value=0.0, temp=0.0, mole=0.0;
  double tot_tmp=0.0, err_tmp=0.0;
  char line[LINELEN];
  char Tname[50];
  unsigned fmode=F_CONFOR;
  FILE *TF;
  int Natoms;
  char s[30][9];
  Vector V={0.0,0.0,0.0};
  long fpos=0;

  if(M==NULL) return(0);
  if((at=First_Atom(M->atoms))==NULL) return(0);
  strcpy(Tname,M->Fname);
  New_Extension(Tname,".trj");
  if((TF=Open_Input_File(Tname,&fmode))==NULL)
  {
    sprintf(line,"Enter name of trajectory file : ");
    Get_Input_Line(line,set->input);
    if((TF=Open_Input_File(Tname,&fmode))==NULL)
    {
      Error_Message(E_CONFOP,"Calc Conformer Average");
      return(0);
    }
  }
  if(fmode!=F_CONFOR)
  {
    fclose(TF);
    Error_Message(E_NOTCNF,"Calc Conformer Average");
    return(0);
  }

  if(ster->conf!=NULL) ster->conf=Close_All_Conformers(ster->conf);
  get_next_line(TF,"Totmovatm:",line);
  sscanf(line,"%s%d",s[0],&Natoms);
  if(Natoms!=M->num_atoms)
  {
    Error_Message(E_CONFAT,"Calc Conformer Average");
    fclose(TF);
    return(0);
  }

  sprintf(line,"Calculating the individual conformer %ss",ster->name);
  Out_Message(line,O_NEWLN);

  tot_tmp=ster->tot_val;
  err_tmp=ster->err_val;
  if(M->origin!=NULL) V=VequalV(M->origin->v);
  for(i=1;(fpos=get_next_line(TF,"Time:",line))!=0;fseek(TF,fpos,0),i++)
  {
	if(fgets(line,148,TF)==NULL) break;
	if(fgets(line,148,TF)==NULL) break;
    if(sscanf(line,"%lf%lf",&temp,&Ec)<2) break;
    if(!get_next_line(TF,"Cordinat:",line)) break;

    for(num=0;num<Natoms;) /* get the next conformers atomic positions in */
    {
      if(Natoms-num>1) sscanf(line+9,"%d%lf%lf%lf%d%lf%lf%lf"
	    ,&num,&(Goto_Atom(at,num+1)->v.x),&(Goto_Atom(at,num+1)->v.y),&(Goto_Atom(at,num+1)->v.z)
	    ,&num,&(Goto_Atom(at,num+2)->v.x),&(Goto_Atom(at,num+2)->v.y),&(Goto_Atom(at,num+2)->v.z));
      else sscanf(line+9,"%d%lf%lf%lf"
	    ,&num,&(Goto_Atom(at,num+1)->v.x),&(Goto_Atom(at,num+1)->v.y),&(Goto_Atom(at,num+1)->v.z));
	  if((fgets(line,148,TF)==NULL)||(num>Natoms)) break;
    }
                         /* clear the atomic count for counting algorithm */

    for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next)
    {
      n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Calc Conformer Average");
      at->count=0;
    }
    if((at=First_Atom(M->atoms))==NULL) break;

                         /* calculate the relevant steric parameter value */

    if(M->origin) V=VequalV(M->origin->v);
    Calculate_Parameters(M,V,set);
    Set_Total_SVAngles(M);
    value=0.0;
    switch(mode)
    {
      case CONE : value=Cone(M,ster); break;
      case TOLM : value=Tolman(M,ster); break;
      case OLDL : value=Solid_Leach(M,ster,set,(MOVG_BIT|SCOR_BIT)&set->mode); break;
      case RYAN : value=Solid_Ryan(M,ster,set,0); break;
      case CRAI : value=Solid_Craig(M,ster,set,0); break;
      case NUME : value=Solid_Numerical(M,ster,set,0.0); break;
      case VAOV : value=VA_Overlap(M,ster,set,0); break;
      case SAOV : value=SA_Overlap(M,ster,set,0); break;
      default   : break;
    }
                  /* store steric value, conformer energy and temperature */

    if((current=New_Conformer(current))==NULL) break;
    if(ster->conf==NULL) ster->conf=current;
    current->value=value;
    current->Ec=Ec;
    current->temperature=temp;
    current->mol=0.0;
    sprintf(line,"%-3d, %-8s: value= %f  (Ec=%7.4f)"
	            ,i,M->name,value,Ec);
    Out_Message(line,O_NEWLN);
  }
  ster->tot_val=tot_tmp;
  ster->err_val=err_tmp;
  value=0.0;
  temp=0.0;
  sprintf(line,"Calculating the current molar ratios for %s ligand",M->name);
  Out_Message(line,O_NEWLN);
  if((first=First_Conformer(current))==NULL)
  {
    fclose(TF);
    return(0);
  }
  for(current=first,i=1;current!=NULL;current=current->next,i++)
  {
    mole=0.0;
    for(prev=first;prev!=NULL;prev=prev->next)
    {
      mole+=pow(E,(current->Ec-prev->Ec)/(GAS*current->temperature));
    }
    mole=1/(mole);
    temp+=mole;
    current->mol=mole;
    value+=current->value*mole;
    sprintf(line,"%-3d, %-8s: (mol=%8.6f)*(value=%8.6f)  =>  %8.6f sr "
	            ,i,M->name,mole,current->value,current->value*mole);
    Out_Message(line,O_NEWLN);
  }
  Out_Message("",O_NEWLN);
  ster->conf=first;
  ster->tot_con=value/temp;
  ster->err_con=0.0;

  sprintf(line,"Tot          : (value=%8.6f)/(mol=%8.6f)  =>  %8.6f sr"
	          ,value,temp,ster->tot_con);
  Out_Message(line,O_NEWLN);

  if((mode==CONE)||(mode==TOLM)||(mode==VAOV))
    sprintf(line,"Conformer Averaged %s for %s is %f radians (%f *2pi)\n"
	            ,ster->name,M->name,ster->tot_con,ster->tot_con/(2*PI));
  else
    sprintf(line,"Conformer Averaged %s for %s is %f steradians (%f *4pi)\n"
	            ,ster->name,M->name,ster->tot_con,ster->tot_con/(4*PI));
  Out_Message(line,O_NEWLN);
  return(1);
}


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