program snlite
      include 'implno.dek'
      include 'const.dek'

c..makes the radioactive portion of a supernova light curve
c..the numbers here are set for supernova 1987a

c..declare
      integer          i,nstep
      double precision xm56,xm57,xm44,xm22,xm60,bignum,
     1                 xln2,tau56,tau57,tau22,tau44,tau60,
     2                 dt,tmin,tmax,tday,t,rhor,xn56,xn57,xn44,xn22,
     3                 xn60,xk56,xk57,xk44,xk22,xk60,xl56,xl57,xl44,
     4                 xl22,xl60,xltot,yn56,yn57,yn60,yn22,yn44,
     5                 xx,w,ye60max,enu,ye60ave,ye44max,ye44ave,
     6                 ye22max,ye22ave,yg60,ye60,yg44,ye44,yg22,ye22,
     7                 yg56,yg57


c..mass of 56ni = co56 = fe56 produced in solar masses
c..mass of 57ni = co57 = fe57 produced
c..mass of ti44 produced in solar masses 
c..mass of na22 produced in solar masses 
c..mass of co60 produced in solar masses 

      xm56 = 0.069d0
      xm57 = 0.0033d0
      xm44 = 1.0d-4
      xm22 = 2.0e-6
      xm60 = 2.0e-5


c..number of atoms at time = 0
      bignum = avo * msol
      yn56   = xm56 * bignum / 56.0d0
      yn57   = xm57 * bignum / 57.0d0
      yn44   = xm44 * bignum / 44.0d0
      yn22   = xm22 * bignum / 22.0d0
      yn60   = xm60 * bignum / 60.0d0


c..half lives of the nuclei in seconds
      xln2  = log(2.0d0)
      tau56 = 77.226d0 * 24.0d0 * 3600.0d0 / xln2
      tau57 = 271.0d0 * 24.0d0 * 3600.0d0 / xln2
      tau22 = 2.60d0 * 3.16d7 / xln2
      tau44 = 60.3d0 * 3.16d7 / xln2
      tau60 = 5.27d0 * 3.16d7 / xln2


c..escape fractions
      xk56 = 0.033d0
      xk57 = 2.4d0 * xk56
      xk44 = 0.04d0
      xk60 = 0.04d0



c..here are the time=0 gamma ray and positron deposition energies

c..gamma rays from 56co 
c..includes positron kinetic energy .125 mev and gamma ray energy 3.57

      yg56     = yn56/tau56 * (3.700d0) * ev2erg * 1.0d6 


c..gamma rays from co57, there is no positron
      yg57     = yn57/tau57 * (0.136d0) * ev2erg * 1.0d6


c..gamma rays from na22
      yg22    = yn22/tau22 * (1.275d0 + 2.0d0*0.511d0) * ev2erg * 1.0d6

c..electron kinetic energy for na22
      ye22max = 0.49d0
      xx      = me * clight * clight 
      w       = 1.0d0 + ye22max/xx
      enu     = 0.5d0 * w * xx * (1.0d0 - 1.0d0/(w*w)) 
     1          * (1.0d0 - 1.0d0/(4.0d0*w) - 1.0d0/(9.0d0*w*w))
      ye22ave = ye22max - enu      
      ye22    = yn22/tau22 * ye22ave * ev2erg * 1.0d6


c..gamma rays from ti44 
      yg44    = yn44/tau44 * 
     1          (1.159d0 + 2.0d0*0.511d0 + 0.068d0 + 0.078d0) *
     2          ev2erg * 1.0d6

c..electron kinetic energy for ti44
      ye44max = 1.47d0
      w       = 1.0d0 + ye44max/xx
      enu     = 0.5d0 * w * xx * (1.0d0 - 1.0d0/(w*w)) 
     1          * (1.0d0 - 1.0d0/(4.0d0*w) - 1.0d0/(9.0d0*w*w))
      ye44ave = ye44max - enu      
      ye44    = yn44/tau44 * ye44ave * ev2erg * 1.0d6


c..gamma rays from co60 
      yg60    = yn60/tau60 * (1.173d0 + 1.332d0) * ev2erg * 1.0d6


c..electron kinetic energy for co60
      ye60max = 0.3178d0
      w       = 1.0d0 + ye60max/xx
      enu     = 0.5d0 * w * xx * (1.0d0 - 1.0d0/(w*w)) 
     1          * (1.0d0 - 1.0d0/(4.0d0*w) - 1.0d0/(9.0d0*w*w))
      ye60ave = ye60max - enu      
      ye60    = yn60/tau60 * ye60ave * ev2erg * 1.0d6



c..write a header
      write(6,10) 'day','co56','co57','lti44','na22','co60','ltot'
10    format(1x,t6,a,t17,a,t28,a,t39,a,t50,a,t61,a,t72,a)




c..dt is the time step in days
      dt     = 10.0d0
      tmin   = 500.0d0
      tmax   = 3500.0d0
      nstep  = int((tmax - tmin)/dt) + 1

      do 100  i=1,nstep
       tday = tmin + float(i-1) * dt
       t    = tday * 24.d0 * 3600.d0


c..column depth from which the radioactivity leaks out 
       rhor = 7.0d4 * (1.0d6/t)**2


c..decay the number of nuclei 
       xn56 = yn56 * exp(-t/tau56)
       xn57 = yn57 * exp(-t/tau57)
       xn44 = yn44 * exp(-t/tau44)
       xn22 = yn22 * exp(-t/tau22)
       xn60 = yn60 * exp(-t/tau60)


c..decay the luminosities
       xl56  = yg56 * (1.0d0 - exp(-xk56*rhor)) * exp(-t/tau56)
       xl57  = yg57 * (1.0d0 - exp(-xk57*rhor)) * exp(-t/tau57)
       xl44  = (yg44 * (1.0d0 - exp(-xk44*rhor)) + ye44) * exp(-t/tau44)
       xl60  = (yg60 * (1.0d0 - exp(-xk60*rhor)) + ye60) * exp(-t/tau60)
       xl22  = (yg22 * (1.0d0 - exp(-xk22*rhor)) + ye22) * exp(-t/tau22)
       xltot = xl56 + xl57 + xl22 + xl44 + xl60


c..log em and write out what we got
       xl56  = log10(xl56)
       xl57  = log10(xl57)
       xl44  = log10(xl44)
       xl60  = log10(xl60)
       xl22  = log10(xl22)
       xltot = log10(xltot)


       write(6,20) tday,xl56,xl57,xl44,xl22,xl60,xltot
20     format(1x,1p7e11.3)
100   continue

      stop 'normal termination'
      end