C++*********************************************************************
C
C CROSRNG_E.F  MERGED CROSRMG_DS.F & CROSRNG_DS    AUG 04 ARDEAN LEITH
C              REWRITE FOR FFTW3                   MAR 08 ARDEAN LEITH
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2008  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_E(CIRC1,CIRC2,LCIRC,NRING,MAXRIN,NUMR,QMAX,POS_MAX,TT,NEG)
C
C PURPOSE: CROSS CORRELATION OF RADIAL RINGS FOR USE IN ROTATIONAL
C          ALIGNMENT. CHECKS ONLY UN-MIRRORED POSITION
C
C PARAMETERS:
C    CIRC1  - FT OF RINGS MULTIPLIED BY WEIGHTS           (SENT)
C    CIRC2  - FT OF RINGS MULTIPLIED BY WEIGHTS           (SENT)
C    LCIRC  - SIZE OF CIRCS ARRAYS                        (SENT)
C    NRING  - NUMBER OF RINGS                             (SENT)
C    MAXRIN - LONGEST RING                                (SENT)
C    NUMR   - RING LOCATION POINTERS                      (SENT)
C    QMAX   - CC MAX                                      (RETURNED)
C    POS_MAX - POSITION OF CC MAX                          (RETURNED)
C    TT     - USED FOR SGI FFT (UNUSED NOW)               (SENT)
C    NEG    - FLAG FOR CONJUGATE (MIRROR) OF 1'ST RING    (SENT)
C
C     U = 11*21+ 12*22  (e.g.)
C
C     ABNORMAL LOCATION   <MAXRIN       NVAL=MAXRIN   (U = USUAL)
C     E.g. NVAL:            256             512
C        1:                 11x21          11x21         (same)
C        2:                   0            12x22 
C        NVAL-1               U              U           (same)
C        NVAL                 U              U (maxrin)  (same)
C        NVAL+1             12x22            0 (unused)
C        MAXRIN               0 (unused)     U               
C        MAXRIN + 1:          0              0           (same) (unused) 
C        MAXRIN + 2:          0              0           (same) (unused)    
C        summation          0..257         0..514
C
C  NOTES:     COMPLEX CONJUGATE OF COMPLEX NUMBER = a - bi
C             USUAL COMPLEX MULTIPLICATION: (ac-bd),(ad+bc)
C             OUR   COMPLEX MULTIPLICATION: (ac+bd),(-ad+bc) 
C             (If loaded by fmrs can use usual method)
C              SPIDER_SIGN flag when data originally transformed))
C
C    Typical radii:  32, 64, 128, 256, 512  (only)
c  ip:            -9    lcirc:          23168
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*********************************************************************



C--************************  CROSRNG_EP_NEW ******************************

        SUBROUTINE CROSRNG_EP_NEW(CIRC1,CIRC2,LCIRC, NRING,MAXRIN,NUMR,
     &                           QMAX,POS_MAX, TT,NEG, FFTW3PLAN)

C       USES NUMR TABLE FOR MAPPING INTO Q ARRAY 
C       USES SIMPLIFIED LOGIC FOR BOUNDARY VALUES, FLOATING PT. ARITH.

        INTEGER, INTENT(IN)           :: NUMR(3,NRING)
        REAL, INTENT(IN)              :: CIRC1(LCIRC), CIRC2(LCIRC)
        DOUBLE PRECISION, INTENT(OUT) :: QMAX
        REAL, INTENT(OUT)             :: POS_MAX
	DOUBLE PRECISION, INTENT(IN)  :: TT(*)
        LOGICAL, INTENT(IN)           :: NEG

        LOGICAL                       :: NEG_USED

C       AUTOMATIC ARRAYS
        REAL                          :: Q(MAXRIN+2)

        NEG_USED = NEG

C       ZERO WHOLE Q ARRAY,  STRAIGHT  = CIRC1 * CONJG(CIRC2)
        Q = 0.0
     
        DO I=1,NRING
           IGO   = NUMR(2,I)
           IGOM1 = IGO - 1
           NVAL  = NUMR(3,I)       
           J1    = 1
           
#ifndef SP_LIBFFTW3
C          NATIVE OR SGI FFT
           Q(1) = Q(1) + (CIRC1(IGO)) * CIRC2(IGO)
           IF (NVAL .NE. MAXRIN) THEN
              Q(NVAL+1) = Q(NVAL+1) + (CIRC1(IGO+1)) * CIRC2(IGO+1)
           ELSE
              Q(2)      = Q(2)      + (CIRC1(IGO+1)) * CIRC2(IGO+1)            
           ENDIF
           J1 = 3
#endif

           IF (NEG_USED) THEN
C             FIRST RING SET IS CONJUGATED (MIRRORED)
              DO J=J1,NVAL,2
                 JC     = J + IGOM1

                 Q(J)   = Q(J)   + (CIRC1(JC))   * CIRC2(JC)   -
     &                             (CIRC1(JC+1)) * CIRC2(JC+1)
                 Q(J+1) = Q(J+1) - (CIRC1(JC))   * CIRC2(JC+1) -
     &                             (CIRC1(JC+1)) * CIRC2(JC)
              ENDDO
           ELSE
C             FIRST RING SET IS NON-CONJUGATED (NOT MIRRORED)
              DO J=J1,NVAL,2
                 JC     = J + IGOM1

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

        CALL CROSRNG_COM_NEW(Q,LCIRC,MAXRIN,QMAX,POS_MAX,TT,FFTW3PLAN)
 
        END


#ifdef NEVER
           J  = NVAL  + 1
           JC = IGOM1 + NVAL 
           IF (NVAL .NE. MAXRIN) THEN
              Q(J)   = Q(J)   + ((CIRC1(JC))   * CIRC2(JC)   + IPOS *
     &                           (CIRC1(JC+1)) * CIRC2(JC+1) / 2.0)
              Q(J+1) = Q(J+1) - ((CIRC1(JC))   * CIRC2(JC+1) + IPOS *
     &                           (CIRC1(JC+1)) * CIRC2(JC) / 2.0)
           ELSE
              Q(J)   = Q(J)   + (CIRC1(JC))   * CIRC2(JC)   + IPOS *
     &                          (CIRC1(JC+1)) * CIRC2(JC+1)
              Q(J+1) = Q(J+1) - (CIRC1(JC))   * CIRC2(JC+1) + IPOS *
     &                          (CIRC1(JC+1)) * CIRC2(JC)
           ENDIF
#endif


C--************************  CROSRNG_COM_NEW ******************************

        SUBROUTINE CROSRNG_COM_NEW(QS,LCIRC,MAXRIN,QMAX,POS_MAX,
     &                             TT,FFTW3PLAN)

C       COMMON CODE FOR REVERSE TRANSFORM AND MAX LOCATION DETERMINATION

        REAL                          :: QS(MAXRIN + 2)
        INTEGER, INTENT(IN)           :: LCIRC,MAXRIN
        DOUBLE PRECISION, INTENT(OUT) :: QMAX
        REAL, INTENT(OUT)             :: POS_MAX
        INTEGER*8, INTENT(IN)         :: FFTW3PLAN  ! STRUCTURE POINTER 
        DOUBLE PRECISION              :: TT(*)      ! FOR SGI LIBFFT

C       AUTOMATIC ARRAYS
        DOUBLE PRECISION              :: T7(-3:3)
        INTEGER                       :: MAXL_ARRRAY(1)

        LOGICAL, PARAMETER            :: SPIDER_SIGN = .FALSE.

        INV = -1

#ifdef SP_LIBFFTW3        
C       REVERSE FOURIER TRANSFORM ON QS   
        CALL FMRS(QS,MAXRIN, 1,1,FFTW3PLAN, SPIDER_SIGN,.FALSE., 
     &            INV,IRTFLG)
#else
#if defined (SP_LIBFFT)
C       SGI FFT USED
        LDA         = 1
        QS(MAXRIN+1) = QS(2)
        QS(2)        = 0.0
        QS(MAXRIN+2) = 0.0

C       REVERSE FOURIER TRANSFORM ON QS   
        CALL ZDFFT1DU(INV,MAXRIN,QS,LDA,TT)
#else
C       NATIVE SPIDER FFT USED
        IP = -LOG2(MAXRIN)

C       REVERSE FOURIER TRANSFORM ON QS   
        CALL FFTR_D(QS,IP)
#endif
#endif
 
C       FIND MAXIMUM AND ITS LOCATION INSIDE QS
        MAXL_ARRRAY = MAXLOC(QS)      ! RETURNS ARRAY OF LENGTH: 1
        MAXL        = MAXL_ARRRAY(1)  ! LOCATION OF MAXIMUM
        QMAX        = QS(MAXL)        ! MAXIMUM VALUE

C       INTERPOLATION OVER 3x3 NEIGHBORHOOD
        DO K=-3,3
           J     = MOD(MAXL+K+MAXRIN-1,MAXRIN) + 1
           T7(K) = QS(J)
        ENDDO

        CALL PRB1D(T7,7,POS)
        POS_MAX = FLOAT(MAXL) + POS    ! SUB-PIXEL LOCATION

#ifdef SP_LIBFFTW3        
        QMAX = 1.00048 * QMAX / (MAXRIN)
#endif

#ifdef DEBUGNEVER
        write(6,901)maxl,qmax,f_maxl,qs(15)
 901    format(' crosrng_com_s; maxl: ',i6,' qmax: ',f16.3,
     &         ' f_maxl: ',f10.3,' qs(15):  ',f16.3)
#endif

        END

C       --------------- CROSRNG_E ------------------------------------

        SUBROUTINE CROSRNG_E(CIRC1,CIRC2,LCIRC, NRING,MAXRIN,NUMR,
     &                           QMAX,POS_MAX,TT,NEG)

C       USES NUMR TABLE FOR MAPPING INTO Q ARRAY 
C       USES SIMPLIFIED LOGIC FOR BOUNDARY VALUES, FLOATING PT. ARITH.

        INTEGER, INTENT(IN) :: NUMR(3,NRING)
        REAL, INTENT(IN)    :: CIRC1(LCIRC), CIRC2(LCIRC)
        DOUBLE PRECISION    :: QMAX
        REAL, INTENT(OUT)   :: POS_MAX
	DOUBLE PRECISION    :: TT(*)
        LOGICAL, INTENT(IN) :: NEG

C       AUTOMATIC ARRAYS
        DOUBLE PRECISION    :: Q(MAXRIN+2)

C       ZERO WHOLE Q ARRAY
        Q = 0.0
     
        DO I=1,NRING
           IGO  = NUMR(2,I)
           NVAL = NUMR(3,I)

           Q(1) = Q(1) + (CIRC1(IGO)) * CIRC2(IGO)

           IF (NVAL .NE. MAXRIN) THEN
              Q(NVAL+1) = Q(NVAL+1) + (CIRC1(IGO+1)) * CIRC2(IGO+1)
           ELSE
              Q(2)      = Q(2)      + (CIRC1(IGO+1)) * CIRC2(IGO+1)            
           ENDIF

           IGOM1 = IGO - 1

           IF (NEG) THEN
C             FIRST RING SET IS CONJUGATED (MIRRORED)
              DO J=3,NVAL,2
                 JC     = J + IGOM1

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

                 Q(J+1) = Q(J+1) - (CIRC1(JC))   * CIRC2(JC+1) -
     &                             (CIRC1(JC+1)) * CIRC2(JC)
              ENDDO
           ELSE
C             FIRST RING SET IS NON-CONJUGATED (NOT MIRRORED)
              DO J=3,NVAL,2
                 JC     = J + IGOM1

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

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

        CALL CROSRNG_COM(Q,LCIRC,MAXRIN,QMAX,POS_MAX,TT)
 
        END

C--************************  CROSRNG_COM ******************************

        SUBROUTINE CROSRNG_COM(QS,LCIRC,MAXRIN,QMAX,POS_MAX,TT)

C       COMMON CODE FOR REVERSE TRANSFORM AND MAX LOCATION DETERMINATION

        DOUBLE PRECISION              :: QS(MAXRIN + 2)
        INTEGER, INTENT(IN)           :: LCIRC,MAXRIN
        DOUBLE PRECISION, INTENT(OUT) :: QMAX
        REAL, INTENT(OUT)             :: POS_MAX
        DOUBLE PRECISION, INTENT(IN)  :: TT(*)

        DOUBLE PRECISION              :: T7(-3:3)
        INTEGER                       :: MAXL_ARRRAY(1)
        
C       REVERSE FOURIER TRANSFORM ON QS   

        IP = MAXRIN

#ifdef SP_LIBFFT
        INV         = -1
        LDA         = 1
        QS(MAXRIN+1) = QS(2)
        QS(2)        = 0.0
        QS(MAXRIN+2) = 0.0
        CALL ZDFFT1DU(INV,IP,QS,LDA,TT)
#else
        IP = -LOG2(IP)
        CALL FFTR_D(QS,IP)
#endif

C       FIND MAXIMUM AND ITS LOCATION INSIDE QS
        MAXL_ARRRAY = MAXLOC(QS)      ! RETURNS ARRAY OF LENGTH: 1
        MAXL        = MAXL_ARRRAY(1)  ! LOCATION OF MAXIMUM
        QMAX        = QS(MAXL)        ! MAXIMUM VALUE

#ifdef SP_LIBFFT
        QMAX = QMAX/MAXRIN
#endif

C       INTERPOLATION OVER 3x3 NEIGHBORHOOD
        DO K=-3,3
           J     = MOD(MAXL+K+MAXRIN-1,MAXRIN) + 1
           T7(K) = QS(J)
        ENDDO

        CALL PRB1D(T7,7,POS)
        POS_MAX = FLOAT(MAXL) + POS    ! SUB-PIXEL LOCATION

#ifdef DEBUGNEVER
        write(6,901)maxl,qmax,f_maxl,qs(15)
 901    format(' crosrng_com_s; maxl: ',i6,' qmax: ',f16.3,
     &         ' f_maxl: ',f10.3,' qs(15):  ',f16.3)
#endif

        END