C++************************************************************************
C
C PR3DB.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 PURPOSE: CALCULATE THE 3-D PHASE RESIDUE OUTSIDE MISSING CONE, 
C          PR OF FOURIER RINGS(RADIUS, DIRECTION RELATIVE TO Z) AND 
C          OF FOURIER SHELLS(RADIUS) IS CALCULATED. 
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*********************************************************************

        SUBROUTINE PR3DB(A,B,PR,AMP,CSUM1,LR,CSUM,CSUM2,
     &     AVSUM,LSD,NSAM,NROW,NSLICE,DSCALE,NSCALE,SCALE1,SSANG,
     &     INC,Y1,WI,SER)

        DIMENSION A(LSD,NROW,NSLICE),B(LSD,NROW,NSLICE)
        DIMENSION PR(NSCALE,INC),AMP(NSCALE,INC),AVSUM(NSCALE,INC)
        DIMENSION CSUM1(INC),CSUM2(INC),LR(INC)
        DIMENSION CSUM(INC)
        CHARACTER*1 SER

        PARAMETER (QUADPI = 3.141592653589793238462643383279502884197)
        PARAMETER (RAD_TO_DGR = (180.0/QUADPI))
        PARAMETER (DGR_TO_RAD = (QUADPI/180))

        ZANG=(90.0-SSANG)*DGR_TO_RAD
        ND2=NSAM/2
        NR2=NROW/2

C       for 2D case set NS2 to one
        NS2 = MAX0(1,NSLICE/2)

        DO  L=1,INC
           LR(L)    = 0
           CSUM(L)  = 0.0
           CSUM1(L) = 0.0
           CSUM2(L) = 0.0
           DO  NSC=1,NSCALE
              PR(NSC,L)=0.0
              AVSUM(NSC,L)=0.0
              AMP(NSC,L) =0.0
           ENDDO
        ENDDO

        DO K=1,NSLICE
           IIK=(K-1)
           IF(IIK.GT.NS2)  IIK=IIK-NSLICE
           PK=(FLOAT(IIK)/FLOAT(NS2))**2
           DO J=1,NROW
              IIJ=(J-1)
              IF(IIJ.GT.NR2)  IIJ=IIJ-NROW
              PJ=(FLOAT(IIJ)/FLOAT(NR2))**2

              DO I=1,LSD,2
                III=(I-1)/2
C               skip Hermitian related values
                IF (III.GT.0 .OR.(IIK.GE.0 .AND. 
     &             (IIJ.GE.0 .OR. IIK.NE.0))) THEN
                PII = (FLOAT(III)/FLOAT(ND2))**2
                R   = SQRT(PII+PJ+PK)*0.5

		IF (R .GT. 0.0)  THEN
                   R1 = SQRT(PII+PK)*0.5          ! r1=(x^2+z^2)^0.5
                   !if(r1.gt.0.05) goto 3
                   IF (SER .EQ. 'C')  THEN
                      IF (ACOS(AMIN1(1.0,SQRT(PK)/R)).LT.ZANG) GO TO 3

                   ELSEIF (SER .EQ. 'W')  THEN
                      IF (R1 .EQ. 0.0)  THEN
                         FI = ACOS(1.0)
                      ELSE
                         FI = ACOS(AMIN1(1.0,SQRT(PK)/R1))
                      ENDIF
                      IF (FI .LT. ZANG) GO TO 3
                   ENDIF

                   L = NINT(R*2*(INC-1))+1
                   IF (L .GT. INC) GO TO 3
                   LR(L) = LR(L) + 2
                   IF (A(I,J,K) .NE. 0.0)THEN
                      PHA = ATAN2(A(I+1,J,K),A(I,J,K))*RAD_TO_DGR
                   ELSE
                      PHA = 0
                   ENDIF 
                   IF (B(I,J,K) .NE. 0.0)THEN
                      PHB = ATAN2(B(I+1,J,K),B(I,J,K))*RAD_TO_DGR
                   ELSE
                      PHB = 0.0
                   ENDIF 
                   DPH = PHA-PHB
                   IF (DPH .GT. 180.0)  DPH=360.0-DPH
                   IF (DPH .LT. -180.0) DPH=360.0+DPH

                   QA = SQRT(A(I,J,K)**2+A(I+1,J,K)**2)
                   QB = SQRT(B(I,J,K)**2+B(I+1,J,K)**2)

                   CSUM(L)=
     &                CSUM(L)+A(I,J,K)*B(I,J,K)+A(I+1,J,K)*B(I+1,J,K)
                   CSUM1(L) = CSUM1(L)+QA*QA
                   CSUM2(L) = CSUM2(L)+QB*QB

C                  SCALES
                   DPH=DPH*DPH
                   DO NSC=1,NSCALE
                      SCALE        = SCALE1+(NSC-1)*DSCALE
                      QBS          = QB*SCALE
                      AVSUM(NSC,L) = AVSUM(NSC,L)+ABS(QA-QBS)
                      PR(NSC,L)    = PR(NSC,L)+((QA+QBS)/2.)*DPH
                      AMP(NSC,L)   = AMP(NSC,L)+((QA+QBS)/2.)
                   ENDDO
		ENDIF
             ENDIF
3               CONTINUE
              ENDDO
           ENDDO
        ENDDO
        END