88 lines
2.6 KiB
Fortran
88 lines
2.6 KiB
Fortran
subroutine rtesol(dtau,st0,rup,rdown,amu0,ri,ali)
|
|
C ================================================
|
|
C
|
|
C formal solver of the radiative transfer equation
|
|
C for one frequency, angle, and for completely known source function;
|
|
c by the Discontinuous Finite Element method
|
|
c Castor, Dykema, Klein, 1992, ApJ 387, 561.
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
parameter (one=1.d0)
|
|
dimension dtau(mdepth),st0(mdepth),ri(mdepth),ali(mdepth),
|
|
* rim(mdepth),rip(mdepth),aim(mdepth),aip(mdepth)
|
|
c
|
|
c incoming intensity
|
|
c
|
|
if(amu0.lt.0) then
|
|
c
|
|
id=1
|
|
rip(id)=rup
|
|
dt0=dtau(id)
|
|
dtaup1=dt0+one
|
|
dtau2=dt0*dt0
|
|
bb=two*dtaup1
|
|
cc=dt0*dtaup1
|
|
aa=dtau2+bb
|
|
rim(id)=(aa*rip(id)-cc*st0(id)+dt0*st0(id+1))/bb
|
|
do id=1,nd-1
|
|
dt0=dtau(id)
|
|
dtaup1=dt0+one
|
|
dtau2=dt0*dt0
|
|
bb=two*dtaup1
|
|
cc=dt0*dtaup1
|
|
aa=dtau2+bb
|
|
rim(id+1)=(two*rim(id)+dt0*st0(id)+cc*st0(id+1))/aa
|
|
rip(id)=(bb*rim(id)+cc*st0(id)-dt0*st0(id+1))/aa
|
|
aim(id+1)=cc/aa
|
|
aip(id)=(cc+bb*aim(id))/aa
|
|
enddo
|
|
do id=2,nd-1
|
|
dtt=un/(dtau(id-1)+dtau(id))
|
|
ri(id)=(rim(id)*dtau(id)+rip(id)*dtau(id-1))*dtt
|
|
ali(id)=(aim(id)*dtau(id)+aip(id)*dtau(id-1))*dtt
|
|
enddo
|
|
ri(1)=rip(1)
|
|
ri(nd)=rim(nd)
|
|
ali(1)=aim(1)
|
|
ali(nd)=aim(nd)
|
|
C
|
|
c outgoing intensity
|
|
c
|
|
else
|
|
c
|
|
rip(nd)=rdown
|
|
id=nd-1
|
|
dt0=dtau(id)
|
|
dtaup1=dt0+one
|
|
dtau2=dt0*dt0
|
|
bb=two*dtaup1
|
|
cc=dt0*dtaup1
|
|
aa=dtau2+bb
|
|
rim(id+1)=(aa*rip(id+1)-cc*st0(id+1)+dt0*st0(id))/bb
|
|
do id=nd-1,1,-1
|
|
dt0=dtau(id)
|
|
dtaup1=dt0+one
|
|
dtau2=dt0*dt0
|
|
bb=two*dtaup1
|
|
cc=dt0*dtaup1
|
|
aa=dtau2+bb
|
|
rim(id)=(two*rim(id+1)+dt0*st0(id+1)+cc*st0(id))/aa
|
|
rip(id+1)=(bb*rim(id+1)+cc*st0(id+1)-dt0*st0(id))/aa
|
|
aim(id)=cc/aa
|
|
aip(id+1)=(cc+bb*aim(id+1))/aa
|
|
enddo
|
|
do id=2,nd-1
|
|
dtt=un/(dtau(id-1)+dtau(id))
|
|
ri(id)=(rim(id)*dtau(id-1)+rip(id)*dtau(id))*dtt
|
|
ali(id)=(aim(id)*dtau(id-1)+aip(id)*dtau(id))*dtt
|
|
enddo
|
|
ri(1)=rim(1)
|
|
ri(nd)=rip(nd)
|
|
ali(1)=aim(1)
|
|
ali(nd)=aim(nd)
|
|
end if
|
|
c
|
|
return
|
|
end
|