c  THIS PROGRAM WILL CALCULATE THE OPTIMUM TIME TO SPEND COUNTING XRD
c  PATTERNS IN EQUAL BLOCKS OF 2THETA - LONGER COUNT TIMES AT HIGH ANGLE
c  TO ALLOW FOR LOWER INTENSITIES
c  THIS IS THE NEW VERSION OF THE PROGRAM  MODIFIED FEBRUARY 1992      
c  CHANGES INCLUDE THE CALCULATION ON A POINT BY POINT BASIS INSTEAD OF THE 
c  ORIGINAL BLOCK STRUCTURE
c  ROUTINES HAVE BEEN ADDED TO ALLOW FOR DEBYE-SCHERRER GEOMETRY
C  THE PROGRAM IS DIMENSIONED TO ALLOW CALCULATION FOR 17001 DATA POINTS
C  THIS REPRESENTS A SCAN FROM 0 to 171 deg 2THETA IN STEPS OF 0.01deg
C
C  NOTE THAT SOME CALCULATIONS ARE BASED ON RESTRICTIONS IMPOSED BY THE 
C  Philips PW1710 and PW3710 DIFFRACTION CONTROLLERS - IN PARTICULAR
C    the minimum step width for a step scan is 0.005 deg
C    the minimum counting time is 0.05 sec
c  VALUES ARE ROUNDED TO THESE INCREMENTS

      character*1 ans
      character*2 tube
      character*4 symb,symbol(15)
      character*12 system(6)
c
      logical exists
c
      parameter(lowlim = 1)
      parameter(upplim = 160)
      parameter(nsize = 17001)
      parameter(ibk = 170)

      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
      common/sym/symbol

      integer*2 ihr, imin, isec, i100th, iyr, imon, iday
      common/datetime/ihr, imin, isec, i100th, iyr, imon, iday

      common/absorp/ab(100,19)
c
      data system/'Cubic       ','Tetragonal  ','Orthorhombic',
     $            'Hexagonal   ','Monoclinic  ','Triclinic   '/
c
      inquire(file='XRDTIME.OUT',exist=exists)
      if(exists) then
        open(unit=8,file='XRDTIME.OUT',status='old')
        close(unit=8,status='delete')
      end if
c
      open(unit=8,file='XRDTIME.OUT',status='NEW')


1     write(*,10)
10    format(' Enter  START angle (deg)            [D:   10   > '\)
      read(*,20) start
20    format(f10.0)
      if (start.eq. 0.0) start = 10.0
      write(*,25) int(upplim)
25    format(' Enter  STOP  angle (deg)            [D:  ',i3,'   > '\)
      read(*,20) stop
      if (stop. eq. 0.0) stop  = 160.0
      if(stop.lt.start.or.start.lt.lowlim.or.
     $   stop.gt.upplim) then
        write(*,30)
30      format(' Incorrect angle limits - TRY AGAIN '/)
        go to 1
      end if
      
      write(*,40)
40    format(' Enter  STEP WIDTH  (deg)            [D: 0.020  > '\)
      read(*,50) step
50    format(f10.0)
      if(step.eq.0.0) step=0.020
      stpp=(int(0.5+200.0*step))/200.0
      if(abs(step-stpp).gt.0.0005) then
        step=stpp
        write(*,60) step
60      format(' Step width not a multiple of 0.005 - changed to',f6.3)
      end if
      nstp=int(1.5+(stop-start)/step)
      if(nstp.gt.nsize) then
        write(*,70) nstp
70      format(i7,' data points is too many - decrease scan range or ',
     $            'increase step width',/)
        go to 1
      end if
c
      blktmp=1.0
2     write(*,80) int(blktmp)
80    format(' Enter BLOCK WIDTH  (deg)            [D:  ',i3,'   > '\)
      read(*,90) blkwdt
90    format(f10.0)
      if(blkwdt.lt.1.0) blkwdt=1.0
      blkwdt=int(blkwdt+0.5)
      
      write(*,110)
110   format(' Enter  TOTAL COUNTING TIME (sec)    [D: 14400  > '\)
      read(*,120) tt
120   format(f10.0)
      if(tt.eq.0.0) tt=14400
c
      timavg=tt/nstp
      timavg=(int(20.0*timavg+0.5))/20.0

C  TIMMIN IS PART OF THE PATTERN WHICH HAS CONSTANT COUNTING TIME
C  CURRENTLY 1/10th (0.10) OF THE TOTAL TIME IS DISTRIBUTED ACROSS THE 
C  PATTERN AS FIXED COUNT TIME  THE REMAINDER (0.90) IS ALLOWED TO BE VARIABLE
C  THIS ENSURES THAT THE COUNTING TIMES AT THE LOW 2THETA END OF THE PATTERN
C  ARE NOT SO SHORT THAT THERE ARE POOR COUNTING STATISTICS IN THE BACKGROUND

      timmin=0.10*timavg
      timmin=(int(20.0*timmin+0.5))/20.0
      write(*,130) timmin
130   format(' Enter  MINIMUM step count time      [D:',f5.2,'   > '\)
      read(*,140) ttemp
140   format(f10.0)
      if(ttemp.ne.0.0) then
         ttemp=(int(20.0*ttemp+0.5))/20.0
         timmin=ttemp
      end if
c
3     write(*,150)
150   format(' Enter  TUBE TARGET  Fe Co Cu Mo Ag  [D:  Cu    > '\)
      read(*,160) tube
160   format(a)
      wave=0.0
      if(tube.eq.'  ') tube='CU'
      if(tube.eq.'cu'.or.tube.eq.'CU') wave=1.54056              
      if(tube.eq.'fe'.or.tube.eq.'FE') wave=1.93604
      if(tube.eq.'co'.or.tube.eq.'CO') wave=1.78896
      if(tube.eq.'mo'.or.tube.eq.'MO') wave=0.70930
      if(tube.eq.'ag'.or.tube.eq.'AG') wave=0.55941
      if(wave.eq.0.0) then
        write(*,170) tube
170     format(' Unknown tube type ',a,' - TRY AGAIN')
        go to 3
      end if
c
4     write(*,180)
180   format(' Is a monochromator fitted  Y/N      [D:   Y    > '\)
      read(*,190) ans
190   format(a)
      monoflag=0
      if(ans.eq.' ') ans='Y'
      if(ans.eq.'n'.or.ans.eq.'N') monoflag=1
      if(ans.eq.'y'.or.ans.eq.'Y') monoflag=2
      if(monoflag.eq.0) then
        write(*,200) ans
200     format(1x,a,'is an inappropriate response  - TRY AGAIN')
        go to 4
      end if
c      
      write(*,210)
210   format(' Enter  THERMAL VIBRATION (B) value  [D:  0.5   > '\)
      read(*,220) avtemp
220   format(f10.0)
      if(avtemp.eq.0.0) avtemp=0.5
c        
      write(*,230)
230   format(/' Enter  uR  for DEBYE-SCHERRER'/
     $       ' geometry   0.0 for flat plate       [D:  0.0   > '\)
      read(*,240) tmv
240   format(f10.0)
c
5     write(*,250)
250   format(/' Enter code for CRYSTAL SYSTEM'/
     $' 1 = Cubic       2 = Tetragonal   3 = Orthorhombic'/
     $' 4 = Hexagonal   5 = Monoclinic   6 = Triclinic     [D:3 > '\)
      read(*,260) isys
260   format(i5)
      if(isys.eq.0) isys=3
      if(isys.lt.1.or.isys.gt.6) then
        write(*,270) isys
270     format(i5,' is an invalid entry - TRY AGAIN')
        go to 5
      end if
c
      write(*,280)
280   format(/' Enter  CHEMICAL FORMULA ?   Y/N    [D:   N    > '\)
      read(*,290) ans
290   format(a)
      if(ans.eq.' ') ans='N'
      if(ans.ne.'y'.and.ans.ne.'Y') then
        noel=1
        symbol(1)='FE  '
        nat(1)=1.0
      end if
c
      if(ans.eq.'y'.or.ans.eq.'Y') then
        write(*,300)
300     format(/)
        noel=0
6       noel=noel+1
7       write(*,310) noel
310     format(' Enter SYMBOL for element #',i2,
     $         ' (Return to end list) > '\)
        read(*,320) symb
320     format(a)
        if(symb.eq.'    ') go to 111
        iflag=0

C CONVERT ALL ELEMENT SYMBOLS TO UPPER CASE        
        do i=1,4
          itmp=ichar(symb(i:i))
          if(itmp.ge.97.and.itmp.le.122) symb(i:i) = char(itmp - 32)
        end do        
        
        do i=1,256
        if(symb.eq.tbxc(i)) then
          symbol(noel)=tbxc(i)
          iflag=1
        end if
        end do
        if(iflag.eq.0) then
          write(*,330) symb
330       format(' Element symbol ',a4,' not found - TRY AGAIN ')
          go to 7
        end if
        go to 6
111     continue
        noel=noel-1
        write(*,340)
340     format(/)
      do i=1,noel
        write(*,350) symbol(i)
350     format(' Enter number of ',a4,' atoms in formula [D:1> '\)
        read(*,360) nat(i)
360     format(f10.0)
        if(nat(i).eq.0.0) nat(i)=1.0
      end do
c
      write(*,460) (symbol(i),i=1,noel)
      write(*,470) (nat(i),i=1,noel)
      end if
c
      call gettim(ihr, imin, isec, i100th)
      call getdat(iyr, imon, iday)
c
      write(8,370) ihr,imin,iday,imon,iyr,start,step,stop,blkwdt
370   format('XRDTIME - calculation of variable count time'
     $' data collection conditions'/
     $'Calculation run at  ',i2,':',i2.2,'  on  ',2(i2.2,'/'),i4,/
     $'All published use should refer to',/
     $' # Madsen and Hill, (1994), J.Appl.Cryst., 27, 385-392 #',
     $//' Start =',f7.3,'    Step =',f7.3,
     $'    Stop =',f8.3,'     Block width (deg)=',f4.0/)
      write(8,380) tt,timmin,timavg
380   format(' Total count time =',f9.1,' sec',
     $'      Minimum step count time =',f6.2,' sec'//
     $' Step count time for fixed time data collection =',f6.2/)
      write(8,390) tube,wave
390   format(1x,a,' target tube - wavelength =',f8.5/)
C
      if(monoflag.eq.1) write(8,400)
400   format(' Monochromator correction NOT INCLUDED in LP'/)
      if(monoflag.eq.2) write(8,410)
410   format(' Monochromator correction INCLUDED in LP'/)
C
      write(8,420) avtemp
420   format(' Average temperature factor of',f6.2,' assumed.'/)
      write(8,430) tmv
430   format(' mu-R used for absorption correction =',f6.3/)
      write(8,440) system(isys)
440   format(' Crystal system used for multiplicity effect is ',a/)
      write(8,450) 
450   format(' Formula used to calculate scattering factor curves =')
      write(8,460) (symbol(i),i=1,noel)
460   format(' Atom type  ',15(2x,a4))
      write(8,470) (nat(i),i=1,noel)
470   format(' # of atoms ',15(f5.2,1x))
      write(8,480)
480   format(/)
c
      nblk=int(0.9+(step*nstp)/blkwdt)

      write(*,485) blkwdt, nblk      
485   format(//' Block width = ', f4.1,'deg   Number of blocks = ',i3//)

      do i=1,nblk
      ttlo(i)=start+(i-1)*blkwdt
      tthi(i)=start-step+i*blkwdt
      npts(i)=int(1.5+(tthi(i)-ttlo(i))/step)
      end do
      tthi(nblk)=stop
      npts(nblk)=int(1.5+(tthi(nblk)-ttlo(nblk))/step)
c
      npoint=0
      do i=1,nblk
      npoint=npoint+npts(i)
      end do
      write(*,500) npoint
500   format(' Total number of points in pattern =',i6)
      write(*,510)
510   format(//'***  Now performing calculation - PLEASE WAIT  ***'/)
c
      call temfac     
      call lpfact(monoflag)
      call scatfact
      call absfact
      call multi(isys)
      call average(avint)
c
      timfix=timmin*nstp
      timres=tt-timfix
      sumi=0.0
c
      do i=1,nblk
      sumi=sumi+(npts(i)/avint(i))
      end do
c      
      const=timres/sumi
c
      do i=1,nblk
      ttime(i)=npts(i)*timmin  +  npts(i)*const/avint(i)
      end do
c
      sumt=0.0
      do i=1,nblk
      time(i)=ttime(i)/npts(i)
      time(i)=(int(20.0*time(i)+0.5))/20.0
c  THE MINIMUM COUNTING TIME IN THE PW1710/PW3710 CONTROLLERS IS 0.05sec
      if(time(i).lt.0.05) time(i)=0.05
      sumt=sumt+npts(i)*time(i)
      end do    
c
      write(8,600)
600   format(/' Block  #      Angle - Limits    Average     Total',
     $'    Step Time'/
     $        '   #   pts     Low       High   Int./1000    Time',   
     $'     sec/step'/)
      do i=1,nblk
      write(8,610) i,npts(i),ttlo(i),tthi(i),avint(i),ttime(i),time(i)
610   format(i4,i6,2f10.3,3f10.2)
      end do
      write(8,620) sumt
620   format(/28x,'Total time =',f10.2)
      close(unit=8)
c
      write(*,700)
700   format(/' Do you wish to write these conditions to',
     $' PW1710 compatible files Y/N [D:n> '\)
      read(*,710) ans
710   format(a)
      if(ans.eq.' ') ans='N'
      if(ans.eq.'y'.or.ans.eq.'Y') call write
      stop
      end
c
c------------------------------------------------------------------------------
      subroutine temfac
C  THIS SUBROUTINE CALCULATES THE TEMPERATURE AT EACH POINT IN THE PATTERN
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c     
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
c
      common/absorp/ab(100,19)
c
      do i=1,nstp
        ang=start+(i-1)*step
        temp(i)=exp(-2.0*avtemp*((sin(0.017453*ang/2.0))/wave)**2.0)
C
C      write(*,500) i,ang,temp(i)
C500   format(i5,f10.3,f10.3)
C
        end do
c
      return
      end
c
c------------------------------------------------------------------------------
      subroutine lpfact(monoflag)
C  THIS SUBROUTINE CALCULATES THE VALUE OF THE Lp CORRECTION FACTOR AT
C  EACH POINT IN THE PATTERN.
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
c
      common/absorp/ab(100,19)
c
c  MONOCHROMATOR CORRECTION FACTOR (CTHM) = 0.8 IF FLAG=2
C  CTHM = 1.0 IF NO MONOCHROMATOR PRESENT
C
      cthm=1.0
      if(monoflag.eq.2) cthm=0.80
c
      do i=1,nstp
      ang = start+(i-1)*step
      the = 0.01745329 * ang / 2.0
c
c  SEE CULLITY PAGES 139 and 183 FOR LP FORMULA DISCUSSION
C  SEE ALSO PAGE 463 FOR DISCUSSION re MODIFIED LP EXPRESSION
C  APPLICABLE TO STEP INTENSITIES
c
         lp(i) = (1.0 +   cthm   * (cos(2.0*the))**2.0) / 
     $           (2.0 * cos(the) * (sin(the))**2.0) 
C
C      write(*,500) i,ang,lp(i)
C500   format(i5,f10.3,f20.3)
C
      end do
c
      return
      end
c
c------------------------------------------------------------------------------
      subroutine scatfact
C  THIS SUBROUTINE CALCULATES THE VALUE OF THE AVERAGE SCATTERING FACTOR 
C  AT EACH POINT IN THE PATTERN.
c
      character*4 symb,symbol(15)
      dimension ai(15,4),bi(15,4),c(15)
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
      common/sym/symbol
c
      common/absorp/ab(100,19)
c
c  PUT THE SCATTERING FACTORS IN THE Ai, Bi AND C ARRAYS
c
      do i=1,noel
         do j=1,256
           if(symbol(i).eq.tbxc(j)) then
           
             ai(i,1)=tbx(j,1)
             ai(i,2)=tbx(j,3)
             ai(i,3)=tbx(j,5)
             ai(i,4)=tbx(j,7)
             bi(i,1)=tbx(j,2)
             bi(i,2)=tbx(j,4)
             bi(i,3)=tbx(j,6)
             bi(i,4)=tbx(j,8)
             c(i)=tbx(j,9)
           
           end if
         end do
      end do
c
      do i=1,nstp
      ang=start+(i-1)*step
      sin2=(sin(0.017453*ang/2.0))**2.0
      scat(i)=0.0
c
          do jj=1,noel
            tem=0.0
C
            do j=1,4
              tem=tem+ai(jj,j)*exp(-1.0*bi(jj,j)*sin2/(wave*wave))
            end do
C
            tem = tem + c(jj)
            tem = tem * nat(jj)
            scat(i)=scat(i) +  tem  
C
          end do
      scat(i) = (scat(i))**2.0
C
C      write(*,500) i,ang,scat(i)
C500   format(i5,f10.3,f20.3)
c
       end do
c
      return
      end
C
c------------------------------------------------------------------------------
      subroutine absfact
C  THIS SUBROUTINE CALCULATES THE VALUE OF uR WHEN THE DEBYE-SCHERRER 
C  GEOMETRY IS USED.  THE EFFECT OF THE CYLINDRICAL GEOMETRY IS TO INCREASE
C  THE INTENSITY AT HIGH 2THETA RELATIVE TO LOW 2THETA, SO SOME COMPENSATION
C  IS REQUIRED
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
c
      common/absorp/ab(100,19)
c
      do i=1,nstp
      ang=start+(i-1)*step
      the=ang/2.0
c
c  TMV = 0.0  FLAT PLATE GEOMETRY  NO CORRECTION NEEDED
c
      if(tmv.lt.0.000001) then
        absorb(i)=1.0
      end if
c
c  TMV BETWEEN 0.0 AND 1.0 - USE ROUSE et al (1970) FORMULA
C
      if(tmv.gt.0.000001.and.tmv.lt.1.0) then
        sst=(sin(0.017453*the))**2.0
        rouse=-(1.7133-0.0368*sst)*tmv + (0.0927+0.375*sst)*tmv*tmv
        absorb(i)=exp(rouse)
      end if
c
c  TMV GREATER THAN OR EQUAL TO 1.0 - USE WEBER (1967) FORMULA
c
      if(tmv.ge.1.0) then
        itmv1=ifix(tmv*10.0)+1
        itmv2=itmv1+1
        iang1=ifix(the/5.0)+1
        iang2=iang1+1
        frtmv=amod(tmv,1.0)
        iiang1=(iang1-1)*5
        frang=(the-iiang1)/5.0
        trans1=ab(itmv1,iang1)+(ab(itmv1,iang2)-ab(itmv1,iang1))*frang
        trans2=ab(itmv2,iang1)+(ab(itmv2,iang2)-ab(itmv2,iang1))*frang
        trans=trans1+(trans2-trans1)*frtmv
        absorb(i)=1.0/trans
      end if
c
C        write(*,500) ang,absorb(i)
C500     format(f12.3,f12.5)
c
      end do
c
      return
      end
c------------------------------------------------------------------------------      
      subroutine multi(isys)
C  THIS ROUTINE CALCULATES THE APPROXIMATE REFLECTION MULTIPLICITY AS 
C  A FUNCTION OF ANGLE AND CRYSTAL SYSTEM.  THE EQUATIONS ARE DERIVED
C  FROM AVERAGE REFLECTION MULTIPLICITIES GENERATED BY RIETVELD RUNS 
C  ASSUMING A UNIT CELL VOLUME OF 500A^3.  THE EQUATIONS HAVE THE FORM
C  mult = A + B(ln(2theta))  WHERE  A & B  ARE VARIABLES
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
c
      common/absorp/ab(100,19)
c
      
      dimension ai(6),bi(6)
c
c               cubic    tetrag    ortho     hexag    monoc  triclin
c       isys =    1         2        3         4        5        6
      data ai/-5.52353, -1.21295, 0.48112, -1.97655, 1.58523, 2.00000/
      data bi/ 8.01637,  2.93023, 1.35769,  4.21472, 0.42562, 0.00000/
c
      if(isys.lt.1.or.isys.gt.6) then
        write(*,10) isys
10      format(/' Unknown crystal system')
        go to 999
      end if
c
      do i=1,nstp
        ang=start+(i-1)*step
        mult(i)=ai(isys) + bi(isys)*log(ang)
C
C      write(*,500) i,ang,mult(i)
C500   format(i5,f10.3,f10.3)
C
        end do
c
999   return
      end
c
c------------------------------------------------------------------------------
      subroutine average(avint)
c  THIS SUBROUTINE TAKES THE AVERAGE OF THE VALUES OF Lp, SCATTERING FACTOR
C  THERMAL PARAMETER AND ABSORPTION WITHIN THE BLOCK OF 2THETA SPECIFIED.
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv
c
      common/absorp/ab(100,19)
c
      write(8,10)
10    format(' Block  Angle - Limits    Avg.     Avg.', 
     $'     Avg.   Avg.    Avg.    Average'/
     $       '   #    Low      High      LP      Scatt',
     $'    Temp   Trans   Mult.   Int./1000'/)
c
      do i=1,nblk
        tmplp=0.0
        tmptem=0.0
        tmpsca=0.0
        tmpabs=0.0
        tmpmul=0.0
        ibegin=int(1.5+(ttlo(i)-start)/step)
        iend=int(1.5+(tthi(i)-start)/step)
c        write(*,3000) i,ibegin,iend
c3000    format(' Start, End for block',i3,' are ',2i6)
c
        do j=ibegin,iend
          tmplp  = tmplp  + lp(j)
          tmptem = tmptem + temp(j)
          tmpsca = tmpsca + scat(j)
          tmpabs = tmpabs + absorb(j)
          tmpmul = tmpmul + mult(j)
        end do
        tmplp=tmplp/npts(i)
        tmptem=tmptem/npts(i)
        tmpsca=tmpsca/npts(i)
        tmpabs=tmpabs/npts(i)
        tmpmul=tmpmul/npts(i)
c        
        avint(i)=0.001* tmplp * tmptem * tmpsca * tmpabs * tmpmul
        write(8,100) i,ttlo(i),tthi(i),tmplp,tmpsca,
     $  tmptem,tmpabs,tmpmul,avint(i)
100     format(i4,2f9.3,2f9.2,3f8.3,f12.2)  
c
      end do
c
      return
      end
c
c------------------------------------------------------------------------------
      subroutine write
C  THIS SUBROUTINE WRITES OUT THE COMMANDS NEEDED TO CONTROL THE 
C  DIFFRACTOMETER.  THE CODE IS APPROPRIATE FOR USE WITH A PHILIPS PW1710
C  DIFFRACTION CONTROLLER.  THE COMMANDS ARE SENT TO THE PW1710 
C  AND THE RESPONSE IS LOGGED TO AN ASCII FILE
c
      logical exists
c
      parameter(nsize=17001)
      parameter(ibk = 170)
      
      real lp(nsize),temp(nsize),scat(nsize),absorb(nsize),mult(nsize)
      real nat(15)
      dimension time(ibk),ttime(ibk),ttlo(ibk),tthi(ibk),avint(ibk)
      dimension npts(ibk)
c
      character*4 tabnc(128),tbxc(256)
      common/tabchr/tabnc,tbxc
      common/tables/tbx(256,10)
c
      common/factors/lp,temp,scat,absorb,mult
      common/d1/ttlo,tthi,nblk,wave,avtemp,nat,noel,npts,nstp
      common/d2/start,step,stop,tt,timmin,time,timavg,tmv

      integer*2 ihr, imin, isec, i100th, iyr, imon, iday
      common/datetime/ihr, imin, isec, i100th, iyr, imon, iday

      common/absorp/ab(100,19)
c
        character*1 ans
        character*30 fname
C
        fname='COMMAND.LIS'
C
        inquire(file=fname,exist=exists)
        if(exists) then
           open(unit=6,file=fname,status='old')
           close(unit=6,status='delete')
        end if
c
        open(unit=6,file=fname,STATUS='NEW')
        write(6,660) step
660     format('ROF'/'DOF'/'STP ',f6.3)
        write(6,670)
670     format('#*HED'/'#VARIABLE COUNT TIME DATA COLLECTION')
        write(6,680) timmin,tt
680     format('#Min step time = ',f6.2,'   Total Time =',f10.2)
        write(6,690) ihr,imin,iday,imon,iyr
690     format('#XRDTIME calculation run at ',
     $  i2,':',i2.2,'  on  ',2(i2.2,'/'),i4/'# ')

        write(6,700)
700     format('#*TOP')
c
        do i=1,nblk
          write(6,710) ttlo(i),time(i),time(i),tthi(i)-ttlo(i)
710       format('SAN ',f7.3,/,'INT ',f6.2,/,'#*NEW TIME ',f6.2,/
     $     'SSC ',f7.3)
           
           end do
c
        write(6,730) ttlo(1)
730     format('#*END',/,'#  ',/,'SAN ',f6.3,/,'# ')
c
        write(*,800)
800     format(/' Do you want to collect  FIXED TIME  data Y/N [D:Y> '\)
        read(*,810) ans
810     format(a)
        if(ans.eq.' ') ans='Y'
        if(ans.eq.'y'.or.ans.eq.'Y') then
          write(6,820) step,timavg
820       format('ROF'/'DOF'/'STP',f6.3/'INT',f6.2)
          write(6,830) timavg,tt
830       format('#*HED'/'#Fixed Count Time Data Collection Time = ',
     $    f6.2,/,'#Total Time = ',f10.2)
          write(6,690) ihr,imin,iday,imon,iyr
          write(6,840) start,stop-start,start
840       format('SAN',f7.3,/,'#*TOP',/,'SSC',f7.3,/,'#*END',/,
     $    '# ',/,'SAN',f7.3)
        end if
c
        close(unit=6)
c
      return
      end
