!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