/* testfit.inc  
  Routines used to check fit.
  Just used for debugging - try to find procedures for
   cutting bad fits.
*/

#if 0// find minimum distance of far side atoms from mainchain
int i,ires;
float mindist;
   mindist = 1000.0;
   for(i=1;i<=fitnum[0];i++){
     printf("Testing fit %d\n",i);
     for(ires=1;ires<=fit[i].nres-2;ires++)      
        findside_maindistfit(fit[i],ires,&mindist);
   }
   printf("mindist = %5.3f\n",mindist);
#endif
#if 0 // find distance between long sides and O for same residue
   int i,j,k;
   float mag;
   for(i=1;i<=fitnum[0];i++)
     for(j=1;j<=fit[i].nres-2;j++)
       if(fit[i].itype[j]==2)
        if(fit[i].nside[j]>=3)
           for(k=3;k<=fit[i].nside[j];k++){
              mag = distv(fit[i].pos[j][4+k],fit[i].pos[j][3]);
              if(mag<=2.5)
                 printf(" For fit %d res %d side-0 dist = %5.3f\n",i,j,mag);
           }
#endif


#if 0 // Find min dist of atoms from all bones
   verbose=2;
//   savefitcheckbone(map[1],geo,bone[1]); //N0 LONGER USED -Only save fits that have good fit to bones
   bonedist(map[1],geo,bone[1]);//N0 LONGER USED - (part of savefitcheckbone)
#endif
    
#if 0 // routine to look for bad CA density If 2 or more bad, it returns -1
   int ii;
   for(ii=1;ii<=fitnum[0];ii++){
     copyfittogeo(geo,fit[ii]);
     printf("For fit %d: testden = %d\n",ii,testden(map[1],geo) );
   }
#endif

#if 0 //NOTE this not used any more - use checkoneendbranch - debug the  checkendbranch routine
  int i;  
   checkendbranch(map[1],geo,bone[1],1); //Checks if end on branched main bone
  for(i=0;i<=2;i++){ // This overwrites colors in bonefile
              bone[1].mcolor[i]=magentav[i];
              bone[1].scolor[i]=cyanv[i];
  }
#endif

#if 0  // debug testconnect routine
   int pt1,pt2,iresult,maxdist;
   maxdist = 25;
   pt1 = 11184;
   pt2 = 5111;
   starttestconnect(bone[1],pt1,pt2,maxdist,&iresult);
   printf("At end of startestconnect, condist = %d\n",iresult);
#endif

#if 0 // check the routine that check Ca dist at helix end
int ii,ifit;
float dist;
   ifit = 1;
   for(ii=1;ii<=fitnum[0];ii++){
//   for(ii=ifit;ii<=ifit;ii++){
     copyfittogeo(geo,fit[ii]);
     if( (fit[ii].itype[2] ==1)&& ((geo.nend[1]-geo.nend[0])>3) ){//If helix - check Ca[i]-Ca[i+4] dist
         dist=distv(geo.pos[geo.nend[0]][1],geo.pos[geo.nend[0]+4][1]);
         printf("  fit %d  Ca[1]-Ca[5] dist = %5.3f:\n",ii,dist);
         if(dist < geo.ca1_4min){
            printf("IN checkres bad  because Ca[1]-Ca[5] dist = %5.3f\n",dist); 
//            return 0;
         }
         dist=distv(geo.pos[geo.nend[1]-4][1],geo.pos[geo.nend[1]][1]);
         printf("  fit %d  Ca[end-4]-Ca[end] dist = %5.3f:\n",ii,dist);
         if(dist < geo.ca1_4min){
            printf("IN checkres bad  because Ca[end-4]-Ca[end] dist = %5.3f\n",dist); 
//            return 0;
         }
    } // helix constraint 
   }
#endif

#if 0 // debug the checkextendenfit, checkoneendbranch or bonedistgeo routine
int ii,*quittrace,origsavlen,ifit,origverbose,endvalue,resbad[200];
   origsavlen = geo.savlen;
//   geo.savlen = geo.savlen-1;
   origverbose = verbose;
//   verbose = 2;
   for(ii=1;ii<=fitnum[0];ii++){
#if 0
   ifit = 36;
    for(ii=ifit;ii<=ifit;ii++){
#endif
//if(fit[ii].nres<=5){ // only check short fits
      printf(" Checking fit %d    nres = %d:\n",ii,fit[ii].nres);
      copyfittogeo(geo,fit[ii]);
      geo.itype = geo.rtype[geo.nend[0]+1];//NOTE assume all res. same type
// if(geo.itype==2)
     *quittrace=0; // default
     bonedistgeo(map[1],geo,bone[1],quittrace);
     if(geo.frcon<0.999)
        bonedistgeo2(map[1],geo,bone[1],quittrace);
     if(*quittrace==2)
        printf("BAD FIT RETURNED FROM BONEDISTGEO\n");
//      checkextendedfit(1,geo,map[1],resbad,quittrace); //Recursive call
//      checkoneendbranch(map[1],geo,bone[1],quittrace);
//      checkterm(map[1],geo,bone[1],&endvalue);
//}//if(fit[ii].nres<=5){ // only check short fits
   }
   geo.savlen=origsavlen;
   verbose  = origverbose;
#endif

#if 0
int i,nii[2],nbad,nres[100],quitcond;
float origrfor2,lastden;
   for(i=1;i<=fitnum[0];i++)
      nres[i] = geo.nend[1]-geo.nend[0];
   checkend(map[1],geo,bone[1],geo.savlen,1); // When all finished, check newfits for bad ends and remove
   for(i=1;i<=fitnum[0];i++){ // Refine each fit that has been shorted before saving
    if( (fit[i].nres-1)>=(geo.savlen-1) ){ // if long enough 
     nii[0]=fit[i].val1;
      nii[1]=fit[i].val2;
      copyfittogeo(geo,fit[i]);
      if( (geo.nend[1]-geo.nend[0]) < nres[i]){ // ends have been cut, refine
          if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
               tordata.adata[1] = 2; 
          geo.fixend=0; // Refine all
          quitcond = 30;
          origrfor2 = tordata.rfor;
          tordata.rfor = 0;
          printf("Refining shortened fit %d before saving  tordata.rfor = %5.3f\n",i,tordata.rfor);
          tordynbat(quitcond,map[1],geo,tordata,&lastden,&nbad);// "1" means use N term density as criteria for q;uiting
          tordata.rfor = origrfor2;
      } // if shortened
      sidenumreg(geo.nend[0],geo.nend[1],map[1],geo); // Not used anymore
      copygeotofit(geo,fit[i],nii);
    }
   }// i  
    savefitcheckbone(map[1],geo,bone[1]); //Only save fits that have good fit to bones
#endif

#if 0 // compare nterm end O versus Ca den
   int ii;
   float oden,caden;
   for(ii=1;ii<=fitnum[0];ii++){
     copyfittogeo(geo,fit[ii]);
     caden= aveatom(map[1],geo,1,geo.nend[0],geo.pos);
     oden= aveatom(map[1],geo,3,geo.nend[0],geo.pos);
     printf("fit %d:  caden = %5.3f  oden = %5.3f  caden/oden = %5.3f\n",
          ii,caden,oden,caden/oden);
   }
#endif

#if 0 // Look at ratio of Cb/Cb2
float avecb2,aveco,finalden;
int i,ii,resbad[NAFIT],nbad,iremove;
   for(ii=1;ii<=fitnum[0];ii++){
     if((fit[ii].itype[2]==2)&&( (fit[ii].nres-1)<=(geo.savlen+1) )){ // if short sheet
       copyfittogeo(geo,fit[ii]);
       aveden(map[1],geo,geo.pos,&finalden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
       printf("Fit %d:  avecb/avecb2 = %5.3f\n",ii,geo.avecb/avecb2);
    }
   }
#endif

#if 0 // compare aveden and avemiden
    avemaxbatch(map[1],geo);
#endif 

#if 0 // run some kind of averaging
float avecb2,aveco,finalden;
int i,ii,resbad[NAFIT],nbad,iremove;
#if 0 // Used after  30may00 - smaller allowed range seems to work better  - used only while tracing - then go back to original
float origgeohelix[10];
for(i=0;i<=7;i++)
   origgeohelix[i] = geo.sheet[1][i]; 
   geo.sheet[1][0]= -DEGTORAD*110.0; /*phi min*/
   geo.sheet[1][1]= -DEGTORAD*50.0;  /* phi max*/
   geo.sheet[1][2]= -DEGTORAD*45.0; /*psi min*/
   geo.sheet[1][3]= DEGTORAD*10.0; /*psi max*/
   geo.sheet[1][7] = DEGTORAD*30.0; //New 30may00 - try smaller range +/- range 
#endif
   for(ii=1;ii<=fitnum[0];ii++){
     copyfittogeo(geo,fit[ii]);
#if 0
       iremove = checkendden(map[1],geo,bone[1],geo.pos,ii);
       printf("iremove = %d\n",iremove);
#endif 
#if 0 // Use aveden
    if(fit[ii].itype[2] ==2){ // if sheet
       aveden(map[1],geo,geo.pos,&finalden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
       printf("In  fit %d:  aveC/aveN = %5.3f\n",ii,geo.aveatom[2]/geo.aveatom[4]);
#if 0
       printf("In  fit %d: Length =%d   finalden = %5.3f avecb = %5.3f avecb2 = %5.3f aveco = %5.3f nbad = %d npoor = %d obad = %d\n",
         ii,geo.nend[1]-geo.nend[0],finalden,geo.avecb,avecb2,aveco,nbad,geo.npoor,geo.obadn);
       for(i=1;i<=6;i++)
         printf("  ave[%d]= %5.3f",i,geo.aveatom[i]);
       printf("\n");
#endif
#if 0 // print respoor (=resbad) and resbad (=res0)
       if((float)geo.npoor/(geo.nend[1]-geo.nend[0]) > 0.39)
           printf("Fit %d:  npoor ratio = %5.3f > 0.39\n",ii,(float)geo.npoor/(geo.nend[1]-geo.nend[0]));
       printf("Fit %d:  coave/(ca + n ave) = %5.3f\n",ii, 2.0*geo.aveatom[3]/(geo.aveatom[1]+geo.aveatom[4]));
#endif
#if 0 // print respoor (=resbad) and resbad (=res0)
      for(i=geo.nend[0];i<geo.nend[1];i++)
         printf("respoor[%d] = %d resbad=%d \n",
              i,resbad[i],geo.res0[i]);
#endif // print respoor
     } // if sheet
#endif // use aveden

#if 0 // test avedenfit0
float aveden;
   avedenfit0(map[1],geo,fit[13],&aveden);
printf("Fit 13 aveden = %5.3f\n",aveden);
#endif


#if 0 // find and test phi/psi angles
float phi,psi,psimin,psimax,phimin,phimax,psim;
int rtype;
    printf("Fit %d phi/psi:",ii);
    for(i=geo.nend[0];i<=geo.nend[1]-2;i++){
       rtype = fit[ii].itype[i-geo.nend[0]];
       if(rtype ==2){
          findphipsi(geo.pos,i,&phi,&psi);
          printf("  %d %5.3f %5.3f",i-geo.nend[0]+1,phi*180.0/PI,psi*180.0/PI);

         psim = geo.sheet[1][5]*phi + geo.sheet[1][6]; // correalted ideal psi
         psimin = psim -geo.sheet[1][7];// correlated psimin
         psimax = psim +geo.sheet[1][7];
         if(psimin < geo.sheet[1][2]) // absolute allowed min.
             psimin = geo.sheet[1][2];
         if(psimax > geo.sheet[1][3]) // absolute allowed max.
             psimax = geo.sheet[1][3];

// printf("\n i= %d: PHI= %6.3f/PSI=%6.3f psim =%5.3f  psimin = %5.3f psimax = %5.3f)\n",
//                i-geo.nend[0]+1,phi*180.0/PI,psi*180.0/PI,psim*180.0/PI,psimin*180.0/PI,psimax*180.0/PI);
          if(checkangrange(geo,phi,psi,rtype,&phimin,&phimax,&psimin,&psimax)==0){
//if(verbose==1)
             printf("\n ERROR IN RES. # %d: PHI= %6.3f/PSI=%6.3f NOT IN RANGE (type = %d: %5.3f %5.3f %5.3f %5.3f)\n",
                i-geo.nend[0]+1,phi*180.0/PI,psi*180.0/PI,rtype,phimin*180.0/PI,phimax*180.0/PI,psimin*180.0/PI,psimax*180.0/PI);
           }
           if(psi*180.0/PI < -55.0)
               printf("SMALLPSI");
        } // rtype
     }// i
     printf("\n");
#if 0
for(i=0;i<=7;i++)
   geo.sheet[1][i]=origgeohelix[i]; 
#endif
#endif // test phi/psi
#if 0 // use aveden for end region
   avedenreg(map[1],geo,geo.pos,&finalden,&avecb2,&nbad,
                  geo.nend[1]-1,geo.nend[1]); // Find average density for last two additions
   printf("In  fit %d: Length =%d   endfinalden = %5.3f  nbad = %d\n",
         ii,geo.nend[1]-geo.nend[0],finalden,nbad);
#endif // ave for end region
   }//ii
#endif //some kind of averaging


#if 0 // find angles between successive Ca-Cb vectors and between co vectors if helix or sheet (select below)
   float sideang,coang,avecorot;
int rtype,countco;
int i,ii;
   for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==2){
//     if(rtype==1){
#if 0//   print sideang rotation angle
       printf("Side rotation angles for fit %d type = %d\n",ii,rtype);
       for(i=1;i<=fit[ii].nres-3;i++){
           sideang =siderot(fit[ii].pos[i][1],fit[ii].pos[i+1][1],
                        fit[ii].pos[i][5],fit[ii].pos[i+1][5]);
           coang =corot(fit[ii].pos[i][2],fit[ii].pos[i+1][2],
                        fit[ii].pos[i][3],fit[ii].pos[i+1][3]);
           printf("Residue %d: sideang (%5.3f) coang (%5.3f)\n",i,sideang,coang);
           if( (sideang>215.0)&&(coang<115) )
             printf("Fit %d Res %d:Both bad: sideang =%5.3f>%5.3f  coang =%5.3f< %5.3f\n",
                      ii,i,sideang,215.0,coang,115.0);
       }
       printf("\n");
#endif
#if 0  // print Co rotation angle for i to i+2 CO
       printf(" Fit %d  max that can be tested = %d:\n",ii,fit[ii].nres-4);
       for(i=0;i<=fit[ii].nres-4;i++){
           coang =corot(fit[ii].pos[i][2],fit[ii].pos[i+2][2],
                        fit[ii].pos[i][3],fit[ii].pos[i+2][3]);
            printf("Res %d co2ang =  %5.3f",i,coang);
           if(coang<25.0)
              printf("BAD Res %d co2ang =  %5.3f",i,coang);
       }
       printf("\n");
#endif
#if 0  // print Co rotation angle
       printf("CO rotation angles for fit %d type = %d\n",ii,rtype);
       countco=0;
       avecorot=0;
       for(i=0;i<=fit[ii].nres-3;i++){
           countco=countco+1;
           coang =corot(fit[ii].pos[i][2],fit[ii].pos[i+1][2],
                        fit[ii].pos[i][3],fit[ii].pos[i+1][3]);
           avecorot = avecorot+coang;
//           if(coang <= 115.0)
//              printf("  BAD:");
           printf(" %5.3f",coang);
              
       }
       avecorot = avecorot/countco;
       printf(" aveco = %5.3f",avecorot);
       printf("\n");
#endif
#if 0 // print Co_side rotation angle  -- ang between Ca-Cb (i)and  co(i) and 
       printf(" Ca-Cb (i)and  co(i) rotation angles for fit %d type = %d\n",ii,rtype);
       countco=0;
       avecorot=0; // ang between co(i) and Ca-Cb (i+1)
       for(i=0;i<=fit[ii].nres-2;i++){
           countco=countco+1;
           coang =corot(fit[ii].pos[i][1],fit[ii].pos[i][2],
                        fit[ii].pos[i][5],fit[ii].pos[i][3]);
           avecorot = avecorot+coang;
//           if(coang <= 65.0)
//              printf("For fit %d res %d (fit.nres = %d): bad co_side ang = %5.3f < %5.3f\n",
//                          ii,i,fit[ii].nres,coang,65.0);
           printf(" %5.3f",coang);
              
       }
       avecorot = avecorot/countco;
//       printf(" aveco = %5.3f",avecorot);
       printf("\n");
#endif
#if 0  // print Co_side rotation angle  -- ang between co(i) and Ca-Cb (i+1)
//       printf("CO_side rotation angles for fit %d type = %d\n",ii,rtype);
       countco=0;
       avecorot=0; // ang between co(i) and Ca-Cb (i+1)
       for(i=0;i<=fit[ii].nres-3;i++){
           countco=countco+1;
           coang =corot(fit[ii].pos[i][2],fit[ii].pos[i+1][1],
                        fit[ii].pos[i][3],fit[ii].pos[i+1][5]);
           avecorot = avecorot+coang;
           if(coang <= 65.0)
              printf("For fit %d res %d (fit.nres = %d): bad co_side ang = %5.3f < %5.3f\n",
                          ii,i,fit[ii].nres,coang,65.0);
//           printf(" %5.3f",coang);
              
       }
       avecorot = avecorot/countco;
//       printf(" aveco = %5.3f",avecorot);
//       printf("\n");
#endif
     }  // rtype
   }//ii        
#endif // find co and side rot angles

#if 0  // Find angle between successive Ca-Ca vectors
int ii,rtype;
//   verbose=2;
   for(ii=1;ii<=fitnum[0];ii++){
     printf("for fit %d\n",ii);
     rtype = fit[ii].itype[2];
//     if(rtype==2){
     if(rtype==1){
      copyfittogeo(geo,fit[ii]);
        if(ca_caang3(geo)==1)
           printf(" BAD ca_aange FIT FOR FIT %d\n",ii);
     }
   }
#endif

#if 0  // Find distance plus/minus one Ca
int ii,rtype;
   for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==1){
        copyfittogeo(geo,fit[ii]);
        if(cadel(geo)==1)
            printf(" BAD cadel FIT FOR FIT %d\n",ii);
     }
   }
#endif

#if 0  // find angle between CO vector and axis
int rtype;
int i,ii;
float coaxisang,nexttolastang,lastang;
   for(ii=1;ii<=fitnum[0];ii++){
     copyfittogeo(geo,fit[ii]);
     if(geo.nend[1]-geo.nend[0]<=3){
        printf(" fit %d is too short to check C0 - axis angle\n",ii);
        goto skipcoaxis;
     }
     rtype = fit[ii].itype[1];
     if(rtype==1){
//     if(rtype==2){
       printf("CO versus axis angles for fit %d type = %d\n",ii,rtype);
       for(i=geo.nend[0];i<=geo.nend[1]-1;i++){
            coaxisang = co_axisang(geo,i);
            printf("  %5.3f",coaxisang);
            if(i==geo.nend[1]-2)
               nexttolastang = coaxisang;
            else if(i==geo.nend[1]-1)
               lastang = coaxisang;
       }
       printf("\n");
       if( (nexttolastang - lastang)>30.0)
          printf("FIT %d Bad last residue because nexttolastang - lastang =%5.3f\n",
                     ii,nexttolastang - lastang);
     }  // rtype
     skipcoaxis:;        
   }//ii
#endif

#if 0  // find distance from ca[i] to ca[i+3] for hellix
int rtype,i,ii;
   float c4dist;
   for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==1){
        printf("Distance from ca[i] to ca[i+3] for fit %d type = %d\n",ii,rtype);
        for(i=3;i<=fit[ii].nres-1;i++){
           c4dist = distv(fit[ii].pos[i][1],fit[ii].pos[i-3][1]);
           printf("  %5.3f",c4dist);
        }
        printf("\n");
     }  // rtype
   }//ii        
#endif


#if 0  // find distance from o[i] to ca[i+2] for hellix
int rtype,i,ii;
  float c4dist;
   for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==1){
        printf("Distance from o[i] to ca[i+2] for fit %d type = %d\n",ii,rtype);
        for(i=2;i<=fit[ii].nres-1;i++){
           c4dist = distv(fit[ii].pos[i][1],fit[ii].pos[i-2][3]);
           printf("  %5.3f",c4dist);
        }
        printf("\n");
     }  // rtype
   }//ii        
#endif


#if 0  // find distance from c[i] to c[i+1] for sheet
   float ccdist;
   int rtype,ii,i;
   for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==2){
        printf("Distance from c[i] to c[i+1] for fit %d type = %d\n",ii,rtype);
        for(i=1;i<=fit[ii].nres-2;i++){
           ccdist = distv(fit[ii].pos[i][2],fit[ii].pos[i-1][2]);
           printf("  %5.3f",ccdist);
        }
        printf("\n");
     }  // rtype
   }//ii        
#endif


#if 0  // find distance from o[i] to o[i+1] for sheet
   float oodist;
   int rtype,ii,i;
   for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==2){
        printf("Distance from o[i] to o[i+1] for fit %d type = %d\n",ii,rtype);
        for(i=1;i<=fit[ii].nres-2;i++){
           oodist = distv(fit[ii].pos[i][3],fit[ii].pos[i-1][3]);
           printf("  %5.3f",oodist);
        }
        printf("\n");
     }  // rtype
   }//ii        
#endif
 
#if 0 // run helixca_caang
  int ii,i,k,rtype;
  float ang,minang;
  
  minang = 1.70;
  for(ii=1;ii<=fitnum[0];ii++){
     rtype = fit[ii].itype[1];
     if(rtype==1){
       for(i=1;i<=fit[ii].nres-3;i++){       
         ang = helixca_caang(fit[ii].pos[i-1][1],fit[ii].pos[i][1],fit[ii].pos[i+2][1]);
         if(ang<minang)
            printf("ERROR:fit %d res %d helixca_caang =%5.3f < %5.3f\n",ii,i,ang,minang); 
       }
     } // rtype
   }//ii
#endif
