!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
! Copyright D.J. Jeffery, 2006jan01.
!
!  lc.f is a code for lotting supernova light curves.
!
!   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 lc 
      use numerical
      use physical 
      implicit none
!
      integer, parameter :: nlightcurve=600000   ! Number of data points. 
      integer, parameter :: nlightcurve_num=10   ! Number of light 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. 
!
      character :: afile(nlightcurve_num)*180
      character :: aputlabel(nputlabel)*80          ! Maybe some could get long.  Use trim.
      character :: bfile*180
      character :: caption(ntable)*180 
      integer :: i,j,k,l,m,n
      integer :: idash(nlightcurve_num)=(/ (i-1,i=1,nlightcurve_num) /) ! Just default.
      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 :: ilightcurve(nlightcurve_num)
      integer :: ilightcurve_num
      integer :: ipost=1
      integer :: iputlabel=0
      integer :: itable=0
      integer :: itmp1,itmp2
      real (kind=nsngl) :: lightcurve(2,nlightcurve,nlightcurve_num)
      real (kind=nsngl) :: putlabel(2,nputlabel)
      real (kind=nsngl) :: scale 
      real (kind=nsngl) :: t
      real (kind=nsngl) :: t_corr
      real (kind=nsngl) :: tmp 
      character :: title*180 
      real (kind=nsngl) :: xcorr=1.
      character :: xlabel*180=                                             &
     & '\rm Time Relative to {\it B} Maximum light (days)'
      real (kind=nsngl) :: xlightcurve(2,nlightcurve,0:nlightcurve_num)
      real (kind=nsngl) :: xmin,xmax
      real (kind=nsngl) :: xtick(2)=0.
      real (kind=nsngl) :: xytable(ntable)
      real (kind=nsngl) :: ycorr=1.
      character :: ylabel*180='\rm Log(Scaled  {\it L})'
      real (kind=nsngl) :: ymin,ymax
      real (kind=nsngl) :: ytick(2)=0.
!
      namelist/lcpar/afile,                                                &
     &              bfile,                                                 &
     &              ilightcurve_num,                                       &
     &              ipost,                                                 &
     &              t_corr,                                                &
     &              idash,                                                 &
     &              idash_table,                                           &
     &              igrid,                                                 &
     &              xtick,ytick,                                           &
     &              scale,                                                 &
     &              title,xlabel,xmin,xmax,ylabel,ymin,ymax,               &
     &              itable,xytable,caption,xcorr,ycorr,                    &
     &              iputlabel,aputlabel,putlabel
!
!-----------------------------------------------------------------------
!
       open(unit=1,file='./lc.par',status='old',action='read')
         read(1,lcpar)
!         write(*,lcpar)
       close(unit=1)
!
       ymin=log10(ymin)
       ymax=log10(ymax)
       putlabel(2,:iputlabel)=log10(putlabel(2,:iputlabel))
!
!-----------------------------------------------------------------------
      write(*,*)
      write(*,*) 'Before reading the lightcurve data.'
!-----------------------------------------------------------------------
!
      do410 : do i=1,ilightcurve_num
        open(unit=2,file=afile(i),status='old',action='read')
         call readcopy(2,0,0,' ',ilightcurve(i))
         write(*,*) 'i,ilightcurve(i)=',i,ilightcurve(i)
         do420 : do j=1,ilightcurve(i)
           read(2,*) xlightcurve(1,j,i),xlightcurve(2,j,i)
!           write(*,*) xlightcurve(1,j,i),xlightcurve(2,j,i)
         end do do420 
         close(unit=2)
         if(trim(afile(i)) .eq. 'lc_dogIa.dat') then
!             write(*,*) 'i,ilightcurve(i)=',i,ilightcurve(i)
             xlightcurve(2,:ilightcurve(i),i)=                                &
     &          -.4d0*xlightcurve(2,:ilightcurve(i),i)
             xlightcurve(1,:ilightcurve(i),i)=                                &
     &        xlightcurve(1,:ilightcurve(i),i)+t_corr
           else
!             write(*,*) 'i,ilightcurve(i)=',i,ilightcurve(i)
             xlightcurve(2,:ilightcurve(i),i)=                                &
     &          log10(xlightcurve(2,:ilightcurve(i),i))
         end if
      end do do 410
!
      do422 : do i=1,ilightcurve_num,-1
        do424 : do j=1,ilightcurve(i)
          write(*,*) xlightcurve(1,j,i),xlightcurve(2,j,i)
        end do do424
      end do do422
!
!-----------------------------------------------------------------------
      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,ilightcurve_num
        write(*,*) 'Connecting the lightcurve ',i,' idash(i) ',idash(i)
        call sm_ltype(idash(i))
        call sm_conn(xlightcurve(1,:ilightcurve(i),i),                    &
     &               xlightcurve(2,:ilightcurve(i),i),ilightcurve(i) )
!
! Note ltype screws up if the the density of points is too high on
! the plot.k
!        call sm_draw(1000.,0.)
      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 lc 
