C ++********************************************************************
C                                                                      *
C EPUR4     REMOVED FROM HPLAN.FOR FILE     APRIL 89 al                    *
C           ID(IDIM,NID)                    12/7/93 ML                                  *             
C           UPDATED SOME                    AUG 1999 ARDEAN LEITH                           *
C           SIMPLIFIED LOGIC FOR OUTLIERS   DEC 2005 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   EPUR4(IDIM,NPTS, X,Y,ID, MOD, PEX, KP,KLIC,IRTFLG,NDAT) 
C 
C   POINTS WHICH FALL BEYOND PEX STANDARD DEVIATIONS ARE COLLAPSED ONTO         
C   THE FRAME OF THE GRAPH.
C                                                     
C   OUTPUT: THE KP MODIFIED POINTS HAVE THEIR INDICES IN KLIC()                 
C   IF THERE ARE MORE THAN 264 POINTS ON THE FRAME THE GRAPH IS ABORTED         
C   IF MOD=1 LABELS HAVE A1 FORMAT - IF MOD=4 LABELS HAVE A7 FORMAT 
C            
C   WARNING: X(*) AND Y(*) ARE DESTROYED IF KP .NE. 0                           
C
C   CALLED BY: HPLAN, HISMAP,  HISMAP4
C
C **********************************************************************

        SUBROUTINE EPUR4(IDIM,NPTS,X,Y,ID,MOD,PEX,KP,KLIC,IRTFLG,NDAT)

        CHARACTER(LEN=7) :: ID(IDIM)
        REAL             :: X(IDIM), Y(IDIM)
        INTEGER          :: KLIC(371)

        DATA  SEUIL / 1.0 E-07 /

        IRTFLG  = 1
        IF (PEX .LE. 0.0) RETURN

        NPTSP1  = NPTS + 1
        IF (NPTSP1 .GT. IDIM) RETURN
                               
        X(NPTSP1) = 0.0
        Y(NPTSP1) = 0.0

        NPTSP1    = NPTS
        SX        = 0.0
        SY        = 0.0
        DO I = 1,NPTS                                                      
           SX = SX + X(I) * X(I)
           SY = SY + Y(I) * Y(I)
        END DO

        SX    = SQRT(SX/FLOAT(NPTS))
        SY    = SQRT(SY/FLOAT(NPTS))
        PX    = PEX * SX 
        PY    = PEX * SY
        KP    = 0 
        PIN   = PEX + 20.0

        DO  I = 1, NPTS                                                      
          IF (ABS(X(I)) .GT. PX .OR. ABS(Y(I)) .GT. PY) THEN
C           POINT IS BEYOND SD BOUNDARY REQUESTED 
            KP = KP + 1 
            IF (MOD.NE.1 .AND. KP .EQ. 1) THEN
C              WRITE HEADER FOR OUTLIER LISTING
               WRITE(NDAT,120) PEX 
120            FORMAT (//15X,'WARNING',//, 
     &          ' THE FOLLOWING POINTS LOCATED BEYOND ',F5.1,
     &          ' STANDARD DEVIATIONS FROM',/,
     &          ' THE ORIGIN ARE PLACED ON THE EDGES OF THE MAP',/)

               WRITE(NDAT,130)
130            FORMAT(' ',29('*'))
            ENDIF
            IF (KP .GT. 264) RETURN
                               
            KLIC(KP) = I 
            IF (MOD .NE. 1) WRITE(NDAT,110) ID(I),X(I),Y(I)
110         FORMAT ('  *',A7,1X,'*',2(F13.5,8X,'*') )
c            X(I)  = X(I) / PIN
c            Y(I)  = Y(I) / PIN 

             IF (ABS(X(I)) .GT. PX)  X(I)  = SIGN(PX, X(I))
             IF (ABS(Y(I)) .GT. PY)  Y(I)  = SIGN(PY, Y(I))
          ENDIF
        ENDDO
        IF (KP .NE. 0) WRITE(NDAT,130)

        IRTFLG  = 0
        END 
                                                                    
#ifdef NEVER
C       REMOVED DEC 05 AS HAD BUG WITH SOME DISTRIBUTIONS
        IF (KP .NE. 0)  THEN
C         RESCALE PERIPHERY VALUES 
C         (INEFFICIENT!!! al)

          IF (MOD .EQ. 4) WRITE(NDAT,(' ',29('*')))
          CALL BORNS(NPTSP1,Y,YMIN,YMAX)
          CALL BORNS(NPTSP1,X,XMIN,XMAX)
          IF (XMIN .EQ. 0.) XMIN = SEUIL
          IF (XMAX .EQ. 0.) XMAX = SEUIL 
          T1    = YMAX / XMAX
          T2    = YMAX / XMIN 
          T3    = YMIN / XMIN
          T4    = YMIN / XMAX

          DO  K = 1, KP 
            KLK = KLIC(K)
            U   = X(KLK)  
            V   = Y(KLK)  
            IF (U .EQ. 0.) U = SEUIL
            T = V / U 
            A = YMAX 
            IF (U .GT. 0.0 .AND. (T4.LT.T .AND. T.LT.T1)) THEN
              A = XMAX  
              U = A 
              V = A * T
            ELSEIF (U .LT. 0.0 .AND.(T2.LT.T .AND. T.LT.T3)) THEN
              A = XMIN 
              U = A 
              V = A * T
            ELSE
              IF (V .LT. 0.) A = YMIN 
              V = A
              IF (T .EQ. 0.) T = SEUIL
              U = A / T
            END IF
                   
            IF (V .GT. YMAX) V = YMAX 
            IF (U .GT. XMAX) U = XMAX  
            X(KLK) = U  
            Y(KLK) = V 
          END DO

          WRITE(NDAT,140) KP, PEX
140       FORMAT ('  ',I4,' POINTS LOCATED BEYOND ',F5.1,
     &            ' STANDARD DEVIATIONS ARE PLACED ON ',
     &            'EDGES OF THE MAP')
        ENDIF
#endif