SUBROUTINE RTEWIN(IU) C ===================== C C Solution of the radiative transfer equation - frequency by C frequency - for the known source function. C C The numerical method used: c Discontinuous Finite Element (DFE) method c Castor, Dykema, Klein, 1992, ApJ 387, 561. C C Input through blank COMMON block: C AB - two-dimensional array absorption coefficient (frequency, C depth) C STH - Thermal source function C C Version including velocity field and extension C radiative transfer along ray IU C INCLUDE 'PARAMS.FOR' INCLUDE 'MODELP.FOR' INCLUDE 'SYNTHP.FOR' INCLUDE 'WINCOM.FOR' PARAMETER (UN=1., TWO=2., HALF=0.5) PARAMETER (TAUREF = 0.6666666666667) DIMENSION ST0(2*MDEPF ),TAU(2*MDEPF ),AB0(2*MDEPF ), * rip(2*MDEPF ),rim(2*MDEPF ) c dimension sc0(2*mdepf) dimension sctd(2*mdepf) COMMON/COPAC/AB(MOPAC,MDEPF),STH(MOPAC,MDEPF),SCH(MFREQC,MDEPF) COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC) COMMON/CONSCV/SCCF(MFREQC,mdepf) COMMON/REFDEP/IREFD(MFREQ) C IUD=NUDF(IU) IF(IU.LE.NREXT) IUD=2*NUDF(IU)-1 IF(IUD.EQ.1) RETURN IF(NFREQ.GT.1) dlama0=(wlobs(nfrobs)-wlobs(1))/(nfrobs-1) C C overall loop over frequencies (observer's frame) C DO 500 IJ=1,NFROBS FR=FRQOBS(IJ) wl0=wlobs(ij) C C Opacity and total source function c interpolation in opacity table C IVK=NOPAC-2 DO ID=1,IUD KY=KRAY(IU,ID) YDR=DRAY(IU,ID) YDR1=UN-YDR dwlcom=wl0*DFRQF(IU,ID) wlcom=wl0+dwlcom if(wlcom.le.wlam(3)) then abd1=ab(1,ky-1) std1=sth(1,ky-1) abd0=ab(1,ky) std0=sth(1,ky) ij1=1 else if(wlcom.ge.wlam(nfreq)) then abd1=ab(nfreq,ky-1) std1=sth(nfreq,ky-1) abd0=ab(nfreq,ky) std0=sth(nfreq,ky) ij1=nfreq else xijap=(wlcom-wlam(3))/dlama0 ijap=int(xijap) ijap=max(ijap,1) ijap=min(ijap,nfreq) wlap=wlam(ijap) if(wlcom.lt.wlap) then ij1=ijap-1 do iji=ijap-1,1,-1 if(wlcom.ge.wlam(iji)) go to 20 end do 20 continue ij1=iji else ij1=ijap+1 do iji=ijap+1,nfreq if(wlcom.lt.wlam(iji)) go to 30 end do 30 continue ij1=iji-1 end if xfa=(wlam(ij1+1)-wlcom)/(wlam(ij1+1)-wlam(ij1)) abd1=xfa*ab(ij1,ky-1)+(1.-xfa)*ab(ij1+1,ky-1) std1=xfa*sth(ij1,ky-1)+(1.-xfa)*sth(ij1+1,ky-1) abd0=xfa*ab(ij1,ky)+(1.-xfa)*ab(ij1+1,ky) std0=xfa*sth(ij1,ky)+(1.-xfa)*sth(ij1+1,ky) end if AB0(ID)=YDR1*Abd1+YDR*abd0 ST0(ID)=YDR1*Std1+YDR*Std0 C C Add scattering C IJC=IJCINT(IJ1) IF(IFREQ.NE.17) THEN SC1=YDR1*SCCF(ijc,KY-1)+YDR*SCCF(ijc,KY) SC2=YDR1*SCCF(ijc+1,KY-1)+YDR*SCCF(ijc+1,KY) SCT=FRX1(ij1)*SC1+(1.-FRX1(ij1))*SC2 sctd(id)=sct/ab0(id) ST0(ID)=ST0(ID)+SCT/AB0(ID) END IF ENDDO C C Optical depth scale C TAU(1)=0. IREF=1 IF(IU.LE.NFIRY) THEN DO ID=1,IUD-1 JD=ID IF(ID.GT.NUDF(IU)) JD=2*NUDF(IU)-ID-1 DT=HALF*(AB0(ID+1)+AB0(ID))*DELZF(IU,JD) TAU(ID+1)=TAU(ID)+DT END DO ELSE DO ID=1,IUD-1 JD=ID IF(ID.GT.NUD(IU)) JD=2*NUD(IU)-ID-1 DT=HALF*(AB0(ID+1)+AB0(ID))*DELZ(IU,JD) TAU(ID+1)=TAU(ID)+DT END DO END IF if(iu.eq.kmu) then DO ID=1,IUD-1 IF(TAU(ID).LE.TAUREF.AND.TAU(ID+1).GT.TAUREF) IREF=ID END DO irefd(ij)=iref end if C C Outgoing intensity C IF(IU.LE.NREXT) THEN C C 1. External rays C ndt=iud rip(ndt)=0. dt0=tau(ndt)-tau(ndt-1) dtaup1=dt0+un dtau2=dt0*dt0 bb=two*dtaup1 cc=dt0*dtaup1 aa=dtau2+bb rim(ndt)=(aa*rip(ndt)-cc*st0(ndt)+dt0*st0(ndt-1))/bb do id=1,iud-1 jd=iud-id dt0=tau(jd+1)-tau(jd) dtaup1=dt0+un dtau2=dt0*dt0 bb=two*dtaup1 cc=dt0*dtaup1 aa=un/(dtau2+bb) rim(jd)=(two*rim(jd+1)+dt0*st0(jd+1)+cc*st0(jd))*aa enddo ELSE C C 2. core rays C NDT=IUD FR15=FR*1.D-15 BNU=BN*FR15*FR15*FR15 PLAND=BNU/(EXP(HK*FR/TEMP(ND))-UN) DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-UN) DPLAN=(PLAND-DPLAN)/(TAU(IUD)-TAU(IUD-1)) RIP(NDT)=PLAND+DPLAN dt0=tau(ndt)-tau(ndt-1) dtaup1=dt0+un dtau2=dt0*dt0 bb=two*dtaup1 cc=dt0*dtaup1 aa=dtau2+bb rim(ndt)=(aa*rip(ndt)-cc*st0(ndt)+dt0*st0(ndt-1))/bb do id=iud-1,1,-1 dt0=tau(id+1)-tau(id) dtaup1=dt0+un dtau2=dt0*dt0 bb=two*dtaup1 cc=dt0*dtaup1 aa=un/(dtau2+bb) rim(id)=(two*rim(id+1)+dt0*st0(id+1)+cc*st0(id))*aa enddo ENDIF FLUX(IJ)=FLUX(IJ)+WMUH(IU)*RIM(1) c c if(ij.eq.1.or.ij.eq.3.or.ij.eq.5.or.ij.eq.9.or.ij.eq.83) then c if(iu.eq.2.or.iu.eq.20.or.iu.eq.60.or.iu.eq.80) then c do id=1,iud c write(79,679) ij,iu,id,ab0(id),st0(id),sctd(id), c * tau(id),rim(id), c * flux(ij) c end do c end if c end if c 679 format(3i5,1p6e12.4) C c CFX=WMUH(IU)*RIM(1) c write(78,780) ij,iu,wlobs(ij),cfx,RIM(1) c 780 format(2i4,f10.3,1p2e16.8) C c if(iflux.ge.1) then C C output of emergent specific intensities to Unit 10 (line points) C or 18 (two continuum points) C c IF(IJ.GT.2) THEN c WRITE(10,618) WLAM(IJ),FLUX(IJ),RIM(1),IU c ELSE c WRITE(18,618) WLAM(IJ),FLUX(IJ),RIM(1),IU c END IF c end if c 618 FORMAT(1H ,f10.3,2pe15.5,i5) C C if needed (if iprin.ge.3), output of interesting physical C quantities at the monochromatic optical depth tau(nu)=2/3 C c IF(IPRIN.GE.3) THEN c T0=LOG(TAU(IREF+1)/TAU(IREF)) c X0=LOG(TAU(IREF+1)/TAUREF)/T0 c X1=LOG(TAUREF/TAU(IREF))/T0 c DMREF=EXP(LOG(DM(IREF))*X0+LOG(DM(IREF+1))*X1) c TREF=EXP(LOG(TEMP(IREF))*X0+LOG(TEMP(IREF+1))*X1) c STREF=EXP(LOG(ST0(IREF))*X0+LOG(ST0(IREF+1))*X1) c SSREF=EXP(LOG(-SS0(IREF))*X0+LOG(-SS0(IREF+1))*X1) c SREF=STREF+SSREF c ALM=2.997925E18/FREQ(IJ) c WRITE(36,636) IJ,ALM,IREF,DMREF,TREF,STREF,SSREF,SREF c 636 FORMAT(1H ,I3,F10.3,I4,1PE10.3,0PF10.1,1X,1P3E10.3) c END IF C C Contribution to J and H C c do id=1,nud(iu) c rad1(id)=rad1(id)+wmuj(iu,id)*uf(id) c ali1(id)=ali1(id)+wmuj(iu,id)*af(id) c end do c FLUXc(IJ)=FLUXc(IJ)+WMUH(IU)*RIM(1) C C C end of the loop over frequencies C 500 CONTINUE RETURN END