	PROGRAM FIVMET
	REAL*4 OILDENSITY, METALDENSITY
	COMMON/HISTO/XPOS(5,500),DIAMETER,VOLTS,V0Z,SIGVX,SIGVZ
	COMMON/RGA/IT,F1,C,R
	INCLUDE GRAPH.INC
	INCLUDE EXEC.INC
	INCLUDE INTUIT.INC
	INCLUDE GCBIG.PRM
	DATA SPACING/1.4/,ALENGTH/500./
     1 	     ,PI/3.14159265/,METALDENSITY/6.1/
     2       ,G/981./,Q/1.602E-19/,OILDENSITY/0.91/
	SIGVX = 0.01
	SIGVZ = 1.0
	IT = 1
c	DENSITY = METALDENSITY
c	WRITE(9,10)
c10	FORMAT('ENTER DROP DIAMTER IN MICRONS')
c	ACCEPT DIAMETER
c	WRITE(9,15)
c15	FORMAT('ENTER DROP VELOCITY IN CM/SEC')
c	ACCEPT V0Z
c	WRITE(9,20)
c20	FORMAT('ENTER PLATE POTENTIAL IN KILOVOLTS')
c	ACCEPT VOLTS
c	WRITE(9,25)
c25	FORMAT('OIL (0) OR METAL (1)')
c	ACCEPT DENSITY
c	IF(DENSITY.EQ.0) DENSITY = OILDENSITY
c	WRITE(9,30)
c30	FORMAT('NUMBER OF ELECTRONS?')
c	ACCEPT ELECTRONS
	DENSITY = OILDENSITY
	DIAMETER = 90.
	V0Z = 1000.
	VOLTS = 30.
	ELECTRONS = 0.
	VOLTS = VOLTS * 1000.
	RADIUS = DIAMETER/2.E4
	DO 40 J = 1,5
	DO 40 I = 1,500
	XPOS(J,I) = 0.
40	CONTINUE
	VDROP = 4*PI*(RADIUS)**3/3.
	DROPMASS = VDROP*DENSITY
	DO 50 I = 1,5
	DO 50 J = 1,500
	CALL RGAUSS(0.,SIGVX,VX)
	CALL RGAUSS(V0Z,SIGVZ,VZ)
	T = (-VZ+SQRT(VZ**2+2.*G*ALENGTH))/G
	AX = Q*(ELECTRON+FLOAT(3-I))*VOLTS*1.E7/(DROPMASS*SPACING)
	XPOS(I,J) = 0.5*AX*T**2
	XPOS(I,J) = XPOS(I,J) + VX*T
50	CONTINUE
	V0Z = V0Z/100.
	CALL HIST
	END

	SUBROUTINE RGAUSS(X0,SIGMA,V)
C
	COMMON/RGA/IT,F1,C,R
	DATA PI/3.14159265/
11	IF(IT.NE.1) GO TO 2
	IT = 2
	R = RANFL(0)
	C = COS(2.0*PI*R)
	F1 = SQRT(-2.*ALOG(RANFL(0)))
	V = F1*C*SIGMA+X0
	GO TO 3
2	IT = 1
	V = F1*SIN(2.*PI*R)*SIGMA+X0
3	RETURN
	END

	SUBROUTINE HIST

	CALL DEVSEL(2,4,IERR)
	IF(IERR.NE.0) STOP 'UNABLE TO SELECT DEVICE 1'
	CALL BAR('X DEFLECTION (CM)'//CHAR(0),'NUMBER'//CHAR(0)
     1		,'DEFLECTION HISTOGRAM'//CHAR(0),0)
	CALL ENDPLT
	CALL RLSDEV
	RETURN
	END

	SUBROUTINE BAR(SXLAB,SYLAB,STITLE,TYPE)

C
C
	COMMON/HISTO/X(5,500),DIAMETER,VOLTS,V0Z,SIGVX,SIGVZ
	INTEGER TYPE,COUNT(6,512)
        REAL*4    XLOW,XHIGH,STEP,FBAR,YLOW,YHIGH,X0,Y0,VX0,VY0,VX1,VY1
        CHARACTER*1 SXLAB(20),SYLAB(20),STITLE(20)
	CHARACTER*20 LABEL1,LABEL2,LABEL3,LABEL4,LABEL5
C
C
        YLOW   = 0.0
        YHIGH  = 1.0
        FBAR   = 512.
C
        XLOW = X(1,1)
	XHIGH = X(1,1)
C
	DO 20 I = 1,5
	DO 10 J = 1,500
	XLOW = AMIN1(XLOW,X(I,J))
	XHIGH = AMAX1(XHIGH,X(I,J))
10	CONTINUE
20	CONTINUE
        STEP   = (XHIGH - XLOW) / FBAR
C
        DO 100 I = 1,512
	DO 100 J = 1,6
	COUNT(J,I) = 0
100	CONTINUE
C
C
        DO 200 I = 1,5
	DO 200 J = 1,500
	IF(X(I,J).GE.XLOW) GO TO 101
	JJ = 1
	GO TO 103
101	IF(X(I,J).LE.XHIGH) GO TO 102
	JJ = 512
	GO TO 103
C
102         JJ      = INT((X(I,J)-XLOW)/STEP) + 1
            IF (JJ .GT. 512) JJ = 512
103         COUNT(I,JJ) = COUNT(I,JJ) + 1
	    COUNT(6,JJ) = COUNT(6,JJ) + 1
C
200     CONTINUE
C
	IHI1 = COUNT(1,1)
	IHI2 = COUNT(6,1)
	DO 210 J = 1,512
	IHI2 = MAX0(IHI2,COUNT(6,J))
	DO 210 I = 1,5
	IHI1 = MAX0(IHI1,COUNT(I,J))
210	CONTINUE
C
	YLOW = 0.0
	YHIGH1 = FLOAT(IHI1) + 0.1 * FLOAT(IHI1)
	YHIGH2 = FLOAT(IHI2) + 0.1 * FLOAT(IHI2)
C
	CALL BGNPLT
	CALL GSCOLR(1,IERR)
	CALL MAPSIZ(0.0,100.0,14.0,50.0,0.0)
	CALL MAPSML(XLOW,XHIGH,YLOW,YHIGH1,SXLAB,SYLAB,STITLE,TYPE)

C
        DO 300 I = 1,5
	CALL GSCOLR(1+I,IERR)
	X0 = XLOW
	Y0 = 0.0
	CALL SCALE(X0,Y0,VX0,VY0)
	CALL GSMOVE(VX0,VY0)
C
	DO 400 J = 1,512
C
	X0   = XLOW + J * STEP
	Y0   = FLOAT(COUNT(I,J))
	CALL SCALE(X0,Y0,VX1,VY1)
	IF(Y0.GT.0.) GO TO 398
	CALL GSMOVE(VX1,VY0)
	GO TO 399
398	CALL GSDRAW(VX0,VY1)
	CALL GSDRAW(VX1,VY1)
	CALL GSDRAW(VX1,VY0)
C
399	VX0 = VX1
400	CONTINUE
300	CONTINUE
C
	CALL GSCOLR(1,IERR)
	CALL MAPSIZ(0.0,100.0,54.0,85.0,0.0)
	CALL MAPSML(XLOW,XHIGH,YLOW,YHIGH2,SXLAB,SYLAB,STITLE,TYPE)
	CALL GSCOLR(7,IERR)
C
	X0 = XLOW
	Y0 = 0.0
	CALL SCALE(X0,Y0,VX0,VY0)
	CALL GSMOVE(VX0,VY0)
C
	DO 500 J = 1,512
C
C
	X0   = XLOW + J * STEP
	Y0   = FLOAT(COUNT(6,J))
	CALL SCALE(X0,Y0,VX1,VY1)
	IF(Y0.GT.0.) GO TO 498
	CALL GSMOVE(VX1,VY0)
	GO TO 499
498	CALL GSDRAW(VX0,VY1)
	CALL GSDRAW(VX1,VY1)
	CALL GSDRAW(VX1,VY0)
C
499	VX0 = VX1
500	CONTINUE
999	VOLTS = VOLTS/1000.
	SIGVX = SIGVX * 10000.
	SIGVZ = SIGVZ * 10000.
	WRITE(LABEL1,510)INT(DIAMETER)
	WRITE(LABEL2,510)INT(VOLTS)
	WRITE(LABEL3,510)INT(V0Z)
	WRITE(LABEL4,510)INT(SIGVX)
	WRITE(LABEL5,510)INT(SIGVZ)
510	FORMAT(I6)
	VLEN1 = GSLENS('DROPLET DIAMETER (MICRONS) = '//CHAR(0))
	VLEN2 = GSLENS('PLATE KILOVOLTAGE = '//CHAR(0))
	VLEN3 = GSLENS('MEAN Z VELOCITY (METERS/SECOND) = '//CHAR(0))
	VLEN4 = GSLENS('X VELOCITY JITTER (MICRONS/SECOND) = '//CHAR(0))
	VLEN5 = GSLENS('Z VELOCITY JITTER (MICRONS/SECOND) = '//CHAR(0))
	CALL MAPSIZE( 0.0,100.0,0.0,12.0,0.0)
	CALL GSCOLR(1,IERR)
	CALL GSMOVE(0.,2.2)
	CALL GSPSTR('Droplet diameter (microns) = '//CHAR(0))
	CALL GSMOVE(0.,1.7)
	CALL GSPSTR('Plate Kilovoltage = '//CHAR(0))
	CALL GSMOVE(0.,1.2)
	CALL GSPSTR('Mean z Velocity (meters/second) = '//CHAR(0))
	CALL GSMOVE(0.,0.7)
	CALL GSPSTR('x Velocity Jitter (microns/second) = '//CHAR(0))
	CALL GSMOVE(0.,0.2)
	CALL GSPSTR('z Velocity Jitter (microns/second) = '//CHAR(0))
	CALL GSMOVE(VLEN1,2.2)
	CALL GSPSTR(LABEL1//CHAR(0))
	CALL GSMOVE(VLEN2,1.7)
	CALL GSPSTR(LABEL2//CHAR(0))
	CALL GSMOVE(VLEN3,1.2)
	CALL GSPSTR(LABEL3//CHAR(0))
	CALL GSMOVE(VLEN4,0.7)
	CALL GSPSTR(LABEL4//CHAR(0))
	CALL GSMOVE(VLEN5,0.2)
	CALL GSPSTR(LABEL5//CHAR(0))
C
C       BYE
C
        RETURN
        END
