
C ++********************************************************************
C                                                                      *
C RFACTSD2                                                             *
C                      REGISTER OUTPUT ADDED    JAN 2005 ARDEAN LEITH  *
C                      VERBOSE                  FEB 2006 ARDEAN LEITH                                                            *
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 RFACTSD2(PR,AMP,CSUM1,LR,CSUM,CSUM2,AVSUM, NSCALE,INC,WI,FACT,NOUT)
C
C IMAGE_PROCESSING_ROUTINE
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************

         SUBROUTINE RFACTSD2(PR,AMP,CSUM1,LR,CSUM,CSUM2,AVSUM,
     &                        NSCALE,INC,WI,FACT,NOUT)

         PARAMETER  (NLIST=6)
         DIMENSION  DLIST(NLIST)
         DIMENSION  PR(NSCALE,INC),AMP(NSCALE,INC),CSUM1(INC),
     &              CSUM(INC),CSUM2(INC),AVSUM(NSCALE,INC)
         INTEGER    LR(INC)
         LOGICAL  :: IFOUNDIT

C        VERBOSE    USER SET CONTROL VARIABLE FOR PRINTING OF INFO
         LOGICAL  :: SILENT,VERBOSE,DELAY_FREE,USE_SPIRE
         COMMON   /IPRTT/IDUM245,NTRACE,NALPH,VERBOSE,USE_SPIRE,SILENT

         DATA  NDOC/88/

#ifdef USE_MPI
        INCLUDE 'mpif.h'

        ICOMM  = MPI_COMM_WORLD
        MPIERR = 0
        CALL MPI_COMM_RANK(ICOMM, MYPID, MPIERR)
#else
        MYPID = -1
#endif

         CALL REG_GET_USED(NSEL_USED)
         IF (NSEL_USED .GT. 0) THEN
C           OUTPUT TO SPIDER'S REGISTERS NEEDED LATER
            XPREV  = 0
            FSCCUT = 0.5       ! FSC CUTOFF (J.f.'s)
            DLIST  = HUGE(FSCLAST)
         ENDIF

         DO  L=1,INC
            JDUM = LR(L)
            IF (JDUM .NE. 0) THEN
               DLIST(1) = L

               SPFLAST  = DLIST(2)
               DLIST(2) = FLOAT(L-1)/FLOAT(INC-1)*0.5
               DLIST(5) = AMIN1(1.0,FACT/SQRT(FLOAT(JDUM)))
               DLIST(6) = JDUM

               RFMIN    = -HUGE(RFMIN)
               NSCM     = 1
               IFOUNDIT = .FALSE.
               DO  NSC=1,NSCALE
		 IF (AMP(NSC,L).GT.TINY(RFMIN))  THEN
                    RFM = AVSUM(NSC,L) / AMAX1(1.0,AMP(NSC,L))
                    IF (RFM .LT. RFMIN) THEN
                       NSCM     = NSC
                       RFMIN    = RFM
                       IFOUNDIT = .TRUE.
                    ENDIF
		 ENDIF
               ENDDO


C              NSCM IS THE NUMBER OF THE ELEMENT IN EACH ARRAY WITH THE
C              CORRECT SCALING. SCALE IS THE CORRECT SCALING.

               BK1 = AMP(NSCM,L)
               BK2 = PR(NSCM,L)
               IF (BK1 .GT. TINY(BK3)) THEN
		  DLIST(3) = SQRT(BK2/BK1)
               ELSE
		  DLIST(3) = 0.0
               ENDIF
               BK3     = CSUM2(L)
               BK4     = CSUM1(L)

               FSCLAST = DLIST(4)

               IF (BK3.GT.TINY(BK3) .AND. BK4.GT.TINY(BK3)) THEN
		  DLIST(4) = CSUM(L) / SQRT(BK4*BK3)
               ELSE
		  DLIST(4) = 0.0
               ENDIF

               CALL SAVD(NDOC,DLIST,NLIST,IRTFLG)

               IF (VERBOSE .AND. IFOUNDIT .AND. MYPID .LE. 0) THEN
                  WRITE(NOUT,6100) L,(DLIST(K),K=2,5),JDUM
6100              FORMAT (1X,I4,4(2X,F12.5),4X,I6)
               ELSEIF (VERBOSE .AND. MYPID .LE. 0) THEN
                  WRITE(NOUT,6101) L,(DLIST(K),K=2,5),JDUM
6101              FORMAT (1X,I4,4(2X,F12.5),4X,I6,'  LACKS MINIMUM!')
               ENDIF

               IF (NSEL_USED .GT. 0) THEN
C                 OUTPUT TO SPIDER'S REGISTERS NEEDED LATER

                  IF (L .GE. 3 .AND. 
     &                FSCLAST  .GE. FSCCUT .AND.
     &                DLIST(4) .LT. FSCCUT) THEN

C                     CROSSED FSCCUT GOING DOWN
                      XPREV   = L - 1     ! LAST INDEX ABOVE CUTOFF

                      FSCPREV = FSCLAST
                      FSCNOW  = DLIST(4)

                      SPFPREV = SPFLAST
                      SPFNOW  = DLIST(2)
                  ENDIF
               ENDIF
            ENDIF
         ENDDO

         CALL SAVDC
         CLOSE(NDOC)

         IF (NSEL_USED .GT. 0) THEN
C           OUTPUT TO SPIDER'S REGISTERS NEEDED
            IF (XPREV .GT. 0) THEN
               FINTERP   = (FSCCUT - FSCPREV) / (FSCNOW - FSCPREV)  

               XINTERP   = XPREV   + FINTERP * (1)
               SPFINTERP = SPFPREV + FINTERP * (SPFNOW - SPFPREV)          
            ELSE
               XINTERP   = INC
               SPFINTERP = DLIST(4)
            ENDIF

            CALL REG_SET_NSEL(1,2,XINTERP,SPFINTERP, 0.0,0.0,0.0,IRTFLG)
         ENDIF

         END
