PROGRAM pgcont_6loc29q
*
* This is a fortran program which calls pgplot subroutines to plot up 
*  to 100 contours of emission line EWs or ratios.
* 
* It reads in the 3 column fort files of the Quasar emission 
*  line grids of Korista et al. (1997).

*
* Set up here for 6 grids of 29 n(H) by 29 Phi(H).
*  ...modify idim (density index), jdim (flux index) if required:

       parameter (cmax=100,idim=29,jdim=29)

* The emission lines plotted are specified by the names of the fort.xx 
*  files in the OPEN statements below.

* User is queried whether to 'plot ratio' or 'plot EW'. The latter are 
*  shown in the Korista et al. paper. 'ratio' refers to the normalized 
*  equivalent widths described in the README.txt file, or if the user
*  has previously taken the ratio of two fort.xx files (which results
*  in a ratio of two line strengths in the density-flux plane).

* User is then asked for min contour and interval. 
*  If using 'ratio' option, try -2, 0.2 
*  If using 'EW' option, try 0.0, 0.2

* User then asked if time tag desired at bottom of plot (answer t or f).
* User then asked if peaks should be marked on plots (answer t or f). 
*  If 't', user must know in advance the x,y coordinate of the peak in 
*  each plot (from looking through the fort.xx files), and is queried for 
*  each position.

* Finally, enter /VPS as the plot option, then look at pgplot.ps with 
*  ghostview etc after program runs.


*************************************************************************

       REAL Z(idim,jdim)
       REAL CO(cmax),TR(6)
       integer jcont(cmax) 
       logical ok, peaks
*       character*32 file
*       character*64 label
*
*************************************************************************
* starting and increment values for density and flux; MODIFY AS NECESSARY
*
      hdst = 7.0
      hdinc = 0.25
      flst = 17.0
      flinc = 0.25
*
*************************************************************************
*
*
* open data files; modify these and their labels, below, as needed.

* upper left
      OPEN(UNIT=31,FILE='fort.59',status='old')  ! Lya 1216
* upper right
      OPEN(UNIT=32,FILE='fort.45',status='old')  ! O VI 1035
* mid left
      OPEN(UNIT=33,FILE='fort.124',status='old')  ! H-beta 4861 
* mid right
      OPEN(UNIT=34,FILE='fort.82',status='old')  ! C IV 1549
* lower left 
      OPEN(UNIT=35,FILE='fort.84',status='old')  ! He II 1640
* lower right
      OPEN(UNIT=36,FILE='fort.108',status='old')  ! Mg II 2798

*
* set maximum number of contours plotted equal to max

      ncont = cmax

* get some info

      print *
      write(6,30)
30    format('(1) plot ratio; (2) plot EW(1215)')
      read(5,*) iopt

*
*************************************************************************
*
      print *

* Begin contour stuff --
* some contouring information for pgplot x,y arrays:
* array to set up x-y axes.

       hdend = float(idim - 1.)*hdinc + hdst
       flend = float(jdim - 1.)*flinc + flst

       TR(1)=hdst - hdinc
       TR(2)=hdinc
       TR(3)=0.0
       TR(4)=flst - flinc
       TR(5)=0.0
       TR(6)=flinc
*
      write(6,*)'enter log value of smallest contour'
      read(5,*) co(1)
      print *

* set increment size for contours

*      finc = 0.2

      write(6,*)'choose log contour increment: e.g., 0.1, 0.2 or 0.5'
      read(5,*) finc
      print *

      jskip = int(1./finc)

* load up contour values

      icont = ncont - 1
      do i=1,icont
        co(i+1)=co(i)+finc
      end do
*
*
* inquire about name/time tag label on bottom of plot

      write(6,'('' Put on name/time tag (true or false)?'')')
      read(5,*) ok
*
* inquire about IDing peaks of EW distributions
      print *
      write(6,'('' Identify peaks in EW distributions? t or f'')')
      read(5,*) peaks
      print *
*
*
       print *
       write(6,*)'Choose /vps pgplot option for creating ps file'
       print *
*
*
****************************
*      Begin FILE # 1      *
****************************
* skip first 3 lines in file

      READ(31,*)
      READ(31,*)
      READ(31,*)

* read data in file

       I=1
4      DO 6 J=1,jdim

         READ(31,*) temp1, temp2, z(i,j)

* reset very small numbers
         if(z(i,j).lt.1.0e-29) then
           z(i,j) = 1.0e-29
         endif  
 
* check to rescale normalized EW to real EWs, referred
* to continuum at 1215A.

         if(iopt.eq.2) then
           z(i,j) = z(i,j) * 1215.
         endif

         z(i,j) = log10(z(i,j))

6      CONTINUE
       I=I+1
       IF (I.LE.idim) GOTO 4
*
*
* call pgplot programs
*
       CALL PGBEGIN(0,'?',1,1)

       CALL PGSLW(3)
       call pgscf(2)

* define viewport for multiwindows
       call pgvsize(0.75, 3.75, 7.0, 10.0)

       CALL PGSCH(1.0)
       CALL PGSLS(1)
       CALL PGWINDOW(hdst+0.001,hdend-0.001,flst+0.001,flend)
       CALL PGBOX('BCST',1.0,2,'BCNST',1.0,2)
       CALL PGLABEL('','Log \\gF(H)',
     &              ' ')

*
* initialize counter of logrithmic decade contours (e.g. 1,6,11,...
* for contours of 0.2 dex intervals).
* decade pointers will have value 1, others 0 in counter.

       do j=1,ncont
         jcont(j) = 0
       end do

       do j=1,ncont,jskip
         jcont(j) = 1
       end do

* plot the contours, solid or dotted:

       do i=1,ncont
         if(jcont(i).eq.1) then
           CALL PGSLS(1)
         else
           CALL PGSLS(4)
         endif
         CALL PGCONS(z,idim,jdim,1,idim,1,jdim,co(i),1,tr)
       end do

* apply labels

       CALL PGSCH(1.0)
       call pgtext(7.50,23.0,'Ly\\ga\ 1216')
*
* apply reference coordinate to grid, corresponding to 
* log U = -2.00, log n = 10.
       call pgsch(2.0)
       call pgpoint(1,10.00,18.47682,18)
*
* if plotting the "peak EW" symbol at coordinate of peak in
* EW distribution:
       if(peaks) then
         print *
         write(6,*)'enter coordinates (n,f) of peak EW for upper left'
         read(5,*) px, py
         call pgsch(2.0)
         call pgpoint(1,px,py,13)
       endif

*
* note: if the peak EW of an emission line appears on a boundary, then
* give it the actual coordinate +/- delta, so that the coordinate is
* found within the x-y plane. otherwise, pgplot will not plot the symbol.
*
 
****************************
*       End FILE # 1       *
****************************
*
*
****************************
*      Begin FILE # 2      *
****************************
* skip first 3 lines in file

      READ(32,*)
      READ(32,*)
      READ(32,*)

* read data in file

       I=1
7      DO 9 J=1,jdim

         READ(32,*) temp1, temp2, z(i,j)

* reset very small numbers
         if(z(i,j).lt.1.0e-29) then
           z(i,j) = 1.0e-29
         endif  
 
* check to rescale normalized EW to real EWs, referred
* to continuum at 1215A.

         if(iopt.eq.2) then
           z(i,j) = z(i,j) * 1215.
         endif

         z(i,j) = log10(z(i,j))

9      CONTINUE
       I=I+1
       IF (I.LE.idim) GOTO 7
*
*
* call pgplot programs
*
*       CALL PGBEGIN(0,'?',1,1)
       CALL PGSCH(1.0)
       CALL PGSLS(1)
       CALL PGSLW(3)
       call pgscf(2)

* define viewport for multiwindows
       call pgvsize(3.75, 6.75, 7.0, 10.0)

       CALL PGBOX('BCST',1.0,2,'BCST',1.0,2)
       CALL PGWINDOW(hdst+0.001,hdend-0.001,flst+0.001,flend-0.001)
       CALL PGLABEL('','',
     &              ' ')

*
* initialize counter of logrithmic decade contours (e.g. 1,6,11,...
* for contours of 0.2 dex intervals).
* decade pointers will have value 1, others 0 in counter.

       do j=1,ncont
         jcont(j) = 0
       end do

       do j=1,ncont,jskip
         jcont(j) = 1
       end do

* plot the contours 

       do i=1,ncont
         if(jcont(i).eq.1) then
           CALL PGSLS(1)
         else
           CALL PGSLS(4)
         endif
         CALL PGCONS(z,idim,jdim,1,idim,1,jdim,co(i),1,tr)
       end do

* apply labels

       CALL PGSCH(1.0)
       call pgtext(7.50,23.,'O VI 1034')
       call pgsch(2.0)
       call pgpoint(1,10.0,18.47682,18)
       if(peaks) then
         print *
         write(6,*)'enter coordinates (n,f) of peak EW for upper right'
         read(5,*) px, py
         call pgsch(2.0)
         call pgpoint(1,px,py,13)
       endif
 
****************************
*       End FILE # 2       *
****************************
*
****************************
*      Begin FILE # 3      *
****************************
* skip first 3 lines in file

      READ(33,*)
      READ(33,*)
      READ(33,*)

* read data in file

       I=1
10       DO 12 J=1,jdim

          READ(33,*) temp1, temp2, z(i,j)

* reset very small numbers
         if(z(i,j).lt.1.0e-29) then
           z(i,j) = 1.0e-29
         endif  
 
* check to rescale normalized EW to real EWs, referred
* to continuum at 1215A.

         if(iopt.eq.2) then
           z(i,j) = z(i,j) * 1215.
         endif

         z(i,j) = log10(z(i,j))

12      CONTINUE
       I=I+1
       IF (I.LE.idim) GOTO 10
*
*
* call pgplot programs
*
*       CALL PGBEGIN(0,'?',1,1)
       CALL PGSCH(1.0)
       CALL PGSLS(1)
       CALL PGSLW(3)
       call pgscf(2)

* define viewport for multiwindows
       call pgvsize(0.75, 3.75, 4.0, 7.0)

       CALL PGWINDOW(hdst+0.001,hdend-0.001,flst+0.001,flend-0.001)
       CALL PGBOX('BCST',1.0,2,'BCNST',1.0,2)
       CALL PGLABEL('','Log \\gF(H)',
     &              ' ')

*
* initialize counter of logrithmic decade contours (e.g. 1,6,11,...
* for contours of 0.2 dex intervals).
* decade pointers will have value 1, others 0 in counter.

       do j=1,ncont
         jcont(j) = 0
       end do

       do j=1,ncont,jskip
         jcont(j) = 1
       end do

* plot the contours 

       do i=1,ncont
         if(jcont(i).eq.1) then
           CALL PGSLS(1)
         else
           CALL PGSLS(4)
         endif
         CALL PGCONS(z,idim,jdim,1,idim,1,jdim,co(i),1,tr)
       end do

* apply labels

       CALL PGSCH(1.0)
       call pgtext(7.50,23.,'H\\gb 4861')
       call pgsch(2.0)
       call pgpoint(1,10.0,18.47682,18)
       if(peaks) then
         print *
         write(6,*)'enter coordinates (n,f) of peak EW for middle left'
         read(5,*) px, py
         call pgsch(2.0)
         call pgpoint(1,px,py,13)
       endif

 
****************************
*       End FILE # 3       *
****************************
*
****************************
*      Begin FILE # 4      *
****************************
* skip first 3 lines in file

      READ(34,*)
      READ(34,*)
      READ(34,*)

* read data in file

       I=1
13      DO 15 J=1,jdim

         READ(34,*) temp1, temp2, z(i,j)

* reset very small numbers
         if(z(i,j).lt.1.0e-29) then
           z(i,j) = 1.0e-29
         endif  
 
* check to rescale normalized EW to real EWs, referred
* to continuum at 1215A.

         if(iopt.eq.2) then
           z(i,j) = z(i,j) * 1215.
         endif

         z(i,j) = log10(z(i,j))

15      CONTINUE
       I=I+1
       IF (I.LE.idim) GOTO 13
*
*
* call pgplot programs
*
*       CALL PGBEGIN(0,'?',1,1)
       CALL PGSCH(1.0)
       CALL PGSLS(1)
       CALL PGSLW(3)
       call pgscf(2)

* define viewport for multiwindows
       call pgvsize(3.75, 6.75, 4.0, 7.0)

       CALL PGWINDOW(hdst+0.001,hdend-0.001,flst+0.001,flend-0.001)
       CALL PGBOX('BCST',1.0,2,'BCST',1.0,2)
*       CALL PGLABEL('Log n(H)','',
*     &              ' ')

*
* initialize counter of logrithmic decade contours (e.g. 1,6,11,...
* for contours of 0.2 dex intervals).
* decade pointers will have value 1, others 0 in counter.

       do j=1,ncont
         jcont(j) = 0
       end do

       do j=1,ncont,jskip
         jcont(j) = 1
       end do

* plot the contours 

       do i=1,ncont
         if(jcont(i).eq.1) then
           CALL PGSLS(1)
         else
           CALL PGSLS(4)
         endif
         CALL PGCONS(z,idim,jdim,1,idim,1,jdim,co(i),1,tr)
       end do

* apply labels

       CALL PGSCH(1.0)
       call pgtext(7.50,23.,'C IV 1549')
       call pgsch(2.0)
       call pgpoint(1,10.0,18.47682,18)
       if(peaks) then
         print *
         write(6,*)'enter coordinates (n,f) of peak EW for middle right'
         read(5,*) px, py
         call pgsch(2.0)
         call pgpoint(1,px,py,13)
       endif
 
****************************
*       End FILE # 4       *
****************************
*
****************************
*      Begin FILE # 5      *
****************************
* skip first 3 lines in file

      READ(35,*)
      READ(35,*)
      READ(35,*)

* read data in file

       I=1
16       DO 18 J=1,jdim

          READ(35,*) temp1, temp2, z(i,j)

* reset very small numbers
         if(z(i,j).lt.1.0e-29) then
           z(i,j) = 1.0e-29
         endif  
 
* check to rescale normalized EW to real EWs, referred
* to continuum at 1215A.

         if(iopt.eq.2) then
           z(i,j) = z(i,j) * 1215.
         endif

         z(i,j) = log10(z(i,j))

18     CONTINUE
       I=I+1
       IF (I.LE.idim) GOTO 16
*
*
* call pgplot programs
*
*       CALL PGBEGIN(0,'?',1,1)
       CALL PGSCH(1.0)
       CALL PGSLS(1)
       CALL PGSLW(3)
       call pgscf(2)

* define viewport for multiwindows
       call pgvsize(0.75, 3.75, 1.0, 4.0)

       CALL PGWINDOW(hdst,hdend-0.001,flst,flend-0.001)
       CALL PGBOX('BCNST',1.0,2,'BCNST',1.0,2)
       CALL PGLABEL('Log n(H)','Log \\gF(H)',
     &              ' ')

*
* initialize counter of logrithmic decade contours (e.g. 1,6,11,...
* for contours of 0.2 dex intervals).
* decade pointers will have value 1, others 0 in counter.

       do j=1,ncont
         jcont(j) = 0
       end do

       do j=1,ncont,jskip
         jcont(j) = 1
       end do

* plot the contours 

       do i=1,ncont
         if(jcont(i).eq.1) then
           CALL PGSLS(1)
         else
           CALL PGSLS(4)
         endif
         CALL PGCONS(z,idim,jdim,1,idim,1,jdim,co(i),1,tr)
       end do

* apply labels

       CALL PGSCH(1.0)
       call pgtext(7.50,23.,'He II 1640')
       call pgsch(2.0)
       call pgpoint(1,10.0,18.47682,18)
       if(peaks) then
         print *
         write(6,*)'enter coordinates (n,f) of peak EW for lower left'
         read(5,*) px, py
         call pgsch(2.0)
         call pgpoint(1,px,py,13)
       endif
 
****************************
*       End FILE # 5       *
****************************
*
****************************
*      Begin FILE # 6      *
****************************
* skip first 3 lines in file

      READ(36,*)
      READ(36,*)
      READ(36,*)

* read data in file

       I=1
19      DO 21 J=1,jdim

         READ(36,*) temp1, temp2, z(i,j)

* reset very small numbers
         if(z(i,j).lt.1.0e-29) then
           z(i,j) = 1.0e-29
         endif  
 
* check to rescale normalized EW to real EWs, referred
* to continuum at 1215A.

         if(iopt.eq.2) then
           z(i,j) = z(i,j) * 1215.
         endif

         z(i,j) = log10(z(i,j))

21      CONTINUE
       I=I+1
       IF (I.LE.idim) GOTO 19
*
*
* call pgplot programs
*
*       CALL PGBEGIN(0,'?',1,1)
       CALL PGSCH(1.0)
       CALL PGSLS(1)
       CALL PGSLW(3)
       call pgscf(2)

* define viewport for multiwindows
       call pgvsize(3.75, 6.75, 1.0, 4.0)

       CALL PGSCH(1.0)
       CALL PGWINDOW(hdst+0.001,hdend,flst+0.001,flend-0.001)
       CALL PGBOX('BCNST',1.0,2,'BCST',1.0,2)
       CALL PGLABEL('Log n(H)','',
     &              ' ')

*
* initialize counter of logrithmic decade contours (e.g. 1,6,11,...
* for contours of 0.2 dex intervals).
* decade pointers will have value 1, others 0 in counter.

       do j=1,ncont
         jcont(j) = 0
       end do

       do j=1,ncont,jskip
         jcont(j) = 1
       end do

* plot the contours 

       do i=1,ncont
         if(jcont(i).eq.1) then
           CALL PGSLS(1)
         else
           CALL PGSLS(4)
         endif
         CALL PGCONS(z,idim,jdim,1,idim,1,jdim,co(i),1,tr)
       end do

* apply labels

       CALL PGSCH(1.0)
       call pgtext(7.50,23.,'Mg II 2798')
       call pgsch(2.0)
       call pgpoint(1,10.0,18.47682,18)
       if(peaks) then
         print *
         write(6,*)'enter coordinates (n,f) of peak EW for lower right'
         read(5,*) px, py
         call pgsch(2.0)
         call pgpoint(1,px,py,13)
       endif
 
****************************
*       End FILE # 6       *
****************************
*
      if( ok ) then
       call pgiden
      endif

       CALL PGEND
       END