C @(#)cloudabs.for	19.1 (ESO-DMD) 02/25/03 13:23:29
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

       PROGRAM CLABS
C*************************************************************
C        M.PIERRE 1-88                                       *
C                                                            *
C       CALCULATION OF A SYNTHTETIC ABSORPTION SPECTRUM      *
C                                                            * 
C              50 ABSORPTION LINES MAX                       *       
C           VOIGTIAN OR MAXWELLIAN PROFILES                  *
C                                                            *
C       MAIN PROGRAM :ABS                                    *
C       SUBROUTINES  :ABSCOF absorption coefficient          *
C                   :IMAGE  final image                      *
C                   :CONVL1 1D convolution                   * 
C       Modified - START, STEP in/out in double precision    *
C*************************************************************
C
       IMPLICIT NONE
C 
       CHARACTER*60 TBLP,TBLC,IMA,NOMP,NOMI,NOMO,NOMOC
       CHARACTER*72 IDENTI,IDENTO,IDENTP
       CHARACTER*48 CUNIT
C 
       INTEGER    I,NPIXI(1),NPIXO(1),IST,NCARA
       INTEGER    NPIXP
       INTEGER    NPARAM,NCLD,NAXIS,IPROF
       INTEGER    KNUL,KUN,IMNP,IMNI,IMNO,IMNOC,NAC,NAR
       INTEGER    MADRID, TIDP, TIDC, I1
       INTEGER*8  PNTRO,PNTRI,PNTROC,PNTRP
C 
       REAL Z
       REAL ABSA(500),ABSBL(500),ABSC(500)
       REAL START,STEPI,STEPO,ALAM(5000)
       REAL STARTP,STEPP,X1,X2,X3,X4
       REAL PROFPS(5000),YO(5000)
C 
       DOUBLE PRECISION DSTA,DSTP
C 
       LOGICAL NULL,NULL1,NULL2,NULL3,NULL4,ISEL
C
C DEFINITION DES PARAMETRES ATOMIQUES
C       50 transitions max , multiplicite max 6
C       valeur actuelle du nb de transitions :NTR
       INTEGER Q1(500,6),QK(500,6),MUL(500),NAT(500)
       REAL    WAVEO(500,6),SAKI(500,6),F1K(500,6),XNAT(500)
       INCLUDE 'MID_INCLUDE:TABLES.INC'
       COMMON/VMR/MADRID(1)
       INCLUDE 'MID_INCLUDE:TABLED.INC'
C
C       
       CALL STSPRO('CLABS')
C
C******** ENTREE DES DONNEES  ***********
C
       CALL STKRDC('IN_A',1,1,60,NCARA,NOMI, KUN,KNUL,IST)
       CALL STKRDC('IN_B',1,1,60,NCARA,IMA,  KUN,KNUL,IST)
       CALL STKRDC('IN_C',1,1,60,NCARA,TBLC, KUN,KNUL,IST)
       CALL STKRDC('IN_D',1,1,60,NCARA,TBLP, KUN,KNUL,IST)
       CALL STKRDC('IN_E',1,1,60,NCARA,NOMP, KUN,KNUL,IST)
       CALL STKRDR('CLDZ',  1,1, NCARA,Z,    KUN,KNUL,IST)
       CALL STKRDI('CLDOP', 1,1, NCARA,IPROF,KUN,KNUL,IST)
C
       CALL STIGET(NOMP,10,0,1,1,NAXIS,NPIXP,DSTA,DSTP,
     :  IDENTP,CUNIT,PNTRP,IMNP,IST)
       STARTP = DSTA
       STEPP  = DSTP
C
C******** INITIALISATION DE L'ENVIRONEMENT MIDAS *******
C
       CALL TBTOPN(TBLP,F_I_MODE,TIDP,IST)
       CALL TBTOPN(TBLC,F_U_MODE,TIDC,IST)
       CALL TBIGET(TIDP,I,NPARAM,I,NAC,NAR,IST)
       CALL TBIGET(TIDC,I,NCLD,I,NAC,NAR,IST)
C
C******** LECTURE DES FICHIERS ********
C
C       IMAGE:
       CALL STIGET(NOMI,10,0,1,1,NAXIS,NPIXI,DSTA,DSTP,
     :       IDENTI,CUNIT,PNTRI,IMNI,IST)
       START = DSTA
       STEPI = DSTP
C
       IF((STEPI.GT.STEPP).OR.(STEPI.LT.STEPP)) THEN
          CALL STTPUT(
     :    'ERROR :IMAGE AND PSF HAVE DIFFERENT STEPSIZE',IST)
          GO TO 1111
       ENDIF
C
C       Absorption:
       I1 = 0
       DO 10 I=1,NCLD
          CALL TBSGET(TIDC,I,ISEL,IST)
          IF (ISEL) THEN
             CALL TBERDR(TIDC,I,2,X1,NULL1,IST)
             CALL TBERDR(TIDC,I,3,X2,NULL2,IST)
             CALL TBERDR(TIDC,I,4,X3,NULL3,IST)
             CALL TBERDR(TIDC,I,5,X4,NULL4,IST)
             IF (NULL1 .OR. NULL2 .OR. NULL3 .OR. NULL4) GOTO 5
               I1        = I1 + 1
               ABSA(I1)  = X1
               ABSBL(I1) = X2
               ABSC(I1)  = X3
               XNAT(I1)  = X4
               NAT(I1)   = INT(XNAT(I1))
5            CONTINUE
          ENDIF
10     CONTINUE
       NCLD = I1
C
C       TABLE DES PARAMETRES PHYSIQUES
       DO 20 I=1,NPARAM
          CALL TBERDI(TIDP,I,2, MUL(I),    NULL,IST)
          CALL TBERDR(TIDP,I,3, WAVEO(I,1),NULL,IST)
          CALL TBERDR(TIDP,I,4, WAVEO(I,2),NULL,IST)
          CALL TBERDI(TIDP,I,5, Q1(I,1),   NULL,IST)
          CALL TBERDI(TIDP,I,6, Q1(I,2),   NULL,IST)
          CALL TBERDI(TIDP,I,7, QK(I,1),   NULL,IST)
          CALL TBERDI(TIDP,I,8, QK(I,2),   NULL,IST)
          CALL TBERDR(TIDP,I,9, SAKI(I,1), NULL,IST)
          CALL TBERDR(TIDP,I,10,SAKI(I,2), NULL,IST)
          CALL TBERDR(TIDP,I,11,F1K(I,1),  NULL,IST)
          CALL TBERDR(TIDP,I,12,F1K(I,2),  NULL,IST)
20     CONTINUE
C 
C******* CALCUL DE L'IMAGE DE SORTIE *********
C
       STEPO=STEPI
       NPIXO(1)=NPIXI(1)
       DO 30 I=1,NPIXO(1)
          PROFPS(I)=START+(FLOAT(I)-1.)*STEPO
30     CONTINUE
       I=INDEX(IMA,' ')-1
       IF (IPROF.GT.0.) THEN 
         NOMO=IMA(1:I)//'1'
         NOMOC=IMA(1:I)//'1c'
C        NOMO='SYN1'
C        NOMOC='SYN1C'
         IDENTO='MAXWELL + DAMPING'
       ELSE
         NOMO=IMA(1:I)//'0'
         NOMOC=IMA(1:I)//'0c'
         IDENTO='MAXWELL SEUL'
       ENDIF
C
       DSTA = START
       DSTP = STEPO
       CALL STIPUT(NOMO,10,1,1,NAXIS,NPIXO,DSTA,DSTP,IDENTO,CUNIT,
     :        PNTRO,IMNO,IST)
       CALL STIPUT(NOMOC,10,1,1,NAXIS,NPIXO,DSTA,DSTP,IDENTO,CUNIT,
     :        PNTROC,IMNOC,IST)
C
C       Coefficient d'absorption
       CALL ABSCOF(TIDC,ABSC,ABSBL,ABSA,
     :       Z,NCLD,IPROF,NAT,MUL, 
     :       NPIXO(1),PROFPS,ALAM,
     :       WAVEO,Q1,QK,F1K,SAKI)       
C
C       Image finale
       CALL IMAGE(MADRID(PNTRI),MADRID(PNTRO),ALAM,NPIXO(1))
       CALL IMAGE(MADRID(PNTRI),YO,ALAM,NPIXO(1))
       CALL CONVL1(YO,NPIXO(1),MADRID(PNTRP),NPIXP,MADRID(PNTROC))
C      CALL PUTKEY_ST('NOM','C',NOMO,1,60,IST)
C
C****** LIBERATION DES FICHIERS ET SORTIE DE L'ENVIRONEMENT MIDAS
       CALL STTPUT('   ',IST)
       CALL STTPUT('RESULTING ABSORPTION SPECTRA :',IST)
       CALL STTPUT(NOMO,IST)
       CALL STTPUT(NOMOC,IST)
       CALL TBTCLO(TIDC,IST)
       CALL TBTCLO(TIDP,IST)
1111   CALL STSEPI
       STOP
       END
C
C********************************
C******* SOUS PROGRAMMES ********
C********************************
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
       SUBROUTINE ABSCOF (TIDC,B,VR,WEI,
     :  Z,NCLD,IPROF,NAT,MUL,
     :  NPTS,PROFPS,ALAM,
     :  WLO,Q1,QK,F1K,SAKI)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C+
C        Subroutine to compute the absorption
C        coefficient for an array of given
C        wavelengths.
C
C        Given;
C        B - Array of "B" values for cloud model.
C         (VR - Array of radial velcities for cloud model.)
C  ICI VR  position observee de l'absorption (en Angstroms)
C        WEI - Array of weights for cloud model.
C        WLO  :position centrale  des raies
C        F1K - F value.
C        QK
C        Q1
C        NAT :nature de la transition
C        MUL :multiplicite de la transition
C        AKI :valeurs des Aki
C        SAKI:some des Aki
C        NPTS - No. of pts. in profile.
C        PROFPS - Array of wavelength positions.
C        MAXPROF - Size of arrays PROFPS & ALAM.
C
C        Returned;
C        ALAM - Absorption coeff. evaluated at each wavelength.
C
C        Function called: VOIGT
C
C        A C Davenhall. /ROE/                            28/10/81.
C
C        Based on a routine by A C Davenhall /UCL/      Spring 1977.
C
C-
       IMPLICIT NONE
       INTEGER  TIDC
       REAL     PROFPS(5000),ALAM(5000)
       REAL     WAVE1,WAVE4,WAVEI,WAVE
       INTEGER  NPTS,I,J,K,N
       REAL     B(500),VR(500),WEI(500),ZA(500)
       INTEGER  NCLD
       REAL     CO1,CO2,C,DUMMY,Z
       REAL     DELTAM,VNU,BB,WW,AA
C
       INTEGER IST,IPROF
       INTEGER QK(500,6),Q1(500,6),MUL(500),NAT(500)
       REAL WLO(500,6),SAKI(500,6),F1K(500,6)
       REAL AK1(500,6),CONST(500,6),DELTA(500,6)
       REAL LAM,BETA,R,VITES,XX1,XX2,CONV
       REAL VOIGT
C
C        Speed of light in Km/sec.
C
        DATA C/2.9979251E5/
C
       DATA CO1/2.654418269E-20/
       DATA CO2/7.487975686E-37/
       DATA CONV/1.499E-16/
C
C
       DO 20 J=1,NCLD
          N=NAT(J)
          LAM=VR(J) !voir definition de VR 
          ZA(J)=LAM/WLO(N,1)-1.
          R=(1.+Z)/(1.+ZA(J))
          BETA=(R**2-1.)/(R**2+1.)
          VITES=BETA*C
          CALL TBEWRR(TIDC,J,6,VITES,IST)       
          CALL TBEWRR(TIDC,J,7,ZA(J),IST)       
          DO 10 K=1,MUL(N)                    ! constantes atomiques
            AK1(N,K)=F1K(N,K)*FLOAT(Q1(N,K))
     :       /(FLOAT(QK(N,K))*(WLO(N,K)**2)*CONV)
            CONST(N,K)=CO2*AK1(N,K)*FLOAT(QK(N,K))/FLOAT(Q1(N,K))
            DELTA(N,K)=CO1*SAKI(N,K)
10        CONTINUE
20     CONTINUE
C*        

       DO 70 I=1,NPTS             !calcul du coefficient d'absorption
         DUMMY=0.0E0
         WAVEI=PROFPS(I)

C
         IF (IPROF.EQ.1) THEN
C CAS MAXWELL+ DAMPING WINGS
            DO 40  J=1,NCLD
               N=NAT(J)
               WW=WEI(J)
               WAVE=WAVEI/(1.+ZA(J))
               WAVE4=WAVE**4
               DO 30 K=1,MUL(N)
                   WAVE1=C/WLO(N,K)
                   BB=WAVE1/B(J) 
                   DELTAM=DELTA(N,K)*(WAVE**2)
                   AA=DELTAM*BB
                   VNU=(WAVE-WLO(N,K))*BB
                   DUMMY=DUMMY+
     :                   (WW*CONST(N,K)*WAVE4*VOIGT(VNU,AA)*BB)
30             CONTINUE
40          CONTINUE
         ELSE
C CAS MAXWELL SEUL
            DO 60 J=1,NCLD
               WW=WEI(J)
               BB=B(J)
               N=NAT(J)       
               WAVE=WAVEI/(1.+ZA(J))
               DO 50 K=1,MUL(N) 
                  XX1=WW*WLO(N,K)/BB*(1.495E-15)*F1K(N,K)
                  XX2=C*(WLO(N,K)-WAVE)/WAVE/BB
                  DUMMY=DUMMY+XX1*EXP(-XX2*XX2)
50             CONTINUE
60          CONTINUE
         ENDIF
C
         ALAM(I)=DUMMY
70     CONTINUE
C
       RETURN
       END
C
C***************************************************************
C
       REAL FUNCTION VOIGT(V,A)
C
C        Subroutine to evaluate a Voigtian.
C
C        Obtained from; S L Wright /UCL/              Spring 1977.
C
      IMPLICIT NONE
C
      REAL   V,A
      INTEGER N,N1
      REAL    V0,V1,V2
      REAL   H0(41),H1(81),H2(41)
      DATA H0/
     1 1.0000000, 0.9900500, 0.9607890, 0.9139310, 0.8521440, 0.7788010,
     2 0.6976760, 0.6126260, 0.5272920, 0.4448580, 0.3678790, 0.2984970,
     3 0.2369280, 0.1845200, 0.1408580, 0.1053990, 0.0773050, 0.0555760,
     4 0.0391640, 0.0270520, 0.0183156, 0.0121552, 0.0079071, 0.0050418,
     5 0.0031511, 0.0019305, 0.0011592, 0.0006823, 0.0003937, 0.0002226,
     6 0.0001234, 0.0000671, 0.0000357, 0.0000186, 0.0000095, 0.0000048,
     7 0.0000024, 0.0000011, 0.0000005, 0.0000002, 0.0000001/
      DATA H1/
     1-1.1283800,-1.1059600,-1.0404800,-0.9370300,-0.8034600,-0.6494500,
     2-0.4855200,-0.3219200,-0.1677200,-0.0301200, 0.0859400, 0.1778900,
     3 0.2453700, 0.2898100, 0.3139400, 0.3213000, 0.3157300, 0.3009400,
     4 0.2802700, 0.2564800, 0.2317260, 0.2075280, 0.1848820, 0.1643410,
     5 0.1461280, 0.1302360, 0.1165150, 0.1047390, 0.0946530, 0.0860050,
     6 0.0785650, 0.0721290, 0.0665260, 0.0616150, 0.0572810, 0.0534300,
     7 0.0499880, 0.0468940, 0.0440980, 0.0415610, 0.0392500, 0.0351950,
     8 0.0317620, 0.0288240, 0.0262880, 0.0240810, 0.0221460, 0.0204410,
     9 0.0189290, 0.0175820, 0.0163750, 0.0152910, 0.0143120, 0.0134260,
     A 0.0126200, 0.0118860, 0.0112145, 0.0105990, 0.0100332, 0.0095119,
     B 0.0090306, 0.0085852, 0.0081722, 0.0077885, 0.0074314, 0.0070985,
     C 0.0067875, 0.0064967, 0.0062243, 0.0059688, 0.0057287, 0.0055030,
     D 0.0052903, 0.0050898, 0.0049006, 0.0047217, 0.0045526, 0.0043924,
     E 0.0042405, 0.0040964, 0.0039595/
      DATA H2/
     1 1.0000000, 0.9702000, 0.8839000, 0.7494000, 0.5795000, 0.3894000,
     2 0.1953000, 0.0123000,-0.1476000,-0.2758000,-0.3679000,-0.4234000,
     3-0.4454000,-0.4392000,-0.4113000,-0.3689000,-0.3185000,-0.2657000,
     4-0.2146000,-0.1683000,-0.1282100,-0.0950500,-0.0686300,-0.0483000,
     5-0.0331500,-0.0222000,-0.0145100,-0.0092700,-0.0057800,-0.0035200,
     6-0.0021000,-0.0012200,-0.0007000,-0.0003900,-0.0002100,-0.0001100,
     7-0.0000600,-0.0000300,-0.0000100,-0.0000100, 0.0000000/
      V=ABS(V)
      V0=V*10.E0
      N=V0
      IF(N.LT.40) GO TO 1
      IF(N.LT.120)GO TO 2
      VOIGT=(0.56419E0+0.846E0/(V*V))/(V*V)*A
      RETURN
    1 V1=N
      N=N+1
      V2=V0-V1
      N1=N+1
      VOIGT=V2*(H0(N1)-H0(N)+A*(H1(N1)-H1(N)+A*(H2(N1)-H2(N))))+
     1 H0(N)+A*(H1(N)+A*H2(N))
      RETURN
    2 N=N/2+20
      V1=(N-20)*2
      N=N+1
      V2=(V0-V1)/2.E0
      N1=N+1
      VOIGT=A*((H1(N1)-H1(N))*V2+H1(N))
      RETURN
      END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
       SUBROUTINE IMAGE(YI,YO,ALAM,NPIXO)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C       Calcul de l'image synthetique + absorption
C
       IMPLICIT NONE
C
       INTEGER NPIXO
       REAL YO(1),YI(1),ALAM(5000)
       INTEGER I
C
       DO 10 I=1,NPIXO
         YO(I)=YI(I)*EXP(-ALAM(I))   
10     CONTINUE
       RETURN 
       END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE CONVL1(F,NF,PSF,NP,C)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION:
C  subroutine CONVL1        version 2.00       840116
C  P. GROSBOL               ESO - GARCHING
C  K. Banse (streamlined for MIDAS)
C
C.KEYWORDS:
C  CONVOLUTION
C
C.PURPOSE:
C  CONVOLES TWO ONE DIMENSIONAL ARRAYS WITH EACH OTHER
C
C.ALGORITHM:
C  A ONE DIMENSIONAL ARRAY IS CONVOLED WITH A ONE DIMENSIONAL
C  POINT SPREAD FUNCTION. OUTSIDE THE LIMITS OF THE FUNCTION
C  ARRAY THE FIRST/LAST VALUE IS TAKEN.
C
C.INPUT/OUTPUT:
C  THE SUBROUTINE IS CALLED AS
C      CALL CONVL1(F,NF,PSF,NP,C)
C  WHERE THE PARAMETERS ARE
C  F(NF)   :  FUNCTION ARRAY                          (INPUT)
C  NF      :  DIMENSION OF FUNCTION ARRAY             (INPUT)
C  PSF(NP) :  POINT SPREAD FUNCTION ARRAY             (INPUT)
C  NP      :  DIMENSION OF POINT SPREAD FUNCTON       (INPUT)
C  C(NF)   :  CONVOLVED RESULT ARRAY                  (OUTPUT)
C
C  NOTE :  THE RESULT IS ADDED TO THE OUTPUT ARRAY 'C'.
C
C---------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER NF, NP
      REAL F(NF),C(NF),PSF(NP)
C
      INTEGER NOS, NL1, NL2, NL3, I, NR, K
      REAL    R, FACT
C
C  initialize
      NOS = NP/2
      NL1 = NOS + 1
      NL2 = NF - NOS
      NL3 = NL2 + 1
C
C  if NP = 1, we only use middle section
       IF (NOS.EQ.0) THEN
          DO 20 I = NL1,NL2
             NR = I - NOS
             R = 0.
             DO 10 K = 1,NP
               R = R + PSF(K)*F(NR)
               NR = NR + 1
10           CONTINUE
             C(I) = R
20        CONTINUE
          RETURN
       ENDIF
C       
C  for NP > 1, handle 1. section, where values before F(1) would be needed
       FACT = F(1)
       DO 50 I = 1,NOS
           R = 0.
          DO 30 K=1,1+NOS-I
             R = R + PSF(K)*FACT
30        CONTINUE
          NR = 1
           DO 40 K=2+NOS-I,NP
              R = R + PSF(K)*F(NR)
              NR = NR + 1
40         CONTINUE
           C(I) = R
50     CONTINUE
C
C  in the middle no problems
       DO 70 I = NL1,NL2
          NR = I - NOS
          R = 0.
          DO 60 K = 1,NP
             R = R + PSF(K)*F(NR)
             NR = NR + 1
60        CONTINUE
          C(I) = R
70     CONTINUE
C
C  last section, values after F(NF) would be needed
       FACT = F(NF)
       DO 100 I = NL3,NF
          NR = I - NOS
           R = 0.
          DO 80 K=1,NF+NOS-I
             R = R + PSF(K)*F(NR)
              NR = NR + 1
80        CONTINUE
          DO 90 K=NF+NOS-I+1,NP
              R = R + PSF(K)*FACT
90        CONTINUE
          C(I) = R
100    CONTINUE
C
       RETURN
       END
