      
      subroutine thom (u,v,w,mySeed)
c***********************************************************************
c     selects a scattering angle for a thomson collision and then
c     generates new direction cosines.
c***********************************************************************
      implicit double precision (a-h,o-z)
      integer mySeed(*)
      real*8 ranf
      external ranf
      dimension csn (128)
      
      data(csn(i),i=1,67)/
     |     1.       d0,.98952871d0,.9789466 d0,.96825135d0,.95744053d0,
     x .94651168d0,
     |     .93546224d0,.92428964d0,.91299119d0,.90156417d0,.89000577d0,
     x .87831314d0,
     |     .86648335d0,.8545134d0, .84240023d0,.83014072d0,.81773168d0,
     x .80516986d0,
     |     .79245197d0,.77957465d0,.76653448d0,.75332801d0,.73995175d0,
     x .72640218d0,
     |     .71267574d0,.69876886d0,.68467798d0,.67039954d0,.65592999d0,
     x .64126583d0,
     |     .62640361d0,.61133997d0,.59607164d0,.58059548d0,.5649085 d0,
     x .5490079d0,
     |     .53289112d0,.51655582d0,.5d0,       .48322199d0,.46622053d0,
     x .44899477d0,
     |     .43154441d0,.41386968d0,.39597143d0,.3778512 d0,.35951125d0,
     x .34095466d0,
     |     .32218536d0,.30320817d0,.2840289 d0,.26465437d0,.24509241d0,
     x .22535195d0,
     |     .20544297d0,.18537655d0,.1651648 d0,.14482089d0,.12435893d0,
     x .10379394d0,
     |     .083141761d0,.062418936d0,.041642596d0,.020830322d0,.0d0,
     |     -.020830322d0,-.041642596d0/
      
      data(csn(i),i=68,128)/
     |     -.062418936d0,-.083141761d0,-.10379394d0,-.12435893d0,
     x -.14482089d0,
     |     -.1651648d0,-.18537655d0,-.20544297d0,-.22535195d0,
     x -.24509241d0,
     |     -.26465437d0,-.2840289d0,-.30320817d0,-.32218536d0,
     x -.34095466d0,-.35951125d0,
     |     -.3778512d0,-.39597143d0,-.41386968d0,-.43154441d0,
     x -.44899477d0,-.46622053d0,
     |     -.48322199d0,-.5d0,-.51655582d0,-.53289112d0,-.5490079d0,
     x -.5649085d0,
     |     -.58059548d0,-.59607164d0,-.61133997d0,-.62640361d0,
     x -.64126583d0,-.65592999d0,
     |     -.67039954d0,-.68467798d0,-.69876886d0,-.71267574d0,
     x -.72640218d0,-.73995175d0,
     |     -.75332801d0,-.76653448d0,-.77957465d0,-.79245197d0,
     x -.80516986d0,-.81773168d0,
     |     -.83014072d0,-.84240023d0,-.8545134d0,-.86648335d0,
     x -.87831314d0,
     |     -.89000577d0,-.90156417d0,-.91299119d0,
     |     -.92428964d0,-.93546224d0,
     |     -.94651168d0,-.95744053d0,-.96825135d0,-.9789466d0,
     x -.98952871d0/
      
      index  = ranf(mySeed)*128 + 1
      sctang = csn (index)
      t3     = 6.2831853d0 * ranf(mySeed)
      t4     = sqrt(1.0d0 - sctang*sctang)
      t5     = sqrt(1.0d0 - w*w)
      
c--------------------update scattered photon direction
c--------------------(from tid-26607)
      temp1 = t4/t5
      temp2 = sin (t3)
      temp3 = cos (t3)
      temp4 = w*temp3
      t1    = u
      
c--------------------new u, v, w
      
      u = u*sctang + temp1*(u*temp4 - v*temp2)
      v = v*sctang + temp1*(v*temp4 + t1*temp2)
      w = w*sctang - (t4*t5)*temp3
      
      return
      end
