!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 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=10            ! I'll probably never need more than 10. 
!
      character :: afile(nlightcurve_num)*180
      character :: aputlabel(nputlabel)*80          ! Maybe some could get long.  Use trim.
      character :: bfile*180
      character :: caption(ntable)*180 
      real (kind=nsngl) :: decaycon
      real (kind=nsngl) :: decayf 
      integer :: i,j,k,l,m,n
      integer :: idash(0:nlightcurve_num)=(/ (i,i=0,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(0:nlightcurve_num)
      integer :: ilightcurve_num
      integer :: ipost
      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) :: qval1=1.729   ! Q_deposited for Ni-56.
      real (kind=nsngl) :: qval2=3.752   ! Q_deposited for Co-56.
      real (kind=nsngl) :: qvalco57=0.1590  ! Q_deposited for Co-57.
      character :: record*180 
      real (kind=nsngl) :: rw7co57=.0225d0/.627d0
      real (kind=nsngl) :: scale 
      real (kind=nsngl) :: t
      real (kind=nsngl) :: t_corr
      real (kind=nsngl) :: t_del
      real (kind=nsngl) :: t_e1
      real (kind=nsngl) :: t_e2
      real (kind=nsngl) :: te_co57
      real (kind=nsngl) :: t_half1=6.077  ! Half-life of Ni-56.
      real (kind=nsngl) :: t_half2=77.27  ! Half-life of Co-56.
      real (kind=nsngl) :: t_half_co57=271.79  ! Half-life of Co-57.
      real (kind=nsngl) :: t_lcmax
      real (kind=nsngl) :: t_max
      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,                                       &
     &              t_corr,                                                &
     &              t_del,                                                 &
     &              t_lcmax,                                               &
     &              t_max,                                                 &
     &              idash,                                                 &
     &              idash_table,                                           &
     &              igrid,                                                 &
     &              xtick,ytick,                                           &
     &              scale,                                                 &
     &              title,xlabel,xmin,xmax,ylabel,ymin,ymax,               &
     &              itable,xytable,caption,xcorr,ycorr,
     &              iputlabel,aputlabel,putlabel
!
       decayf(t)=decaycon*(  (qval1/(56.d0*t_e1))*exp(-t/t_e1)             &
     &                      +(qval2/(56.d0*(t_e2-t_e1)))                   &
     &                        *(exp(-t/t_e2)-exp(-t/t_e1))                 &
     &                      +rw7co57*                                      &
     &                       (qvalco57/(57.d0*te_co57))*exp(-t/te_co57)    &
     &                    )
!
!-----------------------------------------------------------------------
!
       t_e1=t_half1/log(2.d0)
       t_e2=t_half2/log(2.d0)
       te_co57=t_half_co57/log(2.d0)
!
!-----------------------------------------------------------------------
      write(*,*) 't_half1,t_e1,t_half2,t_e2'
      write(*,*)  t_half1,t_e1,t_half2,t_e2
!-----------------------------------------------------------------------
!
       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 creating the decay lightcurve data.'
!-----------------------------------------------------------------------
!
       decaycon=1.d0
       decaycon=1.d0/decayf(t_lcmax) 
!       do t=0.,10.,.1
!         write(*,*) t,decayf(t) 
!       end do
       i=0
       ilightcurve(i)=ceiling(t_max/t_del)+1 
       do402 : do j=1,ilightcurve(i)
         xlightcurve(1,j,i)=t_del*real(j-1,nsngl)
!         write(*,*) j,xlightcurve(1,j,i),xlightcurve(2,j,i)
         t=xlightcurve(1,j,i)           ! It's a mystery why there is
         xlightcurve(2,j,i)=decayf(t)   ! a seg-fault without this.
!         xlightcurve(2,j,i)=decayf(xlightcurve(1,j,i))
!         write(*,*) j,xlightcurve(1,j,i),xlightcurve(2,j,i)
!         if(j .eq. 23) stop
       end do do402
!
       xlightcurve(2,:ilightcurve(i),i)=                                   &
     &   log10( xlightcurve(2,:ilightcurve(i),i) )
!
!-----------------------------------------------------------------------
      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.'
!-----------------------------------------------------------------------
!
      ipost=1
      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
        call sm_ltype(idash(i))
        call sm_conn(xlightcurve(1,:ilightcurve(i),i),                    &
     &               xlightcurve(2,:ilightcurve(i),i),ilightcurve(i) )
      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 
