'$DYNAMIC
'************************************************
'***  MXD-HIST(ory)                           ***
'***                                          ***
'***  Plot T-P-E-D-A-B-C... - TIME            ***
'***                                          ***
'***  1993-Dec-27  version (kats)             ***
'***  1994-Feb-01  version for LIPS3 (kats)   ***
'***  1996-May-10  5000 steps (kats)          ***
'***       Oct-14  f-basic95                  ***
'***  2000-Oct-18  width function (kats)      ***
'***  2005-Aug-14  FreeBASIC version          ***
'***  2005-Oct-31  Average and standard devi. ***
'***  2008-Nov-22  Minor change               ***
'***  2008-Nov-22  New format of file09v.dat  ***
'***  2009-Jan-27  bald lines                 ***
'************************************************
'
 deflng I-N, X-X
 defdbl a-h,o-W,y-z
'
 ' screen 19 :  xscreen= 800: yscreen=600
   screen 20 :  xscreen=1024: yscreen=768 : lcx=xscreen/8 : lcy=yscreen/16
                                            MAXLINE = LCY

  LVA = 44: NLIM = 50001:   '  nlim = 51
                            '  NLIM = 101
  dim D1(NLIM),D4(NLIM),D7(NLIM), S$(LVA), WA(LVA), DMAX(LVA)
  dim D2(NLIM),D5(NLIM),D8(NLIM), WD(LVA), DMIN(LVA)
  dim D3(NLIM),D6(NLIM),D9(NLIM), FACT(LVA), IVAL(9), JVAL(LVA)
  dim EGI(LVA),D(NLIM), kval(9,11)
  dim atom$(10), nion(10), ions(2,10), spres(6),al(6), kcol(12)

  dim ncompo,npstep,npstepi,no,itps,ivl,njob1,njob2
  dim dtime
  dim title$,ss$,s2$,s3$

  for I= 1 to LVA: EGI(I) =0: next
'
   s$(1) = "T/K    ":   s$(17) = "D/gcm-3"
   s$(2) = "P/GPa  ":   s$(18) = "Vm/cm3 "
   s$(3) = "Pxx/GPa":   s$(19) = "a/A    ":   s$(22) = "alpha  "
   s$(4) = "Pyy/GPa":   s$(20) = "b/A    ":   s$(23) = "beta   "
   s$(5) = "Pzz/GPa":   s$(21) = "c/A    ":   s$(24) = "gamma  "
   s$(6) = "Pxy/GPa":
   s$(7) = "Pxz/GPa":   s$(25) = "T(1)/K ":   s$(35) = "Msd(1) "
   s$(8) = "Pyz/GPa":   s$(26) = "T(2)/K ":   s$(36) = "Msd(2) "
   s$(9) = "U(C)/kJ":   s$(27) = "T(3)/K ":   s$(37) = "Msd(3) "
  s$(10) = "U(s)/kJ":   s$(28) = "T(4)/K ":   s$(38) = "Msd(4) "
  s$(11) = "U(3)/kJ":   s$(29) = "T(5)/K ":   s$(39) = "Msd(5) "
  s$(12) = "U(t)/kJ":   s$(30) = "T(6)/K ":   s$(40) = "Msd(6) "
  s$(13) = "E(k)/kJ":   s$(31) = "T(7)/K ":   s$(41) = "Msd(7) "
  s$(14) = "E(t)/kJ":   s$(32) = "T(8)/K ":   s$(42) = "Msd(8) "
  s$(15) = "PV/kJ  ":   s$(33) = "T(9)/K ":   s$(43) = "Msd(9) "
  s$(16) = "H/kJ   ":   s$(34) = "T(10)/K":   s$(44) = "Msd(10)"
'
   kcol(1)=1: kcol(2)= 2: kcol(3)= 3: kcol( 4)= 4: kcol( 5)= 5: kcol( 6)= 6
   kcol(7)=9: kcol(8)=10: kcol(9)=11: kcol(10)=12: kcol(11)=13: kcol(12)=34
'
                  FOR j = 1 TO lva: jval(j) = 0: NEXT j
'                 ================================================
   NVAL = 5:      IVAL(1) = 1:  JVAL(1) = 1  ' 1: Temperature (1)
                  IVAL(2) = 2:  JVAL(2) = 1  ' 2: Pressure    (2)
                  IVAL(3) = 14: JVAL(14) = 1 ' 4: E(total)    (14)
                  IVAL(4) = 17: JVAL(17) = 1 ' 5: Density     (17)
                  IVAL(5) = 19: JVAL(19) = 1 ' 6: A           (19)
                  IVAL(6) = 20: JVAL(20) = 1 ' 7: B           (20)
                  IVAL(7) = 21: JVAL(21) = 1 ' 8: C           (21)
'                 ================================================
'
kval(1,1)= 1:kval(2,1)= 2:kval(3,1)=14:kval(4,1)=17:kval(5,1)=19:kval(6,1)=20:kval(7,1)=21:kval(8,1)=0:kval(9,1)=0
kval(1,2)= 2:kval(2,2)= 3:kval(3,2)= 4:kval(4,2)= 5:kval(5,2)= 6:kval(6,2)= 7:kval(7,2)= 8:kval(8,2)=0:kval(9,2)=0
kval(1,3)=17:kval(2,3)=19:kval(3,3)=20:kval(4,3)=21:kval(5,3)=22:kval(6,3)=23:kval(7,3)=24:kval(8,3)=0:kval(9,3)=0
'
   GOSUB Menu
'
      YINTV=6
'
 '--------------------------------- Read data from file07.dat and file09v.dat
 '------------------------------------------------------------ Data directory
     d$=""
     '
     CLS 0
     GOSUB readf07 '------------------------------------ Read from file07.dat
     for I=1 to NCOMPO
        S$(24+I) = "T("   + left$(ATOM$(I),2) + ")/K"
        S$(34+I) = "Msd(" + left$(ATOM$(I),2) + ")"
     next I
     '
     PRINT : COLOR 14 
     PRINT "No. of positions recorded:"; npstep&, "   Limit of plot:"; nlim
     PRINT : ISTEP = 1
     input "Interval of data plotted (1-m, default:1) "; ISTEP
     PRINT : COLOR 7 : IF ISTEP < 1 THEN  ISTEP = 1
     DT = DTIME * npstepi * ISTEP
     '
     GOSUB DATAREAD09v  '>---------------------> Go to subroutine DATAREAD09v
     COLOR 14 
     PRINT : PRINT "Push [return] key for plotting": COLOR 7 
L1:  a$ = INKEY$: IF a$ = "" THEN  GOTO L1
'
     GOSUB Scale '-------------------------------------------- Calculate scale
'
LG:  GOSUB GRAPH  '>----------------------------------> Go to subroutine GRAPH
'
     locate MAXLINE,1: color 0: print "Type ";
                   color 1: print "[Q]"; : color 2: print " to quit ";
                   color 1: print "[S]"; : color 2: print " to change scale";
L3:  A$ = inkey$: if A$ = "" then  goto L3
     if A$ = "S" or A$ = "s" then  gosub SET_SCALE : goto LG
     if A$ = "Q" or A$ = "q" then  cls 0: end
     IF a$ = "E" OR a$ = "e" THEN  CLS 0: END
     GOTO L3
     end
'
'
'================================================================ Subroutines
MENU: '------------------------------------------------------------------ Menu
     CLS 0
     itm=0
     for I = 1 to LVA
        IX = (i-1) mod 6 
        IY = (I-1)\6+2
        color 15
        if JVAL(I) = 1 then  color 14 : itm=itm+1
        locate IY,IX*16+2
        print using "## &     &"; I; S$(I);
     next I: color 15
           LOCATE 11,50 : COLOR 14
           PRINT USING "Number of selected items is ##";itm;
     '
     color 15
     LOCATE 13,2: PRINT "The yellow items will be plotted. "
     locate 14,2: print "       (Maximum number of items is 7)"
     LOCATE 16,6: PRINT "1-44   : Type number (1-44) to switch  [ON] <-> [OFF]"
     locate 18,6: print "g,G    : group selection"
     LOCATE 20,6: PRINT "r,R    : read file09v and plot"
     LOCATE 22,6: PRINT "q,Q    : quit to system "

     for k=1 to 4 ' ---------------------------------------- Groups
         l=26 : locate l+k, 5 : print using "##:";k;
         for i=1 to 7: locate l+k,i*12 : print using "##:";kval(i,k); 
                                        print s$(kval(i,k)); : next i
         print
     next k

A:  locate 24,1 : color 15
    color 12: input "Type-in ? ",a$: color 15 : if A$="" then  goto A
              jj=val(a$)
              if jj >0 and jj < 45 then
                     if jval(jj)=1 then
                            jval(jj)=0
                     else
                            jval(jj)=1
                     end if
              END IF
     if a$ = "g" or a$ = "G" then '---------------------- Group select
           locate 26+5,5 : color 12 : input "Typein ?";n
           for i=1 to 44 : jval(i)=0 : next i
           for i=1 to 9 : j=kval(i,n) : if j>0 then jval(j)=1 
           next : goto menu
     end if
     IF a$ = "q" OR a$ = "Q" THEN STOP
     IF a$ = "r" or a$ = "R" THEN GOTO sel
     goto menu

sel: nval = 0
     FOR i = 1 TO lva
        IF jval(i) = 1 THEN
              nval = nval + 1
              ival(nval) = i
        END IF
     NEXT i: LOCATE 21,1
     '
     return
     end
'
'
GRAPH: '--------------------------------------------------------------- graph
       cls 0
       line (0,0)-(xscreen,yscreen),15,bf
       color 0,15: locate 1,12 : print TITLE$;
       NPLT = 0: NOPLOTS = NVAL
       ICOL = 1: YL = NVAL*YINTV*16+11
                 XO = 76: YYOO=18
                 XL = 1000-XO - 46
                 FX = XL / NO
       YINCRE = YINTV*16: YO = YINCRE / 2 + 23
       line (XO,   YYOO  )-(XO+XL,   Yyoo+yl  ), 0, b    
       line (XO-1, YYOO+1)-(XO+XL+1, Yyoo+yl-1), 0, b    
       '
       TPS = DT*NO /1E-12
       M = 0: X = XO: fPS = xl/tps
       JS =1 : if TPS >  9 then JS =  2
               if TPS > 19 then JS =  5
               if TPS > 49 then JS = 10
               if TPS > 99 then JS = 20
               if TPS >199 then JS = 50
               if TPS >501 then JS = 100
       ky=(yyoo+yl)/16+1
       for IT = 0 to TPS step JS '------------------------------- Time marker
           IX = fPS*IT +XO
           if IX<=xscreen then L=6
           if ITPS <21 and IT mod 5 = 0 then L = 9
           if IT mod 10 =0 then L=7
                  locate ky,(IX -8)/8-2
                  print using "####";IT;
           line (IX,YYOO)-(IX+1,YYOO+L), 1,b
           line (IX,Yyoo+yl)-(IX+1,Yyoo+yl-l), 1,b
       next IT
       '
       color 0,15:' locate ky,7 : print "0";
                  ' locate ky,71: print using "####"; M - 1;
                    locate ky+1,55  : print "Time / ps";
'
'              =====================
               for IIII = 1 to NVAL
                     IVL = IVAL(IIII)
                     gosub PLOT
               next IIII
'              =====================
   return
   end
'
'
PLOT: '----------------------------------------------------------------- plot
      nplt = nplt + 1
      FACTJ = FACT(IVL) :  FACTK = FACT(IVL) / 12 
      if IIII = 1 then for I = 1 to NO: D(I) = D1(I): next
      if IIII = 2 then for I = 1 to NO: D(I) = D2(I): next
      if IIII = 3 then for I = 1 to NO: D(I) = D3(I): next
      if IIII = 4 then for I = 1 to NO: D(I) = D4(I): next
      if IIII = 5 then for I = 1 to NO: D(I) = D5(I): next
      if IIII = 6 then for I = 1 to NO: D(I) = D6(I): next
      if IIII = 7 then for I = 1 to NO: D(I) = D7(I): next
      if IIII = 8 then for I = 1 to NO: D(I) = D8(I): next
      if IIII = 9 then for I = 1 to NO: D(I) = D9(I): next
      ICOL = kcol(nplt): if ICOL > 13 then ICOL = 1
      color ICOL,15

      line (XO, YO)-(XO + XL, YO), 1
       for IT = 0 to TPS step JS '------------------------------- Time marker
           IX = fPS*IT +XO
           if IX<=xscreen then L=5
           if ITPS <21 and IT mod 5 = 0 then L = 8
           line (IX,YO-L)-(IX,YO+L), 1
       next IT
  
      DAVJ = WA(IVL)
      if IVL>=2 and IVL<=8 then
            if S$(2)="P/MPa  " then DAVJ=DAVJ/1000: FACTJ=FACTJ*1000
      end if
      I0 = 0:       
      IX = XO: IY = YO - (D(1) - DAVJ) / FACTK
      for I = 2 to NO
          A = I*FX
          JX = XO + A
          JY = YO - (D(I) - DAVJ) / FACTK
          if iy>=yyoo and IY<=YSCREEN and jy>=yyoo and JY<=YSCREEN then
                line (IX, IY)-(JX, JY), iCOL
          end if
          IX = JX: IY = JY:' i0 = d(i)
      NEXT i
      line (XO, YO)-(XO + XL, YO), 1
       for IT = 0 to TPS step JS '------------------------------- Time marker
           IX = fPS*IT +XO
           if IX<=xscreen then L=5
           if ITPS <21 and IT mod 5 = 0 then L = 8
           line (IX,YO-L)-(IX,YO+L), 1
       next IT
      line (XO+XL, YO)-(XO+XL+5, YO), 1
      line (XO+XL, YO+12)-(XO+XL-7, YO+12), 1
      line (XO+XL, YO-12)-(XO+XL-7, YO-12), 1
      LL = NPLT *yintv-2
      SS$=S$(IVL): S2$=left$(SS$,2)
      locate LL,2   : print SS$;
      locate LL+1,2 : AD = WA(IVL)
           if S2$="T(" then
                  print using "####.###";AD
           elseif S2$="P/" or S2$="Px" or S2$="Py" or S3$="Pz" then
                  if S$(2)="P/MPa  " then
                      print using "####.###";WA(IVL);
                  else
                      if AD>=100            then print using "###.####";WA(IVL);
                      if AD<100 and AD>= 10 then print using "##.#####";WA(IVL);
                      if AD< 10 and AD>=  0 then print using "#.######";WA(IVL);
                      if AD<  0 and AD> -10 then print using "##.#####";WA(IVL);
                      if AD<-10 and AD>-100 then print using "###.####";WA(IVL);
                      if AD<-100            then print using "####.###";WA(IVL);
                  end if
           else
                              if Ad >=100000 then print using "########"; WA(IVL);
              if AD < 100000 and AD >= 10000 then print using "######.#"; WA(IVL);
              if AD <  10000 and AD >=  1000 then print using "#####.##"; WA(IVL);
              if AD <   1000 and AD >=   100 then print using "####.###"; WA(IVL);
              if AD <    100 and AD >=    10 then print using "###.####"; WA(IVL);
              if AD <     10 and AD >=     1 then print using "##.#####"; WA(IVL);
              if AD <      1 and AD >=     0 then print using "#.######"; WA(IVL);
              if AD <    0.0 and AD >    -10 then print using "##.#####"; WA(IVL);
              if AD <=   -10 and AD >   -100 then print using "###.####"; WA(IVL);
              if AD <=  -100 and AD >  -1000 then print using "####.###"; WA(IVL);
              if AD <= -1000 and AD > -10000 then print using "#####.##"; WA(IVL);
              if AD <=-10000 and AD >-100000 then print using "######.#"; WA(IVL);
              if AD <-100000                 then print using "########"; WA(IVL);
           end if
      locate ll,(xo+xl)/8+2
      IF factj > 9.9 THEN  PRINT USING "#######"; factj; : GOTO pl1
      IF factj > .99 THEN  PRINT USING "###.###"; factj; : GOTO pl1
      IF factj > .09 THEN  PRINT USING "##.####"; factj; : GOTO pl1
      if factj >= .01 then PRINT USING "#.#####"; factj; : goto pl1
      if factj <  .01 then print using "#.#####"; factj; : goto pl1
                           print using "#.#####";factj;
      
pl1: yo = yo + yincre
     RETURN
     end
'
'
READF07: '--------------------------------------------------------------------
         OPEN d$ + "FILE07.DAT" FOR INPUT AS #2: '------- READ FROM FILE07.DAT
         LINE INPUT #2, a0$: title$ = LEFT$(a0$, 60)
                             NHIST = VAL(MID$(a0$, 66, 5))
         LINE INPUT #2, a1$: PRINT a1$
                             NTION   = val(left$(A1$, 7))
                             NCOMPO  = val(mid$(A1$, 8, 3))
                             NSTEP&  = val(mid$(A1$,11,10))
                             NPSTEP& = val(mid$(A1$,51,10))
                             NPSTEPI = NSTEP&*1.01 / NPSTEP&
         LINE INPUT #2, a2$: PRINT a2$
         LINE INPUT #2, a3$: PRINT a3$
         LINE INPUT #2, a4$: PRINT a4$
         LINE INPUT #2, a5$: PRINT a5$
                FOR i = 1 TO ncompo: i6 = 6 * (i - 1) + 1
                    atom$(i) = MID$(a2$, i6 + 2, 6)
                    nion(i) = VAL(MID$(a3$, i6, 6))
                    ions(1, i) = VAL(MID$(a4$, i6, 6))
                    ions(2, i) = VAL(MID$(a5$, i6, 6))
                NEXT i
'
         LINE INPUT #2, a6$: PRINT a6$
                        t = VAL(MID$(a6$, 1, 10))
                        delta = VAL(MID$(a6$, 11, 10))
                        tget = VAL(MID$(a6$, 21, 10))
                        spres(1) = VAL(MID$(a6$, 31, 10))
                        spres(2) = VAL(MID$(a6$, 41, 10))
                        spres(3) = VAL(MID$(a6$, 51, 10))
         LINE INPUT #2, a7$: PRINT a7$
                        DTIME = VAL(MID$(a7$, 1, 10))
                        al(1) = VAL(MID$(a7$, 21, 10))
                        al(2) = VAL(MID$(a7$, 31, 10))
                        al(3) = VAL(MID$(a7$, 41, 10))
         LINE INPUT #2, a8$: PRINT a8$
                        DNSTY = val(mid$(A8$, 1, 10))
                        HTENSOR$ = mid$(A8$, 11, 10)
                    IF Htensor$ = "H-TENSOR  " THEN
                          LINE INPUT #2, a81$
                          LINE INPUT #2, a82$
                          LINE INPUT #2, a83$
                    END IF

         PRINT title$
         PRINT "NTION="; ntion, "NSTEP="; nstep&,
         PRINT "JOB NO.="; njob1; njob2, "NHIST="; NHIST
         FOR i = 1 TO 5: i10 = 10 * (i - 1)
            PRINT "     "; i; atom$(i), nion(i), ions(1, i); ions(2, i)
         NEXT i
         print "T="; T, "DT="; DTIME, "AL(1-3)="; AL(1); AL(2); AL(3)
         PRINT "DELTA="; delta, "TGET ="; tget
         PRINT "DNSTY="; dnsty, "SPRES(1-3)="; spres(1); spres(2); spres(3)
         PRINT
         '
        'PRINT "READING ATOMIC COORDINATES AT "; : ly = CSRLIN: lx = POS(0)
        'FOR i = 1 TO ntion: LINE INPUT #1, a$
        'FOR j = 1 TO 3: p(j, i) = VAL(MID$(a$, j * 9 - 8)): NEXT j
        'FOR j = 1 TO 3: v(j, i) = VAL(MID$(a$, 21 + j * 8)) * .0625: NEXT j
        'FOR j = 1 TO 3: p0(j, i) = VAL(MID$(a$, 45 + j * 9)): NEXT j
        '   LOCATE ly,lx : PRINT USING "#### "; i;
        'NEXT i
         '
        'mhist = INT(NHIST / 4): PRINT
        'FOR i = 1 TO mhist: LINE INPUT #1, a$: PRINT a$
        '  INPUT #1,(7800)  ((JHIST(J,I),J=1,4),I=1,NHIST)
        'NEXT i
         CLOSE #2: PRINT
         '
         RETURN
         end
'
Scale: ' --------------------------------------------- Calculate scale factors
       scl = 1
       fact(1) = INT(wd(1) * scl / 2) * 2
                  IF fact(1) = 0 THEN fact(1) = 1.0:     ' Temperature
       FOR i = 25 TO 34: fact(i) = fact(1): NEXT         ' Temperature of ions
       '
       iii = 0
       FOR i = 1 TO nval
           IF ival(i) >= 3 AND ival(i) <= 8 THEN iii = 1
       NEXT i
                       fact2 = INT(wd(2) * scl * 50) / 50
       IF iii = 1 THEN fact2 = INT(wd(3) * scl * 50) / 50
                    IF fact2 = 0   THEN fact2 = 0.02 :        'Pressure
                    if wd(2)<0.001 then fact2 = 0.001
       PMAX = -99999
       for I = 2 to 8: FACT(I) = FACT2                    'Pressure tensol
           if PMAX < WA(I) then PMAX = WA(I)
       next I
       if PMAX < 0.9 then
             S$(2) = "P/MPa  "
             S$(3) = "Pxx/MPa" :  S$(6) = "Pyz/MPa"
             S$(4) = "Pyy/MPa" :  S$(7) = "Pxz/MPa"
             S$(5) = "Pzz/MPa" :  S$(8) = "Pxy/MPa"
             for I = 2 to 8
                 WA(I) = WA(I)*1000: WD(I)=WD(I)*1000
             next I
       end if
       '
                              wwwww = wd(12)
       IF wd(13) > wwwww THEN wwwww = wd(13)
       IF wd(14) > wwwww THEN wwwww = wd(14)
       fact9 = INT(wwwww * scl * 50) / 50
                        IF fact9 = 0 THEN fact9 = .01
       FACT( 9) = FACT9:                               'U(coulomb)
       fact(10) = fact9:                               'U(short)
       fact(11) = fact9:                               'U(3-body)
       fact(12) = fact9:                               'U(total)
       fact(13) = fact9:                               'U(kinetic)
       fact(14) = fact9:                               'E(total)
       '          if WWWWW > WD(14)*10 then FACT(14) = FACT9 * 10
       fact(15) = fact9:                               'PV
       fact(16) = fact9:                               'H(enthalpy)
                            'fact(14) = fact9 / 2
       '
       fact(17) = INT(wd(17) * scl * 1000) / 1000
                  IF fact(17) = 0 THEN fact(17) = 0.001  :  ' Density
                  if wd(17)<0.001 then fact(17) = 0.0001 
                  if wd(17)<0.0001 then fact(17)=0.00001
       '
       fact(18) = INT(wd(18) * scl * 1000) / 1000
                  IF fact(18) = 0 THEN fact(18) = 0.0001:   ' Molar volume
       '
       WWWW = WD(19): if WWWW < WD(20) then  WWWW = WD(20)
                      if WWWW < WD(21) then  WWWW = WD(21)
       FACT19 = int(WWWW * SCL * 1000) / 1000
                IF fact19 = 0 THEN fact19 = .001:
       fact(19) = fact19                                 'Cell a
       fact(20) = fact19                                 'Cell b
       fact(21) = fact19                                 'Cell c
       fact(22) = INT(wd(22) * scl * 1000) / 1000        'Cell alpha
       fact(23) = fact(22)                               'Cell beta
       fact(24) = fact(22)                               'Cell gamma
       if FACT(22)<0.00001 then
          FACT(23)=int(WD(23)*SCL*1000) /1000
          FACT(22)=FACT(23)
          FACT(24)=FACT(23)
          if FACT(23)<0.00001 then
             FACT(24)=int(WD(24)*SCL*1000) /1000
             FACT(22)=FACT(24)
             FACT(23)=FACT(24)
          end if
       end if
       if FACT(22)<0.00001 then
          FACT(22)=0.001
          FACT(23)=0.001
          FACT(24)=0.001
       end if
       '
       TMAX = 0.0: TMIN = 99999
       FMAX = 0.0: FMIN = 99999
       for IO=1 to NCOMPO
          JJ = IO+24
          TT = int(WD(JJ) * SCL * 1) / 1
          if TT = 0 then TT = 1
          if TMAX < TT then TMAX = TT
          if TMIN > TT then TMIN = TT
          II = IO+34
          F = int(WD(II) * SCL * 100) / 100
          if F = 0 then F = .01
          if FMAX < F then FMAX = F
          if FMIN > F then FMIN = F
       next IO
       TSCALE = TMIN
       if TMAX > TMIN*5    then TSCALE = TMIN *2
       if TMAX > TMIN*10   then TSCALE = TMIN *4
       if TMAX > TMIN*20   then TSCALE = TMIN *8
       if TMAX > TMIN*50   then TSCALE = TMIN *20
       if TMAX > TMIN*100  then TSCALE = TMIN *40
       if TMAX > TMIN*200  then TSCALE = TMIN *80
       if TMAX > TMIN*500  then TSCALE = TMIN *200
       if TMAX > TMIN*1000 then TSCALE = TMIN *400
      
       FSCALE = FMIN
       if FMAX > FMIN*5    then FSCALE = FMIN *2
       if FMAX > FMIN*10   then FSCALE = FMIN *4
       if FMAX > FMIN*20   then FSCALE = FMIN *8
       if FMAX > FMIN*50   then FSCALE = FMIN *20
       if FMAX > FMIN*100  then FSCALE = FMIN *40
       if FMAX > FMIN*200  then FSCALE = FMIN *80
       if FMAX > FMIN*500  then FSCALE = FMIN *200
       if FMAX > FMIN*1000 then FSCALE = FMIN *400
       for IO = 1 to 10:
            J=IO+24: FACT(J)=TSCALE             'Temp  (1-10)
            I=IO+34: FACT(I)=FSCALE: next       'M.s.d.(1-10)
       return
       end
'
'
'=================================================================== St_scale
SET_SCALE: cls: color 7: print
            print       "  0:         quit to graph draw"
            print using "  1: 1       Temperature            ####.#####";FACT(1)
            print using "  2: 2-8     P Px Py Pz Pyz Pxz Pxy ####.#####";FACT(2)
            print using "  3: 9-16    Uc Us U3 Ut Ek Et PV H ####.#####";FACT(9)
            print using "  4: 17      d g/cm3                ####.#####";FACT(17)
            print using "  5: 18      V cm3/mol              ####.#####";FACT(18)
            print using "  6: 19-21   a b c                  ####.#####";FACT(19)
            print using "  7: 22-24   alpfa beta gamma       ####.#####";FACT(22)
            print using "  8: 25-34   Temperature of atoms   ####.#####";FACT(25)
            print using "  9: 35-44   msd of atoms           ####.#####";FACT(35)
            print : input "  No 1-9    ";II
            if II>0 and II<10 then
                    input "  New scale "; FF
                 if II=1 then FACT(1)=FF
                 if II=2 then for I=2 to 8: FACT(I)=FF: next
                 if II=3 then for I=9 to 16: FACT(I)=FF: next
                 if II=4 then FACT(17)=FF
                 if II=5 then FACT(18)=FF
                 if II=6 then for I=19 to 21: FACT(I)=FF: next
                 if II=7 then for I=22 to 24: FACT(I)=FF: next
                 if II=8 then for I=25 to 34: FACT(I)=FF: next
                 if II=9 then for I=35 to 44: FACT(I)=FF: next
                 goto set_scale
            end if
            return
            end
'
'====================================================== Read from file09v.dat
DATAREAD09V:  for I = 1 to LVA: WA(I) = 0: DMAX(I) = -9.9E+19
                                WD(I) = 0: DMIN(I) =  9.9E+19:  next I
              MO& = 0: INTV =10: if NPSTEP&>=  200 then INTV = 20
              NO  = 0:           if NPSTEP&>= 1000 then INTV = 50
                                 if NPSTEP&>= 2000 then INTV = 100
                                 if npstep&>=10000 then intv = 200
                                 if npstep&>=20000 then intv = 500
'
           open D$ + "FILE09V.DAT" for input as #1: ' Read from file09v.dat
'
DATAREAD: line input #1, D1$
           for I =  1 to  8: EGI(I) = val(mid$(D1$,(I- 1)*10+1,10)): next
           line input #1, D2$
           for I =  9 to 16: EGI(I) = val(mid$(D2$,(I- 9)*10+1,10)): next
           line input #1, D3$
           for I = 17 to 24: EGI(I) = val(mid$(D3$,(I-17)*10+1,10)): next
           line input #1, D4$
           for I = 25 to 34: EGI(I) = val(mid$(D4$,(I-25)*9+1,8)): next
           line input #1, D5$
           for I = 35 to 44: EGI(I) = val(mid$(D5$,(I-35)*9+1,8)): next
           '
           for I = 1 to LVA
                  WE = EGI(I)
                  if WE < DMIN(I) then DMIN(I) = WE
                  if WE > DMAX(I) then DMAX(I) = WE
                  WA(I) = WA(I) + WE
                  WD(I) = WD(I) + WE^2
           next I
           MO& = MO& + 1
  if egi(1)<0.1 then
     print mo&
  end if
           if MO& mod ISTEP = 0 then
                NO = NO + 1
                D1(NO) = EGI(IVAL(1))
                D2(NO) = EGI(IVAL(2))
                D3(NO) = EGI(IVAL(3))
                D4(NO) = EGI(IVAL(4))
                D5(NO) = EGI(IVAL(5))
                D6(NO) = EGI(IVAL(6))
                D7(NO) = EGI(IVAL(7))
                D8(NO) = EGI(IVAL(8))
                D9(NO) = EGI(IVAL(9))
                if NO mod INTV = 0 then
                      print using "###### #####  "; MO&; NO;
                      print using "#####.#  ###.####  "; EGI(1); EGI(2);
                      print using "######.#  ##.##### "; EGI(14); EGI(17);
                      for I = 35 to 37: print using " ######.###"; EGI(I); : next I
                      print
                end if
'                if egi(1)<100.0 then sleep
           end if
           if mo& >= NPSTEP& then NO = NO: MO& = MO&: goto READEND
           if eof(1) < 0    then NO = NO: MO& = MO&: goto READEND
           if NO >= NLIM    then NO = NO: MO& = MO&: goto READEND
           goto DATAREAD
           end
'
READEND: close #1
         '------------------------------- Calculate standard deviations, etc.
         print string$(79, "=")
         for I = 1 to LVA: WWW = WA(I) ^ 2 / MO&
              if WD(I) > WWW then
                     WD(I) = sqr((WD(I) - WWW) / MO&)
              else
                     WD(I) = 0
              end if
                     WA(I) = WA(I) / MO&
          next I
         print using "         AVR  #####.#  ###.####  ######.#  ##.##### "; WA(1); WA(2); WA(14); WA(17);
         for I = 35 to 37: print using "######.####"; WA(I); : next I: print
         print using "         STD  #####.#  ###.####  ######.#  ##.##### "; WD(1); WD(2); WD(14); WD(17);
         for I = 35 to 37: print using "######.####"; WD(I); : next I: print
         print string$(79, "-")
         print using "         MAX  #####.#  ###.####  ######.#  ##.##### "; DMAX(1); DMAX(2); DMAX(14); DMAX(17);
         for I = 35 to 37: print using "######.####"; DMAX(I); : next I: print
         print using "         MIN  #####.#  ###.####  ######.#  ##.##### "; DMIN(1); DMIN(2); DMIN(14); DMIN(17);
         for I = 35 to 37: print using "######.####"; DMIN(I); : next I: print
         print string$(79, "=")
         print : color 15
         print using "####### steps computed";nstep&
         print       "I---------------------------------------------------------------------------I"
         print using "I  U =######.####(##.####)kJ/mol   E = UpK  = #######.####(##.####)kJ/mol   I";wa(12);wd(12);wa(14);wd(14)
         print using "I  K =######.####(##.####)kJ/mol   H = EpPV = #######.####(##.####)kJ/mol   I";wa(13);wd(13);wa(16);wd(16)
         print using "I  PV=######.####(##.####)kJ/mol   Molar volume=######.####(###.####)kJ/mol I";wa(15);wd(15);wa(18);wd(18)
         print       "I---------------------------------------------------------------------------I"
         print
         print       "=============================================================================="
         print       "  T/K    P/GPa        A         B         C      Alpha     Beta      Gamma    "
         print       "------------------------------------------------------------------------------"
         print using " ####.# ##.####  ###.##### ###.##### ###.##### ####.#### ####.#### ####.####  ";wa(1);wa(2);wa(19);wa(20);wa(21);wa(22);wa(23);wa(24)
         print using " ####.# ##.####  ###.##### ###.##### ###.##### ####.#### ####.#### ####.####  ";wd(1);wd(2);wd(19);wd(20);wd(21);wd(22);wd(23);wd(24)
         print       "=============================================================================="
         return
'
'============================================================================