C++*********************************************************************
C
C CROSRNG_DS.F 
C
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
C=*                                                                    *
C=* This program is free software; you can redistribute it and/or      *
C=* modify it under the terms of the GNU General Public License as     *
C=* published by the Free Software Foundation; either version 2 of the *
C=* License, or (at your option) any later version.                    *
C=*                                                                    *
C=* This program is distributed in the hope that it will be useful,    *
C=* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C
C CROSRNG_DS(CIRC1,CIRC2,LCIRC,NRING,MAXRIN,NUMR,QN,TOT,TT)
C
C  INPUT - FOURIER TRANSFORMS OF RINGS!!!
C  CIRC1 ALREADY MULTIPLIED BY WEIGHTS!!
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*********************************************************************

	SUBROUTINE CROSRNG_DS(CIRC1,CIRC2,LCIRC,NRING,MAXRIN,NUMR,
     &                        QN,TOT,TT)

	INTEGER           NUMR(3,NRING)
	DIMENSION         CIRC1(LCIRC),CIRC2(LCIRC)
	DOUBLE PRECISION  TT(*),QN

C       AUTOMATIC ARRAYS
	DIMENSION         T(MAXRIN+2)
	DOUBLE PRECISION  Q(MAXRIN+2)
	DOUBLE PRECISION  T7(-3:3)

        DATA  ICALL/0/

	IP = MAXRIN
#ifndef SP_LIBFFT
	IP=-LOG2(IP)
#endif


C       ZERO Q ARRAY
	Q = 0.0D0
     
	T(MAXRIN+1) = 0.0
	T(MAXRIN+2) = 0.0

 	DO I=1,NRING
	   NUMR3I = NUMR(3,I)
	   NUMR2I = NUMR(2,I)

	   T(1) = (CIRC1(NUMR2I)) * CIRC2(NUMR2I)

	   IF (NUMR3I .NE. MAXRIN)  THEN
	      T(NUMR3I+1) = (CIRC1(NUMR2I+1))*CIRC2(NUMR2I+1)
	      T(2)        = 0.0

	      DO J=3,NUMR3I,2
	         JC     = J+NUMR2I-1

	         T(J)   =  (CIRC1(JC))*CIRC2(JC) +
     &                     (CIRC1(JC+1))*CIRC2(JC+1)

	         T(J+1) = -(CIRC1(JC))*CIRC2(JC+1) +
     &                     (CIRC1(JC+1))*CIRC2(JC)
	      ENDDO

	      Q(1:NUMR3I+1) = Q(1:NUMR3I+1) + T(1:NUMR3I+1)

	   ELSE
	      T(2) = CIRC1(NUMR2I+1) * CIRC2(NUMR2I+1)

	      DO J=3,MAXRIN,2
	         JC     =  J+NUMR2I-1

	         T(J)   =  (CIRC1(JC))*CIRC2(JC) +
     &                     (CIRC1(JC+1))*CIRC2(JC+1)

	         T(J+1) = -(CIRC1(JC))*CIRC2(JC+1) +
     &                     (CIRC1(JC+1))*CIRC2(JC)
	      ENDDO
	      Q = Q + T
           ENDIF
	ENDDO
 
#ifdef SP_LIBFFT
	INV         = -1
	LDA         = 1
	Q(MAXRIN+1) = Q(2)
	Q(2)        = 0.0
	Q(MAXRIN+2) = 0.0
	CALL  ZDFFT1DU(INV,IP,Q,LDA,TT)

C       SKIP THE NORMALIZATION, DIVIDE THE MAXIMUM INSTEAD.
C	CALL DSCAL1D(IP,1.0D0/DBLE(FLOAT(IP)),Q,LDA)
#else
	CALL FFTR_D(Q,IP)
#endif

	QN = -1.0D20
	DO J=1,MAXRIN
	   IF (Q(J) .GE. QN)  THEN
	      QN   = Q(J)
	      JTOT = J
	   ENDIF
	ENDDO
#ifdef SP_LIBFFT
	QN=QN/MAXRIN
#endif

	DO K=-3,3
           J     = MOD(JTOT+K+MAXRIN-1,MAXRIN) + 1
	   T7(K) = Q(J)
	ENDDO

	CALL  PRB1D(T7,7,POS)

	TOT = FLOAT(JTOT) + POS

	END