C++*********************************************************************
C
C POLAR_CC.F
C              OPFILEC                               FEB 03 ARDEAN LEITH
C              DE-OMP PARRALEL LOOP FOR FFTW3 USE    MAR 08 ARDEAN LEITH
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 **********************************************************************

        SUBROUTINE POLAR_CC

        INCLUDE 'CMBLOCK.INC'
        INCLUDE 'CMLIMIT.INC' 
 
        CHARACTER(LEN=MAXNAM)   ::   FINPIC,FINPAT
        COMMON  /F_SPEC/  FINPAT,NLET,FINPIC

        REAL, ALLOCATABLE, DIMENSION(:,:) :: X,REF
        CHARACTER(LEN=1) ::   MODE,ASK,NULL
        DOUBLE PRECISION  T7(-3:3)

        DATA  INPIC/77/,INREF/76/,LUNO/75/

        NULL = CHAR(0)

C       ASK FOR DATA FILES
        MAXIM = 0
        CALL OPFILEC(0,.TRUE.,FINPIC,INPIC,'O',ITYPE,NSAM,NROW,
     &             NSLICE,MAXIM,'INPUT',.FALSE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

        MAXIM = 0
        CALL OPFILEC(0,.TRUE.,FINPIC,INREF,'O',ITYPE,NSAMP,NROWP,
     &              NSLICE,MAXIM,'REFERENCE',.FALSE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

        IF (NSAM.NE.NSAMP .OR. NROW.NE.NROWP)  THEN
           CALL ERRT(1,'POLAR_CC',IER)
           RETURN
        ENDIF

C       ASK FOR OUTPUT FILE
        MAXIM = 0
	ITYPE=1
	NROWP=1
	NSLICEP=1
        CALL OPFILEC(INPIC,.TRUE.,FINPIC,LUNO,'N',ITYPE,NSAM,NROWP,
     &             NSLICEP,MAXIM,'OUTUT',.FALSE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

        NA=1
        CALL  RDPRMC(ASK,NA,.TRUE.,'(F)ULL OR (H)ALF CIRCLE',NULL,IRT)
        IF (ASK.EQ.'H')  THEN
           MODE='H'
        ELSE
           MODE='F'
        ENDIF

        PI   = 4*DATAN(1.0D0)
        LSAM = NSAM+2-MOD(NSAM,2)
        ALLOCATE (X(LSAM,NROW),REF(LSAM,NROW), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'POLAR_CC, X,REF',IER)
           CLOSE(INPIC)
           CLOSE(INREF)
           RETURN
        ENDIF

        DO J=1,NROW
           CALL REDLIN(INPIC,X(1,J),NSAM,J)
        ENDDO
        CLOSE(INPIC)
        DO J=1,NROW
           CALL REDLIN(INREF,REF(1,J),NSAM,J)
        ENDDO
        CLOSE(INREF)

        INV = +1
C       parallel do private(j),shared(invt) REMOVED FOR FFTW3 al
        DO  J=1,NROW
           INVT = INV
           CALL FMRS_1(X(1,J),NSAM,INVT)
           CALL FMRS_1(REF(1,J),NSAM,INVT)
        ENDDO
        IF (INV .EQ. 0) THEN
           CALL ERRT(101,'*** ERROR IN POLAR_CC, INV = 0',NDUM)
           RETURN
        ENDIF

        CALL  MF1D(X,REF,LSAM/2,NROW)
        INV=-1
        CALL FMRS_1(X(1,1),NSAM,INV)
        QM=-1.0E20
        DO    I=1,NSAM
           IF (X(I,1) .GE. QM)  THEN
              QM   = X(I,1)
              JTOT = I
           ENDIF
        ENDDO
        DO    K=-3,3
           J     = MOD(JTOT+K+NSAM-1,NSAM)+1
           T7(K) = X(J,1)
        ENDDO
        CALL  PRB1D(T7,7,POS)
        RANG=REAL(JTOT)+POS
        IF (MODE.EQ.'F')  THEN
           RANG=(RANG-1.0)/NSAM*360.0
        ELSE
           RANG=(RANG-1.0)/NSAM*180.0
        ENDIF   
        WRITE(NOUT,2700)  RANG,QM
2700    FORMAT('    Angle = ',F10.4,'    Peak = ',G12.5)

        CALL REG_SET_NSEL(1,2,RANG,QM,0.0,0.0,0.0,IRTFLG)
C       Shift 1D CCF to have the origin at INT(NSAM/2)+1
	DO I=1,NSAM
	  REF(MOD(I+INT(NSAM/2)-1,NSAM)+1,1)=X(I,1)
	ENDDO
	CALL  WRTLIN(LUNO,REF,NSAM,1)
	CLOSE(LUNO)

        DEALLOCATE (X,REF)

        END


C       -------------------- MF1D -------------------------------------

        SUBROUTINE  MF1D(X,REF,L,NROW)

        COMPLEX  X(L,NROW),REF(L,NROW)

        X(:,1) = X(:,1)*CONJG(REF(:,1))
        IF (NROW .EQ. 1)  RETURN

        DO J=2,NROW
           X(:,1) = X(:,1) + X(:,J) * CONJG(REF(:,J))
        ENDDO
        END