!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!
!
!----------------------------------------------------------------------
!
      program lc_model 
      use numerical
      use physical 
      use astronomical 
      implicit none
!
      real (kind=nprecision) :: alpha=1.d0      ! Default.
      real (kind=nprecision) :: beta=1.d0       ! Default.
      integer :: i,j,k,l,m,n
      integer :: ilightcurve
      integer :: irad
      integer :: istart
      integer :: itrap
      real (kind=nprecision) :: t 
      real (kind=nprecision) :: t_del=0.1d0     ! Default.
      real (kind=nprecision) :: t_max=500.d0    ! Default.
      real (kind=nprecision) :: t_lcmax=19.5d0  ! Bmax according to conley2006.
      real (kind=nprecision) :: x 
      real (kind=nprecision) :: xlumf 
!
      namelist/lc_modelpar/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='lc_model.dat',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,t_lcmax,alpha,beta,0,irad,itrap)  ! Do do the initializations in xlumf.
       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,t_lcmax,alpha,beta,1,irad,itrap) ! irad is irrelevant with istart=1.
         write(1,'(f10.4,e25.15)') t,x
       end do do410
      close(unit=1)
!
      end program lc_model 
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
!
      function xlumf(t,t_lcmax,alpha,beta,istart,irad,itrap)
      use numerical
      use physical
      use astronomical
      implicit none
!
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_sc44=0.21854115145076516,          &
     &                         fphoton_sc44=0.7814588485492349
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_ti44=0.06928166558653383,          &
     &                         fphoton_ti44=0.9307183344134662
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_fe55=0.6895724212202615,           &
     &                         fphoton_fe55=0.3104275787797385      
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_co56=0.034567238826442844,         & 
     &                         fphoton_co56=0.9654327611735571
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_ni56=3.989706506630354E-3,         &
     &                         fphoton_ni56=0.9960102934933697
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_co57=0.11068628132212158,          &
     &                         fphoton_co57=0.8893137186778785
      real (kind=nprecision),  parameter ::                               &
     &                         flepton_ni57=0.07575296055923067,          &
     &                         fphoton_ni57=0.9242470394407694
!
      real (kind=nprecision),  parameter :: qval_nu_omit_sc44=2.736967
      real (kind=nprecision),  parameter :: qval_nu_omit_ti44=0.157329
      real (kind=nprecision),  parameter :: qval_nu_omit_fe55=5.5107e-3
      real (kind=nprecision),  parameter :: qval_nu_omit_co56=3.7521
      real (kind=nprecision),  parameter :: qval_nu_omit_co57=0.159008
      real (kind=nprecision),  parameter :: qval_nu_omit_ni56=1.7529
      real (kind=nprecision),  parameter :: qval_nu_omit_ni57=2.100248
      real (kind=nprecision),  parameter :: 
!     &                                     rw7_sc44=.0,                  & Insign.
     &                                      rw7_ti44=.181d-4,             & 
     &                                      rw7_fe55=.614d-2,             & sum of co55 & fe55 
!     &                                     rw7_co56=.0,                  & Insign.
     &                                      rw7_ni56=.627d0,              & 
     &                                      rw7_co57=.903d-3,             & 
     &                                      rw7_ni57=.216d-1 
      real (kind=nprecision),  parameter :: t_half_sc44=3.927/24.0
      real (kind=nprecision),  parameter :: t_half_ti44=63.*yearjday
      real (kind=nprecision),  parameter :: t_half_fe55=2.73*yearjday
      real (kind=nprecision),  parameter :: t_half_co56=77.27d0
      real (kind=nprecision),  parameter :: t_half_co57=271.79
      real (kind=nprecision),  parameter :: t_half_ni56=6.077d0
      real (kind=nprecision),  parameter :: t_half_ni57=35.60/24.0
!
      real (kind=nprecision) :: alpha
      real (kind=nprecision) :: beta 
      real (kind=nprecision) :: coef_max=1.d0
      real (kind=nprecision) :: coef_late 
      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) :: func_trap
      real (kind=nprecision) :: func_trap1
!     real (kind=nprecision) :: gsc44   ! gsc44 unneeded.
      real (kind=nprecision) :: hsc44 
      real (kind=nprecision) :: gti44
      real (kind=nprecision) :: gfe55
!     real (kind=nprecision) :: gco56   ! gco56 unneeded.
      real (kind=nprecision) :: hco56 
      real (kind=nprecision) :: gco57
      real (kind=nprecision) :: hco57
      real (kind=nprecision) :: gni56
      real (kind=nprecision) :: gni57
      real (kind=nprecision) :: heaviside 
      integer :: istart
      integer :: irad
      integer :: itrap
      real (kind=nprecision) :: small
      real (kind=nprecision) :: smallest 
      real (kind=nprecision) :: t
      real (kind=nprecision) :: te_sc44
      real (kind=nprecision) :: te_ti44
      real (kind=nprecision) :: te_fe55
      real (kind=nprecision) :: te_co56
      real (kind=nprecision) :: te_co57
      real (kind=nprecision) :: te_ni56
      real (kind=nprecision) :: te_ni57
      real (kind=nprecision) :: t_late
      real (kind=nprecision) :: t_lcmax
      real (kind=nprecision) :: xlumf
!
      func_trap1(t)=1.d0-exp(-(t_late/max(t,small))**2) 
      func_trap(t)=max(func_trap1(t),real(itrap,nprecision))
      func_late(t)=coef_late*(                                           &
     &     gni56*exp(-t/te_ni56)                                         &
     &      *(flepton_ni56+fphoton_ni56*func_trap(t))                    & ! Unneeded almost
     &    +hco56*(exp(-t/te_co56)-exp(-t/te_ni56))                       &
     &      *(flepton_co56+fphoton_co56*func_trap(t))                    &
!
     &    +gni57*exp(-t/te_ni57)                                         &
     &      *(flepton_ni57+fphoton_ni57*func_trap(t))                    & ! Unneeded
     &    +hco57*(exp(-t/te_co57)-exp(-t/te_ni57))                       &
     &      *(flepton_co57+fphoton_co57*func_trap(t))                    &
     &    +gco57*exp(-t/te_co57)                                         &
     &      *(flepton_co57+fphoton_co57*func_trap(t))                    &
!
     &    +gfe55*exp(-t/te_fe55)                                         &
     &      *(flepton_fe55+fphoton_fe55*func_trap(t))                    &
!
     &    +gti44*exp(-t/te_ti44)                                         &
     &      *(flepton_ti44+fphoton_ti44*func_trap(t))                    &
     &    +hsc44*(exp(-t/te_sc44)-exp(-t/te_ti44))                       &
     &      *(flepton_sc44+fphoton_sc44*func_trap(t))                    &
     &                       )
      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
          small=1.d-10
          smallest=10.d0**(-precision(0.d0))
          te_sc44=t_half_sc44/log(2.d0)
          te_ti44=t_half_ti44/log(2.d0)
          te_fe55=t_half_fe55/log(2.d0)
          te_co56=t_half_co56/log(2.d0)
          te_co57=t_half_co57/log(2.d0)
          te_ni56=t_half_ni56/log(2.d0)
          te_ni57=t_half_ni57/log(2.d0)
!
          if(irad .ge. 4) then
              gti44=( qval_nu_omit_ti44/(44.d0*te_ti44) )*rw7_ti44
              hsc44=( qval_nu_omit_sc44/(44.d0*(te_sc44-te_ti44)) )       &
     &                                                   *rw7_ti44
!             gsc44 insignificant
            else
              gti44=0.d0
              hsc44=0.d0
          end if
!
          if(irad .ge. 3) then
              gfe55=( qval_nu_omit_fe55/(55.d0*te_fe55) )*rw7_fe55
            else
              gfe55=0.d0
          end if
!
          gni56=( qval_nu_omit_ni56/(56.d0*te_ni56) )*rw7_ni56
          hco56=( qval_nu_omit_co56/(56.d0*(te_co56-te_ni56)) )*rw7_ni56
!          gco56 insignificant
!
          if(irad .ge. 2) then
              gni57=( qval_nu_omit_ni57/(57.d0*te_ni57) )*rw7_ni57
              hco57=( qval_nu_omit_co57/(57.d0*(te_co57-te_ni57)) )       &
     &                                                   *rw7_ni57
              gco57=( qval_nu_omit_co57/(57.d0*te_co57) )*rw7_co57
            else
              gni57=0.d0
              hco57=0.d0
              gco57=0.d0
          end if
!
          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))
          coef_rise=1.d0
          coef_rise=1.d0/func_rise(t_lcmax)
!          write(*,*) '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
