C @(#)gra.for	10.2 (ESO-IPG) 2/9/96 17:22:13

C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C       Internet e-mail: midas@eso.org
C       Postal address: European Southern Observatory
C                       Data Management Division 
C                       Karl-Schwarzschild-Strasse 2
C                       D 85748 Garching bei Muenchen 
C                       GERMANY
C===========================================================================
C

      
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE DefGrW(ireg)
C     
C     Permette di plottare diverse regioni dello spettro
C     
C     Lavora con le variabili di default XSPETTRO
C     
C     Ritorna il numero di regioni selezionate.
C     
C     Version 1 - 26/X/93  at ESO
C     Version 2 - 07/XII/95  at OAR:
C                + Velocity space
C                + returns meaningful ireg
C                + global variables updated only when input is complete
C     
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      
      implicit none

      integer ireg      
      include 'MID_REL_INCL:fit_var.inc'

      character space
      character*72 stt_ap
      character*65 Commnt
      character*60 lbl(NMXLIN)
      character*40 domand
      integer istt,iredo,i,icurer,ij,nj,ik,j
      logical iexist,i_red
      real WinCen(NMXLIN),ydum,vel,c,lam0
      double precision ranvel,WinSiz(NMXLIN),z

      parameter (c=2.997e5)

      INTEGER MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      
      COMMON/VMR/MADRID
      
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'

     
C&&&&&&&
C&&   Menu con loop se l'input non e' corretto (la flag e' "iscelta")
C&&&&&&&
      
 100  call sttdis(' ',0,istt)
      call sttdis(' ',0,istt)
      call sttdis('****** Define Graphic Window  ******',0,istt)
      write(Commnt,'(A30,F10.2,A5,F10.2)')
     1     '---> Spectrum ranging from: ',XLAMBD(1),' to :',
     2     XLAMBD(NPUNTI)
      call sttdis(Commnt,0,istt)
      call sttdis(' ',0,istt)

      ireg = -1                 !default return

      if (I_Z) space='V'
      if (.not.I_Z) space='L'
      i_red = I_Z
 99   call AskC('(L)ambda.or.(V)elocity space?',space,iredo)
      if (iredo.eq.-1) return
      if (space.eq.'L'.or.space.eq.'l') I_red=.false. 
      if (space.eq.'V'.or.space.eq.'v') I_red=.true. 

                                ! NUMERO  DELLE REGIONI
      call AskI('### Number of regions:',NREGIO,iredo)
      if (iredo.eq.-1) goto 100 !torna al menu


                                ! LIMITI DELLE REGIONI
      if (I_red) then

         call sttdis('Plot will be performed in the VELOCITY space',
     1        0,icurer)
         z = REDSH
         call AskD('Input central redshift or wavelength',z,iredo)
         if (iredo.eq.-1) goto 100 !torna al menu
         if (I_GRAP.and.iredo.eq.-98) then ! legge cursore
           call WrnMsg('I will use the WAVELENGTH from the cursor:')
           call GetXCr(WinCen(1),ydum,z,vel,icurer)  
           z = WinCen(1)
         end if

      else

         call sttdis('Plot will be performed in the WAVELENGTH space',
     1        0,icurer)
         do i  = 1,NREGIO
            write(Commnt,'(a18,i2,a10)')' Input CENTER of ',i
     1           ,' ^ region '
            if (REGMIN(i).ge.0) WinCen(i) = (REGMIN(i) + REGMAX(i))/2.
            call AskR(Commnt,WinCen(i),iredo)
            if (iredo.eq.-1) goto 100 !torna al menu
            if (I_GRAP.and.iredo.eq.-98) then ! legge cursore
               call GetXCr(WinCen(i),ydum,z,vel,icurer) 
            end if
         end do

      endif

 200  continue                  ! SIZES
      if (i_red) then
         ranvel = VELRAN
         call AskD('Input VELOCITY range',ranvel,iredo)
         if (iredo.eq.-1) goto 100
      else
         do i  = 1,NREGIO
            write(Commnt,'(a18,i2,a10)')' Input SIZE of ',i
     1           ,' ^ region'
            winsiz(i) = WINSTE(i)
            if (winsiz(i).le.0) winsiz(i)=winsiz(i-1)
            call AskD(Commnt,winsiz(i),iredo)
            if (iredo.eq.-1) goto 100
         end do
      endif


                                ! LABELS
      do i  = 1,NREGIO
         write(Commnt,'(a,i2,a)')' Input LABEL for ',i
     1        ,' ^ region (NULL to clear)'
         iredo=0
         lbl(i) = LLABEL(i)
         domand=lbl(i)
 1002    call AskC(Commnt,domand,iredo)
         lbl(i) = domand
         if (lbl(i).eq.'null'.or.lbl(i).eq.'NULL') lbl(i)=' '
         if (iredo.eq.-1) goto 100 !torna al menu
                                ! help
         if (Lbl(i).eq.'help'.or.Lbl(i).eq.'HELP') then
            call sttdis('*** List of elements: **',0,istt)
            nj=nint(AT_N/5.)+1
            do ij=1,nj
               write(stt_ap,'(1x,5A14)')(AT_NAM(ij+ik*nj),ik=0,4)
               call sttdis(stt_ap,0,istt)
            end do
            goto 1002
         end if

         if (i_red) then          ! Checks whether it exists.
            iexist = .false.    
            do ij=1,AT_N
               if (Lbl(i).eq.AT_NAM(ij))  iexist = .true.
            end do  
            if (.not.iexist) goto 1002
         endif
      end do        



                                ! computes limits..
      do i = 1 , NREGIO
         LLABEL(i) = lbl(i)
      end do

      if (i_red) then

         REDSH = z
         VELRAN = ranvel
         
         do i = 1 , NREGIO       ! trova le lam0
           lam0 = -1.
           do j = 1 ,AT_N
              if (LLABEL(i).eq.AT_NAM(j)) lam0=AT_LAM(j)
           end do
           if (i.eq.1.and.REDSH.gt.100) REDSH = REDSH / lam0 -1.
           REGMIN(i)=(1.d0-VELRAN/c)/(1.d0+VELRAN/c)*lam0*(REDSH+1.D0)
           REGMAX(i)=(1.d0+VELRAN/c)/(1.d0-VELRAN/c)*lam0*(REDSH+1.D0)
           WINSTE(i)= REGMAX(i) - REGMIN(i)           
        end do

      else
        
         REDSH=0.
         VELRAN=0.
         do i  = 1,NREGIO
            WINSTE(i) = winsiz(i)
            REGMIN(i)= WinCen(i) - (WINSTE(i) / 2.)
            REGMAX(i)= WinCen(i) + (WINSTE(i) / 2.)
         end do        
      
      endif
      
      I_GRAP = .TRUE.
      I_Z = i_red
      ireg=0

      return
      end



C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        SUBROUTINE   GraMai(ierr)
C
C     Main plot routine.
C     Calls appropriate graph routine
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        IMPLICIT NONE
        integer ierr

        integer ist
        include 'MID_REL_INCL:fit_var.inc'

C%%%%%%
C       Is graph active? Are regions defined?
C%%%%%%
        if (.not.I_GRAP) then
           ierr=-1
           return
        end if

        if (NREGIO .eq. 0) then
           NREGIO = 1
           REGMIN(1) = XLAMBD(1)
           REGMAX(1) = XLAMBD(NPUNTI)
           I_Z = .false.
        end if

        if (I_Z) then
           call GraRed(ist)     ! z plot
           ierr = ist
        else
           call GraLam(ist)     ! lambda plot 
           ierr = ist
        endif
        
        return
        end


C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        SUBROUTINE   GraRed(ierr)
C
C     Plot routine in the redshift space.
C     Calls appropriate graph routine
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        IMPLICIT NONE
        integer ierr

        include 'MID_REL_INCL:fit_var.inc'


        real rx,ry,rgroff
        real xdef(4),ydef(4)
        real   dummy(2)
        real grf_w(NMXWRK),grf_f(NMXWRK),grf_r(NMXWRK),lin_ww(NMXLIN)
        real grf_ww(NMXWRK),grf_fw(NMXWRK),grf_rw(NMXWRK)
        real lim_x(4), lim_y(4),rgra
        real lam0(20),z,dv
        integer colore(4),ipt,i,imode,k,nptgra,j,m

C        include 'MID_REL_INCL:fit_mid.inc'

        double precision c
        parameter (c=2.997D5)
        INTEGER MADRID(1)
        INCLUDE 'MID_INCLUDE:ST_DEF.INC'

        COMMON/VMR/MADRID

        INCLUDE 'MID_INCLUDE:ST_DAT.INC'

        if (NGRALI .gt. 0) then
           do i = 1,NGRALI
              YLIN(i) = GR_TIC
           end do
        end if

        do i = 1 , NREGIO       ! trova le lam0
           lam0(i) = -1.
           do j = 1 ,AT_N
              if (LLABEL(i).eq.AT_NAM(j)) lam0(i)=AT_LAM(j)
           end do
           if (lam0(i).lt.0) then
              call ErrMsg('Transition or element unknown')
              ierr = -1
              return
           end if
        end do

        if (REDSH.gt. 100)  REDSH = REDSH/ lam0(1) - 1.
        z = REDSH + 1

                                ! ricalcola REGMIN e REGMAX
        do i = 1, NREGIO
           REGMIN(i) = (1.d0-VELRAN/2./c) / (1.d0+VELRAN/2./c) * 
     1          lam0(i) * z
           REGMAX(i) = (1.d0+VELRAN/2./c) / (1.d0-VELRAN/2./c) * 
     +          lam0(i) * z
           WINSTE(i) = REGMAX(i) - REGMIN(i)
        end do

c%%%%%%% 
C        inserisco limiti e set up
C%%%%%%%

        GR_YMN = - GR_BOT
        GR_DEL = 1 + GR_BOT + GR_TOP
        GR_YMX = NREGIO * GR_DEL - GR_BOT


        call PTKWRC('BINMODE','ON')
        CALL PTKWRR('SCALES',0,dummy)
        CALL PTKWRR('OFFS',0,dummy)
        if (NREGIO .ne. 1) then
           call PTKWRC('YFORMAT','NONE') 
        end if

        xdef(1)=-VELRAN
        xdef(2)=VELRAN
        xdef(3)=2*VELRAN/5.
        xdef(4)=xdef(3)/5.
        
        ydef(1)=GR_YMN 
        ydef(2)=GR_YMX 
        ydef(3) = 1.            !(GR_YMX - GR_YMN) * 100.
        ydef(4)=  ydef(3) / 5


        do i = 1,NREGIO

C%%%%%% 
C     creo i vettori grafici
C%%%%%%
          nptgra=1
          do j = 1,NPUNTI
             dv = 2.*c*(XLAMBD(j) - lam0(i) * z )/(XLAMBD(j)+lam0(i)*z)
             if(abs(dv).le.VELRAN)then 
                GraWav(nptgra)=dv
                GraSpe(nptgra)=SPECTR(j)
                GraSig(nptgra)=SIGNOI(j)
                nptgra=nptgra+1
             end if
          end do   
          nptgra=nptgra-1



C%%%%%%%
C       OPEN WINDOW
C%%%%%%%
          call PTKWRR('XWNDL',4,xdef)
          call PTKWRR('YWNDL',4,ydef)
          imode = -1
          if (i.eq.1)then
            call PTOPEN(' ','fitlyman.plt',0,imode) ! plot
            call AGVERS()
          else
             call PTOPEN(' ','fitlyman.plt',1,imode)    ! overplot     
          end if
          if (i.eq.1) call PTFRAM(xdef,ydef,'Velocity (km/s)',
     1         'Normalized flux') 

          GR_OFF = i * GR_DEL - GR_BOT - GR_TOP - 1.
          rgroff = GR_OFF

          colore(1)=GCOLOR(2)
          call PTKWRI('COLOUR',1,colore)
          call PTDATA(5,1,0,GraWav,CONTIN,rgroff,nptgra)
          call PTDATA(5,3,0,GraWav,SKYLEV,rgroff,nptgra)   

          colore(1)=GCOLOR(1)
          call PTKWRI('COLOUR',1,colore)
          call PTDATA(5,1,1,GraWav,GraSpe,rgroff,nptgra)
          ry = rgroff +1 + GR_LAB
          rx = - 4 * VELRAN / 5
          colore(1)=GCOLOR(6)
          call PTKWRI('COLOUR',1,colore)
          call PTTEXT(LLABEL(i),rx,ry,0.,.7,1)
          colore(1)=GCOLOR(1)
          call PTKWRI('COLOUR',1,colore)

        end do


C!!!!!!!
C       Draws the fit
C!!!!!!!

C     reads fdummy.grf
        open(25,file='fdummy.grf',status='old',iostat=ipt,err=29)
        ipt=1
 122    read(25,*,end=123,err=123)grf_w(ipt),grf_f(ipt),grf_r(ipt)
        ipt=ipt+1
        goto 122
 123    continue
        ipt=ipt -1
        if (ipt.lt.1) goto29
        close(25)

        lim_x(1)=-VELRAN
        lim_x(2)=VELRAN
        lim_y(1)= GR_YMN
        lim_y(2)= GR_YMX
        imode = -1
        
        do i = 1,NREGIO
            do k=1,NMXWRK
               grf_ww(k)=0.
               grf_rw(k)=0.
               grf_fw(k)=0.
            end do


            j = 1
            do k = 1,ipt
             dv = 2.*c*(grf_w(k) - lam0(i) * z )/(grf_w(k)+lam0(i)*z)
               if(abs(dv).le.VELRAN)then
                  grf_ww(j)=dv
                  grf_fw(j)=grf_f(k)
                  grf_rw(j)=grf_r(k)
                  j = j+1
               end if
            end do
            j = j-1

            call PTKWRR('XWNDL',2,lim_x)        
            CALL AGWDEF(lim_x(1),lim_x(2),lim_y(1),lim_y(2))

                                ! position of fitted lines
            if (NGRALI .gt. 0) then
               m = 1
               do k = 1,NGRALI
                  dv=2.*c*(Lamcen(k) - lam0(i)*z)/(lamcen(k)+lam0(i)*z)
                  if(abs(dv).le.VELRAN)then
                     lin_ww(m)=dv
                     m = m+1
                  end if
               end do
               m = m-1

               colore(1)=GCOLOR(7)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               dummy(1)=2.
               CALL PTKWRR('SSIZE',1,dummy)
               call PTDATA(13,0,0,Lin_ww,YLIN,rGRA,NGRALI)
               dummy(1)=1.
               CALL PTKWRR('SSIZE',1,dummy)
            end if

                                ! residuals
            if (I_RESD.and.j.gt.0) then
               colore(1)=GCOLOR(4)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               call PTDATA(5,1,1,grf_ww,grf_rw,rGRA,j)
            end if

                                ! fitted profile
            if (j.gt.0) then
               colore(1)=GCOLOR(3)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               call PTDATA(5,1,0,grf_ww,grf_fw,rGRA,j)
            end if
                                ! variance
            if (I_VAR) then
               colore(1)=GCOLOR(5)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               call PTDATA(5,1,0,GraWav,GraSig,rGRA,nptgra)
            end if

          end do

          colore(1)=GCOLOR(1)
          call PTKWRI('COLOUR',1,colore)

          call PTKWRR('XWNDL',2,lim_x)        
          CALL AGWDEF(lim_x(1),lim_x(2),lim_y(1),lim_y(2))


 29     continue


c        call PTKWRC('BINMODE','OFF')
        call PTKWRC('YFORMAT','AUTO') 


        return
        end

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        SUBROUTINE   GraLam(ierr)
C
C     Plot routine in the lambda space.
C     Calls appropriate graph routine
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        IMPLICIT NONE
        integer ierr

        include 'MID_REL_INCL:fit_var.inc'
        real rx,ry,rgroff
        real xdef(4),ydef(4)
c        real xcont(10),ycont(10)
        character*15 cwave
        real   dummy(2)
        real grf_w(NMXWRK),grf_f(NMXWRK),grf_r(NMXWRK)
        real grf_ww(NMXWRK),grf_fw(NMXWRK),grf_rw(NMXWRK)
        real lim_x(4), lim_y(4),rgra
        integer colore(4),ipt,i,imode,k,nptgra,j

C        include 'MID_REL_INCL:fit_mid.inc'
        INTEGER MADRID(1)
        INCLUDE 'MID_INCLUDE:ST_DEF.INC'

        COMMON/VMR/MADRID

        INCLUDE 'MID_INCLUDE:ST_DAT.INC'



C%%%%%% 
C     creo i vettori grafici
C%%%%%%
        nptgra=1
        do i = 1,NPUNTI
          do k = 1 , NREGIO
             if(XLAMBD(i).gt.REGMIN(k).and.XLAMBD(i).lt.REGMAX(k))then 
                GraWav(nptgra)=XLAMBD(i)
                GraSpe(nptgra)=SPECTR(i)
                GraSig(nptgra)=SIGNOI(i)
                nptgra=nptgra+1
                goto 10
             end if
 10          continue
          end do   
        end do
        nptgra=nptgra-1

        if (NGRALI .gt. 0) then
           do i = 1,NGRALI
              YLIN(i) = GR_TIC
           end do
        end if

c%%%%%%% 
C        inserisco limiti e set up
C%%%%%%%

        GR_YMN = - GR_BOT
        GR_DEL = 1 + GR_BOT + GR_TOP
        GR_YMX = NREGIO * GR_DEL - GR_BOT


        call PTKWRC('BINMODE','ON')
        CALL PTKWRR('SCALES',0,dummy)
        CALL PTKWRR('OFFS',0,dummy)
        if (NREGIO .ne. 1) then
           call PTKWRC('YFORMAT','NONE') 
        end if


        do i = 1,NREGIO
C%%%%%%%
C       OPEN WINDOW
C%%%%%%%
          xdef(1)=REGMIN(i)
          xdef(2)=REGMAX(i)
          xdef(3)=-(REGMIN(i)-REGMAX(i))/5.
          xdef(4)=xdef(3)/5.

          ydef(1)=GR_YMN 
          ydef(2)=GR_YMX 

          ydef(3) = 1. !(GR_YMX - GR_YMN) * 100.
          ydef(4)=  ydef(3) / 5

          call PTKWRR('XWNDL',4,xdef)
          call PTKWRR('YWNDL',4,ydef)
          imode = -1
          if (i.eq.1)then

            call PTOPEN(' ','fitlyman.plt',0,imode) ! plot
            call AGVERS()
          else
             call PTOPEN(' ','fitlyman.plt',1,imode)    ! overplot     

                               ! Draws wavelengths of upper spectra
            ry= i * GR_DEL - GR_BOT - GR_TOP - 1
            rx=REGMIN(i)
            write(cwave,'(f5.0)')rx
            call PTTEXT(cwave,rx,ry,0.,.9,2)
            rx=REGMAX(i)
            write(cwave,'(f5.0)')rx
            call PTTEXT(cwave,rx,ry,0.,.9,1)
          end if
          if (NREGIO .eq. 1) then
             call PTFRAM(xdef,ydef,'Wavelength (A)','Normalized Flux') 
          else
             if (i.eq.1) call PTFRAM(xdef,ydef,'Wavelength (A)',' ') 
          end if

          GR_OFF = i * GR_DEL - GR_BOT - GR_TOP - 1.
          rgroff = GR_OFF

          colore(1)=GCOLOR(2)
          call PTKWRI('COLOUR',1,colore)
          call PTDATA(5,1,0,GraWav,CONTIN,rgroff,nptgra)
          call PTDATA(5,3,0,GraWav,SKYLEV,rgroff,nptgra)   

          colore(1)=GCOLOR(1)
          call PTKWRI('COLOUR',1,colore)
          call PTDATA(5,1,1,GraWav,GraSpe,rgroff,nptgra)
          ry = rgroff +1 + GR_LAB
          rx = REGMIN(i) + (REGMAX(i) - REGMIN(i)) / 10
          colore(1)=GCOLOR(6)
          call PTKWRI('COLOUR',1,colore)
          call PTTEXT(LLABEL(i),rx,ry,0.,.7,1)
          colore(1)=GCOLOR(1)
          call PTKWRI('COLOUR',1,colore)

        end do


C!!!!!!!
C       Draws the fit
C!!!!!!!

C     reads fdummy.grf
          open(25,file='fdummy.grf',status='old',iostat=ipt,err=29)
          ipt=1
 122      read(25,*,end=123,err=123)grf_w(ipt),grf_f(ipt),grf_r(ipt)
          ipt=ipt+1
          goto 122
 123      continue
          ipt=ipt -1
          if (ipt.lt.1) goto29
          close(25)

          imode = -1
          do i = 1,NREGIO

            lim_x(1)=REGMIN(i)
            lim_x(2)=REGMAX(i)
            lim_y(1)= GR_YMN
            lim_y(2)= GR_YMX

            do k=1,NMXWRK
               grf_ww(k)=0.
               grf_rw(k)=0.
               grf_fw(k)=0.
            end do

            j = 1
            do k = 1,ipt
               if(grf_w(k).ge.REGMIN(i).and.grf_w(k).le.REGMAX(i))then
                  grf_ww(j)=grf_w(k)
                  grf_fw(j)=grf_f(k)
                  grf_rw(j)=grf_r(k)
                  j = j+1
               end if
            end do
            j = j-1
            call PTKWRR('XWNDL',2,lim_x)        
            CALL AGWDEF(lim_x(1),lim_x(2),lim_y(1),lim_y(2))

                                ! position of fitted lines
            if (NGRALI .gt. 0) then
               colore(1)=GCOLOR(7)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               dummy(1)=2.
               CALL PTKWRR('SSIZE',1,dummy)
               call PTDATA(13,0,0,LamCen,YLIN,rGRA,NGRALI)
               dummy(1)=1.
               CALL PTKWRR('SSIZE',1,dummy)
            end if

                                ! residuals
            if (I_RESD.and.j.gt.0) then
               colore(1)=GCOLOR(4)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               call PTDATA(5,1,1,grf_ww,grf_rw,rGRA,j)
            end if

                                ! fitted profile
            if (j.gt.0) then
               colore(1)=GCOLOR(3)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               call PTDATA(5,1,0,grf_ww,grf_fw,rGRA,j)
            end if
                                ! variance
            if (I_VAR) then
               colore(1)=GCOLOR(5)
               call PTKWRI('COLOUR',1,colore)
               call PTOPEN(' ','fitlyman.plt',1,imode) ! overplot
               GR_OFF = i * GR_DEL - GR_BOT -GR_TOP - 1.
               rgra = GR_OFF
               call PTDATA(5,1,0,GraWav,GraSig,rGRA,nptgra)
            end if

          end do

          colore(1)=GCOLOR(1)
          call PTKWRI('COLOUR',1,colore)

          lim_x(1)=REGMIN(1)
          lim_x(2)=REGMAX(1)
          lim_y(1)=GR_YMN
          lim_y(2)=GR_YMX
          
          call PTKWRR('XWNDL',2,lim_x)        
          CALL AGWDEF(lim_x(1),lim_x(2),lim_y(1),lim_y(2))


 29     continue


c        call PTKWRC('BINMODE','OFF')
        call PTKWRC('YFORMAT','AUTO') 

        return
        end




C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        SUBROUTINE GetXCr(wave,flux,z,vel,ierr)
C
C
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

        IMPLICIT NONE
        real WAVE,flux,riplot,vel
        integer ierr

        include 'MID_REL_INCL:fit_var.inc'

        integer ikey,ier,icur,iplot
        real lim_x(2),xcur,ycur
        character*40 strapp
        DOUBLE PRECISION scafac,z

        DOUBLE PRECISION c,zr,dcur
        parameter (c=2.997D5)

C        include 'MID_REL_INCL:fit_mid.inc'
        INTEGER MADRID(1)
        INCLUDE 'MID_INCLUDE:ST_DEF.INC'

        COMMON/VMR/MADRID

        INCLUDE 'MID_INCLUDE:ST_DAT.INC'


          if (I_Z) then
             xcur = 0.
             ycur=1.
             lim_x(1)=-VELRAN
             lim_x(2)=VELRAN
          else
             xcur = (REGMIN(1) + REGMAX(1)) / 2.
             ycur = 1.
             lim_x(1)=REGMIN(1)
             lim_x(2)=REGMAX(1)
          endif

         call PTKWRR('XWNDL',2,lim_x)        

         call PTOPEN(' ',' ',-1,1) 

         z=0.
         vel=0.
         if (I_Z) then
           call sttdis(' VELOCITY  REDSHIFT  WAVELENGTH    FLUX',0,ier)
         else
           call sttdis(' WAVELENGTH        FLUX',0,ier)
         ENDIF
         strapp = ' '
          call sttdis(strapp,0,ier)
         
         call PTGCUR(xcur,ycur,ikey,icur)

 100    continue

         CALL AGGPLM(xcur,ycur,1,4) ! Overplots a cross

          ierr=icur
          GR_DEL = 1 + GR_BOT + GR_TOP
          riplot = (ycur + GR_BOT) / GR_DEL
          iplot = riplot
          iplot = iplot + 1

          scafac=(REGMAX(iplot)-REGMIN(iplot))/(REGMAX(1)-REGMIN(1))
          flux = ycur - (iplot - 1)*GR_DEL 
          
          strapp = ' '
          if (I_Z) then
             dcur = xcur
             wave = (1.D0 + dcur/2./c)/(1.D0-dcur/2./c)*
     1            (REGMAX(iplot)+REGMIN(iplot))/2.D0
             zr = (2.*dcur + REDSH * dcur + 2.*c*REDSH)/(2.*c - dcur)
             write(strapp,'(f9.2,f10.6,2f10.2)')xcur,zr,wave,flux
             vel=xcur
             z=zr
          else
             wave = (xcur - REGMIN(1))*scafac + REGMIN(iplot)
             write(strapp,'(f9.2,5x,f9.2)')wave,flux
          endif

          call sttdis(strapp,0,ier)
          call PTGCUR(xcur,ycur,ikey,icur)

          if (ikey.ne.32) goto 100

        return
        end