/* editbones.c++  Routines to edit bones at ends of fits
14jul99: changed so that first point on main bone has mtype = 0
Problem with endo - fit5 den 1.0
    ERROR: jj = 3 > 2  setting = 2
  Try being smarter in which one is picked
PROCEDURE:
23jun99: Prodedure for meating other fits:
   Set side pts and mainpt = -1,
   Set mainpts at ends of fits = ifit.
   If traced bonepoint hits a -1, treat as endpoint and quit
   If traced bonepoint hits an ifit, keep tracing along until 
     no more bone points hit this ifit, and then assume that this is
     the only correct bonetrace, and remove the rest.   
*/
#include "drawwindow.h"
#include "editbones.h"
#include "math.h"
#include "trace.h" // includes a number of "extern .... " lines
#include "peptide.h"
#include "tordyn.h"
#include "extendfit.h"
#include "initialize.h"
#include "bone.h"
#include "checkfit.h"
#include <string.h>

extern struct extbones ext; // All the data for edited bones used to extend one end of a given fit
extern struct maiddatafile maidfile;
extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern int fitnum[NFIT];// see trace.c++ for definition
extern int verbose; //Defined in trace.c++  if = 1, print a lot of stuff
extern FILE *fpout; // = stdout or logfile, depending whether batch or screen version



void startnewmain(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
     struct mainbone& main,int mainnum,int ipos,int npos,int startnum,int origmainnum);
void nextbr(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
    struct mainbone& main,int mainnum,int ipos,int npos,int ncon);
void tracetobr(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       int *ipos,int *npos,int *ncon,int *sidenum,int tpos[],
       int startnum,int *metend,float *aveden);
void tracetobrrec(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       int *ipos,int *npos,int *ncon,int *sidenum,int tpos[],
       int startnum,int *metend,float *aveden1,float *aveden2);
void nearestbonept(struct griddata& den1,float pos[3],struct geometry& geo,struct bonedata& bone,
    int *nearnum,float *neardist);
void nearestbrbonept(struct griddata& den1,float pos[3],struct geometry& geo,struct bonedata& bone,
    int *nearnum,float *neardist);
int extendmainonfitcheck(struct geometry& geo,int endmain);

#if 0  // Never used or debugged
void checkobranch(int checknear,struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,float posc[3],float posca[3],
            float poso[3],int *ibranch)
/*
 Check if branch seems to be going along c-O
Input positions C, Ca and O
Find nearest bonept to C:
The see if there is branch moving in direction of O or Ca
If find a branch:
   Trace out on branch and find bone pts that bracket O and Ca and
     find distance perpendicular distance from O and Ca to this line

*/
{
   int i,j,k,nearnum,ncon,metend;
   int icon,mcon,sidenum,tpos[MAXSIDEPTS],ipos,npos;
   float fvec[3],distc,neardist,aveden;  
   int imino,iminca;
   float mino,minca,minox[3],mincax[3],perpdisto,perpdistca;
   float dist,dist1,dist2,minperpdisto,minperpdistca,fvec1[3],fvec2[3];
   float linedisto,linedistca;

   *ibranch = 0;
   metend = -1;  // This should not be needed here but is needed for tracetobr
   nearestbonept(den1,posc,geo,bone1,&nearnum,&neardist);
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,fvec);// convert pos to world fvec[3]
   distc=distv(posc,fvec);
   if(distc>=2.0){ // If nearest point far away - assume there is not a main branch
       *ibranch= 0;
       goto quitrout;
   }
   ncon = bone1.ptv[nearnum].num; // number of connections
   if(ncon >4){ 
    if(verbose>=2)
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
   minperpdistca = minperpdisto = 100000.0;
   for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
       ipos= nearnum;
       npos= icon;
       tracetobr(den1,geo,bone1, &ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden);// non -recursive version
       if(mcon>2){// then found a branch   find perpendic. dist to O and Ca
          mino = minca = 10000.0;
          for(i=1;i<=sidenum;i++){ // find closest point to O and ca
             d1tow(den1,bone1.cmax,bone1.ptv[ tpos[i] ].pos,fvec);// convert pos to world fvec[3]
             dist = distv(fvec,poso);
             if(dist < mino){
                mino = dist;
                imino = i;
                for(k=0;k<=2;k++)
                  minox[k] = fvec[k];
             }
             dist = distv(fvec,posca);
             if(dist < minca){
                minca = dist;
                iminca = i;
                for(k=0;k<=2;k++)
                  mincax[k] = fvec[k];
             }
          }//i
       } // if mcon
       if(mino < 1.5){ // find perp dist
           d1tow(den1,bone1.cmax,bone1.ptv[imino-1].pos,fvec1);// convert pos to world fvec[3]
           dist1 = distv(fvec1,poso);
           d1tow(den1,bone1.cmax,bone1.ptv[imino+1].pos,fvec2);// convert pos to world fvec[3]
           dist2 = distv(fvec2,poso);
           if(dist1<dist2)
              closedist2(minox,fvec1,poso,&perpdisto,&linedisto);
           else
              closedist2(minox,fvec2,poso,&perpdisto,&linedisto);
        }
        else
           perpdisto = mino;
        if(minca < 1.5){ // find perp dist
           d1tow(den1,bone1.cmax,bone1.ptv[iminca-1].pos,fvec1);// convert pos to world fvec[3]
           dist1 = distv(fvec1,posca);
           d1tow(den1,bone1.cmax,bone1.ptv[iminca+1].pos,fvec2);// convert pos to world fvec[3]
           dist2 = distv(fvec2,posca);
           if(dist1<dist2)
              closedist2(mincax,fvec1,posca,&perpdistca,&linedistca);
           else
              closedist2(mincax,fvec2,poso,&perpdistca,&linedistca);
        }
        else
           perpdistca = minca;
        if(perpdisto < minperpdisto)
           minperpdisto = perpdisto;
        if(perpdistca < minperpdistca)
           minperpdistca = perpdistca;                  
    } // j
    quitrout:;
} 

void checkside_ca_branch(int checknear,struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,float posca1[3],float posca2[3],
           float poscb[3],int *ibranch)
/*
return 1 if want to cut end, 0 otherwise.
Checks if side branch starting at posca2 comes close to posca1
Procedure:  
   Finds nearest branch point to Ca2. If point > 1.5 Ang. for Ca2, then
      quit and return 0
   Then from this branch there should be one branch that move toward Ca1
     and one that move toward Cb. Trace the the one that goes to Cb to
     end or branch and  and checks dist of bone pts from Ca1 - 
          if < minside_cadist (=2.0) Angstromgs - return 1 
*/
{

   int i,j,nearnum,ncon,metend;
   int nposcb,nposca;
   int icon,mcon,sidenum,tpos[MAXSIDEPTS],ipos,npos;
   float fvec[3],neardist,aveden;
   float dist,distca,distca2,distcb,distcb2,delca,delcb,minside_cadist;

   *ibranch = 0;
   minside_cadist= 2.0;
   metend = -1;  // This should not be needed here but is needed for tracetobr
   nearestbrbonept(den1,posca2,geo,bone1,&nearnum,&neardist);
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,fvec);// convert pos to world fvec[3]
   distca=distv(posca2,fvec); // dist ca2 to nearpt
   if(distca>=1.5){ // If nearest point far away - assume there is not a side branch
       *ibranch= 0;
       goto quitrout;
   }
   distca = distv(posca1,fvec); // distance from Ca1 to nearpt
   distcb = distv(poscb,fvec); 
   ncon = bone1.ptv[nearnum].num; // number of connections - should be >2
   if(ncon >4){ 
    if(verbose>=2)
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
   nposca = nposcb = -1;  // indicates did not find appropriate branch
   for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
       distca2 = distv(posca1,fvec); // distance from Ca1 to connect point
       distcb2 = distv(posca1,fvec); // distance from Cb to connect point
       delca = distca-distca2; // change in distance
       delcb = distcb-distcb2;
       if( (delcb>0) && (delcb>delca) ){ // If moving toward Cb more than toward Ca
           nposcb = icon;
       }
       else if( (delca>0) && (delca>delcb) ){ // If moving toward Ca more than toward Cb
           nposca = icon;
       }
    }// j
    if( (nposca>=0) && (nposcb>=0) ){ // Found the correct set of branch
           ipos= nearnum;
           npos= nposcb;
           tracetobr(den1,geo,bone1, &ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden);// non -recursive version
           for(i=1;i<=sidenum;i++){// see if side comes close to Ca
              d1tow(den1,bone1.cmax,bone1.ptv[ tpos[i] ].pos,fvec);// convert pos to world fvec[3]
              dist = distv(posca1,fvec);
              if(dist< minside_cadist){
                 *ibranch = 1;
                 if(verbose>=2)
                    printf("Bad end because sidebranch come close to Ca\n");
                 goto quitrout;
              }
           }// i
    } // if OK branch
    *ibranch = 0;  // default - end OK
    quitrout:;
} 
#endif

void copymainpt(int main0,int pt0,int main1,int pt1)
/* Copies mainpt from main[main1].mainpt[pt1] to main0, pt0
  Does not copy mtype, stype or ctype since they are either not needed
    or will be different
*/
{
   int k;
   ext.main[main0].mainpt[pt0].num = ext.main[main1].mainpt[pt1].num;
   for(k=0;k<=2;k++)
       ext.main[main0].mainpt[pt0].pos[k] =
              ext.main[main1].mainpt[pt1].pos[k];
}

void trace2(struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,int near1,int near2,float pos2[3],
            int *inum,float *neardist)
/*
    Start with near 1 and find all points connected
         going down 2 connections -
      If one of these points = near2, then trace3 return inum = -1;
         else trace returns the inum that is closest to  near2 and
           and the neardist of inum.
*/
{
   int j,j2,ncon,icon,icon2,ncon2,imin;
   float dist,mindist,fvec[3];
     
   *inum = 0; // default
   d1tow(den1,bone1.cmax,bone1.ptv[near1].pos,fvec);// convert pos to world fvec[3]
   mindist= distv(pos2,fvec); // distance from pos2 to initial pt on entry
   ncon = bone1.ptv[near1].num; // number of connections
   for(j=0;j<ncon;j++){
     icon = bone1.ptv[near1].con[j]; // the pt num of the connect point
     if(icon ==near2){ // then connected
        *inum = -1; // indicates connected
        goto quitrout;
     }
     ncon2 = bone1.ptv[icon].num; // number of connections
     d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
     dist = distv(pos2,fvec); // distance from C to connect point
     if(dist< mindist){
        mindist = dist;
        imin = icon;
     }
     for(j2=0;j2<ncon2;j2++){
        icon2 = bone1.ptv[icon].con[j2]; // the pt num of the connect point
        if(icon2 ==near2){ // then connected
          *inum = -1; // indicates connected
           goto quitrout;
        }
        d1tow(den1,bone1.cmax,bone1.ptv[icon2].pos,fvec);// convert pos to world fvec[3]
        dist = distv(pos2,fvec); // distance from C to connect point
        if(dist< mindist){ 
            mindist = dist;
            imin = icon2;
        }        
     }// j2
   }// j
   *inum = imin;
   *neardist = mindist;
   quitrout:;
}

void trace3(struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,int near1,int near2,float pos2[3],
            int *inum,float *neardist)
/*
FOUND BAD BAD ERROR - (NOT REPRODUCABLE!!!)
  If dist > great than  mindist for entire routine -
   retuns  *inum = imin   but imin NEVER SET.
 FIXED:  ON INPUT SET  imin = near1 - default value
   Same as trace2 but goes 3 connections down
*/
{
   int j,j2,j3,ncon,icon,icon2,ncon2,ncon3,icon3,imin;
   float dist,mindist,fvec[3];
     
   *inum = 0; // default
   imin = near1; // Best point on input
   d1tow(den1,bone1.cmax,bone1.ptv[near1].pos,fvec);// convert pos to world fvec[3]
   mindist= distv(pos2,fvec); // distance from pos2 to initial pt on entry
   ncon = bone1.ptv[near1].num; // number of connections
   for(j=0;j<ncon;j++){
     icon = bone1.ptv[near1].con[j]; // the pt num of the connect point
     if(icon ==near2){ // then connected
        *inum = -1; // indicates connected
        goto quitrout;
     }
     ncon2 = bone1.ptv[icon].num; // number of connections
     d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
     dist = distv(pos2,fvec); // distance from C to connect point
     if(dist< mindist){
        mindist = dist;
        imin = icon;
     }
     for(j2=0;j2<ncon2;j2++){
        icon2 = bone1.ptv[icon].con[j2]; // the pt num of the connect point
        if(icon2 ==near2){ // then connected
          *inum = -1; // indicates connected
           goto quitrout;
        }
        ncon3 = bone1.ptv[icon2].num; // number of connections
        d1tow(den1,bone1.cmax,bone1.ptv[icon2].pos,fvec);// convert pos to world fvec[3]
        dist = distv(pos2,fvec); // distance from C to connect point
        if(dist< mindist){ 
            mindist = dist;
            imin = icon2;
        }        
        for(j3=0;j3<ncon3;j3++){
            icon3 = bone1.ptv[icon2].con[j3]; // the pt num of the connect point
            if(icon3 ==near2){ // then connected
              *inum = -1; // indicates connected
              goto quitrout;
            }
            d1tow(den1,bone1.cmax,bone1.ptv[icon3].pos,fvec);// convert pos to world fvec[3]
            dist = distv(pos2,fvec); // distance from C to connect point
            if(dist< mindist){ 
               mindist = dist;
               imin = icon3;
            }        
        }// j3
     }//j2
   }// j
   *inum = imin;
   *neardist = mindist;
   quitrout:;
}


int connected2(struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,int near1,int near2,
            float pos1[3],float pos2[3])
/* 20feb01 - new routine - Starting at point near 1, Call tracetobr twice
     and sees if any points on along fit to branch = near2
  Just returns first connection -
    May not be shortest if there is a loop
*/
{
  int ii,j,j2,j3,ncon1,conpts;
  int icon,mcon,sidenum,tpos[MAXSIDEPTS],ipos,npos,metend;
  int ipos2,npos2,sidenum2,mcon2,tpos2[MAXSIDEPTS];
  int ipos3,npos3,sidenum3,mcon3,tpos3[MAXSIDEPTS];
  float aveden;

   conpts = -1; // default - indicates not connected
   ncon1 = bone1.ptv[near1].num; // number of connections
   metend = -1; //NEW -  7nov01 NOT PRESENT PREVIOUS - SEROUS BUG IN ALPHA MACHINE
   for(j=0;j<ncon1;j++){
       icon = bone1.ptv[near1].con[j]; // the pt num of the connect point
       ipos = near1;
       npos = icon;
       tracetobr(den1,geo,bone1,&ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden);// non -recursive version -NOTE RETURN NEW  ipos, npos
       for(ii=1;ii<=sidenum;ii++){
          if(tpos[ii]==near2){
             conpts = ii;
             return conpts;
          }
       }
       if(mcon>1){// then call again
         for(j2=0;j2<mcon;j2++){
           ipos2 = tpos[sidenum];
           npos2 = bone1.ptv[ipos2 ].con[j2]; // the pt num of the connect point
           if(npos2!= tpos[ sidenum-1 ]){ // do not go back
               tracetobr(den1,geo,bone1,&ipos2,&npos2,&mcon2,&sidenum2,tpos2,0,
                  &metend,&aveden);// non -recursive version
               for(ii=1;ii<=sidenum2;ii++){
                   if(tpos2[ii]==near2){
                     conpts =sidenum+ ii;
                     return conpts;
                   }
               }
               if(mcon2>1){// then call again
                  for(j3=0;j3<mcon2;j3++){
                   ipos3 = tpos2[sidenum2];
                   npos3 = bone1.ptv[ipos3 ].con[j3]; // the pt num of the connect point
                   if(npos3!= tpos2[ sidenum2-1 ]){ // do not go back
                     tracetobr(den1,geo,bone1,&ipos3,&npos3,&mcon3,&sidenum3,tpos3,0,
                       &metend,&aveden);// non -recursive version
                     for(ii=1;ii<=sidenum3;ii++){
                        if(tpos3[ii]==near2){
                          conpts =sidenum+ sidenum2+ ii;
                          return conpts;
                        }
                     }
                   }// npos3
                 }//j3
               }// mcon2     
            }// if npos2
         }//j2      
       }// if mcon
    }// j
    return conpts;
}              

int connected(struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,int near1,int near2,
            float pos1[3],float pos2[3])
/* checks if point near1 and nearnum (which are assumed to be fairly close)
  are connected.
  Return 1 if connected, 0 if not connected
Procedure:
  Starting with near1, 
    1)Call trace2(near1,near2,..), which finds all points connected
         going down 2 connections -
      If one of these points = near2, then trace3 return inum = -1;
         else trace returns the inum that is closest to  near2 and
           and the neardist of inum.
    2) If inum!= -1 (not connected)
         if(inum = near1, then assume points not connected
         else - call trace2(inum,near2)
            and repeat.
    3) Finally - call trace2 a third time  or fourth- this is 2*3 = 6 steps through
        bones and should find a connection if there is one.
*/   
{ 
     int i,icall,initnum,inum;
     float neardist;

//     icall = 4; // number of times call trace2
     icall = 3; // number of times call trace3
     inum = near1; // for first call trace2
     for(i=1;i<=icall;i++){
//printf(" i= %d in connected()\n",i);
        initnum = inum;
//        trace2(den1,geo,bone1,inum,near2,pos2,&inum,&neardist);
        trace3(den1,geo,bone1,inum,near2,pos2,&inum,&neardist);
//printf(" inum = %d near2 = %d neardist = %5.3f\n",inum,near2,neardist);
        if(inum == -1) // connected
           return 1;
        else if(inum == initnum) // if did not get any closer after 2 steps, assume not connected
           return 0;
     } // each time reach here - call trace2 with new inum 
     return 0; // Should not reach here       
}

void checkifbranchrec(int checknear,struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,float pos1[3],float pos2[3],int *ibranch,int pt[])
/*
NEW ADDITION - 6nov00 - Check if point closest to pos1 and pos2 are connected
     if not connected - return ibranch = -1 - indicates bad end 
Return either 0, no branch, or 1 = branch or  2 = branch on branch
New version  called by checkoneendbr - after 5oct00 - uses recursive form of tracetobr
  Also traces to branch even if nearnum is a branch itself.
New: aug2500: if all points connected to nearpt are moving closer
   to pos1 - this is indication of something wrong - return 3
If checknear = 1 then
 If nearnum is a branch pt or is connected to one, return 2
 If there is a branch pt with 2 pts of nearnum and this 
    pt is <1.5 Ang. from pos2, then return 2
  Check if branch along bones in direction from pos1 to pos2.
  Used to see if end residue is good or along side chain.
 Procedure:
    Find nearest point to pos2 = pt.  Check all connections to
        pt = ptcon. If dist from connected ptcon to pos1 is greater 
        than dist from pt to pos, then trace out and see if find branch.
  Return ibranch = 0 if no branch,,(NOT ANH MOREibranch = 2 if there is a)
     a branch point within 1 of pos2, and ibranch = 1 if connected to
     branch in direction pos1 to pos2
*/
{

   int j,nearnum,nearnum1,ncon,metend,imin,near1;
   int ncon2;
   int icon,mcon,sidenum,tpos[MAXSIDEPTS],ipos,npos;
   int maxsidenum,maxmcon;
   float nvec[3],fvec[3],dist,dist2,neardist,aveden,aveden2;
   float mindist;
   int nextpt;

   *ibranch = 0;
   metend = -1;  // This should not be needed here but is needed for tracetobr
   nearestbonept(den1,pos2,geo,bone1,&nearnum,&neardist);  // bone pt nearest to end Ca
   nearnum1 = nearnum;
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,fvec);// convert pos to world fvec[3]
   dist=distv(pos2,fvec);
   if(dist>=2.0){ // If nearest point far away - assume there is not a main branch
       *ibranch= 0;
       goto quitrout;
   }
   ncon = bone1.ptv[nearnum].num; // number of connections
#if 0
   if(ncon>=3){ // if near branch, return 2
       *ibranch = 2;
       goto quitrout;
   }
#endif
   nearestbonept(den1,pos1,geo,bone1,&near1,&neardist);
#if 1
// if(connected(den1,geo,bone1,near1,nearnum,pos1,pos2)==0){ // check if these 2 points are connected
 if(connected2(den1,geo,bone1,near1,nearnum,pos1,pos2)== -1){ // check if these 2 points are connected
    *ibranch= 0; // incidates the 2 pts are not connected, bad end
    if(verbose>=1)
       printf(" Bad end because Ca end pt not connected to second atoms\n");
     goto quitrout;
 }    
#endif
 if(checknear==1){
   if(ncon<=2){ // If not branch - see if there is close branch point conneected to nearnum
     for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       if(ncon2>2){
          nearnum = icon;
          ncon = ncon2;
          goto testnear;
       }
     }
   }// ncon
 }// checknear
   testnear:;
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,nvec);// convert pos to world fvec[3]
   if(ncon >4){ 
    if(verbose>=2)
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
//   mindist = 1000.0; // find the point closest to pos1
   imin= -1; // default - nearnum is closest to pos1
   mindist = distv(pos1,nvec);
   for(j=0;j<ncon;j++){// Do not use the branch that is closest to pos1
           icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
           d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
           dist2 = distv(pos1,fvec); // distance from C to connect point
           if(dist2< mindist){ // if connect point moving away from C
              mindist = dist2;
              imin = j;
           }
   }// j
   maxsidenum= maxmcon = -1;
   if(imin== -1){ // nearnum closest to pos1 (= N) then just check path from nearnum to earnum1
           ipos= nearnum;
           npos= nearnum1;
           tracetobrrec(den1,geo,bone1, &ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden,&aveden2);// recursive version (if second branch, then aveden2>0)
           if(mcon>maxmcon){
              nextpt = icon;
              maxmcon = mcon;
           }
           if(sidenum>maxsidenum){
              maxsidenum = sidenum;
              nextpt = icon;
           }
   }
   else{ // check all paths moving away from pos1
    for(j=0;j<ncon;j++){
     if(j!=imin){
           ipos= nearnum;
           icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
           npos= icon;
           tracetobrrec(den1,geo,bone1, &ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden,&aveden2);// recursive version (if second branch, then aveden2>0)
           if(mcon>maxmcon){
              nextpt = icon;
              maxmcon = mcon;
           }
           if(sidenum>maxsidenum){
              maxsidenum = sidenum;
              nextpt = icon;
           }
       }
     } // j
    }//else
    if((maxmcon>2)||(maxsidenum>ext.sidelim) ){ // Then branch along this chain
       if(verbose>=2)
          printf(" Branch on chain - mcon = %d sidenum = %d\n",mcon,sidenum);
       *ibranch = 1;
       pt[0] = nearnum;
       pt[1] = nextpt;
       goto quitrout;
    }
    quitrout:;
} 

void checkifbranch(int checknear,struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,float pos1[3],float pos2[3],int *ibranch)
/*
NEW  18oct00:  if not on a branch (old - returned 0) - now return ibranch = -sidenum
Does not use recursive form of tracetobr
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
New: aug2500: if all points connected to nearpt are moving closer
   to pos1 - this is indication of something wrong - return 3
If checknear = 1 then
 If nearnum is a branch pt or is connected to one, return 2
 If there is a branch pt with 2 pts of nearnum and this 
    pt is <1.5 Ang. from pos2, then return 2
  Check if branch along bones in direction from pos1 to pos2.
  Used to see if end residue is good or along side chain.
 Procedure:
    Find nearest point to pos2 = pt.  Check all connections to
        pt = ptcon. If dist from connected ptcon to pos1 is greater 
        than dist from pt to pos, then trace out and see if find branch.
  Return ibranch = 0 if no branch,,ibranch = 2 if there is a
     a branch point within 1 of pos2, and ibranch = 1 if connected to
     branch in direction pos1 to pos2
*/
{

   int j,nearnum,ncon,metend;
   int icon2,ncon2,ncon3,j2;
   int icon,mcon,sidenum,tpos[MAXSIDEPTS],ipos,npos;
   int maxsidenum,maxmcon;
   float nvec[3],fvec[3],dist,dist1,dist2,neardist,aveden;

   *ibranch = 0;
   sidenum = 0;
   metend = -1;  // This should not be needed here but is needed for tracetobr
   nearestbonept(den1,pos2,geo,bone1,&nearnum,&neardist);
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,fvec);// convert pos to world fvec[3]
   dist=distv(pos2,fvec);
   if(dist>=2.0){ // If nearest point far away - assume there is not a main branch
       *ibranch= 0;
       goto quitrout;
   }
   ncon = bone1.ptv[nearnum].num; // number of connections
   if(ncon>=3){ // if near branch, return 2
       *ibranch = 2;
       goto quitrout;
   }
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,nvec);// convert pos to world fvec[3]
   dist1 = distv(pos1,nvec); // distance from pt to pos1       
   if(ncon >4){ 
    if(verbose>=2)
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
   maxsidenum= maxmcon = -1;
   for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
       dist2 = distv(pos1,fvec); // distance from C to connect point
       if(dist2>dist1){ // if connect point moving away from C
           ipos= nearnum;
           npos= icon;
           tracetobr(den1,geo,bone1, &ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden);// non -recursive version
           if(mcon>maxmcon)
              maxmcon = mcon;
           if(sidenum>maxsidenum)
              maxsidenum = sidenum;
       }
    } // j
    if( (ncon>1)&&(dist<1.0)&&(maxmcon == -1) ){// all points moving closer to pos1 - something wrong
          *ibranch = 3;
          goto quitrout;
    }
    if((maxmcon>2)||(maxsidenum>ext.sidelim) ){ // Then branch along this chain
       if(verbose>=2)
          printf(" Branch on chain - mcon = %d sidenum = %d\n",mcon,sidenum);
       *ibranch = 1;
       goto quitrout;
    }
 if(checknear==1){
   if(ncon<=2){ // If not branch - see if there is close branch point conneected to nearnum
     for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       if(ncon2>2){
          d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
          dist2=distv(pos2,fvec);
          dist1=distv(pos1,fvec);       
//             if (dist2<= 2.0) { //If not too close to pos 1 then assume a brach at end
         if( (dist2<= 1.7)&&(dist1>1.0) ){ //If not too close to pos 1 then assume a brach at end
             *ibranch = 2;
             goto quitrout;
          }
       }
     }
   }// ncon
   if(ncon<=2){ // If not branch - see if there is close  branch point within 2 poionts
     for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       for(j2=0;j2<ncon2;j2++){
          icon2 = bone1.ptv[icon].con[j2]; // the pt num of the connect point
          ncon3 = bone1.ptv[icon2].num; // number of connections
          if(ncon3>2){
             d1tow(den1,bone1.cmax,bone1.ptv[icon2].pos,fvec);// convert pos to world fvec[3]
             dist2=distv(pos2,fvec);
             dist1=distv(pos1,fvec);       
//             if (dist2<= 2.0) { //If not too close to pos 1 then assume a brach at end
             if( (dist2<= 1.7)&&(dist1>1.0) ){ //If not too close to pos 1 then assume a brach at end
               *ibranch = 2;
                goto quitrout;
             }
          }// ncon2
        }//j2
      } // j
   } // ncon 
 }// checknear 
    quitrout:;
    if(*ibranch ==0) // new return -sidenum instead of 0
        *ibranch = -sidenum;
} 


void checksidebranch(int checknear,struct griddata& den1,struct geometry& geo,
            struct bonedata& bone1,float pos1[3],float pos2[3],int *ibranch)
/*
19feb01 - identical to checkifbranch but checks for loops
  Also if closest bonept to pos2  = CB2 is < 1 Angstrgom from pos1 = C, then assume no branch.
  Also if nearest point to Cb2 is a branch pt, then use this as starting point
     looking for another branch (previosly - just quit and return a branch)  
NEW  18oct00:  if not on a branch (old - returned 0) - now return ibranch = -sidenum
Does not use recursive form of tracetobr
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
New: aug2500: if all points connected to nearpt are moving closer
   to pos1 - this is indication of something wrong - return 3
If checknear = 1 then
 If nearnum is a branch pt or is connected to one, return 2
 If there is a branch pt with 2 pts of nearnum and this 
    pt is <1.5 Ang. from pos2, then return 2
  Check if branch along bones in direction from pos1 to pos2.
  Used to see if end residue is good or along side chain.
 Procedure:
    Find nearest point to pos2 = pt.  Check all connections to
        pt = ptcon. If dist from connected ptcon to pos1 is greater 
        than dist from pt to pos, then trace out and see if find branch.
  Return ibranch = 0 if no branch,,ibranch = 2 if there is a
     a branch point within 1 of pos2, and ibranch = 1 if connected to
     branch in direction pos1 to pos2
*/
{

   int ii,j,nearnum,ncon,metend;
   int icon2,ncon2,ncon3,j2;
   int icon,mcon,sidenum,tpos[MAXSIDEPTS],ipos,npos;
   int maxsidenum,maxmcon;
   float nvec[3],fvec[3],fvec2[3],dist,dist1,dist2,linedist,dist4,neardist,aveden;

   *ibranch = 0;
   sidenum = 0;
   metend = -1;  // This should not be needed here but is needed for tracetobr
   nearestbonept(den1,pos2,geo,bone1,&nearnum,&neardist);
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,fvec);// convert pos to world fvec[3]
   dist=distv(pos2,fvec);// distance for Cb2 to near pt
   if(dist>=2.0){ // If nearest point far away - assume there is not a main branch
       *ibranch= 0;
       goto quitrout;
   }
   dist=distv(pos1,fvec);//distance from near pt to C
   if(dist< 1.0){// if very close to C assum no branch 
       *ibranch= 0;
       goto quitrout;
   }
   ncon = bone1.ptv[nearnum].num; // number of connections
#if 0 //19feb01 - if Cb2 near branch then can be just a normal residue branch - look for another branch
   if(ncon>=3){ // if near branch, return 2
       *ibranch = 2;
       goto quitrout;
   }
#endif
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,nvec);// convert pos to world fvec[3]
   dist1 = distv(pos1,nvec); // distance from pt to pos1       
   if(ncon >4){ 
    if(verbose>=2)
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
   maxsidenum= maxmcon = -1;
   for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
       dist2 = distv(pos1,fvec); // distance from C to connect point
       if(dist2>dist1){ // if connect point moving away from C
           ipos= nearnum;
           npos= icon;
           tracetobr(den1,geo,bone1, &ipos,&npos,&mcon,&sidenum,tpos,0,
                  &metend,&aveden);// non -recursive version
           if(sidenum> MAXSIDEPTS-3){ // Assume this indicates a continous loop  - not a bracnh
              sidenum=mcon=1;
           }  
           else if( (mcon>2)&&(sidenum>2) ){ //NEW 7may01 only check if sidenum>2 - if reached branch pt, check if just looped back to main -linedist = along bones to branch, deldist = increase in dist of branch versus icon from C
              linedist=0.0;
              for(ii=2;ii<=sidenum;ii++){              
                 d1tow(den1,bone1.cmax,bone1.ptv[ tpos[ii] ].pos,fvec);// convert pos to world fvec[3]
                 d1tow(den1,bone1.cmax,bone1.ptv[ tpos[ii-1] ].pos,fvec2);// convert pos to world fvec[3]
                 linedist = linedist + distv(fvec,fvec2); // distance from C to connect point
              }
              d1tow(den1,bone1.cmax,bone1.ptv[ tpos[sidenum] ].pos,fvec);// tpos[sidenum] = branch pt
              dist4 = distv(pos1,fvec); // distance from C to connect point
//  printf(" icon = %d tpos[1] = %d  sidenum = %d linedist = %5.3f deldist = %5.3f\n",
//             icon,tpos[1],sidenum,linedist,dist4-dist2); 
             if(dist4-dist2 < 0.3*linedist){
                 if(verbose>=1)
                    printf("Side branch looped back to main - assuming not a branch -sidenum=%d  deldist = %5.3f linedist = %5.3f\n",
                                 sidenum,dist4-dist2,linedist);
                  mcon = sidenum = 1; // set mcon = 1 as if not a branch
             } 
           }// mcon>2
           if(mcon>maxmcon)
              maxmcon = mcon;
           if(sidenum>maxsidenum)
              maxsidenum = sidenum;
       }
    } // j
    if( (ncon>1)&&(dist<1.0)&&(maxmcon == -1) ){// all points moving closer to pos1 - something wrong
          *ibranch = 3;
          goto quitrout;
    }
    if((maxmcon>2)||(maxsidenum>ext.sidelim) ){ // Then branch along this chain
       if(verbose>=2)
          printf(" Branch on chain - mcon = %d sidenum = %d\n",mcon,sidenum);
       *ibranch = 1;
       goto quitrout;
    }
 if(checknear==1){
   if(ncon<=2){ // If not branch - see if there is close branch point conneected to nearnum
     for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       if(ncon2>2){
          d1tow(den1,bone1.cmax,bone1.ptv[icon].pos,fvec);// convert pos to world fvec[3]
          dist2=distv(pos2,fvec);
          dist1=distv(pos1,fvec);       
//             if (dist2<= 2.0) { //If not too close to pos 1 then assume a brach at end
         if( (dist2<= 1.7)&&(dist1>1.0) ){ //If not too close to pos 1 then assume a brach at end
             *ibranch = 2;
             goto quitrout;
          }
       }
     }
   }// ncon
   if(ncon<=2){ // If not branch - see if there is close  branch point within 2 poionts
     for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       for(j2=0;j2<ncon2;j2++){
          icon2 = bone1.ptv[icon].con[j2]; // the pt num of the connect point
          ncon3 = bone1.ptv[icon2].num; // number of connections
          if(ncon3>2){
             d1tow(den1,bone1.cmax,bone1.ptv[icon2].pos,fvec);// convert pos to world fvec[3]
             dist2=distv(pos2,fvec);
             dist1=distv(pos1,fvec);       
//             if (dist2<= 2.0) { //If not too close to pos 1 then assume a brach at end
             if( (dist2<= 1.7)&&(dist1>1.0) ){ //If not too close to pos 1 then assume a brach at end
               *ibranch = 2;
                goto quitrout;
             }
          }// ncon2
        }//j2
      } // j
   } // ncon 
 }// checknear 
    quitrout:;
    if(*ibranch ==0) // new return -sidenum instead of 0
        *ibranch = -sidenum;
} 

void nearestbonept(struct griddata& den1,float pos[3],struct geometry& geo,struct bonedata& bone,
    int *nearnum,float *neardist)
/* finds the nearest bone pt to pos)
*/
{
    int i,k,ivec[3],mini,totnum;
    float fvec[3],dist,mindist,dvec[3];

    mindist = 10000.0;
    totnum = bone.ptv.end()-bone.ptv.begin();
    for(i=0;i<totnum;++i){
     d1tod3(bone.cmax,bone.ptv[i].pos,ivec);
     imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
     for(k=0;k<=2;k++)
        dvec[k]= fvec[k]-pos[k];
     dist = dot(dvec,dvec);
     if(dist<mindist){
        mindist = dist;
        mini = i;
     }
   }
   if(verbose>=2)
     printf(" Minimum dist to nearest bone pt = %5.3f for pt # %d\n",mindist,mini);
   *nearnum = mini;
   *neardist = sqrt(mindist);
    d1tod3(bone.cmax,bone.ptv[mini].pos,ivec);
    imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
    for(k=0;k<=2;k++)// world co0rd of near pos
      ext.stpos[k]= fvec[k];
}

void nearestbrbonept(struct griddata& den1,float pos[3],struct geometry& geo,struct bonedata& bone,
    int *nearnum,float *neardist)
/* finds the nearest bone pt to pos with ncon>2
*/
{
    int i,k,ivec[3],mini,totnum,ncon;
    float fvec[3],dist,mindist,dvec[3];

    mindist = 10000.0;
    totnum = bone.ptv.end()-bone.ptv.begin();
    for(i=0;i<totnum;++i){
     ncon = bone.ptv[i].num; // number of connections
     if(ncon>2){
       d1tod3(bone.cmax,bone.ptv[i].pos,ivec);
       imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
       for(k=0;k<=2;k++)
          dvec[k]= fvec[k]-pos[k];
       dist = dot(dvec,dvec);
       if(dist<mindist){
         mindist = dist;
         mini = i;
       }
     }// ncon
   }
   if(verbose>=2)
     printf(" Minimum dist to nearest branch bone pt = %5.3f for pt # %d\n",mindist,mini);
   *nearnum = mini;
   *neardist = sqrt(mindist);
}

void addmainpt(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       struct mainbone& main,int mainnum,int ipos,int mtype,int stype)
/* Assign values to ext.mainbone[ext.nmain]
   mtype = main chain number if new mainchain starts with this point, else = 0;
   stype = 0 or side branch number if side starts at this point
*/
{
   int i,k,ivec[3];
   float vec[3];

   i=main.num;
   main.mainpt[i].num =ipos;
   d1tod3(bone1.cmax,bone1.ptv[ipos].pos,ivec);
   imcv(den1,ivec,vec);// convert int grid ivec to float world co-ord fvec
   geo.set[ivec[0]][ivec[1]][ivec[2]] = -1; // set this point
   for(k=0;k<=2;k++)
      main.mainpt[i].pos[k]=vec[k];
   main.mainpt[i].mtype = mtype;
   main.mainpt[i].stype = stype;
   main.num = main.num+1; // Add one point on main - NOTE: points run from 0,1,2,3 ...
   main.mainpt[i].ctype[0] = 0; //Default - no connection to other mains
}

void addbrtomain(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       struct mainbone& main,int mainnum,int tpos[],int sidenum,int ncon)
/* Add point up to brach to main and continue with new branch
  NOTE: set mtype = stype = 0, this is changed later
*/
{
   int j,ivec[3];
   float origdist,vec[3];

 if(verbose>=2)
   printf("In addbrtomain:  main.num = %d\n",main.num);
#if 0
   if(ext.metend>0)
       main.metend = ext.metend;
#endif
   for(j=1;j<=sidenum;j++){ // add  points  to initial main
      if (main.num >= MAXMAINPTS-1){
        if(verbose>=1)
         printf("REACHED LIMIT in addbrtomain: main.num = %d >= MAXMAINPTS-1\n",main.num);
         *reachedlim = 10; // indicates reacehd limit in addbrtomain
         goto quitrout;
      }
      addmainpt(den1,geo,bone1,main,mainnum,tpos[j],0,0);// NOTE: set mtype = stype = 0 for branch point - this is changed later
  }
   d1tod3(bone1.cmax,bone1.ptv[ tpos[sidenum] ].pos,ivec);
   imcv(den1,ivec,vec);// convert int grid ivec to float world co-ord fvec
   origdist = distv(vec,ext.stpos);
   if(ncon>1){
       if (origdist>=ext.distlim){
        if(verbose>=2)
          printf("Quitting this branch because reached DISTANCE LIMIT\n");
       }
       else{
        if(verbose>=2)
          printf("Starting nextbr in addbromain: dist from origin =%d\n",
          main.num+main.startnum);
          nextbr(reachedlim,den1,geo,bone1,main,mainnum,tpos[sidenum-1],
          tpos[sidenum],ncon);
       }
  }
  
  quitrout:;
}

void addsidept(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       struct mainbone& main,int mainnum,int ipos)
/* Assign values to main.side
*/
{
   int i,k,num,ivec[3];
   float vec[3];

   i=main.nside; 
   num = main.side[i].num;
   main.side[i].ipos[num] =ipos;
   main.side[i].ctype[0]=0;  //Default, side does not connect to other main
   d1tod3(bone1.cmax,bone1.ptv[ipos].pos,ivec);
   imcv(den1,ivec,vec);// convert int grid ivec to float world co-ord fvec
   geo.set[ivec[0]][ivec[1]][ivec[2]] = -1; // set this point
   for(k=0;k<=2;k++)
      main.side[i].pos[num][k]=vec[k];
   main.side[i].num =main.side[i].num+1; // add one point
}

void addsidechain(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       struct mainbone& main,int mainnum,int tpos[],int sidenum,int origmainnum)
/* create a new side chain
*/
{
   int j;
   if(sidenum>0){
     main.nside = main.nside+1;
    if(verbose>=2)
     printf("In add side chain nside = %d\n",main.nside);
     if(main.nside>= MAXSIDE){
       if(verbose>=2)
         printf("REACHED LIMIT: main.side = %d\n",main.nside);
         *reachedlim = 5;
         goto quitrout;
     }
     main.side[main.nside].num = 0;// initial number of points = 0 point = bracnh
     main.side[main.nside].mnum = origmainnum; // main number side starts at
     for(j=0;j<=sidenum;j++){ // First side point is branch point = tpos[0]
        addsidept(den1,geo,bone1,main,mainnum,tpos[j]);// first point is branch point
     }
     main.mainpt[origmainnum].stype = main.nside; // Set stype for main pt
  }
  quitrout:;
}

void addpt(struct bonedata& bone1,int icon,int *sidenum,int tpos[],
         int *ipos,int *npos,int *ncon)
{
    *sidenum = *sidenum+1;
    if(*sidenum > MAXSIDEPTS-3){
          printf("ERROR: num. pts in tracetobr > MAXSIDEPTS-3  Setting ncon = 3\n");
          *ncon = 3;
          tpos[*sidenum] =icon;
          *ipos = *npos;
          *npos = icon;
          goto quitrout;
    }
    tpos[*sidenum] =icon;
    *ipos = *npos;
    *npos = icon;
    *ncon = bone1.ptv[*npos].num; // number of connections
    quitrout:;
#if 0
    if(verbose>=2)
      printf("In addpt: ipos = %d npos = %d ncon = %d\n",*ipos,*npos,*ncon);
#endif
}

void tracetobrrec(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       int *ipos,int *npos,int *ncon,int *sidenum,int tpos[],
       int startnum,int *metend,float *aveden1,float *aveden2)
/* 
 This version of tracetobrrec - calls the original tracetobr
    for each branch to find next branch with highest density
 aveden2 <0 if there is not another branch after first branch
   else aveden2 = density on second branch.
Traces out simple points until either reach branch or end
  NEW NEW aveden is highest ave den on next branch 
  New returns aveden which is either ave on minimum den in next 2 0r 3 pts.
  1) ipos is previous point
  2)  npos is next point
 Returns:
  3) ncon of last point = either 1 or >2 
  4) sidenum = number of side points up to end
  5) tpos[1,2,...  upto to sidelim] the pos of each point
  6 and 7) ipos and npos are the last points before exiting
  8) If hit end of another fit, set ext.metend = to geo.set[] = fitnum
       and then only continue along points along this fit
*/
{
  int icon,j,ivec[3],iset,icapos,nden;
  float fvec[3],dist,fden;
  float min1,minden2,min2ave;

   tpos[0]= *ipos;
//   *aveden = 0; // if return aveden
//   *aveden1 = *aveden2 =  1000.0; // now is minimum den
   if(verbose>=2)
     printf("Enter tracetobr_rec:  ipos = %d  npos = %d\n",*ipos,*npos);
   *sidenum =1;
   tpos[*sidenum] = *npos; // first point is branch point 
   *ncon = bone1.ptv[*npos].num; // number of connections
   *aveden1 = 0.0;
   *aveden2 = -10000.0;
   nden =0;
   if( *ncon==1)
     *aveden1 = -1000.0;
   minden2 = 10000.0;
   while(*ncon==2){
      for(j=0;j< *ncon;j++){
            icon = bone1.ptv[*npos].con[j]; // the pt num of the connect point
            if(icon!= *ipos)// do not go backward
              goto foundfor;
      }// j
      foundfor:;
      d1tod3(bone1.cmax,bone1.ptv[icon].pos,ivec);
#if 1 // NEW retrun minimum of two consequtive points
       nden = nden +1;     
       fden = (float)den1.den[ivec[0]][ivec[1]][ivec[2]];
       if(nden>1){
          min2ave = (min1+fden)/2.0;
          min1 = fden;
          if(min2ave< minden2){
            *aveden1 = min2ave;
            minden2 = min2ave;
          }
       }
       else{
           *aveden1 = fden;
           min1=fden;
       }
#endif
#if 0 //Used up to 2sep00: NEW aveden1 is is aveden for next 3 pts.
      if(*sidenum<=3){ // only ave den for first 2 or 3 pts
          nden = nden +1;     
          fden = (float)den1.den[ivec[0]][ivec[1]][ivec[2]];
          *aveden1 = *aveden1 + fden;
      }
      *aveden1 = *aveden1/nden;
#endif
      iset = geo.set[ivec[0]][ivec[1]][ivec[2]];
      if(geo.useset==0) //This is only used for checkifbranch()
         iset = -10;
      if(iset!= -1) { //side and central main pts set = -1 if not previous set
       if(iset >0){ // hit an end of a fit num = geo.set
        if(iset ==ext.ifit){ //  the same fit as started at - check if just begnning
    if(verbose>=2)
     printf(" ISET = IFIT, SIDENUM = %d STARTNUM = %d\n",*sidenum,startnum);
             if(*sidenum + startnum < 5){ // if just beginning
                   addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
                   goto repeat;
             }
          }
          else{ // if(iset !=ext.ifit)  Not the same fit as started at
                 imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
                 if(geo.fixend ==2) // if extending from Cterm, must hit nterm end of fit
                    icapos = 0;
                 else // if extending from nterm, must hit cterm end (Ca) of fit
                    icapos= fit[iset].nres -1;
                 dist = distv(fvec,fit[iset].pos[icapos][1]); // Dist from cterm Ca to icon
    if(verbose>=2)
      printf(" In tracetobrrec dist = %f  icapos =%d icon = %d\n",dist,icapos,icon);
                 if(dist<=2.3){ //If close enough to end Ca then has made correct contact
//printf("at dist<2  ivec[0] = %d\n",ivec[0]);
//                      ext.metend = *metend = geo.set[ivec[0]][ivec[1]][ivec[2]]; // the fit number it hit
                      *metend = geo.set[ivec[0]][ivec[1]][ivec[2]]; // the fit number it hit
                      addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
                      *ncon = 1; //stop  Treat as simple end point so leaves loop
                      goto repeat;                     
                 }
            } //iset !=ext.ifit)
        } //iset>0
        else{ // set != -1 and <0 - i.e. = -10
               addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
               goto repeat;
        }
      } //iset != -1
      else{// iset == -1  addpt and quit
          addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
          *ncon = -1; // Indicates met a previously fit point - force to be side
          goto quitrout;
      }
     *ncon = 1; //If reach here - stop  Treat as simple end point
      repeat:;
  } // while
int jj,npos2,ipos2,metend2,sidenum2[4];
int mcon[4],tpos2[4][MAXSIDEPTS],iposm[4],nposm[4];
float avedenr,minden; // aveden for brach after nextbr
  if(*metend>0)
      *ncon = 3;
  else if(*ncon==3){ // Then find density for next branch 
    jj=0; // counts forward branches
    metend2=0;
    npos2 = *npos;
    ipos2 = *ipos;
    minden = -10000.0;
    for(j=0;j< *ncon;j++){
       icon = bone1.ptv[npos2].con[j]; // the pt num of the connect point
       if(icon!= ipos2){ // go forward
           d1tod3(bone1.cmax,bone1.ptv[icon].pos,ivec);
           iset = geo.set[ivec[0]][ivec[1]][ivec[2]];
           if(iset == -1){ // Then quit this trace
              avedenr = - 100.0; // give it a bad density
           } 
           else{ 
              iposm[jj]= npos2;
              nposm[jj]= icon;
              metend2=  0;
              tracetobr(den1,geo,bone1, &iposm[jj],&nposm[jj],&mcon[jj],
                 &sidenum2[jj],tpos2[jj],startnum,
                 &metend2,&avedenr);  // The old, non-recursive, version
              if(metend2>0){// 8mar01 - if meets another end - make this branch hav e best density
                *aveden1 = 10000.0;
                goto quitrout;
              }
           }          
           if( (mcon[jj]>1)&&(avedenr> minden) ){//if reach another branch - return this best den.
              *aveden2 = avedenr;
              minden = avedenr;
           }
           jj = jj+1;
       } // icon !=ipos
     } // j
  } // ncon = 3
  else{// ncon ==1) // give aveden a bad density
    *aveden2 = -1000.0;
  }
  if(verbose>=2)
    printf("sidenum = %d\n",*sidenum);
  quitrout:;
#if 0 // use this if aveden is actual average
  if(nden >=1)
    *aveden = *aveden/nden;
#endif
}


void tracetobr(struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
       int *ipos,int *npos,int *ncon,int *sidenum,int tpos[],
       int startnum,int *metend,float *aveden)
/* Traces out simple points until either reach branch or end
  NEW NEW aveden is highest ave den on next branch 
  New returns aveden which is either ave on minimum den in next 2 0r 3 pts.
  1) ipos is previous point
  2)  npos is next point
 Returns:
  3) ncon of last point = either 1 or >2 
  4) sidenum = number of side points up to end
  5) tpos[1,2,...  upto to sidelim] the pos of each point
  6 and 7) ipos and npos are the last points before exiting
  8) If hit end of another fit, set ext.metend = to geo.set[] = fitnum
       and then only continue along points along this fit
*/
{
  int icon,j,ivec[3],iset,icapos,nden;
  float fvec[3],dist,fden;

   tpos[0]= *ipos;
   if(verbose>=2)
     printf("Enter tracetobr:  ipos = %d  npos = %d\n",*ipos,*npos);
#if 0  // for degugging
if(*ipos==10592)
  printf("\n");
#endif
#if 0  //15dec99
   if(*metend>0){ //THIS SHOULD NO LONGER HAPPEN: if have met another fit - only extend along this fit
       d1tod3(bone1.cmax,bone1.ptv[*npos].pos,ivec);
       if(geo.set[ivec[0]][ivec[1]][ivec[2]] !=  *metend){ // If not set, treat as sidechain of 0 length
          *sidenum = 0;
          *ncon = 1;
          goto quitrout;
       }
   }
#endif
   *sidenum =1;
   tpos[*sidenum] = *npos; // first point is branch point 
   *ncon = bone1.ptv[*npos].num; // number of connections
   *aveden = 0.0; 
   nden =0;
   if( *ncon==1)
     *aveden = -1000.0;
   while(*ncon==2){
      for(j=0;j< *ncon;j++){
            icon = bone1.ptv[*npos].con[j]; // the pt num of the connect point
            if(icon!= *ipos)// do not go backward
              goto foundfor;
      }// j
      foundfor:;
      d1tod3(bone1.cmax,bone1.ptv[icon].pos,ivec);
#if 0
if(icon==15376){
   printf("ivec = %d %d %d\n",ivec[0],ivec[1],ivec[2]);
   exit(0);
}
#endif
#if 1 // NEW aveden1 is is aveden for next 3 pts.
      if(*sidenum<=3){ // only ave den for first 2 or 3 pts
          nden = nden +1;     
          fden = (float)den1.den[ivec[0]][ivec[1]][ivec[2]];
          *aveden = *aveden + fden;
      }
      *aveden = *aveden/nden;
#endif
      iset = geo.set[ivec[0]][ivec[1]][ivec[2]];
      if(geo.useset==0) //This is only used for checkifbranch()
         iset = -10;
#if 0
      if(*metend>0){ // if have met another fit - only extend along this fit
         if(iset == *metend){ 
                addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
                goto repeat;
         }
      }
      else if(iset != -1){ //side and central main pts set = -1 if not previous set
#endif
      if(iset!= -1) { //side and central main pts set = -1 if not previous set
       if(iset >0){ // hit an end of a fit num = geo.set
        if(iset ==ext.ifit){ //  the same fit as started at - check if just begnning
    if(verbose>=2)
     printf(" ISET = IFIT, SIDENUM = %d STARTNUM = %d\n",*sidenum,startnum);
             if(*sidenum + startnum < 5){ // if just beginning
                   addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
                   goto repeat;
             }
          }
          else{ // if(iset !=ext.ifit)  Not the same fit as started at
                 imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
                 if(geo.fixend ==2) // if extending from Cterm, must hit nterm end of fit
                    icapos = 0;
                 else // if extending from nterm, must hit cterm end (Ca) of fit
                    icapos= fit[iset].nres -1;
                 dist = distv(fvec,fit[iset].pos[icapos][1]); // Dist from cterm Ca to icon
    if(verbose>=2)
      printf(" In tracetobr dist = %f  icapos =%d icon = %d\n",dist,icapos,icon);
//                 if(dist<=2.0){ //If close enough to end Ca then has made correct contact
                 if(dist<=2.3){ //New 17jul00 - increase from 2.0 to 2.5  If close enough to end Ca then has made correct contact
                      if(verbose>=2)
                          printf(" Setting *metend = %d in tracetobr\n",geo.set[ivec[0]][ivec[1]][ivec[2]]);
                      *metend = geo.set[ivec[0]][ivec[1]][ivec[2]]; // the fit number it hit
                      addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
                      *ncon = 1; //stop  Treat as simple end point so leaves loop
                      goto repeat;                     
                 }
            } //iset !=ext.ifit)
        } //iset>0
        else{ // set != -1 and <0 - i.e. = -10
               addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
               goto repeat;
        }
      } //iset != -1
      else{// iset == -1  addpt and quit
          addpt(bone1,icon,sidenum,tpos,ipos,npos,ncon);
          *ncon = -1; // Indicates met a previously fit point - force to be side
          goto quitrout;
      }
     *ncon = 1; //If reach here - stop  Treat as simple end point
      repeat:;
  } // while
  if(*metend>0)
      *ncon = 3;
  if(verbose>=2)
    printf("sidenum = %d\n",*sidenum);
  quitrout:;
#if 0 // use this if aveden is actual average
  if(nden >=1)
    *aveden = *aveden/nden;
#endif
}


#if 1
void add2branchs(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
    struct mainbone& main,int mainnum,int tpos[4][MAXSIDEPTS],int mcon[],
    int sidenum[],float aveden1[],float aveden2[],int metoneend[]) 
/* 
New NEW: if 2 main branchs :
    1) go to branch and run tracetobr to find next branchs
    2) For each of these next branches - find the average density
         Record the highest average density
    3) Extend first on the branch that has highest average density on 
         branch         
      (avegabed over first 3 points in tracetobr)
    If same aveden - then extendf first on short one
Extend along 2 branches  - either as side or main chains
New: if 2 main branchs - extend first on one with higher ave density
      (avegabed over first 3 points in tracetobr)
    If same aveden - then extendf first on short one
Extend along 2 branches  - either as side or main chains
*/
{
   int i,j,jmain,jside,startnum,ivec[3],iset,endnum,origmainnum;
  if(verbose>=2)
   printf("In add2branch:  tpos[0][0,1] = %d %d tpos[1][0,1] = %d %d\n",
     tpos[0][0], tpos[0][1], tpos[1][0], tpos[1][1]);
//if(tpos[0][0]==17288)
// printf("\n");
   startnum = main.startnum+main.num; // initial value
   origmainnum = main.num-1;
   if( (mcon[0]<2)&&(mcon[1]<2) ){ // Case I - both side chains -
      if(main.metend>0){ // if metend, then j= 0 is main
         jmain = 0;
         jside = 1;
      }
      else if( sidenum[0] >= sidenum[1] ){
         jmain = 0;
         jside = 1;
      }
      else{
         jmain =1;
         jside = 0;
      }
      if( mcon[jmain] == -1){ // must be side chain
          addsidechain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jmain],sidenum[jmain],origmainnum);
          if(mcon[jside] == -1) // must be side chain
              addsidechain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside],sidenum[jside],origmainnum);
          else // add to mainchain
               addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside],sidenum[jside],mcon[jside]);
      }
      else{ // mcon = 1
         addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jmain],sidenum[jmain],mcon[jmain]);
         if( (sidenum[jside] <= ext.sidelim)|| (mcon[jside] == -1) )  // create sidechain with these points 
             addsidechain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside],sidenum[jside],origmainnum);
         else  // sidenum > sidelim - make a simple mainchain
            startnewmain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside][0],
               tpos[jside][1],startnum,origmainnum);
      }
      if(*reachedlim>0)
           goto quitrout;
   } // Case I
#if 1  // new
   else if ( (mcon[0]>2) && (mcon[1]>2) ){ // Case II both main chains (May have met end)- assign br. with higher minden to orig. main chain and add new main chain    
#if 1 // New - 30jun00
      if(metoneend[0]>0){
           jmain=0;
           jside=1; 
      }
      else if(metoneend[1]>0){
           jmain=1;
           jside=0; 
      }
#endif
      else if( (aveden2[1]< -100.0)&&(aveden2[0]> -100) ){//< -100 indicates branch ends after next branch
            jmain =0;
            jside = 1;
      }
      else if( (aveden2[0]< -100.0)&&(aveden2[1]> -100) ){//< -100 indicates branch ends after next branch
            jmain =1;
            jside = 0;
      }
#if 1   
      else if(aveden1[0]> 1.2*aveden1[1]){  // First look for large differencein aveden1
         jmain = 0;
         jside = 1;
      }
      else if (aveden1[1]>1.2*aveden1[0]){ 
         jmain = 1;
         jside = 0;
      }
#endif
      else if( ((sidenum[0]>3)||(sidenum[1]>3))
//      else if( ((sidenum[0]>4)||(sidenum[1]>4))
               &&(sidenum[0]!=sidenum[1]) ){ // if large number of pts before branch use shortest
         if(sidenum[0]<sidenum[1]){
            jmain =0;
            jside = 1;
         }
         else{
            jmain =1;
            jside = 0;
         }
      } 
      else if(aveden2[0]> aveden2[1]){  //Then look at aveden2
         jmain = 0;
         jside = 1;
      }
      else if (aveden2[1]>aveden2[0]){ 
         jmain = 1;
         jside = 0;
      }
      else{ // both branches same density - assigne shorter to main
         if(sidenum[0]<sidenum[1]){
            jmain =0;
            jside = 1;
         }
         else{
            jmain =1;
            jside = 0;
         }
      } // same den
#if 1  // New 16dec99
      if(metoneend[jmain]>0){ // make jmet
        ext.metend = main.metend = metoneend[jmain];
        mcon[jmain]=1; 
        addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jmain],
             sidenum[jmain],mcon[jmain]);
        goto quitrout;
      } // metend
#endif
      addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jmain],sidenum[jmain],mcon[jmain]);
      if(*reachedlim>0){ //New: used all mains, throw away new ones and got to start new main on jside
//          goto quitrout;
#if 0
            if (ext.nmain >= MAXMAIN-1){
              if(verbose>=2)
                  printf(" RESETTING ext.main from %d to %d\n",ext.nmain,mainnum);
               ext.nmain = mainnum;
               ext.metend = -1;  //Reset in case metend in previus addbrtomain
               for(i=0;i<4;i++)
               metoneend[i] = -1;
                goto startnew;
            }
#endif
            if (ext.nmain >= MAXMAIN-1){ //New -30jun00
                goto quitrout;
            }
      }
      d1tod3(bone1.cmax,bone1.ptv[ tpos[jside][1] ].pos,ivec);
      iset = geo.set[ivec[0]][ivec[1]][ivec[2]];
      if(iset == -1){ // Check if PREVIOUS BRANCH HAS BECOME SIDE OF ANOTHER
        if(verbose>=2)
         printf(" PREVIOUS BRANCH HAS BECOME SIDE OF ANOTHER\n");
         for(i=1;i<=ext.nmain;i++){
            for(j=1;j<=ext.main[i].nside;j++){ // each side for this main
              endnum = ext.main[i].side[j].ipos[ ext.main[i].side[j].num-1 ]; // bone num of end of side
              if(endnum == tpos[0][0]){
                  main.mainpt[origmainnum].ctype[0] = 1;
                  main.mainpt[origmainnum].ctype[1] = i; // main bone side originates from 
                  main.mainpt[origmainnum].ctype[2] = 
                         ext.main[i].side[j].mnum;// main pt num side starts at 
                  main.mainpt[origmainnum].ctype[3] = j; // sidenum that connects
                  ext.main[i].side[j].ctype[0] = 1;
                  ext.main[i].side[j].ctype[1] = mainnum;
                  ext.main[i].side[j].ctype[2] = origmainnum;
  if(verbose>=2)
    printf(" This connection from main %d side %d\n",i,j);
                  goto quitrout;
              }
            }//j
         } // i
         goto quitrout;
      }
//      main.mainpt[origmainnum].mtype = ext.nmain+1; //Indicates new main starts at this oint
      startnew:;
      if(startnum<MAINONEPTS)
         startnewmain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside][0],
             tpos[jside][1],startnum,origmainnum);
   } // CAse II
#endif
   else{ // Case III - one main,mainnum, one side
#if 0
      if(main.metend>0){ // if metend, then j= 0 is main
         jmain = 0;
         jside = 1;
      }
      else if( mcon[0]>2){
#endif
      if( mcon[0]>2){
         jmain = 0;
         jside = 1;
      }
      else{
         jmain = 1;
         jside = 0;
      }
#if 1  // New 16dec99
   if(metoneend[jmain]>0){ // make jmet
     ext.metend = main.metend = metoneend[jmain]; 
     mcon[jmain]=1;
     addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jmain],
             sidenum[jmain],mcon[jmain]);
     goto quitrout;
   } // metend
#endif
     addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jmain],sidenum[jmain],mcon[jmain]);
      if(*reachedlim>0)
           goto quitrout;
      if( (sidenum[jside] <= ext.sidelim)|| (mcon[jside] == -1) ) { // create sidechain with these points 
         addsidechain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside],sidenum[jside],origmainnum);
         if(*reachedlim>0)
           goto quitrout;
      }
      else{  // sidenum > sidelim - make a simple mainchain
//        main.mainpt[origmainnum].mtype = ext.nmain+1; //Indicates new main starts at this oint
        startnewmain(reachedlim,den1,geo,bone1,main,mainnum,tpos[jside][0],
              tpos[jside][1],startnum,origmainnum);
      }
   }
   quitrout:;
}
#endif

void nextbr(int *reachedlim,struct griddata& den1,struct geometry& geo,
    struct bonedata& bone1,
    struct mainbone& main,int mainnum,int ipos,int npos,int ncon)
/* Recursive relation which moves from branch to branch
NEW - do not set metend if tracetobr finds 2mains
  ASSUME MAX OF 2 BRANCHS - (3 CONNECTIONS)
  ALSO - cut out metend and nmet, no longer used.
  Case I) If both branches do not have another branch, then
      assign longer side to main,mainnum, shorter to side,
      Later will try to continue this using lower density 
  Case II) Both branches have branch- assign one to main,mainnum, and for
     other start new mainchain by calling firstpt
  Case III) If one branch has branch, other doesent,
      assign branched branch to main - other to side and call nextbr
  Case IV) If side is > ext.sidelim in length, then, treat as a 
     another mainchain
*/
{
   int icon,mcon[4],sidenum[4],tpos[4][MAXSIDEPTS],iposm[4],nposm[4];
   int ii,j,jj,ivec[3],iset,metend,metoneend[4],jmet;
   float aveden1[4],aveden2[4];

if(verbose>=2){
  printf(" Enter nextbr: npos = %d\n",npos);
  if (npos==1718)
     printf("DEbug 4\n");
}
   if(ncon >4){ 
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
   jj=0; // counts forward branches
   metend = 0;
   for(j=0;j<ncon;j++){
       icon = bone1.ptv[npos].con[j]; // the pt num of the connect point
       if(icon!=ipos){ // go forward
           d1tod3(bone1.cmax,bone1.ptv[icon].pos,ivec);
           iset = geo.set[ivec[0]][ivec[1]][ivec[2]];
           if(iset == -1){ // Then quit this trace
              sidenum[jj]=1; 
              tpos[jj][0] = npos;
              tpos[jj][1] = icon;
              mcon[jj] = -1; // Indicates met prevous set point - force to be side
           } 
           else{ 
              iposm[jj]= npos;
              nposm[jj]= icon;
              metoneend[jj] =metend = 0;
              tracetobrrec(den1,geo,bone1, &iposm[jj],&nposm[jj],&mcon[jj],
                 &sidenum[jj],tpos[jj],main.startnum+main.num,
                 &metend,&aveden1[jj],&aveden2[jj]); // New recursive version
              metoneend[jj]=metend;
#if 0 // comment out 16dec99
              if(metoneend[jj]>0){
                 jmet = jj;
                 metend = metoneend[jj];
                 jj = jj+1; // because skip line below
                 goto quitfor;
              }
#endif
           }
           jj = jj+1;
       }
   } // j
   quitfor:;
   float pos[3],posden,lowestden[3];
   int endpt,imain,kk;
#if 0// comment out 16dec99
   if(metend>0){ // make jmet
     main.metend = metend; 
     addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,
            tpos[jmet],sidenum[jmet],mcon[jmet]);
     goto quitrout;
   } // metend
#endif
#if 1
   for(ii=0;ii<jj;ii++){// check and see if end of branch is on previous man - if it is 
                              //quit this main because previous main will extend into it
      if(mcon[ii] == -1){ // indicates met a point on a previous main - delete the current main
         endpt = tpos[ii][sidenum[ii]];
         for(imain = 1;imain< ext.nmain;imain++){ // check and see if end of branch is on previous man
            for(j=0;j<ext.main[imain].num;j++){ // each saved othemain pt
               if(endpt == ext.main[imain].mainpt[j].num){ // then metmain - throw away current main
                   goto foundmain;
               }
            } //j
         } // imain
         foundmain:;
         if(verbose>=2)
            printf("Main %d MET PREVIOUS on MAIN %d\n",mainnum,imain);
#if 0
         if(mainnum!=imain){
//            *reachedlim = 1; // HAVe commented this out - do not know why needed
             goto quitrout;
         }
         else { //mainnum = imain: Came back on itself,  just continue on forward branch 
#endif
         for(kk=0;kk<jj;kk++){
              if(mcon[kk] != -1){
#if 1  // New 16dec99
                if(metoneend[kk]>0){ 
                  ext.metend = main.metend = metoneend[kk]; 
                  mcon[kk] = 1;
                }
#endif
                addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,
                      tpos[kk],sidenum[kk],mcon[kk]);
                goto quitrout;
              }// mcon
            } //kk
//         } // else
      } // mcon= -1              
   } // ii 
#endif  
   if(tpos[0][sidenum[0]] == tpos[1][sidenum[1]]){ // If a loop, only use branch with best minden
      for(ii=0;ii<=1;ii++){
        lowestden[ii] = 1000.0;
        for(j=1;j<=sidenum[ii];j++){
           d1tod3(bone1.cmax,bone1.ptv[ tpos[ii][j] ].pos,ivec);
           imcv(den1,ivec,pos);// convert int grid ivec to float world co-ord fvec
           posden = avepos(geo,den1,pos);
           if(posden < lowestden[ii])
              lowestden[ii] = posden;
        }
      } // ii
      if(verbose>=2)
         printf("LOOP  to pos = %d  lowestden[0,1]=%5.3f %5.3f \n",
                  tpos[0][sidenum[0]],lowestden[0],lowestden[1]);
      if( lowestden[0]<lowestden[1])
          jmet =1;
      else
          jmet = 0;
      addbrtomain(reachedlim,den1,geo,bone1,main,mainnum,
            tpos[jmet],sidenum[jmet],mcon[jmet]);
   } // if loop
   else { // did not meet end
      if(jj>2){ // assume max of 2 branches  (jj = 0,1)
         printf(" ERROR: jj = %d > 2  setting = 2\n",jj);
         jj=2;
      }
      add2branchs(reachedlim,den1,geo,bone1,main,mainnum,tpos,mcon,
         sidenum,aveden1,aveden2,metoneend);
   }
   quitrout:;
}    
      
    
void findclose(struct geometry& geo,int na[],int ifit,float endpos[3],
     int *closeres,int *closeatom,float *closedist)
/* 
Finds the closest fit point to endpos
  Starts with end residue and moves towards center
  Checks all atoms in residue, keeps track of minimum distance
     and corresponding atoms #
  Continues until min. dist increases.
  Then returns close atom num ii  (closeres, closeatom)
  Note: returned number is ii, actual atom num = na[ii]
NOTE: HAVE ONLY DEBUGGED CTERM END
       na[1] = 4; // order atoms in residue = N - Ca-  C
       na[2] = 1;
       na[3] = 2;
       na[0] = 3; // # of atoms
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int ii,resatom,ires,iresmin;
   float resmin,dist;

#if 0
int istep;
   if(geo.fixend==2){ // extending at C term end--then must have met nterm end
      ires= -1;  // first residue tested = ires+istep
      istep = 1;
   }
   else {  // extending at N term end--then must have met cterm end
   }
#endif
   resmin =10000.0;
   for(ires = 0;ires<=3;ires++) { // try first 3 residuews - NOTE ONLY CORRECT FOR CTERM EXTENSION
     for(ii= 1;ii<=na[0];ii++){
       if( (ires ==0) && (na[ii]==4) ) // N not defined for first res.
           goto nextii;
       dist = distv(endpos,fit[ifit].pos[ires][na[ii]]);
       if(dist<resmin){
            resmin = dist;
            resatom = ii;
            iresmin = ires;
       }
       nextii:;
     }// ii
   }// ires
   *closeres = iresmin;
   *closeatom = resatom;
   *closedist = resmin;
}


void firstpt(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
    struct mainbone& main,int mainnum,int ipos,int npos)
/* Recursive routine that given first 2 points on main (ipos,npos)
   of main - then traces out end -
  If it finds a branch point - either one branch is main and other side
    or both main,mainnum, in which case, it starts another main trace by calling
  firstpt recursively
*/
{
    int j,ncon,sidenum,tpos[MAXSIDEPTS];
    float vec[3],aveden1,aveden2;

    main.nside = 0; // No side branches yet
    d1tow(den1,bone1.cmax,bone1.ptv[ipos].pos,vec);// convert pos to world nvec[3]
    main.num = 0; // First point on main = 0
    addmainpt(den1,geo,bone1,main,mainnum,ipos,0,0);// add first point - set mtype = ext.nmain
    if(verbose>1)
       printf("Calling tracebrrec from firstpt()");
    tracetobrrec(den1,geo,bone1, &ipos,&npos,&ncon,&sidenum,tpos,
            main.startnum,&main.metend,&aveden1,&aveden2); // new recursive version
    for(j=1;j<=sidenum;j++){ // add  points  to initial main - starting with npos
       if (main.num >= MAXMAINPTS-1){
         printf("REACHED LIMIT: main.num = %d >= MAXMAINPTS-1\n",main.num);
         *reachedlim = 1;
         goto quitrout;
       }
       addmainpt(den1,geo,bone1,main,mainnum,tpos[j],0,0);// NOTE: set mtype = stype = 0 for branch point - this is changed later
    }
//    if(ncon==1 ){ // reached end without reaching branch (or met another main) - just add this point and quit

#if 1  // New - added on 10may00
    if(main.metend>0){
        ext.metend = main.metend;
        goto quitrout;
    }
#endif
    if( (ncon==1)||(ncon == -1) ){ // reached end without reaching branch (or met another main) - just add this point and quit
       goto quitrout;
    }
    else{ 
      nextbr(reachedlim,den1,geo,bone1,main,mainnum,ipos,npos,ncon);
      if(*reachedlim >0)
          goto quitrout;
    }
    quitrout:;
}
         
void startnewmain(int *reachedlim,struct griddata& den1,struct geometry& geo,struct bonedata& bone1,
     struct mainbone& main,int mainnum,int ipos,int npos,int startnum,int origmainnum)
// Start another main trace off of previous main
{
    if (ext.nmain >= MAXMAIN-1){
       if(verbose>=1)
          printf("REACHED LIMIT: ext.nmain = %d>= MAXMAIN-1\n",ext.nmain);
       *reachedlim = 1;
       goto quitrout;
    }
    ext.nmain = ext.nmain+1; // add another main chain
    ext.main[ext.nmain].zeronum = 0; // default start of main
    ext.main[ext.nmain].mainnum = origmainnum;
    ext.main[ext.nmain].startmain = mainnum;
    ext.main[ext.nmain].startnum = startnum;// startnum of this main
 if(verbose>=2)
   printf("Enter startnewmain:  nmain = %d  startnum = %d\n",ext.nmain,startnum);
    firstpt(reachedlim,den1,geo,bone1,ext.main[ext.nmain],ext.nmain,ipos,npos);
    if(*reachedlim>0)
       goto quitrout;
    main.mainpt[origmainnum].mtype = ext.nmain; //Now set in calling routine - Set mtype for main pt
 if(verbose>=2)
   printf("Leave startnewmain:  nmain = %d  startnum = %d\n",ext.nmain,startnum);
   quitrout:;
}

void copyoneeditbones(struct mainonebone& mainone1,struct mainonebone& mainone2)
// Just copies mainone1 to mainone2 
{
   int j,kk;

if(verbose>=2)
  printf("Copying onemain1 to onemain2\n");
  mainone2.zeronum = 0; // default value
  mainone2.num = mainone1.num;
  for(j=0;j<mainone1.num;j++){ //  the main pts
          mainone2.mainpt[j].num=mainone1.mainpt[j].num;
          mainone2.mainpt[j].mtype=mainone1.mainpt[j].mtype;
          for(kk=0;kk<=2;kk++)
            mainone2.mainpt[j].pos[kk]=mainone1.mainpt[j].pos[kk];
  }
}

void copyeditbones(struct extbones& ext1,struct extbones& ext2)
// Copy ext1 to ext2
{
   int i,j,kk,k;

if(verbose>=2)
  printf("Copying ext1 to ext2\n");
    ext2.length = ext1.length;
    ext2.metend = ext1.metend;
    ext2.nmain = ext1.nmain;
    for(i=1;i<=ext1.nmain;i++){ // each main
        ext2.main[i].zeronum = 0; // default value
        ext2.main[i].num = ext1.main[i].num;
        for(j=0;j<ext1.main[i].num;j++){ //  the main pts
          ext2.main[i].mainpt[j].num=ext1.main[i].mainpt[j].num;
          ext2.main[i].mainpt[j].mtype=ext1.main[i].mainpt[j].mtype;
          for(kk=0;kk<=2;kk++)
            ext2.main[i].mainpt[j].pos[kk]=ext1.main[i].mainpt[j].pos[kk];
        }
        ext2.main[i].nside=ext1.main[i].nside;
        for(j=1;j<=ext1.main[i].nside;j++){ // each side for this main
          ext2.main[i].side[j].num=ext1.main[i].side[j].num;
          ext2.main[i].side[j].mnum=ext1.main[i].side[j].mnum;
          for(k=0;k<ext1.main[i].side[j].num;k++){ 
             ext2.main[i].side[j].ipos[k]=ext1.main[i].side[j].ipos[k];  
             for(kk=0;kk<=2;kk++)
                ext2.main[i].side[j].pos[k][kk]=ext1.main[i].side[j].pos[k][kk];
          }
        }// j  
   }// i
}

void readeditbones(const char editbone[])
{
   int i,j,k,kk;
   FILE *fp;

 if(verbose>=2)
   printf("Reading edited bones to file 'editbones'\n");
   if( (fp=fopen(editbone,"r")) == NULL)
      fprintf(fpout,"ERROR: There is no editbones file\n" );
   fscanf(fp,"%d",&ext.metend);
   fscanf(fp,"%d",&ext.nmain);
  if(verbose>=2)
     printf("Reading %d edited main bones from file 'editbones'\n",ext.nmain);
 if(ext.nmain == -150){// this is indicator that just save onemain
    ext.mainone.zeronum = 0; // default value
    fscanf(fp,"%d",&ext.mainone.num);
    for(j=0;j<ext.mainone.num;j++){ // each  mainpt
        fscanf(fp,"%d",&ext.mainone.mainpt[j].num);
    }
    for(j=0;j<ext.mainone.num;j++){ // each remaining mainpt
        for(kk=0;kk<=2;kk++)
           fscanf(fp,"%f",&ext.mainone.mainpt[j].pos[kk]);
    }
  }// ext.nmain == -150
  else{ // standard procedure
   for(i=1;i<=ext.nmain;i++){ // each main
        ext.main[i].zeronum = 0; // default value
        fscanf(fp,"%d",&ext.main[i].num);
        fscanf(fp,"%d",&ext.main[i].startnum);
 if(verbose>=2)
   printf("main[%d] has %d bone point  startlenght = %d totlength= %d\n",
           i,ext.main[i].num,ext.main[i].startnum,ext.main[i].startnum+ext.main[i].num);
        for(j=0;j<ext.main[i].num;j++){ //  the main pts
           fscanf(fp,"%d",&ext.main[i].mainpt[j].num);
           fscanf(fp,"%d",&ext.main[i].mainpt[j].mtype);
        }
        for(j=0;j<ext.main[i].num;j++) //  the main pts
          for(kk=0;kk<=2;kk++)
             fscanf(fp,"%f",&ext.main[i].mainpt[j].pos[kk]);
        fscanf(fp,"%d",&ext.main[i].nside); // number of sides
        for(j=1;j<=ext.main[i].nside;j++){ // each side for this main
          fscanf(fp,"%d",&ext.main[i].side[j].num);
          fscanf(fp,"%d",&ext.main[i].side[j].mnum);
          for(k=0;k<ext.main[i].side[j].num;k++)   
             fscanf(fp," %d",&ext.main[i].side[j].ipos[k]);
          for(k=0;k<ext.main[i].side[j].num;k++)   
             for(kk=0;kk<=2;kk++)
                fscanf(fp,"%f",&ext.main[i].side[j].pos[k][kk]);
        }// j  
   }// i
 } // standard procedure
 fclose(fp);
if(verbose>=2)
 printf("ext.main[1].nside = %d ext.main[1].side[1].mnum = %d  ext.main[1].num = %d\n",
            ext.main[1].nside,ext.main[1].side[1].mnum,ext.main[1].num);
}


void saveeditbones(const char editbone[],int testmetend)
/*
New: if testmetend !=1, then just save complete set
NOTE: FOUND ERROR: (12oct99) - .mtype is number of main that starts
    at mainpt - must be sure to give it correct value when renumber mains
    oldmain[i] = 0 = mype if this mtype is no longer included 
      else = new main number
 NOTE: 
   if metend>0, a side chain may come after main is terminated
     That is, not all these side chains will touch mains
     This should be corrected by recalculating side chains
   NEED TO DO THIS - NOT DONE
*/
{
   int count,i,j,k,kk,oldmain[MAXMAIN],sidecount,sidebr;
   FILE *fp;
int ii;

if(verbose>=2)
   printf("Writing edited bones to file 'editbones'\n");
   fp=fopen(editbone,"w");
   fprintf(fp,"%d\n",ext.metend);
   if(ext.nmain == -150){// this is indicator that just save onemain
       fprintf(fp,"%d\n",ext.nmain);
       fprintf(fp,"%d\n",ext.mainone.num );
       for(j=0;j<ext.mainone.num;j++){ // each  mainpt
         fprintf(fp," %d",ext.mainone.mainpt[j].num);
       }
       fprintf(fp,"\n");
       for(j=0;j<ext.mainone.num;j++){ // each remaining mainpt
          for(kk=0;kk<=2;kk++)
            fprintf(fp," %6.3f",ext.mainone.mainpt[j].pos[kk]);
          fprintf(fp,"\n");
       }
    }// ext.nmain == -150
#if 0 // OLD VERSion before 3nov00
    else if( (ext.metend>0)&&(testmetend==1) ){// Then just write the main bones on path to fit
     if(verbose>=2)
        printf("ext.metend = %d\n",ext.metend);
      count = 0;
      for(i=1;i<=ext.nmain;i++){ // count the metend mains
        oldmain[i] = 0; 
        if(ext.main[i].metend>0){
           count= count+1;
           oldmain[i] = count; // keep track of relation of new number to old
        }
      }
      fprintf(fp,"%d\n",count);
      for(i=1;i<=ext.nmain;i++){ // each main
        if(ext.main[i].metend>0){
           sidecount = 0;
           for(j=ext.main[i].zeronum;j<ext.main[i].num-ext.main[i].zeronum;j++){ //  count the number sides left after editing
              if(ext.main[i].mainpt[j].stype>0)
                 sidecount = sidecount+1;
           }
           fprintf(fp,"%d\n",(ext.main[i].num-ext.main[i].zeronum) );
           fprintf(fp,"%d\n",ext.main[i].startnum);
            for(j=ext.main[i].zeronum;j<ext.main[i].num-ext.main[i].zeronum;j++){ // each remaining mainpt
             fprintf(fp," %d",ext.main[i].mainpt[j].num);
             if(ext.main[i].mainpt[j].mtype>0) // if main starts here - use its new  number
                 fprintf(fp," %d",oldmain[ext.main[i].mainpt[j].mtype]);
             else // mtype = 0;
                 fprintf(fp," %d",ext.main[i].mainpt[j].mtype);
           }
           fprintf(fp,"\n");
           for(j=ext.main[i].zeronum;j<ext.main[i].num-ext.main[i].zeronum;j++){ // each remaining mainpt
             for(kk=0;kk<=2;kk++)
                fprintf(fp," %6.3f",ext.main[i].mainpt[j].pos[kk]);
             fprintf(fp,"\n");
           }
           fprintf(fp,"%d\n",sidecount); // number of sides
           for(j=ext.main[i].zeronum;j<ext.main[i].num-ext.main[i].zeronum;j++){ // each remaining mainpt
               sidebr = ext.main[i].mainpt[j].stype; // either 0 or side branch number
               if(sidebr>0){ // then write data for this side
                  fprintf(fp,"%d\n",ext.main[i].side[sidebr].num);
                  fprintf(fp,"%d\n",ext.main[i].side[sidebr].mnum);
                  for(k=0;k<ext.main[i].side[sidebr].num;k++)   
                      fprintf(fp," %d",ext.main[i].side[sidebr].ipos[k]);
                  fprintf(fp,"\n");
                  for(k=0;k<ext.main[i].side[sidebr].num;k++){   
                     for(kk=0;kk<=2;kk++)
                       fprintf(fp," %6.3f",ext.main[i].side[sidebr].pos[k][kk]);
                     fprintf(fp,"\n");
                  } // k
                } // if sidebr
           }// j 
         }// if metend 
       }// i
   } // standard procedure
#endif
#if 1 // New VERSion after 3nov00 - use passed variable metend
    else if( (ext.metend>0)&&(testmetend==1) ){// Then just write the main bones on path to fit
     if(verbose>=2)
        printf("ext.metend = %d\n",ext.metend);
      count = 0;
      for(i=1;i<=ext.nmain;i++){ // count the metend mains
        oldmain[i] = 0; 
        if(ext.main[i].metend>0){
           count= count+1;
           oldmain[i] = count; // keep track of relation of new number to old
        }
      }
      fprintf(fp,"%d\n",ext.metmain[0]);
      for(ii=ext.metmain[0];ii>=1;ii--){ // Mains in order of path to metend
        i=ext.metmain[ii];
           sidecount = 0;
           for(j=ext.main[i].zeronum;j<ext.main[i].num;j++){ //  count the number sides left after editing
              if(ext.main[i].mainpt[j].stype>0)
                 sidecount = sidecount+1;
           }
           fprintf(fp,"%d\n",(ext.main[i].num-ext.main[i].zeronum) );
           fprintf(fp,"%d\n",ext.main[i].startnum);
            for(j=ext.main[i].zeronum;j<ext.main[i].num;j++){ // each remaining mainpt
             fprintf(fp," %d",ext.main[i].mainpt[j].num);
             if(ext.main[i].mainpt[j].mtype>0) // if main starts here - use its new  number
                 fprintf(fp," %d",oldmain[ext.main[i].mainpt[j].mtype]);
             else // mtype = 0;
                 fprintf(fp," %d",ext.main[i].mainpt[j].mtype);
           }
           fprintf(fp,"\n");
           for(j=ext.main[i].zeronum;j<ext.main[i].num;j++){ // each remaining mainpt
             for(kk=0;kk<=2;kk++)
                fprintf(fp," %6.3f",ext.main[i].mainpt[j].pos[kk]);
             fprintf(fp,"\n");
           }
           fprintf(fp,"%d\n",sidecount); // number of sides
           for(j=ext.main[i].zeronum;j<ext.main[i].num;j++){ // each remaining mainpt
               sidebr = ext.main[i].mainpt[j].stype; // either 0 or side branch number
               if(sidebr>0){ // then write data for this side
                  fprintf(fp,"%d\n",ext.main[i].side[sidebr].num);
                  fprintf(fp,"%d\n",ext.main[i].side[sidebr].mnum);
                  for(k=0;k<ext.main[i].side[sidebr].num;k++)   
                      fprintf(fp," %d",ext.main[i].side[sidebr].ipos[k]);
                  fprintf(fp,"\n");
                  for(k=0;k<ext.main[i].side[sidebr].num;k++){   
                     for(kk=0;kk<=2;kk++)
                       fprintf(fp," %6.3f",ext.main[i].side[sidebr].pos[k][kk]);
                     fprintf(fp,"\n");
                  } // k
                } // if sidebr
           }// j 
       }// ii
   } // standard procedure
#endif
   else{ // ext.metend = -1 or testemetend = -1
     fprintf(fp,"%d\n",ext.nmain);
     for(i=1;i<=ext.nmain;i++){ // each main
        fprintf(fp,"%d\n",ext.main[i].num);
        fprintf(fp,"%d\n",ext.main[i].startnum);
        for(j=0;j<ext.main[i].num;j++){ //  the main pts
           fprintf(fp," %d",ext.main[i].mainpt[j].num);
           fprintf(fp," %d",ext.main[i].mainpt[j].mtype);
        }
        fprintf(fp,"\n");
        for(j=0;j<ext.main[i].num;j++){ //  the main pts
             for(kk=0;kk<=2;kk++)
                fprintf(fp," %6.3f",ext.main[i].mainpt[j].pos[kk]);
             fprintf(fp,"\n");
        }
        fprintf(fp,"%d\n",ext.main[i].nside); // number of sides
        for(j=1;j<=ext.main[i].nside;j++){ // each side for this main
//printf("for main %d side %d num = %d mnum = %d\n",
//    i,j,ext.main[i].side[j].num, ext.main[i].side[j].mnum);
           fprintf(fp,"%d\n",ext.main[i].side[j].num);
           fprintf(fp,"%d\n",ext.main[i].side[j].mnum);
           for(k=0;k<ext.main[i].side[j].num;k++)   
             fprintf(fp," %d",ext.main[i].side[j].ipos[k]);
           fprintf(fp,"\n");
           for(k=0;k<ext.main[i].side[j].num;k++){   
               for(kk=0;kk<=2;kk++)
                  fprintf(fp," %6.3f",ext.main[i].side[j].pos[k][kk]);
               fprintf(fp,"\n");
           }
        }// j  
     }
   } // else
   fclose(fp);
}

void increasemainden(struct griddata& den1,struct bonedata& bone1)
/*
 If metend>0:
   move main to neighbor pt if it has higher den
  Go through the main pts (not ends) -
  Example:  pt i:
    Check all neighbors j of i:
      if j not = pt i-1 or pt i+1 and j with +-1 of i-1 and i and
        if den[j] > den [i] - then move pos of i to j.
*/
{
   int i,num[3],ivec[3][3],j,jj,j1,i0,i1,i2,k,kk,itest[3];
   int inum;

 if(ext.metend>0){
   for(i=1;i<=ext.nmain;i++){ // each main
     for(j=1;j<ext.main[i].num-1;j++){ //  the main pts - excluding ends
        for(jj= j-1;jj<=j+1;jj++){
           inum =   ext.main[i].mainpt[jj].num;    
           j1= jj- (j-1);
           num[j1] =  bone1.ptv[inum].pos;
           d1tod3(bone1.cmax,num[j1],ivec[j1]); // convert from 1d to 3d vec[i,j,k]
        }//jj
        for(i0 = -1;i0<=1;i0++){// test all neighbors
          for(i1 = -1;i1<=1;i1++)
            for(i2 = -1;i2<=1;i2++)
              if( (i0!=0)||(i1!=0)||(i2!=0) ){
                itest[0] = ivec[1][0]+i0;
                itest[1] = ivec[1][1]+i1;
                itest[2] = ivec[1][2]+i2;
                for(kk = 0;kk<=2;kk = kk+2){
                   if( (itest[0]==ivec[kk][0])&&(itest[1]==ivec[kk][1])&&
                        (itest[2]==ivec[kk][2]) ) // then = other main pt
                          goto nextpt;
                   for(k=0;k<=2;k++)
                     if(abs(itest[k]-ivec[kk][k]) >1) // point step too large
                        goto nextpt; 
                }// kk
                // if reach here - test point is OK
if(verbose>=2)
  printf("denneigh[%d][%d][%d] = %d   den[%d][%d][%d] = %d\n",
   itest[0],itest[1],itest[2],(int)den1.den[itest[0]][itest[1]][itest[2]],
   ivec[1][0],ivec[1][1],ivec[1][2],(int)den1.den[ivec[1][0]][ivec[1][1]][ivec[1][2]]);  
                if((int)den1.den[itest[0]][itest[1]][itest[2]]>
                   (int)den1.den[ivec[1][0]][ivec[1][1]][ivec[1][2]]){ // if OK neigbor has higher den                     
if(verbose>=2)
  printf("denneigh >den\n");
                      ext.main[i].mainpt[j].num = 
                            d3tod1(bone1.cmax,itest[0],itest[1],itest[2]);
                     d1tow(den1,bone1.cmax,num[1],ext.main[i].mainpt[j].pos);// convert pos to world fvec[3]
                     goto nextpt; 
                } // if den>
                nextpt:;
              }//if( (i0!=0)||(i1!=0)||(i2!=0) 
        } // i0               
      } // j
    } // i
  } // ext.metend
}

int extendmainonfitcheck(struct griddata& den1,struct bonedata& bone1,
             struct geometry& geo,int nearnum)
/*
Similar to extendmainonfit  -
  But is called intially in editmain if ext.metend >0.
  If closedist > 2.0; then it assumes that
   cannot be extended and return -1;
  else returns 1;
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ii,k,iend,ifit,imain,closeres,closeatom;
   int maxext,extmax,na[6],istep,ires,iatom,newbone;
   float closedist;
 
   
   maxext = 4; 
   imain =1;
   ext.main[1].num = 0;  // firsst pt
   addmainpt(den1,geo,bone1,ext.main[1],1,nearnum,0,0);// add first point - set mtype = ext.nmain

   ifit = ext.metend;  // the fit it met
   iend= ext.main[imain].num-1; //the main num of the end pt of main
   if(geo.fixend==2){ // if extending at cterm end - then must have met nterm end
       na[1] = 4; // order atoms in residue = N - Ca-  C
       na[2] = 1;
       na[3] = 2;
       na[0] = 3; // # of atoms
       istep = 1;
       if(maxext<fit[ifit].nres-2) // do not ext. into last residue
         extmax = maxext; // max # of residue extend main bone on
       else
         extmax = fit[ifit].nres-2;
   }
   else {  //if extending at nterm end 
   }
   findclose(geo,na,ifit,ext.main[imain].mainpt[iend].pos,
                         &closeres,&closeatom,&closedist);
   if(closedist>2.0)
       return -1;
   if(closeatom  == 3){
       ires = closeres+istep;
       iatom = 1;
   }
   else{
       ires = closeres;
       iatom = closeatom+1;
   }
#if 0 
   for(i=1;i<= na[0];i++){ // find first extension point on main atoms
     if(closeatom == na[i]){
       if(i==na[0]){
          ires = closeres+istep;
          iatom = 1;
       }
       else{
          ires = closeres;
          iatom = i+1;
       }
     } // if closeartom
   }// i
#endif
if(verbose>=2)
  printf("Continue extension at fit: at iatom = %d, ires = %d\n",iatom,ires);
   for(i=1;i<=extmax;i++){ // extend out extmax residues
      for(ii= iatom;ii<= na[0];ii++){
         newbone = ext.main[imain].num; // the new bonept num
         ext.main[imain].num = newbone+1; // new total number
         ext.main[imain].mainpt[newbone].mtype = 0;
         ext.main[imain].mainpt[newbone].stype = 0;
         ext.main[imain].mainpt[newbone].ctype[0] = 0;
         for(k=0;k<=2;k++)
            ext.main[imain].mainpt[newbone].pos[k] =
                   fit[ifit].pos[ires][ na[ii] ][k];
      }
      ires = ires+istep;
      iatom = 1; // start next res with first atom
   }
   return 1;
}   
     

void extendmainonfit(struct geometry& geo,int endmain)
/* Called only if metend>0.
NEW: 28mar00: Only remove last point if next to last point is closer 
      to Ca then it is to last point.
    Never remove 2 points.
New: routinely remove last 1 point before continuing
  Remove 2 points if dist from end to C is less then 1.8
ERROR FOUND ON 14OCT99:  the main that meets fit is not necessarily the
        highest numbered main 
    Fix by keeping track of main that meets fit (= endmain) and passing it when called
   Extends main along mainchain atoms the fit it met.
      for distance of 4 residues. This should quarantee that
      extendfit() which goes in steps of 3 residues, will overlap the
      fit by 1 residue, which can then be annealed.
This is wrong
  NOTE:  assume last main is one that meets fit  -
    I am not sure this is always true
Instead - find the main that has metend>=0
    mains run from i= 1,,2 .. ext,nmain;
    mainpts run from 0,1,2, .. ext.main[i].num-1
If there is a branch at next to last pt of main that meets
   end, then this branch may be correct direction.
   As cure for this, just remove last point of main 
     and extend from there.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
{
   int i,ii,k,ifit,imain,iend,iremove;
   int maxext,extmax,na[6],istep,ires,iatom,newbone;
   float dist,dist2;
 
   
   imain = 0;
   maxext = 4;
   imain = endmain;
   geo.quit = 0; // default - no problems in extension
   if(imain ==0){
      printf("ERROR: No main has metend>0\n");
      goto quitrout;
   }
   ifit = ext.main[imain].metend;  // the fit it met
//goto quitrout;
   na[0] = 3; // # of atoms
   if(maxext<fit[ifit].nres-2) // do not ext. into last residue
         extmax = maxext; // max # of residue extend main bone on
   else
         extmax = fit[ifit].nres-2;
   if(geo.fixend==2){ // if extending at cterm end - then must have met nterm end
       na[1] = 4; // order atoms in residue = N - Ca-  C
       na[2] = 1;
       na[3] = 2;
       istep = 1;
#if 1 // New 7dec99 assume close at isterm  Ca  
       ires = 0;
       iatom = 2; //na[2] = Calpha
#endif  
#if 0
       findclose(geo,na,ifit,ext.main[imain].mainpt[iend].pos,
                         &closeres,&closeatom,&closedist);
       if(closeatom  == 3){
          ires = closeres+istep;
          iatom = 1;
       }
       else{
          ires = closeres;
          iatom = closeatom+1;
       }
#endif
   }
   else {  //if extending at nterm end - must have met cterm end
       na[1] = 2; // order atoms in residue =  Ca - N  - C
       na[2] = 1;
       na[3] = 4;
       istep = -1;
       ires = fit[ifit].nres-1;
       iatom = 2; //na[2] = Calpha
   }
     iend =ext.main[imain].num;
     dist = distv(ext.main[imain].mainpt[iend-1].pos,
                   fit[ifit].pos[ires][ 2 ]); // iend-1 to Ca dist
     dist2 = distv(ext.main[imain].mainpt[iend-1].pos,
                   ext.main[imain].mainpt[iend-1].pos); //NOTE ERROR - ALWAYS 0 iend-1 to iend dist
     if(dist<dist2){// if iend-1 is closer to Ca then to iend, remove iend
       iremove = 1;
       ext.main[imain].num = ext.main[imain].num-1;
    }
    else
      iremove = 0;
#if 0
   iremove = 1;
   ext.main[imain].num = ext.main[imain].num-1;
   if(ext.main[imain].num>1){
     iend =ext.main[imain].num;
     dist = distv(ext.main[imain].mainpt[iend].pos,
                   fit[ifit].pos[ires][ 2 ]); // iend to C dist
//     if(dist < 1.8){ // remove another
     if(dist < 0.5){ // remove another
       ext.main[imain].num = ext.main[imain].num-1; 
       iremove = iremove+1;
     }
    }
#endif
if(verbose>=2)
  printf("Removing %d end atoms and Continue extension at fit: at iatom = %d, ires = %d\n",iremove,iatom,ires);
  ext.length = ext.length-iremove;  //this is the bone pt that equal nterm Ca
//  ext.length = ext.length+1; //NO! - I think have done this previouslyNew 11jul00 - add last pt = nterm Ca.
   for(i=1;i<=extmax;i++){ // extend out extmax residues
      for(ii= iatom;ii<= na[0];ii++){
         newbone = ext.main[imain].num; // the new bonept num
         ext.main[imain].num = newbone+1; // new total number
         if(ext.main[imain].num>MAXMAINPTS-2){
            printf("ERROR: too many mainpts in extendmainonfit - quitting\n");
            geo.quit = 1; // indicates error found
            goto quitrout;
         }
         ext.main[imain].mainpt[newbone].mtype = 0;
         ext.main[imain].mainpt[newbone].stype = 0;
         ext.main[imain].mainpt[newbone].ctype[0] = 0;
         for(k=0;k<=2;k++)
            ext.main[imain].mainpt[newbone].pos[k] =
                   fit[ifit].pos[ires][ na[ii] ][k];
      }
      ires = ires+istep;
      iatom = 1; // start next res with first atom
   }
   quitrout:;
}   
    

void editmain(int *endmain)
/*
 MAJOR MODIFICATION OF 6nov00:    
 This was commented out on 6nov00 - because if e.g. metmain[1]=1 before entering, never get out of while()
         for(k=1;k<=j;k++){
            if(i== ext.metmain[k])
               goto nextmain;
         }
#endif
 MAJOR MODIFICATION OF 3nov00:    
  1)Find the highest number main that ends on metmain[j] 
     and returns the mains on path to metend
     ext.metmain[0] = number of mains
     arranged in reverse order:  metmain[1] is main the meets end
      and metmain[ ext.metmain[0] ] should = main #1
  2)Presumably - the highest numbered one should have branched earlier and
     should be shorter branch
  3) If no fit meets metmain[j] - then just trace back to beginning
       as done in previous version.
 Called if ext.metend>0 - edits mainbranch so only learve
     main path to other fit
  Procedure:
    1) Find mains that are intitially have metend>0
       Probably will be only one, but may be 2
          Choose the shortest as the correct one
    2) Go to start point of this main,mainnum, 
        If startpt = startpt of main #1, then simple delete other mand
        Else, find the main which has this startpt number,
              Cut out all higher number pts from this main
           ALSO (NEW) CUT OUT  SIDEs THAT START AT POINT > THAN END MAIN
              This requires copying higher num sides to lower num
            Go back to (2) and repeat. until startpt = pt of main #1
    3) Must also allow for case where there are two paths from 1 main to the metmain
       a) First run step (2); Find the direct path and cut off ends of mains not on direct path
       b) Then for each main on this path (ext.main[i].metend = ext.metend)
            Go through all points (after cut in (2) above) and see if there is
              a connect type point (ctype[0]=1). 
           If so, trace this connection back though the mains until it hits
             one of the connect mains.  If it hits it at a point that has not
             been cut in (2) - then must include all mains on this branch - and must
             cut parts of these mains that are not on the direct branch
         FOR NOW - ASSSUME ONLY ONE CONNECT IN THIS BRANCH, THUS, ONLY NEED
            TO TRACE EACH MAIN BACK TO ITS START POINT
        I think connect can only come from higher numbered main
       As go though all points, may find several connect branches
         keep track of the maximum number of points needed on each main
           in this branch as follows:
         For each potential connect branch -
             tmain[i] = -1 if main not on branch, else = pt num where branch starts
         Then if branch meets a metend main bone, copy tmain to 
             pmain[i][1,2,3,4=max]  = pt num of up to 4 times this main is in a branch
      Then, cut from maxnum to end of main bones - HOWEVER - it is possible one point on
         this region that plan to cut, contacts point on path - must check this
      
*/
{
   int ii,count,i,j,k,kk,kj,sidecount;
   int metstart,istart,imin,length,minlen;

   j=0; // counts mains on path to met fit
   *endmain = 0; // default value - indicates not yet assigned
   for(i=1;i<=ext.nmain;i++){
      if(ext.main[i].metend>0){
if(verbose>=2)
  printf(" Main[%d] metend %d\n",i,ext.main[i].metend);
         j=j+1;
         ext.metmain[j] =i;
      }
   }
   if(j>1){  // If more than one main that met fit, use the shortest
    if(verbose>=2)
       printf(" %d INITIAL MAINS with ext.metmain>0\n",j);
     minlen = 10000;
     for(i=1;i<=j;i++){ // find shortest and return value of length
        length = ext.main[ ext.metmain[i] ].startnum + ext.main[ ext.metmain[i] ].num;
        if(length < minlen){
           minlen = length;
           imin = i;
        }
     }
     for(i=1;i<=j;i++){ // set longer ones back to -1
       if(i != imin)
         ext.main[ ext.metmain[i] ].metend = -1;
     }
     ext.metend = ext.main[ ext.metmain[imin] ].metend;
     j=1; // just use the shortest one
     ext.metmain[j]= ext.metmain[imin];
     *endmain = ext.metmain[j]; // this is passed to extendmainonfit
   }
   else{
     *endmain = ext.metmain[1];  
     minlen = ext.main[ ext.metmain[1] ].startnum + ext.main[ ext.metmain[1] ].num;
   }
   ext.length = minlen;  // use global value
   if(verbose>=2)
     printf("minlen = %d\n",minlen);
   istart = ext.main[1].mainpt[0].num;
   j=1; // Assumed number of initial met mains - counts # of mains used
   metstart = ext.main[ ext.metmain[j] ].mainpt[0].num;
// NOTE - comment the following section out if want to find edited main on higest density branch
#if 1 // NEW 3nov00 - find the shortest path to metend
int imet,endnum,metnum;
   imet=ext.metmain[j];
   for(i=ext.nmain;i>=1;i--){// go through all mains, find highest number main that ends on metmain[j] 
      endnum = ext.main[i].mainpt[ ext.main[i].num-1 ].num;
      for(k=0;k<ext.main[imet].num-1;k++){
        metnum=ext.main[imet].mainpt[k].num;
        if(metnum==endnum) { // Then trace back on this main
          j= j+1;
          ext.metmain[j]=i;
          metstart = ext.main[ ext.metmain[j] ].mainpt[0].num; // Use beginning of this main
          ext.main[i].metend = ext.metend;    
          ext.main[imet].zeronum =k;
          goto foundend;
        }
      }//k
   }//i;
   metstart = ext.main[ ext.metmain[j] ].mainpt[0].num; // If reach here, no mains end on metmain, just use begnning pt
   foundend:;
#endif
   while(metstart!=istart){ // Find direct path of mains to to metmain[1]
      for(i=1;i<=ext.nmain;i++){ // Go through other mains, look for start pt of mainn[j]
#if 0 // This was commented out on 6nov00 - because if e.g. metmain[1]=1 before entering, never get out of while()
         if(i== ext.metmain[1])
               goto nextmain;
#endif
         count=0;
         for(k=0;k<ext.main[i].num;k++){
            count = count+1;
            if(ext.main[i].mainpt[k].num==metstart){
                ext.main[i].num = count; // terminate this main at this point
                sidecount =0;
                for(ii=1;ii<=ext.main[i].nside;ii++){// only keep side that start on point less than this end
                   if(ext.main[i].side[ii].mnum<=count){ // starts on point on edited main
                     sidecount= sidecount +1;
                     if(sidecount!=ii){ // then must copy from ii to sidecouunt
                        for(kk=0;kk<ext.main[i].side[ii].num;kk++){   
                           ext.main[i].side[sidecount].ipos[kk]=ext.main[i].side[ii].ipos[kk];
                           for(kj=0;kj<=2;kj++)
                               ext.main[i].side[sidecount].pos[kk][kj]=ext.main[i].side[ii].pos[kk][kj];
                        }
                        ext.main[i].side[sidecount].mnum=ext.main[i].side[ii].mnum;     
                        ext.main[i].side[sidecount].num= ext.main[i].side[ii].num;
                     }
                   }// if side OK
                }// ii
                ext.main[i].nside = sidecount; // use new value for side number
                j= j+1;
                ext.metmain[j]=i;
                metstart = ext.main[ ext.metmain[j] ].mainpt[0].num; 
                ext.main[i].metend = ext.metend;    
                goto foundmet;  
            } 
          }// k 
          nextmain:;
       }//i
       foundmet:;
   } // while
  int pmain[MAXMAIN][MAXPATH],tmain[MAXMAIN],iside;
  int cmain,imain,nmet,mnum,maxnum[MAXMAIN];
   for(i=0;i<=ext.nmain;i++)  
       pmain[i][0] = 0; //counts number of time main i in a connect branch
   nmet = j; // counts number of metend mains
   for(i=1;i<=nmet;i++){ // Go through all edited mains, 
     imain = ext.metmain[i];
     for(k=0;k<ext.main[imain].num;k++){ // each point on main
        if(ext.main[imain].mainpt[k].ctype[0]!=0){// a side connect to this main
           for(ii= 1;ii<=ext.nmain;ii++) // reset tmain
              tmain[ii] = -1;
           cmain = ext.main[imain].mainpt[k].ctype[1];// main chain num that side starts on 
           mnum = ext.main[imain].mainpt[k].ctype[2]; // main num. that side stats a
           repeat:;
           for(kk=1;kk<=nmet;kk++){ // test if cmain is part of metmain
             if(cmain == ext.metmain[kk]){ // If meets a metend
                if (mnum<= ext.main[cmain].num){// if connect point still on edited metmain
                   for(ii= 1;ii<=ext.nmain;ii++){ // add mains on path to pmain
                     if(tmain[ii]>=0){
                       if(pmain[ii][0]>= MAXPATH-1){
                          printf("WARNING: pmain[ii][0]>= MAXPATH\n");
                          goto quitpath;
                        }
                        pmain[ii][0] =pmain[ii][0]  +1; // add another path to pmain
                        pmain[ii][ pmain[ii][0] ] = tmain[ii];
                     }
                   }
                }// if successfull path
              } // if cmain == metend
            }//kk - New 24aug00 -Now look for good path for all fits - previously never allowed more than k=1 to be tried
            if( cmain ==1)
                  goto nextk;
            else{ //- New 24aug00 - If reach here Add main to tmain and Go back to bone that cmain came from and repeat
                 tmain[cmain]=mnum;
                 mnum = ext.main[cmain].mainnum; // main num. that side stats a
                 cmain= ext.main[cmain].startmain;
                 goto repeat;
             }
//          }
        } // if ctype
        nextk:;
      } // k
    }// i
    quitpath:;
    for(ii= 1;ii<=ext.nmain;ii++){ // Find max num and add this bone to metend list
       maxnum[ii] = -1;
       for(i=1;i<=pmain[ii][0];i++){  // Find max. of num for this main
          if(pmain[ii][i]>maxnum[ii]){
             maxnum[ii] = pmain[ii][i];
          }
          nmet = nmet+1;             
          ext.main[ii].metend = ext.metend;
       }
    }
    for(ii=1;ii<=ext.nmain;ii++){ // Before cutting from maxnum to end, check if double back on itself - Start with lowest numberand work back
       if(pmain[ii][0]>=1) {  //indicate new main to be added to list
          for(k= ext.main[ii].num-1;k>maxnum[ii];k--){ // check if points plan to cut are connected to lower number point on bone
             if(ext.main[ii].mainpt[k].stype>0){// a side branch begins at this point
               iside = ext.main[ii].mainpt[k].stype;
               if(ext.main[ii].side[iside].ctype[0]==1){ // this side connects to another main
                  cmain = ext.main[ii].side[iside].ctype[1];// main chain num that side starts on 
                  mnum = ext.main[ii].side[iside].ctype[2]; // main num. that side stats a
                  if(ext.main[cmain].metend >0){ // if bone connects to is part of main
                    if (mnum<= ext.main[cmain].num){// if connect point still on edited metmain - then end here  
                       ext.main[ii].num = k+1; //Note: main pts run from 0,1,2 .. num-1 terminate this main at this point
                       goto nextii;
                    }
                  } // metnend>0
               } // ctype
             }// stype
          } // k
          ext.main[ii].num = maxnum[ii]+1; //If reach here - use the original maxnum
      } // if pmain
      nextii:;
   } // ii
#if 0   // Check if there is a side that connects 2 mains that metend - assume that this side is shortest distance - cut out rest of main
// This works for simple connection e.g. shad, fit 8, bone 1.4, but bombs for more 
//    complicated, e.g.shad, fit 6, bone 1.0 or  endo fit 9, bone 1.0
   for(i=1;i<=nmet;i++){ // Go through all edited mains, 
     imain = ext.metmain[i];
     for(k=0;k<ext.main[imain].num;k++){ // each point on main
        if(ext.main[imain].mainpt[k].ctype[0]>0){// a side connect to this main 
          cmain = ext.main[imain].mainpt[k].ctype[1]; // the main it connects to
          mnum = ext.main[imain].mainpt[k].ctype[2]; // the main pt it connects to
          if(ext.main[cmain].metend == ext.metend){ // check if main met end   
             if(mnum<= ext.main[cmain].num-1){ //if mnum still part of edited main then cut lower numberedmain
                ext.main[imain].num = k+1; // cut rest of imain;
//                ext.main[cmain].zeronum = mnum-1; //THIS CAUSES TROUBLE  cut the beginning of cmain
             }  // k
          } // if metend
        }// if side
     } // k
   }// i      
#endif 
   ext.metmain[0] = nmet;    
}


void copytoonemainpt(int pt0,int main1,int pt1)
/* Copies mainpt from main[main1].mainpt[pt1] to mainone, pt0
  Does not copy mtype, stype or ctype since they are either not needed
    or will be different
*/
{
   int k;
   ext.mainone.mainpt[pt0].num = ext.main[main1].mainpt[pt1].num;
//printf("main1=%d  pt1=%d  num= %d\n",main1,pt0,ext.mainone.mainpt[pt0].num);
   for(k=0;k<=2;k++)
       ext.mainone.mainpt[pt0].pos[k] =
              ext.main[main1].mainpt[pt1].pos[k];
//printf("   pos[0]= %5.3f\n",ext.mainone.mainpt[pt0].pos[0]);
}

#if 0 // old version copies all pt to main[1]
void makeonemain(void)
/* 
Copies all the main pts on path to one main (=1) with no side pts, etc.
  I think this should be easy because after saving and reading edited mains,
  the main pts are just in simple order on main 1,2,3..
 First mainpts will be identical to those in main[1], so that start
   adding pts with main 2,3, etc.
*/
{
   int i,jbest,j;

   jbest = ext.main[1].num-1; // this is number of point in best main -adding after first main
   for(i=2;i<=ext.nmain;i++){ // each additional main
      for(j=0;j<ext.main[i].num;j++){ //  the main pts
         jbest = jbest+1;
         if(jbest> MAXMAINPTS-2) { // this should not happen
           printf("ERROR: jbest> MAXMAINPTS-2 in makeonemain\n");
           ext.metend = -1; // do not try to extend
           jbest = jbest -1;
           goto quitcopy;
         }
         copymainpt(1,jbest,i,j);
      }
   }
   quitcopy:;
   ext.nmain = 1;
   ext.main[1].num = jbest+1;
   ext.main[1].zeronum = 0;  // start number
   ext.main[1].startnum = 0;
   ext.main[1].nside = 0;  // no sides
   for(j=0;j<ext.main[1].num;j++){ // side side data etc to zero
      ext.main[1].mainpt[j].mtype = 0;
      ext.main[1].mainpt[j].stype = 0;
      ext.main[1].mainpt[j].ctype[0] = 0;
   }
   ext.main[1].mainpt[0].mtype = 0; // Need this for first point in maincadist
   ext.main[1].metend = ext.metend;
}
#endif

#if 1 // NEW 3july00 version copies all pt to mainone
void makeonemain(void)
/* 
Copies all the main pts on path to one main (=1) with no side pts, etc.
  I think this should be easy because after saving and reading edited mains,
  the main pts are just in simple order on main 1,2,3..
 First mainpts will be identical to those in main[1], so that start
   adding pts with main 2,3, etc.
*/
{
   int i,jbest,j;

   jbest = ext.main[1].num-1; // this is number of point in best main -adding after first main
   for(j=0;j<=jbest;j++){ // copy main[1] to mainone
      copytoonemainpt(j,1,j);
   }
   for(i=2;i<=ext.nmain;i++){ // each additional main
      for(j=0;j<ext.main[i].num;j++){ //  the main pts
         jbest = jbest+1;
         if(jbest> MAINONEPTS-2) { // this should not happen
           printf("ERROR: jbest> MAINONEPTS-2 in makeonemain\n");
           ext.metend = -1; // do not try to extend
           jbest = jbest -1;
           goto quitcopy;
         }
         copytoonemainpt(jbest,i,j);
      }
   }
   quitcopy:;
   ext.nmain = -150; // this is indicator have made one mainone
   ext.mainone.num = jbest+1; 
   ext.mainone.zeronum = 0;  // start number
   ext.mainone.startnum = 0;
   ext.mainone.nside = 0;  // no sides
   ext.mainone.mainpt[0].mtype = 0; // Need this for first point in maincadist
   ext.mainone.metend = ext.metend;
}
#endif
int testclose(struct griddata& den1,struct bonedata& bone1,
      struct geometry& geo,int nearnum)
/*
Need special test to see if initial pt is close
  to another Ca because if two Ca's are close, set
  may overlap.
   returns 0 if not close
   returns iset = the other fit, if is close
*/
{
    int i,ivec[3],icapos;
    float fvec[3],dist;

    if (verbose>=2)
       printf("In routine testclose\n");
    d1tod3(bone1.cmax,bone1.ptv[nearnum].pos,ivec);
    imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
    for(i=1;i<=fitnum[0];i++){ // set points for each fit - set side and end separately
       if(i!=ext.ifit){ // not itself
         if(geo.fixend ==2) // if extending from Cterm, must hit nterm end of fit
             icapos = 0;
         else // if extending from nterm, must hit cterm end (Ca) of fit
             icapos= fit[i].nres -1;
         dist = distv(fvec,fit[i].pos[icapos][1]); // Dist from cterm Ca to icon
         if(dist<=2.0){//If close enough to end Ca then has made correct contact
            ext.ifit = nearnum;
            ext.metend=i; 
            ext.nmain = 1; 
            ext.main[ext.nmain].num = 1;
            ext.main[ext.nmain].num = 0;
            ext.main[ext.nmain].zeronum = 0; // default start of main
            ext.main[ext.nmain].nside = 0; // default start of main
            ext.main[ext.nmain].metend = i; // default start of main
            ext.main[ext.nmain].mainnum = 1;
            ext.main[ext.nmain].startmain = 0;
            ext.main[ext.nmain].startnum = 0;// startnum of this main
            addmainpt(den1,geo,bone1,ext.main[ext.nmain],ext.nmain,
                nearnum,0,0);
            return i;
         }
       } // not itself
    }  // if iset
    return 0;  // defaut - not close enough to another Ca 
}


void editbones(struct griddata& den1,float pos1[3],float pos2[3],
      struct geometry& geo,struct bonedata& bone1)
/* Routine for editing bones for extenstion:
 NOTE: Only cterm extension has been written:
  Distinguish between cterm and nterm exttension by geo.fixend:
         if(geo.fixend ==2) // if extending from Cterm, must hit nterm end of fit
                    icapos = 0;
             else // if extending from nterm, must hit cterm end (Ca) of fit
                    icapos= fit[iset].nres -1;

  For cterm extension on input: pos1 = N, pos2 = CA  or C and Ca residues
  Note: bone2 not used now    (bone1 = default (1.4) bone; bone2 = 0.8 bones
  Edited bones are stored in ext (see extendfit.h for definition
-NOTE  DO NOT ALLOW FOR 2 MAIN CHAINS STARTING AT FIRST POINT
*/
{
   int i,k,j,jj,ii,nearnum,ncon,ncon2,icon2,icon,ivec[3];
   float neardist;
   float dist1,dist2[5];
   float fvec[3],nvec[3];
   int mcon[4],sidenum[4],tpos[4][MAXSIDEPTS],iposm[4],nposm[4],metoneend[4];
   int reachedlim,startnum,endmain;
   int testmetend;
   float aveden1[4],aveden2[4];
   char editbone[300];
   float maxdist1,maxdist2;
   int imin1,imin2;
//   float ang;

   if (verbose>=2)
       printf("In routine editbones\n");
   reachedlim = 0;
   testmetend = 1; // This means that when saving bones, sort and save only those on path to fit
   for(j=0;j<4;j++)
     metoneend[j] = -1;
   tpos[0][0]=0; // Default value - should not be needed
   nearestbonept(den1,pos2,geo,bone1,&nearnum,&neardist);
#if 1 // This was only added on 6nov00 - but should obviously have been used before
   if(neardist>2.0){
      if(verbose>=1)
        printf("Nearest point dist = %5.3f - too far, quitting\n",neardist);
      ext.metend = -1;
      goto quitrout;
   }
#endif
   strcpy(editbone,maidfile.editbones);
   for(j=0;j<MAXMAIN;j++)
      ext.main[j].metend = -1; // default value - has not met other fit
   if(testclose(den1,bone1,geo,nearnum)>0) // Quit with just one bone pt
      goto onept;
   ncon = bone1.ptv[nearnum].num; // number of connections
   if(ncon<=2){ // If not branch - see if there is branch point conneected to nearnum
     for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       if(ncon2>2){
          nearnum = icon;
          goto newnear;
       }
     }
   }// ncon
   newnear:;
   if(testclose(den1,bone1,geo,nearnum)>0) // Quit with just one bone pt
      goto onept;
   d1tow(den1,bone1.cmax,bone1.ptv[nearnum].pos,nvec);// convert pos to world fvec[3]
   dist1 = distv(pos1,nvec); // distance from C to nearnum 
   ext.metend = -1;  // Indicates have not yet met end of another fit      
   ext.nmain = 0; // Number of main traces
   ext.main[0].num=1; // This is defined only so that startnewmain works OK
   ext.main[0].startnum =0; // initializes totnum
   ncon = bone1.ptv[nearnum].num; // number of connections
if(verbose>=2)
 printf("ncon = %d  nearnum = %d\n",ncon,nearnum);
   if(ncon >4){ 
      printf(" ERROR: ncon = %d > 4  setting = 4\n",ncon);
      ncon = 4;
   }
   imin1 = imin2 = -1; // set to values that will not be used unless repeat
findjjagain:;  // Come back and repeat if jj>2
   jj=0; // counts forward branches
   for(j=0;j<ncon;j++){
       icon = bone1.ptv[nearnum].con[j]; // the pt num of the connect point
       ncon2 = bone1.ptv[icon].num; // number of connections
       if(ncon2==2){ // If simple point - go to next point before checking dist
         for(ii=0;ii< ncon2;ii++){
           icon2 = bone1.ptv[icon].con[ii]; // the pt num of the connect point
           if(icon2!= nearnum){ // do not go backward
              goto testdist;
           }
         }//ii
       }// ncon2
       icon2= icon; // If ncon2 !=2, then test dist with first point
       testdist:;
       d1tow(den1,bone1.cmax,bone1.ptv[icon2].pos,fvec);// convert pos to world fvec[3]
       dist2[j] = distv(pos1,fvec); // distance from C to connect point
       if(dist2[j]>(dist1-1.0)){ //New - allow some slack for perpendicualar branch if connect point moving away from C
           if( (j==imin1)|| (j==imin2) ) // if repeat, do not use these
              goto nextj;
           iposm[jj]= nearnum;
           nposm[jj]= icon;
           tracetobrrec(den1,geo,bone1, &iposm[jj],&nposm[jj],&mcon[jj],
              &sidenum[jj],tpos[jj],ext.main[0].startnum,
              &metoneend[jj],&aveden1[jj],&aveden2[jj]);
//           ext.main[0].metend=metoneend[jj];
           jj = jj+1;
#if 0// 16dec99 - comment this out 
           if(ext.metend>0){ // Then met end, just use this one branch
              if(jj==2){ // copy 1 to jj=0
                iposm[0]=iposm[1];
                nposm[0]=nposm[1];
                mcon[0]=mcon[1];
                for(ii=1;ii<=sidenum[1];ii++)
                   tpos[0][ii]=tpos[1][ii];
                sidenum[0]=sidenum[1];
              }
              jj=1;
              ext.metend = -1; // must resetn because startnewmain call tracetobr again 
              goto quitloop;
          } // ext.metend>0 
#endif
      } //dist2>dist1
      nextj:;
   } // j
   quitloop:;
   if(jj>2){ // find the two with largest dist2 and repeat above
      maxdist1=maxdist2 = 10000.0;
      for(i=0;i<jj;i++){// find shortest
         if(dist2[i]<maxdist1){
              imin1 = i;
              maxdist1 = dist2[i];
         }
      }
      if(jj==4){ // find next shortest
         for(i=0;i<jj;i++){
            if(i!=imin1){
               if(dist2[i]<maxdist2){
                  imin2 = i;
                  maxdist2 = dist2[i];
               }
             }
         }
      } // jj=4;
      if(jj>4)
         printf("BIG BIG ERROR:  jj=%d >4\n",jj);
      goto findjjagain;
   }     
   if(jj==0)
     goto quitrout;
   d1tod3(bone1.cmax,bone1.ptv[ tpos[0][0] ].pos,ivec);
   geo.set[ivec[0]][ivec[1]][ivec[2]] = -1; // set this point
   startnum = 0; // count points from start
   if(jj>2){ // assume max of 2 branches  (jj = 0,1)
      printf(" ERROR: jj = %d > 2  setting = 2\n",jj);
      jj=2;
   }
   if(jj==1){ // just start main[1] (main[0] is dummy)
//        ext.main[1].mainpt[0].mtype = 1; //Indicates new main starts at this oint
        startnewmain(&reachedlim,den1,geo,bone1,ext.main[0],0,
             tpos[jj-1][0],tpos[jj-1][1],startnum,1); // 1 = origmainnum
   }
   else if(jj==2){
      ext.nmain=1;
      ext.main[1].num = 0; // First point on main = 0
      addmainpt(den1,geo,bone1,ext.main[1],1,nearnum,0,0);// add first point - set mtype = ext.nmain
      ext.main[1].num = 1;
      ext.main[1].nside=0;
      ext.main[1].startnum =0; // initializes totnum
      add2branchs(&reachedlim,den1,geo,bone1,ext.main[1],1,tpos,
             mcon,sidenum,aveden1,aveden2,metoneend);

   }
   quitrout:; 
//ext.metend = -1;
onept:; // jump to here if if just meet end at first pt.
ext.length = 100; // for degubbing
#if 0 // only used for debugging - save all bones 
   saveeditbones(editbone,-1); //this writes only the main chains that are on path to met end      
   readeditbones(editbone); // read the saved edited bones
exit(0);
#endif
#if 1 // for debugging - do not edit bones or make onemain    
   ext.length = -100; // Value if does not meet end
   if(ext.metend>0){ // then some main has met other fit - edit mains
     editmain(&endmain); // return endmain = main num that contacted fit
//     increasemainden(den1,bone1);  //NOT USED  move main to neighbor pt if it has higher den.
     if(geo.fixend==2){ // NEW 29aug00 - used to be done in extendfit.c++ extending at cterm end - set 0 point = term Ca
          for(k=0;k<=2;k++) 
              ext.main[1].mainpt[0].pos[k] = pos2[k];
     }
     else{ // extending at nterm
          for(k=0;k<=2;k++) // New - set first pt = Nterm Ca
            ext.main[1].mainpt[0].pos[k] = pos2[k];
     }
     extendmainonfit(geo,endmain);
     if(geo.quit==1){
        printf("ERROR FOUND IN  extendmainonfit() - CANNOT MAKE ONE MAIN\n");
        ext.metend = -1;
        goto skipsave;
     }
#if 0 // 9 nov01:  This should not be here - must first convert to onemain using makeonemain before sav
     saveeditbones(editbone,testmetend); //this writes only the main chains that are on path to met end      
     readeditbones(editbone); // read the saved edited bones
#endif
     if(ext.length< MAINONEPTS -2 ){
        makeonemain();  // write all mainpts to one main
        saveeditbones(editbone,testmetend); // just save the mainone main -metmain not used     
        readeditbones(editbone); // read the saved edited bones
     }
     else{ 
         printf(" ERROR: CANNOT MAKE ONE MAIN BECAUSE ext.length (=%d) is too long\n",ext.length);
         ext.metend = -1;
     }
  } // ext.metend>0
  skipsave:;
  #if 0 //for debugging  SET = 1 IF WANT TO SEE BONES EVEN IF DO NOT MEET FIT
   saveeditbones(editbone,testmetend); //this writes only the main chains that are on path to met end      
   readeditbones(editbone); // read the saved edited bones
  #endif
#endif
#if 0 // Debug, check if main and side start in right place
int i;
for(i=1;i<=ext.nmain;i++){
     for(j=0;j<ext.main[i].num;j++){ // draw the main pts
        if(ext.main[i].mainpt[j].mtype!=0){
           printf("main %d starts at main %d point %d ipos %d\n",
              ext.main[i].mainpt[j].mtype,i,j,ext.main[i].mainpt[j].num);
        } 
        if(ext.main[i].mainpt[j].stype!=0){
           printf("side %d starts at main %d point %d ipos %d\n",
              ext.main[i].mainpt[j].stype,i,j,ext.main[i].mainpt[j].num);
        } 
        if(ext.main[i].mainpt[j].ctype[0]!=0){
           printf("connect from main %d side %d starts at main %d point %d ipos %d\n",
              ext.main[i].mainpt[j].ctype[1],ext.main[i].mainpt[j].ctype[2],
               i,j,ext.main[i].mainpt[j].num);
        }
     }//j
}// i 
#endif
}


#if 1 

void testconnect(struct bonedata& bone1,int spt,int npt,int pt2,int maxdist,
         int condist,int *iresult)
/* Recursive function - 
  Trace path in direction of spt to npt until either find branch or
    pt = pt2 when return condist = number of point up to branch.
  Finally returns for 3 cases:
    1) reaches pt2  condist = total number of pts from origin 
    2) reaches end of trace
    3) condist =  maximum number of pts from origin
*/
{
   int i,ncon,ppt,icon0[20],count,icon,i1;
  
   ncon = bone1.ptv[npt].num; // number of connections to pt1
   ppt = spt;
   condist = condist +1;
   while(ncon==2){
      for(i=0;i<2;i++){
         icon =  bone1.ptv[npt].con[i];
         if(icon!=ppt){
            ppt = npt;
            npt = icon;
            ncon = bone1.ptv[icon].num; // number of connections to pt1
            condist = condist +1;
            if(npt==pt2){
               *iresult = condist;
                goto quitrout;
            }     
            goto foundcon;
         } 
      }
      foundcon:;
   }
   if( (condist>maxdist)||(ncon==1) )
      *iresult = -1;
   else{ // reached another branch - call testconnect again
      count = 0;
      for(i=0;i<ncon;i++){
        i1 =  bone1.ptv[icon].con[i];
        if(i1!= ppt){
          if(i1==pt2){
             *iresult = condist+1;
             goto quitrout;
          }
          count = count+1;
          icon0[count-1] = i1;
          
        }
      }
      for(i=0;i<count;i++){
         testconnect(bone1,icon,icon0[i],pt2,maxdist,condist,iresult);
         if(*iresult>0)
           goto quitrout;
      }
   }
   quitrout:;
}

void starttestconnect(struct bonedata& bone1,int pt1,int pt2,int maxdist,
        int *iresult)
/* Routine to test whether pt1 and pt2 are connected within
   and max separation of cmax pts, return condist = -1 if not 
   connnected or condist = number of points between pt1 and pt2
 Procedure:  Find all pts connected to pt1 and call
    testconnect foreach pair of pts (go in all directions)
    testconnect is recursive function
*/
{
   int ncon,i,icon0[20],condist;

   condist = 0;
   ncon = bone1.ptv[pt1].num; // number of connections to pt1
   for(i=0;i<ncon;i++)
      icon0[i] =  bone1.ptv[pt1].con[i];
   for(i=0;i<ncon;i++){
      testconnect(bone1,pt1,icon0[i],pt2,maxdist,condist,iresult);
      if(*iresult>0)
        goto quitrout;
   }
   *iresult = -1; // if reach here, never found connection
   quitrout:;
}
#endif
