/* routines to find if assigned region has a non-crystallographic
   symmetry related region among the other fits.
TENTATIVE PROCEDURE:
   1) Only apply it to assigned regions that have been found by
        searching entire PDB sequence.  (the other regions are found
        extension from known regions).
   2)Then - starting with first Ca in Assigned region, find rotation and
        translation matrix that provides the best supperpositon for
       (say) Ca1-Ca3 of assigned to Ca1-Ca3 of uassgined.(Only compare sheets
        with sheets, and helixes with helices.
        This is done by use of quaternions (see notation in paper
                by A. L. Mackay (Acta Ctyst, 1984 A40, 165-166)
   3)Then using this translation- rotation matrix as starting point, use Powell
       minimization to find best denstity of the entire Assigned region.
       If this fits unknown region well, then assume that they are symmetry mates.

   4) Repeat step 2) for the other Ca's in assigned and unassigned.
1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in pos
*/
#include "drawwindow.h"
#include "math.h"
#include "trace.h"
#include "peptide.h"
#include "pdbtolev.h"
#include "symmate.h"
#include "assignseq.h"
#define MAXRES  20
#define MAXV   5*MAXRES+1 // Assume 4 atoms per res
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; // if = 1,print some stuff;  = 2 print a lot of stuff
extern int side[AA][NBRN][NPTS];  // defined globally (originally in pdbtolev.c++)  because needed in assignseg.c++ 
int maxside = 3; // maximum number of sidechain atoms rotated

void eigsrt(float d[],float v[4][4],int n)
/* Given the eigenvalues d[1,,,n] and eigenvectors v[1..n][1..n]- routine sorts
    eigenvalues into descending order, and rearanges the columns of v.
*/
{
    int k,j,i;
    float p;

    for(i=1;i<n;i++){
      p=d[k=i];
      for(j=i+1;j<=n;j++)
        if(d[j]>=p) p=d[k=j];
      if(k!=i){
        d[k]=d[i];
        d[i]=p;
        for(j=1;j<=n;j++){
          p=v[j][i];
          v[j][i]=v[j][k];
          v[j][k]=p;
        }
      }
    }
}


#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
	a[k][l]=h+s*(g-h*tau);

void jacobi(float a[4][4],int n,float d[],float v[4][4],int *nrot)
/* Computes all eigenvalues and eigenvectors of a real symmetric matrix 
    a[1..n][1..n]. On output, elements of a avove the diagonal are destroyed.
   d[1..n] returns the eigenvalues of a.  v[1..n][1..n] is a matrix whose 
   columns conatin, on output, the normalized eigenvectors of a.  nrot returns
   the number of Jacobi rotations that were required.
*/
{
	int j,iq,ip,i;
	float tresh,theta,tau,t,sm,s,h,g,c,*b,*z,*vector();

	b= new float[n+1];
	z= new float[n+1];
	for (ip=1;ip<=n;ip++) {
		for (iq=1;iq<=n;iq++) v[ip][iq]=0.0;
		v[ip][ip]=1.0;
	}
	for (ip=1;ip<=n;ip++) {
		b[ip]=d[ip]=a[ip][ip];
		z[ip]=0.0;
	}
	*nrot=0;
	for (i=1;i<=50;i++) {
		sm=0.0;
		for (ip=1;ip<=n-1;ip++) {
			for (iq=ip+1;iq<=n;iq++)
				sm += fabs(a[ip][iq]);
		}
		if (sm == 0.0) {
                        delete []z;
                        delete []b;
			return;
		}
		if (i < 4)
			tresh=0.2*sm/(n*n);
		else
			tresh=0.0;
		for (ip=1;ip<=n-1;ip++) {
			for (iq=ip+1;iq<=n;iq++) {
				g=100.0*fabs(a[ip][iq]);
				if (i > 4 && fabs(d[ip])+g == fabs(d[ip])
					&& fabs(d[iq])+g == fabs(d[iq]))
					a[ip][iq]=0.0;
				else if (fabs(a[ip][iq]) > tresh) {
					h=d[iq]-d[ip];
					if (fabs(h)+g == fabs(h))
						t=(a[ip][iq])/h;
					else {
						theta=0.5*h/(a[ip][iq]);
						t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
						if (theta < 0.0) t = -t;
					}
					c=1.0/sqrt(1+t*t);
					s=t*c;
					tau=s/(1.0+c);
					h=t*a[ip][iq];
					z[ip] -= h;
					z[iq] += h;
					d[ip] -= h;
					d[iq] += h;
					a[ip][iq]=0.0;
					for (j=1;j<=ip-1;j++) {
						ROTATE(a,j,ip,j,iq)
					}
					for (j=ip+1;j<=iq-1;j++) {
						ROTATE(a,ip,j,j,iq)
					}
					for (j=iq+1;j<=n;j++) {
						ROTATE(a,ip,j,iq,j)
					}
					for (j=1;j<=n;j++) {
						ROTATE(v,j,ip,j,iq)
					}
					++(*nrot);
				}
			}
		}
	}
	for (ip=1;ip<=n;ip++) {
		b[ip] += z[ip];
		d[ip]=b[ip];
		z[ip]=0.0;
	}
        printf("ERROR: Too many iterations in routine JACOBI\n");
}

#undef ROTATE


void findrotpar(struct pepfit fitu,struct pepfit sfit,int ures0,int ares0,int nres,
          float acm[3],float ucm[3],float rotm[3][3])
/* 
New version - use Kabsch routine (acta cryst A34, 827, 1978
x[i=1,2,3,4,..][i,j,k] is assigned fit vector (normalized to center)  // start with 1 so can use Recipe directly
  want to rotate this to position of y[] - unassigned fit
Rotate sfit onto ufit:   ufit = U X sfit
*/
{
   int i,j,k,nvec,count,ares,ures,im,nrot;
   float x[MAXV][4],y[MAXV][4],w[MAXV],r[4][4],rt[4][4],rr[4][4],rp[3][3],b[3][3];
   float d[4],v[4][4],a[3][3],mag;

   nvec=0;
   count = 0;
   for(k=0;k<=2;k++){
      acm[k]=0.0;
      ucm[k]=0.0;
   }
   for(i=0;i<=nres;i++){//center vectors about center of mass
        ares = ares0+i;
        ures = ures0+i;      
        for(im = 1;im<=4;im++){ // the 4 main chain atoms - 1=ca,2=c,3=o;4=n- Generate the matrix used for quaternions.
           count = count+1;
           for(k=0;k<=2;k++){
              acm[k] = acm[k]+ sfit.pos[ares][im][k];
              ucm[k] = ucm[k]+ fitu.pos[ures][im][k];
           }
        }
   }
   for(k=0;k<=2;k++){
      acm[k] = acm[k]/count;
      ucm[k] = ucm[k]/count;
   }
   count = 0;
   for(i=0;i<=nres;i++){//center vectors about center of mass
        ares = ares0+i;
        ures = ures0+i;      
        for(im = 1;im<=4;im++){ // the 4 main chain atoms - 1=ca,2=c,3=o;4=n- Generate the matrix used for quaternions.
           count = count+1;
           for(k=0;k<=2;k++){ // y = Ux
              x[count][k+1] = sfit.pos[ares][im][k]-acm[k];
              y[count][k+1] = fitu.pos[ures][im][k]-ucm[k];
           }
        }
   }
   nvec = count;  // number of vectors
   for(k=1;k<=nvec;k++)  // assign weights
       w[k] = 1.0;
   for(i=1;i<=3;i++)
     for(j=1;j<=3;j++)  
        r[i][j] = 0;
   for(i=1;i<=3;i++)
     for(j=1;j<=3;j++)  
       for(k=1;k<=nvec;k++)  // generate R matrix 
          r[i][j] = r[i][j]+w[k]*y[k][i]*x[k][j];
   for(i=1;i<=3;i++)// make rT matrix
     for(j=1;j<=3;j++)  
         rt[i][j] = r[j][i]; 
   for(i=1;i<=3;i++)
     for(j=1;j<=3;j++)
        rr[i][j]=0;
   for(i=1;i<=3;i++)  // rr = rt X rr
     for(j=1;j<=3;j++)
        for(k=1;k<=3;k++)
          rr[i][j]=rr[i][j]+rt[i][k]*r[k][j];
   jacobi(rr,3,d,v,&nrot);
   eigsrt(d,v,3);
   for(i=1;i<=3;i++) // Find eigenvectors (a[0,1,2][0,1,2]  which are columns of v
     for(k=1;k<=3;k++)
        a[i-1][k-1]= v[k][i];
#if 0 // for debugging
   float c[3];
   printf(" Eigenvectors:\n");
   for(i=0;i<=2;i++) // eigenvectors are columns of v
      printf("    %5.3f %5.3f %5.3f\n",a[i][0],a[i][1],a[i][2]);
   crossprod(a[0],a[1],c);
   printf(" a0 X a1 = %5.3f %5.3f %5.3f\n",c[0],c[1],c[2]);
   crossprod(a[0],a[2],c);
   printf(" a0 X a2 = %5.3f %5.3f %5.3f\n",c[0],c[1],c[2]);
#endif
   crossprod(a[0],a[1],a[2]);  // set a[2] = a0 X a1  for right handed system
   for(i=1;i<=3;i++) // rp[0,1,2][] = r[1,2,3][]
      for(k=1;k<=3;k++)
        rp[i-1][k-1] = r[i][k];
   for(i=0;i<=1;i++){
     multvec(rp,a[i],b[i]); /* returns vector v = m*v1*/
     mag = sqrt( dot(b[i],b[i]) );
     for(k=0;k<=2;k++)
       b[i][k]=b[i][k]/mag;
   }
   crossprod(b[0],b[1],b[2]);
   for(i=0;i<=2;i++) // Finally  form rotation matrix
     for(j=0;j<=2;j++)
        rotm[i][j] = 0.0;
   for(i=0;i<=2;i++) // Finally  form rotation matrix
     for(j=0;j<=2;j++)
       for(k=0;k<=2;k++)
         rotm[i][j] = rotm[i][j]+ b[k][i]*a[k][j];
}

void movefitaa(struct pepfit fit1,struct pepfit& fit2,float rotm[3][3],float tran[3],
    int maxside)
/*
  Similar to movefit, but rotates entire amino acid residue side
*/
{
   int i,ii,jj,rottype,ires,im,count,atomname; 
   
    allocatefit(fit2,fit1.nres +2); // make sure big enough allocation
    fit2.nres =fit1.nres;
    fit2.val1= fit1.val1;
    fit2.val2= fit1.val2;
    fit2.irank[0]=0;
    for(i=0;i<fit2.nres;i++){
          fit2.itype[i] = fit1.itype[i];
    }
    for(ires =0;ires<fit1.nres;ires++){ // main chain atoms
      for(im=1;im<=4;im++){ // 4 main chain atoms = 1=ca,2=c,3=o;4=n-
       if( (ires==0)&&(im==4) )
         goto nextim;
       if( (ires==(fit1.nres-1) ) && ( (im==2)||(im==3) ) )
         goto nextim;
       rottranvec(rotm,tran,fit1.pos[ires][im],fit2.pos[ires][im]);
       nextim:;
     } // im
    }
    fit2.nside[0]=fit2.nside[fit1.nres-1] = 0; 
    for(ires =1;ires<fit1.nres-1;ires++){ // side chain atoms
       if(fit1.nside[ires]>= 0){// side series of C atoms
          if(fit1.nside[ires]>= maxside)
              fit2.nside[ires] = maxside;
          else
              fit2.nside[ires] = fit1.nside[ires];
          for(i=1;i<=fit2.nside[ires];i++){
             rottranvec(rotm,tran,fit1.pos[ires][4+i],fit2.pos[ires][4+i]);
          }
       }
       else { // side is a amino acid
          fit2.irot[ires]=fit1.irot[ires];
          fit2.nside[ires] = fit1.nside[ires];
          fit2.cis[ires] = fit1.cis[ires];
          rottype = -fit1.nside[ires];
          count = 0; // counts side chain atoms copied
          for(ii=1;ii<=side[rottype][0][0]+1;ii++){ // copy each side chain branch
             for(jj=0;jj<side[rottype][ii][0];jj++){ 
               atomname=side[rottype][ii][jj+1];
               if(atomname>5){ // not part of mainchain or Cb
                     count = count+1;
                     rottranvec(rotm,tran,fit1.pos[ires][5+count],fit2.pos[ires][5+count]);
               }
               else // = 5 = Cb
                  rottranvec(rotm,tran,fit1.pos[ires][atomname],fit2.pos[ires][atomname]);
            }//jj
         }// ii
       } // else is amino acid
       nextires:;
    } // ires 
}


void movefit(struct pepfit fit1,struct pepfit fit2,float rotm[3][3],float tran[3],
    int maxside)
/*
  Similar to rottranfit - except uses rot matrix and tran vector
    F' = RxF + T
*/
{
   int i,jj,rottype,ires,im,count; 
   
    for(ires =0;ires<fit1.nres;ires++){ // main chain atoms
      for(im=1;im<=4;im++){ // 4 main chain atoms = 1=ca,2=c,3=o;4=n-
       if( (ires==0)&&(im==4) )
         goto nextim;
       if( (ires==(fit1.nres-1) ) && ( (im==2)||(im==3) ) )
         goto nextim;
       rottranvec(rotm,tran,fit1.pos[ires][im],fit2.pos[ires][im]);
       nextim:;
     } // im
    }
    fit2.nside[0]=fit2.nside[fit1.nres-1] = 0; 
    for(ires =1;ires<fit1.nres-1;ires++){ // side chain atoms
       if(fit1.nside[ires]>= 0){// side series of C atoms
          if(fit1.nside[ires]>= maxside)
              fit2.nside[ires] = maxside;
          else
              fit2.nside[ires] = fit1.nside[ires];
          for(i=1;i<=fit2.nside[ires];i++){
             rottranvec(rotm,tran,fit1.pos[ires][4+i],fit2.pos[ires][4+i]);
          }
       }
       else { // side is a amino acid
          rottype = -fit1.nside[ires];
          if(rottype==GLY)
              fit2.nside[ires] = 0;
          else{ // copy maxside atoms from first branch
             count = 0; // counts side chain atoms copied
             for(jj=1;jj<side[rottype][1][0];jj++){ //start with jj=1 (CB) copy atoms on first branch
               if(count<maxside){
                  if(rottype==PRO){ // special case, do not include 4th atom on branch = N
                      if(count==3){
                         fit2.nside[ires] = 3;
                         goto nextires;
                      }
                  }
                  else{
                     count = count+1;
                     rottranvec(rotm,tran,fit1.pos[ires][4+count],fit2.pos[ires][4+count]);
                  }
               } // count<
            }//jj
            fit2.nside[ires] = count;
         }// copy atoms
       } // else is amino acid
       nextires:;
    } // ires 
}


void rottranfit(struct pepfit fit1,struct pepfit fit2,float initca[3],float finalca[3],float rotm[3][3],
    int maxside)
/*
  Rotates/translates fit1 - returns new fit2
  IF FIT 1 HAS AMINO ACID SIDE CHAINS, CONVERTS SO THAT FIT2 JUST HAS SERIES OF C ATOMS
  Rotates and translates fit1 and returns the new fit2:  
  USES THE rotation paramters rotpar[0,1,2,3] =
     0,1,2 = direction  cosines of rotation axis (l.m.n) and rotpar[3] = theta = angle of rotation
  initca and finalca are the postitions of the Ca1 atom before and after translation
     and the rotation assume that Ca1 is at origin.  
   maxside = maximum number of side atoms to rotate (usually = 1).
Procedure:
  First translate ifit so Ca1 is at origin.
  Then use rotpar and rotate all the atoms
  then translate, to put Ca1 at new position.
*/
{
   int i,jj,rottype,ires,im,count; 
   float vnew[3];
   
    for(ires =0;ires<fit1.nres;ires++){ // main chain atoms
      for(im=1;im<=4;im++){ // 4 main chain atoms = 1=ca,2=c,3=o;4=n-
       if( (ires==0)&&(im==4) )
         goto nextim;
       if( (ires==(fit1.nres-1) ) && ( (im==2)||(im==3) ) )
         goto nextim;
       trannewm(fit1.pos[ires][im],fit2.pos[ires][im],initca); // translate to origin
       multvec(rotm,fit2.pos[ires][im],vnew); /*mult rot matrix times vector rpt */
       trannew(vnew,fit2.pos[ires][im],finalca); // translate to final Ca pos
       nextim:;
     } // im
    }
    fit2.nside[0]=fit2.nside[fit1.nres-1] = 0; 
    for(ires =1;ires<fit1.nres-1;ires++){ // side chain atoms
       if(fit1.nside[ires]>= 0){// side series of C atoms
          if(fit1.nside[ires]>= maxside)
              fit2.nside[ires] = maxside;
          else
              fit2.nside[ires] = fit1.nside[ires];
          for(i=1;i<=fit2.nside[ires];i++){
             trannewm(fit1.pos[ires][4+i],fit2.pos[ires][4+i],initca); // translate to origin
             multvec(rotm,fit2.pos[ires][4+i],vnew); /*mult rot matrix times vector rpt */
             trannew(vnew,fit2.pos[ires][4+i],finalca); // translate to final Ca pos
          }
       }
       else { // side is a amino acid
          rottype = -fit1.nside[ires];
          if(rottype==GLY)
              fit2.nside[ires] = 0;
          else{ // copy maxside atoms from first branch
             count = 0; // counts side chain atoms copied
             for(jj=1;jj<side[rottype][1][0];jj++){ //start with jj=1 (CB) copy atoms on first branch
               if(count<maxside){
                  if(rottype==PRO){ // special case, do not include 4th atom on branch = N
                      if(count==3){
                         fit2.nside[ires] = 3;
                         goto nextires;
                      }
                  }
                  else{
                     count = count+1;
                     trannewm(fit1.pos[ires][4+count],fit2.pos[ires][4+count],initca); // translate to origin
                     multvec(rotm,fit2.pos[ires][4+count],vnew); /*mult rot matrix times vector rpt */
                     trannew(vnew,fit2.pos[ires][4+count],finalca); // translate to final Ca pos
                  }
               } // count<
            }//jj
            fit2.nside[ires] = count;
         }// copy atoms
       } // else is amino acid
       nextires:;
    } // ires 
}


float rottranfunc(float p[],phead)
/* Used for powell minimization to find rotation/translation that
   optimizes fit to density.
Uses a funny set or parameters that is obtained from paper
  by MacKay (Acta Cryst A40, 165) that uses quaternions.
It seems to work OK - I think it sort of linearizes the problem.
   p[3,4,5]  correspod to  lt,mt,nt of MacKay.
     where l,m,n are directon cosines of axis of rotation vector
     and t = tan(theta/2) where theta = angle or rotation.
   p[1,2,3] = position of center of rotation.
The initial step size for p is 1.0, this is too large for rotation\
  Thus, scale pp down by factor of 10
*/
{
   int i,j;
   float nca[3],mag,theta,rotm[3][3],finalca[3],aveden,dist,frout;

   for(i=1;i<=6;i++)
       p[i] = p[i]/10.0;
   mag = sqr(p[3]) + sqr(p[4]) + sqr(p[5]);
   mag = sqrt(mag);
   dist = sqrt(sqr(p[1]) + sqr(p[2]) + sqr(p[3]) );// Total translation dist
#if 1
   if(dist > 1.5) // if dist too big return bad value
      return 100.0;
#endif
   if(mag<0.001){// then use unit rotm
      for(i=0;i<=2;i++)
        for(j=0;j<=2;j++)
           rotm[i][j] = 0.0;
      for(i=0;i<=2;i++)
         rotm[i][i] = 1.0;
   }
   else{
     theta = 2.0*atan(mag); // = angle of rotation about l,m,n axis  0<=theta<= PI, NOTE PROBLEM IF theta - PI
     for(i=0;i<3;i++) // Finally, solve for l, m, and n
        nca[i] = p[i+3]/mag;
     rotaxis(geo.de,geo.ep,theta,nca,rotm); //theta = angle of rotation
   }
   for(i=0;i<=2;i++)
        finalca[i] = origpsi[i]+p[i+1];
   rottranfit(fit1,fit2,origpsi,finalca,rotm,maxside);
   avedenfit0(den1,geo,fit2,&aveden,&frout);
#if 0
for(i=1;i<=6;i++)
   printf(" p[%d]= %5.3f",i,p[i]);
printf(" aveden = %5.3f\n",aveden);
#endif
    aveden = -aveden;  //change sign becuse want maximum value, but Powell finds minimum    
   return aveden;
}

void makenewrotfit(struct geometry& geo,struct griddata& den1,struct pepfit fit1,
       float brot[3][3],float btran[3],int maxside,int *newfit,float *aveden,float *frout)
/*  Rotates and translates fit[afit] to create fit[newfit]
*/
{
        *newfit = fitnum[0]+1; // Need to make the newfit in order to find fraction out of range
        fitnum[0] = fitnum[0]+1;
        fitnum[ fitnum[0] ] = fitnum[0];
        movefitaa(fit1,fit[*newfit],brot,btran,maxside);
        avedenfit0(den1,geo,fit[*newfit],aveden,frout); // frout = fraction out of range
}

void findbestsym3(int maxrot,struct geometry& geo,struct griddata& den1,
     int afit,int ufit,int ures0,int ares0,
     float *initden,float *finalden,float rotm[3][3],float tran[3],
     float minden)
/* 
maxrot is max number of residues to rotate in assigned fit
  if <0, rotate all.
Identical to findbestsym but inputs struct pepfit sfit in place of afit (fit[afit])
*/
{
   int i,nvar,ibeg,iend,isres,nres;
   float nca[3],frout;
   float ucm[3],acm[3],rotm2[3][3],rotm3[3][3];
   float pp[NDIM];
   float mag,theta;
   struct pepfit fit1,fit2,sfit;//NOTE fit2 is used in powellmin
// Following are not used but are passed in pfunc
     int cafix =0;
     float origphi[3],origpsi[3],ipos[NAFIT][NAA][3],fpos[NAFIT][NAA][3];
     float zaxis[3],xaxis[3]; // this is needed for pfunc definition but is not used.

   if(verbose>1)
      printf("Running findbestsym for ufit=%d ures0=%d ares0=%d\n",ufit,ures0,ares0);
   nvar = 6; // Number of adjustable variables
   fit1.nalloc = fit2.nalloc = sfit.nalloc = -10;
   if(maxrot<0){ // use the entire afit
       allocatefit(fit1,fit[afit].nres +2); // make sure big enough allocation
       allocatefit(fit2,fit[afit].nres +2); // make sure big enough allocation
       allocatefit(sfit,fit[afit].nres +2); // make sure big enough allocation
   }
   else{// copy afit to smaller sfit
      allocatefit(fit1,maxrot +3); // make sure big enough allocation
      allocatefit(fit2,maxrot +3); // make sure big enough allocation
      allocatefit(sfit,maxrot +3); // make sure big enough allocation
   }
    if((maxrot>0)&& (fit[afit].nres-2>maxrot) ){ // only check limited regin of afit
         if(ares0+maxrot>fit[afit].nres-1){
            iend = fit[afit].nres-1;
            ibeg = iend -maxrot; 
         }
         else{                      
          ibeg = ares0-1;
          iend = ibeg+maxrot;
        }
        isres = ares0-ibeg;
        copyfittofit(fit[afit],sfit,ibeg,iend);
     }// use limited region
     else{
       isres = ares0;
       copyfittofit(fit[afit],sfit,0,fit[afit].nres);
     }
     fit1.nres = fit2.nres= sfit.nres;
    if((sfit.nres-isres)<(fit[ufit].nres-ures0) ) 
       nres = sfit.nres-isres-2;
    else
       nres = fit[ufit].nres-ures0 -2;
#if 0 // Compares as many atoms as possible
    if(nres>MAXRES-1)
      nres = MAXRES-1;
#endif
#if 1 // Compares maximum of 3 - will not be influences as much by bad fit
    if(nres>3)
      nres = 3;
    else if (nres<1){ // skip the rest
       *finalden = minden -10; // set finalden < mindein - signal that bad nrs
       goto quitrout;
    }
#endif
    findrotpar(fit[ufit],sfit,ures0,isres,nres,acm,ucm,rotm2);// this returns the rotation matrix rotm
    rottranfit(sfit,fit1,acm,ucm,rotm2,maxside);
#if 0 // for debugging - find rmsd difference between atom positions
  int count;
  float avedist;
      count = 0;
      avedist = 0.0;      
      for(i=0;i<=nres;i++){// all the residues that are used to find rotpar
        for(k = 1;k<=4;k++){ // the 4 main chain atoms - 1=ca,2=c,3=o;4=n- Generate the matrix used for quaternions.
            count = count+1;
            avedist = avedist+ distv(fit1.pos[isres+i][k],fit[ufit].pos[ures0+i][k]);  
        }
      }//i
      avedist = avedist/count;
      printf(" AVE DIST: fit1-ufit = %5.3f\n",avedist);		
#endif     
#if 0 // Save the initial fit rotated
 int newfit; 
      newfit = fitnum[0]+1;
      fitnum[0] = fitnum[0]+1;
      fitnum[ fitnum[0] ] = fitnum[0];
      allocatefit(fit[newfit],sfit.nres +2); // make sure big enough allocation
      fit[newfit].nres = sfit.nres;
      fit[newfit].val1= sfit.val1;
      fit[newfit].val2= sfit.val2;
      fit[newfit].irank[0]=0;
      for(i=0;i<fit[newfit].nres;i++)
         fit[newfit].itype[i] = sfit.itype[i];
       rottranfit(sfit,fit[newfit],acm,ucm,rotm2,maxside);
      savefit();
      goto quitrout;
//exit(0);
#endif
      avedenfit0(den1,geo,fit1,initden,&frout); // frout = fraction out of range
      if(verbose>=2)
         printf("Initial den = %5.3f\n",*initden);
      for(i=0;i<=2;i++) // copy finalca to origpsi - passed to powellmin)
         origpsi[i]= ucm[i];
      powellmin(finalden,pfunc,rottranfunc); //(returns best pos = fpos; finalden is pointer Rigid refinemnent returns best tran, rot of ipos
      if(verbose>=2)
//      if(verbose>=1)
         printf("After powell, finalden = %5.3f\n",*finalden);
      if(*finalden>minden){  // return the overall best rotation and translation matrixs
        for(i=1;i<=6;i++) // pp is scaled by 10
          pp[i] = pp[i]/10.0;
        mag = sqr(pp[3]) + sqr(pp[4]) + sqr(pp[5]);
        mag = sqrt(mag);
        theta = 2.0*atan(mag); // = angle of rotation about l,m,n axis  0<=theta<= PI, NOTE PROBLEM IF theta - PI
        for(i=0;i<3;i++) // Finally, solve for l, m, and n
           nca[i] = pp[i+3]/mag;
        rotaxis(geo.de,geo.ep,theta,nca,rotm3); //theta = angle of rotation
        mult33_33(rotm3,rotm2,rotm); // return rotm
        multvec(rotm,acm,tran); /* returns vector v = m*v1*/
        for(i=0;i<=2;i++)
           acm[i] = origpsi[i]+pp[i+1];
        for(i=0;i<=2;i++) // return overall translation vector
           tran[i] = acm[i]-tran[i];
     }// if finalden>
#if 0 // Save the rotated
 int newfit; 
      newfit = fitnum[0]+1;
      fitnum[0] = fitnum[0]+1;
      fitnum[ fitnum[0] ] = fitnum[0];
      allocatefit(fit[newfit],sfit.nres +2); // make sure big enough allocation
      fit[newfit].nres = sfit.nres;
      fit[newfit].val1= sfit.val1;
      fit[newfit].val2= sfit.val2;
      fit[newfit].irank[0]=0;
      for(i=0;i<fit[newfit].nres;i++)
         fit[newfit].itype[i] = sfit.itype[i];
      rottranfit(fit1,fit[newfit],ucm,acm,rotm3,maxside);
      savefit();
//exit(0);
#endif
quitrout:;
   freepepfit(fit1);
   freepepfit(fit2);
   freepepfit(sfit);
}

void printsymoperator(float tran[],float rot[3][3])
// Prints to standard output the translation and rotation matrix
{
   printf("TRAN/ROT OPERATOR:  (Q' = ROT x Q  + TRAN):\n");
   printf("   Translation Vector =  %5.3f %5.3f %5.3f\n",tran[0],tran[1],tran[2]);
   printf("   Rotation Matrix  =    %5.3f %5.3f %5.3f\n",rot[0][0],rot[0][1],rot[0][2]);    
   printf("                         %5.3f %5.3f %5.3f\n",rot[1][0],rot[1][1],rot[1][2]);    
   printf("                         %5.3f %5.3f %5.3f\n",rot[2][0],rot[2][1],rot[2][2]);    
}

void  addrottran(int *nsym,float btran[NSYM][3],float brot[NSYM][3][3],
                float tran[3],float rot[3][3])
{
   int i,j;

   *nsym = *nsym+1;
   for(i=0;i<=2;i++){
      btran[*nsym][i]=tran[i];
      for(j=0;j<=2;j++)
         brot[*nsym][i][j]=rot[i][j];
   }
} 

void analyzesym(int *ifind,struct geometry& geo,struct griddata& den1,
    int *nsym,float brot[NSYM][3][3],float btran[NSYM][3],int del[],int afit,
    int ufit,int ures,int ares,float minden)
/*
Analyses the symmetery fit that is found:
  Checks if finalden is high enoguth and if atoms in density
     Then checks to see if rot/tran function is new, and saves
     it if it is
*nsym = # good and different rot/tran functions found
*/
{
   int i,j,k,ilim;
   float initden,finalden;
   float rot[3][3],tran[3];
   int newfit;
   float frout,trandiff,rotdiff;

   if(*nsym==NSYM-1){
      printf("ERROR: Reached allowed limited of number of rot/tran operators\n");
      printf("    Increase NSYM and recompile\n");
      goto quitrout;
   }
   trandiff = 3.0;  // difference in Angstroms required for new operator
   rotdiff =  0.1; // difference in rot matrix for new operator  
   ilim = -1; // if ilim = maxrot <0 rotate all the residues in assigned fit
   findbestsym3(ilim,geo,den1,afit,ufit,ures,ares,&initden,&finalden,
                    rot,tran,minden);
   if(finalden>minden){ 
        makenewrotfit(geo,den1,fit[afit],rot,tran,maxside,&newfit,&initden,&frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
        fitnum[0]= fitnum[0]-1; //this removes the new sym. mate fit
        if(verbose>=1){
            printf("FOUND SYM MATE - finalden for entire fit = %5.3f fraction of range = %5.3f\n",
                 finalden,frout);
            printf("   fit %d (res %d)is good sym. mate to fit %d  (res %d)\n",
                             ufit,ures,afit,ares);
            if(verbose>=2)
                printsymoperator(tran,rot);
         }
         if(frout> 0.5){ // do not use newfit and delete fit # = isym[1][0]
#if 0 // NEW 14nov00 - do not want to delete this fit
               del[0]=1; //delete 1 fit
               del[1] = ufit;
#endif
               if(verbose>=1){
                 printf("    Rot/tran operator is not used because > 50 percent atoms out of map range\n");
//                 printf("         and fit %d is markded for deletion\n",ufit);
               }
          }
          else{ // save this rot/tran - if new
             if(*nsym>0){ // check if unique
               for(i=1;i<= *nsym;i++){ // check if new operator
                 for(k=0;k<=2;k++){
                    if(fabs(tran[k]-btran[i][k])> trandiff)// then different 
                       goto nexti;
                    for(j=0;j<=2;j++){
                        if(fabs(rot[k][j]-brot[i][k][j])> rotdiff)// then different 
                            goto nexti;
                    }//j
                  }// k
                  if(verbose>=1)// if reach here new operator identical to i
                       printf("This rot/tran operator is NOT saved, because not new\n");
                  goto skipadd; 
                  nexti:;
               }// i
             } // nsym>0
             addrottran(nsym,btran,brot,tran,rot); 
             if(verbose>=1)
                  printf("Saving this operator as number = %d\n",*nsym);
             skipadd:;
          } // else check if new
     }// finalden>minden
     else{
          if(verbose>=1)
             printf("  Poor density for symmate\n");
     }
     quitrout:;
}


void analyzeonesym(int *ifind,struct geometry& geo,struct griddata& den1,
    float rot[3][3],float tran[3],int afit,
    int ufit,int ures,int ares,float *minden)
/*
NOTE: this is only called with minden = -10000 so that findsym3 always returns at rotm
Similar to analysesym - but does not check frout -
  Returns minden, and the rot/tran matrix

*/
{
   int ilim;
   float initden,finalden;
   int newfit;
   float frout;

   ilim = -1; // if ilim = maxrot <0 rotate all the residues in assigned fit
   findbestsym3(ilim,geo,den1,afit,ufit,ures,ares,&initden,&finalden,
                    rot,tran,*minden);//NOTE: minden = -10000.0, so always returns rotm
   if(finalden> *minden){//Only not true if nres <1
     makenewrotfit(geo,den1,fit[afit],rot,tran,maxside,&newfit,&initden,&frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
     fitnum[0]= fitnum[0]-1; //this removes the new sym. mate fit
     *minden = finalden;
     if(verbose>=1){
            printf("FOUND SYM MATE - finalden for entire fit = %5.3f fraction of range = %5.3f\n",
                 finalden,frout);
            printf("   fit %d (res %d)is good sym. mate to fit %d  (res %d)\n",
                             ufit,ures,afit,ares);
//            printsymoperator(tran,rot);
     }
   }
   quitrout:;
}

  
void symfit(int iplus,int *ifind,int icheck,struct geometry& geo,struct griddata& den1,
     int afit,int isym[NSYM][2],int *nsym,float brot[NSYM][3][3],float btran[NSYM][3],
     int del[])
/*
if iplus = 1 - then only check ufits>afit
if icheck = 1, only check fits that are from different files,
   i.e. have different absolute values of val2
   ALSO - icheck  IS ONLY SET = 1 FOR CASE WHERE CALLING AFTER MERGING TWO FILES.
    IN THIS CASE MAY HAVE VERY LARGE ASSIGNED FITS, ONLY CHECK SYM. FOR
    LIMITED REGION OF afit.
   NOTE: Criteria for sym fit is that second best den < 0.7 * best.
    (typically, second best = 20, best = 36 - but depends on lenght of
        assigned sequence that is rotated.  Longer the sequence, the
        more discriminating).
      THIS CRITERIA ASSUMEs THAT AT MOST 1 UNASSIGNED FIT IS SYM. RELATED TO
         afit.  If more than 1, then both will have high density, and
         symmetry will not be found. 
Tests assigned fit fit[afit] to see if there other unassigned fits with same
   helix/sheet type that it can be rigidly rotated and translated
   to get a high density.
For Ca1 to Ca(nmax = fit[afit].nres-2-nres) - compare
   with unassigned fits using Ca1 to Ca(nmax= fit[?].nres-2-nres
  where (nres+1) (=3) is the number of residues that are initial compared
  before finding optimal density.
NOTE: Do not have to use each assigned and unassigned as starting points,
  since, for example,  comparing assigned residues 1,2,3, with unassigned
  1,2,3 should give the same result as comparint assinged 2,3,4 with
  unassigned 2,3,4.
  THUS, test first assigned residue in helix/sheet string with all allowed
     unassigned residues (1,2,....N, since if unassigned, all must be same 
     constraint type) , and then need only test each subsequent assigned residue
     with unassinged res #1.
For supperimosing Cai to Cai+2 is assigned and unassgined
   and find best density of assigned - when trans/rot to unatsigned position
As test, Val 212 in Metb (fit4, res6)
     corresponds to helix fit2, res1 in Meta.
ares = assigned res number (runs from ares0 to ares0+nres)
ures = unasighed res number
afit = assigned fit number
ufit = unassigned fit number
nres = number of residues tested (0,1,..nres)
The assigned fit corresonds to r coord, and unassinged fit the r' cord which
    r will be rotated to by the l,m,n and theta
FOR NOW Assume Ca's that are compared do not include ends -so N and C O are always defined.
Checks all possible unassigned fits with same constraint and the
   ranks the qualality in rankf[1,2,..]
if(rankf[2]< 0.7 rankf[1]) then assumed region is a symtry mate.
RETURNS isym[0,1]
  isym[0] = -1 if no symmetry mate found
  else isym[1] = fit number of sym. mate.
    and isym[2] = sequence number that should be assigned to res1 in this sym. mate  
*/
{
   int i,j,nres0,nmax,atype,ires,tres,ilim;
   int ares0,ures0,ufit;
   int intprop[4],ranki[4][MAXRANK];
   float initden,finalden,rankf[MAXRANK],minden,mindentemp;
   float rotm[3][3],tran[3]; // rotation and best rotation
   int newfit,uresbeg,uresend,aresbeg,ierror,itotal;
   float frout,dist;

   if(verbose>=1)
     printf(" CHECKING fit %d for unassigned symetry mates (fitnum=%d)\n",afit,fitnum[0]);
   isym[0][0] = 0;
   minden = 20.0;  // Minimum average density required for potential sym. fit
//minden = 5.0;  // Minimum average density required for potential sym. fit
   nres0=2; // mininum number of residues to compare to find rotparameters
   ilim = 12; // if icheck = 1, copy region of afit of ilim residues to sfit
   ilim = 25; // if icheck = 1, copy region of afit of ilim residues to sfit
 //ilim = 30; // if icheck = 1, copy region of afit of ilim residues to sfit

#if 0 // debug - save rotate fit and look at it
   ures0=1;
   ires=1; // The correct res to compare with ufit = 1
   ufit=40; //use with shad_testsym
   findbestsym3(ilim,geo,den1,afit,ufit,ures0,ires,&initden,&finalden,
           rotm,tran,minden);
   atype= fit[afit].itype[ires];
   printf("Finalden  = %5.3f  Initden = %5.3f\n",finalden,initden);
   printf("   For type %d symfit for assign fit %d res %d, unassign fit %d res %d\n",
                          atype,afit,ires,ufit,ures0);
#if 0 // Save the rotated
      makenewrotfit(geo,den1,fit[afit],rotm,tran,maxside,&newfit,&initden,&frout);
      savefit();
#endif
   goto quitrout;
#endif

   nmax = fit[afit].nres -2-nres0;
   ares0 = 1;
   do{
     if(  (fit[afit].itype[ares0]==1)||(fit[afit].itype[ares0]==2) ){ // helix or sheet
         atype= fit[afit].itype[ares0];
         for(ires=ares0+1;ires<fit[afit].nres-1;ires++){ // check next residues, see if smae type
            if(fit[afit].itype[ires]!=atype){
               if( (ires-ares0)> nres0){ // then acceptabale string
                  tres = ires-1; // ares0 to tres is string of type atype
                  goto foundsequence;
               }
               else{ // not acceptable sequence
                  ares0 = ires+1;
                  goto nextares;
               }
             }  // != atype
             else if(ires==fit[afit].nres-2){ // reached end and still same type
                  if( (ires-ares0)>= nres0){ // then acceptabale string
                     tres = ires; // ares0 to tres is string of type atype
                     goto foundsequence;
                  }      
             } // if at end
          } // ires
      } // helix or sheet
      else{ // not helix or sheet
         ares0 = ares0+1;
         goto nextares;
      }
      foundsequence:; // if reach here - have sequence from ares0 to tres of type atype
//printf("Testing sequence %d to %d of assigned fit %d\n",ares0,tres,afit);
      ranki[0][0]=0;
      rankf[1]=rankf[2]=1.0;  //Default -corresponds to no rankings found
      for(ufit=1;ufit<=fitnum[0];ufit++){ 
//      for(ufit=40;ufit<=40;ufit++){ 
	    if( (icheck==1)&&( abs(fit[afit].val2)==abs(fit[ufit].val2) ) ) //If from same file, do not check
            goto nextufit; 
        if( (fit[ufit].val1>=0)&& (fit[ufit].itype[3]==atype) ){ // unassigned and of type atype
           for(ires = ares0;ires<=tres-nres0;ires++){// test assigned sequence starting with ares0 with all allowed unassigned
             for(ures0=1;ures0<=fit[ufit].nres-2-nres0;ures0++){
                if( (ires > ares0)&&(ures0>1) ){ // ONLY TEST ares0 with all ures0, and other ires just with first ures0 (SEE ABOVE)
                   goto nextires;
                }
                findbestsym3(ilim,geo,den1,afit,ufit,ures0,ires,&initden,&finalden,
                               rotm,tran,minden);
                if(finalden<=minden)//NEW 11sep01 - rotm only returned if finalden>minden
                   goto nextires;
                makenewrotfit(geo,den1,fit[afit],rotm,tran,maxside,&newfit,&initden,
                      &frout);//NEW 14nov00 - need to find frout
                fitnum[0]=fitnum[0]-1; // remove the new fit made by makenewrotfit
                if(frout>0.1) // do not use
                    finalden = -1000.0;
#if 0
                if(finalden>minden){ 
                    if(finalden > rankf[1]){   // store the best rot, tran matrix
                        for(j=0;j<=2;j++){
                          btran[j]= tran[j];
                          for(k=0;k<=2;k++)
                             brot[j][k] = rotm[j][k];
                        }
                    } 
                    if(verbose>=2){
//                    if(verbose>=1){
                       printf("Finalden  = %5.3f  Initden = %5.3f\n",finalden,initden);
                       printf("   For type %d symfit for assign fit %d res %d, unassign fit %d res %d\n",
                          atype,afit,ires,ufit,ures0);
                    }
                } // if finalden
#endif
                intprop[0]=ires;
                intprop[1]=ufit;
                intprop[2]=ures0;
                intprop[3] = 0;
                newrankqual(finalden,ranki,intprop,rankf);
                nextires:;
             }// ures0
           }// ires
         } // if atype
         nextufit:;
      }// ufit
      if(verbose>=2)
//      if(verbose>=1)
//        for(i=1;i<=ranki[0][0];i++)
        for(i=1;i<=3;i++)
            printf("  qual = %5.3f  ares0=%d  ufit=%d  ures0=%d\n",
                       rankf[i],ranki[0][i],ranki[1][i],ranki[2][i]);
      mindentemp= 0.0;//set minden = 0 here so findbestsym3 should always return a rotm.
      if( (rankf[2]<= 0.7*rankf[1])&& (rankf[1]>minden) ){ // Indicate found potential sym mate     
         ufit = ranki[1][1];
         findbestsym3(-1,geo,den1,afit,ranki[1][1],ranki[2][1],
                      ranki[0][1],&initden,&finalden,rotm,tran,mindentemp);// set ilim = -1 and minden = 0.0
         if(finalden<=mindentemp){//New 11sep01 - case where no rotm returned NOTE- this should not happen, but just in case!
            goto badden;
         }
         makenewrotfit(geo,den1,fit[afit],rotm,tran,maxside,&newfit,&initden,&frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
         fitnum[0]=fitnum[0]-1; // remove the new fit made by makenewrotfit
         if(ranki[0][1]>=ranki[2][1])//Find region of overlap between afit and ufit -ranki[2][1]=ures, ranki[0][1]=ares
            uresbeg = 0;
         else
            uresbeg = ranki[2][1]-ranki[0][1]; 
         if( (fit[ufit].nres -ranki[2][1])<=(fit[afit].nres -ranki[0][1]) )
            uresend = fit[ufit].nres-1;
         else
            uresend = fit[ufit].nres -1 -
                ( (fit[ufit].nres -ranki[2][1]) - (fit[afit].nres -ranki[0][1]) );
         aresbeg = ranki[0][1] - (ranki[2][1] - uresbeg);
         ierror = itotal = 0;
         for(i=uresbeg;i<=uresend;i++){
            for(j=1;j<=4;j++){
              if( (i==uresbeg)&& (j==4) )
                 goto nextj;
              if( (i==uresend)&&( (j==2)||(j==3) ) )
                 goto nextj;
              dist = distv( fit[ufit].pos[i][j],fit[newfit].pos[aresbeg+i-uresbeg][j]);
              if(dist> 1.0)
                 ierror = ierror+1;
              itotal = itotal+1;
              nextj:;
            } // j
         } //i
         if( (initden > minden) &&  ((float)ierror/itotal < 0.15) ){       
            isym[0][0]=isym[0][0]+1;
            if(isym[0][0]>= NSYM){
               printf("ERROR: too many unassigned fits related to assigned fit %d\n",afit);
                goto foundsym;
            }
            isym[ isym[0][0] ][0] = ranki[1][1]; // the sym fit #
            isym[ isym[0][0] ][1] = (-fit[afit].val1 + ranki[0][1] - ranki[2][1]); // the sequence number of res1 in this fit
            if(verbose>=1)
                printf("Found NCS fit # %d: aveden = %5.3f,ierror =%d fit #=%d seq #= %d\n",
                             ufit,initden,ierror,isym[ isym[0][0] ][0],isym[ isym[0][0] ][1]);
#if 0 // for degugging - save rotated fit
  fitnum[0]=fitnum[0]+1;
  savefit();
#endif 
            if(*ifind==1)// then quit after finding sym mate - return best rot/tran for entire fit
               goto foundsym;
         }
         else{
           if(verbose>=1)
              printf("Found NCS symfit but NOT used because bad density = %5.3f, and/or %d bad atoms/%d total\n",initden,ierror,itotal);
         }
      }// if rank
      badden:;
      ares0 = tres+1; // ONLY reach here if found testable sequence
      nextares:;
//printf(" ares0 = %d\n",ares0);
   } while(ares0<nmax);
   foundsym:;
   if( (*ifind==1)&&(isym[0][0]>0) )
      analyzesym(ifind,geo,den1,nsym,brot,btran,del,afit,ranki[1][1],ranki[2][1],
                ranki[0][1],minden);
   quitrout:;
}

void symfitass(int iplus,int *ifind,int icheck,struct geometry& geo,struct griddata& den1,
     int afit,int isym[NSYM][2],int *nsym,float brot[NSYM][3][3],float btran[NSYM][3],int del[])
/*
if iplus = 1 - then only check ufits>afit
 Similar to symfit - but only for assigned fits-
   Symmetry related fits are defined as having same amino acid seq. 
      assignment
*/
{
   int ufit,ustart,astart,ures,ares;
   float minden;
    
   if(verbose>=1)
        printf(" CHECKING fit %d for assigned symetry mates (fitnum=%d)\n",afit,fitnum[0]);
    minden = 20.0;  // Minimum average density required for potential sym. fit
    astart = -fit[afit].val1;
    for(ufit=1;ufit<=fitnum[0];ufit++){ 
//    for(ufit=40;ufit<=40;ufit++){ 
	if( (icheck==1)&&( abs(fit[afit].val2)==abs(fit[ufit].val2) ) ) //If from same file, do not check
             goto nextufit; 
        if( (ufit!=afit)&&(fit[ufit].val1<0) ){ // assigned 
            ustart = -fit[ufit].val1;
            ures = -1; // idicates not sym fit found
            if( ustart <= astart){
                if( ustart+fit[ufit].nres-1 >= astart){
                   ures = astart-ustart+1; 
                   ares = 1;
                }
            }
            else {// ustart>astart
                if( astart+fit[afit].nres-1 >= ustart){
                   ures = 1; 
                   ares = ustart -astart+1;
                }
            }
            if(ures>=0){ // indicates found sym fit
               if(verbose>=1)
                  printf(" Assigned Fits %d and %d overlap\n",afit,ufit); 
               analyzesym(ifind,geo,den1,nsym,brot,btran,del,afit,
                   ufit,ures,ares,minden);
            }
         }
         nextufit:;
    }//ufit
}

void symonefitass(int *ifind,struct geometry& geo,struct griddata& den1,
     int afit,int isym[NSYM][2],int *nsym,float rot[3][3],float tran[3],
     int del[],int ufit,float *minden)
/*
 Similar to symfitass - but only check one fit
 Returns ifind =1 if find density >minden, - no conditons of frout
*/
{
   int ustart,astart,ures,ares;
    
   *minden = -10000.0;
   if(verbose>=1)
        printf(" CHECKING fit %d for assigned symetry mates (fitnum=%d)\n",afit,fitnum[0]);
    astart = -fit[afit].val1;
        if( (ufit!=afit)&&(fit[ufit].val1<0) ){ // assigned 
            ustart = -fit[ufit].val1;
            ures = -1; // idicates not sym fit found
            if( ustart <= astart){
                if( ustart+fit[ufit].nres-1 >= astart){
                   ures = astart-ustart+1; 
                   ares = 1;
                }
            }
            else {// ustart>astart
                if( astart+fit[afit].nres-1 >= ustart){
                   ures = 1; 
                   ares = ustart -astart+1;
                }
            }
            if(ures>=0){ // indicates found sym fit
               if(verbose>=1)
                  printf(" Assigned Fits %d and %d overlap\n",afit,ufit); 
               analyzeonesym(ifind,geo,den1,rot,tran,afit,ufit,ures,
                     ares,minden);
            }
         }
         nextufit:;
}

void refinerottran(struct geometry& geo,struct griddata& den1,
        float brot[3][3],float btran[3],int maxside,float rotm[3][3],
        float tran[3],int ifit,float minden,float *finalden,float *frout)
/* Uses brot and btran as starting point and return refined rot/tran vector
     and refined fit as fitnum[0] (increased by 1)
*/
{
   int ii,k,im,nvar,ares,newfit,count;
   float nca[3],initden;
   float acm[3];
   float pp[NDIM];
   float mag,theta;
   struct pepfit fit1,fit2;//NOTE fit2 is used in powellmin
// Following are not used but are passed in pfunc
     int cafix =0;
     float origphi[3],origpsi[3],ipos[NAFIT][NAA][3],fpos[NAFIT][NAA][3];
     float zaxis[3],xaxis[3]; // this is needed for pfunc definition but is not used.

   nvar = 6; // Number of adjustable variables
   fit1.nalloc = fit2.nalloc = -10;
          allocatefit(fit1,fit[ifit].nres +2); // make sure big enough allocation
          allocatefit(fit2,fit[ifit].nres +2); // make sure big enough allocation
          fit1.nres = fit2.nres= fit[ifit].nres;
          makenewrotfit(geo,den1,fit[ifit],brot,btran,maxside,
              &newfit,&initden,frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
          copyfittofit(fit[newfit],fit1,0,fit[ifit].nres);
          fitnum[0]=fitnum[0]-1; // get rid of newfit
          count = 0; // find center of mass to refine rot/tran about
          for(k=0;k<=2;k++)
              acm[k]=0.0;
          for(ares=1;ares<=fit[ifit].nres-2;ares++){//center vectors about center of mass
             for(im = 1;im<=4;im++){ // the 4 main chain atoms - 1=ca,2=c,3=o;4=n- Generate the matrix used for quaternions.
                count = count+1;
                for(k=0;k<=2;k++)
                   acm[k] = acm[k]+ fit1.pos[ares][im][k];
             }
          }
          for(k=0;k<=2;k++)
            acm[k] = acm[k]/count;
          for(k=0;k<=2;k++) 
             origpsi[k]= acm[k];
          powellmin(finalden,pfunc,rottranfunc); //NOTE - finalden is for max. of 3 side chain atoms(returns best pos = fpos; finalden is pointer Rigid refinemnent returns best tran, rot of ipos
          if(*finalden>minden){  // return the overall best rotation and translation matrixs
            for(k=1;k<=6;k++) // pp is scaled by 10
              pp[k] = pp[k]/10.0;
            mag = sqr(pp[3]) + sqr(pp[4]) + sqr(pp[5]);
            mag = sqrt(mag);
            if(mag<0.001){// then use unit rotm
              for(ii=0;ii<=2;ii++)
               for(k=0;k<=2;k++)
                 rotm[ii][k] = 0.0;
              for(ii=0;ii<=2;ii++)
                rotm[ii][ii] = 1.0;
            }
            else{
                theta = 2.0*atan(mag); // = angle of rotation about l,m,n axis  0<=theta<= PI, NOTE PROBLEM IF theta - PI
                for(ii=0;ii<3;ii++) // Finally, solve for l, m, and n
                    nca[ii] = pp[ii+3]/mag;
                rotaxis(geo.de,geo.ep,theta,nca,rotm); //theta = angle of rotation
            }
//            mult33_33(rotm3,rotm2,rotm); // return rotm
            multvec(rotm,acm,tran); /* returns vector v = m*v1*/
            for(ii=0;ii<=2;ii++)
               acm[ii] = origpsi[ii]+pp[ii+1];
            for(ii=0;ii<=2;ii++) // return overall translation vector
               tran[ii] = acm[ii]-tran[ii];
            *finalden = -100.0; // default value
            makenewrotfit(geo,den1,fit1,rotm,tran,maxside,
                  &newfit,finalden,frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
          } // finalden>minden
          else{ // set to very bad value
             *finalden = -10000.0;
          }
    freepepfit(fit1);
    freepepfit(fit2);
}


void applyrottran(struct geometry& geo,struct griddata& den1,int nsym,
             float brot[NSYM][3][3],float btran[NSYM][3],int del[MAXDEL])
/* applies rot/tran operatory to each fit - then uses rigid body Powell refinement,
     Checks finalden
    If finalden> minden and frout < 0.5  then add fit.
*/
{
   int ifit,i,id,origfitnum,iselect;
   float minden,frout,finalden;
   float rotm[3][3],tran[3];
// Following are not used but are passed in pfunc
 
   if(verbose>=1)
        printf("\n\n APPLYING ROT/TRAN TO EACH FIT AND REMOVING OVERLAPING FITS\n");
   iselect = 0; // indicates use all checks in removebad
   origfitnum=fitnum[0];
   minden = 20.0;  // Minimum average density required for potential sym. fit
   for(i=1;i<=nsym;i++){
       for(ifit=1;ifit<=origfitnum;ifit++){// ue origfitnum because fitnum[0] changes
          for(id= 1;id<=del[0];id++) // marked to be deleted
            if(ifit ==  del[id] )
              goto nextifit;
          refinerottran(geo,den1,brot[i],btran[i],maxside,rotm,tran,ifit,minden,
                 &finalden,&frout);
#if 0
          allocatefit(fit1,fit[ifit].nres +2); // make sure big enough allocation
          allocatefit(fit2,fit[ifit].nres +2); // make sure big enough allocation
          fit1.nres = fit2.nres= fit[ifit].nres;
          makenewrotfit(geo,den1,fit[ifit],brot[i],btran[i],maxside,
              &newfit,&initden,&frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
          copyfittofit(fit[newfit],fit1,0,fit[ifit].nres);
          fitnum[0]=fitnum[0]-1; // get rid of newfit
          count = 0; // find center of mass to refine rot/tran about
          for(k=0;k<=2;k++)
              acm[k]=0.0;
          for(ares=1;ares<=fit[ifit].nres-2;ares++){//center vectors about center of mass
             for(im = 1;im<=4;im++){ // the 4 main chain atoms - 1=ca,2=c,3=o;4=n- Generate the matrix used for quaternions.
                count = count+1;
                for(k=0;k<=2;k++)
                   acm[k] = acm[k]+ fit1.pos[ares][im][k];
             }
          }
          for(k=0;k<=2;k++)
            acm[k] = acm[k]/count;
          for(k=0;k<=2;k++) 
             origpsi[k]= acm[k];
          powellmin(&finalden,pfunc,rottranfunc); //NOTE - finalden is for max. of 3 side chain atoms(returns best pos = fpos; finalden is pointer Rigid refinemnent returns best tran, rot of ipos
          if(finalden>minden){  // return the overall best rotation and translation matrixs
            for(k=1;k<=6;k++) // pp is scaled by 10
              pp[k] = pp[k]/10.0;
            mag = sqr(pp[3]) + sqr(pp[4]) + sqr(pp[5]);
            mag = sqrt(mag);
            if(mag<0.001){// then use unit rotm
              for(ii=0;ii<=2;ii++)
               for(k=0;k<=2;k++)
                 rotm[ii][k] = 0.0;
              for(ii=0;ii<=2;ii++)
                rotm[ii][ii] = 1.0;
            }
            else{
                theta = 2.0*atan(mag); // = angle of rotation about l,m,n axis  0<=theta<= PI, NOTE PROBLEM IF theta - PI
                for(ii=0;ii<3;ii++) // Finally, solve for l, m, and n
                    nca[ii] = pp[ii+3]/mag;
                rotaxis(geo.de,geo.ep,theta,nca,rotm); //theta = angle of rotation
            }
//            mult33_33(rotm3,rotm2,rotm); // return rotm
            multvec(rotm,acm,tran); /* returns vector v = m*v1*/
            for(ii=0;ii<=2;ii++)
               acm[ii] = origpsi[ii]+pp[ii+1];
            for(ii=0;ii<=2;ii++) // return overall translation vector
               tran[ii] = acm[ii]-tran[ii];
          if(finalden>minden){  // return the overall best rotation and translation matrixs
            finalden = -100.0; // default value
            makenewrotfit(geo,den1,fit1,rotm,tran,maxside,
                  &newfit,&finalden,&frout);//Make newwfit = fitnum[0]+1 and increae fitnum[0] Need to make in order to find frout
#endif
            if( finalden <= -900.0) // indicates never called last makenewrotfit( 
               goto nextifit;
            if( (finalden < minden)||(frout >0.5) ){ // delete new fit
               fitnum[0] = fitnum[0]-1;
               if(verbose>=1)
                   printf("isym=%d ifit = %d: NOT USED because  finalden = %5.3f frout = %5.3f\n",i,ifit,finalden,frout);
            }
            else{  // check new fit for identical
              removeonebadfit(fitnum[0],iselect,den1,geo,del,ifit);//ifit = fit # that new fit was made from f by symop 
              if(del[ del[0] ] == fitnum[0]){ // for this case newfit completely overlapped or bad assignment- do not add
                 if(verbose>=1)
                   printf("isym=%d ifit = %d: NOT USED because  completely overlaped by another fit\n",i,ifit);
                 del[0] = del[0]-1;
                 fitnum[0]=fitnum[0]-1;
              }
              else{              
                 if(verbose>=1)
                   printf("isym=%d ifit = %d: used\n",i,ifit);
              }
            }// check for identical
          nextifit:;
       }//ifit
    }// isym
}       

