!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
! Copyright D.J. Jeffery, 2006jan01.
!
!  giant.f is for investigating the giant steps procedure in the
!  case of the plane-parallel grey atmosphere.
!
!   References:
!
!     \bibitem[Metcalf et al.(2004)]{metcalf} Metcalf, M., Reid, J., \&~Cohen, M. 2004,
!          Fortran 95/2003 Explained
!          (Oxford:  Oxford University Press) (MRC)
!          % Pretty good reference book.
!
!     \bibitem[Press et al.(1992)]{press1992} Press, W. H., Teukolsky, S. A.,
!          Vetterling, W. T., \& Flannery, B. P. 1992,
!            Numerical Recipes in Fortran
!            (Cambridge:  Cambridge University Press)
!            % Bill should have given me a complimentary copy when
!            % I co-taught with him in 1993!
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
!
      module late_mod
      use numerical
      use physical
      use astronomical
      use rde_mod2
      implicit none
!
      real (kind=nprecision) :: alpha=1.d0      ! Default.
      real (kind=nprecision) :: beta=1.d0       ! Default.
      real (kind=nprecision) :: coef_late
      integer :: iirad
      integer :: irad
      integer :: istart=0  ! Starting value.
      integer :: itrap
      integer :: jrad(16)=(/10,13,11,14,7,9,3,4,1,2,15,16,8,12,5,6/)
      real (kind=nprecision) :: small
      real (kind=nprecision) :: t_lcmax=19.5d0  ! Default Bmax according 
!                                               !   to conley2006.
      real (kind=nprecision) :: t_late
!
      end module late_mod
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
!
      program lc_model 
      use numerical
      use physical 
      use astronomical 
      use late_mod 
      implicit none
!
      character :: afile*180='lc_model.dat'     ! Default.
      integer :: i,j,k,l,m,n
      integer :: ilightcurve
      real (kind=nprecision) :: t 
      real (kind=nprecision) :: t_del=0.1d0     ! Default.
      real (kind=nprecision) :: t_max=500.d0    ! Default.
      real (kind=nprecision) :: x 
      real (kind=nprecision) :: xlumf 
!
      namelist/lc_modelpar/afile,alpha,beta,irad,itrap,                  &
     & t_del,t_lcmax,t_max
!
      open(unit=1,file='lc_model.par',status='unknown',action='read')
       read(1,lc_modelpar)
       write(*,lc_modelpar)
      close(unit=1) 
!
      open(unit=1,file=afile,status='unknown',action='write')
       write(1,*) 
       write(1,'(a)') 'Free-parameter-free SN Ia light curve model'
       write(1,*) 
       write(1,lc_modelpar) 
       write(1,*) 
       write(1,'(a)') 'END:'
       ilightcurve=ceiling(t_max/t_del)+1
       t=0.d0
       x=xlumf(t)
       do410 : do i=1,ilightcurve   ! xlum(t=0)=0 can't be plotted on a log plot.
         t=t_del*real(i,nprecision)
         x=xlumf(t)
         write(1,'(f10.4,e25.15)') t,x
       end do do410
      close(unit=1)
!
      end program lc_model 
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
!
      function xlumf(t)
      use numerical
      use physical
      use astronomical
      use rde_mod2
      use late_mod
      implicit none
!
      real (kind=nprecision) :: coef_max=1.d0
      real (kind=nprecision) :: coef_rise
      real (kind=nprecision) :: func_late 
      real (kind=nprecision) :: func_max
      real (kind=nprecision) :: func_max1
      real (kind=nprecision) :: func_rise
      real (kind=nprecision) :: heaviside 
      integer :: i,j,k,l,m,n
      real (kind=nprecision) :: smallest 
      real (kind=nprecision) :: t
      real (kind=nprecision) :: xlumf
!
      func_max1(t)=exp(-( (t-t_lcmax)/t_lcmax )**2 )  ! This must come before func_max.
      func_max(t)=coef_max*func_max1(t)
      func_rise(t)=coef_rise*(1.d0-exp(-(t/t_lcmax)**2))
!
      if(istart .eq. 0) then
          istart=1
          small=1.d-10
          smallest=10.d0**(-precision(0.d0))
!
          write(*,*) 'In xlumf before reading radioactive decay data.' 
          open(unit=10,file='rde_out2.dat',status='old',action='read')
           i=1
  110      continue
           read(10,*,end=120) radio(i)
           write(*,*) radio(i)
           i=i+1
           go to 110
  120      continue
           iirad=i-1
          close(unit=10)
!
          write(*,*) 'In xlumf before calculating relevant quantities.' 
          t_late=beta*t_lcmax*sqrt(third*10.d0)  ! Compton down by 3, but free electrons up
!                                          !  up by about 10.   ! Fiducial.
          coef_late=1.d0 
          coef_late=coef_max/(alpha*func_late(t_lcmax))
          write(*,*) 't_late,coef_late'
          write(*,*) t_late,coef_late
          coef_rise=1.d0
          coef_rise=1.d0/func_rise(t_lcmax)
          write(*,*) 't_late,coef_late,coef_rise'
          write(*,*) t_late,coef_late,coef_rise
          write(*,*) 'In xlumf, after istart set-up.'
!          stop
      end if 
!
      if(itrap .eq. 0) then
          xlumf=(func_rise(t)*heaviside(t_lcmax-t)+heaviside(t-t_lcmax))      &
     &          *(func_max(t)+func_late(t)                                    &
     &                    *heaviside(t-t_lcmax)*(1.d0-func_max1(t)) )
        else
          xlumf=func_late(t)
      end if
!
      end function xlumf
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
!
      function func_late(t)
      use numerical
      use physical
      use astronomical
      use rde_mod2
      use late_mod 
      implicit none
!
      integer :: i,j,k,l,m,n
!
      real (kind=nprecision) :: func_late
      real (kind=nprecision) :: func_trap
      real (kind=nprecision) :: func_trap1
      real (kind=nprecision) :: t 
      real (kind=nprecision) :: tmp 
!
      func_trap1(t)=1.d0-exp(-(t_late/max(t,small))**2)
      func_trap(t)=max(func_trap1(t),real(itrap,nprecision))
!
      func_late=0.d0
      do410 : do i=1,irad
        j=jrad(i)
        if(radio(j)%ia_parent .eq. 0) then
            tmp=0.d0
          else
            tmp=t/radio(j)%t_efold_parent
        end if
        func_late=func_late + (                                          &
     &    radio(j)%gfactor  *exp(-t/radio(j)%t_efold)                    &
     &   +radio(j)%hfactor*( exp(-t/radio(j)%t_efold)-exp(-tmp) ) )      &
     &      *(radio(j)%flepton+radio(j)%fphoton*func_trap(t))   
!        write(*,*) 'i,j,func_late,radio(j)'
!        write(*,*) i,j,func_late,radio(j)
!        write(*,*) 'i,j,func_late,t_late,t,func_trap(t)'
!        write(*,*) i,j,func_late,t_late,t,func_trap(t)
!        write(*,*) 'exp(-t/radio(j)%t_efold)-exp(-tmp)'
!        write(*,*) exp(-t/radio(j)%t_efold)-exp(-tmp)
!        write(*,*) 't/radio(j)%t_efold,tmp'
!        write(*,*) t/radio(j)%t_efold,tmp
!        write(*,*) radio(j)
      end do do410
!
      func_late=func_late*coef_late
!
      end function func_late
