108 lines
3.3 KiB
Fortran
108 lines
3.3 KiB
Fortran
SUBROUTINE COMPT0(IJ,ID,ab,compa,compb,compc,compe,comps,compd)
|
|
C ===============================================================
|
|
C
|
|
c auxiliary quantities for the Compton scattering source function
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
|
|
common/auxcbc/cden1m(mdepth),cden10(mdepth),
|
|
* cden2m(mdepth),cden20(mdepth)
|
|
c
|
|
IJI=NFREQ-KIJ(IJ)+1
|
|
if(iji.eq.1) then
|
|
compa=0.
|
|
compb=0.
|
|
compc=0.
|
|
compd=0.
|
|
compe=0.
|
|
comps=0.
|
|
return
|
|
end if
|
|
c
|
|
FR=FREQ(IJ)
|
|
frp=freq(ijorig(iji+1))
|
|
frm=freq(ijorig(iji-1))
|
|
xcomp=fr*xcon
|
|
e2=ycon*temp(id)
|
|
e1=xcomp-3.*e2
|
|
c
|
|
del0=two/(dlnfr(iji)+dlnfr(iji-1))
|
|
cder1p(iji)=(un-delj(iji,id))*del0
|
|
cder1m(iji)=-delj(iji-1,id)*del0
|
|
cder10(iji)=-del0*(un-delj(iji-1,id)-delj(iji,id))
|
|
ss0=elec(id)*sige/ab
|
|
if(ichcoo.eq.0) then
|
|
cder10(iji)=-cder1m(iji)-cder1p(iji)
|
|
compa=ss0*(e1*cder1m(iji)+e2*cder2m(iji))
|
|
compb=ss0*(un-xcomp-sigec(ij)/sige+e1*cder10(iji)+
|
|
* e2*cder20(iji))
|
|
compc=ss0*(e1*cder1p(iji)+e2*cder2p(iji))
|
|
else
|
|
epsnu=(ab-elec(id)*sigec(ij))/ab
|
|
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
|
|
compa=ss0*(-delj(iji-1,id)*zxxm12+e2*cder2m(iji))
|
|
compc=ss0*((un-delj(iji,id))*zxxp12+e2*cder2p(iji))
|
|
compb=ss0*(delj(iji,id)*zxxp12-(un-delj(iji-1,id))*zxxm12+
|
|
* e2*cder20(iji)-sigec(ij)/sige)-epsnu+1.
|
|
compe=0.
|
|
end if
|
|
compd=(-3.*cder10(iji)+cder20(iji))*rad(iji,id)
|
|
c
|
|
IF(ICOMDE.EQ.0) THEN
|
|
COMPA=0.
|
|
COMPC=0.
|
|
COMPB=0.
|
|
END IF
|
|
c
|
|
x0=ss0*bnus(iji)
|
|
if(icomst.eq.0) x0=0.
|
|
if(ichcoo.eq.0) then
|
|
compe=x0*(cder10(iji)-un)
|
|
compu=x0*cder1m(iji)
|
|
compv=x0*cder1p(iji)
|
|
cbs=compe*rad(iji,id)
|
|
compe=cbs
|
|
end if
|
|
comps=compb*rad(iji,id)
|
|
if(iji.gt.1) then
|
|
if(ichcoo.eq.0) cbs=cbs+compu*rad(iji-1,id)
|
|
comps=comps+compa*rad(iji-1,id)
|
|
compd=compd+(-3.*cder1m(iji)+cder2m(iji))*rad(iji-1,id)
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
if(ichcoo.eq.0) cbs=cbs+compv*rad(iji+1,id)
|
|
comps=comps+compc*rad(iji+1,id)
|
|
compd=compd+(-3.*cder1p(iji)+cder2p(iji))*rad(iji+1,id)
|
|
end if
|
|
if(ichcoo.eq.0) then
|
|
compb=compb+cbs
|
|
compa=compa+compu*rad(iji,id)
|
|
compc=compc+compv*rad(iji,id)
|
|
comps=comps+cbs*rad(iji,id)
|
|
end if
|
|
compd=compd*ss0*ycon
|
|
IF(ICOMDE.EQ.0) COMPD=0.
|
|
c
|
|
c a variant with ICOMPT=2 - no off-diagonal terms in intensity
|
|
c
|
|
if(icompt.eq.2) then
|
|
if(iji.gt.1) compb=compb+compa*rad(iji-1,id)
|
|
if(iji.lt.nfreq) compb=compb+compc*rad(iji+1,id)
|
|
compa=0.
|
|
compc=0.
|
|
else if(icompt.eq.3) then
|
|
compa=0.
|
|
compb=0.
|
|
compc=0.
|
|
end if
|
|
return
|
|
end
|