!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 ni56 use numerical use physical use astronomical implicit none ! integer, parameter :: ncurve=4 ! Number of curves. integer, parameter :: nputlabel=50 ! I'll probably never need more than 50. integer, parameter :: ntable=20 ! I'll probably never need more than 20. integer, parameter :: nx=10000 ! Number of data points real (kind=nprecision), parameter :: & & t_efold_co56=111.47704580949019 ! More digits than significant. real (kind=nprecision), parameter :: & & t_efold_ni56=8.767257763482231 ! More digits than significant. ! character :: aputlabel(nputlabel)*80 ! Maybe some could get long. Use trim. character :: bfile*80 character :: caption(ntable)*180 real (kind=nprecision) :: fcof real (kind=nprecision) :: fcof2 real (kind=nprecision) :: fco0=0.d0 real (kind=nprecision) :: fco02=1.d0 real (kind=nprecision) :: ffef real (kind=nprecision) :: ffe0=0.d0 real (kind=nprecision) :: fnif real (kind=nprecision) :: fni0=1.d0 integer :: i,j,k,l,m,n integer :: idash(ncurve)=(/0,0,0,1/) integer :: idash_table(ntable) integer :: igrid=-1 ! For grid on plot: -1 none; 0 for major tick marks; ! ! 1 for major and minor tick marks. integer :: ipost=1 integer :: iputlabel=0 integer :: itable=0 integer :: ix real (kind=nsngl) :: putlabel(2,nputlabel) real (kind=nsngl) :: scale real (kind=nsngl) :: tmp character :: title*180 real (kind=nsngl) :: xcorr=1. character :: xlabel*180= & & '\rm Time Relative to Time Zero (days)' real (kind=nsngl) :: x(nx) real (kind=nprecision) :: xx real (kind=nprecision) :: x_del=.1d0 real (kind=nsngl) :: xmin,xmax real (kind=nsngl) :: xtick(2)=0. real (kind=nsngl) :: xytable(ntable) real (kind=nsngl) :: y(nx,ncurve) real (kind=nsngl) :: ycorr=1. character :: ylabel*180='Nuclide Fraction' real (kind=nsngl) :: ymin,ymax real (kind=nsngl) :: ytick(2)=0. ! namelist/param/ & bfile, & & fco0, & & fco02, & & ffe0, & & fni0, & & ipost, & & idash, & & idash_table, & & igrid, & & xtick,ytick, & & scale, & & x_del, & & title,xlabel,xmin,xmax,ylabel,ymin,ymax, & & itable,xytable,caption,xcorr,ycorr, & & iputlabel,aputlabel,putlabel & & ! fnif(xx)=fni0*exp(-xx/t_efold_ni56) fcof(xx)=fco0*exp(-xx/t_efold_co56) & & +fni0*(t_efold_co56/(t_efold_co56-t_efold_ni56)) & & *(exp(-xx/t_efold_co56)-exp(-xx/t_efold_ni56)) ffef(xx)=ffe0+(fni0+fco0-fnif(xx)-fcof(xx)) fcof2(xx)=fco02*exp(-xx/t_efold_co56) ! !----------------------------------------------------------------------- ! open(unit=1,file='./ni56.par',status='old',action='read') read(1,param) ! write(*,param) close(unit=1) ! ix=ceiling( (xmax-xmin)/x_del ) x_del=(xmax-xmin)/real(ix) ix=ix+1 !----------------------------------------------------------------------- write(*,*) write(*,*) 'Before calculating the fractional abundance curves.' !----------------------------------------------------------------------- ! do410 : do i=1,ix xx=x_del*real(i-1,nprecision) ! Prevent error creep. x(i)=real(xx,nsngl) y(i,1)=real(fnif(xx),nsngl) y(i,2)=real(fcof(xx),nsngl) y(i,3)=real(ffef(xx),nsngl) y(i,4)=real(fcof2(xx),nsngl) ! write(*,'(i5,3f15.10)') i,xx,sum(y(i,1:3)),sum(y(i,1:2)) write(*,'(i5,5f15.10)') i,xx,y(i,1),y(i,2),y(i,3),y(i,4) end do do 410 ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'Before calling the SM.' !----------------------------------------------------------------------- ! if(ipost .eq. 1) then call sm_device( & & 'postencap '//bfile) ! Square plot else if(ipost .eq. 2) then call sm_device( & & 'postportencap '//bfile) ! Portrait plot else call sm_device( & & 'postlandencap '//bfile) ! Landscape plot end if call sm_graphics call sm_defvar('TeX_strings','1') ! call sm_ticksize(xtick(1),xtick(2),ytick(1),ytick(2)) ! print*,'Before sm_limits.' write(*,*) 'xmin,xmax,ymin,ymax' write(*,*) xmin,xmax,ymin,ymax call sm_limits(xmin,xmax,ymin,ymax) ! print*,'Before sm_box.' call sm_box(1,2,0,0) ! print*,'Before sm_expand.' call sm_expand(1.) call sm_xlabel(trim(xlabel)) ! print*,'Before sm_label.' call sm_ylabel(trim(ylabel)) ! call sm_expand(1.0) ! print*,'Before sm_grelocate.' call sm_grelocate(2000,0) ! Location in screen coordinates: they run 0--32787 call sm_putlabel(9,trim(title)) ! write(*,*) 'Before connecting the points.' do480 : do i=1,ncurve write(*,*) 'Connecting the nuclide fraction',i, & & ' idash(i) ',idash(i) call sm_ltype(idash(i)) call sm_conn(x,y(:,i),ix) end do do480 ! write(*,*) 'After connecting the points.' ! if(itable .ne. 0) then call sim_scale(xmin,ymin,xmax,ymax) ! xytable(2)=log10(xytable(2)) call sim_table(xytable,itable,idash_table,caption, & & xcorr,ycorr) end if ! do490 : do i=1,iputlabel call sm_relocate(putlabel(1,i),putlabel(2,i)) call sm_putlabel(5, trim(aputlabel(i)) ) end do do490 ! if(igrid .ge. 0) call sm_grid(igrid,0) ! The second parameter is ! ! to prevent a bus error. ! call sm_alpha call sm_hardcopy print*,'Completed is spectrum plot. ' ! end program ni56