subroutine allard(xl,t,hneutr,hcharg,prof,iq,jq) c ================================================ c c quasi-molecular opacity for Lyman alpha, beta, and Balmer alpha c modified routine provided originally by D. Koester c c Input: xl: wavelength in [A] c hneutr: neutral H particle density [cm-3] c hcharg: ionized H particle density [cm-3] c iq: quantum number of the lower level c jq: quantum number of the upper level; c =2 - Lyman alpha c =3 - Lyman beta c Output: prof: Lyman alpha line profile, normalized to 1.0e8 c if integrated over A; c It then renormalized by multiplying by c 8.853e-29*lambda_0^2*f_ij c INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' parameter (NXMAX=1400,NNMAX=5,NTAMAX=6) parameter (xnorma=8.8528e-29*1215.6*1215.6*0.41618, * xnormb=8.8528e-29*1025.73*1025.7*0.0791, * xnormg=8.8528e-29*972.53*972.53*0.0290, * xnormc=8.8528e-29*6562.*6562.*0.6407) common /callarda/xlalp(NXMAX),plalp(NXMAX,NNMAX),stnnea,stncha, * vneua,vchaa,nxalp,iwarna common /callardb/xlbet(NXMAX),plbet(NXMAX,NNMAX),stnneb,stnchb, * vneub,vchab,nxbet,iwarnb common /callardg/xlgam(NXMAX),plgam(NXMAX,NNMAX),stnneg,stnchg, * vneug,vchag,nxgam,iwarng common /callardc/xlbal(NXMAX),plbal(NXMAX,NNMAX),stnnec,stnchc, * vneuc,vchac,nxbal,iwarnc common /calphatd/xlalpd(NXMAX,NTAMAX),plalpd(NXMAX,NNMAX,NTAMAX), * stnead(ntamax),stnchd(ntamax), * vneuad(ntamax),vchaad(ntamax), * talpd(ntamax),nxalpd(ntamax),ntalpd common/quasun/tqmprf,iquasi,nunalp,nunbet,nungam,nunbal c prof=0. c c Lyman alpha c if(iq.eq.1.and.jq.eq.2) then if(nunalp.lt.0) then call allardt(xl,t,hneutr,hcharg,prof) else if(xl.lt.xlalp(1).or.xl.gt.xlalp(nxalp)) return vn1=hneutr/stnnea vn2=hcharg/stncha vns=vn1*vneua+vn2*vchaa if(iwarna.eq.0) then if(vn1*vneua.gt.0.3.or.vn2*vchaa.gt.0.3) then write(*,*) ' warning: density too high for', * ' Lyman alpha expansion' iwarna=1 endif endif vn11=vn1*vn1 vn22=vn2*vn2 vn12=vn1*vn2 xnorm=1.0/(1.0+vns+0.5*vns*vns) c jl=0 ju=nxalp+1 10 if(ju-jl.gt.1) then jm=(ju+jl)/2 if((xlalp(nxalp).gt.xlalp(1)).eqv.(xl.gt.xlalp(jm))) then jl=jm else ju=jm endif go to 10 endif j=jl c if(j.eq.0) j=1 if(j.eq.nxalp) j=j-1 a1=(xl-xlalp(j))/(xlalp(j+1)-xlalp(j)) p1= vn1*((1.0-a1)*plalp(j,1)+a1*plalp(j+1,1)) p11=vn11*((1.0-a1)*plalp(j,2)+a1*plalp(j+1,2)) p2= vn2*((1.0-a1)*plalp(j,3)+a1*plalp(j+1,3)) p22=vn22*((1.0-a1)*plalp(j,4)+a1*plalp(j+1,4)) p12=vn12*((1.0-a1)*plalp(j,5)+a1*plalp(j+1,5)) prof=(p1+p2+p11+p22+p12)*xnorm*xnorma return end if end if c c Lyman beta c if(iq.eq.1.and.jq.eq.3) then if(nxbet.eq.0) return if(xl.lt.xlbet(1).or.xl.gt.xlbet(nxbet)) return vn1=hneutr/stnneb vn2=hcharg/stnchb vns=vn1*vneub+vn2*vchab if(iwarnb.eq.0) then if(vn1*vneub.gt.0.3.or.vn2*vchab.gt.0.3) then write(*,*) ' warning: density too high for', * ' Lyman beta expansion' iwarnb=1 endif endif vn11=vn1*vn1 vn22=vn2*vn2 vn12=vn1*vn2 xnorm=1.0/(1.0+vns+0.5*vns*vns) c jl=0 ju=nxbet+1 20 if(ju-jl.gt.1) then jm=(ju+jl)/2 if((xlbet(nxbet).gt.xlbet(1)).eqv.(xl.gt.xlbet(jm))) then jl=jm else ju=jm endif go to 20 endif j=jl c if(j.eq.0) j=1 if(j.eq.nxbet) j=j-1 a1=(xl-xlbet(j))/(xlbet(j+1)-xlbet(j)) p1= vn1*((1.0-a1)*plbet(j,1)+a1*plbet(j+1,1)) p11=vn11*((1.0-a1)*plbet(j,2)+a1*plbet(j+1,2)) p2= vn2*((1.0-a1)*plbet(j,3)+a1*plbet(j+1,3)) p22=vn22*((1.0-a1)*plbet(j,4)+a1*plbet(j+1,4)) p12=vn12*((1.0-a1)*plbet(j,5)+a1*plbet(j+1,5)) prof=(p1+p2+p11+p22+p12)*xnorm*xnormb return end if c c Lyman gamma c if(iq.eq.1.and.jq.eq.4) then if(nxgam.eq.0) return if(xl.lt.xlgam(1).or.xl.gt.xlgam(nxgam)) return vn1=hneutr/stnneg vn2=hcharg/stnchg vns=vn1*vneug+vn2*vchag if(iwarng.eq.0) then if(vn1*vneug.gt.0.3.or.vn2*vchag.gt.0.3) then write(*,*) ' warning: density too high for', * ' Lyman gamma expansion' iwarng=1 endif endif vn11=vn1*vn1 vn22=vn2*vn2 vn12=vn1*vn2 xnorm=1.0/(1.0+vns+0.5*vns*vns) c jl=0 ju=nxgam+1 30 if(ju-jl.gt.1) then jm=(ju+jl)/2 if((xlgam(nxgam).gt.xlgam(1)).eqv.(xl.gt.xlgam(jm))) then jl=jm else ju=jm endif go to 30 endif j=jl c if(j.eq.0) j=1 if(j.eq.nxgam) j=j-1 a1=(xl-xlgam(j))/(xlgam(j+1)-xlgam(j)) p1= vn1*((1.0-a1)*plgam(j,1)+a1*plgam(j+1,1)) p11=vn11*((1.0-a1)*plgam(j,2)+a1*plgam(j+1,2)) p2= vn2*((1.0-a1)*plgam(j,3)+a1*plgam(j+1,3)) p22=vn22*((1.0-a1)*plgam(j,4)+a1*plgam(j+1,4)) p12=vn12*((1.0-a1)*plgam(j,5)+a1*plgam(j+1,5)) prof=(p1+p2+p11+p22+p12)*xnorm*xnormg return end if c c Balmer alpha c if(iq.eq.2.and.jq.eq.3) then if(xl.lt.xlbal(1).or.xl.gt.xlbal(nxbal)) return vn1=0. vn2=hcharg/stnchc vns=vn1*vneuc+vn2*vchac vn11=vn1*vn1 vn22=vn2*vn2 vn12=vn1*vn2 xnorm=1.0/(1.0+vns+0.5*vns*vns) c jl=0 ju=nxbal+1 40 if(ju-jl.gt.1) then jm=(ju+jl)/2 if((xlbal(nxbal).gt.xlbal(1)).eqv.(xl.gt.xlbal(jm))) then jl=jm else ju=jm endif go to 40 endif j=jl c if(j.eq.0) j=1 if(j.eq.nxbal) j=j-1 a1=(xl-xlbal(j))/(xlbal(j+1)-xlbal(j)) p1= vn1*((1.0-a1)*plbal(j,1)+a1*plbal(j+1,1)) p11=vn11*((1.0-a1)*plbal(j,2)+a1*plbal(j+1,2)) p2= vn2*((1.0-a1)*plbal(j,3)+a1*plbal(j+1,3)) p22=vn22*((1.0-a1)*plbal(j,4)+a1*plbal(j+1,4)) p12=vn12*((1.0-a1)*plbal(j,5)+a1*plbal(j+1,5)) prof=(p1+p2+p11+p22+p12)*xnorm*xnormc end if c return end