PROGRAM CALLTARGET
C
C Program CALLTARGET is used to generate target arrays using target
C generation routines employed by DDSCAT.
C
C B.T.Draine, Princeton Univ. Obs., 91/1/1
C History:
C 91/05/24 (BTD): Modified to allow entry target names in lower case
C 91/07/16 (BTD): Replaced DO...ENDDO with DO #...# CONTINUE
C 92/12/16 (BTD): Added IDVSHP to argument list of TARGET
C 93/01/07 (BTD): Added options TWOELL and TWOAEL
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C Adjustable parameters:
C    MXNAT=max.number of dipoles in target
      INTEGER MXNAT
      PARAMETER(MXNAT=200000)
C Derived parameters:
      INTEGER MXN3
      PARAMETER(MXN3=3*MXNAT)
C Variables:
      CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6
      INTEGER IDVSHP,IDVOUT,IOSHP,J,JJ,NAT,NAT3
      INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT),ICOMP(MXNAT,3)
      REAL ERR,F,PI,RGYR2,RXY,RXZ,RYX,RYZ,RZX,RZY,
     &     X,XCM,Y,YCM,Z,ZCM
      REAL A1(3),A2(3),SHPAR(6)
      DOUBLE PRECISION DX,DY,DZ,X2,Y2,Z2
C
C***********************************************************************
      DATA IDVOUT/0/,IOSHP/15/
      PI=4.*ATAN(1.)
C Define variables CFLSHP and IDVSHP to please flint:
      CFLSHP='none'
      IDVSHP=999
 1000 WRITE(IDVOUT,*)'What shape? Choices:'
      WRITE(IDVOUT,*)'  TETRAH = regular tetrahedron'
      WRITE(IDVOUT,*)'  SPHERE = sphere'
      WRITE(IDVOUT,*)'  ELLIPS = ellipsoid'
      WRITE(IDVOUT,*)'  CYLNDR = cylinder'
      WRITE(IDVOUT,*)'  RCTNGL = rectangular solid'
      WRITE(IDVOUT,*)'  HEXGON = hexagonal prism'
      WRITE(IDVOUT,*)'  TWOELL = two touching isotropic ellipsoids'
      WRITE(IDVOUT,*)'  TWOAEL = two touching anisotropic ellipsoids'
      WRITE(IDVOUT,*)'  THRELL = three collinear touching ellipsoids'
      WRITE(IDVOUT,*)'  THRAEL = three collinear touching anisotropic ',
     &               'ellipsoids'
      READ(*,*)CSHAPE
      IF(CSHAPE.EQ.'TETRAH'.OR.CSHAPE.EQ.'tetrah')THEN
         CSHAPE='TETRAH'
         DO 1100 J=1,100
            WRITE(IDVOUT,*)'Enter tetrahedron side length ',
     &                     ' (0 to stop)'
            READ(*,*)SHPAR(1)
            SHPAR(2)=SHPAR(1)
            SHPAR(3)=SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(*,6000)SHPAR(1),NAT
 1100    CONTINUE
         STOP
      ELSEIF(CSHAPE.EQ.'SPHERE'.OR.CSHAPE.EQ.'sphere')THEN
         CSHAPE='SPHERE'
         DO 1200 J=1,100
            WRITE(IDVOUT,*)'Enter sphere diameter (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)stop
            SHPAR(2)=SHPAR(1)
            SHPAR(3)=SHPAR(1)
            CSHAPE='ELLIPS'
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            RGYR2=0.
            XCM=0.
            YCM=0.
            ZCM=0.
            DO 1150 JJ=1,NAT
               X=IX(JJ)
               Y=IY(JJ)
               Z=IZ(JJ)
               XCM=XCM+X
               YCM=YCM+Y
               ZCM=ZCM+Z
               RGYR2=RGYR2+X*X+Y*Y+Z*Z+0.25
 1150       CONTINUE
            XCM=XCM/REAL(NAT)
            YCM=YCM/REAL(NAT)
            ZCM=ZCM/REAL(NAT)
            RGYR2=RGYR2/REAL(NAT)-XCM**2-YCM**2-ZCM**2
            F=(5./3.)*(4.*PI/REAL(3.*NAT))**(2./3.)*RGYR2
            WRITE(*,6100)SHPAR(1),SHPAR(2),SHPAR(3),NAT,F
 1200       CONTINUE
         STOP
      ELSEIF(CSHAPE.EQ.'ELLIPS'.OR.CSHAPE.EQ.'ellips')THEN
         CSHAPE='ELLIPS'
         DO 1300 J=1,100
            WRITE(IDVOUT,*)'Enter ellipsoid diameters xv,yv,zv',
     &                     ' (0,0,0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
            IF(SHPAR(1).LE.0.)stop
            CSHAPE='ELLIPS'
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            XCM=0.
            YCM=0.
            ZCM=0.
            DO 1250 JJ=1,NAT
               X=IX(JJ)
               Y=IY(JJ)
               Z=IZ(JJ)
               XCM=XCM+X
               YCM=YCM+Y
               ZCM=ZCM+Z
 1250       CONTINUE
            XCM=XCM/REAL(NAT)
            YCM=YCM/REAL(NAT)
            ZCM=ZCM/REAL(NAT)
            X2=0.
            Y2=0.
            Z2=0.
            DO 1260 JJ=1,NAT
               DX=IX(JJ)-XCM
               DY=IY(JJ)-YCM
               DZ=IZ(JJ)-ZCM
               X2=X2+DX*DX
               Y2=Y2+DY*DY
               Z2=Z2+DZ*DZ
 1260       ENDDO
            X2=X2/NAT+1./12.
            Y2=Y2/NAT+1./12.
            Z2=Z2/NAT+1./12.
            ERR=20.*PI*SQRT(5.*X2*Y2*Z2)/(3.*NAT)-1.
            RXY=SQRT(X2/Y2)
            RYZ=SQRT(Y2/Z2)
            RZX=SQRT(Z2/X2)
            RYX=1./RXY
            RZY=1./RYZ
            RXZ=1./RZX
            X2=2.*SQRT(5.*X2)
            Y2=2.*SQRT(5.*Y2)
            Z2=2.*SQRT(5.*Z2)
            WRITE(*,6200)SHPAR(1),SHPAR(2),SHPAR(3),NAT,X2,Y2,Z2,ERR,
     &                   RYX,RZX,RYZ,RXY,RXZ,RZY
 1300    CONTINUE
         STOP
      ELSEIF(CSHAPE.EQ.'CYLNDR'.OR.CSHAPE.EQ.'cylndr')THEN
         CSHAPE='CYLNDR'
         DO 1400 J=1,100
            WRITE(IDVOUT,*)'Enter cylinder length (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            WRITE(IDVOUT,*)'Enter cylinder diameter'
            READ(*,*)SHPAR(2)
            WRITE(IDVOUT,*)'Enter 1,2,3 for axis along x,y,z'
            READ(*,*)SHPAR(3)
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 1400    CONTINUE
         STOP
      ELSEIF(CSHAPE.EQ.'RCTNGL'.OR.CSHAPE.EQ.'rctngl')THEN
         CSHAPE='RCTNGL'
         DO 1500 J=1,100
            WRITE(IDVOUT,*)'Enter x-length (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            WRITE(IDVOUT,*)'Enter y-length'
            READ(*,*)SHPAR(2)
            WRITE(IDVOUT,*)'Enter z-length'
            READ(*,*)SHPAR(3)
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 1500    CONTINUE
      ELSEIF(CSHAPE.EQ.'HEXGON'.OR.CSHAPE.EQ.'hexgon')THEN
         CSHAPE='HEXGON'
         WRITE(IDVOUT,*)'Select orientation (1-6):'
         WRITE(IDVOUT,*)'1 or 2 for axis in x direction'
         WRITE(IDVOUT,*)'3 or 4             y'
         WRITE(IDVOUT,*)'5 or 6             z'
         READ(*,*)SHPAR(3)
         DO 1600 J=1,100
            WRITE(IDVOUT,*)'Enter length (along symmetry axis)',
     &                     ' (0 to stop)'
            READ(*,*)SHPAR(1)
            IF(SHPAR(1).LE.0.)STOP
            WRITE(IDVOUT,*)'Enter hexagon diameter = 2*side width'
            READ(*,*)SHPAR(2)
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 1600    CONTINUE
      ELSEIF(CSHAPE.EQ.'TWOELL'.OR.CSHAPE.EQ.'twoell')THEN
         CSHAPE='TWOELL'
         DO 1700 J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
               IF(SHPAR(1).EQ.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 1700    CONTINUE
      ELSEIF(CSHAPE.EQ.'TWOAEL'.OR.CSHAPE.EQ.'twoael')THEN
         CSHAPE='TWOAEL'
         DO 1800 J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
               IF(SHPAR(1).EQ.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 1800    CONTINUE
      ELSEIF(CSHAPE.EQ.'THRELL'.OR.CSHAPE.EQ.'threll')THEN
         CSHAPE='THRELL'
         DO 1900 J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
               IF(SHPAR(1).EQ.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 1900    CONTINUE
      ELSEIF(CSHAPE.EQ.'THRAEL'.OR.CSHAPE.EQ.'thrael')THEN
         CSHAPE='THRAEL'
         DO 2000 J=1,100
            WRITE(IDVOUT,*)'Enter length of one ellipsoid in x, y, z',
     &                     ' directions (0 0 0 to stop)'
            READ(*,*)SHPAR(1),SHPAR(2),SHPAR(3)
               IF(SHPAR(1).EQ.0)STOP
            CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT,MXN3,NAT3)
            WRITE(IDVOUT,*)'NAT=',NAT
 2000    CONTINUE
      ELSE
         WRITE(IDVOUT,*)'Do not know how to make shape = ',CSHAPE
         GOTO 1000
      ENDIF
 6000 FORMAT(' AX=',F7.4,' NAT=',I7)
 6100 FORMAT(' AX=',F7.4,' AY=',F7.4,' AZ=',F7.4,' NAT=',I7,' F=',F10.6)
 6200 FORMAT('    XV =',F7.4,'    YV =',F7.4,'    ZV =',F7.4,
     &       ' NAT=',I8,/,
     &       ' 2a_eff=',F7.4,' 2b_eff=',F7.4,' 2c_eff=',F7.4,
     &       ' err=',F8.7,/,
     &       '   y:x =',F7.4,'   z:x =',F7.4,'   y:z =',F7.4,/,
     &       '   x:y =',F7.4,'   x:z =',F7.4,'   z:y =',F7.4)
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C***********************************************************************
C This version of timeit uses the system call "etime" available under
C SunOS and some other bsd-like Unix systems (e.g., Convex unix).
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C 94/06/20 (PJF) add dtime to the formal parameters
C Arguments:
      CHARACTER CMSGTM*(*)
C
C External system calls:
      REAL ETIME
      EXTERNAL ETIME
C
C Local variables:
      REAL DTIME,OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C
C External subroutines:
      EXTERNAL WRIMSG
C
C Local variables to be saved:
      SAVE CSTA,OLDTIM
C
C Data statements:
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
         WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
         CALL WRIMSG('TIMEIT',CMSGNM)
         CSTA='ONE'
         DTIME=ETIME(TARRAY)-OLDTIM
         WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
         CALL WRIMSG('TIMEIT',CMSGNM)
      ENDIF
      RETURN
      END

      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_convex
C
C This version of timeit uses the system call "etime" available under
C Convex unix and other bsd-like Unix systems.
C
C 94/06/20 (PJF) add dtime to the formal parameters
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      CHARACTER CMSGTM* (*)
C Local variables:
      REAL DTIME,OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C External Subroutines:
      EXTERNAL WRIMSG
C External Functions:
      REAL ETIME
      EXTERNAL ETIME
C***********************************************************************
      DATA CSTA/'ONE'/
      IF(CSTA.EQ.'ONE')THEN
              CSTA='TWO'
              OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
          WRITE(CMSGNM,FMT='(A, A)')'Timing results for: ', CMSGTM
          CALL WRIMSG('TIMEIT',CMSGNM)
          CSTA='ONE'
          DTIME=ETIME(TARRAY)-OLDTIM
          WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
          CALL WRIMSG('TIMEIT',CMSGNM)
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_cray
C
C CRAY version. Call to second is executed.
C 94/06/20 (PJF) add dtime to the formal parameters
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C     .. Scalar Arguments ..
      CHARACTER CMSGTM* (*)
C     ..
C     .. Local Scalars ..
      REAL DTIME, OLDTIM
      CHARACTER CSTA*3, CMSGNM*70
C
C     ..
C
      save csta, time1
C     .. Data statements ..
      DATA CSTA/'ONE'/
C     ..
      IF (CSTA .EQ. 'ONE') THEN
              CSTA = 'TWO'
              call second(time1)
      ELSE IF (CSTA .EQ. 'TWO') THEN
          WRITE (CMSGNM,FMT='(A, A)') 'Timing results for: ', CMSGTM
          CALL WRIMSG('TIMEIT',CMSGNM)
          CSTA = 'ONE'
         call second(time2)
         dtime = time2-time1 
          WRITE (CMSGNM,FMT='(G15.6, A)') DTIME,
     $          ' CPU time (in secs) '
          CALL WRIMSG('TIMEIT',CMSGNM)
      END IF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C timeit_hp
C
C This version uses the HP-AUX system calls SYSCONF and TIMES
C 94/06/20 (PJF) add dtime to the formal parameters
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C External system calls:
$ALIAS sysconf = 'sysconf' (%val)
      INTEGER SYSCONF,TIMES
      EXTERNAL SYSCONF,TIMES
C
C Arguments:
      CHARACTER CMSGTM*(*)
C
C Local variables:
      INTEGER*2 NAME
      INTEGER CLK_TCK,I
      INTEGER TMS(4)
      REAL DTIME, OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
C
C External subroutines:
      EXTERNAL WRIMSG

C Data statements:
      DATA CSTA/'ONE'/
C     ..
      IF(CSTA.EQ.'ONE')THEN
         NAME=2
         CLK_TCK=SYSCONF(NAME)
         OLDTIM=REAL(TMS(1)+TMS(2))/REAL(CLK_TCK)
      ELSEIF(CSTA.EQ.'TWO')THEN
         WRITE(CMSGNM,FMT='(A, A)')'Timing results for: ', CMSGTM
         CALL WRIMSG('TIMEIT',CMSGNM)
         CSTA='ONE'
         I=TIMES(TMS)
         DTIME=REAL(TMS(1)+TMS(2))/REAL(CLK_TCK)-OLDTIM
         WRITE(CMSGNM,FMT='(G15.6,A)')DTIME,' CPU time (in secs) '
         CALL WRIMSG('TIMEIT',CMSGNM)
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_ibm6000
C 94/06/20 (PJF) add dtime to the formal parameters
C
C This version of timeit uses the system call "mclock" available under
C RISC6000/320 IBM AIX3.1
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      CHARACTER CMSGTM* (*)
C
C Local variables:
      REAL DTIME, OLDTIM
      CHARACTER CSTA*3, CMSGNM*70
      REAL TARRAY(2)
C
C External Subroutines ..
      EXTERNAL WRIMSG
C
C External Functions ..
      REAL ETIME
      EXTERNAL ETIME
C
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         I1=MCLOCK()
      ELSEIF(CSTA.EQ.'TWO')THEN
         WRITE(CMSGNM,FMT='(A,A)')'Timing results for: ',CMSGTM
         CALL WRIMSG('TIMEIT',CMSGNM)
         CSTA='ONE'
         I2=MCLOCK()
         IALL=I2-I1
         DTIME=FLOAT(IALL)/100.
         WRITE(CMSGNM,FMT='(G15.6,A)')DTIME,
     $         ' CPU time (in secs) '
         CALL WRIMSG('TIMEIT',CMSGNM)
      ENDIF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C***********************************************************************
C This version of timeit uses the SGI IRIX system call "etime"
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C 94/06/20 (PJF) add dtime to the formal parameters
C Arguments:
      CHARACTER CMSGTM*(*)
C
C External system calls:
      REAL ETIME
      EXTERNAL ETIME
C
C Local variables:
      REAL DTIME,OLDTIM
      CHARACTER CSTA*3,CMSGNM*70
      REAL TARRAY(2)
C
C External subroutines:
      EXTERNAL WRIMSG
C
C Local variables to be saved:
      SAVE CSTA,OLDTIM
C
C Data statements:
      DATA CSTA/'ONE'/
C
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
         OLDTIM=ETIME(TARRAY)
      ELSEIF(CSTA.EQ.'TWO')THEN
         WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
         CALL WRIMSG('TIMEIT',CMSGNM)
         CSTA='ONE'
         DTIME=ETIME(TARRAY)-OLDTIM
         WRITE(CMSGNM,FMT='(G15.6, A)')DTIME,' CPU time (in secs) '
         CALL WRIMSG('TIMEIT',CMSGNM)
      ENDIF
      RETURN
      END


      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_titan
C
C This version of TIMEIT is for Titan computers.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C
C 94/06/20 (PJF) add dtime to the formal parameters
C***********************************************************************
C     .. Scalar Arguments ..
      CHARACTER CMSGTM* (*)
C     ..
C     .. Local Scalars ..
      REAL DTIME, OLDTIM, TIME
      INTEGER IAD
      CHARACTER CSTA*3, CMSGNM*70
C     ..
C     .. External Subroutines ..
      EXTERNAL WRIMSG
C     ..
C     .. External Functions ..
      REAL CPUTIM
      EXTERNAL CPUTIM
C     ..
C     .. Data statements ..
      DATA CSTA/'ONE'/
C     ..
      IF (CSTA .EQ. 'ONE') THEN
          CSTA = 'TWO'
          OLDTIM = CPUTIM(0.)
      ELSE IF (CSTA .EQ. 'TWO') THEN
          WRITE (CMSGNM,FMT='(A, A)') 'TIMING RESULTS OF: ', CMSGTM
          CALL WRIMSG('TIMEIT',CMSGNM)
          CSTA = 'ONE'
          TIME = CPUTIM(0.)
          DTIME = TIME - OLDTIM
          WRITE (CMSGNM,FMT='(G15.6, A)') DTIME,
     $          ' CPU TIME (IN SECS) '
              CALL WRIMSG('TIMEIT',CMSGNM)
      END IF
      RETURN
      END
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_null
C
C This version of timeit is a dummy which does not provide any
C timing information.
C 94/06/20 (PJF) add dtime to the formal parameters
C
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      CHARACTER CMSGTM*(*)
C Local variables:
      CHARACTER CSTA*3,CMSGNM*70
C***********************************************************************
      DATA CSTA/'ONE'/
      IF(CSTA.EQ.'ONE')THEN
          CSTA='TWO'
      ELSEIF(CSTA.EQ.'TWO')THEN
          CSTA='ONE'
          WRITE(CMSGNM,FMT='(A,A)')'Timing results for: ',CMSGTM
          CALL WRIMSG('TIMEIT',CMSGNM)
          WRITE(CMSGNM,FMT='(A)')
     &         'Dummy timing routine: no timings available'
          CALL WRIMSG('TIMEIT',CMSGNM)
      ENDIF
      RETURN
      END
      PROGRAM TSTFFT
C Parameters:
      INTEGER MX,MY,MZ
        PARAMETER (MX=100, MY=100, MZ=100)
C Local variables:
      INTEGER IN,N,NX,NY,NZ
      INTEGER NF235(45)
      COMPLEX C(MX,MY,MZ),CW(MX,MY,MZ)
C***********************************************************************
C Purpose:
C Tests different 3D FFT implementations (CPU time).
C Required routines:
C cfft99f.o  cxfft3.o  cxfft3n.o  
C errmsg.o  fourx.o  getset.o  
C gpfa.o  timeit_sgi.o  tst.o  wrimsg.o gen.o
C
C History:
C 94.06.20 (PJF): Original version
C 94.12.18 (BTD): Minor changes in output formatting.
C                 Explicit declaration of all variables.
C
C Copyright (C) 1994 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C
C***********************************************************************
C Elements of NF235 are numbers.le.4096  which can be factored by 2,3,5:
       DATA NF235/2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,
     &    45,48,50,54,60,64,72,75,80,81,90,96,100,108,120,125,128,135,
     &    144,150,160,162,180,192,200/
       PRINT *,"Calculation time for different 3-D FFT's methods"
       PRINT *,"================================================ " 
       PRINT *,"MACHINE         LVR=??? "
       PRINT *,
     &  "                    BRENNER   TEMPERTON(OLD)  TEMPERTON(NEW)"
       PRINT *," "
       open(unit=20,file='err.dat',status='unknown')
       open(unit=21,file='res.dat',status='unknown')
       WRITE(21,6021)
 6021 FORMAT(' CPU time (sec) for different 3-D FFT methods',/,
     &   '================================================',/,
     &   'MACHINE         LVR=???',/,/,
     &   '    NX  NY  NZ          BRENNER   TEMPERTON(OLD)  ',
     &   'TEMPERTON(NEW)')
      rewind(20)
      DO 10 IN=1,31
         N=NF235(IN)
         NX=N
         NY=N
         NZ=N
         CALL CASE(C,CW,NX,NY,NZ)
10    CONTINUE
      STOP
      END
      SUBROUTINE CASE(C,CW,NX,NY,NZ)
C Arguments:
      INTEGER NX,NY,NZ
      COMPLEX C(NX,NY,NZ),CW(NX,NY,NZ)
C Local variables:
      INTEGER NN(3)
      REAL T0,T1,T2,T3
      COMPLEX C1, C2, C3, C4, C5, C6
C***********************************************************************
C Purpose: interface to 3D FFT codes (for use in FFT tests)
C History:
C 94.06.20 (PJF): Original version
C 94.12.19 (BTD): Explicit declaration of all variables.
C
C Copyright (C) 1994, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
***********************************************************************C
C--------- BRENNR
         NN(1)=NX
         NN(2)=NY
         NN(3)=NZ
         call  gen(c,nx,ny,nz)
         CALL TIMEIT('BREN',T0)         
         CALL FOURX(C,NN,3,1)
         CALL TIMEIT('BREN',T1)
         c1=c(1,1,1)
         c2=c(2,1,1)

C--------- TEMP-OLD
         call  gen(c,nx,ny,nz)
         CALL TIMEIT('OLD',T0)
         CALL CXFFT3(C,NX,NY,NZ,1,CW)
         CALL TIMEIT('OLD',T2)
         c3=c(1,1,1)
         c4=c(2,1,1)

C--------- TEMP-NEW
         call  gen(c,nx,ny,nz)
         CALL TIMEIT('NEW',T0)
         CALL CXFFT3n (C,NX,NY,NZ,1)
         CALL TIMEIT('NEW',T3)
         c5=c(1,1,1)
         c6=c(2,1,1)
C--------- Print
C         WRITE(21, '(1x,3i5,3g15.6)') NX,NY,NZ,T1,T2,T3 
         WRITE(21, '(1x,3i5,3F15.6)') NX,NY,NZ,T1,T2,T3 
         write(20,100) nx, c1,c2
         write(20,100) nx, c3,c4
         write(20,100) nx, c5,c6
 100     format(1x,i5,4g15.4)
         RETURN
         END
      subroutine gen(c,nx,ny,nz)
C Arguments:
      INTEGER NX,NY,NZ
      complex c(nx,ny,nz)
C Local variables:
      INTEGER I1,I2,I3
      REAL F
C Purpose:
C generates a 3D matrix for use in 3D FFT CPU timmings
      do 103 i1=1, nx
         do 102 i2=1, ny
            do 101 i3=1, nz
c              f=i1*i2*i3
               f=sin(float(i1))+sin(float(i2))+sin(float(i3))
               c(i1,i2,i3)=cmplx(f*1.,f*0.5)
 101        CONTINUE
 102     CONTINUE
 103  continue
      return
      end