call random_seed ()
do j=1,Nhist
   x=0
   mu=1
10 j=j+1
   call random_number(gamma)
   la=-((1/sigt)*Log(gamma))
   x=x+mu*lambda
   if (x>d) then
   Nout=Nout+1
   go to 10
   else
   Nref=Nref+1
   go to 10

   print *, Nref, Nout

   else
   call random_number(gamma)
   if (gamma<(sigc/sigt)) then
   Nab=Nab+1
   go to 10
   else call random_number(gamma)
   mu=2*gamma-      end if
end do
end