/* checkfit.c++   Moved routines used to check and cut the final fit
   for trace.c++ to here
*/

#include "drawwindow.h"
#include "math.h"
#include "trace.h"
#include "peptide.h"
#include "tordyn.h"
#include "checkfit.h"
#include "editbones.h"
#include "pdbtolev.h"
#include "assignseq.h"
#include "bone.h"
#include "string.h"
#include "stdlib.h"

#define MAX1D 1000  //maximum number of positions saved in routine avemax
#define MAXDEL 50 // Maximum nuber of fits deleted by avemaxbatch
extern struct tordyn tordata;
extern struct maiddatafile maidfile;
extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern struct extbones ext; // All the data for edited bones used to extend one end of a given fit
extern int fitnum[NFIT];// see trace.c++ for definition
extern int dynquit;  // if = 0, then use stand. quitcond - save best end density; if = N, just run N steps and save last 
extern int verbose; // if = 1,print some stuff;  = 2 print a lot of stuff
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version
extern float fminden[10]; // defined in bone.c++ and  mapmenu.c++


void closedist2(float pos1[3],float pos2[3],float capos[3],float *perpdist,
           float *linedist)
/*
Similar to closedist, but returns, perpdist (perp. dist from line)
    and linedist  (dist along line between to points of closest point
*/
{
   int k;
   float v[3],pv,rv,mag,pos[3],al;
   
   for(k=0;k<=2;k++) // vector from ipt to ipt+1
      v[k] = pos2[k] -pos1[k];
   mag = sqrt( dot(v,v) );
   for(k=0;k<=2;k++) // vector from ipt to ipt+1 
      v[k] = v[k]/mag;
   rv= dot(pos1,v);
   pv = dot(capos,v);
   al = pv-rv; // distance along unit vector v of closest point (should be >0)
   for(k=0;k<=2;k++) 
       pos[k] = pos1[k]+al*v[k];
   *perpdist = distv(pos,capos);
   *linedist = al*mag;
}


void closedist(float pos1[3],float pos2[3],float capos[3],float *dist,
           float pos[3],int *iplus)
/*
Given that close point lies between pos1 and pos2 on nmain,
  find closest dist = dist and pos
close point must be on line between pos1 and pos2
  (allow slightly less than pos1)
if close dist is closer to to top  point, return iplus =1, else = 0
*/
{
   int k;
   float v[3],al,pv,rv,mag;
   
   *iplus = 0;
   for(k=0;k<=2;k++) // vector from ipt to ipt+1
      v[k] = pos2[k] -pos1[k];
   mag = sqrt( dot(v,v) );
   for(k=0;k<=2;k++) // vector from ipt to ipt+1 
      v[k] = v[k]/mag;
   rv= dot(pos1,v);
   pv = dot(capos,v);
   al = pv-rv; // distance along unit vector v of closest point (should be >0)
   if( (al> -0.3)&&( (al-mag)< 0.3) ){ // close point must be with in 0.3 angstroms of  line between pos1 and pos2
      for(k=0;k<=2;k++) 
          pos[k] = pos1[k]+al*v[k];
      *dist = distv(pos,capos);
      if(al> 0.5*mag)// closer to top point
         *iplus = 1;
   }
   else
      *dist = 1000.0;
}

void linedist(float pos1[3],float pos2[3],float capos[3],float *dist,
           float pos[3],int *idist)
/*
New: 11sep01: if pos1=pos2 then just return idist=0
Find dist from capos to the line determined by pos1 and pos2
Define  v = pos2 - pos1.  and u = capos -pos1
Form right triangle  with hypotonuse = |u| base a (along v) and 
  height d (dist of point from line).
  Then u dot v = |u||v|cos(ang) = |v| a, and
   thus: d = sqrt(u*2 - a*2) 
Want closepoint to be with the p1 to p2 interval, 
   (i.e.  0<a< |v})
  Allow some leaway;  -0.5 < a < |v|+0.5
*/
{
   int k;
   float v[3],u[3],magu,magv,a;
   for(k=0;k<=2;k++){ 
      v[k] = pos2[k] -pos1[k];// vector from ipt to ipt+1
      u[k] = capos[k]-pos1[k];
   }
 magv = sqrt( dot(v,v) );
 if(magv<0.001){ //pos1=pos2:  then return dist from capos to pos1
      *idist=0; // Notin interval - since no inverval exist
 }
 else{ 
   magu = sqrt( dot(u,u) );
   a = dot(u,v)/magv;
   if( (a >0.0) && (a< magv ) ){ // in interval
       *idist = 1;
       *dist  = sqrt( magu*magu - a*a);
       for(k=0;k<=2;k++) 
          pos[k] = pos1[k] + a*v[k]/magv;
   }
   else{
      *idist = 0;
   }
 }// magv not = 0
}

#if 1 //  Used to cut bad Ca's
int testcaden(struct griddata& den1,struct geometry& geo)
/* routine to look for poor CA density 
  If 2 or more bad, it returns -1 - indicates fit should be deleted
  If 1 - returns the Ca value to be cut (if long enough) by checkextendedfit
  Returns 0 if OK
*/
{
   int i,j,countbad,ibad;
   float aden;

   j=1; // Ca
   countbad = 0;
   for(i=geo.nend[0]+1;i<=geo.nend[1]-1;i++){ // just look at atoms with sides
       aden= aveatom(den1,geo,j,i,geo.pos);
       if(aden<=geo.poorrigid ){
          countbad = countbad+1;
          ibad = i; // last bad residue
       }
   }// i
   if(countbad>=2)
     return -1;
   else if( countbad ==1)
     return ibad;
   else
     return 0;
}
#endif

#if 0 // for degugging prints out all densities
int testden(struct griddata& den1,struct geometry& geo)
/* routine to look for bad density and too much variation
*/
{
   int i,j;
   float aden;

   for(i=geo.nend[0]+1;i<=geo.nend[1]-1;i++){ // just look at atoms with sides
     printf("   For res %d:",i-geo.nend[0]);
     for(j=1;j<=5;j++){
       aden= aveatom(den1,geo,j,i,geo.pos);
       printf("  %5.3f",aden);
     }
     printf("\n");
   }// i
   return 0;  // dummy
}
#endif

void refine(struct geometry& geo,struct griddata& den1,
      int savlen,int *quittrace)
/* Runs tordynbat and rechecks  - returns quitrace
  New:  returns quittrace = 2 if < savlen -1
  If  = savlen-1 returns quittrace =10 and refines
*/
{
    int quitcond;
    int nbad;
    float lastden,origrfor;
//    float aveco,avecb2;

  if(dynquit==0)
     quitcond = geo.fixend;
  else 
     quitcond = dynquit;
   origrfor = tordata.rfor;
   if ( (geo.nend[1]-geo.nend[0])<savlen ){
//if(verbose >=1)
  fprintf(fpout," Bad fit after recheck because extended fit < savlen = %d\n",savlen);
         *quittrace = 2; // found bad fit
          goto endrout;
     } 
     else{// refine
#if 1 
        if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
            tordata.adata[1] = 2; 
        tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);//  run quitcond rounds of refinement at use last value
        if(tordata.irepeat>=1){ //run 30 more times with small random force
//         tordata.rfor = tordata.rfor2;
#if 1 // New  3 mar00:  run on all atoms for 30 times with rfor = 0
           geo.fixend = 0; // Run on all atoms
           quitcond = 30;
           tordata.rfor = 0;
#endif
           tordata.adata[1] = 2; 
           tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);//  run quitcond rounds of refinement at use last value
           tordata.rfor = origrfor;
       }// irepeat
#endif
#if 0
        tordata.adata[1] = 2;
        quitcond = 30;
        geo.fixend=0; 
        tordata.rfor = 0;
        tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);//  run refinement for entire chain
#endif
     }
     sidenumreg(geo.nend[0],geo.nend[1],den1,geo); //Assign new side lengths
     endrout:;
     tordata.rfor = origrfor;
}

void cutend(struct griddata& den1,struct geometry& geo,
    int *quittrace,int iremove)
{
  int icut,mmin;

      if(iremove==0)
         icut = 0;
      else if(  (iremove == 1)||(iremove==2) )
         icut =1;
      else // cut both ends
         icut = 2;
      if(iremove !=0){
         mmin = geo.savlen -(geo.nend[1]-geo.nend[0] -icut);
         if( mmin>0 )
         if(verbose>=1)
           printf("   Cutting end reduces length %d below savlen\n",mmin);
      }
      if(iremove ==1){
         if(verbose>=1)
              printf(" Remove nterm end because bad checkend\n");
          geo.nafit = geo.nafit- 1;
          geo.nend[0]= geo.nend[0]+1;  
          geo.fixend = 1; //refine N term end
          refine(geo,den1,geo.savlen,quittrace); // quitrace is pointer
          goto quitrout;
      }
      else if (iremove==2){
         if(verbose>=1)
             printf(" Remove cterm end because bad checkend\n");
          geo.nafit = geo.nafit- 1;
          geo.nend[1] = geo.nend[1] -1;
          geo.fixend = 2;
          refine(geo,den1,geo.savlen,quittrace); // quitrace is pointer
          goto quitrout;
     }
       else if (iremove==3){
         if(verbose>=1)
             printf("Remove nterm and cterm end because bad checkend\n");
          geo.nafit = geo.nafit- 2;
          geo.nend[0]= geo.nend[0]+1;  
          geo.nend[1] = geo.nend[1] -1;
          geo.fixend = 0; // refine all residues
          refine(geo,den1,geo.savlen,quittrace); // quitrace is pointer
          goto quitrout;
    }
    quitrout:;
}

void cutco(struct geometry& geo,struct griddata& den1,
          int resbad[],int savlen,int *quittrace,int i, const char angdist[20]) 
/*IDENTICAL TO cutside
 Used for cutting when bad co angle
 cuts ith residue, and then refines - angdist is either angle
*/
{
   int ncut;
   if( i< (geo.nend[1]+geo.nend[0])/2 ){// cut N term end
           ncut = i-geo.nend[0]+1;
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
if(verbose >=1)
  fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
   }
   else{ //cut C term end
           ncut = geo.nend[1]-1 - i;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           geo.fixend = 2;
if(verbose >=1)
  fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
    }
    quitrout:;
}


void cutside(struct geometry& geo,struct griddata& den1,
          int resbad[],int savlen,int *quittrace,int i, const char angdist[20]) 
/* 
IDENTICAL TO cutco
Used for cutting when bad helix or sheet sideang
 cuts ith residue, and then refines - angdist is either angle
*/
{
   int ncut;
   if( i< (geo.nend[1]+geo.nend[0])/2 ){// cut N term end
           ncut = i-geo.nend[0]+1;
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
if(verbose >=1)
  fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
   }
   else{ //cut C term end
           ncut = geo.nend[1] - i;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           geo.fixend = 2;
if(verbose >=1)
  fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
    }
    quitrout:;
}

void cutsidefit(struct geometry& geo,struct griddata& den1,
          int ifit,const char angdist[20]) 
/*  Inputs  fit:  where fit[ifit].itemp2[0][0][0] is the 
   number of bad residues in ifit and fit[ifit].itemp2[ires][0][0] is
   the residue number of the bad fit
   ires runs from 0 to fit.nres-1.
     for(ires=0;ires<fit[ifit].nres;ires++){ // each residue
   i0 = geo.nend[0];
        fit.pos[j][1][i] = geo.pos[i0+j][1][i];
   Returns      fit[ifit].itemp[0]=0; if  bad fit 
*/ 
{
   int ncut,i0,ires,ctermcut,ntermcut,cutgeo[10],nii[2];
   int quittrace;
   
   copyfittogeo(geo,fit[ifit]);
   i0 = geo.nend[0];
   ncut=fit[ifit].itemp2[0][0][0];
   for(ires=1;ires<= ncut;ires++)
      cutgeo[ires] = i0+ fit[ifit].itemp2[ires][0][0];
   ntermcut = -10;
   ctermcut = 10000;
   geo.fixend = -1;
   quittrace=4;  //Used as indicator that testfit was run in this routine Set to good value, then set to bad below
   for(ires=1;ires<= ncut;ires++){// find largest Nterm and C term cuts
      if( (cutgeo[ires]< (geo.nend[1]+geo.nend[0])/2) && 
               (cutgeo[ires]> ntermcut) )
          ntermcut = cutgeo[ires]; 
      else if( (cutgeo[ires]>= (geo.nend[1]+geo.nend[0])/2) &&
               ( cutgeo[ires]< ctermcut) )
          ctermcut = cutgeo[ires]; 
   }
   if(ntermcut > 0){ // cut N term end   
           ncut = ntermcut-geo.nend[0]+1;
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
         if(verbose >=1)
               fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
   }
   if(ctermcut <1000){ // cut C term end   
           ncut = geo.nend[1] - ctermcut;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           if(geo.fixend==1)// indicates cut both end
              geo.fixend = 0; //refine all
           else
              geo.fixend = 2;//refine C term end
          if(verbose >=1)
             fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
  }
  if( (ntermcut>0) || ctermcut<1000)
        refine(geo,den1,geo.savlen,&quittrace); // quitrace is pointer
  if(quittrace==2)
      fit[ifit].itemp[0]=0; //indicates  bad fit  - must be removed
  else{
      nii[0]= fit[ifit].val1;
      nii[1] = fit[ifit].val2;
      copygeotofit(geo,fit[ifit],nii);
  }
}

void cutmultside(struct geometry& geo,struct griddata& den1,
          int resbad[],int savlen,int *quittrace, const char angdist[20]) 
/* 
  may be more than 1 bad res -
  figures out where to cut
*/ 
{
   int j,ncut,ires,ctermcut,ntermcut,cutgeo[10];
   
   ncut = 0;
   for(j=geo.nend[0];j<=geo.nend[1];j++){ 
      if(resbad[j]>0){
          ncut = ncut+1;
          cutgeo[ncut] = j;
      }
   }
   
   ntermcut = -10;
   ctermcut = 10000;
   geo.fixend = -1;
   for(ires=1;ires<= ncut;ires++){// find largest Nterm and C term cuts
      if( (cutgeo[ires]< (geo.nend[1]+geo.nend[0])/2) && 
               (cutgeo[ires]> ntermcut) )
          ntermcut = cutgeo[ires]; 
      else if( (cutgeo[ires]>= (geo.nend[1]+geo.nend[0])/2) &&
               ( cutgeo[ires]< ctermcut) )
          ctermcut = cutgeo[ires]; 
   }
   if(ntermcut > 0){ // cut N term end   
           ncut = ntermcut-geo.nend[0]+1;
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
         if(verbose >=1)
               fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
   }
   if(ctermcut <1000){ // cut C term end   
           ncut = geo.nend[1]+1 - ctermcut;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           if(geo.fixend==1)// indicates cut both end
              geo.fixend = 0; //refine all
           else
              geo.fixend = 2;//refine C term end
          if(verbose >=1)
             fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
  }
  if( (ntermcut>0) || ctermcut<1000)
        refine(geo,den1,geo.savlen,quittrace); // quitrace is pointer
}


void cutca4(struct geometry& geo,struct griddata& den1,
          int resbad[],int savlen,int *quittrace,int i, const char angdist[20]) 
// Used for cutting when bad Ca[i] -Ca[i+4] dist
// cuts ith residue, and then refines - angdist is either angle
{
   int ncut;
   if( (i-2)< (geo.nend[1]+geo.nend[0])/2 ){// cut N term end
           ncut = i-geo.nend[0]-2;
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
if(verbose>=1)
  fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
   }
   else{ //cut C term end
           ncut = geo.nend[1]-i+1;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           geo.fixend = 2;
if(verbose >=1)
  fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
    }
    quitrout:;
}
void cutca3(struct geometry& geo,struct griddata& den1,
          int resbad[],int savlen,int *quittrace,int i, const char angdist[20]) 
// Used for cutting when bad 0[i] ca[i+2] dist
// cuts ith residue, and then refines - angdist is either angle
{
   int ncut;
   if( (i-1)< (geo.nend[1]+geo.nend[0])/2 ){// cut N term end
           ncut = i-geo.nend[0]-1;
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
if(verbose>=1)
  fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
   }
   else{ //cut C term end
           ncut = geo.nend[1]-i+1;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           geo.fixend = 2;
if(verbose >=1)
  fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
    }
    quitrout:;
}

void cutoo(struct geometry& geo,struct griddata& den1,
          int resbad[],int savlen,int *quittrace,int i, const char angdist[20]) 
// Used for cutting when bad oodist in sheet
// cuts ith residue, and then refines - 
{
   int ncut;
   if( i< (geo.nend[1]+geo.nend[0])/2 ){// cut N term end
           ncut = i-geo.nend[0];
           geo.nafit = geo.nafit- ncut;
           geo.nend[0]= geo.nend[0]+ncut;  
           geo.fixend = 1; //refine N term end
if(verbose >=1)
  fprintf(fpout,"Removing %d N term residue in checkfit because bad %s\n",
                   ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
   }
   else{ //cut C term end
           ncut = geo.nend[1] - i;
           geo.nafit = geo.nafit - ncut;
           geo.nend[1] = geo.nend[1] -ncut;
           geo.fixend = 2;
if(verbose >=1)
  fprintf(fpout,"Removing %d C term residue in checkfit because bad %s\n",
                 ncut,angdist);
           refine(geo,den1,savlen,quittrace); // quitrace is pointer
           goto quitrout;
    }
    quitrout:;
}

void checkextendedfit(int icheckhelsh,struct geometry& geo,
   struct griddata& den1,int resbad[],int *quittrace)
/* Recursive function.
NEW: 17feb00:  Add two new criteria for deleting helixes.
New 2dec99: check o-o dist for sheet
NEW: 1dec99:  reduces value of minco
New - check each end and then check for bad angles

if icheckhelsh = 1 then check sheet/helix angles
New: at end, check if residue second from end has nbad 
   indicated by geo.res0[i]>0,  
   Remove both residues at end if it will not make too short.
 
 Check if bad atom at ends (if so, remove) or in center (may need to quit)
 Procedure: 
 New - only check end residues
    1) Check if need to remove residue
    2) If remove, refine with tordyn and then recursively call checkextendedfit
  quittrace = 4 on input, if bad fit, then set quittrace =2
*/
{
    int nbad,initialnum;
    int i,j;
    float aveco,lastden,avecb2;
int savlen;
int checkca_ca3dist=1;  // check ca3dist
//int checkca_ca3dist=0;  
int checko_odist =1; // check oodist in sheet
    
  initialnum = geo.nafit;
  *quittrace=4;  //Used as indicator that testfit was run in this routine Set to good value, then set to bad below
   savlen = geo.savlen;

  if(icheckhelsh==1){ // New  If sheet, find angles between succesive Ca-Cb vectors and delete bad residues and between neighbor CO bonds
    float ang[NAFIT],coang[NAFIT],minco;
    float maxhelco,mincoend,minhelang,mincocen;
#if 0 // Used up to 30may00 
    minhelang = 200.0; // Minimum allowed side angle for helix
#endif
#if 1 // New 30may00 
    minhelang = 190.0; // Minimum allowed side angle for helix
#endif
    maxhelco = 74.0;  // Maximum co angle for helix
#if 0 // Used up to 1 dec99 (seems to remove some good residues)
    mincoend = 100.0;  // minimum allowed angle between neighbor CO's for end
    mincocen = 120.0; // min for internal angle 
#endif
#if 1 // Used after to 1 dec99 
    mincoend = 75.0;  //Used up to 5jun00 - minimum allowed angle between neighbor CO's for end
//    mincoend = 95.0;  //Tried after 5jun00 minimum allowed angle between neighbor CO's for end
    mincocen = 95.0; // min for internal angle 
#endif
    if(geo.itype==2){
       for(i=geo.nend[0];i<=geo.nend[1]-2;i++){ 
         coang[i]=corot(geo.pos[i][2],geo.pos[i+1][2],
                        geo.pos[i][3],geo.pos[i+1][3]);
         if( (i==geo.nend[0])||(i==geo.nend[1]-2) )// end angles for corot
           minco= mincoend;
         else
           minco= mincocen;
//printf(" coang[%d] = %5.3f  minco = %5.3f\n",i,coang[i],minco);
         if(coang[i]< minco ){ // bad angle
             if(verbose>=1)
                printf("Bad coang[%d] = %5.3f   minco = %5.3f\n",i,coang[i],minco);
             cutco(geo,den1,resbad,savlen,quittrace,i,"shcoang"); 
             if(*quittrace==2)
                goto endrout;
          } // bad coang
          else if(i>geo.nend[0]) { // check if both bad coang and siderot
              ang[i] =siderot(geo.pos[i][1],geo.pos[i+1][1],
                        geo.pos[i][5],geo.pos[i+1][5]);
//              if(ang[i]>240.0){
              if(ang[i]>250.0){  // New 1feb01 - increase to 250
                 if(verbose>=1)
                    printf("Residue %d: BAD sideang (%5.3f)\n",i,ang[i]);
                 cutside(geo,den1,resbad,savlen,quittrace,i,"sideangle"); 
                 if(*quittrace==2)
                    goto endrout;
                 else
                    goto testhelside; // must jump out of loop because geo.nend changed
              }
              if( (ang[i]>217.0)&&(coang[i]<115.0) ){
                 if(verbose>=1)
                    printf("Residue %d: BAD sideang (%5.3f) and coang (%5.3f)\n",i,ang[i],coang[i]);
                 cutside(geo,den1,resbad,savlen,quittrace,i,"side_and_coang"); 
                 if(*quittrace==2)
                    goto endrout;
                 else
                    goto testhelside; // must jump out of loop because geo.nend changed
              }
          } // test both angles
       }// i
    } // sheet
//    testsideang:; // NOTE - CHANGED 23MAY00 - DELETED testsideangle - see checkfit.232may00 for old version
    testhelside:;
    if(geo.itype==1){ // for helix check  siderot angle
      for(i=geo.nend[0]+1;i<=geo.nend[1]-2;i++){ 
         ang[i] =siderot(geo.pos[i][1],geo.pos[i+1][1],
                        geo.pos[i][5],geo.pos[i+1][5]);
         if(ang[i] < minhelang){ 
             if(verbose>=2)
                printf(" Bad helsideang[%d] = %5.3f < minhelang = %5.3f\n",i,ang[i],minhelang);
             cutside(geo,den1,resbad,savlen,quittrace,i,"helsideang"); 
             if(*quittrace==2)
                goto endrout;
             else
                goto testcoaxis;
         } // bad angle
      } // i
    } // helix
    testcoaxis:;
#if 0 // 2mar00 - this removes severa; good helix - do not use it
    float maxcoaxisang;
    maxcoaxisang = 65.0;  // Max allowed angle between axis and CO for helix
    if(geo.itype==1){ // for helix check  count number of corot>74 and then cut
     if(geo.nend[1]-geo.nend[0]<=4){
if(verbose>=2)
        printf(" geo is too short to check C0 - axis angle\n");
        goto testcohel;
     }
     for(i=geo.nend[0];i<=geo.nend[1]-1;i++){
        ang[i] = co_axisang(geo,i);
        if(ang[i] > maxcoaxisang){ 
             cutside(geo,den1,resbad,savlen,quittrace,i,"co_axisang"); 
             if(*quittrace==2)
                goto endrout;
             else
                goto testcohel;
        } // bad angle
      } // i
    } // helix
    testcohel:;
#endif
#if 1 // 31may00 - add criteria for sheet
    float maxcoaxisang;
    maxcoaxisang = 110.0;  // Max allowed angle between axis and CO for sheet
    if(geo.itype==2){ // for helix check  count number of corot>74 and then cut
     if(geo.nend[1]-geo.nend[0]<=4){
if(verbose>=2)
        printf(" geo is too short to check C0 - axis angle\n");
        goto testside_Odist;
     }
     for(i=geo.nend[0];i<=geo.nend[1]-1;i++){
        ang[i] = co_axisang(geo,i);
        if(ang[i] > maxcoaxisang){ 
             cutside(geo,den1,resbad,savlen,quittrace,i,"co_axisang"); 
             if(*quittrace==2)
                goto endrout;
             else
                goto testside_Odist;
        } // bad angle
      } // i
    } //sheet
#endif
testside_Odist:;
#if 1 // 31may00 - add NEW criteria for sheet CHECK O atoms to Side atom distnace
   float mag;
   if(geo.itype==2){ // for sheet
      for(i=geo.nend[0]+1;i<=geo.nend[1]-2;i++){ 
        if(geo.iside[i]>=3){
           for(j=3;j<=geo.iside[i];j++){
              mag = distv(geo.pos[i][4+j],geo.pos[i][3]);
              if(mag<=2.0){
                  cutside(geo,den1,resbad,savlen,quittrace,i,"O_sheetsidedist"); 
                  if(*quittrace==2)
                     goto endrout;
                  else
                    goto testcohel;
              } // <2.0
           }// j
        } // >=3
      }//i
   } // itype
#endif
    testcohel:;
    if(geo.itype==1){ // for helix check  count number of corot>74 and then cut
      int count,totcount;
      count = totcount = 0;
      for(i=geo.nend[0];i<=geo.nend[1]-2;i++){ 
         totcount = totcount+1;
         coang[i]=corot(geo.pos[i][2],geo.pos[i+1][2],
                        geo.pos[i][3],geo.pos[i+1][3]);
         if(coang[i]>maxhelco)
            count = count+1;
      }
      if(count > totcount/2.0){
         for(i=geo.nend[0];i<=geo.nend[1]-2;i++){ 
             coang[i]=corot(geo.pos[i][2],geo.pos[i+1][2],
                        geo.pos[i][3],geo.pos[i+1][3]);
             if(coang[i]> maxhelco){ 
                cutco(geo,den1,resbad,savlen,quittrace,i,"helcoang"); 
                if(*quittrace==2)
                   goto endrout;
                else
                   goto testcadist;
             } // bad angle
         } // i
       } // large count
    } // helix
} //icheckhelsh=1
   testcadist:;
if(checkca_ca3dist==1){ // NEW- 11oct00 - add another helix check (helixca_caang) -New 28apr99 - check ca[i] to ca[i+3] and o[i] to ca[i+2] dist in helixes
   float maxca4dist,ca4dist,minang,ang;
   float mincocadist,cocadist;
   maxca4dist = 6.2; // maximum allowed ca[i] to ca[i+3] dist
//   maxca4dist = 7.0; //NEW 22jun00 - increase maximum allowed ca[i] to ca[i+3] dist
   mincocadist = 3.2; // minimum allowed Distance from o[i] to ca[i+2]
   minang = 1.70;
   if(geo.itype==1){ // only check helix
      for(i=geo.nend[0]+1;i<=geo.nend[1]-2;i++){ 
         ang = helixca_caang(geo.pos[i-1][1],geo.pos[i][1],geo.pos[i+2][1]);
         if(ang<minang){
            if(verbose>=1)
               printf("Bad helixca_caang:res %d helixca_caang =%5.3f < %5.3f\n",i,ang,minang); 
            cutco(geo,den1,resbad,savlen,quittrace,i,"helixca_caang"); // cuts at residue i, identical to cutside
            if(*quittrace==2)
               goto endrout;
            else
                goto checkca4cadist;
       } // bad ca4dist
     } // i
     checkca4cadist:;
     for(i=geo.nend[0]+3;i<=geo.nend[1];i++){ 
       ca4dist = distv(geo.pos[i][1],geo.pos[i-3][1]);
       if(ca4dist>maxca4dist){ // then bad
           if(verbose>=2)
               printf("Bad ca4dist[i=%d] = %5.3f > %5.3f\n", i,ca4dist,maxca4dist);
           cutca4(geo,den1,resbad,savlen,quittrace,i,"ca4dist"); 
           if(*quittrace==2)
               goto endrout;
           else
                goto checkcocadist;
       } // bad ca4dist
     } // i
     checkcocadist:;
     for(i=geo.nend[0]+2;i<=geo.nend[1];i++){ // Distance from o[i] to ca[i+2]
        cocadist = distv(geo.pos[i][1],geo.pos[i-2][3]);
        if(cocadist < mincocadist){ // then bad
           cutca3(geo,den1,resbad,savlen,quittrace,i,"cocadist"); 
           if(*quittrace==2)
               goto endrout;
           else
                goto endoca;
       } // bad cocadist
     }
    } // helix
} //icheckca_ca3dsist=1
endoca:;
if(checko_odist ==1){ // New 1dec99 -- check o-o dist for sheet
//    float minoodist = 3.2; // Used upto 30may00
//    float minoodist = 3.0; // tried after 30may00
    float minoodist = 2.85; // tried after 24aug00 - needed for fah protein?
    float oodist;
    if(geo.itype==2){ // only check sheet
      for(i=geo.nend[0]+1;i<geo.nend[1];i++){ 
        oodist = distv(geo.pos[i][3],geo.pos[i-1][3]);
       if(oodist<minoodist){ // then bad
printf(" oodist[i=%d] = %5.3f  minoodist = %5.3f\n", i,oodist,minoodist);
           cutoo(geo,den1,resbad,savlen,quittrace,i,"oodist"); 
           if(*quittrace==2)
               goto endrout;
           else
                goto endoo;
       } // bad oodist
     } // i
    } // itype
} //  checkoodist
endoo:;
 
#if 0 // do not use these last criteria
// Add 3 new criteria for cutting - These routine Simply remove the bad fit - they do not cut out bad  
   aveden(den1,geo,geo.pos,&lastden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
#if 0 // Now called in trace.c++  (testfit() ) after removing end
   if( (geo.itype==2)&&((geo.nend[1]-geo.nend[0]) <= geo.savlen+1) ){ // if short sheet
      if(geo.avecb/avecb2 < 0.8){ // delete fit if cb/cb2 den ratio is small
        *quittrace=2;
        goto endrout;
      }
   }
#endif      
   cutmultside(geo,den1,resbad,savlen,quittrace,"density"); 
   if(*quittrace==2)
          goto endrout;
int testca;
   testca = testcaden(den1,geo);
   if(testca == -1){// routine to look for bad CA density If 2 or more bad, it returns -1
      printf(" DELETING  FIT BECAUSE  2 poor Ca den\n");
      *quittrace=2;
      goto endrout;
   }
   else  if(testca >0){ // incidcates one bad Ca, make cut at i = testcaden
       cutside(geo,den1,resbad,savlen,quittrace,testca,"badcaden"); 
       if(*quittrace==2)
             goto endrout;
   }
#if 0  // Now moved to testfit()
   if(geo.itype==1){// Add 2 new criteria for cutting helix - These routine Simply remove the bad fit - they do not cut out bad
        if(cadel(geo)==1){// NEW (17feb00) cut helix if bad Cai-1 to Cai+1 dist
            printf(" DELETING HELIX FIT BECAUSE BAD CA(i+1)--CA(i-1) dist\n");
            *quittrace=2;
            goto endrout;
        }
        if(ca_caang3(geo)==1){// Find angle between successive Ca-Ca vectors
            printf(" DELETING HELIX FIT BECAUSE BAD angle between successive Ca-Ca vectors\n");
            *quittrace=2;
            goto endrout;
        }
   }
#endif
   if( (resbad[geo.nend[0]]>0) || (resbad[geo.nend[0]+1]>0)  ){ // remove N term end
       geo.nafit=geo.nafit-1; 
       geo.nend[0] = geo.nend[0]+1;
       geo.fixend = 1; //refine N term end
//if(verbose >=1)
  fprintf(fpout,"Removing N term residue in checkextended fit because bad resbad\n");
       refine(geo,den1,savlen,quittrace); // quitrace is pointer
       *quittrace=2;
       goto endrout;
    } // removed N term end
    if( (resbad[geo.nend[1]-1]>0)||(resbad[geo.nend[1]-2]>0) ){ // remove C term end
       geo.nafit=geo.nafit-1; 
       geo.nend[1] = geo.nend[1]-1;
       geo.fixend = 2;// refine C term end
//if(verbose >=1)
  fprintf(fpout,"Removing C term residue in checkextended fit because bad resbad\n");
       refine(geo,den1,savlen,quittrace); // quitrace is pointer
       *quittrace=2;
       goto endrout;
   } // removed C term end
#endif // last 3 criteria

   if(initialnum!=geo.nafit){
      if ( (geo.nend[1]-geo.nend[0])<savlen ){
if(verbose >=1)
  fprintf(fpout," Bad fit in checkext because extended fit = %d not long enough\n",geo.nend[1]-geo.nend[0]);
           *quittrace = 2; // found bad fit
           goto endrout;           
      }
      aveden(den1,geo,geo.pos,&lastden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
      checkextendedfit(icheckhelsh,geo,den1,resbad,quittrace); //Recursive call
      if(*quittrace==2)
              goto endrout;
   }
   endrout:;
}

#if 0 // old version
void bonedistgeo(struct griddata& den1,struct geometry& geo,
         struct bonedata& bone,int *quittrace)
/* same as bonedist - but only applied to one geo and uses geo in place of fit
returns *quittrace = 2 if bad and should be deleted
*/
{
    int i,j,ires,k,k2,ivec[3],ibone,totnum,cbadres;
    int i0,n0,nres;
    int npt1,npt2,ncon,icon,icon2,ncon2,ccount,tcount,nbad,nbadres,nbadres2,atomcount;
    float fvec[3],dist,dvec[3];
//    float mdist1,mdist2;
    float perpdist,perpdist2,linedist,linedist2,pos1[3],pos2[3];
    float aveperp,aveline,fractbad,baddist;
    int ibadres[10];
    float mag1,mag2,vec1[3],vec2[3],ang;
    int iden;
    char origbonename[300];
    char ch[10];
    float ftemp[NAFIT][5][2];
    int  itemp2[NAFIT][5][2],resbad[NAFIT];
    FILE *fp;

    i0 = geo.nend[0]; //geo[i0] corresponds to fit[0]
    n0 = geo.nend[1]-geo.nend[0];
    nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0)       
    baddist = 1.0;  // perp. bad dist
    printf("Checking agreement of geo to bones 1.0 and deleting bad fits\n");
    strcpy(origbonename,bone.bonename);
    iden = 1;  // this corresponds to bone 1.0
    bone.minden = fminden[iden]; 
    strcpy(bone.bonename,den1.mapname);
    strcat(bone.bonename,".bone"); // bone  name is mapname with .bone_minden appended
    sprintf(ch,"%1.1f",bone.minden);
    strcat(bone.bonename,ch);
    fp=fopen(bone.bonename,"r");
    if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone);
    }
    else{
       fclose(fp);
       readbonefile(bone); // This reopens bone file
    }

    for(ires=i0;ires<nres;ires++) // each residue 
       for(i=1;i<=4;i++){
           ftemp[ires][i][0]=ftemp[ires][i][1] = 10000.0;
           for(k=0;k<=1;k++)
              itemp2[ires][i][k]= 1;
       }
    totnum = bone.ptv.end()-bone.ptv.begin();
    for(ibone=0;ibone<totnum;++ibone){
       d1tod3(bone.cmax,bone.ptv[ibone].pos,ivec);
       imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
        if( (nres-1)>=(geo.savlen-1) ){ // if long enough 
         for(ires=i0;ires<nres;ires++) // each residue 
          for(i=1;i<=4;i++){
            if( (i==4)&&(ires==i0) ) // N not defined for first residue
               goto nextatom; 
            if( ((i==3)||(i==2))&&(ires== nres-1) ) // O and C not defined for last atom
               goto nextatom; 
            for(k=0;k<=2;k++)
               dvec[k]= fvec[k]-geo.pos[ires][i][k];
            dist = dot(dvec,dvec); // squared dist from bone
            if(dist< ftemp[ires][i][0]){
               ftemp[ires][i][1]=ftemp[ires][i][0];
               ftemp[ires][i][0] = dist;
               itemp2[ires][i][1]=itemp2[ires][i][0];
               itemp2[ires][i][0]=ibone;
            }
            else if(dist< ftemp[ires][i][1]){
               ftemp[ires][i][1] = dist;
               itemp2[ires][i][1]=ibone;
            }   
            nextatom:;
         } // i
       } // long enoug
    } // bone
    if( (nres-1)>=(geo.savlen-1) ){ // if long enough 
     tcount = ccount= nbad = cbadres = 0; 
     aveline = aveperp =  0.0;
     for(ires=i0;ires<nres;ires++){ // each residue
       nbadres=nbadres2= atomcount = 0;
       if(verbose>=2) 
          printf("res %d:",ires);
       for(i=1;i<=4;i++){
         if(i!=3){ // do not include O
           ang = -1000.0;
           if( (i==4)&&(ires==i0) ) // N not defined for first residue
               goto nextatom2; 
           if( ((i==3)||(i==2))&&(ires== nres-1) ) // O and C not defined for last atom
               goto nextatom2;
           atomcount= atomcount+1;
           tcount = tcount+1;
           npt1= itemp2[ires][i][0];
           npt2= itemp2[ires][i][1];
           ncon = bone.ptv[npt1].num; // number of connections
           for(k=0;k<ncon;k++){ // check if close pts are connected
              icon = bone.ptv[npt1].con[k]; // the pt num of the connect point
              if(icon == npt2){ // close pts are connected
                itemp2[ires][i][1]=1; //indicates connected
                ccount = ccount+1; // count connected
                d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
                imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
                imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                closedist2(pos1,pos2,geo.pos[ires][i],&perpdist,
                                &linedist);
                goto foundconnect;
              }
           } // k
           for(k=0;k<ncon;k++){ // check if close pts are connected by one pt and find closest perpdist
              icon = bone.ptv[npt1].con[k]; // the pt num of the connect point
              ncon2 = bone.ptv[icon].num; // number of connections
              for(k2=0;k2<ncon2;k2++){                
                icon2 = bone.ptv[icon].con[k2]; // the pt num of the connect point
                if(icon2 == npt2){ // close pts are connected through another pt = icon2
                   ccount = ccount+1; // count connected
                   itemp2[ires][i][1]=2; //indicates connected through another pt                   
                   d1tod3(bone.cmax,bone.ptv[icon2].pos,ivec);
                   imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                   d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,geo.pos[ires][i],&perpdist,
                                &linedist);
                   d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,geo.pos[ires][i],&perpdist2,
                                &linedist2);
                   if(perpdist2<perpdist){
                       perpdist = perpdist2;
                       linedist2 = linedist2;
                   }
                   goto foundconnect;
                 } // if connected
               }// k2
            }// k
           itemp2[ires][i][1]=0; //if reach here not connected
           nbad = nbad+1;  // not connected points are considered bad no matter what perpdist
           nbadres2 = nbadres2+1; // counts bad perp and not connected aotms
           perpdist = ftemp[ires][i][0];// set perpdist = dist to closedt point
           if(verbose>=2) 
               printf(" noc %5.3f",ftemp[ires][i][0]);
           goto nextatom2;
           foundconnect:; // skip to here if want to interpolate between close pts
           aveline = aveline +fabs(linedist);
           aveperp = aveperp + perpdist;
           if(perpdist>baddist){
              nbad = nbad+1;
              d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);// find angle between atom and vec. to 2 nearest pts
              imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
              d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
              imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
              for(j=0;j<=2;j++){ 
                 vec1[j] = geo.pos[ires][i][j]-pos1[j];
                 vec2[j] = pos2[j] - -pos1[j];
              }
              mag1 = sqrt( dot(vec1,vec1) );
              mag2 = sqrt( dot(vec2,vec2) );
              ang = (180.0/PI)*acos( dot(vec1,vec2)/(mag1*mag2) );
           }
           if(verbose>=2){
              if(ang<0.0){// default - good perpdist
                 printf(" %5.3f %5.3f  ", perpdist,linedist); 
                 if(fabs(linedist)>1.0)
                    nbadres2 = nbadres2+1; // counts connected, but bad points
              }
              else if(ang <= 110.0){ // then bad dist and points behave as if connected              
                 nbadres = nbadres+1; // counts connected, but bad points
                 nbadres2 = nbadres2+1; // counts connected, but bad points
                 printf("ang: %5.3f %5.3f %5.3f  ",ang, perpdist,linedist);
              }
              else{ // bad perp dist but really not connected
                 nbadres2 = nbadres2+1; // counts connected, but bad points
                 printf(" noc  ang: %5.3f %5.3f",ang,ftemp[ires][i][0]);
              }
           }
           nextatom2:;
         } // not O
       }// i 
       if(verbose>=2)
          printf("\n");
       if( (ires<=1)||(ires>=nres-2) ){
         if( (float)nbadres2/atomcount>0.6 ){//Cut end residue if bad perp or not connected.
           cbadres = cbadres+1;
           ibadres[cbadres] = ires;
           if(verbose>=2) 
              printf(" 60 percent atoms are bad for end res %d\n ",ires);
         }
       }   
       else if( (float)nbadres/atomcount>0.6){//Cut middle residue if bad perp
           cbadres = cbadres+1;
           ibadres[cbadres] = ires;
           if(verbose>=2) 
              printf(" 60 percent atoms are bad for res %d ",ires);
           if(cbadres>5){
              *quittrace= 2; // indicates bad fit 
              printf("Bad fit because more than 5 bad res in bonedist\n");
              goto quitrout;
           } 
       }
      } // ires
      aveperp = aveperp/ccount;
      aveline = aveline/ccount;
      fractbad = (float)nbad/tcount;
      if(verbose>=1){ 
          printf("fract bad = %5.3f fract con = %5.3f  aveperp = %5.3f aveline = %5.3f\n",
               fractbad,(float)ccount/tcount,aveperp,aveline);
          printf("\n");
      }
      if(cbadres> 0){
         itemp2[0][0][0]= cbadres;
         for(i=1;i<=cbadres;i++)
            itemp2[i][0][0]= ibadres[i];
         if(verbose>=1){
            printf("%d bad Residues:",itemp2[0][0][0]);
            for(i=1;i<=itemp2[0][0][0];i++)
                printf("   %d",itemp2[i][0][0]); 
            printf("\n");
         }
//         cutsidefit(geo,den1,ifit," bones"); 
         for(ires=i0;ires<nres;ires++) // Use cutmultside in place of cutsidefit
            resbad[ires]=0;
         for(i=1;i<=itemp2[0][0][0];i++)
            resbad[ itemp2[i][0][0] ]=1; 
         cutmultside(geo,den1,resbad,geo.savlen,quittrace," bad bone res"); 
      }
      if(aveperp >0.7){
         *quittrace= 2; // indicates bad fit 
         printf("Bad geo because aveperp = %5.3f\n",aveperp);
      }
      else if(fractbad > 0.25){
         *quittrace= 2; // indicates bad fit 
         printf("Bad geo because fraction bad = %5.3f\n",fractbad);
      }
     } // long enough
   quitrout:;
   strcpy(bone.bonename,origbonename); // put original bonemane in bone1
   readbonefile(bone);
   bone.iexist = -2; // turn off drawing of bone
}
#endif

void connect3(struct bonedata& bone,int npt1,int npt2,int ipath[])
/* - checks if point pt1 and pt2 are connected, 
    Goes 3 deep
  returns ipath[0]= -1 if not connected
  return ipath[0]= number of points on path
     ipath[i]= point numbers
     ipath[1]= npt1, ......  ipath[ ipath[0] ] = npt2
*/
{
   int k,k1,k2,k3,k4;
   int ncon,ncon2,ncon3,ncon4,ncon5;
   int icon1,icon2,icon3,icon4,icon5;

   ipath[0] = -1; // default not connected
   ncon = bone.ptv[npt1].num; // number of connections
#if 0 // NO even at at end of chain - can be connected
   if(ncon==1){ // then assume not connected
               goto quitrout;
   }
#endif
   for(k=0;k<ncon;k++){ // check if close pts are connected
       icon1 = bone.ptv[npt1].con[k]; // the pt num of the connect point
       if(icon1 == npt2){ // close pts are connected
           ipath[0]=2;
           ipath[1]=npt1;
           ipath[2]=npt2;
           goto quitrout;
       }
       ncon2 = bone.ptv[icon1].num; // number of connections
       for(k1=0;k1<ncon2;k1++){ // check if close pts are connected
          icon2 = bone.ptv[icon1].con[k1]; // the pt num of the connect point
          if(icon2 == npt2){ // close pts are connected
               ipath[0]=3;
               ipath[1]=npt1;
               ipath[2]=icon1;
               ipath[3]=npt2;
               goto quitrout;
          }
          ncon3 = bone.ptv[icon2].num; // number of connections
          for(k2=0;k2<ncon3;k2++){ // check if close pts are connected
             icon3 = bone.ptv[icon2].con[k2]; // the pt num of the connect point
             if(icon3 == npt2){ // close pts are connected
                 ipath[0]=4;
                 ipath[1]=npt1;
                 ipath[2]=icon1;
                 ipath[3]=icon2;
                 ipath[4]=npt2;
                 goto quitrout;
             }
             ncon4 = bone.ptv[icon3].num; // number of connections
             for(k3=0;k3<ncon4;k3++){ // check if close pts are connected
                icon4 = bone.ptv[icon3].con[k3]; // the pt num of the connect point
                if(icon4 == npt2){ // close pts are connected
                   ipath[0]=5;
                   ipath[1]=npt1;
                   ipath[2]=icon1;
                   ipath[3]=icon2;
                   ipath[4]=icon3;
                   ipath[5]=npt2;
                   goto quitrout;
                }

                ncon5 = bone.ptv[icon4].num; // number of connections
                for(k4=0;k4<ncon5;k4++){ // check if close pts are connected
                   icon5 = bone.ptv[icon4].con[k4]; // the pt num of the connect point
                   if(icon5 == npt2){ // close pts are connected
                     ipath[0]=6;
                     ipath[1]=npt1;
                     ipath[2]=icon1;
                     ipath[3]=icon2;
                     ipath[4]=icon3;
                     ipath[5]=icon4;
                     ipath[6]=npt2;
                     goto quitrout;
                  }
                }//k4
             }//k3


          }//k2
       }//k1
    }//k
    quitrout:;
}


#if 1 // new version
void bonedistgeo(struct griddata& den1,struct geometry& geo,
         struct bonedata& bone,int *quittrace)
/* same as bonedist - but only applied to one geo and uses geo in place of fit
returns *quittrace = 2 if bad and should be deleted
nbad - counts total number of bad atoms - either bad perp or not connected.
nbadres - counts atoms per res with bad perp
nbadres2  counts atoms per res that are not connected
cbadres counts residues with 2 bad atoms
*/
{
    int i,ii,kk,ires,k,ivec[3],ibone,totnum,cbadres;
    int i0,n0,nres;
    int npt1,npt2,ccount,tcount,nbad,nbadres,nbadres2,atomcount;
    float fvec[3],dist,dvec[3];
//    float mdist1,mdist2;
    float perpdist,linedist,pos1[3],pos2[3];
    float aveperp,fractbad,baddist;
    int ibadres[10];
    int iden;
    char origbonename[300];
    char ch[10];
    float ftemp[NAFIT][5][2];
    int  itemp2[NAFIT][5][2],resbad[NAFIT];
    int oldnpt1,oldnpt2,iatom[5];
    float minndist,npos1[3],npos2[3],opos1[3],opos2[3];
    int ipath[10];
    float minperp,minline;
    float aveline;
    FILE *fp;

    i0 = geo.nend[0]; //geo[i0] corresponds to fit[0]
    n0 = geo.nend[1]-geo.nend[0];
    nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0)       
    baddist = 1.0;  // perp. bad dist
    printf("Checking agreement of geo to bones 1.0 and deleting bad fits\n");
    strcpy(origbonename,bone.bonename);
    tcount = ccount= nbad = cbadres = 0; 
    aveline = aveperp =  0.0;
#if 0// start with higher bone - NOTE GOOD IDEA FOR TEST CONNECT, BUT NOT FOR PERP DIST WHERE CALL cultmultside
 int idenmin,idenmax;
 idenmax= 3;
 idenmin= 1;
 for(iden = idenmax;iden>=idenmin;iden--){ // Go from idenmax to idenmin, quit if OK
    *quittrace=0;  // default - good fit
#endif
 iden = 1; // corresponds to bone 1.0
    tcount = ccount= nbad = cbadres = 0; 
    aveline = aveperp =  0.0;
    bone.minden = fminden[iden]; 
    strcpy(bone.bonename,den1.mapname);
    strcat(bone.bonename,".bone"); // bone  name is mapname with .bone_minden appended
    sprintf(ch,"%1.1f",bone.minden);
    strcat(bone.bonename,ch);
    fp=fopen(bone.bonename,"r");
    iatom[1]=4; // iatom[1,2,3] = N,Ca,C
    iatom[2]=1;
    iatom[3]=2;
    if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone);
    }
    else{
       fclose(fp);
       readbonefile(bone); // This reopens bone file
    }

    for(ires=i0;ires<nres;ires++) // each residue 
       for(i=1;i<=4;i++){
           ftemp[ires][i][0]=ftemp[ires][i][1] = 10000.0;
           for(k=0;k<=1;k++)
              itemp2[ires][i][k]= 1;
       }
    totnum = bone.ptv.end()-bone.ptv.begin();
    for(ibone=0;ibone<totnum;++ibone){
       d1tod3(bone.cmax,bone.ptv[ibone].pos,ivec);
       imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
        if( (nres-1)>=(geo.savlen-1) ){ // if long enough 
         for(ires=i0;ires<nres;ires++) // each residue 
          for(ii=1;ii<=3;ii++){
            i=iatom[ii];
            if( (i==4)&&(ires==i0) ) // N not defined for first residue
               goto nextatom; 
            if( ((i==3)||(i==2))&&(ires== nres-1) ) // O and C not defined for last atom
               goto nextatom; 
            for(k=0;k<=2;k++)
               dvec[k]= fvec[k]-geo.pos[ires][i][k];
            dist = dot(dvec,dvec); // squared dist from bone
            if(dist< ftemp[ires][i][0]){
               ftemp[ires][i][1]=ftemp[ires][i][0];
               ftemp[ires][i][0] = dist;
               itemp2[ires][i][1]=itemp2[ires][i][0];
               itemp2[ires][i][0]=ibone;
            }
            else if(dist< ftemp[ires][i][1]){
               ftemp[ires][i][1] = dist;
               itemp2[ires][i][1]=ibone;
            }   
            nextatom:;
         } // i
       } // long enoug
    } // bone
    oldnpt1 = oldnpt2 = npt1 = npt2 = -10; // indicates not assigned
    if( (nres-1)>=(geo.savlen-1) ){ // if long enough 
     for(ires=i0;ires<nres;ires++){ // each residue
       nbadres=nbadres2= atomcount = 0;
       if(verbose>=2) 
          printf("res %d:",ires);
       for(ii=1;ii<=3;ii++){
           i=iatom[ii];
           if( (i==4)&&(ires==i0) ) // N not defined for first residue
               goto nextatom2; 
           if( ((i==3)||(i==2))&&(ires== nres-1) ) // O and C not defined for last atom
               goto nextatom2;
           atomcount= atomcount+1;
           tcount = tcount+1; // count total number atoms
           oldnpt1 = npt1;
           oldnpt2 = npt2;
           npt1= itemp2[ires][i][0];
           npt2= itemp2[ires][i][1];
           if(oldnpt1>=0){ // then find min. dist between old and new
              minndist = 1000.0;
              d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
              imcv(den1,ivec,npos1);// convert int grid ivec to float world co-ord fvec
              d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
              imcv(den1,ivec,npos2);// convert int grid ivec to float world co-ord fvec
              d1tod3(bone.cmax,bone.ptv[oldnpt1].pos,ivec);
              imcv(den1,ivec,opos1);// convert int grid ivec to float world co-ord fvec
              d1tod3(bone.cmax,bone.ptv[oldnpt2].pos,ivec);
              imcv(den1,ivec,opos2);// convert int grid ivec to float world co-ord fvec
              dist = distv(npos1,opos1);
              if(dist<minndist)
                 minndist = dist;
              dist = distv(npos1,opos2);
              if(dist<minndist)
                 minndist = dist;
              dist = distv(npos2,opos1);
              if(dist<minndist)
                 minndist = dist;
              dist = distv(npos2,opos2);
              if(dist<minndist)
                 minndist = dist;
              if (minndist>3.0){ // NEW 26oct00 criteria
                  if(verbose>=1)
                      printf("BAD mindist = %5.3f ires = %d iatom = %d  npt1 = %d npt2 = %d\n",minndist,ires,i,npt1,npt2);
                  cbadres = cbadres+1;
                  ibadres[cbadres] = ires;
              }
           }
#if 0 // Old version
           ncon = bone.ptv[npt1].num; // number of connections
           if(ncon==1){ // then assume not connected
               goto notconnect;
           }
           for(k=0;k<ncon;k++){ // check if close pts are connected
              icon = bone.ptv[npt1].con[k]; // the pt num of the connect point
              if(icon == npt2){ // close pts are connected
                d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
                imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
                imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                closedist2(pos1,pos2,geo.pos[ires][i],&perpdist,
                                &linedist);
                goto foundconnect;
              }
           } // k
           for(k=0;k<ncon;k++){ // check if close pts are connected by one pt and find closest perpdist
              icon = bone.ptv[npt1].con[k]; // the pt num of the connect point
              ncon2 = bone.ptv[icon].num; // number of connections
              for(k2=0;k2<ncon2;k2++){                
                icon2 = bone.ptv[icon].con[k2]; // the pt num of the connect point
                if(icon2 == npt2){ // close pts are connected through another pt = icon2
                   d1tod3(bone.cmax,bone.ptv[icon2].pos,ivec);
                   imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                   d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,geo.pos[ires][i],&perpdist,
                                &linedist);
                   d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,geo.pos[ires][i],&perpdist2,
                                &linedist2);
                   if(perpdist2<perpdist){
                       perpdist = perpdist2;
                       linedist = linedist2;
                   }
                   goto foundconnect;
                 } // if connected
               }// k2
            }// k
           notconnect:;
           nbad = nbad+1;  // not connected points are considered bad no matter what perpdist
           nbadres2 = nbadres2+1; // counts  not connected aotms/per residue
           if(verbose>=2) 
               printf(" noc %5.3f",ftemp[ires][i][0]);
           goto nextatom2;
           foundconnect:; // skip to here if want to interpolate between close pts
           if(verbose>=1)
               printf(" for res %d atom %d perpdist = %5.3f  linedist = %5.3f\n",ires,i,perpdist,linedist);
#endif

#if 1 // New  15feb01 - uses connect3 - checks if points connect by total of 4 point and returns path
           connect3(bone,npt1,npt2,ipath); // if ipath[0]= -1 not connected else connected by path - ipath[1]= npt1, ......  ipath[ ipath[0] ] = npt2
           if(ipath[0]== -1){ // not connected - then point is bad independent of perpdist
//              nbad = nbad+1;  // not connected points are considered bad no matter what perpdist
              nbad = nbad+2;  //19feb01 -  not connected points = 2 bad atoms
              nbadres2 = nbadres2+1; // counts  not connected aotms/per residue
              if(verbose>=1) 
                 printf(" res %d atom %d is not connected\n",ires,i);
              goto nextatom2;
           }
           else{ // connected on path ipath[1]= npt1, ......  ipath[ ipath[0] ] = npt2
               minperp = 10000.0;
               for(kk=2;kk<=ipath[0];kk++){
                   d1tod3(bone.cmax,bone.ptv[ ipath[kk-1] ].pos,ivec);
                   imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                   d1tod3(bone.cmax,bone.ptv[ ipath[kk] ].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,geo.pos[ires][i],&perpdist,
                                &linedist);
                   if(perpdist<minperp){
                       minperp = perpdist;
                       minline = linedist;
                   }
                }//kk
            }// connected
           perpdist = minperp;
           linedist = minline; // NOTE - linedist is not used - not meaningful - can have negative value
           if(verbose>=2)
               printf(" for res %d atom %d perpdist = %5.3f  linedist = %5.3f\n",ires,i,perpdist,linedist);
#endif
           ccount = ccount+1; // counts total connected
           aveperp = aveperp + perpdist;
           aveline = aveline +fabs(linedist);
#if 0 // NOT USED - NO GOOD  - look at neighbor perpdists
           perpd[0]=perpd[1]; //keeps track of ordered perdist
           perpd[1]=perpd[2];
           perpd[2]=perpdist;
           if( perpd[1]>baddist){
printf(" Bad perpd[1] = %5.3f, perpd[0 and 2] = %5.3f,%5.3f\n",perpd[1],perpd[0],perpd[2]);
           }
#endif
           if(perpdist>baddist){
//printf("Bad perpdist = %5.3f\n",perpdist);
              nbad = nbad+1;
              nbadres = nbadres+1; // counts aotms with bad perp
           }
#if 0 /* NOT USED - DID NOT SEEM TO BE ANY GOOD FOR DISTINGUISHING BAD SHEETSNew - 23oct00 - find angle between C-N vector and bone vector
    At this point - pos1 and pos2 are two closest bone pts to C if i=2.
*/
float cnvec[3],bonev[3],dist1,dist2,ang;
           for(k=0;k<=2;k++)// vector C to N
               cnvec[k] = geo.pos[ires][4][k]-geo.pos[ires][2][k];
           if(i==2){ // 
              dist1= distv(geo.pos[ires][4],pos1); // dist N to pos1
              dist2= distv(geo.pos[ires][4],pos2); // dist N to pos1
              if(dist1<dist2) // pos1 closer to N
                 for(k=0;k<=2;k++)
                    bonev[k] = pos1[k]-pos2[k];
              else
                 for(k=0;k<=2;k++)
                    bonev[k] = pos2[k]-pos1[k];
              dist1 = sqrt( dot(cnvec,cnvec) );
              dist2 = sqrt( dot(bonev,bonev) );
              ang = acos( dot(cnvec,bonev)/(dist1*dist2) );
if(ang*180.0/PI>50.0)
  printf(" res %d: %5.3f\n",ires,ang*180.0/PI);
           }
#endif
           nextatom2:;
       }// ii 
       if(verbose>=2) 
//   if(verbose>=1) 
         printf(" fraction bad perp = %5.3f  fraction not con = %5.3f",
              (float)nbadres/atomcount,(float)nbadres2/atomcount);
       if(verbose>=2)
//   if(verbose>=1)
          printf("\n");
     if( nbadres+nbadres2>=2){
           if(verbose>=1)
              printf(" Res # %d has %d bad atoms - should cut\n",ires-i0,nbadres+nbadres2);
           cbadres = cbadres+1;
           ibadres[cbadres] = ires;
       }
      } // ires
      aveperp = aveperp/ccount;
      aveline = aveline/ccount;
      fractbad = (float)nbad/tcount;
      geo.frcon = (float)ccount/tcount;
      if(verbose>=1){ 
          printf("fract bad = %5.3f fract con = %5.3f  aveperp = %5.3f aveline=%5.3f\n",
               fractbad,geo.frcon,aveperp,aveline);
          printf("\n");
      }
      if(cbadres> 0){
         if(verbose>=1){
            printf("%d bad Residues:",cbadres);
            for(i=1;i<=cbadres;i++)
                printf("   %d",ibadres[i]); 
            printf("\n");
         }
         for(ires=i0;ires<nres;ires++) // Use cutmultside in place of cutsidefit
            resbad[ires]=0;
         for(i=1;i<=cbadres;i++)
            resbad[ ibadres[i] ]=1; 
         cutmultside(geo,den1,resbad,geo.savlen,quittrace," bad bone res");
      }
      else if(aveperp >0.7){
         *quittrace= 2; // indicates bad fit 
         if(verbose>=1)
            printf("Bad geo because aveperp = %5.3f> 0.7\n",aveperp);
      }
#if 0 // 15feb01 - Do not use linedist criteria - I think it may be meaningless - e.g. depends on which end measure from
      else if(aveline > geo.maxline){ // only used to cut short fits
         *quittrace= 2; // indicates bad fit 
         if(verbose>=1)
            printf("Bad geo because aveline = %5.3f> %5.3f\n",aveline,geo.maxline);
      }
#endif

#if 1 // for debugging  23feb01
      else if(fractbad > 0.13){
         *quittrace= 2; // indicates bad fit 
         if(verbose>=1)
             printf("Bad geo because fraction bad = %5.3f> 0.13\n",fractbad);
      }
#endif
     } // long enough
     if(*quittrace !=2){// If OK - quit, do not try another iden
        goto quitrout;
     }
// }// iden - no longer let iden vary
 quitrout:;
 strcpy(bone.bonename,origbonename); // put original bonemane in bone1
 readbonefile(bone);
 bone.iexist = -2; // turn off drawing of bone
}
#endif


#if 1 // Call again - use den0.8 and cut any atom not connected
void bonedistgeo2(struct griddata& den1,struct geometry& geo,
         struct bonedata& bone,int *quittrace)
/* same as bonedist - but only applied to one geo and uses geo in place of fit
returns *quittrace = 2 if bad and should be deleted
nbad - counts total number of bad atoms - either bad perp or not connected.
nbadres - counts atoms per res with bad perp
nbadres2  counts atoms per res that are not connected
cbadres counts residues with 2 bad atoms
*/
{
    int i,ii,ires,k,ivec[3],ibone,totnum,cbadres;
    int i0,n0,nres;
    int npt1,npt2,tcount,nbad,nbadres,nbadres2,atomcount;
    float fvec[3],dist,dvec[3];
//    float mdist1,mdist2;
    int ibadres[10];
    int iden;
    char origbonename[300];
    char ch[10];
    float ftemp[NAFIT][5][2];
    int  itemp2[NAFIT][5][2],resbad[NAFIT];
    int iatom[5];
    int ipath[10];
    FILE *fp;

    i0 = geo.nend[0]; //geo[i0] corresponds to fit[0]
    n0 = geo.nend[1]-geo.nend[0];
    nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0)       
    printf("Checking connectivity of geo to bones 0.8 and deleting bad fits\n");
    strcpy(origbonename,bone.bonename);
    tcount = nbad = cbadres = 0; 
 iden = 0; // corresponds to bone 0.8
    bone.minden = fminden[iden]; 
    strcpy(bone.bonename,den1.mapname);
    strcat(bone.bonename,".bone"); // bone  name is mapname with .bone_minden appended
    sprintf(ch,"%1.1f",bone.minden);
    strcat(bone.bonename,ch);
    fp=fopen(bone.bonename,"r");
    iatom[1]=4; // iatom[1,2,3] = N,Ca,C
    iatom[2]=1;
    iatom[3]=2;
    if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone);
    }
    else{
       fclose(fp);
       readbonefile(bone); // This reopens bone file
    }

    for(ires=i0;ires<nres;ires++) // each residue 
       for(i=1;i<=4;i++){
           ftemp[ires][i][0]=ftemp[ires][i][1] = 10000.0;
           for(k=0;k<=1;k++)
              itemp2[ires][i][k]= 1;
       }
    totnum = bone.ptv.end()-bone.ptv.begin();
    for(ibone=0;ibone<totnum;++ibone){
       d1tod3(bone.cmax,bone.ptv[ibone].pos,ivec);
       imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
        if( (nres-1)>=(geo.savlen-1) ){ // if long enough 
         for(ires=i0;ires<nres;ires++) // each residue 
          for(ii=1;ii<=3;ii++){
            i=iatom[ii];
            if( (i==4)&&(ires==i0) ) // N not defined for first residue
               goto nextatom; 
            if( ((i==3)||(i==2))&&(ires== nres-1) ) // O and C not defined for last atom
               goto nextatom; 
            for(k=0;k<=2;k++)
               dvec[k]= fvec[k]-geo.pos[ires][i][k];
            dist = dot(dvec,dvec); // squared dist from bone
            if(dist< ftemp[ires][i][0]){
               ftemp[ires][i][1]=ftemp[ires][i][0];
               ftemp[ires][i][0] = dist;
               itemp2[ires][i][1]=itemp2[ires][i][0];
               itemp2[ires][i][0]=ibone;
            }
            else if(dist< ftemp[ires][i][1]){
               ftemp[ires][i][1] = dist;
               itemp2[ires][i][1]=ibone;
            }   
            nextatom:;
         } // i
       } // long enoug
    } // bone
    npt1 = npt2 = -10; // indicates not assigned
    if( (nres-1)>=(geo.savlen-1) ){ // if long enough 
     cbadres = 0;
     for(ires=i0;ires<nres;ires++){ // each residue
       nbadres2 = 0;
       if(verbose>=2) 
          printf("res %d:",ires);
       for(ii=1;ii<=3;ii++){
           i=iatom[ii];
           if( (i==4)&&(ires==i0) ) // N not defined for first residue
               goto nextatom2; 
           if( ((i==3)||(i==2))&&(ires== nres-1) ) // O and C not defined for last atom
               goto nextatom2;
           atomcount= atomcount+1;
           tcount = tcount+1; // count total number atoms
           npt1= itemp2[ires][i][0];
           npt2= itemp2[ires][i][1];
#if 1 // New  15feb01 - uses connect3 - checks if points connect by total of 4 point and returns path
           connect3(bone,npt1,npt2,ipath); // if ipath[0]= -1 not connected else connected by path - ipath[1]= npt1, ......  ipath[ ipath[0] ] = npt2
           if(ipath[0]== -1){ // not connected - then point is bad independent of perpdist
              nbadres2 = nbadres2+1; // counts  not connected aotms/per residue
              if(verbose>=1) 
                 printf(" res %d atom %d is not connected\n",ires,i);
              goto nextatom2;
           }
#endif
           nextatom2:;
       }// ii 
       if( nbadres2>=1){
           if(verbose>=1)
              printf(" Res #%d not connected - should cut\n",ires);
           cbadres = cbadres+1;
           ibadres[cbadres] = ires;
       }
      } // ires
      if(cbadres> 0){
         if(verbose>=1){
            printf("%d bad Residues:",cbadres);
            for(i=1;i<=cbadres;i++)
                printf("   %d",ibadres[i]); 
            printf("\n");
         }
         for(ires=i0;ires<nres;ires++) // Use cutmultside in place of cutsidefit
            resbad[ires]=0;
         for(i=1;i<=cbadres;i++)
            resbad[ ibadres[i] ]=1; 
         cutmultside(geo,den1,resbad,geo.savlen,quittrace," bad bone res");
      }
     } // long enough
     if(*quittrace !=2){// If OK - quit, do not try another iden
        goto quitrout;
     }
// }// iden - no longer let iden vary
 quitrout:;
 strcpy(bone.bonename,origbonename); // put original bonemane in bone1
 readbonefile(bone);
 bone.iexist = -2; // turn off drawing of bone
}
#endif

void bonedist(struct griddata& den1,struct geometry& geo,struct bonedata& bone)
/* 
I DO NOT THINK THIS IS USED.
finds the minimum dist of each fit atom from nearest bone pt
  Goes through each main chain atom and stores min. dist in
    fit[ifit].ftemp[1,2,3,4]
 nbad are pts that are either not connected or have
   perpdist > 1.0
 perpdist - ave dist from pts that are connected (either directly or though 
               another point.
 nbad  - atoms that are either not connected or have big perpdist
 nbadres counts the atoms in each residue that have bad perpdist 
 nbadres2 counts the atoms in each residue that have bad perpdist or are not connected
Cut residue at end if nbadres2 is bad,
Cut residue in middle if nbadres is bad  
*/
{
    int ifit,i,j,ires,k,k2,ivec[3],ibone,totnum,cbadres;
    int npt1,npt2,ncon,icon,icon2,ncon2,ccount,tcount,nbad,nbadres,nbadres2,atomcount;
    float fvec[3],dist,dvec[3];
//    float mdist1,mdist2;
    float perpdist,perpdist2,linedist,linedist2,pos1[3],pos2[3];
    float aveperp,averesperp,aveline,fractbad,baddist;
    int ibadres[10];
    float mag1,mag2,vec1[3],vec2[3],ang;
    int iden;
    char origbonename[300];
    char ch[10];
    FILE *fp;

    baddist = 1.0;  // perp. bad dist
    printf("Checking agreement of fit to bones 1.0 and deleting bad fits\n");
    strcpy(origbonename,bone.bonename);
    iden = 1;  // this corresponds to bone 1.0
    bone.minden = fminden[iden]; 
    strcpy(bone.bonename,den1.mapname);
    strcat(bone.bonename,".bone"); // bone  name is mapname with .bone_minden appended
    sprintf(ch,"%1.1f",bone.minden);
    strcat(bone.bonename,ch);
    fp=fopen(bone.bonename,"r");
    if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone);
    }
    else{
       fclose(fp);
       readbonefile(bone); // This reopens bone file
    }

    for(ifit=1;ifit<=fitnum[0];ifit++)
      for(ires=0;ires<fit[ifit].nres;ires++) // each residue 
       for(i=1;i<=4;i++){
           fit[ifit].ftemp[ires][i][0]=fit[ifit].ftemp[ires][i][1] = 10000.0;
           for(k=0;k<=1;k++)
              fit[ifit].itemp2[ires][i][k]= 1;
       }
    totnum = bone.ptv.end()-bone.ptv.begin();
    for(ibone=0;ibone<totnum;++ibone){
       d1tod3(bone.cmax,bone.ptv[ibone].pos,ivec);
       imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
       for(ifit=1;ifit<=fitnum[0];ifit++){
        if( (fit[ifit].nres-1)>=(geo.savlen-1) ){ // if long enough 
         for(ires=0;ires<fit[ifit].nres;ires++) // each residue 
          for(i=1;i<=4;i++){
            if( (i==4)&&(ires==0) ) // N not defined for first residue
               goto nextatom; 
            if( ((i==3)||(i==2))&&(ires==fit[ifit].nres-1) ) // O and C not defined for last atom
               goto nextatom; 
            for(k=0;k<=2;k++)
               dvec[k]= fvec[k]-fit[ifit].pos[ires][i][k];
            dist = dot(dvec,dvec); // squared dist from bone
            if(dist< fit[ifit].ftemp[ires][i][0]){
               fit[ifit].ftemp[ires][i][1]=fit[ifit].ftemp[ires][i][0];
               fit[ifit].ftemp[ires][i][0] = dist;
               fit[ifit].itemp2[ires][i][1]=fit[ifit].itemp2[ires][i][0];
               fit[ifit].itemp2[ires][i][0]=ibone;
            }
            else if(dist< fit[ifit].ftemp[ires][i][1]){
               fit[ifit].ftemp[ires][i][1] = dist;
               fit[ifit].itemp2[ires][i][1]=ibone;
            }   
            nextatom:;
         } // i
       } // long enoug
      } // ifit
    } // bone
   for(ifit=1;ifit<=fitnum[0];ifit++){ // Use two close pts to analyze data
    if( (fit[ifit].nres-1)>=(geo.savlen-1) ){ // if long enough 
     tcount = ccount= nbad = cbadres = 0; 
     aveline = aveperp =  0.0;
     for(ires=0;ires<fit[ifit].nres;ires++){ // each residue
       nbadres=nbadres2= atomcount = 0;
       if(verbose>=2) 
          printf("fit %d res %d:",ifit,ires);
       for(i=1;i<=4;i++){
         averesperp= 0.0;
         if(i!=3){ // do not include O
           ang = -1000.0;
           if( (i==4)&&(ires==0) ) // N not defined for first residue
               goto nextatom2; 
           if( ((i==3)||(i==2))&&(ires==fit[ifit].nres-1) ) // O and C not defined for last atom
               goto nextatom2;
           atomcount= atomcount+1;
           tcount = tcount+1;
           npt1= fit[ifit].itemp2[ires][i][0];
           npt2= fit[ifit].itemp2[ires][i][1];
           ncon = bone.ptv[npt1].num; // number of connections
           for(k=0;k<ncon;k++){ // check if close pts are connected
              icon = bone.ptv[npt1].con[k]; // the pt num of the connect point
              if(icon == npt2){ // close pts are connected
                fit[ifit].itemp2[ires][i][1]=1; //indicates connected
                ccount = ccount+1; // count connected
                d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
                imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
                imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                closedist2(pos1,pos2,fit[ifit].pos[ires][i],&perpdist,
                                &linedist);
                goto foundconnect;
              }
           } // k
           for(k=0;k<ncon;k++){ // check if close pts are connected by one pt and find closest perpdist
              icon = bone.ptv[npt1].con[k]; // the pt num of the connect point
              ncon2 = bone.ptv[icon].num; // number of connections
              for(k2=0;k2<ncon2;k2++){                
                icon2 = bone.ptv[icon].con[k2]; // the pt num of the connect point
                if(icon2 == npt2){ // close pts are connected through another pt = icon2
                   ccount = ccount+1; // count connected
                   fit[ifit].itemp2[ires][i][1]=2; //indicates connected through another pt                   
                   d1tod3(bone.cmax,bone.ptv[icon2].pos,ivec);
                   imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
                   d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,fit[ifit].pos[ires][i],&perpdist,
                                &linedist);
                   d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
                   imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
                   closedist2(pos1,pos2,fit[ifit].pos[ires][i],&perpdist2,
                                &linedist2);
                   if(perpdist2<perpdist){
                       perpdist = perpdist2;
                       linedist2 = linedist2;
                   }
                   goto foundconnect;
                 } // if connected
               }// k2
            }// k
           fit[ifit].itemp2[ires][i][1]=0; //if reach here not connected
           nbad = nbad+1;  // not connected points are considered bad no matter what perpdist
           nbadres2 = nbadres2+1; // counts bad perp and not connected aotms
 #if 0
printf("Not connected: npt1 %d  npt2 %d ires %d  i %d  ncon %d: %d %d %d\n",
   npt1,npt2,ires,i,ncon,bone.ptv[npt1].con[0],bone.ptv[npt1].con[1],bone.ptv[npt1].con[2]);
#endif
 #if 0
printf("Not connected: npt1 %d  npt2 %d ires %d  i %d  ncon %d: %d %d %d\n",
   npt1,npt2,ires,i,ncon,bone.ptv[npt1].con[0],bone.ptv[npt1].con[1],bone.ptv[npt1].con[2]);
#endif
           perpdist = fit[ifit].ftemp[ires][i][0];// set perpdist = dist to closedt point
           if(verbose>=2) 
               printf(" noc %5.3f",fit[ifit].ftemp[ires][i][0]);
           goto nextatom2;
           foundconnect:; // skip to here if want to interpolate between close pts
           aveline = aveline +fabs(linedist);
           aveperp = aveperp + perpdist;
           averesperp = averesperp+ perpdist;
           if(perpdist>baddist){
              nbad = nbad+1;
              d1tod3(bone.cmax,bone.ptv[npt1].pos,ivec);// find angle between atom and vec. to 2 nearest pts
              imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
              d1tod3(bone.cmax,bone.ptv[npt2].pos,ivec);
              imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
              for(j=0;j<=2;j++){ 
                 vec1[j] = fit[ifit].pos[ires][i][j]-pos1[j];
                 vec2[j] = pos2[j] - -pos1[j];
              }
              mag1 = sqrt( dot(vec1,vec1) );
              mag2 = sqrt( dot(vec2,vec2) );
              ang = (180.0/PI)*acos( dot(vec1,vec2)/(mag1*mag2) );
           }
           if(verbose>=2){
              if(ang<0.0){// default - good perpdist
                 printf(" %5.3f %5.3f  ", perpdist,linedist); 
                 if(fabs(linedist)>1.0)
                    nbadres2 = nbadres2+1; // counts connected, but bad points
              }
              else if(ang <= 110.0){ // then bad dist and points behave as if connected              
                 nbadres = nbadres+1; // counts connected, but bad points
                 nbadres2 = nbadres2+1; // counts connected, but bad points
                 printf("ang: %5.3f %5.3f %5.3f  ",ang, perpdist,linedist);
              }
              else{ // bad perp dist but really not connected
                 nbadres2 = nbadres2+1; // counts connected, but bad points
                 printf(" noc  ang: %5.3f %5.3f",ang,fit[ifit].ftemp[ires][i][0]);
              }
           }
           nextatom2:;
         } // not O
       }// i 
       if(verbose>=2)
          printf("\n");
       if( (ires<=1)||(ires>=fit[ifit].nres-2) ){
printf(" For END  res %d  averesperp = %5.3f\n",ires,averesperp);
         if( (float)nbadres2/atomcount>0.6 ){//Cut end residue if bad perp or not connected.
           cbadres = cbadres+1;
           ibadres[cbadres] = ires;
           if(verbose>=2) 
              printf(" 60 percent atoms are bad for end res %d\n ",ires);
         }
       }   
       else if( (float)nbadres/atomcount>0.6){//Cut middle residue if bad perp
           cbadres = cbadres+1;
           ibadres[cbadres] = ires;
           if(verbose>=2) 
              printf(" 60 percent atoms are bad for res %d ",ires);
           if(cbadres>5){
              fit[ifit].itemp[0]=0; // indicates bad fit 
              printf("Bad fit because more than 5 bad res in bonedist\n");
              goto nextifit;
           } 
       }
      } // ires
       aveperp = aveperp/ccount;
       aveline = aveline/ccount;
       fractbad = (float)nbad/tcount;
       if(verbose>=2){ 
          printf("fract bad = %5.3f fract con = %5.3f  aveperp = %5.3f aveline = %5.3f\n",
               fractbad,(float)ccount/tcount,aveperp,aveline);
          printf("\n");
      }
      if(cbadres> 0){
         fit[ifit].itemp[0]=2; //indicates bad residue that must be removed
         fit[ifit].itemp2[0][0][0]= cbadres;
         for(i=1;i<=cbadres;i++)
            fit[ifit].itemp2[i][0][0]= ibadres[i];
         if(verbose>=2){
            printf("%d bad Residues:",fit[ifit].itemp2[0][0][0]);
            for(i=1;i<=fit[ifit].itemp2[0][0][0];i++)
                 printf("   %d",fit[ifit].itemp2[i][0][0]); 
            printf("\n");
         }
         cutsidefit(geo,den1,ifit," bones"); 
      }
      if(aveperp >0.7){
         fit[ifit].itemp[0]=0; // indicates bad fit 
         printf("Bad fit %d because aveperp = %5.3f\n",ifit,aveperp);
      }
      else if(fractbad > 0.25){
         fit[ifit].itemp[0]=0; // indicates bad fit 
         printf("Bad fit %d because fraction bad = %5.3f\n",ifit,fractbad);
      }
      else 
         fit[ifit].itemp[0]=1; // indicates good fit 
     } // long enough
     nextifit:;
   } //ifit
   strcpy(bone.bonename,origbonename); // put original bonemane in bone1
   readbonefile(bone);
   bone.iexist = -2; // turn off drawing of bone
}


#if 0 // New version  - IS NOT USED - THE NEW checkterm() served same purpose- uses recursive form of checkifbranchrec()
void checkoneendbranch(struct griddata& den1,struct geometry& geo,
        struct bonedata& bone1,int *quittrace)
/*
NOW USED IN PLACE OF checkendbranch
Similar to checkendbranch - but only checks one geo fit
returns quittrace ==2 if fit cut too short
  NOTE: Allow geo.savlen to be one less than in generall routine
   This allows to cut so that length = original geo.savlen -1
Add new condition:  for helix - check if cterm C-O is on main branch and Ca is not
checkifbranchrec RETUrns ibranch =
   0 if not on banch
   1 if on branch, but branch does not go to another branch
   2 if on branch, and branch  goes to another branch
   3 if nearnum is a branch pt
   4  if all points moving closer to pos1 - something wrong - treat as possible branch
 If checknear = 1  (usually used for end Ca pt) then returns ibranch =
   3 if nearnum is close to a branch pt
*/
{
  int iden,mbranch,sbranch,iremove,itemp[2],int pt[3];
  int i0,n0,nres,idenmax,idenmin,testside,endvalue;
  char ch[10];
  FILE *fp;
  char origbonename[300];
  int istore[2];

//  ext.sidelim=7; // If number of side pts > sidelim, then assume it is a main branch
  ext.sidelim=4; // If number of side pts > sidelim, then assume it is a main branch
  strcpy(origbonename,bone1.bonename);
  itemp[0]=itemp[1] = -1;  // indicates not yet assigned
//  for(iden = 2;iden<=5;iden++){// den 1.2 to den 1.8
  idenmax = 4; // corresponds to 1.6
  idenmin = 2; // corresponds to 1.2
//  testside = geo.maxside;
  testside =3; // NOTE - must be >=2
  istore[0]=istore[1] = 0;  // if = 1, indicates possible bad end
  checkterm(den1,geo,bone1,&endvalue); // just checks if end on branched chain
  if(endvalue ==1) // bad nterm 
     itemp[0]=1;
  else if(endvalue==2) // bad cterm
     itemp[1]=1;
  else if(endvalue==3){ // both ends bad
     itemp[0]=itemp[1]=1;
  }
  for(iden = idenmin;iden<=idenmax;iden++){// den 1.6 to den 1.2
     bone1.minden = fminden[iden]; 
     strcpy(bone1.bonename,den1.mapname);
     strcat(bone1.bonename,".bone"); // bone  name is mapname with .bone_minden appended
     sprintf(ch,"%1.1f",bone1.minden);
     strcat(bone1.bonename,ch);
     fp=fopen(bone1.bonename,"r");
     if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone1);
     }
     else{
       fclose(fp);
       readbonefile(bone1); // This reopens bone file
     }
     i0 = geo.nend[0];
     n0 = geo.nend[1]-geo.nend[0];
     nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0) 
     if( (itemp[0] == -1)&&(geo.iside[i0+1]>=testside) ){ // then test nterm end      
#if 0 // Procedure to see if end followed O - written but never used or debugged
        checkside_ca_branch(1,den1,geo,bone1,geo.pos[i0][1],geo.pos[i0+1][1],
           geo.pos[i0+1][6],&asbranch);
        if(asbranch==1){// indicates bad nterm 
           itemp[0]=1;
           goto checkcterm;
        }
         checkobranch(1,den1,geo,bone1,geo.pos[i0][2],geo.pos[i0][1],geo.pos[i0][3],
             &obranch);
#endif
         checkifbranchrec(1,den1,geo,bone1,geo.pos[i0+1][4],
                geo.pos[i0][1],&mbranch,pt); //"1" means see if a branch pt is with 2 of pos2 
         if(mbranch== -1){ // Ca and N not connected - bad end
            if(verbose>=1)
               printf("Ca and N are not connected at nterm - remove end\n");
            itemp[0] = 1; // bad end, remove nterm 
            goto checkcterm;
         }    
         checkifbranchrec(0,den1,geo,bone1,geo.pos[i0+1][1],
                 geo.pos[i0+1][6],&sbranch,pt);
         if( (mbranch==3)||( (sbranch==3)&&(mbranch!=1) ) ){ // Store possibility of bad end
              if(verbose>=1)
                  printf("iden = %d Bad nterm because mbranch =%d  sbranch = %d\n",iden,mbranch,sbranch);
             istore[0] = 1;
         }
         if((mbranch==0)&&(sbranch!=0)){
              if(verbose>=1)
                  printf(" Bad nterm because Side on branch and Ca is not\n");
              itemp[0] = 1; // bad end, remove nterm 
         }    
         else if( (mbranch!=0)&&(sbranch==0) )
              itemp[0] = 0; // good end 
     } // itemp[0]
     checkcterm:;
     if(itemp[1] == -1){ // then test cterm end      
         if((geo.rtype[nres-2]==2)&&(geo.iside[nres-2]>=testside) ){ // check cterm side only for sheet
#if 0 // Procedure to see if end followed O - written but never used or debugged
            checkside_ca_branch(1,den1,geo,bone1,geo.pos[nres-1][1],geo.pos[nres-2][1],
                  geo.pos[nres-2][6],&asbranch);
          if(asbranch==1){// indicates bad nterm 
             itemp[1]=1;
             goto endcterm;
          }
            checkobranch(1,den1,geo,bone1,geo.pos[nres-2][2],geo.pos[nres-1][1],
                 geo.pos[nres-2][3],&obranch);
#endif
            checkifbranchrec(1,den1,geo,bone1,geo.pos[nres-2 ][2],
                  geo.pos[nres-1 ][1],&mbranch,pt);
            if(mbranch== -1){ // Ca and N not connected - bad end
               if(verbose>=1)
                 printf("Ca and C are not connected at cterm - remove end\n");
               itemp[1] = 1; // bad end, remove cterm 
               goto endcterm;
            }    
            checkifbranchrec(0,den1,geo,bone1,geo.pos[nres-2 ][1],
                   geo.pos[nres-2 ][6],&sbranch,pt);
            if( (mbranch==3)||( (sbranch==3)&&(mbranch!=1) ) ){
              if(verbose>=1)
                  printf("iden = %d Bad cterm because mbranch =%d  sbranch = %d\n",iden,mbranch,sbranch);
              istore[1] = 1;
            }
            if( (mbranch==0)&&(sbranch!=0) ){
                 if(verbose>=1)
                    printf(" Bad  cterm because side  on branch, Ca is not\n");
                 itemp[1] = 1; // bad end, remove nterm  
            }   
            else if( (mbranch!=0)&&(sbranch==0) )
                  itemp[1] = 0; // good end 
         } // testside
         if( (geo.rtype[nres-2]==1)&&(iden ==idenmin) ){  // for helix -  check C-O branch at c term end at lowest den
            checkifbranchrec(1,den1,geo,bone1,geo.pos[nres-2 ][2],
                  geo.pos[nres-1 ][1],&mbranch);
            checkifbranchrec(0,den1,geo,bone1,geo.pos[nres-2 ][2],
                    geo.pos[nres-2 ][3],&sbranch,pt); // C-O branch
            if( (mbranch==0)&&(sbranch!=0) ){
                if(verbose>=1)
                    printf(" Bad  cterm becausehelix CO on branch, Ca is not\n");
                itemp[1] = 1; // bad end, remove nterm 
            } 
        } // if helix
        endcterm:;
     } // test cterm
   } //iden
   if(itemp[0] == -1){ // if still not set (i.e. both side and end  have braanch - assume OK
      if(istore[0]==0)
          itemp[0]=0;
      else
          itemp[0]=1;
   }
   if(itemp[1] == -1){ // if still not set (i.e. both side and end  have braanch - assume OK
      if(istore[1]==0)
          itemp[1]=0;
      else
          itemp[1]=1;
   }
   if( (itemp[0]==1)&&(itemp[1]==1) ){// remove both ends
       printf(" Cutting nterm and cterm end because bad checkifbranch\n");      
       iremove = 3; 
   }
   else if(itemp[0]==1){// remove nterm
       printf(" Cutting nterm end because bad checkifbranch\n");
       iremove = 1;
   }
   else if (itemp[1] ==1){ // remove cterm 
       printf(" Cutting cterm end because bad checkifbranch\n");
       iremove = 2;
   }
   else // do not remove
      iremove = 0;
   cutend(den1,geo,quittrace,iremove);
   if( (*quittrace!=2)&&(iremove!=0) ){ //if not bad fit and  end was cut, call again
       checkoneendbranch(den1,geo,bone1,quittrace);
   }
   strcpy(bone1.bonename,origbonename); // put original bonemane in bone1
   readbonefile(bone1);
   bone1.iexist = -2; // turn off drawing of bone
}     
#endif

#if 1 // HALF NEW  version - does not use recursive form of checkifbranch itself
      // but calls new version of checkterm() which calls checkifbranchrec
void checkoneendbranch(struct griddata& den1,struct geometry& geo,
        struct bonedata& bone1,int *quittrace)
/*
NOW USED IN PLACE OF checkendbranch
Similar to checkendbranch - but only checks one geo fit
returns quittrace ==2 if fit cut too short
  NOTE: Allow geo.savlen to be one less than in generall routine
   This allows to cut so that length = original geo.savlen -1
Add new condition:  for helix - check if cterm C-O is on main branch and Ca is not
checkifbranch RETUrns ibranch =
   0 if not on banch
   1 if on branch
   2 if nearnum is a branch pt
   3  if all points moving closer to pos1 - something wrong
 If checknear = 1  (usually used for end Ca pt) then returns ibranch =
   2 if nearnum is close to a branch pt
*/
{
  int iden,mbranch,sbranch,iremove,itemp[2];
  int i0,n0,nres,idenmax,idenmin,testside,endvalue;
  char ch[10];
  FILE *fp;
  char origbonename[300];
  int istore[2];

//  ext.sidelim=7; // If number of side pts > sidelim, then assume it is a main branch
  ext.sidelim=4; // If number of side pts > sidelim, then assume it is a main branch
  strcpy(origbonename,bone1.bonename);
  itemp[0]=itemp[1] = -1;  // indicates not yet assigned
//  for(iden = 2;iden<=5;iden++){// den 1.2 to den 1.8
//  idenmax = 4; // corresponds to 1.6
  idenmax = 6; //19feb01  corresponds to 2.0
  idenmin = 2; // corresponds to 1.2
//idenmin=idenmax;
//  testside = geo.maxside;
  testside =3; // NOTE - must be >=2 so that geo.pos[i0+1][6] is defined below
  istore[0]=istore[1] = 0;  // if = 1, indicates possible bad end
  checkterm(den1,geo,bone1,&endvalue); // just checks if end on branched chain
  if(endvalue ==1) // bad nterm 
     itemp[0]=1;
  else if(endvalue==2) // bad cterm
     itemp[1]=1;
  else if(endvalue==3){ // both ends bad
     itemp[0]=itemp[1]=1;
  }
  for(iden = idenmin;iden<=idenmax;iden++){// den 2.0 to den 1.2
     bone1.minden = fminden[iden]; 
     strcpy(bone1.bonename,den1.mapname);
     strcat(bone1.bonename,".bone"); // bone  name is mapname with .bone_minden appended
     sprintf(ch,"%1.1f",bone1.minden);
     strcat(bone1.bonename,ch);
     fp=fopen(bone1.bonename,"r");
     if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone1);
     }
     else{
       fclose(fp);
       readbonefile(bone1); // This reopens bone file
     }
     i0 = geo.nend[0];
     n0 = geo.nend[1]-geo.nend[0];
     nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0) 
     if( (itemp[0] == -1)&&(geo.iside[i0+1]>=testside) ){ // then test nterm end      
//     if(itemp[0] == -1){ //18feb01 -CANNOT USE - SEE NOTE 19feb01 - in trace.notes   == always check side branches then test nterm end      
#if 0 // Procedure to see if end followed O - written but never used or debugged
        checkside_ca_branch(1,den1,geo,bone1,geo.pos[i0][1],geo.pos[i0+1][1],
           geo.pos[i0+1][6],&asbranch);
        if(asbranch==1){// indicates bad nterm 
           itemp[0]=1;
           goto checkcterm;
        }
         checkobranch(1,den1,geo,bone1,geo.pos[i0][2],geo.pos[i0][1],geo.pos[i0][3],
             &obranch);
#endif
         checkifbranch(1,den1,geo,bone1,geo.pos[i0+1][4],
                geo.pos[i0][1],&mbranch); //"1" means see if a branch pt is with 2 of pos2
         checksidebranch(0,den1,geo,bone1,geo.pos[i0+1][1],
                 geo.pos[i0+1][6],&sbranch);
         if( (mbranch==3)||( (sbranch==3)&&(mbranch!=1) ) ){ // Store possibility of bad end
              if(verbose>=2)
                  printf("iden = %d  Possible Bad nterm because mbranch =%d  sbranch = %d\n",iden,mbranch,sbranch);
             istore[0] = 1;
         }
         if(mbranch<=0){ // main not on branch
            if( (sbranch>0) || ((mbranch>= -1)&&(sbranch<= -3)) ){ // if main on short and side on long side 
                if(verbose>=1)
                     printf(" Bad nterm because sbranch=%d and mbranch = %d iden = %d\n",sbranch,mbranch,iden);
              itemp[0] = 1; // bad end, remove nterm 
            }
         } 
#if 0 // old before 18oct00
         if((mbranch<=0)&&(sbranch>0)){
              if(verbose>=1)
                  printf(" Bad nterm because Side on branch and Ca is not\n");
              itemp[0] = 1; // bad end, remove nterm 
         } 
#endif   
         else if( (mbranch>0)&&(sbranch<=0) )
              itemp[0] = 0; // good end
         else if( (sbranch==1)&&(geo.trmetfit[0]==1) ){
            printf("Bad nterm end because side on branch and end met another fit\n");
            itemp[0]=1;// bad end, remove nterm 
         }
     } // itemp[0] - check n term side

     if( (itemp[0] == -1)&&(geo.rtype[i0+1]==2) ){ //9mar01  If sheet check if nterm C-O on branch     
         checkifbranch(1,den1,geo,bone1,geo.pos[i0+1][4],
                geo.pos[i0][1],&mbranch); //"1" means see if a branch pt is with 2 of pos2
         checksidebranch(0,den1,geo,bone1,geo.pos[i0][2],
                 geo.pos[i0][3],&sbranch);
         if( (mbranch!=1)&&(sbranch==1) ){
                if(verbose>=1)
                    printf(" Bad  nterm sheet nterm CO on branch, Ca is not iden = %d\n",iden);
                itemp[0] = 1; // bad end, remove nterm 
         } 
      } // sheet

     checkcterm:;
     if(itemp[1] == -1){ // then test cterm end      
         if((geo.rtype[nres-2]==2)&&(geo.iside[nres-2]>=testside) ){ // check cterm side only for sheet
//         if(geo.rtype[nres-2]==2){ //19feb01 - check no matter lenght of side chain check cterm side only for sheet
#if 0 // Procedure to see if end followed O - written but never used or debugged
            checkside_ca_branch(1,den1,geo,bone1,geo.pos[nres-1][1],geo.pos[nres-2][1],
                  geo.pos[nres-2][6],&asbranch);
          if(asbranch==1){// indicates bad nterm 
             itemp[1]=1;
             goto endcterm;
          }
            checkobranch(1,den1,geo,bone1,geo.pos[nres-2][2],geo.pos[nres-1][1],
                 geo.pos[nres-2][3],&obranch);
#endif
            checkifbranch(1,den1,geo,bone1,geo.pos[nres-2 ][2],
                  geo.pos[nres-1 ][1],&mbranch);
            checksidebranch(0,den1,geo,bone1,geo.pos[nres-2 ][1],
                   geo.pos[nres-2 ][6],&sbranch);
            if( (mbranch==3)||( (sbranch==3)&&(mbranch!=1) ) ){
              if(verbose>=2)
                  printf("iden = %d  Possible Bad cterm because mbranch =%d  sbranch = %d\n",iden,mbranch,sbranch);
              istore[1] = 1;
            }
            if(mbranch<=0){ // main not on branch
              if( (sbranch>0) || ((mbranch>= -1)&&(sbranch<= -3)) ){ // if main on short and side on long side 
                if(verbose>=1)
                     printf(" Bad cterm because sbranch=%d and mbranch = %d iden = %d\n",sbranch,mbranch,iden);
                itemp[1] = 1; // bad end, remove cterm 
              }
           } 
#if 0 // Old before 18oct00
            if( (mbranch<=0)&&(sbranch>0) ){
                 if(verbose>=1)
                    printf(" Bad  cterm because side  on branch, Ca is not\n");
                 itemp[1] = 1; // bad end, remove nterm  
            } 
#endif  
            else if( (mbranch>0)&&(sbranch<=0) )
                  itemp[1] = 0; // good end 
#if 1
            else if( (sbranch==1)&&(geo.trmetfit[1]==1) ){
               printf("Bad cterm end because side on branch and end met another fit\n");
              itemp[1]=1;// bad end, remove cterm 
            }
#endif
         } // testside
         if( (geo.rtype[nres-2]==1)&&(iden ==idenmin) ){  // for helix -  check C-O branch at c term end at lowest den
            checkifbranch(1,den1,geo,bone1,geo.pos[nres-2 ][2],
                  geo.pos[nres-1 ][1],&mbranch);
            checksidebranch(0,den1,geo,bone1,geo.pos[nres-2 ][2],
                    geo.pos[nres-2 ][3],&sbranch); // C-O branch
            if( (mbranch<=0)&&(sbranch>0) ){
                if(verbose>=1)
                    printf(" Bad  cterm becausehelix CO on branch, Ca is not iden = %d\n",iden);
                itemp[1] = 1; // bad end, remove nterm 
            } 
        } // if helix
        endcterm:;
     } // test cterm
   } //iden
   if(itemp[0] == -1){ // if still not set (i.e. both side and end  have braanch - assume OK
      if(istore[0]==0)
          itemp[0]=0;
      else
          itemp[0]=1;
   }
   if(itemp[1] == -1){ // if still not set (i.e. both side and end  have braanch - assume OK
      if(istore[1]==0)
          itemp[1]=0;
      else
          itemp[1]=1;
   }
   if( (itemp[0]==1)&&(itemp[1]==1) ){// remove both ends
       printf(" Cutting nterm and cterm end because bad checkifbranch\n");      
       iremove = 3; 
   }
   else if(itemp[0]==1){// remove nterm
       printf(" Cutting nterm end because bad checkifbranch\n");
       iremove = 1;
   }
   else if (itemp[1] ==1){ // remove cterm 
       printf(" Cutting cterm end because bad checkifbranch\n");
       iremove = 2;
   }
   else // do not remove
      iremove = 0;
   cutend(den1,geo,quittrace,iremove);
   if( (*quittrace!=2)&&(iremove!=0) ){ //if not bad fit and  end was cut, call again
       checkoneendbranch(den1,geo,bone1,quittrace);
   }
   strcpy(bone1.bonename,origbonename); // put original bonemane in bone1
   readbonefile(bone1);
   bone1.iexist = -2; // turn off drawing of bone
}     
#endif

#if 0 // weak version - requires branch !=0
void checkterm(struct griddata& den1,struct geometry& geo,
        struct bonedata& bone1,int *endvalue)
/*
Checks if end lies on main branch
Returns endvalue = 0 if both ends good, 1 if bad nterm 2 bad cterm 3 bothends bad
*/
{
  int iden,nbranch,cbranch;
  int i0,n0,nres;
  char ch[10];
  FILE *fp;
  char origbonename[300];

  *endvalue = 0;  // default - indicates both ends good.
//  ext.sidelim=7; // If number of side pts > sidelim, then assume it is a main branch
  ext.sidelim=4; // If number of side pts > sidelim, then assume it is a main branch
  strcpy(origbonename,bone1.bonename);
//  for(iden = 2;iden<=5;iden++){// den 1.2 to den 1.8
//  iden = 0; // den0.8
  iden = 1; // den1.0  - this is strong criteria - removes a few good ends- but removes many bad ends
     bone1.minden = fminden[iden]; 
     strcpy(bone1.bonename,den1.mapname);
     strcat(bone1.bonename,".bone"); // bone  name is mapname with .bone_minden appended
     sprintf(ch,"%1.1f",bone1.minden);
     strcat(bone1.bonename,ch);
     fp=fopen(bone1.bonename,"r");
     if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone1);
     }
     else{
       fclose(fp);
       readbonefile(bone1); // This reopens bone file
     }
     i0 = geo.nend[0];
     n0 = geo.nend[1]-geo.nend[0];
     nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0) 
     checkifbranch(1,den1,geo,bone1,geo.pos[i0+1][4],
                geo.pos[i0][1],&nbranch); //"1" means see if a branch pt is with 2 of pos2 
     checkifbranch(1,den1,geo,bone1,geo.pos[nres-2 ][2],
                  geo.pos[nres-1 ][1],&cbranch);
     if(nbranch<=0){
         *endvalue = 1; //nterm bad
         if(verbose>=1)
             printf("BAD NTERM\n");
     }
     if(cbranch<=0){
         if(*endvalue==1)
            *endvalue = 3; // both ends bad
         else
            *endvalue = 2; // cterm bad
         if(verbose>=1)
             printf("BAD CTERM\n");
     }
     if(verbose>=2) 
         printf("At nterm: branch = %d   cterm:  branch = %d\n",nbranch,cbranch);
}
#endif

#if 1 // New  version  - calls checkifbranchrec
void checkterm(struct griddata& den1,struct geometry& geo,
        struct bonedata& bone1,int *endvalue)
/*
NEW - checkifbranchrec return the point numbers of nearpoint and the nextpt on main branch
   Use this to check if this branch is parallel to C-O at nterm end.
Checks if end lies on main branch
Returns endvalue = 0 if both ends good, 1 if bad nterm 2 bad cterm 3 bothends bad
*/
{
  int iden,idenmin,nbranch,cbranch,pt[3];
  int i0,n0,nres;
  char ch[10];
  FILE *fp;
  char origbonename[300];
//  float pos1[3],pos2[3],mag2,mag1,vec[3],cavec[3],ang;

  *endvalue = 0;  // default - indicates both ends good.
  ext.sidelim=7; // If number of side pts > sidelim, then assume it is a main branch
//  ext.sidelim=4; // If number of side pts > sidelim, then assume it is a main branch
  strcpy(origbonename,bone1.bonename);
//  for(iden = 2;iden<=5;iden++){// den 1.2 to den 1.8
//  iden = 0; // den0.8
  idenmin = 1; // Used up to 20feb01 = den1.0  - this is strong criteria - removes a few good ends- but removes many bad ends
 cbranch=nbranch= 0; // default - bad end
 for(iden=1;iden>=idenmin;iden--){
     bone1.minden = fminden[iden]; 
     strcpy(bone1.bonename,den1.mapname);
     strcat(bone1.bonename,".bone"); // bone  name is mapname with .bone_minden appended
     sprintf(ch,"%1.1f",bone1.minden);
     strcat(bone1.bonename,ch);
     fp=fopen(bone1.bonename,"r");
     if(fp == NULL){ // file does not exist NOTE: In linux - closing a NULL file causes crash
//         fclose(fp);
         findbone(den1,bone1);
     }
     else{
       fclose(fp); 
       readbonefile(bone1); // This reopens bone file
     }
     i0 = geo.nend[0];
     n0 = geo.nend[1]-geo.nend[0];
     nres = i0+n0+1;// corresponds to fit.nres   total number of res = n0+1  (0,1,....n0) 
     if(nbranch==0)
        checkifbranchrec(1,den1,geo,bone1,geo.pos[i0+1][4],
                geo.pos[i0][1],&nbranch,pt); //"1" means see if a branch pt is with 2 of pos2 
     if(cbranch==0)
         checkifbranchrec(1,den1,geo,bone1,geo.pos[nres-2 ][2],
                  geo.pos[nres-1 ][1],&cbranch,pt);
  }
     if(nbranch==0){
         *endvalue = 1; //nterm bad
         if(verbose>=1)
             printf("BAD NTERM\n");
     }
#if 0 // NOT used - does not work well
     else{ // NEW 13oct00 - if branch at nterm end - see if antiparallel to C-Ca
        d1tod3(bone1.cmax,bone1.ptv[ pt[0] ].pos,ivec);
        imcv(den1,ivec,pos1);// convert int grid ivec to float world co-ord fvec
        d1tod3(bone1.cmax,bone1.ptv[ pt[1] ].pos,ivec);
        imcv(den1,ivec,pos2);// convert int grid ivec to float world co-ord fvec
        for(k=0;k<=2;k++){
           vec[k] = pos2[k]-pos1[k];
           cavec[k] = geo.pos[i0][1][k]-geo.pos[i0][2][k]; // C-Ca vector
        } 
        mag1 = sqrt( dot(vec,vec) );
        mag2 = sqrt( dot(cavec,cavec) );
        ang = acos( dot(vec,cavec)/(mag1*mag2) );
        printf(" Nterm C-Ca ang (degrees) = %5.3f\n",ang*180.0/PI);
     }
#endif
     if(cbranch==0){
         if(*endvalue==1)
            *endvalue = 3; // both ends bad
         else
            *endvalue = 2; // cterm bad
         if(verbose>=1)
             printf("BAD CTERM\n");
     }
     if(verbose>=2) 
         printf("At nterm: branch = %d   cterm:  branch = %d\n",nbranch,cbranch);
}
#endif

void checkendbranch(struct griddata& den1,struct geometry& geo,
        struct bonedata& bone1,int startfitnum)
/*  NOTE NOT USED ANYMORE  (CALLED BY CHECKEND WHICH IS NOT USED)
Checks if nterm end in on main branch:
 NOTE: this assume that geo.maxside>=2  
 returns fit[i].itemp[0,1]  (0 for nterm, 1 for cterm)
 Initially itemp set = -1
   Only check if = -1.
  If main on branch bone and side not - then it is OK return itemp = 0
  If neither main and side on branch point - or side on branch and main not
     then then bad return itemp = 1
 Starts with den1.2 - do above check on all fits
 then goto den 1.4, 1.6,1.8 
   If itemp = -1, then carry out above check
 If after 1.8, itemp still = -1, (i.e. both main and side on branch asume it
    is OK and set itemp = 0;
*/
{
  int iden,ifit,mbranch,sbranch;
  char ch[10];
  FILE *fp;
  struct pepfit fit1;
  char origbonename[300];

  fit1.nalloc = -10; // initial default - indicates have not yet allocated 
  ext.sidelim=7; // If number of side pts > sidelim, then assume it is a main branch
  for(ifit=startfitnum;ifit<=fitnum[0];ifit++) // set itemp to default value
     fit[ifit].itemp[0] = fit[ifit].itemp[1] = -1;
  strcpy(origbonename,bone1.bonename);
  for(iden = 2;iden<=5;iden++){// den 1.2 to den 1.8
     bone1.minden = fminden[iden]; 
     strcpy(bone1.bonename,den1.mapname);
     strcat(bone1.bonename,".bone"); // bone  name is mapname with .bone_minden appended
     sprintf(ch,"%1.1f",bone1.minden);
     strcat(bone1.bonename,ch);
     fp=fopen(bone1.bonename,"r");
     if(fp == NULL){ // file does not exist
//         fclose(fp);
         findbone(den1,bone1);
     }
     else{
       fclose(fp);
       readbonefile(bone1); // This reopens bone file
     }
     for(ifit=startfitnum;ifit<=fitnum[0];ifit++){ // check end residues and write out new fit file
      if( (fit[ifit].nres-1)>=(geo.savlen-1) ){ // if long enough 
       copyfittogeo(geo,fit[ifit]);
       if( (fit[ifit].itype[1]==2)&&(fit[ifit].nside[1]>=geo.maxside) ){ // check if nterm end Ca on branched mainbone
         if(fit[ifit].itemp[0]== -1){ // itemp not yet assigned
           checkifbranch(1,den1,geo,bone1,fit[ifit].pos[0][2],
                fit[ifit].pos[0][1],&mbranch); //"1" means see if a branch pt is with 2 of pos2  
           checkifbranch(0,den1,geo,bone1,fit[ifit].pos[1][1],
                fit[ifit].pos[1][6],&sbranch);
//printf("sheet nterm: ifit = %d iden = %d mbranch = %d sbranch = %d\n",ifit,iden,mbranch,sbranch);
           if( (mbranch<=0)&&(sbranch>0) ){
              fit[ifit].itemp[0] = 1; // bad end, remove     
           }
           else if( (mbranch>0)&&(sbranch<=0) ){
              fit[ifit].itemp[0] = 0; // good end    
           }
         }
       }
       if( (fit[ifit].itype[1]==1)&&(fit[ifit].nside[1]>=geo.maxside) ){// helix
         if(fit[ifit].itemp[0]== -1){ // itemp not yet assigned
           checkifbranch(1,den1,geo,bone1,fit[ifit].pos[0][2],
             fit[ifit].pos[0][1],&mbranch);
           checkifbranch(0,den1,geo,bone1,fit[ifit].pos[1][1],
             fit[ifit].pos[1][6],&sbranch);
//printf("helix nterm ifit = %d iden = %d mbranch = %d sbranch = %d\n",ifit,iden,mbranch,sbranch);   
           if( (mbranch<=0)&&(sbranch>0) ){
              fit[ifit].itemp[0] = 1; // bad end, remove     
           }
           else if( (mbranch>0)&&(sbranch<=0) ){
              fit[ifit].itemp[0] = 0; // good end    
           }
         }
       }
       if( (fit[ifit].itype[1]==2)&&(fit[ifit].nside[fit[ifit].nres-2 ]>=geo.maxside) ){ //for sheet check if cterm end Ca on branched mainbone
         if(fit[ifit].itemp[1]== -1){ // itemp not yet assigned
           checkifbranch(1,den1,geo,bone1,fit[ifit].pos[fit[ifit].nres-1 ][4],
              fit[ifit].pos[ fit[ifit].nres-1 ][1],&mbranch);
           checkifbranch(0,den1,geo,bone1,fit[ifit].pos[fit[ifit].nres-2 ][1],
              fit[ifit].pos[ fit[ifit].nres-2 ][6],&sbranch);
//printf("sheet cterm ifit = %d iden = %d mbranch = %d sbranch = %d\n",ifit,iden,mbranch,sbranch);
           if( (mbranch<=0)&&(sbranch>0) ){
              fit[ifit].itemp[1] = 1; // bad end, remove     
           }
           else if( (mbranch>0)&&(sbranch<=0) ){
              fit[ifit].itemp[1] = 0; // good end    
           }
         }//fit[ifit].itemp[1]== -1
       }//if(fit[ifit].itype[1]==2)
      } // if long enough
     }// ifit
   } //iden
   
#if 1 // now removed in checkend
   for(ifit=startfitnum;ifit<=fitnum[0];ifit++){ // remove end residues 
      if(fit[ifit].itemp[0] == -1) // if still not set (i.e. both side and end  have braanch - assume OK
          fit[ifit].itemp[0]=0;
      if(fit[ifit].itemp[1] == -1) // if still not set (i.e. both side and end  have braanch - assume OK
          fit[ifit].itemp[1]=0;
      if(fit[ifit].itemp[0]==1){ // remove nterm end
         printf("Bad nterm end of fit %d because not on brached main\n",ifit);   
      }
      if(fit[ifit].itemp[1]==1){ // remove cterm end
          printf("Bad cterm end of fit %d because not on brached main\n",ifit);   
      }
   }
#endif
   if(fit1.nalloc!= -10)
        freepepfit(fit1);
   strcpy(bone1.bonename,origbonename); // put original bonemane in bone1
   readbonefile(bone1);
   bone1.iexist = -2; // turn off drawing of bone
}     


int checkendden(struct griddata& den1,struct geometry& geo,
    float pos[NAFIT][NAA][3])
/*
Returns 0 if OK, 1 if bad nterm, 2 if bad cterm 3 if both bad
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ires,ireturn;
   float badden;
   float aden;
   float aveden[3],ca1den[3],ca2den[3],oden[3],nden[3];
   float cden[3],cbden[3],noden[3];

   i=0;
   ireturn = 0;
   aveden[i] = 0;
   badden = 0.25*maidfile.denscale;
   ires = geo.nend[0]; // n term end
   ca1den[i] = aden= aveatom(den1,geo,1,ires,pos);
   if(aden< badden){
      if(verbose>=1)
        printf("Bad nterm end because bad Ca den = %5.3f at res %d\n",aden,ires);
      ireturn = 1; 
      goto cend;
   }     
   aveden[i] = aveden[i] +aden;
   oden[i] = aden = aveatom(den1,geo,3,ires,pos);
   if(aden<  badden){
      if(verbose>=1)
        printf("Bad nterm end because bad O den = %5.3f at res %d\n",aden,ires);
      ireturn = 1;      
      goto cend;
   }
   if(geo.rtype[ geo.nend[0]+1 ]==2){// if sheet
      if(ca1den[0]/oden[0] < 0.47){ 
         printf(" Cutting sheet nterm end because bad caden/oden ratio = %5.3f\n",ca1den[0]/oden[0]);    
//printf("WARNING FOR DEBUGGING - NOT USING CA1DEN/ODEN RATIO\n");
#if 1
         ireturn = 1;      
         goto cend;
#endif
      }
   }
   aveden[i] = aveden[i] +aden;
   cden[i]=aden= aveatom(den1,geo,2,ires,pos);
   if(aden<  badden){
      if(verbose>=1)
        printf("Bad nterm end because bad C den = %5.3f at res %d\n",aden,ires);
      ireturn = 1;      
      goto cend;
   }     
   aveden[i] = aveden[i] +aden;
   ires = ires+1;
   cbden[i]=aden= aveatom(den1,geo,5,ires,pos);
#if 0 // New 30may00 - do not include Cb
   if(aden<  badden){
      ireturn = 1;      
      goto cend;
   }
#endif     
   nden[i]=aden= aveatom(den1,geo,4,ires,pos);
   if(aden<  badden){
      if(verbose>=1)
        printf("Bad nterm end because bad N den = %5.3f at res %d\n",aden,ires);
      ireturn = 1;      
      goto cend;
   }     
   aveden[i] = aveden[i] +aden;
   ca2den[i]=aden= aveatom(den1,geo,1,ires,pos);
   if(aden<  badden){
      if(verbose>=1)
        printf("Bad nterm end because bad Ca2 den = %5.3f at res %d\n",aden,ires);
      ireturn = 1;      
      goto cend;
   } 
#if 0 // New 30may00 - Do not use rations    
//printf("NTERM: ca1den=%5.3f  ca2den = %5.3f  cden/c2den = %5.3f\n",ca1den[i],ca2den[i],ca1den[i]/ca2den[i]);
   if(ca1den[i]/ca2den[i] < 0.5){
      ireturn = 1;      
      goto cend;
   }     
   c2den[i]=aveatom(den1,geo,2,ires,pos);
//printf(" cden=%5.3f  c2den = %5.3f  cden/c2den = %5.3f\n",cden[i],c2den[i],cden[i]/c2den[i]);
   ratio = (cden[i]+ca1den[i])/(c2den[i]+ca2den[i]);
   if(cden[i]/c2den[i] < 0.5){
      ireturn = 1;      
      goto cend;
   }
   if(ratio < 0.55){ // NEW
      ireturn = 1;      
      goto cend;
   }
#endif
   aveden[i] = aveden[i] +aden;
   aveden[i]=aveden[i]/5.0;
   noden[i] = (nden[i]+oden[i])/2.0;

   cend:;
   i=1;
   aveden[i] = 0;
   ires = geo.nend[1]; // c term end
   ca1den[i] = aden= aveatom(den1,geo,1,ires,pos);
   if(aden< 0.0){
      if(verbose>=1)
        printf("Bad cterm end because bad Ca1 den = %5.3f at res %d\n",aden,ires);
      return (ireturn+2); 
   }     
   aveden[i] = aveden[i] +aden;
   nden[1]=aden= aveatom(den1,geo,4,ires,pos);
   if(aden< badden){
      if(verbose>=1)
        printf("Bad cterm end because bad N den = %5.3f at res %d\n",aden,ires);
      return (ireturn+2); 
   }     
   aveden[i] = aveden[i] +aden;
   ires = ires-1;
   cbden[i]=aden= aveatom(den1,geo,5,ires,pos);
   oden[i] = aden = aveatom(den1,geo,3,ires,pos);
   if(aden< badden){
      printf("Bad O den =%5.3f for cterm end res=%d\n",aden,ires);
      return (ireturn+2);
   }     
   aveden[i] = aveden[i] +aden;
   ca2den[i]=aden= aveatom(den1,geo,1,ires,pos);
   if(aden< badden){
      if(verbose>=1)
        printf("Bad cterm end because bad Ca2 den = %5.3f at res %d\n",aden,ires);
      return (ireturn+2); 
   }     
   aveden[i] = aveden[i] +aden;
   cden[i]=aden= aveatom(den1,geo,2,ires,pos);
   if(aden< badden){
      if(verbose>=1)
        printf("Bad cterm end because bad C den = %5.3f at res %d\n",aden,ires);
      return (ireturn+2); 
   }     
//printf("CTERM: ca1den=%5.3f  ca2den = %5.3f  cden/c2den = %5.3f\n",ca1den[i],ca2den[i],ca1den[i]/ca2den[i]);
   if(ca1den[i]/ca2den[i] < 0.5){
      if(verbose>=1)
         printf("Bad cterm end because bad Caden/ca2den = %5.3f at res %d\n",ca1den[i]/ca2den[i],ires);
      return (ireturn+2);
  }
#if 0 // New 30may00 - Do not use rations    
   n2den[1]=aveatom(den1,geo,4,ires,pos);
   ratio = (nden[i]+ca1den[i])/(n2den[i]+ca2den[i]);
//printf(" nden=%5.3f  n2den = %5.3f  nden/n2den = %5.3f\n",nden[i],n2den[i],nden[i]/n2den[i]);
   if(nden[i]/n2den[i] < 0.5)
      return (ireturn+2);      
   if(ratio < 0.55)
      return (ireturn+2);
#endif 
//printf(" cterm Cb den = %5.3f\n",aden); 
   aden = aveatom(den1,geo,5,geo.nend[1]-1,pos);
//printf(" nterm Cb den = %5.3f\n",aden);  
   if(aden < badden){
      if(verbose>=1)
        printf("Bad cterm end because bad Cb den = %5.3f at res %d\n",aden,ires);
      return (ireturn+2); 
   }

   aveden[i] = aveden[i] +aden;
   aveden[i]=aveden[i]/5.0;
   noden[i] = (nden[i]+oden[i])/2.0;
if(verbose>=2){
//if(verbose>=1){
   i=0;
   int itype =  geo.rtype[ geo.nend[0]+1 ];
   printf("For  type %d: end %d: ca1 = %5.3f n=%5.3f c=%5.3f o = %5.3f ca2=%5.3f\n",
                itype,i,ca1den[i],nden[i],cden[i],oden[i],ca2den[i]);
   printf("  ca1/ca2=%5.3f  cb = %5.3f  n+0 = %5.3f  ave = %5.3f\n",
                          ca1den[i]/ca2den[i],cbden[i],noden[i],aveden[i]); 
   i=1;
   printf("For  type %d: end %d: ca1 = %5.3f n=%5.3f c=%5.3f o = %5.3f ca2=%5.3f\n",
                itype,i,ca1den[i],nden[i],cden[i],oden[i],ca2den[i]);
   printf("  ca1/ca2=%5.3f  cb = %5.3f  n+0 = %5.3f  ave = %5.3f\n",
                          ca1den[i]/ca2den[i],cbden[i],noden[i],aveden[i]); 
}
   return ireturn;
}

void checkoneend(struct griddata& den1,struct geometry& geo,int *icut)
/* checks ends of one geo
   returns icut = number of end residues rmoved
*/
{
   int iremove,ii;
   float c4dist;

      iremove = checkendden(den1,geo,geo.pos);      
      if(geo.rtype[ geo.nend[0]+1 ]==1){ // helix
           ii= geo.nend[0]+2;
           c4dist = distv(geo.pos[ii][1],geo.pos[ii-2][3]);
           if(c4dist> 5.75){
               printf("Cutting nterm end of helix because bad c4dist = %5.3f\n",c4dist);
               if( iremove ==0)
                  iremove = 1;
               if(iremove ==2)
                  iremove = 2;
           }
      }
      if(iremove==0)
         *icut = 0;
      else if(  (iremove == 1)||(iremove==2) )
         *icut =1;
      else // cut both ends
         *icut = 2;
}


void checkoneenden(struct griddata& den1,struct geometry& geo,
    int *quittrace)
/* 
  Adapted from checkend - applied to one fit
  Returns quitrace ==2 if bad
  NOTE: Allow geo.savlen to be one less than in generall routine
   This allows to cut so that length = original geo.savlen -1
*/
{
   int ii,iremove;
   float c4dist;

     iremove = checkendden(den1,geo,geo.pos);
     if(geo.rtype[ geo.nend[0]+1 ]==1){ // helix
           ii= geo.nend[0]+2;
           c4dist = distv(geo.pos[ii][1],geo.pos[ii-2][3]);
           if(c4dist> 5.75){
               printf("Cutting nterm end of helix because bad c4dist = %5.3f\n",c4dist);
               if( iremove ==0)
                  iremove = 1;
               if(iremove ==2)
                  iremove = 2;
           }
      }
    cutend(den1,geo,quittrace,iremove);
    if( (iremove!=0)&&(*quittrace!=2) ){ // if end removed and OK, call recursively
       checkoneenden(den1,geo,quittrace);
    }
}

void checkend(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
        int savlen,int startfitnum)
/*
NOTE: This is not used anymore!!
NOTE: THIS IS ONLY CALLED BY TRACE - ASSSUMES THAT FITS ARE SHORT 
   ENOUGH THAT CAN USE copyfittogeo
NOTE: IN FUTURE - THIS WILL BE CALLED IN trace so that does not need to be
   done here. 
checks end residues to see if questionable
NEW: 8mar00 - Add 2 new conditons
  1) c4dist for nterm end of helix.
  2) caden/oden ratio for nterm end of sheet.
NEW 1dec99 - only cut if long enough
Two checks:
  1) Finds den. for end residues - if any atom has bad density -remove
  2) Finds den of Ca(end) if less than 0.5  - remove
  3) Takes ratio of Ca(end)/Ca(end-1).  If less than <0.5 remove
  4) Very special case -  nterm end of helix: if nside = geo.maxside and mainchain
       does not lie on mainbranch (using 1.4 density) - then remove
  checkendden returns 0 if OK, 1 if bad nterm, 2 if bad cterm 3 if both bad
*/
{
   int i,iremove,icut,ii,mmin;
   struct pepfit fit1;
   float c4dist;

   fit1.nalloc = -10; // initial default - indicates have not yet allocated 
   for(i=startfitnum;i<=fitnum[0];i++){ // check end residues and write out new fit file
      copyfittogeo(geo,fit[i]);
      iremove = checkendden(den1,geo,geo.pos);
      if(fit[i].itype[1]==1){ // helix
           ii=2;
           c4dist = distv(fit[i].pos[ii][1],fit[i].pos[ii-2][3]);
           if(c4dist> 5.75){
               printf("Cutting nterm end of helix because bad c4dist = %5.3f\n",c4dist);
               if( iremove ==0)
                  iremove = 1;
               if(iremove ==2)
                  iremove = 2;
           }
      }
      if(iremove==0)
         icut = 0;
      else if(  (iremove == 1)||(iremove==2) )
         icut =1;
      else // cut both ends
         icut = 2;
      if(iremove !=0){
         mmin = savlen -(geo.nend[1]-geo.nend[0] -icut);
         if( mmin>0 )
           printf("   Cutting end reduces length %d below savlen\n",mmin);
      }
      if(iremove ==1){
          printf("Fit %d: Remove nterm end because bad checkend\n",i);
          copyfittofit(fit[i],fit1,1,fit[i].nres); // remove 1 nterm res. 
          copyfittofit(fit1,fit[i],0,fit1.nres);
      }
      else if (iremove==2){
          printf("Fit %d: Remove cterm end because bad checkend\n",i);
          copyfittofit(fit[i],fit1,0,fit[i].nres-1); // remove 1 cterm res. 
          copyfittofit(fit1,fit[i],0,fit1.nres);
      }
       else if (iremove==3){
          printf("Fit %d: Remove nterm and cterm end because bad checkend\n",i);
          copyfittofit(fit[i],fit1,1,fit[i].nres-1); 
          copyfittofit(fit1,fit[i],0,fit1.nres);
      }
    }
    checkendbranch(den1,geo,bone1,startfitnum); //Checks if end on branched main bone
int ifit;
    for(ifit=startfitnum;ifit<=fitnum[0];ifit++){ // remove end residues after checkendbranch 
      if(fit[ifit].itemp[0]==1){ // remove nterm end
         printf("Removing nterm end of fit %d because not on brached main\n",ifit);   
         copyfittofit(fit[ifit],fit1,1,fit[ifit].nres); // remove 1 nterm res. 
         copyfittofit(fit1,fit[ifit],0,fit1.nres);
      }
      if(fit[ifit].itemp[1]==1){ // remove cterm end
          printf("Removing cterm end of fit %d because not on brached main\n",ifit);   
          copyfittofit(fit[ifit],fit1,0,fit[ifit].nres-1); // remove 1 cterm res. 
          copyfittofit(fit1,fit[ifit],0,fit1.nres);
      }
    }
//    savefit();
    quitrout:;
    if(fit1.nalloc!= -10)
        freepepfit(fit1);
}   

void maxdenpos(struct griddata& den1,int cmax[3],
      struct geometry& geo,float pos[3],float radius,float *maxden,int max1d[])
/* Given the positon of atoms center (=pos), tests all grid points with
  radius (=1 Angstrom) of center and returns the value of maximum density
NOTE - assume that only need to check 2 grid steps from center
*/
{
   int i,ir[3],ivec[3],k,ii,jj,kk,pos1d;
   float fd[3],dist,den,fvec[3];

   *maxden = -1000.0;
   xtof(den1,pos,fd); /* convert from double word to float sccreen*/
   for(k=0;k<=2;k++)
       ir[k]=fd[k]+0.5; /*round  to integer*/
   for(ii=ir[0]-2;ii<=ir[0]+2;ii++){
     if( (ii>=1)&& (ii<=den1.amax[0]-den1.amin[0]-1) )
        for(jj=ir[1]-2;jj<=ir[1]+2;jj++)
          if( (jj>=1)&& (jj<=den1.amax[1]-den1.amin[1]-1) )
            for(kk=ir[2]-2;kk<=ir[2]+2;kk++){
              if( (kk>=1)&& (kk<=den1.amax[2]-den1.amin[2]-1) ){
                  ivec[0]=ii;
                  ivec[1]=jj;
                  ivec[2]=kk;
                  imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
                  dist = distv(fvec,pos);
                  if(dist<=radius){
                    den= den1.den[ii][jj][kk];
                    if(den> *maxden){
                       pos1d = d3tod1(cmax,ii,jj,kk);
                       for(i=1;i<=max1d[0];i++){
                          if(pos1d ==max1d[i])
                             goto nextpt;
                       }
                       *maxden = den;// if reach here, maxpt is a new pt
                       max1d[0] = max1d[0]+1;
                       if(max1d[0]>=MAX1D){
                          printf("ERROR: TOO MANY ATOMS IN fit for maxdenpos- quitting\n");
                          goto quitrout;
                       }
                       max1d[ max1d[0] ] = pos1d;
                     }
                  }
               }// kk>=1
               nextpt:;
             }//kk
   }//ii
   quitrout:;
}


void avemax(struct griddata& den1,struct geometry& geo,
       float pos[NAFIT][NAA][3],float *avemaxden)
/*
For each atom find the grid point with highest density within
   distance of radius of atom center.
  Saves 1d pos of max point and does not allow same pt to be used for more
    than one atom
  First finds max pt for Cb2 - and excludes this point for main chain max points
Returns the average of these densities.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int j,count;
   int max1d[MAX1D],cmax[3];
   float radius,maxden;

   *avemaxden = 0.0;
   count=0;
   max1d[0]=0; // counts saved max pts
   radius = 1.5; //check all points with radius of atoms center
   for(j=0;j<=2;j++)
      cmax[j] = den1.amax[j]-den1.amin[j]-1;
#if 1
   for(j=geo.nend[0]+1;j<=geo.nend[1]-1;j++){  // Exclude the points with radiud of cb2
      if(geo.iside[j]>1) // If cb2 exists - add max pt to max1d[] - but do not include in avemadden
         maxdenpos(den1,cmax,geo,pos[j][6],radius,&maxden,max1d);
         if(max1d[0]>=MAX1D){ // indicates reached limit - too many atoms
            *avemaxden = -100.0;
            goto quitrout;
         }
   } 
#endif    
   for(j=geo.nend[0];j<=geo.nend[1];j++){  // Ca
      maxdenpos(den1,cmax,geo,pos[j][1],radius,&maxden,max1d);
      if(max1d[0]>=MAX1D){ // indicates reached limit - too many atoms
            *avemaxden = -100.0;
            goto quitrout;
      }
      if(maxden>0){  // if some point with in radius of center
         *avemaxden = *avemaxden + maxden;
         count = count+1;
      }
   } //Ca
   for(j=geo.nend[0]+1;j<=geo.nend[1];j++){ // N
     maxdenpos(den1,cmax,geo,pos[j][4],radius,&maxden,max1d); 
     if(max1d[0]>=MAX1D){ // indicates reached limit - too many atoms
            *avemaxden = -100.0;
            goto quitrout;
     }
      if(maxden>0){  // if some point with in radius of center
         *avemaxden = *avemaxden + maxden;
         count = count+1;
      }
   }
   for(j=geo.nend[0];j<=geo.nend[1]-1;j++){ 
     maxdenpos(den1,cmax,geo,pos[j][2],radius,&maxden,max1d); 
     if(max1d[0]>=MAX1D){ // indicates reached limit - too many atoms
            *avemaxden = -100.0;
            goto quitrout;
     }
      if(maxden>0){  // if some point with in radius of center
         *avemaxden = *avemaxden + maxden;
         count = count+1;
      }
     maxdenpos(den1,cmax,geo,pos[j][3],radius,&maxden,max1d); 
     if(max1d[0]>=MAX1D){ // indicates reached limit - too many atoms
            *avemaxden = -100.0;
            goto quitrout;
     }
     if(maxden>0){  // if some point with in radius of center
         *avemaxden = *avemaxden + maxden;
         count = count+1;
      }
   }// j
   *avemaxden = *avemaxden/count;
   quitrout:;
}

void avemaxbatch(struct griddata& den1,struct geometry& geo)
/* First calls avemax - finds ave value of ration ave/max
  Then finds cut off value defining bad fits and repeats,
   indicating bad fits
  Use standard cutoff of 0.6 for fah which has averatio of 0.706
*/
{
   int ii,count,nbad,resbad[NAFIT],del[MAXDEL];
   float ratio,finalden,avemaxden,averatio,aveco,avecb2,cutoff;

   count = 0;
   averatio = 0.0;
   if(verbose>=1)
      printf("FINDING BAD FITS BASED ON RATIO OF AVEDEN/AVEMAXDEN\n");
   del[0]= 0; // counts fit that are to be deleted
   for(ii=1;ii<=fitnum[0];ii++){ // first find ave value of ratio for long fits
     if(fit[ii].nres>= geo.savlen+1){ // long fits - presumably the better ones
        copyfittogeo(geo,fit[ii]);     
        aveden(den1,geo,geo.pos,&finalden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
        avemax(den1,geo,geo.pos,&avemaxden);
        if(avemaxden < -1.0)// indicates too many atoms
          goto quitrout;
        if(verbose>=2)
           printf("Fit %d: nres=%d  aveden = %5.3f,,avemaxden =%5.3f  ave/max = %5.3f\n",
               ii,fit[ii].nres,finalden,avemaxden,finalden/avemaxden);
        ratio = finalden/avemaxden;
        averatio = averatio+ratio;
        count = count+1;
      } // lng
   }// ii
   averatio = averatio/count;
   if(verbose>=2)
      printf(" ave of finalden/avemaxden = %5.3f\n",averatio);
   cutoff = 0.6 * averatio/0.706;  // Use standard cut off of 0.6 for fah which has ratio of 0.706
   for(ii=1;ii<=fitnum[0];ii++){ // first find ave value of ratio for long fits
        copyfittogeo(geo,fit[ii]);     
        aveden(den1,geo,geo.pos,&finalden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
        avemax(den1,geo,geo.pos,&avemaxden);
        if(avemaxden < -1.0)// indicates too many atoms
          goto quitrout;
        if(verbose>=2)
           printf("Fit %d: nres=%d  aveden = %5.3f,,avemaxden =%5.3f  ave/max = %5.3f\n",
              ii,fit[ii].nres,finalden,avemaxden,finalden/avemaxden);
        ratio = finalden/avemaxden;
        if( (ratio<cutoff)&&(fit[ii].itype[1]==2) ){ // 8mar01 - only delete if sheet
          if(verbose>=1)
             printf("BAD SHEET FIT %d becasue aveden/maxden = %5.3f < cutoff = %5.3f\n",ii,ratio,cutoff);
          del[0]= del[0]+1; // add to list to be removed
          if( del[0]>= MAXDEL){
             printf("ERROR: avemaxbatch() tried to remove %d fits - assume some error, do not remove any\n",del[0]);
             goto quitrout;
          }
          del[ del[0] ] = ii;
        }
   }// ii
   removenfits(fit,del);  // Removes fits ordered with smallest first -This always saves current fit
   quitrout:;
}
