      include 'ase_frdm_composition_module.f'

      program example

      use ase_frdm_composition_module

      implicit none
      
!.....naz(1:kmax): array with number densities of nuclei [fm^-3]
!.....xaz(1:kmax): array with mass fractions of nuclei []
!.....nn: number density of unbound neutrons [fm^-3]
!.....xn: mass fraction of unbound neutrons []
!.....np: number density of unbound protons [fm^-3]
!.....xp: mass fraction of unbound protons []
!.....az(1:kmax,1): array with mass numbers of nuclei
!.....az(1:kmax,2): array with charge numbers of nuclei
!.....kmax: number of considered nuclei
!.....sflag: success flag for composition routine with interpolation

      double precision, dimension(kmax) :: naz, xaz
      double precision :: xn, xp, nn, np
      integer :: sflag
      
      integer :: inuc, i, j, k, ife
      double precision :: t, nb, ye

!.....read in the table with the data for the calculation of the 
!.....composition
      call compdata_readin
      
      write(*,*)'Number of considered nuclei, kmax:'
      write(*,*)kmax
      write(*,*)

!.....EXAMPLES:
!.....examples where the distribution is calculated for given 
!.....EOS grid-point i,j,k, with 
!.....temperature index i, electron fraction index j,
!.....and baryon number density index k:

!.....Example 1: neutron matter at the lowest temperature and density:
      i=1
      j=0
      k=1
!.....find the index for 56Fe
      do inuc=1,kmax
        if (az(inuc,1).eq.56.and.az(inuc,2).eq.26) then
          ife=inuc
        end if
      end do
      write(*,*)'Example 1 for i, j, k:'
      write(*,'(3i5)')i,j,k
      write(*,*)
      call sub_dist_grid(i,j,k,xaz,xn,xp,naz,nn,np)
      write(*,*)'   A    Z  n(A,Z)        X(A,Z)'
      write(*,'(2i5,2e14.6)')1,0,nn,xn
      write(*,'(2i5,2e14.6)')1,1,np,xp
      write(*,'(2i5,2e14.6)')az(ife,1),az(ife,2),naz(ife),xaz(ife)  
      write(*,*)
      write(*,*)
!.....uncomment this part if you want to see the abundances of all nuclei
!       do inuc=1,kmax
!         write(*,'(2i5,2e14.6)')az(inuc,1),az(inuc,2),naz(inuc),xaz(inuc)
!       end do

!.....Example 2: some arbitrary low density grid-point:
      i=24
      j=46
      k=136
      write(*,*)'Example 2 for i, j, k:'
      write(*,'(3i5)')i,j,k
      write(*,*)
      call sub_dist_grid(i,j,k,xaz,xn,xp,naz,nn,np)
      write(*,*)'   A    Z  n(A,Z)        X(A,Z)'
      write(*,'(2i5,2e14.6)')1,0,nn,xn
      write(*,'(2i5,2e14.6)')1,1,np,xp
      write(*,'(2i5,2e14.6)')az(ife,1),az(ife,2),naz(ife),xaz(ife)      
      write(*,*)
      write(*,*)
      
!.....Example 3: the largest t, ye, nb of the table:
      i=81
      j=60
      k=326
      write(*,*)'Example 3 for i, j, k:'
      write(*,'(3i5)')i,j,k
      write(*,*)
      call sub_dist_grid(i,j,k,xaz,xn,xp,naz,nn,np)
      write(*,*)'   A    Z  n(A,Z)        X(A,Z)'
      write(*,'(2i5,2e14.6)')1,0,nn,xn
      write(*,'(2i5,2e14.6)')1,1,np,xp
      write(*,'(2i5,2e14.6)')az(ife,1),az(ife,2),naz(ife),xaz(ife)      
      write(*,*)
      write(*,*)

!.....Examples where the distribution is calculated for given 
!.....temperature t [MeV], electron fraction ye [], 
!.....and baryon number density [fm^-3]

!.....Example 4: neutron matter at the lowest temperature and density
!.....(as example 1):
      t=0.1d0
      ye=0
      nb=1.0d-12
      write(*,*)'Example 4 for t, ye, nb:'
      write(*,'(3e14.6)')t,ye,nb
      write(*,*)
      call sub_dist_interpol(t,ye,nb,xaz,xn,xp,naz,nn,np,sflag)
      write(*,*)'sflag: ',sflag
      write(*,*)'   A    Z  n(A,Z)        X(A,Z)'
      write(*,'(2i5,2e14.6)')1,0,nn,xn
      write(*,'(2i5,2e14.6)')1,1,np,xp
      write(*,'(2i5,2e14.6)')az(ife,1),az(ife,2),naz(ife),xaz(ife)  
      write(*,*)
      write(*,*)
      
!.....Example 5: some arbitrary point inside the table
      t=4.0d0*datan(1.0d0)
      ye=dexp(1.0d0)*0.1d0
      nb=1.0d-3*datan(1.0d0)
      write(*,*)'Example 5 for t, ye, nb:'
      write(*,'(3e14.6)')t,ye,nb
      write(*,*)
      call sub_dist_interpol(t,ye,nb,xaz,xn,xp,naz,nn,np,sflag)
      write(*,*)'sflag: ',sflag
      write(*,*)'   A    Z  n(A,Z)        X(A,Z)'
      write(*,'(2i5,2e14.6)')1,0,nn,xn
      write(*,'(2i5,2e14.6)')1,1,np,xp
      write(*,'(2i5,2e14.6)')az(ife,1),az(ife,2),naz(ife),xaz(ife)  
      write(*,*)
      write(*,*)
      
!.....Example 6: some arbitrary point outside the table
      t=1.0d3
      ye=0.3d0
      nb=1.0d-5
      write(*,*)'Example 6 for t, ye, nb:'
      write(*,'(3e14.6)')t,ye,nb
      write(*,*)      
      call sub_dist_interpol(t,ye,nb,xaz,xn,xp,naz,nn,np,sflag)
      write(*,*)'sflag: ',sflag
      write(*,*)'   A    Z  n(A,Z)        X(A,Z)'
      write(*,'(2i5,2e14.6)')1,0,nn,xn
      write(*,'(2i5,2e14.6)')1,1,np,xp
      write(*,'(2i5,2e14.6)')az(ife,1),az(ife,2),naz(ife),xaz(ife)  
     
      end program example