!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
! Copyright D.J. Jeffery, 2006jan01.
!
!  solar.f is for plotting the solar composition.
!
!   References:
!
!     \bibitem[Metcalf et al.(2004)]{metcalf2004} Metcalf, M., Reid, J., \&~Cohen, M. 2004,
!          Fortran 95/2003 Explained
!          (Oxford:  Oxford University Press) (MRC)
!          % Pretty good reference book.
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
      program solar 
      use numerical
      use physical
      use astronomical
      use atomic_mod
      implicit none
!
      real (kind=nsngl) :: abund(natom) 
      integer :: i,j,k,l,m,n
      integer :: idata
      integer :: ipost 
      real (kind=nprecision) :: sum1
      real (kind=nsngl) :: xmax=95.
      real (kind=nsngl) :: xmin=0.
      real (kind=nsngl) :: ylabel(natom) 
      real (kind=nsngl) :: ymax=0.    ! 0 for everything.
      real (kind=nsngl) :: ymin=-11.  ! 11 for everything.
      real (kind=nsngl) :: z_metal 
!
      open(unit=1,file='../../../../aalib/atomic.dat',                   &
     &     status='old',action='read')
        read(1,atomic_data)
      close(unit=1)
!
      abund_solar(:)=10.d0**abund_solar(:)
      sum1=sum(abund_solar(:)*awt(:))                  ! Converting from number fraction 
      abund(:)=real(abund_solar(:)*awt(:)/sum1,nsngl)  !   to mass fraction.
!
      z_metal=sum(abund(3:))/sum(abund(:))
!
!-----------------------------------------------------------------------
      write(*,*)
      write(*,*) 'Solar metalicity ',z_metal
!-----------------------------------------------------------------------
!
      open(unit=1,file='solar.html',status='old',action='read')
      open(unit=2,file='tmp',status='unknown',action='write')
        call readcopy(1,2,1,'<pre>',idata)  
        write(2,*)
        write(2,"(' Element',2x,'At.No.',2x,'Fiducial At.Wt.',           &
     &            5x,'Solar At.Wt.',                                     &
     &            4x,'No. Frac.',3x,'Mass Frac.')")
        do402 : do i=1,natom
           write(2,'(6x,a2,f8.0,f17.0,f17.6,1p,2e13.3)')                 &
     &      aname(i),anumber(i),awt(i),awt2(i),                          &
     &      abund_solar(i),abund(i) 
        end do do402
        write(2,*)
        write(2,"('</pre>')")
        write(2,"('</html>')")
      close(unit=2)
      close(unit=1)
!
      abund(:)=log10(max(abund(:),1.e-30))      ! Remember abund is single precision. 
      ylabel(:)=abund(:)+(ymax-ymin)*.0
!
!-----------------------------------------------------------------------
      write(*,*)
      write(*,*) 'Before composition plot.'
!-----------------------------------------------------------------------
!
      ipost=1
      if(ipost .eq. 1) then
          call sm_device(                                                  &
     &         'postencap '//'solar.ps')        ! Square plot
        else if(ipost .eq. 2) then
          call sm_device(                                                  &
     &         'postportencap '//'solar.ps')    ! Portrait plot
        else
          call sm_device(                                                  &
     &         'postlandencap '//'solar.ps')    ! Landscape plot
      end if
      call sm_graphics
      call sm_defvar('TeX_strings','1')
!
      call sm_ticksize(1.,0.,-1.,10.)
!
      print*,'Before sm_limits.'
      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('\rm Atomic Number')
!      print*,'Before sm_label.'
      call sm_ylabel('\rm Mass Fraction')
!      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,                                                  &
     &  '\rm  Solar Composition (Asplund et al. 2005)')
!
      write(*,*) 'Before connecting the curves.'
      call sm_ltype(0)
      call sm_conn(real(anumber,nsngl),abund,natom)
      do410 : do i=1,natom 
        if(ylabel(i) .lt. ymin) cycle do410
        call sm_relocate(real(anumber(i),nsngl),ylabel(i))
        call sm_putlabel(5,trim(aname(i)))
      end do do410
!
      call sm_alpha
      call sm_hardcopy
!
      write(*,*) 'The abundance plot is completed.'
!
      end program solar