C @(#)necompix.for	17.1.1.1 (ESO-DMD) 01/25/02 17:51:30
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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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  Program COMPIXEL.FOR
C
C  Computes the pixelsize of an Echelle session
C
C
C  By A.Fontana
C  Version 1.0   Sun&IBM 01-01-95
C  Version 1.1   Sun&IBM 20-02-95 YFORMAT ok.
C  Version 1.2   S.Wolf  05-10-99 return minimum and maximum
c                                 pixelsize in OUTPUTD(1..2),
c                                 plot flag added (inputi(3)),
c                                 use {LINTAB} instead of line.tbl.
C!!!!!!!
C     MIDAS SET-UP
C!!!!!!!

       implicit none
       integer NMXORD,NMXPIX
       parameter (NMXORD=100,NMXPIX=5000)
       double precision wave(NMXORD,NMXPIX)
       double precision pixsiz(NMXORD,NMXPIX)
       double precision coefd(7*NMXORD)
       double precision wavmin,wavmax,pixmin,pixmax,xj

       integer ord(2),TabID,coenum,xsize

       integer iread,inul,istat,i,j,k

       real xdef(4),ydef(4)                 ! graphic
       real   dummy(2),XLAMBD(NMXPIX),SPECTR(NMXPIX)
       integer colore(4),imode,stype,ltype,plt
       character*60 LINTAB
c       character*70 stringa


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

C      Initialize midas
       call STSPRO('compixel')

C      Reads ord1,ord2 from INPUTI and x-axis size from IMSIZE
       call STKRDI('INPUTI',1,2,iread,ord,inul,inul,istat)
       call STKRDI('INPUTI',3,1,iread,plt,inul,inul,istat)
       call STKRDI('IMSIZE',1,1,iread,xsize,inul,inul,istat)
       call STKRDC('LINTAB',1,1,60,iread,LINTAB,inul,inul,istat)
 
C      Opens line.tbl and reads coefd
       call TBTOPN(LINTAB,F_D_MODE,TabID,istat)
       call STDRDD(TabID,'COEFD',1,7*NMXORD,coenum,coefd,
     1       inul,inul,istat)
       call TBTCLO(TabID,istat)

C      Computes wavelengths
       do i = ord(1), ord(2)
         do j = 1,xsize
             wave(i,j) = coefd(7*i -6)
             xj = 1
             do k = 2,7
               xj = xj * dfloat(j)
               wave(i,j)=wave(i,j)+coefd(7*i-7+k)*xj
             end do
             wave(i,j)=wave(i,j) * 1.D4
c debug
c             if (j.gt.1250.and.j.lt.1350) then
c                write(stringa,*)j,wave(i,j)
c                call sttdis(stringa,o,istat)
c             end if
c  endebug
         end do
       end do

C      Computes pixelsizes
       do i = ord(1), ord(2)
         do j = 1,xsize-1
             pixsiz(i,j) = wave(i,j+1) - wave(i,j)
         end do
         pixsiz(i,xsize) = pixsiz(i,xsize-1)
       end do

C      Looks for maxima and minima
       wavmin = wave(Ord(1),1)
       wavmax = wave(Ord(2),xsize)

       pixmin = pixsiz(Ord(1),xsize)  ! should be the right one...
       pixmax = pixsiz(Ord(2),1)      ! should be the right one...
                                      ! but I check....
       do i = ord(1), ord(2)
         do j = 1,xsize
           if (pixmin .gt. pixsiz(i,j)) pixmin = pixsiz(i,j)
           if (pixmax .lt. pixsiz(i,j)) pixmax = pixsiz(i,j)
         end do
       end do

       if (plt.eq.0) goto 10
c%%%%%%% 
C        Graphic set-up
C%%%%%%%


       CALL PTKWRR('SCALES',0,dummy)
       CALL PTKWRR('OFFS',0,dummy)


       xdef(1)=wavmin - (wavmax - wavmin) / 10.
       xdef(2)=wavmax + (wavmax - wavmin) / 10.
       xdef(3)=(xdef(2) - xdef(1))/5.
       xdef(4)=xdef(3)/5.


       ydef(1)= pixmin - (pixmax - pixmin ) / 10.
       ydef(2)= pixmax + (pixmax - pixmin ) / 10.
       ydef(3) = (ydef(2) - ydef(1))/5
       ydef(4)=  ydef(3) / 5.

       call PTKWRR('XWNDL',4,xdef)
       call PTKWRR('YWNDL',4,ydef)
       call PTKRDI('STYPE',1,inul,stype)
       call PTKRDI('LTYPE',1,inul,ltype)
       imode = -1
       colore(1)=1
       call PTKWRI('COLOUR',1,colore)

C     Makes the plots

       do i = Ord(1),Ord(2)

         if (i.eq.Ord(1))then
           call PTOPEN(' ','compixel.plt',0,imode) ! plot
           call AGVERS()
           call PTFRAM(xdef,ydef,'Wavelength (A)','Pixel Size (A)') 
         else
           call PTOPEN(' ','fitlyman.plt',1,imode)    ! overplot     
         end if
  
         do j = 1,xsize
           XLAMBD(j) = wave(i,j)
           SPECTR(j) = pixsiz(i,j)
         end do

         if (colore(1) .eq. 2 ) then 
           colore(1)=1
         else
           colore(1)=2
         end if
         call PTKWRI('COLOUR',1,colore)
         call PTDATA(stype,ltype,0,XLAMBD,SPECTR,0.,xsize)
         
       end do

       colore(1)=1
       call PTKWRI('COLOUR',1,colore)
C
10     coefd(1) = pixmin
       coefd(2) = pixmax
       CALL STKWRD('OUTPUTD',coefd,1,2,inul,istat)

       call STSEPI() 
       end