/* Copyright (C) 2001-2003 Damir Zucic */

/*=============================================================================

			    is_part_of_ring.c

Purpose:
	Check is the specified bond part of a ring or of some other closed
	structural element.  The array indices of two atoms which form the
	bond are used to specify the chosen bond.

Input:
	(1) Pointer to MolComplexS structure.
	(2) The macromolecular complex index.
	(3) The array index of the first atom of the chosen bond.
	(4) The array index of the second atom of the chosen bond.

Output:
	(1) Return value.

Return value:
	(1) Positive if bond is part of a ring.
	(2) Zero if  the specified bond is  not involved in formation of a
	    ring or some other closed structural element.
	(3) Negative,  if the macromolecular complex contains no atoms, if
	    two specified atoms  do not belong to  the same residue  or if
	    this residue does not contain any atoms (rather strange).

Notes:
	(1) Imagine that chemical bonds form  a system of tunnels dug in a
	    mountain. Atoms are treated as places where two, three or more
	    tunnels are connected. The shock wave is initiated starting at
	    the position of the first atom,  propagating in a direction of
	    the second atom. One or more shock waves will leave the second
	    atom, but there will be no reflected wave.  In each step, each
	    of the waves will travel  the distance between two neighboring 
	    atoms, i.e. the length of a single bond. The waves may branch,
	    but no reflection is allowed.  After the shock wave passes any
	    tunnel section (half of a single bond), this section collapses
	    so that  shock waves  cannot traverse it.  If any of the waves
	    reaches the atom  where  the whole process started  (the first
	    atom),  both the first atom and the second atom belong to some
	    closed structure (an aromatic ring or something like that). If
	    this is the case  no rotation  is allowed  about the specified
	    bond.

	(2) Step indices are assigned to  the atoms  which were reached by
	    shock waves.  The first  step index  (one) is assigned  to the
	    first input atom.  The second step index  (two) is assigned to
	    the second atom.  After that,  the whole residue is repeatedly
	    scanned, atom by atom.  For each atom,  all bonds are scanned.
	    The atom  is reached by the  shock wave in the current step if
	    the step index assigned to at least one  bond partner is equal
	    to the previous step index.

	(3) The auxiliaryI index should be reinitialized before every use.

========includes:============================================================*/

#include <stdio.h>

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <X11/Xatom.h>

#include "defines.h"
#include "typedefs.h"

/*======check is bond part of a closed structure:============================*/

int IsPartOfRing_ (MolComplexS *mol_complexSP, int  mol_complexI,
		   size_t atom1_arrayI, size_t atom2_arrayI)
{
MolComplexS		*curr_mol_complexSP;
size_t			atomsN;
int			residueI;
int			start_atomI, end_atomI;
int			stepI, previous_stepI;
AtomS			*curr_atomSP;
int			atoms_in_residueN;
int			propagationF;
int			atomI;
int			atom_reachedF;
int			bondI;
TrueBondS		*curr_bondSP;
AtomS			*curr_neighborSP;

/* Pointer to the chosen macromolecular complex: */
curr_mol_complexSP = mol_complexSP + mol_complexI;

/* Prepare and check the number of atoms: */
atomsN = curr_mol_complexSP->atomsN;
if (atomsN == 0) return -1;

/* Two specified atoms should belong to the same residue: */
residueI = (curr_mol_complexSP->atomSP + atom1_arrayI)->residue_arrayI;
if ((int) (curr_mol_complexSP->atomSP + atom2_arrayI)->residue_arrayI !=
								residueI)
	{
	return -2;
	}

/* If this point is reached, both atoms belong to the same residue. */

/* Prepare and check the atomic array indices which */
/* define the residue  to which  both atoms belong. */
/* This check might be redundant but I am paranoic. */
start_atomI = (curr_mol_complexSP->residueSP + residueI)->residue_startI;
end_atomI   = (curr_mol_complexSP->residueSP + residueI)->residue_endI;
if (((int) atom1_arrayI < start_atomI) || ((int) atom1_arrayI > end_atomI))
	{
	return -3;
	}
if (((int) atom2_arrayI < start_atomI) || ((int) atom2_arrayI > end_atomI))
	{
	return -4;
	}

/* Initialize the step index: */
stepI = 1;

/*------reset auxiliaryI for each atom of a chosen residue:------------------*/

/* Scan the chosen residue, atom by atom: */
for (atomI = start_atomI; atomI <= end_atomI; atomI++)
	{
	/* Pointer to the current atom: */
	curr_atomSP = curr_mol_complexSP->atomSP + atomI;

	/* Reset the auxiliaryI: */
	curr_atomSP->auxiliaryI = 0;
	}

/*------the shock wave starts from the first atom:---------------------------*/

/* Associate the first step index with the first atom: */
curr_atomSP = curr_mol_complexSP->atomSP + atom1_arrayI;
curr_atomSP->auxiliaryI = 1;

/*------the shock wave reaches the second atom:------------------------------*/

/* Associate the second step index with the second atom: */
curr_atomSP = curr_mol_complexSP->atomSP + atom2_arrayI;
curr_atomSP->auxiliaryI = 2;

/*------further propagation of shock waves:----------------------------------*/

/* Now scan the residue repeatedly. The number of scans should not */
/* exceed the number of atoms  in the residue,  because the length */
/* of the fully extended residue is equal to  number_of_atoms - 1. */
atoms_in_residueN = end_atomI - start_atomI + 1;

/* Initialize the index of the previous step: */
previous_stepI = 2;

/* The loop which counts steps should start from the third step: */
for (stepI = 3; stepI < atoms_in_residueN; stepI++)
	{
	/* Reset the wave propagation flag: */
	propagationF = 0;

	/* In each step, scan the entire residue, atom by atom: */
	for (atomI = start_atomI; atomI <= end_atomI; atomI++)
		{
		/* In the third step,  the first atom of */
		/* the specified bond should be skipped: */ 
		if ((stepI == 3) && (atomI == (int) atom1_arrayI)) continue;

		/* Pointer to the current atom: */
		curr_atomSP = curr_mol_complexSP->atomSP + atomI;

		/* If the current atom was reached */
		/* in the previous step,  skip it: */
		if (curr_atomSP->auxiliaryI == previous_stepI) continue;

		/* Reset the flag which says that the */
		/* atom is reached by the shock wave: */
		atom_reachedF = 0;

		/* Scan all bonds of the current atom: */
		for (bondI = 0; bondI < curr_atomSP->bondsN; bondI++)
			{
			/* Pointer to the current bond: */
			curr_bondSP = curr_atomSP->true_bondSA + bondI;

			/* Only the covalent bonds  are taken  into account. */
			/* Hydrogen, disulfide and pseudo-bonds are ignored. */
			if (curr_bondSP->bond_typeI != 1) continue;

			/* Bonds formed with atoms from another */
			/* macromolecular complex  are ignored: */
			if (curr_bondSP->neighbor_mol_complexI != mol_complexI)
				{
				continue;
				}

			/* If the bond partner was reached by the shock */
			/* wave in the previous step,  the current atom */
			/* is reached  by the shock wave  in this step: */
			curr_neighborSP = mol_complexSP->atomSP +
					  curr_bondSP->neighbor_arrayI;

			if (curr_neighborSP->auxiliaryI == previous_stepI)
				{
				/* Set the flag which says  that the */
				/* shock wave continues propagation: */
				propagationF = 1;

				/* Set the flag which says  that this */
				/* atom is reached by the shock wave: */
				atom_reachedF = 1;

				/* Do not inspect the remaining */
				/* bonds,  it is not necessary: */
				break;
				}
			}

		/* If the current atom is reached in this step: */
		if (atom_reachedF)
			{
			/* If the current atom is  the first atom of the */
			/* specified bond and if it was reached  in this */
			/* step, the bond is involved in ring formation: */
			if (atomI == (int) atom1_arrayI) return 1;

			/* If the shock wave reached the current atom */
			/* for the first time,  store the step index: */
			if (curr_atomSP->auxiliaryI == 0)
				{
				curr_atomSP->auxiliaryI = stepI;
				}
			}
		}

	/* If the shock wave is not propagating any more, break: */
	if (!propagationF) break;

	/* Update the value of the previous step index: */
	previous_stepI = stepI;
	}

/*---------------------------------------------------------------------------*/

/* Return zero if this point is reached: */
return 0;
}

/*===========================================================================*/