!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