C     DEMO2
C
C     DEMONSTRATION MODEL FOR THE MICHIGAN EXPERIMENTAL SIMULATION
C     SUPERVISOR VERSION 3
C     ROBERT L. STOUT, SIMULATION DESIGNER AND PROGRAMMER
C     SEPTEMBER, 1973
C
      SUBROUTINE MODEL(N,NC1,*)
C
C
      LOGICAL*4 FLGS(12), VPT(6), XINC, XLAZY, MISS, INITLS
      LOGICAL*4 NOTIME, NOQUAL, TBLS
      INTEGER*4 SI(12), SF(12), CONDS(16), TMM, TMF, QM
      INTEGER*4 IV(12), TASK, WORKCD, HLDRS, ODEV1
      REAL*4 FV(12), DMISS(6), DV(6), INCOMP, INCTST, LTEST, INCME
      REAL*4 INC, LTDIST(4), LAZME, PMD(3), LCONST, QUALIN(5,2)
      REAL*4 FTIME(5,3), LAMBDA, FLOAF(2), FERR(2), TC(4), INCLEV(5,2)
      COMMON /VARS1/IV,FV,SI,SF,FLGS
      COMMON /CNDTNS/NCD,CONDS
      COMMON /IO/IDEV1,IDEV2, IDEV3, IDEV4, ODEV1
C
!!      DATA DMISS/6*-1./, VPT/6*.FALSE./
C
C     GIVE IDENTIFIABLE NAMES TO THE INDEPENDENT VARIABLES
      EQUIVALENCE (INCOMP, FV(1)), (LAZY, IV(1)), (WORKCD, IV(2)),
     1 (TASK, IV(3)), (XINC, FLGS(1)), (XLAZY, FLGS(2)), (NOTIME,
     2 FLGS(3)), (NOQUAL, FLGS(4))
C
C     DEPENDENT AND CONCOMITANT VARIABLES GO IN ARRAY DV
      EQUIVALENCE (INCTST, DV(1)), (LTEST, DV(2)), (TIME, DV(3)),
     1 (QUAL, DV(4))
C
C     THE FOLLOWING MODEL PARAMETERS ARE CONTROLLED BY FLOATING POINT
C     VARIABLES WHICH ARE USUALLY DECLARED INT
      EQUIVALENCE (INCME, FV(2)), (LAZME, FV(3)), (C1, FV(4)),
     1 (LCONST, FV(5)), (FSD, FV(6)), (QSD, FV(7)), (T1, FV(8)),
     2 (T2, FV(9)), (T3, FV(10)), (T4, FV(11)), (TBLS, FLGS(12))
C
C     TABLE PMD CONTROLS THE EFFECT OF WORKING CONDITIONS ON THE PROBA-
C     BILITY OF MISSING DATA
!!      DATA PMD/0.02, 0.08, 0.15/
C
C     PINIT = PROB. THAT A GIVEN NAPKIN HOLDER WILL HAVE TO HAVE
C     INITIALS STAMPED ON IT.  PREJ = PROB. THAT A GIVEN HOLDER WILL
C     BE FOUND DEFECTIVE ON FINAL INSPECTION.
!!      DATA PINIT/ .4/, PREJ/ .02/
C
C     COST PARAMETERS
      INTEGER COST5
!!      DATA COST1M/5./, COST1F/2./, COST2/2.5/, COST3/10./, COST4/.10/,
!!     1 COST5/2/, COST6/5./
C
C	MODEL NAME
	DOUBLE PRECISION MODNM1,MODNM2,MODNM3
	COMMON /DATANM/ MODNM1,MODNM2,MODNM3
	DATA MODNM1/'DEMO.DAT'/,MODNM2/'DEMO3.DAT'/,MODNM3/'DEMO.BIN'/
      DATA DMISS/6*-1./, VPT/6*.FALSE./
      DATA PMD/0.02, 0.08, 0.15/
      DATA PINIT/ .4/, PREJ/ .02/
      DATA COST1M/5./, COST1F/2./, COST2/2.5/, COST3/10./, COST4/.10/,
     1 COST5/2/, COST6/5./
C********************
C
C     SET UP OUTPUT FLAGS FOR FPOUT
      VPT(1)=.FALSE.
      VPT(2)=.FALSE.
      VPT(3)=.FALSE.
      VPT(4)=.FALSE.
      DO 1 I=1,NC1
      J=CONDS(I)
      CALL CVALS(J)
C
C     INCOMPETENCE TEST SCORES ARE PRINTED ONLY IF 'XINC' IS SPECIFIED
      VPT(1)=VPT(1) .OR. XINC
C     SAME FOR LAZINESS
      VPT(2)=VPT(2) .OR. XLAZY
C
C     TIME IS TO BE PRINTED UNLESS 'NOTIME' IS SPECIFIED
      VPT(3)=VPT(3) .OR. .NOT.NOTIME
C     SAME FOR QUALITY
 1    VPT(4)=VPT(4) .OR. .NOT.NOQUAL
      CALL FPOUT1(VPT, DMISS, &99)
C
      TMM=0
      TMF=0
      QM=0
C
C********************
C
C     SIMULATE THE REQUISITE NUMBER OF SUBJECTS
      DO 10 I=1,N
      K=CONDS(1)
      CALL CVALS(K)
C
C--------------------
C
C     SELECT INDIVIDUAL DIFFERENCES PARAMETERS AS NECESSARY
C
C     INCOMPETENCE:
      INCTST=INCOMP
C
C     EVEN IF INC IS CONTROLLED, IT IS NOT PERFECTLY CONTROLLED
      CALL URAND1( -.5, .5, X)
      INCTST=INCTST+X
      IF(INCTST .LT. 0.)INCTST=0.
      IF(INCTST .GT. 10.)INCTST=10.
C
C     IF INCOMPETENCE IS RANDOM, SELECT A TEST SCORE FROM AN UNIFORM
C     DISTRIBUTION
      IF(SF(1) .NE. 0)CALL URAND1(0., 10., INCTST)
C
C     ROUND INC SCORE TO TWO DECIMAL PLACES
      CALL ROUND(INCTST, .01)
C
C     ADD MEASUREMENT ERROR TO INCOMPETENCE TEST SCORE TO DETERMINE
C     'TRUE' INCOMPETENCE LEVEL
      CALL NRAND(INCTST, INCME, INC)
C
C     RESCALE AND TRUNCATE
      IF(INC .GT. 10.)INC=10.
      INC=.05*INC+.0025*INC*INC
      IF(INC .LT. 0.02)INC=0.02
C
C
C     LAZINESS:
C     TAKE CARE OF NUMBER-KEYWORD EQUIVALENCE
      LAZY=LAZY+SI(1)
      LTEST=LAZY
      IF(LAZY .LT. 5)GO TO 12
C
C     LAZINESS IS RANDOM--DECIDE TEST SCORE ACCORDING TO THE MULTINOMIAL
C     DISTRIBUTION SPECIFIED IN LTDIST.
      CALL MULNOM(LTDIST, 4, 1, 1, 1, K, J, J, J)
      LTEST=K
C
C     GENERATE REAL LAZINESS PARAMETER GIVEN TEST SCORE
C     MEAN UNDERLYING LAZINESS PARAMETER IS LINEAR WITH TEST SCORE
C     EXCEPT FOR THE HIGHEST LEVEL OF LAZINESS
 12   X=LTEST
      IF(X .GE. 3.99)X=8.
C
C     ADD RANDOM ERROR
      CALL NRAND(X, LAZME, QLAZ)
C
C     RESCALE AND TRUNCATE
      QLAZ=0.1*QLAZ
      IF(QLAZ.LT.0.1)QLAZ=0.1
      IF(QLAZ.GT.0.9)QLAZ=0.9
C
C
      MISS=.FALSE.
      FTLAST=0.
C
C********************
C
C     SIMULATE CONDITIONS WITHIN SUBJECTS
      DO 20 J=1,NC1
C
C     SET UP CURRENT VALUES FOR ALL IVS
      K=CONDS(J)
      CALL CVALS(K)
C
C     SET UP VARIABLE PRINT INDICATORS FOR FPOUT
      NOTIME=.NOT.NOTIME
      NOQUAL=.NOT.NOQUAL
C
C     DO NOTHING IF SUBJECT HAS QUIT OR BEEN FIRED ON A PREVIOUS
C     TRIAL
      TIME=-1.
      QUAL=-1.
      IF(MISS)GO TO 19
C
C     WHAT TO DO NEXT DEPENDS ON WHAT TASK IS TO BE DONE
      IF(TASK.EQ.2)GO TO 31
C
C--------------------
C
C     THE TASK IS MANUFACTURING NAPKIN HOLDERS
C     MUST FIRST DETERMINE WHETHER SUBJECT QUITS OR IS FIRED
C     ON THIS TRIAL
C
C     THE PROBABILITY OF HAVING MISSING DATA ON THIS TRIAL DEPENDS
C     ON WORKING CONDITIONS PLUS THE CHANCE THAT THE SUBJECT IS FIRED
C     FOR INCOMPETENCE AND/OR LAZINESS
C
      P=PMD(WORKCD)+C1*INC*QLAZ
      CALL BINOM(1,P,K)
      MISS=K .GE. 1
      IF(MISS)GO TO 19
C
C     KEEP TRACK OF NO. OF TIMES TIME AND QUALITY ARE MEASURED
C     FOR USE IN COST CALCULATIONS
C     NOTE:NOTIME AND NOQUAL WERE INVERTED ABOVE
      IF(NOTIME)TMM=TMM+1
      IF(NOQUAL)QM=QM+1
C
C--------------------
C
C     SIMULATE THE PRODUCTION OF 10 PORNOGRAPHIC NAPKIN HOLDERS
      HLDRS=0
      TIME=0.
C
C
C     DECIDE IF THE NEXT NAPKIN HOLDER MUST BE 'PERSONALIZED' BY
C     HAVING INITIALS STAMPED ON IT
 21   CALL BINOM(1, PINIT, K)
      INITLS=K.GT.0
C
C
C     INITIAL CUTTING AND STAMPING OF PARTS REQUIRES A CONSTANT MIN-
C     IMUM OF TIME PLUS A RANDOM PART DISTRIBUTED AS A GAMMA
C     RANDOM VARIABLE WITH SHAPE PARAMETER ALPHA=2, AND SCALE
C     PARAMETER BETA=TC(1), PLUS OPPORTUNITY FOR LOAFING
 22   TIME=TIME+T1
C
      CALL RGAMMA(2, TC(1), X)
      TIME=TIME+X
C
C     ADD TIME FOR LOAFING AS APPROPRIATE
      CALL LOAF(QLAZ, LCONST, TIME)
C
C     PROBABILITY OF HAVING TO REDO THIS STEP DEPENDS ON INCOMPETENCE
      CALL CHANCE(INC, &22)
C
C
C     DO INITIALS HAVE TO BE STAMPED ON THIS HOLDER?
      IF(.NOT.INITLS)GO TO 23
C
C     YES--TIME REQUIRED IS A CONSTANT MINIMUM PLUS A GAMMA RV
C     WITH ALPHA=1, BETA=TC(2), PLUS A SHORT OPPORTUNITY FOR LOAFING
      TIME=TIME+T2
      CALL RGAMMA(1, TC(2), X)
      TIME=TIME+X
      X=LCONST*0.5
      CALL LOAF( QLAZ, X, TIME )
C
C     A TRULY INCOMPETENT WORKER CAN RUIN THE PRODUCT AT THIS STEP
      X=.25*INC
      CALL CHANCE(X, &22)
C
C
C     ASSEMBLING PARTS REQUIRES A CONSTANT AMOUNT OF TIME PLUS A
C     GAMMA RV WITH ALPHA=1 AND BETA=TC(3), PLUS ALLOWANCE FOR
C     LOAFING
 23   TIME=TIME+T3
      CALL RGAMMA(1, TC(3), X)
      TIME=TIME+X
      CALL LOAF( QLAZ, LCONST, TIME)
C
C     THERE IS A SMALL CHANCE DEPENDING SLIGHTLY ON INCOMPETENCE THAT
C     THE ITEM DOES NOT PASS INSPECTION AT THIS POINT
      P=PREJ+.05*INC
      CALL CHANCE(P, &22)
C
C
C     TIME REQUIRED FOR PACKAGING IS A CONSTANT PLUS A GAMMA RV
C     WITH ALPHA=1 AND BETA=TC(4), PLUS OPPORTUNITY FOR LOAFING
      TIME=TIME+T4
      CALL RGAMMA(1, TC(4), X)
      TIME=TIME+X
      CALL LOAF(QLAZ, LCONST, TIME)
C
C
C     ONE MORE PORNOGRAPHIC NAPKIN HOLDER IS AVAILABLE TO THE WORLD
      HLDRS=HLDRS+1
      IF(HLDRS .LT. 10)GO TO 21
C
C--------------------
C
C     CALCULATE QUALITY INDEX:
C     AVERAGE QUALITY DEPENDS ONLY ON INCOMPETENCE--CALCULATE MEAN BY
C     INTERPOLATION
      CALL PLFN(QUALIN, 5, INC, X)
C
C     ADD RANDOM ERROR, CHOP, AND ROUND TO ONE DECIMAL PLACE
      CALL NRAND(X, QSD, QUAL)
      IF(QUAL .LT. 1.)QUAL=1.
      IF(QUAL .GT. 3.)QUAL=3.
C
C     TASK COMPLETE
      GO TO 50
C
C--------------------
C
C     FILLING-OUT-FORMS TASK
C     MEAN TIME TO FILL OUT FORMS INVOLVES AN INTERACTION OF INCOM-
C     PETENCE AND WORKING CONDITIONS AS DESCRIBED IN TABLE FTIME,
C     PLUS NUMEROUS OPPORTUNITIES FOR LOAFING, PLUS AN INCREASE IN
C     TIME ACROSS TRIALS.
C
C     RESCALE INC FOR INTPOL
 31   CALL PLFN(INCLEV, 5, INC, XI)
      XW=WORKCD
      CALL INTPOL(FTIME, 5, 3, 1, 1, XI, XW, 1., 1., X)
C
C     ADD NORMALLY DISTRIBUTED ERROR
      CALL NRAND(X, FSD, TIME)
      IF(TIME .LT. 5.)TIME=5.
C
C
C     FILLING OUT FORMS ALLOWS MANY OPPORTUNITIES FOR LOAFING--
C     NUMBER OF LOAFING INCIDENTS IS DETERMINED BY A POISSON PROCESS
C     WITH LAMBDA DEPENDENT ON LAZINESS
C
      LAMBDA=FLOAF(1)+FLOAF(2)*QLAZ
      CALL POISSN(LAMBDA, K)
C
C     SIMULATE THE REQUIRED NUMBER OF LOAFING EPISODES
      IF(K .LE. 0)GO TO 34
      DO 35 L=1,K
 35   CALL LOAF(1., LCONST, TIME)
C
C     TO SIMULATE CUMULATIVE EFFECT OF DISGUST, ADD 25% OF TIME ON LAST
C     TRIAL TO TIME FOR THIS TRIAL
      TIME=TIME+.25*FTLAST
      FTLAST=TIME
C
C
C     'QUALITY' IN THIS CASE MEANS THE NUMBER OF ERRORS MADE IN
C     FILLING OUT THE FORMS--ANOTHER POISSON RV, THIS ONE CON-
C     TROLLED BY INCOMPETENCE
 34   LAMBDA=FERR(1)+FERR(2)*INC
      CALL POISSN(LAMBDA, K)
      QUAL=K
C
C     KEEP TRACK OF NO. OF TIMES TIME AND QUALITY ARE MEASURED
C     FOR USE IN COST CALCULATIONS
C     NOTE:NOTIME AND NOQUAL WERE INVERTED ABOVE
      IF(NOTIME)TMF=TMF+1
      IF(NOQUAL)QM=QM+1
C
C--------------------
C
C     END OF CONDITIONS LOOP
C     ROUND TIME AND QUALITY AND PRINT
 50   CALL ROUND(TIME, .25)
      CALL ROUND(QUAL, .1)
19	ITMP=I
	JTMP=J
      CALL FPOUT2(ITMP, JTMP, FLGS, DV)
20	CONTINUE
C     END OF SUBJECTS LOOP
 10   CONTINUE
C
C********************
C
C     COST CALCULATION:
C     COST FOR MEASURING TIME IN THE MFG TASK IS COST1M DOLLARS
C     PER SUBJECT PER CONDITION, WITH ALLOWANCE FOR MISSING DATA.
C     IN THE FORMS TASK, THE COST IS COST1F DOLLARS PER SUBJECT
C     PER CONDITION.  QUALITY MEASUREMENT COSTS COST2 DOLLARS PER
C     SUBJECT PER CONDITION IN EACH TASK.  THE COST OF SETUP AND
C     DATA ANALYSIS IS COST3 + COST4*NO. OF SS.  MEASURING INCOM-
C     PETENCE AND LAZINESS COSTS NOTHING.  A LITTLE RANDOM VARIATION
C     IS ADDED TO THE FINAL COST IN THE FORM OF A GAMMA RV WITH
C     ALPHA=COST5 AND BETA=COST6 + .1*AMOUNT OF TIME MEASUREMENT.
      X=TMM
      COST=X*COST1M
      X=TMF
      COST=COST+X*COST1F
      X=QM
      COST=COST+X*COST2
      X=N
      COST=COST+COST3+X*COST4
      Y=TMM+TMF
      Y=.1*Y+COST6
      CALL RGAMMA(COST5,Y,X)
      COST=COST+X
C
C     ROUND OFF TO DOLLARS AND CENTS
      CALL ROUND(COST, .01)
      CALL FPOUT3(COST)
C
      IF(.NOT.TBLS)RETURN
C
C     PRINT TABLES FOR DEBUGGING PURPOSES
      WRITE(ODEV1,1002)LTDIST
      WRITE(ODEV1,1002)TC
      WRITE(ODEV1,1002)QUALIN
      WRITE(ODEV1,1002)INCLEV
      WRITE(ODEV1,1002)FTIME
      WRITE(ODEV1,1002)FLOAF
      WRITE(ODEV1,1002)FERR
      RETURN
C
C********************
C
C     MODEL INITIALIZATION:
      ENTRY MINIT
      OPEN(UNIT=IDEV3,ACCESS='SEQIN',FILE=MODNM2,DIRECTORY='4030,11')
C  
      READ(IDEV3,1001)LTDIST
      READ(IDEV3,1001)TC
      READ(IDEV3,1001)QUALIN
      READ(IDEV3,1001)INCLEV
      READ(IDEV3,1001)FTIME
      READ(IDEV3,1001)FLOAF
      READ(IDEV3,1001)FERR
      RETURN
C
C
 1001 FORMAT(6X,5G5.3)
 1002 FORMAT(1H0/(1X,5G11.3))
C
C     ERROR DETECTED BY FPOUT--PROBABLY NO VARIABLES
C     TO BE PRINTED (NOTIME AND NOQUAL SPECIFIED, BUT NOT XINC
C     OR XLAZY)
 99   RETURN 1
      END
      SUBROUTINE LOAF(QLAZ, LCONST, TIME)
      REAL*4 LCONST
C
C     THIS SUBROUTINE DECIDES IF A LOAFING INCIDENT IS TO OCCUR
C     AND, IF SO, SIMULATES THE TIME IT TAKES BY GENERATING
C     A GAMMA RANDOM VARIABLE WITH SHAPE PARAMETER ALPHA=2 AND
C     SCALE PARAMETER BETA=LCONST, AND ADDING IT TO TIME
C
C
      CALL CHANCE(QLAZ, &1)
      RETURN
C
C     A LOAFING INCIDENT IS TO OCCUR--CALCULATE TIME
 1    CALL RGAMMA(2, LCONST, X)
      TIME=TIME+X
      RETURN
      END

