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

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

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

double Calculate_philimit(double a, double b, double c
             ,double ap, double bp, double cp)
{
  double xnp, ynp, xn, yn, xp, yp, deltan;
  double a2  = a*a;
  double ap2  = ap*ap;
  double b2  = b*b;
  double bp2 = bp*bp;
  double c2  = c*c;
  double cp2 = cp*cp;
  double HGQ;

  xn = 0;
  if (b < bp) yn = b;
  else yn = bp;

  do
  {
    yp = yn;
    xp = xn;
    deltan = 4*yn*((xn+c)/(a2*bp2) - (xn-cp)/(ap2*b2));

    xnp = (2*yn/deltan) *
        ((xn*xn - c2)/(a2*bp2) - (xn*xn - cp2)/(ap2*b2) + 1/bp2 - 1/b2);
    ynp = (2/deltan) *
        (((xn + c)/a2)*((yn*yn)/bp2 + 1) - ((yn*yn)/b2 + 1)*((xn-cp)/ap2) +
         ((xn+c)*(xn-cp)*(c+cp))/(a2*ap2));
    xn = xnp;
    yn = ynp;
    if (yn < 0) return(69);  /* total overlap */
  } while (!AlmostZero(xn-xp) || !AlmostZero(yn-yp));

  HGQ = arcTan(yn,xn);
  if (HGQ >= PI)
  {
    Error_Message(E_BDHGQ,"Calculate philimit");
    exit(0);
  }
  return(HGQ);
}

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

double Calculate_Ellipses(double *Aa, double *Ab, double *Ac
             ,double *Ba, double *Bb, double *Bc
             ,double alpha, double beta
             ,double Adelta, double Bdelta)
{
  double As, Bs;

  As=cos(2*alpha)+cos(2*Adelta);
  *Aa=sin(2*alpha)/As;
  *Ab=msqrt(2*pow(sin(alpha),2)/As,"msqrt -> Calculate Ellipses (Ab)");
  *Ac=sin(2*Adelta)/As;
  Bs=cos(2*beta)+cos(2*Bdelta);
  *Ba=sin(2*beta)/Bs;
  *Bb=msqrt(2*pow(sin(beta),2)/Bs,"msqrt -> Calculate Ellipses (Bb)");
  *Bc=sin(2*Bdelta)/Bs;
  return(Calculate_philimit(*Aa,*Ab,*Ac,*Ba,*Bb,*Bc));
}

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

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

void Setup_integrand(double a, double b, double c, int sign, int sign2)
{
  iA=a;
  iB=b;
  iC=c;
  isign=sign;
  isign2=sign2;
}
/**************************************************************************/

double integrand(double phi)
{
  double p,q;
  double s;
  double eta;

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

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

double Old_Leach_Double_Cone(Atms *A, Atms *B, double chi
                        ,double eps, unsigned short mode)
{
  double philimit=PI_2;
  double alpha, beta, gamma, Adelta, Bdelta;
  double Aa,Ab,Ac,Ba,Bb,Bc;
  double integral=0;
  double temp;
  double phibranch=PI_2;
  char line[100];
  Vector V[2]={{0.0,0.0,0.0},{0.0,0.0,0.0}};
  Vector I[2]={{0.0,0.0,0.0},{0.0,0.0,0.0}};

  alpha=A->SVangle;
  beta=B->SVangle;
  gamma=(alpha+beta+chi)/2;
  Adelta=gamma-alpha;
  Bdelta=gamma-beta;
  if (gamma==PI_2)
  {
    printf("[%f] ",2*PI);
    return(2*PI);
  }
  if (gamma>PI_2)
  {
    printf("[%f] ",4*PI);
    return(4*PI);
  }
  if(mode&MOVG_BIT)
  {
    V[0]=unit_vector(A->v);
    V[1]=unit_vector(B->v);
    Find_2atom_Intersepts(alpha, beta, V, I);
    I[0]=unit_vector(Vsum(I[0],I[1]));
    Adelta=VangleV(I[0],V[0]);
    Bdelta=VangleV(I[0],V[1]);
  }

  philimit=Calculate_Ellipses(&Aa,&Ab,&Ac,&Ba,&Bb,&Bc,alpha,beta,Adelta,Bdelta);

  Setup_integrand(Ba,Bb,Bc,POSITIVE,POSITIVE);
  integral+=qsimp(integrand,0.0,philimit,eps);
  if ((mode&SCOR_BIT)&&(Bc>Ba))
  {
    phibranch=arcTan(Bb,msqrt(Bc*Bc-Ba*Ba,"msqrt -> Old Leach Double Cone (B.branch)"));
    Setup_integrand(Ba,Bb,Bc,POSITIVE,POSITIVE);
    temp=qsimp(integrand,philimit,phibranch,eps);
    integral+=temp;
    Setup_integrand(Ba,Bb,Bc,POSITIVE,NEGATIVE);
    temp=-qsimp(integrand,philimit,phibranch,eps);
    integral+=temp;
  }

  Setup_integrand(Aa,Ab,Ac,NEGATIVE,POSITIVE);
  integral+=qsimp(integrand,philimit,PI,eps);
  if ((mode&SCOR_BIT)&&(Ac>Aa))
  {
    phibranch=PI-arcTan(Ab,msqrt(Ac*Ac-Aa*Aa,"msqrt -> Old Leach Double Cone (A.branch)"));
    Setup_integrand(Aa,Ab,Ac,NEGATIVE,POSITIVE);
    integral+=qsimp(integrand,phibranch,philimit,eps);
    Setup_integrand(Aa,Ab,Ac,NEGATIVE,NEGATIVE);
    integral-=qsimp(integrand,phibranch,philimit,eps);
  }
  integral=2*PI-2*integral;
  if(mode&VIS_BIT)
  {
    sprintf(line,"double(%s[%d]-%s[%d]):[%f] ",A->name,Get_Atom_Number(A),B->name,Get_Atom_Number(B),integral);
    Out_Message(line,O_BLANK);
  }
  return(integral);
}

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

double Craig_Counting(Mol *M, double eps, unsigned short mode)
{
  int count, i;
  double solid=0;
  double chi=0.0;
  Atms *A=NULL;
  Atms *B=NULL;
  char line[100];

  M->atoms=First_Atom(M->atoms);
  for(A=M->atoms,count=0;A!=NULL;A=A->next)  /* zero all atomic counts    */
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Craig Counting");
    A->count=0;
  }

  for(A=M->atoms,count=0;A!=NULL;A=A->next) /* count A though whole list  */
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Craig Counting");
    if((A->count==0)&&(A->SVangle!=0.0)&&(A->stat&1))
    {
      for(B=M->atoms,i=0;B!=NULL;B=B->next) /* count B though whole list  */
      {
        i++; if(i>M->num_atoms) Error_Message(E_ATMNUM,"Craig Counting");
        if ((A!=B)&&(B->SVangle!=0.0)&&(B->stat&1))
        {
          chi=VangleV(A->v,B->v);
          if (chi<(A->SVangle+B->SVangle))  /* A-B steric overlap found   */
          {
            if (chi<=fabs(A->SVangle-B->SVangle))
            {                                     /* total overlap found  */
              if(mode&VIS_BIT) Out_Message("+ ",O_BLANK);
              if (A->SVangle>=B->SVangle)                 /* A overlaps B */
              solid+=Single_Atom_Solid_Angle(A,mode);
              else solid+=Single_Atom_Solid_Angle(B,mode);/* B overlaps A */
              if(mode&VIS_BIT)
              {
                sprintf(line,"= %f",solid);
                Out_Message(line,O_NEWLN);
              }
            }
            else                                 /* partial overlap found */
            {
              if(mode&VIS_BIT) Out_Message("+ ",O_BLANK);
              solid+=Old_Leach_Double_Cone(A,B,chi,eps,mode);/*leach algorithm*/
              if(mode&VIS_BIT)
              {
                sprintf(line,"= %f",solid);
                Out_Message(line,O_NEWLN);
              }
            }
            A->count++;                    /* indicate A has been counted */
            B->count++;                    /* indicate B has been counted */
          }
        }
      }
    }
    if ((A->count==0)&&(A->SVangle!=0)&&(A->stat&1))
    {                                           /* no overlaps were found */
      if(mode&VIS_BIT) Out_Message("+ ",O_BLANK);
      solid+=Single_Atom_Solid_Angle(A,mode);
      if(mode&VIS_BIT)
      {
        sprintf(line,"= %f",solid);
        Out_Message(line,O_NEWLN);
      }
      A->count++;                               /* count only atom A once */
    }
  }

  for(A=M->atoms,count=0;A!=NULL;A=A->next) /* count A though whole list  */
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"Craig Counting");
    if ((A->count>1)&&(A->SVangle!=0))
    {
      if(mode&VIS_BIT)
      {
        sprintf(line,"- %d*",A->count-1);
        Out_Message(line,O_BLANK);
      }
      solid-=(A->count-1)*Single_Atom_Solid_Angle(A,mode);
      if(mode&VIS_BIT)
      {
        sprintf(line,"= %f",solid);
        Out_Message(line,O_NEWLN);
      }
      A->count=1;                           /* remove extra single atoms  */
    }
  }
  if (solid>4*PI) { Error_Message(E_STOBIG,"Craig Counting"); solid=4*PI; }
  if (solid<0.0) { Error_Message(E_STOSML,"Craig Counting"); solid=0.0; }
  return(solid);
}

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