
C++*********************************************************************
C
C MODEL.FOR                                      REWRITTEN APRIL 97 al
C           RDPRAF REMOVED                       DEC 05 ARDEAN LEITH 
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   MODEL(LUN1,NSAM,NROW)
C
C   PURPOSE:  PREPARES TEST PICTURES
C
C   PARAMETERS:
C        LUN1        LOGICAL UNIT NUMBER OF FILE
C        NSAM,NROW   DIMENSIONS OF FILE
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************

	SUBROUTINE MODEL(LUN1,NSAM,NROW)
#ifdef SP_NT
	use dfport
#endif
	INCLUDE 'CMBLOCK.INC'

	COMMON NX(20),NY(20),CX(20)

	INTEGER    IRTFLG,NVAL
	DIMENSION  PH(20),FWA(4)
	DIMENSION  A0(NSAM)
	PARAMETER  (QUADPI = 3.1415926535897932)
	PARAMETER  (TWOPI=2*QUADPI)

	CHARACTER    ANS,NULL,GA

	DATA  MAXSIN/20/

	NULL = CHAR(0)

	WRITE(NOUT,100)
 100  FORMAT(
     &      ' .MENU: (B)LANK    -- CONSTANT DENSITY IMAGE'/
     &      '        (C)IRCLE   -- FILLED CIRCLE'/
     &      '        (G)AUSSIAN -- GAUSSIAN FUNCTION'/
     &      '        (R)ANDOM   -- RANDOM STATISTICS'/
     &      '        (S)INE     -- SET OF SINE WAVES'/
     &      '        (T)EST     -- 2D SINE WAVE'/
     &      '        (W)EDGE    -- DENSITY WEDGE'/)

1010  CALL RDPRMC(ANS,NC,.TRUE.,'B/C/G/R/S/T/W',NULL,IRTFLG)
      IF (IRTFLG .NE. 0) RETURN

      IF (ANS .EQ. 'B') THEN
C       BLANK -------------------------------------------------- BLANK
        CALL RDPRM(BACK,NOT_USED,'BACKGROUND CONSTANT')

        CALL BLANK(LUN1,NSAM,NROW,1,BACK)

      ELSEIF (ANS .EQ. 'C' .OR. ANS .EQ. 'W') THEN
C       CIRCLE OR WEDGE ------------------------------- CIRCLE OR WEDGE
        IF (ANS .EQ. 'C') THEN
           CALL RDPRM(RAD,NOT_USED,'RADIUS (FLOATING POINT)')
           RAD2 = RAD**2
        ENDIF

        SCF = 2.0/FLOAT(NROW+NSAM)
	DO I=1,NROW
          IF (ANS .EQ. 'C') THEN
C           CIRCLE WANTED
            FI2 = (FLOAT(I-NROW/2-1))**2
            DO  J=1,NSAM
               A0(J) = 0.0
               IF (FI2+FLOAT(J-NSAM/2-1)**2.LT.RAD2) A0(J) = 1.0
            ENDDO 
          ELSE
C           WEDGE WANTED
            FI = FLOAT(I)*SCF
            DO  J = 1,NSAM
               A0(J) = FI+ SCF * FLOAT(J)
            ENDDO
          ENDIF
C         SEND RECORD TO FILE
          CALL WRTLIN(LUN1,A0,NSAM,I)
        ENDDO

      ELSEIF (ANS .EQ. 'G') THEN
C       GAUSSIAN   ------------------------------------------ GAUSSIAN
	CALL RDPRM2(DX,DY,NOT_USED,'COORDINATES OF THE CENTER')
	CALL RDPRM2(SX,SY,NOT_USED,'STD. DEVIATIONS X,Y')
	IF(SX.LE.0.0.OR.SY.LE.0.0)  THEN
	  CALL ERRT(31,'MO',NE)
	  RETURN
	ENDIF
	GNM=1.0/SX/SY/2/QUADPI
	TNM=ALOG(1.0/TINY(GNM))
	SX=SX*SX
	SY=SY*SY
        DO I = 1,NROW
          DO K = 1,NSAM
	       EEE=0.5*((K-DX)**2/SX+(I-DY)**2/SY)
	       IF (EEE.GE.TNM)  THEN
	          A0(K)=0.0
	       ELSE
                  A0(K)=GNM*EXP(-EEE)
	       ENDIF
          ENDDO
          CALL WRTLIN(LUN1,A0,NSAM,I)
        ENDDO
C
      ELSEIF (ANS .EQ. 'R') THEN
C       PUT RANDOM NUMBERS IN THE IMAGE ------------------------- 
        CALL RDPRMC(GA,NC,.TRUE.,'GAUSSIAN DISTRIBUTION? (Y/N)',
     &             NULL,IRT)
        IF (IRT .NE. 0) GOTO 1010

        IF (GA .EQ. 'Y') THEN
          CALL  RDPRM2(XM,SD,NOT_USED,
     &        'MEAN AND STANDARD DEVIATION OF GAUSSIAN DIST.')
        ENDIF
        DO I = 1,NROW
         IF (GA.EQ.'Y') THEN
          DO  K = 1,NSAM
               A0(K) = RANN(XM,SD)
          ENDDO
         ELSE
	   CALL  RANDOM_NUMBER(HARVEST=A0)
         ENDIF
          CALL WRTLIN(LUN1,A0,NSAM,I)
        ENDDO

      ELSEIF (ANS .EQ. 'S' .OR. ANS .EQ. 'T') THEN

C       PUT SINE WAVES OF INTENSITY IN THE IMAGE ----------------- SINE
        IF (ANS .EQ. 'T') THEN
C          ONLY ONE WAVE IS WANTED
           NS    = 1
           CX(1) = 1.
           NX(1) = 2
           NY(1) = 2
           PH(1) = 0.0

        ELSEIF (ANS .EQ. 'S') THEN

C         FIND NUMBER OF SINE WAVES WANTED IN THE IMAGE ---------------

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

          DO I = 1,NS
9006        IRTFLG = 0
            CALL RDPRA('AMPLITUDE, PHASE, SP. FREQUENCY(KX,KY)',
     &         4,0,.FALSE.,FWA,NVAL,IRTFLG)
            IF (IRTFLG .NE. 0) GOTO 9006

            CX(I) = FWA(1)
            PH(I) = FWA(2)
            ANX   = FWA(3)
            ANY   = FWA(4)
	    
            NX(I) = ANX+0.5
            IF (ANX .LT .0.0) NX(I) = ANX-0.5
            NY(I) = ANY+0.5
            IF (ANY .LT. 0.0) NY(I) = ANY-0.5
            WRITE(NOUT,9008)I,CX(I),PH(I),NX(I),NY(I)
9008        FORMAT(1X,I5,2F8.2,2I6)
            PH(I) = PH(I)*TWOPI/360.0
          ENDDO
        ENDIF

C       PLACE SINE WAVE(S) OF INTENSITY IN THE IMAGE
        DO I=1,NROW
           PHASE = FLOAT(I-1)*TWOPI/FLOAT(NROW)
           DO J=1,NSAM
              A0(J) = 0.0
              DO K = 1,NS
              A0(J)=CX(K)*SIN(FLOAT(J-1)*TWOPI*FLOAT(NX(K))/FLOAT(NSAM)
     1              +PHASE*FLOAT(NY(K))+PH(K))+A0(J)
              ENDDO
           ENDDO
           CALL WRTLIN(LUN1,A0,NSAM,I)
         ENDDO
      ENDIF

      END
