C ++********************************************************************
C                                                                      *
C                                                                      *
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                                                                      *
C                                                                      *
C  PURPOSE:                                                            *
C                                                                      *
C INPUT:
C     P3DREF(3,NTPT)= POINTS IN REFERENCE IMAGE
C     P3D1(3,NTPT)= POINTS IN VIEWS TO BE ALIGNED, on output
C                   aligned points and the errors
C     ERRORP(NTPT) - error per point (output)
C     PHI= ANGLE TO ROTATE VIEW ABOUT Z (IN PLANE)
C     THETA= TILT ANGLE OF VIEW (ABOUT Y AXIS)
C     PSI= ANGLE TO TURN ABOUT NEW Z AXIS
C     SCAL= MULTIPLICATIVE SCALING FACTOR
C  PARAMETERS:                                                         *
C                                                                      *
C        0         2         3         4         5         6         7 *
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************

      SUBROUTINE MRERROR
     &	(P3DREF,P3D1,ERRORP,NTPT,PHI,THETA,PSI,SHIFT,SCAL,ERG,NOUT)

      DIMENSION  P3DREF(3,NTPT),P3D1(3,NTPT),ERRORP(NTPT)
      DIMENSION  SHIFT(3),CNEW(3)

      CZ=COS(PHI)
      SZ=SIN(PHI)
      CY=COS(THETA)
      SY=SIN(THETA)
      CNZ=COS(PSI)
      SNZ=SIN(PSI)

      WRITE(NOUT,*)  '  Reference point coordinates'
      WRITE(NOUT,401)  (M,(P3DREF(I,M),I=1,3),M=1,NTPT)
401   FORMAT(1X,I5,2X,3F12.3)
      WRITE(NOUT,*)  '  Aligned point coordinates and the errors'
      ERG = 0.0
      DO M=1,NTPT
         CNEW(1) =
     &    SCAL*((CNZ*CZ*CY - SZ*SNZ)*P3D1(1,M)
     &    + (SZ*CY*CNZ + SNZ*CZ)*P3D1(2,M) - CNZ*SY*P3D1(3,M))

         CNEW(2) =
     &    SCAL*(-(SZ*CNZ + CZ*SNZ*CY)*P3D1(1,M)
     &    + (CZ*CNZ - SNZ*CY*SZ)*P3D1(2,M) +SNZ*SY*P3D1(3,M))

         CNEW(3) =
     &     SCAL*(SY*CZ*P3D1(1,M) + SY*SZ*P3D1(2,M) + CY*P3D1(3,M))

         ERT = 0.0
         DO J=1,3
            ERT=ERT+(CNEW(J)-P3DREF(J,M))**2
C           OVERWRITE INPUT POINTS WITH ALIGNED POINTS
            P3D1(J,M)=CNEW(J)
         ENDDO
         ERRORP(M)=SQRT(ERT)
         WRITE(NOUT,402)  M,CNEW,ERRORP(M)
402      FORMAT(1X,I5,2X,4F12.3)
         ERG=ERG+ERT
      ENDDO
      ERG=SQRT(ERG/NTPT)
      WRITE(NOUT,*)  '  Total error=',ERG
      WRITE(NOUT,*)  ' '

      END