216 lines
6.4 KiB
Fortran
216 lines
6.4 KiB
Fortran
SUBROUTINE LINPRO(ITR,ID,PRF)
|
|
C =============================
|
|
C
|
|
C Line profile coefficient for the line ITR;
|
|
C for "classical" lines (i.e. those which are not represented by ODF's).
|
|
C It is either calculated for each depth (if LCOMP=true), or is
|
|
C just taken as constant, depth-independent quantity, already
|
|
C calculated in START and stored in array PROF
|
|
C
|
|
C Input: ITR - index of transition
|
|
C ID - depth index
|
|
C Output: PRF - array of absorption profile
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ODFPAR.FOR'
|
|
common/quasun/tqmprf,iquasi,nunalp,nunbet,nungam,nunbal
|
|
PARAMETER (BOL2=2.76108D-16, CIN=UN/2.997925D10, OS0=0.02654,
|
|
* F0C1=1.25D-9, TTW=2./3., PISQ1=UN/1.77245385090551D0,
|
|
* C18IN=UN/CAS,F0C2=3.906D-11,cca=2.997925e18,
|
|
* CINV=UN/cca,AL10=2.3025851)
|
|
DIMENSION PRF(MFREQL),PRF0(MHWL),prfb(mfreql),prfr(mfreql)
|
|
C
|
|
MODE=IABS(INDEXP(ITR))
|
|
IJ0=IFR0(ITR)
|
|
IJ1=IFR1(ITR)
|
|
INTM0=INTMOD(ITR)
|
|
INTM=IABS(INTM0)
|
|
IP=IABS(IPROF(ITR))
|
|
C
|
|
C Doppler width
|
|
C
|
|
IF(LCOMP(ITR).OR.IP.EQ.2) THEN
|
|
IAT=IATM(ILOW(ITR))
|
|
AM=BOL2/AMASS(IAT)*TEMP(ID)
|
|
DOP=FR0(ITR)*CIN*SQRT(AM+VTURBS(ID)*VTURBS(ID))
|
|
DOP1=UN/DOP
|
|
END IF
|
|
IF(LCOMP(ITR).AND.IP.EQ.1.OR.IP.LT.0) THEN
|
|
CALL DOPGAM(ITR,ID,TEMP(ID),DOP,AGAM)
|
|
DOP1=UN/DOP
|
|
END IF
|
|
S=OSC0(ITR)*OS0
|
|
XNORM=PISQ1*S*DOP1
|
|
C
|
|
C Depth-independent profile
|
|
C
|
|
IF(ISPODF.EQ.0) THEN
|
|
DO IJ=IJ0,IJ1
|
|
PRF(IJ-IJ0+1)=PROF(IJ)
|
|
END DO
|
|
ELSE
|
|
DO IJ=KFR0(ITR),KFR1(ITR)
|
|
PRF(IJ-KFR0(ITR)+1)=PROF(IJ)
|
|
END DO
|
|
END IF
|
|
IF(.NOT.LCOMP(ITR)) RETURN
|
|
c
|
|
c evaluation of the depth-dependent profile
|
|
C
|
|
IF(IP.EQ.0) THEN
|
|
DO IJ=IJ0,IJ1
|
|
V=(FREQ(IJ)-FR0(ITR))*DOP1
|
|
IF(ABS(V).LE.13.) PRF(IJ-IJ0+1)=EXP(-V*V)*XNORM
|
|
END DO
|
|
ELSE IF(IP.EQ.1) THEN
|
|
DO IJ=IJ0,IJ1
|
|
V=(FREQ(IJ)-FR0(ITR))*DOP1
|
|
PRF(IJ-IJ0+1)=VOIGT(V,AGAM)*XNORM
|
|
END DO
|
|
ELSE IF(IP.EQ.2) THEN
|
|
C
|
|
C for IP=2 - approximate Stark profile for hydrogen lines,
|
|
C evaluate necessary frequency-independent quantities
|
|
C
|
|
ANE=ELEC(ID)
|
|
F000=EXP(TTW*LOG(ANE))
|
|
II=NQUANT(ILOW(ITR))
|
|
JJ=NQUANT(IUP(ITR))
|
|
IZZ=IZ(IEL(ILOW(ITR)))
|
|
FAC=TWO
|
|
F00=F000*F0C1
|
|
if(iquasi.gt.0.and.ii.eq.1) then
|
|
if(jj.eq.2) fac=un
|
|
if(jj.eq.3.and.iquasi.gt.1) fac=un
|
|
end if
|
|
IF(IZZ.EQ.2) THEN
|
|
FAC=UN
|
|
F00=F000*F0C2
|
|
END IF
|
|
CALL STARK0(II,JJ,IZZ,XKIJ0,WL00,FIJ0)
|
|
FXK=F00*XKIJ0
|
|
DBETA=WL00*WL00*C18IN/FXK
|
|
BETAD=DOP*DBETA
|
|
FID=OS0*FIJ0*DBETA
|
|
FID0=OS0*(OSC0(ITR)-FIJ0)*DOP1*PISQ1
|
|
CALL DIVSTR(IZZ)
|
|
C
|
|
C loop over frequencies
|
|
C
|
|
DO IJ=IJ0,IJ1
|
|
BETA=DBETA*ABS(FREQ(IJ)-FR0(ITR))
|
|
SG=STARKA(BETA,fac)*FID
|
|
SG0=0.
|
|
V=(FREQ(IJ)-FR0(ITR))*DOP1
|
|
IF(ABS(V).LE.13.) SG0=EXP(-V*V)*FID0
|
|
PRF(IJ-IJ0+1)=SG+SG0
|
|
sgmax=max(sgmax,sg+sg0)
|
|
END DO
|
|
ELSE IF(IP.EQ.3.or.ip.eq.4) THEN
|
|
C
|
|
C for IP=3 - exact Stark profile for hydrogen lines (Lemke tables)
|
|
C
|
|
II=NQUANT(ILOW(ITR))
|
|
JJ=NQUANT(IUP(ITR))
|
|
ANE=ELEC(ID)
|
|
F00=F0C1*EXP(TTW*LOG(ANE))
|
|
FID=OS0*OSC0(ITR)
|
|
anel=log10(ane)
|
|
tl=log10(temp(id))
|
|
c
|
|
c switch to either original Lemke/Tremblay of Xenomorph
|
|
c
|
|
if(ilxen(ii,jj).eq.0.or.anel.lt.xnemin) then
|
|
c
|
|
c original Lemke/Tremblay
|
|
c
|
|
WLI0=cca/fr0(itr)
|
|
fac=un
|
|
if(iquasi.gt.0.and.ii.eq.1) then
|
|
if(jj.eq.2) fac=half
|
|
if(jj.eq.3.and.iquasi.gt.1) fac=half
|
|
end if
|
|
iline=0
|
|
if(ip.eq.3) then
|
|
IF(II.LE.4.AND.JJ.LE.22) iline=ilinh(ii,jj)
|
|
else
|
|
IF(II.LE.2.AND.JJ.LE.10) iline=ilinh(ii,jj)
|
|
end if
|
|
c
|
|
if(iline.gt.0) then
|
|
CALL INTLEM(PRF0,WLI0,ILINE,ID)
|
|
NWL=NWLHYD(ILINE)
|
|
DO IJ=IJ0,IJ1
|
|
al=abs(wli0-cca/freq(ij))
|
|
IF(AL.LT.1.E-4) AL=1.E-4
|
|
AL=AL/F00
|
|
AL=LOG10(AL)
|
|
DO IWL=1,NWL-1
|
|
IW0=IWL
|
|
IF(AL.LE.WLHYD(ILINE,IWL+1)) GO TO 40
|
|
END DO
|
|
40 IW1=IW0+1
|
|
PRFF=(PRF0(IW0)*(WLHYD(ILINE,IW1)-AL)+PRF0(IW1)*
|
|
* (AL-WLHYD(ILINE,IW0)))/
|
|
* (WLHYD(ILINE,IW1)-WLHYD(ILINE,IW0))
|
|
SG=EXP(PRFF*AL10)*FID
|
|
SG=SG*WLI0**2*CINV/F00
|
|
PRF(IJ-IJ0+1)=SG*fac
|
|
END DO
|
|
end if
|
|
c
|
|
c XENOMORPH data for selected lines
|
|
c
|
|
else
|
|
ixn=ilxen(ii,jj)
|
|
nwl=nwlxen(ixn)
|
|
do iwl=1,nwl
|
|
call intxen(prfb0,prfr0,tl,anel,iwl,ixn,id)
|
|
prfb(iwl)=prfb0
|
|
prfr(iwl)=prfb0
|
|
end do
|
|
do ij=ij0,ij1
|
|
al=(freq(ij)-fr0(itr))/f00
|
|
if(abs(al).lt.1.e-4) al=1.e-4
|
|
all=log10(abs(al))
|
|
do iwl=1,nwl-1
|
|
iw0=iwl
|
|
if(all.le.alxen(ixn,iwl+1)) go to 50
|
|
end do
|
|
50 iw1=iw0+1
|
|
if(al.gt.0.) then
|
|
prff=(prfb(iw0)*(alxen(ixn,iw1)-all)+
|
|
* prfb(iw1)*(all-alxen(ixn,iw0)))/
|
|
* (alxen(ixn,iw1)-alxen(ixn,iw0))
|
|
else
|
|
prff=(prfr(iw0)*(alxen(ixn,iw1)-all)+
|
|
* prfr(iw1)*(all-alxen(ixn,iw0)))/
|
|
* (alxen(ixn,iw1)-alxen(ixn,iw0))
|
|
end if
|
|
sg=exp(prff*al10)*fid/f00
|
|
prf(ij-ij0+1)=sg
|
|
end do
|
|
END IF
|
|
c
|
|
c for IP ge 10 - special line profile
|
|
c
|
|
ELSE IF(IP.GE.10) THEN
|
|
DO IJ=IJ0,IJ1
|
|
PRF(IJ-IJ0+1)=PROFSP(FREQ(IJ),DOP,ITR,ID)
|
|
END DO
|
|
END IF
|
|
C
|
|
C if required, force the profile at the endpoints to zero
|
|
C
|
|
IF(INTM0.EQ.-1) THEN
|
|
PRF(IJ1-IJ0+1)=0.
|
|
ELSE IF(INTM0.EQ.-2) THEN
|
|
PRF(1)=0.
|
|
PRF(IJ1-IJ0+1)=0.
|
|
END IF
|
|
RETURN
|
|
END
|