Frozen Fish 1: Amiga/FrozenFish-Apr94.iso/bbs/alib/d2xx/d267/diglib.lha/Diglib/diglib.zoo/diglib/fivmet.for

Text File | 1989-06-20 | 5.2 KB | 240 lines

    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