C ++********************************************************************
C                                                                      *
C  RPRQ.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:                                                            *
C                                                                      *
C  PARAMETERS:                                                         *
C
C IMAGE_PROCESSING_ROUTINE                                             *
C
C        0         2         3         4         5         6         7 *
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************

        SUBROUTINE RPRQ(N,B,CUBE,IPCUBE,NN,PHI,THETA,PSI,DM,RI)

        INCLUDE 'CMBLOCK.INC'

        DIMENSION  B(N,N),CUBE(*)
        INTEGER           IPCUBE(5,NN)
        DOUBLE PRECISION  CPHI,SPHI,CTHE,STHE,CPSI,SPSI
        DOUBLE PRECISION  QUADPI,DGR_TO_RAD
        DIMENSION         DM(9)
        LOGICAL           ALLOK

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

        COMMON /PAR/  LDP,NM,LDPNM

C IPCUBE: 1 - beginning
C         2 - end
C         3 - ix
C         4 - iy
C         5 - iz

        R=RI*RI

        CPHI=DCOS(DBLE(PHI)*DGR_TO_RAD)
        SPHI=DSIN(DBLE(PHI)*DGR_TO_RAD)
        CTHE=DCOS(DBLE(THETA)*DGR_TO_RAD)
        STHE=DSIN(DBLE(THETA)*DGR_TO_RAD)
        CPSI=DCOS(DBLE(PSI)*DGR_TO_RAD)
        SPSI=DSIN(DBLE(PSI)*DGR_TO_RAD)

        DM(1)=CPHI*CTHE*CPSI-SPHI*SPSI
        DM(2)=SPHI*CTHE*CPSI+CPHI*SPSI
        DM(3)=-STHE*CPSI
        DM(4)=-CPHI*CTHE*SPSI-SPHI*CPSI
        DM(5)=-SPHI*CTHE*SPSI+CPHI*CPSI
        DM(6)=STHE*SPSI
        DM(7)=STHE*CPHI
        DM(8)=STHE*SPHI
        DM(9)=CTHE     

        ALLOK = .TRUE.
        DM1   = DM(1)
        DM4   = DM(4)

c$omp parallel do private(i,j,xb,yb,xbb,ybb,iqx,iqy,dipy),shared(allok)

        DO I=1,NN
           XB = (IPCUBE(3,I)-LDP)*DM(1)+(IPCUBE(4,I)-LDP)*DM(2)+
     &          (IPCUBE(5,I)-LDP)*DM(3)
           YB = (IPCUBE(3,I)-LDP)*DM(4)+(IPCUBE(4,I)-LDP)*DM(5)+
     &          (IPCUBE(5,I)-LDP)*DM(6)

           DO J=IPCUBE(1,I),IPCUBE(2,I)

              XBB = (J - IPCUBE(1,I)) * DM1 + XB
              IQX = IFIX(XBB + FLOAT(LDPNM))

              YBB = (J - IPCUBE(1,I)) * DM4 + YB
              IQY = IFIX(YBB + FLOAT(LDPNM))
              IF (IQX.LT.1 .OR. IQX.GE.N .OR.
     &            IQY.LT.1 .OR. IQY.GE.N) THEN
                 ALLOK = .FALSE.
              ELSE

                 DIPY = YBB + LDPNM - IQY

C                EVEN FASTER VERSION
                 CUBE(J) = CUBE(J) + B(IQX,IQY) + 
     &              DIPY * (B(IQX,IQY+1)   - B(IQX,IQY)) +
     &              (XBB + LDPNM - IQX) * (B(IQX+1,IQY) - B(IQX,IQY) + 
     &              DIPY * (B(IQX+1,IQY+1) - B(IQX+1,IQY) -
     &                      B(IQX,IQY+1)   + B(IQX,IQY)))

C               FASTER VERSION
C               CUBE(J)=CUBE(J)
C    &          +B(IQX,IQY)+DIPY*(B(IQX,IQY+1)-B(IQX,IQY))
C    &          +DIPX*(B(IQX+1,IQY)-B(IQX,IQY)
C    &          +DIPY*(B(IQX+1,IQY+1)-B(IQX+1,IQY)
C    &          -B(IQX,IQY+1)+B(IQX,IQY)))
C
C               ORIGINAL VERSION
C                CUBE(J)=CUBE(J)
C     &                 +(1.0-DIPX)*(1.0-DIPY)*B(MAP(IQX,IQY))
C     &                 +     DIPX *(1.0-DIPY)*B(MAP(IQX+1,IQY))
C     &                 +(1.0-DIPX)*     DIPY *B(MAP(IQX,IQY+1))
C     &                 +     DIPX *     DIPY *B(MAP(IQX+1,IQY+1))

             ENDIF
           ENDDO
        ENDDO

        IF (.NOT. ALLOK) THEN
           WRITE(NOUT,*) '*** IQX or IQY out of range - ',
     &                   ' reduce the radius!'
           CALL ERRT(100,'BP RP',NE)
         ENDIF

         END