/*///////////////////////////////////////////////////////////////////////
//  Sieve - A tool for describing protein backbone geometry            //
//                                                                     //
//  Authors: Peter R{\o}gen & Robert Sinclair                          //
///////////////////////////////////////////////////////////////////////*/

/* 
Copyright (c)   The researchers and professors of the Technical 
University of Denmark. All Rights Reserved
 
Permission to use, copy, modify and distribute any part of this GI 
software for educational, research and non-profit purposes, without fee, 
and without a written agreement is hereby granted, provided that the 
above copyright notice, this paragraph and the following four paragraphs 
appear in all copies.
 
Those desiring to incorporate this Software into commercial products or 
use for commercial purposes should contact the Technology Transfer 
Office, The Secretariat of Research and Innovation, Technical University 
of Denmark, Anker Engelundsvej, Building 101 A, 2800 Kgs. Lyngby, 
DENMARK. Fax (+45) 45888040.
 
IN NO EVENT SHALL THE TECHNICAL UNIVERSITY OF DENMARK BE LIABLE TO ANY 
PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL 
DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE,
EVEN IF THE TECHNICAL UNIVERSITY OF DENMARK  HAS BEEN ADVISED OF THE 
POSSIBILITY OF SUCH DAMAGE.
 
THE  SOFTWARE PROVIDED HEREIN IS ON AN "AS IS" BASIS, AND THE TECHNICAL
UNIVERSITY OF DENMARK HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, 
UPDATES, ENHANCEMENTS OR MODIFICATIONS. THE TECHNICAL UNIVERSITY OF 
DENMARK MAKES NO REPRESENTATIONS AND EXTENDS NO WARRANTIES OF ANY KIND,
EITHER IMPLIED OR EXPRESS, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR 
THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, TRADEMARK OR
OTHER RIGHTS.

CITATIONS
For publication of results, please cite: 
P. R{\o}gen & R. Sinclair, Computing a new Family of Shape Descriptors 
for Protein Structures , J. Chem. Inf. Comp. Sci., 43, 1740-1747, 2003.

*/


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

int no_go;

char farrayo_name[1024];
FILE *farrayo;

int min_length=-1;
int max_length=900;

/* maximum protein length (fairly arbitrary; can be increased) */
#define NMAX 2048

/* maximum block size (fairly arbitrary) */
#define BMAX 16384

double global_invariant[30];
int n_proteins=0;

struct block
{
  int data[BMAX];
  struct block *next;
} *blocks;

struct xyz
{
  double x, y, z;
} protein[NMAX];
int N;

struct triangle
{
  double x[3];
  double y[3];
  double z[3];
  struct block *b, *jb, *kb;
  int i, j, k;
  int n, jn, kn;
  int jj, kk;
  struct triangle *next;
} *triangles;
int n_triangles;

struct patch
{
  double x[4];
  double y[4];
  double z[4];
};

struct area
{
  double x[13];
  double y[13];
  double z[13];
  int n;
};

#define IN 1
#define OUT 0
#define UNKNOWN 3

struct segment
{
  double x1, x2, ix[2];
  double y1, y2, iy[2];
  double z1, z2, iz[2];
  double cx, cy, cz;
  int s1, s2, status, n_intersections;
};

struct polygon
{
  struct segment seg[12];
  int n;
};

double perturbation=0.0;
long n_triplets;
long n_real_triplets;

long int clock_overflows;
clock_t last_clock_count, initial_clock_count;
double max_clock_count=0.0;

void find_max_clock_count()
{
  clock_t cl;

  initial_clock_count=last_clock_count=clock();
  clock_overflows=0L;

  if (sizeof(clock_t)==sizeof(short int))
    max_clock_count=USHRT_MAX+1.0;
  else if (sizeof(clock_t)==sizeof(int))
    max_clock_count=UINT_MAX+1.0;
  else if (sizeof(clock_t)==sizeof(long int))
    max_clock_count=ULONG_MAX+1.0;
  else
  {
    printf("Internal error in find_max_clock_count()\n");
    exit(0);
  }
}

void update_processor_time()
{
  clock_t count;

  count=clock();
  if (count<last_clock_count) clock_overflows++;
  last_clock_count=count;

  return;
}

double read_processor_time()
{
  update_processor_time();

  return ((clock_overflows*max_clock_count
          -initial_clock_count)+clock())/CLOCKS_PER_SEC;
}

void copy_polygon(struct polygon *p1, struct polygon *p2)
{
  struct segment *seg1, *seg2;
  int i, j;

  for (i=0; i<p1->n; i++)
  {
    seg1=p1->seg+i; seg2=p2->seg+i;

    seg2->x1=seg1->x1; seg2->x2=seg1->x2;
    seg2->y1=seg1->y1; seg2->y2=seg1->y2;
    seg2->z1=seg1->z1; seg2->z2=seg1->z2;

    seg2->cx=seg1->cx;
    seg2->cy=seg1->cy;
    seg2->cz=seg1->cz;

    seg2->n_intersections=seg1->n_intersections;
    seg2->status=seg1->status;
    seg2->s1=seg1->s1;
    seg2->s2=seg1->s2;

    for (j=0; j<seg1->n_intersections; j++)
    {
      seg2->ix[j]=seg1->ix[j];
      seg2->iy[j]=seg1->iy[j];
      seg2->iz[j]=seg1->iz[j];
    }
  }

  p2->n=p1->n;
}

void neg_polygon(struct polygon *p1, struct polygon *p2)
{
  struct segment *seg1, *seg2;
  int i, j;

  for (i=0; i<p1->n; i++)
  {
    seg1=p1->seg+i; seg2=p2->seg+i;

    seg2->x1=-seg1->x1; seg2->x2=-seg1->x2;
    seg2->y1=-seg1->y1; seg2->y2=-seg1->y2;
    seg2->z1=-seg1->z1; seg2->z2=-seg1->z2;

    seg2->cx=-seg1->cx;
    seg2->cy=-seg1->cy;
    seg2->cz=-seg1->cz;

    seg2->n_intersections=seg1->n_intersections;
    seg2->status=seg1->status;
    seg2->s1=seg1->s1;
    seg2->s2=seg1->s2;

    for (j=0; j<seg1->n_intersections; j++)
    {
      seg2->ix[j]=-seg1->ix[j];
      seg2->iy[j]=-seg1->iy[j];
      seg2->iz[j]=-seg1->iz[j];
    }
  }

  p2->n=p1->n;
}

void make_polygon(int i, int j, struct polygon *pol)
{
  struct segment *seg;
  double mx, my, mz, ux, uy, uz, vx, vy, vz;
  int k, l, m;
  int c1, c2, c3, c4, t;

  mx =pol->seg[0].x1=protein[j  ].x-protein[i  ].x;
  mx+=pol->seg[1].x1=protein[j+1].x-protein[i  ].x;
  mx+=pol->seg[2].x1=protein[j+1].x-protein[i+1].x;
  mx+=pol->seg[3].x1=protein[j  ].x-protein[i+1].x;

  my =pol->seg[0].y1=protein[j  ].y-protein[i  ].y;
  my+=pol->seg[1].y1=protein[j+1].y-protein[i  ].y;
  my+=pol->seg[2].y1=protein[j+1].y-protein[i+1].y;
  my+=pol->seg[3].y1=protein[j  ].y-protein[i+1].y;

  mz =pol->seg[0].z1=protein[j  ].z-protein[i  ].z;
  mz+=pol->seg[1].z1=protein[j+1].z-protein[i  ].z;
  mz+=pol->seg[2].z1=protein[j+1].z-protein[i+1].z;
  mz+=pol->seg[3].z1=protein[j  ].z-protein[i+1].z;

  for (k=0; k<4; k++)
  {
    seg=pol->seg+k;

    l=(k+1)%4;
    seg->x2=pol->seg[l].x1;
    seg->y2=pol->seg[l].y1;
    seg->z2=pol->seg[l].z1;

    c1=j+(((k>>1)^k)&1); c2=i+((k&2)>>1);
    c3=j+(((l>>1)^l)&1); c4=i+((l&2)>>1);

         if (c1==c2) { c2=c3; c3=c4; }
    else if (c1==c3) { c3=c4; }
    else if (c2==c3) { c3=c4; }

         if (c2<c1 && c2<c3) { t=c1; c1=c2; c2=t; }
    else if (c3<c1 && c3<c2) { t=c1; c1=c3; c3=t; }

    if (c2>c3) { t=c2; c2=c3; c3=t; }

    ux=protein[c2].x-protein[c1].x;
    uy=protein[c2].y-protein[c1].y;
    uz=protein[c2].z-protein[c1].z;

    vx=protein[c3].x-protein[c1].x;
    vy=protein[c3].y-protein[c1].y;
    vz=protein[c3].z-protein[c1].z;

    seg->cx=uy*vz-vy*uz;
    seg->cy=uz*vx-vz*ux;
    seg->cz=ux*vy-vx*uy;

    if (seg->cx*mx+seg->cy*my+seg->cz*mz<0.0)
    {
      seg->cx=-seg->cx;
      seg->cy=-seg->cy;
      seg->cz=-seg->cz;
    }

    seg->n_intersections=0;
    seg->status=UNKNOWN;
    seg->s1=UNKNOWN;
    seg->s2=UNKNOWN;
  }

  pol->n=4;
}

void points_in_convex(struct segment *s, struct polygon *p)
{
  struct segment *seg;
  int i;

  s->s1=IN;
  s->s2=IN;
  for (i=0; i<p->n; i++)
  {
    seg=p->seg+i;

    if (seg->x1==s->x1 && seg->y1==s->y1 && seg->z1==s->z1) s->s1=OUT;
    if (seg->x1==s->x2 && seg->y1==s->y2 && seg->z1==s->z2) s->s2=OUT;

    if (seg->cx*s->x1+seg->cy*s->y1+seg->cz*s->z1<=0.0) s->s1=OUT;
    if (seg->cx*s->x2+seg->cy*s->y2+seg->cz*s->z2<=0.0) s->s2=OUT;
  }

  return;
}

int midpoint_in_convex(struct segment *s, struct polygon *p)
{
  int i;

  for (i=0; i<p->n; i++)
    if (p->seg[i].cx*(s->x1+s->x2)
       +p->seg[i].cy*(s->y1+s->y2)
       +p->seg[i].cz*(s->z1+s->z2)<=0.0) return OUT;

  return IN;
}

int point_in_convex(double x, double y, double z, struct polygon *p)
{
  int i;

  for (i=0; i<p->n; i++)
  {
    if (p->seg[i].x1==x && p->seg[i].y1==y && p->seg[i].z1==z) return OUT;
    if (p->seg[i].cx*x+p->seg[i].cy*y+p->seg[i].cz*z<=0.0) return OUT;
  }

  return IN;
}

void intersect_segments(struct polygon *p1, int i, struct polygon *p2, int j)
{
  double d1, d2, d3, d4, ix, iy, iz, ix2, iy2, iz2, l;
  struct segment *s1, *s2;

  s1=p1->seg+i;
  s2=p2->seg+j;

  if ((s1->x1==s2->x1 && s1->y1==s2->y1 && s1->z1==s2->z1)
   || (s1->x1==s2->x2 && s1->y1==s2->y2 && s1->z1==s2->z2)
   || (s1->x2==s2->x2 && s1->y2==s2->y2 && s1->z2==s2->z2)
   || (s1->x2==s2->x1 && s1->y2==s2->y1 && s1->z2==s2->z1)) return;

  if ((s1->cx==-s2->cx && s1->cy==-s2->cy && s1->cz==-s2->cz)
   || (s1->cx==s2->cx && s1->cy==s2->cy && s1->cz==s2->cz)) return;

  d1=s1->x1*s2->cx+s1->y1*s2->cy+s1->z1*s2->cz;
  d2=s1->x2*s2->cx+s1->y2*s2->cy+s1->z2*s2->cz;
  d3=s2->x1*s1->cx+s2->y1*s1->cy+s2->z1*s1->cz;
  d4=s2->x2*s1->cx+s2->y2*s1->cy+s2->z2*s1->cz;

  if (d1*d2>=0.0 || d3*d4>=0.0) return;

  d1=fabs(d1);
  d2=fabs(d2);
  d3=fabs(d3);
  d4=fabs(d4);

  ix=d2*s1->x1+d1*s1->x2;
  iy=d2*s1->y1+d1*s1->y2;
  iz=d2*s1->z1+d1*s1->z2;

  ix2=d4*s2->x1+d3*s2->x2;
  iy2=d4*s2->y1+d3*s2->y2;
  iz2=d4*s2->z1+d3*s2->z2;

  if (ix*ix2+iy*iy2+iz*iz2<=0.0) return;

  l=1.0/sqrt(ix*ix+iy*iy+iz*iz);
  ix*=l; iy*=l; iz*=l;
  l=1.0/sqrt(ix2*ix2+iy2*iy2+iz2*iz2);
  ix2*=l; iy2*=l; iz2*=l;

  if ((ix-ix2)*(ix-ix2)+(iy-iy2)*(iy-iy2)+(iz-iz2)*(iz-iz2)>1e-20)
  {
    printf("Internal error in intersect_segments(): 1\n");
    no_go=1;
  }

  if (s1->n_intersections>=2)
  {
    printf("Internal error in intersect_segments(): 2\n");
    no_go=1;
  }
  s1->ix[s1->n_intersections  ]=ix;
  s1->iy[s1->n_intersections  ]=iy;
  s1->iz[s1->n_intersections++]=iz;

  if (s2->n_intersections>=2)
  {
    printf("Internal error in intersect_segments(): 3\n");
    no_go=1;
  }
  s2->ix[s2->n_intersections  ]=ix;
  s2->iy[s2->n_intersections  ]=iy;
  s2->iz[s2->n_intersections++]=iz;
}

void update_segment_status(struct polygon *p, struct polygon *p2)
{
  struct segment *seg;
  int i;

  for (i=0; i<p->n; i++)
  {
    seg=p->seg+i;

    if (seg->status==IN) continue;

    switch (seg->n_intersections)
    {
      case 0:
        if (midpoint_in_convex(p->seg+i,p2)) seg->status=IN;
        else seg->status=OUT;
        break;
      case 1:
        if (seg->s1==IN && seg->s2==OUT)
        {
          seg->x2=seg->ix[0];
          seg->y2=seg->iy[0];
          seg->z2=seg->iz[0];
          seg->status=IN;
        }
        else if (seg->s2==IN && seg->s1==OUT)
        {
          seg->x1=seg->ix[0];
          seg->y1=seg->iy[0];
          seg->z1=seg->iz[0];
          seg->status=IN;
        }
        else if (seg->s2==OUT && seg->s1==OUT)
        {
          if (point_in_convex(seg->x1+seg->ix[0],
                              seg->y1+seg->iy[0],
                              seg->z1+seg->iz[0],p2))
          {
            seg->x2=seg->ix[0];
            seg->y2=seg->iy[0];
            seg->z2=seg->iz[0];
            seg->status=IN;
          }
          else if (point_in_convex(seg->x2+seg->ix[0],
                                   seg->y2+seg->iy[0],
                                   seg->z2+seg->iz[0],p2))
          {
            seg->x1=seg->ix[0];
            seg->y1=seg->iy[0];
            seg->z1=seg->iz[0];
            seg->status=IN;
          }
          else seg->status=OUT;
        }
        else
        {
          printf("Internal error in update_segments(): 1\n");
          no_go=1;
        }
        break;
      case 2:
        seg->x1=seg->ix[0];
        seg->y1=seg->iy[0];
        seg->z1=seg->iz[0];
        seg->x2=seg->ix[1];
        seg->y2=seg->iy[1];
        seg->z2=seg->iz[1];
        seg->status=IN;
        break;
      default:
        printf("Internal error in update_segments(): 2\n");
        no_go=1;
    }
  }
}

int intersect_polygons(struct polygon *p1, struct polygon *p2, struct polygon *p3)
{
  struct segment *sp, *sp3;
  double t1, t2;
  int i, j;

  p3->n=0;

  t1=p1->seg[0].x1;
  for (i=1; i<p1->n; i++)
    if (t1*p1->seg[i].x1<=0.0) goto cont_1;
  t2=p2->seg[0].x1;
  for (i=1; i<p2->n; i++)
    if (t2*p2->seg[i].x1<=0.0) goto cont_1;
  if (t1*t2<0.0) return 0;

cont_1:
  t1=p1->seg[0].y1;
  for (i=1; i<p1->n; i++)
    if (t1*p1->seg[i].y1<=0.0) goto cont_2;
  t2=p2->seg[0].y1;
  for (i=1; i<p2->n; i++)
    if (t2*p2->seg[i].y1<=0.0) goto cont_2;
  if (t1*t2<0.0) return 0;

cont_2:
  t1=p1->seg[0].z1;
  for (i=1; i<p1->n; i++)
    if (t1*p1->seg[i].z1<=0.0) goto cont_3;
  t2=p2->seg[0].z1;
  for (i=1; i<p2->n; i++)
    if (t2*p2->seg[i].z1<=0.0) goto cont_3;
  if (t1*t2<0.0) return 0;

cont_3:
  for (i=0; i<p1->n; i++)
    points_in_convex(p1->seg+i,p2);

  for (i=0; i<p2->n; i++)
    points_in_convex(p2->seg+i,p1);

  for (i=0; i<p1->n; i++)
  {
    if (p1->seg[i].s1==IN && p1->seg[i].s2==IN) p1->seg[i].status=IN;
    else p1->seg[i].status=UNKNOWN;
  }

  for (i=0; i<p2->n; i++)
  {
    if (p2->seg[i].s1==IN && p2->seg[i].s2==IN) p2->seg[i].status=IN;
    else p2->seg[i].status=UNKNOWN;
  }

  for (i=0; i<p1->n; i++)
    if (p1->seg[i].status==UNKNOWN)
      for (j=0; j<p2->n; j++)
        if (p2->seg[j].status==UNKNOWN)
          intersect_segments(p1,i,p2,j);

  update_segment_status(p1,p2);
  update_segment_status(p2,p1);

  for (i=0; i<p1->n; i++)
  {
    sp=p1->seg+i;
    if (sp->status!=OUT)
    {
      sp3=p3->seg+p3->n++;
      sp3->x1=sp->x1;
      sp3->y1=sp->y1;
      sp3->z1=sp->z1;

      sp3->x2=sp->x2;
      sp3->y2=sp->y2;
      sp3->z2=sp->z2;

      sp3->cx=sp->cx;
      sp3->cy=sp->cy;
      sp3->cz=sp->cz;

      sp3->status=IN;
      sp3->n_intersections=0;
    }
  }

  for (i=0; i<p2->n; i++)
  {
    sp=p2->seg+i;
    if (sp->status!=OUT)
    {
      sp3=p3->seg+p3->n++;
      sp3->x1=sp->x1;
      sp3->y1=sp->y1;
      sp3->z1=sp->z1;

      sp3->x2=sp->x2;
      sp3->y2=sp->y2;
      sp3->z2=sp->z2;

      sp3->cx=sp->cx;
      sp3->cy=sp->cy;
      sp3->cz=sp->cz;

      sp3->status=IN;
      sp3->n_intersections=0;
    }
  }

  if (p3->n>12)
  {
    printf("Internal error in intersect_polygons()\n");
    no_go=1;
    return 0;
  }

  if (p3->n==0) return 0;
  return 1;
}

int feq(double x1, double y1, double z1, double x2, double y2, double z2)
{
  if ((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)<1e-24) return 1;
  return 0;
}

int make_area(struct polygon *p, struct area *a, int n_areas)
{
  struct segment *seg;
  struct area *an;
  struct polygon cp;
  double x0, y0, z0, x, y, z, l;
  int i, flag, ii_1, ii_2, ii_3;

  if (p->n<3) return 0;

  copy_polygon(p,&cp);

  an=a+n_areas;

  an->n=0;
  flag=0;
  for (i=0; i<p->n; i++)
  {
    seg=p->seg+i;

    if (seg->status==IN)
    {
      x0=an->x[an->n  ]=seg->x1;
      y0=an->y[an->n  ]=seg->y1;
      z0=an->z[an->n++]=seg->z1;
      x =an->x[an->n  ]=seg->x2;
      y =an->y[an->n  ]=seg->y2;
      z =an->z[an->n++]=seg->z2;
      seg->status=OUT;
      flag=1;
      break;
    }
  }
  if (flag==0) return 0;

  while (!(x==x0 && y==y0 && z==z0))
  {
    flag=0;
    for (i=0; i<p->n; i++)
    {
      seg=p->seg+i;

      if (seg->status==IN
       && seg->x1==x && seg->y1==y && seg->z1==z)
      {
        x=an->x[an->n  ]=seg->x2;
        y=an->y[an->n  ]=seg->y2;
        z=an->z[an->n++]=seg->z2;
        seg->status=OUT;
        flag=1;
        break;
      }
      if (seg->status==IN
       && seg->x2==x && seg->y2==y && seg->z2==z)
      {
        x=an->x[an->n  ]=seg->x1;
        y=an->y[an->n  ]=seg->y1;
        z=an->z[an->n++]=seg->z1;
        seg->status=OUT;
        flag=1;
        break;
      }
    }
    if (!flag)
    {
      copy_polygon(&cp,p);

      for (i=0; i<p->n; i++)
      {
        seg=p->seg+i;

        if (seg->status==IN)
        {
          x=seg->x1;
          y=seg->y1;
          z=seg->z1;
          l=1.0/sqrt(x*x+y*y+z*z);
          seg->x1*=l;
          seg->y1*=l;
          seg->z1*=l;

          x=seg->x2;
          y=seg->y2;
          z=seg->z2;
          l=1.0/sqrt(x*x+y*y+z*z);
          seg->x2*=l;
          seg->y2*=l;
          seg->z2*=l;

          if (feq(seg->x1,seg->y1,seg->z1,
                  seg->x2,seg->y2,seg->z2)) seg->status=OUT;
        }
      }

      an->n=0;
      flag=0;
      for (i=0; i<p->n; i++)
      {
        seg=p->seg+i;

        if (seg->status==IN)
        {
          x0=an->x[an->n  ]=seg->x1;
          y0=an->y[an->n  ]=seg->y1;
          z0=an->z[an->n++]=seg->z1;
          x =an->x[an->n  ]=seg->x2;
          y =an->y[an->n  ]=seg->y2;
          z =an->z[an->n++]=seg->z2;
          seg->status=OUT;
          flag=1;
          break;
        }
      }

      while (!feq(x,y,z,x0,y0,z0))
      {
        flag=0;
        for (i=0; i<p->n; i++)
        {
          seg=p->seg+i;

          if (seg->status==IN && feq(seg->x1,seg->y1,seg->z1,x,y,z))
          {
            x=an->x[an->n  ]=seg->x2;
            y=an->y[an->n  ]=seg->y2;
            z=an->z[an->n++]=seg->z2;
            seg->status=OUT;
            flag=1;
            break;
          }
          else if (seg->status==IN && feq(seg->x2,seg->y2,seg->z2,x,y,z))
          {
            x=an->x[an->n  ]=seg->x1;
            y=an->y[an->n  ]=seg->y1;
            z=an->z[an->n++]=seg->z1;
            seg->status=OUT;
            flag=1;
            break;
          }
        }
        if (!flag)
        {
          printf("Internal error in make_area()\n");
          no_go=1;
          return 0;
        }
      }

      if (an->n<4) return 0;

      return 1;
    }
  }

  if (an->n<4) return 0;

  for (i=0; i<an->n; i++)
  {
    x=an->x[i];
    y=an->y[i];
    z=an->z[i];
    l=1.0/sqrt(x*x+y*y+z*z);
    an->x[i]*=l;
    an->y[i]*=l;
    an->z[i]*=l;
  }

  flag=1;
  while (flag)
  {
    flag=0;
    for (ii_1=0; ii_1<an->n-1; ii_1++)
    {
      ii_2=(ii_1+1)%(an->n-1);
      if (feq(an->x[ii_1],an->y[ii_1],an->z[ii_1],
              an->x[ii_2],an->y[ii_2],an->z[ii_2]))
      {
        for (ii_3=ii_2; ii_3<an->n-1; ii_3++)
        {
          an->x[ii_3]=an->x[ii_3+1];
          an->y[ii_3]=an->y[ii_3+1];
          an->z[ii_3]=an->z[ii_3+1];
        }
        an->n--;
        flag=1;
      }
    }
  }

  if (an->n<4) return 0;

  return 1;
}

double spangle(double a[3], double b[3], double c[3], int *scary)
{
  register double ab, bc, ac, det_pos, det_neg, t, y, x;
  
  ab=a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
  bc=c[0]*b[0]+c[1]*b[1]+c[2]*b[2];
  ac=a[0]*c[0]+a[1]*c[1]+a[2]*c[2];
  det_pos=a[0]*b[1]*c[2]+b[0]*c[1]*a[2]+c[0]*a[1]*b[2];
  det_neg=a[0]*c[1]*b[2]+b[0]*a[1]*c[2]+c[0]*b[1]*a[2];

  if ((det_pos!=0.0 && det_neg!=0.0 &&
       fabs((det_pos-det_neg)/(fabs(det_pos)+fabs(det_neg)))<1e-6)
   || (ab*bc!=0.0 && ac!=0.0 && fabs((ab*bc-ac)/(fabs(ab*bc)+fabs(ac)))<1e-6))
    (*scary)=1;
  else
    (*scary)=0;

  t=fabs(det_pos-det_neg)+fabs(ab*bc-ac);
  if (t==0.0)
  {
    (*scary)=1;
    return 0.0;
  }

  return atan2((det_pos-det_neg)/t,(ab*bc-ac)/t);
}

double det(double a[3], double b[3], double c[3])
{
  register double det_pos, det_neg;

  det_pos=a[0]*b[1]*c[2]+b[0]*c[1]*a[2]+c[0]*a[1]*b[2];
  det_neg=a[0]*c[1]*b[2]+b[0]*a[1]*c[2]+c[0]*b[1]*a[2];

  return det_pos-det_neg;
}

double spherical_triangle_area(double da, double db, double dc)
{
  double a, b, c, p;

  a=2.0*asin(0.5*da);
  b=2.0*asin(0.5*db);
  c=2.0*asin(0.5*dc);

  p=0.5*(a+b+c);

  return 4.0*atan(sqrt(fabs(tan(0.5*p)*tan(0.5*(p-a))
                           *tan(0.5*(p-b))*tan(0.5*(p-c)))));
}

double compute_an_area(struct area *ar)
{
  int ii, ii_1, ii_2, scary;
  double a[3], b[3], c[3], sum, first_angle, tmp_angle, il;
  double area_1, mx, my, mz, x1, y1, z1, x2, y2, z2, da, db, dc, area_2;

  if (ar->n-1<3) return 0.0;

  sum=0.0;
  for (ii=0; ii<(ar->n-1); ii++)
  {
    ii_1=(ii+1)%(ar->n-1);
    ii_2=(ii+2)%(ar->n-1);

    a[0]=ar->x[ii  ]; a[1]=ar->y[ii  ]; a[2]=ar->z[ii  ];
    b[0]=ar->x[ii_1]; b[1]=ar->y[ii_1]; b[2]=ar->z[ii_1];
    c[0]=ar->x[ii_2]; c[1]=ar->y[ii_2]; c[2]=ar->z[ii_2];

    tmp_angle=spangle(a,b,c,&scary);
    if (ii==0) first_angle=tmp_angle;

    if (tmp_angle*first_angle<=0.0 || scary==1) goto difficult;

    sum+=tmp_angle;
  }

  if (sum>=0.0) area_1=fabs( 2.0*M_PI-sum);
           else area_1=fabs(-2.0*M_PI-sum);

  if (area_1<1e-6) goto difficult;

  return area_1;

difficult:

  mx=my=mz=0.0;
  for (ii=0; ii<(ar->n-1); ii++)
  {
    mx+=ar->x[ii]; my+=ar->y[ii]; mz+=ar->z[ii];
  }
  il=1.0/sqrt(mx*mx+my*my+mz*mz);
  mx*=il;
  my*=il;
  mz*=il;

  area_2=0.0;
  for (ii=0; ii<(ar->n-1); ii++)
  {
    ii_1=(ii+1)%(ar->n-1);

    x1=ar->x[ii  ]; y1=ar->y[ii  ]; z1=ar->z[ii  ];
    x2=ar->x[ii_1]; y2=ar->y[ii_1]; z2=ar->z[ii_1];
    da=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
    db=sqrt((mx-x2)*(mx-x2)+(my-y2)*(my-y2)+(mz-z2)*(mz-z2));
    dc=sqrt((x1-mx)*(x1-mx)+(y1-my)*(y1-my)+(z1-mz)*(z1-mz));

    area_2+=spherical_triangle_area(da,db,dc);
  }

  return fabs(area_2);
}

void compute_areas(int i, int j, int k, int l, int m, int n,
                   struct area *areas, int n_areas)
{
  double a;
  int jj;

  for (jj=0; jj<n_areas; jj++)
  {
    a=compute_an_area(areas+jj);

    global_invariant[third_order_inv_index(i,j,k,l,m,n)]+=sign(i,j)*sign(k,l)*sign(m,n)*a;
  } 
}

void compute_areas_2(int i, int j, int k, int l,
                     struct area *areas, int n_areas)
{
  double tmp;
  int jj, sij, skl;

  for (jj=0; jj<n_areas; jj++)
  {
    tmp=compute_an_area(areas+jj);

    sij=sign(i,j);
    skl=sign(k,l);
    if (j<k)
    {
      global_invariant[3]+=sij*skl*tmp;
      global_invariant[4]+=skl*tmp;
      global_invariant[5]+=sij*tmp;
      global_invariant[6]+=tmp;
    }
    else
    {
      if (j<l)
      {
	global_invariant[ 7]+=sij*skl*tmp;
	global_invariant[ 8]+=skl*tmp;
	global_invariant[ 9]+=sij*tmp;
	global_invariant[10]+=tmp;
      }
      else
      {
	global_invariant[11]+=sij*skl*tmp;
	global_invariant[12]+=skl*tmp;
	global_invariant[13]+=sij*tmp;
	global_invariant[14]+=tmp;
      }
    }
  } 
}

void compute_areas_1(int i, int j, struct area *areas, int n_areas)
{
  double a;
  int jj;
  
  for (jj=0; jj<n_areas; jj++)
  {
    a=compute_an_area(areas+jj);

    global_invariant[1]+=sign(i,j)*a;
    global_invariant[2]+=a;
  } 
}

int third_order_inv_index(int i, int j, int k, int l, int m, int n)
{
  int num_index_less_than_j=0, l_index;
 
  if (k<j)
  {
    num_index_less_than_j++;
    if (l<j)
      num_index_less_than_j++;
    if (m<j)
    {
      num_index_less_than_j++;
      if (n<j)
        num_index_less_than_j++;
    }
  }

  l_index=0;

  if (m<l)
  {
    if (l<n) l_index=1;
        else l_index=2;
  }

  return 3*num_index_less_than_j+l_index+15;
}

int sign(int i, int j)
{
  double a[3], b[3], c[3];
  
  a[0]=protein[i+1].x-protein[i].x;
  a[1]=protein[i+1].y-protein[i].y;
  a[2]=protein[i+1].z-protein[i].z;
  b[0]=protein[j+1].x-protein[j].x;
  b[1]=protein[j+1].y-protein[j].y;
  b[2]=protein[j+1].z-protein[j].z;
  c[0]=protein[i+1].x+protein[i].x-protein[j+1].x-protein[j].x;
  c[1]=protein[i+1].y+protein[i].y-protein[j+1].y-protein[j].y;
  c[2]=protein[i+1].z+protein[i].z-protein[j+1].z-protein[j].z;

  if (det(a,b,c)>=0.0) return 1;
  return -1;
}

void singlet(int i, int j)
{
  struct polygon p;
  struct area areas[1];
  int m;

  make_polygon(i,j,&p);
  for (m=0; m<4; m++) p.seg[m].status=IN;
  make_area(&p,areas,0);

  compute_areas_1(i,j,areas,1);
}

void doublet(int i, int j, int k, int l)
{
  struct polygon pij, tmp_pij, pkl_p, pkl_m, out_pol;
  struct area areas[2];
  int n_areas;

  n_areas=0;
  make_polygon(i,j,&pij);
  copy_polygon(&pij,&tmp_pij);
  make_polygon(k,l,&pkl_p);
  neg_polygon(&pkl_p,&pkl_m);

  intersect_polygons(&pij,&pkl_p,&out_pol);
  n_areas+=make_area(&out_pol,areas,n_areas);

  intersect_polygons(&tmp_pij,&pkl_m,&out_pol);
  n_areas+=make_area(&out_pol,areas,n_areas);

  if (n_areas==0) return;

  compute_areas_2(i,j,k,l,areas,n_areas);

  update_processor_time();
}

void triplet(int i, int j, int k, int l, int m, int n, int new)
{
  static struct polygon sta_pol_pp, sta_pol_pm, sta_pol_mp, sta_pol_mm;
  struct polygon tmp_pol_pp, tmp_pol_pm, tmp_pol_mp, tmp_pol_mm;
  struct polygon pol_1_p, pol_2_p, pol_2_m, pol_3, pol_tmp;
  struct polygon out_pol_pp, out_pol_pm, out_pol_mp, out_pol_mm;
  struct area areas[4];
  int n_areas;

  n_areas=0;
  n_triplets++;

  if (new)
  {
    update_processor_time();

    make_polygon(i,j,&pol_1_p);
    make_polygon(k,l,&pol_2_p);
    copy_polygon(&pol_1_p,&pol_tmp);
    neg_polygon(&pol_2_p,&pol_2_m);

    if (intersect_polygons(&pol_1_p,&pol_2_p,&sta_pol_pp))
      neg_polygon(&sta_pol_pp,&sta_pol_mm);
    else
      sta_pol_mm.n=0;

    if (intersect_polygons(&pol_tmp,&pol_2_m,&sta_pol_pm))
      neg_polygon(&sta_pol_pm,&sta_pol_mp);
    else
      sta_pol_mp.n=0;
  }
  if (sta_pol_pp.n==0 && sta_pol_pm.n==0) return;
  if (sta_pol_pp.n>0)
  {
    copy_polygon(&sta_pol_pp,&tmp_pol_pp);
    copy_polygon(&sta_pol_mm,&tmp_pol_mm);
  }
  else
    tmp_pol_pp.n=tmp_pol_mm.n=0;

  if (sta_pol_pm.n>0)
  {
    copy_polygon(&sta_pol_pm,&tmp_pol_pm);
    copy_polygon(&sta_pol_mp,&tmp_pol_mp);
  }
  else
    tmp_pol_pm.n=tmp_pol_mp.n=0;

  make_polygon(m,n,&pol_3);

  if (tmp_pol_pp.n>0)
  {
    copy_polygon(&pol_3,&pol_tmp);
    intersect_polygons(&tmp_pol_pp,&pol_tmp,&out_pol_pp);
    copy_polygon(&pol_3,&pol_tmp);
    intersect_polygons(&tmp_pol_mm,&pol_tmp,&out_pol_mm);

    n_areas+=make_area(&out_pol_pp,areas,n_areas);
    n_areas+=make_area(&out_pol_mm,areas,n_areas);
  }
  if (tmp_pol_pm.n>0)
  {
    copy_polygon(&pol_3,&pol_tmp);
    intersect_polygons(&tmp_pol_pm,&pol_tmp,&out_pol_pm);
    copy_polygon(&pol_3,&pol_tmp);
    intersect_polygons(&tmp_pol_mp,&pol_tmp,&out_pol_mp);

    n_areas+=make_area(&out_pol_pm,areas,n_areas);
    n_areas+=make_area(&out_pol_mp,areas,n_areas);
  }

  if (n_areas==0) return;

  n_real_triplets++;

  compute_areas(i,j,k,l,m,n,areas,n_areas);
}

void create_triangles(int n)
{
  double si, co, si_0, co_0, si_1, co_1, dth_0, dth_1, th_0, th_1;
  int i, j, cnt, nj, nj_0, nj_1, j_0, j_1;
  
  cnt=((int)(4.0*n));
  for (i=1; i<=(n-1); i++)
    cnt+=2*((int)(4.0*n*cos(M_PI*i/(n*2.0))));
  n_triangles=cnt;

  triangles=(struct triangle *)malloc(cnt*sizeof(struct triangle));
  if (triangles==(struct triangle *)NULL)
  {
    printf("Sorry - out of memory in create_triangles()\n");
    exit(0);
  }

  for (i=0; i<=(n-2); i++)
  {
    si_0=sin(M_PI*(i+0.0)/(n*2.0));
    co_0=cos(M_PI*(i+0.0)/(n*2.0));
    nj_0=(int)(4.0*n*cos(M_PI*(i+0.0)/(n*2.0)));
    dth_0=M_PI/(0.5*nj_0);
    th_0=0.0;

    si_1=sin(M_PI*(i+1.0)/(n*2.0));
    co_1=cos(M_PI*(i+1.0)/(n*2.0));
    nj_1=(int)(4.0*n*cos(M_PI*(i+1.0)/(n*2.0)));
    dth_1=M_PI/(0.5*nj_1);
    th_1=0.0;

    for (j_0=0, j_1=0; j_0<nj_0 || j_1<nj_1; )
    {
      if (th_0+dth_0<=th_1+dth_1)
      {
        triangles[--cnt].x[0]=co_0*cos(th_0+dth_0);
        triangles[  cnt].y[0]=co_0*sin(th_0+dth_0);
        triangles[  cnt].z[0]=si_0;
        triangles[  cnt].x[1]=co_0*cos(th_0);
        triangles[  cnt].y[1]=co_0*sin(th_0);
        triangles[  cnt].z[1]=si_0;
        triangles[  cnt].x[2]=co_1*cos(th_1);
        triangles[  cnt].y[2]=co_1*sin(th_1);
        triangles[  cnt].z[2]=si_1;

        th_0+=dth_0;
        j_0++;
      }
      else
      {
        triangles[--cnt].x[0]=co_1*cos(th_1+dth_1);
        triangles[  cnt].y[0]=co_1*sin(th_1+dth_1);
        triangles[  cnt].z[0]=si_1;
        triangles[  cnt].x[1]=co_0*cos(th_0);
        triangles[  cnt].y[1]=co_0*sin(th_0);
        triangles[  cnt].z[1]=si_0;
        triangles[  cnt].x[2]=co_1*cos(th_1);
        triangles[  cnt].y[2]=co_1*sin(th_1);
        triangles[  cnt].z[2]=si_1;

        th_1+=dth_1;
        j_1++;
      }
    }
  }

  si=sin(M_PI*(n-1.0)/(n*2.0));
  co=cos(M_PI*(n-1.0)/(n*2.0));
  nj=(int)(4.0*n*cos(M_PI*(n-1.0)/(n*2.0)));
  for (j=0; j<nj; j++)
  {
    triangles[--cnt].x[0]=0.0;
    triangles[  cnt].y[0]=0.0;
    triangles[  cnt].z[0]=1.0;
    triangles[  cnt].x[1]=co*cos(M_PI*(j+0.0)/(0.5*nj));
    triangles[  cnt].y[1]=co*sin(M_PI*(j+0.0)/(0.5*nj));
    triangles[  cnt].z[1]=si;
    triangles[  cnt].x[2]=co*cos(M_PI*(j+1.0)/(0.5*nj));
    triangles[  cnt].y[2]=co*sin(M_PI*(j+1.0)/(0.5*nj));
    triangles[  cnt].z[2]=si;
  }

  printf("There are %d triangles\n",n_triangles);
}

int proj_point_in_convex(double px, double py, double pz,
                         double *x, double *y, double *z, int n)
{
  register double pic_0, pic;
  register double cx, cy, cz;
  int i, j;

  if (px==x[0] && py==y[0] && pz==z[0]) return 0;

  cx=y[0]*z[1]-y[1]*z[0];
  cy=z[0]*x[1]-z[1]*x[0];
  cz=x[0]*y[1]-x[1]*y[0];
  pic_0=cx*px+cy*py+cz*pz;
  if (pic_0==0.0) return 0;

  for (i=1; i<n; i++)
  {
    if (px==x[i] && py==y[i] && pz==z[i]) return 0;

    j=(i+1)%n;
    cx=y[i]*z[j]-y[j]*z[i];
    cy=z[i]*x[j]-z[j]*x[i];
    cz=x[i]*y[j]-x[j]*y[i];
    pic=cx*px+cy*py+cz*pz;
    if (pic*pic_0<=0.0) return 0;
  }

  return 1;
}

int curves_intersect(double *ax, double *ay, double *az, int na, int i,
                     double *bx, double *by, double *bz, int nb, int j)
{
  register double cx, cy, cz, x1, y1, z1, x2, y2, z2;
  register double x3, y3, z3, x4, y4, z4;
  register int ii, jj;

  ii=(i+1)%na;
  jj=(j+1)%nb;

  x1=ax[i];
  y1=ay[i];
  z1=az[i];
  x2=ax[ii];
  y2=ay[ii];
  z2=az[ii];
  cx=y1*z2-y2*z1;
  cy=z1*x2-z2*x1;
  cz=x1*y2-x2*y1;
  if ((bx[j]*cx+by[j]*cy+bz[j]*cz)*(bx[jj]*cx+by[jj]*cy+bz[jj]*cz)>0.0) return 0;

  x3=bx[j];
  y3=by[j];
  z3=bz[j];
  x4=bx[jj];
  y4=by[jj];
  z4=bz[jj];
  cx=y3*z4-y4*z3;
  cy=z3*x4-z4*x3;
  cz=x3*y4-x4*y3;
  if ((ax[i]*cx+ay[i]*cy+az[i]*cz)*(ax[ii]*cx+ay[ii]*cy+az[ii]*cz)>0.0) return 0;

  if ((x1==x3 && y1==y3 && z1==z3)
   || (x2==x3 && y2==y3 && z2==z3)
   || (x1==x4 && y1==y4 && z1==z4)
   || (x2==x4 && y2==y4 && z2==z4)) return 0;

  return 1;
}

int overlap(double *ax, double *ay, double *az, int na,
            double *bx, double *by, double *bz, int nb)
{
  double t1, t2, t3, t4;
  int i, j, pos, neg;

  t1=ax[0];
  for (i=1; i<na; i++)
    if (t1*ax[i]<=0.0) goto Try_yz;
  t3=bx[0];
  for (i=1; i<nb; i++)
    if (t3*bx[i]<=0.0) goto Try_yz;
  t2=ay[0];
  for (i=1; i<na; i++)
    if (t2*ay[i]<=0.0) goto try_xz;
  t4=by[0];
  for (i=1; i<nb; i++)
    if (t4*by[i]<=0.0) goto try_xz;
  if ((t1*t3<0.0 && t2*t4>0.0)
   || (t1*t3>0.0 && t2*t4<0.0)) return 0;

try_xz:

  t1=ax[0];
  t2=az[0];
  for (i=1; i<na; i++)
    if (t2*az[i]<=0.0) goto further;
  t3=bx[0];
  t4=bz[0];
  for (i=1; i<nb; i++)
    if (t4*bz[i]<=0.0) goto further;
  if ((t1*t3<0.0 && t2*t4>0.0)
   || (t1*t3>0.0 && t2*t4<0.0)) return 0;

Try_yz:

  t1=ay[0];
  for (i=1; i<na; i++)
    if (t1*ay[i]<=0.0) goto further;
  t3=by[0];
  for (i=1; i<nb; i++)
    if (t3*by[i]<=0.0) goto further;

try_yz:

  t1=ay[0];
  t2=az[0];
  for (i=1; i<na; i++)
    if (t2*az[i]<=0.0) goto further;
  t3=by[0];
  t4=bz[0];
  for (i=1; i<nb; i++)
    if (t4*bz[i]<=0.0) goto further;
  if ((t1*t3<0.0 && t2*t4>0.0)
   || (t1*t3>0.0 && t2*t4<0.0)) return 0;

further:

  for (i=0; i<na; i++)
    if (proj_point_in_convex(ax[i],ay[i],az[i],bx,by,bz,nb)) return 1;

  for (i=0; i<nb; i++)
    if (proj_point_in_convex(bx[i],by[i],bz[i],ax,ay,az,na)) return 1;

  for (i=0; i<na; i++)
    for (j=0; j<nb; j++)
      if (curves_intersect(ax,ay,az,na,i,bx,by,bz,nb,j)) return 1;

  t1=t2=t3=0.0;
  for (i=0; i<na; i++) { t1+=ax[i]; t2+=ay[i]; t3+=az[i]; }
  if (proj_point_in_convex(t1,t2,t3,bx,by,bz,nb)) return 1;

  t1=t2=t3=0.0;
  for (i=0; i<nb; i++) { t1+=bx[i]; t2+=by[i]; t3+=bz[i]; }
  if (proj_point_in_convex(t1,t2,t3,ax,ay,az,na)) return 1;

  return 0;
}

int next_patch_j(struct triangle *t)
{
  t->jn--;
  if (t->jn<0) { t->jn=-1; return t->jj=N*32768+N; }
  if (t->j>=BMAX)
  {
    t->jb=t->jb->next;
    t->j=1;
    return t->jj=t->jb->data[0];
  }
  return t->jj=t->jb->data[t->j++];
}

int next_patch_k(struct triangle *t)
{
  t->kn--;
  if (t->kn<0) { t->kn=-1; return t->kk=N*32768+N; }
  if (t->k>=BMAX)
  {
    t->kb=t->kb->next;
    t->k=1;
    return t->kk=t->kb->data[0];
  }
  return t->kk=t->kb->data[t->k++];
}

void sieve()
{
  struct triangle *first_triangle, *last_triangle, *tri;
  struct block *current_block;
  struct patch p;
  int i, ii, j, k, l, m, jj, jk, kk, min_jj, min_kk;
  int tr_i, tr_j, tr_k, tr_l, tr_m, tr_n;
  int jk_old, jj_old, new;

  printf("N=%i\n",N);
  
  n_triplets=0L;
  n_real_triplets=0L;

  if (blocks==(struct block *)NULL)
  {
    current_block=blocks=(struct block *)malloc(sizeof(struct block));
    if (blocks==(struct block *)NULL)
    {
      printf("Sorry - out of memory in sieve(): 1\n");
      exit(0);
    }
    current_block->next=(struct block *)NULL;
  }
  else
    current_block=blocks;

  i=0;

  for (ii=0; ii<n_triangles; ii++)
  {
    triangles[ii].b=current_block;
    triangles[ii].i=i;
    triangles[ii].n=0;
    for (j=0; j<N; j++)
      for (k=j+2; k<N; k++)
      {
        p.x[0]=protein[k  ].x-protein[j  ].x;
        p.x[1]=protein[k+1].x-protein[j  ].x;
        p.x[2]=protein[k+1].x-protein[j+1].x;
        p.x[3]=protein[k  ].x-protein[j+1].x;

        p.y[0]=protein[k  ].y-protein[j  ].y;
        p.y[1]=protein[k+1].y-protein[j  ].y;
        p.y[2]=protein[k+1].y-protein[j+1].y;
        p.y[3]=protein[k  ].y-protein[j+1].y;

        p.z[0]=protein[k  ].z-protein[j  ].z;
        p.z[1]=protein[k+1].z-protein[j  ].z;
        p.z[2]=protein[k+1].z-protein[j+1].z;
        p.z[3]=protein[k  ].z-protein[j+1].z;

        if (overlap(triangles[ii].x,triangles[ii].y,triangles[ii].z,3,
                    p.x,p.y,p.z,4))
        {
          if (i>=BMAX)
          {
            if (current_block->next==(struct block *)NULL)
            {
              current_block->next=(struct block *)malloc(sizeof(struct block));
              if (current_block->next==(struct block *)NULL)
              {
                printf("Sorry - out of memory in sieve(): 2\n");
                exit(0);
              }
              current_block=current_block->next;
              current_block->next=(struct block *)NULL;
            }
            else
              current_block=current_block->next;

            i=0;
          }
          current_block->data[i++]=j*32768+k;
          triangles[ii].n++;
        }
      }
  }

  /* First order */
  for (j=0; j<N; j++)
    for (k=j+2; k<N; k++)
      singlet(j,k);

  /* Second order */
  for (j=0; j<N; j++)
    for (k=j+2; k<N; k++)
    {
      jk=j*32768+k;
      first_triangle=(struct triangle*)0;
      for (ii=0; ii<n_triangles; ii++)
      {
        triangles[ii].jb=triangles[ii].b;
        triangles[ii].j=triangles[ii].i;
        triangles[ii].jn=triangles[ii].n;
        for (i=0; i<triangles[ii].n; i++)
        {
          if (next_patch_j(triangles+ii)==jk)
          {
            if (first_triangle==(struct triangle*)0)
              first_triangle=triangles+ii;
            else
              last_triangle->next=triangles+ii;
            triangles[ii].next=(struct triangle*)0;
            last_triangle=triangles+ii;
            break;
          }
        }
      }

      if (first_triangle!=(struct triangle*)0)
      {
        jj=N*32768+N;
        for (tri=first_triangle; tri!=(struct triangle*)0; tri=tri->next)
        {
          l=next_patch_j(tri);
          if (l<jj) jj=l;
        }
        while (jj!=N*32768+N)
        {
          tr_i=jk/32768; tr_j=jk-tr_i*32768;
          tr_k=jj/32768; tr_l=jj-tr_k*32768;
          if (tr_i!=tr_k && tr_i!=tr_l && tr_j!=tr_k && tr_j!=tr_l)
          {
            doublet(tr_i,tr_j,tr_k,tr_l);
            if (no_go==1) return;
          }

          for (tri=first_triangle; tri!=(struct triangle*)0; tri=tri->next)
            if (tri->jj==jj)
              next_patch_j(tri);

          jj=N*32768+N;
          for (tri=first_triangle; tri!=(struct triangle*)0; tri=tri->next)
          {
            l=tri->jj;
            if (l<jj) jj=l;
          }
        }
      }
    }

  /* Third order */
  jk_old=jj_old=N*32768+N;
  for (j=0; j<N; j++)
    for (k=j+2; k<N; k++)
    {
      jk=j*32768+k;
      first_triangle=(struct triangle*)0;
      for (ii=0; ii<n_triangles; ii++)
      {
        triangles[ii].jb=triangles[ii].b;
        triangles[ii].j=triangles[ii].i;
        triangles[ii].jn=triangles[ii].n;
        for (i=0; i<triangles[ii].n; i++)
        {
          if (next_patch_j(triangles+ii)==jk && triangles[ii].jn>=2)
          {
            if (first_triangle==(struct triangle*)0)
              first_triangle=triangles+ii;
            else
              last_triangle->next=triangles+ii;
            triangles[ii].next=(struct triangle*)0;
            last_triangle=triangles+ii;
            break;
          }
        }
      }

      if (first_triangle!=(struct triangle*)0)
      {
        jj=kk=N*32768+N;
        for (tri=first_triangle; tri!=(struct triangle*)0; tri=tri->next)
        {
          l=next_patch_j(tri);
          tri->kb=tri->jb;
          tri->k=tri->j;
          tri->kn=tri->jn;
          m=next_patch_k(tri);
          if (l<jj) { jj=l; kk=m; }
          else if (l==jj && m<kk) { jj=l; kk=m; }
        }
        while (jj!=N*32768+N && kk!=N*32768+N)
        {
          tr_i=jk/32768; tr_j=jk-tr_i*32768;
          tr_k=jj/32768; tr_l=jj-tr_k*32768;
          tr_m=kk/32768; tr_n=kk-tr_m*32768;
          if (tr_i!=tr_k && tr_i!=tr_l && tr_i!=tr_m && tr_i!=tr_n
           && tr_j!=tr_k && tr_j!=tr_l && tr_j!=tr_m && tr_j!=tr_n
                                       && tr_k!=tr_m && tr_k!=tr_n
                                       && tr_l!=tr_m && tr_l!=tr_n)
          {
            if (jk_old!=jk || jj_old!=jj) new=1; else new=0;
            triplet(tr_i,tr_j,tr_k,tr_l,tr_m,tr_n,new);
            if (no_go==1) return;
            jk_old=jk; jj_old=jj;
          }
          for (tri=first_triangle; tri!=(struct triangle*)0; tri=tri->next)
            if (tri->jj==jj && tri->kk==kk)
              if (next_patch_k(tri)==N*32768+N)
              {
                next_patch_j(tri);
                tri->kb=tri->jb;
                tri->k=tri->j;
                tri->kn=tri->jn;
                if (next_patch_k(tri)==N*32768+N)
                {
                  tri->jj=tri->kk;
                  tri->j=tri->k;
                  tri->jn=tri->kn;
                }
              }
          jj=kk=N*32768+N;
          for (tri=first_triangle; tri!=(struct triangle*)0; tri=tri->next)
          {
            l=tri->jj;
            m=tri->kk;
            if (l<jj) { jj=l; kk=m; }
            else if (l==jj && m<kk) { jj=l; kk=m; }
          }
        }
      }
    }
  if (no_go==1) return;

  printf("There were %li possible triplets\n",n_triplets);
  printf("       and %li actual triplets\n",n_real_triplets);
}


void calculate_third_order_invariants()
{
  int i;
  
  global_invariant[0]=N+1;
  for (i=1;i<30;i++) global_invariant[i]=0.0;
  
  if (N>32700)
  {
    printf("Sorry, this protein is too long\n");
    exit(0);
  }

  sieve();
}

char line[1000];

int getatom(FILE *fi)
{
  char *cp;
  while (cp=fgets(line,999,fi),strstr(line,"ATOM")!=line && cp!=NULL);
  if (cp==(char *)NULL) return 0;
  return 1;
}

double rnd()
{
  return rand()/(0.5*RAND_MAX)-1.0;
}

void look_at_file(char *dir, char *fname)
{
  FILE *fi;
  double x, y, z;
  int nr, no=-1500, link_count=0, start_chain,i;
  char chain_name, new_chain_name;
  char dir_fname[1024], *cp;
  
  strcpy(dir_fname,dir);
  strcat(dir_fname,fname); 

  farrayo=fopen(farrayo_name,"r");
  if (farrayo!=(FILE *)NULL)
  {
    while (cp=fgets(line,999,farrayo),cp!=NULL)
      if (strstr(line,fname)==line)
      {
        printf("There is already an entry for <%s> in file <%s>\n",fname,farrayo_name);
        fclose(farrayo);
        return;
      }
    fclose(farrayo);
  }

  fi=fopen(dir_fname,"r");
  if (fi==(FILE *)NULL)
  {
    printf("Sorry - could not open file <%s>\n",dir_fname);
    return;
  }
  chain_name='a';

loop:
  if (getatom(fi)==0) 
  { 
    fclose(fi);
    if (link_count<max_length && link_count>min_length && link_count>0
     && no-start_chain+1-link_count<=3)
    {
      N=link_count-1;
      global_invariant[0]=link_count;
      no_go=0;
      calculate_third_order_invariants();
      if (no_go==1) return;
      n_proteins++;
      farrayo=fopen(farrayo_name,"a");
      if (farrayo==(FILE *)NULL)
      {
        printf("Sorry - could not open file <%s>\n",farrayo_name);
        return;
      }
      if (chain_name==' ') chain_name='0';
      fprintf(farrayo,"%s %1c %i %i",fname,chain_name,no-start_chain+1-link_count,link_count);
      for (i=1; i<30; i++) fprintf(farrayo," %15.10f",global_invariant[i]/(2.0*M_PI));
      fprintf(farrayo,"\n");
      fclose(farrayo);
      printf("Written new entry in <%s>\n",farrayo_name);
      return;
    }
    return;
  }
  sscanf(line+22,"%4i",&nr);
  sscanf(line+21,"%1c",&new_chain_name);
  if (nr!=no && line[13]=='C' && line[14]=='A')
  {
/*
      if(new_chain_name!=chain_name && chain_name!='a')
	{
          if (link_count<max_length && link_count>min_length && link_count>0)
	    {
	      printing pdb entry  fprintf(farrayo,"%c%c%c%c:%c  %i %i",line[72],line[73],line[74],line[75],chain_name,link_count,no-start_chain+1-link_count);
	      fprintf(farrayo,"%s %i %i",fname,link_count,no-start_chain+1-link_count);
	      N=link_count-1;
	      global_invariant[0]=link_count;
	      n_proteins++;
	      printf("%i\n",link_count);
	      calculate_third_order_invariants();
	    }  
	  link_count=0;
	}
*/
    link_count++;
    if (link_count>max_length-1) { fclose(fi); return; }
    no=nr;
    if (link_count==1) start_chain=no;
    sscanf(line+30,"%lf %lf %lf",&x,&y,&z);
    protein[link_count-1].x=x+perturbation*rnd();
    protein[link_count-1].y=y+perturbation*rnd();
    protein[link_count-1].z=z+perturbation*rnd();
    chain_name=new_chain_name;
  }
  goto loop;
}

void dirwalk(char *dir)
{
  char name[1024];
  struct dirent *dp;
  DIR *dfd;

  if ((dfd=opendir(dir))==NULL)
  {
    printf("Sorry - could not open directory <%s>\n",dir);
    return;
  }
  while ((dp=readdir(dfd))!=NULL)
  {
    if (strcmp(dp->d_name,".")==0 || strcmp(dp->d_name,"..")==0) continue;
    sprintf(name,"%s/%s",dir,dp->d_name);
    if (dp->d_name[strlen(dp->d_name)-4]=='.'
     && dp->d_name[strlen(dp->d_name)-3]=='p'
     && dp->d_name[strlen(dp->d_name)-2]=='d'
     && dp->d_name[strlen(dp->d_name)-1]=='b')
    {
      printf("%i: Now looking at <%s>\n",n_proteins,dp->d_name);
      look_at_file(dir,dp->d_name);
    }
  }
  closedir(dfd);
}

main(int argc, char *argv[])
{
  int n;
  
  if (((32067*32768+30767)/32768)!=32067
   || (32067*32768+30767)-32067*32768!=30767)
  {
    printf("Sorry, this machine's integers are not long enough\n");
    exit(0);
  }

  find_max_clock_count();

  if (argc<6)
  {
    printf("Please give min_length, max_length, n,\n");
    printf("input directory (with trailing \"/\") and output file name as arguments.\n");
    printf("A noise amplitude may follow these.\n");
    return;
  }

  if (argc==7)
  {
    sscanf(argv[6],"%lf",&perturbation);
  }
  printf("Noise amplitude=%f\n",perturbation);
  strncpy(farrayo_name,argv[5],1023);
  sscanf(argv[1],"%i",&min_length);
  printf("min_length=%i\n",min_length);
  sscanf(argv[2],"%i",&max_length);
  printf("max_length=%i\n",max_length);
  sscanf(argv[3],"%i",&n);
  printf("The triangulation parameter is %i\n",n);
  create_triangles(n);

  dirwalk(argv[4]);
  printf("30 measures have been calculated for %i protein strands.\n",n_proteins);
  printf("%20.3f seconds CPU time.\n",read_processor_time());
  fflush(stdout);
}