C @(#)cloudabs.for 17.1.1.1 (ESO-DMD) 01/25/02 17:13:49 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