C ++********************************************************************
C                                                                      *
C CENGR3.F                                                               *
C                                       LONG FILENAMES JAN 89 al
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  CENGR3(LUN3)                                                             *
C                                                                      *
C  PURPOSE:                                                            *
C     CALCULATES CENTER OF GRAVITY AND RADIUS OF GYRATION OF A DENSITY
C    MASS STORED IN A THREE-D FILE. FOR REFERENCE, SEE INTERNATIONAL
C    TABLES OF CRYSTALLOGRAPHY, VOL.III P. 327.     --- JF 6/01/83
C
C  PARAMETERS:                                                         *
C
C  SUPPORT_ROUTINE
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************

        SUBROUTINE CENGR3(LUN3)

        INCLUDE 'CMBLOCK.INC'
        INCLUDE 'CMLIMIT.INC' 
 
        CHARACTER(LEN=MAXNAM)   ::  FILNAM
        COMMON /COMMUN/FILNAM
        REAL, ALLOCATABLE, DIMENSION(:) :: AIMG

        CHARACTER     NULL

        MAXIM = 0
        CALL OPFILEC(0,.TRUE.,FILNAM,LUN3,'O',ITYPE,
     &               NSAM,NROW,NSLICE,MAXIM,
     &               'INPUT',.FALSE.,IRTFLG)
        IF (IRTFLG .NE. 0) RETURN
	
        IF (ITYPE .NE. 3 .AND. ITYPE .NE. 1)  THEN
           CLOSE(LUN3)
           CALL  ERRT(2,'CG 3',NE)
           RETURN
        ENDIF

        CALL RDPRM(TH,NOT_USED,'THRESHOLD')

C       FIRST CALCULATE CENTER OF GRAVITY
        XM=0.0
        YM=0.0
        ZM=0.0
        BT=0.0

        ALLOCATE (AIMG(NSAM), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'CG 3, AIMG',IER)
           RETURN
        ENDIF
        DO  L=1,NSLICE
           I1=(L-1)*NROW
           DO  I=1,NROW
              CALL REDLIN(LUN3,AIMG,NSAM,I+I1)
              DO  K=1,NSAM
                 IF (AIMG(K).GT.TH) THEN
                    B=AIMG(K)
                    BT=B+BT
                    XM=XM+B*K
                    YM=YM+B*I
                    ZM=ZM+B*L
                 ENDIF  
              ENDDO
           ENDDO
        ENDDO
        ZM=ZM/BT
        YM=YM/BT
        XM=XM/BT
        

C       NOW CALCULATE RADIUS OF GYRATION, RG
        IF (ITYPE.EQ.3) THEN
           BT=0.0
           RG=0.0
           DO  L=1,NSLICE
              I1=(L-1)*NROW
              DO  I=1,NROW
                 REST=(YM-I)**2+(ZM-L)**2
                 CALL REDLIN(LUN3,AIMG,NSAM,I+I1)
                 DO K=1,NSAM
                    IF (AIMG(K).GT.TH) THEN
                       RG=RG+AIMG(K)*((XM-K)**2+REST)
                       BT=BT+AIMG(K)
                    ENDIF
                 ENDDO
              ENDDO
           ENDDO
           RG=SQRT(RG/BT)
        ENDIF        

        XM=XM-NSAM/2-1
        YM=YM-NROW/2-1
        ZM=ZM-NSLICE/2-1


        IF (ITYPE.EQ.3) THEN
           WRITE(NOUT,201)XM,YM,ZM,RG
201        FORMAT(' ** COMPONENTS OF CENTER OF GRAVITY:'//,
     &     '           X = ',F7.2,'  Y = ',F7.2,'  Z = ',F7.2,
     &     //'      RADIUS OF GYRATION = ',F12.2/)

           CALL REG_SET_NSEL(1,4,XM,YM,ZM,RG,0.0,IRTFLG)

        ELSE
           WRITE(NOUT,202)XM,YM
202        FORMAT(' ** COMPONENTS OF CENTER OF GRAVITY:'//,
     &     '           X = ',F7.2,'  Y = ',F7.2)

           CALL REG_SET_NSEL(1,2,XM,YM,0.0,0.0,0.0,IRTFLG)
        ENDIF

        DEALLOCATE (AIMG)

        END