!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 
