/*******************************************************************************
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 * 
 * Contact Info:
 * 	Bruce Donald
 * 	Duke University
 * 	Department of Computer Science
 * 	Levine Science Research Center (LSRC)
 * 	Durham
 * 	NC 27708-0129 
 * 	USA
 * 	brd@cs.duke.edu
 * 
 * Copyright (C) 2011 Jeffrey W. Martin and Bruce R. Donald
 * 
 * <signature of Bruce Donald>, April 2011
 * Bruce Donald, Professor of Computer Science
 ******************************************************************************/


package edu.duke.donaldLab.jdshot.compute;

import java.util.Iterator;

import edu.duke.donaldLab.jdshot.analyze.AnalysisContext;
import edu.duke.donaldLab.jdshot.grid.Axes;
import edu.duke.donaldLab.jdshot.grid.dn.DnGridPoint;
import edu.duke.donaldLab.jdshot.search.SearchContext;
import edu.duke.donaldLab.share.geom.Vector3;
import edu.duke.donaldLab.share.protein.Atom;
import edu.duke.donaldLab.share.protein.AtomAddressInternal;
import edu.duke.donaldLab.share.protein.Subunit;

public class AxesFinder
{
	/**************************
	 *   Data Members
	 **************************/
	
	SearchContext m_searchContext;
	AnalysisContext m_analysisContext;
	
	
	/**************************
	 *   Constructors
	 **************************/
	
	public AxesFinder( SearchContext searchContext, AnalysisContext analysisContext )
	{
		m_searchContext = searchContext;
		m_analysisContext = analysisContext;
	}
	
	
	/**************************
	 *   Methods
	 **************************/
	
	public DnGridPoint find( )
	{
		// get our subunits
		Subunit monomerSubunit = m_searchContext.getMonomer();
		Subunit mainSubunit = m_analysisContext.getOligomer().getSubunit( 0 );
		Subunit intradimerSubunit = m_analysisContext.getOligomer().getSubunit( 2 );
		Subunit interdimerSubunit = m_analysisContext.getOligomer().getSubunit( 1 );
		
		// first, find the offset when we center the main subunit
		Atom monomerAtom = monomerSubunit.getAtom( monomerSubunit.getBackboneAtomIndex().get( 0 ) );
		Atom mainSubunitAtom = monomerSubunit.getAtom( mainSubunit.getBackboneAtomIndex().get( 0 ) );
		Vector3 offset = new Vector3( monomerAtom.getPosition() );
		offset.subtract( mainSubunitAtom.getPosition() );
		
		// offset the entire oligomer
		for( int i=0; i<4; i++ )
		{
			Iterator<AtomAddressInternal> iterAtom = m_analysisContext.getOligomer().getSubunit( i ).getAtomIndex().iterator();
			while( iterAtom.hasNext() )
			{
				Atom atom = m_analysisContext.getOligomer().getAtom( iterAtom.next() );
				atom.getPosition().add( offset );
			}
		}
		
		// first, find the geometric center of the oligomer - this will be our t
		Vector3 center = new Vector3();
		int count = 0;
		for( int i=0; i<4; i++ )
		{
			Iterator<AtomAddressInternal> iterAtom = m_analysisContext.getOligomer().getSubunit( i ).getAtomIndex().iterator();
			while( iterAtom.hasNext() )
			{
				Vector3 pos = m_analysisContext.getOligomer().getAtom( iterAtom.next() ).getPosition();
				center.add( pos );
				count++;
			}
		}
		center.scale( 1 / (double)count );
		
		// then, find the geometric center of the dimer
		Vector3 dimerCenter = new Vector3();
		count = 0;
		Iterator<AtomAddressInternal> iterAtom = mainSubunit.getAtomIndex().iterator();
		while( iterAtom.hasNext() )
		{
			Vector3 pos = mainSubunit.getAtom( iterAtom.next() ).getPosition();
			dimerCenter.add( pos );
			count++;
		}
		iterAtom = intradimerSubunit.getAtomIndex().iterator();
		while( iterAtom.hasNext() )
		{
			Vector3 pos = intradimerSubunit.getAtom( iterAtom.next() ).getPosition();
			dimerCenter.add( pos );
			count++;
		}
		dimerCenter.scale( 1 / (double)count );
		
		// now, we have the flip axis
		Vector3 flip = new Vector3( dimerCenter );
		flip.subtract( center );
		flip.normalize();
		
		// to find the rotation vector, we need the center of the bizarro dimer
		Vector3 bizarroDimerCenter = new Vector3();
		count = 0;
		iterAtom = mainSubunit.getAtomIndex().iterator();
		while( iterAtom.hasNext() )
		{
			Vector3 pos = mainSubunit.getAtom( iterAtom.next() ).getPosition();
			bizarroDimerCenter.add( pos );
			count++;
		}
		iterAtom = interdimerSubunit.getAtomIndex().iterator();
		while( iterAtom.hasNext() )
		{
			Vector3 pos = interdimerSubunit.getAtom( iterAtom.next() ).getPosition();
			bizarroDimerCenter.add( pos );
			count++;
		}
		bizarroDimerCenter.scale( 1 / (double)count );
		
		// now, the rotation axis is the flip axis cross the bizarro flip axis
		Vector3 bizarroFlip = new Vector3( bizarroDimerCenter );
		bizarroFlip.subtract( center );
		bizarroFlip.normalize();
		Vector3 rotation = new Vector3();
		flip.getCross( rotation, bizarroFlip );
		
		// compute the grid point from the axes
		DnGridPoint point = new DnGridPoint();
		point.x = center.x;
		point.y = center.y;
		point.z = center.z;
		
		// back-calculate the angles
		double[] angles = Axes.getAngles( rotation, flip );
		point.phi = angles[0];
		point.theta = angles[1];
		point.rho = angles[2];
		
		// TEMP: make sure if we backcompute from this point, we get similar axes
		Vector3 reRotation = new Vector3();
		Vector3 reFlip = new Vector3();
		Axes.getRotationVector( reRotation, point.phi, point.theta );
		Axes.getFlipVector( reFlip, point.phi, point.theta, point.rho );
		System.out.println( "rotation:\t" + rotation );
		System.out.println( "re rotation:\t" + reRotation );
		System.out.println( "flip:\t\t" + flip );
		System.out.println( "re flip:\t" + reFlip );
		
		// TEMP: make sure the two axes are perpendicular
		System.out.println( "dot product: " + rotation.getDot( flip ) );
		System.out.println( "re dot product: " + reRotation.getDot( reFlip ) );
		
		return point;
	}
}
