172 lines
5.1 KiB
Fortran
172 lines
5.1 KiB
Fortran
SUBROUTINE RTECF0(IJ)
|
|
C =====================
|
|
C
|
|
C Setup of the individual matrix elements of matrices A,B,C, E,U,V,
|
|
C and alpha, beta, gamma, for solving the coupled transfer equation
|
|
C with Compton scattering
|
|
C Evaluation for a given frequency point IJ.
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
COMMON/OPTDPT/DT(MDEPTH)
|
|
PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
|
|
COMMON/AUXRTE/
|
|
* COMA(MDEPTH),COMB(MDEPTH),COMC(MDEPTH),VL(MDEPTH),
|
|
* COME(MDEPTH),U(MDEPTH),V(MDEPTH),BS(MDEPTH),
|
|
* AL(MDEPTH),BE(MDEPTH),GA(MDEPTH)
|
|
common/auxcbc/cden1m(mdepth),cden10(mdepth),
|
|
* cden2m(mdepth),cden20(mdepth)
|
|
c
|
|
IJI=NFREQ-KIJ(IJ)+1
|
|
FR=FREQ(IJ)
|
|
frp=freq(ijorig(iji+1))
|
|
frm=freq(ijorig(iji-1))
|
|
xcomp=fr*xcon
|
|
C
|
|
C optical depth scale
|
|
C
|
|
do id=1,nd-1
|
|
dt(id)=deldmz(id)*(absot(id+1)+absot(id))
|
|
end do
|
|
C
|
|
C depth discretization matrices
|
|
C
|
|
c 1. upper boundary
|
|
c
|
|
id=1
|
|
dtp1=dt(id)
|
|
bb0=un/dtp1
|
|
bb1=two*bb0*bb0
|
|
be(id)=bb0*two*fh(ij)+bb1*fak(ij,id)
|
|
ga(id)=bb1*fak(ij,id+1)
|
|
sext=two*bb0*hextrd(ij)
|
|
c
|
|
c 2. normal depth point
|
|
c
|
|
do id=2,nd-1
|
|
dtm1=dtp1
|
|
dtp1=dt(id)
|
|
dt0=two/(dtm1+dtp1)
|
|
al(id)=fak(ij,id-1)/dtm1*dt0
|
|
ga(id)=fak(ij,id+1)/dtp1*dt0
|
|
be(id)=fak(ij,id)*dt0*(un/dtm1+un/dtp1)
|
|
end do
|
|
c
|
|
c 3. lower boundary
|
|
c
|
|
id=nd
|
|
c
|
|
c stellar atmospheric
|
|
c
|
|
IF(IDISK.EQ.0.OR.IFZ0.LT.0) then
|
|
IF(IBC.EQ.0) THEN
|
|
be(ID)=fak(ij,id)/DTP1+HALF
|
|
al(ID)=fak(ij,id-1)/DTP1
|
|
ELSE IF(IBC.LT.4) THEN
|
|
B=UN/DTP1
|
|
A=TWO*B*B
|
|
be(id)=B*TWO*fhd(ij)+a*fak(ij,id)
|
|
al(id)=a*fak(ij,id-1)
|
|
ELSE
|
|
B=UN/DTP1
|
|
A=TWO*B*B
|
|
be(id)=b+a*fak(ij,id)
|
|
al(id)=a*fak(ij,id-1)
|
|
END IF
|
|
c
|
|
c accretion disk - symmetric boundary
|
|
c
|
|
ELSE
|
|
bb0=un/dtp1
|
|
bb1=two*bb0*bb0
|
|
B=TWO/DTP1
|
|
be(id)=bb1*fak(ij,id)
|
|
al(id)=bb1*fak(ij,id-1)
|
|
END IF
|
|
C
|
|
C scattering matrices
|
|
C
|
|
do id=1,nd
|
|
scat0=elec(id)*sige
|
|
sa0=emis1(id)/abso1(id)
|
|
ss0=scat0/abso1(id)
|
|
epsnu=(abso1(id)-scat1(id))/abso1(id)
|
|
x0=ss0
|
|
e2=ycon*temp(id)+0.7*xcomp*xcomp
|
|
e1=xcomp-3.*e2-0.7*xcomp*xcomp
|
|
e0=1.-xcomp-4.2*xcomp*xcomp
|
|
coma(id)=0.
|
|
comc(id)=0.
|
|
u(id)=0.
|
|
v(id)=0.
|
|
vl(id)=sa0
|
|
if(id.eq.1.and.iwinbl.lt.0) vl(id)=sa0+sext
|
|
bs(id)=0.
|
|
if(iji.eq.1) then
|
|
comc(id)=0.
|
|
comb(id)=x0*(un-2.*xcomp)
|
|
else if(iji.lt.nfreq) then
|
|
del0=two/(dlnfr(iji)+dlnfr(iji-1))
|
|
cder1p(iji)=(un-delj(iji,id))*del0
|
|
cder1m(iji)=-delj(iji-1,id)*del0
|
|
if(ichcoo.eq.0) then
|
|
cder10(iji)=-cder1m(iji)-cder1p(iji)
|
|
coma(id)=x0*(e1*cder1m(iji)+e2*cder2m(iji))
|
|
comb(id)=x0*(e0+e1*cder10(iji)+e2*cder20(iji))
|
|
comc(id)=x0*(e1*cder1p(iji)+e2*cder2p(iji))
|
|
x0=ss0*bnus(iji)
|
|
IF(ICOMST.EQ.0) X0=0.
|
|
come(id)=x0*(cder10(iji)-un)
|
|
u(id)=x0*cder1m(iji)
|
|
v(id)=x0*cder1p(iji)
|
|
bs(id)=come(id)*rad(iji,id)+
|
|
* u(id)*rad(iji-1,id)+v(id)*rad(iji+1,id)
|
|
else
|
|
cder10(iji)=-del0*(un-delj(iji-1,id)-delj(iji,id))
|
|
zxxp=xcon*frp+0.5*bnus(iji+1)*rad(iji+1,id)-3.*e2
|
|
zxx0=xcomp+0.5*bnus(iji)*rad(iji,id)-3.*e2
|
|
zxxm=xcon*frm+0.5*bnus(iji-1)*rad(iji-1,id)-3.*e2
|
|
zxxp12=((un-delj(iji,id))*zxxp+delj(iji,id)*zxx0)*del0
|
|
zxxm12=((un-delj(iji-1,id))*zxx0+delj(iji-1,id)*zxxm)*
|
|
* del0
|
|
coma(id)=x0*(-delj(iji-1,id)*zxxm12+e2*cder2m(iji))
|
|
comc(id)=x0*((un-delj(iji,id))*zxxp12+e2*cder2p(iji))
|
|
comb(id)=x0*(delj(iji,id)*zxxp12-(un-delj(iji-1,id))*
|
|
* zxxm12+e2*cder20(iji))-epsnu+1
|
|
end if
|
|
else
|
|
dlt=delj(iji-1,id)
|
|
zj1=exp(-hk*freq(ij)/temp(id))
|
|
zj2=exp(-hk*freq(ij+1)/temp(id))
|
|
if(ichcoo.eq.0) then
|
|
zj0=un/(hk*sqrt(freq(ij)*freq(ij+1))/temp(id))
|
|
zxx=un-3.*zj0+(un-dlt)*zj1+dlt*zj2
|
|
comb(id)=zj0/dlnfr(iji-1)+(un-dlt)*zxx
|
|
coma(id)=-zj0/dlnfr(iji-1)+dlt*zxx
|
|
else
|
|
zxx0=xcomp*(un+zj1)-3.*e2
|
|
zxxm=xcon*frm*(un+zj2)-3.*e2
|
|
zxx=(un-dlt)*zxx0+dlt*zxxm
|
|
comb(id)=e2/dlnfr(iji-1)+(un-dlt)*zxx
|
|
coma(id)=-e2/dlnfr(iji-1)+dlt*zxx
|
|
end if
|
|
vl(id)=0.
|
|
if(icomde.ne.0) then
|
|
al(id)=0.
|
|
be(id)=-un
|
|
ga(id)=0.
|
|
end if
|
|
end if
|
|
if(icomde.eq.0) then
|
|
coma(id)=0.
|
|
comc(id)=0.
|
|
comb(id)=x0*(un-2.*xcomp)
|
|
end if
|
|
end do
|
|
C
|
|
return
|
|
end
|