/* Copyright (C) 2003 Damir Zucic */

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

				beta_membrane.c

Purpose:
	Find the position and orientation of the membrane with respect to
	the macromolecular structure. The structure in question should be
	some membrane protein of beta barrel type. This function wil fail
	to work properly for alpha helix bundle proteins.

Input:
	(1) Pointer to MolComplexS structure,  with the chosen structure.
	(2) Pointer to RuntimeS structure.
	(3) Pointer to ConfigS structure (stereo angle required).
	(4) Pointer to GUIS structure.

Output:
	(1) The membrane position and orientation  will be calculated for
	    the first time or updated.
	(2) Return value.

Return value:
	(1) Positive on success.
	(2) Negative on failure.  
Notes:
	(1) A short description of some tricks used in this function:
	    -------------------------------------------------------------
	    Step 1: Find the geometric center for CA atoms. Ignore hetero
	    atoms (water, lipids, cofactors etc.). Ignore all other atoms
	    too (C, N, O, side chain atoms etc.).
	    -------------------------------------------------------------
	    Step 2: Use the CA geometric center and some arbitrary vector
	    to create an axis in space.  Associate  a cylinder  with this
	    axis. The length of this cylinder should be 120  angstroms and
	    the radius may be arbitrary. Divide the cylinder into patches
	    of equal size.  The height of  a single patch  should be  3.0
	    angstroms  (see the parameter CELL_LINEAR_WIDTH in defines.h)
	    and the patch angular width should  correspond to the angular
	    range of 36 degrees. Note that the radius of the cylinder and
	    the width of  a single patch  are not important:  what really
	    matters is the angular range of a single patch.  Associate an
	    array of cells (data structures) with this cylinder. Scan the
	    protein structure,  residue by residue.  For a given residue,
	    check the CA-CB vector.  If this  vector points  inwards with
	    respect to the  given axis,  ignore this residue.  Otherwise,
	    project the  CB atom position to the surface of the cylinder.
	    Identify the patch to which the  CB atom was projected.  Find
	    the distance  between the  CB  atom and  the axis.  Check the
	    distance which was stored before to  the cell associated with
	    the patch to which the CB atom was projected.  If the current
	    CB is more distant from  the axis than the  CB atom which was
	    the most distant atom projected before to this patch,  occupy
	    the cell with the current  CB atom:  store the hydrophobicity
	    of the current residue and  the distance between  CB  and the
	    axis.
	    -------------------------------------------------------------
	    Step 3: Divide the cylinder surface into stripes, parallel to
	    the cylinder surface. For each cell in each stripe, calculate
	    the average hydrophobicity over a given number of neighboring
	    cells.  These cells  should belong to  the same  stripe.  The
	    sliding window should be comparable to the membrane thickness
	    (or slightly smaller). The unused cells should be ignored. It
	    is possible that average hydrophobicity will not be available
	    for all cells after this step,  but this should not cause any
	    major problems.  Most of the cells in the central part of the
	    cylinder surface will be filled with  average hydrophobicity.
	    -------------------------------------------------------------
	    Step 4: Divide the cylinder surface into rings. In this step,
	    the averaged hydrophobicities  will be used  to calculate the
	    second average,  this time  combining values  from  the cells
	    which belong to the same ring. This procedure is repeated for
	    each ring.  Ignore cells which  were skipped in  the previous
	    step,  i.e.  which  are  missing  the average  hydrophobicity
	    value.
	    -------------------------------------------------------------
	    Step 5:  Find  the most  hydrophobic ring.  Store  the second
	    average of hydrophobicity for later usage. The same procedure
	    is repeated  for a  number of axes,  to cover  the full space
	    angle of  4 PI  steradijans.  Find  the axis which  gives the
	    highest second average of hydrophobicity. Use the unit vector
	    of this axis as the initial unit vector  perpendicular to the
	    membrane. Together with the CA center, this is just enough to
	    define the membrane.
	    -------------------------------------------------------------
	    Step 6: Refine the orientation of the normal vector. Look for
	    two rings of aromatic residues, PHE, TRP and TYR. Ignore HIS.
	    The transmembrane part of  the beta barrel protein has a ring
	    of aromatic residues PHE, TRP and TYR at each end. Ignore all
	    residues pointing outwards with respect to the axis.  Project
	    eacg CG atom to a number of cells, covering few angstroms.
	    -------------------------------------------------------------
	    Step 7: Refine the position of the membrane center. Up to now
	    the position of CA center was used as the membrane center but
	    this is obviously wrong. Project CG atoms of PHE, TRP and TYR
	    over a number of cells, so that each CG covers few angstroms.
	    -------------------------------------------------------------

	(2) The minimal distance between the point (x, y, z) and the axis
	    passing through the point (x0, y0, z0), pointing in direction
	    (l, m, n) is:
                                        
	                      1
	     d^2   =   ---------------  { [ (x - x0) m - (y - y0) l ]^2 +
                       l^2 + m^2 + n^2
                                          [ (y - y0) n - (z - z0) m ]^2 +

                                          [ (z - z0) l - (x - x0) n ]^2 }

	    Using  x - x0 = dx, y - y0 = dy, z - z0 = dz,  l = ex, m = ey
	    and n = ez,  where ex, ey and ez  are three components of the
	    axis unit vector, the equation looks simpler:

	    distance_squared  = (dx * ey - dy * ex)^2 +
				(dy * ez - dz * ey)^2 +
				(dz * ex - dx * ez)^2 .

	    The formula  was taken from  the famous mathematical handbook
	    by I.N. Bronstein and K.A. Semendjajew. However,  it was used
	    only for debugging and later it was removed.  The distance is
	    now calculated using  the absolute value of the perpendicular
	    component of  the  CA_center_CB_vectorS.  The expression from
	    the book was used only to do one extra check.

	(3) Two strange terms are used here: beta cells and simple cells.
	    The auxiliary cylinder  is divided into  "beta cells" and the
	    axis used for refinement is divided into "simple cells".

	(4) Indentation is exceptionally 4 spaces.

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

#include <stdio.h>

#include <string.h>
#include <math.h>

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

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

/*======function prototypes:=================================================*/

int		ExtractCA_ (VectorS *, AtomS *, size_t, size_t);
int		ExtractCACB_ (VectorS *, VectorS *, AtomS *, size_t, size_t);
double		ScalarProduct_ (VectorS *, VectorS *);
void		VectorProduct_ (VectorS *, VectorS *, VectorS *);
double		AbsoluteValue_ (VectorS *);
int		ExtractCG_ (VectorS *, AtomS *, size_t, size_t);

/*======find membrane position (beta barrel version):========================*/

int BetaMembrane_ (MolComplexS *mol_complexSP,
		   RuntimeS *runtimeSP, ConfigS *configSP, GUIS *guiSP)
{
double		membrane_thickness;
int		simple_cells_along_axisN;
double		simple_cell_width, recip_simple_width;
double		projection_width;
int		simple_window_width, half_simple_window_width;
double		scan_half_width;
int		simple_scan_half_width;
int		max_simpleI, simple_thickness, half_simple_thickness;
int		beta_cells_along_axisN, cells_in_ringN, half_cells_in_ringN;
double		beta_cell_linear_width, recip_beta_lin_width;
double		beta_cell_angular_width, recip_beta_ang_width;
int		half_window_width, window_width;
int		beta_index_offset, simple_index_offset;
size_t		atomsN;
int		residuesN, residueI;
ResidueS	*residueSP;
int		n;
VectorS		CA_vectorS;
VectorS		CA_center_vectorS;
int		CA_atomsN;
double		reciprocal_denominator;
int		thetaI, max_thetaI, phiI, max_phiI;
double		theta, theta_step, phi, phi_step;
VectorS		axis_unit_vectorS;
VectorS		aux_vectorS;
double		scalar_product;
double		abs_value, reciprocal_abs_value;
VectorS		unit_vector2S, unit_vector3S;
int		beta_cellI;
BetaCellS	*beta_cellSP;
VectorS		CB_vectorS;
VectorS		CA_CB_vectorS;
VectorS		CA_center_CB_vectorS;
VectorS		parallel_vectorS;
VectorS		perpendicular_vectorS;
double		distance;
double		sin_angle, cos_angle, angle;
int		axialI, radialI;
int		localI, combined_radialI;
AtomS		*first_atomSP;
double		hydrophobicity, average_hydrophobicity;
int		cells_usedN;
int		windowI;
int		combined_axialI;
double		second_average;
int		best_second_average_foundF;
double		best_second_average;
VectorS		best_unit_vectorS;
VectorS		i1_unit_vectorS, j1_unit_vectorS, k1_unit_vectorS;
int		total_phe_trp_tyrN;
char		*pure_residue_nameP;
double		cos_theta, sin_theta, cos_phi, sin_phi;
int		two_rings_count, best_two_rings_count;
int		central_simpleI, start_simpleI, end_simpleI, simpleI;
int		phe_trp_tyrNA[SIMPLE_CELLS_ALONG_AXIS];
VectorS		CG_vectorS;
VectorS		CA_center_CG_vectorS;
int		start_pointI, end_pointI, pointI;
double		shift, best_shift;
VectorS		center_vectorS;
VectorS		radius_vectorS;
double		projection;
int		ring1I, ring2I;

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

/* Copy the membrane thickness: */
membrane_thickness = mol_complexSP->membraneS.thickness;

/* Copy the number of simple cells along the axis: */
simple_cells_along_axisN = SIMPLE_CELLS_ALONG_AXIS;

/* Prepare the linear width of a single simple cell: */
simple_cell_width = SIMPLE_CELL_WIDTH;
recip_simple_width = 1.0 / simple_cell_width;

/* Prepare the projection width. The projection of */
/* CG atom will be smeared over a number of cells. */
projection_width = 5.0;              /* angstroms. */

/* Prepare the width of a simple window.  A single CG atom will */
/* be projected over a number of cells. The simple_window_width */
/* defines how many cells  will be covered by a single CG atom. */
simple_window_width = recip_simple_width * projection_width;
half_simple_window_width = simple_window_width / 2;

/* During the refinement of  the membrane center,  a range of positions */
/* along the refined unit vector should be scanned. It is defined here: */
scan_half_width = 100.0;
simple_scan_half_width = (int) (scan_half_width * recip_simple_width);

/* The initial number of simple cells across  the given membrane thickness. */
/* This value will be changed before the refinement of the membrane center. */
simple_thickness = (int) (membrane_thickness / simple_cell_width);
if (simple_thickness < 3) simple_thickness = 3;
half_simple_thickness = simple_thickness / 2;

/* Maximal value of the simple index. The simple index is */
/* used during the refinement. The maximal value set here */
/* should prevent the array overflow. Note that ten cells */
/* are wasted, to ensure there will be no array overflow. */
max_simpleI = simple_cells_along_axisN - 1 - simple_thickness - 10;
if (max_simpleI < 0) max_simpleI = 0;

/* Copy the number of beta cells along the axis of the auxiliary cylinder: */
beta_cells_along_axisN = BETA_CELLS_ALONG_AXIS;

/* Copy the number of cells in a single ring around the auxiliary cylinder: */
cells_in_ringN = BETA_CELLS_IN_RING;
half_cells_in_ringN = cells_in_ringN / 2;

/* Prepare the linear width of a single beta cell: */
beta_cell_linear_width = BETA_CELL_LINEAR_WIDTH;
recip_beta_lin_width = 1.0 / beta_cell_linear_width;

/* Prepare the angular width of a single beta cell: */
beta_cell_angular_width = 6.2831853 / (double) cells_in_ringN;
recip_beta_ang_width = 1.0 / beta_cell_angular_width;

/* Prepare the sliding window width.  This number should */
/* be odd. It should not reach the hydrophilic residues. */
n = (int) (membrane_thickness / beta_cell_linear_width);
half_window_width = (n - 1) / 2;
if (half_window_width < 1) return -1;
window_width = 2 * half_window_width + 1;

/* The height of the CA center with respect to the base of the */
/* cylinder should reach the middle of the auxiliary cylinder. */
/* To ensure this, the axial index should be shifted properly: */
beta_index_offset = beta_cells_along_axisN / 2;

/* During the refinement, the CA center should be */
/* projected to the central cell  along the axis. */
simple_index_offset = simple_cells_along_axisN / 2;

/* Copy and check the number of atoms: */
atomsN = mol_complexSP->atomsN;
if (atomsN == 0) return -2;

/* Copy and check the number of residues: */
residuesN = mol_complexSP->residuesN;
if (residuesN < window_width) return -3;

/*------(1) find the geometric center for CA atoms:-------------------------*/

/* Initialize the CA center: */
CA_center_vectorS.x = 0.0;
CA_center_vectorS.y = 0.0;
CA_center_vectorS.z = 0.0;

/* Initialize the number of CA atoms: */
CA_atomsN = 0;

/* Scan the macromol. complex, residue by residue, ignoring water molecules: */
for (residueI = 0; residueI < residuesN; residueI++)
    {
    /* Pointer to the current residue: */
    residueSP = mol_complexSP->residueSP + residueI;

    /* Try to extraxt the coordinates of CA atom: */
    n = ExtractCA_ (&CA_vectorS, mol_complexSP->atomSP,
		    residueSP->residue_startI,
		    residueSP->residue_endI);
    if (n < 1) continue;

    /* Update the CA geometric center: */
    CA_center_vectorS.x += CA_vectorS.x;
    CA_center_vectorS.y += CA_vectorS.y;
    CA_center_vectorS.z += CA_vectorS.z;

    /* Update the counter of CA atoms: */
    CA_atomsN++;
    }

/* If there were no CA atoms at all, return error indicator: */
if (CA_atomsN == 0) return -4;

/* Calculate the geometric center of CA atoms: */
reciprocal_denominator = 1.0 / (double) CA_atomsN;
CA_center_vectorS.x *= reciprocal_denominator;
CA_center_vectorS.y *= reciprocal_denominator;
CA_center_vectorS.z *= reciprocal_denominator;

/* Take this point as the initial membrane center: */
mol_complexSP->membraneS.center_x = CA_center_vectorS.x;
mol_complexSP->membraneS.center_y = CA_center_vectorS.y;
mol_complexSP->membraneS.center_z = CA_center_vectorS.z;

/*------scan a number of directions in space:--------------------------------*/

/* Initialize the parameters which define a total number of directions: */
max_thetaI = 18;
max_phiI   = 72;

/* Prepare the angular increments. Note that only half of the full theta */
/* angle is used:  we do not distinguish two sides of the membrane here. */
theta_step = 1.5707963 / (double) max_thetaI;
phi_step   = 6.2831853 / (double) max_phiI;

/* Initialize the flag which will be  set to one if the */
/* best second average of hydrophobicity will be found: */
best_second_average_foundF = 0;

/* Initialize the best second average, just to avoid compiler complaints: */
best_second_average = -99999.00;

/* Initialize the best axis unit vector, for any case: */
best_unit_vectorS.x = 0.0;
best_unit_vectorS.y = 1.0;
best_unit_vectorS.z = 0.0;

/* Scan a number of directions in space. The spherical polar system is used. */

/* Theta scan: */
for (thetaI = 0; thetaI <= max_thetaI; thetaI++)
    {
    /* Prepare the angle theta: */
    theta = (double) thetaI * theta_step;

    /* Phi scan: */
    for (phiI = 0; phiI <= max_phiI; phiI++)
	{
	/* Avoid to check too many directions close */
	/* to the pole. If the theta angle is equal */
	/* to zero, one phi value should be enough: */
	if (thetaI == 0)
	    {
	    if (phiI != 0) continue;
	    }

	/* If theta is too small, skip every second direction: */
	if (thetaI <= 2)
	    {
	    n = phiI / 2;
	    if (n * 2 != phiI) continue;
	    }

	/* Prepare the angle phi: */
	phi = (double) phiI * phi_step;

	/* The unit vector which defines the axis direction: */
	axis_unit_vectorS.x = sin (theta) * cos (phi);
	axis_unit_vectorS.y = sin (theta) * sin (phi);
	axis_unit_vectorS.z = cos (theta);

	/* Two additional vectors  are associated  with each axis: */
	/* unit_vector2S and unit_vector3S. These two vectors will */
	/* be used to define  which cell is the first one  in each */
	/* ring of cells.  Without  these two vectors,  it will be */
	/* unclear to which cell the CB atom  should be projected. */

	/* To prepare the  unit_vector2S,  use the unit vector parallel */
	/* to the x axis to define the plane. If it is too close to the */
	/* cylinder axis use  the unit vector  parallel to the  y axis. */

	/* The first choice of auxiliary vector: */
	aux_vectorS.x = 1.0;
	aux_vectorS.y = 0.0;
	aux_vectorS.z = 0.0;

	/* Check the scalar product with the axis unit vector.  Both */
	/* vectors are normalized,  so in the worst case  the scalar */
	/* product will be equal to one.  A value above 0.9 is taken */
	/* as a threshold to recognize the bad pair of unit vectors. */
	if (ScalarProduct_ (&aux_vectorS, &axis_unit_vectorS) > 0.9)
	    {
	    /* The second choice: */
	    aux_vectorS.x = 0.0;
	    aux_vectorS.y = 1.0;
	    aux_vectorS.z = 0.0;
	    }

	/* Now it is certain that  axis_unit_vectorS  and */
	/* aux_vectorS are  not  parallel.  Calculate and */
	/* normalize the vector product of these vectors. */ 
	VectorProduct_ (&unit_vector2S, &axis_unit_vectorS, &aux_vectorS);
	abs_value = AbsoluteValue_ (&unit_vector2S);
	if (abs_value == 0.0) continue;  /* Yes, I am paranoic! */
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector2S.x *= reciprocal_abs_value;
	unit_vector2S.y *= reciprocal_abs_value;
	unit_vector2S.z *= reciprocal_abs_value;

	/* Prepare  the third unit vector.  As both  axis_unit_vectorS */
	/* and unit_vector2S  are normalized and  mutually orthogonal, */
	/* the vector product of these vectors will be normalized too. */
	VectorProduct_ (&unit_vector3S, &axis_unit_vectorS, &unit_vector2S);

	/* Reset the array of BetaCellS structures: */
	for (beta_cellI = 0; beta_cellI < runtimeSP->beta_cellsN; beta_cellI++)
	    {
	    beta_cellSP = runtimeSP->beta_cellSP + beta_cellI;
	    beta_cellSP->cell_usedF = 0;
	    beta_cellSP->distance = 0.0;
	    beta_cellSP->hydrophobicity = 0.0;
	    beta_cellSP->cells_usedN = 0;
	    beta_cellSP->average_calculatedF = 0;
	    beta_cellSP->average_hydrophobicity = 0.0;
	    }

	/*------(2) project surface CB atoms to cylinder surface:------------*/

	/* Scan residues: */
	for (residueI = 0; residueI < residuesN; residueI++)
	    {
	    /* Pointer to the current residue: */
	    residueSP = mol_complexSP->residueSP + residueI;

	    /* Try to extract the CA and CB coordinates: */
	    n = ExtractCACB_ (&CA_vectorS, &CB_vectorS,
			      mol_complexSP->atomSP,
			      residueSP->residue_startI,
			      residueSP->residue_endI);
	    if (n < 2) continue;

	    /* Prepare the vector from CA to CB: */
	    CA_CB_vectorS.x = CB_vectorS.x - CA_vectorS.x;
	    CA_CB_vectorS.y = CB_vectorS.y - CA_vectorS.y;
	    CA_CB_vectorS.z = CB_vectorS.z - CA_vectorS.z;

	    /* The vector from the CA geometric center */
	    /* to the CB atom of  the current residue: */
	    CA_center_CB_vectorS.x = CB_vectorS.x - CA_center_vectorS.x;
	    CA_center_CB_vectorS.y = CB_vectorS.y - CA_center_vectorS.y;
	    CA_center_CB_vectorS.z = CB_vectorS.z - CA_center_vectorS.z;

	    /* The scalar product  between this */
	    /* vector and the axis unit vector: */
	    scalar_product = ScalarProduct_ (&CA_center_CB_vectorS,
					     &axis_unit_vectorS);

	    /* The  parallel component of  the vector  from the */
	    /* CA center to the CB atom of the current residue: */
	    parallel_vectorS.x = scalar_product * axis_unit_vectorS.x;
	    parallel_vectorS.y = scalar_product * axis_unit_vectorS.y;
	    parallel_vectorS.z = scalar_product * axis_unit_vectorS.z;

	    /* The perpendicular component: */
	    perpendicular_vectorS.x = CA_center_CB_vectorS.x -
				      parallel_vectorS.x;
	    perpendicular_vectorS.y = CA_center_CB_vectorS.y -
				      parallel_vectorS.y;
	    perpendicular_vectorS.z = CA_center_CB_vectorS.z -
				      parallel_vectorS.z;

	    /* Calculate and check the distance  from the axis. */
	    /* Skip residues whose  CB  atoms are  too close to */
	    /* the axis, they do not belong to the barrel wall. */
	    /* A limit of  eight angstroms might be reasonable. */
	    distance = AbsoluteValue_ (&perpendicular_vectorS);
	    if (distance <= 8.0) continue;

	    /* Check the scalar product between the CA-CB vector and */
	    /* the perpendicular component of the vector from the CA */
	    /* center to the  CB atom of the current residue.  If it */
	    /* is negative the CA-CB vector points towards the axis. */
	    scalar_product = ScalarProduct_ (&CA_CB_vectorS,
					     &perpendicular_vectorS);
	    if (scalar_product < 0.0) continue;

	    /* If this point is reached,  the CA-CB bond */
	    /* points outwards with respect to the axis. */

	    /* Project the  CB coordinates to the */
	    /* surface of the auxiliary cylinder. */

	    /* To find the cell to which CB should be projected, */
	    /* two indices are required: axial and radial index. */

	    /* The parallel part of the vector from the CA center */
	    /* to the CB atom will be used to calculate the axial */
	    /* index of the cell  to which the  CB atom should be */
	    /* projected. The perpendicular component of the same */
            /* vector will be used to calculate the radial index. */

	    /* Axial index  is calculated  from  the projection of the */
	    /* parallel component of CA_center_CB_vectorS to the axis: */
	    scalar_product = ScalarProduct_ (&parallel_vectorS,
					     &axis_unit_vectorS);
	    axialI = (int) (recip_beta_lin_width * scalar_product) +
		     beta_index_offset;

	    /* Check the axial index, to avoid the arary overflow: */
	    if (axialI < 0) continue;
	    if (axialI >= beta_cells_along_axisN) continue;

	    /* The angle in radians  is required to calculate the */
	    /* radial index. The angular range if from 0 to 2 PI. */

	    /* Calculate the arc cosine using the scalar product between the */
	    /* perpendicular component of the  CA_center_CB_vectorS  and the */
	    /* unit_vector2S.  The result will be in the range from 0 to PI. */
	    scalar_product = ScalarProduct_ (&perpendicular_vectorS,
					     &unit_vector2S);
	    cos_angle = scalar_product / distance;
	    if      (cos_angle <= -1.0) angle = 3.1415927;
	    else if (cos_angle >=  1.0) angle = 0.0;
	    else angle = acos (cos_angle);

	    /* The  scalar product  between  the perpendicular */
	    /* component  of the  CA_center_CB_vectorS and the */
	    /* unit_vector3S is used to resolve the ambiguity. */
	    /* If it is negative,  the angle  should be fixed. */
	    scalar_product = ScalarProduct_ (&perpendicular_vectorS,
					     &unit_vector3S);
	    if (scalar_product < 0.0) angle = 6.2831853 - angle;

	    /* Now it is possible to calculate the radial index: */
	    radialI = recip_beta_ang_width * angle;

	    /* Check the radial index, to avoid the arary overflow: */
	    if (radialI < 0) continue;
	    if (radialI >= cells_in_ringN) continue;

	    /* Check the cell defined by  the given pair */
	    /* of indices and two  adjacent cells in the */
	    /* plane which is perpendicular to the axis. */

	    /* Scan the window of three cells: */
	    for (localI = -1; localI <= 1; localI++)
		{
		/* Prepare and check the combined index: */
		combined_radialI = radialI + localI;
		if (combined_radialI < 0)
		    {
		    combined_radialI = cells_in_ringN - 1;
		    }
		else if (combined_radialI >= cells_in_ringN)
		    {
		    combined_radialI = 0;
		    }

		/* Now the cell to which the  CB  atom should be */
		/* projected is  identified.  Check was it  used */
		/* before.  If not,  just store  the information */
		/* about this CB atom.  If it was used,  compare */
		/* the distance from  the axis of this  CB  atom */
		/* with  the value  which was  stored  before to */
		/* the cell  defined by  the  given  axialI  and */
		/* combined_radialI.  If this CB atom is not the */
		/* most distant CB atom  projected to this cell, */
		/* ignore it. If it is, change the cell content. */

		/* The index of  BetaCellsS  structure to */
		/* which the CB atom should be projected. */
		/* The array  is indexed  ring  by  ring. */
		beta_cellI = axialI * cells_in_ringN + combined_radialI;

		/* A true paranoic will make some extra checks: */
		if (beta_cellI < 0) continue;
		if (beta_cellI >= runtimeSP->beta_cellsN) continue;

		/* The pointer to this BetaCellS structure: */
		beta_cellSP = runtimeSP->beta_cellSP + beta_cellI;

		/* If this cell was not used before, just fill it: */
		if (beta_cellSP->cell_usedF == 0)
		    {
		    beta_cellSP->cell_usedF = 1;
		    beta_cellSP->distance = distance;
		    first_atomSP = mol_complexSP->atomSP +
				   residueSP->residue_startI;
		    hydrophobicity = first_atomSP->raw_atomS.hydrophobicity;
		    beta_cellSP->hydrophobicity = hydrophobicity;
		    }

		/* If it was used before, compare the old and new distance: */
		else
		    {
		    if (beta_cellSP->distance > distance) continue;

		    /* If this point is reached, the given CB atom is the */
		    /* most distant CB atom  projected to the given cell. */

		    beta_cellSP->distance = distance;
		    first_atomSP = mol_complexSP->atomSP +
				   residueSP->residue_startI;
		    hydrophobicity = first_atomSP->raw_atomS.hydrophobicity;
		    beta_cellSP->hydrophobicity = hydrophobicity;
		    }

		/* End of localI loop: */
		}

	    /* End of residueI loop: */
	    }

	/*------(3) calculate average hydrophob. (axial scan):---------------*/

	/* Now comes the averaging. */

	/* Scan  the cylinder surface,  stripe by stripe.  All  stripes  are */
	/* parallel to the axis. The sliding window width (used to calculate */
	/* the average  hydrophobicity  with a given residue  in the center) */
	/* is set to cover  the entire membrane thickness  (but no more than */
	/* that).  For example,  if membrane thickness is  30  angstroms the */
	/* window width will be about 10, because in axis direction one cell */
	/* covers 3 angstroms. Odd number is recommended, so the final width */
	/* is 9.  If 11 cells were used hydrophilic residues may be reached. */

	/* The outermost loop scans the surface, stripe by stripe: */
	for (radialI = 0; radialI < cells_in_ringN; radialI++)
	    {

	    /* The next loop scans the given stripe, cell by cell. The cells */
	    /* which are too close to either top or bottom edge are skipped. */
	    for (axialI = 0; axialI < beta_cells_along_axisN; axialI++)
		{
		if (axialI < half_window_width) continue;
		if (axialI + half_window_width >= beta_cells_along_axisN)
		    {
		    continue;
		    }

		/* Reset the number of useful cells in the sliding window: */
		cells_usedN = 0;

		/* Reset the average hydrophobicity: */
		average_hydrophobicity = 0.0;

		/* Somewhere above the flag average_calculatedF */
		/* was cleared (set to zero) for all cells.  It */
		/* is not necessary to clear these flags again. */

		/* The innermost loop scans the sliding */
		/* window centered  at  the given cell: */
		for (windowI = 0; windowI < window_width; windowI++)
		    {
		    /* Prepare the combined axial index: */
		    combined_axialI = axialI + windowI - half_window_width;

		    /* The cell array index: */
		    beta_cellI = combined_axialI * cells_in_ringN + radialI;

		    /* Pointer to the current cell: */
		    beta_cellSP = runtimeSP->beta_cellSP + beta_cellI;

		    /* If this cell was unused, skip it: */
		    if (beta_cellSP->cell_usedF == 0) continue;

		    /* If this point is reached, this cell is useful. */

		    /* Add the hydrophobicity of this cell to the total: */
		    average_hydrophobicity += beta_cellSP->hydrophobicity;

		    /* Update the counter: */
		    cells_usedN++;
		    }

		/* If  there were  too many  unused cells in  this sliding */
		/* window,  the cell in the center of  this sliding window */
		/* should be treated as useless.  The cell will be treated */
		/* as useful if  at least  half_window_width + 1  cells in */
		/* the sliding window associated with this cell were used. */
		if (cells_usedN < half_window_width + 1) continue;

		/* Divide the total hydrophobicity by the number of cells: */
		average_hydrophobicity /= (double) cells_usedN;

		/* Store the average hydrophobicity and the number */
		/* of used cells and set  the flag which says that */
		/* average hydrophob. was successfully calculated. */
		/* Store it to  the cell  which was  in the middle */
		/* of  the sliding window  which was just scanned. */
                beta_cellI = axialI * cells_in_ringN + radialI; 
                beta_cellSP = runtimeSP->beta_cellSP + beta_cellI;
		beta_cellSP->average_hydrophobicity = average_hydrophobicity;
		beta_cellSP->cells_usedN = cells_usedN;
                beta_cellSP->average_calculatedF = 1;
		}
	    }

        /*------(4) calculate average hydrophob. (radial scan):--------------*/

	/* The outermost loop scans the surface, ring by ring: */
	for (axialI = 0; axialI < beta_cells_along_axisN; axialI++)
	    {
	    /* Reset the number of useful cells in the ring: */
	    cells_usedN = 0;

	    /* Reset the second average of hydrophobicity: */
	    second_average = 0.0;

	    /* The next loop  scans  the entire ring  at the */
	    /* given height defined by axialI, cell by cell: */
	    for (radialI = 0; radialI < cells_in_ringN; radialI++)
		{
		/* The cell array index: */
		beta_cellI = axialI * cells_in_ringN + radialI;

		/* Pointer to this cell: */
		beta_cellSP = runtimeSP->beta_cellSP + beta_cellI;

		/* If the average hydrophobicity is not */
		/* available  for  this cell,  skip it: */
		if (beta_cellSP->average_calculatedF == 0) continue;

		/* If this point is reached, this cell is useful. */

		/* Add the hydrophobicity of this cell to the total: */
		second_average += beta_cellSP->average_hydrophobicity;

		/* Update the counter: */
		cells_usedN++;
		}

	    /* At least half of the cells in this ring should be */
	    /* useful,  otherwise the second average is useless: */
	    if (cells_usedN < half_cells_in_ringN) continue;

	    /* Divide the total by the number of cells: */
	    second_average /= (double) cells_usedN;

	    /* If this is the first second average at all, just copy it: */
	    if (best_second_average_foundF == 0)
		{
		best_second_average = second_average;
		best_second_average_foundF = 1;
		best_unit_vectorS.x = axis_unit_vectorS.x;
		best_unit_vectorS.y = axis_unit_vectorS.y;
		best_unit_vectorS.z = axis_unit_vectorS.z;
		continue;
		}

	    /*------(5) check is it the most hydrophobic ring:---------------*/

	    /* Now check is there any chance that this value is the */
            /* highest value of the second average up to now. If it */
	    /* is,  update the best value and  store theta and phi: */
	    if (second_average > best_second_average)
		{
		best_second_average = second_average;
		best_second_average_foundF = 1;
		best_unit_vectorS.x = axis_unit_vectorS.x;
		best_unit_vectorS.y = axis_unit_vectorS.y;
		best_unit_vectorS.z = axis_unit_vectorS.z;
		}
	    }

	/* End of phiI loop: */
	}

    /* End of thetaI loop: */
    }

/*------store the normal vector (raw):---------------------------------------*/

/* Store (copy) the components of the normal */
/* vector because  the refinement  may fail. */

/* The normal vector of the first plane points to the same direction: */
mol_complexSP->membraneS.plane1S.normal_x[0] =  best_unit_vectorS.x;
mol_complexSP->membraneS.plane1S.normal_y    =  best_unit_vectorS.y;
mol_complexSP->membraneS.plane1S.normal_z[0] =  best_unit_vectorS.z;

/* The normal vector of the second plane points to the opposite direction: */
mol_complexSP->membraneS.plane2S.normal_x[0] = -best_unit_vectorS.x;
mol_complexSP->membraneS.plane2S.normal_y    = -best_unit_vectorS.y;
mol_complexSP->membraneS.plane2S.normal_z[0] = -best_unit_vectorS.z;

/*------prepare stereo data (raw):-------------------------------------------*/

/* Prepare the stereo data for both normal vectors. */
/* This is useful because  the refinement may fail. */

/* Prepare the sine and cosine of the stereo angle: */
sin_angle = sin (configSP->stereo_angle);
cos_angle = cos (configSP->stereo_angle);

/* Calculate the stereo data: */
mol_complexSP->membraneS.plane1S.normal_x[1] =
                 mol_complexSP->membraneS.plane1S.normal_x[0] * cos_angle +
                 mol_complexSP->membraneS.plane1S.normal_z[0] * sin_angle;
mol_complexSP->membraneS.plane1S.normal_z[1] =
                -mol_complexSP->membraneS.plane1S.normal_x[0] * sin_angle +
                 mol_complexSP->membraneS.plane1S.normal_z[0] * cos_angle;
mol_complexSP->membraneS.plane2S.normal_x[1] =
                 mol_complexSP->membraneS.plane2S.normal_x[0] * cos_angle +
                 mol_complexSP->membraneS.plane2S.normal_z[0] * sin_angle;
mol_complexSP->membraneS.plane2S.normal_z[1] =
                -mol_complexSP->membraneS.plane2S.normal_x[0] * sin_angle +
                 mol_complexSP->membraneS.plane2S.normal_z[0] * cos_angle;

/*------(6) search two rings of PHE, TRP and TYR:----------------------------*/

/* In this procedure,  a new  coordinate system  will be used.  Three axes */
/* of  this system  are defined by  i1_unit_vectorS,  j1_unit_vectorS  and */
/* k1_unit_vectorS. The best unit vector found in the procedure above will */
/* be used to define the z axis,  i.e. it will be used as k1_unit_vectorS. */

/* Prepare the unit vector which define the z axis of a new system: */
k1_unit_vectorS.x = best_unit_vectorS.x;
k1_unit_vectorS.y = best_unit_vectorS.y;
k1_unit_vectorS.z = best_unit_vectorS.z;

/* To prepare  the unit vector  which  defines  the y axis of  a new */
/* system,  use the unit vector parallel to  the x axis of  the main */
/* coordinate system to define the plane.  If it is too close to new */
/* z axis use the unit vector parallel to y axis of the main system. */

/* The first choice of auxiliary vector: */
aux_vectorS.x = 1.0;
aux_vectorS.y = 0.0;
aux_vectorS.z = 0.0;

/* Check the scalar product with k1_unit_vectorS. Both these */
/* vectors are normalized,  so in the worst case  the scalar */
/* product will be equal to one.  A value above 0.9 is taken */
/* as a threshold to recognize the bad pair of unit vectors. */
if (ScalarProduct_ (&aux_vectorS, &k1_unit_vectorS) > 0.9)
    {
    /* The second choice: */
    aux_vectorS.x = 0.0;
    aux_vectorS.y = 1.0;
    aux_vectorS.z = 0.0;
    }

/* Now it is certain that k1_unit_vectorS and aux_vectorS are */
/* not parallel.  Calculate and normalize the vector product. */
VectorProduct_ (&j1_unit_vectorS, &k1_unit_vectorS, &aux_vectorS);
abs_value = AbsoluteValue_ (&j1_unit_vectorS);
if (abs_value == 0.0) return -5;	/* Yes, I am paranoic! */
reciprocal_abs_value = 1.0 / abs_value;
j1_unit_vectorS.x *= reciprocal_abs_value;
j1_unit_vectorS.y *= reciprocal_abs_value;
j1_unit_vectorS.z *= reciprocal_abs_value;

/* The third unit vector, parallel to x axis of a new system: */
VectorProduct_ (&i1_unit_vectorS, &j1_unit_vectorS, &k1_unit_vectorS);

/* Initialize the total number of  PHE, TRP and */
/* TYR residues in this macromolecular complex. */
total_phe_trp_tyrN = 0;

/* Scan the macromolecular complex, residue by residue.  For each */
/* residue, check is it PHE, TRP or TYR. If it is, set auxiliaryI */ 
/* of the first atom to one.  If not, set the auxiliaryI to zero. */
for (residueI = 0; residueI < residuesN; residueI++)
    {
    /* Pointer to the current residue: */
    residueSP = mol_complexSP->residueSP + residueI;

    /* Pointer to the first atom: */
    first_atomSP = mol_complexSP->atomSP + residueSP->residue_startI;

    /* The purified residue name (without spaces): */
    pure_residue_nameP = first_atomSP->raw_atomS.pure_residue_nameA;

    /* Check the purified residue name: */
    if (strcmp (pure_residue_nameP, "PHE") == 0)
	{
	first_atomSP->auxiliaryI = 1;
	total_phe_trp_tyrN++;
	}
    else if (strcmp (pure_residue_nameP, "TRP") == 0)
	{
	first_atomSP->auxiliaryI = 2;
	total_phe_trp_tyrN++;
	}
    else if (strcmp (pure_residue_nameP, "TYR") == 0)
	{
	first_atomSP->auxiliaryI = 3;
	total_phe_trp_tyrN++;
	}
    else first_atomSP->auxiliaryI = 0;
    }

if (total_phe_trp_tyrN <= 1) return -6;

/* Scan a cone about the k1_unit_vectorS. */
/* The  spherical  polar  system is used. */

/* Initialize the parameters which define a total number of directions: */
max_thetaI = 30;
max_phiI   = 72;

/* Prepare the angular increments. The maximal allowed */
/* deviation from  the initial vector is  45  degrees. */
theta_step = 0.7853982 / (double) max_thetaI;
phi_step   = 6.2831853 / (double) max_phiI;

/* Reset  the highest number of  PHE,  TRP */
/* and TYR residues in two selected rings: */
best_two_rings_count = 0;

/* Theta scan: */
for (thetaI = 0; thetaI <= max_thetaI; thetaI++)
    {
    /* Prepare the angle theta: */
    theta = (double) thetaI * theta_step;

    /* Phi scan: */
    for (phiI = 0; phiI <= max_phiI; phiI++)
	{
	/* Avoid to check too many directions close */
	/* to the pole. If the theta angle is equal */
	/* to zero, one phi value should be enough: */
	if (thetaI == 0)
	    {
	    if (phiI != 0) continue;
	    }

	/* If theta is too small, skip every second direction: */
	if (thetaI <= 2)
	    {
	    n = phiI / 2;
	    if (n * 2 != phiI) continue;
	    }

	/* Prepare the angle phi: */
	phi = (double) phiI * phi_step;

	/* The unit vector which defines the axis direction: */
	cos_theta = cos (theta);
	sin_theta = sin (theta);
	cos_phi   = cos (phi);
	sin_phi   = sin (phi);
	axis_unit_vectorS.x = sin_theta * cos_phi * i1_unit_vectorS.x +
			      sin_theta * sin_phi * j1_unit_vectorS.x +
			      cos_theta *           k1_unit_vectorS.x;
	axis_unit_vectorS.y = sin_theta * cos_phi * i1_unit_vectorS.y +
			      sin_theta * sin_phi * j1_unit_vectorS.y +
			      cos_theta *           k1_unit_vectorS.y;
	axis_unit_vectorS.z = sin_theta * cos_phi * i1_unit_vectorS.z +
			      sin_theta * sin_phi * j1_unit_vectorS.z +
			      cos_theta *           k1_unit_vectorS.z;

	/* Reset the array which counts  PHE,  TRP */
	/* and TYR residues projected to the axis: */
	for (simpleI = 0; simpleI < simple_cells_along_axisN; simpleI++)
	    {
	    phe_trp_tyrNA[simpleI] = 0;
	    }

	/* Scan residues: */
	for (residueI = 0; residueI < residuesN; residueI++)
	    {
	    /* Pointer to the current residue: */
	    residueSP = mol_complexSP->residueSP + residueI;

	    /* Pointer to the first atom of this residue: */
	    first_atomSP = mol_complexSP->atomSP + residueSP->residue_startI;

	    /* Check the auxiliaryI of the first atom; if */
	    /* it is equal to zero,  ignore this residue. */
	    if (first_atomSP->auxiliaryI == 0) continue;

	    /* Try to extract the CA and CB coordinates: */
	    n = ExtractCACB_ (&CA_vectorS, &CB_vectorS,
			      mol_complexSP->atomSP,
			      residueSP->residue_startI,
			      residueSP->residue_endI);
	    if (n < 2) continue;

	    /* Prepare the vector from CA to CB: */
	    CA_CB_vectorS.x = CB_vectorS.x - CA_vectorS.x;
	    CA_CB_vectorS.y = CB_vectorS.y - CA_vectorS.y;
	    CA_CB_vectorS.z = CB_vectorS.z - CA_vectorS.z;

	    /* The vector from the CA geometric center */
	    /* to the CB atom of  the current residue: */
	    CA_center_CB_vectorS.x = CB_vectorS.x - CA_center_vectorS.x;
	    CA_center_CB_vectorS.y = CB_vectorS.y - CA_center_vectorS.y;
	    CA_center_CB_vectorS.z = CB_vectorS.z - CA_center_vectorS.z;

	    /* The scalar product  between this */
	    /* vector and the axis unit vector: */
	    scalar_product = ScalarProduct_ (&CA_center_CB_vectorS,
					     &axis_unit_vectorS);

	    /* The  parallel component of  the vector  from the */
	    /* CA center to the CB atom of the current residue: */
	    parallel_vectorS.x = scalar_product * axis_unit_vectorS.x;
	    parallel_vectorS.y = scalar_product * axis_unit_vectorS.y;
	    parallel_vectorS.z = scalar_product * axis_unit_vectorS.z;

	    /* The perpendicular component: */
	    perpendicular_vectorS.x = CA_center_CB_vectorS.x -
				      parallel_vectorS.x;
	    perpendicular_vectorS.y = CA_center_CB_vectorS.y -
				      parallel_vectorS.y;
	    perpendicular_vectorS.z = CA_center_CB_vectorS.z -
				      parallel_vectorS.z;

	    /* Calculate and check the distance  from the axis. */
	    /* Skip residues whose  CB  atoms are  too close to */
	    /* the axis, they do not belong to the barrel wall. */
	    /* A limit of eleven angstroms might be reasonable. */
	    distance = AbsoluteValue_ (&perpendicular_vectorS);
	    if (distance <= 11.0) continue;

	    /* Check the scalar product between the CA-CB vector and */
	    /* the perpendicular component of the vector from the CA */
	    /* center to the  CB atom of the current residue.  If it */
	    /* is negative the CA-CB vector points towards the axis. */
	    scalar_product = ScalarProduct_ (&CA_CB_vectorS,
					     &perpendicular_vectorS);
	    if (scalar_product < 0.0) continue;

	    /* If this point is reached,  the CA-CB bond */
	    /* points outwards with respect to the axis. */

	    /* Extract the coordinates of the CG atom: */
	    n = ExtractCG_ (&CG_vectorS,
			    mol_complexSP->atomSP,
			    residueSP->residue_startI,
			    residueSP->residue_endI);
	    if (n < 1) continue;

	    /* The  CG atom will be projected over a number of cells. */
	    /* The index of the central cell  will be calculated from */
	    /* parallel part of the vector from CA center to CG atom. */

	    /* The vector from the CA geometric center */
	    /* to the CG atom of  the current residue: */
	    CA_center_CG_vectorS.x = CG_vectorS.x - CA_center_vectorS.x;
	    CA_center_CG_vectorS.y = CG_vectorS.y - CA_center_vectorS.y;
	    CA_center_CG_vectorS.z = CG_vectorS.z - CA_center_vectorS.z;

	    /* The scalar product  between this */
	    /* vector and the axis unit vector: */
	    scalar_product = ScalarProduct_ (&CA_center_CG_vectorS,
					     &axis_unit_vectorS);

	    /* The  parallel component of  the vector  from the */
	    /* CA center to the CG atom of the current residue: */
	    parallel_vectorS.x = scalar_product * axis_unit_vectorS.x;
	    parallel_vectorS.y = scalar_product * axis_unit_vectorS.y;
	    parallel_vectorS.z = scalar_product * axis_unit_vectorS.z;

	    /* The index of central  simple cell  is calculated */
	    /* from the projection of the parallel component of */
	    /* CA_center_CB_vectorS  to the  axis_unit_vectorS: */
	    scalar_product = ScalarProduct_ (&parallel_vectorS,
					     &axis_unit_vectorS);
	    central_simpleI = (int) (recip_simple_width * scalar_product) +
			      simple_index_offset;

	    /* The CG atom will be projected over a number of cells. */
	    /* The  side  chains  which  form  the ring of  aromatic */
	    /* residues  are not  perfectly  aligned  and  for  that */
	    /* reason the  CG atom will be smeared over  more cells. */
	    /* This  should improve  recognition of  aromatic rings. */
	    start_simpleI = central_simpleI - half_simple_window_width;
	    end_simpleI   = start_simpleI + simple_window_width;

	    /* Smear CG projection over a number of cells: */
	    for (simpleI = start_simpleI; simpleI <= end_simpleI; simpleI++)
		{
		/* Do not allow array overflow: */
		if (simpleI < 0) continue;
		if (simpleI >= simple_cells_along_axisN) break;

		/* Update the number of PHE, TRP and TYR residues which */
		/* were projected to the cell defined by this  simpleI: */
		phe_trp_tyrNA[simpleI]++;
		}

	    /* End of residueI loop: */
	    }

	/* For this orientation, check the total number of PHE, TRP and TYR */
	/* residues in two rings separated by the given membrane thickness. */
	for (simpleI = 0; simpleI < max_simpleI; simpleI++)
	    {
	    two_rings_count = phe_trp_tyrNA[simpleI] +
			      phe_trp_tyrNA[simpleI + simple_thickness];
	    if (two_rings_count > best_two_rings_count)
		{
		best_two_rings_count = two_rings_count;
		best_unit_vectorS.x = axis_unit_vectorS.x;
		best_unit_vectorS.y = axis_unit_vectorS.y;
		best_unit_vectorS.z = axis_unit_vectorS.z;
		}
	    }

	/* End of phiI loop: */
	}

    /* End of thetaI loop: */
    }

/*------store the normal vector (refined):-----------------------------------*/

/* Store (copy) the components of the normal vector. */

/* The normal vector of the first plane points to the same direction: */
mol_complexSP->membraneS.plane1S.normal_x[0] =  best_unit_vectorS.x;
mol_complexSP->membraneS.plane1S.normal_y    =  best_unit_vectorS.y;
mol_complexSP->membraneS.plane1S.normal_z[0] =  best_unit_vectorS.z;

/* The normal vector of the second plane points to the opposite direction: */
mol_complexSP->membraneS.plane2S.normal_x[0] = -best_unit_vectorS.x;
mol_complexSP->membraneS.plane2S.normal_y    = -best_unit_vectorS.y;
mol_complexSP->membraneS.plane2S.normal_z[0] = -best_unit_vectorS.z;

/*------prepare stereo data (refined):---------------------------------------*/

/* Prepare the stereo data for both normal vectors. */

/* Prepare the sine and cosine of the stereo angle: */
sin_angle = sin (configSP->stereo_angle);
cos_angle = cos (configSP->stereo_angle);

/* Calculate the stereo data: */
mol_complexSP->membraneS.plane1S.normal_x[1] =
		 mol_complexSP->membraneS.plane1S.normal_x[0] * cos_angle +
		 mol_complexSP->membraneS.plane1S.normal_z[0] * sin_angle;
mol_complexSP->membraneS.plane1S.normal_z[1] =
		-mol_complexSP->membraneS.plane1S.normal_x[0] * sin_angle +
		 mol_complexSP->membraneS.plane1S.normal_z[0] * cos_angle;
mol_complexSP->membraneS.plane2S.normal_x[1] =
		 mol_complexSP->membraneS.plane2S.normal_x[0] * cos_angle +
		 mol_complexSP->membraneS.plane2S.normal_z[0] * sin_angle;
mol_complexSP->membraneS.plane2S.normal_z[1] =
		-mol_complexSP->membraneS.plane2S.normal_x[0] * sin_angle +
		 mol_complexSP->membraneS.plane2S.normal_z[0] * cos_angle;

/*------update auxiliaryI:---------------------------------------------------*/

/* In this section,  the auxiliaryI for  PHE,  TRP and TYR will */
/* be set to zero for side  chains which are pointing  inwards. */
/* This should improve performance: the axis orientation during */
/* the membrane center  refinement  is not  subject  to change. */
for (residueI = 0; residueI < residuesN; residueI++)
    {
    /* Pointer to the current residue: */
    residueSP = mol_complexSP->residueSP + residueI;

    /* Pointer to the first atom of this residue: */
    first_atomSP = mol_complexSP->atomSP + residueSP->residue_startI;

    /* Check the auxiliaryI of the first atom; if */
    /* it is equal to zero,  ignore this residue. */
    if (first_atomSP->auxiliaryI == 0) continue;

    /* Try to extract the CA and CB coordinates: */
    n = ExtractCACB_ (&CA_vectorS, &CB_vectorS,
		      mol_complexSP->atomSP,
		      residueSP->residue_startI,
		      residueSP->residue_endI);
    if (n < 2) continue;

    /* Prepare the vector from CA to CB: */
    CA_CB_vectorS.x = CB_vectorS.x - CA_vectorS.x;
    CA_CB_vectorS.y = CB_vectorS.y - CA_vectorS.y;
    CA_CB_vectorS.z = CB_vectorS.z - CA_vectorS.z;

    /* The vector from the CA geometric center */
    /* to the CB atom of  the current residue: */
    CA_center_CB_vectorS.x = CB_vectorS.x - CA_center_vectorS.x;
    CA_center_CB_vectorS.y = CB_vectorS.y - CA_center_vectorS.y;
    CA_center_CB_vectorS.z = CB_vectorS.z - CA_center_vectorS.z;

    /* The scalar product  between this */
    /* vector and the axis unit vector: */
    scalar_product = ScalarProduct_ (&CA_center_CB_vectorS,
				     &axis_unit_vectorS);

    /* The  parallel component of  the vector  from the */
    /* CA center to the CB atom of the current residue: */
    parallel_vectorS.x = scalar_product * axis_unit_vectorS.x;
    parallel_vectorS.y = scalar_product * axis_unit_vectorS.y;
    parallel_vectorS.z = scalar_product * axis_unit_vectorS.z;

    /* The perpendicular component: */
    perpendicular_vectorS.x = CA_center_CB_vectorS.x - parallel_vectorS.x;
    perpendicular_vectorS.y = CA_center_CB_vectorS.y - parallel_vectorS.y;
    perpendicular_vectorS.z = CA_center_CB_vectorS.z - parallel_vectorS.z;

    /* Calculate and check the distance  from the axis. */
    /* Skip residues whose  CB  atoms are  too close to */
    /* the axis, they do not belong to the barrel wall. */
    /* A limit of eleven angstroms might be reasonable. */
    distance = AbsoluteValue_ (&perpendicular_vectorS);
    if (distance <= 11.0) continue;

    /* Check the scalar product between the CA-CB vector and */
    /* the perpendicular component of the vector from the CA */
    /* center to the  CB atom of the current residue.  If it */
    /* is negative the CA-CB vector points towards the axis. */
    scalar_product = ScalarProduct_ (&CA_CB_vectorS, &perpendicular_vectorS);
    if (scalar_product < 0.0) first_atomSP->auxiliaryI = 0;
    }

/*------(7) refine the membrane center:--------------------------------------*/

/* Initialize the best shift: */
best_shift = 0.0;

/* Reset  the highest number of  PHE,  TRP */
/* and TYR residues in two selected rings: */
best_two_rings_count = 0;

/* Scan a number of positions along the refined axis.  The old CA */
/* center should fall in the middle of the array of simple cells. */
start_pointI = -simple_scan_half_width;
end_pointI   =  simple_scan_half_width;
for (pointI = start_pointI; pointI <= end_pointI; pointI++)
    {
    /* Shift, in angstroms: */
    shift = (double) pointI * simple_cell_width;

    /* Presumed membrane center position: */
    center_vectorS.x = CA_center_vectorS.x + shift * best_unit_vectorS.x;
    center_vectorS.y = CA_center_vectorS.y + shift * best_unit_vectorS.y;
    center_vectorS.z = CA_center_vectorS.z + shift * best_unit_vectorS.z;

    /* Reset the array which counts  PHE,  TRP */
    /* and TYR residues projected to the axis: */
    for (simpleI = 0; simpleI < simple_cells_along_axisN; simpleI++)
	{
	phe_trp_tyrNA[simpleI] = 0;
	}

    /* Scan the macromolecular complex, residue by residue: */
    for (residueI = 0; residueI < residuesN; residueI++)
	{
	/* Pointer to the current residue: */
	residueSP = mol_complexSP->residueSP + residueI;

	/* Pointer to the first atom of this residue: */
	first_atomSP = mol_complexSP->atomSP + residueSP->residue_startI;

	/* Check the  auxiliaryI of the first atom;  if */ 
	/* it is equal to zero, ignore this residue. It */
	/* is either of  wrong type or  points inwards. */
	if (first_atomSP->auxiliaryI == 0) continue;

	/* Extract the coordinates of the CG atom: */
	n = ExtractCG_ (&CG_vectorS, mol_complexSP->atomSP,
			residueSP->residue_startI, residueSP->residue_endI);
	if (n < 1) continue;

	/* The parallel part of the vector from the presumed membrane */
	/* center to the CG atom will be used to calculate  the index */
	/* of the cell  to which  the  CG atom  should  be projected. */

	/* The vector from the current center to */
	/* the  CG atom of  the current residue: */
	radius_vectorS.x = CG_vectorS.x - center_vectorS.x;
	radius_vectorS.y = CG_vectorS.y - center_vectorS.y;
	radius_vectorS.z = CG_vectorS.z - center_vectorS.z;

	/* Projection of this vector to the axis: */
	projection = ScalarProduct_ (&radius_vectorS, &best_unit_vectorS);

	/* The index of the central simple cell: */
	central_simpleI = (int) (recip_simple_width * projection) +
			  simple_index_offset;

	/* The CG atom will be projected over a number of cells. */
	/* The  side  chains  which  form  the ring of  aromatic */
	/* residues  are not  perfectly  aligned  and  for  that */
	/* reason the  CG atom will be smeared over  more cells. */
	start_simpleI = central_simpleI - half_simple_window_width;
	end_simpleI   = start_simpleI + simple_window_width;

	/* Smear CG projection over a number of cells: */
	for (simpleI = start_simpleI; simpleI <= end_simpleI; simpleI++)
	    {
	    /* Do not allow array overflow: */
	    if (simpleI < 0) continue;
	    if (simpleI >= simple_cells_along_axisN) break;

	    /* Update the number of PHE, TRP and TYR residues which */
	    /* were projected to the cell defined by this  simpleI: */
	    phe_trp_tyrNA[simpleI]++;
	    }

	/* End of residueI loop: */
	}

    /* Check the count of PHE, TRP and TYR residues in two rings: */
    ring1I = simple_index_offset - half_simple_thickness;
    if (ring1I < 0) continue;
    if (ring1I >= simple_cells_along_axisN) continue;
    ring2I = ring1I + simple_thickness;
    if (ring2I < 0) continue;
    if (ring2I >= simple_cells_along_axisN) continue;
    two_rings_count = phe_trp_tyrNA[ring1I] + phe_trp_tyrNA[ring2I];
    if (two_rings_count > best_two_rings_count)
	{
	best_two_rings_count = two_rings_count;
	best_shift = shift;
	}

    /* End of pointI loop: */
    }

/*------update the membrane center:------------------------------------------*/

/* Update the membrane center: */

/* Update the membrane center: */
mol_complexSP->membraneS.center_x = CA_center_vectorS.x +
				    best_shift * best_unit_vectorS.x;
mol_complexSP->membraneS.center_y = CA_center_vectorS.y +
				    best_shift * best_unit_vectorS.y;
mol_complexSP->membraneS.center_z = CA_center_vectorS.z +
				    best_shift * best_unit_vectorS.z;

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

/* Return positive value on success: */
return 1;
}

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