C++*********************************************************************
C
C MODEL3.F                         FILENAMES LENGTHENED DEC 7 88 al
C                                  CHAR. VARIABLES AUG 89 al
C                                  REWRITTEN SOME APRIL 97 al
C                                  FUNCTIONS ADDED ??  by ?
C                                  HELIX BUGS FIXED AUG 02 al
C                                  RDPRAF REMOVED   DEC 05 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    MODEL3(LUN1,NDOC,FILNAM,NSAM,NROW,NSLICE)
C
C    PURPOSE:   PREPARES 3-D MODEL DENSITY DISTRIBUTIONS
C
C    PARAMETERS:       LUN1        LOGICAL UNIT NUMBER OF VOLUME FILE
C                      NDOC        LOGICAL UNIT NUMBER OF DOCUMENT FILE
C                      FILNAM      NAME OF VOLUME FILE
C                      NSAM,NROW,NSLICE   DIMENSIONS OF FILE
C
C UPDATED 9/17/85 J.F. 2/9/88 DOCUMENT FILE OPTION
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C--********************************************************************

	SUBROUTINE MODEL3(LUN1,NDOC,FILNAM,NSAM,NROW,NSLICE)

        INCLUDE 'CMBLOCK.INC'
        INCLUDE 'CMLIMIT.INC'

	PARAMETER (MAXSIN=50)
	PARAMETER (MAXSPH=300)

	COMMON  NX(MAXSIN),NY(MAXSIN),CX(MAXSIN), NZ(MAXSIN),
     &    ZZ(MAXSPH),YY(MAXSPH),XX(MAXSPH),FL2(MAXSPH),FL22(MAXSPH),
     &    RAD1(MAXSPH), RAD2(MAXSPH), DENS(MAXSPH)

        REAL                  :: A0(NSAM)
        REAL                  :: FWA(4),ADUM(3)
	LOGICAL               :: LDOC,GETOLD,GFLAG
        CHARACTER(LEN=*)      :: FILNAM
        CHARACTER(LEN=MAXNAM) :: DOCNAM
        CHARACTER(LEN=4)      :: ANS
        CHARACTER(LEN=1)      :: NULL,GA,YN

	PARAMETER     (QUADPI = 3.1415926535897932)
	PARAMETER     (TWOPI=2*QUADPI)
      
        NULL = CHAR(0)

	XOFF = 0.
	YOFF = 0.
	ZOFF = 0.

       WRITE(NOUT,100)
 100   FORMAT(
     &      ' .MENU: (B)LANK    -- CONSTANT DENSITY VOLUME'/
     &      '        (C)YLINDER -- SET OF CYLINDERS'/
     &      '        (G)AUSSIAN -- GAUSSIAN FUNCTION'/
     &      '        (H)ELIX    -- HELIX OF SPHERES'/
     &      '        (HA)       -- HELIX OF SPHERES, ADD DENSITIES'/
     &      '        (NUM)BERS  -- LINE NUMBERS'/
     &      '        (R)ANDOM   -- RANDOM STATISTICS'/
     &      '        (S)INE     -- SET OF SINE WAVES'/
     &      '        (SP)HERE   -- SET OF SPHERES'/
     &      '        (SPA)      -- SET OF SPHERES, ADD DENSITIES'/
     &      '        (SPV)      -- SET OF SPHERES, VARIABLE DENSITIES'/
     &      '        (T)EST     -- 3D SINE WAVE'/
     &      '        (W)EDGE    -- DENSITY WEDGE'/)

	CALL RDPRMC(ANS,NCHAR,.TRUE.,
     &      'B/C/G/H/HA/NUM/R/S/SP/SPA/SPV/T/W', NULL,IRT)

	CALL CHKINPQ('T,S,SP,SPA,SPV,W,B,H,HA,C,NUM,R,G$',
     &              ANS(:NCHAR),ICALL)

	GOTO(5500,200,300,2000,2000,2000,1500,4100,2600,2600,
     &       4600,6000,7000,8000) ,ICALL

C     TEST & SINE WAVES *********************************** TEST & SINE

  200 NS    = 1
      CX(1) = 1.0
      NX(1) = 2
      NY(1) = 2
      NZ(1) = 2
      GOTO 1100

  300 CALL RDPRMI(NS,IDUM,NOT_USED,'NUMBER OF SINE WAVES')
      IF (NS .GT. MAXSIN) THEN
         WRITE(NOUT,400)MAXSIN
  400    FORMAT(' *** BUFFER SPACE RESTRICTED TO ',I3,' SINE WAVES')
         NS = MAXSIN
      ENDIF

      AV = 0.0
      DO I = 1,NS
         CALL RDPRA('REL. AMPL., SPATIAL FREQUENCY(KX,KY,KZ)',
     &         4,0,.FALSE.,FWA,NVAL,IRTFLG)
         IF (IRTFLG .NE. 0) RETURN

         CX(I) = FWA(1)
         NX(I) = FWA(2)
         NY(I) = FWA(3)
         NZ(I) = FWA(4)
         AV    = AV + CX(I)
      ENDDO

      DO I = 1,NS
         CX(I) = CX(I) * 2.0 / AV
         WRITE(NOUT, 900)I,CX(I),NX(I),NY(I),NZ(I)
  900    FORMAT(1X,I5,F5.2,3I6)
      ENDDO

 1100 CONTINUE

      DO L = 1,NSLICE
         PHASEZ = FLOAT(L-1)*TWOPI/FLOAT(NSLICE)                 
         DO I=1,NROW
           PHASE = FLOAT(I-1)*TWOPI/FLOAT(NROW)
           A0    = 0.
           DO J=1,NSAM
             DO K = 1,NS
             A0(J)=CX(K)*SIN(FLOAT(J-1)*TWOPI*FLOAT(NX(K))/FLOAT(NSAM)
     &              +PHASE*FLOAT(NY(K))+PHASEZ*FLOAT(NZ(K)))+A0(J)
             ENDDO
           ENDDO
           CALL WRTLIN(LUN1,A0,NSAM,I+(L-1)*NROW)
        ENDDO
      ENDDO
      RETURN

C     WEDGE **************************************************** WEDGE

 1500  SCF = 2./FLOAT(NROW+NSAM+NSLICE)
       DO  L= 1,NSLICE
         FI = FLOAT(L) *SCF
         DO  I=1,NROW
           FI = FI + SCF*FLOAT(I)
           DO  J = 1,NSAM
             A0(J) = FI + SCF*FLOAT(J)
 	   ENDDO
           CALL WRTLIN(LUN1,A0,NSAM,(L-1)*NROW+I)
 	 ENDDO
       ENDDO

      RETURN

C     SPHERES ************************************************* SPHERES

 2000   GETOLD = .FALSE.
        CALL RDPRMC(YN,NCHAR,.TRUE.,
     &        'GET COORDINATES FROM DOCUMENT FILE (Y/N)',NULL,IRTFLG)
        IF (IRTFLG .EQ. -1) RETURN

        LDOC = .FALSE.
        IF (YN .NE. 'N') THEN
C          USE DOC FILE FOR COORDINATES
           LDOC  = .TRUE.
           NOPEN = 0
           CALL FILERD(DOCNAM,NLET,NULL,'DOCUMENT',IRTFLG)
           IF (IRTFLG .EQ. -1) RETURN
           CALL RDPRMI(NNSPH,IREG,NOT_USED,
     &          'NUMBER OF SPHERES, STARTING COL. FOR X,Y, & Z')

           CALL RDPRM3S(XOFF,YOFF,ZOFF,NOT_USED,
     &          'X, Y, & Z OFFSETS',IRTFLG)
           IF (IRTFLG .NE. 0) RETURN
	ENDIF


        GETOLD = .FALSE.
        IF (ANS(1:3) .EQ. 'SPV') THEN
C          VARIABLE SPHERE DENSITIES READ FROM DOC FILE
           IDEN = 1
           IF (ANS(4:4) .EQ. 'A') GETOLD = .TRUE.
        ELSE
C          ALL SPHERES WILL HAVE SAME DENSITIES
           IDEN = 0
           CALL RDPRM(AMQ, NOT_USED, 
     &        'ENTER MASS VALUE PER VOXEL (DEFAULT = 2.0)')
           IF (AMQ .EQ. 0.0) AMQ = 2.0     
c          DEFAULT VALUE 2.0 
        ENDIF
        NSPH = 0
C       THE DEFAULT COORDINATES OF THE ALL THE SPHERES ARE CENTER.

        DO  L = 1, MAXSPH
           RAD1(L)  = 0.0
           RAD2(L)  = 0.0
           XX(L) = NSAM/2 + 1
           YY(L) = NROW/2 + 1
           ZZ(L) = NSLICE/2 + 1
        ENDDO
		
 	CALL RDPRM2(RAD2(1),RAD1(1),
     &       NOT_USED,'OUTER AND INNER RADII OF SPHERES')

C       CHANGED 10/12/88 MR
C       FIND HOW MANY RADII HAVE TO BE READ FROM DOCUMENT FILE:
	NSPH = 1
        IRAD = 0
        IF (RAD2(1) .LT. 0) IRAD = 1
        IF (RAD1(1) .LT. 0) IRAD = 2

	IF (.NOT. LDOC) WRITE(NOUT,*) ' TYPE: 0, 0, 0  TO STOP INPUT'
2200	RAD2(NSPH)   = RAD2(1)
	RAD1(NSPH)   = RAD1(1)

C       'SPA' OPTION -- ADD MASSES WHERE SPHERES OVERLAP
C       'SP'  OPTION -- LET SPHERES INTERPENETRATE

	SWITCH = 0.0
	IF (ANS(3:3) .EQ. 'A') SWITCH = 1.0

	IF (LDOC) THEN
           CALL UNSAV(DOCNAM,NOPEN,NDOC,NSPH,CX,IREG+2+IRAD+IDEN,
     &                LERR,1)  
c          SEQUENTIAL READ SPECIFIED FOR DOC FILE
           NOPEN   = 1
           ADUM(1) = CX(IREG)
           ADUM(2) = CX(IREG+1)
           ADUM(3) = CX(IREG+2)
           IF (IRAD.GT.0) RAD2(NSPH) = CX(IREG+3)
           IF (IRAD.GT.1) RAD1(NSPH) = CX(IREG+4)
           IF (IDEN.GT.0) DENS(NSPH) = CX(IREG+3+IRAD)
           WRITE(NOUT,2400)NSPH,(ADUM(K),K=1,3)
	ELSE
	
           CALL RDPRM3S(ADUM(1),ADUM(2),ADUM(3),NOT_USED,
     &                 'X,Y,Z',IRTFLG)

C*         ALTERED SO DOES NOT ACCEPT LAST 0.0,0.0,0,0 al june 88
           IF (ADUM(1).EQ.0..AND.ADUM(2).EQ.0..AND.ADUM(3).EQ.0.)THEN
              NSPH = NSPH - 1
              IF (NSPH .EQ. 0) RETURN
              GOTO 3000
           ENDIF
           WRITE(NOUT,2400) NSPH,(ADUM(K),K=1,3)
 2400      FORMAT(I4,3F12.6)
	ENDIF

	XX(NSPH) = ADUM(1) + XOFF
	YY(NSPH) = ADUM(2) + YOFF
	ZZ(NSPH) = ADUM(3) + ZOFF
 
        IF (ZZ(NSPH)+RAD2(NSPH) .GT. NSLICE .OR.
     &      ZZ(NSPH)-RAD2(NSPH) .LT. 0      .OR.
     &      YY(NSPH)+RAD2(NSPH) .GT. NROW   .OR.
     &      YY(NSPH)-RAD2(NSPH) .LT. 0      .OR.
     &      XX(NSPH)+RAD2(NSPH) .GT. NSAM   .OR.
     &      XX(NSPH)-RAD2(NSPH) .LT. 0)      WRITE(NOUT,2450)NSPH

 2450   FORMAT(' *** SPHERE ',I3,' WILL BE TRUNCATED')
	IF (LDOC .AND. NSPH .EQ. NNSPH) GOTO 3000       
	IF (NSPH .EQ. MAXSPH) THEN
	   WRITE(NOUT,2500)MAXSPH
 2500	   FORMAT(' ** NUMBER OF SPHERES LIMITED TO',I3)
	   GOTO 3000
	ENDIF
	NSPH = NSPH+1
	GOTO 2200

C       SPHERES ARRANGED ON HELIX ******************************* HELIX

 2600   CONTINUE
        GETOLD = .FALSE.
        IDEN   = 0
        AMQ    = 2.0
        SWITCH = 0.0
        CALL RDPRM1S(AMQ, NOT_USED, 
     &      'ENTER MASS VALUE PER VOXEL (DEFAULT = 2.0)',IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

        RADT = 3
        RHEL = 3
	CALL RDPRM2S(RADT,RHEL,NOT_USED,
     &            'SPHERE RADIUS, HELIX RADIUS',IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

	CALL RDPRMI(NSPH,NTURN,NOT_USED,'NO. OF SPHERES, NO. OF TURNS')
	IF (NSPH .EQ. MAXSPH) THEN
	   CALL ERRT(102,'NUMBER OF SPHERES LIMITED TO',MAXSPH)
	   RETURN
	ENDIF

	DO N=1,NSPH
           RAD1(N) = 0.0
           RAD2(N) = RADT
           XX(N)   = NSAM/2+1    + RHEL * COS(TWOPI * FLOAT(N-1) * 
     &               FLOAT(NTURN) / FLOAT(NSPH))
           ZZ(N)   = NSLICE/2+1  + RHEL * SIN(TWOPI * FLOAT(N-1) *
     &               FLOAT(NTURN) / FLOAT(NSPH))
           YY(N)   = FLOAT(N-1)  * FLOAT(NROW) / FLOAT(NSPH)
        ENDDO

 3000 IF (NSPH .EQ. 0) NSPH = 1
      DO L1 = 1, NSLICE
         DO  N = 1, NSPH
           FL2(N) = (L1-ZZ(N))**2
         ENDDO

         DO L2 = 1, NROW
           DO N  =1, NSPH
              FL22(N) = FL2(N) + (L2-YY(N))**2
           ENDDO

           IF (GETOLD) THEN
C             WANT TO ADD NEW SPERES TO EXISTING FILE
              CALL REDLIN(LUN1,A0,NSAM,(L1-1)*NROW+L2)
           ENDIF

           DO L3 = 1, NSAM
             IF (.NOT. GETOLD) A0(L3) = 0.0

             DO N=1, NSPH
               RAD41 = RAD1(N)**2
               RAD44 = RAD2(N)**2
               IF (IDEN .GT. 0) AMQ = DENS(N)
               XYZ   = FL22(N) + (L3-XX(N))**2
               IF (XYZ .GE. RAD41 .AND.  XYZ .LE. RAD44) THEN
C                NAIK GENERAL AMQ NOT 2
                 A0(L3) = AMQ + A0(L3) * SWITCH    
               ENDIF
             ENDDO
           ENDDO
           CALL WRTLIN(LUN1,A0,NSAM,(L1-1)*NROW+L2)
         ENDDO
      ENDDO
      RETURN

C     RANDOM DISTRIBUTION ************************************* RANDOM

 7000 CALL RDPRMC(GA,NCHAR,.TRUE.,'GAUSSIAN DISTRIBUTION.? (Y/N)',
     &            NULL,IRTFLG)
      IF (IRTFLG .EQ. -1) RETURN
      
      GFLAG = (GA .EQ. 'Y')
      IF (GFLAG) THEN
	CALL RDPRM2(XM,SD,NOT_USED,
     &            'MEAN AND STANDARD DEVIATION OF GAUSSIAN DIST.')
      ENDIF

      DO L = 1,NSLICE
        DO I = 1,NROW
          IF (GFLAG) THEN
             DO K = 1,NSAM
                A0(K) = RANN(XM,SD)
             ENDDO
          ELSE
	     CALL  RANDOM_NUMBER(HARVEST=A0)
          ENDIF
          CALL WRTLIN(LUN1,A0,NSAM,I+(L-1)*NROW)
        ENDDO
      ENDDO

      RETURN

C     BLANK (CONSTANT DENSITY) ********************************** BLANK

 4100 CALL RDPRM(BACK,NOT_USED,'BACKGROUND CONSTANT')
      A0 = BACK
      DO L=1,NSLICE
         DO I=1,NROW
            CALL WRTLIN(LUN1,A0,NSAM,I+(L-1)*NROW)
         ENDDO
      ENDDO

      RETURN

C     CYLINDERS (INTERPENETRATING) ************************** CYLINDERS

C     CLEAR VOLUME FIRST (SET VOLUME BACKGROUND)

 4600 CALL RDPRM(BVAL,NOT_USED,'VALUE OF VOLUME BACKGROUND')
      A0 = BVAL
      DO  L=1,NSLICE
        I1=(L-1)*NROW
        DO  I=1,NROW
          CALL WRTLIN(LUN1,A0,NSAM,I1+I)
        ENDDO
      ENDDO

 5000 WRITE(NOUT,*) 'YOU HAVE A CHOICE OF 3 CYLINDER AXES.'
      CALL RDPRMC(YN,NA,.TRUE.,
     &            'CHOOSE X, Y, Z (OR Q TO END CYLINDER ENTRY)',
     &            NULL,IRTFLG)
      IF (YN .NE. 'X' .AND. YN .NE. 'Y' .AND. YN .NE. 'Z') RETURN 

      CALL RDPRM2(RADC,HT,NOT_USED,'RADIUS, HEIGHT')
      IF (RADC .EQ. 0.0 .OR.  HT   .EQ. 0.0) RETURN
      HT = HT / 2.0
      CALL RDPRM2(XH,YH,NOT_USED,'X & Y COORDINATES OF CENTER')
      CALL RDPRM2(ZH,VALUE,NOT_USED,'Z COORDINATE, DENSITY')

C     ****************** CYLINDER AXIS ALONG Y ***********
      IF (YN .EQ. 'Y') THEN 
         NSL1  = ZH-RADC-0.5
         NSL2  = ZH+RADC+0.5
         IF (NSL1.LT.1) NSL1=1
         IF (NSL2.GT.NSLICE) NSL2=NSLICE
         NX1   = XH-RADC-0.5
         NX2   = XH+RADC+0.5
         IF (NX1.LT.1) NX1=1
         IF (NX2.GT.NSAM) NX2=NSAM
         NY1   = YH-HT-0.5
         NY2   = YH+HT+0.5
         IF (NY1.LT.1) NY1=1
         IF (NY2.GT.NROW) NY2=NROW
         RAD44 = RADC**2
         DO L=NSL1,NSL2
           I1 = (L-1)*NROW
           Z2 =( FLOAT(L)-ZH)**2
           DO I=NY1,NY2
             CALL REDLIN(LUN1,A0,NSAM,I+I1)
             DO K=NX1,NX2
                IF ((FLOAT(K)-XH)**2+Z2.LE.RAD44) A0(K) = VALUE                 
             ENDDO
             CALL WRTLIN(LUN1,A0,NSAM,I1+I)
           ENDDO
         ENDDO
      ENDIF
 
C     ****************** CYLINDER AXIS ALONG Z ***********
      IF (YN .EQ. 'Z') THEN 

         NSL1  = ZH-HT-0.5
         NSL2  = ZH+HT+0.5
         IF (NSL1 .LT.1) NSL1=1
         IF (NSL2 .GT. NSLICE) NSL2=NSLICE
         NX1   = XH-RADC-0.5
         NX2   = XH+RADC+0.5
         IF (NX1 .LT.1) NX1=1
         IF (NX2 .GT. NSAM) NX2=NSAM
         NY1   = YH-RADC-0.5
         NY2   = YH+RADC+0.5
         IF( NY1 .LT.1) NY1=1
         IF (NY2 .GT. NROW) NY2=NROW
         RAD44 = RADC**2

         DO L=NSL1,NSL2
           I1 = (L-1)*NROW
           DO I=NY1,NY2
             CALL REDLIN(LUN1,A0,NSAM,I+I1)
             Y2 = (FLOAT(I)-YH)**2
             DO  K=NX1,NX2
                IF((FLOAT(K)-XH)**2+Y2.LE.RAD44) A0(K) = VALUE 
             ENDDO                
             CALL WRTLIN(LUN1,A0,NSAM,I1+I)
           ENDDO
         ENDDO
      ENDIF

C     ****************** CYLINDER AXIS ALONG X ***********
      IF (YN .EQ. 'X') THEN 
         NSL1  = ZH-RADC-0.5
         NSL2  = ZH+RADC+0.5
         IF (NSL1.LT.1) NSL1=1
         IF (NSL2.GT.NSLICE) NSL2=NSLICE
         NX1   = XH-HT-0.5
         NX2   = XH+HT+0.5
         IF (NX1.LT.1) NX1=1
         IF (NX2.GT.NSAM) NX2=NSAM
         NY1   = YH-RADC-0.5
         NY2   = YH+RADC+0.5
         IF (NY1.LT.1) NY1=1
         IF (NY2.GT.NROW) NY2=NROW
         RAD44 = RADC**2
         DO L=NSL1,NSL2
            I1=(L-1)*NROW
            Z2=(FLOAT(L)-ZH)**2
            DO I=NY1,NY2
               Y2=(FLOAT(I)-YH)**2
               IF (Z2+Y2 .LE. RAD44) THEN
                  CALL REDLIN(LUN1,A0,NSAM,I+I1)
                  DO  K=NX1,NX2
                     A0(K) = VALUE                 
                  ENDDO
                CALL WRTLIN(LUN1,A0,NSAM,I1+I)
              ENDIF
           ENDDO
         ENDDO
      ENDIF
      GOTO 5000

C     PUT NUMBERS INTO VOLUME:******************************** NUMBERS

 6000 CALL RDPRMI(IPOS,IDUM,NOT_USED,'POSITION IN LINE')
      CALL RDPRM(SCALE,NOT_USED,'SCALING FACTOR')
      SUM = 0.0
      DO L=1,NSLICE
        DO I=1,NROW
          DO K=1,NSAM
            A0(K) = 0.0
          ENDDO
          A0(IPOS) = FLOAT((L-1)*NROW+I)*SCALE
          SUM      = SUM+A0(IPOS)
          CALL WRTLIN(LUN1,A0,NSAM,I+(L-1)*NROW)
        ENDDO
      ENDDO

      CALL RDPRMC(ANS,NCHAR,.TRUE.,'NUMBER PIXELS IN LINE? (Y/N)',
     &              NULL,IRT)
      IF (ANS(1:1) .EQ. 'Y') THEN
        CALL RDPRMI(ISLIC,ILINE,NOT_USED,'SLICE NUMBER, LINE NUMBER')
        CALL RDPRM(SCALE2,NOT_USED,'SCALING FACTOR')
        ILL = (ISLIC-1) * NROW + ILINE
        CALL REDLIN(LUN1,A0,NSAM,ILL)
        DO  I=1,NSAM
           A0(I) = FLOAT(I) * SCALE2
           SUM   = SUM + A0(I)
        ENDDO
        SUM = SUM - FLOAT(IPOS)*SCALE
        CALL WRTLIN(LUN1,A0,NSAM,ILL)
      ENDIF
      RETURN

C     ******************************************* 3D GAUSSIAN FUNCTION

8000	CALL RDPRA('COORDINATES OF THE CENTER X,Y,Z',
     &         3,0,.FALSE.,RAD1,NVAL,IRTFLG)

	CALL RDPRA('STD. DEVIATIONS X,Y,Z',
     &         3,0,.FALSE.,RAD2,NVAL,IRTFLG)

	GNM = 1.0/RAD2(1)/RAD2(2)/RAD2(3)/(2*QUADPI)**1.5
	TNM = ALOG(1.0/TINY(GNM))

	DO I=1,3
	   RAD2(I) = RAD2(I) * RAD2(I)
	ENDDO

        DO K = 1,NSLICE
           DO J = 1,NROW
              DO I = 1,NSAM
	        EEE=0.5* ((I-RAD1(1))**2/RAD2(1)+
     &		    (J-RAD1(2))**2/RAD2(2)+
     &		    (K-RAD1(3))**2/RAD2(3))
	        IF (EEE.GE.TNM) THEN
	          A0(I) = 0.0
	        ELSE
                  A0(I)=GNM*EXP(-EEE)
	        ENDIF
             ENDDO
             CALL WRTLIN(LUN1,A0,NSAM,J+(K-1)*NROW)
           ENDDO
        ENDDO
	RETURN


 5500 WRITE(NOUT,*) '*** UNIDENTIFIED OPTION'
      CALL ERRT(100,'MODEL3',NE)

      RETURN
      END