SUBROUTINE RYBMAT(IJ) C ===================== C C evaluation of the complete-linearization matrices in the Rybicki c formalism c INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ALIPAR.FOR' INCLUDE 'ARRAY1.FOR' COMMON/RYBMTX/RA(MDEPTH),RB(MDEPTH),RC(MDEPTH),VR(MDEPTH), * UA(MDEPTH),UB(MDEPTH),UC(MDEPTH), * VA(MDEPTH),VB(MDEPTH),VC(MDEPTH),WR(MDEPTH), * WM(MDEPTH,MDEPTH) common/dsctva/dsct1(mdepth),dscn1(mdepth) PARAMETER (SIXTH=UN/6.D0, * THIRD=UN/3.D0, * TWOTHR=TWO/3.D0) C IJT=IJFR(IJ) C C ============================================================== C 1. components corresponding to the radiative transfer equation C ============================================================== C C ----------------------------------- C ID = 1 - upper boundary condition C ID=1 DDM=(DM(ID+1)-DM(ID))*HALF DTM=UN/((ABSO1(ID)+ABSO1(ID+1))*DDM) DTM2=DTM*DTM ALF=DTM*DDM FD=TWO*FH(IJ) EXTI=EXTRAD(IJ) BET=(EXTI-FD*RAD1(ID))*DTM GAM=(FAK1(ID)*RAD1(ID)-FAK1(ID+1)*RAD1(ID+1))*TWO*DTM2 S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID) C1=ALF*(TWO*GAM-BET) B1=C1-S0/ABSO1(ID) UNQ=UN+TWO*DTM*Q0(IJ) c unq=un RB(ID)=-(UN+DTM*(FD+TWO*FAK1(ID)*DTM))+ * SCAT1(ID)/ABSO1(ID)*UNQ RC(ID)= TWO*FAK1(ID+1)*DTM2 VR(ID)= GAM-BET+RAD1(ID)-S0*UNQ UB(ID)= B1*DABT1(ID)+(DEMT1(ID)+ * DSCT1(ID)*RAD1(ID))/ABSO1(ID)*UNQ- * emis1(id)/abso1(id)**2*dabt1(id)*two*dtm*q0(ij)+ * dtm*s0*(uu0(ij)*dm(1)*dabt1(id)- * two*q0(ij)*dtm*ddm*dabt1(id)) UC(ID)= C1*DABT1(ID+1)-two*dtm*q0(ij)*s0*dtm*ddm*dabt1(id+1) c if(iubc.gt.0) then corf=half/dtm rb(id)=rb(id)*corf rc(id)=rc(id)*corf vr(id)=vr(id)*corf c1=(gam-rad1(id)+s0)*corf*alf b1=c1-corf*s0/abso1(id) UB(ID)= B1*DABT1(ID)+(DEMT1(ID)+ * DSCT1(ID)*RAD1(ID))/ABSO1(ID)*corf UC(ID)= C1*DABT1(ID+1) end if C C C ---------------------------------- C 1 < ID < ND - normal depth point C DO ID=2,ND-1 DDM=(DM(ID)-DM(ID-1))*HALF DDP=(DM(ID+1)-DM(ID))*HALF DZM=ABSO1(ID)+ABSO1(ID-1) DZP=ABSO1(ID)+ABSO1(ID+1) DTAUP=DZP*DDP DTAUM=DZM*DDM DTAU0=HALF *(DTAUP+DTAUM) FRD=FAK1(ID)*RAD1(ID) ALF1=(FRD-FAK1(ID+1)*RAD1(ID+1))/DTAUP/DTAU0 GAM1=(FRD-FAK1(ID-1)*RAD1(ID-1))/DTAUM/DTAU0 BET1=ALF1+GAM1 X1=HALF *BET1/DTAU0 A1=(GAM1+X1*DTAUM)/DZM C1=(ALF1+X1*DTAUP)/DZP B1=A1+C1 BS=UN CHIELM=SCAT1(ID-1) CHIEL0=SCAT1(ID) CHIELP=SCAT1(ID+1) S0=(EMIS1(ID)+CHIEL0*RAD1(ID))/ABSO1(ID) AS=0. CS=0. A2=0. C2=0. A3=0. C3=0. BET2=0. SM=0. SP=0. IF(MOD(ISPLIN,3).EQ.0) GO TO 60 SM=(EMIS1(ID-1)+RAD1(ID-1)*CHIELM)/ABSO1(ID-1) SP=(EMIS1(ID+1)+RAD1(ID+1)*CHIELP)/ABSO1(ID+1) IF(ISPLIN.EQ.1) THEN C C spline collocation (ISPLIN=1) C AS=DTAUM/DTAU0*SIXTH CS=DTAUP/DTAU0*SIXTH BS=0.666666666666667D0 ALF2=AS*(RAD1(ID-1)-SM) GAM2=CS*(RAD1(ID+1)-SP) BET2=ALF2+GAM2 X =HALF *BET2/DTAU0 A2=(GAM2-X*DTAUM)/DZM C2=(ALF2-X*DTAUP)/DZP ELSE C C Hermitian method (ISPLIN=2) C AS=DTAUP*DTAUP/DTAUM/DTAU0 CS=DTAUM*DTAUM/DTAUP/DTAU0 AL3=(RAD1(ID+1)-SP-RAD1(ID)+S0)*SIXTH GA3=(RAD1(ID-1)-SM-RAD1(ID)+S0)*SIXTH AV=AL3*CS CV=GA3*AS AS=(UN-HALF *AS)*SIXTH CS=(UN-HALF *CS)*SIXTH BS=UN-AS-CS X=(AV+CV)/DTAU0/4.D0 A2=(X*DTAUM+HALF *CV-AV)/DZM C2=(X*DTAUP+HALF *AV-CV)/DZP BET2=AS*(RAD1(ID-1)-SM)+CS*(RAD1(ID+1)-SP) END IF C C auxiliary quantities C B1=B1-(A2+C2) A1=A1-A2 C1=C1-C2 A2=AS/ABSO1(ID-1) C2=CS/ABSO1(ID+1) A3=A2*SM C3=C2*SP C 60 CONTINUE B2=BS/ABSO1(ID) B3=B2*S0 A1=A1-A3 B1=B1-B3 C1=C1-C3 C C *** elements of the matrices C RA(ID)= FAK1(ID-1)/DTAUM/DTAU0-AS*(UN-CHIELM/ABSO1(ID-1)) RB(ID)=-FAK1(ID)/DTAU0*(UN/DTAUP+UN/DTAUM)- * BS*(UN-CHIEL0/ABSO1(ID)) RC(ID)= FAK1(ID+1)/DTAUP/DTAU0-CS*(UN-CHIELP/ABSO1(ID+1)) VR(ID)= BET1+BET2+BS*(RAD1(ID)-S0) UA(ID)= A1*DABT1(ID-1)+A2*(DEMT1(ID-1)+DSCT1(ID-1)*RAD1(ID-1)) UB(ID)= B1*DABT1(ID)+B2*(DEMT1(ID)+DSCT1(ID)*RAD1(ID)) UC(ID)= C1*DABT1(ID+1)+C2*(DEMT1(ID+1)+DSCT1(ID+1)*RAD1(ID+1)) END DO C C ---------------------------------- C ID=ND - lower boundary condition C ID=ND DDM=HALF*(DM(ID)-DM(ID-1)) T0=TEMP(ID) TM=TEMP(ID-1) IF(IBC.GT.0.AND.IBC.LT.4.AND.IDISK.EQ.0) THEN DTM=UN/((ABSO1(ID-1)+ABSO1(ID))*DDM) DTM2=DTM*DTM FD=TWO*FHD(IJT) FR=FREQ(IJT) FR15=FR*1.D-15 BNU=BN*FR15*FR15*FR15 X0=HK*FR/T0 XM=HK*FR/TM PLAND=BNU/(EXP(X0)-UN) PLANM=BNU/(EXP(XM)-UN) DPLDT0=PLAND/(UN-EXP(-X0))*X0/T0 DPLDTM=PLANM/(UN-EXP(-XM))*XM/TM DPLAN=(PLAND-PLANM)*DTM ALF=DTM*DDM BET=(PLAND-FD*RAD1(ID))*DTM GAM=(FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)- * THIRD*(PLAND-PLANM))*TWO*DTM2 S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID) A1=ALF*(TWO*GAM-BET) B1=A1-S0/ABSO1(ID) RA(ID)= TWO*FAK1(ID-1)*DTM2 RB(ID)=-(UN+DTM*(FD+TWO*FAK1(ID)*DTM))+ * SCAT1(ID)/ABSO1(ID) VR(ID)= GAM-BET+RAD1(ID)-S0 UA(ID)= B1*DABT1(ID-1) * +(DEMT1(ID-1)+ * DSCT1(ID-1)*RAD1(ID-1))/ABSO1(ID-1) * -DPLDTM*DTM2*TWOTHR UB(ID)= B1*DABT1(ID)+(DEMT1(ID)+ * DSCT1(ID)*RAD1(ID))/ABSO1(ID)+ * DPLDT0*DTM*(UN+TWOTHR*DTM) C if(ifryb.gt.0) then DTM=UN/((ABSO1(ID-1)+ABSO1(ID))*DDM) FR=FREQ(IJT) FR15=FR*1.D-15 BNU=BN*FR15*FR15*FR15 X0=HK*FR/T0 XM=HK*FR/TM PLAND=BNU/(EXP(X0)-UN) PLANM=BNU/(EXP(XM)-UN) DPLDT0=PLAND/(UN-EXP(-X0))*X0/T0 DPLDTM=PLANM/(UN-EXP(-XM))*XM/TM GAM=(FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)- * THIRD*(PLAND-PLANM))*DTM S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID) BS=HALF/DTM BET=BS*(RAD1(ID)-S0) A1=(GAM-BET)*DTM*DDM B1=A1-BS*S0/ABSO1(ID) RA(ID)=FAK1(ID-1)*DTM RB(ID)=-FAK1(ID)*DTM-BS*(UN-SCAT1(ID)/ABSO1(ID))-FHD(IJT) VR(ID)=GAM+BET-HALF*PLAND+FHD(IJT)*RAD1(ID) UA(ID)=A1*DABT1(ID-1)-DPLDTM*DTM*THIRD UB(ID)=B1*DABT1(ID)+ * BS*(DEMT1(ID)+DSCT1(ID)*RAD1(ID))/ABSO1(ID)+ * (HALF+THIRD*DTM)*DPLDT0 end if C ELSE C C for disks - C lower b.c. expresses just I(taumax,-mu,nu)=I(taumax,+mu,nu) C DZM=ABSO1(ID)+ABSO1(ID-1) FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1) DTAUM=DZM*DDM GAM1=FRD/DTAUM A1=GAM1/DZM AS=0. BS=DTAUM*HALF S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID) GAM2=BS*(RAD1(ID)-S0) X1=GAM2/DZM A1=A1-X1 B2=BS/ABSO1(ID) B1=A1-B2*S0 RA(ID)=FAK1(ID-1)/DTAUM-AS*(UN-SCAT1(ID-1)/ABSO1(ID-1)) RB(ID)=-FAK1(ID)/DTAUM-BS*(UN-SCAT1(ID)/ABSO1(ID)) UA(ID)=A1*DABT1(ID-1)+A2*(DEMT1(ID-1)+DSCT1(ID-1)*RAD1(ID-1)) UB(ID)=B1*DABT1(ID)+B2*(DEMT1(ID)+DSCT1(ID)*RAD1(ID)) VR(ID)=GAM1+GAM2 END IF C C C ===================================================== C components corresponding to the radiative equilibrium C ===================================================== C DO ID=1,ND C C ********* integral equation part C IF(REINT(ID).GT.0.) THEN HEAT = ABSO1(ID)-SCAT1(ID) DHEAT = (DABT1(ID)-DSCT1(ID))*RAD1(ID) WDR = W(IJ)*dens(id)*reint(id) VB(ID) = HEAT*WDR WM(ID,ID)= WM(ID,ID)+(DHEAT-DEMT1(ID))*WDR WR(ID) = WR(ID)-(HEAT*RAD1(ID)-EMIS1(ID))*WDR END IF END DO C C ********* differential equation part C ID=1 IF(REDIF(ID).GT.0.) THEN WF=W(IJ)*FH(IJT)*REDIF(ID) VB(ID)=VB(ID)+WF WR(ID)=WR(ID)-WF*RAD1(ID) END IF c c original variant with 1st-order mid-point differences c if(icentr.eq.0) then nd1=nd if(ilbc.ne.0) nd1=nd-1 DO ID=2,nd1 IF(REDIF(ID).GT.0.) THEN DDM=(DM(ID)-DM(ID-1))*HALF OMEG0=ABSO1(ID) OMEGM=ABSO1(ID-1) DTAUM=(OMEG0+OMEGM)*DDM FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1) GAMR=FRD/DTAUM*redif(id) A1=GAMR/(OMEG0+OMEGM) VA(ID) =-W(IJ)*FAK1(ID-1)/DTAUM*redif(id) VB(ID) = VB(ID)+W(IJ)*FAK1(ID)/DTAUM*redif(id) WM(ID,ID-1)= WM(ID,ID-1)-A1*W(IJ)*DABT1(ID-1) WM(ID,ID) = WM(ID,ID)-A1*W(IJ)*DABT1(ID) WR(ID) = WR(ID)-W(IJ)*GAMR END IF END DO c if(ilbc.gt.0) then id=nd IF(REDIF(ID).GT.0.) THEN DDM=(DM(ID)-DM(ID-1))*HALF DTAUM=(ABSO1(ID)+ABSO1(ID-1))*DDM*3.D0 FR=FREQ(IJT) FR15=FR*1.D-15 BNU=BN*FR15*FR15*FR15 X0=HK*FR/TEMP(ID) XM=HK*FR/TEMP(ID-1) PLAND=BNU/(EXP(X0)-UN) PLANM=BNU/(EXP(XM)-UN) DPLDT0=PLAND/(UN-EXP(-X0))*X0/TEMP(ID) DPLDTM=PLANM/(UN-EXP(-XM))*XM/TEMP(ID-1) FLX=(PLAND-PLANM)/DTAUM*redif(id) A1=FLX/(ABSO1(ID)+ABSO1(ID-1)) WM(ID,ID-1)= WM(ID,ID-1)+W(IJ)*(DPLDTM-A1*DABT1(ID-1)) WM(ID,ID) = WM(ID,ID)+W(IJ)*(DPLDT0-A1*DABT1(ID)) WR(ID) = WR(ID)-W(IJ)*FLX END IF end if c c centered difference variant c else do id=2,nd-1 if(redif(id).gt.0) then wwr=w(ij)*redif(id) ddm=half*(dm(id)-dm(id-1)) ddp=half*(dm(id+1)-dm(id)) dtm=(abso1(id-1)+abso1(id))*ddm dtp=(abso1(id+1)+abso1(id))*ddp dt0=dtm+dtp frm=fak1(id)*rad1(id)-fak1(id-1)*rad1(id-1) frp=fak1(id+1)*rad1(id+1)-fak1(id)*rad1(id) alp=dtp/dtm/dt0 gam=dtm/dtp/dt0 flx=alp*frm+gam*frp c c matrix elements c va(id)=-wwr*fak1(id-1)*alp vb(id)=vb(id)+wwr*fak1(id)*(alp-gam) vc(id)=wwr*fak1(id+1)*gam c dmtm=ddm/(dtm*dtm) dmtp=ddp/(dtp*dtp) rm=dtm/dt0 rp=dtp/dt0 wm(id,id-1) = wm(id,id-1)+wwr*dabt1(id-1)*dmtm* * (rm*rm*frp-(un+rm)*frm) wm(id,id) = wm(id,id)+wwr*dabt1(id)* * ((dmtp*rp*rp-dmtm*rp*(un+rm))*frm + * (dmtm*rm*rm-dmtp*rm*(un+rp)*frp)) wr(id)=wr(id)-wwr*flx end if end do c id=nd IF(REDIF(ID).GT.0.) THEN DDM=(DM(ID)-DM(ID-1))*HALF DTAUM=(ABSO1(ID)+ABSO1(ID-1))*DDM FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1) FLX=FRD/DTAUM*redif(id) A1=FLX/(ABSO1(ID)+ABSO1(ID-1)) VA(ID) =-W(IJ)*FAK1(ID-1)/DTAUM*redif(id) VB(ID) = VB(ID)+W(IJ)*FAK1(ID)/DTAUM*redif(id) WM(ID,ID-1)= WM(ID,ID-1)-A1*W(IJ)*DABT1(ID-1) WM(ID,ID) = WM(ID,ID)-A1*W(IJ)*DABT1(ID) WR(ID) = WR(ID)-W(IJ)*FLX END IF end if c if(nretc.gt.0) then do id=1,nretc wr(id)=0. vb(id)=0. va(id)=0. wm(id,id)=1. if(id.gt.1) wm(id,id-1)=0. end do end if RETURN END