# include "centroid.h"

int FileSelectToXYZ(file_records *f, NRVectP X,NRVectP Y,NRVectP Z)
{
  int ii;
  int npts = 0;
  
  for (ii = 0;ii<f->atomnum; ii++) {
    atom_record * atom = f->atoms[ii];
    if (SELECTED_ATOM_P(atom)) {
      npts++ ;
      X[npts] = atom->x ;
      Y[npts] = atom->y ;
      Z[npts] = atom->z ;
    }
  }
  if (npts == 0) 
    STDERR("npts == 0");
  return npts ;
}


void XYZToFileSelect(file_records *f, NRVectP X,NRVectP Y,NRVectP Z, int N)
{
  int ii;
  int npts = 0;
  
  for (ii = 0;ii<f->atomnum; ii++) {
    atom_record * atom = f->atoms[ii];
    if (SELECTED_ATOM_P(atom)) {
      npts++ ;
      atom->x = X[npts] ;
      atom->y = Y[npts] ;
      atom->z = Z[npts] ;
    }
  }
  if (npts > N) 
    STDERR("npts > N");
}


int CentroidFromStructure(file_records *f, NRVectP v)
{
  int 
    an = f->atomnum ;
  NRVectP 
    X = dvector(1,an),
    Y = dvector(1,an),
    Z = dvector(1,an);
  int n = FileSelectToXYZ(f,X,Y,Z);
  v[1] = NRVect_Mean(X,n);
  v[2] = NRVect_Mean(Y,n);
  v[3] = NRVect_Mean(Z,n);
  free_dvector(X,1,an);
  free_dvector(Y,1,an);
  free_dvector(Z,1,an);
  return n;
}

void HelixTipsFromStructure(file_records *f, NRVectP ntip, NRVectP ctip)
{
  int 
    an = f->atomnum ;
  NRVectP 
    X = dvector(1,an),
    Y = dvector(1,an),
    Z = dvector(1,an);
  AllocateNRVect3 m, EigenVals;
  NRMatrix33P C = dmatrix(1,3,1,3);
  NRMatrix33P EigenVecs = dmatrix(1,3,1,3);
  int N = FileSelectToXYZ(f,X,Y,Z);
  m[1] = -NRVect_Mean(X,N);
  m[2] = -NRVect_Mean(Y,N);
  m[3] = -NRVect_Mean(Z,N);
  NRTOE(m);
  TranslateXYZby1(m,X,Y,Z,N);
# ifdef OFF
  XYZToFileSelect(f,X,Y,Z,N);
  write_pdb_file(stdout,f,0);
# endif
  NR_CalcCovarMat33(C,X,Y,Z,N);
  DiagMatrix33(C,EigenVals,EigenVecs);
  NRTOE(EigenVals);
# ifdef OFF
  PermuteEigenValsVecsSoLargestLast(EigenVals,EigenVecs); 
  NRTOE(EigenVals);
  assert(EigenVals[3] >= EigenVals[2]);
  assert(EigenVals[3] >= EigenVals[1]);
# endif
  Transpose(EigenVecs);
  RotateXYZby1(EigenVecs,X,Y,Z,N);
# ifdef OFF
  XYZToFileSelect(f,X,Y,Z,N);
  write_pdb_file(stdout,f,0);
  exit ;
# endif
  {
    AllocateNRVect3 V ;
    Transpose(EigenVecs);
    m[1] *= -1 ; m[2] *= -1; m[3] *= -1;

    if (LargestElementIndex(EigenVals,3) == 2) { 
      V[1]=0; V[2]= Y[1]; V[3] = 0;
    } else if (LargestElementIndex(EigenVals,3) == 1) {
      V[1]=X[1]; V[2]= 0; V[3] = 0;
    } else  {
      V[1]=0; V[2]= 0; V[3] = Z[1];
    }
    VectRotate_1by2(V,EigenVecs);
    VectTranslate_1by2(V,m);
    CopyNRVect3_1to2(V,ntip);
    NRTOE(ntip);

    if (LargestElementIndex(EigenVals,3) == 2) {
      V[1]=0; V[2]= Y[N]; V[3] = 0;
    } else if (LargestElementIndex(EigenVals,3) == 1)  {
      V[1]=X[N]; V[2]= 0; V[3] = 0;
    } else  {
      V[1]=0; V[2]= 0; V[3] = Z[N];
    }
    VectRotate_1by2(V,EigenVecs);
    VectTranslate_1by2(V,m);
    CopyNRVect3_1to2(V,ctip);
    NRTOE(ctip);
  }

  free_dmatrix(C,1,3,1,3);
  free_dmatrix(EigenVecs,1,3,1,3);
  free_dvector(X,1,an);
  free_dvector(Y,1,an);
  free_dvector(Z,1,an);
}



double FindFitTranslationOnly
(double * TVect, file_records *Moving, file_records *Fixed)
{ 
  AllocateNRVect3 vF, vM;
  int nFix = CentroidFromStructure(Fixed,vF) ;
  int nMove = CentroidFromStructure(Moving,vM) ;
  fprintf(stderr,"FindFitTranslationOnly(): Fixed (centroid)=( %7.4f %7.4f %7.4f )\n",
	  vF[1],vF[2],vF[3]);
  fprintf(stderr,"FindFitTranslationOnly(): Moving(centroid)=( %7.4f %7.4f %7.4f )\n",
	  vM[1],vM[2],vM[3]);
  if (nMove != nFix ) {
    fprintf(stderr,
	    "FindFitTranslationOnly(): %6d = nFix != nMove = %6d\n",nMove,nFix);
  }
  VectSub_1is2minus3(TVect,vF,vM);
  return 0;
}