C @(#)complyman.for	19.1 (ESO-IPG) 02/25/03 13:28:03
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

      program complyman
      
      implicit none
      include 'MID_REL_INCL:fit_var.inc'
      
                                ! copies of keywords coming from line
      character*80 SPESYS,LINSYS,SINSYS
      character SELSYS
      real NOISYS(3)

      character*14 Elem(1000)
      double precision lam_lin(1000),nh_lin(1000)
      double precision bk_lin(1000),bt_lin(1000)
      double precision outflx(NMXSPE),noivar(NMXSPE),inppxs(NMXSPE)
      double precision inpwav(NMXSPE),inpfwh(NMXSPE),fluxun(NMXSPE)
      double precision dnoi1, dnoi2
      double precision gasdev

      integer ist,i,ijk,iok,nj,j
      integer TabID,i1,i2,i3,i4,icol(20)
      logical la
      double precision d1
      character*2 des
      character*14 Enam

      
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     MIDAS SET-UP
C!!!!!!!
      call STSPRO('fitlyman')
      
      call STKRDC('SPECTAB',1,1,80,i,SPESYS,ijk,ijk,iok)
      if (iok.ne.0) then
         call ErrMsg('Error in reading keywords: check context..')
         call STSEPI() 
         stop
      end if   
      call STKRDC('OUTTAB',1,1,80,i,LINSYS,ijk,ijk,iok)
      if (iok.ne.0) then
         call ErrMsg('Error in reading keywords: check context..')
         call STSEPI() 
         stop
      end if   
      call STKRDC('SINTTAB',1,1,80,i,SINSYS,ijk,ijk,iok)
      if (iok.ne.0) then
         call ErrMsg('Error in reading keywords: check context..')
         call STSEPI() 
         stop
      end if   
      call STKRDC('INPUTC',1,1,1,i,SELSYS,ijk,ijk,iok)
      if (SELSYS.eq.'l') SELSYS = 'L'
      if (SELSYS.eq.'a') SELSYS = 'A'
      if (SELSYS.eq.'m') SELSYS = 'M'

      call STKRDR('INPUTR',1,3,i,NOISYS,ijk,ijk,iok)
      if (iok.ne.0) then
         call ErrMsg('Error in reading keywords: check context..')
         call STSEPI() 
         stop
      end if   

      

C!!!!!!!
C     Videata iniziale
C!!!!!!!
      call DisMsg('Loading data...')
      call DisMsg('Wait please..')
      
C!!!!!!!
C     Inizializzazione programma
C!!!!!!!
      call AtmRD(ist)            !read atompar.dat
      NPUNTI=NMXSPE
      call reaspe(SPESYS,NPUNTI)    !read spectrum
      if (NPUNTI.lt.1) goto 666
      do i=1,NPUNTI
         inpwav(i)=XLAMBD(i)
         inpfwh(i)=FWHM(i)
         inppxs(i)=PIXSIZ(i)
      end do
      
      call TBTOPN(LINSYS,F_IO_MODE,TabID,ist) !open line table
      call TBLSER(TabID,'Element',icol(1),ist)
      call TBLSER(TabID,'LamCen',icol(2),ist)
      call TBLSER(TabID,'ColDen',icol(3),ist)
      call TBLSER(TabID,'B',icol(4),ist)
      call TBLSER(TabID,'BTur',icol(5),ist)


      i1=1                      !doesn't stop on errors
      i2=0
      i3=0
      call STECNT('PUT',i1,i2,i3)

      call TBIGET(TABID,i1,i2,i3,i4,nj,ist) !n. of rows


      j=0
      do i=1,nj
         call TBSGET(TABID,i,la,ist) ! selected?
         if (la) then
            call TBERDC(TABID,i,icol(1),ENam,la,ist) !reads element
            if (.not.la) then
               la = .false.
               des=Enam(1:2)
               if (SELSYS.eq.'A') then
                  la =.true.
               else
                  if (des.eq.'Ly'.and.SELSYS.eq.'L') la =.true.
                  if (des.ne.'Ly'.and.SELSYS.eq.'M') la =.true.
               end if
               if (la) then
                  j=j+1
                  Elem(j)=ENam
                  call TBERDD(TABID,i,icol(2),d1,la,ist)
                  if (la) then
                     j=j-1
                     goto 1011
                  end if
                  lam_lin(j)=d1

                  call TBERDD(TABID,i,icol(3),d1,la,ist)
                  if (la) then
                     j=j-1
                     goto 1011
                  end if
                  nh_lin(j)=d1

                  call TBERDD(TABID,i,icol(4),d1,la,ist)
                  if (la) then
                     j=j-1
                     goto 1011
                  end if
                  bk_lin(j)=d1

                  call TBERDD(TABID,i,icol(5),d1,la,ist)
                  if (la) then
                     j=j-1
                     goto 1011
                  end if
                  bt_lin(j)=d1

               end if

            end if

         end if
 1011    continue
      end do

      call TBTCLO(TABID,ist)


c     computes new spectrum

      call DisMsg('...done')
      call DisMsg('Computing spectrum....')

      do i=1,j
         bk_lin(i) = dsqrt(bk_lin(i)*bk_lin(i)+bt_lin(i)*bt_lin(i))
      end do

      d1=0.                     ! tau gunn-peterson. = 0
      call CHECK(d1,Elem,lam_lin,bk_lin,NH_lin,inpfwh,j,inpwav,inppxs,
     1    NPUNTI,outflx,fluxun)


c     add noise

                                ! NOISYS(1)= noise of sky, (stdev)
                                ! NOISYS(2)= total noise 
                                ! NOISYS(3)=idum
      ijk=NOISYS(3)
      if (ijk .gt.0) ijk=-ijk
      dnoi2 = abs(NOISYS(2)*NOISYS(2)-NOISYS(1)*NOISYS(1))
      dnoi1 =NOISYS(1)*NOISYS(1)
      do i=1,npunti
         noivar(i)=dsqrt(dnoi2*outflx(i) + dnoi1)
         outflx(i)=outflx(i) +noivar(i)*gasdev(ijk)
      end do

C     Saves new spectrum

      call TBTINI(SINSYS,F_TRANS,F_O_MODE,5,10,TABID,ist)
      call TBCINI(TabID,D_R8_FORMAT,1,'F8.3',' ','WAVE',icol(1),ist)
      call TBCINI(TabID,D_R8_FORMAT,1,'F8.5',' ','NORMFLUX',icol(2),
     1     ist)
      call TBCINI(TabID,D_R8_FORMAT,1,'F6.3',' ','FWHM',icol(3),ist)
      call TBCINI(TabID,D_R8_FORMAT,1,'F8.5',' ','STDEV',icol(4),ist)
      call TBCINI(TabID,D_R8_FORMAT,1,'F8.5',' ','FLUX_UNCONV',icol(5)
     1     ,ist)
      
      do i=1,NPUNTI
         call TBEWRD(TABID,i,icol(1),inpwav(i),ist)
         call TBEWRD(TABID,i,icol(2),outflx(i),ist)
         call TBEWRD(TABID,i,icol(3),inpfwh(i),ist)
         call TBEWRD(TABID,i,icol(4),noivar(i),ist)
         call TBEWRD(TABID,i,icol(5),fluxun(i),ist)
      end do

      call TBTCLO(TABID,ist)


      call DisMsg('...done')

 666  continue
      call STSEPI()
      end









