c   This program determines appropriate kiln
c   schedules given the specific gravity of the species.
c
c
c   Steve Verrill, 6/10/96
c
c
c

      implicit double precision (a-h,o-z)

      integer temp(40),depr(40),mc(40)
      integer number(40)

      integer harsh(40)

      integer iorigord(40)

      integer thick

      real*8 gb(40),sdgb(40)
    
      real*8 logsdgb(40),logprob(40),sdgb22(40)

      real*8 loglik(40),liksort(40)

   
      real temp0,depr0,mc0
      integer temp1,depr1,mc1,harsh1
      integer temp2,depr2


c      character*20 results

      character*80 argv,intwrite
      character*1 t,mcsched,alpha(6)

      data alpha/"A","B","C","D","E","F"/

      t = "T"
      

      temp( 1) =      1  
      temp( 2) =      2  
      temp( 3) =      2  
      temp( 4) =      2  
      temp( 5) =      2  
      temp( 6) =      3  
      temp( 7) =      3  
      temp( 8) =      3  
      temp( 9) =      3  
      temp(10) =      4  
      temp(11) =      5  
      temp(12) =      5  
      temp(13) =      5  
      temp(14) =      6  
      temp(15) =      6  
      temp(16) =      6  
      temp(17) =      6  
      temp(18) =      6  
      temp(19) =      6  
      temp(20) =      6  
      temp(21) =      7  
      temp(22) =      7  
      temp(23) =      8  
      temp(24) =      8  
      temp(25) =      8  
      temp(26) =      8  
      temp(27) =      8  
      temp(28) =      8  
      temp(29) =     10  
      temp(30) =     10  
      temp(31) =     10  
      temp(32) =     10  
      temp(33) =     10  
      temp(34) =     10  
      temp(35) =     11  
      temp(36) =     12  
      temp(37) =     12  
      temp(38) =     12  
      temp(39) =     13  
      temp(40) =     14  

      depr( 1) =      1
      depr( 2) =      2
      depr( 3) =      2
      depr( 4) =      2
      depr( 5) =      4
      depr( 6) =      1
      depr( 7) =      1
      depr( 8) =      2
      depr( 9) =      2
      depr(10) =      2
      depr(11) =      2
      depr(12) =      2
      depr(13) =      3
      depr(14) =      2
      depr(15) =      3
      depr(16) =      3
      depr(17) =      3
      depr(18) =      3
      depr(19) =      4
      depr(20) =      4
      depr(21) =      3
      depr(22) =      6
      depr(23) =      2
      depr(24) =      3
      depr(25) =      3
      depr(26) =      4
      depr(27) =      4
      depr(28) =      4
      depr(29) =      4
      depr(30) =      4
      depr(31) =      4
      depr(32) =      4
      depr(33) =      5
      depr(34) =      5
      depr(35) =      4
      depr(36) =      5
      depr(37) =      5
      depr(38) =      7
      depr(39) =      4
      depr(40) =      6

      mc( 1) =      2
      mc( 2) =      2
      mc( 3) =      3
      mc( 4) =      4
      mc( 5) =      4
      mc( 6) =      3
      mc( 7) =      4
      mc( 8) =      3
      mc( 9) =      4
      mc(10) =      2
      mc(11) =      2
      mc(12) =      4
      mc(13) =      3
      mc(14) =      4
      mc(15) =      1
      mc(16) =      2
      mc(17) =      3
      mc(18) =      6
      mc(19) =      1
      mc(20) =      4
      mc(21) =      2
      mc(22) =      2
      mc(23) =      3
      mc(24) =      2
      mc(25) =      4
      mc(26) =      2
      mc(27) =      3
      mc(28) =      4
      mc(29) =      3
      mc(30) =      4
      mc(31) =      5
      mc(32) =      6
      mc(33) =      4
      mc(34) =      6
      mc(35) =      4
      mc(36) =      4
      mc(37) =      5
      mc(38) =      5
      mc(39) =      3
      mc(40) =      3

      number( 1) =      2
      number( 2) =      1
      number( 3) =     25
      number( 4) =      1
      number( 5) =     19
      number( 6) =      2
      number( 7) =      2
      number( 8) =     40
      number( 9) =      9
      number(10) =      1
      number(11) =      1
      number(12) =      1
      number(13) =      1
      number(14) =     63
      number(15) =      1
      number(16) =      1
      number(17) =      2
      number(18) =      1
      number(19) =      1
      number(20) =     18
      number(21) =      2
      number(22) =      1
      number(23) =      1
      number(24) =      3
      number(25) =      1
      number(26) =      2
      number(27) =      3
      number(28) =      4
      number(29) =      1
      number(30) =     20
      number(31) =      3
      number(32) =      2
      number(33) =     16
      number(34) =      2
      number(35) =      1
      number(36) =      1
      number(37) =      1
      number(38) =      2
      number(39) =      9
      number(40) =      1

      gb( 1) =    0.91450  
      gb( 2) =    0.65900  
      gb( 3) =    0.75412  
      gb( 4) =    0.53200  
      gb( 5) =    0.67453  
      gb( 6) =    0.86250  
      gb( 7) =    0.61950  
      gb( 8) =    0.64773  
      gb( 9) =    0.68100  
      gb(10) =    0.58000  
      gb(11) =    0.74000  
      gb(12) =    0.53000  
      gb(13) =    0.41000  
      gb(14) =    0.57454  
      gb(15) =    0.66000  
      gb(16) =    0.40000  
      gb(17) =    0.62500  
      gb(18) =    0.40000  
      gb(19) =    0.51000  
      gb(20) =    0.50856  
      gb(21) =    0.75000  
      gb(22) =    0.36000  
      gb(23) =    0.56000  
      gb(24) =    0.61000  
      gb(25) =    0.60000  
      gb(26) =    0.51000  
      gb(27) =    0.49333  
      gb(28) =    0.47000  
      gb(29) =    0.48000  
      gb(30) =    0.44905  
      gb(31) =    0.39800  
      gb(32) =    0.34500  
      gb(33) =    0.44806  
      gb(34) =    0.34000  
      gb(35) =    0.40000  
      gb(36) =    0.22600  
      gb(37) =    0.46000  
      gb(38) =    0.33500  
      gb(39) =    0.35011  
      gb(40) =    0.33400  

      sdgb( 1) =    0.08273  
      sdgb( 2) =    0.06921  
      sdgb( 3) =    0.10557  
      sdgb( 4) =    0.06921  
      sdgb( 5) =    0.17456  
      sdgb( 6) =    0.06718  
      sdgb( 7) =    0.07000  
      sdgb( 8) =    0.11637  
      sdgb( 9) =    0.07114  
      sdgb(10) =    0.06921  
      sdgb(11) =    0.06921  
      sdgb(12) =    0.06921  
      sdgb(13) =    0.06921  
      sdgb(14) =    0.10735  
      sdgb(15) =    0.06921  
      sdgb(16) =    0.06921  
      sdgb(17) =    0.02121  
      sdgb(18) =    0.06921  
      sdgb(19) =    0.06921  
      sdgb(20) =    0.10029  
      sdgb(21) =    0.04243  
      sdgb(22) =    0.06921  
      sdgb(23) =    0.06921  
      sdgb(24) =    0.14000  
      sdgb(25) =    0.06921  
      sdgb(26) =    0.05657  
      sdgb(27) =    0.05508  
      sdgb(28) =    0.06782  
      sdgb(29) =    0.06921  
      sdgb(30) =    0.11489  
      sdgb(31) =    0.03704  
      sdgb(32) =    0.02121  
      sdgb(33) =    0.08762  
      sdgb(34) =    0.04243  
      sdgb(35) =    0.06921  
      sdgb(36) =    0.06921  
      sdgb(37) =    0.06921  
      sdgb(38) =    0.02121  
      sdgb(39) =    0.13640  
      sdgb(40) =    0.06921  


      harsh( 1) =     1  
      harsh( 2) =     1  
      harsh( 3) =     1  
      harsh( 4) =     1  
      harsh( 5) =     1  
      harsh( 6) =     2  
      harsh( 7) =     2  
      harsh( 8) =     2  
      harsh( 9) =     2  
      harsh(10) =     2  
      harsh(11) =     3  
      harsh(12) =     3  
      harsh(13) =     3  
      harsh(14) =     3  
      harsh(15) =     3  
      harsh(16) =     3  
      harsh(17) =     3  
      harsh(18) =     3  
      harsh(19) =     3  
      harsh(20) =     3  
      harsh(21) =     4  
      harsh(22) =     4  
      harsh(23) =     4  
      harsh(24) =     4  
      harsh(25) =     4  
      harsh(26) =     4  
      harsh(27) =     4  
      harsh(28) =     4  
      harsh(29) =     5  
      harsh(30) =     5  
      harsh(31) =     5  
      harsh(32) =     5  
      harsh(33) =     5  
      harsh(34) =     5  
      harsh(35) =     6
      harsh(36) =     7  
      harsh(37) =     7  
      harsh(38) =     8  
      harsh(39) =     9  
      harsh(40) =     10  

c 4    write(6,5)
c 5    format(/,1x,'What name do you want for your results file?',/)
c      read(5,7) results
c 7    format(a80)

c 9    open(unit=7,file=results,status="new",err=1990)
c      go to 2000

c 1990 write(6,1992) results
c 1992 format(/,1x,'It was not possible to open file',1x,a20,/,1x,
c     x'for writing.  Please use a different file.',/)
c      go to 4


 2000 do 10 i = 1,40

         logprob(i) = log(number(i)/268.0)
         logsdgb(i) = log(sdgb(i))
         sdgb22(i) = 2.0*sdgb(i)*sdgb(i)

 10   continue

c      write(6,20)
c      write(7,20)
c 20   format(/,1x,'What is specific gravity of the species in',
c     x' which you are interested?',/)
c      read(5,*) sg0
c      write(7,*) sg0

      call getarg(1,argv)
      write(intwrite,7) argv
 7    format(a80)
      read(intwrite,*) sg0

      call getarg(2,argv)
      write(intwrite,7) argv
      read(intwrite,*) thick


      do 40 i = 1,40

         loglik(i) = -logsdgb(i) - (sg0 - gb(i))*
     x          (sg0 - gb(i))/sdgb22(i)
     x          + logprob(i)


 40   continue

c
c   sort the loglik's
c

      call rsort(loglik,liksort,iorigord,40)

c
c   obtain the estimated schedules from the regression approach
c
      temp0 = 13.65 - 13.62*sg0
      depr0 =  5.20 -  3.95*sg0
      mc0   =  4.51 -  1.56*sg0

      temp1 = nint(temp0)
      depr1 = nint(depr0)
      mc1 = nint(mc0)

      write(6,44) temp1,mc1,depr1
c      write(7,44) temp1,depr1,mc1

 44   format('

We give results for both the regression', x' and classification approaches below. If the two', x' approaches yield schedules that differ significantly,', x' it might be wise to choose the mildest of the', x' schedules. In this case we would also appreciate', x' hearing from you (there might be a bug in the program)', x'. Please contact us at ', x'', x'.
', x'

Regression Approach

', x'Rounded to the nearest integer, the temp,' x' mc, and depr
predicted', x' by the regression equations are:

', x1x,3(i3,2x)) if (thick .eq. 50) then if (temp1 .lt. 5) then depr1 = depr1 - 1 else if (temp1 .eq. 5) then temp1 = temp1 - 2 depr1 = depr1 - 1 else if (temp1 .lt. 9) then temp1 = temp1 - 3 depr1 = depr1 - 1 else temp1 = temp1 - 2 depr1 = depr1 - 1 endif write(6,41) temp1,mc1,depr1 41 format('

Because the thickness is 50 mm,', x ' the temp, mc, and depr values are changed', x ' to:

', x 1x,3(i3,2x)) endif mcsched = alpha(mc1) write(6,42) 42 format('

') if (temp1 .lt. 10) then if (temp1 .eq. 0) temp1 = 1 if (depr1 .eq. 0) depr1 = 1 write(6,1045) t,temp1,mcsched,depr1 1045 format('') write(6,45) t,temp1,mcsched,depr1 45 format('This corresponds to the schedule', x 2x,'.
(Click on the box', x ' to obtain the full schedule.)

') else write(6,1046) t,temp1,mcsched,depr1 1046 format('') write(6,46) t,temp1,mcsched,depr1 46 format('This corresponds to the schedule', x 2x,'.
(Click on the box', x ' to obtain the full schedule.)

') endif write(6,48) 48 format('
') c 44 format(1x,'Rounded to the nearest integer, the temp,' c x' depr, and mc',/,1x,'predicted', c x' by the regression equations are',//, c x1x,3(i3,2x),//) c23456789012345678901234567890123456789012345678901234567890123456789012 if (thick .eq. 2538) then write(6,49) 49 format('
Classification Approach

', x 'In the table below we list the schedules with', x ' the 5 largest relative likelihoods. These are the', x ' schedules', x ' (among the 40 currently in use) that are most likely to be', x ' appropriate for the specific gravity under consideration.', x '
') else write(6,1049) 1049 format('
Classification Approach

', x 'In the table below we list the schedules with', x ' the 5 largest relative likelihoods. These are', x ' the schedules', x ' (among the 40 currently in use) that are most likely to be', x ' appropriate for the specific gravity under consideration.', x ' They are then modified', x ' to reflect the fact that the thickness', x ' is 50 mm.
') endif if (thick .eq. 2538) then write(6,50) 50 format('
The temperature, moisture content,', x ' wet bulb depression, harshness, ', x 'relative likelihood, and schedule values', x ' are (click on a box to obtain a detailed schedule)', x ':

') write(6,8010) 8010 format("") write(6,8012) 8012 format(" ") write(6,8014) 8014 format("
temp mc depr harsh rel likeli", x " schedule") do 70 i = 40,36,-1 id = iorigord(i) diff = liksort(i) - liksort(40) if (diff .lt. -100.0d0) then ratio = 0.0d0 else ratio = dexp(liksort(i) - liksort(40)) endif temp1 = temp(id) mc1 = mc(id) depr1 = depr(id) harsh1 = harsh(id) mcsched = alpha(mc1) if (temp1 .lt. 10) then write(6,65) temp1,mc1,depr1,harsh1,ratio, x t,temp1,mcsched,depr1, x t,temp1,mcsched,depr1 65 format('
',i2,' ',i1,' ',i1, x ' ',i2,' ',f5.3,' ', x '
', x'', x'
') else write(6,66) temp1,mc1,depr1,harsh1,ratio, x t,temp1,mcsched,depr1, x t,temp1,mcsched,depr1 66 format('
',i2,' ',i1,' ',i1, x ' ',i2,' ',f5.3,' ', x '
', x'', x'
') endif 70 continue write(6,8040) 8040 format("
") else write(6,1050) 1050 format('
The temperature, moisture content,', x ' wet bulb depression, harshness, ', x 'relative likelihood, schedule, and modified schedule', x ' values', x ' are (click on a box to obtain a detailed modified', x ' schedule)', x ':

') write(6,8010) write(6,8012) write(6,9014) 9014 format(" temp mc depr harsh rel likeli", x " schedule modified schedule") do 1070 i = 40,36,-1 id = iorigord(i) diff = liksort(i) - liksort(40) if (diff .lt. -100.0d0) then ratio = 0.0d0 else ratio = dexp(liksort(i) - liksort(40)) endif temp1 = temp(id) mc1 = mc(id) depr1 = depr(id) harsh1 = harsh(id) mcsched = alpha(mc1) if (temp1 .lt. 5) then temp2 = temp1 depr2 = depr1 - 1 else if (temp1 .eq. 5) then temp2 = temp1 - 2 depr2 = depr1 - 1 else if (temp1 .lt. 9) then temp2 = temp1 - 3 depr2 = depr1 - 1 else temp2 = temp1 - 2 depr2 = depr1 - 1 endif if (temp2 .lt. 1) temp2 = 1 if (depr2 .lt. 1) depr2 = 1 if (temp1 .lt. 10) then write(6,1065) temp1,mc1,depr1,harsh1,ratio, x t,temp1,mcsched,depr1, x t,temp2,mcsched,depr2, x t,temp2,mcsched,depr2 1065 format('',i2,' ',i1,' ',i1, x ' ',i2,' ',f5.3,' ', x a1,i1,'-',a1,i1,' ', x '

', x'', x'
') else if (temp2 .lt. 10) then write(6,1066) temp1,mc1,depr1,harsh1,ratio, x t,temp1,mcsched,depr1, x t,temp2,mcsched,depr2, x t,temp2,mcsched,depr2 1066 format('',i2,' ',i1,' ',i1, x ' ',i2,' ',f5.3,' ', x a1,i2,'-',a1,i1,' ', x '
', x'', x'
') else write(6,1067) temp1,mc1,depr1,harsh1,ratio, x t,temp1,mcsched,depr1, x t,temp2,mcsched,depr2, x t,temp2,mcsched,depr2 1067 format('',i2,' ',i1,' ',i1, x ' ',i2,' ',f5.3,' ', x a1,i2,'-',a1,i1,' ', x '
', x'', x'
') endif 1070 continue write(6,8040) endif stop end c This subroutine attempts to perform an n*log(n) sort c rather than an n*n sort. c It sorts reals. c It also keeps track of the original order of the x's. c c c Steve Verrill, 5/2/89 c c c subroutine rsort(x,xsort,iorigord,n) implicit double precision (a-h,o-z) real*8 x(n),xsort(n) integer iorigord(n) do 10 i = 1,n iorigord(i) = i 10 continue num = 2 if (x(1) .lt. x(2)) then xsort(1) = x(1) iorigord(1) = 1 xsort(2) = x(2) iorigord(2) = 2 else xsort(1) = x(2) iorigord(1) = 2 xsort(2) = x(1) iorigord(2) = 1 endif do 1000 j = 3,n c c insert the next value c if (x(j) .ge. xsort(num)) go to 200 if (x(j) .le. xsort(1)) go to 160 il = 1 iu = num 145 numb = iu - il if (numb .eq. 1) go to 170 numbd2 = numb/2 itest = il + numbd2 if (x(j) .le. xsort(itest)) go to 150 il = itest go to 145 150 iu = itest go to 145 160 il = 0 170 num = num + 1 do 180 k = num,il + 2,-1 km1 = k - 1 xsort(k) = xsort(km1) iorigord(k) = iorigord(km1) 180 continue 190 ilp1 = il + 1 xsort(ilp1) = x(j) iorigord(ilp1) = j go to 1000 200 num = num + 1 xsort(num) = x(j) iorigord(num) = j 1000 continue return end