Title                  :TWODISANL
Keywords               :CTEM, DISOLOCATIONS
Computer               :DEC VAX 11/730-785, DEC PDP 11/2-11/73
Operating System       :VAXVMS, RT-11
Programming Language   :Fortran IV
Hardware Requirements  :None
Author(s)              :R. J. Holton 
Correspondence Address :Argonne Nat. Lab, Electron Microscopy Center,Bldg 212
                       :Materials Science Division, Argonne, Illinois 60439,
Abstract:

TWODISANL is a modified version of the original TWODIS program developed
by the group of A.K. Head, P. Humble, L.M. Clarebrough, A.J. Morton and
C.T. Forwood at CSIRO Division of Tribophysics, University of Melbourne, 
Autralia for calculating disolcation images for the TEM. This version
modifies only their output algorithms for output to a serial data file rather
than overprinting upon a line printer. Using the included program (IMOUT) 
this can then be displayed as a randomized dither pattern on a DEC LA100 or 
LA50 printer
-------------------------------------------------------------------------------

                       
Title                  :TWODISANL
Keywords               :CTEM, DISOLOCATIONS
Computer               :DEC VAX 11/730-785, DEC PDP 11/2-11/73
Operating System       :VAXVMS, RT-11
Programming Language   :Fortran IV
Hardware Requirements  :None
Author(s)              :R. J. Holton 
Correspondence Address :Argonne Nat. Lab, Electron Microscopy Center,Bldg 212
                       :Materials Science Division, Argonne, Illinois 60439,
DOCUMENTATION:

TWODISANL is a modified version of the original TWODIS program developed
by the group of A.K. Head, P. Humble, L.M. Clarebrough, A.J. Morton and
C.T. Forwood at CSIRO Division of Tribophysics, University of Melbourne, 
Autralia for calculating disolcation images for the TEM. This version
modifies only their output algorithms for output to a serial data file rather
than overprinting upon a line printer. Using the included program (IMOUT) 
this can then be displayed as a randomized dither pattern on a DEC LA100 or 
LA50 printer.


Extensive documentation for the ONEDIS program can be found in the book:

Computed Electron Micrographs and Defect Identification
A.K. Head, P. Humble, L.M. Clarebrough, A.J. Morton and C.T. Forwood
CSIRO Division of Tribophysics
University of Melbourne, Autralia

1973 Noth Holland Publishing Co.


Compilation Procedure: RT11 

FORTRAN/NOSWAP/NOLINENUMBERS  TWODIS,ANCALC,NEWTON,DERIV,RKM,HALFTN
FORTRAN/NOSWAP/NOLINENUMBERS  PUTBLK,IRANDM

R LINK
WD:TWODIS,TM:=RD:TWODIS,RD:ANCALC,RD:NEWTON,RD:DERIV/C
RD:RKM,RD:HALFTN,RD:PUTBLK,RD:IRANDM
^C

Test Data:

SEE REFERENCE BOOK DOCUMENTED ABOVE 

-------------------------------------------------------------------------------

               
Title                  :TWODISANL
Keywords               :CTEM, DISOLOCATIONS
Computer               :DEC VAX 11/730-785, DEC PDP 11/2-11/73
Operating System       :VAXVMS, RT-11
Programming Language   :Fortran IV
Hardware Requirements  :None
Author(s)              :R. J. Holton
Correspondence Address :Argonne Nat. Lab, Electron Microscopy Center,Bldg 212
                       :Materials Science Division, Argonne, Illinois 60439,
SOURCE CODE:

THIS IS A COMBINATION OF THE ORIGINAL TWODIS SOURCE CODE AND THE MODIFIED
OUTPUT ROUTINE UTILIZED AT ANL

810       PROGRAM TWODIS
C
      INTEGER SYMBOL(3)
C
      DIMENSION LB(3),LU(3),LG(3),LBM(3),LFN(3),IY(15),IZ(26),BD(3),
     1 GD(3),BM(3),FN(3),FNX(3),DCX(3,3),DR(3),DI(3),UR(3,3),UI(3,3),
     2 VR(3,3),VI(3,3),CB(3),CU(3),CG(3),CBM(3),CFN(3),TB(160),
     3 TEMPY(8),FX(80,4)
     4 ,LB2(3),LFP(3,3),LS1(3),LS2(3),LS3(3),LS(3,3),LQ(3),CQ(3),FP(3)
     5 ,CB2(3),FP1X(3),FP2X(3),FP3X(3),B2D(3),CFP(3,3),FPX(3,3),VEC1(3)
     6 ,VEC2(3),ALPHA(3),SINA(3),COSA(3),POS(2,4),ITYPE(4)
      EQUIVALENCE (FP1X(1),FPX(1)),(FP2X(1),FPX(4)),(FP3X(1),FPX(7)),
     1            ( LS1(1), LS(1)),( LS2(1), LS(4)),( LS3(1), LS(7))
      COMMON/RKMDRV/CN(30),X,X1,Y(8),ERROR,SKIP,Q,D(8),ANO
      COMMON/ANCNEW/NEW,ZR,ZI,QR(7),QI(7),KRASH,C11,C12,C44,DC(3,3),
     1 C(6,6),PR(3),PI(3),AR(3,3),AI(3,3),ELR(3,3),ELI(3,3),EMR(3,3),
     2 EMI(3,3),B(3,3),H(3,3)
      COMMON/DATA/NP(3),NQ(3),MM(3),NN(3),L1(6),L2(6),L3(3,3),PY,MR,MW
C     DATA (NP=2,3,1),(NQ=3,1,2),(MM=1,6,5),(NN=6,2,4),(L1=1,2,3,2,3,1)
C     DATA (L2=1,2,3,3,1,2),(L3=1,6,5,6,2,4,5,4,3),(PY=3.1415926536)
C     DATA (MR=60),(MW=61)
      DATA NP/2,3,1/, NQ/3,1,2/, MM/1,6,5/, NN/6,2,4/, L1/1,2,3,2,3,1/
      DATA L2/1,2,3,3,1,2/, L3/1,6,5,6,2,4,5,4,3/, PY/3.1415926536/
      DATA MR/2/, MW/7/
      DATA SYMBOL /' ', 'M', 'D'/

7     WRITE (7,1)
1     FORMAT (1X,'Enter input data file name: ',$)
      CALL ASSIGN(MR,DUMMY,-1,'OLD')
10    READ(MR,20)C11,C12,C44,IY,ANO
C     WRITE (6,*) C11,C12,C44,ANO
20    FORMAT(3F10.0,15A1,25X,F9.0)
      IF(C11)30,30,40
30    STOP
40    CONTINUE
1733  WRITE(7,1735)
1735  FORMAT(' ','Output to a (F)ile, (L)A-100, or (P)eritek? ',$)
      READ(5,1736) IANS
1736  FORMAT(A1)
      IF (IANS .EQ. 'F') GOTO 1799
      IF (IANS .EQ. 'L') GOTO 1755
      IF (IANS .NE. 'P') GOTO 1733
      IDEV = 0
      ITV = 0
      GOTO 1800
1755  NESC = 27
      WRITE(6,1734) NESC
1734  FORMAT(' ',A1,'P1q')
      IDEV = -1
      GOTO 1800
1799  WRITE (7,11)
11    FORMAT(' Enter output data file name: ',$)
      IDEV = 3
      CALL ASSIGN(IDEV,DUMMY,-1,'NEW','CC')
1800  CONTINUE
      READ(MR,50)LB,LD,LU,LG,LBM,LFN,W,THICK,START,FINISH,LPR,IZ
50    FORMAT(3I2,X,I1,X,3I2,X,3I2,X,3I2,X,3I2,4F4.0,I2,26A1)
      READ(MR,60)LB2,LD2,SEP,LFP,LS1,LQ(1),LS2,LQ(2),LS3,LQ(3)
60    FORMAT(3I2,X,I1,F4.0,3I2,X,3I2,X,3I2,X,3I2,X,I1,X,3I2,X,I1,X,
     &       3I2,X,I1)
C     TIME=TIMEF(X)
      CN(14)=2.0*W
      IF(LD)80,70,80
70    LD=1
80    IF(ANO)100,90,100
90    ANO=0.1
100   ANO=-ANO
      IF(FINISH)130,110,130
110   IF(START)130,120,130
120   START=0.0
      FINISH=THICK
130   IF(FINISH-START)140,140,160
140   WRITE(MW,150)
      GO TO 1350
150   FORMAT(/,/,20H START AFTER FINISH   ,/,/)
160   IF(LPR)180,170,180
170   LPR=1
180   IF(LFN(1)**2+LFN(2)**2+LFN(3)**2)210,190,210
190   DO 200 J=1,3
200   LFN(J)=LBM(J)
210   IF(LD2)220,220,230
220   LD2=1
230   DO 270 J=1,3
      IF(LQ(J))240,240,250
240   LQ(J)=1
250   CQ(J)=LQ(J)
      DO 260  K=1,3
260   CFP(J,K)=LFP(J,K)
      CB2(J)=FLOAT(LB2(J))/FLOAT(LD2)
      CB(J)=FLOAT(LB(J))/FLOAT(LD)
      CU(J)=LU(J)
      CG(J)=LG(J)
      CBM(J)=LBM(J)
270   CFN(J)=LFN(J)
      DO 280 J=1,3
      DC(3,J)=CU(J)
      K=NP(J)
      L=NQ(J)
280   DC(1,J)=CBM(K)*CU(L)-CBM(L)*CU(K)
      DO 290 J=1,3
      K=NP(J)
      L=NQ(J)
290   DC(2,J)=DC(3,K)*DC(1,L)-DC(3,L)*DC(1,K)
      DO 340 J=1,3
      Z=0.0
      DO 300 K=1,3
300   Z=Z+DC(J,K)**2
      IF(Z-0.0001)310,310,330
310   WRITE(MW,320)
      GO TO 1350
320   FORMAT(/,/,16H BEAM PARALLEL U  ,/,/)
330   Z=1.0/SQRT(Z)
      DO 340 K=1,3
340   DC(J,K)=DC(J,K)*Z
      DO 350 J=1,3
      DCX(1,J)=-DC(1,J)
350   DCX(2,J)=-CBM(J)
      DO 360 J=1,3
      K=NP(J)
      L=NQ(J)
360   DCX(3,J)=DCX(1,K)*DCX(2,L)-DCX(1,L)*DCX(2,K)
      DO 380 J=1,3
      Z=0.0
      DO 370 K=1,3
370   Z=Z+DCX(J,K)**2
      Z=1.0/SQRT(Z)
      DO 380 K=1,3
380   DCX(J,K)=DCX(J,K)*Z
      DO 410 J=1,3
      B2D(J)=0.0
      FP(J)=0.0
      DO 390  K=1,3
390   FPX(J,K)=0.0
      BD(J)=0.0
      GD(J)=0.0
      BM(J)=0.0
      FN(J)=0.0
      FNX(J)=0.0
      DO 410 K=1,3
      B2D(J)=B2D(J)+DC(J,K)*CB2(K)
      FP(J)=FP(J)+DC(J,K)*CFP(K,2)
      DO 400  L=1,3
400   FPX(J,L)=FPX(J,L)+DCX(J,K)*CFP(K,L)
      BD(J)=BD(J)+DC(J,K)*CB(K)
      BM(J)=BM(J)+DC(J,K)*CBM(K)
      FN(J)=FN(J)+DC(J,K)*CFN(K)
      FNX(J)=FNX(J)+DCX(J,K)*CFN(K)
410   GD(J)=GD(J)+DC(J,K)*CG(K)
      IF(LBM(1)*LG(1)+LBM(2)*LG(2)+LBM(3)*LG(3))420,440,420
420   WRITE(MW,430)
430   FORMAT(/,/,19H BEAM NOT PERP TO G   ,/,/)
      GO TO 1350
440   Z=SQRT(FN(1)**2+FN(2)**2+FN(3)**2)
      X=SQRT(BM(1)**2+BM(2)**2+BM(3)**2)
      DO 450 J=1,3
      BM(J)=BM(J)/X
450   FN(J)=FN(J)/Z
      FNBM=0.0
      DO 460 J=1,3
460   FNBM=FNBM+FN(J)*BM(J)
      IF(FN(3))470,470,490
470   WRITE(MW,480)
      GO TO 1350
480   FORMAT(/,/,29H U AND FOIL NORMAL NOT ACUTE   ,/,/)
490   IF(FNBM)500,500,520
500   WRITE(MW,510)
      GO TO 1350
510   FORMAT(/,/,32H FOIL NORMAL AND BEAM NOT ACUTE  ,/,/)
520   DO 550  J=1,3
      Z=LS2(J)
      IF(FP2X(2))530 ,530 ,540
530   Z=-Z
540   VEC1(J)=FLOAT(LS1(J))/CQ(1)+Z/CQ(2)+CB(J)
550   VEC2(J)=-Z/CQ(2)-FLOAT(LS3(J))/CQ(3)+CB2(J)
      DO 580  J=1,3
      IF(LU(1)*LFP(1,J)+LU(2)*LFP(2,J)+LU(3)*LFP(3,J))560 ,580 ,560
560   WRITE(MW,570 ) J
      GO TO 1350
570   FORMAT(/,/,22H U NOT IN FAULT PLANE ,I2,/,/)
580   CONTINUE
      DO 590  J=1,3
      M=LG(1)*LS(1,J)*LG(2)*LS(2,J)+LG(3)*LS(3,J)
      Z=M/LQ(J)
      ALPHA(J)=2.0*PY*(FLOAT(M)/CQ(J)-Z)
      SINA(J)=SIN(ALPHA(J))
590   COSA(J)=COS(ALPHA(J))
      CALL ANCALC
      IF(KRASH)1350,600,1350
600   CONTINUE
      DO 610 J=1,3
      DR(J)=0.0
      DI(J)=0.0
      DO 610 K=1,3
      DR(J)=DR(J)+GD(K)*AR(K,J)
610   DI(J)=DI(J)+GD(K)*AI(K,J)
      DO 620 J=1,3
      Z=DR(J)
      DR(J)=Z*PR(J)-DI(J)*PI(J)
620   DI(J)=Z*PI(J)+DI(J)*PR(J)
      DO 630 JA=1,3
      DO 630 L=1,3
      UR(JA,L)=0.0
      UI(JA,L)=0.0
      DO 630 J=1,3
      UR(JA,L)=UR(JA,L)+EMR(JA,J)*H(J,L)
630   UI(JA,L)=UI(JA,L)+EMI(JA,J)*H(J,L)
      DO 640 JA=1,3
      DO 640 L=1,3
      VR(JA,L)=DR(JA)*UR(JA,L)-DI(JA)*UI(JA,L)
640   VI(JA,L)=DR(JA)*UI(JA,L)+DI(JA)*UR(JA,L)
      DO 650 JA=1,3
      DO 650 L=1,3
650   UR(JA,L)=VR(JA,L)*PR(JA)+VI(JA,L)*PI(JA)
      DO 660 J=1,3
      CN(J+6)=PR(J)
      CN(J+9)=PI(J)**2
      CN(J)=0.0
      CN(J+3)=0.0
      CN(J+20)=0.0
      CN(J+23)=0.0
      DO 660 L=1,3
      CN(J)=CN(J)+VR(J,L)*BD(L)
      CN(J+3)=CN(J+3)+UR(J,L)*BD(L)
      CN(J+20)=CN(J+20)+VR(J,L)*B2D(L)
660   CN(J+23)=CN(J+23)+UR(J,L)*B2D(L)
      GO TO (670,680,690,700,710,720,730,740,750),LPR
670   BLACK=0.313
      WHITE=1.154
      GO TO 760
680   BLACK=0.129
      WHITE=1.253
      GO TO 760
690   BLACK=0.058
      WHITE=1.340
      GO TO 760
700   BLACK=0.397
      WHITE=1.462
      GO TO 760
710   BLACK=0.191
      WHITE=1.866
      GO TO 760
720   BLACK=0.098
      WHITE=2.280
      GO TO 760
730   BLACK=0.259
      WHITE=0.959
      GO TO 760
740   BLACK=0.089
      WHITE=0.869
      GO TO 760
750   BLACK=0.035
      WHITE=0.809
760   NROWS = 64	!128	!MUST BE EVEN
      NCOLS = 155		!MUST BE ODD
      HCPI = 1.
      VCPI = 1.
      ROWS = NROWS
      COLS = NCOLS
      NEND = (NCOLS - 1) / 2
      NEND2 = 2 * NEND
      CALL HALFTN(NCOLS,TB,MARK,BLACK,WHITE,-1,IDEV)
      CALL HALFTN(NCOLS,TB,MARK,BLACK,WHITE,0,IDEV)
C      CALL HALFTN(ITV,129,TB,MARK,BLACK,WHITE,-1)
C      CALL HALFTN(ITV,129,TB,MARK,BLACK,WHITE,0)
C      WRITE(MW,770)
C770   FORMAT(1H1)
      Z=SQRT(FP(1)**2+FP(2)**2)
      IF(Z)810 ,780 ,810
780   Z=1.0
      FP2X(1)=1.0
      IF(SEP)790,810 ,790
790   WRITE(MW,800)
      GO TO 1350
800   FORMAT(/,/,36H FAULT PLANE 2 ZERO WITH SEP NONZERO,/,/)
810   PT=SEP*FP(2)/Z
      SL=-SEP*FP(1)/(Z*BM(2))
      CN(16)=PY*PT/2.0
      CN(17)=PY*SL/2.0
      CN(30)=CN(16)/BM(2)
      EXTRA=ABS(SL+PT*FNX(1)/FNX(2))
      DIVISR =BM(3)/BM(2)-FNX(3)/FNX(2)
      THBM=THICK/FNBM
      TBP=PY*THBM
      FRACTN=(FINISH-START)/THICK
C     DELT=PY*FRACTN*(THBM+EXTRA)/64.0
C     DELT=FRACTN*TBP/64.0
      DELT=PY*FRACTN*(THBM+EXTRA)/((COLS-1.0)/2.0)
      WL=((THICK*BM(2)/FN(3))+EXTRA/DIVISR)*FRACTN
C     WW=59.0*10.0*WL/(6.0*128.0)
C     WW=59.0*10.0*WL/(6.0*128.0)
      WW = (ROWS-1.0)*HCPI*WL / (VCPI*(COLS-1.0))
C     DELW=PY*WW/59.0
C     DELW = PY * WW / 59.0
      DELW = PY * WW / (ROWS - 1.0)
C     SHIFT=DELW/2.0
      SHIFT = DELW * .5
      IF(SHIFT-0.01)830 ,820 ,820
820   SHIFT=0.01
830   CN(29)=1000.0
      X=0.0
      Q=0.0
      ERROR=0.0001
      DO 840 JK=1,8
840   Y(JK)=0.0
      Y(1)=1.0
      X1=DELT
      CALL RKM
      X1=TBP
      CALL RKM
      BACK=Y(1)**2+Y(2)**2
C     DO 1340 JC=1,60
      DO 1340 JC=1,NROWS
C     CN(15)=(FLOAT(JC)-30.5)*DELW
      CN(15) = (FLOAT(JC)-((ROWS*.5)+.5))*DELW
      MOVE=0
      Z=CN(15)-CN(16)
      IF(SHIFT-ABS(Z))870 ,870 ,850
850   CN(15)=CN(16)+SHIFT
      IF(Z)860 ,900 ,900
860   CN(15)=CN(16)-SHIFT
      GO TO 900
870   Z=CN(15)+CN(16)
      IF(SHIFT-ABS(Z))910 ,910 ,880
880   CN(15)=-CN(16)+SHIFT
      IF(Z)890 ,900 ,900
890   CN(15)=-CN(16)-SHIFT
900   MOVE=1
910   XXX=CN(15)+CN(16)
      YYY=CN(15)-CN(16)
      CN(29)=CN(15)/BM(2)
      DO 920 J=1,4
      POS(1,J)=-10000.0
920   ITYPE(J)=J
      IF(YYY)950,950,930
930   IF(ALPHA(1)*FP1X(2))940,950,940
940   POS(1,1)=CN(17)-(CN(15)-CN(16))*FP1X(1)/FP1X(2)
950   IF(XXX*YYY)960,980,980
960   IF(ALPHA(2)*FP2X(2))970,980,970
970   POS(1,2)=-CN(15)*FP2X(1)/FP2X(2)
980   IF(XXX)990,1010,1010
990   IF(ALPHA(3)*FP3X(2))1000,1010,1000
1000  POS(1,3)=-CN(17)-(CN(15)+CN(16))*FP3X(1)/FP3X(2)
1010  STARTA=PY*(EXTRA/2.0-(THBM+EXTRA)*FINISH/THICK)-CN(15)*
     1  FNX(1)/FNX(2)
      SURFAC=STARTA+TBP
      POS(1,4)=SURFAC
      DO 1040 J=1,3
      LUCK=0
      DO 1030 K=1,3
      IF(POS(1,K)-POS(1,K+1))1030,1030,1020
1020  Z=POS(1,K+1)
      POS(1,K+1)=POS(1,K)
      POS(1,K)=Z
      ISTORE=ITYPE(K+1)
      ITYPE(K+1)=ITYPE(K)
      ITYPE(K)=ISTORE
      LUCK=-1
1030  CONTINUE
      IF(LUCK)1040,1050,1040
1040  CONTINUE
1050  LSWTCH=0
      DO 1100 J=1,4
      IF(ITYPE(J)-4)1060,1080,1060
1060  IF(LSWTCH)1070,1090,1070
1070  POS(2,J)=POS(1,J)
      GO TO 1100
1080  LSWTCH=-1
1090  POS(2,J)=-10050.0
1100  CONTINUE
      X=STARTA
      X1=X
      IFLAG=0
      DO 1110 JK=1,8
1110  Y(JK)=0.0
      Y(1)=1.0
      Y(7)=1.0
      KOUNTF=1
      N=1
      LINK=1
C     DO 1130 JT=1,64
      DO 1130 JT=1,NEND
      X1=X1+DELT
      GO TO 1390
1120  CALL RKM
      DNR=Y(1)*Y(7)-Y(2)*Y(8)-Y(3)*Y(5)+Y(4)*Y(6)
      DNI=Y(1)*Y(8)+Y(2)*Y(7)-Y(3)*Y(6)-Y(4)*Y(5)
      DNN=1.0/(DNR**2+DNI**2)
      FX(JT,1)=DNN*(Y(7)*DNR+Y(8)*DNI)
      FX(JT,2)=DNN*(Y(8)*DNR-Y(7)*DNI)
      FX(JT,3)=-DNN*(Y(3)*DNR+Y(4)*DNI)
      FX(JT,4)=DNN*(Y(3)*DNI-Y(4)*DNR)
1130  CONTINUE
      IF(IFLAG)1150,1140,1150
1140  X1=SURFAC
      LINK=3
      GO TO 1390
1150  X=SURFAC
      X1=X
      DO 1160 JK=1,8
1160  Y(JK)=TEMPY(JK)
      KOUNTF=1
      N=2
      LINK=2
C     DO 1180 JM=1,64
      DO 1180 JM=1,NEND
      X1=X1+DELT
      GO TO 1390
1170  CALL RKM
      TT=(FX(JM,1)*Y(1)-FX(JM,2)*Y(2)+FX(JM,3)*Y(5)-FX(JM,4)*Y(6))**2
     1 +(FX(JM,1)*Y(2)+FX(JM,2)*Y(1)+FX(JM,3)*Y(6)+FX(JM,4)*Y(5))**2
1180  TB(2*JM+1)=TT/BACK
      TB(1)=(TEMPY(1)**2+TEMPY(2)**2)/BACK
C     DO 1190 JZ=2,128,2
      DO 1190 JZ=2,NEND2,2
1190  TB(JZ)=0.5*(TB(JZ-1)+TB(JZ+1))
      MARK=SYMBOL(1)
      STAGR=EXTRA/DIVISR
C     LSTAG=STAGR*129.0/WL
      LSTAG=STAGR*COLS/WL
C      DELL=DELW/2.0+0.00000001
      DELL = DELW*.5+0.00000001
      EFP2X1=FP2X(1)/SQRT(FP2X(1)**2+FP2X(2)**2)
      EFNX1=FNX(1)/SQRT(FNX(1)**2+FNX(2)**2)
      IF(ABS(XXX)-DELL) 1210,1210,1200
1200  IF(ABS(YYY)-DELL) 1210,1210,1330
1210  IF(FP2X(2))1220,1230,1240
1220  IF(EFP2X1-EFNX1) 1250,1270,1260
1230  IF(SEP) 1270,1270,1280
1240  IF(EFP2X1+EFNX1) 1260,1270,1250
1250  IF(CN(15)) 1290,1330,1300
1260  IF(CN(15)) 1300,1330,1290
1270  TB(1)=-1.0
C     TB(129)=-1.0
      TB(NCOLS)=-1.0
      GO TO 1310
1280  TB(1)=-1.0
C     TB(129-LSTAG)=-1.0
      TB(NCOLS-LSTAG)=-1.0
1290  TB(LSTAG+1)=-1.0
C     TB(129)=-1.0
      TB(NCOLS)=-1.0
      GO TO 1310
1300  TB(1)=-1.0
C     TB(129-LSTAG)=-1.0
      TB(NCOLS-LSTAG)=-1.0
1310  MARK=SYMBOL(3)
      IF(MOVE)1320,1330,1320
1320  MARK=SYMBOL(2)
C1330  CALL HALFTN(ITV,129,TB,MARK,BLACK,WHITE,1)
1330  CALL HALFTN(NCOLS,TB,MARK,BLACK,WHILE,1,IDEV)
1340  CONTINUE
C     TIME=(TIMEF(X)-TIME)/1000.0
      TIME = 0.0
1350  WRITE(MW,1360) C11,C12,C44,IY,TIME,WL,WW,START,FINISH,THICK,THBM,
     1VEC1
1360  FORMAT(1H ,F5.2,4HC11 F5.2,4HC12 F5.2,4HC44 15A1
     1  ,F6.1,5H SECSF6.2,3H WLF6.2,3H WWF5.2,6H STRT F5.2,5H FIN 
     2F7.3,2HTHF7.3,4HTHBM3F6.2,2HV1 )
      WRITE(MW,1370)LB,LD,LU,LG,LBM,LFN,W,ANO,BACK,BLACK,WHITE,VEC2
1370  FORMAT(1H 3I2,1H/I1,5HB    3I2,5HU    3I2,5HG    3I2,5HBM   3I2,
     12HFNF7.3,4HW   F6.3,5HANO  F6.3,6HBACK  F5.3,7HBLACK  F5.3,7HWHITE 
     2  3F6.2,2HV2 )
      WRITE(MW,1380)LB2,LD2,SEP,LFP,LS1,LQ(1),LS2,LQ(2),LS3,LQ(3),IZ
1380  FORMAT(1H 3I2,1H/I1,4HB2  F5.2,5HSEP  3I2,5HFP1  3I2,4HFP2 3I2,
     1 5HFP3  3I2,1H/I1,5HSH1  3I2,1H/I1,5HSH2  3I2,1H/I1,5HSH3  26A1,
     2 12HTWODIS TRIBO )
      IF (IDEV .LT. 0) GOTO 1885
      IF (IDEV .GT. 0) GOTO 1887
      GOTO 10
1885  WRITE(6,1835) NESC
1835  FORMAT(' ','-',A1,'\')
      GO TO 10
1887  CALL CLOSE(3)
      GOTO 10
1390  IF(KOUNTF-5)1410,1400,1400
1400  GO TO (1120,1170,1150),LINK
1410  IF(X1-POS(N,KOUNTF))1400,1420,1420
1420  IF(X-POS(N,KOUNTF))1430,1430,1460
1430  XX1=X1
      X1=POS(N,KOUNTF)
      CALL RKM
      IP=ITYPE(KOUNTF)
      GO TO (1440,1440,1440,1470),IP
1440  Z=Y(3)
      ZZ=Y(7)
      Y(3)=Y(3)*COSA(IP)-Y(4)*SINA(IP)
      Y(7)=Y(7)*COSA(IP)-Y(8)*SINA(IP)
      Y(4)=Y(4)*COSA(IP)+Z*SINA(IP)
      Y(8)=Y(8)*COSA(IP)+ZZ*SINA(IP)
1450  X1=XX1
      POS(N,KOUNTF)=-9000.0
1460  KOUNTF=KOUNTF+1
      GO TO 1390
1470  DO 1480 JK=1,8
1480  TEMPY(JK)=Y(JK)
      IFLAG=1
      GO TO 1450
      END
      SUBROUTINE ANCALC
      DIMENSION D(6,6), DR(3,3),DI(3,3),G(9),E(9)
      COMMON/DATA/NP(3),NQ(3),MM(3),NN(3),L1(6),L2(6),L3(3,3),PY,MR,MW
      COMMON/ANCNEW/NEW,ZR,ZI,QR(7),QI(7),KRASH,C11,C12,C44,DC(3,3),
     1 C(6,6),PR(3),PI(3),AR(3,3),AI(3,3),ELR(3,3),ELI(3,3),EMR(3,3),
     2 EMI(3,3),B(3,3),H(3,3)
      DO 30  JA=1,6
      DO 30  JB=1,6
30    D(JA,JB)=0.0
      D(1,1)=C11/C44
      D(2,2)=D(1,1)
      D(3,3)=D(1,1)
      D(1,2)=C12/C44
      D(2,1)=D(1,2)
      D(2,3)=D(1,2)
      D(3,2)=D(1,2)
      D(1,3)=D(1,2)
      D(3,1)=D(1,2)
      D(4,4)=1.0
      D(5,5)=1.0
      D(6,6)=1.0
      DO 34 M=1,6
      I=L1(M)
      J=L2(M)
      DO 34 N=1,M
      K=L1(N)
      L=L2(N)
      X=0.0
      DO 33 LP=1,3
      Y=0.0
      DO 32 LQ=1,3
      LT=L3(LP,LQ)
32    Y=Y+DC(J,LQ)*
     1   (DC(K,1)*(DC(L,1)*D(LT,1)+DC(L,2)*D(LT,6)+DC(L,3)*D(LT,5))
     2   +DC(K,2)*(DC(L,1)*D(LT,6)+DC(L,2)*D(LT,2)+DC(L,3)*D(LT,4))
     3   +DC(K,3)*(DC(L,1)*D(LT,5)+DC(L,2)*D(LT,4)+DC(L,3)*D(LT,3)))
33    X=X+DC(I,LP)*Y
      C(M,N)=X
34    C(N,M)=X
      G(1)=C(5,5)
      G(2)=2.0*C(4,5)
      G(3)=C(4,4)
      G(4)=C(6,6)
      G(5)=2.0*C(2,6)
      G(6)=C(2,2)
      G(7)=C(1,1)
      G(8)=2.0*C(1,6)
      G(9)=C(6,6)
      E(1)=C(5,6)
      E(2)=C(2,5)+C(4,6)
      E(3)=C(2,4)
      E(4)=C(1,5)
      E(5)=C(5,6)+C(1,4)
      E(6)=C(4,6)
      E(7)=C(1,6)
      E(8)=C(6,6)+C(1,2)
      E(9)=C(2,6)
      DO 50 KP=1,7
      QR(KP)=0.0
50    QI(KP)=0.0
      DO 51 KQ=1,3
      DO 51 KR=1,3
      DO 51 KS=1,3
      KT=KQ+KR+KS-2
51    QR(KT)=QR(KT)+G(KQ)*G(KR+3)*G(KS+6)+2.0*E(KQ)*E(KR+3)*E(KS+6)-
     1 E(KQ)*E(KR)*G(KS+6)-E(KQ+3)*E(KR+3)*G(KS+3)-E(KQ+6)*E(KR+6)*G(KS)
      DO 52 KP=1,7
52    QR(KP)=QR(KP)/QR(7)
      KRASH=0
      NEW=7
      ZR=0.1
      ZI=1.0
      CALL NEWTON
      IF(KRASH.NE.0)GOTO 69
61    PR(1)=ZR
      PI(1)=ABS(ZI)
      ZI=-ZI
      CALL NEWTON
      IF(KRASH.NE.0)GOTO 69
62    ZR=0.5
      ZI=0.9
      CALL NEWTON
      IF(KRASH.NE.0)GOTO 69
63    PR(2)=ZR
      PI(2)=ABS(ZI)
      ZI=-ZI
      CALL NEWTON
      IF(KRASH.NE.0)GOTO 69
64    ZR=-ZR
      CALL NEWTON
      IF(KRASH.NE.0)GOTO 69
65    PR(3)=ZR
      PI(3)=ABS(ZI)
      ZR=-C(4,5)/C(4,4)
      ZI=SQRT(ABS(C(4,4)*C(5,5)-C(4,5)**2))/C(4,4)
      DO 67 N=1,2
      IF((ZR-PR(N))**2*(ZI-PI(N))**2.GE.(ZR-PR(N+1))**2-(ZI-PI(N+1))**2)
     &  GOTO 67
66    Z=PR(N)
      PR(N)=PR(N+1)
      PR(N+1)=Z
      Z=PI(N)
      PI(N)=PI(N+1)
      PI(N+1)=Z
67    CONTINUE
      GO TO 74
69    IF (KRASH .LT. 0) GOTO 70
      WRITE(MW,73)
      GO TO 99
70    WRITE(MW,71)
      GO TO 99
71    FORMAT(/,/,11H NOCONVERGE,/,/)
73    FORMAT(/,/,9H REALROOT,/,/)
74    CONTINUE
      DO 80 K=1,3
      I=NP(K)
      L=NQ(K)
      PRK=PR(K)
      PIK=PI(K)
      SQR=PRK**2-PIK**2
      SQI=2.0*PRK*PIK
      DR(1,1)=C(1,1)+PRK*2.0*C(1,6)+SQR*C(6,6)
      DR(2,2)=C(6,6)+PRK*2.0*C(2,6)+SQR*C(2,2)
      DR(3,3)=C(5,5)+PRK*2.0*C(4,5)+SQR*C(4,4)
      DR(1,2)=C(1,6)+PRK*(C(1,2)+C(6,6))+SQR*C(2,6)
      DR(2,3)=C(5,6)+PRK*(C(4,6)+C(2,5))+SQR*C(2,4)
      DR(3,1)=C(1,5)+PRK*(C(1,4)+C(5,6))+SQR*C(4,6)
      DR(2,1)=DR(1,2)
      DR(3,2)=DR(2,3)
      DR(1,3)=DR(3,1)
      DI(1,1)=PIK*2.0*C(1,6)+SQI*C(6,6)
      DI(2,2)=PIK*2.0*C(2,6)+SQI*C(2,2)
      DI(3,3)=PIK*2.0*C(4,5)+SQI*C(4,4)
      DI(1,2)=PIK*(C(1,2)+C(6,6))+SQI*C(2,6)
      DI(2,3)=PIK*(C(4,6)+C(2,5))+SQI*C(2,4)
      DI(3,1)=PIK*(C(1,4)+C(5,6))+SQI*C(4,6)
      DI(2,1)=DI(1,2)
      DI(3,2)=DI(2,3)
      DI(1,3)=DI(3,1)
      DO 80 J=1,3
      M=NP(J)
      N=NQ(J)
      AR(J,K)=DR(I,M)*DR(L,N)-DI(I,M)*DI(L,N)-DR(I,N)*DR(L,M)+DI(I,N)*
     1 DI(L,M)
80    AI(J,K)=DR(I,M)*DI(L,N)+DI(I,M)*DR(L,N)-DR(I,N)*DI(L,M)-DI(I,N)*
     1 DR(L,M)
      DO 82 J=1,3
      NJ=NN(J)
      DO 82 K=1,3
      XR=0.0
      XI=0.0
      DO 81 L=1,3
      NL=NN(L)
      ML=MM(L)
      YR=C(NJ,ML)+C(NJ,NL)*PR(K)
      YI=C(NJ,NL)*PI(K)
      XR=XR+YR*AR(L,K)-YI*AI(L,K)
81    XI=XI+YI*AR(L,K)+YR*AI(L,K)
      ELR(J,K)=XR
82    ELI(J,K)=XI
      DO 83 J=1,3
      J1=NP(J)
      J2=NQ(J)
      DO 83 K=1,3
      K1=NP(K)
      K2=NQ(K)
      EMR(K,J) = ELR(J1,K1)*ELR(J2,K2) -ELI(J1,K1)*ELI(J2,K2)
     1           -ELR(J1,K2)*ELR(J2,K1)+ELI(J1,K2)*ELI(J2,K1)
83    EMI(K,J) = ELR(J1,K1)*ELI(J2,K2) +ELI(J1,K1)*ELR(J2,K2)
     1           -ELR(J1,K2)*ELI(J2,K1) -ELI(J1,K2)*ELR(J2,K1)
      DELR=0.0
      DELI=0.0
      DO 84 J=1,3
      DELR = DELR  +ELR(3,J)*EMR(J,3)  -ELI(3,J)*EMI(J,3)
84    DELI = DELI  +ELR(3,J)*EMI(J,3)  +ELI(3,J)*EMR(J,3)
      AUMR = DELR/(DELR**2+DELI**2)
      AUMI=-DELI/(DELR**2+DELI**2)
      DO 85 J=1,3
      DO 85 K=1,3
      X = EMR(J,K)*AUMR -EMI(J,K)*AUMI
      EMI(J,K) = EMR(J,K)*AUMI +EMI(J,K)*AUMR
85    EMR(J,K) = X
      DO 86 I=1,3
      DO 86 J=1,3
      B(I,J)=0.0
      DO 86 K=1,3
86    B(I,J)=B(I,J) -AR(I,K)*EMI(K,J) -AI(I,K)*EMR(K,J)
      DO 87 I=1,3
      I1=NP(I)
      I2=NQ(I)
      DO 87 J=1,3
      J1=NP(J)
      J2=NQ(J)
87    H(I,J) =  B(I1,J1)*B(I2,J2) -B(I1,J2)*B(I2,J1)
      DEL=B(3,1)*H(3,1)+B(3,2)*H(3,2)+B(3,3)*H(3,3)
      DO 88 I=1,3
      DO 88 J=1,3
88    H(I,J)=H(I,J)/DEL
99    RETURN
      END
      SUBROUTINE DERIV
      COMMON/RKMDRV/CN(30),X,X1,Y(8),ERROR,SKIP,Q,D(8),ANO
      IF (SKIP .NE. 0) GOTO 3
      IF (X .EQ. 0) X = .0000000001
2     R=CN(29)/X
      BETA=CN(14)+((R*CN(1)+CN(4))/((R+CN(7))**2+CN(10))+(R*CN(2)+
     1  CN(5))/((R+CN(8))**2+CN(11))+(R*CN(3)+CN(6))/((R+CN(9))**2+
     2  CN(12)))/X
3     Z=ANO*(Y(1)+Y(3))
      D(1)=Z-Y(4)
      D(3)=-BETA*Y(4)+Z-Y(2)
      Z=ANO*(Y(2)+Y(4))
      D(2)=Z+Y(3)
      D(4)=BETA*Y(3)+Z+Y(1)
      Z=ANO*(Y(5)+Y(7))
      D(5)=Z-Y(8)
      D(7)=-BETA*Y(8)+Z-Y(6)
      Z=ANO*(Y(6)+Y(8))
      D(6)=Z+Y(7)
      D(8)=BETA*Y(7)+Z+Y(5)
      RETURN
      END
      SUBROUTINE DERIV
      COMMON/RKMDRV/CN(30),X,X1,Y(8),ERROR,SKIP,Q,D(8),ANO
      IF(SKIP)7,1,7
1     X11=X-CN(17)
      IF(X11)3,2,3
2     X11=0.000000001
3     R1=(CN(29)-CN(30))/X11
      X22=X+CN(17)
      IF(X22)6,5,6
5     X22=0.000000001
6     R2=(CN(29)+CN(30))/X22
      BETA=CN(14)
     1+(((R1*CN(1)+CN(4))/((R1+CN(7))**2+CN(10)))+((R1*CN(2)+CN(5))/((
     2R1+CN(8))**2+CN(11)))+((R1*CN(3)+CN(6))/((R1+CN(9))**2+CN(12))))
     3 /X11
     4+(((R2*CN(21)+CN(24))/((R2+CN(7))**2+CN(10)))+((R2*CN(22)+CN(25))/
     5((R2+CN(8))**2+CN(11)))+((R2*CN(23)+CN(26))/((R2+CN(9))**2+CN(12))
     6))/X22
7     Z=ANO*(Y(1)+Y(3))
      D(1)=Z-Y(4)
      D(3)=-BETA*Y(4)+Z-Y(2)
      Z=ANO*(Y(2)+Y(4))
      D(2)=Z+Y(3)
      D(4)=BETA*Y(3)+Z+Y(1)
      Z=ANO*(Y(5)+Y(7))
      D(5)=Z-Y(8)
      D(7)=-BETA*Y(8)+Z-Y(6)
      Z=ANO*(Y(6)+Y(8))
      D(6)=Z+Y(7)
      D(8)=BETA*Y(7)+Z+Y(5)
      RETURN
      END
	FUNCTION IRANDM (LOW, HIGH)
C
	INTEGER LOW, HIGH
C
C	THIS FUNCTION RETURNS A RANDOM INTEGER UNIFORMLY DISTRIBUTED
C	FROM LOW TO HIGH.
C
C	NOTE THAT THIS PROGRAM MAKE A CALL TO A RANDOM NUMBER GENERATOR.
C
	INTEGER RANX1, RANX2
	COMMON /IRANX /RANX1, RANX2
C	
	IF (LOW .GE. HIGH) GOTO 99
C
	IRANDM = IFIX((HIGH - LOW + 1) * RAN(RANX1,RANX2)) + LOW
	RETURN
C
99	IRANDM = LOW
	RETURN
	END
      SUBROUTINE RKM
      DIMENSION YT(8),DT(8,3)
      COMMON/RKMDRV/CN(30),X,X1,Y(8),ERROR,SKIP,Q,D(8),ANO
      LAST=0
      IF (Q .NE. 0.0) GOTO 16
1     Q=X1-X
      SKIP=0.0
      Q15=0.0
      ERHIGH=ERROR*5.0
      ERLOW=ERHIGH*0.03125
2     LAST=1
      Q1=Q
      Q=X1-X
3     XT=X
      DO 4 M=1,8
4     YT(M)=Y(M)
5     IF (Q*1.5 .EQ. Q15) GOTO 7
6     Q2 = Q * .5
      Q3 = Q * .3333333333333
      Q4=4.0*Q3
      Q6 = Q * .1666666666667
      Q15=Q*1.5
7     CALL DERIV
      DO 8 M=1,8
      DT(M,1)=Q3*D(M)
8     Y(M)=DT(M,1)+YT(M)
      X=X+Q3
      CALL DERIV
      DO 9 M=1,8
9     Y(M)=Q6*D(M)+0.5*DT(M,1)+YT(M)
      SKIP=1.0
      CALL DERIV
      SKIP=0.0
      DO 10 M=1,8
      DT(M,2)=Q15*D(M)
10    Y(M)=0.375*DT(M,1)+0.25*DT(M,2)+YT(M)
      X=XT+Q2
      CALL DERIV
      DO 11 M=1,8
      DT(M,1)=Q4*D(M)+DT(M,1)
11    Y(M)=1.5*DT(M,1)-DT(M,2)+YT(M)
      X=XT+Q
      CALL DERIV
      DOUBLE=2.0
      DO 14 M=1,8
      DT(M,3)=Q6*D(M)
      TEST=ABS(DT(M,1)-DT(M,2)-DT(M,3))
      IF (TEST .GT. ERHIGH) GOTO 18
      IF (TEST .GE. ERLOW) DOUBLE = 1.0
14    CONTINUE
      DO 15 M=1,8
15    Y(M)=0.5*DT(M,1)+DT(M,3)+YT(M)
      Q=DOUBLE*Q
      IF (LAST .NE. 0) GOTO 20
16    IF ( (X1-X) .GT. Q ) GOTO 3
      GOTO 2
18    Q = Q * .5
      LAST=0
      DO 19 M=1,8
19    Y(M)=YT(M)
      X=XT
      GO TO 6
20    IF (Q .LT. Q1) Q = Q1
22    RETURN
      END
      SUBROUTINE NEWTON
      COMMON/ANCNEW/NEW,ZR,ZI,QR(7),QI(7),KRASH
      KONVRG=0
      DO 6 KOUNT =1,70
      XR=0.0
      XI=0.0
      YR=0.0
      YI=0.0
      DO 3 J=1,NEW
      TR=ZR*YR-ZI*YI+XR
      YI=ZR*YI+ZI*YR+XI
      YR=TR
      M=NEW+1-J
      TR=ZR*XR-ZI*XI+QR(M)
      TI=ZR*XI+ZI*XR+QI(M)
      IF (KONVRG .EQ. 0) GOTO 2
1     QR(M)=XR
      QI(M)=XI
2     XR=TR
3     XI=TI
      IF (KONVRG .NE. 0) GOTO 7
4     F=1.0/(YR**2+YI**2)
      TR=F*(XR*YR+XI*YI)
      TI=F*(XI*YR-XR*YI)
      ZR=ZR-TR
      ZI=ZI-TI
C      IF((TR**2+TI**2)/(ZR**2+ZI**2)-0.1E-11)5,5,6
      IF ( (TR**2+TI**2)/(ZR**2+ZI**2) .LE. 0.1E-11 ) KONVRG = 1
C5     KONVRG=1
6     CONTINUE
      KRASH=-70
      GO TO 10
7     IF(ABS(ZR) .GT. 100000.0*ABS(ZI)) KRASH = NEW
C8     KRASH=NEW
9     NEW=NEW-1
10    RETURN
      END
      SUBROUTINE HALFTN(N,ARRAY,MARK,BLACK,WHITE,JUMP,IFILE)
C  THIS SUBROUTINE OUTPUTS ONE LINE OF AN IMAGE.
C IBM VERSION FOR SUPERPROG.FROM MIKE O'K IN U.K.,1980 (J.C.H.S)
      INTEGER FILE
      DIMENSION ARRAY(157),JRAY(157)
      INTEGER RANX1, RANX2
      ISCLE = 36
      MW=7
C  JUMP TO SET RANGE, WRITE SCALE, OUTPUT IMAGE LINE
      IF(JUMP)190,160,10
C  OUTPUT IMAGE LINE
C  CHECK VALUE OF N
   10 IF(N*(N-157))40,20,20
   20 WRITE(MW,1200)N
      GO TO 230
C  BRANCH TO OUTPUT TO LA100, PERITEC, OR FILE
40    IF (IFILE) 42,230,900 
42    WRITE (6,41)
41    FORMAT (' ','-')
C  FIND STARTING AND FINISHING COLUMNS TO CENTRE IMAGE LINE
      IF (JS .LE. 1) GOTO 31
C      DO 31 J = 1, (JS - 1)
C      CALL PUTBLK (0)
31    CONTINUE
C  LOOP OVER INTENSITIES IN LINE
      DO 130 J=JS,JF
      K=J-JS+1
C  SCALE INTENSITIES TO MATCH SYMBOL TABLE
      IOUT= ISCLE - ((ARRAY(K)-DARK)*SCALE)  
C  CHECK FOR VALUES OUT-OF-RANGE
      IF(IOUT.LT.0)IOUT=0
      IF(IOUT.GT.ISCLE)IOUT=ISCLE
      CALL PUTBLK( IOUT )
130   CONTINUE
160   CONTINUE
      GO TO 230
C  OUTPUT TO FILE
900   IF (JS .LE. 1) GOTO 931
      DO 931 J = 1, (JS - 1)
      JRAY(J) = 0
931    CONTINUE
C  LOOP OVER INTENSITIES IN LINE
      DO 9130 J=JS,JF
      K=J-JS+1
C  SCALE INTENSITIES TO MATCH SYMBOL TABLE
      IOUT= ISCLE - ((ARRAY(K)-DARK)*SCALE)  
C  CHECK FOR VALUES OUT-OF-RANGE
      IF(IOUT.LT.0)IOUT=0
      IF(IOUT.GT.ISCLE)IOUT=ISCLE
      JRAY(J) = IOUT
9130   CONTINUE
      WRITE(IFILE,9140) (JRAY(J),J=1,55)
      WRITE(IFILE,9140) (JRAY(J),J=56,110)
      WRITE(IFILE,9140) (JRAY(J),J=111,JF)
9140  FORMAT(' ',55I2)
      GO TO 230
C  SET RANGE
  190 DARK=BLACK
      BRITE=WHITE
      RANGE =BRITE-DARK
      IF(ABS(RANGE).GT.1E-08) GO TO 50
      WRITE(MW,1500) RANGE
      GO TO 230
50    SCALE = FLOAT(ISCLE)/RANGE
      JS=(158-N)/2
      JF=(156+N)/2
      RANX1 = 0
      RANX2 = 0
      DO 200 NN = 1, 200
      I = IRANDM(1,2)
200   CONTINUE
      IF (IFILE .GT. 0) WRITE(IFILE,201) (JF - 110)
201   FORMAT (' ','035555',I2)
 1200 FORMAT(23H N IS WRONG, AND EQUALS,I15,/)
 1500 FORMAT(39H0INTENSITY RANGE LT 1.0E-8 AND EQUAL TO,E8.1,/)
  230 RETURN
      END
      SUBROUTINE HALFTN(N,ARRAY,MARK,BLACK,WHITE,JUMP,IFILE)
C  THIS SUBROUTINE OUTPUTS ONE LINE OF AN IMAGE.
C IBM VERSION FOR SUPERPROG.FROM MIKE O'K IN U.K.,1980 (J.C.H.S)
C MODIFIED FOR USE WITH ASU MULTISLICE AND LA100 GRAPHICS MODE
C BY R. HOLTON, ARGONNE NATIONAL LAB, MST, 1984
      INTEGER FILE
      DIMENSION ARRAY(157),JRAY(157)
      INTEGER RANX1, RANX2
      COMMON /TV/ ITV,NCOLS,XLPI,XCPI,ISCLE
      DATA MW/7/
C  JUMP TO SET RANGE, WRITE SCALE, OUTPUT IMAGE LINE
      IF(JUMP)190,160,10
C  OUTPUT IMAGE LINE
C  CHECK VALUE OF N
   10 IF(N*(N-MAXWD))900,20,20
   20 WRITE(MW,1200)N
      GO TO 230
C  OUTPUT TO FILE
900   IF (JS .LE. 1) GOTO 931
      DO 931 J = 1, (JS - 1)
      JRAY(J) = ISCLE
931    CONTINUE
C  LOOP OVER INTENSITIES IN LINE
      DO 9130 J=JS,JF
      K=J-JS+1
C  SCALE INTENSITIES TO MATCH SYMBOL TABLE
      IOUT= ISCLE - ((ARRAY(K)-DARK)*SCALE)  
C  CHECK FOR VALUES OUT-OF-RANGE
      IF(IOUT.LT.0)IOUT=0
      IF(IOUT.GT.66)IOUT=ISCLE
      JRAY(J) = IOUT
9130   CONTINUE
      WRITE(IFILE,9140) (JRAY(J),J=1,55)
      WRITE(IFILE,9140) (JRAY(J),J=56,110)
      WRITE(IFILE,9140) (JRAY(J),J=111,JF)
9140  FORMAT(' ',55I2)
      GO TO 230
C  SET RANGE
  190 DARK=BLACK
      BRITE=WHITE
      RANGE =BRITE-DARK
      IF(ABS(RANGE).GT.1E-08) GO TO 50
      WRITE(MW,1500) RANGE
      GO TO 230
50    SCALE = FLOAT(ISCLE)/RANGE
      MAXWD = NCOLS + 2
      JS=(MAXWD+1-N)/2
      JF=(MAXWD-1+N)/2
      RANX1 = 0
      RANX2 = 0
      DO 200 NN = 1, 200
      I = IRANDM(1,2)
200   CONTINUE
      IF (IFILE .GT. 0) WRITE(IFILE,201) (JF - 110)
201   FORMAT (' ','035555',I2)
 1200 FORMAT(23H N IS WRONG, AND EQUALS,I15,/)
 1500 FORMAT(39H0INTENSITY RANGE LT 1.0E-8 AND EQUAL TO,E8.1,/)
  230 RETURN
      END
	SUBROUTINE SETBLK (TONE, BLOCK)
C
	INTEGER TONE
	BYTE BLOCK(6)
C
C	THIS SUBROUTINE RETURNS ONE 6 BY 6 BLOCK IN LA100 FORMAT 
C	FOR A HALFTONE IMAGE.  TONE IS A INTEGER FROM 0 TO 36 (THUS
C	THERE ARE 37 DIFFERENT TONES TO THE GREYSCALE), WHERE 0 IS  
C	WHITE AND 36 IS BLACK.  THE ROUTINE RANDOMLY SELECTS WHICH
C	POINTS WITHIN THE BLOCK ARE TO BE PRINTED.
C
C	AUTHOR:  R. HOLTON,  ARGONNE NATIONAL LABORATORY, 8-84
C
	INTEGER
     &	CLEAR,		! ORIGINAL SETTING OF ALL POINTS IN THE BLOCK
     &	FILL,		! SETTING OF THOSE POINTS RANDOMLY SELECTED
     &	POINTS,		! NUMBER OF CLEAR POINTS LEFT
     &	HALF,		! HALF OF THE TOTAL NUMBER OF POINTS
     &	SELECT,		! WHICH CLEAR POINT TO FILL
     &	LOC,		! POINT TO CHECK IF FILL OR CLEAR
     &	BOX(6,6),	! THE BLOCK TO BE OUTPUT (6 COLUMNS, 6 ROWS)
     &	ARRAY(36),	! LINEAR ARRAY EQUIVALENCED WITH BOX
     &	NCOLS,		! NUMBER OF COLUMNS (CURRENTLY 6)
     &	NROWS,		! NUMBER OF ROWS (CURRENTLY 6...LIMITED BY LA100)
     &	FINISH		! TERMINATING VALUE FOR LOOP
C
	EQUIVALENCE (ARRAY,BOX)
C
	NCOLS = 6
	NROWS = 6
	POINTS = NCOLS * NROWS
	HALF = POINTS / 2
C
	IF (TONE .GE. 0) GOTO 20
		TONE = 0
		GOTO 40
C
20	IF (TONE .LE. POINTS) GOTO 30
		TONE = POINTS
		GOTO 50
C
30	IF (TONE .GE. HALF) GOTO 50
40		CLEAR = 0
		FILL = 1
		FINISH = (POINTS - TONE) + 1
		GOTO 60
C	ELSE
50		CLEAR = 1
		FILL = 0
		FINISH = TONE + 1
C
60	DO 70  J = 1,POINTS
		ARRAY(J) = CLEAR
70		CONTINUE
C
	IF (FINISH .GT. POINTS) GOTO 140
C
	DO 100  J = POINTS ,FINISH, -1
		SELECT = IRANDM(1,J)
		LOC = 0
		DO 90  K = 1, SELECT
80			LOC = LOC + 1
			IF (ARRAY(LOC) .EQ. FILL) GOTO 80
90			CONTINUE
		ARRAY(LOC) = FILL
100		CONTINUE
C
	DO 130  J = 1, NCOLS
		BLOCK(J) = 63
		DO 110  K = 1, NROWS
			BLOCK(J) = BLOCK(J) + BOX(J,K)*2**(K-1)
110			CONTINUE
C
130		CONTINUE
C
140	RETURN
	END
	SUBROUTINE PUTBLK (TONE)
C
	INTEGER TONE
C
C	THIS SUBROUTINE OUTPUTS TO LOGICAL DEVICE PRNT ONE BLOCK
C	FOR A HALFTONE IMAGE.  TONE IS A INTEGER FROM 0 TO 36 (THUS
C	THERE ARE 37 DIFFERENT TONES TO THE GREYSCALE), WHERE 0 IS  
C	WHITE AND 36 IS BLACK.  THE ROUTINE RANDOMLY SELECTS WHICH
C	POINTS WITHIN THE BLOCK ARE TO BE PRINTED.
C	PRNT SHOULD EITHER BE AN LA100 PRINTER, OR A SEQUENTIAL ACCESS
C	FILE THAT WILL LATER BE OUTPUT TO AN LA100 PRINTER.
C
C	AUTHOR:  R. HOLTON,  ARGONNE NATIONAL LABORATORY, 8-84
C
	INTEGER
     &	CLEAR,		!ORIGINAL SETTING OF ALL POINTS IN THE BLOCK
     &	FILL,		!SETTING OF THOSE POINTS RANDOMLY SELECTED
     &	POINTS,		!NUMBER OF CLEAR POINTS LEFT
     &	HALF,		!HALF OF THE TOTAL NUMBER OF POINTS
     &	SELECT,		!WHICH CLEAR POINT TO FILL
     &	LOC,		!POINT TO CHECK IF FILL OR CLEAR
     &	BOX(6,6),	!THE BLOCK TO BE OUTPUT (6 COLUMNS, 6 ROWS)
     &	ARRAY(36),	!LINEAR ARRAY EQUIVALENCED WITH BOX
     &	NCOLS,		!NUMBER OF COLUMNS (CURRENTLY 6)
     &	NROWS,		!NUMBER OF ROWS (CURRENTLY 6...LIMITED BY LA100)
     &	FINISH,		!TERMINATING VALUE FOR LOOP
     &	PRNT,		!LOGICAL UNIT FOR OUTPUT...CURRENTLY LOGICAL UNIT 6
     &	WIRES		!ONE VERTICAL COLUMN IN OUTPUT FORMAT
C
	EQUIVALENCE (ARRAY,BOX)
C
	PRNT = 6
	NCOLS = 6
	NROWS = 6
	POINTS = NROWS * NCOLS
	HALF = POINTS / 2
C
	IF (TONE .GT. 0) GOTO 20

	THI	WRITE(PRNT,10)
10		FORMAT(' ',6('?'))
		GOTO 140
C
20	IF (TONE .LT. POINTS) GOTO 40
		WRITE(PRNT,30)
30		FORMAT(' ',6('~')) 
		GOTO 140
C
40	IF (TONE .GE. HALF) GOTO 50
		CLEAR = 0
		FILL = 1
		FINISH = (POINTS - TONE) + 1
		GOTO 60
C	ELSE
C
	DO 100  J = POINTS ,FINISH, -1
		SELECT = IRANDM(1,J)
		LOC = 0
		DO 90  K = 1, SELECT
80			LOC = LOC + 1
			IF (ARRAY(LOC) .EQ. FILL) GOTO 80
90			CONTINUE
		ARRAY(LOC) = FILL
100		CONTINUE
C
	DO 130  J = 1, NCOLS
		WIRES = 63
		DO 110  K = 1, NROWS
			WIRES = WIRES + BOX(J,K)*2**(K-1)
110			CONTINUE
C
		WRITE(PRNT,120) WIRES
120		FORMAT(' ',A1)
130		CONTINUE
C
140	RETURN
	END
	
	
	=============================================================
	
	 THIS THE OUTPUT PROGRAM WHICH WILL PRINT THE RESULTING DATA
	FILE OUT ONTO AN LA100 OR LA50 PRINTER
	
	
	        PROGRAM IMOUT
C
        INTEGER NCHAR(9), BUFFER(132)
C
        DATA NESC /27/
C
        WRITE(7,10)
10      FORMAT(' Enter the name of the file to be displayed: ',$)
        CALL ASSIGN (3,DUMMY,-1,'OLD',80)
        DO 15 J = 1,150
        K = IRANDM(1,2)
15      CONTINUE
        READ(3,20) NRECS,(NCHAR(J),J=1,NRECS)
20      FORMAT (1X,10I2)
C       WRITE(7,22) NRECS, (NCHAR(J),J=1,NRECS)
C22     FORMAT(' ',10(2X,I2))
        WRITE(6,25) NESC
25      FORMAT(' ',A1,'P1q')
30      WRITE(6,35)
35      FORMAT (' ','-')
        DO 60 K = 1,NRECS
        READ(3,40,END=90) (BUFFER(J),J=1,NCHAR(K))
40      FORMAT(1X,132I2)
        DO 50 L = 1,NCHAR(K)
        CALL PUTBLK(BUFFER(L))
50      CONTINUE
60      CONTINUE
        GOTO 30
C
90      WRITE (6,100) NESC
100     FORMAT (' ','-',A1,'\')
        CALL EXIT
        END

        SUBROUTINE PUTBLK (TONE)
C
        INTEGER TONE
C
C       THIS SUBROUTINE OUTPUTS TO LOGICAL DEVICE PRNT ONE BLOCK
C       FOR A HALFTONE IMAGE.  TONE IS A INTEGER FROM 0 TO 36 (THUS
C       THERE ARE 37 DIFFERENT TONES TO THE GREYSCALE), WHERE 0 IS  
C       WHITE AND 36 IS BLACK.  THE ROUTINE RANDOMLY SELECTS WHICH
C       POINTS WITHIN THE BLOCK ARE TO BE PRINTED.
C       PRNT SHOULD EITHER BE AN LA100 PRINTER, OR A SEQUENTIAL ACCESS
C       FILE THAT WILL LATER BE OUTPUT TO AN LA100 PRINTER.
C
C       AUTHOR:  R. HOLTON,  ARGONNE NATIONAL LABORATORY, 8-84
C
        INTEGER
     &  CLEAR,          !ORIGINAL SETTING OF ALL POINTS IN THE BLOCK
     &  FILL,           !SETTING OF THOSE POINTS RANDOMLY SELECTED
     &  POINTS,         !NUMBER OF CLEAR POINTS LEFT
     &  HALF,           !HALF OF THE TOTAL NUMBER OF POINTS
     &  SELECT,         !WHICH CLEAR POINT TO FILL
     &  LOC,            !POINT TO CHECK IF FILL OR CLEAR
     &  BOX(6,6),       !THE BLOCK TO BE OUTPUT (6 COLUMNS, 6 ROWS)
     &  ARRAY(36),      !LINEAR ARRAY EQUIVALENCED WITH BOX
     &  NCOLS,          !NUMBER OF COLUMNS (CURRENTLY 6)
     &  NROWS,          !NUMBER OF ROWS (CURRENTLY 6...LIMITED BY LA100)
     &  FINISH,         !TERMINATING VALUE FOR LOOP
     &  PRNT,           !LOGICAL UNIT FOR OUTPUT...CURRENTLY LOGICAL UNIT 6
     &  WIRES           !ONE VERTICAL COLUMN IN OUTPUT FORMAT
C
        EQUIVALENCE (ARRAY,BOX)
C
        PRNT = 6
        NCOLS = 6
        NROWS = 6
        POINTS = NROWS * NCOLS
        HALF = POINTS / 2
C
        IF (TONE .GT. 0) GOTO 20
                WRITE(PRNT,10)
10              FORMAT(' ',6('?'))
                GOTO 140
C
20      IF (TONE .LT. POINTS) GOTO 40
                WRITE(PRNT,30)
30              FORMAT(' ',6('~')) 
                GOTO 140
C
40      IF (TONE .GE. HALF) GOTO 50
                CLEAR = 0
                FILL = 1
                FINISH = (POINTS - TONE) + 1
                GOTO 60
C       ELSE
50              CLEAR = 1
                FILL = 0
                FINISH = TONE + 1
C
60      DO 70  J = 1,POINTS
                ARRAY(J) = CLEAR
70              CONTINUE
C
        DO 100  J = POINTS ,FINISH, -1
                SELECT = IRANDM(1,J)
                LOC = 0
                DO 90  K = 1, SELECT
80                      LOC = LOC + 1
                        IF (ARRAY(LOC) .EQ. FILL) GOTO 80
90                      CONTINUE
                ARRAY(LOC) = FILL
100             CONTINUE
C
        DO 130  J = 1, NCOLS
                WIRES = 63
                DO 110  K = 1, NROWS
                        WIRES = WIRES + BOX(J,K)*2**(K-1)
110                     CONTINUE
C
                WRITE(PRNT,120) WIRES
120             FORMAT(' ',A1)
130             CONTINUE
C
140     RETURN
        END

        FUNCTION IRANDM (LOW, HIGH)
C
        INTEGER LOW, HIGH
C
C       THIS FUNCTION RETURNS A RANDOM INTEGER UNIFORMLY DISTRIBUTED
C       FROM LOW TO HIGH.
C
C       NOTE THAT THIS PROGRAM MAKE A CALL TO A RANDOM NUMBER GENERATOR.
C
        INTEGER RANX1, RANX2
        COMMON /IRANX /RANX1, RANX2
C       
        IF (LOW .GE. HIGH) GOTO 99
C
        IRANDM = IFIX((HIGH - LOW + 1) * RAN(RANX1,RANX2)) + LOW
        RETURN
C
99      IRANDM = LOW
        RETURN
        END
        
        
        ==============END OF CODE =================